Combining Progressive Hedging with a Frank-Wolfe Method to Compute Lagrangian Dual Bounds in Stochastic Mixed-Integer Programming

TL;DR

This paper introduces a novel primal-dual algorithm combining progressive hedging with a Frank-Wolfe inspired linearization to efficiently compute Lagrangian dual bounds in stochastic mixed-integer programming, with proven convergence and superior empirical performance.

Contribution

The paper presents a new algorithm that guarantees convergence to the optimal dual value and improves bound computation in SMIP by integrating Frank-Wolfe linearization steps into progressive hedging.

Findings

Algorithm converges to the optimal Lagrangian dual value.

Empirical results show improved bounds over standard progressive hedging.

Demonstrates efficiency in solving stochastic mixed-integer programs.

Abstract

We present a new primal-dual algorithm for computing the value of the Lagrangian dual of a stochastic mixed-integer program (SMIP) formed by relaxing its nonanticipativity constraints. This dual is widely used in decomposition methods for the solution of SMIPs. The algorithm relies on the well-known progressive hedging method, but unlike previous progressive hedging approaches for SMIP, our algorithm can be shown to converge to the optimal Lagrangian dual value. The key improvement in the new algorithm is an inner loop of optimized linearization steps, similar to those taken in the classical Frank-Wolfe method. Numerical results demonstrate that our new algorithm empirically outperforms the standard implementation of progressive hedging for obtaining bounds in SMIP.

Tables

Table 1. Table 1: Result summary for CAP-101-250, with the absolute percentage gap based on the known optimal value 733827.3

	Percentage gap			# Iterations			Termination
Penalty	PH	FW-PH		PH	FW-PH		PH	FW-PH
		$\alpha =0$	$\alpha =1$		$\alpha =0$	$\alpha =1$		$\alpha =0$	$\alpha =1$
20	0.08%	0.10%	0.11%	466	439	430	T	T	T
100	0.01%	0.00%	0.00%	178	406	437	C	T	T
500	0.07%	0.00%	0.00%	468	92	93	T	C	C
1000	0.15%	0.00%	0.00%	516	127	130	T	C	C
2500	0.34%	0.00%	0.00%	469	259	274	T	C	C
5000	0.66%	0.00%	0.00%	33	431	464	C	T	T
7500	0.99%	0.00%	0.00%	28	18	19	C	C	C
15000	1.59%	0.00%	0.00%	567	28	33	T	C	C

Table 2. Table 2: Result summary for DCAP-233-500, with the absolute percentage gap based on the best known lower bound 1737.7.

	Percentage gap			# Iterations			Termination
Penalty	PH	FW-PH		PH	FW-PH		PH	FW-PH
		$\alpha =0$	$\alpha =1$		$\alpha =0$	$\alpha =1$		$\alpha =0$	$\alpha =1$
2	0.13%	0.12%	0.12%	1717	574	600	T	T	T
5	0.22%	0.09%	0.09%	2074	589	574	T	T	T
10	0.23%	0.07%	0.07%	2598	592	587	T	T	T
20	0.35%	0.07%	0.07%	1942	590	599	T	T	T
50	1.25%	0.06%	0.06%	2718	597	533	T	T	T
100	1.29%	0.06%	0.06%	2772	428	438	T	C	C
200	2.58%	0.06%	0.06%	2695	256	262	T	C	C
500	2.58%	0.07%	0.07%	2871	244	246	T	C	C

Table 3. Table 3: Result summary for SSLP-5-25-50, with the absolute percentage gap based on the known optimal value -121.6

	Percentage gap			# Iterations			Termination
Penalty	PH	FW-PH		PH	FW-PH		PH	FW-PH
		$\alpha =0$	$\alpha =1$		$\alpha =0$	$\alpha =1$		$\alpha =0$	$\alpha =1$
1	0.30%	0.00%	0.00%	105	115	116	C	C	C
2	0.73%	0.00%	0.00%	51	56	56	C	C	C
5	0.91%	0.00%	0.00%	25	26	27	C	C	C
15	3.15%	0.00%	0.00%	12	16	17	C	C	C
30	6.45%	0.00%	0.00%	12	18	18	C	C	C
50	9.48%	0.00%	0.00%	18	25	26	C	C	C
100	9.48%	0.00%	0.00%	8	45	45	C	C	C

Table 4. Table 4: Result summary for SSLP-10-50-100, with the absolute percentage gap based on the known optimal value -354.2

	Percentage gap			# Iterations			Termination
Penalty	PH	FW-PH		PH	FW-PH		PH	FW-PH
		$\alpha =0$	$\alpha =1$		$\alpha =0$	$\alpha =1$		$\alpha =0$	$\alpha =1$
1	0.57%	0.22%	0.23%	126	234	233	T	T	T
2	0.63%	0.03%	0.03%	127	226	228	T	T	T
5	1.00%	0.00%	0.00%	104	219	220	C	T	T
15	2.92%	0.00%	0.00%	33	45	118	C	C	C
30	4.63%	0.00%	0.00%	18	21	22	C	C	C
50	4.63%	0.00%	0.00%	11	26	27	C	C	C
100	4.63%	0.00%	0.00%	9	43	45	C	C	C

Equations

$${\zeta }^{SMIP}:=\underset{x}{\min }\left\{c^{\top }x+\mathcal{Q}(x):x\in X\right\},$$

$$\mathcal{Q}(x):={\mathbb{E}}_{\xi }\left[\underset{y}{\min }\left\{q(\xi {)}^{\top }y:W(\xi )y=h(\xi )-T(\xi )x,y\in Y(\xi )\right\}\right].$$

$$\zeta^{SMIP}=\min_{x,y}\left\{\begin{array}[]{c}c^{\top}x+\displaystyle\sum_{s\in\mathcal{S}}p_{s}q_{s}^{\top}y_{s}\;:\;(x,y_{s})\in K_{s},\forall s\in\mathcal{S}\end{array}\right\},$$

$$\zeta^{SMIP}=\min_{x,y,{z}}\left\{\begin{array}[]{c}\!\!\!\displaystyle\sum_{s\in\mathcal{S}}p_{s}(c^{\top}x_{s}+q_{s}^{\top}y_{s})\;:\;(x_{s},y_{s})\in K_{s},x_{s}={z},\;\forall s\in\mathcal{S},z\in\mathbb{R}^{n_{x}}\end{array}\!\!\!\right\}.$$

$$\phi(\mu):=\min_{x,y,{z}}\left\{\begin{array}[]{c}\sum_{s\in\mathcal{S}}\left[p_{s}(c^{\top}x_{s}+q_{s}^{\top}y_{s})+\mu_{s}^{\top}(x_{s}-{z})\right]\;:\\ (x_{s},y_{s})\in K_{s},\ \forall s\in\mathcal{S},z\in\mathbb{R}^{n_{x}}\end{array}\right\},$$

$$\varphi (\omega ):=\underset{x,y,z}{\min }\left\{\sum _{s\in \mathcal{S}}{p}_s{L}_s(x_s,y_s,z,{\omega }_s):(x_s,y_s)\in K_s,\forall s\in \mathcal{S},z\in {\mathbb{R}}^{n_x}\right\},$$

