Implementation

Master model

The master model reads in the data, defines the solution algorithm, coordinates the communication between the submodels, and prints out the solution at the end. For step 2 of the algorithm (solving the dual problem with fixed integer variables) we have the choice to solve either the primal problem and retrieve the dual solution values from the Optimizer or to define the dual problem ourselves and solve it. Model parameter ALG lets the user choose between these two options.

The main part of the master model looks as follows. Prior to the start of the solution algorithm itself all submodels are compiled, loaded, and started so that in each step of the algorithm we simply need to trigger the (re)solving of the corresponding submodel.

model "Benders (master model)"
 uses "mmxprs", "mmjobs"

 parameters
  NCTVAR = 3
  NINTVAR = 3
  NC = 4
  BIGM = 1000
  ALG = 1                             ! 1: Use Benders decomposition (dual)
                                      ! 2: Use Benders decomposition (primal)
  DATAFILE = "bprob33.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: 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 then exit(1); end-if
 create(stepmod(1)); load(stepmod(1), "benders_int.bim")
 if compile("benders_dual.mos")<>0 then exit(2); end-if

 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 then exit(3); end-if
  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

The subroutines starting the different submodels send a `start solving' event and retrieve the solution once the submodel solving has finished. For the generation of the start solution we need to choose the right submodel, according to the settings of the parameter ALG. If this problem is found to be infeasible, then the whole problem is infeasible and we stop the execution of the model.

! 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

The master model also tests whether the termination criterion is fulfilled (function eval_solution) and prints out the final solution (procedure print_solution). The latter procedure also stops all submodels:

 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 

Submodel 1: fixed continuous variables

In the first step of the decomposition algorithm we need to solve a pure integer problem. When the execution of this model is started it reads in the invariant data and sets up the variables. It then halts at the wait statement (first line of the repeat-until loop) until the master model sends it a (solving) event. At each invocation of solving this problem, the solution values of the previous run of the continuous problem—read from memory—are used to define a new constraint MC(k) for the integer problem. The whole model, and with it the endless loop into which the solving is embedded, will be terminated only by the `stop model' command from the master model. The complete source of this submodel (file benders_int.mos) is listed below.

model "Benders (integer problem)"
 uses "mmxprs", "mmjobs"

 parameters
  NINTVAR = 3
  NC = 4
  BIGM = 1000
 end-parameters

 declarations
  STEP_0=2                            ! Event codes sent to submodels
  STEP_1=3
  EVENT_SOLVED=6                      ! Event codes sent by submodels
  EVENT_READY=8

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

  B: array(Ctrs,IntVars) of integer   ! Coeff.s of discrete variables
  b: array(Ctrs) of integer           ! RHS values
  D: array(IntVars) of integer        ! Obj. coeff.s of discrete variables
  MC: array(range) of linctr          ! Constraints generated by alg.

  sol_u: array(Ctrs) of real          ! Solution of dual problem
  sol_y: array(IntVars) of real       ! Solution of primal prob.

  y: array(IntVars) of mpvar          ! Discrete variables
  z: mpvar                            ! Objective function variable
 end-declarations

 initializations from "bin:shmem:probdata"
  B  b  D
 end-initializations

 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

 send(EVENT_READY,0)                  ! Model is ready (= running)

 repeat
  wait
  ev:= getnextevent
  ct:= integer(getvalue(ev))

  initializations from "bin:shmem:sol"
   sol_u
  end-initializations

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

  initializations to "bin:shmem:sol"
   sol_y
  end-initializations

  send(EVENT_SOLVED, getobjval)

  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

 until false

end-model 

Since the problems we are solving differ only by a single constraint from one iteration to the next, it may be worthwhile to save the basis of the solution to the root LP-relaxation (not the basis to the MIP solution) and reload it for the next optimization run. However, for our small test case we did not observe any improvements in terms of execution speed. For saving and re-reading the basis, the call to minimize needs to be replaced by the following sequence of statements:

declarations
 bas: basis
