The Backpropagation Artificial Neural Network Based on Elite Particle Swam Optimization Algorithm for Stochastic Linear Bilevel Programming Problem

Abstract

For a class of stochastic linear bilevel programming problem, we firstly transform it into a deterministic linear bilevel covariance programming problem. Then, the deterministic bilevel covariance programming problem is solved by backpropagation artificial neural network based on elite particle swam optimization algorithm (BPANN-PSO). Finally, we perform the simulation experiments and the results show that the computational efficiency of the proposed algorithm has a potential upside compared with the classical algorithm.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The bilevel programming problem (BLP) is a nested optimizations problem in which the feasible region of the upper level problem is determined implicitly by the solution set of the lower level problem. As an optimization tool, the BLP has been widely used in variety practical problems, for example, in homeland security [1–3], model production processes [4], the optimal tax policies formulation [5–7], the strategic for deregulating markets [8], and the optimization of retail channel structures [9]. In addition, the optimization theory of the BLP has been integrated in many other disciplines, such as in management [10, 11], chemical engineering [12, 13], structural optimization [14, 15], optimal control problems [16, 17], facility location [10, 18, 19], and transportation [20–22].

Therefore, many researchers are devoted to develop the algorithms for BLP and propose many efficient algorithms. Traditional methods commonly used to handle BLP include Karus-Kuhn-Tucker approach [23–26], Branch-and-bound method [27], and penalty function approach [28–31]. Despite a significant progress made in traditional optimization towards solving BLPPs, the properties such as differentiation and continuity are necessary for these algorithms.

Due to the limitation of the traditional algorithms, the heuristics such as evolutionary algorithms are recognized as potent tools for solving BLPPs. Mathieu et al. [32] firstly developed the genetic algorithm (GA) for bilevel linear programming problem. Motivated by the same reason, other kinds of GAs for BLPPs were also presented in [33–36]. Owing to its high speed of convergence and relative simplicity, the particle swam optimization (PSO) algorithm has been employed for solving BLP problems recently [37–42].

However, it is worth noting that the papers mentioned above only focus on deterministic bilevel programming problem and the stochastic bilevel programming has seldom been studied so far. In 1999, Patriksson and Wynter [43] firstly proposed the stochastic mathematical programs with equilibrium constraints and introduced a framework for hierarchical decision-making problem under uncertainty. However, they did not give a numerical experiment. Gao, Liu, and Gen [44] presented a hybrid intelligent algorithm for a decentralized multilevel decision-making problem in stochastic environment in 2004. For the contracting arrangements of the long-term contracts and the spot markets transactions under uncertain electricity spot market, Wan et al. [45] proposed a stochastic bilevel programming model for the optimal bidding strategies between power seller and buyer. They solved the model by Monte Carlo approximation method. It is worth noting that the decision variable in both levels are one-dimensional variable. Soon after, Wan, Fan, and Wang [46] proposed an interactive fuzzy decision-making method for the model in [45]. In 2013, He and Feng [47] presented an approximation algorithm for the compensated stochastic bilevel programming problem. In addition, the stability analysis and convergence analysis for bilevel stochastic programming problem can be seen in [48–50]. Obviously, they only researched the simple stochastic bilevel model and few of them have studied the numerical performance of the algorithm.

In this paper, we consider the general stochastic linear bilevel programming problem in which the coefficient of objective functions and the coefficient of constraint functions are random variables. For the problem, we firstly transformed it into a deterministic linear bilevel covariance programming problem with expected constraints, and then the deterministic bilevel covariance programming model is solved by the BPANN-PSO algorithm. Finally, we perform the simulation experiments and the results suggest that the variance obtained by our algorithm is better than the results in reference when the means of the upper objective function value is same. Furthermore, the computational efficiency of our algorithm performs better with the dimension increasing.

The rest of this paper is organized as follows. Section 2 introduces the definitions and properties of stochastic linear bilevel programming problem. Section 3 proposes the BPANN-PSO algorithm for stochastic linear bilevel programming problem. We use three test problems from the reference to measure and evaluate the proposed algorithm in Section 4, while the conclusion is reached in Section 5.

2. Stochastic Linear Bilevel Programming

Let $x\in R^{n_1}$, $y\in R^{n_2}$, $F:R^{n_1}\times R^{n_2}\to R$, $f:{\mathrm{R}}^{n_1}\times R^{n_2}⟶R,\hspace{0.17em}\hspace{0.17em}A\in m\times n_{\mathrm{1}}$, $B\in m\times n_{\mathrm{2}}$, we consider stochastic linear bilevel programming problem as follows: $$\begin{array}{c}\text{(1)}& \underset{x,y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}F\left(x,y\right)=c_{\mathrm{1}}^{T}x+d_{\mathrm{2}}^{T}y\underset{y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}f\left(x,y\right)=c_{\mathrm{2}}^{T}x+d_{\mathrm{2}}^{T}y\\ \text{s.t.}\hspace{1em}Ax+By\le b\\ \hspace{1em}x\ge \mathrm{0},\hspace{0.17em}\hspace{0.17em}y\ge \mathrm{0}.\end{array}$$where $F(x,y)$ and $f(x,y)$ are the upper level and the lower level objective functions, respectively. The coefficient $c_{\mathrm{1}},c_{\mathrm{2}}\in R^{n_{\mathrm{1}}}$, $d_{\mathrm{1}},d_{\mathrm{2}}\in R^{n_{\mathrm{2}}}$ are the random variable.

Let ${\alpha }_i$ be a probability of the extent to which the $ith$ constraint violation is admitted; $P\left(·\right)$ means a probability measure. The $ith$ constraint of problem (1) is interpreted as follows:$$\begin{array}{}\text{(2)}& P\left(A_{i}x+B_{i}y\le b_i\right)\ge a_i\hspace{1em}\left(i=\mathrm{1,2},\dots ,m\right)\end{array}$$

Inequality (2) means that the $ith$ constraint may be violated, but at most $\mathrm{1}-{\alpha }_i$ proportion of the time. Let $F_i(\tau )$ be the distribution function of random variables $b_i$, then, inequality (2) is presented as$$\begin{array}{}\text{(3)}& P\left(A_{i}x+B_{i}y\le b_i\right)=\mathrm{1}-F\left(A_{i}x+B_{i}y\right)\end{array}$$According to (2) and (3), we can obtain inequality (4),$$\begin{array}{}\text{(4)}& F\left(A_{i}x+B_{i}y\right)\le \mathrm{1}-{\alpha }_i\end{array}$$

Let $\tau =F^{-\mathrm{1}}(\mathrm{1}-{\alpha }_i),\hspace{0.17em}\hspace{0.17em}{K}_{\mathrm{1}-{\alpha }_i}=\max \tau$. Then, from the monotonicity of the distribution function $F_i(\tau )$, inequality (2) is rewritten as$$\begin{array}{}\text{(5)}& A_{i}x+B_{i}y\le K_{\mathrm{1}-\alpha }\hspace{1em}\left(i=\mathrm{1,2},\dots ,m\right)\end{array}$$That is, $$\begin{array}{}\text{(6)}& Ax+By\le K_{\mathrm{1}-\alpha }\end{array}$$where $K_{\mathrm{1}-\alpha }=(K_{\mathrm{1}-{\alpha }_{\mathrm{1}}},K_{\mathrm{1}-{\alpha }_{\mathrm{2}}},\dots ,K_{\mathrm{1}-{\alpha }_m})$. Based on (6), problem (1) is rewritten as the following problem with deterministic constraints:$$\begin{array}{c}\text{(7)}& \underset{x,y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}F\left(x,y\right)=c_{\mathrm{1}}^{T}x+d_{\mathrm{1}}^{T}y\underset{y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}f\left(x,y\right)=c_{\mathrm{2}}^{T}x+d_{\mathrm{2}}^{T}y\\ \text{s.t.}\hspace{1em}Ax+By\le K_{\mathrm{1}-\alpha }\\ \hspace{1em}x\ge \mathrm{0},\hspace{0.17em}\hspace{0.17em}y\ge \mathrm{0}.\end{array}$$

