Evolutionary algorithm using surrogate models for solving bilevel multiobjective programming problems

A bilevel programming problem with multiple objectives at the leader’s and/or follower’s levels, known as a bilevel multiobjective programming problem (BMPP), is extraordinarily hard as this problem accumulates the computational complexity of both hierarchical structures and multiobjective optimisation. As a strongly NP-hard problem, the BMPP incurs a significant computational cost in obtaining non-dominated solutions at both levels, and few studies have addressed this issue. In this study, an evolutionary algorithm is developed using surrogate optimisation models to solve such problems. First, a dynamic weighted sum method is adopted to address the follower’s multiple objective cases, in which the follower’s problem is categorised into several single-objective ones. Next, for each the leader’s variable values, the optimal solutions to the transformed follower’s programs can be approximated by adaptively improved surrogate models instead of solving the follower’s problems. Finally, these techniques are embedded in MOEA/D, by which the leader’s non-dominated solutions can be obtained. In addition, a heuristic crossover operator is designed using gradient information in the evolutionary procedure. The proposed algorithm is executed on some computational examples including linear and nonlinear cases, and the simulation results demonstrate the efficiency of the approach.


Problem models
The bilevel programming problem (BLPP), a hierarchical optimisation problem, has a nested optimisation structure because of a leader and a follower. In a BLPP, optimisation procedures are executed successively at both the leader's and follower's levels. The upper-level optimisation is provided by the leader, whereas the lower-level optimisation is controlled by the follower. The related problems are called the upper-level/leader's problem and lower-level/ follower's problem. Both the leader's and follower's problems have their own decision Where x = (x 1 , . . ., x n ) are the leader's variables, y = (y 1 , . . ., y m ) the follower's variables,F : R nþm ! R the leader's objective function,f : R nþm ! R the follower's objective function, G k : R n+m ! R, k = 1, 2, � � �, K, the leader's constraints, and g j : R n+m ! R, j = 1, 2, � � �, J, the follower's constraints. When the number of objectives in a leader's and/or follower's problems exceeds one, the problem is called a bilevel multiobjective programming problem (BMPP). The BMPP is one of the most difficult optimisation problems, because it accumulates all computational costs of hierarchical and multiobjective optimisations. Consequently, it is strongly NP-hard. The BMPP can be formulated as follows: Here, a � x � b, where a, and b are the lower and upper bounds of x, respectively. In general, the optimisation procedure of problem (2) can be described as follows. The leader first selects a strategy x to optimise his objective, and then the follower reacts to the leader's selection by providing a group of non-dominated solutions y i to the follower's problem. The pairs (x, y i ) are then used as bilevel feasible solutions to (2). When the leader selects another strategy, the related feasible point pairs can be obtained. To solve (2) is analogous to determining a set of non-dominated solutions that is distributed well for the leader's objective in all bilevel feasible solutions. Unlike BLPPs with a single objective at each level, the BMPP incurs an additional computational cost for the selection of the follower's non-dominated solutions. In addition, the selection of Pareto solutions in the leader's problem renders the BMPP much harder than the single objective case.