end-declarations

loadprob(z)
loadbasis(bas)
minimize(XPRS_LPSTOP, z)
savebasis(bas)
minimize(XPRS_CONT, z)

Submodel 2: fixed integer variables

The second step of our decomposition algorithm consists in solving a subproblem where all integer variables are fixed to their solution values found in the first step. The structure of the model implementing this step is quite similar to the previous submodel. When the model is run, it reads the invariant data from memory and sets up the objective function. It then halts at the line wait at the beginning of the loop to wait for a step 2 solving event sent by the master model. At every solving iteration the constraints CTR are redefined using the coefficients read from memory and the solution is written back to memory. Below follows the source of this model (file benders_cont.mos).

model "Benders (continuous problem)"
 uses "mmxprs", "mmjobs"

 parameters
  NCTVAR = 3
  NINTVAR = 3
  NC = 4
  BIGM = 1000
 end-parameters

 declarations
  STEP_0=2                            ! Event codes sent to submodels
  STEP_2=4
  STEP_3=5
  EVENT_SOLVED=6                      ! Event codes sent by submodels
  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
  Ctr: array(Ctrs) of linctr          ! Constraints of orig. problem

  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)

  x: array(CtVars) of mpvar           ! Continuous variables
 end-declarations

 initializations from "bin:shmem:probdata"
  A  B  b  C
 end-initializations

 Obj:= sum(i in CtVars) C(i)*x(i)

 send(EVENT_READY,0)                  ! Model is ready (= running)

! (Re)solve this model until it is stopped by event "STEP_3"
 repeat
  wait
  dropnextevent

  initializations from "bin:shmem:sol"
   sol_y
  end-initializations

  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))

  initializations to "bin:shmem:sol"
   sol_u  sol_x
  end-initializations

  send(EVENT_SOLVED, getobjval)

  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

 until false

end-model 

Submodel 0: start solution

To start the decomposition algorithm we need to generate an initial set of values for the continuous variables. This can be done by solving the dual problem in the continuous variables with a dummy objective function. A second use of the dual problem is for Step 2 of the algorithm, replacing the primal model we have seen in the previous section. The implementation of this submodel takes into account these two cases: within the solving loop we test for the type (class) of event that has been sent by the master problem and choose the problem to be solved accordingly.

The main part of this model is implemented by the following Mosel code (file benders_dual.mos).

model "Benders (dual problem)"
 uses "mmxprs", "mmjobs"

 parameters
  NCTVAR = 3
  NINTVAR = 3
  NC = 4
  BIGM = 1000
 end-parameters

 forward procedure save_solution

 declarations
  STEP_0=2                            ! Event codes sent to submodels
  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

  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)

  u: array(Ctrs) of mpvar             ! Dual variables
 end-declarations

 initializations from "bin:shmem:probdata"
  A  B  b  C
 end-initializations

 forall(i in CtVars) CtrD(i):= sum(j in Ctrs) u(j)*A(j,i) <= C(i)

 send(EVENT_READY,0)                  ! Model is ready (= running)

! (Re)solve this model until it is stopped by event "STEP_3"
 repeat
  wait
  ev:= getnextevent
  Alg:= getclass(ev)

  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")
    send(EVENT_INFEAS,0)              ! Problem is infeasible
   else
    write_("**** Start solution: ")
    save_solution
   end-if

  else                                ! STEP 2: Solve the dual problem for
                                      ! given solution values of y
   initializations from "bin:shmem:sol"
    sol_y
   end-initializations

   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_solution                      ! Write solution to memory
   BigM:= 0                           ! Reset the 'BigM' constraint
  end-if

 until false 

This model is completed by the definition of the subroutine save_solution that writes the solution to memory and informs the master model of it being available by sending the EVENT_SOLVED message.

! Process solution data
 procedure save_solution
 ! 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))

  initializations to "bin:shmem:sol"
   sol_u  sol_x
  end-initializations

  send(EVENT_SOLVED, getobjval)

  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

end-model