Chaotic Multi-Objective Particle Swarm Optimization Algorithm Incorporating Clone Immunity

Abstract

It is generally known that the balance between convergence and diversity is a key issue for solving multi-objective optimization problems. Thus, a chaotic multi-objective particle swarm optimization approach incorporating clone immunity (CICMOPSO) is proposed in this paper. First, points in a non-dominated solution set are mapped to a parallel-cell coordinate system. Then, the status of the particles is evaluated by the Pareto entropy and difference entropy. At the same time, the algorithm parameters are adjusted by feedback information. At the late stage of the algorithm, the local-search ability of the particle swarm still needs to be improved. Logistic mapping and the neighboring immune operator are used to maintain and change the external archive. Experimental test results show that the convergence and diversity of the algorithm are improved.

1. Introduction

Most problems [1,2,3,4,5,6] in engineering and science are multi-objective optimization problems, and these objectives usually conflict with each other. The focus of academics and engineers is on finding an optimal solution for these problems. For multi-objective optimization problems, however, when a certain objective is improved, it often simultaneously impacts other objectives. Usually, people make a trade-off and compromise to solve this problem. A very effective method to solve such problems is by using an intelligent algorithm. In recent years, intelligent algorithms have been paid increasing attention in many areas [7,8,9,10] and have achieved good results. However, when solving multi-objective optimization problems using intelligent algorithms, there are three goals that need to be satisfied: (1) particles must be as close to the Pareto optimal front as possible; (2) there must be a maximal number of particles in the Pareto optimal front; and (3) particles must be distributed as evenly as possible in the Pareto optimal front.

The particle swarm optimization (PSO) algorithm [11,12] is one of the most important and studied paradigms in computational swarm intelligence. It was put forward by Dr. Kennedy and Professor Eberhart in 1995. The simple form, simple parameters, and rapid convergence of PSO have contributed to its rapid development and it becoming a research hotspot over the last 20 years. For multi-objective PSO algorithms, there are six research directions [13]: (1) aggregating approaches, which combine all the objectives of the problem into a single one. In other words, the multi-objective problem is transformed into a single-objective problem. It is not a new idea, since aggregating functions can be derived from the well-known K-T conditions for non-dominated solutions. (2) Lexicographic ordering, in which the objectives are ranked in order of importance. The optimal solution is obtained by separately minimizing the objective functions, starting with the most important one and proceeding according to the order of importance of the objectives. Lexicographic ordering tends to be useful only when few objective functions are used, and it maybe sensitive to the ordering of the objectives. (3) Subpopulation approaches, which use several subpopulations as single-objective optimizers. In order to balance the objectives, the subpopulations exchange information with each other. Approaches: The information exchange with Sub-Population is not well controlled, and it is difficult to ensure concurrent evolution of each objective. (4) Pareto-based approaches, which use leader selection techniques based on Pareto dominance. This method is the mainstream method, such as multi-objective particle swarm (MOPSO) [14]. It is useful for the multi-objective optimization problems with few objective functions, but when the objective functions increase, the selection pressure increases. (5) Combined approaches, that combine the PSO algorithm with other algorithms. Some examples include the genetic algorithm (GA) [15], the cultural algorithm [16], and the differential evolution (DE) algorithm [17]. (6) Other approaches [18,19] that do not fall into the above five types, such as designing a threshold for the particle updating strategy [20].

The remainder of this paper is organized as follows. In Section 2, we describe the multi-objective optimization and PSO algorithm. Thereafter, in Section 3, we explain a computational method to improve the multi-objective PSO algorithm. Section 4 outlines the strategy and flow of this algorithm. Test problems, performance measures, and results are provided in Section 5, and conclusions are presented in Section 6.

2. Propaedeutics

2.1. Multi-Objective Optimization Problem

In general, a multi-objective optimization problem with n decision variables and M objective functions can be described as [21]:

$$\begin{array}{cc}min\hfill & y=f\left(x\right)=[f_1\left(x\right),f_2\left(x\right),\cdots ,f_M\left(x\right)]\hfill \\ s.t.\hfill & \begin{array}{cccc}{g}_i\left(x\right)\hfill & \le & \hfill 0,& \hfill i=1,2,\cdots ,p\\ h_j\left(x\right)\hfill & =& \hfill 0,& \hfill j=1,2,\cdots ,q\end{array}\hfill \end{array}$$

(1)

where, $x=(x_1,x_2,\cdots ,x_n)\in X$ is a single decision variable, X is the decision variable space, $y=(f_1,f_2,\cdots ,f_M)\in Y$ is the objective function vector, Y is the objective function space, $g_i\left(x\right)$ is the $i^{\mathrm{th}}$ inequality constraint, and $h_j\left(x\right)$ is the $j^{\mathrm{th}}$ equality constraint.

The following definitions are useful for the conceptual framework of a multi-objective optimization problem:

Definition 1.

(Pareto dominance) Vector $x^{\prime }=(x_1^{\prime },x_2^{\prime },\cdots ,x_n^{\prime })$ dominates vector $x=(x_1,x_2,\cdots ,x_n)$ if and only if the next statement is true:

$$\forall i=(1,2,\cdots ,n):f_i\left(x^{\prime }\right)\le f_i\left(x\right)and\exists i\in (1,2,\cdots ,n):f_i\left(x^{\prime }\right)<f_i\left(x\right).$$

We denote this dominance as $x^{\prime }\prec x$.

Definition 2.

(Pareto optimality) A solution $x^{*}\in D$ is a Pareto optimal solution if there is no other $x\in D$ that satisfies $f\left(x\right)\prec f\left(x^{*}\right)$.

Definition 3.

(Pareto optimal set) The Pareto optimal set is defined as the set of all Pareto optimal solutions.

Definition 4.

(Pareto optimal front) The Pareto optimal front consists of the values of the objective functions at the solutions in the Pareto optimal set.

2.2. Particle Swarm Optimization Algorithm

