# A Novel Algorithm of Stochastic Chance-Constrained Linear Programming and Its Application

Authors: Xiaodong Ding; Chengliang Wang
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

The computation problem is discussed for the stochastic chance-constrained linear programming, and a novel direct algorithm, that is, simplex algorithm based on stochastic simulation, is proposed. The considered programming problem in this paper is linear programming with chance constraints and random coefficients, and
therefore the stochastic simulation is an important implement of the proposed algorithm. By theoretical analysis, the theory basis of the proposed algorithm is obtained and, by numerical examples, the feasibility and validness of this algorithm are illustrated. The detailed algorithm procedure is given, which is easily converted into the executable codes of software tools. Then, we compare it with some algorithms to verify its superiority. Finally, a practical example is presented to show its practicability.


---

## Body

## 1. Introduction

In the late 1950s, stochastic linear programming (SLP) appeared with the further application of linear programming. SLP is a special kind of linear programming problem in which a part or all of coefficients are random variables with joint probability distribution. Generally speaking, there are two sorts of SLP models, one of which can be described as “wait-and-see” model based on the hypothesis that the decision maker can wait until the random variables come true and another one is called “here-and-now” model in which the decision maker must make decision before the random variables come true.Stochastic chance-constrained programming (SCP), firstly proposed by Charnes and Cooper [1], offers a powerful means of modeling stochastic decision and control systems (see, e.g., [2, 3]). SCP is mainly concerned with the problem that the decision maker must give his solution before the random variables come true. In this problem, the made decision may not satisfy the constraints in some degree, but the probability of decision satisfying the constraints cannot be less than some given confidence level α. Stochastic chance-constrained linear programming (SCLP) is an important part of SCP.As we know the traditional method to solve SCLP is converting it into an equivalent deterministic linear programming and then obtaining the optimal solution by some deterministic algorithms [4]. However, this method is only effective for some special cases. Generally, SCLP cannot be converted into some deterministic linear programming and convex programming, see the work of Kall and Wallace in [5]. However, for those which can be converted to deterministic problem, it is usually a difficult work to convert and one can always obtain a complicated nonlinear programming which is traditional hard problem. Therefore, it is necessary and urgent to find the direct and effective algorithm to solve SCLP. Luckily, with the rapid development of computer, genetic algorithms based on stochastic simulation were designed for SCP, see, for example, [6–8]. These intelligent algorithms are more direct and effective than the method of converting SCP to deterministic programming. However, the disadvantages of these genetic algorithms are obvious. That is, they are always designed aiming to solve every specified problem, depending on experiment and lacking the common theory basis. Therefore, it is a necessary and urgent task to find a better algorithm, which gives rise to the motivation for the present study in this paper.Summarizing the above discussion, we aim to develop a novel direct and universal algorithm to solve the computation problem of the stochastic chance-constrained linear programming. The main contribution of this paper can be given as follows: (i) several simple approaches are illustrated to be ineffective or limited to SCLP by numerical example; (ii) a novel and direct algorithm, that is, simplex algorithm based on stochastic simulation, is proposed; (iii) the theory basis of the proposed algorithm is proved; (iv) the detailed procedures of the simplex algorithm are given. By numerical examples, our algorithm is compared with the traditional methods and genetic algorithm based on stochastic simulation and verified to be feasible, valid, and better than the traditional ones and intelligent technique. Finally, a real-world example is given to illustrate the practicability of the developed algorithm.


## 2. Research Model and Computation Problem of SCLP

In this paper, unless otherwise specified,(Ω,ℱ,P) denotes a complete probability space, where Ω is a nonempty sample space, ℱ is the power set of Ω, and P is a probability measure on ℱ, and a random variable is defined as a function from the probability space (Ω,ℱ,P) to the set of ℱ.


### 2.1. Stochastic Linear Programming (SLP) Model

Firstly, introduce the following SLP model(2.1)min⁡CTXs.t.AX≤bX≥0,
where X is an n dimensional vector to be determined, and A=(aij(ω))m×n,b=(bi(ω))m×1,C=(cj(ω))n×1, with ai,j(ω),bi(ω),cj(ω)(i=1,…,m,j=1,…,n) being random variables and ω∈Ω.


### 2.2. Stochastic Chance-Constrained Linear Programming (SCLP) Model

SCLP usually includes two sorts of models which can be formulated as follows:(2.2)min⁡CTXs.t.P(AX≤b)≥αX≥0,(2.3)min⁡CTXs.t.P(AiX≤bi)≥αi,i=1,2,…,mX≥0,
where Ai is the ith row of A, bi is the ith element of b, and αi is the ith confidence level of the constraints.


### 2.3. Computation Problem of SCLP