$$L_s(x_s,y_s,z,{\omega }_s):=c^{\top }{x}_s+q_s^{\top }{y}_s+{\omega }_s^{\top }(x_s-z).$$

$$\varphi (\omega )=\sum _{s\in \mathcal{S}}{p}_s{\varphi }_s({\omega }_s),$$

$${\varphi }_s({\omega }_s):=\underset{x,y}{\min }\left\{(c+{\omega }_s{)}^{\top }x+q_s^{\top }y:(x,y)\in K_s\right\}.$$

$${\zeta }^{LD}:=\underset{\omega }{\sup }\left\{\varphi (\omega ):\sum _{s\in \mathcal{S}}{p}_s{\omega }_s=0\right\}.$$

$${\zeta }^{LD}=\underset{x,y,z}{\min }\left\{\sum _{s\in \mathcal{S}}{p}_s(c^{\top }{x}_s+q_s^{\top }{y}_s):(x_s,y_s)\in \text{conv}(K_s),x_s=z,\forall s\in \mathcal{S}\right\},$$

$$L^{\rho }(x,y,z,\omega ):=\sum _{s\in \mathcal{S}}{p}_s{L}_s^{\rho }(x_s,y_s,z,{\omega }_s),$$

$$L_s^{\rho }(x_s,y_s,z,{\omega }_s):=c^{\top }{x}_s+q_s^{\top }{y}_s+{\omega }_s^{\top }(x_s-z)+\frac{\rho }{2}{\Vert x_s-z\Vert }_2^2$$

$$\sum _{s\in \mathcal{S}}{p}_s\left[{\Vert x_s^k-z^k\Vert }_2^2+{\Vert z^k-z^{k-1}\Vert }_2^2\right]=\sum _{s\in \mathcal{S}}{p}_s{\Vert x_s^k-z^{k-1}\Vert }_2^2$$

$$(x_s^k,y_s^k)\in \underset{x,y}{\mathrm{argmin}}\left\{L_s^{\rho }(x,y,z^{k-1},{\omega }_s^k):(x,y)\in \text{conv}(K_s)\right\}$$

$${\zeta }^{FW}:=\underset{x}{\min }\left\{f(x):x\in D\right\},$$

$$\hat{x}\in \underset{x}{\mathrm{argmin}}\left\{{\nabla }_{x}f(x^{t-1}{)}^{\top }x:x\in D\right\}$$

$$x^t\in \underset{x}{\mathrm{argmin}}\left\{f(x):x\in D^t\right\}.$$

$${\Gamma }^t:=-{\nabla }_{x}f(x^{t-1}{)}^{\top }(\hat{x}-x^{t-1})<\tau ,$$

$$\nabla_{(x,y)}L_{s}^{\rho}(x,y,{z},\omega_{s})|_{(x,y)=\left(x^{t-1}_{s},y^{t-1}_{s}\right)}=\left[\begin{array}[]{c}c+\omega_{s}+\rho(x^{t-1}_{s}-z)\\ q_{s}\end{array}\right]=\left[\begin{array}[]{c}c+\widehat{\omega}_{s}\\ q_{s}\end{array}\right],$$

$$(x^{t}_{s},y^{t}_{s},a)\in\operatorname*{argmin}_{x,y,a}\left\{\begin{array}[]{c}L_{s}^{\rho}(x,y,z,\omega_{s})\;:\;(x,y)=\sum_{(\widehat{x}^{i},\widehat{y}^{i})\in V^{t}_{s}}a_{i}(\widehat{x}^{i},\widehat{y}^{i}),\\ \sum_{i=1,\dots,|V^{t}_{s}|}a_{i}=1,\;\text{and}\;a_{i}\geq 0\;\text{for}\;i=1,\dots,|V^{t}_{s}|\end{array}\right\}.$$

$$(x_s^t,y_s^t)\notin \underset{x,y}{\mathrm{argmin}}\left\{L_s^{\rho }(x,y,z^{k-1},{\omega }_s^k):(x,y)\in \text{conv}(K_s)\right\}$$

$$(x_s^t,y_s^t)\in \underset{x,y}{\mathrm{argmin}}\left\{L_{\rho }^s(x,y,z^{k-1},{\omega }_s^k):(x,y)\in \text{conv}(V_s^t)\right\}$$

$$\nabla_{(x,y)}L_{\rho}^{s}(x_{s}^{t},y_{s}^{t},{z}^{k-1},\omega_{s}^{k})\left[\begin{array}[]{c}x-x_{s}^{t}\\ y-y_{s}^{t}\end{array}\right]\geq 0\quad\text{for all}\;(x,y)\in\text{conv}(V^{t}_{s})$$

$$\min_{x,y}\left\{\nabla_{(x,y)}L_{\rho}^{s}(x_{s}^{t},y_{s}^{t},{z}^{k-1},\omega_{s}^{k})\left[\begin{array}[]{c}x-x_{s}^{t}\\ y-y_{s}^{t}\end{array}\right]\;:\;(x,y)\in\text{conv}(K_{s})\right\}<0.$$

$${\tilde{\omega }}_s^k$$

$$\sum _{s\in \mathcal{S}}{p}_s(x_s^{k-1}-z^{k-1})=0\text{and}\sum _{s\in \mathcal{S}}{p}_s{\omega }_s^k=0,$$

$$\bigcap _{s\in \mathcal{S}}{\text{Proj}}_x(\text{conv}(V_s^0))\ne \varnothing$$

$$\text{conv}(V_s^k)=:{\overline{D}}_s\subseteq \text{conv}(K_s).$$

$$\begin{gathered} (x_s^k,y_s^k)\in \underset{x,y}{\mathrm{argmin}}\left\{L_s^{\rho }(x,y,z^{k-1},{\omega }_s^k): \\ (x,y)\in \text{conv}(K_s)\right\}. \end{gathered}$$

Full text

Keywords: mixed-integer stochastic programming, Lagrangian duality, progressive hedging, Frank-Wolfe method

1 Introduction

Stochastic programming with recourse provides a framework for modeling problems where decisions are made in stages. Between stages, some uncertainty in the problem parameters is unveiled, and decisions in subsequent stages may depend on the outcome of this uncertainty. When some decisions are modeled using discrete variables, the problem is known as a Stochastic Mixed-Integer Programming (SMIP) problem. The ability to simultaneously model uncertainty and discrete decisions make SMIP a powerful modeling paradigm for applications. Important applications employing SMIP models include unit commitment and hydro-thermal generation scheduling [25, 34], military operations [32], vaccination planning [28, 35], air traffic flow management [4], forestry management and forest fire response [5, 26], and supply chain and logistics planning [19, 21]. However, the combination of uncertainty and discreteness makes this class of problems extremely challenging from a computational perspective. In this paper, we present a new and effective algorithm for computing lower bounds that arise from a Lagrangian-relaxation approach.

The mathematical statement of a two-stage SMIP is

[TABLE]

where the vector $c\in\mathbb{R}^{n_{x}}$ is known, and $X$ is a mixed-integer linear set consisting of linear constraints and integer restrictions on some components of $x$. The function $\mathcal{Q}:\mathbb{R}^{n_{x}}\mapsto\mathbb{R}$ is the expected recourse value

