Combining Penalty-based and Gauss-Seidel Methods for solving Stochastic Mixed-Integer Problems

TL;DR

This paper introduces a new decomposition method combining penalty-based Lagrangian and Gauss-Seidel techniques to efficiently solve stochastic mixed-integer programming problems, demonstrating superior computational performance and solution quality.

Contribution

The paper presents a novel PBGS approach that leverages problem structure for improved efficiency over existing methods like Progressive Hedging.

Findings

PBGS outperforms PH in computational speed

PBGS achieves higher solution quality

Numerical experiments validate efficiency and effectiveness

Abstract

In this paper, we propose a novel decomposition approach for mixed-integer stochastic programming (SMIP) problems that is inspired by the combination of penalty-based Lagrangian and block Gauss-Seidel methods (PBGS). In this sense, PBGS is developed such that the inherent decomposable structure that SMIPs present can be exploited in a computationally efficient manner. The performance of the proposed method is compared with the Progressive Hedging method (PH), which also can be viewed as a Lagrangian-based method for obtaining solutions for SMIP. Numerical experiments performed using instances from the literature illustrate the efficiency of the proposed method in terms of computational performance and solution quality.

Tables

Table 1. Table 1 : CAP101

	$\rho$	500				2500				7500
	$\beta$	1.25		1.11		1.25		1.11		1.25		1.11
	$\gamma$	0.5	1	0.5	1	0.5	1	0.5	1	0.5	1	0.5	1
Obj. diff.	Average	0.12%	0.12%	0.13%	0.09%	0.06%	0.12%	0.12%	0.20%	0.16%	0.08%	0.17%	0.04%
	St. dev.	0.13%	0.15%	0.18%	0.14%	0.10%	0.14%	0.14%	0.20%	0.14%	0.10%	0.14%	0.07%
Speed-up	Average	2.02	1.12	2.14	1.24	1.63	0.93	1.57	1.01	2.12	1.44	2.28	1.60
	St. dev.	0.47	0.29	0.62	0.31	2.31	1.04	1.91	1.12	0.82	0.55	0.86	0.53
PH conv. Fraction.		96.0%				92.0%				94.0%

Table 2. Table 2 : CAP111

	$\rho$	500				2500				7500
	$\beta$	1.25		1.11		1.25		1.11		1.25		1.11
	$\gamma$	0.5	1	0.5	1	0.5	1	0.5	1	0.5	1	0.5	1
Obj. diff.	Average	0.04%	0.10%	0.01%	0.01%	0.02%	0.12%	0.00%	0.02%	0.07%	0.07%	0.06%	0.04%
	St. dev.	0.06%	0.11%	0.03%	0.01%	0.16%	0.31%	0.17%	0.16%	0.05%	0.05%	0.06%	0.05%
Speed-up	Average	3.95	3.67	2.06	1.87	2.98	2.93	1.60	1.73	1.97	2.06	1.19	1.25
	St. dev.	0.27	0.01	0.11	0.37	1.28	1.31	0.51	0.61	0.76	0.81	0.48	0.51
PH conv. Fraction.		4.0%				86.0%				92.0%

Table 3. Table 3 : DCAP233

	$\rho$	5				10				50
	$\beta$	1.25		1.11		1.25		1.11		1.25		1.11
	$\gamma$	0.5	1	0.5	1	0.5	1	0.5	1	0.5	1	0.5	1
Obj. diff.	Average	13.41%	9.56%	7.78%	6.16%	14.22%	13.28%	11.04%	8.10%	18.75%	19.73%	16.86%	14.80%
	St. dev.	2.18%	2.44%	2.26%	2.53%	2.35%	2.88%	2.78%	2.65%	2.06%	2.33%	2.21%	2.57%
Speed-up	Average	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A
	St. dev.	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A
PH conv. Fraction.		0.0%				0.0%				0.0%

Table 4. Table 4 : DCAP342

	$\rho$	5				10				50
	$\beta$	1.25		1.11		1.25		1.11		1.25		1.11
	$\gamma$	0.5	1	0.5	1	0.5	1	0.5	1	0.5	1	0.5	1
Obj. diff.	Average	6.77%	8.36%	7.08%	8.51%	6.01%	9.61%	6.31%	7.59%	6.69%	5.24%	7.27%	6.92%
	St. dev.	4.78%	4.45%	4.47%	4.39%	3.26%	5.36%	5.07%	3.70%	3.82%	2.59%	2.90%	3.67%
Speed-up	Average	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A
	St. dev.	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A
PH conv. Fraction.		0.0%				0.0%				0.0%

Table 5. Table 5 : SSLP5-50

	$\rho$	5				50				100
	$\beta$	1.25		1.11		1.25		1.11		1.25		1.11
	$\gamma$	0.5	1	0.5	1	0.5	1	0.5	1	0.5	1	0.5	1
Obj. diff.	Average	0.23%	0.46%	0.11%	0.10%	0.45%	0.45%	0.21%	0.16%	0.40%	0.51%	0.40%	0.40%
	St. dev.	0.76%	1.14%	0.49%	0.32%	1.08%	1.08%	0.64%	0.53%	1.07%	1.25%	1.07%	1.07%
Speed-up	Average	1.29	1.29	0.76	0.81	1.32	1.40	0.93	1.03	1.12	1.21	0.82	0.86
	St. dev.	0.60	0.56	0.30	0.31	0.91	0.96	0.61	0.67	0.68	0.82	0.58	0.49
PH conv. Fraction.		100.0%				98.0%				98.0%

Table 6. Table 6 : SSLP10-50

	$\rho$	5				50				100
	$\beta$	1.25		1.11		1.25		1.11		1.25		1.11
	$\gamma$	0.5	1	0.5	1	0.5	1	0.5	1	0.5	1	0.5	1
		0.5	1	0.5	1	0.5	1	0.5	1	0.5	1	0.5	1
Obj. diff.	Average	0.00%	0.00%	0.47%	0.00%	4.80%	5.41%	4.19%	2.84%	3.50%	4.54%	2.85%	3.50%
	St. dev.	0.0%	0.0%	0.5%	0.0%	12.5%	12.9%	12.0%	10.5%	8.8%	9.9%	8.0%	8.8%
Speed-up	Average	7.95	7.86	2.95	4.07	2.27	2.35	1.35	1.41	1.69	1.72	1.03	1.10
	St. dev.	2.14	2.08	0.36	1.13	1.11	1.29	0.78	0.74	1.04	1.05	0.69	0.75
PH conv. Fraction.		4.0%				96.0%				90.0%

Equations

$$\displaystyle\zeta^{SIP}:=~{}\min_{x,y}\$$

$$x\in X$$

$$y_s\in Y_s(x),\forall s\in S,$$

$$\displaystyle\zeta^{SIP}:=~{}\min_{x,y,z}\$$

$$x_s-z=0,\forall s\in S$$

$$x_s\in X,\forall s\in S$$

$$y_s\in Y_s(x_s),\forall s\in S.$$

$$\displaystyle\zeta^{LR}(\omega):=~{}\min_{x,y,z}\$$

$$x_s\in X,\forall s\in S$$

$$y_s\in Y_s(x_s),\forall s\in S,$$

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

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

$${\zeta }^{LD}:=\underset{\omega \in \Omega }{\max }{\zeta }^{LR}(\omega ).$$

$$\displaystyle\zeta^{LR+}_{\rho}(\omega):=~{}\min_{x,y,z}\$$