In this subsection, by solving a numerical example of SLP, the computation problem of SCLP is analyzed and several possible approaches are tried to obtain the optimal solution to the example.Example 2.1.
Let(a,b) be the random variables with uniform distribution in rectangle {1≤a≤4,1/3≤b≤1} and consider the following SLP problem:
(2.4)min⁡x1+x2s.t.ax1+x2≥7bx1+x2≥4x1≥0,x2≥0.Solution 1.
In (2.4), the expectations of the random variables a and b are easily derived. Therefore, a very simple idea is to replace the stochastic parameters a and b by their expectations, respectively, and solve the corresponding deterministic linear programming problem described as follows:
(2.5)min⁡x1+x2s.t.52x1+x2≥723x1+x2≥4x1≥0,x2≥0.
From (2.5), by deterministic linear programming algorithms, it is easy to derive the unique optimal solution x1*=18/11,x2*=32/11.
Now, in order to analyze the feasibility of this approach, we assumeθ*=(18/11,32/11), and, from the example, we can easily know the feasible region D={(x1,x2)∣ax1+x2≥7,bx1+x2≥4}. Then
(2.6)P(θ*∈D)=P(1811a+3211≥7,1811b+3211≥4)=P(a≥52,b≥23)=0.25.
That is, the probability of the solution θ* taking value in the feasible region D is only 0.25, which verifies that the method of directly replacing random parameters by their expectations is invaluable.Solution 2.
It is well known that the samples of random variables are also easily obtained; therefore, another simple technique is to produce some random samples of these random parameters and solve all the deterministic linear programming problems corresponding to the samples, and then we choose the best solution as the optimal solution. Now, produce 10 samples ofa and b in (2.4):
(2.7)a:(1.99,2.02,2.47,3.23,2.56,3.78,1.68,3.95,2.33,3.65),b:(0.63,0.62,0.41,0.58,0.85,0.59,0.47,0.80,0.57,0.62).
From Table 1, it is easy to see that the largest probability P(θ*∈D) is no more than 0.5. Consequently, this approach is almost useless for any practical problems.
To further testify the above conclusion about the technique in Solution2, we consider the following stochastic chance-constrained programming, which relaxes the constraints in (2.4):
(2.8)min⁡x1+x2s.t.P{ax1+x2≥7,bx1+x2≥4}≥0.95×0.95xi≥0,i=1,2.
Generate 10000 samples of(a,b), replace (a,b) by the samples, and calculate all the deterministic programming problems. We can find that there is almost no solution satisfying the constraints of (2.8) among 10000 solutions.Table 1
The solutions and their probabilitiesθ*∈D.
NumberSolution(x1,x2)x1+x2P(θ*∈D)1(2.20, 2.62)4.820.372(2.13, 2.68)4.810.383(1.46, 3.40)4.860.454(1.13, 3.34)4.480.165(1.75, 2.50)4.260.116(0.95, 3.44)4.390.057(3.18, 1.66)4.840.378(0.95, 3.23)4.190.009(1.70, 3.03)4.740.3610(0.99, 3.39)4.810.38Solution 3.
In the following, we analyze a traditional but indirect approach, which is to convert SLP into an equivalent deterministic programming problem and obtain the optimal solution by some deterministic programming algorithms.
From the constraints in (2.8), we have
(2.9)P{ax1+x2≥7,bx1+x2≥4}=P{a≥7-x2x1,b≥4-x2x1}=P{a≥a*,b≥b*}=[(4-a*)(1-b*)]÷(3×23),
where a*=(7-x2)/x1,b*=(4-x2)/x1.
Therefore, the constraints in (2.8) are equivalent to
(2.10)[(4-7-x2x1)×(1-4-x2x1)]÷(3×23)≥0.9025,(4-7-x2x1)÷3≥0.9025,(1-4-x2x1)÷23≥0.9025,
and then the SCLP model (2.8) is equivalent to the following deterministic nonlinear programming:
(2.11)min⁡x1+x2s.t.1.805x12-(4x1+x2-7)(x1+x2-4)≤00.9025x1-4x1-(7-x2)3≤00.9025x1-x1-(4-x2)2/3≤0x1≥0,x2≥0.
Setting the initial valueX0=(0,0), then the optimal solution θ* and the optimal value f* of this nonlinear programming model are
(2.12)θ*=(3.36,2.84)T,f*=6.19,
and we can have the probability P(θ*∈D)=0.9025.From the above calculation and analysis of the three methods, we can conclude that the former two are not feasible regardless of their simpleness and the third one is valid. However, the third measure is indirect and the obtained equivalent deterministic programming is usually a nonlinear programming whose calculation is difficult and sometimes impossible. Therefore, it is necessary to find some novel direct computation method. Recently, a direct and effective algorithm, that is, genetic algorithms, has been put forward and rapidly developed to solve the stochastic chance-constrained linear programming, see, for example, the work of Liu in [9], of Ding et al. in [10], and of Ding and Sun in [11]. However, they are always designed aiming to solve every specified problem, depending on experiment and lacking the common theory basis.In next section, we propose a direct and universal algorithm of SCLP, simplex algorithm based on stochastic simulation, and build its theory basis. Then we design the detailed procedures of the algorithm, which are easily changed into the executable codes of software.


## 2.1. Stochastic Linear Programming (SLP) Model

Firstly, introduce the following SLP model(2.1)min⁡CTXs.t.AX≤bX≥0,
where X is an n dimensional vector to be determined, and A=(aij(ω))m×n,b=(bi(ω))m×1,C=(cj(ω))n×1, with ai,j(ω),bi(ω),cj(ω)(i=1,…,m,j=1,…,n) being random variables and ω∈Ω.


## 2.2. Stochastic Chance-Constrained Linear Programming (SCLP) Model

SCLP usually includes two sorts of models which can be formulated as follows:(2.2)min⁡CTXs.t.P(AX≤b)≥αX≥0,(2.3)min⁡CTXs.t.P(AiX≤bi)≥αi,i=1,2,…,mX≥0,
where Ai is the ith row of A, bi is the ith element of b, and αi is the ith confidence level of the constraints.


## 2.3. Computation Problem of SCLP

