Smoothing Techniques and Augmented Lagrangian Method for Recourse Problem of Two-Stage Stochastic Linear Programming

Abstract

The augmented Lagrangian method can be used for solving recourse problems and obtaining their normal solution in solving two-stage stochastic linear programming problems. The augmented Lagrangian objective function of a stochastic linear problem is not twice differentiable which precludes the use of a Newton method. In this paper, we apply the smoothing techniques and a fast Newton-Armijo algorithm for solving an unconstrained smooth reformulation of this problem. Computational results and comparisons are given to show the effectiveness and speed of the algorithm.

1. Introduction

In stochastic programming, some data are random variables with specific possibility distribution [1], which was first introduced by the designer of linear programming problems, Dantzig, in [2].

In this paper, we consider the following two-stage stochastic linear program (slp) with recourse which involves the calculation of an expectation over a discrete set of scenarios: min_𝑥∈𝑋𝑓(𝑥)=𝑐^𝑇𝑥+𝜙(𝑥),𝑋=𝑥∈ℝ^𝑛∶𝐴𝑥=𝑏,𝑥≥0, (1) where 𝜙(𝑥)=𝐸(𝑄(𝑥,𝜔))=^𝑁_𝑖=1𝑄𝑥,𝜔^𝑖𝜌𝜔^𝑖 (2) and 𝐸 shows the expectation of function 𝑄(𝑥,𝜔) which depend on the random variable 𝜔. The function 𝑄 is defined as follows: 𝑄(𝑥,𝜔)=min_{𝑦∈ℝ𝑛2}𝑞(𝜔)^𝑇𝑦∣𝑊(𝜔)^𝑇𝑦≥ℎ(𝜔)−𝑇(𝜔)𝑥, (3) where 𝐴∈ℝ^𝑚×𝑛, 𝑐∈ℝ^𝑛, and 𝑏∈ℝ^𝑚. Also, in the problem (3) vector of coefficients 𝑞(⋅)∈ℝ^𝑛2, matrix of coefficients 𝑊^𝑇(⋅)∈ℝ^𝑚2×𝑛2, demand vector ℎ(⋅)∈ℝ^𝑚2, and matrix 𝑇(⋅)∈ℝ^𝑚2×𝑛 depend on the random vector 𝜔 with support space Ω. The problems (1) and (3) are called master and recourse problems of stochastic programming, respectively.

We assume that the problem (3) has a solution for each 𝑥∈𝑋 and 𝜔∈Ω.

In general, the recourse function 𝜙(𝑥) is not differentiable everywhere. Therefore, the traditional methods use nonsmooth optimization techniques [3–5]. However, in the last decade, it is proposed smoothing method for recourse function in standard form of recourse problem [6–11]. In this paper, we apply a smooth approximation technique to smooth recourse function that the recourse problem has inequality linear constrained. For more explanation see Section 2. The approximated problem is based on the least two-norm solution of recourse problem. This paper considers the augmented Lagrangian method to obtain least two-norm solution (Section 3). For convenience, Euclidean least two-norm solution of linear programming problem is named normal solution. This effective method contains solving an unconstrained quadratic problem which its objective function is not twice differentiable. To apply a fast Newton method we use the soothing technique and replace plus function by an accurate smooth approximation [12, 13]. In Section 4, the smoothing algorithm and the numerical results are presented. Also, concluding remarks are given in Section 5.

We now describe our notation. Let 𝑎=[𝑎_𝑖] be a vector in ℝ^𝑛. By 𝑎₊ we mean a vector in ℝ^𝑛 whose 𝑖th entry is 0 if 𝑎_𝑖<0 and equals 𝑎_𝑖 if 𝑎_𝑖≥0. By 𝐴^𝑇 we mean the transpose of matrix 𝐴, and ∇𝑓(𝑥₀) is the gradient of 𝑓 at 𝑥₀. For 𝑥∈ℝ^𝑛, ‖𝑥‖ and ‖𝑥‖_∞ denote 2-norm and infinity norm, respectively.

2. Approximation of Recourse Function

As mentioned the objective function of (1) is nondifferentiable. This disadvantage property occurs on the recourse function. In this section, there is an attempt to approximate it to a differentiable function.

