A Hybrid Estimation of Distribution Algorithm and Nelder-Mead Simplex Method for Solving a Class of Nonlinear Bilevel Programming Problems

We propose a hybrid algorithm based on estimation of distribution algorithm (EDA) and Nelder-Mead simplex method (NM) to solve a class of nonlinear bilevel programming problems where the follower’s problem is linear with respect to the lower level variable. The bilevel programming is an NP-hard optimization problem, for which EDA-NM is applied as a new tool aiming at obtaining global optimal solutions of such a problem. In fact, EDA-NM is very easy to be implemented since it does not require gradients information. Moreover, the hybrid algorithm intends to produce faster and more accurate convergence. In the proposed approach, for fixed upper level variable, we make use of the optimality conditions of linear programming to deal with the follower’s problem and obtain its optimal solution. Further, the leader’s objective function is taken as the fitness function. Based on these schemes, the hybrid algorithm is designed by combining EDA with NM. To verify the performance of EDA-NM, simulations on some test problems are made, and the results demonstrate that the proposed algorithm has a better performance than the compared algorithms. Finally, the proposed approach is used to solve a practical example about pollution charges problem.


Introduction
The bilevel programming problem (BLP) is a nested optimization problem with two levels (viz., the upper and lower levels) in a hierarchy.The upper level decision maker (the leader) and the lower level decision maker (the follower) control their own sets of decision variables, respectively, and have their own objective functions and constraints.The leader makes a decision first, and thereafter the follower chooses his/her strategy according to the leader's action.Anticipating the reaction of the follower, the leader selects the parameters so as to optimize his/her own objective function.The leader can influence, but cannot control, the decision of the follower.As optimization problems with hierarchical structure, the bilevel programming problem can be widely used in such areas as resource allocation, decentralized control, network design, and so forth [1].Over the past decades, the bilevel programming problem has been increasingly addressed in the literature including some useful reviews [1,2], surveys [3,4], and good textbooks [5,6].
However, it is extremely difficult to solve the bilevel programming problem due to its nonconvexity and noncontinuity, especially the nonlinear bilevel programming problem.For this context, many researchers have devoted themselves into developing efficient algorithms for solving the problem over the years.To date, there already have been some traditional approaches for solving such a problem, such as vertex enumeration method, approach based on the Kuhn-Tucker condition, descent approach, and penalty function approach.Most of these algorithms can accurately work out the optimal solution, but they are computationally costly even when solving small sized problems.As the scale of the problem increases, these exact algorithms could no longer afford the computation efficiently within reasonable time duration.
To overcome the above shortcomings, intelligent algorithms (such as evolutionary algorithm, genetic algorithm, tuba search approach, and simulated annealing) have been widely used to solve different problems in optimal areas and have been extended to deal with the bilevel programming for their good characteristics.For the linear bilevel programming problem, Mathieu et al. [7] firstly proposed a genetic algorithm for solving the problem.Subsequently, other kinds of intelligent algorithms [8][9][10][11] have been appearing for the same problem.For the nonlinear case, Wang et al. [12] presented an evolutionary algorithm for solving nonlinear bilevel programming problems with nonconvex objective functions.Deb and Sinha [13] designed a hybrid evolutionary-cum-local-search-based algorithm for dealing with bilevel multiobjective programming problems.In addition, a novel particle swarm optimization based on CHKS smoothing function was proposed in [14].By using duality conditions, an evolutionary algorithm was developed to cope with a class of bilevel programming with a linear lower level problem [15].Wan et al. [16] addressed a hybrid intelligent algorithm by combining the particle swarm optimization with chaos searching technique for solving nonlinear bilevel programming problems.
In recent years, an increasing interest has been concentrated on a class of probabilistic and graphical model based evolutionary computational methods, commonly called as the estimation of distribution algorithm (EDA) [17,18].Compared with other evolutionary algorithms, the EDA do not use crossover or mutation.Instead, it selects the best solutions from the current population and explicitly extracts global statistical information from the selected solutions.The EDA has some advantages to solve complex optimization problems, especially a high convergent reliability and low time consumption, and has been used widely in many real world problems.Nevertheless, to the best of our knowledge, there is no research work about the EDA for solving the bilevel programming problem.To enhance the performance of the EDA in our work, a local search algorithm, the Nelder-Mead simplex method, is integrated with the EDA to solve the bilevel programming.
In this paper, we consider a class of nonlinear bilevel programming problems where the follower's problem is linear only with respect to the lower level variable.To deal with such problems, based on the optimality results of linear programming, a hybrid algorithm is proposed by combing the EDA with a local simplex search technique.As a mater of fact, it can effectively integrate the characteristic of global search in the EDA and the capability of local search in simplex search technique to avoid converging ahead of time and to raise the accuracy of problem solving, as well as to make the search converge faster than pure EDA.Numerical simulations on some test problems are carried out, and the results demonstrate that the proposed algorithm has a better performance than the compared algorithms.
The remainder of the paper is organized as follows.Section 2 describes the formulation of the nonlinear bilevel programming problem and introduces the related definitions.Section 3 proposes a hybrid algorithm by combining the estimation of distribution algorithm and the Nelder-Mead simplex method for solving nonlinear bilevel programming problems.We conduct the simulations on the proposed algorithm and compare the results with some other algorithms for some test problems in Section 4. Finally, a conclusion is given in Section 5.