$$x_s\in X,\forall s\in S$$

$$y_s\in Y_s(x_s),\forall s\in S,$$

$$\displaystyle\zeta^{LD+}_{\rho}:=\max_{\omega\in\Omega}\$$

$$\displaystyle\zeta^{SIP}:~{}\min_{x,y,z}\$$

$$x_s-z\le 0,\forall s\in S$$

$$-(x_s-z)\le 0,\forall s\in S$$

$$x_s\in X,\forall s\in S$$

$$y_s\in Y_s(x_s),\forall s\in S.$$

$${\psi }_{\rho }^s(u_s):={\underline{\rho }}_s^{\top }[u_s{]}^{-}+{\overline{\rho }}_s^{\top }[-u_s{]}^{-},$$

$${\psi }_{\rho }(u):=\sum _{s\in S}{\psi }_{\rho }^s(u_s)=\left(\sum _{s\in S}{\underline{\rho }}_s^{\top }[u_s{]}^{-}+\sum _{s\in S}{\overline{\rho }}_s^{\top }[-u_s{]}^{-}\right).$$

$${\psi }_{\infty }(u):=\Vert u{\Vert }_{\infty }$$

$${\psi }_1(u):=\Vert u{\Vert }_1$$

$$\text{ or equivalently }\Vert u{\Vert }_1$$

$${\psi }_{\infty }^N(u)$$

$${\psi }_1^N(u)$$

$${\psi }_1^{N_{\rho }}(u)$$

$${\psi }_A(u)\ge \gamma {\psi }_B(u)\text{ for all }u\in {\mathbb{R}}^m$$

$${\psi }_A(u)=||u||{\psi }_A\left(\frac{u}{||u||}\right)=||u||\frac{{\psi }_A\left(\frac{u}{||u||}\right)}{{\psi }_B\left(\frac{u}{||u||}\right)}{\psi }_B\left(\frac{u}{||u||}\right)\ge \gamma ||u||{\psi }_B\left(\frac{u}{||u||}\right)=\gamma {\psi }_B\left(u\right).$$

$$\delta$$

$$\frac{\epsilon }{\Vert u{\Vert }_{\infty }}{\psi }_{\infty }^N(u)\ge \delta >0$$

$${\psi }_{\infty }^N(u)\ge \delta \frac{\Vert u{\Vert }_{\infty }}{\epsilon }\ge \delta >0,$$

$${\psi }_1^N(u)\ge \eta {\psi }_{\infty }^N(u)\ge \eta \delta >0\text{ for all }u\notin V.$$

$$\frac{\epsilon }{\Vert u{\Vert }_{\infty }}{\psi }_{\infty }^N(u)\ge \delta >0$$

$${\psi }_{\infty }^N(u)\ge \delta \frac{\Vert u{\Vert }_{\infty }}{\epsilon }\ge \frac{\delta }{\epsilon }\Vert u{\Vert }_{\infty }>0,$$

$${\psi }_1^N(u)\ge \eta {\psi }_{\infty }^N(u)\ge \eta \gamma \Vert u{\Vert }_{\infty }>0.$$

Full text

††journal: International Transactions of Operational Research

1 Introduction

Inspired by recent advances and the increased availability of parallel computation resources, there has been a recent surge of methods which are capable of exploiting the structure of large-scale mathematical programming problems to achieve increased efficiency.

One relevant class of problems that can benefit from this paradigm is stochastic mixed-integer programming (SMIP) problems. The modelling framework for this class of problems is versatile as it simultaneously allows for representation of integer-valued decisions and uncertainty in the input data. However, they are frequently challenging in terms of computational tractability due to their inherent NP-hard nature and their large-scale that arise from their scenario-based representations. Parallel computation is particularly appropriate for solving SMIPs, since the special structure of these problems (i.e. their partial separation by decision stage and by outcome scenario) makes them easier to decompose into smaller subproblems which may then be solved simultaneously.

The opportunities arising from approaching SMIP problems by means of decomposition has prompted the development of several different theoretical and algorithmic approaches. For example, the Integer L-Shaped method [22] employs Benders’ decomposition to achieve stage-wise decomposition of SMIP problems. Other algorithms employ Lagrangian duality to achieve scenario-wise decomposition, such as the Dual Decomposition algorithm [10] which uses Lagrangian dual bounds in a branch-and-bound framework, or Progressive Hedging (PH) [29, 23, 33, 32] which applies an Alternating-Direction-type method to the augmented Lagrangian dual problem. Recent studies and applications of these methods include [2, 19, 24, 15] and references therein.

All of the above methods are based on the concept of duality, and therefore they must consider the duality gap that may exist between the optimal solution values of the original (primal) problem and the dual problem. This duality gap is frequently nonzero in the context of non-convex problems, such as those with integer decision variables. If the duality gap for a particular problem is large, any algorithm based on that dual is unlikely to be effective [11].

Several possible approaches to modifying Lagrangian duality to deal with the duality gap which arises in non-convex problems have been considered in the literature. These approaches include $l_{1}$-like penalty functions [11], indicator augmenting functions [21], nonlinear Lagrangian functions [35], and semi-Lagrangian duality [3]. With the exception of nonlinear Lagrangian functions, these approaches have not yet been widely exploited in terms of experimental investigation and practical applications despite numerous theoretical developments available in literature.

In this paper, we propose an alternative approach to deal with SMIP problems that builds upon recent theoretical results from [8] and [13] showing that duality gaps can be diminished with the use of finite-valued penalties for specific class of penalty functions. We show that one can obtain reasonable penalty functions using positive bases and that parallelisation can be obtained by the application of a block Gauss-Seidel approach. The combination of these two frameworks allows us to develop an efficient heuristics that is capable of providing solutions for large-scale SMIP problems. In terms of objective value quality and computational time, the developed approach is shown to be competitive with existing approaches such as PH, which, despite its heuristic nature in the context of SMIPs, has been relied upon as an efficient solution method (see, for example [30, 27, 32]). Furthermore, the theoretical basis for PH does not apply for problems containing integer variables. On the other hand there exists some supporting theory for the PBGS approach which we present in this paper. This partial theory has the potential to inform directions for the further improvement of heuristic methods of this kind.

This paper is structured as follows. In Section 2 we cover the technical background which will be used in the development of the proposed method. Section 3 describes the development of the penalty-based block Gauss-Seidel method, while in Section 4 we discuss computational aspects of the algorithm. Section 5 provides results of the numerical experiments performed. Finally, in Section Acknowledgements we provide conclusions and directions for further development of this research.

2 Technical Background

In the following developments, we consider two-stage stochastic mixed-integer programming problems of the form:

[TABLE]

where $x\in\mathbb{R}^{n_{x}},y\in\mathbb{R}^{n_{y}\times|S|}$ are decision variables, $c\in\mathbb{R}^{n_{x}}$ and $q\in\mathbb{R}^{n_{y}\times|S|}$ are input parameters, and $p_{s}\in\mathbb{R}^{|S|}$ represents the scenario probabilities. Sets $X\subset\mathbb{R}^{n_{x}}$ and $Y_{s}(x)\subset\mathbb{R}^{n_{y}\times|S|}$ define the feasible decision set and consist of linear constraints and integrality restrictions on $x$ and $y$.

To obtain a formulation that is amenable to decomposition, i.e., that can be exploited in terms of its block-angular structure, $\zeta^{SIP}$ may be equivalently rewritten as:

