# A Smoothing Algorithm for a New Two-Stage Stochastic Model of Supply Chain Based on Sample Average Approximation

Authors: Liu Yang; Yao Xiong; Xiao-jiao Tong
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

We construct a new two-stage stochastic model of supply chain with multiple factories and distributors for perishable product. By introducing a second-order stochastic dominance (SSD) constraint, we can describe the preference consistency of the risk taker while minimizing the expected cost of company. To solve this problem, we convert it into a one-stage stochastic model equivalently; then we use sample average approximation (SAA) method to approximate the expected values of the underlying random functions. A smoothing approach is proposed with which we can get the global solution and avoid introducing new variables and constraints. Meanwhile, we investigate the convergence of an optimal value from solving the transformed model and show that, with probability approaching one at exponential rate, the optimal value converges to its counterpart as the sample size increases. Numerical results show the effectiveness of the proposed algorithm and analysis.


---

## Body

## 1. Introduction

A supply chain planning (SCP) is a network of suppliers, manufacturing plants, warehouses, and distribution channels organized to acquire raw materials, convert these raw materials into specified final products, and distribute these products to customers. In fact, the SCP involves the commitment of substantial capital resources over long periods of time and it makes the SCP problem extremely important in practice.Beginning with the seminal work of Geoffrion and Graves [1] on multicommodity distribution system design, a large number of optimization-based approaches have been proposed for the problem of SCP; see, for example, [2–4]. However, most real SCP problems are characterized by numerous sources of technical and commercial uncertainty; critical parameters such as customer demands, suppliers, and resource capacities are quite uncertain. There are a few research works addressing comprehensive SCP with two-stage stochastic models. As far as we know, Tsiakis et al. considered a two-stage stochastic programming model for SCP under demand uncertainty, where the authors developed a large-scale mixed-integer linear programming model; see [5]. Santoso et al. integrated a sampling strategy with an accelerated Benders decomposition to solve SCP problems with continuous distributions for the uncertain parameters; see [6]. Azaron et al. developed a multiobjective stochastic programming approach for SCP under uncertainty, which also used the goal attainment technique to obtain the Pareto-optimal solutions for decision-making; see [7]. However, very few people considered applying the thought of stochastic dominant to SCP. Recently, optimization models with stochastic dominance constraints have increasingly been the subject of theoretical considerations and have been applied in many areas such as finance, energy, and transportation; see [8–11]. Stochastic dominance is one of the popular established ways to formalize preferences among random values. We can see the concept of the first and the second-order stochastic dominance in detail in [12]. Dentcheva and Ruszczyński have done much work in the applications of stochastic dominance. They firstly provided a stochastic optimal problem with SSD constraint and then discussed the optimality, properties, algorithms, and the applications; see [13–16]. This is an attractive approach for managing risks in an optimization setting. While pursuing expected profits, one can avoid high risks by choosing options that are preferable to a random benchmark. Dentcheva and Martinez also considered a two-stage stochastic optimization problem with a random sequence constraint based on recourse; they proposed two decomposition algorithms and applied them in a two-stage SCP problem; see [17]. Sun discussed the properties, algorithms, and the applications of the stochastic optimization model with SSD constraints. He proposed an improved cutting plane algorithm and applied it to solve a portfolio problem and a SCP problem; see [18].In this paper, we introduce aSSD constraint to a two-stage SCP optimal model in which we design an acceptable benchmark according to the risk preferences consistency. Our model is a new model to use the SSD constraint in the two-stage SCP model which is a special model of [17, 18]. We do not only minimize the expected cost but also prevent possibly very large cost, which may occur with relatively small probability. For solving such a two-stage SCP optimal model, firstly we use the sample average values to approximate the expected values of the underlying random functions and reformulate them as an ordinary nonlinear programming model with finite number of constraints. Then we provide a smoothing penalty function algorithm to deal with the nonsmoothness resulting from the plus function and penalty function. There are four advantages in our models and methods: (i) we introduce a SSD constraint to describe the preference consistency of the risk taker while minimizing the expected cost of company; (ii) we transfer the two-stage problem into an equivalent one-stage problem which can be solved more easily; (iii) the smoothing SAA method can avoid the infinite constraints in the transformed models and the size of the smoothing algorithm model will not increase as the sample grows; (iv) the smoothing algorithm can get the global optimal solution because the smoothing function maintains the original convexity. Meanwhile, we investigate the convergence of the optimal value of the transformed model and show that the optimal value converges to its counterpart with probability approaching one at exponential rate as the sample size increases.The remainder of this paper is organized as follows: in Section2, we set a two-stage stochastic model for SCP by introducing SSD constraint and reformulate it into a one-stage model equivalently. In Section 3, we use a SAA scheme to tackle the expected values. In Section 4, we provide a smoothing penalty function algorithm to solve the approximated SAA model and discuss the convergence of the optimal value. Finally, we give numerical results of a specific small scale SCP to illustrate the proposed methodology in Section 5.