Related work
For problem (1), BLPPs have been extensively applied in economic management and engineering, urban traffic and transportation, supply chain planning, resource assignment, engineering design, structural optimisation, and game-playing strategies. For example, Zhu and Guo [2] proposed a BLPP with a maxmin or minmax operator in the follower's levels for a manufacturer, where the manufacturer plans to produce multiple short life cycle products using one-shot decision theory. The model was solved using typically used optimisation methods. Based on a bilevel complementary model, Nasrolahpouret et al. [3] developed an energy storage system decision tool for merchant price-making, that can determine the most advantageous trading behaviour in a pool-based market. More examples in real-world applications are summarised in [4][5][6][7][8].
As the applications of BLPPs are becoming increasingly extensive and diverse, researchers have focused on developing efficient solution strategies for such problems. However, owing to the intrinsic non-convexity and non-differentiability of BLPPs, a general BLPP is always complex and extremely challenging to solve using canonical optimisation methods that involve gradient information. Thus far, only a limited number of BLPPs can be solved effectively, such as linear [9][10][11][12][13] and convex quadratic BLPPs [14][15][16]. Liu et al. [10] used a new variant of the penalty method and Karush-Kuhn-Tucker(K-K-T) conditions to solve weak linear BLPPs. Franke et al. [16] used K-K-T conditions and the optimal value approach to transform a bilevel convex programming problem into a single-level optimisation problem, and introduced Mstationarity for mathematical programs with complementarities in Banach spaces. For most general BLPPs, only local optima can be obtained using these gradient approaches. To effectively solve BLPPs, another algorithm framework, i.e. the swarm optimisation method, has been widely adopted. This method is based on population search technology and affords good global convergence. Evolutionary algorithms [17][18][19], as representative swarm optimisation techniques, have been widely adopted to develop various bilevel optimisation algorithms in the past decades. For example, Sinha et al. [19] presented a single-level reduction of BLPPs using approximate K-K-T conditions and embedded the neighbourhood measure idea into an evolutionary algorithm. Based on a new constraint-handling scheme, Wang, et al. [20] proposed an evolutionary algorithm using K-K-T conditions. In Wang's method, the new constraint-handling scheme enables individuals to satisfy linear constraints.
For problem (2), even though it is extremely difficult to design efficient solution approaches for the BMPP, a large number of practical applications stimulate the researches on theoretical results and algorithmic design. Guo and Xu [21] developed a BMPP model to study the seismic risk of transportation system reconstruction in large construction projects, and fuzzy random variable transformation and fuzzy variable decomposition methods were proposed to solve the model. Brian et al. [22] proposed a BMPP model for coordinating multiple design problems according to conflicting criteria. The design of a hybrid vehicle layout was expressed as a twostage decomposition problem including vehicle class and battery class, and a multiobjective decomposition algorithm was developed. Alves et al. [23] developed a particle swarm optimisation algorithm to solve linear BMPPs with simple multiple objectives at the leader's level. Chakraborti et al. [24] investigated environmental-economic power generation and despatch problems using the BMPP.
For linear BMPPs, Calice et al. [34] first converted the linear BMPP into two multiobjective problems, and then used some specific optimisation conditions to prove that the common solutions to each multiobjective problem are optimal to the original problem. Kirlik et al. [35] proposed the use of a global optimisation method for discrete linear BMPPs. However, this method is computationally intensive. Liu et al. [36] investigated a class of pessimistic semi-rectorial BLPPs, and used secularisation to, transform the original problem into a scalar objective optimisation problem. Subsequently, the authors used the generalised differentiation calculus of Mordukhovich to establish optimality conditions for linear multiobjective problems. Lv and Wan [37] used the duality gap as a penalty on the leader's objective, and the follower's problem was transformed into a single objective case by adopting a weighted sum scheme. Alves [25] used a multiobjective particle swarm optimisation algorithm to solve linear BMPPs with multiple objective functions at the leader's level.
Only a few studies have been reported regarding nonlinear BMPPs. Deb and Sinha [29] presented an evolutionary-local-search algorithm for solving nonlinear BMPPs. In that study, the mapping method was used for the first time to solve the follower's programming problems. Jia et al. [30] used genetic algorithms to solve BMPPs, in which a uniform design was used to generate weights, crossover and mutation operators. In addition, the follower's objectives were converted into a single objective function by the weight sum method. In this method, MOEA/ D was used to solve the leader's multiobjective convex programming problem. Besides, based on the sensitivity theorem, an iterative method was used to solve a class of nonlinear BMPPs [33]. Zhang et al. [38] presented an improved particle swarm optimisation for solving BMPPs.