To deal with the objective functions with random variables in both levels, the minimum covariance model [51] is applied in this paper. Then, problem (7) is rewritten as follows:$$\begin{array}{c}\text{(8)}& \underset{x,y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}\mathrm{var}\left(F\left(x,y\right)\right)=\left(x^T,y^T\right)V_{\mathrm{1}}\left(\begin{array}{c}x\\ y\end{array}\right)\underset{y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}\mathrm{var}\left(f\left(x,y\right)\right)=\left(x^T,y^T\right)V_{\mathrm{2}}\left(\begin{array}{c}x\\ y\end{array}\right)\\ \text{s.t.}\hspace{1em}Ax+By\le K_{\mathrm{1}-\alpha }\\ \hspace{1em}x\ge \mathrm{0},\hspace{0.17em}\hspace{0.17em}y\ge \mathrm{0}.\end{array}$$where $\mathrm{var}(·)$ denotes the variance of objective function; $V_{\mathrm{1}}$ and $V_{\mathrm{2}}$ are the covariance matrices of $(c_{\mathrm{1}},d_{\mathrm{1}})$ and $(c_{\mathrm{2}},d_{\mathrm{2}})$, respectively. In this paper, we assume that $V_{\mathrm{1}}$ and $V_{\mathrm{2}}$ are positive-definite.

Though problem (8) can guarantee uniform distribution of the objective function values of both levels with a small drop, the decision makers in both levels often have their own expectation. Let ${\beta }_{\mathrm{1}}$ and ${\beta }_{\mathrm{2}}$ denote the means of the objective function values for the leader and the follower, respectively. Then, problem (8) is rewritten as the follows:$$\begin{array}{c}\text{(9)}& \underset{x,y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}\mathrm{var}\left(F\left(\overline{x},y\right)\right)=\left({\overline{x}}^T,y^T\right)V_{\mathrm{1}}\left(\begin{array}{c}\overline{x}\\ y\end{array}\right)\underset{y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}\mathrm{var}\left(f\left(\overline{x},y\right)\right)=\left({\overline{x}}^T,y^T\right)V_{\mathrm{2}}\left(\begin{array}{c}\overline{x}\\ y\end{array}\right)\\ \text{s.t.}\hspace{1em}A\overline{x}+By\le K_{\mathrm{1}-\alpha }\\ \hspace{1em}{m}_{\mathrm{1}}^c\overline{x}+m_{\mathrm{1}}^{d}y\le {\beta }_{\mathrm{1}}\\ \hspace{1em}{m}_{\mathrm{2}}^c\overline{x}+m_{\mathrm{2}}^{d}y\le {\beta }_{\mathrm{2}}\\ \hspace{1em}x\ge \mathrm{0},\hspace{0.17em}\hspace{0.17em}y\ge \mathrm{0}.\end{array}$$where $m_{\mathrm{1}}^c$ and $m_{\mathrm{1}}^d$ are the means of $c_{\mathrm{1}}$ and $d_{\mathrm{1}}$, respectively. $m_{\mathrm{2}}^c$ and $m_{\mathrm{2}}^d$ are the means of $c_{\mathrm{2}}$ and $d_{\mathrm{2}}$, respectively. For the fixed upper level decision variable $\stackrel{-}{x}$, the rational reaction set $R\left(x\right)$ can be solved by problem (10).$$\begin{array}{c}\text{(10)}& \underset{y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}\mathrm{var}\left(f\left(\overline{x},y\right)\right)=\left({\overline{x}}^T,y^T\right)V_{\mathrm{2}}\left(\begin{array}{c}\overline{x}\\ y\end{array}\right)\text{s.t.}\hspace{1em}A\overline{x}+By\le K_{1-\alpha }\\ \hspace{1em}{m}_{\mathrm{1}}^c\overline{x}+m_{\mathrm{1}}^{d}y\le {\beta }_{\mathrm{1}}\\ \hspace{1em}{m}_{\mathrm{2}}^c\overline{x}+m_{\mathrm{2}}^{d}y\le {\beta }_{\mathrm{2}}\\ \hspace{1em}y\ge 0.\end{array}$$Let $A_1^{\prime }={\left(A,m_1^c,m_2^c\right)}^T$, $A_2^{\prime }={\left(B,m_1^d,m_2^d\right)}^T$, $b^{\prime }={\left(K_{1-\alpha },{\beta }_1,{\beta }_2\right)}^T$. Then, problem (9) is formulated as$$\begin{array}{c}\text{(11)}& \underset{x,y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}\mathrm{var}\left(F\left(\overline{x},y\right)\right)=\left({\overline{x}}^T,y^T\right)V_{\mathrm{1}}\left(\begin{array}{c}\overline{x}\\ y\end{array}\right)\underset{y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}\mathrm{var}\left(f\left(\overline{x},y\right)\right)=\left({\overline{x}}^T,y^T\right)V_{\mathrm{2}}\left(\begin{array}{c}\overline{x}\\ y\end{array}\right)\\ \text{s.t.}\hspace{1em}{A}_1^{\prime }\overline{x}+A_2^{\prime }y\le b^{\prime }\end{array}$$For problem (11), the basic notions of BLP are recalled as follows:

(a)

Constraint region of the BLP:$$\begin{array}{}\text{(12)}& S=\left\{\left(x,y\right)\mid A_{\mathrm{1}}^{\prime }x+A_{\mathrm{2}}^{\prime }y\le b\right\}\end{array}$$

(b)

Feasible set for the lower level problem for each fixed $x$:$$\begin{array}{}\text{(13)}& S\left(x\right)=\left\{\left(x,y\right)\mid A_{\mathrm{1}}^{\prime }x+A_{\mathrm{2}}^{\prime }y\le b\right\}\end{array}$$

(c)

Projection of onto the upper level maker’s decision space:$$\begin{array}{}\text{(14)}& S\left(X\right)=\left\{x\mid \exists y,A_{\mathrm{1}}^{\prime }x+A_{\mathrm{2}}^{\prime }y\le b\right\}\end{array}$$

(d)

The lower level maker’s rational reaction set for each fixed $x\in S(X)$:$$\begin{array}{}\text{(15)}& P\left(x\right)=\left\{y:y\in \arg \min \left[f\left(x,y\right):y\in S\left(x\right)\right]\right\}\end{array}$$

(e)

Inducible region: $$\begin{array}{}\text{(16)}& IR=\left\{\left(x,y\right):\left(x,y\right)\in S,\hspace{0.17em}\hspace{0.17em}y\in P\left(x\right)\right\}\end{array}$$

Definition 1.

A point $(x,y)$ is feasible if $(x,y)\in IR$.

Definition 2.

A feasible point $(x^{\ast },y^{\ast })$ is an optimal solution if $(x^{\ast },y^{\ast })\in IR$ and $F(x^{\ast },y^{\ast })\le F(\overline{x},\overline{y})$, $\forall (\overline{x},\overline{y})\in IR$.