A Class of Nonlinear Bilevel Programming Problems
In this paper, we restrict our attention to the following nonlinear bilevel programming problem (NBLP): where  ∈  ⊂  Next we give the following definitions of the nonlinear bilevel programming problem: (1) constraint region: (2) for fixed , the feasible region of the lower level problem: (3) projection of  onto the upper level's decision space: (4) for each fixed  ∈  1 , the rational reaction set of the lower level problem: (5) inducible region: In order to ensure that the nonlinear bilevel programming problem is well posed, assume that in our work the constraint region  is nonempty and compact, and the lower level decision maker has some room to respond for all decisions taken by the upper level decision maker; that is () ̸ = 0.It is worth mentioning that the solution set () of the lower level problem is not always a singleton for a given  ∈  1 .In this case, the upper level problem is not a well-defined optimization problem.To overcome such an unpleasant situation, there are some strategies available for the leader, such as the optimistic approach and the pessimistic approach [6].In this paper we will avoid this unpleasant situation by assuming that there is a unique solution to the lower level problem for each fixed  ∈  1 .
Observing that the follower's problem is linear with respect to the lower level variable.For fixed , it is obvious that the follower's problem is a linear programming problem.Now the optimality conditions for the follower's problem are discussed according to the optimality conditions for the standard linear programming problem.For fixed , the term () is constant in the follower's problem, and we can ignore and delete it.As a result, the follower's problem can be written as the following problem by adding slack variables: where () = ((), ),  = (, 0 ) ∈  + , () = ((), 0), and  0 is a slack vector.For given , the term () is an  × ( + ) constant matrix, and its rank is .Assume that () ≥ 0. Given a set of indices  ⊂ {1, . . ., (+)}, let  be a basis of problem (7), and let  = {1, . . ., ( + )} −  be a nonbasis.A variable   with index  is called basic if  ∈ , nonbasic otherwise.Then we split (), (), and  into basis and nonbasis parts by means of  and  as index sets.Denoted by () = (()  , ()  ), ()  = (()   , ()   ), and  = (  ,   ).According to the theory of the simplex method, an optimal basis of problem (7) is characterized by the following inequalities: As a matter of fact, the inequalities (8) can be considered as the optimality conditions for the follower's problem.If the inequalities hold, then an optimal solution () = ((),  0 ()) is provided to the follower's problem; that is to say, the basic variable values of the solution are taken as () −1  (), while the values of nonbasic variables are 0. In essence, the presence of such a basis guarantees that there exists an optimal solution to the follower's problem.

Hybrid EDA-NM for Solving Nonlinear
Bilevel Programming Problems

Overview of Estimation of Distribution Algorithm.
The estimation of distribution algorithm (EDA), first proposed by Mühlenbein and Paaß [27] and later developed by Pelikan [18], is a new class of evolutionary optimization techniques.
Essentially, the EDA is also population-based algorithms similar to other evolutionary algorithms.The main difference between the EDA and other evolutionary algorithms is the way in which new offspring is generated in each generation.
There is neither crossover nor mutation operator in the EDA, while the algorithm selects the most promising individuals from the current population, constructs a probability model based on statistical information from the selected individuals, and then generates new offspring through sampling from the constructed probability model.The EDA works generally as follows.
Step 2. Select population of promising individuals from the current population and build the probability model of the selected individuals.
Step 3. Sample new offspring from the probability model.
Step 4. Replace some or all of the individuals in the previous population by the new offspring.
Step 5.If the stopping condition is met, stop.Otherwise, go to Step 2.
The EDA is able to overcome some drawbacks exhibited by traditional evolutionary algorithms.One of the biggest advantages of the EDA over other evolutionary algorithms is its ability to adapt their operators to the structure of the problem.Furthermore, this important difference allows the EDA to solve some problems for which other algorithms scale poorly.Nevertheless, the EDA pays more attention to global exploration, while its exploitation capability is relatively limited.So, an effective EDA should balance the exploration and the exploitation abilities.A very natural way to improve the performance of the EDA is to hybridize local search with the EDA.