Research motivation
As mentioned above, because BMPPs have high computational complexity, only a few efficient solution methods exist for general BMPPs, and most of the existing approaches have been developed to solve cases where only the leader's problem is multiobjective. When the follower's problem involves multiple objectives, the non-dominated procedure can cause a large amount of computation. In addition, recall that in BLPPs, the optimisation procedure of the follower's problem is frequently executed to obtain bilevel feasible points, which can increase the computational cost. Consequently, if one intends to develop an efficient approach, two key issues must be addressed. One is to reduce the computational cost caused by the follower's optimisation procedure, whereas the other is to avoid excessive dominant comparisons among individuals. Hence, the weighted sum method is always adopted to delete non-dominated sorting procedures [40], and K-K-T conditions are applied to transform BLPPs into single-level ones. Furthermore, either MOEA/D [41,42] or NSGA-II [43,44] can be used to solve multiobjective optimisation problems. However, weight sum and K-K-T transformation techniques cannot perform the computation effectively when addressing the follower's problems. In the literature [26] Li et al used evolutionary algorithms to solve linear BMPPs in which the leader's objective is multiobjective. In Li's method, the programming problem was first transformed into a single-level multiobjective problem by K-K-T conditions; subsequently. The multiobjective problem was solved via the weighted sum method, Tchebycheff approach and NSGA-II method, separately. Finally, the computational results obtained using the three approaches were compared. Sinha et al. [31] presented an approximation of the set-mapping method to solve the follower's problem. This method can effectively reduce the amount of computation caused by the follower's optimisation. However, it is complicated to establish multiple quadratic fibers between the variables of the leader's and follower's problems. The methods used in the literature [26] and [31] requires a high computational cost.
Motivated by both weighted aggregation and approximate solution methods, in this study, an evolutionary algorithm using surrogate models was developed to solve BMPPs. Unlike the existing algorithms, the proposed algorithm is characterised as follows. First, a weighted aggregation method using a uniform design is adopted to convert the follower's problem into several single objective ones. Second, surrogate optimisation models and interpolation functions are used to replace the solution functions of the follower's problem and updated periodically using new sample points. Finally, as heuristic information, the gradient direction is embedded into genetic operators to produce potentially high-quality offspring.
Definition 0.1 (Dominated relations): . . . ; q, and there exists at least j 2 {1, 2, . . ., q}, such that F j 2 < F j 1 . If the objective value � Fðx 1 Þ of x 1 is dominated by the objective value � Fðx 2 Þ of x 2 , we also denote the relation by x 2 � x 1 . After a decision x 2 R n is made by the leader, the optimal solution set to the follower's problem is defined as follows: Assumptions: (A1) For each x taken in B 1 \B 2 by the leader, the follower has room to react, that is to say, O x 6 ¼ Ø.
(A2) F(x, y) is differential in x, and f(x, y) and g(x, y) is differentiable and convex in y for x fixed.
Here, A1 is a common assumption in BLPPs, which makes BLPPs well-posed. A2 is presented to conveniently utilize gradient information.

Transformation of follower problem
BMPPs will incur a high computational cost if the follower's problem is addressed as a general multiobjective optimisation. In the proposed approach, the solution of the follower's problem is categorised into two steps. First, multiobjective problem is converted into several singleobjective problems by the weighted sum of the objectives. Subsequently, these converted problems are solved using simplified surrogate models that can efficiently decrease the computational cost; this will be presented in the next Section. 1). The uniform design of weights: Uniform design and spherical coordinates are used to generate weights uniformly distrib- As a classical experimental design method, uniform design has become popular since the 1980's. It was originally developed to obtain designed points that are uniformly distributed over the experimental domain [45,46]. A uniform design array for k factors with q levels and H combinations is often denoted by L H (K q ). For example, L 16 (4 5 ) involves four factors, and each factor has five levels. Therefore, 1024(= 4 5 ) combinations exist. However, in a uniform design, one must only select 16 combinations from those 1024 possible combinations. H selected combinations can be represented by a uniform matrix U(n, H) = [U iℓ ] H × n , where U iℓ is the level of the ℓ th factor in the i th combination. In this subsection, we use the concept of uniform design to generate uniformly distributed points in 0; p 2 � � pÀ 1 . For a closed and bounded set φ � R p−1 , the main purpose of the uniform design is to sample a small group of points from set φ, such that the sampled points are uniformly scattered in set φ. Herein, we consider only sample points in the set 0; p 2 � � pÀ 1 . A brief description of the uniform design in the set is presented as follows.
For any point θ = (θ 1 , θ 2 , . . ., θ p−1 ) in φ, a hyper-rectangle between 0 and θ, denoted by φ (θ), can be defined as [40]: For a set of v points on φ, we assume that v(θ) points exist in the hyper-rectangle φ(θ). Therefore, the ratio of points in the hyper-rectangle φ(θ) is vðyÞ v , the volume of the hyper-rectangle is p 2 À � pÀ 1 , and the percentage of the hyper-rectangle volume is Subsequently, the uniform design on φ is defined as determining v points on φ such that the following discrepancy is minimised.
To determine uniformly scattered points on φ, we used the good-lattice-point method [47] to generate a set of v uniformly scattered points on φ, denoted by c(p − 1, v). The procedure is as follows: a v × (p − 1) uniform array is first generated: where Hence, we have For example, when v = 7 and p = 5, it can be seen that σ = 3 from Table 1, Thus: 2). Weight vector We can obtain v weight vectors as follows: where if p = 2, we have otherwise, if p > 2, we also have 3). Transformed problems We used the weight vectors designed in the above to address the follower's multiple objectives, transforming the follower's problem of the original problem into several single-objective ones. Subsequently, then (2) can be transformed into v BMPPs with a single objective in their follower's problems: The follower's problems of (14) are:

