# Solving Bilevel Multiobjective Programming Problem by Elite Quantum Behaved Particle Swarm Optimization

Authors: Tao Zhang; Tiesong Hu; Jia-wei Chen; Zhongping Wan; Xuning Guo
Journal: Abstract and Applied Analysis
 (2012)

---

## Abstract

An elite quantum behaved particle swarm optimization (EQPSO) algorithm is proposed, in which an elite strategy is exerted for the global best particle to prevent premature convergence of the swarm. The EQPSO algorithm is employed for solving bilevel multiobjective programming problem (BLMPP) in this study, which has never been reported in other literatures. Finally, we use eight different test problems to measure and evaluate the proposed algorithm, including low dimension and high dimension BLMPPs, as well as attempt to solve the BLMPPs whose theoretical Pareto optimal front is not known. The experimental results show that the proposed algorithm is a feasible and efficient method for solving BLMPPs.


---

## Body

## 1. Introduction

Bilevel programming problem (BLPP) arises in a wide variety of scientific and engineering applications including optimal control, process optimization, game-playing strategy development, transportation problem, and so on. Thus the BLPP has been developed and researched by many scholars. The reviews, monographs, and surveys on the BLPP can refer to [1–11]. Moreover, the evolutionary algorithms (EAs) have been employed to address BLPP in papers [12–16].For the multiobjective characteristics widely existing in the BLPP, the bilevel multiobjective programming problem (BLMPP) has attracted many researchers to study it. For example, Shi and Xia [17, 18], Abo-Sinna and Baky [19], Nishizaki and Sakawa [20], Zheng et al. [21] presented an interactive algorithm for BLMPP. Eichfelder [22] presented a method for solving nonlinear bilevel multiobjective optimization problems with coupled upper level constraints. Thereafter, Eichfelder [23] developed a numerical method for solving nonlinear nonconvex bilevel multiobjective optimization problems. In recent years, the metaheuristic has attracted considerable attention as an alternative method for BLMPP. For example, Deb and Sinha [24–26], as well as Sinha and Deb [27] discussed BLMPP based on evolutionary multiobjective optimization principles. Based on those studies, Deb and Sinha [28] proposed a viable and hybrid evolutionary-local-search based algorithm and presented challenging test problems. Sinha [29]   presented a progressively interactive evolutionary multiobjective optimization method for BLMPP.Particle swarm optimization (PSO) is a relatively novel heuristic algorithm inspired by the choreography of a bird flock, which has been found to be quite successful in a wide variety of optimization tasks [30]. Due to its high speed of convergence and relative simplicity, the PSO algorithm has been employed by many researchers for solving BLPPs. For example, Li et al. [31] proposed a hierarchical PSO for solving BLPP. Kuo and Huang [32] applied the PSO algorithm for solving bilevel linear programming problem. Gao et al. [33] presented a method to solve bilevel pricing problems in supply chains using PSO. However, it is worth noting that the papers mentioned above only for bilevel single objective problems and the BLMPP have seldom been studied using PSO so far. There are probably two reasons for this situation. One reason is that the added complexities associated with solving each level, and the other reason is that the global convergence of the PSO cannot be guaranteed [34].In this paper, a global convergence guaranteed method called as EQPSO is proposed, in which an elite strategy is exerted for global best particle to prevent premature convergence of the swarm. The EQPSO is employed for solving the BLMPP in this study, which has not been reported in other literatures. For such problems, the proposed algorithm directly simulates the decision process of the bilevel programming, which is different from most traditional algorithms designed for specific versions or based on specific assumptions. The BLMPP is transformed to solve multiobjective optimization problems in the upper level and the lower level interactively by the EQPSO. And a set of approximate Pareto optimal solutions for BLMPP are obtained using the elite strategy. This interactive procedure is repeated until the accurate Pareto optimal solutions of the original problem are found. The rest of the paper is organized as follows. In Section2, the problem formulation is provided. The proposed algorithm for solving bilevel multiobjective problem is presented in Section 3. In Section 4, some numerical examples are given to demonstrate the feasibility and efficiency of the proposed algorithm.


## 2. Problem Formulation

Letx∈Rn1,y∈Rn2, F:Rn1×Rn2→Rm1, f:Rn1×Rn2→Rm2, G:Rn1×Rn2→Rp, g:Rn1×Rn2→Rq. The general model of the BLMPP can be written as follows:
(2.1)minxF(x,y)s.t.G(x,y)≥0,minyf(x,y)s.t.g(x,y)≥0,
where F(x,y) and f(x,y) are the upper level and the lower level objective functions, respectively. G(x,y) and g(x,y) denote the upper level and the lower level constraints, respectively.LetS={(x,y)∣G(x,y)≥0,g(x,y)≥0}, X={x∣∃y,G(x,y)≥0,g(x,y)≥0}, S(x)={y∣g(x,y)≥0}, and for the fixed x∈X, let S¯(X) denote the weak efficiency set of solutions to the lower level problem, the feasible solution set of problem (2.1) is denoted as: IR={(x,y)∣(x,y)∈S,y∈S¯(X)}.Definition 2.1.
For a fixedx∈X, if y is a Pareto optimal solution to the lower level problem, then (x,y) is a feasible solution to the problem (2.1).Definition 2.2.
If(x*,y*) is a feasible solution to the problem (2.1), and there are no (x,y)∈IR, such that F(x,y)≺F(x*,y*), then (x*,y*) is a Pareto optimal solution to the problem (2.1), where “≺” denotes Pareto preference.For problem (2.1), it is noted that a solution (x*,y*) is feasible for the upper level problem if and only if y* is an optimal solution for the lower level problem with x=x*. In practice, we often make the approximate Pareto optimal solutions of the lower level problem as the optimal response feed back to the upper level problem, and this point of view is accepted usually. Based on this fact, the EQPSO algorithm may have a great potential for solving BLMPP. On the other hand, unlike the traditional point-by-point approach mentioned in Section 1, the EQPSO algorithm uses a group of points in its operation, thus the EQPSO can be developed as a new way for solving BLMPP. We next present the algorithm based on the EQPSO is presented for (2.1).


## 3. The Algorithm

### 3.1. The EQPSO

The quantum behaved particle swarm optimization (QPSO) is the integration of PSO and quantum computing theory developed by [35–38]. Compared with PSO, it needs no velocity vectors for particles and has fewer parameters to adjust. Moreover, its global convergence can be guaranteed [39]. Due to its global convergence and relative simplicity, it has been found to be quite successful in a wide variety of optimization tasks. For example, a wide range of continuous optimization problems [40–45] are solved by QPSO and the experiment results show that the QPSO works better than standard PSO. Some improved QPSO algorithms can refer to [46–48]. In this paper, the EQPSO algorithm is proposed, in which an elite strategy is exerted for global best particle to prevent premature convergence of the swarm, and it makes the proposed algorithm has good performance for solving the high dimension BLMPPS. The EQPSO has the same design principle with the QPSO except for the global optimal particle selection criterion, so the global convergence proof of the EQPSO can refer to [39]. In the EQPSO, the particles move according to the following iterative equation:
(3.1)zt+1=pt-αt(mBestt-zt)*ln(1u)ifk≥0.5,zt+1=pt+αt(mBestt-zt)*ln(1u)ifk<0.5,
where
(3.2)pt=φ*ppBestt+(1-φ)*pgBestt,mBestt=1N∑i=1NppBestit,αt=m-(m-n)*tT,pgBest∈rand(At),
where the z denotes the particle’s position. mBest denotes the mean best position of all the particles’ best positions. The k, u, and φ are random numbers distributed uniformly on (0,1), respectively. α(t) is the expansion-contraction coefficient. In general, m=1,n=0.5, t is the current iteration number, and T is the maximum number of iterations. The ppBest and pgBest are the particle’s personal best position and the global best position, respectively. At is the elite set which is introduced in following parts (see Algorithm: Step 3).


