Abstract
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 implementedsince 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.
1. 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–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.
2. A Class of Nonlinear Bilevel Programming Problems
In this paper, we restrict our attention to the following nonlinear bilevel programming problem (NBLP): where , are the decision variables under the control of the upper and lower level problems, respectively; the sets and represent box sets of vectors and ; , are objective functions of the upper and lower level problems, respectively. , , , and whose rank is .
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 , 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 .
It is worth mentioning that the solution set of the lower level problem is not always a singleton for a given . 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 .
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 , , , and is a slack vector.
For given , the term is an constant matrix, and its rank is . Assume that . Given a set of indices , let be a basis of problem (7), and let 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 is provided to the follower’s problem; that is to say, the basic variable values of the solution are taken as , 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.
3. Hybrid EDA-NM for Solving Nonlinear Bilevel Programming Problems
3.1. 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 1. Randomly generate initial population.
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.
3.2. 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 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:
reflection: ,
expansion: ,
contraction: ,
shrink: ,
where , , , and are constants, is the worst vector corresponding to the maximum function value, and is the center of the simplex excluding 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.
3.3. 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.
3.3.1. 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 from the search space of the upper level variable. To solve the follower’s problem for given , , 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 , where the top individuals are updated by the estimation of distribution algorithm and the last 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 roulette-wheel 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 distinct vectors. For this purpose, the last individuals from the current population based on rank-based 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 new individuals.
The search process is outlined below.
Step 1. From ranked population, select the last 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 by individuals satisfying the inequalities (8). Set .
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 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 .
Step 5 (termination). An algorithm terminates when it satisfies a stopping criterion; otherwise, set and go to Step 2.
4. 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–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.
4.1. 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 , the value of is from range of to , 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 . 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.
4.2. 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.
In the analysis of the convergence reliability, we focus on the best solutions obtained by EDA-NM and EDA to achieve predefined convergence results. The convergence property of EDA-NM and EDA for different types of test problems is shown in Figures 1, 2, 3, 4, 5, and 6 in which the superiority of EDA-NM over EDA is evident. Figures 1–6 illustrate the performances of EDA-NM and EDA for the test problems 1, 2, 4, 5, 8, and 14 to find the global optimum by plotting the best fitness versus the number of iterations for a single run. From these figures, we can see that our hybrid EDA-NM algorithm converges more quickly and less likely to be trapped in local optima than EDA. Furthermore, EDA-NM method’s fitness drops quickly in maximized models and is raised rapidly in minimized BLPP at the evolution process. As a result, EDA-NM performs better than EDA from the aspects of computation efficiency.
From the above discussions, it is obvious that our algorithm is effective, is stable, and performs better than the compared algorithms.
4.3. Practical Application: Pollution Charges Problem
In this section, pollution charges decision-making problem based on bilevel programming is used to test the feasibility of the proposed algorithm. Pollution charges problem can be divided into two levels including the pollution tax collected by the government development and the pollutant discharge fees paid by the sewage enterprise, and each of them has its own objective function and decision variable. The government development first establishes pollution charges policy, and thereafter the sewage enterprise chooses its strategy according to the decision made by the government development. It is worth mentioning that the government development and the sewage enterprise independently seek their own interest, but influence each other. Therefore, pollution charges problem fits Stackelberg game model and can be solved by the bilevel programming model.
Assume that the sewage enterprise produces kinds of pollutant substance during the manufacturing process. The following notions will be used henceforth. , : cost prices of pollutant substances set by the government development; , : the pollutant discharge from the sewage enterprise; : the tax collected by the government development according to the ratio of the severity of the effects of the pollutant discharge in the environment; : the tax collected by the government development according to a certain percentage of the amount of the pollutant discharge from the sewage enterprise. Therefore, the total tax of the government department is . Subsequently, the sewage enterprise pays for pollutant discharge fees, that is, . In addition, assume that the variable and the variable have a certain constraint: , where and are matrix, , and is the shadow price caused by pollutants. To satisfy needs of maximizing the pollution tax collected by the government development and the cost of pollutant discharge of the sewage enterprise, the nonlinear bilevel programming is formulated as follows:
In the following, the numerical values of the parameters in above nonlinear bilevel problem are offered:
We use EDA-NM to deal with the bilevel programming problem. The parameters are selected as Section 4.1. We execute the algorithm in 20 independent runs and record the best solutions and the upper objective value as well as the lower objective value. The best optimal solution obtained by the algorithm is , and the optimal value is .
5. 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.
Acknowledgments
This work was supported by National Natural Science Foundation of China (no. 61272119), the Science Foundation of the Education, Department of Shaanxi Province of China (no. 2013JK0593), the Scientific Research Foundation of Xi’an Polytechnic University (BS1014), China Postdoctoral Science Foundation (20110491668), and National Natural Science Foundations of China (Grant nos. 11201362 and 110791).