Surrogate models
In BLPPs/BMPPs, for each the leader's variable values, the follower's problem must be optimised. The procedure result in can a significant amount of computation in solving BLPPs/ BMPPs, particularly when the problem is large. According to the optimisation procedure, the optimal solutions to the follower's problem are always determined by the leader's variables. This means that the optimal solution of the follower's problem is a function of the leader's variables. However, the function is often implicit and can not be obtained analytically. In the proposed approach, we used the interpolation function as surrogate models to estimate the optimal solutions to the follower's problems. Therefore, some values x i of the leader's variables are first used, and then the follower's problems are optimised individually. The optimal solutions are denoted as y i , point pairs (x i , y i ) are regarded as interpolation sample points and used to generate interpolation functions. The interpolation function demonstrates better performance in fitting unknown functions [48] and can efficiently decrease the computational times of the follower's problems. In the proposed algorithm, the cubic spline interpolation method is adopted when the interpolation function is one-dimensional, whereas the linear interpolation method is used for other cases.
As an example, on the one-dimensional coordinate plane, for w sample points x i , i = 1, 2, . . ., w, we constructed an interpolation function y = ψ(x) to approximate the optimal solution function of the follower's problem. Therefore, y = ψ(x � ) can be regarded as the approximate solution to the follower's problem when x is fixed at x � . The relationship between the approximate and real solutions is shown in Fig 1. For interpolation functions, denser sample interpolation points resulted in a better approximation effect. It is noteworthy that for these interpolation points, each follower's variable value must be optimal when the leader's components are fixed. In the proposed algorithm, the interpolation function can be obtained as follows. First, an initial population of N points is randomly generated for the leader's variable space, and these points are denoted by x i , i = 1, 2, . . ., N. Subsequently, for each leader's value among x i , i = 1, 2, . . ., N, the optimal solutions to the follower's problem (15) are denoted by y ℓi , i = 1, 2, . . ., N, ℓ = 1, 2, . . ., v. Consequently, v × N point pairs of (x i , y ℓi ), i = 1, . . ., N, ℓ = 1, 2, . . ., v can be obtained. These point pairs are used as interpolation nodes to generate v interpolation functions. Where i.e. each of y j ' ; j ¼ 1; . . . ; m; ' ¼ 1; . . . ; v, is a function of x and y l ðxÞ ¼ ðy 1 ' ; y 2 ' ; . . . ; y m ' Þ; ' ¼ 1; . . . ; v. According to the above-mentioned interpolation method, we can obtain the approximate optimal solutions to the follower's problem (15). In the proposed algorithm, we periodically update both the interpolation sample points and to improve interpolation functions.

