A Review on the Performance of Linear and Mixed Integer Two-Stage Stochastic Programming Software

Abstract

This paper presents a tutorial on the state-of-the-art software for the solution of two-stage (mixed-integer) linear stochastic programs and provides a list of software designed for this purpose. The methodologies are classified according to the decomposition alternatives and the types of the variables in the problem. We review the fundamentals of Benders decomposition, dual decomposition and progressive hedging, as well as possible improvements and variants. We also present extensive numerical results to underline the properties and performance of each algorithm using software implementations, including DECIS, FORTSP, PySP, and DSP. Finally, we discuss the strengths and weaknesses of each methodology and propose future research directions.

1. Introduction

In the modeling and optimization of real-world problems, there is usually a level of uncertainty associated with the input parameters and their future outcomes. Stochastic programming (SP) models have been widely studied to solve optimization problems under uncertainty over the past decades [1,2]. SP is acknowledged for providing superior results compared to the corresponding deterministic model with nominal values for the uncertain parameters, which can lead to suboptimal or infeasible solutions. SP applications in process systems engineering include manufacturing networks and supply chain optimization [3,4], production scheduling [5], synthesis of process networks [6].

Two-stage stochastic programming is a commonly applied framework for cases where parameter uncertainties are decision independent (exogenous uncertainties). In stochastic programming, uncertain parameters are explicitly represented by a set of scenarios. Each scenario corresponds to one possible realization of the uncertain parameters, according to a discretized probability distribution. The goal is to optimize the expected value of the objective function over the full set of scenarios, subject to the implementation of common decisions at the beginning of the planning horizon.

Stochastic programs are often difficult to solve due to their large size and complexity that grows with the number of scenarios. To overcome these difficulties, decomposition algorithms, such as Benders decomposition [7], Lagrangian decomposition [8], and progressive hedging [9], have been developed to solve linear programming (LP) and mixed-integer linear programming (MILP) stochastic problems. For a comprehensive review of recent algorithmic advances in two-stage stochastic MIP, we refer the readers to the tutorial [10].

On the software development front, several modeling systems and optimization platforms have included extensions for a convenient algebraic representation of stochastic problems, as offered by major software vendors, such as GAMS, LINDO, XpressMP, AIMMS, and Maximal.

In recent years, various commercial and open-source applications have been developed specifically to represent and solve two-stage and multistage SP problems. Some of them include capabilities to read and build stochastic MPS (SMPS) files, the standard exchange format for SP applications. However, despite advances in the field and proven benefits, SP has not been widely used in industrial applications. Existing applications include power systems planning and scheduling [11,12,13], and process systems engineering [14,15].

The scope of this paper is to provide a tutorial for beginners on linear and mixed-integer linear stochastic programming, introducing the fundamentals on solution methodologies with a focus on available software. To accomplish this objective, we present a general description of techniques based on Benders decomposition and dual decomposition, which are the fundamentals of most known solvers designed specifically to tackle special problem structures. We then review the currently available state-of-the-art software for the solution of two-stage stochastic programming problems and evaluate their performance, using large-scale test libraries in SMPS format.

The remainder of this paper is organized as follows. In Section 2, we explain the mathematical formulation of (mixed-integer) linear stochastic problems. Section 3 describes the classical L-shaped algorithm. Section 4 summarizes the enhancement strategies to improve the performance of Benders decomposition. Section 5 describes scenario decomposition methods and algorithmic enhancements. In Section 6, we describe algorithmic innovations in software packages for dual decomposition and show some computational results. Finally, in Section 7, we summarize the conclusions.

2. Problem Statement

We consider a two-stage stochastic mixed-integer linear problem (P) in the following form:

$$\left(P\right)\underset{x,y}{min}TC=c^{T}x+\sum _{\omega \in \mathrm{\Omega }}{\tau }_{\omega }{d}_{\omega }^T{y}_{\omega }$$

(1a)

$$\mathrm{s}.\mathrm{t}.Ax\le b$$

(1b)

$$x\in X,X=\left\{x:x_i\in \{0,1\}\forall i\in I_1,x_i\ge 0\forall i\in I\backslash I_1\right\}$$

(1c)

$$W_{\omega }{y}_{\omega }\le h_{\omega }-T_{\omega }x\forall \omega \in \mathrm{\Omega }$$

(1d)

$$y_{\omega }\in Y\forall \omega \in \mathrm{\Omega }$$

(1e)

where x denotes the ‘here and now’ decisions, taken at the beginning of the planning horizon before the uncertainties unfold, and $\mathrm{\Omega }$ is the set of scenarios. Vector $y_{\omega }$ represents the recourse or corrective actions (wait and see), applied after the realization of the uncertainty. Matrix $A\in {\mathbb{R}}^{m_1\times n_1}$ and vector $b\in {\mathbb{R}}^{m_1}$ represent the first-stage constraints. Matrices $T_{\omega }$ and $W_{\omega }$, and vector $h_{\omega }\in {\mathbb{R}}^{m_2}$ represent the second-stage problem. Matrices $T_{\omega }\in {\mathbb{R}}^{m_2\times n_1}$ and $W_{\omega }\in {\mathbb{R}}^{m_2\times n_2}$ are called technological and recourse matrices, respectively. Let $I=\{1,2,\dots ,n_1\}$ be the index set of all first-stage variables. Set $I_1\subseteq I$ as the subset of indices for binary first-stage variables. Let $J=\{1,2,\dots ,n_2\}$ be the index set of all second-stage variables. If the second-stage variables are mixed integer, $Y=\left\{y:y_j\in \{0,1\},\forall j\in J_1,y_j\ge 0\forall j\in J\backslash J_1\right\}$, where set $J_1\subseteq J$ is the subset of indices for binary second-stage variables. If all the second-stage variables are continuous, set $J_1=\varnothing$ and $Y=\left\{y:y_j\ge 0\forall j\in J\right\}$. The objective function ($TC$) minimizes the total expected cost with the scenario probability ${\tau }_{\omega }$, and the cost vectors c and $d_{\omega }$. It should be noted that the approaches discussed in this paper can be applied to mixed-integer variables. Here, we restrict the variables to be mixed binary for simplicity. Equation (1) is often referred to as the deterministic equivalent, or extensive form of the SP.

Formulation (P) can be rewritten in an equivalent form (PNAC) with nonanticipativity constraints (NACs), where the first-stage variables are no longer shared, and each scenario represents an instance of a deterministic problem with a specific realization outcome [16,17].

$$\left(PNAC\right)\underset{x_{\omega },y_{\omega }}{min}TC=\sum _{\omega \in \mathrm{\Omega }}{\tau }_{\omega }(c^T{x}_{\omega }+d_{\omega }^T{y}_{\omega })$$

(2a)

$$\mathrm{s}.\mathrm{t}.\sum _{\omega \in \mathrm{\Omega }}{H}_{\omega }{x}_{\omega }=0$$

(2b)

$$(x_{\omega },y_{\omega })\in G_{\omega }\forall \omega \in \mathrm{\Omega }$$

(2c)

In Equation (2c), $G_{\omega }$ represents the feasible region for scenario $\omega$, which is defined by constraints (1b)–(1e). Nonanticipativity constraints (2b) are added to ensure that the first-stage decisions are the same across all scenarios. Nonanticipativity constraints are represented by a suitable sequence of matrices $H_{\omega }\in {\mathrm{R}}^{n_1·\left(\right|\mathrm{\Omega }|-1)\times n_1}$. One example of such constraints is the following:

$$\begin{array}{c}{x}_{\omega }=x_{\omega -1}\forall \omega =2,3,\dots ,\left|\mathrm{\Omega }\right|\end{array}$$

(3)

