AMPL Development Tutorial 5/6 – Parallelizing Subproblem Solves in Benders Decomposition
Description: In the fifth installment of our six-part series, we delve deeper by showing how to evolve our Benders decomposition Python script from a serial execution to one that solves subproblems in parallel.
This transition has the potential to significantly reduce solution times, especially when dealing with complex subproblems or a large number of scenarios that generate many subproblems.
In our concluding notebook, we will introduce the AMPLS library, a robust resource designed to further enhance Benders decomposition’s efficiency.
It achieves this by enabling the export of AMPL model instances as persistent solver representations, thus streamlining the solution process.
Tags: amplpy, ampl, mip, stochastic, facility location, benders, decomposition, parallel solves
Notebook author: 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)
Introduction
To ensure a cohesive understanding and for ease of reference, we will revisit the foundational concepts of stochastic programming and the theory behind Benders decomposition in this notebook.
This repetition will help reinforce these critical concepts and provide a seamless transition into the PYTHON script implementing Benders decomposition of our stochastic facility location problem.
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
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}
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}
Let
\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$.
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).
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}
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$:
\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}
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:
\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.
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
Given:
Task:
Variables
$x_i \in {0, 1} \quad \forall i \in I$
$y_{ij}^s \geq 0 \quad \forall i \in I, \forall j \in J, \forall s \in 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 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}
$
AMPL implementation
Converting the mathematical formulation of our decomposition scheme into an AMPL model is straightforward.
Similar to the extensive form, AMPL enables the modeling of problems in a way that closely aligns with their conceptual structure.
In AMPL, sets and parameters are globally accessible, allowing for seamless integration across different problem components.
AMPL facilitates the assignment of variables, objective functions, and constraints to specific named problems using the syntax problem <problemName>: <list of model entities>.
In the model file presented below, we’ve organized the content to clearly indicate which elements pertain to the master problem and which to the subproblems.
This organization is primarily for clarity and ease of understanding; the AMPL interpreter treats the entire file as a single problem until it encounters explicit problem statements.
Beyond the entities familiar from the extensive form, we introduce additional parameters like sub_facility_open, which convey the master problem’s solutions to the subproblems.
We also incorporate parameters for the dual values from the subproblems, namely customer_price and facility_price.
New constraints, often referred to as cuts, are added to the master problem based on the subproblems’ feasibility and optimality relative to the master problem’s objectives.
Furthermore, a new set of variables, sub_variable_cost, is added to the master problem. These variables account for the variable costs identified in each scenario’s subproblem, contributing to the overall total cost calculation.
Overwriting floc_bend.mod
Data as Python objects
In this notebook, we continue to work with the same dataset as in the previous installments.
However, the format in which we store our data differs; we now utilize Python lists and Pandas data frames. Initially, we create these Python objects and subsequently employ amplpy methods to seamlessly transfer the data into AMPL.
Solving subproblems in parallel
Leveraging multiple processor cores enables us to solve subproblems in parallel, enhancing efficiency by allowing each subproblem to be addressed independently.
To facilitate this, we’ll develop a dedicated worker function specifically tailored for solving these subproblems.
Each subproblem will be processed in its own AMPL process, with the solutions integrated back into the master problem to update its parameters.
Noteworthy Changes
In contrast to our approach in previous notebooks, where a single AMPL process managed both the master problem and all subproblems, our new setup diverges significantly.
We will initiate separate AMPL processes for each subproblem, in addition to one for the master problem, resulting in a total of “number of scenarios plus one” AMPL processes.
This change precludes the possibility of direct data transfers between subproblems and the master problem using AMPL’s let commands, such as:
let {j in CUSTOMERS} customer_price[j, sub_scenario, nITER] := satisfying_customer_demand[j].dual;
let {i in FACILITIES} facility_price[i, sub_scenario, nITER] := facility_capacity_limits[i].dual;
Instead, data transfer between the subproblems and the master problem will be managed through Python.
This requires the master problem’s Python process to extract data from the subproblems and then load it into the master problem’s AMPL process, and vice versa.
To manage these data transformations efficiently, additional Python code is necessary.
A practical solution is to develop a Python module comprised of utility and worker functions tailored for these tasks.
Given our Jupyter Notebook environment, we will define these functions in a separate code cell for ease of access and organization.
Utility and worker Functions
To streamline our approach before tackling the core of our task, we will establish a series of utility functions designed to handle data transformations and to reduce redundancy in our code.
Alongside these, we’ll introduce a pivotal worker function, solve_sub(), designated to be run by individual threads in our parallel setup.
Worker functions typically have access to global symbols and data.
However, in the realm of parallel computing, it’s advisable to directly supply data to the worker function.
This practice promotes clarity and ensures that each operation is distinct and self-contained.
Adhering to this principle, the solve_sub() function will use a specific AMPL process for each scenario, handling the necessary computations for each subproblem independently.
Parallel Benders
Key elements of our parallel Benders implementation include:
Data Management: We employ dictionaries for each data type (Python lists, Pandas data frames) to streamline and clarify data loading operations, reducing code redundancy.
Core Availability Check: We assess the number of available processor cores. If the number of scenarios exceeds the available cores, we issue a warning to the user.
AMPL Object Initialization: For each scenario, we initialize an ampl_dict object, each containing an AMPL object configured to run as an independent process.
Subproblem Executor: We establish sub_executor, a dedicated object for managing the execution of subproblems.
While Loop Execution: The iterative process begins similarly to previous implementations, but with the solve operations are conducted outside the for loop in parallel.
Violation Checks and Cuts: Within the for loop, we simply identify violations and append the necessary cuts to the master problem.
Cleanup Operations: Upon exiting the while loop, we undertake cleanup operations to prevent the persistence of unused processes, ensuring system resources are efficiently managed.
This structured approach allows us to leverage parallel computing capabilities, potentially enhancing the efficiency and speed of the Benders decomposition process.
option solver highs;
ITERATION 1
HiGHS 1.6.0: pre:solve = off
alg:simplex = 1
HiGHS 1.6.0: HiGHS 1.6.0: pre:solve = off
pre:solve = off
alg:simplex = 1
HiGHS 1.6.0: optimal solution; objective 11621793.46
4 simplex iterations
0 barrier iterations
alg:simplex = 1
HiGHS 1.6.0: optimal solution; objective 21050711.2
6 simplex iterations
0 barrier iterations
HiGHS 1.6.0: optimal solution; objective 14779784.86
5 simplex iterations
0 barrier iterations
1: Optimality cut added for scenario Low
1: Optimality cut added for scenario Medium
1: Optimality cut added for scenario High
SOLVING THE MASTER PROBLEM
HiGHS 1.6.0: pre:solve = off
alg:simplex = 1
HiGHS 1.6.0: optimal solution; objective 16688018.6
0 simplex iterations
1 branching nodes
ITERATION 2
HiGHS 1.6.0: HiGHS 1.6.0: HiGHS 1.6.0: pre:solve = off
alg:simplex = 1
pre:solve = off
pre:solve = off
alg:simplex = 1
alg:simplex = 1
HiGHS 1.6.0: optimal solution; objective 11621793.46
0 simplex iterations
0 barrier iterations
HiGHS 1.6.0: optimal solution; objective 14966984.87
1 simplex iterations
0 barrier iterations
HiGHS 1.6.0: optimal solution; objective 21570711.2
0 simplex iterations
0 barrier iterations
2: No cut needed for scenario Low
2: Optimality cut added for scenario Medium
2: No cut needed for scenario High
SOLVING THE MASTER PROBLEM
HiGHS 1.6.0: pre:solve = off
alg:simplex = 1
HiGHS 1.6.0: optimal solution; objective 16758018.6
1 simplex iterations
1 branching nodes
ITERATION 3
HiGHS 1.6.0: pre:solve = off
HiGHS 1.6.0: alg:simplex = 1
HiGHS 1.6.0: pre:solve = off
pre:solve = off
alg:simplex = 1
alg:simplex = 1
HiGHS 1.6.0: optimal solution; objective 11621793.46
0 simplex iterations
0 barrier iterations
HiGHS 1.6.0: optimal solution; objective 14779784.87
1 simplex iterations
0 barrier iterations
HiGHS 1.6.0: optimal solution; objective 21050711.2
0 simplex iterations
0 barrier iterations
3: No cut needed for scenario Low
3: No cut needed for scenario Medium
3: No cut needed for scenario High
OPTIMAL SOLUTION FOUND
Expected cost = 16758018.596250001
facility_open [*] :=
Baton_Rouge_LA 1
Baytown_TX 1
Beaumont_TX 1
;
Conclusion
In this section, we’ve illustrated the process of integrating parallel computations into the Benders decomposition algorithm, with only minimal modifications required from our original serial Python implementation.
As the complexity of our code increases, we’ve adopted certain practices to enhance the modularity and robustness of our application.
Typically, this stage would involve the development of Python modules and the implementation of comprehensive exception handling to prepare the application for deployment.
However, a detailed exploration of these advanced topics falls beyond the scope of this series.
While our current implementation efficiently addresses subproblems, each solver instance still begins its computation from scratch.
In the concluding notebook of this series, we will explore the use of the AMPLS library to maintain persistent solver instances, thereby further improving the efficiency of our solution process.