Solution of Stochastic Quadratic Programming with Imperfect Probability Distribution Using Nelder-Mead Simplex Method

Show more

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 $\left(\Omega ,\Sigma ,P\right)$ be a probability space, $h=h\left(\omega \right)\in {R}^{m}$ is a stochastic vector in this space, $C\in {R}^{n\times n}$ be symmetric positive semi-definite and $H\in {R}^{m\times m}$ be symmetric positive definite. We consider following problem:

$\begin{array}{l}\underset{x\in R{}^{n},y\in {R}^{m}}{\mathrm{min}}\text{}\frac{1}{2}{x}^{\text{T}}Cx+{c}^{\text{T}}x+f\left(y\right)\\ \text{s}\text{.t}\text{.}Ax\le b\\ \text{}Tx=y\end{array}$ (1)

where

$f\left(y\right)=\underset{P\in \pi}{\mathrm{max}}{E}_{P}\left(g\left(y,\omega \right)\right)$ (2)

$\begin{array}{l}g\left(y,{\omega}_{i}\right)=\underset{z\in {R}^{m}}{\mathrm{max}}\text{}-\frac{1}{2}{z}^{\text{T}}Hz+{z}^{\text{T}}\left(h\left({\omega}_{i}\right)-y\right)\\ \text{s}\text{.t}\text{.}Wz\le q\end{array}$ (3)

Here $x\in {R}^{n},y\in {R}^{m}$ are decision variables, $c\in {R}^{n},A\in R{}^{r\times n},b\in {R}^{r},T\in {R}^{m\times n}$ , $q\in {R}^{{m}_{1}},W\in {R}^{{m}_{1}\times m}$ are fixed matrices, $\omega \in {R}^{{m}_{2}}$ is a random vector with support $\Omega \subseteq {R}^{{m}_{2}}$ , $P={\left({p}_{1},{p}_{2},\cdots ,{p}_{N}\right)}^{\text{T}}$ is the probability distributions of limited sample set $\Omega =\left({\omega}_{1},\cdots ,{\omega}_{N}\right)$ , that is ${p}_{i}=prob\left(\left\{\omega ={\omega}_{i}\right\}\right)$ . Assumed that the probability distribution of random variable has the following linear partial information:

$\pi =\left\{P\in {R}^{N}|BP\le d,{\displaystyle \underset{i=1}{\overset{N}{\sum}}{p}_{i}=1;{p}_{i}}\ge 0,i=1,\cdots ,N\right\}$

Here $d\in {R}^{s},B\in {R}^{s\times N}$ are fixed matrices, assumed the set of probability distributions $\pi $ is a polyhedron. Thus the two-stage function can be written as

