import gurobipy as gp
from gurobipy import GRB
import numpy as np
class LpSolverGurobi:
def __init__(self, obj_coeff, constr_mat, rhs, verbose=False):
# initialize environment and model
self._env = gp.Env('GurobiEnv', empty=True)
# self._env.setParam('LogToConsole', 1 if verbose else 0)
self._env.setParam('OutputFlag', 1 if verbose else 0)
self._env.start()
self._model = gp.Model(env=self._env, name='GurobiLpSolver')
# prepare data
self._obj_coeff = obj_coeff
# print(self._obj_coeff)
self._constr_mat = constr_mat
# print(self._constr_mat)
self._rhs = rhs
self._num_vars = len(self._obj_coeff)
self._num_constrs = len(self._rhs)
# create decision variables
self._vars = self._model.addMVar(self._num_vars,
vtype=GRB.CONTINUOUS,
lb=0)
# create constraints
self._constrs = self._model.addConstr(
self._constr_mat@self._vars == self._rhs
)
# create objective
self._model.setObjective(self._obj_coeff@self._vars,
GRB.MINIMIZE)
def optimize(self, verbose=False):
self._model.optimize()
if self._model.status == GRB.OPTIMAL:
print(f'Optimal solution found!')
print(f'Optimal objective = {self._model.objVal:.2f}')
elif self._model.status == GRB.UNBOUNDED:
print(f'Model is unbounded!')
elif self._model.status == GRB.INFEASIBLE:
print(f'Model is infeasible!')
else:
print(f'Unknown error occurred!')
def save_model(self, filename):
self._model.write(filename)
def clean_up(self):
self._model.dispose()
self._env.dispose()3 Solving Linear Programming Problems with Benders Decomposition
In this chapter, we use Benders decomposition to solve several linear programming (LP) problems in order to demonstrate the decomposition logic, especially how the restricted Benders master problem interacts with the subproblem in an iterative approach to reach final optimality. Most linear programs could be solved efficiently nowadays by either open source or commercial solvers without resorting to any decomposition approaches. However, by working through the example problems in the following sections, we aim to showcase the implementation details when applying Benders decomposition algorithm on real problems, which helps solidify our understanding of Benders decomposition. Hopefully, by the end of this chapter, we will build up enough intuition as well as hands-on experience such that we are ready to tackle most involved problems in the following chapters.
In the following sections, we employ several steps to illustrate the problem solving process of Benders decomposition.
- We will first create two linear programming solvers based on Gurobi and SCIP that can solve any linear programs defined in the standard form. They are used in later section to validate the correctness of the solutions produced by Benders decomposition.
- Next, we use a specific linear program and give the corresponding RBMP and DSP to prepare for the implementations.
- Then, we will solve the example linear program step by step by examining the outputs of the RBMP and DSP to decided the next set of actions.
- Futhermore, a holistic Benders decomposition implementation is then developed to solve the example linear program.
- Following the previous step, a more generic Benders decomposition implementation is created.
- Then, we will examine an alternative implementation using Gurobi callback functions.
- We will also provide an implementation based on SCIP.
- In the final section, we will do several benchmarking testing.
3.0.1 LP solvers based on Gurobi and SCIP
We aim to use Benders decomposition to solve several linear programming problems in the following sections. To do that, we intentionally decompose the LP problem under consideration into two sets, one set of complicating variables and the other set containing the remaining variables. In order to validate the correctness of the results obtained by Benders decomposition, we implement two additional ways of solving the target linear programming problems directly. The first option is based on the Gurobi API in python and the other is based on the open source solve SCIP. The two implementations defined here assume the LP problems under consideration follow the below standard form:
\[ \begin{aligned} \text{min.} &\quad \mathbf{c}^T \mathbf{x} \\ \text{s.t.} &\quad \mathbf{A} \mathbf{x} = \mathbf{b} \\ &\quad \mathbf{x} \geq 0 \end{aligned} \]
Listing 3.1 defines a solver for LP problems using Gurobi. It takes three constructor parameters:
obj_coeff: this corresponds to the objective coefficients \(\mathbf{c}\).constr_mat: this refers to the constraint matrix \(\mathbf{A}\).rhs: this is the right-hand side \(\mathbf{b}\).
Inside the constructor __init__(), a solver environment _env is first created and then used to initialize a model object _model. The input parameters are then used to create decision variables _vars, constraints _constrs and objective function respectively. The optimize() function simply solves the problem and shows the solving status. Finally, the clean_up() function frees up the computing resources.
Listing 3.2 presents an LP solver implementation in class LpSolverSCIP using SCIP. The constructor requires the same of parameters as defined in LpSolverGurobi. The model building process is similar with minor changes when required to create decision variables, constraints and the objective function.
import pyscipopt as scip
from pyscipopt import SCIP_PARAMSETTING
class LpSolverSCIP:
def __init__(self, obj_coeff, constr_mat, rhs, verbose=False):
self._model = scip.Model('LpModel')
if not verbose:
self._model.hideOutput()
# create variables
self._vars = {
i: self._model.addVar(lb=0, vtype='C')
for i in range(len(obj_coeff))
}
# create constraints
for c in range(len(rhs)):
expr = [
constr_mat[c][j] * self._vars.get(j)
for j in range(len(obj_coeff))
]
self._model.addCons(scip.quicksum(expr) == rhs[c])
# create objective
obj_expr = [
obj_coeff[i] * self._vars.get(i)
for i in range(len(obj_coeff))
]
self._model.setObjective(scip.quicksum(obj_expr), "minimize")
def optimize(self):
self._model.optimize()
status = self._model.getStatus()
if status == "optimal":
print(f'Optimal solution found!')
print(f'Optimal objective = {self._model.getObjVal():.2f}')
elif status == "unbounded":
print(f'Model is unbounded!')
elif status == "infeasible":
print(f'Model is infeasible!')
else:
print(f'Unknown error occurred!')Listing 3.3 generates a LP problem with 20 decision variables and 5 constraints.
import numpy as np
np.random.seed(42)
c = np.random.randint(1, 6, size=20)
A = np.random.randint(-10, 12, size=(5, 20))
b = np.random.randint(20, 100, size=5)Listing 3.4 solves the generated LP using LpSolverGurobi and the solver output shows that an optimal solution was found with objective value of 36.90.
lpsolver_gurobi = LpSolverGurobi(obj_coeff=c, constr_mat=A, rhs=b)
lpsolver_gurobi.optimize()Optimal solution found!
Optimal objective = 36.90Listing 3.5 solves the same LP problem using SCIP and not surprisingly, the same optimal objective value was found. This is not exciting, as it only indicates that the two solvers agree on the optimal solution on such a small LP problem as expected. However, they will become more useful in the following sections when we use them to validate our Benders decomposition results.
lpsolver_scip = LpSolverSCIP(obj_coeff=c, constr_mat=A, rhs=b)
lpsolver_scip.optimize()Optimal solution found!
Optimal objective = 36.903.0.2 A serious LP problem that cannot wait to be decomposed!
With the validation tools available for use, we are ready to solve some serious LP problems using Benders decomposition! What we have below is a LP problem with five decision variables, three of which are denoted by \(\mathbf{x} = (x_1, x_2, x_3)\) and the remaining two variables are denoted by \(\mathbf{y} = (y_1, y_2)\). We assume that \(\mathbf{y}\) are the complicating variables.
\[ \begin{aligned*} \text{min.} &\quad 8x_1 + 12x_2 +10x_3 + 15y_1 + 18y_2 \\ \text{s.t.} &\quad 2x_1 + 3x_2 + 2x_3 + 4y_1 + 5y_2 = 300 \\ &\quad 4x_1 + 2x_2 + 3x_3 + 2y_1 + 3y_2 = 220 \\ &\quad x_i \geq 0, \ \forall i = 1, \cdots, 3 \\ &\quad y_i \geq 0, \ \forall j = 1, 2 \end{aligned*} \]
According to the standard LP form presented in the previous section, \(\mathbf{c}^T = (8, 12, 10)\), \(\mathbf{f}^T = (15, 18)\) and \(\mathbf{b}^T = (300, 220)\). In addition,
\[\begin{equation*} \mathbf{A} = \begin{bmatrix} 2 & 3 & 2 \\ 4 & 2 & 3 \\ \end{bmatrix} \qquad \mathbf{B} = \begin{bmatrix} 4 & 5 \\ 2 & 3 \\ \end{bmatrix} \end{equation*}\]
Listing 3.6 solves the LP problem using the two solvers define above. Both solvers confirm the optimal objective value is 1091.43. With this information in mind, we will apply Benders decomposition to see if the same optimal solution could be identified or not.
c = np.array([8, 12, 10])
f = np.array([15, 18])
obj_coeff = np.concatenate([c, f])
A = np.array([[2, 3, 2],
[4, 2, 3]])
B = np.array([[4, 5],
[2, 3]])
constr_mat = np.concatenate([A, B], axis=1)
rhs = np.array([300, 220])
lpsolver_gurobi = LpSolverGurobi(obj_coeff, constr_mat, rhs, verbose=False)
lpsolver_gurobi.optimize()
# Optimal objective = 1091.43
lpsolver_scip = LpSolverSCIP(obj_coeff, constr_mat, rhs, verbose=False)
lpsolver_scip.optimize()
# Optimal objective = 1091.43Optimal solution found!
Optimal objective = 1091.43
Optimal solution found!
Optimal objective = 1091.433.0.3 Benders decomposition formulations
With \(\mathbf{y}\) being the complicating variable, we state the Benders subproblem (SP) below for the \(\mathbf{y}\) assuming fixed values \(\mathbf{\bar{y}} = (\bar{y_1}, \bar{y_2})\):
\[ \begin{aligned*} \text{min.} &\quad 8x_1 + 12x_2 +10x_3 \\ \text{s.t.} &\quad 2x_1 + 3x_2 + 2x_3 = 300 - 4\bar{y_1} - 5\bar{y_2} \\ &\quad 4x_1 + 2x_2 + 3x_3 = 220 - 2\bar{y_1} - 3\bar{y_2} \\ &\quad x_i \geq 0, \ \forall i = 1, \cdots, 3 \end{aligned*} \]
We define the dual variable \(\mathbf{u} = (u_1, u_2)\) to associate with the constraints in the (SP). The dual subproblem (DSP) could then be stated as follows:
\[ \begin{aligned*} \text{max.} &\quad (300 - 4\bar{y_1} - 5\bar{y_2}) u_1 + (220 - 2\bar{y_1} - 3\bar{y_2}) u_2 \\ \text{s.t.} &\quad 2u_1 + 4u_2 \leq 8\\ &\quad 3u_1 + 2u_2 \leq 12 \\ &\quad 2u_1 + 3u_2 \leq 10 \\ &\quad u_1, u_2\ \text{unrestricted} \end{aligned*} \]
The (RBMP) can be stated as below. Note that \(\mathbf{u} = (0, 0)\) is a feasible solution to (DSP) and the corresponding objective value is 0, which is the reason we restrict the variable \(g\) to be nonnegative.
\[ \begin{aligned*} \text{min.} &\quad 15 y_1 + 18 y_2 + g \\ \text{s.t.} &\quad y_1, y_2 \geq 0 \\ &\quad g \geq 0 \end{aligned*} \]
3.0.4 Benders decomposition step by step
Benders decomposition defines a problem solving process in which the restricted Benders master problem and the dual subproblem interact iteratively to identify the optimal solution or conclude infeasibility/unboundedness. To facilitate our understanding of the process, we demonstrate in this section the workings of Benders decomposition by solving the target LP problem step by step.
Listing 3.7 shows the codes that initialize the lower bound lb, upper bound ub and threshold value eps. Furthermore, the restricted Benders master problem rbmp is created with three variables and a minimizing objective function. Note that no constraints are yet included in the model at this moment.
import numpy as np
import gurobipy as gp
from gurobipy import GRB
# initialize lower/upper bounds and threshold value
lb = -GRB.INFINITY
ub = GRB.INFINITY
eps = 1.0e-5
# create restricted Benders master problem
env = gp.Env('benders', empty=True)
env.setParam('OutputFlag', 0)
env.start()
rbmp = gp.Model(env=env, name='RBMP')
# create decision variables
y1 = rbmp.addVar(vtype=GRB.CONTINUOUS, lb=0, name='y1')
y2 = rbmp.addVar(vtype=GRB.CONTINUOUS, lb=0, name='y2')
g = rbmp.addVar(vtype=GRB.CONTINUOUS, lb=0, name='g')
# create objective
rbmp.setObjective(15*y1 + 18*y2 + g, GRB.MINIMIZE)Listing 3.8 initializes the Benders subproblem with two decision variables and three constraints.
# create dual subproblem
dsp = gp.Model(env=env, name='DSP')
# create decision variables
u1 = dsp.addVar(vtype=GRB.CONTINUOUS,
lb=-GRB.INFINITY,
ub=GRB.INFINITY,
name='u1')
u2 = dsp.addVar(vtype=GRB.CONTINUOUS,
lb=-GRB.INFINITY,
ub=GRB.INFINITY,
name='u2')
# create objective function
dsp.setObjective(300*u1 + 220*u2)
# create constraints
dsp.addConstr(2*u1 + 4*u2 <= 8, name='c1')
dsp.addConstr(3*u1 + 2*u2 <= 12, name='c2')
dsp.addConstr(2*u1 + 3*u2 <= 10, name='c3')
dsp.update()In Listing 3.9, we solve the (RBMP) for the first time. It has an optimal solution with \((\bar{y_1}, \bar{y_2}, \bar{g}) = (0, 0, 0)\) and optimal objective value of 0. This is expected as all the variables assume their minimal possible values in order to minimize the objective function. This objective value also serves as the new lower bound.
rbmp.optimize()
if rbmp.status == GRB.OPTIMAL:
print(f'Optimal solution found! Objective value = {rbmp.objVal:.2f}')
y1_opt, y2_opt, g_opt = y1.X, y2.X, g.X
lb = np.max([lb, rbmp.objVal])
print(f'(y1, y2, g) = ({y1_opt:.2f}, {y2_opt:.2f}, {g_opt:.2f})')
print(f'lb={lb}, ub={ub}')
elif rbmp.status == GRB.INFEASIBLE:
print(f'original problem is infeasible!')Optimal solution found! Objective value = 0.00
(y1, y2, g) = (0.00, 0.00, 0.00)
lb=0.0, ub=1e+100Given that \((\bar{y_1}, \bar{y_2}, \bar{g}) = (0, 0, 0)\), we now feed the values of \(\bar{y_1}\) and \(\bar{y_2}\) into the Benders dual subproblem (DSP) by updating its objective function, as shown in Listing 3.10:
# update objective function
dsp.setObjective((300-4*y1_opt-5*y2_opt)*u1 +
(220-2*y1_opt-3*y2_opt)*u2,
GRB.MAXIMIZE)
dsp.update()
dsp.optimize()
if dsp.status == GRB.OPTIMAL:
u1_opt, u2_opt = u1.X, u2.X
print(f'Optimal objective = {dsp.objVal:.2f}')
print(f'(u1, u2) = ({u1_opt:.2f}, {u2_opt:.2f})')
ub = np.min([ub, 15*y1_opt + 18*y2_opt + dsp.objVal])
print(f'lb={lb}, ub={ub}')Optimal objective = 1200.00
(u1, u2) = (4.00, 0.00)
lb=0.0, ub=1200.0We see that the dual subproblem has an optimal solution. The upper bound ub is also updated. Since the optimal objective value of the subproblem turns out to be 1200 and is greater than \(\bar{g} = 0\), which implies that an optimality cut is needed to make sure that the variable \(g\) in the restricted Benders master problem reflects this newly obtained information from the subproblem.
In Listing 3.11, the new optimality cut is added to the (RBMP), which is then solved to optimality.
# add optimality cut
rbmp.addConstr((300-4*y1-5*y2)*u1_opt
+ (220-2*y1-3*y2)*u2_opt <= g,
name='c1')
rbmp.update()
rbmp.optimize()
if rbmp.status == GRB.OPTIMAL:
print(f'Optimal solution found! Objective value = {rbmp.objVal:.2f}')
y1_opt, y2_opt, g_opt = y1.X, y2.X, g.X
lb = np.max([lb, rbmp.objVal])
print(f'(y1, y2, g) = ({y1_opt:.2f}, {y2_opt:.2f}, {g_opt:.2f})')
print(f'lb={lb}, ub={ub}')
elif rbmp.status == GRB.INFEASIBLE:
print(f'original problem is infeasible!')Optimal solution found! Objective value = 1080.00
(y1, y2, g) = (0.00, 60.00, 0.00)
lb=1080.0, ub=1200.0Armed with the optimal solution \((\bar{y_1}, \bar{y_2}, \bar{g}) = (0, 60, 0)\), Listing 3.12 updates the objective function of the (DSP) and obtains its optimal solution.
# update objective function
dsp.setObjective((300-4*y1_opt-5*y2_opt)*u1
+ (220-2*y1_opt-3*y2_opt)*u2,
GRB.MAXIMIZE)
dsp.update()
dsp.optimize()
if dsp.status == GRB.OPTIMAL:
u1_opt, u2_opt = u1.X, u2.X
print(f'Optimal objective = {dsp.objVal:.2f}')
print(f'(u1, u2) = ({u1_opt:.2f}, {u2_opt:.2f})')
ub = np.min([ub, 15*y1_opt + 18*y2_opt + dsp.objVal])
print(f'lb={lb}, ub={ub:.2f}')
elif dsp.Status == GRB.UNBOUNDED:
print(f'DSP is unbounded!')
u1_ray = u1.UnbdRay
u2_ray = u2.UnbdRay
print(f'retrieve extreme ray (u1, u2) = ({u1_ray}, {u2_ray})')
else:
print(f'DSP solve error')DSP is unbounded!
retrieve extreme ray (u1, u2) = (-2.0, 1.0)Since the dual subproblem is unbounded, a feasibility cut is further needed. In Listing 3.13, we add the new cut and solve the restricted Benders master problem again.
# add optimality cut
rbmp.addConstr((300-4*y1-5*y2)*u1_ray
+ (220-2*y1-3*y2)*u2_ray <= 0,
name='c2')
rbmp.update()
rbmp.optimize()
if rbmp.status == GRB.OPTIMAL:
print(f'Optimal solution found! Objective value = {rbmp.objVal:.2f}')
y1_opt, y2_opt, g_opt = y1.X, y2.X, g.X
lb = np.max([lb, rbmp.objVal])
print(f'(y1, y2, g) = ({y1_opt:.2f}, {y2_opt:.2f}, {g_opt:.2f})')
print(f'lb={lb:.2f}, ub={ub:.2f}')
elif rbmp.status == GRB.INFEASIBLE:
print(f'original problem is infeasible!')Optimal solution found! Objective value = 1091.43
(y1, y2, g) = (0.00, 54.29, 114.29)
lb=1091.43, ub=1200.00Note that a new lower bound is obtained after solving the master problem. Since there is still a large gap between the lower bound and upper bound, we continue solving the subproblem in Listing 3.14.
# update objective function
dsp.setObjective((300-4*y1_opt-5*y2_opt)*u1
+ (220-2*y1_opt-3*y2_opt)*u2,
GRB.MAXIMIZE)
dsp.update()
dsp.optimize()
if dsp.status == GRB.OPTIMAL:
u1_opt, u2_opt = u1.X, u2.X
print(f'Optimal objective = {dsp.objVal:.2f}')
print(f'(u1, u2) = ({u1_opt:.2f}, {u2_opt:.2f})')
ub = np.min([ub, 15*y1_opt + 18*y2_opt + dsp.objVal])
print(f'lb={lb:.2f}, ub={ub:.2f}')Optimal objective = 114.29
(u1, u2) = (4.00, 0.00)
lb=1091.43, ub=1091.43Now that the difference between lb and ub is less than the preset threshold eps, we conclude that an optimal solution is reached and the computation resources are freed up.
# release resources
rbmp.dispose()
dsp.dispose()
env.dispose()3.0.5 Benders decomposition automated
It typically takes Benders decomposition many iterations to reach optimality or conclude infeasibility/unboundedness. In this section, we put everything we have learned from the manual approach above into an automatic workflow.
Listing 3.15 defines an Enum class that specifies four possible optimization statuses. The meanings of these statuses are self-explanatory from their corresponding names and further explanations are omitted here.
import gurobipy as gp
from gurobipy import GRB
import numpy as np
from enum import Enum
class OptStatus(Enum):
OPTIMAL = 0
UNBOUNDED = 1
INFEASIBLE = 2
ERROR = 3Listing 3.16 defines a ManualBendersMasterSolver class that models the (RBMP). Its constructor contains the variable and objective function definitions. The ability to take in either feasibility or optimality cuts is implemented in separate functions add_feasibility_cut() and add_optimality_cut(), respectively. The solve() function is responsible for optimizing the model and retrieve the optimal solution if any.
class ManualBendersMasterSolver:
def __init__(self, env):
self._model = gp.Model(env=env, name='RBMP')
# create decision variables
self._y1 = self._model.addVar(vtype=GRB.CONTINUOUS,
lb=0, name='y1')
self._y2 = self._model.addVar(vtype=GRB.CONTINUOUS,
lb=0, name='y2')
self._g = self._model.addVar(vtype=GRB.CONTINUOUS,
lb=0, name='g')
# create objective
self._model.setObjective(15*self._y1+18*self._y2+self._g,
GRB.MINIMIZE)
self._opt_obj = None
self._opt_y1 = None
self._opt_y2 = None
self._opt_g = None
def solve(self) -> OptStatus:
print('-' * 50)
print(f'Start solving master problem.')
self._model.optimize()
opt_status = None
if self._model.status == GRB.OPTIMAL:
opt_status = OptStatus.OPTIMAL
self._opt_obj = self._model.objVal
self._opt_y1 = self._y1.X
self._opt_y2 = self._y2.X
self._opt_g = self._g.X
print(f'\tmaster problem is optimal.')
print(f'\topt_obj={self._opt_obj:.2f}')
print(f'\topt_y1={self._opt_y1:.2f}, opt_y2={self._opt_y2:.2f}')
print(f'\topt_g={self._opt_g:.2f}')
elif self._model.status == GRB.INFEASIBLE:
print(f'\tmaster problem is infeasible.')
opt_status = OptStatus.INFEASIBLE
else:
print(f'\tmaster problem encountered error.')
opt_status = OptStatus.ERROR
print(f'Finish solving master problem.')
print('-' * 50)
return opt_status
def add_feasibility_cut(self, ray_u1, ray_u2) -> None:
self._model.addConstr((300-4*self._y1-5*self._y2)*ray_u1 +
(220-2*self._y1-3*self._y2)*ray_u2 <=
0)
print(f'Benders feasibility cut added!')
def add_optimality_cut(self, opt_u1, opt_u2) -> None:
self._model.addConstr((300-4*self._y1-5*self._y2)*opt_u1 +
(220-2*self._y1-3*self._y2)*opt_u2 <=
self._g)
print(f'Benders optimality cut added!')
def clean_up(self):
self._model.dispose()
@property
def opt_obj(self):
return self._opt_obj
@property
def opt_y1(self):
return self._opt_y1
@property
def opt_y2(self):
return self._opt_y2
@property
def opt_g(self):
return self._gListing 3.17 defines the Benders dual subproblem in a similar fashion. Notice that the update_objective() function is used to set an updated objective function based on the optimal solution identified in the restricted Benders master problem.
class ManualBendersSubprobSolver:
def __init__(self, env):
self._model = gp.Model(env=env, name='DSP')
# create decision variables
self._u1 = self._model.addVar(vtype=GRB.CONTINUOUS,
lb=-GRB.INFINITY,
ub=GRB.INFINITY,
name='u1')
self._u2 = self._model.addVar(vtype=GRB.CONTINUOUS,
lb=-GRB.INFINITY,
ub=GRB.INFINITY,
name='u2')
# create constraints
self._model.addConstr(2*self._u1+4*self._u2 <= 8, name='c1')
self._model.addConstr(3*self._u1+2*self._u2 <= 12, name='c2')
self._model.addConstr(2*self._u1+3*self._u2 <= 10, name='c3')
self._model.setObjective(1, GRB.MAXIMIZE)
self._model.update()
self._opt_obj = None
self._opt_u1 = None
self._opt_u2 = None
self._ray_u1 = None
self._ray_u2 = None
def solve(self):
print('-' * 50)
print(f'Start solving dual subproblem.')
self._model.optimize()
status = None
if self._model.status == GRB.OPTIMAL:
self._opt_obj = self._model.objVal
self._opt_u1 = self._u1.X
self._opt_u2 = self._u2.X
status = OptStatus.OPTIMAL
print(f'\tdual subproblem is optimal.')
print(f'\topt_obj={self._opt_obj:.2f}')
print(f'\topt_y1={self._opt_u1:.2f}, opt_y2={self._opt_u2:.2f}')
elif self._model.status == GRB.UNBOUNDED:
status = OptStatus.UNBOUNDED
print(f'\tdual subproblem is unbounded!')
self._ray_u1 = self._u1.UnbdRay
self._ray_u2 = self._u2.UnbdRay
print(f'\textreme ray (u1, u2) = ({self._ray_u1}, {self._ray_u2})')
else:
status = OptStatus.ERROR
print(f'Finish solving dual subproblem.')
print('-' * 50)
return status
def update_objective(self, opt_y1, opt_y2):
self._model.setObjective(
(300-4*opt_y1-5*opt_y2)*self._u1 +
(220-2*opt_y1-3*opt_y2)*self._u2,
GRB.MAXIMIZE)
print(f'dual subproblem objective updated!')
def clean_up(self):
self._model.dispose()
@property
def opt_obj(self):
return self._opt_obj
@property
def opt_u1(self):
return self._opt_u1
@property
def opt_u2(self):
return self._opt_u2
@property
def ray_u1(self):
return self._ray_u1
@property
def ray_u2(self):
return self._ray_u2Listing 3.18 shows the control flow of Benders decomposition. The main logic is stated as a while loop, in which the master problem and dual subproblem are solved sequentially within each iteration. Depending on whether the subproblem is optimal or unbounded, an optimality or feasibility cut is added to the master problem. The process continues until the gap between the lower bound and the upper bound is within a certain threshold.
class ManualBendersDecomposition:
def __init__(self, master_solver, dual_subprob_solver):
self._master_solver = master_solver
self._dual_subprob_solver = dual_subprob_solver
def optimize(self) -> OptStatus:
eps = 1.0e-5
lb = -np.inf
ub = np.inf
iter = 1
while True:
print(f"\nIteration: {iter}")
iter += 1
# solve master problem
master_status = self._master_solver.solve()
if master_status == OptStatus.INFEASIBLE:
return OptStatus.INFEASIBLE
# update lower bound
lb = np.max([lb, self._master_solver.opt_obj])
print(f'Bounds: lb={lb:.2f}, ub={ub:.2f}')
opt_y1 = self._master_solver.opt_y1
opt_y2 = self._master_solver.opt_y2
# solve subproblem
self._dual_subprob_solver.update_objective(opt_y1, opt_y2)
dsp_status = self._dual_subprob_solver.solve()
if dsp_status == OptStatus.OPTIMAL:
# update upper bound
opt_obj = self._dual_subprob_solver.opt_obj
ub = np.min([ub, 15*opt_y1 + 18*opt_y2 + opt_obj])
print(f'Bounds: lb={lb:.2f}, ub={ub:.2f}')
if ub - lb <= eps:
break
opt_u1 = self._dual_subprob_solver.opt_u1
opt_u2 = self._dual_subprob_solver.opt_u2
self._master_solver.add_optimality_cut(opt_u1, opt_u2)
elif dsp_status == OptStatus.UNBOUNDED:
ray_u1 = self._dual_subprob_solver.ray_u1
ray_u2 = self._dual_subprob_solver.ray_u2
self._master_solver.add_feasibility_cut(ray_u1, ray_u2)
Listing 3.19 solves the LP problem using the wholesome Benders solver. The optimal solution agrees with the solution obtained in the manual approach.
env = gp.Env('benders', empty=True)
env.setParam("OutputFlag",0)
env.start()
master_solver = ManualBendersMasterSolver(env)
dual_subprob_solver = ManualBendersSubprobSolver(env)
benders_decomposition = ManualBendersDecomposition(master_solver, dual_subprob_solver)
benders_decomposition.optimize()
Iteration: 1
--------------------------------------------------
Start solving master problem.
master problem is optimal.
opt_obj=0.00
opt_y1=0.00, opt_y2=0.00
opt_g=0.00
Finish solving master problem.
--------------------------------------------------
Bounds: lb=0.00, ub=inf
dual subproblem objective updated!
--------------------------------------------------
Start solving dual subproblem.
dual subproblem is optimal.
opt_obj=1200.00
opt_y1=4.00, opt_y2=0.00
Finish solving dual subproblem.
--------------------------------------------------
Bounds: lb=0.00, ub=1200.00
Benders optimality cut added!
Iteration: 2
--------------------------------------------------
Start solving master problem.
master problem is optimal.
opt_obj=1080.00
opt_y1=0.00, opt_y2=60.00
opt_g=0.00
Finish solving master problem.
--------------------------------------------------
Bounds: lb=1080.00, ub=1200.00
dual subproblem objective updated!
--------------------------------------------------
Start solving dual subproblem.
dual subproblem is unbounded!
extreme ray (u1, u2) = (-2.0, 1.0)
Finish solving dual subproblem.
--------------------------------------------------
Benders feasibility cut added!
Iteration: 3
--------------------------------------------------
Start solving master problem.
master problem is optimal.
opt_obj=1091.43
opt_y1=0.00, opt_y2=54.29
opt_g=114.29
Finish solving master problem.
--------------------------------------------------
Bounds: lb=1091.43, ub=1200.00
dual subproblem objective updated!
--------------------------------------------------
Start solving dual subproblem.
dual subproblem is optimal.
opt_obj=114.29
opt_y1=4.00, opt_y2=0.00
Finish solving dual subproblem.
--------------------------------------------------
Bounds: lb=1091.43, ub=1091.433.0.6 Benders decomposition - design and implementation
The Benders decomposition algorithm we put together in the last section is a big step towards automating the interaction between the master problem optimizer and the dual subproblem optimizer. However, it is still limited in that both optimizers are tied to a specific problem instance. Even a slight change in the instance data requires an updated optimizer implementation. Ideally, we would like to to have an algorithm that could solve various LP problems as long as they follow the same form.
In addition, we might want to have the flexibility to switch to different solvers within the master and dual subproblem optimizers. A good algorithm design should allow alternative implementations with minimal impact on the overall algorithm structure.
With these needs in mind, in this section, we aim to design a Benders decomposition algorithm that can solve any LP problems following a standard form and allow various implementations in both the master problem and dual subproblem optimizers.
We develop the algorithm in three steps. First, we present a Benders decomposition algorithm design that focuses on the overall algorithm workflow and only the abstract implementations of the master problem and dual subproblem optimizers. Next, an implementation of the Benders decomposition based on Gurobi is demonstrated. Lastly, an alternative implementation based on the open source solver SCIP is illustrated.
3.0.6.1 Benders decomposition algorithm design
We could see from the previous sections that Benders decomposition involves three key components:
- An optimizer that solves the master problem
- An optimizer that solves the dual subproblem
- An orchestrator that dictates the interaction between the two optimizers
To be able to solve LP problems of various sizes, both the master and dual subproblem optimizers must take a generic form and not be tied to a fixed number of various and/or constraints. Furthermore, to be able to switch to different solvers, both optimizers need to conform to a common interface.
To this end, Listing 3.21 defines a base class for master problem optimizers. The constructor __init__() takes three inputs, namely, objective coefficients, constraint matrix and the right-hand side. Any concrete implementation of the base class is responsible for constructing a model from these three input values using solver-specific modeling languages. The base class also defines several other key functions:
solve(): this is the function that invokes solver-specific optimization process and saves the corresponding optimization results.add_feasibility_cut(): this function takes an extreme ray identified by the dual subproblem optimizer and adds a feasibility cut to the master problem.add_optimality_cut(): this function takes an optimal solution identified by the dual subproblem optimizer and adds an optimality cut to the master problem.opt_obj_val(): this function returns the optimal objective value obtained.opt_val_for_complicating_vars(): this function returns the optimal solution values for all the complicating decision variables.opt_val_for_surrogate_var(): this function returns the optimal value of the surrogate variable.
from abc import ABC, abstractmethod
from typing import Dict
import numpy as np
class BendersMasterOptimizer(ABC):
"""base class for master solver implementation
Args:
ABC (ABCMeta): helper class
"""
@abstractmethod
def __init__(self, objective_coefficients: np.array,
constraint_matrix: np.array,
right_hand_side: np.array):
"""constructor for master problem model initialization.
Note that two types of variables will be created:
- the complicating variables
- the surrogate variable
Args:
objective_coefficients (np.array): an array of size n that
defines the coefficient value for each variable in the
master problem objective function.
constraint_matrix (np.array): a matrix of size m * n that
defines the constraint coefficients.
right_hand_side (np.array): an array of size n that
defines the right hand side for each constraint.
"""
raise NotImplementedError
@abstractmethod
def solve(self) -> OptStatus:
"""solve the problem and return optimization status
Returns:
OptStatus: an enum object
"""
raise NotImplementedError
@abstractmethod
def add_feasibility_cut(self, extreme_ray: Dict[int, float]) -> None:
"""add feasibility cut.
Args:
extreme_ray (Dict[int, float]): a mapping between variable key
and extreme ray element. For example:
extreme_ray = {
0: 1.0,
1: 2.0
}
It is assumed that the key type is integer.
"""
raise NotImplementedError
@abstractmethod
def add_optimality_cut(self, opt_sol: Dict[int, float]) -> None:
"""add optimality cut.
Args:
opt_sol (Dict[int, float]): a mapping between variable key
and optimal solution element. For example:
opt_sol = {
0: 1.0,
1: 2.0
}
It is assumed that the key type is integer.
"""
raise NotImplementedError
@abstractmethod
def opt_obj_val(self) -> float:
"""return the optimal objective value if exists
Returns:
float: objective value
"""
raise NotImplementedError
@abstractmethod
def opt_obj_val_comp(self) -> float:
"""return the optimal objective value for the complicating variables
Returns:
float: objective value for the complicating variables
"""
raise NotImplementedError
@abstractmethod
def opt_val_for_complicating_vars(self) -> Dict[int, float]:
"""return the optimal solution values for complicating variables
Returns:
Dict[int, float]: mapping between variable key and optimal value.
For example, opt_val = {
0: 2.0,
1: 3.0
}
"""
raise NotImplementedError
@abstractmethod
def opt_val_for_surrogate_var(self) -> float:
"""return the optimal value for the surrogate variable
Returns:
float: optimal value for the surrogate variable
"""
raise NotImplementedErrorSimilarly, Listing 3.21 presents a base class for dual subproblem optimizers. The constructor __init__() takes the objective coefficients and constraint matrix that correspond to the variables in the subproblem. These information are required to create variables and constraints for the dual subproblem. However, the objective function will not be instantiated without the optimal solution of the restricted master problem. This base class also contains several other key functions:
solve(): this is the function that calls solver-specific procedure to solve the underlying dual subproblem.update_objective(): this function takes the optimal solution obtained by the master problem optimizer and updates the objective function for the dual subproblem.opt_obj_val(): this function returns the optimal objective values identified by the optimizer.opt_sol(): this function returns the optimal solution identified by the optimizer.extreme_ray(): this function returns the extreme ray identified by the optimizer.
from abc import ABC, abstractmethod
from typing import Dict
import numpy as np
class BendersDspOptimizer(ABC):
"""base class for dual subproblem implementation
Args:
ABC (ABCMeta): helper class
"""
@abstractmethod
def __init__(self, objective_coefficients: np.array,
constraint_matrix: np.array):
"""constructor for dual subproblem model initialization.
Args:
objective_coefficients (np.array): an array of size m that
represents the coefficient value for each variable in the
original problem objective function.
constraint_matrix (np.array): a matrix of size m * n that
defines the constraint coefficients in the original problem.
"""
raise NotImplementedError
@abstractmethod
def solve(self) -> OptStatus:
"""solve the problem and return optimization status
Returns:
OptStatus: an enum object
"""
raise NotImplementedError
@abstractmethod
def update_objective(self, opt_sol_rbmp: Dict[int, float],
constraint_matrix_rbmp: np.array,
right_hand_side_rbmp: np.array) -> None:
"""update the objective function for the dual subproblem.
Args:
opt_sol_rbmp (Dict[int, float]): optimal solution of the
restricted master problem
constraint_matrix_rbmp (np.array): constraint matrix
associated with the complicating variables in the
master problem
right_hand_side_rbmp (np.array): right-hand side of
the original problem
"""
raise NotImplementedError
@abstractmethod
def opt_obj_val(self) -> float:
"""return the optimal objective value if exists
Returns:
float: objective value
"""
raise NotImplementedError
@abstractmethod
def opt_sol(self) -> Dict[int, float]:
"""return the optimal solution values
Returns:
Dict[int, float]: mapping between variable key and optimal value.
For example, opt_sol = {
0: 2.0,
1: 4.0
}
"""
raise NotImplementedError
@abstractmethod
def extreme_ray(self) -> Dict[int, float]:
"""return the identified extreme ray
Returns:
Dict[int, float]: mapping between variable key and ray element.
For example, extreme_ray = {
0: 3.5,
1: 6.0
}
"""
raise NotImplementedErrorclass BendersOptimizerConsole:
def __init__(self, master_optimizer: BendersMasterOptimizer,
dsp_optimizer: BendersDspOptimizer):
self._master_optimizer = master_optimizer
self._dsp_optimizer = dsp_optimizer
def optimize(self, verbose=False) -> OptStatus:
eps = 1.0e-5
lb = -np.inf
ub = np.inf
iter = 0
while True:
iter += 1
if verbose:
print(f'\nIteration: {iter}')
# solve master problem
master_status = self._master_optimizer.solve()
if master_status == OptStatus.INFEASIBLE:
if verbose:
print(f'Model is infeasible!')
return OptStatus.INFEASIBLE
# update lower bound
lb = np.max([lb, self._master_optimizer.opt_obj_val])
if verbose:
print(f'Bounds: lb={lb:.2f}, ub={ub:.2f}')
# solve subproblem
opt_val_comp = self._master_optimizer.opt_val_for_complicating_vars
self._dsp_optimizer.update_objective(opt_val_comp)
dsp_status = self._dsp_optimizer.solve()
if dsp_status == OptStatus.OPTIMAL:
if verbose:
print(f'DSP is optimal!')
# update upper bound
opt_obj = self._dsp_optimizer.opt_obj_val
opt_obj_val_comp = self._master_optimizer.opt_obj_val_comp
ub = np.min([ub, opt_obj_val_comp + opt_obj])
if verbose:
print(f'Bounds: lb={lb:.2f}, ub={ub:.2f}')
if ub - lb <= eps:
return OptStatus.OPTIMAL
opt_sol = self._dsp_optimizer.opt_sol
self._master_optimizer.add_optimality_cut(opt_sol)
elif dsp_status == OptStatus.UNBOUNDED:
if verbose:
print(f'DSP is unbounded!!!')
extreme_ray = self._dsp_optimizer.extreme_ray
self._master_optimizer.add_feasibility_cut(extreme_ray)
else:
if verbose:
print(f"DSP solve ERROR!!!")
return OptStatus.ERROR3.0.6.2 Master and subproblem solvers based on Gurobi
Listing 3.23 implements a master problem solver class named GenericMasterSolverGurobi for (RBMP) characterized by the cost coefficient f, the constraint matrix B and right-hand side b. The constructor initializes the complicating variables _y and the dummy variable _g. The solve() function solves the problem and retrieves the optimal solution if exists. Moreover, the solver provides functions to add optimality and feasibility cuts to the existing model.
class BendersMasterOptimizerGurobi(BendersMasterOptimizer):
def __init__(self, f: np.array, B: np.array, b: np.array):
# save data
self._f = f
self._B = B
self._b = b
# env and model
self._env = gp.Env('MasterEnv', empty=True)
self._env.setParam("OutputFlag",0)
self._env.start()
self._model = gp.Model(env=self._env, name='MasterSolver')
# create variables
self._num_y_vars = len(f)
self._y = self._model.addVars(self._num_y_vars,
lb=0,
vtype=GRB.CONTINUOUS,
name='y')
self._g = self._model.addVar(vtype=GRB.CONTINUOUS,
lb=0,
name='g')
# create objective
self._model.setObjective(
gp.quicksum(self._f[i] * self._y.get(i)
for i in range(self._num_y_vars)) +
self._g,
GRB.MINIMIZE)
self._model.update()
self._opt_obj = None
self._opt_obj_y = None
self._opt_y = None
self._opt_g = None
def solve(self) -> OptStatus:
print('-' * 50)
print(f'Start solving master problem.')
self._model.optimize()
opt_status = None
if self._model.status == GRB.OPTIMAL:
opt_status = OptStatus.OPTIMAL
self._opt_obj = self._model.objVal
self._opt_y = {
i: self._y.get(i).X
for i in range(self._num_y_vars)
}
self._opt_g = self._g.X
self._opt_obj_y = self._opt_obj - self._opt_g
print(f'\tmaster problem is optimal.')
print(f'\topt_obj={self._opt_obj:.2f}')
print(f'\topt_g={self._opt_g:.2f}')
for j in range(self._num_y_vars):
print(f'\topt_y{j}={self._y.get(j).X}')
elif self._model.status == GRB.INFEASIBLE:
print(f'\tmaster problem is infeasible.')
opt_status = OptStatus.INFEASIBLE
else:
print(f'\tmaster problem encountered error.')
opt_status = OptStatus.ERROR
print(f'Finish solving master problem.')
print('-' * 50)
return opt_status
def add_feasibility_cut(self, ray_u: dict) -> None:
constr_expr = [
ray_u.get(u_idx) *
(
self._b[u_idx] -
gp.quicksum(
self._B[u_idx][j] * self._y.get(j)
for j in range(self._num_y_vars)
)
)
for u_idx in ray_u.keys()
]
self._model.addConstr(gp.quicksum(constr_expr) <= 0)
print(f'Benders feasibility cut added!')
def add_optimality_cut(self, opt_u: dict) -> None:
constr_expr = [
opt_u.get(u_idx) *
(
self._b[u_idx] -
gp.quicksum(
self._B[u_idx][j] * self._y.get(j)
for j in range(self._num_y_vars)
)
)
for u_idx in opt_u.keys()
]
self._model.addConstr(gp.quicksum(constr_expr) <= self._g)
self._model.update()
print(f'Benders optimality cut added!')
def clean_up(self):
self._model.dispose()
self._env.dispose()
def save_model(self, filename):
self._model.write(filename)
@property
def f(self):
return self._f
@property
def opt_obj(self):
return self._opt_obj
@property
def opt_obj_y(self):
return self._opt_obj_y
@property
def opt_y(self):
return self._opt_y
@property
def opt_g(self):
return self._gListing 3.24 presents a solver for (DSP) that’s defined by the objective function coefficient c and constraint matrix A. The model objective function could be updated by update_objective() with the latest value of y. Notice that, the optimal solution is saved to _opt_u if the underlying problem is optimal. Otherwise, an extreme ray is retrieved and stored in _extreme_ray.
class GenericSubprobSolverGurobi:
def __init__(self, A: np.array, c: np.array, B: np.array, b: np.array):
# save data
self._A = A
self._c = c
self._b = b
self._B = B
# env and model
self._env = gp.Env('SubprobEnv', empty=True)
self._env.setParam("OutputFlag",0)
self._env.start()
self._model = gp.Model(env=self._env, name='SubprobSolver')
# create variables
self._num_vars = len(b)
self._u = self._model.addVars(self._num_vars,
vtype=GRB.CONTINUOUS,
lb=-GRB.INFINITY,
ub=GRB.INFINITY,
name='u')
# create constraints
for c_idx in range(len(c)):
self._model.addConstr(
gp.quicksum(A[:,c_idx][i] * self._u.get(i)
for i in range(len(b))) <= c[c_idx]
)
self._opt_obj = None
self._opt_u = None
self._ray_u = None
def solve(self):
print('-' * 50)
print(f'Start solving dual subproblem.')
self._model.setParam(GRB.Param.DualReductions, 0)
self._model.setParam(GRB.Param.InfUnbdInfo, 1)
self._model.optimize()
status = None
if self._model.status == GRB.OPTIMAL:
self._opt_obj = self._model.objVal
self._opt_u = {
i: self._u.get(i).X
for i in range(self._num_vars)
}
status = OptStatus.OPTIMAL
print(f'\tdual subproblem is optimal.')
print(f'\topt_obj={self._opt_obj:.2f}')
for i in range(self._num_vars):
print(f'\topt_u{i}={self._u.get(i).X}')
elif self._model.status == GRB.UNBOUNDED:
status = OptStatus.UNBOUNDED
self._ray_u = {
i: self._u.get(i).UnbdRay
for i in range(self._num_vars)
}
print(f'dual subproblem is unbounded')
else:
print(f'dual subproblem solve ERROR!')
status = OptStatus.ERROR
print(f'Finish solving dual subproblem.')
print('-' * 50)
return status
def update_objective(self, opt_y: dict):
obj_expr = [
self._u.get(u_idx) *
(
self._b[u_idx] -
sum(self._B[u_idx][j] * opt_y.get(j)
for j in range(len(opt_y))
)
)
for u_idx in range(self._num_vars)
]
self._model.setObjective(gp.quicksum(obj_expr), GRB.MAXIMIZE)
print(f'dual subproblem objective updated!')
def clean_up(self):
self._model.dispose()
self._env.dispose()
def save_model(self, filename):
self._model.write(filename)
@property
def opt_obj(self):
return self._opt_obj
@property
def opt_u(self):
return self._opt_u
@property
def ray_u(self):
return self._ray_uIn Listing 3.26, we utilize the freshly baked solver to tackle the serious LP problem presented in the previous sections. The output shows that the same optimal objective value of 1091.43 was obtained in the end.
import numpy as np
c = np.array([8, 12, 10])
f = np.array([15, 18])
A = np.array([[2, 3, 2], [4, 2, 3]])
B = np.array([[4, 5], [2, 3]])
b = np.array([300, 220])
master_solver = GenericMasterSolverGurobi(f, B, b)
dual_subprob_solver = GenericSubprobSolverGurobi(A, c, B, b)
benders_solver = GenericBendersDecomposition(master_solver, dual_subprob_solver)
status = benders_solver.optimize()
Iteration: 1
--------------------------------------------------
Start solving master problem.
master problem is optimal.
opt_obj=0.00
opt_g=0.00
opt_y0=0.0
opt_y1=0.0
Finish solving master problem.
--------------------------------------------------
Bounds: lb=0.00, ub=inf
dual subproblem objective updated!
--------------------------------------------------
Start solving dual subproblem.
dual subproblem is optimal.
opt_obj=1200.00
opt_u0=4.0
opt_u1=0.0
Finish solving dual subproblem.
--------------------------------------------------
DSP is optimal!
Bounds: lb=0.00, ub=1200.00
Benders optimality cut added!
Iteration: 2
--------------------------------------------------
Start solving master problem.
master problem is optimal.
opt_obj=1080.00
opt_g=0.00
opt_y0=0.0
opt_y1=60.0
Finish solving master problem.
--------------------------------------------------
Bounds: lb=1080.00, ub=1200.00
dual subproblem objective updated!
--------------------------------------------------
Start solving dual subproblem.
dual subproblem is unbounded
Finish solving dual subproblem.
--------------------------------------------------
DSP is unbounded!!!
Benders feasibility cut added!
Iteration: 3
--------------------------------------------------
Start solving master problem.
master problem is optimal.
opt_obj=1091.43
opt_g=114.29
opt_y0=0.0
opt_y1=54.285714285714285
Finish solving master problem.
--------------------------------------------------
Bounds: lb=1091.43, ub=1200.00
dual subproblem objective updated!
--------------------------------------------------
Start solving dual subproblem.
dual subproblem is optimal.
opt_obj=114.29
opt_u0=4.0
opt_u1=0.0
Finish solving dual subproblem.
--------------------------------------------------
DSP is optimal!
Bounds: lb=1091.43, ub=1091.433.0.6.3 Master and subproblem solvers based on SCIP
import pyscipopt
class GenericMasterSolverSCIP:
def __init__(self, f: np.array, B: np.array, b: np.array):
# save data
self._f = f
self._B = B
self._b = b
# env and model
self._model = pyscipopt.Model('MasterSolver')
# create variables
self._num_y_vars = len(f)
self._y = {
i: self._model.addVar(lb=0, vtype='C', name=f'y{i}')
for i in range(self._num_y_vars)
}
self._g = self._model.addVar(lb=0,
vtype='C',
name='g')
# create objective
self._model.setObjective(
pyscipopt.quicksum(self._f[i] * self._y.get(i)
for i in range(self._num_y_vars)) +
self._g,
'maximize')
self._opt_obj = None
self._opt_obj_y = None
self._opt_y = None
self._opt_g = None
def solve(self) -> OptStatus:
print('-' * 50)
print(f'Start solving master problem.')
self._model.setPresolve(pyscipopt.SCIP_PARAMSETTING.OFF)
self._model.setHeuristics(pyscipopt.SCIP_PARAMSETTING.OFF)
self._model.disablePropagation()
self._model.optimize()
opt_status = None
print(self._model.getStatus())
if self._model.getStatus() == 'optimal':
opt_status = OptStatus.OPTIMAL
self._opt_obj = self._model.getObjVal()
self._opt_y = {
i: self._model.getVal(self._y.get(i))
for i in range(self._num_y_vars)
}
self._opt_g = self._model.getVal(self._g)
self._opt_obj_y = self._opt_obj - self._opt_g
print(f'\tmaster problem is optimal.')
print(f'\topt_obj={self._opt_obj:.2f}')
print(f'\topt_g={self._opt_g:.2f}')
for j in range(self._num_y_vars):
print(f'\topt_y{j}={self._opt_y.get(j)}')
elif self._model.getStatus() == 'infeasible':
print(f'\tmaster problem is infeasible.')
opt_status = OptStatus.INFEASIBLE
else:
print(f'\tmaster problem encountered error.')
opt_status = OptStatus.ERROR
print(f'Finish solving master problem.')
print('-' * 50)
return opt_status
def add_feasibility_cut(self, ray_u: dict) -> None:
constr_expr = [
ray_u.get(u_idx) *
(
self._b[u_idx] -
pyscipopt.quicksum(
self._B[u_idx][j] * self._y.get(j)
for j in range(self._num_y_vars)
)
)
for u_idx in ray_u.keys()
]
self._model.addCons(pyscipopt.quicksum(constr_expr) <= 0)
print(f'Benders feasibility cut added!')
def add_optimality_cut(self, opt_u: dict) -> None:
constr_expr = [
opt_u.get(u_idx) *
(
self._b[u_idx] -
pyscipopt.quicksum(
self._B[u_idx][j] * self._y.get(j)
for j in range(self._num_y_vars)
)
)
for u_idx in opt_u.keys()
]
self._model.addCons(pyscipopt.quicksum(constr_expr) <= self._g)
self._model.update()
print(f'Benders optimality cut added!')
@property
def f(self):
return self._f
@property
def opt_obj(self):
return self._opt_obj
@property
def opt_obj_y(self):
return self._opt_obj_y
@property
def opt_y(self):
return self._opt_y
@property
def opt_g(self):
return self._gimport pyscipopt
class GenericSubprobSolverSCIP:
def __init__(self, A: np.array, c: np.array, B: np.array, b: np.array):
# save data
self._A = A
self._c = c
self._b = b
self._B = B
# env and model
self._model = pyscipopt.Model('SubprobSolver')
# create variables
self._num_vars = len(b)
self._u = {
i: self._model.addVar(vtype='C',
lb=-self._model.infinity(),
ub=self._model.infinity(),
name=f'u{i}')
for i in range(self._num_vars)
}
# create constraints
for c_idx in range(len(c)):
self._model.addCons(
pyscipopt.quicksum(A[:,c_idx][i] * self._u.get(i)
for i in range(len(b))) <= c[c_idx]
)
self._opt_obj = None
self._opt_u = None
self._ray_u = None
def solve(self):
print('-' * 50)
print(f'Start solving dual subproblem.')
self._model.setPresolve(pyscipopt.SCIP_PARAMSETTING.OFF)
self._model.setHeuristics(pyscipopt.SCIP_PARAMSETTING.OFF)
self._model.disablePropagation()
self._model.optimize()
status = None
if self._model.getStatus() == 'optimal':
self._opt_obj = self._model.getObjVal()
self._opt_u = {
i: self._model.getVal(self._u.get(i))
for i in range(self._num_vars)
}
status = OptStatus.OPTIMAL
print(f'\tdual subproblem is optimal.')
print(f'\topt_obj={self._opt_obj:.2f}')
for i in range(self._num_vars):
print(f'\topt_u{i}={self._opt_u.get(i)}')
elif self._model.getStatus() == 'unbounded':
status = OptStatus.UNBOUNDED
hasRay = self._model.hasPrimalRay()
print(f'primal ray exists!')
ray = self._model.getPrimalRay()
self._ray_u = {
i: ray[i]
for i in range(self._num_vars)
}
print(f'dual subproblem is unbounded')
else:
print(f'dual subproblem solve ERROR!')
status = OptStatus.ERROR
print(f'Finish solving dual subproblem.')
print('-' * 50)
return status
def update_objective(self, opt_y: dict):
obj_expr = [
self._u.get(u_idx) *
(
self._b[u_idx] -
sum(self._B[u_idx][j] * opt_y.get(j)
for j in range(len(opt_y))
)
)
for u_idx in range(self._num_vars)
]
self._model.setObjective(pyscipopt.quicksum(obj_expr), GRB.MAXIMIZE)
print(f'dual subproblem objective updated!')
@property
def opt_obj(self):
return self._opt_obj
@property
def opt_u(self):
return self._opt_u
@property
def ray_u(self):
return self._ray_uListing 3.26 solves same the LP problem that we have been tackling in the previous sections. The output confirms that the same optimal solution is identified as in the last section.
import numpy as np
c = np.array([8, 12, 10])
f = np.array([15, 18])
A = np.array([[2, 3, 2], [4, 2, 3]])
B = np.array([[4, 5], [2, 3]])
b = np.array([300, 220])
master_solver = GenericMasterSolverSCIP(f, B, b)
dual_subprob_solver = GenericSubprobSolverSCIP(A, c, B, b)
benders_solver = GenericBendersDecomposition(master_solver, dual_subprob_solver)
benders_solver.optimize()
Iteration: 1
--------------------------------------------------
Start solving master problem.
unbounded
master problem encountered error.
Finish solving master problem.
--------------------------------------------------
presolving:
(0.0s) symmetry computation started: requiring (bin +, int +, cont +), (fixed: bin -, int -, cont -)
(0.0s) no symmetry present (symcode time: 0.00)
presolving (0 rounds: 0 fast, 0 medium, 0 exhaustive):
0 deleted vars, 0 deleted constraints, 0 added constraints, 0 tightened bounds, 0 added holes, 0 changed sides, 0 changed coefficients
0 implications, 0 cliques
presolved problem has 3 variables (0 bin, 0 int, 0 impl, 3 cont) and 0 constraints
Presolving Time: 0.00
time | node | left |LP iter|LP it/n|mem/heur|mdpt |vars |cons |rows |cuts |sepa|confs|strbr| dualbound | primalbound | gap | compl.
* 0.0s| 1 | 0 | 0 | - | LP | 0 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | -- | -- | 0.00%| unknown
0.0s| 1 | 0 | 0 | - | 568k | 0 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | -- | -- | 0.00%| unknown
SCIP Status : problem is solved [unbounded]
Solving Time (sec) : 0.00
Solving Nodes : 1
Primal Bound : +1.00000000000000e+20 (1 solutions)
Dual Bound : +1.00000000000000e+20
Gap : 0.00 %TypeError: '>=' not supported between instances of 'float' and 'NoneType'3.0.6.4 Detect infeasibility
import numpy as np
np.random.seed(142)
c = np.random.randint(2, 6, size=20)
f = np.random.randint(1, 15, size=10)
A = np.random.randint(2, 6, size=(20, 20))
B = np.random.randint(2, 26, size=(20, 10))
b = np.random.randint(20, 50, size=20)
obj_coeff = np.concatenate([c, f])
constr_mat = np.concatenate([A, B], axis=1)
rhs = bgurobi_solver = LpSolverGurobi(obj_coeff, constr_mat, rhs)
gurobi_solver.save_model('problem2.lp')
gurobi_solver.optimize()Model is infeasible!scip_solver = LpSolverSCIP(obj_coeff, constr_mat, rhs)
scip_solver.optimize()Model is infeasible!master_solver = GenericLpMasterSolver(f, B, b)
dual_subprob_solver = GenericLpSubprobSolver(A, c, B, b)
benders_solver = GenericBendersSolver(master_solver, dual_subprob_solver)
benders_solver.optimize()
Iteration: 1
--------------------------------------------------
Start solving master problem.
master problem is optimal.
opt_obj=0.00
opt_g=0.00
opt_y0=0.0
opt_y1=0.0
Finish solving master problem.
--------------------------------------------------
Bounds: lb=0.00, ub=inf
dual subproblem objective updated!
--------------------------------------------------
Start solving dual subproblem.
dual subproblem is optimal.
opt_obj=1200.00
opt_u0=4.0
opt_u1=0.0
Finish solving dual subproblem.
--------------------------------------------------
DSP is optimal!
Bounds: lb=0.00, ub=1200.00
Benders optimality cut added!
Iteration: 2
--------------------------------------------------
Start solving master problem.
master problem is optimal.
opt_obj=1080.00
opt_g=0.00
opt_y0=0.0
opt_y1=60.0
Finish solving master problem.
--------------------------------------------------
Bounds: lb=1080.00, ub=1200.00
dual subproblem objective updated!
--------------------------------------------------
Start solving dual subproblem.
dual subproblem is unbounded
Finish solving dual subproblem.
--------------------------------------------------
DSP is unbounded!!!
Benders feasibility cut added!
Iteration: 3
--------------------------------------------------
Start solving master problem.
master problem is optimal.
opt_obj=1091.43
opt_g=114.29
opt_y0=0.0
opt_y1=54.285714285714285
Finish solving master problem.
--------------------------------------------------
Bounds: lb=1091.43, ub=1200.00
dual subproblem objective updated!
--------------------------------------------------
Start solving dual subproblem.
dual subproblem is optimal.
opt_obj=114.29
opt_u0=4.0
opt_u1=0.0
Finish solving dual subproblem.
--------------------------------------------------
DSP is optimal!
Bounds: lb=1091.43, ub=1091.43