Overview of Nelder-Mead Simplex Search
Method.The Nelder-Mead simplex search method (NM) [28] is a local search method developed for nonlinear and deterministic optimization without the need for gradient information.The NM algorithm can be implemented based on four basic procedures: reflection, expansion, contraction, and shrinkage depending on the values at the vertices and center of the simplex.To be more specific, an -dimensional simplex is defined as the convex hull of +1 vertices.If any vertex of an nondegenerate simplex is taken as the origin, then the rest  vertices define vector directions that span the -dimensional vector space.The extreme point of the simplex with the worst function value is moved, and the reflected point is generated.On one hand, if the reflected point is better than the best point, the expansion of the simplex is performed in one or another direction to take larger steps.On the other hand, a contraction step will be taken when the reflected point is the worst point in the new simplex, restricting the search on a smaller region.If the earlier worst point is better than the contracted one, the shrinkage step is performed.Through these operations, the simplex can improve itself and come closer and closer to a local optimum point sequentially.In the following, the above fundamental operations are given: where , , , and  are constants,  high is the worst vector corresponding to the maximum function value, and  cent is the center of the simplex excluding  high in the minimization case.
In essence, The NM can be regarded as a direct line-search method of steepest descent kind.Therefore, it is extremely flexible and very useful for exploring difficult problems.Moreover, the NM usually offers less computational costs when compared with evolutionary computational methods.However, the convergence properties of the NM are in general poor.

The Proposed Algorithm.
As is stated in previous section, the estimation of distribution algorithm is one of the global search algorithms.The algorithm is less likely to be entrapped in local optima but requires much computational effort.While the Nelder-Mead simplex method is a very efficient local search procedure, its convergence is extremely sensitive to the selected starting point.Therefore, the main idea of integrating the estimation of distribution algorithm with the Nelder-Mead simplex method in our work is to combine their advantages and avoid disadvantages.Furthermore, an effective algorithm should balance the exploration and the exploitation abilities.Consequently, to enhance the performance of the estimation of distribution algorithm, we use the Nelder-Mead simplex method as the local search mechanism to improve the efficiency of the search in exploitation.