Given the mathematical structure of the deterministic equivalent formulations (P) and (PNAC), (mixed-integer) linear stochastic problems can be solved using decomposition methods derived from duality theory [18]. Such methods split the deterministic equivalent into a master problem and a series of smaller subproblems to decentralize the overall computational burden. Decomposition methodologies are classified in two groups: (i) node-based or vertical decomposition, which includes Benders decomposition and variants where the problem is decomposed by the nodes in the scenario tree, and (ii) scenario-based or horizontal decomposition, where the problem is decomposed by scenarios. In the following Section 3, we provide a tutorial overview of Benders decomposition. In Section 4, we also provide a tutorial review of scenario decomposition methods, including the dual decomposition algorithm and the progressive hedging algorithm.

3. L-Shaped Algorithm/Benders Decomposition

If the second-stage variables are all continuous (i.e., $Y = \left\{y:y_j\ge 0\forall j\in J\right\}$), problem (P) can be solved with Benders decomposition. Benders decomposition (BD) was originally developed in 1962 by Benders [19] to solve large-scale mixed-integer linear problems (MILP) with complicating variables. This concept has been extended to solve a broader range of optimization problems [20], including multistage, bilevel, and nonlinear programming. When applied to stochastic problems, it is commonly referred to as the L-shaped algorithm [7].

The L-shaped algorithm partitions the deterministic formulation (P) into multiple problems: (i) a master problem (MP) that contains all the first-stage constraints and variables, which can be mixed integer; and (ii) a collection of subproblems that include corrective future actions for the given first-stage solution. The master problem (MP) is derived from the projection of (P) on the variables x:

$$\left(MP\right)\underset{x}{min}TC=c^{T}x+Q\left(x\right)$$

(4a)

$$\mathrm{s}.\mathrm{t}.Ax\le b$$

(4b)

$$x\in X$$

(4c)

where $Q\left(x\right)={\sum }_{\omega \in \mathrm{\Omega }}{\tau }_{\omega }{\theta }_{\omega }\left(x\right)$ is defined as the recourse function (or expected second-stage value function); and ${\theta }_{\omega }\left(x\right)$ is defined by the primal second-stage program for scenario $\omega$, $\left(BSP{p}_{\omega }\right)$,

$$\left(BSP{p}_{\omega }\right){\theta }_{\omega }\left(x\right)=\underset{y_{\omega }}{min}{d}_{\omega }^T{y}_{\omega }$$

(5a)

$$\mathrm{s}.\mathrm{t}{W}_{\omega }{y}_{\omega }\le h_{\omega }-T_{\omega }x$$

(5b)

$$y_{\omega }\ge 0$$

(5c)

The recourse functions ${\theta }_{\omega }\left(x\right)$ and $Q\left(x\right)$ are convex, differentiable, and piece-wise linear, characteristics that are exploited in the BD method [1]. These conditions do not hold when integer variables are included in the second-stage program. For the case of integer recourse, a logic-based Benders framework [21], second-stage convexification techniques [22,23,24], specialized branch-and-bound schemes [25,26] or dual decomposition methods [16] may be applied to solve large stochastic problems. In this section, we only focus on Benders decomposition for SP with continuous second-stage variables.

Formulation ($BSP{p}_{\omega }$) is a linear program for any given feasible value of x. By the strong duality theorem, the second-stage program is equivalent to its dual ($BSP{d}_{\omega }$) if ($BSP{p}_{\omega }$) is feasible and bounded. Vector ${\pi }_{\omega }$ represents the Lagrangian multipliers associated with the second-stage constraints given by Equation (5b):

$$\left(BSP{d}_{\omega }\right){\theta }_{\omega }\left(x\right)=\underset{{\pi }_{\omega }}{max}{(h_{\omega }-T_{\omega }x)}^T{\pi }_{\omega }$$

(6a)

$$\mathrm{s}.\mathrm{t}.W_{\omega }^T{\pi }_{\omega }\le d_{\omega }$$

(6b)

$${\pi }_{\omega }\ge 0$$

(6c)

BD introduces a set of piece-wise linear approximations of the recourse function in the problem MP, known as optimality cuts, which are built from dual solutions of the second-stage program. It is important to highlight that the dual feasible region does not depend on the value of x. Thus, the exact representation of the expected cost consists of the computation of all the extreme points of problems ($BSP{d}_{\omega }$).

However, the second-stage program may not be feasible for some values of x. BD enforces second-stage constraints (5b) by adding feasibility cuts, which are valid inequalities that exclude infeasible first-stage solutions from the MP. Subproblem feasibility is evaluated by solving the following recourse reformulation for scenario $\omega$:

$$V_{\omega }\left(x\right)=\underset{y_{\omega },v^{-}}{min}{e}^T{v}^{-}$$

(7a)

$$\mathrm{s}.\mathrm{t}.W_{\omega }{y}_{\omega }-v^{-}\le h_{\omega }-T_{\omega }x$$

(7b)

$$v^{-}\ge 0,y_{\omega }\ge 0$$

(7c)

where $e\in {\mathbb{R}}^{m_2}$ is a vector with all-1 entries, and $v^{-}\in {\mathbb{R}}^{m_2}$ is the slack of constraints (5b). The objective function $V_{\omega }\left(x\right)$ measures the amount by which these constraints are violated; thus, if $V_{\omega }\left(x\right)$ equals zero, it implies that the original subproblem (5) is feasible. To derive feasibility cuts in terms of x, BD considers the dual of problem (7) to generate an expression equivalent to Equation (7a). The optimal solution $\mu \in {\mathbb{R}}^{m_2}$ of the dual feasibility problem (8) corresponds to one of the extreme rays (or directions) of the recourse subproblem (6):

$$V_{\omega }\left(x\right)=\underset{\mu }{max}{(h_{\omega }-T_{\omega }x)}^T\mu$$

(8a)

$$\mathrm{s}.\mathrm{t}.W_{\omega }^T\mu \le 0$$

(8b)

$$0\le \mu \le e$$

(8c)

The master problem (4) is linearized by (i) substituting function $Q\left(x\right)$ with the weighted sum of the future cost estimation (6a), and (ii) applying feasibility cuts as needed. This reformulation is referred to as the multi-cut Benders master problem (BMP),

$$\left(BMP\right)T{C}_d=\underset{x,{\theta }_{\omega }}{min}{c}^{T}x+\sum _{\omega \in \mathrm{\Omega }}{\tau }_{\omega }{\theta }_{\omega }$$

(9a)

$$\mathrm{s}.\mathrm{t}.Ax\le b,x\in X$$

(9b)

$${(h_j-T_{j}x)}^T{\overline{\mu }}_j\le 0\forall j\in E$$

(9c)

$${(h_{\omega }-T_{\omega }x)}^T{\overline{\pi }}_{\omega }^k\le {\theta }_{\omega }\forall \omega \in \mathrm{\Omega },k\in K$$

(9d)

where variables ${\theta }_{\omega }\in {\mathbb{R}}^{\left|\mathrm{\Omega }\right|}$ represent the outer linearization of the second-stage cost ${\theta }_{\omega }\left(x\right)$. Parameters ${\overline{\pi }}_{\omega }^k$ and ${\overline{\mu }}_j$ represent the extreme points and rays from the dual form of the recourse program ($BSP{d}_{\omega }$), which are stored in sets E and K, respectively. Constraints (9c) and (9d) denote the feasibility and optimality cuts, $j\in E$ and $k\in K$ respectively. Matrices $h_j$ and $T_j$ correspond to the matrices $h_{\omega }$ and $T_{\omega }$ for the scenario where a feasibility cut can be found.

The complete enumeration of the extreme points and rays of the dual second-stage program is impractical, if not impossible. Instead, the L-shaped algorithm relaxes the MP by initially considering a subset of the optimality and feasibility cuts. BD iteratively solves the relaxed problem to generate a candidate solution for the first-stage variables ($\overline{x}$), and then solves the collection of scenarios subproblems at fixed $\overline{x}$ to generate a new group of optimality or feasibility cuts. This process is repeated until the optimal solution is found [26].