In this subsection, by solving a numerical example of SLP, the computation problem of SCLP is analyzed and several possible approaches are tried to obtain the optimal solution to the example.Example 2.1.
Let(a,b) be the random variables with uniform distribution in rectangle {1≤a≤4,1/3≤b≤1} and consider the following SLP problem:
(2.4)min⁡x1+x2s.t.ax1+x2≥7bx1+x2≥4x1≥0,x2≥0.Solution 1.
In (2.4), the expectations of the random variables a and b are easily derived. Therefore, a very simple idea is to replace the stochastic parameters a and b by their expectations, respectively, and solve the corresponding deterministic linear programming problem described as follows:
(2.5)min⁡x1+x2s.t.52x1+x2≥723x1+x2≥4x1≥0,x2≥0.
From (2.5), by deterministic linear programming algorithms, it is easy to derive the unique optimal solution x1*=18/11,x2*=32/11.
Now, in order to analyze the feasibility of this approach, we assumeθ*=(18/11,32/11), and, from the example, we can easily know the feasible region D={(x1,x2)∣ax1+x2≥7,bx1+x2≥4}. Then
(2.6)P(θ*∈D)=P(1811a+3211≥7,1811b+3211≥4)=P(a≥52,b≥23)=0.25.
That is, the probability of the solution θ* taking value in the feasible region D is only 0.25, which verifies that the method of directly replacing random parameters by their expectations is invaluable.Solution 2.
It is well known that the samples of random variables are also easily obtained; therefore, another simple technique is to produce some random samples of these random parameters and solve all the deterministic linear programming problems corresponding to the samples, and then we choose the best solution as the optimal solution. Now, produce 10 samples ofa and b in (2.4):
(2.7)a:(1.99,2.02,2.47,3.23,2.56,3.78,1.68,3.95,2.33,3.65),b:(0.63,0.62,0.41,0.58,0.85,0.59,0.47,0.80,0.57,0.62).
From Table 1, it is easy to see that the largest probability P(θ*∈D) is no more than 0.5. Consequently, this approach is almost useless for any practical problems.
To further testify the above conclusion about the technique in Solution2, we consider the following stochastic chance-constrained programming, which relaxes the constraints in (2.4):
(2.8)min⁡x1+x2s.t.P{ax1+x2≥7,bx1+x2≥4}≥0.95×0.95xi≥0,i=1,2.
Generate 10000 samples of(a,b), replace (a,b) by the samples, and calculate all the deterministic programming problems. We can find that there is almost no solution satisfying the constraints of (2.8) among 10000 solutions.Table 1
The solutions and their probabilitiesθ*∈D.
NumberSolution(x1,x2)x1+x2P(θ*∈D)1(2.20, 2.62)4.820.372(2.13, 2.68)4.810.383(1.46, 3.40)4.860.454(1.13, 3.34)4.480.165(1.75, 2.50)4.260.116(0.95, 3.44)4.390.057(3.18, 1.66)4.840.378(0.95, 3.23)4.190.009(1.70, 3.03)4.740.3610(0.99, 3.39)4.810.38Solution 3.
In the following, we analyze a traditional but indirect approach, which is to convert SLP into an equivalent deterministic programming problem and obtain the optimal solution by some deterministic programming algorithms.
From the constraints in (2.8), we have
(2.9)P{ax1+x2≥7,bx1+x2≥4}=P{a≥7-x2x1,b≥4-x2x1}=P{a≥a*,b≥b*}=[(4-a*)(1-b*)]÷(3×23),
where a*=(7-x2)/x1,b*=(4-x2)/x1.
Therefore, the constraints in (2.8) are equivalent to
(2.10)[(4-7-x2x1)×(1-4-x2x1)]÷(3×23)≥0.9025,(4-7-x2x1)÷3≥0.9025,(1-4-x2x1)÷23≥0.9025,
and then the SCLP model (2.8) is equivalent to the following deterministic nonlinear programming:
(2.11)min⁡x1+x2s.t.1.805x12-(4x1+x2-7)(x1+x2-4)≤00.9025x1-4x1-(7-x2)3≤00.9025x1-x1-(4-x2)2/3≤0x1≥0,x2≥0.
Setting the initial valueX0=(0,0), then the optimal solution θ* and the optimal value f* of this nonlinear programming model are
(2.12)θ*=(3.36,2.84)T,f*=6.19,
and we can have the probability P(θ*∈D)=0.9025.From the above calculation and analysis of the three methods, we can conclude that the former two are not feasible regardless of their simpleness and the third one is valid. However, the third measure is indirect and the obtained equivalent deterministic programming is usually a nonlinear programming whose calculation is difficult and sometimes impossible. Therefore, it is necessary to find some novel direct computation method. Recently, a direct and effective algorithm, that is, genetic algorithms, has been put forward and rapidly developed to solve the stochastic chance-constrained linear programming, see, for example, the work of Liu in [9], of Ding et al. in [10], and of Ding and Sun in [11]. However, they are always designed aiming to solve every specified problem, depending on experiment and lacking the common theory basis.In next section, we propose a direct and universal algorithm of SCLP, simplex algorithm based on stochastic simulation, and build its theory basis. Then we design the detailed procedures of the algorithm, which are easily changed into the executable codes of software.


## 3. Simplex Algorithm Based on Stochastic Simulation

In SCLP model, on one hand, the meaning of “min” is not clear because of C being random vector; on the other hand, it is difficult to judge the convexity of SCLP which is required by optimization theory. Therefore, the computation problem of the SCLP is very difficult. In this section, we propose a satisfying algorithm, that is, simplex algorithm based on stochastic simulation, to overcome this difficulty. Firstly, several problems in SCLP are handled by stochastic simulation. Then, we build the theory basis of this algorithm by theoretical analysis. Finally, the detailed procedures of this algorithm are designed.According to the theory of stochastic chance-constrained linear programming, we can transform (2.2) into the following programming model:(3.1)min⁡fs.t.P(CTX≤f)≥βP(AX≤b)≥αX≥0,
where f is the target function and α,β are the confidence levels of the target function and constraint, respectively.Similarly, model (2.3) can be transformed into(3.2)min⁡fs.t.P(CTX≤f)≥βP(AiX≤bi)≥αi,i=1,2,…,mX≥0.


### 3.1. Stochastic Simulation

#### 3.1.1. Judging the Chance Constraint

Consider the chance constraints in (3.1): (3.3)P{AX≤b}≥α,
where the random matrix A and the random vector b have a known joint probability distribution Φ(A,b). Then we check whether the chance constraint holds or not if we have a solution X by applying the stochastic simulation method (the Monte Carlo simulation). The algorithm is as follows.Algorithm 3.1.
Chance constraint judging algorithm.
Step 1.
SetN′=0.Step 2.
Sample a random vector{(ai,jk,bik),i=1,2,…,m,j=1,2,…,n} according to the joint probability distribution of Φ(A,b).Step 3.
CalculateAX; if AX≤b, then N′++.Step 4.
Repeat Step2 to Step 3  N times.Step 5.
IfN′/N≥α, return FEASIBLE, or else return INFEASIBLE.


#### 3.1.2. Handling the Target Function

Consider the target function with the random parameter vectorC: (3.4)P(CTX≤f)≥β.For any given vectorX, the minimum objective function f can always be found by stochastic simulation, and the algorithm is as follows.Algorithm 3.2.
Minimum target function searching algorithm.
Step 1.
Sample a random vector{Ck,k=1,2,…,N} according to the probability distribution of Φ(C).Step 2.
Compute{CkTX,k=1,2,…,N} and arrange the results according to the ascending order.Step 3.
SetN′ as the integer part of βN.Step 4.
Return theN′th largest element in the set {CkTX,k=1,2,…,N} as the estimation of f.