Definition 3.

If $(x^o,y^o)$ is the optimistic solution for problem (11), $(x^o,y^o)$ is given by ${\min }_{x,y}\{F(x,y)\mid y\in P(x),A_{\mathrm{1}}^{\prime }x+A_{\mathrm{2}}^{\prime }y\le b\}$.

3. The Algorithm

For problem (11), it is noted that a solution $(x^{\ast },y^{\ast })$ is feasible for the upper level problem if and only if $y^{\ast }$ is an optimal solution for the lower level problem with $x=x^{\ast }$. In practice, we often make the approximate optimal solutions of the lower level problem as the optimal response feedback to the upper level problem, and this point of view is accepted usually. Based on this fact, the BPANN-PSO algorithm may have a great potential for solving problem (11). In the following, the BPANN-PSO is presented for solving problem (11).

3.1. The Structure of BPANN

In this paper, the BPANN includes three layers: the input layer, the hidden layer, and the output layer. The number of nodes per layer is $n_{\mathrm{1}}$, $m$, and $n_{\mathrm{2}}$ ($m=\sqrt{n_{\mathrm{1}}+n_{\mathrm{2}}}$). The connectivity weight from the $i-th$ node of the input to the $j-th$ node of the hidden layer is ${\omega }_{ij}$. $w_{jk}$ denotes the connectivity weight from the $j-th$ node of the hidden layer to the $k-th$ node of the output layer. Let the Sigmoid function $f(\mathrm{x})=\mathrm{1}/\left(\mathrm{1}+e^{-x}\right)$. Figure 1 shows the diagram of the BPANN.

[figure omitted; refer to PDF]

Let $x_{ni}$ represent the $i-th$ input value of the $n-th$ sample. The output of the hidden layer is given by $$\begin{array}{}\text{(17)}& s_{nj}=f\left({\displaystyle \sum }_{i=\mathrm{0}}^{n_{\mathrm{1}}}{w}_{ij}{x}_{ni}\right),\hspace{1em}j=\mathrm{1,2},\dots ,m.\end{array}$$

For the output layer, the output is given by $$\begin{array}{}\text{(18)}& y_{nl}=f\left({\displaystyle \sum }_{j=\mathrm{0}}^m{w}_{jk}{s}_{nj}\right),\hspace{1em}l=\mathrm{1,2},\dots ,n_{\mathrm{2}}.\end{array}$$

$y_{nk}^{optimal}$ denotes the $k-th$ node’s expected output of the output layer. The squared differences between the computed and desired outputs are given by (19) and the error sum of all $N$ training pairs is given by (20):$$\begin{array}{}\text{(19)}& E_n=\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum }_{l=\mathrm{1}}^{n_{\mathrm{2}}}\left(t_{nl}-y_{nl}\right)\text{(20)}& E={\displaystyle \sum }_{n=\mathrm{1}}^N{E}_n=\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum }_{n=\mathrm{1}}^N\hspace{0.17em}{\displaystyle \sum }_{l=\mathrm{1}}^{n_{\mathrm{2}}}{\left(t_{nl}-y_{nl}\right)}^{\mathrm{2}}\end{array}$$

3.2. The BPANN-PSO Algorithm

In this algorithm, the position of the particle represents the connectivity weight of the BPANN and the particles in the elite set $A_t$ are the batch training samples. Based on this idea, the BPANN-PSO algorithm can be designed as follows.

Suppose $W_{ij}=\{w_{ij}\}=\{w_{\mathrm{01}},w_{\mathrm{11}},w_{\mathrm{21}},\dots ,w_{n_{\mathrm{1}}\mathrm{1}};w_{\mathrm{02}},w_{\mathrm{12}},w_{\mathrm{22}},\dots ,w_{n_{\mathrm{1}}\mathrm{2}};\dots ;w_{\mathrm{0}m},w_{\mathrm{1}m}$, $w_{\mathrm{2}m},\dots ,w_{n_{\mathrm{1}}m}\}$, $W_{jk}=\{w_{jk}\}=\{w_{\mathrm{11}},w_{\mathrm{21}},\dots ,w_{m\mathrm{1}};w_{\mathrm{12}},w_{\mathrm{22}},\dots ,w_{m\mathrm{2}};\dots ;w_{\mathrm{1}{n}_{\mathrm{2}}},w_{\mathrm{2}{n}_{\mathrm{2}}},\dots$, $w_{m{n}_{\mathrm{2}}}\}$. The $l-th$ particle in $A_t$ can be presented as $A^l=\{X^l,Y^{l,optimal}\}=\{x_{\mathrm{1}}^l,x_{\mathrm{2}}^l,\dots ,x_{n_{\mathrm{1}}}^l;y_{\mathrm{1}}^{l,optimal},y_{\mathrm{2}}^{l,optimal},\dots ,y_{n_{\mathrm{2}}}^{l,optimal}\}\hspace{0.17em}\hspace{0.17em}(l=\mathrm{1,2},\dots ,N(A_t))$.

Algorithm 4.

 

Step 1. Initialize the population $P$ with $N$ particles and each particle can be presented as $P_p=\left\{W_{ij}^p;W_{jk}^p\right\}\hspace{0.17em}\hspace{0.17em}(p=\mathrm{1,2},\dots ,N)$. Let the $N(A_t)$ be the number of the particles in the elite set $A_t$.

Step 2. Train the BPANN using the PSO.

Step 2.1. Calculate its fitness value according to (20).

Step 2.2. Determine its personal best particle and global best particle.

Step 2.3. Update particle’s position and velocity.

Step 3 (stopping criterion). If $gbest(p^t)\le \epsilon =\mathrm{1}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}$, stop. Otherwise, go to step 2.

In step 2.3, the inertia weight $\omega =\mathrm{0.5}+r/\mathrm{2.0},\hspace{0.17em}\hspace{0.17em}r=\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}(\mathrm{0,1})$, the acceleration constants $c_{\mathrm{1}}=c_{\mathrm{2}}=\mathrm{1.49618}$, and $r_{\mathrm{1}},r_{\mathrm{2}}\in \mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{m}(\mathrm{0,1})$.

3.3. The Algorithm for BLP

We proposed the BPANN-PSO for problem (11) and we update the upper level decision variables using the method in [52]. However, the way to solve the lower level problems is completely different. In reference [52], for each updated upper level decision variable, we need to reexecute the lower level optimization problem by cooperative coevolution PSO (CCPSO). In other words, the historical optimal solutions do not contribute to the current optimal solution. In our algorithm, for the updated upper level decision variable, we can predict the corresponding lower level optimal solution directly by BPANN. From the above analysis, we can see that the computational efficiency of our proposed algorithm performs well than the method in the reference [52]. Furthermore, we also give the basic working framework for these two algorithms (see Figure 2)

[figure omitted; refer to PDF]

Algorithm 5.

 

Step 1. Initialization scheme. Initialize a random population ($N_u$) of the upper level variable $X$.

Step 2. For the fixed $X^i$ ($i=\mathrm{1,2},\dots ,N_u$), initialize a random population ($N_l$) of the lower level variable $Y$. Then, we perform a lower level optimization procedure to determine the corresponding lower level optimal solution $Y_i^{optimal}$ using PSO.

Step 3. Evaluate the fitness value of the complete solutions $Z_i=(X_i,Y_i^{optimal})$ ($i=\mathrm{1,2},\dots ,N_u$) by the upper level function $F$ and constraints conditions. Copy all the members in an elite set $A_t$.