$f\left(y\right)=\underset{P\in \pi}{\mathrm{max}}{E}_{P}\left(g\left(y,\omega \right)\right)=\underset{P\in \pi}{\mathrm{max}}{\displaystyle \underset{i=1}{\overset{N}{\sum}}{p}_{i}g\left(y,{\omega}_{i}\right)}$ (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 $g\left(y,{\omega}_{i}\right)$ is the convex function about y ( [8] ), $f\left(y\right)$ 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 $\lambda $ , convex problems (1)-(3) can be written as an unconstrained problem:

$\theta \left(\mu \right)=\mathrm{min}\frac{1}{2}{x}^{\text{T}}Px+{c}^{\text{T}}x+f\left(y\right)+{\lambda}^{\text{T}}\left(Ax-b\right)$ (5)

$\mu ={\left(x,\lambda \right)}^{\text{T}}$ , let ${\mu}^{1},{\mu}^{2},\cdots ,{\mu}^{n+1}$ be the $n+1$ points of n-dimensional space of ${R}^{n}$ , which are not on the same plane. Let ${\mu}^{h},{\mu}^{s},{\mu}^{l}$ represent the points that have the highest, second highest, and lowest estimates of function values, ${\mu}^{c}$ is

the centroid of all vertices other than ${\mu}^{h}$ , ${\mu}^{c}=\frac{1}{n}{\displaystyle \underset{i\ne h}{\sum}{\mu}^{i}}.$

Reflection: since ${\mu}^{h}$ is the vertex with the higher value among the vertices, we can expect to find a lower value at the reflection of ${\mu}^{h}$ in the opposite face formed by all vertices ${\mu}^{i}$ except ${\mu}^{h}$ . Generate a new point ${\mu}^{r}$ , by reflecting ${\mu}^{h}$ through ${\mu}^{c}$ according to

${\mu}^{r}={\mu}^{c}+\alpha \left({\mu}^{c}-{\mu}^{h}\right)\text{with}\left(\alpha >0\right)$ .

If the reflected point is better than the second worst, but not better than the best, i.e. $\theta \left({\mu}^{l}\right)<\theta \left({\mu}^{r}\right)<\theta \left({\mu}^{s}\right)$ , then obtain a new simplex by replacing the worst point ${\mu}^{h}$ with the reflected point ${\mu}^{r}$ .

Expansion: if the reflected point is the best point so far, we can expect to find interesting values along the direction from ${\mu}^{c}$ to ${\mu}^{r}$ , that is if $\theta \left({\mu}^{r}\right)<\theta \left({\mu}^{l}\right)$ , then search direction is favorable, compute the expansion point using

${\mu}^{e}={\mu}^{c}+\beta \left({\mu}^{c}-{\mu}^{h}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}with\text{\hspace{0.17em}}\left(\beta >1\right)$ .

If the expanded point is better than the reflected point, that is $\theta \left({\mu}^{e}\right)<\theta \left({\mu}^{r}\right)$ , then replace ${\mu}^{h}$ by ${\mu}^{e}$ , otherwise, obtain a new simplex by replacing the worst point ${\mu}^{h}$ with the reflected point ${\mu}^{r}$ .

Contraction: here it is certain that $\theta \left({\mu}^{r}\right)>\theta \left({\mu}^{s}\right)$ , 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 $\theta \left({\mu}^{s}\right)<\theta \left({\mu}^{r}\right)<\theta \left({\mu}^{h}\right)$ , the contracted point is determined by ${\mu}^{p}={\mu}^{c}+\gamma \left({\mu}^{r}-{\mu}^{c}\right)$ with $0\le \gamma <1$ , if $\theta \left({\mu}^{p}\right)<\theta \left({\mu}^{h}\right)$ , the contraction is accepted, Replaced ${\mu}^{h}$ by ${\mu}^{p}$ .

2) If $\theta \left({\mu}^{r}\right)\ge \theta \left({\mu}^{h}\right)$ , the contracted point is determined by ${\mu}^{p}={\mu}^{c}+\gamma \left({\mu}^{h}-{\mu}^{c}\right)$ , if $\theta \left({\mu}^{p}\right)<\theta \left({\mu}^{h}\right)$ , the contraction is accepted. Replaced ${\mu}^{h}$ by ${\mu}^{p}$ .

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 ${\mu}^{l}$ with ${\mu}^{i}={\mu}^{l}+\delta \left({\mu}^{i}-{\mu}^{l}\right)$ . 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 $F\left({\mu}^{i}\right)=-\theta \left({\mu}^{i}\right)+M$ , where M is a fully large number, calculate the probability of obtaining ${\mu}^{i}$ by $F\left({\mu}^{i}\right)/{\displaystyle {\sum}_{i=1}^{n+1}F\left({\mu}^{i}\right)}$ . 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 $\delta \left({\mu}^{i}\right)$ of ${\mu}^{i}$ is defined as

$\delta \left({\mu}^{i}\right)=\left\{\mu :\Vert \mu -{\mu}^{i}\Vert \le \mathrm{min}\left\{\Vert {\mu}^{i}-{\mu}^{j}\Vert ,\forall j\ne i\right\}\right\}$

Algorithm termination condition: There are different criteria to determine the termination conditions of the NM algorithm in practice, in this paper, we use

$\left|\frac{\theta \left({\mu}^{h}\right)-\theta \left({\mu}^{l}\right)}{\theta \left({\mu}^{l}\right)}\right|\le \epsilon ,$

as our convergence criterion.

Parameters choice: The polyhedron transform in the NM algorithm mainly includes four parameters, $\alpha $ for reflection, $\beta $ for expansionand $\gamma $ for contraction, assumed they satisfy the following constraints:

$\alpha >0,\beta >1,\beta >\alpha ,0<\gamma <1.$

Algorithm

Choose $\alpha >0,\beta >1,\beta >\alpha ,0<\gamma <1$ , Convergence criterion $\epsilon >0$ .

Step 1 calculate the function values of n+1 points, rank all points and identify ${\mu}_{k}^{h},{\mu}_{k}^{s},{\mu}_{k}^{l}$ , find ${\mu}_{k}^{c}$ , the centroid of all vertices other than ${\mu}_{k}^{h}$ , generating a new point ${\mu}_{k}^{r}$ by reflecting point ${\mu}_{k}^{h}$ through ${\mu}_{k}^{c}$ according to reflection rule ${\mu}_{k}^{r}={\mu}_{k}^{c}+\alpha \left({\mu}_{k}^{c}-{\mu}_{k}^{h}\right)$ ;