The Structure of the Proposed Algorithm.
As far as the performance of evolutionary algorithms is concerned, a good initialization method is very important for improving the convergence speed and the quality of the final solution.For this purpose, uniform design, proposed by Fang [29] and Wang and Fang [30], is applied to design starting experiments, so that initial points could be scattered over the feasible region as uniformly as possible and the number of uniform design points should be required as small as possible.For more details on uniform design please refer to [31].
The initial population consists of  individuals associated with the upper level variable.To generate the initial population, uniform design is first used to generate  points  1 ,  2 , . . .,   from the search space of the upper level variable.To solve the follower's problem for given   ,  = 1, 2, . . ., , the inequalities (8) are checked for feasibility.If it is feasible, the initial point can be put into the initial population.If it is not feasible, the initial point can be discarded.Then, a new sampled point from the search space of the upper level variable is randomly generated, and we repeat the procedure until the inequalities (8) are satisfied, then select this sampled point to join the initial population.It is worth mentioning that the initial population is constructed by  different individuals to guarantee the initial population with certain quality and enough diversity.
In our work, the fitness function of the individual is defined as the leader's objective function.A total of  individuals are sorted according to the fitness function.Taking into account the balance between global exploration and local exploitation, the current population is divided into two parts, and population size is set  =  + ( + 1), where the top  individuals are updated by the estimation of distribution algorithm and the last  + 1 individuals are adjusted by the Nelder-Mead simplex method.
In the process of the EDA, a subset of the promising individuals from the top  individuals is selected with a selection method based on the fitness function.As a matter of fact, any selection method biased towards better fitness can be used, and in this paper, the traditional roulettewheel selection is applied to select  ( < ) individuals from the current  individuals.It is worth mentioning that the EDA produces offspring by sampling according to a probability model rather than crossover and mutation operators like other evolutionary algorithms.Therefore, it is crucial to choose a good probabilistic model in the EDA.
Here we employ a Gaussian model as a probabilistic model with very low computational cost.To be more specific, based on the statistical information extracted from these  selected individuals, a Gaussian model is constructed, and new individuals are generated according to the model.Then the inequalities (8) are checked for feasibility for each new individual, and the set of feasible individuals is formed.By this means,  new individuals are generated.The searching procedure is presented as follows.
Step 1. From ranked population, and select the top  individuals.
Step 2. Use roulette wheel to select a set of the promising  individuals from the current  individuals.
Step 3. Build a Gaussian model based on the statistical information extracted from these  selected individuals, and generate new individuals.
Step 4. For each new individual, check the inequalities (8).
At the local exploitation part, the algorithm works like the basic NM stated in previous section.In order to start up the NM, one has to define the initial simplex, in principle composed of  + 1 distinct vectors.For this purpose, the last  + 1 individuals from the current population based on rankbased fitness form the initial simplex thus avoiding extra function evaluations.Meanwhile the best point, the worst point and the second worst point are determined, and the center of the simplex excluding the worst point is computed.Then the operations of this method are to rescale the simplex based on the local behavior of the function by using four basic procedures: reflection, expansion, contraction, and shrinkage.For newly generated points at each step in these procedures, it is noted that the follower's problem is solved to get its optimal solution, so that we can compute its fitness value and continue the next process.To this end, we test the inequalities (8) for feasibility for each new individual.If the inequalities (8) are feasible, we continue the next process; otherwise, the new individual is replaced by the old individual and go to the next process.Through these procedures, the simplex is updated and replaced by the new simplex.In this way, we generate  + 1 new individuals.
The search process is outlined below.
Step 1. From ranked population, select the last  + 1 individuals to form a simplex.
Step 2. Apply four basic operators: reflection, expansion, contraction, and shrinkage.In each operator, the inequalities (8) are tested for feasibility for each new individual.
Step 3. Update the simplex.
The result is sorted in preparation for repeating the entire run.To guarantee the global convergence, we select the best  individuals from the current population and all the offspring generated by EDA and NM.These  individuals constitute the next population.If the algorithm is executed to the maximal number of generations and the best solution in the population has not been improved in successive generations, then stop.The best solution found in the last population is then taken as the approximate global optimal solution.

3.3.2.
Steps of the Proposed Algorithm.Based on the above procedure, the steps of the proposed algorithm are listed in details as follows.
Step 1 (initialization).Generate  individuals with respect to the upper level variable by means of the uniform design technique.Check the inequalities (8) for feasibility for each individual.Generate the initial population pop(0) by  individuals satisfying the inequalities (8).Set  = 0.
Step 2 (evaluation and ranking).Evaluate the fitness of each individual in current population and rank them based on the fitness results.
Step 3 (update).The first  individuals are updated by the estimation of distribution algorithm, and the +1 individuals in the end of the list are updated by the Nelder-Mead simplex method.
Step 4 (selection).Select the best  individuals among the current population and all the offspring.These selected individuals form the next population pop( + 1).
Step 5 (termination).An algorithm terminates when it satisfies a stopping criterion; otherwise, set  =  + 1 and go to Step 2.