Step 4. Update the upper level decision variables using the simulated binary crossover operator (SBX) and the polynomial variation method (PM).

Step 5. For the new upper level variable $\overline{X^i}$. Predict the lower level optimal solutions using Algorithm 4.

Step 6. Update the elite set $A_t$. Copy all the offspring from the previous step to the elite set $A_t$ and choose the best $\left|A_t\right|$ members as the new members of $A_t$.

Step 7. Perform a termination check. If the termination check is false, go to step 4.

In the Algorithm 5, the PSO parameters are set as follows: the inertia weight $w=\mathrm{0.7298}$, acceleration coefficients $c_{\mathrm{1}}=c_{\mathrm{2}}=\mathrm{1.49618}$, and $r_{\mathrm{1}},r_{\mathrm{2}}\in \mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{m}(\mathrm{0,1})$, $N=\mathrm{20}$. At step 7, the termination criteria based on variance is described as follows. For each upper level variable $i$, the variance at the $j-th$ generation is $v_j^i$. Then, the stop metrics $\alpha$ is computed as$$\begin{array}{}\text{(21)}& {\alpha }_j={{\displaystyle \sum }}_{i=\mathrm{1}}^{N_u}\frac{v_j^i}{v_{\mathrm{0}}^i}\end{array}$$

In this paper, the value of ${\alpha }_{stop}$ is set as 10⁻⁴ for the upper level.

3.4. Algorithm Complexity Analysis

By the algorithm in [51], for the vector $u_i$ with dimension $(m+n_{\mathrm{2}}+\mathrm{2})$, each component of the vector has two states, zero or nonzero, therefore, the vector has a total of ${\mathrm{2}}^{m+n_{\mathrm{2}}+\mathrm{2}}$ forms. We examine the case where the calculation is the most expensive here, that is, until the last combination satisfies the complementary relaxation condition, suppose the number of vertices is $K$ and let the time be $t_{\mathrm{1}}$ for computing the $u_i{g}_i(x^k,y^k)$, and the total time required to complete a check is as follows: $T_{\mathrm{1}}=K\times {\mathrm{2}}^{m+n_{\mathrm{2}}+\mathrm{2}}\times t_{\mathrm{1}}$.

In Algorithm 4 of our paper, the number of nodes of the input layer, the hidden layer, and the output layer are $n_{\mathrm{1}}$, $m$, and $n_{\mathrm{2}}$ ($m=\sqrt{n_{\mathrm{1}}+n_{\mathrm{2}}}$), respectively. The batch training sample is $N(A_t)$ and the time for computing multiplication is $t_{\mathrm{2}}$; for one training sample, the computing time is as follows: $T_{\mathrm{2}}=(n_{\mathrm{1}}\times m+m\times n_{\mathrm{2}})\times t_{\mathrm{2}}=(n_{\mathrm{1}}\times \sqrt{n_{\mathrm{1}}+n_{\mathrm{2}}}+n_{\mathrm{2}}\times \sqrt{n_{\mathrm{1}}+n_{\mathrm{2}}})\times t_{\mathrm{2}}=(n_{\mathrm{1}}+n_{\mathrm{2}}{)}^{\mathrm{3}/\mathrm{2}}\times t_{\mathrm{2}}$. The total computing time of all samples is as follows: $T_{\mathrm{3}}=N(A_t)\times (n_{\mathrm{1}}+n_{\mathrm{2}}{)}^{\mathrm{3}/\mathrm{2}}\times t_{\mathrm{2}}$. In Algorithm 5, suppose the number of iterations is $G$, the computing time required to solve problem (11) is as follows: $T_{\mathrm{4}}=G\times N(A_t)\times (n_{\mathrm{1}}+n_{\mathrm{2}}{)}^{\mathrm{3}/\mathrm{2}}\times t_{\mathrm{2}}$. $t_{\mathrm{1}}$ and $t_{\mathrm{2}}$ are almost the same when they have the same operating platform. Comparing $T_{\mathrm{1}}$ and $T_{\mathrm{4}}$, our algorithm has obvious advantages when the dimension of the decision variable increases.

4. Numerical Experiment

In this section, we test the algorithm using three examples. All results presented in this paper have been obtained on a personal computer (CPU: AMD Phenon(tm)IIX6 1055T 2.80 GHz; RAM: 3.25 GB) using a C# implementation of the proposed algorithm.

Example 1.

Let $x\in R^{\mathrm{1}},y\in R^{\mathrm{1}}$; we consider the problem from [51]$$\begin{array}{c}\text{(22)}& \underset{x,y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}F\left(x,y\right)=c_{\mathrm{1}}x+d_{\mathrm{1}}y\underset{y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}f\left(x,y\right)=c_{\mathrm{2}}x+d_{\mathrm{2}}y\\ \text{s.t.}\hspace{1em}-x+\mathrm{3}y\le b_{\mathrm{1}},\\ \hspace{1em}\mathrm{10}x-y\le b_{\mathrm{2}}\\ \hspace{1em}\mathrm{3}x+y\ge b_{\mathrm{3}},\\ \hspace{1em}x+\mathrm{2}y\ge b_{\mathrm{4}}\\ \hspace{1em}\mathrm{3}x+\mathrm{2}y\ge b_{\mathrm{5}},\\ \hspace{1em}x\ge \mathrm{0},\hspace{0.17em}\hspace{0.17em}y\ge \mathrm{0}.\end{array}$$where $m_{\mathrm{1}}^c=-\mathrm{2.0},\hspace{0.17em}\hspace{0.17em}{m}_{\mathrm{1}}^d=-\mathrm{3.0},\hspace{0.17em}\hspace{0.17em}{m}_{\mathrm{2}}^c=\mathrm{2.0},\hspace{0.17em}\hspace{0.17em}{m}_{\mathrm{2}}^d=\mathrm{1.0}$. $V_{\mathrm{1}}=\left(\begin{array}{cc}\mathrm{2}& \mathrm{1}\\ \mathrm{1}& \mathrm{3}\end{array}\right),\hspace{0.17em}\hspace{0.17em}{V}_{\mathrm{2}}=\left(\begin{array}{cc}\mathrm{1}& -\mathrm{1}\\ -\mathrm{1}& \mathrm{6}\end{array}\right)$. All the right-hand side values $b=(b_{\mathrm{1}},b_{\mathrm{2}},\dots ,b_{\mathrm{5}})$ are normal random variables. The mean, variance, and probability of $b$ are shown in Table 1.

According to the above conditions, problem (22) can be rewritten as follows:$$\begin{array}{c}\text{(23)}& \underset{x,y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}F\left(x,y\right)=\left(x,y\right)\left(\begin{array}{cc}\mathrm{2}& \mathrm{1}\\ \mathrm{1}& \mathrm{3}\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)\underset{y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}f\left(x,y\right)=\left(x,y\right)\left(\begin{array}{cc}\mathrm{1}& -\mathrm{1}\\ -\mathrm{1}& \mathrm{6}\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)\\ \text{s.t.}\hspace{1em}\begin{array}{c}-x\end{array}+\mathrm{3}y\le \mathrm{47},\\ \hspace{1em}\mathrm{10}x-y\le \mathrm{110}\\ \hspace{1em}\mathrm{3}x+y\ge -\mathrm{19},\\ \hspace{1em}x+\mathrm{2}y\ge -\mathrm{15}\\ \hspace{1em}\mathrm{3}x+\mathrm{2}y\ge -\mathrm{29},\\ \hspace{1em}\begin{array}{c}-\mathrm{2}x\end{array}-\mathrm{3}y\le -\mathrm{31}\\ \hspace{1em}\mathrm{2}x+y\le \mathrm{33},\\ \hspace{1em}x\ge \mathrm{0},\hspace{0.17em}\hspace{0.17em}y\ge \mathrm{0}.\end{array}$$