Let n be the dimension of the search space, $x_i=(x_{i1},x_{i2},\cdots ,x_{in})$ be the current position of the $i^{\mathrm{th}}$ particle in the swarm, $p_{i,\mathrm{best}}=(p_{i1,\mathrm{best}},p_{i2,\mathrm{best}},\cdots ,p_{in,\mathrm{best}})$ be the best position of the $i^{\mathrm{th}}$ particle at that time, and $g_{\mathrm{best}}=(g_{1,\mathrm{best}},g_{2,\mathrm{best}},\cdots ,g_{n,\mathrm{best}})$ be the best position that any particle in the entire swarm has visited. The rate of the velocity of the $i^{\mathrm{th}}$ particle is denoted as $v_i=(v_{i1},v_{i2},\cdots ,v_{in})$. The position of each particle $x_{ij}^t$, and its velocity $v_{ij}^t$ updated according to the following:

$$v_{ij}^{t+1}=w{v}_{ij}^t+c_1{r}_1(p_{ij,\mathrm{best}}^t-x_{ij}^t)+c_2{r}_2(g_{ij,\mathrm{best}}^t-x_{ij}^t)$$

(2)

$$x_{ij}^{t+1}=x_{ij}^t+v_{ij}^{t+1},$$

(3)

where $c_1,c_2$ are the learning factors, and $r_1,r_2$ are random numbers in $[0,1]$. w is the inertia weight, which is defined as follows [22]:

$$w=w_{max}-t*(w_{max}-w_{min})/T_{max},$$

(4)

where, $w_{max},w_{min}$ are the maximum and minimum of inertia weight w, respectively, t is the time of the present iteration, and $T_{max}$ is the maximum number of iterations.

3. Pareto Entropy and Difference Entropy-Based Improvements of PSO

3.1. Parallel Cell Coordinate System

The parallel cell coordinate system (PCCS) [23] maps the target vector of non-dominated solutions in the external archive to a two-dimensional plane by parallel coordinates, and then rounds these coordinate values. Its mathematical formula is as follows:

$$L_{k,m}=⌈K\frac{f_{k,m}-f_m^{min}}{f_m^{max}-f_m^{min}}⌉,$$

(5)

where $\lceil x\rceil$ is a top integral function; $k=1,2,\cdots ,K$ where K is the external archive size in the current iteration; $m=1,2,\cdots ,M$ where M is the number of objective functions in the optimization problem; $f_m^{max}=\underset{k}{max}{f}_{k,m}$, and $f_m^{min}=\underset{k}{min}{f}_{k,m}$. When $f_m^{min}=f_m^{max}$, $L_{k,m}=1$.

Consider

Table 1. Parallel grid coordinate system of particles.

tab	F1	F2	F3	${\mathit{L}}_{\mathit{k},1}$	${\mathit{L}}_{\mathit{k},2}$	${\mathit{L}}_{\mathit{k},3}$
1	0.5377	1.8339	−2.2588	4	8	1
2	0.8622	0.3188	−1.3077	5	5	2
3	−0.4336	0.3426	3.5784	2	5	8
4	2.7694	−1.3499	3.0349	8	1	8
5	0.725	−0.0631	0.7147	4	4	5
6	−0.2050	−0.1241	1.4897	3	4	6
7	1.4090	1.4172	0.6715	6	7	5
8	−1.2075	0.7172	1.6302	1	6	6

as an example, in which eight particles were included for three objectives. According to the formula (5), we can calculate $L_{k,m}$ and draw the Figure 1.

In this paper, we use the Pareto entropy to estimate the distribution uniformity of the Pareto front, which is calculated with the following formula:

$$Entropy\left(t\right)=-\sum _{k=1}^K\sum _{m=1}^M\frac{Cel{l}_{k,m}\left(t\right)}{KM}log\frac{Cel{l}_{k,m}\left(t\right)}{KM}$$

(6)

When approximating the Pareto front map to the PCCS, $Cel{l}_{k,m}\left(t\right)$ is the number of cell coordinate components in the kth row and mth column cell.

According to Pareto entropy, we know the distribution of particles in the current iteration. However, in order to judge the change between this iteration and the last iteration, a difference entropy $\Delta Entropy$ is proposed [23]:

$$\Delta Entropy=Entropy\left(t\right)-Entropy(t-1)$$

(7)

Through the difference entropy, information about the external archive can be obtained, and we can dynamically adjust population using this information.

3.2. Difference Entropy Discussion

Previous work has proved that the maximum Pareto entropy is $logM$ and the minimum Pareto entropy is $logKM$ [23]. Here, we discuss a few special situations for the difference entropy:

(1) In the $t^{th}$ iteration, the corresponding coordinate components in the PCCS of all objective vectors on the Pareto front are the most nonuniform, and after maintaining the external archive, the corresponding coordinate components of all objective vectors are the most uniform.

Proof. 

Assume that, in the $t^{th}$ iteration, the distribution of parallel cell coordinates is $Cel{l}_{k,m}\left(t\right)$ and the external archive size is $K_1$. In the ${(t+1)}^{th}$ iteration, the distribution of parallel cell coordinates is $Cel{l}_{k,m}(t+1)$, and the external archive size is $K_2$. Next, make all components of the nth column collect in the nth column of c-row, so that $Cel{l}_{c,n}^{(}t)=K_1$ and $Cel{l}_{k\ne c,n}\left(t\right)=0$. From the L’Hospital’s rule, we can know ${lim}_{x\to 0}xlog\left(x\right)=0$, so we stipulate that $0log\left(0\right)=0$. In the ${(t+1)}^{th}$ iteration, $Cel{l}_{k,m}(t+1)=1$, that is to say, it always holds for kth row and mth column. We obtain Difference Entropy through the following formula:

$$\begin{array}{lll}\hfill \Delta Entropy& =& Entropy(t+1)-Entropy\left(t\right)\hfill \\ & =& -{\displaystyle \sum _{k=1}^{K_2}}{\displaystyle \sum _{m=1}^M}\frac{Cel{l}_{k,m}(t+1)}{K_{2}M}log\frac{Cel{l}_{k,m}(t+1)}{K_{2}M}-\left(-{\displaystyle \sum _{k=1}^{K_1}}{\displaystyle \sum _{m=1}^M}\frac{Cel{l}_{km}\left(t\right)}{K_{1}M}log\frac{Cel{l}_{k,m}\left(t\right)}{K_{1}M}\right)\hfill \\ & =& {\displaystyle \sum _{k=1}^{K_1}}{\displaystyle \sum _{m=1}^M}\frac{Cel{l}_{k,m}\left(t\right)}{K_{1}M}log\frac{Cel{l}_{k,m}\left(t\right)}{K_{1}M}-{\displaystyle \sum _{k=1}^{K_2}}{\displaystyle \sum _{m=1}^M}\frac{Cel{l}_{k,m}(t+1)}{K_{2}M}log\frac{Cel{l}_{k,m}(t+1)}{K_{2}M}\hfill \\ & =& {\displaystyle \sum _{m=1}^M}\frac{K_1}{K_{1}M}log\frac{K_1}{K_{1}M}-{\displaystyle \sum _{k=1}^{K_2}}{\displaystyle \sum _{m=1}^M}\frac{1}{K_{2}M}log\frac{1}{K_{2}M}\hfill \\ & =& log\frac{1}{M}-log\frac{1}{K_{2}M}\hfill \\ & =& log{K}_2\hfill \end{array}$$

(8)

 □

(2) In the $t^{th}$ iteration, the corresponding coordinate components in the PCCS of all objective vectors on the Pareto front are the most nonuniform and after maintaining the external archive, the corresponding coordinate components of all objective vectors remain the most nonuniform.

Proof. 

From Case (1), we know $Cel{l}_{c_1,n_1}\left(t\right)=K_1$ and $Cel{l}_{k\ne c_1,n_1}\left(t\right)=0$. Furthermore, $Cel{l}_{c_2,n_2}(t+1)=K_2$ and $Cel{l}_{k\ne c_2,n_2}(t+1)=0$. From Equation (7):

$$\begin{array}{lll}\hfill \Delta Entropy& =& Entropy(t+1)-Entropy\left(t\right)\hfill \\ & =& -{\displaystyle \sum _{k=1}^{K_2}}{\displaystyle \sum _{m=1}^M}\frac{Cel{l}_{k,m}(t+1)}{K_{2}M}log\frac{Cel{l}_{k,m}(t+1)}{K_{2}M}-\left(-{\displaystyle \sum _{k=1}^{K_1}}{\displaystyle \sum _{m=1}^M}\frac{Cel{l}_{k,m}\left(t\right)}{K_{1}M}log\frac{Cel{l}_{km}\left(t\right)}{K_{1}M}\right)\hfill \\ & =& {\displaystyle \sum _{m=1}^M}\frac{K_1}{K_{1}M}log\frac{K_1}{K_{1}M}-{\displaystyle \sum _{m=1}^M}\frac{K_2}{K_{2}M}log\frac{K_2}{K_{2}M}\hfill \\ & =& log\frac{1}{M}-log\frac{1}{M}\hfill \\ & =& 0\hfill \end{array}$$

(9)

 □

(3) In the $t^{th}$ iteration, the corresponding coordinate components in the PCCS of all objective vectors on the Pareto front are the most uniform and after maintaining the external archive, the corresponding coordinate components of all objective vectors remain the most uniform.

Proof. 

From Cases (1) and (2), we know that $Cel{l}_{k,m}\left(t\right)=1$ and $Cel{l}_{k,m}(t+1)=1$. So:

$$\begin{array}{lll}\hfill \Delta Entropy& =& Entropy(t+1)-Entropy\left(t\right)\hfill \\ & =& -{\displaystyle \sum _{k=1}^{K_2}}{\displaystyle \sum _{m=1}^M}\frac{Cel{l}_{k,m}(t+1)}{K_{2}M}log\frac{Cel{l}_{k,m}(t+1)}{K_{2}M}-\left(-{\displaystyle \sum _{k=1}^{K_1}}{\displaystyle \sum _{m=1}^M}\frac{Cel{l}_{k,m}\left(t\right)}{K_{1}M}log\frac{Cel{l}_{km}\left(t\right)}{K_{1}M}\right)\hfill \\ & =& {\displaystyle \sum _{k=1}^{K_1}}{\displaystyle \sum _{m=1}^M}\frac{1}{K_{1}M}log\frac{1}{K_{1}M}-{\displaystyle \sum _{k=1}^{K_2}}{\displaystyle \sum _{m=1}^M}\frac{1}{K_{2}M}log\frac{1}{K_{2}M}\hfill \\ & =& log\frac{1}{K_{1}M}-log\frac{1}{K_{2}M}\hfill \\ & =& log\frac{K_2}{K_1}\hfill \end{array}$$

(10)

 □

(4) In the $t^{th}$ iteration, the corresponding coordinate components in the PCCS of all objective vectors on the Pareto front are the most uniform, but after maintaining the external archive, the corresponding coordinate components of all objective vectors are the most nonuniform.

Proof. 

From Case (1), we know $Cel{l}_{k,m}\left(t\right)=1$, $Cel{l}_{c_2,n_2}(t+1)=K_2$, and $Cel{l}_{k\ne c_2,n_2}(t+1)=0$.

$$\begin{array}{lll}\hfill \Delta Entropy& =& Entropy(t+1)-Entropy\left(t\right)\hfill \\ & =& -{\displaystyle \sum _{k=1}^{K_2}}{\displaystyle \sum _{m=1}^M}\frac{Cel{l}_{k,m}(t+1)}{K_{2}M}log\frac{Cel{l}_{k,m}(t+1)}{K_{2}M}-\left(-{\displaystyle \sum _{k=1}^{K_1}}{\displaystyle \sum _{m=1}^M}\frac{Cel{l}_{k,m}\left(t\right)}{K_{1}M}log\frac{Cel{l}_{km}\left(t\right)}{K_{1}M}\right)\hfill \\ & =& {\displaystyle \sum _{k=1}^{K_1}}{\displaystyle \sum _{m=1}^M}\frac{1}{K_{1}M}log\frac{1}{K_{1}M}-{\displaystyle \sum _{m=1}^M}\frac{K_2}{K_{2}M}log\frac{K_2}{K_{2}M}\hfill \\ & =& -log\frac{1}{K_{1}M}-log\frac{1}{M}\hfill \\ & =& log\frac{1}{K_1}\hfill \end{array}$$

(11)

 □

In general, there are only four cases for maintaining the external archive, that is, nonuniform to more nonuniform, nonuniform to relatively uniform, uniform to relatively nonuniform, and uniform to more uniform. It is certain that difference entropy has a maximum and minimum, that is to say, the difference entropy is bounded.