[TABLE]

The set of constraints represented by (5) are referred to in the relevant literature as the non-anticipativity constraints (NAC), as they enforce consensus over the decision made before the observation of a given scenario $s\in S$. It is straightforward to see that these constraints prevent the problem from being solved by means of a decomposition approach that can exploit the otherwise block-angular structure of the problem. Hence, a natural way to approach this class of problems consists of relying on frameworks that are capable of considering relaxations for $\zeta^{SIP}$ which does not include constraint (5).

2.1 Lagrangian Relaxation

A natural approach to solve this problem is to relax the NAC by means of Lagrangian relaxation. To achieve this, let $\lambda=(\lambda_{s})_{s\in S}$ be the Lagrangian multipliers associated with the constraints (5). Then, our Lagrangian relaxation can be formulated as the following problem:

[TABLE]

where $\omega:=(\omega_{s})_{s\in S}=\left(\frac{\lambda_{s}}{p_{s}}\right)_{s\in S}$ and

[TABLE]

In order to guarantee that $\zeta^{LR}(\omega)$ has a bounded optimal solution, one must enforce that the dual feasibility condition $\omega\in\Omega:=\{\omega\mid\sum_{s\in S}p_{s}^{\top}\omega_{s}=0\}$ holds. Under this assumption, (8) may be rewritten as

[TABLE]

which forces the $z$ variable to vanish (note that the term $-\omega_{s}^{\top}z$ was in fact the potential cause of the unboundedness of the dual relaxation $\zeta^{LR}(\omega)$) and its removal yields complete separability for each $s\in S$.

It is well known that $\zeta^{LR}(\omega)\leq\zeta^{SIP}$ for any $\omega\in\Omega$. The Lagrangian dual problem consists of finding the $\omega$ which causes $\zeta^{LR}(\omega)$ to most closely approximate or bound $\zeta^{SIP}$ from below, which in practice means solving the problem

[TABLE]

In general only the weak duality condition $\zeta^{LR}(\omega)\leq\zeta^{SIP}$ holds. When equality $\zeta^{LD}=\zeta^{SIP}$ holds, we have strong duality. Due to the presence of integer restricted variables, the primal problem (4)-(7) is not convex, and strong duality (that is $\zeta^{LD}=\zeta^{SIP}$) cannot be guaranteed. Instead, we typically have a duality gap i.e. $\zeta^{LD}<\zeta^{SIP}$.

2.2 Augmented Lagrangian Approach

In the particular domain of mixed-integer problems such as SMIP problems, there has been renewed interest in the use of augmented Lagrangian approaches [9, 14, 16]. The augmented Lagrangian relaxation of $\zeta^{SIP}$ which relaxes the NAC (5) is:

[TABLE]

where $\omega=(\omega_{s})_{s\in S}\in\Omega$ and $\psi_{\rho}^{s}:\mathbb{R}^{n_{x}}\mapsto\mathbb{R}$ is an appropriate penalty function specific to scenario $s$ that depends on the penalty parameter $\rho$. As in the ordinary Lagrangian relaxation, $\omega\in\Omega$ implies that $\sum_{s\in S}p_{s}^{\top}\omega_{s}=0$ so as to ensure that the dual problem has a finite optimal value. The augmented Lagrangian dual problem is:

[TABLE]

A common choice for the penalty function, in this context, is $\psi_{\rho}^{s}(u_{s}):=\frac{\rho}{2}||u_{s}||_{2}^{2}$ for each $s\in S$, which provides smoothness to the original scenario-wise augmented Lagrangian dual function [28, 5, 6].

Recent results have shown that the augmented Lagrangian dual is capable of asymptotically achieving zero duality gap when the weight $\rho$ associated with the penalty function is allowed to go to infinity [8, Prop. 3], [13, Prop. 2 ]. However, despite the theoretical relevance of this observation, it is not practically meaningful to deal with large-valued penalty parameters, in large part due to the associated numerical issues that arise.

Furthermore, [8, Cor. 1], [13, Thm. 4 ] demonstrates that it is possible to circumvent this drawback if the augmentation of the Lagrangian dual is made using a norm as the penalty function. In this case, the theory suggests that it is possible to attain strong duality for a finite value of $\rho$. This result is one of the major motivations for the developments to be presented next.

2.3 Semi-Lagrangian Duality

Semi-Lagrangian duality [3] is a variant of Lagrangian duality in which ”difficult” equality constraints (e.g. $Ax=b$) are reformulated as pairs of inequality constraints ($Ax\leq b$ and $Ax\geq b$). Lagrangian relaxation is then applied to one of the two sets of inequality constraints. Surrogate semi-Lagrangian duality [25, 20] is a variant which replaces one set of inequalities with its weighted sum ($\lambda^{\top}Ax\leq\lambda^{\top}b$, where $\lambda\geq 0$ is a non-negative vector of the weights applied to each inequality) and then applies Lagrangian relaxation to the resulting single inequality. Special classes of problems can exhibit zero duality gap when utilising the semi-Lagrangian dual problem.

The semi-Lagrangian approach is effective when the semi-Lagrangian dual problem is more tractable than the original problem, even though some inequalities remain in the dual problem as explicit constraints. For problems to which semi-Lagrangian relaxation has been previously applied, such as the p-median problem [3] and the uncapacitated facility location problem [4], the semi-Lagrangian dual problem may be simplified by choosing appropriate dual variable values.

The method presented in this paper is similar to semi-Lagrangian duality in that it can be interpreted as a penalty-method analogue of applying Lagrangian duality to a reformulation of the original problem in which equalities are rewritten as inequalities. To attain this objective, let us first reformulate the problem (1)-(3) into the following equivalent form:

[TABLE]

Unfortunately, in the context of the SMIP problems studied in this paper, the inequalities left as explicit constraints by semi-Lagrangian duality would not allow us to achieve our goal of scenario-wise decomposition.

The method presented here instead relaxes inequalities (14) and (15), so that penalty-based approaches can treat the deviations from the separate inequalities differently. Furthermore, the resulting dual problem is separable by scenario within the block Gauss-Seidel framework. As will be observed later, choosing penalty functions based on some positive bases can result in penalty terms analogous to the objective terms obtained through surrogate semi-Lagrangian duality.

2.4 Desirable Properties of Penalty Functions

Our primary objective is to compute $\zeta^{LD+}_{\rho}$ in a decomposed manner, which will require: $i)$ the definition of a suitable penalty functions $\psi_{\rho}^{s}$ for each $s\in S$; and $ii)$ the application of a block Gauss-Seidel (GS)-based approach based on a decomposable structure.

One important result, originally proven for a general mixed-integer programming (MIP) problems, which can be used in this case is Theorem 5 of [13], which is reproduced below (adapted to the context of SMIP problems).

Theorem 1

[13, Thm. 5 ]** Consider a feasible MIP problem given in (1)-(3) whose problem data is formed from rational entries and with its optimal value bounded. If $\psi:\prod_{s\in S}\mathbb{R}^{n_{x}}\mapsto\mathbb{R}$ is a summed augmenting function $\psi(u):=\sum_{s\in S}\psi_{\rho}^{s}(u_{s})$ for problem (11)–(13) such that

$\psi(0)=0$** 2. 2.

$\psi(u)\geq\delta>0,\forall u\not\in V$** 3. 3.

$\psi(u)\geq\gamma||u||_{\infty},\forall u\in V$**