#### 3.1.3. Checking the Estimation Number

In order to check the estimate of probability, we test and return the random variabletj, such that (3.5)P{ω∈Ω∣tj-ĉj>0}≥γ,
where tj is a random variable, ĉj is a deterministic number, and γ is a prescribed confidence level.Algorithm 3.3.
Estimation number checking algorithm.
Step 1.
SetN′=0.Step 2.
Sample random vectors{tj(k),k=1,2,…,N} according to the probability distribution of Φ(tj).Step 3.
Iftj-ĉj>0, then N′++.Step 4.
Repeat Step2 to Step 3  N times.Step 5.
IfN′/N≥γ, return FEASIBLE and execute Step 6.Step 6.
LettingK′ be the integer part of γN, array {tj-ĉj,k=1,2,…,N} according to ascending order, return the (N-K′+1)′th element.Step 7.
Otherwise return INFEASIBLE.


### 3.2. Theoretical Analysis

In this subsection, we aim to build the theory basis of simplex algorithm based on stochastic simulation by theoretical analysis. To begin with, we recall the basic principle of deterministic linear programming. According to the convexity of deterministic linear programming, the programming has optimal solution if the feasible solutions are finite, and the optimal solution must be in the range of the feasible solutions (see the work of Zhang and Xu in [12]). Based on this theory, the simplex method is designed. Specifically, the basic principles of simplex method can be formulated as follows. First of all, find a feasible solution stochastically and check whether it is optimal; if it is not optimal, then we find another feasible solution which can improve the target function and check this solution again; repeat the above process until we find the optimal solution and the corresponding target value, or we can confirm that the programming does not have the optimal solution. According to these principles, we further consider the stochastic chance-constrained linear programming (3.1).In (3.1), replace AX≤b by AX+Y=b with Y≥0 according to the slack variable method and assume H=(A,EY), where EY is an m×m unit matrix. Then the constraint is equivalent to (3.6)P(HZ=b)≥α,
where Z=(XTYT)T. Still mark (CT0YT)T as CT, where 0Y is a m×1 zero vector, and then we have CTX=CTZ.Based on the above assumption, the main principle of simplex algorithm based on stochastic simulation can be described as follows. Firstly, find a base matrix fromH by random sample or Big M method (Big M method is the most direct method and the base matrix is deterministic unit matrix). Secondly, search a feasible basic solution satisfying the chance constraint according to stochastic simulation, calculate the corresponding b (denoted by b̂), then check whether it is the optimal solution by the sample value of C (denoted by Ĉ), which is calculated by stochastic simulation method; if it is not optimal, according to the improvement principle, we can find the sample value of H (denoted by Ĥ). Thirdly, solve the deterministic programming defined by Ĥ,b̂ and Ĉ by simplex algorithm, and we can obtain a solution X which can improve the value of target function. Fourthly, check whether the solution satisfies the constraint. If it does not, change Ĥ,b̂, and Ĉ and check once more until we derive an improved solution satisfying the constraint. Repeat the above steps until we find the optimal solution or are sure that there is no optimal solution (infinite optimal solutions).AssumeZ¯, B, b̂, U (H=(B,U)) are the initial feasible solution, the feasible base, a sample value of b (b̂ is determined by checking whether B-1b̂ satisfies the constraint through stochastic simulation), and the nonbase matrix, respectively, and then we have (3.7)Z¯=(B-1b̂0)=(b¯0),CTZ¯=(CBTCUT)(b¯0)=CBTb¯,
with CB and CU being corresponding to the base variable and the nonbase variable of C, respectively.Computing the value of target function by stochastic simulation and chance constraint(3.8)P(CT(b¯0)≤f¯)≥β,
it is easy to obtain the sample value Ĉ of C and the target value f¯=ĈBTb¯.SetZ=(ZBTZUT)T as being any feasible solution and consider the following stochastic chance constrained programming:(3.9)min⁡f=ĈTZs.t.P(HZ=b̂)≥αZ=(XTYT)T≥0.From (3.9) and the above discussion, we have (3.10)P(ZB=B-1b̂-B-1UZU)≥α,f=ĈTZ=ĈBTZB+ĈUTZU.
therefore,(3.11)P(f=ĈBTb¯-(ĈBTB-1U-ĈUT)ZU)≥α.
Let H=(h1,h2,…,hn) and JU be the lower label set of nonbase variables and (3.12)tj-ĉj=ĈBTB-1hj-ĉj,j∈JU,
where tj-ĉj is defined as the estimation number. Then we can derive (3.13)P(f=f¯-∑j∈JU(tj-ĉj)zj)≥α.Then, the programming (3.9) is converted into(3.14)min⁡f=ĈTZs.t.P(f=f¯-∑j∈JU(tj-ĉj)zj)≥αP(ZB=B-1b̂-B-1UZU)≥αZ=(XTYT)T≥0,
where f¯=ĈBTB-1b̂ with b̂ and Ĉ being obtained by stochastic simulation in (3.1).If the feasible base solutions in programming (3.14) are nondegenerative, from the above discussion, we can have the following theorem.Theorem 3.4.
IfP(tj-ĉj≤0,∀j∈JU)≥α, then Z¯ is the optimal solution of SCLP (3.14) and is denoted as Z*.Proof.
SinceP(tj-ĉj≤0)≥α holds for all j∈JU, there is no new feasible solution satisfying the constraint and reducing target function value. Therefore Z¯ must be the optimal solution.Given a probability space(Ω,ℱ,P), assume (3.15)Ω1={ω∈Ω∣tj-ĉj≤0,∀j∈JU},Ω2={ω∈Ω∣tj-ĉj>0,∀j∈JU}.
Then it is easy to know that Ω=Ω1⋃Ω2,Ω1⋂Ω2=∅, and P(Ω1)+P(Ω2)=1. According to Theorem 3.4, if P(Ω1)≥α does not hold, there may be other new feasible solutions satisfying the constraint and reducing the target function value. Since P(Ω1)<α is equivalent to P(Ω2)≥α, we can compute estimation number P(Ω2)≥1-α by stochastic simulation and return some sample values of tj and also obtain some sample values of hj which is denoted as Ĥ. Then we can derive a feasible solution by solving the following deterministic linear programming by simplex method (see the work of Zhang and Xu in [12]).(3.16)min⁡ĈTZs.t.ĤZ=b̂Z≥0.
Now, we can begin solve the SCLP problem (3.1) by simplex algorithm based on stochastic simulation, and the steps are as follows.Firstly, check the initial feasible solutionZ¯. If all of the estimation numbers tj-ĉj≤0 hold, then Z¯ is the optimal solution; if tj-ĉj>0 and B-1hk≤0, there is no optimal solution; otherwise, if tj-ĉj>0 and some of elements of B-1hk are positive, there must be a new feasible solution Ẑ to reduce the target function value.Secondly, check whetherẐ satisfies the constraint of (3.1) by stochastic simulation. If it does, continue to check whether it is the optimal one; if it does not, change a new Ĥ and the corresponding new b̂ and repeat the above checking.Finally, repeat the above two steps. The number of basic feasible solutions is finite; therefore, it is sure that we can find an optimal solution or the programming problem has no optimal solution.In order to find a new basic feasible solutionẐ, let the vector hk (k is the biggest estimation number of tk-ĉk) enter the base vector and change the nonbase vector into a base vector. Set (3.17)zk=min⁡{b¯ih¯ik∣h¯ik>0}=∧b¯ih¯ik,i=1,2,…,m
and let zBr=0. Then we obtain a new basic feasible solution (3.18)Ẑ=(ẑB1,…,ẑBr-1,0,ẑBr+1,…,ẑBm,0,…,zk,…,0)T,
where ẑBi=b¯i-h¯ikzk≥0,(i=1,…,r-1,r+1,…,m) and ẑBr=zk>0. This is to say ẐB≥0. Then, check whether Ẑ satisfies the constraints by stochastic simulation. If it does, we figure out a new feasible solution.


