Benders Decomposition

🌐 日本語 | English

Overview

Benders decomposition is a method for solving optimization problems that mix binary and continuous variables (mixed-integer programming problems) by splitting them into two smaller problems that are solved alternately.

• Master problem: decides only the binary variables

• Subproblem: decides only the continuous variables, with the binary variables fixed to the master problem’s choice

Each time the subproblem is solved, a cut (a constraint added back to the master problem) is generated from its result. As cuts accumulate, the master problem approximates the original problem better and better, and the lower bound (a value the optimum cannot go below) and the upper bound (the objective value of the best solution found so far) close in from both sides. When they meet, the solution is optimal.

The method works especially well for problems where a small number of binary decisions controls many continuous variables (for example: deciding which facilities to open, and how much to produce at the opened facilities).

The BendersDecomposition class in JijZept Tools automates the whole procedure — splitting the problem, generating cuts, and iterating. You only need to provide two things:

1. a problem written in JijModeling

2. a solver for the master problem (a function in the OMMX-adapter style)

The continuous subproblem is a linear program, so the class solves it internally using HiGHS (an open-source linear programming solver).

Example problem

Consider a small production-planning problem. For two sites A and B, we decide whether to open each site (binary variables \(x, y\)) and how much to produce at the opened sites (continuous variables \(z_A, z_B\)). Opening a site costs a fixed fee (20, 30), production has a unit cost (2, 3), and a total demand of 5 must be met.

\[\begin{split} \begin{aligned} \min \quad & 20x + 30y + 2 z_A + 3 z_B \\ \mathrm{s.t.} \quad & x + y \le 1 \\ & z_A \le 10x, \quad z_B \le 10y \\ & z_A + z_B \ge 5 \\ & x, y \in \{0,1\}, \quad 0 \le z_A, z_B \le 10 \end{aligned} \end{split}\]

BendersDecomposition automatically classifies the constraints in the problem by the kinds of variables they contain:

Constraint	Variables involved	Classification
choose_one (\(x + y \le 1\))	binary only	master problem
capacity_A, capacity_B (\(z_A \le 10x\), …)	binary + continuous	coupling constraints (the constraints linking master and subproblem)
demand (\(z_A + z_B \ge 5\))	continuous only	subproblem

import jijmodeling as jm

# Decision variables
x = jm.BinaryVar("x")  # open site A?
y = jm.BinaryVar("y")  # open site B?
z_a = jm.ContinuousVar("z_a", lower_bound=0, upper_bound=10)  # production at site A
z_b = jm.ContinuousVar("z_b", lower_bound=0, upper_bound=10)  # production at site B

problem = jm.Problem("facility_production")

# Objective: minimize fixed opening costs + production costs
problem += 20 * x + 30 * y + 2 * z_a + 3 * z_b

# Constraints
problem += jm.Constraint("choose_one", x + y <= 1)     # open at most one site (binary only)
problem += jm.Constraint("capacity_A", z_a <= 10 * x)  # production only at opened sites (coupling)
problem += jm.Constraint("capacity_B", z_b <= 10 * y)  # same as above (coupling)
problem += jm.Constraint("demand", z_a + z_b >= 5)     # meet the demand (continuous only)

problem

\[\begin{split}\begin{array}{cccc}\text{Problem:} & \text{facility\_production} & & \\& & \min \quad \displaystyle 20 \cdot x + 30 \cdot y + 2 \cdot z_a + 3 \cdot z_b & \\\text{{s.t.}} & & & \\ & \text{capacity\_A} & \displaystyle z_a \leq 10 \cdot x & \\ & \text{capacity\_B} & \displaystyle z_b \leq 10 \cdot y & \\ & \text{choose\_one} & \displaystyle x + y \leq 1 & \\ & \text{demand} & \displaystyle z_a + z_b \geq 5 & \\\text{{where}} & & & \\& x & 0\text{-dim binary variable}\\& y & 0\text{-dim binary variable}\\& z_a & 0\text{-dim continuous variable}\\ & & \text{lower bound: }0 & \\ & & \text{upper bound: }10 & \\& z_b & 0\text{-dim continuous variable}\\ & & \text{lower bound: }0 & \\ & & \text{upper bound: }10 & \\\end{array}\end{split}\]

Running the algorithm

Instantiate BendersDecomposition with the problem and the instance data (the dict of values for jm.Placeholders; empty here because this example uses no placeholders), then simply call solve().

The first argument of solve() is the solver for the master problem — any function that takes an OMMX Instance and returns a Solution (or a SampleSet) will do. Here we pass the HiGHS OMMX adapter as-is.

