AMPL Development Tutorial 6/6 – Implementing Benders Decomposition with ampls

6_benders_ampls_stoch_floc.ipynb Open In Colab Kaggle Gradient Open In SageMaker Studio Lab Hits

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 <>, Gyorgy Matyasfalvi <>


  • 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:

\begin{equation} \begin{array}{crlllllll} & A x & = & {b} \ & x & \geq & 0. \end{array} \tag{1} \end{equation}

Furthermore, the combination of first and second-stage decisions must satisfy the scenario-specific constraints:

\begin{equation} \begin{array}{crcllllll} & {B}^s x + {D}^s y^s & = & {d}^s \ & y^s & \geq & 0, \end{array} \tag{2} \end{equation}

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:

\begin{equation} \begin{array}{cllllllll} & f^T x & + & \alpha_1 {(q^1)^T} y^1 & + & \alpha^2 {(q^2)^T} y^2 & \ldots & + & \alpha^S {(q^S)^T} y^S \end{array} \tag{3} \end{equation}

This formulation leads to the following two-stage stochastic optimization problem:

\begin{equation} \begin{array}{crcrcrccrcl} \min & f^T x & + & \alpha^1 {(q^1)^T} y_1 & + & \alpha^2 {(q^1)^T} y^2 & \ldots & + & \alpha^s {(q^s)^T} y^S & & \ \textrm{s.t.} & A x & & & & & & & & &= {b} \ & {B^1} x & + & {D^1} y^1 & & & & & & &= {d^1} \ & {B^2} x & + & & & {D^2} y^2 & & & & &= {d^2} \ & \vdots & & & & & \ddots & & & &\vdots \ & {B^S} x & + & & & & & & {D^S} y^S & &= {d^S} \ & x \geq 0, & & y^1 \geq 0, & \ldots & & & & y^S \geq 0 & & \end{array} \tag{4} \end{equation}

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

  1. Given a vector $\hat{x}$, using the dual simplex method solve the subproblem (for all $s \in S$):

    \begin{equation} \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} \tag{5} \end{equation}

  2. The dual of (5) takes the following form:

    \begin{equation} \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} \tag{6} \end{equation}


    \begin{equation} \begin{array}{crrl} & P ;; = & { p & | ;; p^T D \leq q^T } & \ \end{array} \tag{7} \end{equation}

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

  3. 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{equation} \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} \tag{8} \end{equation}

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

  4. 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{equation} \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} \tag{9} \end{equation}

  5. In our algorithm we will only generate constraints that we find to be violated by the current solution.

The algorithm

  1. Initialization step: initialize $(\hat{x}, \hat{z}^1, \ldots, \hat{z}^S)$ to zero.

  2. Sub-problem step: using the dual-simplex method solve the sub-problems individually for every $s \in S$:

\begin{equation} \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} \tag{5} \end{equation}

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

  2. 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:
    \begin{equation} (p^s)^T (d^s - B^s x) \leq z^s, \tag{10} \end{equation} where $p^s$ is an optimal basic feasible solution to the dual of the sub-problem.

  3. 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:
    \begin{equation} (w^s)^T (d^s - B^s x) \leq 0, \tag{11} \end{equation} where $w^s$ is a positive cost extreme ray.

Stochastic facility location a concrete example


  • A set of facilities: $I$.

  • A set of customers: $J$.

  • Set of scenarios: $S$ (representing different customer demands).


  • Find the minimum cost facilities to open such that the customer demand can be satisfied in all scenarios.


  • $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$.


  • $\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:

$ \begin{equation} \begin{array}{rrlrrll} \min \quad & \sum_{i \in I} f_i x_i & + & \sum_{s \in S} \sum_{i \in I} \sum_{j \in J} \alpha^s q_{ij} y_{ij}^s & & & \ & & & & & & \ \textrm{s.t.} \quad & & & \sum_{i \in I} y_{ij}^s & \geq & \lambda_j^s & \forall j \in J, \forall s \in S \ & & & \sum_{j \in J} y_{ij}^s & \leq & k_i x_i & \forall i \in I, \forall s \in S \ & \sum_{i \in I} k_i x_i & & & \geq & \max_{s \in S} \sum_{j \in J} \lambda_j^s & \ & & & & & & \ & x_i & \in {0, 1} & & & & \forall i \in I \ & y_{ij}^s & \geq 0 & & & & \forall i \in I, \forall j \in J, \forall s \in S
\end{array} \tag{12} \end{equation} $