The optimal solution of the relaxed Benders master problem provides a valid lower estimation ($T{C}_d$) of the optimal total cost, TC. On the other hand, the solution of the second-stage programs ($BSP{d}_{\omega }$) at feasible $\overline{x}$ yields an upper bound of the original objective function ($T{C}_p$), given by Equation (10). The solution procedure terminates when the difference between the bounds is closed, as implied by Equation (11). Algorithm 1 summarizes the procedure.

$$T{C}_p\left(\overline{x}\right)=c^T\overline{x}+\sum _{\omega \in \mathrm{\Omega }}{\tau }_{\omega }{\theta }_{\omega }\left(\overline{x}\right)$$

(10)

$$T{C}_d\le TC\le T{C}_p$$

(11)

Algorithm 1: Multi-cut Benders decomposition

The L-shaped method is summarized in Algorithm 1. It is initialized with a guess of the first-stage solution $x_o$ and considers two stopping criteria: (i) the optimality tolerance $ϵ$ that limits the relative gap between the dual $\left(z_{LB}\right)$ and primal $\left(z_{UB}\right)$ bounds of the objective function $\left(TC\right)$, and (ii) the maximum number of allowed iterations $\left(k_{max}\right)$.

4. Benders Decomposition Enhancement Strategies

The application of BD often leads to slow convergence, long computational times, and excessive use of memory resources, particularly for the case when the MILP master problem has poor LP relaxation [27,28,29]. Major BD disadvantages include time-consuming iterations, poor feasibility and optimality cuts, ineffective initial iterations, primal solutions that behave erratically, slow convergence at the end of the algorithm (tailing-off effect), and upper bounds that remain stuck in successive iterations due to second-stage degeneracy [26,30].

Various strategies have been proposed to accelerate the convergence of the standard BD method. Enhancement strategies are mainly split into two categories: reducing the cost of each iteration, or reducing the number of iterations [26,31].

4.1. Reducing the Cost of Each Iteration

Cheaper iterations are achieved by reducing the time spent solving the MP and subproblems. The MP is often the most time-consuming part of the BD algorithm (more than 90% of the execution time) in the case of mixed-integer problems [27]. The overall process can be accelerated by relaxing the integrality of the MP in most of the iterations to rapidly compute a large number of cuts [32]. A variation of this method was proposed by Geoffrion and Graves [28], in which the MP is solved to a non-zero optimality gap. The integrality gap is continuously reduced to ensure global convergence.

Alternatively, the MP might be solved via (meta) heuristics [33,34], which provide good approximate solutions in short time; however, it is still required to solve the MP to optimality to guarantee convergence. The application of heuristics or the LP relaxation of the MP often yields worse bounds and lack of controllability, reducing the ability of BD to generate the necessary cuts [35].

Similarly, suboptimal solutions of the dual subproblems yield valid cuts, known as inexact cuts. Algorithm convergence can still be guaranteed under the conditions described by Zakeri et al. [36]. Additional subproblem acceleration schemes comprise synchronous parallelization and re-optimization. The latter exploits structural similarities between scenarios to solve the subproblems in fewer solver iterations.

4.2. Reducing the Number of Iterations

The number of iterations of the L-shaped algorithm is closely related to the tightness of the LP relaxation of the first-stage problem, as well as the strength of the optimality and feasibility cuts [27]. Better candidates are computed from improvements in the MP problem, especially in the strengthening of the representation of the recourse functions. Tighter formulations can be obtained by adding multiple cuts per iteration (multi-cut reformulation [37]); as well as through the use of heuristics to eliminate inactive cuts and to select the fittest dual variables to be inserted in the MP (size management techniques).

Complementary strategies have been developed to generate cuts that are more efficient. One alternative is the reformulation of the subproblems to select non-dominant dual solutions from the set of optimal multipliers, known as Pareto-optimal cuts [27]. Recently, Ref. [38] proposed a methodology to compute bundles of covering cuts, designed to involve most of the first-stage variables and to carry as much information as possible.

Alternative methods tighten the MP to alleviate some of the drawbacks of BD: cross decomposition, for instance, avoids the generation of low-quality solutions, while quadratic stabilization methods provide a solution for the tailing-off effect. Cross decomposition [39] combines and alternates between iterations of BD and Lagrangian decomposition to provide an additional valid lower bound of the objective function and a set of feasible deterministic solutions $(x_{\omega },y_{\omega })\in G_{\omega }$, which are used to compute Lagrangian-based cuts to strengthen the MP.

Quadratic methods have been proposed to stabilize BD, aiming to improve the quality of the initial iterations and reduce the oscillation that occurs when the algorithm is close to the optimal solution [31]. These methods penalize the distance of the first-stage candidates to the current best solution. Popular variants include regularized decomposition [16,40], the trust-region method [41] and level decomposition [42], which are summarized below. It should be noted that Equations (12)–(14) are all variants of the large family of bundle method, which can be seen from chapters XIV and XV of [43].

• Regularized decomposition (also known as proximal bundle method)

  $$x_{k+1}=\underset{x,{\theta }_{\omega }}{arg\; min}\{c^{T}x+\sum _{\omega \in \mathrm{\Omega }}{\tau }_{\omega }{\theta }_{\omega }+\frac{1}{2t_k}{\parallel x-{\widehat{x}}_k\parallel }_2^2\mathrm{s}.\mathrm{t}.(x,{\theta }_{\omega })\in V_k\}$$

  (12)
• Trust-region method

  $$x_{k+1}=\underset{x,{\theta }_{\omega }}{arg\; min}\{c^{T}x+\sum _{\omega \in \mathrm{\Omega }}{\tau }_{\omega }{\theta }_{\omega }\mathrm{s}.\mathrm{t}.{\parallel x-{\widehat{x}}_k\parallel }_2^2\le R_k,(x,{\theta }_{\omega })\in V_k\}$$

  (13)
• Level decomposition method

  $$x_{k+1}=\underset{x,{\theta }_{\omega }}{arg\; min}\{{\parallel x-{\widehat{x}}_k\parallel }_2^2\mathrm{s}.\mathrm{t}.c^{T}x+\sum _{\omega \in \mathrm{\Omega }}{\tau }_{\omega }{\theta }_{\omega }\le L_k,(x,{\theta }_{\omega })\in V_k\}$$

  (14)

  where $t_k$, $R_k$ and $L_k$ are real-valued, iteration-dependent parameters that balance the minimization of the relaxed MP and the distance to the stability center ${\widehat{x}}_k$. $V_k$ represents the feasible region of the Benders master problem at each iteration, which is defined by the optimality (9d) and feasibility cuts (9c), as well by the first-stage constraints (9b). Stabilization methods were initially introduced for BD with no integer variables; nonetheless, recent improvements have adapted the method to mixed-integer problems [31].

Two software packages using the Benders decomposition algorithm, GAMS– DECIS [44,45] and FORTSP [46] are described in the benchmark in Appendix C.

5. Scenario Decomposition Methods

Scenario decomposition is a popular approach to solve two-stage SP formulations with mixed-integer recourse, i.e., $Y=\left\{y:y_j\in \{0,1\},\forall j\in J_1,y_j\ge 0\forall j\in J\backslash J_1\right\}$ in (PNAC). In contrast to the BD algorithm, scenario decomposition methods dualize the non-anticipativity constraints (NACs) to obtain lower bounds of the original formulation. Scenario-based decomposition addresses the computational difficulties associated with the solution of large stochastic problems by considering each scenario independently and solving the set of subproblems in parallel. Moreover, feasible solutions to the original problem (P) can be obtained by heuristics based on the optimal solutions of the subproblems. In this section, we describe the dual decomposition (DD) algorithm and the progressive hedging (PH) algorithm.

5.1. Dual Decomposition (DD) Method