for some open neighbourhood $V$ of [math], and positive scalars $\delta,\gamma>0$, then there exists a finite $\rho$ such that $\zeta^{LD+}_{\rho}=\zeta^{LR+}_{\rho}(\omega_{LP})=\zeta^{SIP}$, for $\omega_{LP}$ (an optimal multiplier of the linear programming relaxation of the NACs (5)).

Proof. Apply the general theorem [13, Thm. 5 ] to our problem (1)-(3).

Remark 2

In [8] other conditions that do not require the assumption of rationality of the data defining the problem are given that also ensure a limiting zero duality gap.

Remark 3

One may see with little difficulty that the proof of [13] does not rely on the setting of $\omega=\omega_{LP}$. Indeed one can show that for any $\omega$ there still exists a finite (possibly larger) penalty parameter such that Theorem 1 holds true.

In [8] and later in [13] it has been observed that a zero duality gap is achievable for dual problems based on an augmented Lagrangian in MIP problems. In both papers, very general classes of augmenting functions were studied and consequently very little can be inferred as to what would be a practical penalty that one could use on a given problem. It is observed in [13] that the usual quadratic (squared norm) penalty is probably not a practical choice for MIP. One would hope that an augmenting function would lead to a reformulation of the MIP that is not significantly worse to solve than the original problem, which would mean that augmenting functions should lead to a MIP reformulation.

Motivated by the aforementioned facts, we propose a class of augmenting functions based on the use of positive basis [12]. One special case of this class of penalty functions is given as follows. Given discrepancy vector $u:=(u_{s})_{s\in S}\in\prod_{s\in S}\mathbb{R}^{n_{x}}$, we define for each scenario $s$ the penalty function

[TABLE]

where $\rho=(\underline{\rho}_{s},\overline{\rho}_{s})_{s\in S}\in\mathbb{R}^{2n_{x}|S|}_{>0}$ and $[v]^{-}:=-\min\{0,v\}$ (performed component wise), where in this case $v\in\mathbb{R}^{n_{x}}$. Then we define

[TABLE]

In the following developments we will demonstrate that (16) satisfies the conditions of Theorem 1 and indeed lies in a special class of augmenting functions that form a practical set from which one can tailor make an augmenting function for a given problem.

2.5 Positive Basis

A subset of reasonable augmenting functions may be defined by using a positive basis $\{\text{\bf n}_{1},\dots,\text{\bf n}_{l}\}$ to scalarise the deviations $u\in\mathbb{R}^{m}$. (Note that for our purposes, $m=n_{x}\,|S|$.) Such deviations can be associated, for example, with the satisfaction of linear inequalities, where we might have $u=b-Ax$ given a constraint $Ax\leq b$.

Definition 4

We say a set of vectors $\left\{\text{\bf n}_{1},\dots,\ \text{\bf n}_{l}\right\}$ where ${m}+1\leq l\leq 2{m}$ is a positive basis for $\mathbb{R}^{m}$ if and only if every $u\in\mathbb{R}^{{m}}$ can be expressed as a positive combination of these vectors, i.e., there exists $\alpha_{i}\geq 0$ for $i=1,\dots,l$ for which $u=\sum_{i=1}^{l}\alpha_{i}\textbf{\emph{n}}_{i}.$

The following property of positive bases will be useful in the developments which follow.

Theorem 5

([12] Theorem 3.1) $\left\{\textbf{\emph{n}}_{1},\dots,\textbf{\emph{n}}_{l}\right\}$ positively spans $\mathbb{R}^{m}$ if and only if for every non-zero $u$ there exists an index $i$ such that $u\cdot\textbf{\emph{n}}_{i}>0$.

Let $e_{i}$, $i=1,\dots,m$ represent the elementary unit vectors of $\mathbb{R}^{m}$ with entry $i$ set to one and all other entries set to zero. Examples of positive bases on $\mathbb{R}^{m}$ include:

The vertices of a ${m}$-simplex (generalised tetrahedron), centred at the origin.

• 2.

The set of vectors $\{+e_{i}\}^{m}_{i=1}\cup\{\sum_{i=1}^{m}-e_{i}\}$

• 3.

The set of vectors $\{\pm e_{i}\}^{m}_{i=1}$

2.6 Norm-Like Augmenting Functions

As noted in both [8] and [13], norms are viable augmenting functions with appealing theoretical support for overcoming duality gaps. However, norms have the disadvantage that they uniformly penalise constraint violations, which results in a loss of flexibility and precision in the fine-tuning of the penalisation. This limitation motivates the following development of asymmetrical but norm-like augmenting functions. The polyhedral norms $\|\ {\cdot}\ \|_{\infty}$ and $\|\ {\cdot}\ \|_{1}$ may be represented using the positive basis $\{\pm e_{i}\}^{m}_{i=1}$ in the following ways:

[TABLE]

The representation in (19) relies on each vector in the basis having a negative multiple which is also in the basis. The representations in (17) and (18) do not have this limitation, and may be generalised to any positive basis $N:=\{\text{\bf n}_{1},\dots,\text{\bf n}_{l}\}$ as follows:

[TABLE]

The functions $\psi_{\infty}^{N}$ and $\psi_{1}^{N}$ are not necessarily norms, but do share some useful properties with norms. Specifically, these functions are positive homogeneous (which implies that they vanish at zero), strictly positive for all $u\in\mathbb{R}^{m}\setminus\{0\}$, finite valued, sub-additive and coercive.

The proposed augmenting function $\psi_{\rho}(u)$ given in (16) may be represented in the form of (21), using the positive basis $N_{\rho}=\{\overline{\rho}_{s,i}e_{i+(s-1)n_{x}}\mid s\in S,i\in\{1,\dots,{n_{x}}\}\}\cup\{-\underline{\rho}_{s,i}e_{i+(s-1)n_{x}}\mid s\in S,i\in\{1,\dots,{n_{x}}\}\}$, as follows:

[TABLE]

Lemma 6

If two functions $\psi_{A}$ and $\psi_{B}$ are positive homogeneous, strictly positive for all $u\neq 0$, and are finite valued then there exists a finite $\gamma>0$ such that

[TABLE]

Proof. Since they are positive homogeneous, $\psi_{A}$ and $\psi_{B}$ vanish at zero and so (23) trivially holds with equality at $u=0$. To obtain the required inequality in (23) for nonzero $u$, set $V=\{u:||u||=1\}$ (where $||\cdot||$ is any norm) and take $\gamma=\min_{u\in V}\frac{\psi_{A}(u)}{\psi_{B}(u)}$. Since $\psi_{A}(u)$ and $\psi_{B}(u)$ are strictly positive and finite for all $u\neq 0$, and $V$ is compact, this minimum exists and $\gamma$ is strictly positive and finite. For any point $u\in\mathbb{R}^{m}\setminus\{0\}$, $||u||$ is strictly positive and the point $\frac{u}{||u||}$ is in $V$. The required inequality follows from the positive homogeneity of $\psi_{A}$ and $\psi_{B}$:

[TABLE]

Proposition 7

For any positive basis $N$, the augmenting functions $\psi_{\infty}^{N}$ and $\psi_{1}^{N}$ given in (20) and (21) respectively satisfy the conditions given in Theorem 1.

Proof.

Let $V=B^{\infty}_{\varepsilon}(0)$ be an open ball in the infinity norm with radius $\varepsilon>0$ centred at the origin. This is an appropriate open neighbourhood of [math] for the purposes of Conditions 2 and 3 of Theorem 1.