With verbose=True, you can watch the lower bound (LB) and upper bound (UB) close in on each other at every iteration.

from ommx_highs_adapter import OMMXHighsAdapter

from jijzepttools.modeling.algorithm.benders_decomposition import BendersDecomposition

benders = BendersDecomposition(problem, instance_data={})

result = benders.solve(
    OMMXHighsAdapter.solve,  # solver for the master problem
    max_iterations=50,
    tolerance=1e-6,
    verbose=True,
)

Starting Benders decomposition with 2 binary vars, 2 continuous vars
Coupling constraints: 2
------------------------------------------------------------

Iteration   1: LB=-100000.000, UB=5010.000, Gap=105010.000, Feasible=True

Iteration   2: LB=-4970.000, UB=  30.000, Gap=5000.000, Feasible=True

Iteration   3: LB=  30.000, UB=  30.000, Gap=  -0.000, Feasible=True

Converged after 3 iterations!
Optimal value: 30.000000

Reading the result

solve() returns a BendersResult. Its main fields are:

• converged: whether the algorithm converged (the gap between lower and upper bound fell below tolerance)

• iteration: the number of iterations

• lower_bound / upper_bound: the final bounds

• is_feasible: whether the subproblem was feasible

• solution: the final solution (an OMMX Solution object)

As expected, the result is “open only the cheaper site A (\(x=1\)) and produce the full demand of 5 there”, with objective value \(20 + 2 \times 5 = 30\).

print("Converged:", result.converged)
print("Iterations:", result.iteration)
print("Objective value (upper bound):", result.upper_bound)

# Variable values of the final solution (from the OMMX Solution)
result.solution.decision_variables_df[["name", "value"]]

Converged: True
Iterations: 3
Objective value (upper bound): 30.0

	name	value
id
0	x	1.0
1	y	0.0
2	z_a	5.0
3	z_b	0.0

A look inside

At construction time, the variables and constraints are automatically classified into the master side, the subproblem side, and the coupling constraints. The cuts generated at each iteration accumulate in cuts.

If you want to advance one iteration at a time, call iterate() (runs a single iteration and returns a BendersResult) repeatedly instead of solve().

print("Binary variables (master side):", [v.name for v in benders.binary_vars])
print("Continuous variables (subproblem side):", [v.name for v in benders.continuous_vars])
print("Coupling constraints:", [c.name for c in benders.coupling_constraints])
print("Number of cuts added:", len(benders.cuts))

Binary variables (master side): ['x', 'y']
Continuous variables (subproblem side): ['z_a', 'z_b']
Coupling constraints: ['capacity_A', 'capacity_B']
Number of cuts added: 2

Swapping the master-problem solver

The master-problem solver is just “a function that takes an OMMX Instance and solves it”, so switching to another solver is a one-line change. For example, to use the OMMX adapter for SCIP (an open-source mixed-integer programming solver):

from ommx_pyscipopt_adapter import OMMXPySCIPOptAdapter

benders_scip = BendersDecomposition(problem, instance_data={})
result_scip = benders_scip.solve(OMMXPySCIPOptAdapter.solve, max_iterations=50, tolerance=1e-6)

print("Converged:", result_scip.converged)
print("Objective value:", result_scip.upper_bound)

Converged: True
Objective value: 30.0

Because the master problem contains only binary variables, you can also use QUBO-style samplers such as annealing (e.g. OpenJij) as the master-problem solver. When the sampler returns a SampleSet (a collection of multiple solutions), multi-cut mode kicks in: several cuts are generated from the samples in a single iteration, which can reduce the number of iterations.

Key options

The BendersDecomposition constructor accepts the following options, among others (see each argument’s docstring for details):

• slack_M: penalty coefficient for the slack variables (auxiliary variables representing constraint violations) that keep the computation going even when the subproblem is infeasible (default 1000)

• use_two_phase: whether to use a two-phase scheme that first secures feasibility (Phase 1) and then optimizes (Phase 2)

• normalize_cuts: whether to normalize cut coefficients (useful when coefficient magnitudes differ widely)

• max_cuts: maximum number of cuts to keep (least-active cuts are dropped when exceeded; default 100)

• max_samples: when a SampleSet is received, the maximum number of samples used for cut generation (picked in order of objective value)

Related features

• Decomposing a problem along a “series” such as the time axis is covered in the Problem Decomposition tutorial.

• Low-level building blocks for Dantzig-Wolfe decomposition (functions that generate the master problem and subproblems used in column generation) are also available as jijzepttools.modeling.algorithm.dantzig_wolfe_decomposition (no high-level API yet).