An Improved Particle Swarm Optimization for Solving Bilevel Multiobjective Programming Problem

An improved particle swarm optimization PSO algorithm is proposed for solving bilevel multiobjective programming problem BLMPP . For such problems, the proposed algorithm directly simulates the decision process of bilevel programming, which is different from most traditional algorithms designed for specific versions or based on specific assumptions. The BLMPP is transformed to solvemultiobjective optimization problems in the upper level and the lower level interactively by an improved PSO. And a set of approximate Pareto optimal solutions for BLMPP is obtained using the elite strategy. This interactive procedure is repeated until the accurate Pareto optimal solutions of the original problem are found. Finally, some numerical examples are given to illustrate the feasibility of the proposed algorithm.


Introduction
Bilevel programming problem BLPP arises in a wide variety of scientific and engineering applications including optimal control, process optimization, game-playing strategy development, and transportation problem 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 EA have been employed to address BLPP in papers 12-16 .However, the bilevel multiobjective programming problem BLMPP has seldom been studied.Shi and Xia 17, 18 , Abo-Sinna and Baky 19 , Nishizaki and Sakawa 20 , and

Problem Formulation
R p , and g : R n 1 × R n 2 → R q .The general model of the BLMPP can be written as follows: 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.Let S { 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 fixed x ∈ 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 feedback to the upper level problem, and this point of view is accepted usually.Based on this fact, the PSO 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 PSO algorithm uses a group of points in its operation thus, the PSO can be developed as a new way for solving BLMPP.In the following, we present an improved PSO algorithm for solving problem 2.1 .

The Algorithm
The process of the proposed algorithm is an interactive coevolutionary process for both the upper level and the lower level.We first initialize population and then solve multiobjective optimization problems in the upper level and the lower level interactively using an improved PSO.Afterwards, a set of approximate Pareto optimal solutions for problem 1 is obtained by the elite strategy which was adopted in Deb et al. 34 .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. Initialize.Substep 1.1.Initialize the population P 0 with N u particles which is composed by n s N u /N l subswarms of size N l each.The particle's position of the kth k 1, 2, . . ., n s subswarm is presented as z j x j , y j j 1, 2, . . ., n l , and the corresponding velocity is presented as: v j v x j , v y j j 1, 2, . . ., n l , z j and v j are sampled randomly in the feasible space, respectively.Substep 1.2.Initialize the external loop counter t : 0.
Step 2. For the kth subswarm k 1, 2, . . ., n s , each particle is assigned a nondomination rank ND l and a crowding value CD l in f space.Then, all resulting subswarms are combined into one population which is named as the P t .Afterwards, each particle is assigned a nondomination rank ND u and a crowding value CD u in F space.
Step 3. The nondomination particles assigned both ND u 1 and ND l 1 from P t are saved in the elite set A t .
Step 4. For the kth subswarm k 1, 2, . . ., n s , update the lower level decision variables.Substep 4.1.Initialize the lower level loop counter t l : 0. Substep 4.2.Update the jth j 1, 2, . . ., N l particle's position and velocity with the fixed x j and the fixed v j using 3 Substep 4.5.Each particle of the ith subswarm is reassigned a nondomination rank ND l and a crowding value CD l in F space.Then, all resulting subswarms are combined into one population which is renamed as the Q t .Afterwards, each particle is reassigned a nondomination rank ND u and a crowding value CD u in F space.
Step 5. Combine population P t and Q t to form R t .The combined population R t is reassigned a nondomination rank ND u , and the particles within an identical nondomination rank are assigned a crowding distance value CD u in the F space.
Step 6. Choose half particles from R t .The particles of rank ND u 1 are considered first.From the particles of rank ND u 1, the particles with ND l 1 are noted one by one in the order of reducing crowding distance CD u , for each such particle the corresponding subswarm from its source population either P t or Q t is copied in an intermediate population S t .If a subswarm is already copied in S t and a future particle from the same subswarm is found to have ND u ND l 1, the subswarm is not copied again.When all particles of ND u 1 are considered, a similar consideration is continued with ND u 2 and so on till exactly n s subswarms are copied in S t .
Step 7. Update the elite set A t .The nondomination particles assigned both ND u 1 and ND l 1 from S t are saved in the elite set A t .
Step 8. Update the upper level decision variables in S t .Substep 8.1.Initiate the upper level loop counter t u : 0. Substep 8.2.Update the ith i 1, 2, . . ., N u particle's position and velocity with the fixed y i and the fixed v i using Step 9. Consider t : t 1.
Step 10.If t ≥ T , output the elite set A t .Otherwise, go to Step 2.
In Steps 4 and 8, the global best position is chosen at random from the elite set A t .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.

Numerical Examples
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.

Generational Distance (GD)
This metric used by Deb 35 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: where n is the number of the obtained Pareto optimal solutions by the proposed algorithm and d i is the Euclidean distance between each obtained Pareto optimal solution and the nearest member of the theoretical Pareto optimal set.x i The ith particle's position of the upper level problem.

v xi
The velocity of x i .

y j
The jth particle's position of the lower level problem.

v yj
The velocity of y j .

z j
The jth particle's position of BLMPP.

p pbest yj
The jth particle's personal best position for the lower level problem.

p pbest xi
The ith particle's personal best position for the upper level problem.

p gbest l
The particle's global best position for the lower level problem.

p gbest u
The particle's global best position for the upper level problem.

N u
The population size of the upper level problem.

N l
The subswarm size of the lower level problem.
t Current iteration number for the overall problem.

T
The predefined max iteration number for t.

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 for t u .

T l
The predefined max iteration number for t l .
w u Inertia weights for the upper level problem.

w l
Inertia weights the lower level problem.

c 1u
The cognitive learning rate for the upper level problem.

c 2u
The social learning rate for the upper level problem.

c 1l
The cognitive learning rate for the lower level problem.

c 2l
The social learning rate for the lower level problem.

ND u
Nondomination sorting rank of the upper level problem.

CD u
Crowding distance value of the upper level problem.

ND l
Nondomination sorting rank of the lower level problem.

CD l
Crowding distance value of the lower level problem.

P t
The tth iteration population.

Q t
The offspring of P t .

S t
Intermediate population.

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 35 : where number of the upper level objective function, n is the number of the obtained solutions by the proposed algorithm.The PSO parameters are set as follows: r 1u , r 2u , r 1l , r 2l ∈ random 0, 1 , the inertia weight w u w l 0.7298, and acceleration coefficients with c 1u c 2u c 1l c 2l 1.49618.All results presented in this paper have been obtained on a personal computer CPU: AMD Phenom II X6 1055T 2.80 GHz; RAM: 3.25 GB using a C# implementation of the proposed algorithm, and the figures were obtained using the origin 8.0.
Example 4.1.Example 4.1 is taken from 22 .Here x ∈ R 1 , y ∈ R 2 .In this example, the population size and iteration times are set as follows: N u 200, T u 200, N l 40, T l 40, and T 40:

4.3
Figure 1 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.00026, that is, GD 0.00026 see Table 2 .Moreover, the lower SP value SP 0.17569, 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, It is also obvious that all obtained solutions are close to being on the upper level constraint G x boundary 1 y 1 y 2 0 .Example 4.2.Example 4.2 is taken from 36 .Here x ∈ R 1 , y ∈ R 2 .In this example, the population size and iteration times are set as follows: N u 200, T u 50, N l 40, T l 20, and T 40.

4.4
Figure 3 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.00004 see Table 2 .On the other hand, the The obtained Pareto optimal front The theoretical Pareto optimal front obtained Pareto optimal solutions can be distributed uniformly on entire range of theoretical Pareto optimal front based on the fact that the SP value is lower SP 0.00173, see Table 2 .
Figure 4 shows the obtained Pareto optimal solutions; they follow the relationship, that is, x y 1 , y 1 ∈ 0.5, 1 and y 2 0.

4.5
Figure 5 shows the obtained Pareto optimal front of Example 4.3 by the proposed algorithm.
Figure 6 shows all five constrains for all obtained Pareto optimal solutions and it can be seen that the G 1 , g 2 and g 3 are active constrains.Note that, Zhang et al. 37 only obtained a single optimal solution x 146.2955, 28.9394 , and y 0, 67.9318, 0 which lies on the maximum of the F 2 using weighted sum method.In contrast, a set of Pareto optimal solutions is obtained by the proposed algorithm.However, the fact that the single optimal solution in 37 is included in the obtained Pareto optimal solutions illustrates the feasibility of proposed algorithm.

Conclusion
In this paper, an improved PSO is presented for BLMPP.The BLMPP is transformed to solve the multiobjective optimization problems in the upper level and the lower level interactively using the proposed algorithm for a predefined count.And a set of accurate Pareto optimal solutions for BLMPP is obtained by the elite strategy.The experimental results illustrate 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.Furthermore, the proposed algorithm is simple and easy to implement.It also provides another appealing method for further study on BLMPP.
obtained Pareto optimal front by the proposed algorithm The single optimal solution in[37]
Substep 8.3.Consider t u : t u 1. Substep 8.4.If t u ≥ T u , go to Substep 8.5.Otherwise, go to Substep 8.2.Every member is then assigned a nondomination rank ND u and a crowding distance value CD u in F space.

Table 1 :
The notations of the algorithm.

Table 2 :
Results of the Generation Distance GD and Spacing SP metrics for Examples 4.1 and 4.2.