**Condition 1: ** $\psi(0)=0$.

If $u=0$ then $\text{\bf{n}}_{i}^{\top}u=0$ and therefore $\psi_{\infty}^{N}(u)=0$ and $\psi_{1}^{N}(u)=0$, as required.

**Condition 2: ** $\psi(u)\geq\delta>0,\forall u\not\in V$ for some positive scalar $\delta$.

Using Theorem 5, for any $u\neq 0$ we have some $i$ such that $\text{\bf{n}}_{i}^{\top}u>0$ and hence $\psi_{\infty}^{N}(u)>0$. Now define

[TABLE]

where $\delta>0$ follows from the compactness of the $\varepsilon$- ball, the continuity of $u\mapsto\max_{i=1,\dots,l}\{\text{\bf{n}}_{i}^{\top}u\}$, and Theorem 5. For any $u\notin V$, the point $v:=\varepsilon\frac{u}{\|u\|_{\infty}}$ is in $V$ and hence $\psi_{\infty}^{N}(v)\geq\delta>0$. Using the positive homogeneity property we have

[TABLE]

using the fact that $u\notin V$ means $\|u\|_{\infty}\geq\varepsilon$. This is the required inequality for $\psi_{\infty}^{N}$.

Apply Lemma 6 to deduce that there exists a $\eta>0$ such that:

[TABLE]

$\eta\delta$ is also a positive scalar and so this is the required inequality for $\psi_{1}^{N}$.

**Condition 3: ** $\psi(u)\geq\gamma||u||_{\infty},\forall u\in V$ for some positive scalar $\gamma$.

The property holds trivially for $u=0$. For any $u\in V\setminus\{0\}$, the point $v:=\varepsilon\frac{u}{\|u\|_{\infty}}$ is in $V$ and using the same $\delta$ as defined in (24) we have

[TABLE]

and so we may place $\gamma:=\frac{\delta}{\varepsilon}>0$. This is the required inequality for $\psi_{\infty}^{N}$.

As above, apply Lemma 6 to deduce that there exists a $\eta>0$ such that:

[TABLE]

$\eta\gamma$ is also a positive scalar and so this is the required inequality for $\psi_{1}^{N}$.

Corollary 8

Assume that $\zeta^{SIP}$ is feasible, its optimal value is finite, and the data which defines it is rational. Then the optimal value of the augmented Lagrangian dual problem $\zeta^{LD+}_{\rho}$ using an augmenting function of the form of (20) or (21) is equal to the optimal value of $\zeta^{SIP}$ for some finite $\rho$; that is,

[TABLE]

where $\omega_{LP}$ is the optimal multiplier of the linear programming relaxation of (5). In particular, this applies to our proposed augmenting function (16).

Proof. The equalities (25) follow directly from Theorem 1 and Proposition 7. The last claim follows from the observation that (16) may be represented as a function of the form of (21), as demonstrated in (2.6).

Remark 9

Consider a positive basis $N=\left\{\textbf{\emph{n}}_{1},\dots,\textbf{\emph{n}}_{l}\right\}$. Each of the functions $g_{i}(u)=\max\{\textbf{\emph{n}}_{i}^{\top}u,0\}$ is non-negative, positive homogeneous and finite valued, and these properties are preserved if multiple $g_{i}$s are summed, or their maximum is taken. By Theorem 5, for any non-zero $u$ there exists an index $i\in\{1,\dots,l\}$ such that $g_{i}(u)$ is strictly positive. Therefore, if every one of the $g_{i}$ functions is combined using a combination of summation and/or maximisation, the resulting function $g(u)$ will be strictly positive for all non-zero $u$. Applying Lemma 6 to bound $g$ below by a positive multiple of $\psi_{\infty}^{N}$ (as $\psi_{1}^{N}$ was treated in Proposition 7) shows that this function $g(u)$ satisfies the conditions of Theorem 1, and as such will close the duality gap if used as an augmenting function.

Remark 9 implies that we can not only construct $\psi_{\infty}^{N}$ and $\psi_{1}^{N}$ but also a wide variety of other augmenting functions from any given positive basis, depending on the order in which the maximisation and summation operations are applied to the $g_{i}$ functions.

Furthermore, the sum or maximum of any two augmenting functions which satisfy the conditions of Theorem 1 will itself satisfy the conditions of Theorem 1, which yields further flexibility.

Remark 10

By using the positive basis $\{+e_{i}\}^{m}_{i=1}\cup\{\sum_{i=1}^{m}-e_{i}\}$ or similar to define an augmenting function, we can obtain penalty terms analogous to the Lagrangian terms obtained through surrogate semi-Lagrangian relaxation.

3 Developing a penalty-based block Gauss-Seidel method

To exploit the potential for decomposability that this formulation presents, we consider a block Gauss-Seidel (GS) method approach. In Section 3.1, we present a classical framework for GS methods as applied to nonlinear optimisation problems, and in Section 3.2, we show how it can be adapted to obtain solutions for SMIP.

3.1 A block Gauss-Seidel method

We consider the general problem given by

[TABLE]

We assume that $f$ is convex, but not necessarily differentiable. The sets $X$ and $Z$ are closed, but not necessarily convex. The assumption that $x\in X$ and $z\in Z$ are taken from disjoint sets is adequate for our purposes, although block GS approaches have been studied in a more general setting where $(x,z)$ is taken from a set $K\subset\mathbb{R}^{n_{x}+n_{z}}$. (See [34] for a treatment of the case where the constraint set is disjoint, and [17] for the case where the constraint set may not be disjoint and developments are based on biconvexity assumptions.) GS methods solve problem (26) by separating it into two simpler problems. Given an iterate $(x^{k},z^{k})$, problem (26) is solved with respect to $x$ for fixed $z=z^{k}$, yielding a new $x$-iterate $x^{k+1}$. Then, problem (26) is solved with respect to $z$ for fixed $x=x^{k+1}$, yielding a new $z$-iterate $z^{k+1}$. In Algorithm 1, a formal listing of a block GS method applied to problem (26) is given.

The sequence $\left\{(x^{k},z^{k})\right\}$ generated by iterations of Algorithm 1 has limit points when $X$ and $Z$ are compact. When $f$ is furthermore continuous and bounded from below over $X\times Z$, the limit points $(x^{*},z^{*})\in X\times Z$ are easily shown to be partial minima [34]; that is, it holds that

[TABLE]

This claim is formally stated in Proposition 11, whose proof is implicit from the developments of [34] and is given here for the sake of completeness.

Proposition 11

For problem (26), let $f$ be continuous and bounded from below, and let $X$ and $Z$ be compact. Then the limit points $(x^{*},z^{*})$ of the sequence $\left\{(x^{k},z^{k})\right\}$ generated by iterations of Algorithm 1 are partial minima.

