Benders decomposition: working with 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).

Consider the following standard form of a mixed-integer programming problem (problem I).

minimize Cx + Dy
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

C_i·x_i + ∑

NINTVAR

i=1

D_i·y_i, where NCTVAR denotes the number of continuous variables and NINTVAR the number of integer variables.

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:

minimize Cx
Ax ≥ b - By'
x ≥ 0

The dual program of problem II is given by problem IId:

maximize u(b - By')
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, 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 u^p (at a vertex p of the feasible region) we have u^p(b - By') ≤ Cx^*.

To obtain the optimum solution of the original MIP problem (problem I) we may use the following partitioning algorithm that is known as Benders decomposition algorithm:

Step 1

Solve the pure integer program:

minimize Dy + z
y ≥ 0, y integer

This integer problem will accumulate Benders cuts over time and it is always a relaxation of the original problem. The artificial variable z is introduced to take the contribution of x to the objective into account.

Step 2

With the solution y' of Step 1, solve the linear program:

minimize Cx
Ax ≥ b - By'
x ≥ 0

If this subproblem is infeasible, then the corresponding dual (note that we can ignore the objective Cx in this case):

maximize u(b - By')
uA ≤ 0
u ≥ 0

is unbounded and any extreme ray u^* gives a violated feasibility Benders cut to be added to the integer problem in Step 1 of the next iteration:

u^*(b - By') ≤ 0

If the subproblem has a finite optimum, then there are two possible cases:

1. The objective of the subproblem is ≤ z^* : the point (y^*, z^*) is optimal and we can recover an optimal solution for the original problem by concatenating the optimal solution of the subproblem x^* with y^*. The algorithm stops here.
2. The objective of the subproblem is > z^*: from the optimal multipliers u^* of the dual we can derive a violated optimality Benders cut to be added to the integer problem in Step 1 of the next iteration:
   u^*(b - By') ≤ z
   

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 of the original problem can be obtained.

A small example problem

Our implementation of Benders decomposition will solve the following example problem with NINTVAR=4 binary variables x_i, NCTVAR=5 continuous variables y_i, and NC=5 inequality constraints.

minimize

x1

+

6·x2

+

5·x3

+

7·x4

+

9·y1

+

3·y2

+

2·y4

+

3·y5

-5·x1

-

4·x2

+

3·x3

-

2·y1

+

3·y2

-

y3

+

2·y4

-

4·y5

≤5

-3·x2

+

4·x3

-

x4

-

5·y1

+

3·y2

-

2·y3

-

3·y4

≤4

x1

+

4·x3

-

2·x4

+

2·y1

+

y3

-

3·y4

-

5·y5

≤1

-5·x1

-

3·x3

+

x4

+

2·y2

+

3·y3

-

2·y4

+

y5

≤-4

-2·x2

+

3·x3

+

2·x4

+

3·y1

+

4·y2

+

2·y3

-

5·y5

≤-3

x binary, y ≥ 0

Implementation

Base model

The base model (file benders_decomp.mos) reads in the data (file bprob455.dat), defines the solution algorithm, coordinates the communication with the submodels, and prints out the solution at the end. Prior to starting the decomposition, the original problem is solved to optimality without decomposing, for validation purposes. Then, step 1 of the algorithm is performed by submodel 1, which solves the integer (main) problem and accumulates Benders cuts until optimality is achieved. For step 2 of the algorithm (solving the subproblem with fixed integer variables) we solve the problem using the primal version and then retrieve the dual information to derive cuts using the corresponding Mosel functions, however solving the dual subproblem directly would also be a valid option.

The implementation of the base model looks as follows. Prior to the start of the solution algorithm itself, the original problem is solved and 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 (base model)"
 uses "mmxprs", "mmjobs", "mmsystem"

 parameters
  NINTVAR = 4               ! Number of integer variables
  NCTVAR = 5                ! Number of continuous variables
  NC = 5                    ! Number of constraints
  DATAFILE = "bprob455.dat"  ! Data file with problem data
 end-parameters

 forward procedure print_solution
 forward procedure solve_original

 declarations
  ! Event codes exhanged with submodels
  EVENT_READY=2
  EVENT_SOLVED=3
  EVENT_INFEAS=4
  EVENT_ADDFEASCUT=5
  EVENT_ADDOPTCUT=6
  EVENT_MAINSOLVED=7

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

  A1: array(Ctrs,IntVars) of integer   ! Coeff.s of integer variables
  A2: array(Ctrs,CtVars) of integer    ! Coeff.s of cont variables
  b: array(Ctrs) of integer            ! RHS values
  c1: array(IntVars) of integer        ! Obj. coeff.s of integer variables
  c2: array(CtVars) of integer         ! Obj. coeff.s of cont variables

  sol_x: array(IntVars) of real        ! Solution of integer vars
  sol_y: array(CtVars) of real         ! Solution of continuous vars
  sol_obj: real                        ! Objective function value

  RM: range                            ! Model indices
  stepmod: dynamic array(RM) of Model  ! Submodels

  x: array(IntVars) of mpvar           ! Integer variables of original problem
  y: array(CtVars) of mpvar            ! Continuous variables of original problem
 end-declarations

 initializations from DATAFILE
  A1 A2 b c1 c2
 end-initializations

 ! Save problem data to be shared with submodels
 initializations to "bin:shmem:probdata"
  A1 A2 b c1 c2
 end-initializations

 ! Solve original problem without decomposing
 solve_original

 ! Compile + load submodels
 if compile("benders_int.mos")<>0: exit(1)
 create(stepmod(1)); load(stepmod(1), "benders_int.bim")
 if compile("benders_cont.mos")<>0: exit(2)
 create(stepmod(2)); load(stepmod(2), "benders_cont.bim")
                                      ! 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

The base model triggers the different submodels via the send function and retrieves the solution once the corresponding submodel solving has finished.

 ! **** Benders decomposition - Solution algorithm ****
 writeln("\n**** Benders decomposition algorithm ****")
 status:= EVENT_READY                       ! Initialize status flag
 ct:= 1                                     ! Iteration counter
 repeat
 writeln("\n**** Iteration ", ct, " ****")
  send(stepmod(1), status, ct)              ! 1. Solve main problem (with int. vars)
  wait                                      ! Wait for next event
  ev:=getnextevent                          ! Get the event data
  status:= getclass(ev)                     ! Store status of main problem
  if status = EVENT_MAINSOLVED then
   send(stepmod(2), status, ct)             ! 2. Solve subproblem (with cont. vars)
   wait                                     ! Wait for next event
   ev:= getnextevent                        ! Get the event data
   status:= getclass(ev)                    ! Store status of main problem
  else
   writeln("Error occurred with solving the main problem")
   exit(4)
  end-if

  ct+=1
 until status = EVENT_SOLVED

The base model prints out the final solution (procedure print_solution) once the decomposition algorithm ends, stops all submodels and cleans up temporary files:

 if status = EVENT_SOLVED then
  print_solution                           ! Retrieve and display the solution
 end-if

 ! **** Cleaning up temporary files ****
 forall(m in RM) stop(stepmod(m))           ! Stop all submodels
 fdelete("benders_int.bim")
 fdelete("benders_cont.bim")
 fdelete("shmem:probdata")
 fdelete("shmem:sol")

!-----------------------------------------------------------
 procedure print_solution
  ! Retrieve results
  initializations from "bin:shmem:sol"
   sol_x sol_y
  end-initializations

  sol_obj := sum(i in IntVars) c1(i)*sol_x(i) + sum(i in CtVars) c2(i)*sol_y(i)
  writeln("\n**** Final solution (Benders) ****")
  writeln("\nSolution with Benders decomposition: obj: ", sol_obj, ", x: ",
    sol_x, ", y: ", sol_y)
 end-procedure

Lastly, the definition of the function solve_original, which solves the original problem without decomposing prior to the decomposition algorithm starts:

!-----------------------------------------------------------
 ! Solve the original problem without decomposing
 procedure solve_original
  writeln("\n**** Solving original problem without decomposing ****")
  forall(i in IntVars) x(i) is_binary    ! Integrality condition vars x

  forall(j in Ctrs)
   Ctr(j):= sum(i in IntVars) A1(j,i) * x(i) + sum(i in CtVars) A2(j,i) * y(i) <= b(j)

  minimize(sum(i in IntVars) c1(i) * x(i) + sum(i in CtVars) c2(i) * y(i))

  ! Get and print solution
  forall(i in IntVars)
   sol_x(i):= getsol(x(i))
  forall(i in CtVars)
   sol_y(i):= getsol(y(i))

  writeln("\nSolution without decomposing: obj: ", getobjval, ", x: ",
    sol_x, ", y: ", sol_y)

 end-procedure
end-model

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 binary decision variables. It then halts at the wait statement (first line of the repeat-until loop) until the base model sends it a (solving) event. At each invocation of solving this problem, the solution values of a previous run of the continuous problem—read from memory—are used to define a new feasibility or optimality cut (constraint) to 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 base model. The complete source of this submodel (file benders_int.mos) is listed below.

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

 parameters
  NINTVAR = 4               ! Number of integer variables
  NCTVAR = 5                ! Number of continuous variables
  NC = 5                    ! Number of constraints
 end-parameters

 declarations
  ! Event codes exhanged with submodels
  EVENT_READY=2
  EVENT_SOLVED=3
  EVENT_INFEAS=4
  EVENT_ADDFEASCUT=5
  EVENT_ADDOPTCUT=6
  EVENT_MAINSOLVED=7

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

  A1: array(Ctrs,IntVars) of integer  ! Coeff.s of discrete variables
  A2: array(Ctrs,CtVars) of integer   ! Coeff.s of cont variables
  b: array(Ctrs) of integer           ! RHS values
  c1: array(IntVars) of integer       ! Obj. coeff.s of discrete variables
  c2: array(CtVars) of integer        ! Obj. coeff.s of cont variables

  sol_u: array(Ctrs) of real          ! Solution of dual subproblem for cuts (cont.)
  sol_x: array(IntVars) of real       ! Solution of main problem (integers)

  x: array(IntVars) of mpvar          ! Discrete variables
  z: mpvar                            ! Objective function variable
  obj_z: real                         ! Objective value Z to compare

 end-declarations

 initializations from "bin:shmem:probdata"
  A1 A2 b c1 c2
 end-initializations

 forall(i in IntVars) x(i) is_binary   ! Binary variables x

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

 repeat
  wait
  ev:= getnextevent
  status:= getclass(ev)
  ct:= getvalue(ev)                    ! Get the event code and value

  if status = EVENT_ADDFEASCUT then    ! If subproblem is infeasible, add feas. cut

   writeln("Adding feasibility cut to the integer problem...")

   initializations from "bin:shmem:sol"  ! Retrieve the subproblem solution data
    sol_u
   end-initializations

   NewFeasCut(ct) := sum(j in Ctrs) sol_u(j) * (b(j) - sum(i in IntVars) A1(j,i)*x(i)) <= 0

  elif status = EVENT_ADDOPTCUT then   ! If subproblem has finite optimum with obj < z*,
                                       ! add optimality cut
   writeln("Adding optimality cut to the integer problem...")

   initializations from "bin:shmem:sol"     ! Retrieve the subproblem solution data
    sol_u
   end-initializations

   NewOptCut(ct) := sum(j in Ctrs) sol_u(j) * (b(j) - sum(i in IntVars) A1(j,i)*x(i)) <= z
  end-if

  writeln("Solving integer problem...")
  minimize(sum(i in IntVars) c1(i) * x(i) + z)

  if getprobstat=XPRS_OPT then
   forall(i in IntVars) sol_x(i):= getsol(x(i))
   obj_z := getsol(z)

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

   writeln("Step 1 solution: obj: ", getobjval, ", x: ", sol_x, ", z: ", obj_z)

   send(EVENT_MAINSOLVED, getobjval)

  elif getprobstat=XPRS_INF then
   writeln("ERROR: main problem is infeasible")
   send(EVENT_INFEAS, 0)
   exit(4)
  elif getprobstat=XPRS_UNB then
   writeln("ERROR: main problem is unbounded")
   exit(4)
  else
   writeln("Main problem has unknown status: ", getprobstat)
  end-if

 until false

end-model

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 base model. At every solving iteration the constraints Ctr are redefined using the sol_x array read from memory. After solving the subproblem, and depending on its status and objective value, dual coefficients are written back to memory using the dual information (methods getdualray and getdual) to derive a feasibility (if infeasible) or optimality (if subproblem has finite optimum with obj > z^*) cut to be added to the integer problem in the next iteration. In case the subproblem has a finite optimum that matches z^*, the optimal solution for the original problem has been found and can be retrieved by concatenating x^* (from the integer problem) and y^* (from the subproblem). Below follows the source of this model (file benders_cont.mos).

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

 parameters
  NINTVAR = 4               ! Number of integer variables
  NCTVAR = 5                ! Number of continuous variables
  NC = 5                    ! Number of constraints
 end-parameters

 declarations
  ! Event codes exhanged with submodels
  EVENT_READY=2
  EVENT_SOLVED=3
  EVENT_INFEAS=4
  EVENT_ADDFEASCUT=5
  EVENT_ADDOPTCUT=6
  EVENT_MAINSOLVED=7

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

  A1: array(Ctrs,IntVars) of integer  ! Coeff.s of continuous variables
  A2: array(Ctrs,CtVars) of integer   ! Coeff.s of discrete variables
  b: array(Ctrs) of integer           ! RHS values
  c1: array(IntVars) of integer       ! Obj. coeff.s of continuous variables
  c2: array(CtVars) of integer        ! Obj. coeff.s of cont variables

  sol_x: array(IntVars) of real       ! Solution of main prob. (integers)
  sol_y: array(CtVars) of real        ! Solution of primal subprob. (cont.)
  sol_u: array(Ctrs) of real          ! Solution of dual subproblem for cuts (cont.)

  y: array(CtVars) of mpvar           ! Decision variables for continuous subproblem

  constr : set of linctr              ! Set of lin constr to associate with dual ray
  dual_ray : array(constr) of real    ! Array to store dual ray, if infeasible

  obj_z: real                         ! Objective value Z to compare
 end-declarations

 initializations from "bin:shmem:probdata"
  A1 A2 b c1 c2
 end-initializations

 Obj := sum(i in CtVars) c2(i) * y(i)

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

! STEP 2: Solve the dual problem for given solution values of x
 repeat
  wait
  ev:= getnextevent
  status:= getclass(ev)

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

  forall(j in Ctrs)
   Ctr(j):= sum(i in CtVars) A2(j,i) * y(i) <= b(j) - sum(i in IntVars) A1(j,i) * sol_x(i)

  writeln("\nSolving continuous problem (primal)...")

  minimize(Obj)

  ! If infeasible, get dual ray to add the feasibility cut
  if getprobstat=XPRS_INF then
   writeln("The primal subproblem is infeasible -> Add feasibility cut!")
   hasRay := getdualray(dual_ray)
   if hasRay then

    ! Extract the components of the dual ray for feasibility cut
    forall(j in Ctrs)
     sol_u(j):= dual_ray(Ctr(j))

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

    send(EVENT_ADDFEASCUT, 0)
    else
     writeln("ERROR: No dual ray was found")
     exit(4)
    end-if

  elif getprobstat=XPRS_OPT then   ! If problem has finite optimum, then there are two cases

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

    forall(i in CtVars) sol_y(i):= getsol(y(i)) ! Get values of original vars y
    writeln("Step 2 solution: obj: ", getobjval, ", sol_y: ", sol_y)

    if getobjval <= obj_z then      ! Case 1) the objective of the subproblem matches z* =>
                                    ! the point (y*,z*) is optimal, we can recover an optimal
                                    ! solution for the original problem by concatenating the
                                    ! optimal solution of the subproblem x* with y*.
     writeln("The primal subproblem has finite optimum with obj <= z*",
      " -> Optimal solution found!")

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

     send(EVENT_SOLVED, getobjval)

    else                            ! Case 2) the objective of the subproblem is > z*.
                                    ! Then from the optimal multipliers u* of the dual
                                    ! we can derive a violated Benders optimality cut:
     writeln("The primal subproblem has finite optimum with obj > z*",
      " -> Add optimality cut!")
     forall(j in Ctrs)
      sol_u(j):= getdual(Ctr(j))  ! Get values of dual vars u (from the primal subproblem)

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

     send(EVENT_ADDOPTCUT, 0)
    end-if
  else
   writeln("Dual subproblem has unknown status: ", getprobstat)
  end-if

 until false

end-model

Results

The optimal solution to our small test problem has the objective function value of 2.8. 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:

**** Solving original problem without decomposing ****
Solution without decomposing: obj: 2.8, x: [1,0,0,0], y: [0,0,0,0,0.6]
**** Benders decomposition algorithm ****
**** Iteration 1 ****
Solving integer problem...
Step 1 solution: obj: 0, x: [0,0,0,0], z: 0
Solving continuous problem (primal)...
Step 2 solution: obj: 6.4, sol_y: [0,0,0,2.3,0.6]
The primal subproblem has finite optimum with obj > z* -> Add optimality cut!
**** Iteration 2 ****
Adding optimality cut to the integer problem...
Solving integer problem...
Step 1 solution: obj: 2.4, x: [1,0,0,0], z: 1.4
Solving continuous problem (primal)...
Step 2 solution: obj: 1.8, sol_y: [0,0,0,0,0.6]
The primal subproblem has finite optimum with obj > z* -> Add optimality cut!
**** Iteration 3 ****
Adding optimality cut to the integer problem...
Solving integer problem...
Step 1 solution: obj: 2.8, x: [1,0,0,0], z: 1.8
Solving continuous problem (primal)...
Step 2 solution: obj: 1.8, sol_y: [0,0,0,0,0.6]
The primal subproblem has finite optimum with obj <= z* -> Optimal solution found!
**** Final solution (Benders) ****
Solution with Benders decomposition: obj: 2.8, x: [1,0,0,0], y: [0,0,0,0,0.6]