### 3.2. The Algorithm for Solving BLMPP

The process of the proposed algorithm for solving BLMPP is an interactive coevolutionary process. We first initialize population and then solve multiobjective optimization problems in the upper level and the lower level interactively using the EQPSO. For one time of iteration, a set of approximate Pareto optimal solutions for problem 1 is obtained by the elite strategy which was adopted in Deb et al. [49]. This interactive procedure is repeated until the accurate Pareto optimal solutions of problem (2.1) are found. The details of the proposed algorithm are given as follows.Algorithm
Step  1. Initializing.
Step  1.1.  Initialize the population P0 with Nu particles which is composed by ns=Nu/Nl subswarms of size Nl each. The particle’s position of the kth (k=1,2,…,ns) subswarm is presented as: zj=(xj,yj) (j=1,2,…,nl) and zj is sampled randomly in the feasible space.
Step  1.2.  Initialize the external loop counter t:=0.
Step  2.  For the kth subswarm, (k=1,2,…,ns), each particle is assigned a nondomination rank NDl and a crowding value CDl in f space. Then, all resulting subswarms are combined into one population which is named as the Pt. Afterwards, each particle is assigned a nondomination rank NDu and a crowding value CDu in F space.
Step  3.  The nondomination particles assigned both NDu=1 and NDl=1 from Pt are saved in the elite set At.
Step  4.  For the kth subswarm, (k=1,2,…,ns), update the lower level decision variables.
Step  4.1.  Initialize the lower level loop counter tl:=0.
Step  4.2.  Update the jth (j=1,2,…,Nl) particle’s position with the fixed xj according to (3.1) and (3.2).
Step  4.3.  tl:=tl+1.
Step  4.4.  If tl≥Tl, go to Step 4.5. Otherwise, go to Step 4.2
Step  4.5.  Each particle of the ith subswarm is reassigned a nondomination rank NDl and a crowding value CDl in f space. Then, all resulting subswarms are combined into one population which is renamed as the Qt. Afterwards, each particle is reassigned a nondomination rank NDu and a crowding value CDu in F space.
Step  5.  Combined population Pt and Qt to form Rt. The combined population Rt is reassigned a nondomination rank NDu, and the particles within an identical nondomination rank are assigned a crowding distance value CDu in the F space.
Step  6.  Choose half particles from Rt. The particles of rank NDu=1 are considered first. From the particles of rank NDu=1, the particles with NDl=1 are noted one by one in the order of reducing crowding distance CDu, for each such particle the corresponding subswarm from its source population (either Pt or Qt) is copied in an intermediate population St. If a subswarm is already copied in St and a future particle from the same subswarm is found to have NDu=NDl=1, the subswarm is not copied again. When all particles of NDu=1 are considered, a similar consideration is continued with NDu=2 and so on till exactly ns subswarms are copied in St.
Step  7.  Update the elite set At. The nondomination particles assigned both NDu=1 and NDl=1 from St are saved in the elite set At.
Step  8.  Update the upper level decision variables in St.
Step  8.1.  Initiate the upper level loop counter tu:=0.
Step  8.2.  Update the ith (i=1,2,…,Nu) particle’s position with the fixed yi according to (3.1) and (3.2).
Step  8.3.  tu:=tu+1.
Step  8.4.  If tu≥Tu, go to Step 8.5. Otherwise, go to Step 8.2.
Step  8.5.  Every member is then assigned a nondomination rank NDu and a crowding distance value CDu in F space.
Step  9.  t:=t+1.
Step  10.  If t≥T, output the elite set At. Otherwise, go to Step 2.In Steps 4 and 8, the global best position is chosen at random from the elite setAt. The criterion of personal best position choice is that, if the current position is dominated by the previous position, then the previous position is kept; otherwise, the current position replaces the previous one; if neither of them is dominated by the other, then we select one of them randomly. A relatively simple scheme is used to handle constraints. Whenever two individuals are compared, their constraints are checked. If both are feasible, nondomination sorting technology is directly applied to decide which one is selected. If one is feasible and the other is infeasible, the feasible dominates. If both are infeasible, then the one with the lowest amount of constraint violation dominates the other. Notations used in the proposed algorithm are detailed in Table 1.Table 1
The notations of the algorithm.
x
i
Theith particle’s position of the upper level problem.
y
j
Thejth particle’s position of the lower level problem.
z
j
Thejth particle’s position of BLMPP.
N
u
The population size of the upper level problem.
N
l
The sub-swarm size of the lower level problem.
t
Current iteration number for the overall problem.
T
The predefined max iteration number fort.
t
u
Current iteration number for the upper level problem.
t
l
Current iteration number for the lower level problem.
T
u
The predefined max iteration number fortu.
T
l
The predefined max iteration number fortl.
ND
u
Non-domination sorting rank of the upper level problem.
CD
u
Crowding distance value of the upper level problem.
ND
l
Non-domination sorting rank of the lower level problem.
CD
l
Crowding distance value of the lower level problem.
P
t
Thetth iteration population.
Q
t
The offspring ofPt.
S
t
Intermediate population.


## 3.1. The EQPSO

The quantum behaved particle swarm optimization (QPSO) is the integration of PSO and quantum computing theory developed by [35–38]. Compared with PSO, it needs no velocity vectors for particles and has fewer parameters to adjust. Moreover, its global convergence can be guaranteed [39]. Due to its global convergence and relative simplicity, it has been found to be quite successful in a wide variety of optimization tasks. For example, a wide range of continuous optimization problems [40–45] are solved by QPSO and the experiment results show that the QPSO works better than standard PSO. Some improved QPSO algorithms can refer to [46–48]. In this paper, the EQPSO algorithm is proposed, in which an elite strategy is exerted for global best particle to prevent premature convergence of the swarm, and it makes the proposed algorithm has good performance for solving the high dimension BLMPPS. The EQPSO has the same design principle with the QPSO except for the global optimal particle selection criterion, so the global convergence proof of the EQPSO can refer to [39]. In the EQPSO, the particles move according to the following iterative equation:
(3.1)zt+1=pt-αt(mBestt-zt)*ln(1u)ifk≥0.5,zt+1=pt+αt(mBestt-zt)*ln(1u)ifk<0.5,
where
(3.2)pt=φ*ppBestt+(1-φ)*pgBestt,mBestt=1N∑i=1NppBestit,αt=m-(m-n)*tT,pgBest∈rand(At),
where the z denotes the particle’s position. mBest denotes the mean best position of all the particles’ best positions. The k, u, and φ are random numbers distributed uniformly on (0,1), respectively. α(t) is the expansion-contraction coefficient. In general, m=1,n=0.5, t is the current iteration number, and T is the maximum number of iterations. The ppBest and pgBest are the particle’s personal best position and the global best position, respectively. At is the elite set which is introduced in following parts (see Algorithm: Step 3).


## 3.2. The Algorithm for Solving BLMPP