Using dual of the problem (3), function 𝑄(𝑥,𝜔) can be written as follows: 𝑄(𝑥,𝜔)=max_{𝑧∈ℝ𝑚2}(ℎ(𝜔)−𝑇(𝜔)𝑥)^𝑇𝑧s.t.𝑊(𝜔)𝑧=𝑞(𝜔),𝑧≥0. (4) Unlike the linear recourse function, the quadratic recours function is differentiable. Thus in this paper, the approximation is based on the following quadratic problem with helpful properties: 𝑄_𝜖(𝑥,𝜔)=max_{𝑧∈ℝ𝑚2}(ℎ(𝜔)−𝑇(𝜔)𝑥)^𝑇𝜖𝑧−2‖𝑧‖²s.t.𝑊(𝜔)𝑧=𝑞(𝜔),𝑧≥0. (5) The next theorem shows that, for the sufficiently small 𝜖>0, the solution of this problem is the normal solution of the problem (4).

Theorem 1. For functions 𝑄(𝑥,𝜔) and 𝑄_𝜖(𝑥,𝜔) introduced in (4) and (5), the following can be presented:(a)∃𝜖>0 such that, for each 𝜖∈(0,𝜖], the solution for the problem (5) is the normal solution for the problem (4).(b) For each 𝜖>0, function 𝑄_𝜖(𝑥,𝜔) is differentiable with respect to 𝑥.(c) The gradient of function 𝑄_𝜖(𝑥,𝜔) at point 𝑥 is ∇𝑄_𝜖(𝑥,𝜔)=−𝑇^𝑇(𝜔)𝑧^∗_𝜖(𝑥,𝜔) (6) in which 𝑧^∗_𝜖(𝑥,𝜔) is the solution of the problem (5).

Proof. To prove (a), refer to [14, 15].
Also, (b) and (c) can be easily proved considering that function 𝑄_𝜖(𝑥,𝜔) is the conjugate of function 𝑝(𝑧)=𝜖2‖𝑧‖²,𝑧∈𝑍,∞,𝑧∉𝑍, (7) where 𝑍=𝑧∈ℝ^𝑚2∶𝑊(𝜔)𝑧=𝑞(𝜔),𝑧≥0, (8) and Theorems (26-3) and (23-5) in [16].

Using the approximated recourse function 𝑄_𝜖(𝑥,𝜔), we can define a differentiable approximation function to the objective function of (1): 𝑓_𝜖(𝑥)=𝑐^𝑇𝑥+^𝑁_𝑖=1𝑄_𝜖𝑥,𝜔^𝑖𝜌𝜔^𝑖. (9) By (6), the gradient of above function exists and is obtained by ∇𝑓_𝜖(𝑥)=𝑐+^𝑁_𝑖=1∇𝑄_𝜖𝑥,𝜔^𝑖𝜌𝜔^𝑖=𝑐−^𝑁_𝑖=1𝑇^𝑇𝜔^𝑖𝑧^∗_𝜖𝑥,𝜔^𝑖𝜌𝜔^𝑖. (10) This approximation has paved the way to use the optimization algorithm for master problem (1) in which the objective function is substituted by 𝑓_𝜖(𝑥)min_𝑥∈𝑋𝑓_𝜖(𝑥). (11) In [7], it is considered slp problem with inequality constrained in master problem and equality constrained in recourse problem. Also, in Theorem 2.3 of [7], it is shown that a solution of the approximated problem is a good approximation to a solution of master problem. Here we can express a similar theorem for the problem (1) by using the similar technique in the proof of Theorem 2.3 in [7].

Theorem 2. Consider the problem (1). Then, for any 𝑥∈𝑋, there exists an 𝜖(𝑥)>0 such that for any 𝜖∈(0,𝜖(𝑥)]𝑓(𝑥)−𝑓_𝜖(𝑥≤𝜖)2𝑀, (12) where 𝑀 is defined as follows: 𝑀=max_{𝑖=1,2,…,𝑁}𝑧^∗_𝜖(𝑥,𝜔^𝑖)². (13) Let 𝑥^∗ be a solution of (1) and 𝑥^∗_𝜖 a solution of (11). Then, there exists an 𝜖>0 such that for any 0<𝜖≤𝜖max𝑓𝑥^∗_𝜖−𝑓𝑥^∗,𝑓_𝜖𝑥^∗−𝑓_𝜖𝑥^∗_𝜖≤𝜖2𝑀. (14) Further, one assumes that 𝑓 or 𝑓_𝜖 are strongly convex on 𝑋 with modulus 𝜇>0. Then, 𝑥^∗−𝑥^∗_𝜖𝜖≤𝑀𝜇. (15)