The dual decomposition algorithm proposed by Carøe and Schultz [17] applies the Lagrangian relaxation to problem (2) and uses a $branch$-$and$-$bound$ procedure to restore the non-anticipativity conditions. The Lagrangian relaxation of the NACs results in the following dual function:

$$D\left(\lambda \right)=\underset{x,y}{min}\sum _{\omega \in \mathrm{\Omega }}{L}_{\omega }(x_{\omega },y_{\omega },{\lambda }_{\omega })$$

(15a)

$$\mathrm{s}.\mathrm{t}.(x_{\omega },y_{\omega })\in G_{\omega }\forall \omega \in \mathrm{\Omega }$$

(15b)

where

$$L_{\omega }(x_{\omega },y_{\omega },{\lambda }_{\omega })={\tau }_{\omega }(c^T{x}_{\omega }+d_{\omega }^T{y}_{\omega })+{\lambda }_{\omega }^T{H}_{\omega }{x}_{\omega }$$

(16)

In Equation (15a), vector $\lambda \in {\mathbb{R}}^{n_1\times \left(\right|\mathrm{\Omega }|-1)}$ represents the dual multipliers associated with the NACs (2b). ${\lambda }_{\omega }\in {\mathbb{R}}^{n_1}$ represents the Lagrangian multipliers for the NACs associated with scenario $\omega$, as given by Equation (3). Given the independence of the variables and constraints in each scenario, function D can be split into separate subproblems $D_{\omega }\left({\lambda }_{\omega }\right)$:

$$D\left(\lambda \right)=\sum _{\omega \in \mathrm{\Omega }}{D}_{\omega }\left({\lambda }_{\omega }\right)$$

(17a)

$$D_{\omega }\left({\lambda }_{\omega }\right)=\{\underset{x_{\omega },y_{\omega }}{min}{L}_{\omega }(x_{\omega },y_{\omega },{\lambda }_{\omega })\mathrm{s}.\mathrm{t}.(x_{\omega },y_{\omega })\in G_{\omega }\}$$

(17b)

According to the weak duality theorem, the relaxation (17) is always less than or equal to the optimal objective value of problem (2). The best lower bound of (PNAC) is computed by solving the following maximization problem, referred to as the Lagrangian dual problem:

$$Z_{LD}=\underset{\lambda }{max}D\left(\lambda \right)$$

(18)

The Lagrangian dual is a concave non-smooth program and can be solved by subgradient methods, cutting-plane methods, or column-generation methods. The details of these methods can be found in Guignard [47]. We illustrate the dual search approaches by describing the standard cutting-plane algorithm.

5.1.1. Cutting-Plane Method

The cutting-plane algorithm solves the Lagrangian problem iteratively by implementing the outer approximation on (18) and solving the Lagrangian subproblems (17b) to improve the formulation of the relaxed dual function ($R{Z}_{LD}$) in Equation (19a). The outer approximation is given by the Lagrangian master problem (LMP):

$$\left(LMP\right)R{Z}_{LD}=\underset{{\lambda }_{\omega },{\varphi }_{\omega }}{max}\sum _{\omega \in \mathrm{\Omega }}{\varphi }_{\omega }$$

(19a)

$$\mathrm{s}.\mathrm{t}.{\varphi }_{\omega }\le \overline{D_{\omega }^k}\left({\lambda }_{\omega }^k\right)+{\left(H_{\omega }{x}_{\omega }^k\right)}^T({\lambda }_{\omega }-{\lambda }_{\omega }^k)\forall k\in K,\omega \in \mathrm{\Omega }$$

(19b)

where parameters for iteration k and scenario $\omega$, $x_{\omega }^k$ and $\overline{D_{\omega }^k}\left({\lambda }_{\omega }^k\right)$, represent the previous solution of subproblem (17b), and parameter ${\lambda }_{\omega }^k$ represents the vector of the previously considered dual multipliers. The dual search is outlined in Algorithm 2. It should be noted that (19) can be unbounded, especially in the first few iterations. An approach to avoid the unboundedness is to add bounds for the ${\lambda }_{\omega }$ variables. An alternative way is to use the bundle method [48] instead of Algorithm 2.    

Algorithm 2: Cutting-plane dual search

Cutting-plane methods present similar drawbacks as the BD algorithm, such as slow convergence and strong oscillations of the dual variables. This can be explained by the fact that BD is a special case of the cutting-plane method from the view of nonsmooth optimization. Various alternatives have been proposed to accelerate this technique, including the bundle method and the volume algorithm [47]. Additional strategies consider the suboptimal solution of the master problem, using an interior-point method (IPM) in combination with Benders-like cuts to tighten the Lagrangian subproblems (17b) and exclude infeasible first-stage solutions (see [29,49]). Other methodologies such as cross decomposition, exchange information with BD to compute additional cuts derived from feasible first-stage candidates [50].

5.1.2. Branch-and-Bound Method

The DD algorithm proposed by Carøe and Schultz [17] uses the bound $Z_{LD}$ to discard nodes from the first-stage search domain. Algorithm 3 summarizes the branch-and-bound procedure. The set $\mathcal{P}$ denotes the group of active problems and $T{C}^i$, the lower bound associated with program ${\mathcal{P}}_i\in \mathcal{P}$. Commonly, the Lagrangian dual problem yields first-stage solutions that differ in value from one scenario to another. For those instances, a candidate $\widehat{x}$ is estimated by applying a rounding heuristic on the average solution ${\sum }_{\omega \in \mathrm{\Omega }}{\tau }_{\omega }{x}_{\omega }^i$. Note that Algorithm 3 can be applied not only to problems with mixed-binary variables, but to problems with general mixed-integer variables as well. The branching steps assume that the integer variables can be nonbinary. More, recently Kim and Dandurand [51] proposed a new branching strategy based on Dantzig–Wolfe decomposition that allows reduction in the number of nodes and decreases the overall solution time.    

Algorithm 3: DD Branch and Bound method

5.2. Progressive Hedging (PH) Algorithm

The progressive hedging (PH) algorithm [9] is a popular approximation for solving multi-stage stochastic programs. Although it was initially proposed for convex stochastic problems, it has been successfully applied as a heuristic to solve mixed-integer stochastic programs [52,53].

To find a solution of problem (2), PH aggregates a new set of variables $\widehat{x}$ (also known as a first-stage policy) that replaces the NACs (2b). Then, it solves the reformulated program (20) using a specialized variant of the alternating direction method of multipliers (ADMM) [54,55]:

$$\underset{x_{\omega },y_{\omega },\widehat{x}}{min}TC=\sum _{\omega \in \mathrm{\Omega }}(c^T{x}_{\omega }+d_{\omega }^T{y}_{\omega })$$

(20a)

$$\mathrm{s}.\mathrm{t}.(x_{\omega },y_{\omega })\in G_{\omega },x_{\omega }=\widehat{x},\forall \omega \in \mathrm{\Omega },\widehat{x}\in X$$

(20b)

Related to dual decomposition, PH relaxes the non-anticipativity restrictions on the first stage. The augmented Lagrangian relaxation $L^{\rho }$ of constraints $x_{\omega }=\widehat{x}, \forall \omega \in \mathrm{\Omega }$ yields a lower bound $D\left(\lambda \right)$ of the original deterministic formulation (20). The best lower bound is estimated by solving the following problem:

$$TC\ge \underset{\lambda }{max}\{D\left(\lambda \right)\mathrm{s}.\mathrm{t}.\sum _{\omega \in \mathrm{\Omega }}{\tau }_{\omega }{\lambda }_{\omega }=0\}$$

(21a)

where

$$D\left(\lambda \right)=\underset{x,\widehat{x},y}{min}{L}^{\rho }(x,\widehat{x},y,\lambda )\mathrm{s}.\mathrm{t}.(x_{\omega },y_{\omega })\in G_{\omega }\forall \omega \in \mathrm{\Omega },\widehat{x}\in X$$