Step 2 if $\theta \left({\mu}_{k}^{r}\right)<\theta \left({\mu}_{k}^{l}\right)$ , then the reflection point is expanded using the expansion rule ${\mu}_{k}^{e}={\mu}_{k}^{c}+\beta \left({\mu}_{k}^{c}-{\mu}_{k}^{h}\right)$ , if $\theta \left({\mu}_{k}^{e}\right)<\theta \left({\mu}_{k}^{r}\right)$ , then replace ${\mu}_{k}^{h}$ by ${\mu}_{k}^{e}$ , otherwise, let ${\mu}_{k}^{r}$ replaced ${\mu}_{k}^{h}$ , return to Step 4;

Step 3 if $\theta \left({\mu}_{k}^{l}\right)<\theta \left({\mu}_{k}^{r}\right)<\theta \left({\mu}_{k}^{s}\right)$ , then let ${\mu}_{k}^{r}$ replaced ${\mu}_{k}^{h}$ , 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 $\theta \left({\mu}_{k}^{r}\right)\ge \theta \left({\mu}_{k}^{s}\right)$ , then the simplex contracts.

(i) if $\theta \left({\mu}_{k}^{s}\right)<\theta \left({\mu}_{k}^{r}\right)<\theta \left({\mu}_{k}^{h}\right)$ , the contraction point is determined by calculate ${\mu}_{k}^{p}={\mu}_{k}^{c}+\gamma \left({\mu}_{k}^{r}-{\mu}_{k}^{c}\right)$ ,

(ii) if $\theta \left({\mu}_{k}^{r}\right)\ge \theta \left({\mu}_{k}^{h}\right)$ , the contraction point is determined by calculate ${\mu}_{k}^{p}={\mu}_{k}^{c}+\gamma \left({\mu}_{k}^{h}-{\mu}_{k}^{c}\right)$ .

In these case, if $\theta \left({\mu}_{k}^{p}\right)<\theta \left({\mu}_{k}^{r}\right)$ , the contraction accepted, If contraction is accepted, let ${\mu}_{k}^{p}$ replaced ${\mu}_{k}^{h}$ , 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 $X\in {R}^{3},H=I\in {R}^{3\times 3},T=-H$ , and

$C=\left(\begin{array}{ccc}2& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right),W=\left(\begin{array}{ccc}1& -2& 3\\ 0& 1& -1\\ 1& 0& 1\end{array}\right),q=\left(\begin{array}{c}6\\ 3\\ 4\end{array}\right)$

$c=\left(\begin{array}{c}-4\\ -3\\ -2\end{array}\right),A=\left(\begin{array}{ccc}3& 2& 1\\ 1& 2& 1\end{array}\right),b=(128)$

Assumed stochastic vector ${\omega}_{i}={\left({\omega}_{i}^{1},{\omega}_{i}^{2},{\omega}_{i}^{3}\right)}^{\text{T}}$ $\left(i=1,2,\cdots ,N\right)$ is three-dimensional vector, ${\omega}_{i}^{1},{\omega}_{i}^{2}$ and ${\omega}_{i}^{3}$ are independent of each other. We use MATLAB to randomly generate one hundred values for each other. Let $h\left({\omega}_{i}\right)={\omega}_{i}$ , we experiment with the effectiveness of the algorithm under the following conditions.

Case 1: let $N=10$ , the parameter in the linear part information of the probability distribution $\pi $ is set to

$\begin{array}{l}{D}_{8}={\left(1,\cdots ,1\right)}^{\text{T}},B=\left({E}_{8}{\text{,O}}_{\text{8}\times 1}\text{,}{D}_{8}\right),\\ d={\left(1/6,1/8,1/5,1/4,1/7,1/6,1/7,2/7\right)}^{\text{T}}\end{array}$

This means the incomplete information probability distribution is satisfied

${p}_{1}+{p}_{10}\le 1/6,{p}_{2}+{p}_{10}\le 1/8,\cdots ,{p}_{1}+{p}_{2}+\cdots +{p}_{10}=1$

The parameters of modified NM method are given, reflection factor $\alpha =1$ , expansion factor $\beta =2$ , contraction factor $\gamma =0.5$ , convergence criterion $\epsilon =1\times {10}^{-8}$ .

Use MATLAB R2008a to achieve the above problems, Table 1 gives the decision variable x and the function value $\theta $ , the calculation result retains four decimal places.

The actual running results show that the algorithm terminates after 61 times, the optimal solution is $x={\left(0.8634,0.4194,0.2721\right)}^{\text{T}}$ , the optimal function value is $\theta =4.4630$ . The running time is 28.515197 seconds.