According to Theorem 1, it can be found that for obtaining the gradient of function 𝑓_𝜖(𝑥) in each iteration, we need the normal solution of 𝑁 linear programming problems (4). In this paper, the augmented Lagrangian method [17] is used for this purpose.

3. Smooth Approximation and Augmented Lagrangian Method

In the augmented Lagrangian method, the unconstrained maximization problem is solved which gives the project of a point on the solution set of the problem (4).

Assume that 𝑧 is an arbitrary vector. Consider the problem of finding the least 2-norm projection 𝑧_∗ of 𝑧 on the solution set 𝑍_∗ of the problem (4) 12𝑧_∗−𝑧²=min_{𝑧∈𝑍∗}12‖𝑧−𝑧‖²,𝑍_∗=𝑧∈ℝ^𝑚2∶𝑊(𝜔)𝑧=𝑞(𝜔),𝜉^𝑇𝑧=𝑄(𝑥,𝜔),𝑧≥0. (16) In this problem, vector 𝑥 and random variable 𝜔 are constants; therefore, for simplicity, this is assumed to be 𝜉=ℎ(𝜔)−𝑇(𝜔)𝑥, and function 𝑄(𝜉) is defined in a way that 𝑄(𝜉)=𝑄(𝑥,𝜔).

Considering that the objective function of the problem (16) is strictly convex, its solution is unique. Let us introduce the Lagrangian function for the problem (16) as follow: 𝐿=1𝑧,𝑝,𝛽,𝑧,𝜉,𝜔2‖𝑧−𝑧‖²+𝑝^𝑇𝑊(𝜔)(𝜔𝜉𝑧−𝑞)+𝛽^𝑇𝑧−𝑄(𝜉), (17) where 𝑝∈ℝ^𝑛2 and 𝛽∈ℝ are Lagrangian multipliers and 𝜉, 𝑧 are constant values. Therefore, the dual problem of (16) becomes max_𝛽∈ℝmax_{𝑝∈ℝ𝑛2}min_{𝑧∈ℝ𝑚2+}𝐿𝑧,𝑝,𝛽,𝑧,𝜉,𝜔. (18) By solving the inner minimization of the problem (18), duality of the problem (16) is obtained: max_𝛽∈ℝmax_{𝑝∈ℝ𝑛2}𝐿𝑝,𝛽,𝑧,𝜉, (19) where duality function is 𝐿𝑝,𝛽,𝑧,𝜉=𝑞^𝑇(𝜔)1𝑝−2𝑧+𝑊^𝑇(𝜔)𝑝+𝛽𝜉₊²𝑄(𝜉)+1+𝛽2‖𝑧‖². (20) The following theorem states that if 𝛽 is sufficiently large, solving the inner maximization of (19) gives the solution of the problem (16).

Theorem 3 (see [17]). Consider the following maximization problem max_{𝑝∈ℝ𝑛2}𝑆𝑝,𝛽,𝑧,𝜉,𝜔 (21) in which 𝛽, 𝑧, and 𝜉 are constants, and function 𝑆(𝑝,𝛽,𝑧,𝜉) is introduced as follows: 𝑆𝑝,𝛽,𝑧,𝜉,𝜔=𝑞^𝑇(𝜔)1𝑝−2𝑧+𝑊^𝑇(𝜔)𝑝+𝛽𝜉₊². (22) Also, assume that the set 𝑍_∗ is nonempty, and the rank of submatrix 𝑊_𝑙 of 𝑊 corresponding to nonzero components of 𝑧_∗ is 𝑛₂. In such a case, there is 𝛽^∗ which for all 𝛽≥𝛽^∗, 𝑧_∗=(𝑧+𝑊^𝑇𝑝(𝛽)+𝛽𝜉)₊ is the unique and exact solution for the problem (16), where 𝑝(𝛽) is the point obtained from solving the problem (21).

Also, in special conditions, the solution for the problem (3) can be also obtained and the following theorem expresses this issue.