[TABLE]

We assume that the random variable ${\mathbf{\xi}}$ is taken from a discrete distribution indexed by the finite set $\mathcal{S}$, consisting of the realizations, $\xi_{1},\dots,\xi_{|\mathcal{S}|}$, with strictly positive corresponding probabilities of realization, $p_{1},\dots,p_{|\mathcal{S}|}$. When ${\mathbf{\xi}}$ is not discrete, a finite scenario approximation can be obtained via Monte Carlo sampling [18, 23] or other methods [10, 9]. Each realization $\xi_{s}$ of ${\mathbf{\xi}}$, is called a scenario and encodes the realizations observed for each of the random elements $\left(q(\xi_{s}),h(\xi_{s}),W(\xi_{s}),T(\xi_{s}),Y(\xi_{s})\right)$. For notational brevity, we refer to this collection of random elements respectively as $\left(q_{s},h_{s},W_{s},T_{s},Y_{s}\right)$. For each $s\in\mathcal{S}$, the set $Y_{s}\subset\mathbb{R}^{n_{y}}$ is a mixed-integer set containing both linear constraints and integrality constraints on a subset of the variables, $y_{s}$.

The problem (1) may be reformulated as its deterministic equivalent

[TABLE]

where $K_{s}:=\left\{(x,y_{s})\;:\;W_{s}y_{s}=h_{s}-T_{s}x,x\in X,y_{s}\in Y_{s}\right\}$. Problem (2) has a special structure that can be algorithmically exploited by decomposition methods. To induce a decomposable structure, scenario-dependent copies $x_{s}$ for each $s\in S$ of the first-stage variable $x$ are introduced to create the following reformulation of (2):

[TABLE]

The constraints $x_{s}={z}$, $s\in\mathcal{S}$, enforce nonanticipativity for first-stage decisions; the first-stage decisions $x_{s}$ must be the same $({z})$ for each scenario $s\in\mathcal{S}$. Applying Lagrangian relaxation to the nonanticipativity constraints in problem (3) yields the nonanticipative Lagrangian dual function

[TABLE]

where $\mu=(\mu_{1},\dots,\mu_{|\mathcal{S}|})\in\prod_{s\in\mathcal{S}}\mathbb{R}^{n_{x}}$ is the vector of multipliers associated with the relaxed constraints $x_{s}={z}$, $s\in\mathcal{S}$. By setting $\omega_{s}:=\frac{1}{p_{s}}\mu_{s}$, (4) may be rewritten as

[TABLE]

where

[TABLE]

Since ${z}$ is unconstrained in the optimization problem in the definition (5), in order for the Lagrangian function $\phi(\omega)$ to be bounded from below, we require as a condition of dual feasibility that $\sum_{s\in\mathcal{S}}p_{s}\omega_{s}=0$. Under this assumption, the ${z}$ term vanishes, and the Lagrangian dual function (5) decomposes into separable functions,

[TABLE]

where for each $s\in\mathcal{S}$,

[TABLE]

The reformulation (6) is the basis for parallelizable approaches for computing dual bounds that are used, for example, in the dual decomposition methods developed in [8, 22].

For any choice of $\omega=\left(\omega_{1},\dots,\omega_{|\mathcal{S}|}\right)$, it is well-known that the value of the Lagrangian provides a lower bound on the optimal solution to (1): $\phi(\omega)\leq\zeta^{SMIP}$. The problem of finding the best such lower bound is the Lagrangian dual problem:

[TABLE]

The primary contribution of this work is a new and effective method for solving (8), thus enabling a practical and efficient computation of high-quality lower bounds for $\zeta^{SMIP}$.

The function $\phi(\omega)$ is a piecewise-affine concave function, and many methods are known for maximizing such functions. These methods include the subgradient method [33], the augmented Lagrangian (AL) method [15, 29], and the alternating direction method of multipliers (ADMM) [13, 11, 7]. The subgradient method has mainly theoretical significance, since it is difficult to develop reliable and efficient step-size rules for the dual variables $\omega$ (see, e.g., Section 7.1.1 of [31]). As iterative primal-dual approaches, methods based on the AL method or ADMM are more effective in practice. However, in the context of SMIP, both methods require convexification of the constraints $K_{s}$, $s\in\mathcal{S}$ to have a meaningful theoretical support for convergence to the best lower bound value $\zeta^{LD}$. Furthermore, both methods require the solution of additional mixed-integer linear programming (MILP) subproblems in order to recover the Lagrangian lower bounds associated with the dual values, $\omega$ [14]. ADMM has a more straightforward potential for decomposability and parallelization than the AL method, and so in this work, we develop a theoretically-supported modification of a method based on ADMM.

When specialized to the deterministic equivalent problem (2) in the context of stochastic programming, ADMM is referred to as Progressive Hedging (PH) [30, 37]. When the sets $K_{s}$, $s\in\mathcal{S}$, are convex, the limit points of the sequence of solution-multiplier pairs $\left\{((x^{k},y^{k},{z}^{k}),\omega^{k})\right\}_{k=1}^{\infty}$ generated by PH are saddle points of the deterministic equivalent problem (2), whenever such saddle points exist. When the constraints $(x_{s},y_{s})\in K_{s}$, $s\in\mathcal{S}$, enforce nontrivial mixed-integer restrictions, the set $K_{s}$ is not convex and PH becomes a heuristic approach with no guarantees of convergence [20]. Nevertheless, some measure of success, in practice, has been observed in [37] while applying PH to problems of the form (3). More recently, [14] showed that valid Lagrangian lower bounds can be calculated from the iterates of the PH algorithm when the sets $K_{s}$ are not convex. However, their implementation of the algorithm does not offer any guarantee that the lower bounds will converge to the optimal value $\zeta^{LD}$. Moreover, additional computational effort, in solving additional MILP subproblems, must be expended, in order to compute the lower bound. Our contribution is to extend the PH-based approach in [14], creating an algorithm whose lower bound values converge to $\zeta^{LD}$, in theory, and for which lower bound calculations do not require additional computational effort. Computational results in Section 4 demonstrate that the new method outperforms the existing PH-based method, in terms of both quality of bound and efficiency of computation.

To motivate our approach, we first consider the application of PH to the following well-known primal characterization of $\zeta^{LD}$:

[TABLE]

where $\text{conv}(K_{s})$ denotes the convex hull of $K_{s}$ for each $s\in\mathcal{S}$. (See, for example, Theorem 6.2 of [24].) The sequence of Lagrangian bounds $\left\{\phi(\omega^{k})\right\}$ generated by the application of PH to (9) is known to be convergent. Thus, the value of the Lagrangian dual, $\zeta^{LD}$, may, in theory, be computed by applying PH to (9). However, in practice, an explicit polyhedral description of $\text{conv}(K_{s})$, $s\in\mathcal{S}$ is, generally, not available, thus raising the issue of implementability.