3.3. Sketch of State Inspection

The distribution of particles in an external archive can directly react with entropy, and the status of the external archive can be indirectly speculated by the difference entropy. This reaction may be a convergence condition, diverse condition, or stagnant; then, we can change the parameters of the algorithm through the feedback information of these particles. So, the primary question is how to test the state of particles. We directly use the threshold values from the literature [21]:

$${\delta }_c=\frac{2}{H}log2$$

(12)

$${\delta }_s=\frac{2}{MK}log2,$$

(13)

where H is the number of elements in the external archive, K is the maximum capacity of the external archive, and M is the number of objective functions.

Now, the determining conditions of the environment are given as follow:

• Convergence condition: $|\Delta Entropy|>{\delta }_c$ or $\left|H\right(t)-H(t-1\left)\right|>0$
• Divergence condition: ${\delta }_s<|\Delta Entropy|<{\delta }_c$ and $H\left(t\right)=H(t-1)=K$
• Stagnant condition: $|\Delta Entropy|<{\delta }_s$ and $H\left(t\right)=H(t-1)$

Obviously, the determining conditions are sensitive to changes in H. Whenever H increases or decreases, it is seen as a convergence condition because, when H changes, at least a non-dominated solution is introduced from the external archive. The external archive changes, and the PCCS may also change. This is because an H increase or decrease may lead to the maximum or minimum of one dimension in the objective vector change.

4. Multi-Objective Particle Swarm Optimization Algorithm Based on Clone Immunity

4.1. External Archive Update

In a multi-objective particle swarm optimization algorithm, the maintenance of the external archives plays a key role. First, the calculation method of the individual density is introduced [23]. When the particle is mapped to the PCCS, the individual density of $P_i(i,1=2,\cdots ,K,K$ is the external archive size) is obtained according to the following formula:

$$Density\left(P_i\right)=\sum _{j=1,j\ne i}^K\frac{1}{{PCD(P_i,P_j)}^2},$$

(14)