### 3.3. Computation Procedure

In this subsection, we design the detailed steps for the simplex algorithm based on stochastic simulation according to the above algorithm analysis. These procedures can easily be converted into the executable codes of some software tools.Algorithm 3.5.
Simplex algorithm based on stochastic simulation.
Step 1.
Find an initial feasible baseB.Step 2.
Find a basic feasible solution satisfying the chance constraints and the corresponding sample value ofb (denoted by b̂).Step 3.
ComputingZB=B-1b̂=∧b¯ and the target function value f¯=ĈBTb¯ by stochastic simulation, return the sample value Ĉ.Step 4.
CheckP(tj-ĉj>0)>1-α by stochastic simulation, and return sample value tj which satisfies the constraint. Produce a group of Ĥ by stochastic simulation.Step 5.
Select anĤ, and calculate the estimation number by tj-ĉj=ĈBTB-1ĥj-ĉj. Determine the lower label k by tk-ĉk=max{tj-ĉj∣j=1,2,…,n}, and then let zk enter the base vector.Step 6.
Iftk-ĉk≤0, end the procedure. Then, the basic feasible solution is the optimal one and the target function can be calculated as f¯=ĈBTb¯, or else, go to Step 7.Step 7.
Calculateh¯k=B-1ĥk; if h¯k≤0, end the steps. We can get the programming is infinite, or else, go to Step 8.Step 8.
Calculate the ratio(3.19)zk=min⁡{b¯ih¯ik∣h¯ik>0}=∧b¯ih¯ik
to find the lower label r, and set ẑBr=0.Step 9.
ReplaceĥBr by hk, and we have a new base. Then, compute the new basic feasible solution Ẑ.Step 10.
Check whetherẐ satisfies the chance constraints by stochastic simulation. If it does, go to Step 2; otherwise, go to Step 5. If, for all of Ĥ obtained in Step 4, Ẑ does not satisfy the chance constraint, go to Step 2 to find a new b̂.


## 3.1. Stochastic Simulation

### 3.1.1. Judging the Chance Constraint

Consider the chance constraints in (3.1): (3.3)P{AX≤b}≥α,
where the random matrix A and the random vector b have a known joint probability distribution Φ(A,b). Then we check whether the chance constraint holds or not if we have a solution X by applying the stochastic simulation method (the Monte Carlo simulation). The algorithm is as follows.Algorithm 3.1.
Chance constraint judging algorithm.
Step 1.
SetN′=0.Step 2.
Sample a random vector{(ai,jk,bik),i=1,2,…,m,j=1,2,…,n} according to the joint probability distribution of Φ(A,b).Step 3.
CalculateAX; if AX≤b, then N′++.Step 4.
Repeat Step2 to Step 3  N times.Step 5.
IfN′/N≥α, return FEASIBLE, or else return INFEASIBLE.


### 3.1.2. Handling the Target Function

Consider the target function with the random parameter vectorC: (3.4)P(CTX≤f)≥β.For any given vectorX, the minimum objective function f can always be found by stochastic simulation, and the algorithm is as follows.Algorithm 3.2.
Minimum target function searching algorithm.
Step 1.
Sample a random vector{Ck,k=1,2,…,N} according to the probability distribution of Φ(C).Step 2.
Compute{CkTX,k=1,2,…,N} and arrange the results according to the ascending order.Step 3.
SetN′ as the integer part of βN.Step 4.
Return theN′th largest element in the set {CkTX,k=1,2,…,N} as the estimation of f.


### 3.1.3. Checking the Estimation Number

In order to check the estimate of probability, we test and return the random variabletj, such that (3.5)P{ω∈Ω∣tj-ĉj>0}≥γ,
where tj is a random variable, ĉj is a deterministic number, and γ is a prescribed confidence level.Algorithm 3.3.
Estimation number checking algorithm.
Step 1.
SetN′=0.Step 2.
Sample random vectors{tj(k),k=1,2,…,N} according to the probability distribution of Φ(tj).Step 3.
Iftj-ĉj>0, then N′++.Step 4.
Repeat Step2 to Step 3  N times.Step 5.
IfN′/N≥γ, return FEASIBLE and execute Step 6.Step 6.
LettingK′ be the integer part of γN, array {tj-ĉj,k=1,2,…,N} according to ascending order, return the (N-K′+1)′th element.Step 7.
Otherwise return INFEASIBLE.