## 2. A Two-Stage SCP Optimal Model with the Stochastic Dominant Constraint

Consider aSCP has a set F of factories that produce and supply perishable products to a set O of distributors. Obviously, the demand of distributors is a random vector and we assume that the goods are supplied before the demand is observed. Let xij(i=1,…,F;j=1,…,O) denote the quantity of products delivered by factory i to distributor j; the demand quantity of distributor j is Dj,j=1,…,O. Then, if ∑i=1Fxij<Dj, the customers will buy the product in another place and the sales are lost; if ∑i=1Fxij>Dj, the remaining products have to be disposed owing to its corruption. Here we assume that each factory has a limited capacity to produce products and the disposal cost is a deterministic quantity. The quantity produced by factory i is wi; the production cost of one unit of product at factory i is ai; zj is the number of sales at distributor j; yj is the amount disposed at distributor j; dj is the disposal cost on site j; bj is the price distributor j sets for the product; the shipping cost from factory i to distributor j is denoted by cij; the capacity of distributor j is Cj.We adopt the point of view that a random benchmark is available. For the random demand of the distributors, we design an acceptable benchmark strategy depending on preference consistency of risk taker while minimizing the expected cost of the company. Our aim is develop a production and shipping plan to minimize the expected costs and maximize the expected return. The mathematical model for thisSCP problem is formulated as follows: (1)minEPFx,D≔EP∑i=1F∑j=1Ocijxij+∑i=1Faiwi+Qx,Ds.t.-Qx,D⪰2-Z,wi=∑j=1Oxij,xij≥0,i=1,…,F;j=1,…,O,where D≔{Dj,j=1,…,O}. The first and second item in the objective function represent the total cost of shipping and production, respectively. The SSD constraint -Q(x,D)⪰2-Z shows that the cost should dominate the benchmark strategy Z at any time. SSD is a popular order when large outcomes are preferred due to its consistency with risk-averse preferences. The equality constraint denotes the limitation of production and the functional Q(x,D) stands for the optimal value of the following second-stage problem: (2)min∑j=1Odjyj-bjzjs.t.∑i=1Fxij=zj+yj,∑i=1Fxij≤Cj,0≤zj≤Dj,yj≥0,j=1,…,O.Here the objective function denotes the total disposal cost minus the total return of the distributors. The equality constraint comprises the relationship between zj and yj. And the inequality constraint represents the limitation of capacity and sales, respectively.In order to facilitate solving theSCP, we consider transforming it into a one-stage stochastic model firstly. Note that, by using the plus function, we can rewrite second-stage variables—the amount disposed and sold as yj=xij-Dj+ and zj=Dj-Dj-∑i=1Fxij+ according to the relation of the original problem. So with the special construction of this SCP problem, we can obtain the following equivalence one-stage stochastic model easily: (3)minfx≔EPFx,Ds.t.-Qx,D≔-∑j=1Odj∑i=1Fxij-Dj++bjDj-Dj-∑i=1Fxij+⪰2-Z,wi=∑j=1Oxij,∑i=1Fxij≤Cj,xij≥0,i=1,…,F;j=1,…,O,where ·+=max⁡(·,0). Obviously, problem (3) is a convex stochastic optimal model. The solution set X is a closed convex subset of RF×O, D:Ω→Ξ is a vector of random variables defined on probability (Ω,F,P) with support set Ξ⊂RO, F,Q:RF×O→R are continuous functions, and EP[·] denotes the expected value. Mathematically, the second-order dominance constraint in (3) can be reformulated as [13] (4)EPQx,D-η+≤EPZ-η+,∀η∈R,It is well known that above constraint does not satisfySlater’s constraint qualification, which could not obtain the necessary conditions of the optimal solution. In general, a so-called relaxed form is proposed: (5)EPQx,D-η+≤EPZ-η+,∀η∈a,b.