where $P_j(j,1=2,\cdots ,K,j\ne i$, $PCD(P_i,P_j)$ is parallel cell distance between $P_i$ and $P_j$, which can be calculated by the following formula:

$$PCD(P_i,P_j)=\left\{\begin{array}{ll}0.5\hfill & \mathrm{if}\forall m,L_{i,m}=L_{j,m}\hfill \\ {\displaystyle \sum _{m=1}^M}|L_{i,m}-L_{j,m}|\hfill & \mathrm{if}\exists m,L_{i,m}\ne L_{j,m}\hfill \end{array}\right.$$

(15)

Then, the specific update steps in Algorithm 1 [23].

Algorithm 1: Improved external archive updating algorithm.

Input: (i) External archive A.     (ii) Maximum size of the external archive K.     (iii) New solution P obtained by the algorithm.  Output: Updated external archive $A_{new}$.  Step 1: Determine $B=A\cup P$ and calculate the objective vector value of B.  Step 2: Find the non-dominated solution set $B^{{}^{\prime }}$ for B.  Step 3: If $|B^{{}^{\prime }}|\le K$, then $A_{new}=B^{{}^{\prime }}$; if not, turn to Step 4.  Step 4: Calculate the individual density of all particles in $B^{{}^{\prime }}$.  Step 5: Place particles in ascending order by their individual density, and take the first K particles to compose $A_{new}$.

4.2. Update Strategy of the Global Best Position and Personal Best Position

Particles are selected by the lattice dominance strength [23], the specific steps in Algorithm 2.

Algorithm 2: Update strategy of the global best position.

Input: (i) External archive A and the number of objective functions M.     (ii) Status of the external archive.  Output: Global best position.  Step 1: Calculate the lattice coordinate vector $L_{k,m}$ of particles in the external archive.  Step 2: Calculate the lattice dominant strength and individual density $Density\left(i\right)$.  Step3: Determine a and b by the status of the external archive.     (i) convergence: $a=M-1,b=M+1$.     (ii) diversification: $a=M+1,b=M-1$.     (iii) stagnation: $a=M,b=M$.  Step 4: Sort the particles in descending order by the lattice dominant strength and in ascending order by the individual density.  Step 5: Take the first a particles by the individual density and the first b particles by the lattice dominant strength to store in C.  Step 6: Randomly select one particle from C as the global best position.  For the update strategy of the personal best position, we adopted the following method: [24]: $$p_{i,best}^t=\left\{\begin{array}{ll}{p}_{i,best}^{t-1}\hfill & \mathrm{if}f(p_{i,best}^t)\prec f(x_i^t)\hfill \\ x_i^t\hfill & \mathrm{if}f(x_i^t)\prec f(p_{i,best}^t)\hfill \\ randselet(p_{i,best}^{t-1},x_i^t)\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$ (16)

4.3. Parameter Selection Strategy

The PSO algorithm has a simple form and can be easily implemented, but always approaches the global and personal best position, and tends to a local best position after several iterations. In order to overcome the early-maturing phenomenon, the adaptive inertia weight based on the feedback information from status evaluation is used in this paper [25].

$$w\left(t\right)=\left\{\begin{array}{ll}w(t-1)\hfill & \mathrm{stagnation}\hfill \\ w(t-1)-2\times Ste{p}_w(1+|\Delta Entropy\left(t\right)\left|\right)\hfill & \mathrm{convergence}\hfill \\ w(t-1)+2\times Ste{p}_w\left(\right|\Delta Entropy\left(t\right)\left|\right)\hfill & \mathrm{diversificat}ion\hfill \end{array}\right.$$

(17)

where $Ste{p}_w$ is the adjusting stepsize of w, and its formula is as follows:

$$Ste{p}_w=\frac{w_{max}-w_{min}}{T_{max}},$$

(18)

where $w_{max}=0.9$, $w_{min}=0.4$, and $T_{max}$ is the maximum number of iterations. The initial value of w is 0.9.

For $c_1$ and $c_2$, nonlinear functions of the inertia weight w are used [23]:

$$\left\{\begin{array}{l}{c}_1\left(t\right)=1.167\times w^2-0.1167\times w+0.66\hfill \\ c_2\left(t\right)=3-c_1\left(t\right)\hfill \end{array}\right..$$

(19)

Thus, the learning factors $c_1$ and $c_2$ also adjust dynamically with the adjustment of the inertia weight.

According to Figure 2, the inertia weight decreased approximately linearly. With the inertia weight change, the learning factor changed.

4.4. Clone Immune Strategy

With each iteration of the algorithm, the diversity of the population decreases. In our paper, the clone, recombination, and mutation operations are used to solve this problem.

(1) Clone: For the external archive A, the population was cloned according to the crowding distance [26]. The larger the crowding distance, the larger the clone population. Then, the specific steps of clone in Algorithm 3:

Algorithm 3: Clone algorithm.

Input: (i) External archive A and number of objective functions M.     (ii) Maximum size of the clone population $NC$.  Output: Clone population $A^{{}^{\prime }}$.  Step 1: Calculate the crowding degree of the active population: $padis$.  Step 2: Calculate the number of clones of each particle with the following formula: $$q_i=⌈NC\times \frac{padi{s}_i}{{\displaystyle \sum _{j=1}^A}padi{s}_i}⌉.$$ (20)  Step 3: If the total number of clone particles exceeds $NC$, add the first $NC$ particles to $A^{{}^{\prime }}$ and return $A^{{}^{\prime }}$.

Sparse areas get more clone particles by the cloning operation, which can be performed several times.

(2) Recombination: Clone population $A^{{}^{\prime }}$ is recombined [26], the specific steps in Algorithm 4:

Algorithm 4: Recombination algorithm.

Input: (i) Clone population $A^{{}^{\prime }}$ and External archive A.     (ii) Lower bound of particles: $x_{min}$, upper bound of particles: $x_{max}$, the size of the clone population: $N^{{}^{\prime }}$, and parameter $\epsilon ,r_1$.  Output: Recombination population $A^{{}^{″}}$.  Step 1: For each particle $x^{{}^{\prime }}\left(i\right)\in A^{{}^{\prime }}$, randomly select a particle $x\left(t\right)\in A$.  Step 2: Generate a random number $rand\in [0,1]$. If $rand<r_1$, calculate $\Delta =|x^{{}^{\prime }}{\left(i\right)}_j-x{\left(t\right)}_j|$, $y_1=min(x^{{}^{\prime }}{\left(i\right)}_j,x{\left(t\right)}_j)$, $y_2=max(x^{{}^{\prime }}{\left(i\right)}_j,x{\left(t\right)}_j)$.  Step 3: If $\Delta >\epsilon$,    Step 3.1 If $y_1-x_{minj}>x_{maxj}-y_2$, then $Nbeta=1+\frac{2\times (x_{maxj}-y_2)}{y_2-y_1}$;    Step 3.2 Else $Nbeta=1+\frac{2\times (y_1-x_{minj})}{y_2-y_1}$.  Step 4: Let $\beta =\frac{1}{Nbeta}$ and $\alpha =2-{\beta }^{e^p}$.  Step 5: Generate a random number $r\in [0,1]$, if $r<\frac{1}{\alpha }$, then $\alpha =\alpha \times r$; else, $\alpha =2-\alpha \times r$.  Step 6: Let $e^p=\frac{1}{16}$ and $Nbeta={\alpha }^{e^p}$.  Step 7: Calculate $chl{d}_1=0.5\times [(y_1+y_2)-Nbeta\times (y_2-y_1)],chl{d}_2=0.5\times [(y_1+y_2)+Nbeta\times (y_2-y_1)]$.  Step 8: Generate a random number $rand\in [0,1]$, if $rand<r_2$, then $x^{{}^{\prime }}{\left(i\right)}_j=min(max(chl{d}_1,x_{maxj}),x_{minj})$; else, $x^{{}^{\prime }}{\left(i\right)}_j=min(max(chl{d}_2,x_{maxj}),x_{minj})$.

(3) Mutation: The main way to produce new particles is through a mutation that can enhance the diversity of particles [26], the specific steps in Algorithm 5.

Algorithm 5: Mutation algorithm.

Input: (i) Recombination population $A^{{}^{″}}$.     (ii) Mutation parameter $p_m$ and parameters $et{a}_m,r_1$.     (iii) Lower bound of particles: $x_{min}$, upper bound of particles: $x_{max}$, and size of clone population: $N^{{}^{\prime }}$.   Output: Mutation population $A^{{}^{‴}}$.  Step 1: Generate a random number $rand\in [0,1]$, if $rand<p_m$,    Step 1.1 If $x^{{}^{″}}{\left(i\right)}_j-x_{minj}<x_{maxj}-x^{{}^{″}}{\left(i\right)}_j$, then $\delta =\frac{x^{{}^{″}}{\left(i\right)}_j-x_{minj}}{x_{maxj}-x_{minj}}$;    Step 1.2 Else $\delta =\frac{x_{maxj}-x^{{}^{″}}{\left(i\right)}_j}{x_{maxj}-x_{minj}}$.  Step 2: Calculate $xy=1-\delta ,$ and $val=2r_1+(1-2r_1)\times (x{y}^{{\eta }_m+1})$.  Step 3: Generate a random number $rand\in [0,1]$, if $rand<r_1$, ${\delta }_q=va{l}^{\frac{1}{{\eta }_m+1}}-1$; else, ${\delta }_q=1-va{l}^{\frac{1}{{\eta }_m+1}}$.  Step 4: Calculate $x^{{}^{‴}}{\left(i\right)}_j=x^{{}^{″}}{\left(i\right)}_j+{\delta }_q\ast (x_{maxj}-x_{minj})$.

Then, The external archive will be updated by $A^{{}^{‴}}$ and Algorithm 1. We can find that the clone operation can make particles inherit good information, the recombination operation can make good genes be inherited better, and the mutation operation can increase the diversity of particles. So, all operations make the algorithm more effective.

4.5. Chaotic Strategy

In order to enhance the diversity and local search ability of the algorithm in the late iterative stages, a local chaotic searching approach is used for the external archive, the specific steps in Algorithm 6.

Algorithm 6: Local chaotic algorithm.

Input: (i) External archive A.     (ii) Maximum number of chaotic searching agents $M^{\prime }$.  Output: Clone population $A_c$.  Step 1: Randomly generate a chaotic initial point $y_0=rand$, and chaotic searching number $m=1$.  Step 2: According to the logistic map equation in chaotic systems, the following chaotic series is generated: $y_j=4y_j(1-y_j)$.  Step 3: Renew the positions based on the following equation: $x{\left(i\right)}_j=x{\left(i\right)}_j+rand{(-1)}^j{y}_j$.  Step 4: If the termination criterion is satisfied, then output the $A_c$. Otherwise, let $m=m+1$ and return to Step 2.

4.6. Chaotic Multi-Objective Particle Swarm Optimization Algorithm Incorporating Clone Immunity (CICMOPSO)

Step 1: Initialize a population of particles $X_N$, and each particle with random position vector $x_i$ and velocity vector $v_i$. Further initialize the maximum size of the external archive K, maximum size of the active population $NA$, maximum size of the clone population $NC$, lower bound of particles $x_{min}$, upper bound of particles $x_{max}$, maximum number of generations $T_{max}$, and the current generation number $T=0$.

Step 2: Calculate the fitness of all particles in $X_N\left(T\right)$.

Step 3: Initialize an update to the external archive A.

Step 4: Calculate the parallel cell coordinates of particles with Equation (5).

Step 5: Calculate the Pareto entropy and difference entropy of particles with Equations (6) and (7).

Step 6: Evaluate the states of particles by the judgement condition.

Step 7: Adjust the inertia weight by the states of particles and calculate the learning factors.

Step 8: Select the global best solution with Algorithm 2.

Step 9: Renew the positions and velocities of the particles based on Equations (2) and (3).

Step 10: For the external archive, perform the clone operation (Algorithm 3), the recombination operation (Algorithm 4), and the mutation operation (Algorithm 5), then update to the external archive A.

Step 11: Conduct chaotic search for the external archive A with Algorithm 6, obtain the chaotic population, then update to the external archive A.

Step 12: Calculate the fitness of the particles and renew every optimal position of the particles.

Step 13: Renew the external archive A with Algorithm 1.

Step 14: If the termination criterion is satisfied, then output the external archive A. Otherwise, let $T=T+1$ and return to Step 4.

4.7. Computational Complexity

In this section, the computational complexity of proposed algorithm is discussed. We consider the main steps in one generation in the main loop of algorithm. M is the objective number, N indicates the population size, K is external archive size (in our algorithm, $K=N$)

(1) The parallel cell coordinate of particles computing Process requires $O(N\times M)$ time;

(2) The clone operation requires $O\left(N\right)$ time;

(3) The recombination operation requires $O\left(NC\right)$ time, $NC$ indicates the maximum size of Clone population;

(4) The mutation operation requires $O\left(NC\right)$ time;

(5) The chaotic operation requires $O\left(N\right)$ time;

(6) The external archive update process requires $O(N^2\times M)$ time.

To summarize, the worst case overall computational complexity of CICMOPSO whin one generation is $O(N^2\times M)$, which indicates that CICMOPSO is computationally efficient.

5. Numerical Experimentation

5.1. Benchmark Function and Parameter Setting

In order to test the performance of the CICMOPSO algorithm, the ZDT [27] (

Table 2. Benchmark function—ZDT.

Name	Objective Functions	D	Variable Bounds
ZDT1	$\left\{\begin{array}{l}{f}_1\left(x\right)=x_1\hfill \\ f_2\left(x\right)=g\left(x\right)(1-\sqrt{\frac{x_1}{g\left(x\right)}})\hfill \\ g\left(x\right)=1+9\frac{{\displaystyle {\displaystyle \sum _{i=2}^n}}{x}_i}{n-1}\hfill \end{array}\right.$	30	$x_i\in [0,1]$
ZDT2	$\left\{\begin{array}{l}{f}_1\left(x\right)=x_1\hfill \\ f_2\left(x\right)=g\left(x\right)[1-{(x_1/g\left(x\right))}^2]\hfill \\ g\left(x\right)=1+9({\displaystyle {\displaystyle \sum _{i=2}^n}}{x}_i)/(n-1)\hfill \end{array}\right.$	30	$x_i\in [0,1]$
ZDT3	$\left\{\begin{array}{l}{f}_1\left(x\right)=x_1\hfill \\ f_2\left(x\right)=g\left(x\right)(1-\sqrt{x_1/g\left(x\right)}-x_{1}sin\left(10\pi x_1\right)/g\left(x\right))\hfill \\ g\left(x\right)=1+9{\displaystyle {\displaystyle \sum _{i=2}^n}}{x}_i/(n-1)\hfill \end{array}\right.$	30	$x_i\in [0,1]$
ZDT4	$\left\{\begin{array}{l}{f}_1\left(x\right)=x_1\hfill \\ f_2\left(x\right)=g\left(x\right)(1-\sqrt{x_1/g\left(x\right)})\hfill \\ g\left(x\right)=1+10(n-1)+{\displaystyle {\displaystyle \sum _{i=2}^n}}[{x_i}^2-10cos\left(4\pi x_i\right)]\hfill \end{array}\right.$	10	$x_1\in [0,1]$$x_i\in [-5,5]$$i=2,\cdots ,n$
ZDT6	$\left\{\begin{array}{l}{f}_1\left(x\right)=1-exp(-4x_1){sin}^6\left(6\pi x_1\right)\hfill \\ f_2\left(x\right)=g\left(x\right)[1-{(f_1\left(x\right)/g\left(x\right))}^2]\hfill \\ g\left(x\right)=1+9{\left[{\displaystyle {\displaystyle \sum _{i=2}^n}}{x}_i/(n-1)\right]}^{0.25}\hfill \end{array}\right.$	10	$x_i\in [0,1]$

) and DTLZ [28] (

Table 3. Benchmark function—DTLZ.

Name	Objective Functions	D	Variable Bounds
DTLZ1	$\left\{\begin{array}{l}{f}_1\left(x\right)=0.5(1+g\left(x_M\right)){\prod }_{i=1}^{M-1}{x}_i\hfill \\ f_{m=2,M-1}\left(x\right)=0.5(1+g\left(x_M\right))(1-x_{M-m+1}){\prod }_{i=1}^{M-1}{x}_i\hfill \\ f_M\left(x\right)=(1+g\left(x_M\right))sin\left(\frac{x_i\pi }{2}\right)\hfill \\ g\left(x\right)=100(|x_M|+{\displaystyle \sum _{x_i\in M}}{(x-0.5)}^2-cos\left(20\pi (x-0.5)\right))\hfill \end{array}\right.$	30	$x_i\in [0,1]$
DTLZ2	$\left\{\begin{array}{l}{f}_1\left(x\right)=0.5(1+g\left(x_M\right)){\prod }_{i=1}^M-1cos\left(\frac{x_i\pi }{2}\right)\hfill \\ f_{m=2,M-1}\left(x\right)=0.5(1+g\left(x_M\right))sin\left(\frac{x_{M-m+1}\pi }{2}\right){\prod }_{i=1}^M-1cos\left(\frac{x_i\pi }{2}\right)\hfill \\ f_M\left(x\right)=(1+g\left(x_M\right))sin\left(\frac{x_i\pi }{2}\right)\hfill \\ g\left(x\right)={\displaystyle \sum _{x_i\in M}}{(x-0.5)}^2\hfill \end{array}\right.$	30	$x_i\in [0,1]$
DTLZ3	$\left\{\begin{array}{l}{f}_1\left(x\right)=0.5(1+g\left(x_M\right)){\prod }_{i=1}^M-1cos\left(\frac{x_i\pi }{2}\right)\hfill \\ f_{m=2,M-1}\left(x\right)=0.5(1+g\left(x_M\right))sin\left(\frac{x_{M-m+1}\pi }{2}\right){\prod }_{i=1}^M-1cos\left(\frac{x_i\pi }{2}\right)\hfill \\ f_M\left(x\right)=100(|x_M|+{\displaystyle \sum _{x_i\in M}}{(x-0.5)}^2-cos\left(20\pi (x_i-0.5)\right))\hfill \\ g\left(x\right)={\displaystyle \sum _{x_i\in M}}{(x-0.5)}^2\hfill \end{array}\right.$	30	$x_i\in [0,1]$
DTLZ4	$\left\{\begin{array}{l}{f}_1\left(x\right)=0.5(1+g\left(x_M\right)){\prod }_{i=1}^M-1cos\left(\frac{x_i\pi }{2}\right)\hfill \\ f_{m=2,M-1}\left(x\right)=0.5(1+g\left(x_M\right))sin\left(\frac{x_{M-m+1}\pi }{2}\right){\prod }_{i=1}^M-1cos\left(\frac{x_i\pi }{2}\right)\hfill \\ f_M\left(x\right)=(1+g\left(x_M\right))sin\left(x_1^a\pi \right)\hfill \\ g\left(x\right)={\displaystyle \sum _{x_i\in M}}{(x-0.5)}^2,a=100\hfill \end{array}\right.$	30	$x_i\in [0,1]$
DTLZ5	$\left\{\begin{array}{l}{f}_1\left(x\right)=0.5(1+g\left(x_M\right)){\prod }_{i=1}^M-1cos(\frac{\pi (1+2g\left(x_M\right)x_1)}{2(2+(1+g\left(x_M\right)))})\hfill \\ f_{m=2,M-1}\left(x\right)=0.5(1+g\left(x_M\right))sin\left(\frac{x_{M-m+1}\pi }{2}\right){\prod }_{i=1}^M-mcos(\frac{\pi (1+2g\left(x_M\right)x_1)}{2(2+(1+g\left(x_M\right)))})\hfill \\ f_M\left(x\right)=(1+g\left(x_M\right))sin\left(\frac{x_1\pi }{2}\right)\hfill \\ g\left(x\right)={\displaystyle \sum _{x_i\in M}}{(x-0.5)}^2\hfill \end{array}\right.$	30	$x_i\in [0,1]$
DTLZ6	$\left\{\begin{array}{l}{f}_1\left(x\right)=0.5(1+g\left(x_M\right)){\prod }_{i=1}^M-1cos(\frac{\pi (1+2g\left(x_M\right)x_1)}{2(2+(1+g\left(x_M\right)))})\hfill \\ f_{m=2,M-1}\left(x\right)=0.5(1+g\left(x_M\right))sin(\frac{x_{M-m+1}\pi }{2}){\prod }_{i=1}^M-mcos(\frac{\pi (1+2g\left(x_M\right)x_1)}{2(2+(1+g\left(x_M\right)))})\hfill \\ f_M\left(x\right)=(1+g\left(x_M\right))sin\left(\frac{x_1\pi }{2}\right)\hfill \\ g\left(x\right)={\displaystyle \sum _{x_i\in M}}{x}_i^0.1\hfill \end{array}\right.$	30	$x_i\in [0,1]$

) test functions were used in the experiments in this paper. We then compared the algorithm to NSGA-II [21], SPAE2 [29], and NICPSO [26]. For the parameters in this paper, the initial values of the learning factors were $c_1=c_2=1.5$, the initial value of w was $0.9$, the recombination possibility $P_c=0.8$, the mutation probability $P_m=1/n$, which was the reciprocal of the decision variable dimension, the maximum number of the external archive was 100, the population size was 100, and maximum number of generations $T_{max}=300$. Experiments were independently repeated 30 times. The optimization problems are described in Table 2.

5.2. Performance Indicators

(1) Convergence index [26]: Generational distance (GD) indicates the average distance between the obtained Pareto front and the true Pareto front. GD is defined as follow:

$$GD=\frac{{\displaystyle {\displaystyle \sum _{i=1}^k}}{d}_i}{n}$$

(21)

where $d_i=\sqrt{{\displaystyle \sum _{j=1}^n}{(P{F}_j-P{F}_j^t)}^2}$, kis the number of Pareto Solutions, n is the number of objective functions, $P{F}_j$ indicate the $j^{th}$ obtained Pareto solution, $P{F}_j^t$ indicate the nearest point on true Pareto front from $P{F}_j$.

(2) Spacing index [26]: Spacing (S) specify the spread of the obtained Pareto front. It is defined as follow:

$$S=\sqrt{\frac{1}{k-1}{\displaystyle {\displaystyle \sum _{p=2}^k}}{(D_p-\overline{D})}^2}$$

(22)

where $D_p={\displaystyle {\displaystyle \sum _{i=1}^n}}|P{F}_{i,(p-1)}-P{F}_{i,p}|$, is the absolute difference between the two consecutive solutions in the obtained Pareto front $PF$, $\overline{D}$ is the average of all $D_p$.

5.3. Numerical Results

The convergence results are shown in

Table 4. Comparison of results for the convergence index (GD).

		NSGA-II	SPEA2	NICPSO	CICMOPSO
ZDT1	Aver	$1.14\times {10}^{-3}$	$3.82\times {10}^{-3}$	$1.08\times {10}^{-3}$	$2.95\times {10}^{-4}$
	Var	$1.41\times {10}^{-4}$	$4.91\times {10}^{-3}$	$4.26\times {10}^{-5}$	$4.23\times {10}^{-5}$
ZDT2	Aver	$4.88\times {10}^{-1}$	$8.61\times {10}^{-3}$	$7.74\times {10}^{-4}$	$8.29\times {10}^{-5}$
	Var	$2.77\times {10}^{-2}$	$2.60\times {10}^{-3}$	$2.99\times {10}^{-5}$	$2.86\times {10}^{-5}$
ZDT3	Aver	$2.48\times {10}^{-3}$	$9.72\times {10}^{-3}$	$8.93\times {10}^{-4}$	$6.75\times {10}^{-5}$
	Var	$1.27\times {10}^{-4}$	$5.23\times {10}^{-3}$	$4.02\times {10}^{-5}$	$4.03\times {10}^{-5}$
ZDT4	Aver	$5.13\times {10}^{-1}$	$9.98$	$1.37\times {10}^{-3}$	$8.55\times {10}^{-2}$
	Var	$1.18\times {10}^{-1}$	$2.01\times {10}^1$	$7.42\times {10}^{-3}$	$5.72\times {10}^{-2}$
ZDT6	Aver	$7.58\times {10}^{-2}$	$1.93\times {10}^{-2}$	$3.08\times {10}^{-3}$	$1.19\times {10}^{-4}$
	Var	$6.08\times {10}^{-3}$	$1.40\times {10}^{-3}$	$1.33\times {10}^{-4}$	$8.33\times {10}^{-6}$

. The CICMOPSO algorithm had better results than the NSGA-II and SPEA2 algorithms for all benchmark functions. For ZDT4, CICMOPSO was worse than NICPSO, but for other functions, CICMOPSO had the better performance. This finding indicates that the resulting Pareto fronts obtained by CICMOPSO are closer to the true Pareto fronts, and CICMOPSO can effectively improve convergence.

From

Table 5. Comparison of results for the Spacing index(S).

		NSGA-II	SPEA2	NICPSO	CICMOPSO
ZDT1	Aver	$5.04\times {10}^{-1}$	$2.96\times {10}^{-1}$	$8.78\times {10}^{-2}$	$3.33\times {10}^{-2}$
	Var	$3.93\times {10}^{-2}$	$1.09\times {10}^{-1}$	$4.83\times {10}^{-2}$	$1.76\times {10}^{-2}$
ZDT2	Aver	$4.73\times {10}^{-1}$	$5.05\times {10}^{-1}$	$1.22\times {10}^{-2}$	$1.37\times {10}^{-2}$
	Var	$2.99\times {10}^{-2}$	$1.84\times {10}^{-1}$	$1.87\times {10}^{-2}$	$2.09\times {10}^{-2}$
ZDT3	Aver	$5.90\times {10}^{-1}$	$5.03\times {10}^{-1}$	$1.88\times {10}^{-2}$	$1.74\times {10}^{-2}$
	Var	$3.04\times {10}^{-2}$	$9.73\times {10}^{-2}$	$1.66\times {10}^{-2}$	$3.39\times {10}^{-2}$
ZDT4	Aver	$7.03\times {10}^{-1}$	$8.70\times {10}^{-1}$	$1.17\times {10}^{-2}$	$2.23\times {10}^{-2}$
	Var	$6.47\times {10}^{-2}$	$1.01\times {10}^{-1}$	$7.37\times {10}^{-3}$	$5.36\times {10}^{-2}$
ZDT6	Aver	$4.86\times {10}^{-1}$	$2.49\times {10}^{-1}$	$3.27\times {10}^{-3}$	$2.87\times {10}^{-2}$
	Var	$3.61\times {10}^{-2}$	$4.97\times {10}^{-2}$	$1.96\times {10}^{-4}$	$3.59\times {10}^{-2}$

, we can see that NICPSO and CICMOPSO produced better results for the Spacing index (S) than NSGA-II and SPEA2. For ZDT1 and ZDT2, CICMOPSO has better performance than NICPSO, which has minimum mean and variance, but for other functions, NICPSO performed even better.

It is well known that multi-objective optimization problems, especially when the number of objectives ≥3, are more difficult to solve. Simulation results in

Table 6. Convergence index and Spacing index for the DTLZ model.

	Convergence		Uniformity
	Aver	Var	Aver	Var
DTlZ1	$3.54\times {10}^{-1}$	$2.88\times {10}^{-1}$	$3.22\times {10}^{-1}$	$2.76\times {10}^{-1}$
DTlZ2	$1.93\times {10}^{-1}$	$1.62\times {10}^{-1}$	$7.40\times {10}^{-1}$	$5.43\times {10}^{-1}$
DTlZ3	$5.42$	$8.84$	$1.74$	$5.00$
DTlZ4	$1.95$	$6.15$	$9.17$	$1.21$
DTlZ5	$1.37\times {10}^{-1}$	$8.29\times {10}^{-2}$	$4.51\times {10}^{-1}$	$9.62\times {10}^{-2}$
DTlZ6	$1.45$	$1.07$	$9.29$	$6.69$

show that CICMOPSO can solve most objective function problems.

6. Conclusions

A chaotic multi-objective particle swarm optimization algorithm incorporating clone immunity (CICMOPSO), which is used to solve multi-objective problems, was proposed in this paper. CICMOPSO uses the clone immunity strategy to maintain an external file and avoid the algorithm falling into a local optimal solution. Additionally, the Pareto entropy is used to dynamically adjust the algorithm parameters.The experimental results showed that CICMOPSO significantly outperforms all other algorithms based on the test problems with respect to two metrics. For problems with two-objective functions, CICMOPSO performed well, but for more-objective functions, the results are not satisfactory. we will improve CICMOPSO to make it suitable for solving more problems in the near future.