The absence of such an explicit description motivates an application of a solution approach to the PH primal update step that iteratively constructs an improved inner approximation of each $\text{conv}(K_{s})$, $s\in\mathcal{S}$. For this purpose, we apply a solution approach to the PH primal update problem that is based on the Frank-Wolfe (FW) method [12]. Our approach has the additional benefit of providing Lagrangian bounds at no additional computational cost.

One simple, theoretically-supported integration of a FW-like method and PH is realized by having the PH primal updates computed using a method called the Simplicial Decomposition Method (SDM) [16, 36]. SDM is an extension of the FW method that makes use of progressively-improving inner approximations to each set $\text{conv}(K_{s})$, $s\in\mathcal{S}$. The finite optimal convergence of each application of SDM follows directly from the polyhedral structure $\text{conv}(K_{s})$, and the (practically reasonable) assumption that $\text{conv}(K_{s})$ is bounded for each $s\in\mathcal{S}$. inline]Where do we assume/state that $\text{conv}(K_{s})$ is bounded? Proof of convergence in Proposition 1 relies on boundedness. (I thought we just needed it to be polyhedral, but I am happy to assume bounded for this paper)

For computing improvements in the Lagrangian bound efficiently, convergence of SDM to the optimal solution of the subproblem is too costly and not necessary. We thus develop a modified integration whose theoretically-supported convergence analysis is based not on the optimal convergence of SDM, but rather on its ability to adequately extend the inner approximations of each $\text{conv}(K_{s})$, $s\in\mathcal{S}$.

The main contribution of this paper is the development, convergence analysis, and application of a new algorithm, called FW-PH, which is used to compute high-quality Lagrangian bounds for SMIPs efficiently and with a high potential for parallelization. FW-PH is efficient in that, under mild assumptions, each dual update and Lagrangian bound computation may be obtained by solving, for each $s\in\mathcal{S}$, just one MILP problem and one continuous convex quadratic problem. In contrast, each dual update of PH requires the solution of a mixed-integer quadratic programming (MIQP) subproblem for each $s\in\mathcal{S}$, and each PH Lagrangian bound computation requires the solution of one MILP subproblem for each $s\in\mathcal{S}$. In our convergence analysis, conditions are provided under which the sequence of Lagrangian bounds generated by FW-PH converges to the optimal Lagrangian bound $\zeta^{LD}$. To the best of our knowledge, the combination of PH and FW in a manner that is theoretically supported, computationally efficient, and parallelizable is new, in spite of the convergence analyses of both PH and FW being well-developed.

This paper is organized as follows. In Section 2, we present the theoretical background of PH and a brief technical lemma regarding the inner approximations generated by SDM; this background is foundational for the proposed FW-PH method. In Section 3, we present the FW-PH method and a convergence analysis. The results of numerical experiments comparing the Lagrangian bounds computed with PH and those with FW-PH are presented in Section 4. We conclude in Section 5 with a discussion of the results obtained and with suggested directions for further research.

2 Progressive Hedging and Frank-Wolfe-Based Methods

The Augmented Lagrangian (AL) function based on the relaxation of the nonanticipativity constraints $x_{s}={z}$, $s\in\mathcal{S}$, is

[TABLE]

where

[TABLE]

and $\rho>0$ is a penalty parameter. By changing the feasible region, denoted here by $D_{s}$, $s\in\mathcal{S}$, the Augmented Lagrangian can be used in a Progressive Hedging approach to solve either problem (3) or problem (9). Pseudocode for the PH algorithm is given in Algorithm 1.

In Algorithm 1, $k_{max}>0$ is the maximum number of iterations, and $\epsilon>0$ parameterized the convergence tolerance. The initialization of Lines 3–8 provides an initial target primal value $z^{0}$ and dual values $\omega_{s}^{1}$, $s\in\mathcal{S}$, for the main iterations $k\geq 1$. Also, an initial Lagrangian bound $\phi^{0}$ can be computed from this initialization. For $\epsilon>0$, the Algorithm 1 termination criterion $\sqrt{\sum_{s\in\mathcal{S}}p_{s}\left\|x_{s}^{k}-{z}^{k-1}\right\|_{2}^{2}}<\epsilon$ is motivated by the addition of the squared norms of the primal and dual residuals associated with problem (3). (See Section 3.3 of [7].) In summing the squared norm primal residuals $p_{s}\left\|x_{s}^{k}-{z}^{k}\right\|_{2}^{2}$, $s\in\mathcal{S}$, and the squared norm dual residual $\left\|{z}^{k}-{z}^{k-1}\right\|_{2}^{2}$, we have

[TABLE]

The equality in (10) follows since, for each $s\in\mathcal{S}$, the cross term resulting from the expansion of the squared norm $\left\|(x_{s}^{k}-{z}^{k})+({z}^{k}-{z}^{k-1})\right\|_{2}^{2}$ vanishes; this is seen in the equality $\sum_{s\in\mathcal{S}}p_{s}(x_{s}^{k}-{z}^{k})=0$ due to the construction of $z^{k}$.

The Line 11 subproblem of Algorithm 1 is an addition to the original PH algorithm. Its purpose is to compute Lagrangian bounds (Line 14) from the current dual solution $\omega^{k}$ [14]. Thus, the bulk of computational effort in Algorithm 1 applied to problem (3) (the case with $D_{s}=K_{s}$) resides in computing solutions to the MILP (Line 11) and MIQP (Line 12) subproblems. Note that Line 11 may be omitted if the corresponding Lagrangian bound for $\omega^{k}$ is not desired.

2.1 Convergence of PH

The following proposition addresses the convergence of PH applied to problem (9).

Proposition 1**.**

Assume that problem (9) is feasible with $\text{conv}(K_{s})$ bounded for each $s\in\mathcal{S}$, and let Algorithm 1 be applied to problem (9) (so that $D_{s}=\text{conv}(K_{s})$ for each $s\in\mathcal{S}$) with tolerance $\epsilon=0$ for each $k\geq 1$. Then, the limit $\lim_{k\to\infty}\omega^{k}=\omega^{*}$ exists, and furthermore,

$\lim_{k\to\infty}\sum_{s\in\mathcal{S}}p_{s}(c^{\top}x_{s}^{k}+q_{s}^{\top}y_{s}^{k})=\zeta^{LD}$, 2. 2.

$\lim_{k\to\infty}\phi(\omega^{k})=\zeta^{LD}$, 3. 3.

$\lim_{k\to\infty}(x_{s}^{k}-z^{k})=0$* for each $s\in\mathcal{S}$,*