(21b)

$$L^{\rho }(x,y,\widehat{x},\lambda )=\sum _{\omega \in \mathrm{\Omega }}{\tau }_{\omega }{L}_{\omega }(x_{\omega },y_{\omega },\widehat{x},{\lambda }_{\omega })$$

(21c)

$$L_{\omega }(x_{\omega },y_{\omega },\widehat{x},{\lambda }_{\omega })=c^T{x}_{\omega }+d_{\omega }^T{y}_{\omega }+{\lambda }^T(x_{\omega }-\widehat{x})+\rho /2{\parallel x-\widehat{x}\parallel }_2^2$$

(21d)

and $\rho >0$ is a penalty parameter. Constraint ${\sum }_{\omega \in \mathrm{\Omega }}{\tau }_{\omega }{\lambda }_{\omega }=0$ is required to make $L^{\rho }$ bounded from below. To mitigate the computational difficulties of minimizing the augmented Lagrangian dual function (21b), PH decomposes the problem by scenarios. To achieve the complete separability of subproblems, Rockafellar and Wets [9] proposed to fix the first-stage policy temporarily, and repeatedly solve the program (22) to update the multipliers and the value of $\widehat{x}$:

$$\underset{x_{\omega },y_{\omega }}{min}\{c^T{x}_{\omega }+d_{\omega }^T{y}_{\omega }+{\lambda }^T{x}_{\omega }+\rho /2{\parallel x_{\omega }-\widehat{x}\parallel }_2^2\}$$

(22)

Algorithm 4 summarizes the procedure to solve the dual problem (21).    

Algorithm 4: Two-stage progressive hedging algorithm

The termination of the algorithm is achieved when the first-stage policy is non-anticipative. In the case of convex instances, ${\widehat{x}}^{k⟶\infty }$ is equivalent to the optimal solution of the deterministic formulation (2), and the convergence is guaranteed. These conditions do not hold for mixed-integer programs; however, a high-quality solution and upper bound can be computed from a non-convergent value of ${\{{\widehat{x}}^k>\}}_{k=k_{max}}$ and $T{C}_p{>({\widehat{x}}^k>)}_{k=k_{max}}$, respectively [52].

Recent investigations have focused on the improvement and acceleration of PH. Various studies identify the penalty term as a critical factor in the quality of the solution and the convergence rate: larger values of $\rho$ can accelerate the convergence but can lead to suboptimal solutions. On the other hand, lower values can improve the quality of the solutions and lower bounds, although with a very slow convergence rate [56]. Numerous alternatives have been developed to circumvent those problems, from $per$-$component$ and $cost$-$proportional$ heuristics [52], to the dynamic update of the penalty parameter [42,57].

A limitation in applying PH to stochastic mixed-integer programs is the lack of a lower bound to assess the quality of the computed solution. This disadvantage can be alleviated by estimating a valid lower bound from the non-convergent set of Lagrangian weights ${\lambda }_k$ [56], or by combining the Frank–Wolfe and PH methods [58]. These methodologies establish a relationship between the dual decomposition and progressive hedging, which has motivated the development of hybrid solution strategies (see [59]).

6. Software Packages for Scenario Decomposition

In this section, we review two software packages, PySP [52,60] and DSP [49], for scenario decomposition. The two software packages are benchmarked based on the problems in SIPLIB [61], including the SSLP [62], SSLPR [63], and DCAP [64] test problems.

The SSLP test set consists of 12 two-stage stochastic mixed-integer programs arising in stochastic server location problems (SSLPs). The base deterministic server location problem considers building servers in some potential locations to serve clients in given locations. The stochastic version of the server location problem considers different realizations of client locations. Each scenario represents a set of potential clients that do materialize. The decisions in SSLP are all binary variables. In the first stage, we decide whether a server is located at each given location. The second-stage (recourse) actions decide whether any given client is served by any given server. SSLPR (stochastic server location problem with random recourse) is an extension of SSLP. While SSLP assumes fixed demands for the clients, SSLPR considers the demands of the clients as uncertain parameters.

DCAP consists of 12 two-stage stochastic integer programs arising in dynamic capacity acquisition and allocation applications. The deterministic model considers a capacity expansion schedule over T time periods. In each time period, the amount of capacity expansion for each resource needs to be decided. There is a fixed and a variable cost for each capacity expansion. In each time period, each task must be assigned to one of the existing resources, which is represented by binary variables that decide whether a given task is assigned to a given resource. Since there are multiple periods, the stochastic version of this problem should, in principle, be formulated as a multi-stage stochastic program, which is difficult to solve. Ahmed and Garcia [64] proposed to approximate the multi-stage problem with a two-stage stochastic program in which the first-stage decisions are the capacity expansions. The second-stage decisions are the assignment decisions. The uncertainties include the processing requirement for each task and the cost of processing each task.

The sizes of all the test problems are shown in

Table 1. The sizes of the problems tested.

Name	Scenarios	First Stage			Second Stage
		Rows	Cols	Ints	Rows	Cols	Ints
sslp_5_25	50, 100	1	5	5	30	130	125
sslp_10_50	50, 100, 500, 1000	1	10	10	60	510	500
sslp_15_45	5, 10, 15	1	15	15	60	690	675
sslprf_5_25	100	1	5	5	30	130	125
sslprf_5_50	100	1	10	10	60	510	500
dcap 233	200, 300, 500	6	12	6	15	27	27
dcap 243	200, 300, 500	6	12	6	18	36	36
dcap 332	200, 300, 500	6	12	6	12	24	24
dcap 342	200, 300, 500	6	12	6	14	32	32

. The names of the SSLP and SSLPR instances are expressed in the form sslp(rf)_m_n, where m is the number of potential server locations, and n is the number of potential clients. Each instance is tested with a different number of scenarios. The size is expressed as the number of constraints (Rows), variables (Cols), and integer variables (Ints) in the first stage and the second stage per scenario. Note that the SSLP problems have pure binary first-stage variables, and the DCAP problems have mixed-binary first-stage variables. This difference affects the PH algorithm, which is discussed in detail later.

All of the test problems are available in the SMPS format; however, we implement an interpreter to make the format compatible with PySP. All of the tests are run on a server with an Intel Xeon CPU (24 cores) at 2.67 GHz and 128 GB of RAM. The whole set of instances is solved in a synchronous parallel manner to reduce the time of each iteration.

6.1. PySP: Pyomo Stochastic Programming

PySP is a software package implemented in the Python programming language using Pyomo [65] as the optimization modeling framework. PySP enables the user to solve stochastic programs with a specialized progressive hedging algorithm for stochastic mixed-integer programs. In order to use PySP, the user only needs to write a deterministic base model and define the scenario tree structure in Pyomo. With these inputs, PySP is able to apply the progressive hedging algorithm as an effective heuristic for obtaining feasible solutions to multi-stage stochastic programs. We highlight that starting from distribution 6.0, PySP will not be part of Pyomo’s source code; however, it will still be available and maintained. Furthermore, the PySP project will be continued as mpi-sppy, available in https://github.com/Pyomo/mpi-sppy, accessed on 20 January 2022.

6.1.1. Algorithmic Innovations in PySP

The innovations in PySP for multi-stage mixed-integer stochastic programs were described by Watson and Woodruff [52]. Here, we briefly paraphrase those innovations. First, instead of keeping a fixed $\rho$ value for all first-stage decisions in Algorithm 4, the authors proposed several variable-dependent $\rho$ strategies. The cost-proportional (CP) strategy sets $\rho \left(i\right)$ to be proportional to the cost parameter $c\left(i\right)$, i.e., $\rho \left(i\right)=\alpha c\left(i\right)$, where $\alpha$ is a constant multiplier for all first-stage variables i. The other variable-dependent $\rho$ strategy was denoted by SEP in [52], where the $\rho \left(i\right)$ for integer variables is calculated by