## 3.1.1. Judging the Chance Constraint

Consider the chance constraints in (3.1): (3.3)P{AX≤b}≥α,
where the random matrix A and the random vector b have a known joint probability distribution Φ(A,b). Then we check whether the chance constraint holds or not if we have a solution X by applying the stochastic simulation method (the Monte Carlo simulation). The algorithm is as follows.Algorithm 3.1.
Chance constraint judging algorithm.
Step 1.
SetN′=0.Step 2.
Sample a random vector{(ai,jk,bik),i=1,2,…,m,j=1,2,…,n} according to the joint probability distribution of Φ(A,b).Step 3.
CalculateAX; if AX≤b, then N′++.Step 4.
Repeat Step2 to Step 3  N times.Step 5.
IfN′/N≥α, return FEASIBLE, or else return INFEASIBLE.


## 3.1.2. Handling the Target Function

Consider the target function with the random parameter vectorC: (3.4)P(CTX≤f)≥β.For any given vectorX, the minimum objective function f can always be found by stochastic simulation, and the algorithm is as follows.Algorithm 3.2.
Minimum target function searching algorithm.
Step 1.
Sample a random vector{Ck,k=1,2,…,N} according to the probability distribution of Φ(C).Step 2.
Compute{CkTX,k=1,2,…,N} and arrange the results according to the ascending order.Step 3.
SetN′ as the integer part of βN.Step 4.
Return theN′th largest element in the set {CkTX,k=1,2,…,N} as the estimation of f.


## 3.1.3. Checking the Estimation Number

In order to check the estimate of probability, we test and return the random variabletj, such that (3.5)P{ω∈Ω∣tj-ĉj>0}≥γ,
where tj is a random variable, ĉj is a deterministic number, and γ is a prescribed confidence level.Algorithm 3.3.
Estimation number checking algorithm.
Step 1.
SetN′=0.Step 2.
Sample random vectors{tj(k),k=1,2,…,N} according to the probability distribution of Φ(tj).Step 3.
Iftj-ĉj>0, then N′++.Step 4.
Repeat Step2 to Step 3  N times.Step 5.
IfN′/N≥γ, return FEASIBLE and execute Step 6.Step 6.
LettingK′ be the integer part of γN, array {tj-ĉj,k=1,2,…,N} according to ascending order, return the (N-K′+1)′th element.Step 7.
Otherwise return INFEASIBLE.


## 3.2. Theoretical Analysis