We solve problem (23) with constrain of means of the objective function values and we also solved problem (23) without this constrain. From Table 2, we can see that the variance obtained by our algorithm is better than the results in [51] when the means of the upper objective function value is the same. Moreover, the computational efficiency of our algorithm performs better.

Table 1

The mean, variance, and acceptable probability of $b$.

	$$b_1$$	$$b_2$$	$$b_3$$	$$b_4$$	$b_5$
mean	50.11	113.15	15.16	13.16	25.63
variance	9.0	36	9.0	4.0	16
acceptable probability	0.85	0.7	0.9	0.7	0.8

Table 2

The comparison results of Example 1 from 11 runs.

Algorithm	means constrain	optimal solutions	means		variance		CPU time
			UL	LL	UL	LL	(s)
BPANN-PSO	Y	(5.15, 6.90)	-31.00	17.20	266.95	241.11	0.0335
Method in [51]	Y	(5.17, 6.89)	-31.00	17.23	267.12	240.32	0.0879
BPANN-PSO	N	(6.99, 4.00)	-25.98	17.99	201.64	88.94	0.0173
Method in [51]	N	(7.00, 4.00)	-26.00	18.00	202.00	89.00	0.0208

Example 2.

Let $x\in R^{\mathrm{5}},y\in R^{\mathrm{3}}$, $c_{\mathrm{1}},c_{\mathrm{2}}\in R^{\mathrm{5}\times \mathrm{1}}$, $d_{\mathrm{1}},d_{\mathrm{2}}\in R^{\mathrm{3}\times \mathrm{1}}\hspace{0.17em}\hspace{0.17em}A\in R^{\mathrm{3}\times \mathrm{5}}$, $B\in R^{\mathrm{3}\times \mathrm{3}}$, $b\in R^{\mathrm{3}\times \mathrm{1}}$; we consider problem from [51]$$\begin{array}{c}\text{(24)}& \underset{x,y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}F\left(x,y\right)=c_{\mathrm{1}}^{T}x+d_{\mathrm{1}}^{T}y\underset{y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}f\left(x,y\right)=c_{\mathrm{2}}^{T}x+d_{\mathrm{2}}^{T}y\\ \text{s.t.}\hspace{1em}Ax+By\le b,\\ \hspace{1em}x\ge \mathrm{0},\hspace{0.17em}\hspace{0.17em}y\ge \mathrm{0}\end{array}$$where $b_{\mathrm{1}}~N(\mathrm{220},{\mathrm{4}}^{\mathrm{2}})$, $b_{\mathrm{2}}~N(\mathrm{145},{\mathrm{3}}^{\mathrm{2}})$, $b_{\mathrm{3}}~N(-\mathrm{18},{\mathrm{5}}^{\mathrm{2}})$, ${\alpha }_{\mathrm{1}}=\mathrm{0.85}$, ${\alpha }_{\mathrm{2}}=\mathrm{0.95}$, ${\alpha }_{\mathrm{3}}=\mathrm{0.90}$, $K_{\mathrm{1}-{\alpha }_{\mathrm{1}}}=\mathrm{215.85}$, $K_{\mathrm{1}-{\alpha }_{\mathrm{2}}}=\mathrm{140.07}$, $K_{\mathrm{1}-{\alpha }_{\mathrm{3}}}=\mathrm{24.41}$, ${\beta }_{\mathrm{1}}=\mathrm{60.0}$, and ${\beta }_{\mathrm{2}}=\mathrm{35.0}$. The expectation of random variable $c_{\mathrm{1}},c_{\mathrm{2}},d_{\mathrm{1}},d_{\mathrm{2}}$ is shown in Table 3. The coefficient matrices of constraints are as follows:$$\begin{array}{c}\text{(25)}& A=\left(\begin{array}{ccccc}\mathrm{7}& \mathrm{2}& \mathrm{6}& \mathrm{9}& \mathrm{11}\\ \mathrm{5}& \mathrm{6}& -\mathrm{4}& \mathrm{3}& \mathrm{3}\\ -\mathrm{4}& -\mathrm{7}& -\mathrm{2}& -\mathrm{6}& -\mathrm{8}\end{array}\right),B=\left(\begin{array}{ccc}\mathrm{4}& \mathrm{3}& \mathrm{8}\\ -\mathrm{7}& -\mathrm{1}& -\mathrm{3}\\ -\mathrm{3}& -\mathrm{5}& -\mathrm{6}\end{array}\right).\\ V_{\mathrm{1}}=\left(\begin{array}{cccccccc}\mathrm{16.0}& -\mathrm{1.6}& \mathrm{1.8}& -\mathrm{3.5}& \mathrm{1.3}& -\mathrm{2.0}& \mathrm{4.0}& -\mathrm{1.4}\\ -\mathrm{1.6}& \mathrm{25.0}& -\mathrm{2.2}& \mathrm{1.6}& -\mathrm{0.7}& \mathrm{0.5}& -\mathrm{1.3}& \mathrm{2.0}\\ \mathrm{1.8}& -\mathrm{2.2}& \mathrm{25.0}& -\mathrm{2.0}& \mathrm{5.0}& -\mathrm{2.4}& \mathrm{1.2}& -\mathrm{2.1}\\ -\mathrm{3.5}& \mathrm{1.6}& -\mathrm{2.0}& \mathrm{16.0}& -\mathrm{2.0}& \mathrm{3.0}& \mathrm{2.2}& \mathrm{2.8}\\ \mathrm{1.3}& -\mathrm{0.7}& \mathrm{5.0}& -\mathrm{2.0}& \mathrm{4.0}& -\mathrm{1.0}& \mathrm{0.8}& -\mathrm{2.0}\\ -\mathrm{2.0}& \mathrm{0.5}& -\mathrm{2.4}& \mathrm{3.0}& -\mathrm{1.0}& \mathrm{1.0}& -\mathrm{1.5}& \mathrm{0.6}\\ \mathrm{4.0}& -\mathrm{1.3}& \mathrm{1.2}& \mathrm{2.2}& \mathrm{0.8}& -\mathrm{1.5}& \mathrm{4.0}& -\mathrm{2.3}\\ -\mathrm{1.4}& \mathrm{2.0}& -\mathrm{2.1}& \mathrm{2.8}& -\mathrm{2.0}& \mathrm{0.6}& -\mathrm{2.3}& \mathrm{4.0}\end{array}\right)\\ V_{\mathrm{2}}=\left(\begin{array}{cccccccc}\mathrm{4.0}& -\mathrm{1.4}& \mathrm{0.8}& \mathrm{0.2}& \mathrm{1.6}& \mathrm{1.0}& \mathrm{1.2}& \mathrm{2.0}\\ -\mathrm{1.4}& \mathrm{4.0}& \mathrm{0.2}& -\mathrm{1.0}& -\mathrm{2.2}& \mathrm{0.8}& \mathrm{0.9}& \mathrm{1.8}\\ \mathrm{0.8}& \mathrm{0.2}& \mathrm{9.0}& \mathrm{0.2}& -\mathrm{1.5}& \mathrm{1.5}& \mathrm{1.0}& \mathrm{0.6}\\ \mathrm{0.2}& -\mathrm{1.0}& \mathrm{0.2}& \mathrm{36.0}& \mathrm{0.8}& \mathrm{0.4}& -\mathrm{1.5}& \mathrm{0.7}\\ \mathrm{1.6}& -\mathrm{2.2}& -\mathrm{1.5}& \mathrm{0.8}& \mathrm{25.0}& \mathrm{1.2}& -\mathrm{0.2}& \mathrm{2.0}\\ \mathrm{1.0}& \mathrm{0.8}& \mathrm{1.5}& \mathrm{0.4}& \mathrm{1.2}& \mathrm{25.0}& \mathrm{0.5}& \mathrm{1.4}\\ \mathrm{1.2}& \mathrm{0.9}& \mathrm{1.0}& -\mathrm{1.5}& -\mathrm{0.2}& \mathrm{0.5}& \mathrm{9.0}& \mathrm{0.8}\\ \mathrm{2.0}& \mathrm{1.8}& \mathrm{0.6}& \mathrm{0.7}& \mathrm{2.0}& \mathrm{1.4}& \mathrm{0.8}& \mathrm{16.0}\end{array}\right)\end{array}$$