## 3. Sample Average Approximation

Because problem (3) involves the random vector D, it is very difficult to obtain an approximate form of the expected value of the underlying random functions in general. But if the random vector satisfies a finite discrete distribution, we can easily convert (3) into a deterministic nonlinear programming problem. Recently, Sun et al. applied the well-known sample average approximation (SAA) method to solve a stochastic optimization problem with second-order stochastic dominance constraints successfully; see [19]. So we can also construct a SAA problem of (3) in this paper.LetDj1,…,DjS,j=1,…,O be an independent and identically distributed sampling of each component Dj. By Proposition 3.2 in [14], we know that (5) is equivalent to (6)gx,Zs≔EPQx,D-Zs+-EPZ-Zs+≤0,Zs∈a,b,where Zs≔Z(Djs),s=1,2,…,S, are benchmark strategies with sample of size S. Thus, we can get the SAA problem of (3) as follows: (7)minfSx≔1S∑s=1SFx,DSs.t.gSx,Zs≔1S∑s=1SQx,Ds-Zs+-1S∑s=1SZ-Zs+≤0,∀Zs∈a,b,wi=∑j=1Oxij,∑i=1Fxij≤Cj,xij≥0,i=1,…,F;j=1,…,O,where (8)Qx,Ds≔∑j=1Odj∑i=1Fxij-Djs++bjDjs-Djs-∑i=1Fxij+.


## 4. A Smoothing Penalty Function Algorithm