In this subsection, we aim to build the theory basis of simplex algorithm based on stochastic simulation by theoretical analysis. To begin with, we recall the basic principle of deterministic linear programming. According to the convexity of deterministic linear programming, the programming has optimal solution if the feasible solutions are finite, and the optimal solution must be in the range of the feasible solutions (see the work of Zhang and Xu in [12]). Based on this theory, the simplex method is designed. Specifically, the basic principles of simplex method can be formulated as follows. First of all, find a feasible solution stochastically and check whether it is optimal; if it is not optimal, then we find another feasible solution which can improve the target function and check this solution again; repeat the above process until we find the optimal solution and the corresponding target value, or we can confirm that the programming does not have the optimal solution. According to these principles, we further consider the stochastic chance-constrained linear programming (3.1).In (3.1), replace AX≤b by AX+Y=b with Y≥0 according to the slack variable method and assume H=(A,EY), where EY is an m×m unit matrix. Then the constraint is equivalent to (3.6)P(HZ=b)≥α,
where Z=(XTYT)T. Still mark (CT0YT)T as CT, where 0Y is a m×1 zero vector, and then we have CTX=CTZ.Based on the above assumption, the main principle of simplex algorithm based on stochastic simulation can be described as follows. Firstly, find a base matrix fromH by random sample or Big M method (Big M method is the most direct method and the base matrix is deterministic unit matrix). Secondly, search a feasible basic solution satisfying the chance constraint according to stochastic simulation, calculate the corresponding b (denoted by b̂), then check whether it is the optimal solution by the sample value of C (denoted by Ĉ), which is calculated by stochastic simulation method; if it is not optimal, according to the improvement principle, we can find the sample value of H (denoted by Ĥ). Thirdly, solve the deterministic programming defined by Ĥ,b̂ and Ĉ by simplex algorithm, and we can obtain a solution X which can improve the value of target function. Fourthly, check whether the solution satisfies the constraint. If it does not, change Ĥ,b̂, and Ĉ and check once more until we derive an improved solution satisfying the constraint. Repeat the above steps until we find the optimal solution or are sure that there is no optimal solution (infinite optimal solutions).AssumeZ¯, B, b̂, U (H=(B,U)) are the initial feasible solution, the feasible base, a sample value of b (b̂ is determined by checking whether B-1b̂ satisfies the constraint through stochastic simulation), and the nonbase matrix, respectively, and then we have (3.7)Z¯=(B-1b̂0)=(b¯0),CTZ¯=(CBTCUT)(b¯0)=CBTb¯,
with CB and CU being corresponding to the base variable and the nonbase variable of C, respectively.Computing the value of target function by stochastic simulation and chance constraint(3.8)P(CT(b¯0)≤f¯)≥β,
it is easy to obtain the sample value Ĉ of C and the target value f¯=ĈBTb¯.SetZ=(ZBTZUT)T as being any feasible solution and consider the following stochastic chance constrained programming:(3.9)min⁡f=ĈTZs.t.P(HZ=b̂)≥αZ=(XTYT)T≥0.From (3.9) and the above discussion, we have (3.10)P(ZB=B-1b̂-B-1UZU)≥α,f=ĈTZ=ĈBTZB+ĈUTZU.
therefore,(3.11)P(f=ĈBTb¯-(ĈBTB-1U-ĈUT)ZU)≥α.
Let H=(h1,h2,…,hn) and JU be the lower label set of nonbase variables and (3.12)tj-ĉj=ĈBTB-1hj-ĉj,j∈JU,
where tj-ĉj is defined as the estimation number. Then we can derive (3.13)P(f=f¯-∑j∈JU(tj-ĉj)zj)≥α.Then, the programming (3.9) is converted into(3.14)min⁡f=ĈTZs.t.P(f=f¯-∑j∈JU(tj-ĉj)zj)≥αP(ZB=B-1b̂-B-1UZU)≥αZ=(XTYT)T≥0,
where f¯=ĈBTB-1b̂ with b̂ and Ĉ being obtained by stochastic simulation in (3.1).If the feasible base solutions in programming (3.14) are nondegenerative, from the above discussion, we can have the following theorem.Theorem 3.4.
IfP(tj-ĉj≤0,∀j∈JU)≥α, then Z¯ is the optimal solution of SCLP (3.14) and is denoted as Z*.Proof.
SinceP(tj-ĉj≤0)≥α holds for all j∈JU, there is no new feasible solution satisfying the constraint and reducing target function value. Therefore Z¯ must be the optimal solution.Given a probability space(Ω,ℱ,P), assume (3.15)Ω1={ω∈Ω∣tj-ĉj≤0,∀j∈JU},Ω2={ω∈Ω∣tj-ĉj>0,∀j∈JU}.
Then it is easy to know that Ω=Ω1⋃Ω2,Ω1⋂Ω2=∅, and P(Ω1)+P(Ω2)=1. According to Theorem 3.4, if P(Ω1)≥α does not hold, there may be other new feasible solutions satisfying the constraint and reducing the target function value. Since P(Ω1)<α is equivalent to P(Ω2)≥α, we can compute estimation number P(Ω2)≥1-α by stochastic simulation and return some sample values of tj and also obtain some sample values of hj which is denoted as Ĥ. Then we can derive a feasible solution by solving the following deterministic linear programming by simplex method (see the work of Zhang and Xu in [12]).(3.16)min⁡ĈTZs.t.ĤZ=b̂Z≥0.
Now, we can begin solve the SCLP problem (3.1) by simplex algorithm based on stochastic simulation, and the steps are as follows.Firstly, check the initial feasible solutionZ¯. If all of the estimation numbers tj-ĉj≤0 hold, then Z¯ is the optimal solution; if tj-ĉj>0 and B-1hk≤0, there is no optimal solution; otherwise, if tj-ĉj>0 and some of elements of B-1hk are positive, there must be a new feasible solution Ẑ to reduce the target function value.Secondly, check whetherẐ satisfies the constraint of (3.1) by stochastic simulation. If it does, continue to check whether it is the optimal one; if it does not, change a new Ĥ and the corresponding new b̂ and repeat the above checking.Finally, repeat the above two steps. The number of basic feasible solutions is finite; therefore, it is sure that we can find an optimal solution or the programming problem has no optimal solution.In order to find a new basic feasible solutionẐ, let the vector hk (k is the biggest estimation number of tk-ĉk) enter the base vector and change the nonbase vector into a base vector. Set (3.17)zk=min⁡{b¯ih¯ik∣h¯ik>0}=∧b¯ih¯ik,i=1,2,…,m
and let zBr=0. Then we obtain a new basic feasible solution (3.18)Ẑ=(ẑB1,…,ẑBr-1,0,ẑBr+1,…,ẑBm,0,…,zk,…,0)T,
where ẑBi=b¯i-h¯ikzk≥0,(i=1,…,r-1,r+1,…,m) and ẑBr=zk>0. This is to say ẐB≥0. Then, check whether Ẑ satisfies the constraints by stochastic simulation. If it does, we figure out a new feasible solution.


## 3.3. Computation Procedure

In this subsection, we design the detailed steps for the simplex algorithm based on stochastic simulation according to the above algorithm analysis. These procedures can easily be converted into the executable codes of some software tools.Algorithm 3.5.
Simplex algorithm based on stochastic simulation.
Step 1.
Find an initial feasible baseB.Step 2.
Find a basic feasible solution satisfying the chance constraints and the corresponding sample value ofb (denoted by b̂).Step 3.
ComputingZB=B-1b̂=∧b¯ and the target function value f¯=ĈBTb¯ by stochastic simulation, return the sample value Ĉ.Step 4.
CheckP(tj-ĉj>0)>1-α by stochastic simulation, and return sample value tj which satisfies the constraint. Produce a group of Ĥ by stochastic simulation.Step 5.
Select anĤ, and calculate the estimation number by tj-ĉj=ĈBTB-1ĥj-ĉj. Determine the lower label k by tk-ĉk=max{tj-ĉj∣j=1,2,…,n}, and then let zk enter the base vector.Step 6.
Iftk-ĉk≤0, end the procedure. Then, the basic feasible solution is the optimal one and the target function can be calculated as f¯=ĈBTb¯, or else, go to Step 7.Step 7.
Calculateh¯k=B-1ĥk; if h¯k≤0, end the steps. We can get the programming is infinite, or else, go to Step 8.Step 8.
Calculate the ratio(3.19)zk=min⁡{b¯ih¯ik∣h¯ik>0}=∧b¯ih¯ik
to find the lower label r, and set ẑBr=0.Step 9.
ReplaceĥBr by hk, and we have a new base. Then, compute the new basic feasible solution Ẑ.Step 10.
Check whetherẐ satisfies the chance constraints by stochastic simulation. If it does, go to Step 2; otherwise, go to Step 5. If, for all of Ĥ obtained in Step 4, Ẑ does not satisfy the chance constraint, go to Step 2 to find a new b̂.


## 4. A Numerical Example

