(!*******************************************************
   Mosel Example Problems 
   ======================

   file benders_single.mos
   ```````````````````````
   Benders decomposition for solving a simple MIP.
   - Single model containing several subproblems - 

   (c) 2009 Fair Isaac Corporation
       author: S. Heipcke, Jan. 2009, rev. Jul. 2020
*******************************************************!)
        
model "Benders (single model)"
 uses "mmxprs", "mmjobs"

 parameters
  NCTVAR = 3  ! 1
  NINTVAR = 3 ! 2
  NC = 4      ! 3
  BIGM = 1000
  ALG = 1                             ! 1: Use Benders decomposition (dual)
                                      ! 2: Use Benders decomposition (primal)
  DATAFILE = "bprob33.dat"  ! "bprob12.dat"
 end-parameters

 forward procedure start_solution
 forward procedure solve_primal_int(ct: integer)
 forward procedure solve_cont
 forward function eval_solution: boolean
 forward procedure print_solution

 forward procedure define_intprob(prob:mpproblem)
 forward function solve_intprob(prob:mpproblem, ct:integer): real
 forward procedure define_dualprob(prob:mpproblem)
 forward function solve_dualprob(prob:mpproblem, Alg:integer): real
 forward procedure define_contprob(prob:mpproblem)
 forward function solve_contprob(prob:mpproblem): real
 
 declarations
  STEP_0=2                            ! Codes sent to subproblems
  STEP_1=3
  STEP_2=4
  STAT_SOLVED=6                       ! Status codes returned by subproblems
  STAT_INFEAS=7
  STAT_READY=8

  CtVars = 1..NCTVAR                  ! Continuous variables
  IntVars = 1..NINTVAR                ! Discrete variables
  Ctrs = 1..NC                        ! Set of constraints (orig. problem)

  A: array(Ctrs,CtVars) of integer    ! Coeff.s of continuous variables
  B: array(Ctrs,IntVars) of integer   ! Coeff.s of discrete variables
  b: array(Ctrs) of integer           ! RHS values
  C: array(CtVars) of integer         ! Obj. coeff.s of continuous variables
  D: array(IntVars) of integer        ! Obj. coeff.s of discrete variables
  Ctr: array(Ctrs) of linctr          ! Constraints of orig. problem
  CtrD: array(CtVars) of linctr       ! Constraints of dual problem
  MC: array(range) of linctr          ! Constraints generated by alg.
  Obj: linctr                         ! Continuous objective

  y: array(IntVars) of mpvar          ! Discrete variables
  z: mpvar                            ! Objective function variable
  u: array(Ctrs) of mpvar             ! Dual variables
  x: array(CtVars) of mpvar           ! Continuous variables

  sol_u: array(Ctrs) of real          ! Solution of dual problem
  sol_x: array(CtVars) of real        ! Solution of primal prob. (cont.)
  sol_y: array(IntVars) of real       ! Solution of primal prob. (integers)
  sol_obj: real                       ! Objective function value (primal)

  RM: range                           ! Problem indices
  stepprob: dynamic array(RM) of mpproblem    ! Subproblems
  status: dynamic array(mpproblem) of integer ! Subproblem status
 end-declarations

 initializations from DATAFILE
  A B b C D
 end-initializations

! **** Subproblems ****
                                      ! Create and define submodels
 create(stepprob(1)); define_intprob(stepprob(1))
 
 if ALG=1 then
  create(stepprob(2)); define_dualprob(stepprob(2))
 else
  create(stepprob(0)); define_dualprob(stepprob(0))
  create(stepprob(2)); define_contprob(stepprob(2))
 end-if

! **** Solution algorithm ****

 start_solution                       ! 0. Initial solution for cont. var.s
 ct:= 1
 repeat
 writeln("\n**** Iteration: ", ct)
  solve_primal_int(ct)                ! 1. Solve problem with fixed cont.
  solve_cont                          ! 2. Solve problem with fixed int.
  ct+=1
 until eval_solution                  !    Test for optimality
 print_solution                       ! 3. Retrieve and display the solution

!-----------------------------------------------------------
! Produce an initial solution for the dual problem using a dummy objective
 procedure start_solution
  num:= if(ALG=1, 2, 0)
  res:=solve_dualprob(stepprob(num), STEP_0)  ! Start the problem solving

  if status(stepprob(num))=STAT_INFEAS then
   writeln("Problem is infeasible")
   exit(6) 
  end-if
 end-procedure

!-----------------------------------------------------------
! Solve a modified version of the primal problem, replacing continuous
! variables by the solution of the dual
 procedure solve_primal_int(ct: integer)  
  sol_obj:= solve_intprob(stepprob(1), ct)
 end-procedure

!-----------------------------------------------------------
! Solve the Step 2 problem (dual or primal depending on value of ALG)
! for given solution values of y
 procedure solve_cont
  if ALG=1 then                       ! Start the problem solving
   res:= solve_dualprob(stepprob(2), STEP_2)
  else 
   res:= solve_contprob(stepprob(2))
  end-if
 end-procedure

!-----------------------------------------------------------
 function eval_solution: boolean
  write("Test optimality: ", sol_obj - sum(i in IntVars) D(i)*sol_y(i),
          " >= ", sum(j in Ctrs) sol_u(j)* (b(j) - 
	                            sum(i in IntVars) B(j,i)*sol_y(i)) )
  returned:= ( sol_obj - sum(i in IntVars) D(i)*sol_y(i) >= 
      sum(j in Ctrs) sol_u(j)* (b(j) - sum(i in IntVars) B(j,i)*sol_y(i)) )
  writeln(if(returned, " : true", " : false"))
 end-function

!-----------------------------------------------------------
 procedure print_solution  
  ! Reset all submodels
  forall(m in RM) reset(stepprob(m))
   
  write("\n**** Solution (Benders): ", sol_obj, "\n  x: ")
  forall(i in CtVars) write(sol_x(i), " ")
  write("  y: ")
  forall(i in IntVars) write(sol_y(i), " ")
  writeln
 end-procedure

!-----------------------------------------------------------
! Define the integer problem
 procedure define_intprob(prob:mpproblem)
  with prob do
   z is_free                            ! Objective is a free variable
   forall(i in IntVars) y(i) is_integer ! Integrality condition
   forall(i in IntVars) y(i) <= BIGM    ! Set (large) upper bound
  end-do
  status(prob):= STAT_READY
 end-procedure 
 
! Solve the integer problem
 function solve_intprob(prob:mpproblem, ct:integer): real
  with prob do
   status(prob):= STAT_READY
    
   ! Add a new constraint
   MC(ct):= z >= sum(i in IntVars) D(i)*y(i) +  
                 sum(j in Ctrs) sol_u(j)*(b(j) - sum(i in IntVars) B(j,i)*y(i))
    
   minimize(z)
  
   ! Store solution values of y
   forall(i in IntVars) sol_y(i):= getsol(y(i))
  
   returned:=getobjval
   status(prob):= STAT_SOLVED

   write("Step 1: ", getobjval, "\n  y: ")
   forall(i in IntVars) write(sol_y(i), " ")  
   write("\n  Slack: ")
   forall(j in 1..ct) write(getslack(MC(j)), " ")
   writeln
   fflush

  end-do
 end-function
 
!-----------------------------------------------------------
! Define the dual problem
 procedure define_dualprob(prob:mpproblem)
  with prob do
   forall(i in CtVars) CtrD(i):= sum(j in Ctrs) u(j)*A(j,i) <= C(i)
  end-do
  status(prob):= STAT_READY
 end-procedure 

! Process dual solution data
 procedure save_dualsolution(prob:mpproblem)
 ! There is no need for 'with mpproblem do' here since this subroutine
 ! is called from within a subproblem
 
 ! Store values of u and x
  forall(j in Ctrs) sol_u(j):= getsol(u(j))
  forall(i in CtVars) sol_x(i):= getdual(CtrD(i))
  
  status(prob):= STAT_SOLVED

  write(getobjval, "\n  u: ")
  forall(j in Ctrs) write(sol_u(j), " ")
  write("\n  x: ")
  forall(i in CtVars) write(getdual(CtrD(i)), " ")
  writeln
  fflush
 end-procedure

! (Re)solve the dual problem
 function solve_dualprob(prob:mpproblem, Alg:integer): real
  with prob do
   status(prob):= STAT_READY

   if Alg=STEP_0 then                 ! Produce an initial solution for the
                                      ! dual problem using a dummy objective
    maximize(sum(j in Ctrs) u(j))
    if(getprobstat = XPRS_INF) then 
     writeln("Problem is infeasible")
     status(prob):= STAT_INFEAS       ! Problem is infeasible
    else
     write("**** Start solution: ")
     save_dualsolution(prob)
     returned:= getobjval
    end-if

   else                               ! STEP 2: Solve the dual problem for
                                      ! given solution values of y
    Obj:= sum(j in Ctrs) u(j)* (b(j) - sum(i in IntVars) B(j,i)*sol_y(i))
    maximize(XPRS_PRI, Obj)

    if(getprobstat=XPRS_UNB) then
     BigM:= sum(j in Ctrs) u(j) <= BIGM
     maximize(XPRS_PRI, Obj)
    end-if
   
    write("Step 2: ")
    save_dualsolution(prob)           ! Save solution values
    returned:= getobjval
    BigM:= 0                          ! Reset the 'BigM' constraint
   end-if

  end-do
 end-function

!-----------------------------------------------------------
! Define the continuous primal problem
 procedure define_contprob(prob:mpproblem)
  with prob do
   Obj:= sum(i in CtVars) C(i)*x(i)
  end-do
  status(prob):= STAT_READY
 end-procedure 
 
! Solve the continuous problem
 function solve_contprob(prob:mpproblem): real
  with prob do
   status(prob):= STAT_READY
 
  ! (Re)define and solve this problem
   forall(j in Ctrs)
    Ctr(j):= sum(i in CtVars) A(j,i)*x(i) + 
             sum(i in IntVars) B(j,i)*sol_y(i) >= b(j) 

   minimize(Obj)                       ! Solve the problem

  ! Store values of u and x
   forall(j in Ctrs) sol_u(j):= getdual(Ctr(j))
   forall(i in CtVars) sol_x(i):= getsol(x(i))

   returned:=getobjval
   status(prob):= STAT_SOLVED

   write("Step 2: ", getobjval, "\n  u: ")
   forall(j in Ctrs) write(sol_u(j), " ")
   write("\n  x: ")
   forall(i in CtVars) write(getsol(x(i)), " ")
   writeln
   fflush
  
  end-do
 end-function

!-----------------------------------------------------------

end-model