Proof. We have by construction that $f(x^{k},z^{k})\leq f(x^{k},z)$ for all $z\in Z$, and by the continuity of $f$, we have the second requirement $f(x^{*},z^{*})\leq f(x^{*},z),\ \forall z\in Z$ for partial optimality. To establish the first requirement (27), assume for sake of contradiction that there is an $\bar{x}\in X$ for which $f(x^{*},z^{*})>f(\bar{x},z^{*})$. Due to the continuity of $f$, we have, for some infinite subsequence index set $\mathcal{K}$ such that $\lim_{k\to\infty,k\in\mathcal{K}}(x^{k},z^{k})=(x^{*},z^{*})$, the existence of $\gamma>0$ such that $f(x^{k},z^{k})-f(\bar{x},z^{k})>\gamma>0$. Thus, $f(x^{k},z^{k})>f(\bar{x},z^{k})+\gamma\geq f(x^{k+1},z^{k})+\gamma\geq f(x^{k+1},z^{k+1})+\gamma$, which would imply that $\lim_{k\to\infty}f(x^{k},z^{k})=-\infty$ since $\mathcal{K}$ is an infinite index set, so that $f$ is unbounded from below, a contradiction. Therefore, $(x^{*},z^{*})$ must be a partial minimum for problem (26).

Remark 12

For practical purposes, we might approximate the satisfaction of (27) and (28) through the $\epsilon\geq 0$ parameterised termination criterion

[TABLE]

In the setting where $f$ is convex and differentiable, $X$ and $Z$ are nonempty, closed and convex, and $(x,z)\mapsto f(x,z)$ is inf-compact, it is well-known (see, for example, [6, 18, 31]) that the limit points $(x^{*},z^{*})$ are optimal for problem (26). However, in the more general setting where $f$ is non-differentiable and/or $X$ and $Z$ are non-convex, it is well-known that a partial minimum need not be a global, or even a local, minimum. In what follows, we provide a few small examples to illustrate this suboptimal stabilisation, and to motivate heuristic features of our developed algorithm that can mitigate this unfortunate tendency.

Examples:

Let problem (26) be specified so that $f(x,z):\mathbb{R}\times\mathbb{R}\mapsto\mathbb{R}$ is defined to be $f(x,z)=7x^{2}+10xz+7z^{2}$, and let $X=Z=\left\{-2,-1,0,1,2\right\}$. For $(x^{0},z^{0})=(2,-2)$, the application of Algorithm 1 leads immediately to the one limit point $(\bar{x},\bar{z})=(1,-1)$. We have $f(1,-1)=4$, but $f(0,0)=0$, so $(\bar{x},\bar{z})=(1,-1)$ is not optimal. Note here that $f$ is convex and continuously differentiable, but the constraint set $X\times Z$ is nonconvex due to the integer restriction, and this is the reason that the limit point was not guaranteed to be optimal. 2. 2.

Let problem (26) be specified so that $f(x,z):\mathbb{R}\times\mathbb{R}\mapsto\mathbb{R}$ is defined to be $f(x,z)=-2x-z+\rho\left|x-z\right|$, $X=[-2,3]$, and $Z=[0,5]$. For $\rho\in[0,1)$, the optimal solution is $(x^{*},z^{*})=(3,5)$. For $\rho=1$, the optimal solutions are taken from $(x^{*},z^{*})\in\left\{3\right\}\times[3,5]$, and for $\rho>1$, the optimal solution is $(x^{*},z^{*})=(3,3)$.

1. (a)

When applying the GS approach of Algorithm 1 with $\rho\in(0,1)$, the resulting sequence stabilises after one iteration at the optimum $(x^{*},z^{*})=(3,5)$ for any feasible starting point. 2. (b)

For $\rho=1$ with $z^{0}\geq 3$, we have after half an iteration $(x^{1},z^{0})=(3,z^{0})$ which is an optimum solution, and the remaining updates stay at some optimal solution $(x^{*},z^{*})\in\left\{3\right\}\times[3,5]$. For $\rho=1$ with starting point $z^{0}<3$, we have $x^{1}=3$ and $z^{1}\in[3,5]$ and so stabilisation at an optimal solution also occurs. 3. (c)

For $\rho>1$ with $z^{0}\geq 3$, we have $(x^{1},z^{1})=(3,3)$, which is optimal. However, for $\rho>1$ with $z^{0}<3$, we have $x^{1}=z^{0}$ and $z^{1}=z^{0}$, so that stabilisation occurs at $(\bar{x},\bar{z})=(z^{0},z^{0})$, which is not optimal. 3. 3.

Let problem (26) be specified so that $f(x,z):\mathbb{R}^{3}\times\mathbb{R}^{3}\times\mathbb{R}^{3}\mapsto\mathbb{R}$ is defined to be

[TABLE]

and let $X$ and $Z$ be defined so that

[TABLE]

and $Z=\left\{0,1\right\}^{3}$. For $\rho\to\infty$ (simulating the enforcement of constraints $x_{i,j}=z_{j}$ for $i=1,2$ and $j=1,2,3$) we have the optimal solution

[TABLE]

If such constraints are altogether ignored ($\rho=0$), then the optimal $x$-component is $x^{*}=\left((0,0,1),(1,0,0)\right)$. This behaviour would only change at the threshold $\rho=1$. For $\rho>1$, the optimal solution would be $(x^{*},z^{*})=([(0,1,0),(0,1,0)],(0,1,0))$.

1. (a)

Now we consider what happens when the GS approach of Algorithm 1 is applied. Let $z^{0}=(0,0,0)$. Starting with a small penalty such as $\rho=0.5$, we have

[TABLE]

where there is more than one way to choose $z^{1}$. If, for example, we make it a policy to choose $z$ by some bitwise lexicographical rule, then we choose $z^{1}=(0,0,0)$. Keeping this same penalty $\rho=0.5$, we find that stabilisation has occurred, where $x^{k}=x^{1}$ and $z^{k}=z^{1}$ for $k\geq 1$. If we increase the penalty value to $\rho=2$ for iteration $k=2$, then we have the stabilisation $x^{1}=\left((0,0,0),(0,0,0)\right)$ and $z^{1}=(0,0,0)$, which is suboptimal (and $\rho=2$ is the threshold for this change in stabilisation to occur).

If, instead, the $z$ update is chosen by a reverse-lexicographic rule, so that $z^{1}=(1,0,1)$, then we have immediate stabilisation with

[TABLE]

for all $k\geq 1$ for all $\rho>0$. (Notice that no matter how large the penalty is, consensus is not achieved in the GS setting. That is, without additional restriction on how $z$ is updated, the optimal $z$ update may be chosen to always correspond to a consensus solution that is infeasible for both scenarios. In practice, we would need a rule to insure that the $z$ update is chosen to satisfy $\sum_{j=1}^{3}z_{j}\leq 1$ to match with the constraints in the $x$ update subproblems.) 2. (b)

The shortcomings of the above GS approach motivate the introduction of more precision in how the consensus discrepancies are penalised, where $f$ is redefined to be

[TABLE]

That is, instead of one scalar $\rho$, we have term-specific $\rho_{i,j}>0$ for each $i=1,2$ and $j=1,2,3$. We start as before with $z^{0}=(0,0,0)$, and let $\rho_{i,j}=0.5$ for each $i=1,2$ and $j=1,2,3$. Assuming lexicographic rule in choosing $z$, we have as before

[TABLE]

and this is stable if the penalty does not change. Now increase $\rho_{1,3}=\rho_{2,1}=1$, and we have

[TABLE]

and this is stable if the penalty does not change. Increasing $\rho_{1,1}=\rho_{2,3}=1$, we have again

[TABLE]

and this is stable. But once we again increase $\rho_{1,3}=\rho_{2,1}=2$, we have

[TABLE]

which is optimal for the original problem.