The process of the proposed algorithm for solving BLMPP is an interactive coevolutionary process. We first initialize population and then solve multiobjective optimization problems in the upper level and the lower level interactively using the EQPSO. For one time of iteration, a set of approximate Pareto optimal solutions for problem 1 is obtained by the elite strategy which was adopted in Deb et al. [49]. This interactive procedure is repeated until the accurate Pareto optimal solutions of problem (2.1) are found. The details of the proposed algorithm are given as follows.Algorithm
Step  1. Initializing.
Step  1.1.  Initialize the population P0 with Nu particles which is composed by ns=Nu/Nl subswarms of size Nl each. The particle’s position of the kth (k=1,2,…,ns) subswarm is presented as: zj=(xj,yj) (j=1,2,…,nl) and zj is sampled randomly in the feasible space.
Step  1.2.  Initialize the external loop counter t:=0.
Step  2.  For the kth subswarm, (k=1,2,…,ns), each particle is assigned a nondomination rank NDl and a crowding value CDl in f space. Then, all resulting subswarms are combined into one population which is named as the Pt. Afterwards, each particle is assigned a nondomination rank NDu and a crowding value CDu in F space.
Step  3.  The nondomination particles assigned both NDu=1 and NDl=1 from Pt are saved in the elite set At.
Step  4.  For the kth subswarm, (k=1,2,…,ns), update the lower level decision variables.
Step  4.1.  Initialize the lower level loop counter tl:=0.
Step  4.2.  Update the jth (j=1,2,…,Nl) particle’s position with the fixed xj according to (3.1) and (3.2).
Step  4.3.  tl:=tl+1.
Step  4.4.  If tl≥Tl, go to Step 4.5. Otherwise, go to Step 4.2
Step  4.5.  Each particle of the ith subswarm is reassigned a nondomination rank NDl and a crowding value CDl in f space. Then, all resulting subswarms are combined into one population which is renamed as the Qt. Afterwards, each particle is reassigned a nondomination rank NDu and a crowding value CDu in F space.
Step  5.  Combined population Pt and Qt to form Rt. The combined population Rt is reassigned a nondomination rank NDu, and the particles within an identical nondomination rank are assigned a crowding distance value CDu in the F space.
Step  6.  Choose half particles from Rt. The particles of rank NDu=1 are considered first. From the particles of rank NDu=1, the particles with NDl=1 are noted one by one in the order of reducing crowding distance CDu, for each such particle the corresponding subswarm from its source population (either Pt or Qt) is copied in an intermediate population St. If a subswarm is already copied in St and a future particle from the same subswarm is found to have NDu=NDl=1, the subswarm is not copied again. When all particles of NDu=1 are considered, a similar consideration is continued with NDu=2 and so on till exactly ns subswarms are copied in St.
Step  7.  Update the elite set At. The nondomination particles assigned both NDu=1 and NDl=1 from St are saved in the elite set At.
Step  8.  Update the upper level decision variables in St.
Step  8.1.  Initiate the upper level loop counter tu:=0.
Step  8.2.  Update the ith (i=1,2,…,Nu) particle’s position with the fixed yi according to (3.1) and (3.2).
Step  8.3.  tu:=tu+1.
Step  8.4.  If tu≥Tu, go to Step 8.5. Otherwise, go to Step 8.2.
Step  8.5.  Every member is then assigned a nondomination rank NDu and a crowding distance value CDu in F space.
Step  9.  t:=t+1.
Step  10.  If t≥T, output the elite set At. Otherwise, go to Step 2.In Steps 4 and 8, the global best position is chosen at random from the elite setAt. The criterion of personal best position choice is that, if the current position is dominated by the previous position, then the previous position is kept; otherwise, the current position replaces the previous one; if neither of them is dominated by the other, then we select one of them randomly. A relatively simple scheme is used to handle constraints. Whenever two individuals are compared, their constraints are checked. If both are feasible, nondomination sorting technology is directly applied to decide which one is selected. If one is feasible and the other is infeasible, the feasible dominates. If both are infeasible, then the one with the lowest amount of constraint violation dominates the other. Notations used in the proposed algorithm are detailed in Table 1.Table 1
The notations of the algorithm.
x
i
Theith particle’s position of the upper level problem.
y
j
Thejth particle’s position of the lower level problem.
z
j
Thejth particle’s position of BLMPP.
N
u
The population size of the upper level problem.
N
l
The sub-swarm size of the lower level problem.
t
Current iteration number for the overall problem.
T
The predefined max iteration number fort.
t
u
Current iteration number for the upper level problem.
t
l
Current iteration number for the lower level problem.
T
u
The predefined max iteration number fortu.
T
l
The predefined max iteration number fortl.
ND
u
Non-domination sorting rank of the upper level problem.
CD
u
Crowding distance value of the upper level problem.
ND
l
Non-domination sorting rank of the lower level problem.
CD
l
Crowding distance value of the lower level problem.
P
t
Thetth iteration population.
Q
t
The offspring ofPt.
S
t
Intermediate population.


## 4. Numerical Experiment

In this section, three examples will be considered to illustrate the feasibility of the proposed algorithm for problem (2.1). In order to evaluate the closeness between the obtained Pareto optimal front and the theoretical Pareto optimal front, as well as the diversity of the obtained Pareto optimal solutions along the theoretical Pareto optimal front, we adopted the following evaluation metrics.


### 4.1. Performance Evaluation Metrics

(a) Generational Distance (GD): this metric used by Deb [50] is employed in this paper as a way of evaluating the closeness between the obtained Pareto optimal front and the theoretical Pareto optimal front. The GD metric denotes the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front:
(4.1)GD=∑i=1ndi2n,
where n is the number of the obtained Pareto optimal solutions by the proposed algorithm and di is the Euclidean distance between each obtained Pareto optimal solution and the nearest member of the theoretical Pareto optimal set.(b) Spacing (SP): this metric is used to evaluate the diversity of the obtained Pareto optimal solutions by comparing the uniform distribution and the deviation of solutions as described by Deb [50]:
(4.2)SP=∑m=1Mdme+∑i=1n(d¯-di)2∑m=1Mdme+nd¯,
where di=minj(|F1i(x,y)-F1j(x,y)|+|F2i(x,y)-F2j(x,y)|),  i,j=1,2,…,n,  d¯ is the mean of all di,  dme is the Euclidean distance between the extreme solutions in obtained Pareto optimal solution set and the theoretical Pareto optimal solution set on the mth objective, M is the number of the upper level objective function, and n is the number of the obtained solutions by the proposed algorithm.All results presented in this paper have been obtained on a personal computer (CPU: AMD Phenon II X6 1055T 2.80 GHz; RAM: 3.25 GB) using a c# implementation of the proposed algorithm.


### 4.2. Numerical Examples

#### 4.2.1. Low Dimension BLMPPs

