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-2025 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.