For a fixed sample, theSAA problem (7) is an ordinary nonlinear programming problem with finite number of constraints. However, the number of constraints depends on the sample size S heavily [12], and this may lead to difficulty in solving this problem as the number of the samples becomes large. So we consider an exact penalization function method to handle (7); that is, consider the following penalized model: (9)minφSx,ρS≔EPSFx,D+ρSS∑s=1Smax⁡gSx,Zs,0s.t.wi=∑j=1Oxij,∑i=1Fxij≤Cj,xij≥0,i=1,…,F;j=1,…,O,where gSx,Zs=EPSQx,Ds-Zs+-Z-Zs+andρS>0 is a constant depending on S.The objective function andgS(x,Zs) in the penalized problem (9) are nonsmooth because they are containing the plus function. To overcome the difficulty of the nonsmoothness, we can provide a smoothing method which can be differentiable and assure the global convergence of the true model; see [19–22]. In this paper, we propose the smoothing function proposed in [21] to approximate the maximal function (0,·)+. Let t∈R+ be a smoothing parameter; for every t>0, we can get the following smoothing function:(10)p¯t,x≔tln1+exp⁡xt,t>0,0,x+,t=0.It is easy to verify that (11)limt→0+p¯t,x=0,x+,which implies that p¯(t,x) is continuous in t for every x; see [21]. Note that the continuity of p¯(t,x) in t on (0,∞) is obvious. Moreover, p¯(t,x) is uniformly globally Lipschitz continuous with respect to (t,w) with modulus 4ln⁡2; meanwhile, p¯(·,x) and p¯(t,·) are uniformly globally Lipschitz continuous with modulus 1 and modulus 2ln⁡2, respectively; see [21]. So we can construct the smoothing function for the maximal function in (9) as follows:(12)Q¯t,x,Ds=∑j=1Odjtln⁡1+exp⁡∑i=1Fxij-Djst-∑j=1ObjDjs-tln⁡1+exp⁡Djs-∑i=1Fxijt,g¯St,x,Zs=1S∑s=1Stln⁡1+exp⁡Q¯t,x,Ds-Zst-1S∑s=1SZ-Zs+,where Q¯(t,x,Ds) and g¯S(t,x,Zs) are the smoothing function of Q(x,Ds) and gS(x,Zs), respectively. Then we have the following conclusion for any t>0.Theorem 1.
(i)g¯S(t,x,Zs) is increasing with respect to t. Furthermore, we have(13)g¯St,x,Zs-gSx,Zs≤tln⁡t.
(ii)gS(x,Zs) is a convex function with respect to x and g¯S(t,x,Zs) remains the convexity of gS(x,Zs) with respect to x.
(iii)g¯S(t,x,Zs) is continuously differentiable with (14)∇xg¯St,x,Zs=1S∑s=1Se∇xQ¯t,x,Ds-Zs/t1+e∇xQ¯t,x,Ds-Zs/t·∇xQ¯t,x,Ds,where (15)∇xQ¯t,x,Ds=∑j=1Odje∑i=1Fxij-Djs/t1+e∑i=1Fxij-Djs/t+∑j=1ObjeDjs-∑i=1Fxij/t1+eDjs-∑i=1Fxij/t.
(iv) For any fixedx∈RF×O, (16)∇xg¯St,x,Zs-∂gSx,Zs=ot.Proof.
We can obtain conclusions (i), (iii), and (iv) directly from Lemma3.1 in Tong et al. [22]. For conclusion (ii), because g¯S(t,x,Zs) is the sum function of Q¯(t,x,Ds), which is a linear function, Q¯(t,x,Ds) is a convex function, and the sum function of convex functions is still convex. This completes the proof of (ii).From Theorem1, we can obtain the following smoothing penalized problem: (17)minφ¯St,x,ρS≔EPSFx,D+ρsS∑s=1Stln⁡1+exp⁡g¯St,x,Zsts.t.wi=∑j=1Oxij,∑i=1Fxij≤Cj,xij≥0,i=1,…,F;j=1,…,O.In the following, we investigate the convergence of optimal solution of problem (17) as sample size increases.Assumption 2.
There exist a vectoru∈Rn and a constant δ>0 such that (18)g0x,Zs;u≤-δ,∀Zs∈a,b,x∈X,where g0(x,Zs;u) denotes Clarke generalized directional derivative of g(x,Zs) at point x in direction u for given zs∈[a,b].Assumption 3.
There exists a positive integrable functionκ(D) such that (19)max⁡∇xFx,D,∇xQx,D≤κD.Assumption 4.
Problem (3) satisfies the uniform dominance condition (Slater’s constraint qualification); that is, there exists a point x0∈X such that (20)supZs∈a,bgx0,Zs<0.Theorem 5.
(a) Let Assumptions2–4 hold; (b) assume that there exist x0∈X and Zs0∈[a,b] such that EP[F(x0,D)]≤∞, EP[Q(x0,D)-Zs0+]≤∞, and EP[Z-Zs0+]≤∞. X is a compact set; (c) the moment generating function EP[eκ(D)t] of the random variable κ(D) is finite valued as t is close to 0, where κ(D) is defined by Assumption 3; (d) for every x∈X, the moment generating functions EP[e(F(x,D)-EP[F(x,D)])/t],EP[e((Zs-Q(x,D))+-EP[(Zs-Q(x,D))+])/t], and EP[e((Zs-Z)+-EP[(Zs-Z)+])/t] are finite valued as t is close to 0. (e) F(x,D) is convex and Q(x,D) is concave for almost every D∈Ξ. Let ϑ(PS) and ϑ¯(PS,tS), respectively, denote the optimal value function of the SAA problem (7) and smoothing penalized problem (17). Then for any small positive number α, there exist positive constants S∗, C(α), and β(α) such that (21)Probϑ¯PS,ts-ϑP≤α≥1-Cαe-βαS,for S≥S∗.Proof.
Note that(22)ϑ¯PS,tS-ϑP≤ϑ¯PS,tS-ϑPS+ϑPS-ϑP.By Theorem 3.1 [19], there exist positive numbers S∗ and L∗ such that (23)ϑ¯PS-ϑP≤L∗DPS,P,ProbϑPS-ϑP≤α2≥1-Cαe-βαS,for N≥N∗, where DPS,P≔suph∈HEPSh-EPh,H≔gx,Zs:x∈X∪{f(x):x∈X}∪κ·, which is the distance function. By Theorem 4.1 [19], we obtain (24)ϑ¯PS,tS-ϑPS≤supx∈Xφ¯Sx,ρS,tS-φSx,ρS≤supx∈XρSS∑s=1Stln⁡1+exp⁡g¯St,x,Zst-max⁡gSx,Zs,0≤supx∈XρSS∑s=1Sg¯St,x,Zs-gSx,Zs+2ln⁡2≤4ln⁡2tρS;let (4ln⁡2)tρS≤α/2. Then (25)ϑ¯PS,tS-ϑPS≤α2.Combining (22), (23), (24), and (25), we obtain the following formulation for S≥S∗: (26)Probϑ¯PS,ts-ϑP≤α≥Probϑ¯PS,tS-ϑP≤α2≥1-Cαe-βαS.


