Solution of Stochastic Quadratic Programming with Imperfect Probability Distribution Using Nelder-Mead Simplex Method ()
1. Introduction
Stochastic programming is an important method to solve decision problems in random environment. It was proposed by Dantzig, an American economist in 1956 [1] . Currently, the main method to solve the stochastic programming is to transform the stochastic programming into its own deterministic equivalence class and using the existing deterministic planning method to solve it. According to different research problems, stochastic programming mainly consists of three problems: distribution problem, expected value problem, and probabilistic constraint programming problem. Classic stochastic programming with recourse is a type of expected value problem, modeling based on a two-stage decision-making method. It is a method by making decisions before and after observing the value of a variable. With regard to the theory and methods of two-stage stochastic programming, a very systematic study has been conducted and many important solutions have been proposed [2] . In these methods, the dual decomposition L-shape algorithm established in the literature [3] is the most effective algorithm for solving two-stage stochastic programming. It is based on the duality theory, and the algorithm converges to the optimal solution by determining the feasible cutting plane and optimal cutting, and solving the main problem step by step. This method is essentially an external approximation algorithm that can effectively solve the large-scale problems that occur after the stochastic programming is transformed into deterministic mathematical programming. Abaffy and Allevi present a modified version of the L-shaped method in [4] , used to solve two-stage stochastic linear programs with fixed recourse. This method can apply class attributes and special structures to a polyhedron process to solve a certain type of large-scale problems, which greatly reduces the number of arithmetic operations.
While the stochastic programming is transformed into the corresponding equivalence classes, it is generally a nonlinear equation. In recent years, with the introduction of new theories and methods for solving nonlinear equations, especially the infinite dimensional variational inequality theories and the application of smoothness techniques that have received widespread attention in recent years [5] [6] [7] , a stochastic programming solution method based on nonlinear equation theory is proposed. Chen X. expressed the two-stage stochastic programming as a deterministic equivalence problem in the literature [8] , and transformed it into a nonlinear equation problem by introduced Lagrange multiplier. By using the B-differentiable properties of nonlinear functions, a Newton method for solving stochastic programming is proposed. Under certain conditions, the global convergence and local super-linear convergence of the algorithm are proved.
In general, stochastic programming is based on the complete information about probability distribution, but in practical situation, due to the lack of historical data and statistical theory, it is impossible to obtain complete information of the probability distribution, and can only get partial information in fact. In order to solve this problem, literature [9] and [10] based on fuzzy theory, under the condition that the membership function of certain parameters of the probability distribution is known, the method of determining the two-stage recourse function is given , two-stage and multi-stage stochastic programming problems are distributed and discussed. In [11] , based on the linear partial information (LPI) theory of Kofler [12] , a class of two-stage stochastic programming with recourse is established, and an L-shape method based on quadratic programming is given. Based on the literature [8] and literature [11] , this paper establishes a two-stage stochastic programming model under incomplete probability distribution information based on LPI theory, and presents an improved Nelder-Mead solution method. Experiment shows the algorithm is effective.
2. The Model of Stochastic Quadratic Programming with Imperfect Probability Distribution
Let
be a probability space,
is a stochastic vector in this space,
be symmetric positive semi-definite and
be symmetric positive definite. We consider following problem:
(1)
where
(2)
(3)
Here
are decision variables,
,
are fixed matrices,
is a random vector with support
,
is the probability distributions of limited sample set
, that is
. Assumed that the probability distribution of random variable has the following linear partial information:
Here
are fixed matrices, assumed the set of probability distributions
is a polyhedron. Thus the two-stage function can be written as
(4)
We call Equations (1)-(3) stochastic quadratic programming with recourse models under LPI.
Chen established a similar stochastic quadratic programming model in [8] , but assumed that the probability distribution is completely known, that is the “Max” symbol in Equation (2) does not appear. The above model is a new stochastic quadratic programming model based on LPI theory to solve the stochastic programming problem with incomplete information probability distribution. Since
is the convex function about y ( [8] ),
is also the convex function about y (see [13] ), and then the problems (1)-(3) essentially belong to the convex programming problem. Obviously the recourse function is not differential, so the Newton method proposed in [8] is no longer applicable. In order to solve this problem, we design a solution based on the improved Nelder-Mead method. The experimental results show that the method is effective.
3. Modified Nelder-Mead
The Nelder-Mead method (NM) [14] was originally a direct optimization algorithm developed for solving the nonlinear programming, NM algorithm belongs to the modified polyhedron method in nature. It searches for the new solution by reflecting the extreme point with the worst function value through the centroid of the remaining extreme points. Experimental shows, compared to random search, the algorithm can find the optimal solution more efficiently. The NM algorithm does not require any gradient information of the function during the entire optimization procedures, it can handle problems for which the gradient does not exist everywhere. NM allows the simplex to rescale or change its shape based on the local behavior of the response function. When the newly-generated point has good quality, an extension step will be taken in the hope that a better solution can be found. On the other hand, when the newly-generated solution is of poor quality, a contraction step will be taken, restricting the search on a smaller region. Since NM determines its search direction only by comparing the function values, it is insensitive to small inaccuracies in function values.
The classic NM method has several disadvantages in the search process. First, the convergence speed of the algorithm depends too much on the choice of initial polyhedron. Indeed, a too small initial simplex can lead to a local search, consequently the NM may convergent to a local solution. Second, NM might perform the shrink step frequently and in turn reduce the size of simplex to the greatest extent. Consequently, the algorithm can converge prematurely at a non-optimal solution.
Chang [15] propose a new variant of Nelder-Mead, called Stochastic Nelder-Mead simplex method (SNM). This method seeks the optimal solution by gradually increasing the sample size during the iterative process of the algorithm, which not only can effectively save the calculation time, but also can increase the adaptability of the algorithm to prevent premature convergence of the algorithm. This article refers to the design idea of [15] and adds an adaptive random search process to solve problems (1)-(3) in the NM algorithm. The specific process of the algorithm is described as follows:
Firstly, by attaching a Lagrange multiplier vector
, convex problems (1)-(3) can be written as an unconstrained problem:
(5)
, let
be the
points of n-dimensional space of
, which are not on the same plane. Let
represent the points that have the highest, second highest, and lowest estimates of function values,
is
the centroid of all vertices other than
,
Reflection: since
is the vertex with the higher value among the vertices, we can expect to find a lower value at the reflection of
in the opposite face formed by all vertices
except
. Generate a new point
, by reflecting
through
according to
.
If the reflected point is better than the second worst, but not better than the best, i.e.
, then obtain a new simplex by replacing the worst point
with the reflected point
.
Expansion: if the reflected point is the best point so far, we can expect to find interesting values along the direction from
to
, that is if
, then search direction is favorable, compute the expansion point using
.
If the expanded point is better than the reflected point, that is
, then replace
by
, otherwise, obtain a new simplex by replacing the worst point
with the reflected point
.
Contraction: here it is certain that
, in this case, we can expect that a better value will be inside the simplex formed by all the vertices, then the simplex contracts.
1) If
, the contracted point is determined by
with
, if
, the contraction is accepted, Replaced
by
.
2) If
, the contracted point is determined by
, if
, the contraction is accepted. Replaced
by
.
Shrink: although a failed contraction is much rarer, it may happen in some case, In that case, generally we contract towards the lowest point in the expectation of finding a simpler landscape. Replace all points except the best point
with
. This article uses the following process: when contraction fails, using a random search process to generate new points based on fitness of function. Let fitness function be
, where M is a fully large number, calculate the probability of obtaining
by
. According to the roulette algorithm, get a new point by randomly searched in the neighborhood of the point corresponding to the probability interval. The neighborhood
of
is defined as
Algorithm termination condition: There are different criteria to determine the termination conditions of the NM algorithm in practice, in this paper, we use
as our convergence criterion.
Parameters choice: The polyhedron transform in the NM algorithm mainly includes four parameters,
for reflection,
for expansionand
for contraction, assumed they satisfy the following constraints:
Algorithm
Choose
, Convergence criterion
.
Step 1 calculate the function values of n+1 points, rank all points and identify
, find
, the centroid of all vertices other than
, generating a new point
by reflecting point
through
according to reflection rule
;
Step 2 if
, then the reflection point is expanded using the expansion rule
, if
, then replace
by
, otherwise, let
replaced
, return to Step 4;
Step 3 if
, then let
replaced
, go to Step 4, otherwise return to Step 5;
Step 4 if the convergence criterion is met, stop the iteration, otherwise, return to Step 1;
Step 5 if
, then the simplex contracts.
(i) if
, the contraction point is determined by calculate
,
(ii) if
, the contraction point is determined by calculate
.
In these case, if
, the contraction accepted, If contraction is accepted, let
replaced
, return to Step 4;
Step 6 when all previous Step s fail, we use adaptive random search to generate new points, then return to Step 1.
4. Numerical Experiment
Consider the problem (1)-(3) in which
, and
Assumed stochastic vector
is three-dimensional vector,
and
are independent of each other. We use MATLAB to randomly generate one hundred values for each other. Let
, we experiment with the effectiveness of the algorithm under the following conditions.
Case 1: let
, the parameter in the linear part information of the probability distribution
is set to
This means the incomplete information probability distribution is satisfied
The parameters of modified NM method are given, reflection factor
, expansion factor
, contraction factor
, convergence criterion
.
Use MATLAB R2008a to achieve the above problems, Table 1 gives the decision variable x and the function value
, the calculation result retains four decimal places.
The actual running results show that the algorithm terminates after 61 times, the optimal solution is
, the optimal function value is
. The running time is 28.515197 seconds.
Case 2: To verify the effectiveness of the algorithm, we expand the value of N, let
, the parameter in the linear part information of the probability distribution
as follows,
,
,
, that is keep the row of B unchanged, extend the number of random variables to 100, this means
, use MATLAB R2008a to achieve the above problems again, the result is Table 2.
From Table 2, we can see the program stops at 89 times, the optimal solution is
, optimal function value is
. The running time is 382.307942 seconds. Comparing the results, we find that when the value of N is increased by 10 times, the running time increased by 13 times, this is normal, it need more times to calculate the recourse function. However, the number of iterations of the algorithm is only increased by 1/4, which shows that the algorithm in this paper is effective for solving stochastic programming problems.
Case 3: In order to investigate the sensitivity of the linear partial information constraint condition of the probability distribution to the algorithm, we increase
the constraints to observe the change of the optimal solution. Take
, and add the following two constraints to the probability distribution,
,
, following is the running result,
From Table 3 we can see, the program stops at 64 times, the optimal solution is
, optimal function value is
, the running time is 29.490318 seconds. this shows, on the one hand, as the information of the probability distribution changes from incomplete to complete information, the objective function values of the models (1)-(3) constructed tend to have better results. On the other hand, there is no significant increase in the number of iterations of the algorithm during the optimization process.
5. Conclusion
For the case that the probability distribution has incomplete information, this
paper establishes a stochastic quadratic programming model with incomplete probability distribution based on LPI theory, and designs an improved Nelder-Mead algorithm. The numerical examples are used to calculate the results. The results show that the established models and algorithms are reasonable and effective.