Example 4.1.
Example4.1 is taken from [22]. Here x∈R1,y∈R2. In this example, the population size and iteration times are set as follows: Nu=200, Tu=200, Nl=40, Tl=40, and T=40:
(4.3)minxF(x,y)=(y1-x,y2)s.t.G1(y)=1+y1+y2≥0minyf(x,y)=(y1,y2)s.t.g1(x,y)=x2-y12-y22≥0,-1≤y1,y2≤1,0≤x≤1.
Figure1 shows the obtained Pareto front of this example by the proposed algorithm. From Figure 1, it can be seen that the obtained Pareto front is very close to the theoretical Pareto optimal front, and the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front is 0.00024, that is, GD=0.00024 (see Table 2). Moreover, the lower SP value (SP=0.0042, see Table 2) shows that the proposed algorithm is able to obtain a good distribution of solutions on the entire range of the theoretical Pareto optimal front. Figure 2 shows the obtained solutions of this example, which follow the relationship, that is, y1=-1-y2, y2=-1/2±(1/4)8x2-4 and x∈(1/2,1). It is also obvious that all obtained solutions are close to being on the upper level constraint G(x) boundary (1+y1+y2=0).Table 2
Results of the Generation Distance (GD) and Spacing (SP) metrics for the above six examples.
Example
GD
SP
Example4.1
0.00024
0.00442
Example4.2
0.00003
0.00169
Example4.3
0.00027
0.00127
Example4.4
0.00036
0.00235
Example4.5
0.00058
0.00364
Example4.6
0.00039
0.00168Figure 1
The obtained Pareto optimal front of Example4.1.Figure 2
The obtained solutions of Example4.1.Example 4.2.
Example4.2 is taken from [51]. Here x∈R1,y∈R2. In this example, the population size and iteration times are set as follows: Nu=200, Tu=50, Nl=40, Tl=20, and T=40:
(4.4)minxF(x,y)=(x2+(y1-1)2+y22,(x-1)2+(y1-1)2+y22)minyf(x,y)=(y12+y22,(y1-x)2+y22)-1≤x,y1,y2≤2.
Figure3 shows the obtained Pareto optimal front of this example by the proposed algorithm. From Figure 3, it is obvious that the obtained Pareto optimal front is very close to the theoretical Pareto optimal front, the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front is 0.00003 (see Table 2). On the other hand, the obtained Pareto optimal solutions can distribute uniformly on entire range of theoretical Pareto optimal front base on the fact that the SP value is lower (SP=0.00169, see Table 2). Figure 4 shows the obtained Pareto optimal solutions, they follow the relationship, that is, x=y1,y1∈[0.5,1] and y2=0.Figure 3
The obtained Pareto optimal front of Example4.2.Figure 4
The obtained solutions Example4.2.


#### 4.2.2. High Dimension BLMPPs