The sub-problem

Given the above formulation for the extensive form (12) we can express our scenario specific sub-problems as follows

$ \begin{equation} \begin{array}{lcrllll} z_s(\hat{x}) & = & \min \quad & \sum_{i \in I} \sum_{j \in J} q_{ij} y_{ij}^s & & & \ & & & & & & \ & & \textrm{s.t.} \quad & \sum_{i \in I} y_{ij}^s & \geq & \lambda_j^s & \forall j \in J \ & & & \sum_{j \in J} y_{ij}^s & \leq & k_i \hat{x_i} & \forall i \in I \ & & & & & & \ & & & y_{ij}^s \geq 0 & & & \forall ; i \in I, \forall j \in J \end{array} \tag{13} \end{equation} $

The master-problem

The master-problem takes the following form

$ \begin{equation} \begin{array}{rllll} \min \quad & \sum_{i \in I} f_i x_i + \sum_{s \in S} \alpha^s z^s & & & \ & & & & \ \textrm{s.t.} \quad & \sum_{i \in I} k_i x_i & \geq & \max_{s \in S} \sum_{j \in J} \lambda_j^s & \ & \sum_{j \in J} p_j^s \lambda_j^s + \sum_{i \in I} p_i^s k_i x_i & \leq & z^s & \forall s \in S, \forall p \in P^s \ & \sum_{j \in J} w_j^s \lambda_j^s + \sum_{i \in I} w_i^s k_i x_i & \leq & 0 & \forall s \in S, \forall w \in P^s \ & & & & \ & x_i \in {0, 1} & & & \forall i \in I \ \end{array} \tag{12} \end{equation} $

Create amplpy instances of master and sub-problems

  • create_common(): this function initializes a base AMPL object. 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 the AMPL instance for the sub-problem. It begins by invoking create_common(). Following this, it employs the eval() 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 with create_common() to instantiate an AMPL object with the basic shared model elements. It then utilizes the eval() 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"
        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()
        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()
        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


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 (

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 (

Finally, we write a function that adds the Benders cuts using ampls; the optimality and feasibility cuts are very similar to each other: $ \begin{equation} \begin{array}{rllll} & & & & \ \sum_{j \in J} p_j^s \lambda_j^s + \sum_{i \in I} p_i^s k_i x_i & \leq & z^s & \forall s \in S, \forall p \in P^s \ \sum_{j \in J} w_j^s \lambda_j^s + \sum_{i \in I} w_i^s k_i x_i & \leq & 0 & \forall s \in S, \forall w \in P^s \ & & & & \ \end{array} \end{equation} $

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":
    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
    print(f"Added {cutype} cut: {opt.toString()}\n")

Benders decomposition

The main function is, at this point, easy to implement.

  1. 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()


# 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
) = 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
  1. We use the function AMPL.to_ampls to export an instance of the master problem to the ampls library. Note that at this point, solver_master is 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
) = get_maps(solver_master)

Now the master loop, which implements Benders. At every algorithm iteration, we:

  1. Set the scenario index in the subproblem and solve it. We check the solution:

    1. Infeasible -> We extract the dual rays and use them to construct a feasibility cut

    2. Feasible but violating the optimality condition -> We extract the dual values and use them to construct an optimality cut

    3. Feasible and not violating -> We mark this scenario as non-violating

  2. If no scenarios had violations, we reached optimality and quit

  3. Otherwise, we resolve the master problem (with the cuts added in step 1)

  4. We use the solution to assign:

    1. the parameter sub_facility_open in the subproblems

    2. the dictionary sub_variable_cost used 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
        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()
        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()
            n_noviolations += 1

    # If no scenario had a violation, we exit
    if n_noviolations == len(scenarios):


    # 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;

****** Iteration 2 *******
Version identifier: | 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
Version identifier: | 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:


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_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'] >= 

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()
        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_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 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!