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

   file benders_main.mos
   `````````````````````
   Benders decomposition for solving a simple MIP.
   - Main model - 

   *** ATTENTION: This model will return an error if ***
   *** no more than one Xpress licence is available. ***

   (c) 2008 Fair Isaac Corporation
       author: S. Heipcke, Jan. 2006, rev. Sep. 2022
*******************************************************!)
        
model "Benders (main model)"
 uses "mmxprs", "mmjobs", "mmsystem"

 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

 declarations
  STEP_0=2                            ! Event codes sent to submodels
  STEP_1=3
  STEP_2=4
  EVENT_SOLVED=6                      ! Event codes sent by submodels
  EVENT_INFEAS=7
  EVENT_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.

  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                           ! Model indices
  stepmod: dynamic array(RM) of Model ! Submodels
 end-declarations

 initializations from DATAFILE
  A B b C D
 end-initializations

! **** Submodels ****

 initializations to "bin:shmem:probdata"   ! Save data for submodels
  A  B  b  C  D  
 end-initializations

! Compile + load all submodels
 if compile("benders_int.mos")<>0: exit(1)
 create(stepmod(1)); load(stepmod(1), "benders_int.bim")
 if compile("benders_dual.mos")<>0: exit(2)

 if ALG=1 then
  create(stepmod(2)); load(stepmod(2), "benders_dual.bim")
 else
  create(stepmod(0)); load(stepmod(0), "benders_dual.bim")
  if compile("benders_cont.mos")<>0: exit(3)
  create(stepmod(2)); load(stepmod(2), "benders_cont.bim")
  run(stepmod(0), "NCTVAR=" + NCTVAR + ",NINTVAR=" + NINTVAR + ",NC=" + NC)
 end-if
                                      ! Start the execution of the submodels
 run(stepmod(1), "NINTVAR=" + NINTVAR + ",NC=" + NC)
 run(stepmod(2), "NCTVAR=" + NCTVAR + ",NINTVAR=" + NINTVAR + ",NC=" + NC)

 forall(m in RM) do
  wait                                ! Wait for "Ready" messages
  ev:= getnextevent
  if getclass(ev) <> EVENT_READY then
   writeln("Error occurred in a subproblem")
   exit(4)
  end-if
 end-do

! **** 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


! **** Cleaning up temporary files ****
 fdelete("benders_int.bim")
 fdelete("benders_dual.bim")
 if ALG<>1: fdelete("benders_cont.bim")
 fdelete("shmem:probdata")
 fdelete("shmem:sol")
 
!-----------------------------------------------------------
! Produce an initial solution for the dual problem using a dummy objective
 procedure start_solution
  if ALG=1 then                       ! Start the problem solving
   send(stepmod(2), STEP_0, 0)
  else
   send(stepmod(0), STEP_0, 0)
  end-if
  wait                                ! Wait for the solution
  ev:=getnextevent
  if getclass(ev)=EVENT_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)  
  send(stepmod(1), STEP_1, ct)        ! Start the problem solving
  wait                                ! Wait for the solution
  ev:=getnextevent
  sol_obj:= getvalue(ev)              ! Store objective function value

  initializations from "bin:shmem:sol"  ! Retrieve the solution
   sol_y 
  end-initializations
 end-procedure

!-----------------------------------------------------------
! Solve the Step 2 problem (dual or primal depending on value of ALG)
! for given solution values of y
 procedure solve_cont
  send(stepmod(2), STEP_2, 0)         ! Start the problem solving
  wait                                ! Wait for the solution
  dropnextevent

  initializations from "bin:shmem:sol"  ! Retrieve the solution
   sol_u 
  end-initializations
 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  
  ! Retrieve results
  initializations from "bin:shmem:sol"
   sol_x 
  end-initializations

  forall(m in RM) stop(stepmod(m))     ! Stop all submodels
   
  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

end-model