Numerical Simulations
To illustrate the feasibility and efficiency of the proposed algorithm, a series of test problems are selected from the literature [11,12,15,[19][20][21][22][23][24][25][26].These problems are composed of different classes of bilevel programming problems, namely, the leader's problems involving linear, quadratic, fractional, nondifferentiable, and nonconvex objective functions.The variety of dimension and functional forms makes it possible to fairly assess the robustness of the proposed approach within tolerably computational time.In the simulation, the proposed EDA-NM is executed to solve all the test problems on MATLAB (R2011b) platform on Inter(R) Core (TM) i7 CPU 870 at 2.93 GHz having 8 G of RAM under Windows XP.

Comparison of the Proposed EDA-NM with the Existing
Algorithms.In this section, a comparison is made between the proposed EDA-NM with the existing algorithms.These approaches consist of hybrid of genetic algorithm and particle swarm optimization (HGAPSO) of Kuo and Han [11], evolutionary algorithm (EA) of Wang et al. [12], evolutionary algorithm (EA) based on duality conditions of Li and Fang [15], the monotonic approach of Tuy et al. [19], penalty function approach of Calvete and Galé [21], dual-relax penalty function approach of Wan et al. [22], and so forth.It is worth mentioning that some of these algorithms are stochastic, whereas others are deterministic.The algorithm proposed in our work is stochastic.To fairly measure the efficiency and stability of a stochastic algorithm, it is usually required to execute the algorithm in many runs and to compare the best and worst solutions found in these runs.
In our work, the proposed algorithm is executed for 50 independent runs on each of the test problems and record the following data: the best solution and the worst solution found in 50 independent runs, the values of the leader's and the follower's objective functions in the best and worst solutions, and average value (Avg.) and standard deviation (Std.) of the leader's objective function values in 50 independent runs.The parameters are set as follows: for -dimensional problem, the population size is  = 50, the value of  is from range of ( − )/5 to 2( − )/5, reflection coefficient, expansion coefficient, contraction coefficient, and shrinkage coefficient are set 1.0, 2.0, 0.5, and 0.5, and the maximum number of iterations is  = 50.The algorithm stops when the maximum number of iterations is achieved and the best results are unchanged in 10 successive generations.
The results are listed in Table 1.In Table 1, ( * ,  * ), (, ) represent the best and worst solutions found by an algorithm, respectively, and "Reference" stands for the algorithm in the corresponding references.The problems with star " * " are maximized models, whereas others are minimized BLPP.It can be seen from Table 1 that the best results found by EDA-NM are better than or equal to those by the compared algorithms in the references for all test problems except for problems 6, 10, and 14.For problems 4 and 9, the solutions found by EDA-NM are much better than those by the compared algorithms, which indicates that the compared algorithms cannot find the global solutions of the problems.Especially, for problem 9, the algorithm in Wan et al. [22] finds a local solution (0, 0.75, 0, 0.5, and 0) with the objective value of 10.625, while in our algorithm, all results found by EDA-NM in 50 runs are better than the results given by the compared algorithm in the references.Although the best results found by EDA-NM for problems 6, 10, and 14 are a little bit worse than the compared algorithms for these test problems, the precision of solutions found by our algorithm is higher than that provided by the compared algorithm.In addition, Table 1 shows that all of the worst solutions found by EDA-NM in 50 independent runs for all test problems are the same as or very close to the global optimal solutions.This means that the proposed algorithm is stable.

Comparison of EDA-NM with EDA.
In the previous section we have compared EDA-NM with the existing algorithms.In this section, we would compare EDA-NM with EDA.To conduct fair comparisons, we execute both algorithms for 50 independent runs on each of the test problems.The parameters of EDA are the same as those of EDA-NM.All the results are presented in Tables 2 to 3, in which Table 2 provides the comparison of the results found by EDA-NM and EDA for test problems in 50 runs, and Table 3 shows the best solutions found by EDA and EDA-NM in 50 runs and presented in the related references for all test problems.It can be seen from Table 2 that the best solutions found by EDA-NM are better than those found by EDA, and the worst solutions found by EDA-NM are as good as or are very close to the best solutions found by EDA.This means that our algorithm can find high-quality approximate global solutions for test problems.Tables 1-2 indicate that EDA-NM has smaller standard deviation, which means that EDA-NM is a stabler method.

Conclusions
In this paper, a hybrid estimation of distribution algorithm and Nelder-Mead simplex method is proposed to solve a class of nonlinear bilevel programming problems.EDA-NM is very easy to be implemented since it does not require gradients information.Moreover, the hybrid algorithm intends to produce faster and more accurate convergence.The performance of the proposed algorithm is tested on a series of test problems, and the results obtained are compared with those reported in the related references.The comparisons demonstrate that the proposed algorithm has a better performance than the compared algorithms.Moreover, the proposed approach is used to solve a practical example about pollution charges problem.Future work would involve the improvement on the probabilistic model of the EDA operator and the further refinement of the algorithm.

Table 1 :
The results found by EDA-NM and compared algorithms in the corresponding references.

Table 2 :
The results found by EDA with EDA-NM.

Table 3 :
The best solution found by EDA-NM, EDA, and compared algorithms in the corresponding references.