Case 2: To verify the effectiveness of the algorithm, we expand the value of N, let $N=100$ , the parameter in the linear part information of the probability distribution $\pi $ as follows, ${D}_{8}={\left(1,\cdots ,1\right)}^{\text{T}}$ , $B=\left({E}_{8}{\text{,O}}_{\text{8}\times 91}\text{,}{D}_{8}\right)$ , $d={\left(1/6,1/8,1/5,1/4,1/7,1/6,1/7,2/7\right)}^{\text{T}}$ , that is keep the row of B unchanged, extend the number of random variables to 100, this means

${p}_{1}+{p}_{100}\le 1/6,{p}_{2}+{p}_{100}\le 1/8,\cdots ,{p}_{1}+{p}_{2}+\cdots +{p}_{100}=1$ , 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 $x={\left(0.8659,0.3530,0.2682\right)}^{\text{T}}$ , optimal function value is $\theta =4.8989$ . 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

Table 1. Iterative results.

the constraints to observe the change of the optimal solution. Take $N=10$ , and add the following two constraints to the probability distribution, ${p}_{2}+{p}_{3}\le 1/9$ , ${p}_{4}+{p}_{5}\le 1/10$ , following is the running result,

From Table 3 we can see, the program stops at 64 times, the optimal solution is $x={\left(0.8929,0.4113,0.3524\right)}^{\text{T}}$ , optimal function value is $\theta =3.7382$ , 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

Table 2. Iterative results.

Table 3. Iterative results.

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.

References

[1] Dantzig, G.B. and Ferguson, A.R. (1956) The Allocation of Air Craft Routes-An Example of Linear Programming under Uncertainty. Management Science, 3, 23-34.

[2] Birge, J.R. and Louveanx, F. (2003) Introduction to Stochastic Programming. Springer-Verlag, New York.

[3] Slyke, R. and Wets, R. (1969) L-Shaped Linear Programs with Applications to Optimal Control and Stochastic Programming. SIAM Journal on Applied Mathematics, 17, 638-663.

https://doi.org/10.1137/0117061

[4] Abaffy, J. and Allevi, E. (2004) A Modified L-shaped Method. Journal of Optimization Theory and Applications, 11, 255-270.

https://doi.org/10.1007/s10957-004-5148-y

[5] Facchinei, F. and Pang, J.S. (2003) Finite-Dimensional Variational Inequalities and Complementarity Problems. Springer-Verlag, New York.

[6] Qi, L. (1993) Convergence Analysis of Some Algorithms for Solving Non-Smooth Equations. Mathematics and Operations Research, 18, 227-244.

https://doi.org/10.1287/moor.18.1.227

[7] Qi, L. and Sun, J. (1993) A Non-Smooth Version of Newton’s Method. Mathematical Programming, 58, 353-367.

https://doi.org/10.1007/BF01581275

[8] Chen, X., Qi, L, and Womersley, R.S. (1995) Newton’s Method for Quadratic Stochastic Programs with Recourse. Journal of Computational and Applied Mathematics, 60, 29-46.

https://doi.org/10.1016/0377-0427(94)00082-C

[9] Abdelaziz, F.B. and Masri, H. (2005) Stochastic Programming with Fuzzy Linear Partial Information on Probability Distribution. European Journal of Operational Research, 162, 619-629.

https://doi.org/10.1016/j.ejor.2003.10.049

[10] Abdelaziz, F.B. and Masri, H. (2009) Multistage Stochastic Programming with Fuzzy Probability Distribution. Fuzzy Sets and Systems, 160, 3239-3249.

https://doi.org/10.1016/j.fss.2008.10.010

[11] Zhang, Y. and Ma, X. (2016) A Class Recourse Stochastic Programs Algorithm with MaxEM in Evaluation. Operations Research Transactions, 20, 52-60.

[12] Kofler, E. (2001) Linear Partial Information with Applications. Fuzzy Sets and Systems, 118, 167-177.

https://doi.org/10.1016/S0165-0114(99)00088-3

[13] Boyd, S. and Vandenberghe, L. (2004) Convex Optimization. Cambridge University Press, Cambridge.

https://doi.org/10.1017/CBO9780511804441

[14] Nelder, J.A. and Mead, R. (1965) A Simplex Method for Function Minimization. The Computer Journal, 7, 308-313.

https://doi.org/10.1093/comjnl/7.4.308

[15] Chang, K.H. (2012) Stochastic Nelder-Mead Simplex Method—A New Globally Convergent Direct Search Method for Simulation Optimization. European Journal of Operational Research, 220, 684-694.

https://doi.org/10.1016/j.ejor.2012.02.028