## 5. Numerical Results

In this section, we provide a specificSCP for F=2 and O=3. We carry out numerical tests by using solver fmincon built in MATLAB R2012a, which is installed on aLenovo PC with Windows XP operating system and 2 GB of RAM. Assume that each component Dj is set to be gamma distribution with parameters (2,3). The related data is as follows: cij=0.50.80.7;0.60.80.9;ai=10.5;wi=1011;dj=0.50.70.8;bj=4.24.85. Let the sample size S=500, smoothing parameter t=0.001, and penalty parameter ρ=1000; we use the model and algorithm in Sections 2–4 to solve the SCP problem. We obtain that the optimal values are xij=460;506 and E[·]=-254.4191. From the optimal model outcomes, our testing results indeed satisfy the constraints of the SCP model.We consider the effect of different smoothing parameters. Let the sample sizeS=500 and smoothing parameter t=0.001; see Figure 1. It can be observed that the optimal target values do not change, which indicates that the smoothing algorithm is insensitive to the parameter t. It is further shown that the original model based on SSD is equivalent to the smoothing approximation penalty function model when t tends to 0+. In what follows, we discuss the CPU times of new algorithm as the increase of sample size. From Figure 2, an important observation is that the methods scale well with the dimension: increasing the dimension of the true problem (3) does not result in a big increase of the number of CPU times, which may be because an exact penalty function method is subsequently applied to deal with the difficulty that the number of dimensions and constraints heavily depends on sample size.Figure 1
The expected returns for different smoothing parameters.Figure 2
The CPU time as sample size increases.In addition, we can make a comparison of our algorithm and proposed methods in [17]. For the same case of F=3,O=3 and S=100. Running time for our algorithm, which requires only the solution of an ordinary nonlinear convex programming problem with dimension of 2n, is always within one second (0.429 seconds at most). However, the decomposition methods from [17] typically require some minutes (at least 198.0 seconds). Hence, from a computational point of view, our algorithm is clearly more preferable.


## 6. Conclusion

This paper proposes a smoothing algorithm for a new two-stage stochastic optimization model of the supply chain planning (SCP). Our model and method have four advantages: (i) the second-order stochastic dominance constraint can describe the preference consistency of the risk taker while minimizing the expected cost of company; (ii) an equivalent one-stage problem was given to solve the two-stage problem; (iii) a SAA method was used to calculate the expected values and it can avoid the infinite constraints in the original problem; (iv) the smoothing algorithm can get the global solution and avoid introducing new variables and constraints. Meanwhile, we investigate the convergence of the optimal value of the transformed model and show that the optimal value converges to its counterpart with probability approaching one at exponential rate as the sample size increases. Numerical results show the effectiveness of the proposed model and algorithm.



---

*Source: 5681502-2017-10-19.xml*