Theorem 4 (see [17]). Assume that the solution set 𝑍_∗ is nonempty. For each 𝛽>0 and 𝑧∈𝑍_∗, 𝑦_∗=𝑝(𝛽)/𝛽 is one exact solution for the linear programming problem (3), where 𝑝(𝛽) is the solution for the problem (21).

According to the theorems mentioned above, augmented Lagrangian method presents the following iteration process for solving the problem (16): 𝑝_𝑘+1∈argmax_{𝑝∈ℝ𝑛2}𝑞^𝑇(𝜔)1𝑝−2𝑧_𝑘+𝑊(𝜔)^𝑇𝑝+𝛽𝜉₊²,𝑧_𝑘+1=𝑧_𝑘+𝑊^𝑇(𝜔)𝑝_𝑘+1+𝛽𝜉₊, (23) where 𝑧₀ is an arbitrary vector and here we can use zero vector as initial vector for obtaining normal solution of the problem (4).

We note that the problem (23) is a concave problem and its objective function is piecewise quadratic and is not twice differentiable. Applying the smoothing techniques [18, 19] and replacing 𝑥₊ by a smooth approximation, we transform this problem to a twice continuously differentiable problem.

Chen and Mangasarian [19] introduced a family of smoothing functions, which is built as follows. Let 𝜌∶𝑅→[0,∞) be a piecewise continuous density function satisfying ^+∞_−∞𝜌(𝑠)𝑑𝑠=1,^+∞_−∞|𝑠|𝜌(𝑠)𝑑𝑠<∞. (24) It is obvious that the derivative of plus function is step function, that is, (𝑥)₊=∫^𝑥_−∞𝛿(𝑡)𝑑𝑡, where the step function 𝛿(𝑥) is defined 1 if 𝑥>0 and equals 0 if 𝑥≤0. Therefore, a smoothing approximation function of the plus function is defined by 𝜑(𝑥,𝛼)=^𝑥_−∞𝜓(𝑡,𝛼)𝑑𝑡, (25) where 𝜓(𝑥,𝛼) is smoothing approximation function of step function and is defined as 𝜓(𝑥,𝛼)=^𝑥_−∞𝛼𝜌(𝛼𝑡)𝑑𝑡. (26) By choosing 𝜌(𝑠)=𝑒^−𝑠(1+𝑒^−𝑠)², (27) specific cases of these approaches are obtained as follows: 𝜓()=1𝑥,𝛼1+𝑒^−𝛼𝑥(𝑥),𝜑()1≈𝛿𝑥,𝛼=𝑥+𝛼log1+𝑒^−𝛼𝑥≈(𝑥)₊. (28) The function 𝜑 with a smoothing parameter 𝛼 is used here to replace the plus function of (22) to obtain a smooth reformulation of function (22): 𝑆𝑝,𝛽,𝑧,𝜉,𝛽,𝜔,𝛼∶=𝑞^𝑇(𝜔)𝑝−𝜑(𝑧+𝑊^𝑇(𝜔)𝑝+𝛽𝜉,𝛼)². (29) Therefore, we have the following iterative process instead of (23) and (28): 𝑝_𝑘+1∈argmax_{𝑝∈ℝ𝑛2}𝑞^𝑇(𝜔)𝑝−𝜑(𝑧_𝑘+𝑊^𝑇(𝜔)𝑝+𝛽𝜉,𝛼)²,𝑧_𝑘+1=𝜑𝑧_𝑘+𝑊^𝑇(𝜔)𝑝_𝑘+1+𝛽𝜉,𝛼. (30) It can be shown that as the smoothing parameter 𝛼 approaches infinity any solution of smooth problem (29) approaches the solution of the equivalent problem (22) (see [19]).

We begin with a simple lemma that bounds the square difference between the plus function 𝑥₊ and its smooth approximation 𝜑(𝑥,𝛼).

Lemma 5 (see [13]). For 𝑥∈ℝ and |𝑥|<𝜛𝜑²(𝑥,𝛼)−𝑥₊²≤log(2)𝛼²+2𝜛𝛼log(2), (31) where 𝜑(𝑥,𝛼) is the 𝜑 function of (28) with smoothing parameter 𝛼>0.