According to model (8), problem (24) can be rewritten as follows:$$\begin{array}{c}\text{(26)}& \underset{x,y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}F\left(x,y\right)=\left(x^T,y^T\right)V_{\mathrm{1}}\left(\begin{array}{c}x\\ y\end{array}\right)\underset{y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}f\left(x,y\right)=\left(x,y\right)V_{\mathrm{2}}\left(\begin{array}{c}x\\ y\end{array}\right)\\ \text{s.t.}\hspace{1em}Ax+By\le {\left(\mathrm{215.85,140.07,24.41}\right)}^T\\ \hspace{1em}{m}_{c_{\mathrm{1}}}^{T}x+m_{d_{\mathrm{1}}}^{T}y\le \mathrm{60.0}\\ \hspace{1em}{m}_{c_{\mathrm{2}}}^{T}x+m_{d_{\mathrm{2}}}^{T}y\le \mathrm{35.0}\end{array}$$

For problem (26), we solve the problem with constrain of the means of the objective function values and we also solved it without this constrain. From Table 4, we can see that the variance obtained by our algorithm is better than the results in [51] when the mean of the upper objective function value is almost the same. Furthermore, the computational efficiency of our algorithm performs better with the number of dimension increasing.

Examples 1 and 2 are the general linear stochastic bilevel programming problem; that is to say, the coefficient of objective functions and the coefficient of constraint functions are random variables. Furthermore, the algorithm is also effective for the stochastic linear bilevel programming with chance constants. We consider the compensated stochastic linear bilevel programming as follows.

Table 3

The expectation of random variables $c_{\mathrm{1}},c_{\mathrm{2}},d_{\mathrm{1}},d_{\mathrm{2}}$.

$$m_{c_{\mathrm{1}}}$$	5	2	1	2	1	$$m_{d_{\mathrm{1}}}$$	6	1	3
$$m_{c_{\mathrm{2}}}$$	10	-7	1	-2	-5	$$m_{d_{\mathrm{2}}}$$	3	-4	6

Table 4

The comparison results of Example 2 from 11 runs.

Algorithm	means constrain	optimal solutions	means		variance		CPU time
			UL	LL	UL	LL	(s)
BPANN-PSO	Y	(0.0529, 0.2509, 0,0.0072, 0.9501; 0.5987,1.2321,1.1399)	9.9769	-2.2822	5.9689	77.7253	4.3898
Method in [51]	Y	(0.0530,0.2510, 0, 0.0070, 0.9500; 0.5990,1.2320, 1.1400)	9.9770	-2.2820	5.9700	77.7250	21.2172
BPANN-PSO	N	(0.4101,0.6320,0.0671,0.2702,0.8012; 0.2757,0.9731, 0.7488)	9.5979	-3.3727	20.6536	48.3552	0.9654
Method in [51]	N	(0.4100,0.6320,0.0670,0.2700, 0.8010,0.2760,0.9730,0.7490)	9.5980	-3.3720,	20.6540	48.3550	7.1836

Example 3.

Let $x\in R^{\mathrm{1}},y\in R^{\mathrm{1}}$, we consider the problem from [47]$$\begin{array}{c}\text{(27)}& \underset{x,y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}F\left(x,y\right)=-x-\mathrm{2}y\underset{y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}f\left(x,y\right)=-x+\mathrm{11}y\\ \text{s.t.}\hspace{1em}P\left(\begin{array}{c}{A}_{\mathrm{1}}x-B_{\mathrm{1}}y\le \mathrm{4}\\ \mathrm{2}{A}_{\mathrm{2}}x-B_{\mathrm{2}}y\le \mathrm{24}\\ \mathrm{3}{A}_{\mathrm{3}}+\mathrm{4}{B}_{\mathrm{3}}y\le \mathrm{96}\\ A_{\mathrm{4}}x+\mathrm{7}{B}_{\mathrm{4}}y\le \mathrm{126}\end{array}\right)\ge \mathrm{0.8},\\ \hspace{1em}x\ge \mathrm{0},\hspace{0.17em}\hspace{0.17em}y\ge \mathrm{0}.\end{array}$$where $A_{\mathrm{1}},A_{\mathrm{2}},\dots ,A_{\mathrm{4}},B_{\mathrm{1}},B_{\mathrm{2}},\dots ,B_{\mathrm{4}}~U(\mathrm{0,1})$, $p=\mathrm{0.8}$. In this paper, we set the 10 samples model and 30 samples model for problem (27), then, we obtained problem (28) and (29):

10 samples for model:$$\begin{array}{c}\text{(28)}& \underset{x,y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}F\left(x,y\right)=-x-\mathrm{2}y\underset{y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}f\left(x,y\right)=-x+\mathrm{11}y\\ \text{s.t.}\hspace{1em}\mathrm{0.516}x-\mathrm{1.6}y\le \mathrm{3.2}\\ \hspace{1em}\mathrm{1.133}x-\mathrm{0.5441}y\le \mathrm{19.2}\\ \hspace{1em}\mathrm{1.608}x+\mathrm{2.139}y\le \mathrm{76.8}\\ \hspace{1em}\mathrm{0.528}x+\mathrm{3.516}y\le \mathrm{100.8}\\ \hspace{1em}x\ge \mathrm{0},\hspace{0.17em}\hspace{0.17em}y\ge \mathrm{0}.\end{array}$$30 samples for model:$$\begin{array}{c}\text{(29)}& \underset{x,y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}F\left(x,y\right)=-x-\mathrm{2}y\underset{y}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}f\left(x,y\right)=-x+\mathrm{11}y\\ \text{s.t.}\hspace{1em}\mathrm{0.4325}x-\mathrm{1.556}y\le \mathrm{3.2}\\ \hspace{1em}\mathrm{0.793}x-\mathrm{0.4260}y\le \mathrm{19.2}\\ \hspace{1em}\mathrm{1.2186}x+\mathrm{1.6556}y\le \mathrm{76.8}\\ \hspace{1em}\mathrm{0.4865}x+\mathrm{3.4041}y\le \mathrm{100.8}\\ \hspace{1em}x\ge \mathrm{0},\hspace{0.17em}\hspace{0.17em}y\ge \mathrm{0}.\end{array}$$

