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