AMPL Development Tutorial 6/6 – Implementing Benders Decomposition with ampls#
Description: This concluding notebook in our six-part series delves into enhancing the efficiency of our decomposition algorithm by utilizing AMPL Solver Libraries (ampls).
We showcase the use of ampls to implement our Benders decomposition approach for solving our stochastic facility location problem.
AMPL’s intuitive syntax allows for straightforward allocation of variables and constraints to master and subproblems.
However, each iteration typically involves solving the problems from scratch.
To enhance efficiency, we introduce a method using ampls to export the master problem to its solver representation.
This allows us to add cuts directly to the solver’s matrix, enabling solvers to efficiently utilize most of the previously computed solution.
This technique significantly boosts the overall process efficiency allowing to readily embed algorithmic development into production applications.
Tags: amplpy, ampl, ampls, mip, stochastic, facility location, benders
Notebook author: Christian Valente <ccv@ampl.com>, Gyorgy Matyasfalvi <gyorgy@ampl.com>
References:
AMPL a Modeling Language for Mathematical Programming – Robert Fourer et al.
Intro to Linear Optimization – Dimitris Bertsimas et al.
SCIP Optimization Suite example projects – Stephen J. Maher (ZIB)
# Install dependencies
%pip install -q amplpy amplpy-cplex amplpy-gurobi pandas
# Google Colab & Kaggle integration
from amplpy import AMPL, ampl_notebook
ampl = ampl_notebook(
modules=["gurobi", "cplex"], # modules to install
license_uuid="default", # license to use
) # instantiate AMPL object and register magics
# Import all necessary libraries
import amplpy_cplex as ampls
SOLVER = "cplex"
import pandas as pd
import math
Two-stage stochastic programming problem formulation#
Consider a decision-making process occurring over two sequential stages. Initially, a decision vector \(x\) is selected. Following this, new information becomes available, prompting the selection of a second-stage decision vector \(y\). We introduce \(S\) possible scenarios for the newly acquired information, with the actual scenario being revealed only after the selection of \(x\). Each scenario is indexed by \(s\), and its occurrence probability is denoted by \(\alpha^s\), which we assume to be positive. Given the sequential nature of decision-making, the choice of \(y\) can be scenario-dependent, denoted as \(y^s\) to reflect this.
The decision vectors \(x\) and \(y^s\) incur costs represented by vectors \(f\) and \(q\), respectively. The first-stage decisions are constrained by:
Furthermore, the combination of first and second-stage decisions must satisfy the scenario-specific constraints:
for every scenario \(s\). The uncertainty in the model is encapsulated by \(\xi := (q, B, D, d)\), where a constraint is deemed part of the second stage if it involves any nonzero coefficients of variables \({y}\). The matrices \({B}^s\) are referred to as technology matrices, which determine the influence of the decision at the first stage on the second stage. The objective is to determine the choices of \(x\) and \(y^s\) for all \(s=1, \ldots, S\) that minimize the expected total cost:
This formulation leads to the following two-stage stochastic optimization problem:
we will refer to the above problem as the master problem. Notice that even if the number of possible scenarios \(S\) is moderate, this formulation can be a large linear programming problem. However, a decomposition method can help.
Benders decomposition for two-stage problems#
Given a vector \(\hat{x}\), using the dual simplex method solve the subproblem (for all \(s \in S\)):
\[\begin{split} \begin{array}{crl} \min & (q^s)^T y^s & \\ \textrm{s.t.} & D^s y^s & = d^s - B^s \hat{x} \\ & y^s & \geq 0 \end{array} \end{split}\]The dual of (5) takes the following form:
\[\begin{split} \begin{array}{crl} \max & (p^s)^T (d^s - B^s \hat{x}) & \\ \textrm{s.t.} & (p^s)^T D^s & <= (q^s)^T \end{array} \end{split}\]Let
\[\begin{split} \begin{array}{crrl} & P \;\; = & \{ p & | \;\; p^T D \leq q^T \} & \\ \end{array} \end{split}\]We assume that \(P\) is nonempty and has at least one extreme point (so we can apply duality theory). Let \(p_i\), \(i = 1, \ldots, I\), be the extreme points, and let \(w^j\), \(j = 1, \ldots, J\), be a complete set of extreme rays of \(P\).
Define \(z^s(\hat{x})^*\) as the optimal cost obtained from problem (5) described above. Returning to the optimization of \(x\), we encounter the following problem:
\[\begin{split} \begin{array}{crcrcrccrcl} \min & f^T x & + & \sum_{s=1}^S \alpha^s z^s(x) & \\ \textrm{s.t.} & A x & & & = {b} \\ & x & \geq 0 & & \end{array} \end{split}\]Under our assumption that the set \(P\) is nonempty, either the dual problem (6) has an optimal solution and \(z^s(\hat{x})\) is finite, or the optimal dual cost is infinite (indicating that the primal (5) is infeasible).
In particular, \(z^s(x) \leq \infty \leftrightarrow (w^j)^T (d^s - B^s x) \leq 0, \; \forall j\) and whenever \(z^s(x)\) is finite the optimum of (6) must be attained at an extreme point of the set \(P\) meaning \(z^s(x)\) is the smallest number \(z^s\) such that \((p^i)^T (d^s - B^s x) \leq z^s, \; \forall i\), leading us to the following master-problem formulation:
\[\begin{split} \begin{array}{crcrcrccrcl} \min & f^T x & + & \sum_{s=1}^S \alpha^s z^s & & \\ \textrm{s.t.} & A x & & & = b & \\ & (p^i)^T (d^s - B^s x) & & & \leq z^s & \forall \; i, s \\ & (w^j)^T (d^s - B^s x) & & & \leq 0 & \forall \; j, s \\ & x & \geq 0 & & & \end{array} \end{split}\]In our algorithm we will only generate constraints that we find to be violated by the current solution.
The algorithm#
Initialization step: initialize \((\hat{x}, \hat{z}^1, \ldots, \hat{z}^S)\) to zero.
Sub-problem step: using the dual-simplex method solve the sub-problems individually for every \(s \in S\):
Check optimality of master-problem: If for every \(s\), the corresponding sub-problem is feasible and \((q^s)^T \hat{y}^s \leq \hat{z}^s\) then all constraints are satisfied and we have an optimal solution to the master problem.
Add optimality cut: If the sub-problem corresponding to some \(s\) has an optimal solution such that \((q^s)^T \hat{y}^s > \hat{z}^s\) the following cut is added to the master problem:
where \(p^s\) is an optimal basic feasible solution to the dual of the sub-problem.
Add feasibility cut: If the sub-problem corresponding to some \(s\) is infeasible, its dual has infinite cost and the following cut is added to the master problem:
where \(w^s\) is a positive cost extreme ray.
Stochastic facility location a concrete example#
Given:
A set of facilities: \(I\).
A set of customers: \(J\).
Set of scenarios: \(S\) (representing different customer demands).
Task:
Find the minimum cost facilities to open such that the customer demand can be satisfied in all scenarios.
Variables#
\(x_i \in \{0, 1\} \quad \forall i \in I\)
\(x_i = 1\) if facility \(i\) is opened.
\(y_{ij}^s \geq 0 \quad \forall i \in I, \forall j \in J, \forall s \in S\)
\(y_{ij}^s\) is the level of demand for customer \(j\) satisfied by facility \(i\) in scenario \(s\).
Parameters:#
\(\alpha^s\): the probability of scenario \(s\).
\(f_i\): the fixed cost for opening facility \(i\),
\(q_{ij}\): the cost of servicing customer \(j\) from facility \(i\),
\(\lambda_j^s\): the demand of customer \(j\) in scenario \(s\),
\(k_i:\) the capacity of facility \(i\).
The extensive form#
The extensive form of our stochastic program can be formulated as follows:
The sub-problem#
Given the above formulation for the extensive form (12) we can express our scenario specific sub-problems as follows
The master-problem#
The master-problem takes the following form
Create amplpy instances of master and sub-problems#
create_common(): this function initializes a baseAMPLobject. It incorporates the model’s sets and parameters that are common to both the master and the sub-problem.create_sub_problem(): this function is dedicated to setting up theAMPLinstance for the sub-problem. It begins by invokingcreate_common(). Following this, it employs theeval()method to dynamically load the additional model entities and parameters specifically required for the sub-problem.create_master_problem(): similar to the sub-problem setup, this function also starts withcreate_common()to instantiate anAMPLobject with the basic shared model elements. It then utilizes theeval()method to incorporate additional model entities and parameters that are exclusive to the master problem.
def create_common() -> AMPL:
a = AMPL()
a.option["solver"] = SOLVER
a.option["gurobi_options"] = "presolve=0"
a.eval(
r"""
set FACILITIES;
set CUSTOMERS;
set SCENARIOS;
param prob{SCENARIOS} default 1/card(SCENARIOS);
param sub_scenario symbolic in SCENARIOS;
param customer_demand{CUSTOMERS, SCENARIOS} >= 0;
param fixed_cost{ FACILITIES } >= 0;
param facility_capacity{FACILITIES} >= 0;
param variable_cost{ FACILITIES, CUSTOMERS } >= 0;
"""
)
return a
def create_sub_problem() -> AMPL:
a = create_common()
a.eval(
r"""
var production{FACILITIES, CUSTOMERS, SCENARIOS} >= 0;
param sub_facility_open{ FACILITIES } default 1;
minimize operating_cost:
sum{ i in FACILITIES, j in CUSTOMERS }
variable_cost[i, j] * production[i, j, sub_scenario];
s.t.satisfying_customer_demand{ j in CUSTOMERS }:
sum{i in FACILITIES} production[i, j, sub_scenario] >= customer_demand[j, sub_scenario];
s.t.facility_capacity_limits{i in FACILITIES}:
sum{ j in CUSTOMERS } production[i, j, sub_scenario] <= facility_capacity[i] * sub_facility_open[i];
"""
)
return a
def create_master_problem() -> AMPL:
a = create_common()
a.eval(
r"""
var sub_variable_cost{SCENARIOS} >= 0;
var facility_open{ FACILITIES } binary;
minimize TotalCost :
sum{ i in FACILITIES } fixed_cost[i] * facility_open[i] + sum{s in SCENARIOS} prob[s]*sub_variable_cost[s];
s.t. sufficient_production_capacity:
sum{i in FACILITIES} facility_capacity[i]*facility_open[i] >= max{s in SCENARIOS} sum{j in CUSTOMERS} customer_demand[j,s];
"""
)
return a
Data#
We can use the same data file for both the master and the sub-problem.
%%writefile floc_ef.dat
set FACILITIES := Baytown_TX Beaumont_TX Baton_Rouge_LA;
set CUSTOMERS := San_Antonio_TX Dallas_TX Jackson_MS Birmingham_AL;
set SCENARIOS := Low Medium High;
param prob := Low 0.25 Medium 0.5 High 0.25;
param fixed_cost := Baytown_TX 400000 Beaumont_TX 200000 Baton_Rouge_LA 600000;
param facility_capacity := Baytown_TX 1550 Beaumont_TX 650 Baton_Rouge_LA 1750;
param variable_cost: San_Antonio_TX Dallas_TX Jackson_MS Birmingham_AL :=
Baytown_TX 5739.725 6539.725 8650.40 22372.1125
Beaumont_TX 6055.05 6739.055 8050.40 21014.225
Baton_Rouge_LA 8650.40 7539.055 4539.72 15024.325;
param customer_demand: Low Medium High :=
San_Antonio_TX 450 650 887
Dallas_TX 910 1134 1456
Jackson_MS 379 416 673
Birmingham_AL 91 113 207;
Overwriting floc_ef.dat
Utility functions#
We can declare a few python functions that will make the subsequent code more succint.
Firstly, to retain the ability of assigning data from AMPL data files but also allow assigning it directlry from python, we declare get_data_from_ampl that extracts the data from AMPL and creates the python data structures that will be needed when executing the algorithm. Note that in real life usage, the data will mostly come from Python itself.
def get_data_from_ampl(a: AMPL) -> tuple:
"""Get data we need in the algorithm from AMPL"""
facilities = a.set["FACILITIES"].get_values().to_list()
customers = a.set["CUSTOMERS"].get_values().to_list()
scenarios = a.set["SCENARIOS"].get_values().to_list()
facility_capacity = a.param["facility_capacity"].get_values().to_dict()
customer_demand = a.param["customer_demand"].get_values().to_dict()
variable_cost = a.param["variable_cost"].get_values().to_dict()
return (
facilities,
customers,
scenarios,
facility_capacity,
customer_demand,
variable_cost,
)
When exporting the model from AMPL to ampls we will lose the semantics of the model, because ampls stores its representation at solver level. Intuitively, instead of having the AMPL structured model, made of sets, parameters, variables and constraints, we have the Jacobian matrix, a cost vector and a vector of right hand sides.
To be able to add cuts (done via ampls) from values we get by solving the subproblems (via AMPL), we need some helper functions and some maps.
def get_tuple_map(model: ampls.AMPLModel, varname: str):
"""Get a map indexingsetitem -> index for a specified variable"""
var_map = model.get_var_map_filtered(varname)
beg = len(varname) + 2
return {name[beg:-2]: value for name, value in var_map.items()}
def get_maps(model) -> tuple:
"""Get all the maps we will need for the execution"""
index_facility_open = get_tuple_map(model, "facility_open")
index_sub_variable_cost = get_tuple_map(model, "sub_variable_cost")
revindex_facility_open = {
value: name for name, value in index_facility_open.items()
}
revindex_sub_variable_cost = {
value: name for name, value in index_sub_variable_cost.items()
}
return (
index_facility_open,
revindex_facility_open,
index_sub_variable_cost,
revindex_sub_variable_cost,
)
Finally, we write a function that adds the Benders cuts using ampls; the optimality and feasibility cuts are very similar to each other:
which, applied to this problem, correspond to the following (structured) AMPL:
s.t. optimality_cut {k in 1..nITER, s in SCENARIOS}:
sum{j in CUSTOMERS} customer_price[j,s,k]*customer_demand[j,s] +
sum{i in FACILITIES} facility_price[i,s,k]*facility_capacity[i]*facility_open[i] <= sub_variable_cost[s];
s.t. feasibility_cut {k in 1..nITER, s in SCENARIOS}:
sum{j in CUSTOMERS} customer_price[j,s,k]*customer_demand[j,s] +
sum{i in FACILITIES} facility_price[i,s,k]*facility_capacity[i]*facility_open[i] <= 0;
In ampls we lose the structure, so we need to pass the indices (in the Jacobian matrix) and the coefficients of each non-zero element of the cut.
Also, the (numeric) right hand side will have to be calculated explicitly.
Let’s start by reordering the formulation above: the RHS is actually sum{j in CUSTOMERS} customer_price[j,s,k]*customer_demand[j,s], as this is the only numerical data.
Let’s then move the \(z^s\) (which in AMPL is sub_variable_cost[s]) to the left. What we will pass to ampls is akin to:
sub_variable_cost[s] - sum{i in FACILITIES} facility_price[i,s,k]*facility_capacity[i]*facility_open[i]
>= sum{j in CUSTOMERS} customer_price[j,s,k]*customer_demand[j,s];
We therefore need the indices of the variables sub_variable_cost[s] in the master problem, along with the indices of all the facility_open[i]. We can get them from the maps we pass to the function: index_sub_variable_cost and index_facility_open respectively. Also note the negative coefficients for the product facility_price*facility_capacity*facility_open.
def add_benders_cut(
model: ampls.AMPLModel,
cdduals: dict,
fcduals: dict,
index_facility_open: dict,
index_sub_variable_cost: dict,
facility_capacity: dict,
customer_demand: dict,
s: str,
customers: list,
cutype: str,
):
"""Add a Bender's cut to the model."""
# Note that to accomodate for solver view, we formulate it as below
# sub_variable_cost[s]-sum{i in FACILITIES} facility_price[i,s,k]*facility_capacity[i]*facility_open[i]
# >= sum{j in CUSTOMERS} customer_price[j,s,k]*customer_demand[j,s];
indices = [index_facility_open[f] for f in fcduals.keys()]
coeffs = [-fcduals[f] * facility_capacity[f] for f in fcduals.keys()]
if cutype == "optimality":
indices.append(index_sub_variable_cost[s])
coeffs.append(1)
rhs = sum(cdduals[c] * customer_demand[(c, s)] for c in customers)
opt = model.addConstraint(indices, coeffs, ampls.CutDirection.GE, rhs)
# With the following, we keep track of the cuts we add, to then add them back into AMPL
model.record(opt)
print(f"Added {cutype} cut: {opt.toString()}\n")
Benders decomposition#
The main function is, at this point, easy to implement.
We create two instances of AMPL for the master and the sub problem respectively, clearly using the same data and set some options. At this stage, we can also get the data out from one of the AMPL instances to have it readily available in Python for the algorithm implementation. We also initialize some entities that we’ll use in the algorithm.
"""Generic function doing the optimization and reading the results"""
master = create_master_problem()
sub = create_sub_problem()
master.read_data("floc_ef.dat")
sub.read_data("floc_ef.dat")
# Set options
master.option["presolve"] = 0
sub.eval("suffix dunbdd;")
sub.option["presolve"] = 0 # set else AMPL likely presolves the variable facility_open
# Get the data we will use for iterations
(
facilities,
customers,
scenarios,
facility_capacity,
customer_demand,
variable_cost,
) = get_data_from_ampl(master)
# Constraints of subproblems that we will use to extract dual/unbounded rays
# Could do without but this makes the implementation slightly cleaner
satisfying_customer_demand = sub.con["satisfying_customer_demand"]
facility_capacity_limits = sub.con["facility_capacity_limits"]
# Initialize Benders's counter and params
epsilon = 0.00000001
sub_variable_cost = {s: 0 for s in scenarios}
it = 0
We use the function
AMPL.to_amplsto export an instance of the master problem to the ampls library. Note that at this point,solver_masteris a solver-level representation of the master problem, and we lose direct access to AMPL’s structured entities. This is why we also get the maps linking entities names to solver index.
# Export the master problem to ampls
solver_master = master.to_ampls(SOLVER)
# Get maps between solver vars and AMPL entities
(
index_facility_open,
revindex_facility_open,
index_sub_variable_cost,
revindex_sub_variable_cost,
) = get_maps(solver_master)
Now the master loop, which implements Benders. At every algorithm iteration, we:
Set the scenario index in the subproblem and solve it. We check the solution:
Infeasible -> We extract the dual rays and use them to construct a feasibility cut
Feasible but violating the optimality condition -> We extract the dual values and use them to construct an optimality cut
Feasible and not violating -> We mark this scenario as non-violating
If no scenarios had violations, we reached optimality and quit
Otherwise, we resolve the master problem (with the cuts added in step 1)
We use the solution to assign:
the parameter
sub_facility_openin the subproblemsthe dictionary
sub_variable_costused at the next iteration to check for violations
while True:
it += 1
print(f"****** Iteration {it} *******")
n_noviolations = 0
for s in scenarios:
# Solve the subproblem
sub.param["sub_scenario"] = s # Assign scenario
sub.get_output("solve;")
result = sub.get_value("solve_result")
print(f"Scenario {s} Objective={sub.get_current_objective().value()}")
# Decide what cut to add, if any
if result == "infeasible":
cdduals = satisfying_customer_demand.get_values("dunbdd").to_dict()
fcduals = facility_capacity_limits.get_values("dunbdd").to_dict()
add_benders_cut(
solver_master,
cdduals,
fcduals,
index_facility_open,
index_sub_variable_cost,
facility_capacity,
customer_demand,
s,
customers,
"feasibility",
)
elif sub.get_value("operating_cost") > sub_variable_cost[s] + epsilon:
cdduals = satisfying_customer_demand.get_values("dual").to_dict()
fcduals = facility_capacity_limits.get_values("dual").to_dict()
add_benders_cut(
solver_master,
cdduals,
fcduals,
index_facility_open,
index_sub_variable_cost,
facility_capacity,
customer_demand,
s,
customers,
"optimality",
)
else:
n_noviolations += 1
# If no scenario had a violation, we exit
if n_noviolations == len(scenarios):
break
print("SOLVING MASTER PROBLEM")
solver_master.optimize()
# Get the solution vector from ampls
sol = solver_master.get_solution_vector()
# AMPL: let {f in FACILITIES} sub_facility_open[f] := facility_open[f];
sub.param["sub_facility_open"] = {
f: sol[i] for i, f in revindex_facility_open.items()
}
# Assign the costs from the master problem for next iteration
sub_variable_cost = {s: sol[i] for i, s in revindex_sub_variable_cost.items()}
# End of Bender's loop
****** Iteration 1 *******
Scenario Low Objective=11621793.455
Added optimality cut: +1.000000*X[0]>= 11621793.455000;
Scenario Medium Objective=14779784.865
Added optimality cut: +308961.500000*X[3]+1.000000*X[1]>= 15088746.365000;
Scenario High Objective=21050711.200000003
Added optimality cut: +1548961.500000*X[3]+520000.000000*X[4]+1.000000*X[2]>= 23119672.700000;
SOLVING MASTER PROBLEM
****** Iteration 2 *******
Version identifier: 22.1.1.0 | 2022-11-28 | 9160aff4d
Scenario Low Objective=11621793.455
CPXPARAM_Simplex_Display 0
CPXPARAM_MIP_Display 0
CPXPARAM_Barrier_Display 0
Scenario Medium Objective=14966984.865
Added optimality cut: +1548961.500000*X[3]+520000.000000*X[4]+1.000000*X[1]>= 16515946.365000;
Scenario High Objective=21570711.200000003
SOLVING MASTER PROBLEM
Version identifier: 22.1.1.0 | 2022-11-28 | 9160aff4d
****** Iteration 3 *******
CPXPARAM_Simplex_Display 0
CPXPARAM_MIP_Display 0
CPXPARAM_Barrier_Display 0
Scenario Low Objective=11621793.455
Scenario Medium Objective=14779784.865
Scenario High Objective=21050711.200000003
Import solution back to AMPL#
In the following part of code, we have the correct solution in the object solver_master. We could settle at this, but often it is convenient to have the solution back in AMPL.
We therefore use the method AMPL.import_solution to import the solution values back to the AMPL object.
For demonstrational purposes, we also use the line:
master.eval(solver_master.get_recorded_entities())
to show how to import back into AMPL all the entities that have been added via ampls (in our case the cuts) and we print them in AMPL.
print(f"Optimal solution found, cost: {solver_master.get_obj()}\n\n")
# Import the AMPLS model back to AMPL, and show the additional cuts
master.import_solution(solver_master, keep_files=True)
master.eval(solver_master.get_recorded_entities())
master_obj = master.get_current_objective().value()
master.eval("expand {c in 1.._ncons} _con[c];")
Optimal solution found, cost: 16758018.596249998
CPLEX 22.1.1: optimal solution; objective 16758018.59625
0 simplex iterations
CPLEX 22.1.1: optimal solution; objective 16758018.59625
0 simplex iterations
subject to sufficient_production_capacity:
1550*facility_open['Baytown_TX'] + 650*facility_open['Beaumont_TX'] +
1750*facility_open['Baton_Rouge_LA'] >= 3223;
subject to ampls_con_1:
sub_variable_cost['Low'] >= 11621800;
subject to ampls_con_2:
sub_variable_cost['Medium'] + 308962*facility_open['Baytown_TX'] >=
15088700;
subject to ampls_con_3:
sub_variable_cost['High'] + 1548960*facility_open['Baytown_TX'] +
520000*facility_open['Beaumont_TX'] >= 23119700;
subject to ampls_con_4:
sub_variable_cost['Medium'] + 1548960*facility_open['Baytown_TX'] +
520000*facility_open['Beaumont_TX'] >= 16515900;
Check solution#
The following is just to double check the result: we solve the extensive form of the problem and check that its result coincides to the one we get from our implementation of Benders.
def create_ef() -> AMPL:
a = AMPL()
a.eval(
r"""
set FACILITIES;
set CUSTOMERS;
set SCENARIOS;
var facility_open{FACILITIES} binary;
var production{FACILITIES, CUSTOMERS, SCENARIOS} >= 0;
param fixed_cost{FACILITIES} >= 0;
param variable_cost{FACILITIES, CUSTOMERS} >= 0;
param customer_demand{CUSTOMERS, SCENARIOS} >= 0;
param facility_capacity{FACILITIES} >= 0;
param prob{SCENARIOS} default 1/card(SCENARIOS);
minimize TotalCost:
sum{i in FACILITIES} fixed_cost[i] * facility_open[i] +
sum{s in SCENARIOS, i in FACILITIES, j in CUSTOMERS} prob[s] * variable_cost[i,j] * production[i,j,s];
s.t. satisfying_customer_demand{s in SCENARIOS, j in CUSTOMERS}:
sum{i in FACILITIES} production[i,j,s] >= customer_demand[j,s];
s.t. facility_capacity_limits{s in SCENARIOS, i in FACILITIES}:
sum{j in CUSTOMERS} production[i,j,s] <= facility_capacity[i] * facility_open[i];
s.t. sufficient_production_capacity:
sum{i in FACILITIES} facility_capacity[i]*facility_open[i] >= max{s in SCENARIOS} sum{j in CUSTOMERS} customer_demand[j,s];
"""
)
a.option["solver"] = SOLVER
return a
# Check with the extensive form
ef = create_ef()
ef.read_data("floc_ef.dat")
ef.solve()
ef_obj = ef.get_current_objective().value()
print(f"Optimal solution EF form: {ef_obj}")
print(f"Optimal solution Benders: {master_obj}")
assert math.isclose(
ef_obj, master_obj, rel_tol=1e-9
), "EF and Benders objective values do not match!"
print("EF and Benders objective values match!")
CPLEX 22.1.1.0: optimal integer solution; objective 16758018.6
18 MIP simplex iterations
0 branch-and-bound nodes
Optimal solution EF form: 16758018.596250001
Optimal solution Benders: 16758018.596249994
EF and Benders objective values match!