Theorem 6. Consider the problems (21) and max_{𝑝∈ℝ𝑛2}𝑆𝑝,𝛽,𝑧,𝜉,𝜔,𝛼. (32) Then, for any 𝑝∈ℝ^𝑛2 and 𝛼>0𝑆−𝑆≤(2)𝑝,𝛽,𝑧,𝜉,𝜔𝑝,𝛽,𝑧,𝜉,𝜔,𝛼log𝛼²+2𝑀𝑚₂log(2)𝛼, (33) where 𝑀 is defined as follows: 𝑀=max_{1≤𝑖≤𝑚2}𝑧_𝑖+𝑊^𝑇_𝑖(𝜔)𝑝+𝛽𝜉_𝑖. (34) Let 𝑝^∗ be a solution of (21) and 𝑝^∗_𝛼 a solution of (32). Then max𝑆𝑝^∗,𝛽,𝑧,𝜉,𝜔−𝑆𝑝^∗_𝛼,𝑆𝑝,𝛽,𝑧,𝜉,𝜔_∗𝛼,𝛽,𝑧,𝜉,𝜔,𝛼−𝑆𝑝^∗≤(2),𝛽,𝑧,𝜉,𝜔,𝛼log𝛼²+2𝑀𝑚₂log(2)𝛼. (35) Further, one assumes that 𝑊^𝑇 is a full rank matrix. Then, 𝑝^∗−𝑝^∗_𝛼(2)≤𝜇2log𝛼²+8𝑀𝑚₂𝜇log(2)𝛼. (36)

Proof. For any 𝛼>0 and 𝑝∈ℝ^𝑛2𝜑𝑧_𝑖+𝑊^𝑇_𝑖(𝜔)𝑝+𝛽𝜉_𝑖,𝛼≥𝑧_𝑖+𝑊^𝑇_𝑖(𝜔)𝑝+𝛽𝜉_𝑖,𝛼₊. (37) Hence 𝑆−𝑆−𝑆=𝑝,𝛽,𝑧,𝜉,𝜔𝑝,𝛽,𝑧,𝜉,𝜔,𝛼=𝑆𝑝,𝛽,𝑧,𝜉,𝜔𝑝,𝛽,𝑧,𝜉,𝜔,𝛼𝜑(𝑧+𝑊^𝑇(𝜔)𝑝+𝛽𝜉,𝛼)²−𝑧+𝑊^𝑇(𝜔)𝑝+𝛽𝜉,𝛼₊²=^𝑚2_𝑖=1𝜑²𝑧_𝑖+𝑊^𝑇_𝑖(𝜔)𝑝+𝛽𝜉_𝑖−𝑧,𝛼_𝑖+𝑊^𝑇_𝑖(𝜔)𝑝+𝛽𝜉_𝑖,𝛼²₊. (38) By using Lemma 5, we get that 𝑆−𝑆≤𝑝,𝛽,𝑧,𝜉,𝜔𝑝,𝛽,𝑧,𝜉,𝜔,𝛼_𝑚2𝑖=1log(2)𝛼²+2𝑧_𝑖+𝑊^𝑇_𝑖(𝜔)𝑝+𝛽𝜉_𝑖(2),𝛼log𝛼≤log(2)𝛼²(2)+2log𝛼×^𝑚2_𝑖=1max_{1≤𝑖≤𝑚2}𝑧_𝑖+𝑊^𝑇_𝑖(𝜔)𝑝+𝛽𝜉_𝑖=log(2)𝛼²+2𝑀𝑚₂log(2)𝛼. (39) From above inequality, we have 𝑆≤𝑆+(2)𝑝,𝛽,𝑧,𝜉,𝜔,𝛼≤𝑆𝑝,𝛽,𝑧,𝜉,𝜔𝑝,𝛽,𝑧,𝜉,𝜔,𝛼log𝛼²+2𝑀𝑚₂log(2)𝛼. (40) Therefore 𝑆𝑝^∗,𝛽,𝑧,𝜉,𝜔−𝑆𝑝^∗_𝛼≤𝑆𝑝,𝛽,𝑧,𝜉,𝜔_∗,𝛽,𝑧,𝜉,𝜔,𝛼−𝑆𝑝^∗_𝛼+(2),𝛽,𝑧,𝜉,𝜔,𝛼log𝛼²+2𝑀𝑚₂log(2)𝛼≤log(2)𝛼²+2𝑀𝑚₂log(2)𝛼,𝑆𝑝^∗_𝛼,𝛽,𝑧,𝜉,𝜔,𝛼−𝑆𝑝^∗𝑝,𝛽,𝑧,𝜉,𝜔,𝛼≤𝑆^∗_𝛼,𝛽,𝑧,𝜉,𝜔−𝑆𝑝^∗+(2),𝛽,𝑧,𝜉,𝜔log𝛼²+2𝑀𝑚₂log(2)𝛼≤log(2)𝛼²+2𝑀𝑚₂log(2)𝛼. (41) Suppose that 𝑊^𝑇 is full rank. Then the Hessian of 𝑆 is negative definite, and 𝑆 is strongly concave on bounded sets. By the definition of strong concavity, for any 𝛾∈(0,1), 𝑆𝛾𝑝^∗_𝛼+1−𝛾𝑝^∗,𝛽,𝑧,𝜉,𝜔,𝛼−𝛾𝑆𝑝^∗_𝛼−𝑆𝑝,𝛽,𝑧,𝜉,𝜔,𝛼1−𝛾^∗≥1,𝛽,𝑧,𝜉,𝜔,𝛼2𝜇𝛾1−𝛾𝑝^∗_𝛼−𝑝^∗2. (42) Let 𝛾=1/2, then 18𝜇𝑝^∗_𝛼−𝑝^∗2≤𝑆12𝑝^∗_𝛼+𝑝^∗−1,𝛽,𝑧,𝜉,𝜔,𝛼2𝑆𝑝^∗_𝛼,𝛽,𝑧,𝜉,𝜔,𝛼+𝑆𝑝^∗≤𝑆𝑝,𝛽,𝑧,𝜉,𝜔,𝛼_∗𝛼−1,𝛽,𝑧,𝜉,𝜔,𝛼2𝑆𝑝^∗_𝛼,𝛽,𝑧,𝜉,𝜔,𝛼+𝑆𝑝^∗≤1,𝛽,𝑧,𝜉,𝜔,𝛼2𝑆𝑝^∗_𝛼,𝛽,𝑧,𝜉,𝜔,𝛼−𝑆𝑝^∗≤1,𝛽,𝑧,𝜉,𝜔,𝛼2log(2)𝛼²+𝑀𝑚₂log(2)𝛼. (43)