and each limit point $(((x_{s}^{*},y_{s}^{*})_{s\in\mathcal{S}},{z}^{*})$ is an optimal solution for (9).

Proof.

Since the constraint sets $D_{s}=\text{conv}(K_{s})$, $s\in\mathcal{S}$, are bounded, and problem (9) is feasible, problem (9) has an optimal solution $((x_{s}^{*},y_{s}^{*})_{s\in\mathcal{S}},z^{*})$ with optimal value $\zeta^{LD}$. The feasibility of problem (9), the linearity of its objective function, and the bounded polyhedral structure of its constraint set $D_{s}=\text{conv}(K_{s})$, $s\in\mathcal{S}$, imply that the hypotheses for PH convergence to the optimal solution are met (See Theorem 5.1 of [30]). Therefore, $\left\{\omega^{k}\right\}$ converges to some $\omega^{*}$, $\lim_{k\to\infty}\sum_{s\in\mathcal{S}}p_{s}(c^{\top}x_{s}^{k}+q_{s}^{\top}y_{s})=\zeta^{LD}$, $\lim_{k\to\infty}\phi(\omega^{k})=\zeta^{LD}$, and $\lim_{k\to\infty}(x_{s}^{k}-z^{k})=0$ for each $s\in\mathcal{S}$ all hold. The boundedness of each $D_{s}=\text{conv}(K_{s})$, $s\in\mathcal{S}$, furthermore implies the existence of limit points $((x_{s}^{*},y_{s}^{*})_{s\in\mathcal{S}},z^{*})$ of $\left\{((x_{s}^{k},y_{s}^{k})_{s\in\mathcal{S}},z^{k})\right\}$, which are optimal solutions for (9). ∎

Note that the convergence in Proposition 1 applies to the continuous problem (9) but not to the mixed-integer problem (3). In problem (3), the constraint sets $K_{s}$, $s\in\mathcal{S}$, are not convex, so there is no guarantee that Algorithm 1 will converge when applied to (3). However, the application of PH to problem (9) requires, in Line 12, the optimization of the Augmented Lagrangian over the sets $\text{conv}(K_{s})$, $s\in\mathcal{S}$, for which an explicit linear description is unlikely to be known. In the next section, we demonstrate how to circumvent this difficulty by constructing inner approximations of the polyhedral sets $\text{conv}(K_{s})$, $s\in\mathcal{S}$.

2.2 A Frank-Wolfe approach based on Simplicial Decomposition

To use Algorithm 1 to solve (9) requires a method for solving the subproblem

[TABLE]

appearing in Line 12 of the algorithm. Although an explicit description of $\text{conv}(K_{s})$ is not readily available, if we have a linear objective function, then we can replace $\text{conv}(K_{s})$ with $K_{s}$. This motivates the application of a FW algorithm for solving (11), since the FW algorithm solves a sequence of problems in which the nonlinear objective is linearized using a first-order approximation.

The Simplicial Decomposition Method (SDM) is an extension of the FW method, where the line searches of FW are replaced by searches over polyhedral inner approximations. SDM can be applied to solve a feasible, bounded problem of the general form

[TABLE]

with nonempty compact convex set $D$ and continuously differentiable convex function $f$. Generically, given a current solution $x^{t-1}$ and inner approximation $D^{t-1}\subseteq D$, iteration $t$ of the SDM consists of solving

[TABLE]

updating the inner approximation as $D^{t}\leftarrow\text{conv}(D^{t-1}\cup\left\{\widehat{x}\right\})$, and finally choosing

[TABLE]

The algorithm terminates when the bound gap is small, specifically, when

[TABLE]

where $\tau\geq 0$ is a given tolerance.

The application of SDM to solve problem (11), i.e., to minimize $L_{s}^{\rho}(x,y,{z},\omega_{s})$ over $(x,y)\in\text{conv}(K_{s})$, for a given $s\in\mathcal{S}$, is presented in Algorithm 2. Here, $t_{max}$ is the maximum number of iterations and $\tau>0$ is a convergence tolerance. $\Gamma^{t}$ is the bound gap used to measure closeness to optimality, and $\phi_{s}$ is used to compute a Lagrangian bound as described in the next section. The inner approximation to $\text{conv}(K_{s})$ at iteration $t\geq 1$ takes the form $\text{conv}(V^{t}_{s})$, where $V_{s}^{t}$ is a finite set of points, with $V^{t}_{s}\subset\text{conv}(K_{s})$. The points added by Algorithm 2 to the initial set, $V^{0}_{s}$, to form $V^{t}_{s}$, are all in $K_{s}$: here $\mathcal{V}(\text{conv}(K_{s}))$ is the set of extreme points of $\text{conv}(K_{s})$ and, of course, $\mathcal{V}(\text{conv}(K_{s}))\subseteq K_{s}$.

Observe that

[TABLE]

and so the optimization at Line 5 is minimizing the gradient approximation to $L_{s}^{\rho}(x,y,{z},\omega_{s})$ at the point $(x^{t-1}_{s},y^{t-1}_{s})$. Since this is a linear objective function, optimization over $\mathcal{V}(\text{conv}(K_{s}))$ can be accomplished by optimization over $K_{s}$ (see, e.g., [24], Section I.4, Theorem 6.3). Hence Line 5 requires a solution of a single-scenario MILP.

The optimization at Line 11 can be accomplished by expressing $(x,y)$ as a convex combination of the finite set of points, $V^{t}_{s}$, where the weights $a\in\mathbb{R}^{|V^{t}_{s}|}$ in the convex combination are now also decision variables. That is, the Line 11 problem is solved with a solution to the following convex continuous quadratic subproblem

[TABLE]

For implementational purposes, the $x$ and $y$ variables may be substituted out of the objective of problem (13), leaving $a$ as the only decision variable, with the only constraints being nonnegativity of the $a$ components and the requirement that they sum to $1$.

The Simplicial Decomposition Method is known to terminate finitely with an optimal solution when $D$ is polyhedral [16], so the primal update step Line 12, Algorithm 1 with $D_{s}=\text{conv}(K_{s})$ could be accomplished with the SDM, resulting in an algorithm that converges to a solution giving the Lagrangian dual bound $\zeta^{LD}$. However, since each inner iteration of Line 5, Algorithm 2 requires the solution of a MILP, using $t_{\max}$ large enough to ensure SDM terminates optimally is not efficient for our purpose of computing Lagrangian bounds. In the next section, we give an adapation of the algorithm that requires the solution of only one MILP subproblem per scenario at each major iteration of the PH algorithm.

3 The FW-PH method

In order to make the SDM efficient when used with PH to solve the problem (9), the minimization of the Augmented Lagrangian can be done approximately. This insight can greatly reduce the number of MILP subproblems solved at each inner iteration and forms the basis of our algorithm FW-PH. Convergence of FW-PH relies on the following lemma, which states an important expansion property of the inner approximations employed by the SDM.

Lemma 1**.**

For any scenario $s\in\mathcal{S}$ and iteration $k\geq 1$, let Algorithm 2 be applied to the minimization problem (11) for any $t_{max}\geq 2$. For $1\leq t<t_{max}$, if

[TABLE]

*holds, then $\text{conv}(V_{s}^{t+1})\supset\text{conv}(V_{s}^{t})$. *

Proof.

For $s\in\mathcal{S}$ and $k\geq 1$ fixed, we know that by construction

[TABLE]

for $t\geq 1.$ Given the convexity of $(x,y)\mapsto L_{\rho}^{s}(x,y,{z}^{k-1},\omega_{s}^{k})$ and the convexity of $\text{conv}(V^{t}_{s})$, the necessary and sufficient condition for optimality

[TABLE]

is satisfied. By assumption, condition (14) is satisfied, $\text{conv}(K_{s})$ is likewise convex, and so the resulting non-satisfaction of the necessary and sufficient condition of optimality for the problem in (14) takes the form

[TABLE]

In fact, during SDM iteration $t+1$, an optimal solution $(\widehat{x}_{s},\widehat{y}_{s})$ to the problem in condition (16) is computed in Line 5 of Algorithm 2. Therefore, by the satisfaction of condition (15) and the optimality of $(\widehat{x}_{s},\widehat{y}_{s})$ for the problem of condition (16), which is also satisfied, we have $(\widehat{x}_{s},\widehat{y}_{s})\not\in\text{conv}(V^{t}_{s})$. By construction, $V^{t+1}_{s}\leftarrow V^{t}_{s}\cup\left\{(\widehat{x}_{s},\widehat{y}_{s})\right\}$, so that $\text{conv}(V^{t+1}_{s})\supset\text{conv}(V^{t}_{s})$ must hold. ∎

The FW-PH algorithm is stated in pseudocode-form in Algorithm 3. Similar to Algorithm 1, the parameter $\epsilon$ is a convergence tolerance, and $k_{max}$ is the maximum number of (outer) iterations. The parameter $t_{max}$ is the maximum number of (inner) SDM iterations in Algorithm 2.

The parameter $\alpha\in\mathbb{R}$ affects the initial linearization point $\tilde{x}_{s}$ of the SDM method. Any value $\alpha\in\mathbb{R}$ may be used, but the use of $\tilde{x}_{s}=(1-\alpha)z^{k-1}+\alpha x_{s}^{k-1}$ in Line 6 is a crucial component in the efficiency of the FW-PH algorithm, as it enables the computation of a valid dual bound, $\phi^{k}$, at each iteration of FW-PH without the need for additional MILP subproblem solutions. Specifically, we have the following result.

Proposition 2**.**

Assume that the precondition $\sum_{s\in\mathcal{S}}p_{s}\omega^{0}_{s}=0$ holds for Algorithm 3. At each iteration $k\geq 1$ of Algorithm 3, the value, $\phi^{k}$, calculated at Line 9, is the value of the Lagrangian relaxation $\phi(\cdot)$ evaluated at a Lagrangian dual feasible point, and hence provides a finite lower bound on $\zeta^{LD}$.

Proof.

Since $\sum_{s\in\mathcal{S}}p_{s}\omega^{0}_{s}=0$ holds and, by construction, $0=\sum_{s\in\mathcal{S}}p_{s}(x_{s}^{0}-z^{0})$, we have $\sum_{s\in\mathcal{S}}p_{s}\omega^{1}_{s}=0$ also. We proceed by induction on $k\geq 1$. At iteration $k$, the problem solved for each $s\in\mathcal{S}$ at Line 5 in the first iteration ($t=1$) of Algorithm 2 may be solved with the same optimal value by exchanging $\mathcal{V}(\text{conv}(K_{s}))$ for $K_{s}$; this follows from the linearity of the objective function. Thus, an optimal solution computed at Line 5 may be used in the computation of $\phi_{s}(\widetilde{\omega}_{s}^{k})$ carried out in Line 7, where

[TABLE]

By construction, we have at each iteration $k\geq 1$ in Algorithm 3 that

[TABLE]

which establishes that $\sum_{s\in\mathcal{S}}p_{s}\widetilde{\omega}_{s}^{k}=0$. Thus, $\widetilde{\omega}^{k}$ is feasible for the Lagrangian dual problem, so that $\phi(\widetilde{\omega}^{k})=\sum_{s\in\mathcal{S}}p_{s}\phi_{s}^{k}$, and, since each $\phi_{s}^{k}$ is the optimal value of a bounded and feasible mixed-integer linear program, we have $-\infty<\phi(\widetilde{\omega}^{k})<\infty$. ∎

We establish convergence of Algorithm 3 for any $\alpha\in\mathbb{R}$ and $t_{max}\geq 1$. For the special case where we perform only one iteration of SDM for each outer iteration ($t_{max}=1$), we require the additional assumption that the initial scenario vertex sets share a common point. More precisely, we require the assumption

[TABLE]

which can, in practice, be effectively handled through appropriate initialization, under the standard assumption of relatively complete recourse. We describe one initialization procedure in Section 4.

Proposition 3**.**

*Let the convexified separable deterministic equivalent SMIP (9) have an optimal solution, and let Algorithm 3 be applied to (9) with $k_{max}=\infty$, $\epsilon=0$, $\alpha\in\mathbb{R}$, and $t_{max}\geq 1$. If either $t_{max}\geq 2$ or (17) holds, then $\lim_{k\to\infty}\phi^{k}=\zeta^{LD}$. *

Proof.

First note that for any $t_{max}\geq 1$, the sequence of inner approximations $\text{conv}(V_{s}^{k})$, $s\in\mathcal{S}$, will stabilize, in that, for some threshold $0\leq\bar{k}_{s}$, we have for all $k\geq\bar{k}_{s}$

[TABLE]

This follows due to the assumption that each expansion of the inner approximations $\text{conv}(V_{s}^{k})$ take the form $V_{s}^{k}\leftarrow V_{s}^{k-1}\cup\left\{(\widehat{x}_{s},\widehat{y}_{s})\right\}$, where $(\widehat{x}_{s},\widehat{y}_{s})$ is a vertex of $\text{conv}(K_{s})$. Since each polyhedron $\text{conv}(K_{s})$, $s\in\mathcal{S}$ has only a finite number of such vertices, the stabilization (18) must occur at some $\bar{k}_{s}<\infty$.

For $t_{\max}\geq 2$, the stabilizations (18), $s\in\mathcal{S}$, are reached at some iteration $\bar{k}:=\max_{s\in\mathcal{S}}\left\{\bar{k}_{s}\right\}$. Noting that $\overline{D}_{s}=\text{conv}(V_{s}^{k})$ for $k>\overline{k}$ we must have

[TABLE]

Otherwise, due to Lemma 1, the call to SDM on Line 7 must return $V_{s}^{k}\supset V_{s}^{k-1}$, contradicting the finite stabilization (18). Therefore, the $k\geq\bar{k}$ iterations of Algorithm 3 are identical to Algorithm 1 iterations, and so Proposition 1 implies that $\lim_{k\to\infty}x_{s}^{k}-{z}^{k}=0$, $s\in\mathcal{S}$, and $\lim_{k\to\infty}\phi(\omega^{k})=\zeta^{LD}$. By the continuity of $\omega\mapsto\phi_{s}(\omega)$ for each $s\in\mathcal{S}$, we have $\lim_{k\to\infty}\phi^{k}=\lim_{k\to\infty}\sum_{s\in\mathcal{S}}p_{s}\phi_{s}(\omega_{s}^{k}+\alpha(x_{s}^{k-1}-z^{k-1}))=\lim_{k\to\infty}\sum_{s\in\mathcal{S}}p_{s}\phi_{s}(\omega_{s}^{k})=\lim_{k\to\infty}\phi(\omega^{k})=\zeta^{LD}$ for all $\alpha\in\mathbb{R}$.

In the case $t_{max}=1$, we have at each iteration $k\geq 1$ the optimality

[TABLE]

By the stabilization (18), the iterations $k\geq\bar{k}$ of Algorithm 3 are identical to PH iterations applied to the restricted problem

[TABLE]

We have initialized the sets $(V_{s}^{0})_{s\in\mathcal{S}}$ such that $\cap_{s\in\mathcal{S}}\text{Proj}_{x}\text{conv}(V_{s}^{0})\neq\emptyset$, so since the inner approximations to $\text{conv}(K_{s})$ only expand in the algorithm, $\cap_{s\in\mathcal{S}}\text{Proj}_{x}(\overline{D}_{s})\neq\emptyset$. Therefore, problem (20) is a feasible and bounded linear program, and so the PH convergence described in Proposition 1 with $D_{s}=\overline{D}_{s}$, $s\in\mathcal{S}$, holds for its application to problem (20). That is, for each $s\in\mathcal{S}$, we have 1) $\lim_{k\to\infty}\omega_{s}^{k}=\omega_{s}^{*}$ and $\lim_{k\to\infty}(x_{s}^{k}-{z}^{k})=0$; and 2) for all limit points $((x_{s}^{*},y_{s}^{*})_{s\in\mathcal{S}},{z}^{*})$, we have the feasibility and optimality of the limit points, which implies $x_{s}^{*}={z}^{*}$ and

