Benders decomposition: working with several different submodels
Topics covered in this section:
Benders decomposition is a decomposition method for solving large Mixed Integer Programming problems. Instead of solving a MIP problem that may be too large for standard solution methods all-in-one, we work with a sequence of linear and pure integer subproblems (the latter with a smaller number of constraints than the original problem). The description of the decomposition algorithm below is taken from [Hu69] where the interested reader will also find proofs of its finiteness and of the statement that it always results in the optimal solution.
Consider the following standard form of a mixed-integer programming problem (problem I).
Ax + By ≥ b
x, y ≥ 0, y integer
In the above, and all through this section, we use bold letters for vectors and matrices. For instance, Cx + Dy stands for ∑
| NCTVAR |
| i=1 |
| NINTVAR |
| i=1 |
For given values y' of y the problem above is reduced to a linear program (problem II)—we leave out the constant term in the objective:
Ax ≥ b - By'
x ≥ 0
The dual program of problem II is given by problem IId:
uA ≤ C
u ≥ 0
An interesting feature of the dual problem IId is that the feasible region of u is independent of y. Furthermore, from duality theory it follows that if problem IId is infeasible or has no finite optimum solution, then the original problem I has no finite optimum solution. Again from duality theory we know that if problem IId has a finite optimum solution u* then this solution has the same value as the optimum solution to the primal problem (that is, u*(b - By') = Cx*), and for a solution up (at a vertex p of the feasible region) we have up(b - By') ≤ Cx*. Substituting the latter into the objective of the original MIP problem we obtain a constraint of the form
To obtain the optimum solution z* of the original MIP problem (problem I) we may use the following partitioning algorithm that is known as Benders decomposition algorithm:
- Step 0
-
Find a
u' that satisfies
u'A ≤
C
if no such u' exists
then the original problem (I) has no feasible solution
else continue with Step 1
end-if - Step 1
-
Solve the pure integer program
minimize zIf z is unbounded from below, take a z to be any small value z'.
z ≥ u'(b - By) + Dy
y ≥ 0, y integer - Step 2
-
With the solution
y' of Step 1, solve the linear program
maximize u(b - By')If u goes to inifinity with u(b - By')
uA ≤ C
u ≥ 0
then add the constraint ∑iui ≤ M, where M is a large positive constant, and resolve the problem
end-if
Let the solution of this program be u''.
if z' - Dy' ≥ u''(b - By')
then continue with Step 3
else return to Step 1 and add the constraint z' ≥ Dy + u''(b - By)
end-if - Step 3
-
With the solution
y' of Step 1, solve the linear program
minimize Cxx' and y' are the optimum solution z* = Cx' + Dy'
Ax ≥ b - By'
x ≥ 0
This algorithm is provably finite. It results in the optimum solution and at any time during its execution lower and upper bounds on the optimum solution z* can be obtained.
A small example problem
Our implementation of Benders decomposition will solve the following example problem with NCTVAR=3 continuous variables xi, NINTVAR=3 integer variables yi, and NC=4 inequality constraints.
| minimize | 5·x1 | + | 4·x2 | + | 3·x3 | + | 2·y1 | + | 2·y2 | + | 3·y3 | |
| x1 | - | x2 | + | 5·x3 | + | 3·y1 | + | 2·y3 | ≥ 5 | |||
| 4·x2 | + | 3·x3 | - | y1 | + | y2 | ≥ 7 | |||||
| 2·x1 | - | 2·x3 | + | y1 | - | y2 | ≥ 4 | |||||
| 3·x1 | + | 5·y1 | + | 5·y2 | + | 5·y3 | ≥ -2 | |||||
| x, y ≥ 0, y integer | ||||||||||||
Implementation
Main model
The main 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 implementation of the main 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 (main 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: 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 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 main 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 parent 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 parent 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 parent 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 parent 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 parent 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 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.
Main problem
In the (main) 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 main 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 listed 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 Subproblem 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 Subproblem 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 Subproblem 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 Results
The optimal solution to our small test problem has the objective function value 18.1852. Our program produces the following output, showing that the problem is solved to optimality with 3 iterations (looping around steps 1 and 2) of the decomposition algorithm:
**** Start solution: 4.05556 u: 0.740741 1.18519 2.12963 0 x: 0.611111 0.166667 0.111111 **** Iteration: 1 Step 1: -1146.15 y: 1000 0 0 Slack: 0 Step 2: 1007 u: 0 1 0 0 x: 0 251.75 0 Test optimality: -3146.15 = 1007 : false **** Iteration: 2 Step 1: 17.0185 y: 3 0 0 Slack: 0 -1.01852 Step 2: 12.5 u: 0 1 2.5 0 x: 0.5 2.5 0 Test optimality: 11.0185 = 12.5 : false **** Iteration: 3 Step 1: 18.1852 y: 2 0 0 Slack: 0 -5.18519 -0.185185 Step 2: 14.1852 u: 0.740741 1.18519 2.12963 0 x: 1.03704 2.22222 0.037037 Test optimality: 14.1852 = 14.1852 : true **** Solution (Benders): 18.1852 x: 1.03704 2.22222 0.037037 y: 2 0 0
© 2001-2023 Fair Isaac Corporation. All rights reserved. This documentation is the property of Fair Isaac Corporation (“FICO”). Receipt or possession of this documentation does not convey rights to disclose, reproduce, make derivative works, use, or allow others to use it except solely for internal evaluation purposes to determine whether to purchase a license to the software described in this documentation, or as otherwise set forth in a written software license agreement between you and FICO (or a FICO affiliate). Use of this documentation and the software described in it must conform strictly to the foregoing permitted uses, and no other use is permitted.