Considering the advantage of the twice differentiability of the objective function of the problem (32) allows us to use a quadratically convergent Newton algorithm with an Armijo stepsize [20] that makes the algorithm globally convergent.

4. Numerical Results and Algorithm

In each iteration of the process (30), one concave, quadratic, unconstrained maximization problem is solved. For solving it, the fast Newton method can be used.

In the algorithm, the Hessian matrix may be singular, thus we use a modified Newton. The direction in each iteration for solving (30) is obtained through the following relation: 𝑑_𝑠=−∇²_𝑝𝑆𝑝,𝛽,𝑧,𝜉,𝜔,𝛼−𝛿𝐼_𝑛2⁻¹∇_𝑝𝑆,𝑝𝑝,𝛽,𝑧,𝜉,𝜔,𝛼_𝑠+1=𝑝_𝑠+𝜆_𝑠𝑑_𝑠, (44) where 𝛿 is a small positive number, 𝐼_𝑛2 is the identity matrix of order 𝑛₂, and 𝜆_𝑠 is the suitable step length that Armijo algorithm is used for determining it (see Algorithm 1).

Choose a 𝑧0∈𝑅𝑚2, 𝛼>1,𝜄>1, 𝜖>0 be error tolerance and 𝛿 is a small positive number. 𝑖=0; While  𝑧𝑖−𝑧𝑖−1∞≥𝜖 Choose a 𝑝0∈𝑅𝑛2 and set 𝑘=0. While  ∇𝑝𝑆(𝑝𝑘,𝛽,𝑧𝑖,𝜉,𝜔,𝛼)∞≥𝜖 Choose  𝜆𝑘=  max{𝑠,𝑠𝜎,𝑠𝜎2,…} such that 𝑆(𝑝𝑘,𝛽,𝑧𝑖,𝜉,𝜔,𝛼)−𝑆(𝑝𝑘+𝜆𝑘𝑑𝑘,𝛽,𝑧𝑖,𝜉,𝜔,𝛼)≥−𝜆𝑖𝜇∇(𝑆(𝑝𝑘,𝛽,𝑧𝑖,𝜉,𝜔,𝛼))𝑇𝑑𝑘,  where, 𝑑𝑘=−(∇2𝑆(𝑝𝑘,𝛽,𝑧𝑖,𝜉,𝜔,𝛼,)−𝛿𝐼𝑛2)−1∇𝑆(𝑝𝑘,𝛽,𝑧𝑖,𝜉,𝜔,𝛼),𝑠>0  be a constant, 𝜎∈(0,1)  and 𝜇∈(0,1). Put 𝑝𝑘+1=𝑝𝑘+𝜆𝑘𝑑𝑘, 𝑘=𝑘+1  and 𝛼=𝜄𝛼. end Set 𝑝𝑖+1=𝑝𝑘+1, 𝑧𝑖+1=𝜑(𝑧𝑖+𝑊𝑇(𝜔)𝑝𝑖+1+𝛽𝜉,𝛼)  and 𝑖=𝑖+1. end