[TABLE]

Next, for each $s\in\mathcal{S}$ the compactness of $\text{conv}(K_{s})\supseteq\overline{D}_{s}$, the continuity of the minimum value function

[TABLE]

and the limit $\lim_{k\to\infty}\widetilde{\omega}_{s}^{k+1}=\lim_{k\to\infty}{\omega}_{s}^{k+1}+\alpha\rho(x_{s}^{k}-{z}^{k})=\omega_{s}^{*}$, together imply that

[TABLE]

Recall that $\widetilde{\omega}_{s}^{k}=\omega_{s}^{k}+\rho\alpha(x_{s}^{k-1}-z^{k-1})$ is the $t=1$ value of $\widehat{\omega}_{s}^{t}$ defined in Line 4 of Algorithm 2. Thus, for $k+1>\bar{k}$, we have due to the stabilization (18) that

[TABLE]

If equality (23) does not hold, then the inner approximation expansion $\overline{D}_{s}\subset\text{conv}(V_{s}^{k+1})$ must occur, since a point $(\widehat{x}_{s},\widehat{y}_{s})\in\text{conv}(K_{s})$ that can be strictly separated from $\overline{D}_{s}$ would have been discovered during the iteration $k+1$ execution of Algorithm 2, Line 5, $t=1$. The expansion $\overline{D}_{s}\subset\text{conv}(V_{s}^{k+1})$ contradicts the finite stabilization (18), and so (23) holds. Therefore, the equalities (22) and (23) imply that