We solve problem (28) and problem (29) by our algorithm and the method in [47]. From Table 5, we can see that the computational efficiency of our algorithm performs better when they have almost the same optimal solutions.

Table 5

comparison results of Example 3 from 11 runs.

	algorithm	optimal solutions	the optimal value		CPU run time
			UL	LL	(s)
10 sample model	BPANN-PSO	(64.8839, 18.9248)	-102.73335	78.4070	0.00179
	Method in [47]	(64.8840, 18.9250)	-102.7340	78.4070	0.0379
30 sample model	BPANN-PSO	(75.2431, 18.8569)	-112.9569	56.9433	0.00201
	Method in [47]	(75.2432, 18.8570)	-112.9572	56.9395	0.0463

5. Conclusion and Future Works

In this paper, we designed the BPANN-PSO algorithm to solve the general stochastic linear bilevel programming problem and three test problems from the reference are used to measure and evaluate the proposed algorithm. The results suggest that the variance obtained by our algorithm is better than the results in reference when the mean of the upper objective function value is almost the same. Furthermore, the computational efficiency of our algorithm performs better with the number of the dimensions increasing.

In the future works, we will further discuss how to efficiently use the infeasible solution with good performance near the optimal. This kind of discussion could improve the performance of our BPANN-PSO, particularly when the optimal front lies on the boundaries between the feasible and infeasible regions.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (61673006) and China Scholarship Council (201708420111).

References

[1] G. Brown, M. Carlyle, D. Diehl, J. Kline, K. Wood, "A two-sided optimization for theater ballistic missile defense," Operations Research, vol. 53 no. 5, pp. 745-763, DOI: 10.1287/opre.1050.0231, 2005.

[2] L. M. Wein, "Homeland security: From mathematical models to policy implementation: The 2008 Philip McCord Morse lecture," Operations Research, vol. 57 no. 4, pp. 801-811, DOI: 10.1287/opre.1090.0695, 2009.

[3] B. An, F. Ordóñez, M. Tambe, E. Shieh, R. Yang, C. Baldwin, J. DiRenzo, K. Moretti, B. Maule, G. Meyer, "A deployed quantal response-based patrol planning system for the u.s. coast guard," Interfaces, vol. 43 no. 5, pp. 400-420, DOI: 10.1287/inte.2013.0700, 2013.

[4] M. G. Nicholls, "Aluminum production modeling-a nonlinear bilevel programming approach," Operations Research, vol. 43 no. 2, pp. 208-218, DOI: 10.1287/opre.43.2.208, 1995.

[5] M. Labbé, P. Marcotte, G. Savard, "A bilevel model of taxation and its application to optimal highway pricing," Management Science, vol. 44 no. 12, pp. 1608-1622, DOI: 10.1287/mnsc.44.12.1608, 1998.

[6] A. Sinha, P. Malo, A. Frantsev, K. Deb, "Multi-objective Stackelberg game between a regulating authority and a mining company: A case study in environmental economics," Proceedings of the IEEE Congress on Evolutionary Computation, CEC '13, pp. 478-485, .

[7] A. Sinha, P. Malo, K. Deb, "Transportation policy formulation as a multi-objective bilevel optimization problem," Proceedings of the IEEE Congress on Evolutionary Computation, CEC '15, pp. 1651-1658, .

[8] X. Hu, D. Ralph, "Using EPECs to model bilevel games in restructured electricity markets with locational prices," Operations Research, vol. 55 no. 5, pp. 809-827, DOI: 10.1287/opre.1070.0431, 2007.

[9] N. Williams, P. K. Kannan, S. Azarm, "Retail channel structure impact on strategic engineering product design," Management Science, vol. 57 no. 5, pp. 897-914, DOI: 10.1287/mnsc.1110.1326, 2011.

[10] H. Sun, Z. Gao, J. Wu, "A bi-level programming model and solution algorithm for the location of logistics distribution centers," Applied Mathematical Modelling, vol. 32 no. 4, pp. 610-616, DOI: 10.1016/j.apm.2007.02.007, 2008.

[11] J. Bard, "Coordination of multi-divisional firm through two levels of management," Omega, vol. 11 no. 5, pp. 457-465, 1983.

[12] W. Smith, R. Missen, Chemical Reaction Equilibrium Analysis: Theory and Algorithms, 1982.

[13] P. A. Clark, A. W. Westerberg, "Bilevel programming for steady-state chemical process design-I. Fundamentals and algorithms," Computers & Chemical Engineering, vol. 14 no. 1, pp. 87-97, DOI: 10.1016/0098-1354(90)87007-C, 1990.

[14] M. P. Bendsoe, "Optimization of structural topology, shape, and material," Technical report,DOI: 10.1007/978-3-662-03115-5, 1995.

[15] S. Christiansen, M. Patriksson, L. Wynter, "Stochastic bilevel programming in structural optimization," Structural and Multidisciplinary Optimization, vol. 21 no. 5, pp. 361-371, DOI: 10.1007/s001580100115, 2001.

[16] K. Mombaur, A. Truong, J.-P. Laumond, "From human to humanoid locomotion-an inverse optimal control approach," Autonomous Robots, vol. 28 no. 3, pp. 369-383, DOI: 10.1007/s10514-009-9170-7, 2010.

[17] S. Albrecht, K. Ramírez-Amaro, F. Ruiz-Ugalde, D. Weikersdorfer, M. Leibold, M. Ulbrich, M. Beetz, "Imitating human reaching motions using physically inspired optimization principles," Proceedings of the 11th IEEE-RAS International Conference on Humanoid Robots, HUMANOIDS '11, pp. 602-607, .

[18] Q. Jin, S. Feng, "Bi-level simulated annealing algorithm for facility location," Systems Engineering, vol. 2, 2007.

[19] T. Uno, H. Katagiri, K. Kato, "An evolutionary multi-agent based search method for Stackelberg solutions of bilevel facility location problems," International Journal of Innovative Computing, Information and Control, vol. 4 no. 5, pp. 1033-1042, 2008.

[20] A. Migdalas, "Bilevel programming in traffic planning: models, methods and challenge," Journal of Global Optimization, vol. 7 no. 4, pp. 381-405, DOI: 10.1007/BF01099649, 1995.

[21] I. Constantin, M. Florian, "Optimizing frequencies in a transit network: a nonlinear bi-level programming approach," International Transactions in Operational Research, vol. 2 no. 2, pp. 149-164, DOI: 10.1016/0969-6016(94)00023-M, 1995.

[22] L. Brotcorne, M. Labbé, P. Marcotte, G. Savard, "A bilevel model for toll optimization on a multicommodity transportation network," Transportation Science, vol. 35 no. 4, pp. 345-358, DOI: 10.1287/trsc.35.4.345.10433, 2001.

[23] J. F. Bard, "An algorithm for solving the general bilevel programming problem," Mathematics of Operations Research, vol. 8 no. 2, pp. 260-272, DOI: 10.1287/moor.8.2.260, 1983.

[24] T. A. Edmunds, J. F. Bard, "Algorithms for nonlinear bilevel mathematical programs," The Institute of Electrical and Electronics Engineers Systems, Man, and Cybernetics Society, vol. 21 no. 1, pp. 83-89, DOI: 10.1109/21.101139, 1991.

[25] M. A. Amouzegar, "A global optimization method for nonlinear bilevel programming problems," IEEE Transactions on Systems, Man, and Cybernetics, Part B: Cybernetics, vol. 29 no. 6, pp. 771-777, DOI: 10.1109/3477.809031, 1999.