$$\begin{array}{c}\hfill \rho \left(i\right):=\frac{c\left(i\right)}{(x^{max}-x^{min}+1)}\end{array}$$

After PH iteration 0, for each variable x, $x^{max}={max}_{\omega \in \mathrm{\Omega }}{x}_{\omega }^0$ and $x^{min}={min}_{\omega \in \mathrm{\Omega }}{x}_{\omega }^0$. For continuous variables, the $\rho \left(i\right)$ is calculated with

$$\begin{array}{c}\hfill \rho \left(i\right):=\frac{c\left(i\right)}{max\left(\left({\sum }_{\omega \in \mathrm{\Omega }}{\tau }_{\omega }|x_{\omega }^0-{\widehat{x}}^0|\right),1\right)}\end{array}$$

where ${\widehat{x}}^0$ is the weighted average of $x_{\omega }^0$, i.e., ${\widehat{x}}^0={\sum }_{\omega \in \mathrm{\Omega }}{\tau }_{\omega }{x}_{\omega }^0$.

The authors also proposed some heuristics for accelerating convergence. One heuristic is called “variable fixing”. The values of some of the first-stage decisions $x_{\omega ,i}$ are fixed after they stay at a given value $z_i$ for a few iterations for all scenarios. In order to apply this heuristic, the authors introduced a lag parameter $\mu$. At a given PH iteration k, the value of $x_{\omega ,i}^k$ is fixed to $z_i$ for all subsequent iterations $l>k$, if $x_{\omega ,i}^{\left(k\right)}=z_i$ for all $\omega \in \mathrm{\Omega }$ and $m\in \{k-\mu |\mathrm{\Omega }|,\dots ,k\}$, such that $m\ge \mu \left|\mathrm{\Omega }\right|$. Additionally, the authors proposed another more aggressive variable fixing heuristic called “variable slamming”, where the $x_{\omega }^k$ is fixed if they are “sufficiently converged”, i.e., there can be some discrepancies for $x_{\omega }^k$ across all scenarios. In order to decide when variable slamming should be applied, the authors proposed several termination criteria based on the deviations of the first-stage variables.

In solving stochastic mixed-integer programs with PH, the cyclic behavior can be found in some instances. In order to detect the cyclic behavior, the authors proposed a strategy based on the values of the $u_{\omega }$ vectors, i.e., the weights associated with the first-stage decision variable $x_{\omega }$. The authors proposed a simple hashing scheme. Let hash value $h\left(i\right)={\sum }_{\omega \in \mathrm{\Omega }}{z}_{\omega }{u}_{\omega ,i}$, where $z_{\omega }$ is an integer hash weight for each scenario $\omega \in \mathrm{\Omega }$ when the PH is initialized. If equal hash weights are detected, they are interpreted as evidence for a potential cycle. Variable $x_i$ can be fixed if cyclic behaviors are found.

The heuristics, including variable fixing and slamming and cyclic behavior detection, are denoted as WW (Watson–Woodruff) heuristics in the software distribution of PySP.

6.1.2. Computational Results for PySP

We use PySP (Pyomo 5.0.0) to solve the SSLP, SSLPR, and DCAP problems. Each subproblem is solved with the CPLEX (12.7.0) quadratic solver. We use the cost-proportional (CP) heuristic to set the values of $\rho \left(i\right)$. The multipliers $\alpha$ in the CP heuristic are set to 0.1, 0.3, 0.5, and 1.0, respectively. Note that the main results shown in this section are not using WW heuristics, i.e., we do not use the variable fixing and slamming, or cycle-detection heuristics. We will make a comparison of PySP with WW heuristics and PySP without WW heuristics at the end of this section.

The number of iterations and the walltime for different multipliers are shown in Figure 1 and Figure 2, respectively. If the PH algorithm reaches iteration limit, there is an “(i)” label at the top of the column. If the PH algorithm reaches the time limit, there is a “(t)” label on top of the column. From Figure 1 and Figure 2, one can observe that setting the $\alpha$ value to 0.1 makes PH take the largest number of iterations and largest amount of walltime to converge in most of the instances, which may be due to the small step size. On the other hand, setting $\alpha$ to the largest value, i.e., 1.0, takes fewer iterations and less walltime than using other $\alpha$ values in most instances. However, it runs out of the iteration limit in two of the instances. Overall, setting $\alpha$ to 0.3 seems to be a robust choice because cp-0.3 always converges within a reasonable walltime and number of iterations. The details of the SSLP and SSLPR results are shown in

Table A1. Computational results for PySP on SSLP.

Instances	Scenarios	cp Multiplier	Iterations	Time	UB
SSLP_5_25	50	cp 0.1	24	17.59	−121.60
		cp 0.3	11	8.16	−121.60
		cp 0.5	13	8.53	−121.60
		cp 1.0	12	7.80	−121.60
	100	cp 0.1	19	27.38	−127.37
		cp 0.3	11	15.27	−127.37
		cp 0.5	18	21.26	−127.37
		cp 1.0	>300	303.53	−125.59
SSLP_10_50	50	cp 0.1	68	254.15	−364.64
		cp 0.3	22	71.57	−364.64
		cp 0.5	15	49.19	−364.64
		cp 1.0	13	41.91	−364.64
	100	cp 0.1	95	540.21	−354.19
		cp 0.3	31	149.33	−354.19
		cp 0.5	19	88.54	−354.19
		cp 1.0	12	57.98	−354.19
	500	cp 0.1	174	4322.45	−349.14
		cp 0.3	61	1265.18	−349.14
		cp 0.5	35	688.90	−349.14
		cp 1.0	19	379.64	−349.14
	1000	cp 0.1	180	9984.56	−351.71
		cp 0.3	59	2849.90	−357.71
		cp 0.5	35	1604.27	−357.71
		cp 1.0	18	845.16	−351.71
SSLP_15_45	5	cp 0.1	9	12.48	−262.40
		cp 0.3	4	5.85	−261.20
		cp 0.5	>300	197.35	−261.20
		cp 1.0	10	9.65	−261.20
	10	cp 0.1	155	485.38	−260.50
		cp 0.3	55	108.15	−260.50
		cp 0.5	33	58.23	−259.30
		cp 1.0	>300	463.63	−259.30
	15	cp 0.1	186	1416.21	−253.60
		cp 0.3	74	286.78	−253.60
		cp 0.5	38	126.79	−253.60
		cp 1.0	22	65.08	−253.20

and

Table A2. Computational results for PySP on SSLPR.

Instances	cp Multiplier	Iterations	Time	UB
SSLPRF_5_25_100_1	cp 0.1	33	722.61	−74,005.84
	cp 0.3	20	461.1	−74,005.84
	cp 0.5	19	427.04	−74,005.84
	cp 1.0	16	315.67	−74,005.84
SSLPRF_5_25_100_2	cp 0.1	28	852.75	−72,671.95
	cp 0.3	18	517.33	−72,671.95
	cp 0.5	15	436.16	−72,671.95
	cp 1.0	21	539.61	−72,671.95
SSLPRF_5_25_100_3	cp 0.1	34	735.9	−75,664.19
	cp 0.3	21	412.91	−75,664.19
	cp 0.5	20	397.01	−75,664.19
	cp 1.0	18	348.75	−75,664.19
SSLPRF_5_50_100_1	cp 0.1	59	2457.66	138,900.12
	cp 0.3	20	814.04	138,900.12
	cp 0.5	12	429.35	138,900.12
	cp 1.0	6	200.07	138,900.12
SSLPRF_5_50_100_2	cp 0.1	94	>10,800	−
	cp 0.3	95	7382.14	245,424.96
	cp 0.5	114	4315.02	245,424.96
	cp 1.0	75	3008.51	500,144.07