In this section, a simulation example is presented to illustrate the feasibility and effectiveness of the simplex algorithm based on stochastic simulation developed in this paper.Example 4.1.
Consider the stochastic chance-constrained linear programming (2.8) again.
According to Big M method, (2.8) is equivalent to
(4.1)min⁡(x1+x2+Mx3+Nx4)s.t.P{ax1+x2+x3-x5=7bx1+x2+x4-x6=4}≥0.95×0.95xi≥0,i=1,2,…,6,
where x4 and x5 are slack variables.
According to Algorithm3.5, the above SCLP (4.1) can be solved by using MATLAB toolbox and the optimal solution and optimal target function value can be obtained as follows:
(4.2)X*=(3.2010,2.9245)T,f*=6.1255.
Moreover, we can find that, for the above optimal solution,
(4.3)P{ax1+x2≥7,bx1+x2≥4}=0.9035≥0.95×0.95.Remark 4.2.
In SCLP, the coefficients,X and f, are all random variables; therefore, the optimal solution may be different for different computation.  In Table 2, we obtain three couples of optimal solutions by executing Algorithm 3.5 three times.     From the above results, we can know that every optimal solution satisfies the constraints, and there is a different confidence level corresponding to the optimal solution.Table 2
The solutions and their probabilitiesθ*∈D.
NumberSolution(x1,x2)f*P(θ*∈D)1(3.2802, 2.8907)6.17090.90912(3.2708, 2.8899)6.16070.90613(3.4615, 2.8299)6.29140.9252Remark 4.3.
From Solution2 of Example 2.1, it is very difficult to obtain an effective result even if 10000 random samples are utilized. But we can obtain a satisfactory optimal solution using less than 20 samples in our experiment by simplex algorithm based on stochastic simulation. So it is clear that Algorithm 3.5, that is, simplex algorithm based on stochastic simulation, is effective and much better than the method in Solution 2 to Example 2.1.Remark 4.4.
Apparently, simplex algorithm based on stochastic simulation is a direct method and can be applied to any SCLP problem. Thereby, Algorithm3.5 is also better than the traditional approach in Solution 3 to Example 2.1, which is indirect measure and only effective for some special cases.Example 4.5.
In this example, we consider a typical optimal decision problem in oil refinery production (Kall and Wallace [5]). An oil refinery factory refines two kinds of crude oil (denoted by raw1 and raw2) and provides gas (denoted by prod1) for the gas company and burning oil (denoted by prod2) for power company. A plan is needed one week before production. Assume that the yield π(raw1,prod1) of gas by raw1 and yield π(raw2,prod2) of burning oil by raw2 are random and yields of other products are deterministic. Therefore, set
(4.4)π(raw1,prod1)=2+ξ1,π(raw1,prod2)=3,π(raw2,prod1)=6,π(raw2,prod2)=3.4+ξ2,
where ξ1is uniform distribution 𝒰(-0.8,0.8) and ξ2is exponential distribution ℰ𝒳𝒫(0.4). And let the requirement h1 of gas and h2 of burning oil be also random variables
(4.5)h1=180+η1,h2=162+η2,
where η1is normal distribution 𝒩(0,12) and η2is normal distribution 𝒩(0,3).The amount of expendingraw1 and raw2 are denoted as x1 and x2, respectively, and the unit prices of raw1 and raw2 are c1=2 and c2=3, receptively. Therefore, the total cost is 2x1+3x2. Again assume the production capability (the largest amount of raw consumed) per week is 100, that is, the constraint (4.6)x1+x2≤100.
This problem has been dealt with by using a two-stage complement model in [5], and now we solve it by our simplex algorithm based on stochastic simulation developed in this paper.As we know that the decision(x1,x2) is made a week before production and cannot be changed during the next week. Moreover, in these decisions, confidence levels are necessary; these are (4.7)P((2+ξ1)x1+6x2≥180+η1)≥α1,P(3x1+(3.4-ξ2)x2≥162+η2)≥α2.
According to the decision principles: satisfying customer and minimizing loss, we can have the following SCLP problem: (4.8)min⁡f=2x1+3x2s.t.P((2+ξ1)x1+6x2≥180+η1)≥α1P(3x1+(3.4-ξ2)x2≥162+η2)≥α2x1+x2≤100x1,x2≥0.
If the confidence levels α1 and α2 are 0.8 and 0.7, respectively, by our simplex algorithm based on stochastic simulation, we can obtain an optimal solution(4.9)X*=(33.0944,21.7716)T,f*=131.5035,P((2+ξ1)x1+6x2≥180+η1)≥0.8149,P(3x1+(3.4-ξ2)x2≥162+η2)≥0.715.Remark 4.6.
In this SCLP problem, the random variablesξ1,ξ2,η1, and η2 obey uniform, exponential and normal distributions, respectively; therefore the joint distribution of them is so complicated that it is nearly impossible to find Φ-1.  Therefore the programming problem is difficult to be solved by using the approach in Solution 3 to Example 2.1.Remark 4.7.
Genetic algorithm based on stochastic simulation has been put forward by Iwamura and Liu [6] for SCL problems, which is also a direct and effective method. Supposing the scale popsize=30, the cross probability Pc=0.5, the mutation probability Pm=0.05, and the parameter in the rank-based evaluation function α=0.05, after running 500 generations, we obtain an optimal solution
(4.10)X*=(31.95,22.65)T,f*=131.85,P((2+ξ1)x1+6x2≥180+η1)≥0.8859,P(3x1+(3.4-ξ2)x2≥162+η2)≥0.6849.
Comparing (4.9) and (4.10), we can see both of them satisfying all of the constraints. However, (4.9) can obtain a less minimum value; this is a better result than (4.10). We should notice that their confidence levels are different. So simplex algorithm based on stochastic simulation is a practicable method.


## 5. Conclusion

This paper has studied the computation problem of stochastic chance-constrained linear programming and proposed a novel algorithm, simplex algorithm based on stochastic simulation. By a numerical example, several simple approaches have been tried to solve the SCLP problem, and the disadvantages of these methods have been analyzed. By theoretical analysis, the theory basis of the simplex algorithm based on stochastic simulation has been built and a theorem has been proved. Then the detailed procedures of the proposed algorithm have been designed, which are easily executed by some software tools. Finally, by two examples, the introduced algorithm has been verified to be better than the approaches used in Example2.1 and be more effective than genetic algorithm based on stochastic simulation.Based on the algorithm proposed in this paper, some possible further research topics include (i) the direct and universal algorithm for uncertain programming, such as fuzzy programming and nonlinear stochastic programming, see, for example, [10, 13, 14]; (ii) the control and state estimate problems, see, for example, [3, 15–19] and the references therein.



---

*Source: 139271-2011-09-07.xml*