Alternative implementation with multiple problems
Since our decomposition algorithm is formed by a series of sequential optimization runs it is possible to implement the whole approach as a single model defining several problems (instead of several model files, each with a single problem). This alternative implementation requires significantly less effort for data handling (problems simply share data) and coordination of solving (problem solving in a single model always is sequential). As a result, the total Mosel model code is shorter. However, by collecting all problems into a single model, the testing of (sub)problems in isolation is made somewhat more complicated—in the previous version we can simply run one of the (sub)models as a stand-alone model, a feature that certainly is helpful when developing large applications.
Master model
In the master model, all code related to (sub)model handling is replaced by the somewhat simpler handling of (sub)problems that are defined and solved within the same model as the master problem. Here we just display the part of the declarations that are new or modified, followed by the setup of the subproblems and the (unchanged) solution algorithm with the definition of the `start solve' subroutines.
For every submodel we now have two new subroutines, define_*prob and solve_*prob that replace the separate model files. These subroutines are lsited in the following sections.
declarations
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
RM: range ! Problem indices
stepprob: array(RM) of mpproblem ! Subproblems
status: array(mpproblem) of integer ! Subproblem status
end-declarations
! **** 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 Submodel 1: fixed continuous variables
The bounds and integrality conditions for the integer subproblem are stated once (define_intprob). Each time the integer subproblem is solved (solve_intprob) a new constraint gets added. The solving routine also stores the solution values and the problem status into global data structures.
! 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 Submodel 2: fixed integer variables
The setup routine for the continuous problem (define_contprob) only defines the objective function. The constraints are re-defined each time the problem gets solved (solve_contprob).
! 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 Submodel 0: start solution
The setup routine for the dual problem (define_dualprob) states the constraints for this model. The solving subroutines solve_dualprob implements two cases, depending on whether we use this model with a dummy objective for finding a start solution (Step 0) or in the second part of the solution algorithm where the objective is defined with the current solution values (Step 2).
! 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
!-----------------------------------------------------------
! (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 Similarly to the previous implementation as a separate model, we have moved the storing of the solution values into a subroutine (save_dualsolution). Notice that it is not necessary to switch problems using with prob do in this subroutine because it gets called from within a subproblem.
! Process dual solution data
procedure save_dualsolution(prob:mpproblem)
! 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