[TABLE]

Our argument has shown that for all limit points $(x_{s}^{*},y_{s}^{*})$, $s\in\mathcal{S}$, the following stationarity condition is satsfied:

[TABLE]

which together with the feasibility $x_{s}^{*}={z}^{*}$, $s\in\mathcal{S}$ implies that each limit point $((x_{s}^{*},y_{s}^{*})_{s\in\mathcal{S}},{z}^{*})$ is optimal for problem (9) and $\omega^{*}$ is optimal for the dual problem (8).

Thus, for all $t_{max}\geq 1$, we have shown $\lim_{k\to\infty}(x_{s}^{k}-{z}^{k})=0$, $s\in\mathcal{S}$, and $\lim_{k\to\infty}\phi(\omega^{k})=\zeta^{LD}$. By similar reasoning used in the $t_{max}\geq 2$ case, it is straightforward that for all $\alpha\in\mathbb{R}$, we also have $\lim_{k\to\infty}\phi^{k}=\zeta^{LD}$. ∎

While using a large value of $t_{max}$ more closely matches Algorithm 3 to the original PH algorithm as described in Algorithm 1, we are motivated to use a small value of $t_{max}$ since the work per iteration is proportional to $t_{max}$. Specifically, each iteration requires solving $t_{max}|\mathcal{S}|$ MILP subproblems, and $t_{max}|\mathcal{S}|$ continuous convex quadratic subproblems. (For reference, Algorithm 1 applied to problem (3) requires the solution of $|\mathcal{S}|$ MIQP subproblems for each $\omega$ update and $|\mathcal{S}|$ MILP subproblems for each Lagrangian bound $\phi$ computation.)

4 Numerical Experiments

We performed computations using a C++ implementation of Algorithm 1 ($D_{s}=K_{s}$, $s\in\mathcal{S}$) and Algorithm 3 using CPLEX 12.5 [17] as the solver for all subproblems. For reading SMPS files into scenario-specific subproblems and for their interface with CPLEX, we used modified versions of the COIN-OR [3] Smi and Osi libraries. The computing environment is the Raijin cluster maintained by Australia’s National Computing Infrastructure (NCI) and supported by the Australian Government [1]. The Raijin cluster is a high performance computing (HPC) environment which has 3592 nodes (system units), 57472 cores of Intel Xeon E5-2670 processors with up to 8 GB PC1600 memory per core (128 GB per node). All experiments were performed in a serial setting using a single node and one thread per CPLEX solve.

In the experiments with Algorithms 1 and 3, we set the convergence tolerance at $\epsilon=10^{-3}$. For Algorithm 3, we set $t_{max}=1$. Also, for all experiments performed, we set $\omega^{0}=0$. In this case, convergence of our algorithm requires that (17) holds, which can be guaranteed during the initialization of the inner approximations $(V^{0}_{s})_{s\in\mathcal{S}}$. Under the standard assumption of relatively complete resource, i.e., the assumption that for all $x\in X$ and $s\in\mathcal{S}$ there exists $y_{s}$ such that $(x,y_{s})\in K_{s}$, a straightforward mechanism for ensuring this assumption is to solve the recourse problems for any fixed $\widehat{x}\in X$. Specifically, for each $s\in\mathcal{S}$, let

[TABLE]

and initialize $V_{s}^{0}$ for each $s\in\mathcal{S}$ so that $\{(\widehat{x},\widehat{y}_{s})\}\in V_{s}^{0}$. Observe also that this initialization corresponds to a technique for computing a feasible solution to the original problem (2), which is independently useful for obtaining an upper bound on $\zeta^{SMIP}$.

For the computational experiments, we run the following initialization to obtain $(V_{s}^{0})_{s\in\mathcal{S}}$ and $(x_{s}^{0},y_{s}^{0})_{s\in\mathcal{S}}$ that are input into Algorithm 3:

If problem (2) does not have relatively complete recourse, then any means to compute a feasible solution to (2) may be employed to initialize each $V_{s}^{0}$, $s\in\mathcal{S}$, in a way to satisfy (17).