[26] J. B. Etoa Etoa, "Solving quadratic convex bilevel programming problems using a smoothing method," Applied Mathematics and Computation, vol. 217 no. 15, pp. 6680-6690, DOI: 10.1016/j.amc.2011.01.066, 2011.

[27] J. F. Bard, J. E. Falk, "An explicit solution to the multi-level programming problem," Computers & Operations Research, vol. 9 no. 1, pp. 77-100, DOI: 10.1016/0305-0548(82)90007-7, 1982.

[28] K. Shimizu, E. Aiyoshi, "A new computational method for Stackelberg and min-max problems by use of a penalty method," Institute of Electrical and Electronics Engineers Transactions on Automatic Control, vol. 26 no. 2, pp. 460-466, DOI: 10.1109/TAC.1981.1102607, 1981.

[29] E. Aiyoshi, K. Shimizu, "A solution method for the static constrained Stackelberg problem via penalty method," IEEE Transactions on Automatic Control, vol. 29 no. 12, pp. 1111-1114, DOI: 10.1109/TAC.1984.1103455, 1984.

[30] Y. Ishizuka, E. Aiyoshi, "Double penalty method for bilevel optimization problems," Annals of Operations Research, vol. 34 no. 1–4, pp. 73-88, DOI: 10.1007/BF02098173, 1992.

[31] Y. Lv, T. Hu, G. Wang, Z. Wan, "A penalty function method based on Kuhn-Tucker condition for solving linear bilevel programming," Applied Mathematics and Computation, vol. 188 no. 1, pp. 808-813, DOI: 10.1016/j.amc.2006.10.045, 2007.

[32] R. Mathieu, L. Pittard, G. Anandalingam, "Genetic algorithm based approach to bi-level linear programming," Operations Research, vol. 28 no. 1,DOI: 10.1051/ro/1994280100011, 1994.

[33] S. R. Hejazi, A. Memariani, G. Jahanshahloo, M. M. Sepehri, "Linear bilevel programming solution by genetic algorithm," Computers & Operations Research, vol. 29 no. 13, pp. 1913-1925, DOI: 10.1016/S0305-0548(01)00066-1, 2002.

[34] Y.-P. Wang, Y.-C. Jiao, H. Li, "An evolutionary algorithm for solving nonlinear bilevel programming based on a new constraint-handling scheme," IEEE Transactions on Systems, Man, and Cybernetics, Part C: Applications and Reviews, vol. 35 no. 2, pp. 221-232, DOI: 10.1109/TSMCC.2004.841908, 2005.

[35] G.-M. Wang, X.-J. Wang, Z.-P. Wan, S.-H. Jia, "An adaptive genetic algorithm for solving bilevel linear programming problem," Applied Mathematics and Mechanics-English Edition, vol. 28 no. 12, pp. 1605-1612, DOI: 10.1007/s10483-007-1207-1, 2007.

[36] H. I. Calvete, C. Galé, P. M. Mateo, "A new approach for solving linear bilevel problems using genetic algorithms," European Journal of Operational Research, vol. 188 no. 1, pp. 14-28, DOI: 10.1016/j.ejor.2007.03.034, 2008.

[37] X. Li, P. Tian, X. Min, "A hierarchical particle swarm optimization for solving bilevel programming problems," Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics): Preface, vol. 4029, pp. 1169-1178, DOI: 10.1007/11785231_122, 2006.

[38] R. J. Kuo, C. C. Huang, "Application of particle swarm optimization algorithm for solving bi-level linear programming problem," Computers & Mathematics with Applications. An International Journal, vol. 58 no. 4, pp. 678-685, DOI: 10.1016/j.camwa.2009.02.028, 2009.

[39] Y. Jiang, X. Li, C. Huang, X. Wu, "Application of particle swarm optimization based on CHKS smoothing function for solving nonlinear bilevel programming problem," Applied Mathematics and Computation, vol. 219 no. 9, pp. 4332-4339, DOI: 10.1016/j.amc.2012.10.010, 2013.

[40] S. B. Yaakob, J. Watada, "A hybrid intelligent algorithm for solving the bilevel programming models," Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics): Preface, vol. 6277 no. 2, pp. 485-494, DOI: 10.1007/978-3-642-15390-7_50, 2010.

[41] R. J. Kuo, Y. S. Han, "A hybrid of genetic algorithm and particle swarm optimization for solving bi-level linear programming problem—a case study on supply chain model," Applied Mathematical Modelling, vol. 35 no. 8, pp. 3905-3917, DOI: 10.1016/j.apm.2011.02.008, 2011.

[42] Z. Wan, G. Wang, B. Sun, "A hybrid intelligent algorithm by combining particle swarm optimization with chaos searching technique for solving nonlinear bilevel programming problems," Swarm and Evolutionary Computation, vol. 8, pp. 26-32, DOI: 10.1016/j.swevo.2012.08.001, 2013.

[43] M. Patriksson, L. Wynter, "Stochastic mathematical programs with equilibrium constraints," Operations Research Letters, vol. 25 no. 4, pp. 159-167, DOI: 10.1016/S0167-6377(99)00052-8, 1999.

[44] J. Gao, B. Liu, M. Gen, "A Hybrid Intelligent Algorithm for Stochastic Multilevel Programming," IEEJ Transactions on Electronics, Information and Systems, vol. 124 no. 10, pp. 1991-1998, DOI: 10.1541/ieejeiss.124.1991, 2004.

[45] Z. Wan, C. Xiao, X. Wang, K. Xiao, Y. Huang, X. Peng, "Bilevel programming model of optimal bidding strategies under the uncertain electricity markets," Dianli Xitong Zidonghua/Automation of Electric Power Systems, vol. 28 no. 19, pp. 12-16, 2004.

[46] Z. P. Wan, H. Fan, G. M. Wang, X. Y. Peng, Engineering Journal of Wuhan University, vol. 39 no. 5, pp. 85-90, 2006.

[47] Y. He, C. Q. Feng, "An approximating algorithm for compensated bilevel stochastic programming," Basic Sciences Journal of Textile Universities, vol. 26 no. 1, pp. 110-113, 2013.

[48] W. N. Zhou, Y. L. Huo, "Analysis of the Stability for the Solution Set to Stochastic Programming Approximation of Bilevel," Journal of Chongqing Technology Business University, vol. 30 no. 7, pp. 19-23, 2013.

[49] W. N. Zhou, Y. L. Huo, Z. Y. Hu, "The Hausdorff Covergence of the Optimal Solution Set of Approxiation for Bilevel Stochastic Programming," Journal of Natural Science of Hunan Normal University. Hunan Shifan Daxue Ziran Kexue Xuebao, vol. 39 no. 3, pp. 80-83, 2016.

[50] Y. Liu, H. Wang, Y. S. Xu, Y. Li, "Convergence of approximate solutions in bi-level stochastic programming," Pure and Applied Mathematics. Chuncui Shuxue yu Yingyong Shuxue, vol. 24 no. 4, pp. 768-773, 2008.

[51] S. Masatoshi, N. Ichiro, Cooperative and Non-cooperative Multi-Level Programming, 2009.

[52] T. Zhang, Z. Chen, J. Chen, "A cooperative coevolution PSO technique for complex bilevel programming problems and application to watershed water trading decision making problems," Journal of Nonlinear Sciences and Applications. JNSA, vol. 10 no. 4, pp. 2115-2132, DOI: 10.22436/jnsa.010.04.65, 2017.