SSLPRF_5_50_100_3	cp 0.1	88	>10,800	−
	cp 0.3	83	6984.09	258,578.79
	cp 0.5	50	3887.54	258,578.79
	cp 1.0	24	1727.73	258,578.79

in Appendix A.

We also apply PySP to solve DCAP instances. We observe that for all the DCAP instances, PySP is unable to converge within 300 iterations. The details of the results are shown in

Table A3. Computational results for PySP on DCAP.

Instances	Scenarios	cp Multiplier	Iterations	Time	UB
DCAP 233	200	cp 0.1	>300	456.18	2206.7
		cp 0.3	>300	457.5	2431.8
		cp 0.5	>300	458.23	1966.1
		cp 1.0	>300	455.88	1952.6
	300	cp 0.1	>300	734.23	1862.0
		cp 0.3	>300	726.59	1943.3
		cp 0.5	>300	726.77	1831.6
		cp 1.0	>300	741.17	1815.4
	500	cp 0.1	>300	1515.27	2498.1
		cp 0.3	>300	1492.59	2000.0
		cp 0.5	>300	1467.64	1939.2
		cp 1.0	>300	1494.13	1893.8
DCAP 243	200	cp 0.1	>300	481.96	2465.9
		cp 0.3	>300	478.34	2454.2
		cp 0.5	>300	481.51	2369.3
		cp 1.0	>300	466.95	2383.2
	300	cp 0.1	>300	792.01	2825.8
		cp 0.3	>300	756.17	2802.8
		cp 0.5	>300	710.44	2755.0
		cp 1.0	>300	776.97	2743.1
	500	cp 0.1	>300	1690.74	2196.0
		cp 0.3	>300	1622.8	2235.2
		cp 0.5	>300	1674.76	2216.7
		cp 1.0	>300	1536.42	2323.3
DCAP 332	200	cp 0.1	>300	456.58	1362.7
		cp 0.3	>300	427.65	1529.1
		cp 0.5	>300	450.88	1278.0
		cp 1.0	>300	460.35	1171.0
	300	cp 0.1	>300	714.83	1332.1
		cp 0.3	>300	698.52	1948.9
		cp 0.5	>300	709.21	1904.4
		cp 1.0	>300	706.75	1766.3
	500	cp 0.1	>300	1464.11	1768.6
		cp 0.3	>300	1451.73	1822.8
		cp 0.5	>300	1473.66	1846.6
		cp 1.0	>300	1452.64	1861.6
DCAP 342	200	cp 0.1	>300	449.9	1993.4
		cp 0.3	>300	476.77	1990.4
		cp 0.5	>300	445.11	1870.6
		cp 1.0	>300	470.82	1830.8
	300	cp 0.1	>300	722.6	2260.3
		cp 0.3	>300	739.94	2371.1
		cp 0.5	>300	690.76	2497.7
		cp 1.0	>300	702.23	2255.9
	500	cp 0.1	>300	1582.9	2198.1
		cp 0.3	>300	1604.98	2317.5
		cp 0.5	>300	1555.9	2290.6
		cp 1.0	>300	1593.41	2097.5

in Appendix A where the walltime of the upper bound for those instances is reported.We compare the upper bound obtained by PySP with those obtained by DSP in the next subsection. From this experiment, we can see that it is more difficult for PySP to solve problems with mixed-binary first-stage variables than problems with pure binary first-stage variables because it is more difficult for the continuous variables to satisfy the NACs.

Scenario bundling [66,67,68] is a technique that has been used in dual decomposition algorithms. The main idea is to dualize only “some” of the non-anticipativity constraints, rather than dualizing all of them. In other words, the individual scenario subproblems are bundled into larger subproblems in which the NACs are preserved. Ryan et al. [67] used PH with scenario bundling to solve stochastic unit commitment problems. The authors showed that with the use of scenario bundling, PH can obtain solutions with a smaller optimality gap. In order to test the effectiveness of the scenario bundling, we test several instances from the SSLP and DCAP libraries. The computational results are shown in

Table 2. Computational results for PySP with scenario bundling.

Instances	Scenarios	Bundles	Iterations	Time	UB
SSLP_5_25	50	10	4	3.73	−121.60
		50	24	17.59	−121.60
	100	10	2	3.99	−127.37
		50	7	10.03	−127.37
SSLP_10_50	50	10	4	30.55	−364.64
		50	68	254.15	−364.64
	100	10	2	83.67	−354.19
		50	49	263.59	−354.19
		100	95	540.21	−354.19
	500	10	1	476.13	−349.13
		50	2	162.54	−349.14
		500	174	4322.45	−349.14
	1000	10	1	7137.61	−351.71
		50	1	313.07	−351.71
		1000	180	9984.56	−351.71
DCAP233	200	10	>300	342.97	1854.36
		50	>300	232.21	1861.63
		200	>300	456.18	2206.68
	300	10	147	>10,800	−
		50	>300	317.28	1679.80
		300	>300	1515.27	2498.12
	500	10	>300	634.60	1749.87
		50	>300	400.59	1858.98
		500	>300	1494.13	1893.83

. For each instance, we try a different number of bundles. For the SSLP instances, PH with a different number of bundles can obtain the same upper bound. However, the number of bundles has a significant impact on the computational time. For example, for SSLP_10_50 with 1000 scenarios, PH with 50 bundles can reduce the walltime of the original PH with 1000 bundles to 3%. Additionally, it only takes PH with 50 bundles one iteration to converge. For DCAP problems, PH does not converge within 300 iterations for most cases, even with scenario bundling. However, PH is able to obtain better feasible solutions with scenario bundling (see UB in Table 2).

Finally, we evaluate how the use of WW heuristics can affect the performance of PySP on the SSLP and SSLPR libraries. The results on DCAP library are omitted here since PySP does not converge for DCAP instances. The solution time improvements by using WW heuristics for each SSLP and SSLPR instance are shown in Figure 3. Note that there are three cases where the improvements are omitted in the figure: case (1)—neither PH nor PH with WW heuristics can converge in 300 iterations; case (2)—only PH-WW fails to converge in 300 iterations; and case (3)—both PH and PH-WW exceed the time limit of 3 h (10,800 CPU seconds). Using WW heuristics gives significant improvements for small cost-proportional multipliers, i.e., 0.1 and 0.3. As we have described in Table 1, PH with small multipliers usually takes more iterations to converge. Therefore, the WW heuristics can accelerate convergence for those instances, more effectively. However, there are also a few instances where PH can converge, but PH with WW heuristics cannot converge, which are denoted by case (2) in Figure 3.

6.2. DSP: Decompositions for Structured Programming

DSP [49] is an open-source software package implemented in C++ that provides different types of dual decomposition algorithms to solve stochastic mixed-integer programs (SMIPs). DSP can take SMPS files, and JuMP models as input via a Julia package Dsp.jl.

6.2.1. Algorithmic Innovations in DSP

From the description of the dual decomposition algorithm in Section 5, one can observe that the lower bounds of the dual decomposition algorithm are affected by the way the Lagrangian multipliers are updated. One advantage of DSP is that the authors have implemented different dual-search methods including the subgradient method, the cutting plane method, and a novel interior-point cutting-plane method for the Lagrangian master problem. The authors observed that if the simplex algorithm is used, the solutions to the Lagrangian master problem can oscillate significantly, especially when the Lagrangian dual function is not well approximated. Therefore, the authors proposed to solve the Lagrangian master problem suboptimally using an interior point method, which follows from the work of Mitchell [69].

The authors also proposed some tightening inequalities that are valid for the Lagrangian subproblems. These valid inequalities, including feasibility and optimality cuts, are obtained from Benders subproblems, where the integrality constraints are relaxed. Computational results show that the Benders-like cuts can be effective in practice. The latest DSP distribution is able to complete Algorithm 3 using the Dantzig–Wolfe decomposition-based branching strategy described in [51].