Two sets of Algorithm 3 experiments correspond to variants considering $\alpha=0$ and $\alpha=1$. Computations were performed on four problems: the CAP instance 101 with the first 250 scenarios (CAP-101-250) [6], the DCAP instance DCAP233_500 with 500 scenarios, the SSLP instances SSLP5.25.50 with 50 scenarios (SSLP-5-25-50) and SSLP10.50.100 with 100 scenarios (SSLP-10-50-100). The latter three problems are described in detail in [27, 2] and accessible at [2]. For each problem, computations were performed for different penalty values $\rho>0$. The penalty values used in the experiments for the SSLP-5-25-50 instance were chosen to include those penalties that are tested in a computational experiment with PH whose results are depicted in Figure 2 of [14]. For the other problem instances, the set of penalty values $\rho$ tested is chosen to capture a reasonably wide range of performance potential for both PH and FW-PH. All computational experiments were allowed to run for a maximum of two hours in wall clock time.

Tables 1–4 provide a summary indicating the quality of the Lagrangian bounds $\phi$ computed at the end of each experiment for the four problems with varying penalty parameter $\rho$. In each of these tables, the first column lists the values of the penalty parameter $\rho$, while the following are presented for PH and FW-PH (for both $\alpha=0$ and $\alpha=1$) computations in the remaining columns: 1) the absolute percentage gap $\left|\frac{\zeta^{*}-\phi}{\zeta^{*}}\right|*100\%$ between the computed Lagrangian bound $\phi$ and some reference value $\zeta^{*}$ that is either a known optimal value for the problem, or a known best upper bound thereof (column “Percentage Gap”); 2) the total number of dual updates (“# Iterations”); and 3) the indication of whether the algorithm terminated due to the time limit, indicated by letter “T”, or the satisfaction of the convergence criterion $\sqrt{\sum_{s\in\mathcal{S}}p_{s}\left\|x_{s}^{k}-{z}^{k-1}\right\|_{2}^{2}}<\epsilon$, indicated by letter “C” (column “Termination”).

The following observations can be made from the results presented in Tables 1–4. First, for small values of the penalty $\rho$, there is no clear preference between the bounds $\phi$ generated by PH and FW-PH. However, for higher penalties, the bounds $\phi$ obtained by FW-PH are consistently of better quality (i.e., higher) than those obtained by PH, regardless of the variant used (i.e. $\alpha=0$ or $\alpha=1$). This tendency is typically illustrated, for example, in Table 2, where the absolute percentage gap of the Lagrangian lower bound with the known optimal value was 0.06% with $\rho=200$ for FW-PH ($\alpha=0$), while it was 2.58% for the same value of $\rho$ for PH. This improvement is consistently observed for the other problems and the other values of $\rho$ that are not too close to zero. Also, FW-PH did not terminate with suboptimal convergence or display cycling behavior for any of the penalty values $\rho$ in any of the problems considered. For example, all experiments considered in Table 3 terminated due to convergence. The percentage gaps suggest that the convergence for PH was suboptimal, while it was optimal for FW-PH. Moreover, it is possible to see from these tables that the quality of the bounds $\phi$ obtained using FW-PH were not as sensitive to the value of the penalty parameter $\rho$ as obtained using PH.

The FW-PH with $\alpha=0$ versus PH convergence profiles for the experiments performed are given in Figures 1–4, in which we provide plots of wall time versus Lagrangian bound values based on profiling of varying penalty. The times scales for the plots have been set such that trends are meaningfully depicted (1000s for CAP-101-250 and DCAP-233-500, 100 seconds for SSLP-5-25-50, and 3000s for SSLP-10-50-100). The trend of the Lagrangian bounds is depicted with solid lines for FW-PH with $\alpha=0$ and with dashed lines for PH. Plots providing the same comparison for FW-PH with $\alpha=1$ are provided in Appendix A.

As seen in the plots of Figures 1–4, the Lagrangian bounds $\phi$ generated with PH tend to converge suboptimally, often displaying cycling, for large penalty values. In terms of the quality of the bounds obtained, while there is no clear winner when low penalty $\rho$ values are used, for large penalties, the quality of the bounds $\phi$ generated with FW-PH is consistently better than for the bounds generated with PH, regardless of the $\alpha$ value. This last observation is significant because the effective use of large penalty values $\rho$ in methods based on augmented Lagrangian relaxation tends to yield the most rapid early iteration improvement in the Lagrangian bound; this point is most clearly illustrated in the plot of Figure 3.

5 Conclusions and future work

In this paper, we have presented an alternative approach to compute nonanticipativity Lagrangian bounds associated with SMIPs that combines ideas from the Progressive Hedging (PH) and the Frank-Wolfe (FW) based methods. We first note that while Lagrangian bounds can be recovered with PH, this requires—for each iteration and each scenario—the solution of an additional MILP subproblem in addition to the MIQP subproblem. Furthermore, when using the PH method directly, the Lagrangian bounds may converge suboptimally, cycle (for large penalties), or converge very slowly (for small penalties).

To overcome the lack of theoretical support for the above use of PH, we first described a straightforward integration of PH and a FW-like approach such as the Simplicial Decomposition Method (SDM), where SDM is used to compute the primal updates in PH. Its convergence only requires noting that SDM applied to a convex problem with a bounded polyhedral constraint set terminates finitely with optimal convergence. However, for the stated goal of computing high-quality Lagrangian bounds efficiently, the benefits of relying on the optimal convergence of SDM is far outweighed by the computational costs incurred.

As an alternative, we propose the contributed algorithm, FW-PH, that is analyzed under general assumptions on how the Lagrangian bounds are computed and on the number of SDM iterations used for each dual update. Furthermore, under mild assumptions on the initialization of the algorithm, FW-PH only requires the solution of a MILP subproblem and a continuous convex quadratic subproblem for each iteration and each scenario. FW-PH is versatile enough to handle a wide range of SMIPs with integrality restrictions in any stage, while providing rapid improvement in the Lagrangian bound in the early iterations that is consistent across a wide range of penalty parameter values. Although we have opted to focus on two-stage problems with recourse, the generalization of the proposed approach to the multi-stage case is also possible.

Numerical results are encouraging as they suggest that the proposed FW-PH method applied to SMIP problems usually outperforms the traditional PH method with respect to how quickly the quality of the generated Lagrangian bounds improves. This is especially true with the use of larger penalty values. For all problems considered and for all but the smallest penalties considered, the FW-PH method displayed better performance over PH in terms of the quality of the final Lagrangian bounds at the end of the allotted wall clock time.

The improved performance of FW-PH over PH for large penalties is significant because it is the effective use of large penalties enabled by FW-PH that yields the most rapid initial dual improvement. This last feature of FW-PH would be most helpful in its use within a branch-and-bound or branch-and-cut framework for providing strong lower bounds (in the case of minimization). In addition to being another means to compute Lagrangian bounds, PH would still have a role in such frameworks as a heuristic for computing a primal feasible solution to the SMIP, thus providing (in the case of minimization) an upper bound on the optimal value.

Future research on this subject includes the following. First, FW-PH inherits the potential for parallelization from PH. Experiments for exploring the benefit of parallelization are therefore warranted. Second, the theoretical support of FW-PH can be strengthened with a better understanding of the behavior of PH (and its generalization ADMM) applied to infeasible problems. Finally, FW-PH can benefit from a better understanding of how the proximal term penalty coefficient can be varied to improve performance.

Appendix A Additional plots for PH vs. FW-PH for $\alpha=1$