Proposed algorithm
When solving BMPPs, Transformation of follower problem guarantees the reduction of the computational complexity of the follower's problems, surrogate models guarantees the saving of the calculation cost of the follower's problems. Combining the methods in above, in this paper, an evolutionary algorithm [49][50][51]. based on surrogate models, denoted by SMEA, is developed to solve BMPPs. Fig 2 gives the flowchart of SMEA.
The specific procedure is described as follows: Step 0 (Transformation of the follower's problems) The weight vectors uniformly designed in section 3 are used to deal with the follower's multiple objectives, making the follower's problem become v single-objective ones. As a result, we obtain v BMPPs with different follower's objectives.
Step 3 (Crossover) Note that F k (x, y), k = 1, 2, . . ., q are differential in x. To obtain a potential descent direction, we take the negative gradient direction of each leader's function F k (x r , y r ).
Step 6 (Update interpolation function) Use MOEA/D to select non-dominated solutions from set {(x oi , y ℓi ), i = 1, 2, . . ., λ, ℓ = 1, 2, . . ., v}, and put these non-dominated solutions into the set D 2 . For each point in D 2 , update the optimal solutions y ℓi by solving the follower's problems. And update the interpolation function with the exact solutions obtained above.
Step 7 (Update Non-dominated solution set D 1 ) Set D 1 = D 1 S D 2 and D 2 = ϕ. Remove the dominant solutions in D 1 . Once the scale of D 1 exceeds the pre-determined threshold value, then the crowding distance method can be applied to delete some redundant points just as done in NSGA-II.
Step 8 (Selection) Select the best N individuals from set pop 0 (gen) S D 1 to form the next generation of population pop(gen + 1); Step 9 (Termination condition) If the stopping criterion is satisfied, then stop and output the non-dominated solutions set D 1 ; otherwise, set gen = gen + 1, go to Step 3.

Test examples
To demonstrate the feasibility and efficiency of the proposed algorithm, we test SMEA on 20 Examples which are taken from literatures [26,31,32,52] and. In [26] and [53], two EAs have been developed. In spite of the fact that the two algorithms are proposed only for dealing with BMPPs in which the follower's problem is single objective, as a special case, the two approaches can be taken as compared algorithms to demonstrate the performance of the proposed SMEA.

Parameter setting
To compare the experimental results, the parameter settings in this study are consistent with in [26].

Performance measure
To test the effectiveness and practicability of the SMEA algorithm, we tested the non-dominated solutions, HV, IGD, C − metric, CPU time and optimal solution according to the characteristics of the problems. Definitions of several metrics are as follows: 1). HV [54] Select a reference point in the objective space e = (e 1 , e 2 , � � �, e s ) T , denote A = SMEA, NSGA-II, Weighted sum approach or Tchebycheff approach. Compute For reference points in the objective space, the large HV means better quality.
2). IGD [55] Let G be a uniformly distributed subset selected to form the true Pareto Front and G 0 is the approximated set that is obtained by a multiobjective optimisation algorithm. The IGD value of G to G 0 , i.e., IGD(G, G 0 ) is defined as where |G| returns the number of solutions in the set G and d(G i , G 0 ) computes the minnmum Euclidean distance from G i to the solutions of G 0 in objective space. Generally, a lower value of IGD(G, G 0 )is preferred as it indicates that G 0 is distributed more uniformly and closer to the true Pareto front.
3). C − metric [54] We give the following notations: Set SM: non-dominated solution set of SMEA; Set WS: non-dominated solution set of Weighted sum approach; Set TE: non-dominated solution set of Tchebycheff approach; Set NS: non-dominated solution set of NSGA-II. Let:

Simulation results and comparisons
We execute an algorithm on a computer Intel(R) Core(TM)i5-8250U CPU@ 160 GHz 1.80 GHz using Matlab software. The algorithm was independently run ten times on every test instance. For the first 10 Examples, SMEA is compared with the weighted sum method, the Tchebycheff approach and NSGA-II. For each algorithm, we randomly take one among all 10 computational results of a single compared algorithm and show them in Figs 3-12. The analytical(theoretical) solution set of Examples 11 to 13 is provided, as a result, we can compare the computational results of SMEA with these known analytical solutions as shown in Figs 13-15. In Figs 3-12 (example 1-10), for Example 1 (shown in Fig 3), the non-dominated solutions set obtained by SMEA is superior to that obtained by the three compared methods. Example 2 is a maximization problem, from Fig 4 we can see that the non-dominated solutions set obtained by the SMEA is better than those obtained by the three compared methods. In addition, for Example 3, 4, 5, 7 and 9, the non-dominated solutions set obtained by the SMEA is similar to that by the Tchebycheff approach, but better than those by the other two methods. For Example 10, the results presented by the SMEA is almost the same as those obtained by both the Tchebycheff approach and the weighted sum approach, but better than that by the NSGA-II. For Examples 6 and 8, SMEA can also find approximately the non-dominated solutions provided by literature. Figs 13-15 show that the non-dominated solutions sets obtained by SMEA are almost in line with the analytical solutions sets. Table 2 represents the average value of HV obtained by SMEA, Weighted sum approach, Tchebycheff approach and NSGA-II running 10 times independently on Examples 1-10. While Table 3 represents the average value of HV and IGD obtained by SMEA and Analytical points running 10 times independently on Examples 11-13. The symbols "+", "−" and "�" mean the computational result is better than, worse than and almost equal to that obtained by