The last example suggests that there may be no fixed ideal penalty in a GS setting that will lead to both a closing of the duality gap and avoiding the nonoptimal stationarity due to GS iterations. The penalty must vary in a manner that takes the component-wise consensus status into consideration. Any sensible heuristic for varying the penalties would have all penalties start small (but nonzero), and increase carefully, in a ”fine-tuned” manner so as to “suggest” a temporary fixing of certain components of $x$ to the current fixed values of the corresponding components of $z$. The strength of suggestion for each component is always relative to the other components as the magnitude of each component-wise penalty is relative to the magnitude of the other component-wise penalties.

An approach based on such an idea where some subset of variables is subject to “suggested” fixing with strength of suggestion determined by the penalties would be of a “soft” combinatorial nature. This is in contrast with a “hard” combinatorial approach that might be based on the idea of choosing some subset of integer variables at each iteration and simply fixing each one to some constant feasible value while conducting a minimisation over the unfixed variables. The algorithm to be presented later is of a soft combinatorial nature.

3.2 Adapting block Gauss-Seidel method to solve SMIPs using Penalty functions

In this section we present how block GS method can be used to obtain solutions for SMIP problems. The approach will rely on the delayed calculation of variable $z$, which will in turn allow us to obtain a decomposed version of the problem. To do such, let us first explicitly state $\zeta^{LR+}_{\rho}(\omega)$ as

[TABLE]

The following proposition will become useful in the following derivations.

Proposition 13

For problem (30)-(32) with any $\omega=\tilde{\omega}\in\mathbb{R}^{n_{x}\times|S|}$, there exists a finite $\rho^{*}(\tilde{\omega})$ such that $\zeta^{LR+}_{\rho^{*}}(\tilde{\omega})=\zeta^{SIP}$.

Proof. The penalty terms in (30) result from the evaluation of $\psi_{\rho}\left((x_{s}-z)_{s\in S}\right)$ with $\psi_{\rho}$ as defined in (16). Thus, by Corollary 8, the requirements of Theorem 1 are satisfied. Now one can rely on Remark 3 to make a free choice of $\omega$.

Proposition 13 enables us to make the choice of $\omega=0$, which leads to

[TABLE]

The block GS method for solving $\zeta^{LR+}_{\rho}(0)$ proceeds as follows. Let

[TABLE]

where

[TABLE]

and $(\underline{\mu}_{s},\overline{\mu}_{s}):=(\frac{1}{p_{s}}\underline{\rho}_{s},\frac{1}{p_{s}}\overline{\rho}_{s})$ for each $s\in S$. For a given $\rho^{k}_{s}=(\underline{\rho}^{k}_{s},\overline{\rho}^{k}_{s})_{s\in S}$ and an initial $z^{0,0}$, the proposed method will iterate between the solution of following $l=0,1,\dots,l_{\text{max}}$ subproblems:

[TABLE]

and

[TABLE]

followed by $l=l+1$ and successive repetition until partial convergence is approximately achieved in the sense of (29). In this context, partial convergence can be interpreted as having

[TABLE]

given a threshold $\epsilon\geq 0$.

At last, if the current primal infeasibility level, giving by a residual measure such as $||x^{k,l}-z^{k,l}||_{2}^{2}$, is not acceptable for a $\epsilon$ threshold, the set of penalties $\rho^{k}=(\underline{\rho}^{k},\overline{\rho}^{k})$ are then updated to $\rho^{k+1}=(\underline{\rho}^{k+1},\overline{\rho}^{k+1})$ and the process is repeated for iteration $k+1$.

4 Computational Implementation Aspects

Two remarkable features can be exploited in the design of an algorithm based on this idea. First, scenario-wise separability is straightforwardly achieved in the calculation of $(x^{k,l+1},y^{k,l+1})$. This means that instead of solving one large mixed-integer linear programming (MILP) problem in this update step, we can solve several small MILP problems instead, which is typically more efficient due to the exponential nature of the branch-and-cut-based methods used to solve them.

To formulate the $(x^{k,l},y^{k,l})$-update

[TABLE]

it is necessary to explicitly represent the function $[\ \cdot\ ]^{-}$. To do so, we consider an equivalent reformulation of the problem given by

[TABLE]

Second, the calculation of

[TABLE]

may be performed by computing

[TABLE]

where the penalty function $(x,z)\mapsto\zeta^{\rho}(x,z,\rho)$ is defined by

[TABLE]

The last displayed problem can be solved using the following equivalent mathematical programming formulation:

[TABLE]

When the $x$ components are all restricted to take binary values, it is possible to show that the calculation of $z^{k,l+1}$ can be performed in the following closed form where each component of $z^{k,l+1}$ always takes binary value. In that case, its optimal solution is given by

[TABLE]

The cases in which we have a tie might require ”flipping a coin” for deciding on the value for $z^{k,l+1}$, as it becomes a case of multiple minima. The existence of multiple minima can be better understood from the following explicit form of the solution for the general case. In the following proposition, we assume $Z$ is a closed convex set, so that no explicit integrality constraints are enforced.

Proposition 14

Suppose a set of scenario dependent solutions $(x_{s})_{s\in S}$, where $x_{s}=(x_{s,i})_{i=1,\dots,n_{x}}$, are given and $z:=(z_{i})_{i=1,\dots,n_{x}}$. For each $i\in\{1,\dots,n_{x}\}$ define

[TABLE]

Then $z_{i}$ solves problem (37) given fixed $(x_{s})_{s\in S}$ if and only if

[TABLE]

Proof. The index $s$ term of the penalty function $\zeta^{\rho}$ may be written as

[TABLE]

As this is separable in the variables $(z_{1},\dots,z_{n_{x}})$, its subdifferential is defined as the cross product of intervals, one for each component $i$. Thus, the necessary and sufficient condition

[TABLE]

can be equivalently stated as

[TABLE]

for each $i=1,\dots,n_{x}$, which is given by:

[TABLE]

which in turn is equivalent to (42).

We now consider how to update the penalty parameters $\rho^{k}$. A simple strategy is

[TABLE]

By doing so, we are reinforcing the penalties associated with the respective discrepancies. In other words for each $i=1,\dots,n_{x}$:

[TABLE]

[TABLE]

Remark 15

The update in $\rho^{k+1}$ has the effect of changing the left hand side of (42) at the next iteration by the amount:

[TABLE]

*for each $i=1,\dots,n_{x}$. If the addition of this factor ensures the sum in left hand side of (42) at iteration $k+1$ exits the interval $\left[-\sum_{s\in I^{0}(z_{i})}\overline{\rho}_{s,i},\sum_{s\in I^{0}(z_{i})}\underline{\rho}_{s,i}\right]$ associated with the prior choice of $z^{k}_{i}=x^{k}_{s,i}$ then we would be forced to choose new consensus values $z^{k}_{i}$ in order to re-establish the satisfaction of the optimality condition (42). In doing so, a reassignment of the index sets $I^{+}(z_{i}^{k})$, $I^{-}(z_{i}^{k})$, and $I^{0}(z_{i}^{k})$ is effected. As intuition would suggest, the optimality condition (42) is more easily satisfied when $s\in I^{0}(z_{i}^{k})$ for large $\underline{\rho}_{s,i}$ and $\overline{\rho}_{s,i}$, as this makes the target interval $\left[-\sum_{s\in I^{0}(z_{i})}\overline{\rho}_{s,i},\sum_{s\in I^{0}(z_{i})}\underline{\rho}_{s,i}\right]$ larger. *