6.2.2. Computational Results for DSP in Comparison with PySP

We test the dual decomposition algorithm on the SSLP, SSLPR, and DCAP libraries. Each subproblem is solved with CPLEX (12.7.0) mixed-integer linear solver. The interior point method proposed by the authors [49] is used to solve the Lagrangian master problem, which is solved with CPLEX as well. Benders-like cuts are not used because the implementation of Benders cuts in DSP only works with SCIP.

In Figure 4 and Figure 5, we evaluate the best feasible solution (the upper bound) obtained by PySP, and the upper and lower bounds obtained by DSP. For each instance, we include three different gaps. The upper and lower bounds from Ref. [61] are used to evaluate the bounds from PySP and DSP. Note that the bounds from literature are close to the global optimality of each instance. The first column for each instance in both Figure 4 and Figure 5 is the gap between the upper bound from PySP ($PyS{P}_{ub}$) and the lower bound from the literature ($Literatur{e}_{lb}$). The second column represents the gap between the upper bound from DSP ($DS{P}_{ub}$) and the lower bound from the literature ($Literatur{e}_{lb}$). The third column represents the gap between the upper bound from literature ($Literatur{e}_{ub}$) and the lower bound from DSP ($DS{P}_{lb}$).

For the SSLP and SSLPR instances shown in Figure 4, although PySP can converge within the time and iteration limit, the best feasible solution obtained from PySP ($PyS{P}_{ub}$) may not be optimal. There are about 1% gaps for some of the SSLP instances (see the first column of each instance in Figure 4). DSP can solve more instances to optimality than PySP (see the second column of each instance in Figure 4). The lower bounds obtained by DSP are also quite tight, usually less than 0.01% (see the third column of each instance in Figure 4). Note that the literature values for SSLPRF5_50_100(1), SSLPRF5_50_100(2), and SSLPRF5_50_100(3) do not match the values from our experiment. Therefore, we try to solve the deterministic equivalent of these instances to obtain bounds. The literature bounds of SSLPRF5_50_100(3) come from solving the deterministic equivalent. The gaps of SSLPRF5_50_100(1) and SSLPRF5_50_100(2) are omitted since the corresponding deterministic equivalent cannot be solved within 3 h.

For the DCAP instances where we have mixed-integer first-stage variables, the best feasible solutions from PySP ($PyS{P}_{ub}$) and DSP ($DS{P}_{ub}$) are quite far from optimality. The gaps of the first two columns are around 10%. On the other hand, the lower bounds obtained from DSP ($DS{P}_{lb}$) are tight. The gaps between ($Literatur{e}_{ub}$) and ($DS{P}_{lb}$) are around 0.1%. Therefore, in order to improve the relative optimality gap of DSP, the focus should be on designing more advanced heuristics to obtain better feasible solutions.

6.3. Results from Branch-and-Bound Algorithm in DSP

As we were completing this paper, a new version of DSP [51] was developed that uses the branch-and-bound algorithm similar to Algorithm 3. We test the branch-and-bound decomposition on the SLSP, SSLPR, and DCAP test instances. The computational setting is the same as the one described in the previous section.

Table 3. Results of DSP with DW-based branch-and-bound on SSLP instances.

Instances	Scenarios	Time	Solution
SSLP5_25	50	1.70	−121.6
	100	3.16	−127.37
SSLP10_50	50	103.9	−364.64
	100	101.62	−364.191
	500	1110.64	−349.136
	1000	3021.83	−351.711
SSLP15_45	5	11.22	−262.4
	10	151.79	−260.5
	15	377.65	−253.601

,

Table 4. Results of DSP with DW Branch-and- Bound on SSLPRF instances.

Instances	Time	LB	UB	GAP [%]
SSLPRF_5_25_100_1	697.74	−74,005.8	−74,005.8	0
SSLPRF_5_25_100_2	1021.76	−72,675	−72,675	0
SSLPRF_5_25_100_3	673.04	−75,666.9	−75,666.9	0
SSLPRF_5_50_100_1	3071.82	138,900	138,900	0
SSLPRF_5_50_100_2	9083.39	163,931	163,931	0
SSLPRF_5_50_100_3	-	-	-	-

and

Table 5. Results of DSP with DW Branch-and- Bound on DCAP instances.

Instances	Scenarios	Time	LB	UB	GAP [%]
DCAP233	200	417.08	1834.57	1834.57	0
	300	512.97	1644.25	1644.25	0
	500	898.05	1834.57	1834.57	0
DCAP243	200	370.68	2322.5	2322.5	0
	300	525.91	2559.19	2559.19	0
	500	1523.61	2167.4	2167.4	0
DCAP332	200	949.64	1060.75	1060.75	0
	300	1379.98	1252.77	1252.77	0
	500	-	1588.660	1589.203	0.03
DCAP342	200	2732.97	1619.61	1619.61	0
	300	1970.75	2067.5	2067.5	0
	500	3958.25	1904.74	1904.74	0

show the results obtained by DSP. We highlight that almost all the instances were solved to optimality within the walltime limit.

In all the SSLP instances (except SSLP15_45), DSP using the DW-based branch-and-bound (DSP-DW) is able to solve the instances to optimality in times lower than the ones from the DSP implementing only the first iteration of DD. The same trend is replicated on the SSLPRF instances.

On the other hand DSP-DW is remarkably slower in attaining algorithmic convergence in the DCAP test set. This behavior can be explained by looking at the $DS{P}_{lb}$ vs. $Literatur{e}_{ub}$ gap in Figure 4 and Figure 5. The standard implementation of DSP is able to close the optimality gap in the SSLP and SSLPR instances by solving the Lagrangian relaxation problem (15a) and using heuristics to compute a primal bound. Nonetheless, this strategy does not work in the DCAP set for which additional computation is devoted by DSP-DW to close the gap via branch-and-bound.

7. Conclusions

We presented a summary of the state-of-the-art methodologies for the solution of two-stage linear stochastic problems. First, we introduced the mathematical formulation of such programs and highlighted features in their structure, which enable the development of decomposition algorithms. These methodologies are classified in two groups: node-based decomposition and scenario-based decomposition.

For two-stage stochastic programs with continuous recourse, we summarized Benders decomposition, which partitions the problem according to its time structure. BD may present computational problems, which can be alleviated by reducing the cost of each iteration, and/or decreasing the number of iterations.

Scenario decomposition methodologies are popular alternatives in the case of (mixed) integer recourse. The progressive hedging algorithm and dual decomposition relax the nonanticipatity restrictions and provide the user with valid bounds. Our numerical results show that the performance of PH is strongly affected by the constant penalty multiplier. Furthermore, its performance and the quality of the approximate solution may be enhanced by grouping the scenarios in large bundles (or scenario sets). We also tested the dual decomposition algorithm with the DSP package. The computational results show that DSP is able to provide a tight lower bound on the instances that we tested. However, the optimality gaps can be as large as 10%, relative to the upper bound from literature. Therefore, for those tested instances, future effort should be focused on developing more advanced heuristics to improve the best feasible solution. Moreover, we tested the newest DSP implementation, which incorporates a Dantzig–Wolfe decomposition-based branch-and-bound algorithm to close the optimality gap. In most of the instances, it is able to retrieve a global solution in the given walltime limit. We highlight that its computational performance depends on the strength of the dual decomposition.

At the time we completed this paper, several new developments were made in stochastic programming software. These software packages have the capability of solving multi-stage stochastic programs, e.g., SDDP.jl [70], StochasticPrograms.jl [71], and MSPPy [72]. In the terms of algorithmic development, recent works have extended the capability to solve nonlinear problems [4,73,74,75,76,77,78]. Furthermore, it should be noted that the paper focuses on the risk-neutral setting of two-stage stochastic programs. However, the studied approaches can be easily extended to the risk-averse setting [79], e.g., by including some risk measures, such as the conditional value at risk (CVar).