PLOS ONE
Evolutionary algorithm using surrogate models for solving bilevel multiobjective programming problems

PLOS ONE
Evolutionary algorithm using surrogate models for solving bilevel multiobjective programming problems

PLOS ONE
Evolutionary algorithm using surrogate models for solving bilevel multiobjective programming problems

PLOS ONE
Evolutionary algorithm using surrogate models for solving bilevel multiobjective programming problems

PLOS ONE
Evolutionary algorithm using surrogate models for solving bilevel multiobjective programming problems

PLOS ONE
Evolutionary algorithm using surrogate models for solving bilevel multiobjective programming problems our algorithm, respectively. Calculated the standard deviation(std)of HV value in 10 times. Best results are highlighted bold color in Table 2, we can see that the HV obtained by SMEA is much larger than those obtained other three algorithms on test problems 2-5, 7, 9 in Table 2. From Table 2, we can see in SMEA 10 values of HV better than the Weighted sum approach, 10 values of HV better than the Tchebycheff approach and 6 values of HV better than NSGA-II. Which illustrates that the diversity of solutions obtained by SMEA is better than obtained by other algorithms. Meanwhile, it is obvious that the IGD obtained by SMEA is small on Examples 11-13, which indicates for these test problems, the coverage of solutions obtained by SMEA quite well to the analytical solutions.
In terms of the results of the multiple problem Wilcoxon's test in Table 4, all the R + values are bigger than the R − ones, which implies that the performance of our algorithm is superior to that of other competitors. Moreover, the significant difference at α = 0.05 can be found in four cases, i.e., SMEA versus Tchebycheff approach, SMEA versus Weighted sum approach and SMEA versus NSGA-II. Besides, our algorithm ranks first in the Friedmans test(see Table 5 and Fig 16). Table 6 represents the average value of C-metrics running 10 times independently on Examples 1-10. We used SMEA to make a pairwise comparison with three algorithms in the literature.  (NS, SM) on problems 1, 2, 5, 6, 8, 10. It means that SMEA found a better non-dominated set than algorithms in [26].
In addition, CPU time is used to compare the efficiency of the algorithms. For the first 10 Examples, Table 7 shows the CPU times running 10 times independently used by both SMEA and two compared algorithms in [26]. From Table 7 one can see that SMEA can save more computation costs than the compared methods. For other Examples, we also provided their CPU times in Table 7.
It is difficult to solve BMPPs, especially to solve the follower's problem, and there exist only a small number of computational Examples. In the proposed algorithm, the weighted sum approach can transform the original multiobjective problem into single objective ones. Hence, SMEA can also be used to solve BLPPs with a single objective on which the effectiveness of the surrogate models can be verified. When SMEA is used in solving Examples 14-20, the best solution (x � , y � ) in all 10 runs is recorded as well as the corresponding objective values F(x � , y � ) and f(x � , y � ). All results are presented in Table 8 in which the objective values are denoted by Ref-f(x � , y � ) and Ref-F(x � , y � ), respectively. Table 8 shows the average value of optimal results running 10 times independently on Examples 14-20. The best results obtained are highlighted bold in this table. Especially, the optimal solutions of Examples 14-16 and 18 are obviously better than those provided in the existing literature. It follows that the surrogate model technique is effective to solve the problems of this type. Conclusion BMPP is one of the known hardest optimisation models in knowing that this problem always accumulates computational complexity of both the hierarchical structure and multi-objective optimization. In order to reduce the computational cost of the problem, two efficient techniques are embedded in the proposed algorithm. One is the weighted sum method used to deal with multiple objectives in the follower's problem. The other is the surrogate model which can efficiently save the computational cost in obtaining bilevel feasible solutions. In addition, a heuristic crossover operator is designed by making use of the gradient information.

PLOS ONE
Evolutionary algorithm using surrogate models for solving bilevel multiobjective programming problems