To effect a gradual increase in the terms $\Delta^{k}$ in an attempt to improve convergence with the satisfaction of the NAC, we considered an increasing multiplier factor to $\psi_{\rho}$ given by $\beta^{(k-1)}-1$ (where $(k-1)$ represents an exponent and not an iteration index). In other words, we consider the objective at a given iteration $k$ as being

[TABLE]

Combining what have been exposed so far, one first algorithmic approach consists of the following setting presented in Algorithm 2.

5 Experimental setting

In this section we describe the computational experiments performed to assess the performance of the proposed approach. To evaluate the performance of the proposed method, we tested its efficacy on three distinct classes of problems from literature, namely the capacitated facility location problems (CAP) from [7], the dynamic capacity allocation problems (DCAP) available in [1], and the server location under uncertainty problems (SSLP) first introduced in [26]. To provide a more solid base of comparison, 50 random instances of two problems from each class were generated.

The CAP problems are two-stage SMIP problems with pure binary first- and second-stage variables arising in the context of network design problems. We selected the instances coded as 101 and 111 in [7], considering random samples of 100 scenarios from a list of 5000 scenarios available.

The DCAP problems are two-stage SMIP problems arising in dynamic capacity acquisition and allocation under uncertainty. All problem instances have mixed-integer first-stage variables and pure binary second-stage variables. We selected the instances coded as 233 and 342 (which encodes the number of resources, tasks, and periods, respectively), considering random samples of 100 scenarios from the original 500 available.

The SSLP problems are two-stage SMIP problems arising in server location under uncertainty. The problems have pure binary first-stage variables and mixed-binary second-stage variables. We considered the instances coded as 5-50 and 10-50 (which encode the number or servers and the number of clients, respectively) with 100 scenarios that were randomly generated according to the guidelines provided in [26].

To compare and benchmark the performance of the proposed approach against a known quantity we have implemented the Progressive Hedging (PH) algorithm, which was originally proposed by [29] and, as previously discussed, has been widely used as an heuristic approach to solve SMIP problems. The PH algorithm is stated in Algorithm 3 for the sake of completeness. In this algorithm,

[TABLE]

which means that Line 8 comprises the solution of $|S|$ mixed-integer quadratic subproblems at each iteration $k$. The analogous step in PGBS algorithm requires us to solve typically less difficult MIP problems instead.

Another advantage of PGBS in the context of SMIP problems is that $z^{k,l+1}$ tends to (in most cases) satisfy the integrality constraints of the problem. This is in contrast with the consensus value computed with the averaging of PH (Line 10 in Algorithm 3), which tends to steer the consensus value away from integral values. In the case of binary variables, the PH averaging computation of the consensus $z^{k,l+1}$ is especially prone to producing many fractional valued components which can lead to episodic cycling in binary values set in the assignment of scenario specific variables.

In the PGBS experiments, the parameters were chosen from $\beta\in\{1.25,1.11\}$, and $\gamma\in\{0.5\rho^{0},\rho^{0}\}$. Three different initial values for $\rho^{0}$ were used in both the PGBS and PH experiments. In the Progressive Hedging algorithm, dual multipliers were initialised as [math] and the penalty parameter ($\rho$) was set to $\rho=\rho^{0}$. The parameter $z^{0,0}$ has been initialised according to the solution $x^{0}_{s}$ (from Line 3 in Algorithm 2 and Line 3 in Algorithm 3). For Algorithm 2, we initialised $z^{0,0}$ as being

[TABLE]

for all components restricted to be integer variables, where $\lceil\ \cdot\ \rfloor$ denotes rounding operation. For the components without integrality restrictions, we have dropped the rounding operator. In case of PH (Algorithm 3), all components were calculated by dropping the rounding operator.

As both CAP and SSLP problems have pure binary first-stage variables, we have used (41) to perform the step depicted in Line 12 of Algorithm 2. For DCAP, we relied on solving (38)-(40) explicitly.

A time limit of 1000 seconds and termination condition of $\epsilon=10^{-3}$ was used for both methods. A total of 300 ($3\times 2\times 50$) instances were solved with three parameter choices for PH (different choices of $\rho^{0}$) and 12 combinations of parameter choices for PBGS (different choices of $\rho^{0}$, $\beta$ and $\gamma$). The computational experiments were performed on a Intel i7 CPU with 3.40GHz and 8GB of RAM. All methods have been implemented using AIMMS 3.14 and all subproblems have been solved using CPLEX 12.6.3 with its standard configuration.

5.1 Numerical results

A summary of the computational results is presented in Figures 1 to 3, which depicts the average computational time and objective value difference for the 50 instances considered for both PH and PBGS in all parameter settings that have been tested.

The blue bars indicate the average wall clock times for both methods. We highlight that the instances in which PH terminated due to the time limit of 1000s have been removed from the average calculations, these being treated as outliers. The green line shows the average objective value relative difference, which is calculated as

[TABLE]

where $z^{i}_{PBGS}$ and $z^{i}_{PH}$ are the objective function values obtained for the solutions returned by PBGS and PH for instance $i$, respectively, and $N$ is the total number of instances considered for average value calculations. To obtain $z^{i}_{PBGS}$ and $z^{i}_{PH}$, we used the last solution returned by both methods and evaluated it a posteriori. For the cases in which PH returned solutions that were infeasible in regard to integrality restrictions (typically those obtained when the algorithm stopped due to the time criterion), rounding has been performed to recover a feasible solution to be evaluated when applicable.

For the CAP instances, all configurations tested with PBGS and PH presented similar values for the objective function, and in most configurations PBGS presented better performance in terms of computational time. For the DCAP instances, in all cases PH terminated due to the time limit of 1000 seconds. For these problems, a comparison in terms of objective function shows that the differences between the objective function value of the solutions found by PGBS and PH are more pronounced. A similar behaviour can be observed in the SSLP instances, in which PBGS outperforms PH in terms of solution times in most cases while providing solutions that are, in the worst case, 0.5% worse for SSLP5-50 and 5% worse for SSLP10-50. In the Appendix we present a detailed summary of the statistics for each of the problems, including the fraction of the runs in which PH was not able to converge within the specified time limit. Overall, PBGS seems to be able to obtain comparably good solutions however presenting more reliable convergence behaviour.

6 Conclusions

In this paper we have presented an alternative approach for solving stochastic mixed-integer problems based on the combination of penalty-based and block Gauss-Seidel methods. The motivation of such arises from recents theoretical results that motivates the consideration of Lagrangian-based methods under alternative perspectives to approach such problems.

The computational experiments performed suggest that there is potential for exploiting this framework as it allows the development of a competitive approach in terms of computational efficient. It is worth highlighting that the methodology developed is readily amenable to parallelisation, which is a key point for dealing with large-scale SMIPs.

Further developments of this research could be classified under two distinct standpoints. Under a theoretical perspective, suitable alternative extensions of the block Gauss-Seidel approach into non-smooth non-separable problems are worth investigation. A better understanding of how to fine-tune the updates of the penalty coefficients would improve the likelihood (or perhaps even guarantee!) that the block Gauss-Seidel iterations do not display suboptimal stationarity. This would improve the trend of the objective values computed by the main algorithm. In terms of practical considerations, it would be of interest to evaluate the performance of the proposed approach in contexts other than SMIPs and considering its extension to the multi-stage case.

Acknowledgements

The authors would like to acknowledge the support provided by the Australian Research Council (ARC) grant ARC DP140100985.