Algorithm 1 

Newton method with the Armijo rule.

The proposed algorithm was applied to solve some recourse problems. Table 1 compares this algorithm with CPLEX v. 12.1 solver for quadratic convex programming problems (5). As is evident from Table 1, most of recourse problems could be solved more successful by the algorithm which is based on smooth augmented Lagrangian Newton method (SALN) than CPLEX package (for illustration see the problems 21–25 in Table 1). This algorithm gives us high accuracy and the solution with minimum norm in suitable time (see last column of Table 1). Also, we can find that CPLEX is better than the algorithm proposed for some recourse problems in which the matrices are approximately square (Ex. line 5–12).

Table 1 

Comparative between smooth augmented Lagrangian Newton method (SALN) and CPLEX solver.

The test generator generates recourse problems. These problems are generated using the MATLAB code show in Algorithm 2.

%Sgen: Generate random solvable recourse problems: %Input: m,n,d(ensity); Output: W,q,𝜉; m=input(’Enter 𝑛2:’) n=input(’Enter 𝑚2:’) d=input(’Enter d:’) pl=inline(’(abs(x)+x)/2’) W=sprand(𝑛2,𝑚2,d);W=100*(W-0.5*spones(W)); z=sparse(10*pl(rand(𝑚2,1))); q=W*z; y=spdiags((sign(pl(rand(𝑛2,1)-rand(𝑛2,1)))),0,𝑛2,𝑛2) *5*((rand(𝑛2,1)-rand(𝑛2,1))); ξ=W’*y-10*spdiags((ones(𝑚2,1)-sign(z)),0,𝑚2,𝑚2)*ones(𝑚2,1)); format short e; nnz(W)/prod(size(W))

Algorithm 2 

The algorithm considered for solving several recourse problems was run on a computer with 2.5 dual-core CPU and 4 GB memory in MATLAB 7.8 programming environment. Also, in the generated problems, recourse matrix 𝑊 is the Sparse matrix (𝑛₂×𝑚₂) with the density 𝑑. The constants 𝛽 and 𝛿 in the above algorithm in (44) were selected 1 and 10⁻⁸, respectively.

In Table 1, the second column indicates the size and density of matrix 𝑊, the forth column indicates the feasibility of the primal problem (4), and the next column indicates the error norm function of this problem (the MATLAB code of this paper is available from the authors upon request).

5. Conclusion

In this paper, a smooth reformulation process, based on augmented Lagrangian algorithm, was proposed for obtaining the normal solution of recourse problem of a stochastic linear programming. This smooth iterative process allows us to use a quadratically convergent Newton algorithm, which accelerates obtaining the normal solution.

Table 1 shows that the proposed algorithm has appropriate speed in most of the problems. This result, specifically, can be observed in recourse problems with the matrix of coefficients in which the number of constraints is noticeably more than the number of variables. The more challenging is solving the problems which their coefficient matrix is square (the numbers of constraints and variables get closer to each other), and more time is needed by the algorithm for solving the problem.

Acknowledgment

The authors would like to thank the reviewers for their helpful comments.