Example 4.3.
Example4.3 is taken from [28]. Here x∈R10,y∈R10. In this example, the population size and iteration times are set as follows: Nu=400, Tu=50, Nl=40, Tl=20, and T=60:
(4.5)minxF(x,y)=((1+r-cos(απx1))+∑j=2K(xj-j-12)2+τ∑i=2K(yi-xi)2-γcos(γπx12y1),(1+r-sin(απx1))+∑j=2K(xj-j-12)2+τ∑i=2K(yi-xi)2-γsin(γπx12y1)∑j=2K(xj-j-12)2)minyf(x,y)=(y12+∑i=2K(yi-xi)2+∑i=2K10(1-cos(πk(yi-xi))),∑i=1K(yi-xi)2+∑i=2K10|1-sin(πk(yi-xi))|)s.t.-K≤yi≤K,(i=1,2,…,K);1≤x1≤4,-K≤xj≤K,(j=2,3,…,K).α=1,r=0.1,τ=1,γ=1,K=10.Example 4.4.
Example4.4 is taken from [28]. Here x∈R10,y∈R10. In this example, the population size and iteration times are set as follows: Nu=400,Tu=50, Nl=40, Tl=20, and T=80. 											(4.6)minxF(x,y)=v1(x1)+∑j=2Kyj2+10(1-cos(πk)yi)+τ∑i=2K(yi-xi)2-rcos(γπx12y1),v2(x1)+∑j=2Kyj2+10(1-cos(πk)yi)+τ∑i=2K(yi-xi)2-rsin(γπx12y1),minyf(x,y)=(y12+∑i=2K(yi-xi)2,∑i=1Ki(yi-xi)2)s.t.-K≤yi≤K,(i=1,2,…,K);0.001≤x1≤4,-K≤xj≤K,(j=2,3,…,K).α=1,r=0.25,τ=1,γ=4,K=10,
where
(4.7)v1(x1)={cos(0.2π)x1+sin(0.2π)|0.02sin(5πx1)|,for0≤x1≤1,x1-(1-sin(0.2π)),forx1>1,v2(x1)={-sin(0.2π)x1+cos(0.2π)|0.02sin(5πx1)|,for0≤x1≤1,0.1(x1-1)-sin(0.2π),forx1>1.
This problem is more difficult compared to the previous problems (Examples4.1 and 4.2) because the lower level problem of this example has multimodalities, thereby making the lower level problem difficult in finding the upper level Pareto optimal front. From Figure 5, it can be seen that the obtained Pareto front is very close to the theoretical Pareto optimal front, and the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front is 0.00027, that is, GD=0.00027 (see Table 2). Moreover, the lower SP value (SP=0.00127, see Table 2) shows that the proposed algorithm is able to obtain a good distribution of solutions on the entire range of the theoretical Pareto optimal front. Furthermore, two obtained lower level Pareto optimal fronts are given when x1=2 and x1=2.5.Figure 5
The obtained Pareto front of Example4.3.Figure6 shows the obtained Pareto front of Example 4.4 by the proposed algorithm. The upper level problem has multimodalities, thereby causing an algorithm difficulty in finding the upper level Pareto optimal front. From Figure 6, it can be seen that the obtained Pareto front is very close to the theoretical Pareto optimal front, and the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front is 0.00036, that is, GD=0.00036 (see Table 2). Moreover, the lower SP value (SP=0.00235, see Table 2) shows that the proposed algorithm is able to obtain a good distribution of solutions on the entire range of the theoretical Pareto optimal front. Furthermore, all corresponding lower level Pareto optimal fronts are given.Figure 6
The obtained Pareto front of Example4.4.Example 4.5.
Example4.5 is taken from [28]. Here x∈R10,y∈R10. In this example, the population size and iteration times are set as follows: Nu=400,Tu=50, Nl=40, Tl=20, and T=60:
(4.8)minxF(x·y)=(x1+∑j=3K(xj-j2)2+τ∑i=3K(yi-xi)2-cos(4tan-1(x2-y2x1-y1)),x2+∑j=3K(xj-j2)2+τ∑i=3K(yi-xi)2-cos(4tan-1(x2-y2x1-y1)))s.t.G(x)=x2-(1-x12)2≥0minyf(x,y)=(y1+∑i=3K(yi-xi)2,y2+∑i=3K(yi-xi)2)s.t.g1(x,y)=(y1-x1)2+(y2-x2)2≤r2R(x1)=(0.1+0.15|sin(2π(x1-0.1))|),-K≤yi≤K,(i=1,2,…,K);0≤x1≤K,-K≤xj≤K,(j=2,3,…,K).τ=1,r=0.2,K=10.Example 4.6.
Example4.6 is taken from [28]. Here x∈R1,y∈R9. In this example, the population size and iteration times are set as follows: Nu=400,Tu=50, Nl=40, Tl=20, and T=40:
(4.9)minxF(x,y)=((1-y1)(1+∑j=2Kyj2)x1,y1(1+∑j=2Kyj2)x1)s.t.G1(x,y)=-(1-y1)x1-12x1y1≤1minyf(x,y)=((1-y1)(1+∑j=K+1K+Lyj2)x1,y1(1+∑j=K+1K+Lyj2)x1)s.t.1≤x1≤2,-1≤y1≤1,-(K+L)≤yj≤(K+L),(j=2,3,…,K+L),K=5,L=4.
Figure7 shows the obtained Pareto front of Example 4.5 by the proposed algorithm. From Figure 7, it can be seen that the obtained Pareto front is very close to the theoretical Pareto optimal front, and the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front is 0.00058, that is, GD=0.00058 (see Table 2). Moreover, the lower SP value (SP=0.00364, see Table 2) shows that the proposed algorithm is able to obtain a good distribution of solutions on the entire range of the theoretical Pareto optimal front. Furthermore, all obtained lower level Pareto optimal fronts are given. It is also obvious that the Pareto optimal fronts for both the lower and upper level lie on constraint boundaries and every lower level Pareto optimal front has an unequal contribution to the upper level Pareto optimal front.
Figure8 shows the obtained Pareto front of Example 4.6 by the proposed algorithm. From Figure 8, it can be seen that the obtained Pareto front is very close to the theoretical Pareto optimal front, and the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front is 0.00039, that is, GD=0.00039 (see Table 2). Moreover, the lower SP value (SP=0.00168, see Table 2) shows that the proposed algorithm is able to obtain a good distribution of solutions on the entire range of the theoretical Pareto optimal front. Furthermore, three obtained lower level Pareto optimal fronts are given when y1=1, y1=1.5 and y1=2. It can be seen that only one Pareto optimal point from each participating lower level problem qualifies to be on the upper level Pareto optimal front.Figure 7
The obtained Pareto front of Example4.5.Figure 8
The obtained Pareto front of Example4.6.


#### 4.2.3. The BLMPPs with Unknown Theoretical Pareto Optimal Fronts

Example 4.7.
Example4.7 is taken from [52], in which the theoretical Pareto optimal front is not given. Here x∈R2,y∈R3. In this example, the population size and iteration times are set as follows: Nu=100,Tu=50, Nl=20, Tl=10, and T=40:
(4.10)maxxF(x,y)=(x1+9x2+10y1+y2+3x3,9x1+2x2+2y1+7y2+4x3)s.t.G1(x,y)=3x1+9x2+9y1+5y2+3y3≤1039G2(x,y)=-4x1-x2+3y1-3y2+2y3≤94minyf(x,y)=(4x1+6x2+7y1+4y2+8y3,6x1+4x2+8y1+7y2+4y3)s.t.g1(x,y)=3x1-9x2-9y1-4y2≤61g2(x,y)=5x1+9x2+10y1-y2-2y3≤924g3(x,y)=3x1-3x2+y2+5y3≤420x1,x2,y1,y2,y3≥0.Example 4.8.
Example4.8 is taken from [23]. Here x∈R1,y∈R2. In this example, the population size and iteration times are set as follows: Nu=800,Tu=50,Nl=40, Tl=20, and T=40:
(4.11)minxF(x,y)=(x+y12+y22+sin2(x1+y),cos(y2)(0.1+y)(exp(-y10.1+y2)))s.t.G1(x,y)=(x-5)2-(y1-0.5)2-(y2-5)2≤16minyf(x,y)=((y1-2)2+(y2-2)24+xy2+(5-x)216+sin(y210),y12+(y2-6)2-2xy1-(5-x)280)s.t.g1(x,y)=y12-y2≤0g2(x,y)=5y12+y2-10-≤0g3(x,y)=y2-5-x6≤00≤x≤10,0≤y1,y2≤10.
Figure9 shows the obtained Pareto optimal front of Example 4.7 by the proposed algorithm. Note that, Zhang et al. [52] only obtained a single optimal solution x=(146.2955,28.9394) and y=(0,67.9318,0) which lies on the maximum of the F2 using weighted sum method. In contrast, a set of Pareto optimal solutions are obtained by the proposed algorithm. However, the fact that the single optimal solution in [49] is included in the obtained Pareto optimal solutions illustrates the feasibility of proposed algorithm. Figure 10 shows the final archive solutions of the Example 4.8 by the proposed algorithm. For this problem, the exact Pareto optimal front is not known, but the obtained Pareto optimal front by the proposed algorithm is similar to that reported in the previous study [23].Figure 9
The obtained front of Example4.7.Figure 10
The obtained front of Example4.8.


## 4.1. Performance Evaluation Metrics

(a) Generational Distance (GD): this metric used by Deb [50] is employed in this paper as a way of evaluating the closeness between the obtained Pareto optimal front and the theoretical Pareto optimal front. The GD metric denotes the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front:
(4.1)GD=∑i=1ndi2n,
where n is the number of the obtained Pareto optimal solutions by the proposed algorithm and di is the Euclidean distance between each obtained Pareto optimal solution and the nearest member of the theoretical Pareto optimal set.(b) Spacing (SP): this metric is used to evaluate the diversity of the obtained Pareto optimal solutions by comparing the uniform distribution and the deviation of solutions as described by Deb [50]:
(4.2)SP=∑m=1Mdme+∑i=1n(d¯-di)2∑m=1Mdme+nd¯,
where di=minj(|F1i(x,y)-F1j(x,y)|+|F2i(x,y)-F2j(x,y)|),  i,j=1,2,…,n,  d¯ is the mean of all di,  dme is the Euclidean distance between the extreme solutions in obtained Pareto optimal solution set and the theoretical Pareto optimal solution set on the mth objective, M is the number of the upper level objective function, and n is the number of the obtained solutions by the proposed algorithm.All results presented in this paper have been obtained on a personal computer (CPU: AMD Phenon II X6 1055T 2.80 GHz; RAM: 3.25 GB) using a c# implementation of the proposed algorithm.


## 4.2. Numerical Examples

### 4.2.1. Low Dimension BLMPPs

Example 4.1.
Example4.1 is taken from [22]. Here x∈R1,y∈R2. In this example, the population size and iteration times are set as follows: Nu=200, Tu=200, Nl=40, Tl=40, and T=40:
(4.3)minxF(x,y)=(y1-x,y2)s.t.G1(y)=1+y1+y2≥0minyf(x,y)=(y1,y2)s.t.g1(x,y)=x2-y12-y22≥0,-1≤y1,y2≤1,0≤x≤1.
Figure1 shows the obtained Pareto front of this example by the proposed algorithm. From Figure 1, it can be seen that the obtained Pareto front is very close to the theoretical Pareto optimal front, and the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front is 0.00024, that is, GD=0.00024 (see Table 2). Moreover, the lower SP value (SP=0.0042, see Table 2) shows that the proposed algorithm is able to obtain a good distribution of solutions on the entire range of the theoretical Pareto optimal front. Figure 2 shows the obtained solutions of this example, which follow the relationship, that is, y1=-1-y2, y2=-1/2±(1/4)8x2-4 and x∈(1/2,1). It is also obvious that all obtained solutions are close to being on the upper level constraint G(x) boundary (1+y1+y2=0).Table 2
Results of the Generation Distance (GD) and Spacing (SP) metrics for the above six examples.
Example
GD
SP
Example4.1
0.00024
0.00442
Example4.2
0.00003
0.00169
Example4.3
0.00027
0.00127
Example4.4
0.00036
0.00235
Example4.5
0.00058
0.00364
Example4.6
0.00039
0.00168Figure 1
The obtained Pareto optimal front of Example4.1.Figure 2
The obtained solutions of Example4.1.Example 4.2.
Example4.2 is taken from [51]. Here x∈R1,y∈R2. In this example, the population size and iteration times are set as follows: Nu=200, Tu=50, Nl=40, Tl=20, and T=40:
(4.4)minxF(x,y)=(x2+(y1-1)2+y22,(x-1)2+(y1-1)2+y22)minyf(x,y)=(y12+y22,(y1-x)2+y22)-1≤x,y1,y2≤2.
Figure3 shows the obtained Pareto optimal front of this example by the proposed algorithm. From Figure 3, it is obvious that the obtained Pareto optimal front is very close to the theoretical Pareto optimal front, the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front is 0.00003 (see Table 2). On the other hand, the obtained Pareto optimal solutions can distribute uniformly on entire range of theoretical Pareto optimal front base on the fact that the SP value is lower (SP=0.00169, see Table 2). Figure 4 shows the obtained Pareto optimal solutions, they follow the relationship, that is, x=y1,y1∈[0.5,1] and y2=0.Figure 3
The obtained Pareto optimal front of Example4.2.Figure 4
The obtained solutions Example4.2.


### 4.2.2. High Dimension BLMPPs

Example 4.3.
Example4.3 is taken from [28]. Here x∈R10,y∈R10. In this example, the population size and iteration times are set as follows: Nu=400, Tu=50, Nl=40, Tl=20, and T=60:
(4.5)minxF(x,y)=((1+r-cos(απx1))+∑j=2K(xj-j-12)2+τ∑i=2K(yi-xi)2-γcos(γπx12y1),(1+r-sin(απx1))+∑j=2K(xj-j-12)2+τ∑i=2K(yi-xi)2-γsin(γπx12y1)∑j=2K(xj-j-12)2)minyf(x,y)=(y12+∑i=2K(yi-xi)2+∑i=2K10(1-cos(πk(yi-xi))),∑i=1K(yi-xi)2+∑i=2K10|1-sin(πk(yi-xi))|)s.t.-K≤yi≤K,(i=1,2,…,K);1≤x1≤4,-K≤xj≤K,(j=2,3,…,K).α=1,r=0.1,τ=1,γ=1,K=10.Example 4.4.
Example4.4 is taken from [28]. Here x∈R10,y∈R10. In this example, the population size and iteration times are set as follows: Nu=400,Tu=50, Nl=40, Tl=20, and T=80. 											(4.6)minxF(x,y)=v1(x1)+∑j=2Kyj2+10(1-cos(πk)yi)+τ∑i=2K(yi-xi)2-rcos(γπx12y1),v2(x1)+∑j=2Kyj2+10(1-cos(πk)yi)+τ∑i=2K(yi-xi)2-rsin(γπx12y1),minyf(x,y)=(y12+∑i=2K(yi-xi)2,∑i=1Ki(yi-xi)2)s.t.-K≤yi≤K,(i=1,2,…,K);0.001≤x1≤4,-K≤xj≤K,(j=2,3,…,K).α=1,r=0.25,τ=1,γ=4,K=10,
where
(4.7)v1(x1)={cos(0.2π)x1+sin(0.2π)|0.02sin(5πx1)|,for0≤x1≤1,x1-(1-sin(0.2π)),forx1>1,v2(x1)={-sin(0.2π)x1+cos(0.2π)|0.02sin(5πx1)|,for0≤x1≤1,0.1(x1-1)-sin(0.2π),forx1>1.
This problem is more difficult compared to the previous problems (Examples4.1 and 4.2) because the lower level problem of this example has multimodalities, thereby making the lower level problem difficult in finding the upper level Pareto optimal front. From Figure 5, it can be seen that the obtained Pareto front is very close to the theoretical Pareto optimal front, and the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front is 0.00027, that is, GD=0.00027 (see Table 2). Moreover, the lower SP value (SP=0.00127, see Table 2) shows that the proposed algorithm is able to obtain a good distribution of solutions on the entire range of the theoretical Pareto optimal front. Furthermore, two obtained lower level Pareto optimal fronts are given when x1=2 and x1=2.5.Figure 5
The obtained Pareto front of Example4.3.Figure6 shows the obtained Pareto front of Example 4.4 by the proposed algorithm. The upper level problem has multimodalities, thereby causing an algorithm difficulty in finding the upper level Pareto optimal front. From Figure 6, it can be seen that the obtained Pareto front is very close to the theoretical Pareto optimal front, and the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front is 0.00036, that is, GD=0.00036 (see Table 2). Moreover, the lower SP value (SP=0.00235, see Table 2) shows that the proposed algorithm is able to obtain a good distribution of solutions on the entire range of the theoretical Pareto optimal front. Furthermore, all corresponding lower level Pareto optimal fronts are given.Figure 6
The obtained Pareto front of Example4.4.Example 4.5.
Example4.5 is taken from [28]. Here x∈R10,y∈R10. In this example, the population size and iteration times are set as follows: Nu=400,Tu=50, Nl=40, Tl=20, and T=60:
(4.8)minxF(x·y)=(x1+∑j=3K(xj-j2)2+τ∑i=3K(yi-xi)2-cos(4tan-1(x2-y2x1-y1)),x2+∑j=3K(xj-j2)2+τ∑i=3K(yi-xi)2-cos(4tan-1(x2-y2x1-y1)))s.t.G(x)=x2-(1-x12)2≥0minyf(x,y)=(y1+∑i=3K(yi-xi)2,y2+∑i=3K(yi-xi)2)s.t.g1(x,y)=(y1-x1)2+(y2-x2)2≤r2R(x1)=(0.1+0.15|sin(2π(x1-0.1))|),-K≤yi≤K,(i=1,2,…,K);0≤x1≤K,-K≤xj≤K,(j=2,3,…,K).τ=1,r=0.2,K=10.Example 4.6.
Example4.6 is taken from [28]. Here x∈R1,y∈R9. In this example, the population size and iteration times are set as follows: Nu=400,Tu=50, Nl=40, Tl=20, and T=40:
(4.9)minxF(x,y)=((1-y1)(1+∑j=2Kyj2)x1,y1(1+∑j=2Kyj2)x1)s.t.G1(x,y)=-(1-y1)x1-12x1y1≤1minyf(x,y)=((1-y1)(1+∑j=K+1K+Lyj2)x1,y1(1+∑j=K+1K+Lyj2)x1)s.t.1≤x1≤2,-1≤y1≤1,-(K+L)≤yj≤(K+L),(j=2,3,…,K+L),K=5,L=4.
Figure7 shows the obtained Pareto front of Example 4.5 by the proposed algorithm. From Figure 7, it can be seen that the obtained Pareto front is very close to the theoretical Pareto optimal front, and the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front is 0.00058, that is, GD=0.00058 (see Table 2). Moreover, the lower SP value (SP=0.00364, see Table 2) shows that the proposed algorithm is able to obtain a good distribution of solutions on the entire range of the theoretical Pareto optimal front. Furthermore, all obtained lower level Pareto optimal fronts are given. It is also obvious that the Pareto optimal fronts for both the lower and upper level lie on constraint boundaries and every lower level Pareto optimal front has an unequal contribution to the upper level Pareto optimal front.
Figure8 shows the obtained Pareto front of Example 4.6 by the proposed algorithm. From Figure 8, it can be seen that the obtained Pareto front is very close to the theoretical Pareto optimal front, and the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front is 0.00039, that is, GD=0.00039 (see Table 2). Moreover, the lower SP value (SP=0.00168, see Table 2) shows that the proposed algorithm is able to obtain a good distribution of solutions on the entire range of the theoretical Pareto optimal front. Furthermore, three obtained lower level Pareto optimal fronts are given when y1=1, y1=1.5 and y1=2. It can be seen that only one Pareto optimal point from each participating lower level problem qualifies to be on the upper level Pareto optimal front.Figure 7
The obtained Pareto front of Example4.5.Figure 8
The obtained Pareto front of Example4.6.


### 4.2.3. The BLMPPs with Unknown Theoretical Pareto Optimal Fronts

Example 4.7.
Example4.7 is taken from [52], in which the theoretical Pareto optimal front is not given. Here x∈R2,y∈R3. In this example, the population size and iteration times are set as follows: Nu=100,Tu=50, Nl=20, Tl=10, and T=40:
(4.10)maxxF(x,y)=(x1+9x2+10y1+y2+3x3,9x1+2x2+2y1+7y2+4x3)s.t.G1(x,y)=3x1+9x2+9y1+5y2+3y3≤1039G2(x,y)=-4x1-x2+3y1-3y2+2y3≤94minyf(x,y)=(4x1+6x2+7y1+4y2+8y3,6x1+4x2+8y1+7y2+4y3)s.t.g1(x,y)=3x1-9x2-9y1-4y2≤61g2(x,y)=5x1+9x2+10y1-y2-2y3≤924g3(x,y)=3x1-3x2+y2+5y3≤420x1,x2,y1,y2,y3≥0.Example 4.8.
Example4.8 is taken from [23]. Here x∈R1,y∈R2. In this example, the population size and iteration times are set as follows: Nu=800,Tu=50,Nl=40, Tl=20, and T=40:
(4.11)minxF(x,y)=(x+y12+y22+sin2(x1+y),cos(y2)(0.1+y)(exp(-y10.1+y2)))s.t.G1(x,y)=(x-5)2-(y1-0.5)2-(y2-5)2≤16minyf(x,y)=((y1-2)2+(y2-2)24+xy2+(5-x)216+sin(y210),y12+(y2-6)2-2xy1-(5-x)280)s.t.g1(x,y)=y12-y2≤0g2(x,y)=5y12+y2-10-≤0g3(x,y)=y2-5-x6≤00≤x≤10,0≤y1,y2≤10.
Figure9 shows the obtained Pareto optimal front of Example 4.7 by the proposed algorithm. Note that, Zhang et al. [52] only obtained a single optimal solution x=(146.2955,28.9394) and y=(0,67.9318,0) which lies on the maximum of the F2 using weighted sum method. In contrast, a set of Pareto optimal solutions are obtained by the proposed algorithm. However, the fact that the single optimal solution in [49] is included in the obtained Pareto optimal solutions illustrates the feasibility of proposed algorithm. Figure 10 shows the final archive solutions of the Example 4.8 by the proposed algorithm. For this problem, the exact Pareto optimal front is not known, but the obtained Pareto optimal front by the proposed algorithm is similar to that reported in the previous study [23].Figure 9
The obtained front of Example4.7.Figure 10
The obtained front of Example4.8.


## 4.2.1. Low Dimension BLMPPs

Example 4.1.
Example4.1 is taken from [22]. Here x∈R1,y∈R2. In this example, the population size and iteration times are set as follows: Nu=200, Tu=200, Nl=40, Tl=40, and T=40:
(4.3)minxF(x,y)=(y1-x,y2)s.t.G1(y)=1+y1+y2≥0minyf(x,y)=(y1,y2)s.t.g1(x,y)=x2-y12-y22≥0,-1≤y1,y2≤1,0≤x≤1.
Figure1 shows the obtained Pareto front of this example by the proposed algorithm. From Figure 1, it can be seen that the obtained Pareto front is very close to the theoretical Pareto optimal front, and the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front is 0.00024, that is, GD=0.00024 (see Table 2). Moreover, the lower SP value (SP=0.0042, see Table 2) shows that the proposed algorithm is able to obtain a good distribution of solutions on the entire range of the theoretical Pareto optimal front. Figure 2 shows the obtained solutions of this example, which follow the relationship, that is, y1=-1-y2, y2=-1/2±(1/4)8x2-4 and x∈(1/2,1). It is also obvious that all obtained solutions are close to being on the upper level constraint G(x) boundary (1+y1+y2=0).Table 2
Results of the Generation Distance (GD) and Spacing (SP) metrics for the above six examples.
Example
GD
SP
Example4.1
0.00024
0.00442
Example4.2
0.00003
0.00169
Example4.3
0.00027
0.00127
Example4.4
0.00036
0.00235
Example4.5
0.00058
0.00364
Example4.6
0.00039
0.00168Figure 1
The obtained Pareto optimal front of Example4.1.Figure 2
The obtained solutions of Example4.1.Example 4.2.
Example4.2 is taken from [51]. Here x∈R1,y∈R2. In this example, the population size and iteration times are set as follows: Nu=200, Tu=50, Nl=40, Tl=20, and T=40:
(4.4)minxF(x,y)=(x2+(y1-1)2+y22,(x-1)2+(y1-1)2+y22)minyf(x,y)=(y12+y22,(y1-x)2+y22)-1≤x,y1,y2≤2.
Figure3 shows the obtained Pareto optimal front of this example by the proposed algorithm. From Figure 3, it is obvious that the obtained Pareto optimal front is very close to the theoretical Pareto optimal front, the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front is 0.00003 (see Table 2). On the other hand, the obtained Pareto optimal solutions can distribute uniformly on entire range of theoretical Pareto optimal front base on the fact that the SP value is lower (SP=0.00169, see Table 2). Figure 4 shows the obtained Pareto optimal solutions, they follow the relationship, that is, x=y1,y1∈[0.5,1] and y2=0.Figure 3
The obtained Pareto optimal front of Example4.2.Figure 4
The obtained solutions Example4.2.


## 4.2.2. High Dimension BLMPPs

Example 4.3.
Example4.3 is taken from [28]. Here x∈R10,y∈R10. In this example, the population size and iteration times are set as follows: Nu=400, Tu=50, Nl=40, Tl=20, and T=60:
(4.5)minxF(x,y)=((1+r-cos(απx1))+∑j=2K(xj-j-12)2+τ∑i=2K(yi-xi)2-γcos(γπx12y1),(1+r-sin(απx1))+∑j=2K(xj-j-12)2+τ∑i=2K(yi-xi)2-γsin(γπx12y1)∑j=2K(xj-j-12)2)minyf(x,y)=(y12+∑i=2K(yi-xi)2+∑i=2K10(1-cos(πk(yi-xi))),∑i=1K(yi-xi)2+∑i=2K10|1-sin(πk(yi-xi))|)s.t.-K≤yi≤K,(i=1,2,…,K);1≤x1≤4,-K≤xj≤K,(j=2,3,…,K).α=1,r=0.1,τ=1,γ=1,K=10.Example 4.4.
Example4.4 is taken from [28]. Here x∈R10,y∈R10. In this example, the population size and iteration times are set as follows: Nu=400,Tu=50, Nl=40, Tl=20, and T=80. 											(4.6)minxF(x,y)=v1(x1)+∑j=2Kyj2+10(1-cos(πk)yi)+τ∑i=2K(yi-xi)2-rcos(γπx12y1),v2(x1)+∑j=2Kyj2+10(1-cos(πk)yi)+τ∑i=2K(yi-xi)2-rsin(γπx12y1),minyf(x,y)=(y12+∑i=2K(yi-xi)2,∑i=1Ki(yi-xi)2)s.t.-K≤yi≤K,(i=1,2,…,K);0.001≤x1≤4,-K≤xj≤K,(j=2,3,…,K).α=1,r=0.25,τ=1,γ=4,K=10,
where
(4.7)v1(x1)={cos(0.2π)x1+sin(0.2π)|0.02sin(5πx1)|,for0≤x1≤1,x1-(1-sin(0.2π)),forx1>1,v2(x1)={-sin(0.2π)x1+cos(0.2π)|0.02sin(5πx1)|,for0≤x1≤1,0.1(x1-1)-sin(0.2π),forx1>1.
This problem is more difficult compared to the previous problems (Examples4.1 and 4.2) because the lower level problem of this example has multimodalities, thereby making the lower level problem difficult in finding the upper level Pareto optimal front. From Figure 5, it can be seen that the obtained Pareto front is very close to the theoretical Pareto optimal front, and the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front is 0.00027, that is, GD=0.00027 (see Table 2). Moreover, the lower SP value (SP=0.00127, see Table 2) shows that the proposed algorithm is able to obtain a good distribution of solutions on the entire range of the theoretical Pareto optimal front. Furthermore, two obtained lower level Pareto optimal fronts are given when x1=2 and x1=2.5.Figure 5
The obtained Pareto front of Example4.3.Figure6 shows the obtained Pareto front of Example 4.4 by the proposed algorithm. The upper level problem has multimodalities, thereby causing an algorithm difficulty in finding the upper level Pareto optimal front. From Figure 6, it can be seen that the obtained Pareto front is very close to the theoretical Pareto optimal front, and the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front is 0.00036, that is, GD=0.00036 (see Table 2). Moreover, the lower SP value (SP=0.00235, see Table 2) shows that the proposed algorithm is able to obtain a good distribution of solutions on the entire range of the theoretical Pareto optimal front. Furthermore, all corresponding lower level Pareto optimal fronts are given.Figure 6
The obtained Pareto front of Example4.4.Example 4.5.
Example4.5 is taken from [28]. Here x∈R10,y∈R10. In this example, the population size and iteration times are set as follows: Nu=400,Tu=50, Nl=40, Tl=20, and T=60:
(4.8)minxF(x·y)=(x1+∑j=3K(xj-j2)2+τ∑i=3K(yi-xi)2-cos(4tan-1(x2-y2x1-y1)),x2+∑j=3K(xj-j2)2+τ∑i=3K(yi-xi)2-cos(4tan-1(x2-y2x1-y1)))s.t.G(x)=x2-(1-x12)2≥0minyf(x,y)=(y1+∑i=3K(yi-xi)2,y2+∑i=3K(yi-xi)2)s.t.g1(x,y)=(y1-x1)2+(y2-x2)2≤r2R(x1)=(0.1+0.15|sin(2π(x1-0.1))|),-K≤yi≤K,(i=1,2,…,K);0≤x1≤K,-K≤xj≤K,(j=2,3,…,K).τ=1,r=0.2,K=10.Example 4.6.
Example4.6 is taken from [28]. Here x∈R1,y∈R9. In this example, the population size and iteration times are set as follows: Nu=400,Tu=50, Nl=40, Tl=20, and T=40:
(4.9)minxF(x,y)=((1-y1)(1+∑j=2Kyj2)x1,y1(1+∑j=2Kyj2)x1)s.t.G1(x,y)=-(1-y1)x1-12x1y1≤1minyf(x,y)=((1-y1)(1+∑j=K+1K+Lyj2)x1,y1(1+∑j=K+1K+Lyj2)x1)s.t.1≤x1≤2,-1≤y1≤1,-(K+L)≤yj≤(K+L),(j=2,3,…,K+L),K=5,L=4.
Figure7 shows the obtained Pareto front of Example 4.5 by the proposed algorithm. From Figure 7, it can be seen that the obtained Pareto front is very close to the theoretical Pareto optimal front, and the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front is 0.00058, that is, GD=0.00058 (see Table 2). Moreover, the lower SP value (SP=0.00364, see Table 2) shows that the proposed algorithm is able to obtain a good distribution of solutions on the entire range of the theoretical Pareto optimal front. Furthermore, all obtained lower level Pareto optimal fronts are given. It is also obvious that the Pareto optimal fronts for both the lower and upper level lie on constraint boundaries and every lower level Pareto optimal front has an unequal contribution to the upper level Pareto optimal front.
Figure8 shows the obtained Pareto front of Example 4.6 by the proposed algorithm. From Figure 8, it can be seen that the obtained Pareto front is very close to the theoretical Pareto optimal front, and the average distance between the obtained Pareto optimal front and the theoretical Pareto optimal front is 0.00039, that is, GD=0.00039 (see Table 2). Moreover, the lower SP value (SP=0.00168, see Table 2) shows that the proposed algorithm is able to obtain a good distribution of solutions on the entire range of the theoretical Pareto optimal front. Furthermore, three obtained lower level Pareto optimal fronts are given when y1=1, y1=1.5 and y1=2. It can be seen that only one Pareto optimal point from each participating lower level problem qualifies to be on the upper level Pareto optimal front.Figure 7
The obtained Pareto front of Example4.5.Figure 8
The obtained Pareto front of Example4.6.


## 4.2.3. The BLMPPs with Unknown Theoretical Pareto Optimal Fronts

Example 4.7.
Example4.7 is taken from [52], in which the theoretical Pareto optimal front is not given. Here x∈R2,y∈R3. In this example, the population size and iteration times are set as follows: Nu=100,Tu=50, Nl=20, Tl=10, and T=40:
(4.10)maxxF(x,y)=(x1+9x2+10y1+y2+3x3,9x1+2x2+2y1+7y2+4x3)s.t.G1(x,y)=3x1+9x2+9y1+5y2+3y3≤1039G2(x,y)=-4x1-x2+3y1-3y2+2y3≤94minyf(x,y)=(4x1+6x2+7y1+4y2+8y3,6x1+4x2+8y1+7y2+4y3)s.t.g1(x,y)=3x1-9x2-9y1-4y2≤61g2(x,y)=5x1+9x2+10y1-y2-2y3≤924g3(x,y)=3x1-3x2+y2+5y3≤420x1,x2,y1,y2,y3≥0.Example 4.8.
Example4.8 is taken from [23]. Here x∈R1,y∈R2. In this example, the population size and iteration times are set as follows: Nu=800,Tu=50,Nl=40, Tl=20, and T=40:
(4.11)minxF(x,y)=(x+y12+y22+sin2(x1+y),cos(y2)(0.1+y)(exp(-y10.1+y2)))s.t.G1(x,y)=(x-5)2-(y1-0.5)2-(y2-5)2≤16minyf(x,y)=((y1-2)2+(y2-2)24+xy2+(5-x)216+sin(y210),y12+(y2-6)2-2xy1-(5-x)280)s.t.g1(x,y)=y12-y2≤0g2(x,y)=5y12+y2-10-≤0g3(x,y)=y2-5-x6≤00≤x≤10,0≤y1,y2≤10.
Figure9 shows the obtained Pareto optimal front of Example 4.7 by the proposed algorithm. Note that, Zhang et al. [52] only obtained a single optimal solution x=(146.2955,28.9394) and y=(0,67.9318,0) which lies on the maximum of the F2 using weighted sum method. In contrast, a set of Pareto optimal solutions are obtained by the proposed algorithm. However, the fact that the single optimal solution in [49] is included in the obtained Pareto optimal solutions illustrates the feasibility of proposed algorithm. Figure 10 shows the final archive solutions of the Example 4.8 by the proposed algorithm. For this problem, the exact Pareto optimal front is not known, but the obtained Pareto optimal front by the proposed algorithm is similar to that reported in the previous study [23].Figure 9
The obtained front of Example4.7.Figure 10
The obtained front of Example4.8.


## 5. Conclusion

In this paper, an EQPSO is presented, in which an elite strategy is exerted for global best particle to prevent the swarm from clustering, enabling the particle to escape the local optima. The EQPSO algorithm is employed for solving bilevel multiobjective programming problem (BLMPP) for the first time. In this study, some numerical examples are used to explore the feasibility and efficiency of the proposed algorithm. The experimental results indicate that the obtained Pareto front by the proposed algorithm is very close to the theoretical Pareto optimal front, and the solutions are also distributed uniformly on entire range of the theoretical Pareto optimal front. The proposed algorithm is simple and easy to implement, which provides another appealing method for further study on BLMPP.



---

*Source: 102482-2012-11-29.xml*