An integrated cuckoo search optimizer for single and multi-objective optimization problems

Two subpopulation structure

Many species boast complex societies with multiple levels of communities. A common case is when two dominant levels exist, one corresponding to leaders and the other consisting of followers (Ferdinandy, Ozogány & Vicsek, 2017). In this paper we change the structure of population of the CS algorithm. At each iteration, the population is divided into two subpopulations according to the individual fitness. Those who own better fitness form high fitness subpopulation while those who own worse fitness form low fitness subpopulation, seen from Fig. 2. The pentagram symbols represents the individuals with high fitness and the circle symbols represents the individuals with low fitness, as shown in Fig. 2.

The high fitness subpopulation is responsible for exploiting better solutions, while the low fitness subpopulation is responsible for exploring unknown solutions. Exploration plays a major role in the early iterations while exploitation plays a major role in the late iterations. So the size of each subpopulation is dynamically changed as the iteration number increases. The size of the high fitness subpopulation increases with the iteration while the size of the low fitness subpopulation decreases with the iteration. In the whole iteration process, the size of the whole population remains unchanged. The value of $P^L$ indicates the proportion of the low fitness subpopulation in the whole population. The value of $P^L$ decreases linearly from $P_{max}^L$ to $P_{min}^L$ with the increasing of iteration times, as shown in Eq. (4).

(4) $$P^L=P_{max}^L-{\displaystyle \frac{k}{E}\left(P_{max}^L-P_{min}^L\right)}$$

(5) $$P^H=1-P^L$$where $k$ is the current iteration and $E$ is the maximum iterations. $P^L$ is the proportion of the low fitness subpopulation in the whole population. $P^H$ is the proportion of the high fitness subpopulation in the whole population.

The two subpopulations employ different search strategies which can guarantee the diversity. Meanwhile, two subpopulations adopt different search strategies to realize a balanced combination of a local random walk and a global explorative random walk. A lévy flight is employed to exploit the most promising area, so the high fitness population adopts a lévy flight mode. The previous literature has also shown that some of the new solutions should be produced by lévy flight around the best solution obtained so far, which will often speed up the local search (Pavlyukevich, 2007). A random walk method is used to explore the unknown area, so the low fitness population adopts a random walk method. In this paper, a random method is adopted according to Eq. (3).

Adaptive step size strategy

The step scaling factor ( $stepSize$ in Eq. (1)) plays an important role in updating the position. This parameter is related to the scales of the problem. When the variable range increases, the step size should be increased accordingly. When the variable range decreases, the step size should be reduced accordingly.

In addition, it is very important to balance the ability of global exploration and local exploitation in designing a population-based algorithm. Exploration plays a major role at the beginning of the algorithm, and the step size should be larger. In the later stage of the algorithm, exploitation plays a major role and a big step size may make it difficult for the algorithm to converge to the optimal point. So the step size should be smaller in the later stage of the algorithm.

Based on the above two points, in order to make step size adapt to different optimization problems, a self-adaptive step size mechanism for high fitness subpopulation is introduced. The value of the step scaling factor $stepSize$ decreases linearly from $ste{p}_{max}$ to $ste{p}_{min}$ with the increasing of iteration times, as shown in the Eq. (6).

(6) $$stepSize=ste{p}_{max}-{\displaystyle \frac{k}{E}\left(ste{p}_{max}-ste{p}_{min}\right)}$$

(7) $$ste{p}_{max}=0.01\ast range$$

(8) $$ste{p}_{min}=0.001\ast range$$where $k$ is the current iteration and $E$ is the maximum iterations. $range$ is the variable range. For problems with different variable ranges in different dimensions, $range$ is set to the maximum value.

Crossover based on DE operation

Because the two subpopulations adopt different search strategies, the solutions have diversity on the whole, seen from Fig. 2. To take advantage of the diversity of these solutions to find the potential optimal solution, the two subpopulations are combined into a population and a crossover operation is applied to the population.

Differential evolution (DE) proposed by Storn & Price (1997) performs very well in convergence (Gong, Cai & Ling, 2011). Especially, DE has a good performance on searching the local optima and good robustness (Mallipeddi et al., 2011). Based on the fast convergence speed of DE, a crossover phase based on DE is added after the two population are combined into a whole population. From the perspective of biodiversity, this phase can increase the individual information exchange. The steps of the crossover operation are given in Table 2.

(9) $$x_{ij}=x_{rj}+F\left(x_{pj}-x_{qj}\right)$$where $i$, r, $p,q$ are random numbers between 1 and the size of population and mutually different, $j$ is a dimension index between [1, D]. F is the differential weight in the range of [0, 2].

In the crossover operation, the differential mutation operator is the main operation. Figure 3 shows a two-dimensional example for producing new position using Eq. (9). The DE operator is employed to produce a self-adaptive mechanism to change the search range. In the initial stage, individual differences are large, and the algorithm searches in a wide range. At the later stage of algorithm, the population is in a state of convergence, and the individual difference is small. The algorithm searches in a narrow range in the later period of the optimization. Therefore, using DE operator alone will cause the diversity to disappear quickly.

In the proposed ICSO, the initial population is divided into two subpopulations according to the fitness. The two subpopulations update their positions according to different search methods to ensure the diversity of the solutions. Crossover operation is applied after the two subpopulations are merged. After the crossover operation, the whole population is divided into two subpopulations according to their fitness. The above process goes on and on until the end of the algorithm iteration.

The proposed algorithm for single-objective problems

The procedure of ICSO is shown in Table 3 and the flowchart of ICSO is summarized in Fig. 4.

From Table 3, we can see that initialization phase is responsible for initializing the solution. Eq. (10) gives the way to initialize the solution.

(10) $$x_{i,j}=x_j^{min}+rand[0,1](x_j^{max}-x_j^{min})$$where $x_{i,j}$ ( $1\le i\le N,1\le j\le D$) represents the initial value of the jth variable of ith agent. $N$ and D are the size of population and the number of variables, respectively. $x_j^{max}$ and $x_j^{min}$ denote the upper and lower bounds of the jth parameter, respectively.

The proposed algorithm for multi-objective problems

Many problems in real life are multi-objective. Without loss of generality, the mathematical expression of the multi-objective problem is as follows:

(11) $$\mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}f\left(x\right)=\left\{f_i\left(x\right)\right\},i=1,2,\dots ,m$$

(12) $$\mathrm{S}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}:h_j\left(x\right)=0,j=1,2,\dots ,p$$

(13) $$g_k\left(x\right)\le 0,k=1,2,\dots ,q$$

(14) $$x=\left(x_1,x_2,\dots ,x_n\right)$$

(15) $$lo{w}_{\alpha }\le x_{\alpha }\le hig{h}_{\alpha },\alpha =1,2,\dots ,n$$where $m>1$, $f\left(x\right)$ is called multi-objective functions. p is the number of equality constraints. $q$ is the number of inequality constraints. $lo{w}_{\alpha }$ is the lower bound of the $\alpha \mathrm{t}\mathrm{h}$ variable. $hig{h}_{\alpha }$ is the upper bound of the $\alpha \mathrm{t}\mathrm{h}$ variable.

The Pareto front method is used to solve the multi-objective problems through the calculation of Pareto front. Let $U=\left(u_1,u_2,\dots ,u_n\right)$, $V=\left(v_1,v_2,\dots ,v_n\right)$ are two vectors. $U$ is said to dominate $V$ if and only if $U$ is partially less than $V$ in objective space. $x^{\ast }$ is said to be a Pareto optimal solution if and only if any other solutions can’t be detected to dominate $x^{\ast }$. Pareto optimal solutions are also called non-dominated solutions. A set of Pareto optimal solution is called Pareto optimal front.

In order to evaluate the performance of multi-objective optimization algorithm quantitatively, two kinds of performance metrics are given here. One is convergence metric, the other is diversity metric convergence (Deb et al., 2002). The generation distance (GD) measure is defined as the criterion of convergence, which is calculated according to Eq. (16).

(16) $$\mathrm{G}\mathrm{D}=\sqrt{{\displaystyle \frac{1}{n_{pf}}{\sum }_{i=1}^{n_{pf}}{d}_i^2}}$$where $n_{pf}$ is the number of optimal solutions in the obtained Pareto front. $d_i$ is the Euclidean distance from each non dominated solution to the nearest true Pareto front solution.

The computing of GD is illustrated by Fig. 5 where triangles and circles correspond to the solutions on the true Pareto front and the obtained Pareto front, respectively.

The metric of spacing (S) measure is defined as the criterion of diversity, which is calculated according to Eq. (17).

(17) $$S=\sqrt{{\displaystyle \frac{1}{n_{pf}-1}{\sum }_{i=1}^{n_{pf}}{\left(d_i-\overline{d}\right)}^2}}$$where $n_{pf}$ is the number of optimal solutions in the obtained Pareto front. $d_i$ is the crowding distance of the $i$-th solution. $\overline{d}$ is the average of all $d_i$.

The computing of crowding-distance is illustrated by Fig. 6. The crowding distance of the $i$-th solution is the average side length of the rectangle formed by the dotted line (Deb et al., 2002).

For the advantages of the ICSO, this article proposed a multi-objective ICSO (MOICSO) which can deal with multi-objective optimization problems. The main steps of the MOICSO are shown in Table 4.

It is worth noting that compared with ICSO, MOICSO has four more phases. Non-dominated sorting phase is to calculate the non-dominated solutions in the whole population. Crowding distance calculating procedure has been illustrated by Fig. 6.

Benchmark problems

Numerical optimization problem solving capability of ICSO is examined by using 24 classic benchmark problems (Civicioglu & Besdokb, 2019; Karaboga & Akay, 2009), f₁–f₂₄. All benchmark functions used their standard ranges. Detailed mathematical definitions of f₁–f₂₄ are given in Table 5.

Functions $f_{20}$– $f_{24}$ are four rotated functions (Liang et al., 2006). The fitness of the rotated function is calculated with a rotated variable $y$ instead of $x$. The rotated variable $y$ is calculated by the original variable $x$ left multiplied an orthogonal matrix. The orthogonal matrix is calculated using Salomon’s method (Salomon, 1996).

Parameter study

The influences of the parameters including $(f_7)$ and $(f_8)$ are experimented in this section. The dimension of the test functions including Sphere, Powers, Ackley and Griewank is set to 30. The population size of each algorithm is set to 50. The computational cost of each experiment was set to 100,000 function evaluations. Each test problem is solved by assigning different values to the two parameters in 30 independent runs. The numerical results in terms of the mean results (Mean), standard deviation (Std), the best result (Best) and the worst result (Worst) of the optimal function value were given in Tables 6–9. The results demonstrate that ICSO performs best on all test functions when both the parameter $CR$ and $F$ are set to 0.1. Therefore, these two parameters are set to 0.1 in the next experiments.

Comparison with other algorithms

In order to compare the performance of ICSO, DDICS (Wang, Yin & Zhong, 2013), DE/rand (Storn & Price, 1997), CSDE (Zhang, Ding & Jia, 2019) and CSEI (Mareli & Twala, 2018) were employed for comparison. DE/rand is a classical population-based algorithm. Here DE/rand stands for DE/rand/1/bin (Storn & Price, 1997). Detailed descriptions of DDICS, CSDE and CSEI are given in “CS Algorithm and its Variants”.

The population size of ICSO and DDICS was 50. Dimensions of $f_1$– $f_9$ and $f_{11}$– $f_{24}$ are selected as 50. Dimensions of $f_{10}$ is selected as 24. The maximum evaluation count for functions $f_1$– $f_{24}$ is 200,000. $P_{max}^L\mathrm{a}\mathrm{n}\mathrm{d}{P}_{min}^L$ are set to 0.9 and 0.2 respectively. To gather significant statistical results, each algorithm carried out 30 independent runs. For the other specific parameters for involved algorithms, we can follow parameters settings of the original literatures of DDICS (Wang, Yin & Zhong, 2013), CSEI (Mareli & Twala, 2018), DE/rand and CSDE (Zhang, Ding & Jia, 2019).

The results of our experiments are illustrated in Table 10 and the average convergence curves obtained for ICSO, DDICS, CSDE, CSEI and DE/rand are depicted in Fig. 7. M stands for the mean value and S stands for the standard deviation. B stands for the best value and W stands for the worst value. R stands for the performance order of the different algorithms on each standard function. The best mean values obtained on each function are marked as bold. It is obvious that ICSO performed best on most functions.

Table 10 demonstrates that ICSO achieves approving results on most of the twenty-four functions. ICSO has achieved the best mean values except $f_2,f_3,f_7\mathrm{a}\mathrm{n}\mathrm{d}{f}_{13}$. Compared with the other algorithms, ICSO performs very well when solving the functions $f_1$, $f_4$, $f_5$, $f_6$, $f_8$, $f_{12}$, $f_{14}$, $f_{16}$, $f_{18}$, $f_{19}$, $f_{23}$ and $f_{24}$.

Particularly, ICSO has a strong ability to mine the optimal value in the later stage of iteration on $f_4$ and $f_{12}$, seen from Fig. 7. The performance of ICSO on $f_4$ is much similar with it on $f_{12}$. ICSO converged very fast at the beginning and then trapped in local minimum. However, after a certain number of iterations, ICSO jumps out of the local optimum. With multiple search strategies, ICSO keeps sufficient diversity to escape local optimum on functions like $f_4$ and $f_{12}$, even at the later stage of algorithm iteration.

As per the results and discussions of this study, we can conclude that ICSO can enhance both exploration and exploitation capabilities effectively. As can be clearly seen from Fig. 7, ICSO seemed to be able to sustain improvement especially on $f_1$, $f_5$, $f_6$, $f_8$, $f_{14}$, $f_{20}$, $f_{21}$ and $f_{24}$.

Statistical analysis

In this section, two statistical evaluation approaches including Iman-Daveport and Holm tests are considered for validating the performance of the proposed ICSO. A detailed introduction of the two statistical evaluation methods can be found in reference (García et al., 2010).

The results of the Iman–Davenport test and Holm tests are given in Tables 11 and 12 respectively. The critical value listed in the Table 11 and the $\alpha /i$ values listed in the Table 12 are with level of 0.05.

From Table 11, there are significant differences in the ranking of algorithms because the Iman-Davenport values are greater than the critical value. From Table 12, there are significant differences between comparison algorithms and control algorithm (ICSO) because the $p$ values of CSDE, DE/rand, DDICS and CSEI are less than their $\alpha /i$ values.

In this article, the validity and stability of ICSO are also studied by analysis of variance (ANOVA). Figure 8 shows the box plots of the statistical performance of all algorithms. The general characteristics of the distribution can be obtained by looking at these box plots. It can be seen clearly from Fig. 8 that ICSO obtains a good compromise solution variance distribution on most test functions.

Dynamic size of subpopulation

This simulation is conducted to investigate how the size of subpopulation in the ICSO model changes along with the generations on the benchmark environments. In order to clearly see the change of subpopulation size, the initial population size is set to 100. The subpopulation evolution based on ICSO model was simulated on two-dimensional Ackley function. Figure 9 shows the positions of the individuals and the subpopulation variation in different phases, where each red star represents a position with high fitness and each blue circle represents a position with low fitness.

As can be seen from Fig. 9, the initial positions are distributed randomly over the map defined by the two-dimensional Ackley function at the first stage (FEs = 100). We can directly see from Fig. 9 that with the iteration, the size of high fitness subpopulation increases and the size of low fitness subpopulation decreases. From the first phase (FEs = 100) to third phase (FEs = 1,500), the whole population moves toward the optimal position. The two subpopulations adopt different ways to update the position, and exploit the optimal region and explore the unknown region respectively. Under the collaborative effect of the two subpopulations, ICSO can make the whole population move toward the optimal position.

Benchmark problems

The test problems includes seven problems which are denoted by ZDT1, ZDT3, ZDT6, DTLZ1, DTLZ2, DTLZ3 and DTLZ6. ZDT1, ZDT3 and ZDT6 have two objective functions. DTLZ1, DTLZ2, DTLZ3 and DTLZ6 have three objective functions.

The ZDT1 problem (Deb et al., 2002) is described as follows:

(18) $$\mathrm{Z}\mathrm{D}\mathrm{T}1:\{\begin{array}{l}\mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_1\left(x\right)=x_1\\ \mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_2\left(x\right)=g\left(x\right)\left[1-\sqrt{{\displaystyle \frac{x_1}{g\left(x\right)}}}\right]\\ g\left(x\right)=1+{\displaystyle \frac{9\left({\sum }_{i=2}^n{x}_i\right)}{n-1}}\end{array}$$where $0\le x_i\le 1,i=1,2,\dots ,30.$ The Pareto optimal region corresponds to $x_1\in \left[0,1\right]$ and $x_i=0,i=2,3,\dots ,n$.

The ZDT3 problem (Deb et al., 2002) is described as follows:

(19) $$\mathrm{Z}\mathrm{D}\mathrm{T}3:\{\begin{array}{l}\mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_1\left(x\right)=x_1\\ \mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_2\left(x\right)=g\left(x\right)\left[1-\sqrt{x_1/g\left(x\right)}-{\displaystyle \frac{x_1}{g\left(x\right)}\sin \left(10\pi x_1\right)}\right]\\ g\left(x\right)=1+{\displaystyle \frac{9\left({\sum }_{i=2}^n{x}_i\right)}{n-1}}\end{array}$$where $0\le x_i\le 1,i=1,2,\dots ,30.$ The Pareto optimal region corresponds to $x_1\in \left[0,1\right]$ and $x_i=0,i=2,3,\dots ,n$.

The ZDT6 problem (Deb et al., 2002) is described as follows:

(20) $$\mathrm{Z}\mathrm{D}\mathrm{T}6:\{\begin{array}{l}\mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_1\left(x\right)=1-\exp \left(-4x_1\right)si{n}^6\left(6\pi x_1\right)\\ \mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_2\left(x\right)=g\left(x\right)\left[1-{\left({\displaystyle \frac{f_1\left(x\right)}{g\left(x\right)}}\right)}^2\right]\\ g\left(x\right)=1+9{\left[\left(\sum _{i=2}^n{x}_i\right)/\left(n-1\right)\right]}^{0.25}\end{array}$$where $0\le x_i\le 1,i=1,2,\dots ,10.$ The Pareto optimal region corresponds to $x_1\in \left[0,1\right]$ and $x_i=0,i=2,3,\dots ,n$.

The DTLZ1 problem (Deb et al., 2002) is described as follows:

(21) $$\mathrm{D}\mathrm{T}\mathrm{L}\mathrm{Z}1:\{\begin{array}{l}\mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_1\left(x\right)={\displaystyle \frac{1}{2}{x}_1{x}_2\dots x_{M-1}\left(1+g\left(x_M\right)\right)}\\ \mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_2\left(x\right)={\displaystyle \frac{1}{2}{x}_1{x}_2\dots \left(1-x_{M-1}\right)\left(1+g\left(x_M\right)\right)}\\ \begin{array}{ll}⋮& ⋮\end{array}\\ \mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_{M-1}\left(x\right)={\displaystyle \frac{1}{2}{x}_1\left(1-x_2\right)\left(1+g\left(x_M\right)\right)}\\ \mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_M\left(x\right)={\displaystyle \frac{1}{2}\left(1-x_1\right)\left(1+g\left(x_M\right)\right)}\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}0\le x_i\le 1,i=1,2,\dots ,n\\ \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}g\left(x_M\right)=100\left(\left|x_M\right|+\sum _{x_i\in x_M}{\left(x_i-0.5\right)}^2-\cos \left(20\pi \left(x_i-0.5\right)\right)\right)\end{array}$$here, $M$ = 3, $\left|x_M\right|=k$ = 5, the total number of variables is $n=M+k-1$. The Pareto-optimal solution corresponds to $x_i^{\ast }=0.5\left(x_i^{\ast }\in x_M\right)$ and the objective function values lie on the linear hyper plane: ${\sum }_{m=1}^M{f}_m^{\ast }=0.5$.

The DTLZ2 problem (Deb et al., 2002) is described as follows:

(22) $$\mathrm{D}\mathrm{T}\mathrm{L}\mathrm{Z}2:\{\begin{array}{l}\mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_1\left(x\right)=\left(1+g\left(x_M\right)\right)\cos \left({\displaystyle \frac{x_1\pi }{2}}\right)\dots \cos \left({\displaystyle \frac{x_{M-1}\pi }{2}}\right)\\ \mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_2\left(x\right)=\left(1+g\left(x_M\right)\right)\cos \left({\displaystyle \frac{x_1\pi }{2}}\right)\dots \sin \left({\displaystyle \frac{x_{M-1}\pi }{2}}\right)\\ \hfill \begin{array}{cc}⋮& ⋮\end{array}\hfill \\ \mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_M\left(x\right)=\left(1+g\left(x_M\right)\right)\sin \left({\displaystyle \frac{x_1\pi }{2}}\right)\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}0\le x_i\le 1,i=1,2,\dots ,n\\ \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}g\left(x_M\right)=\sum _{x_i\in x_M}{\left(x_i-0.5\right)}^2\end{array}$$here, $M$ = 3, $\left|x_M\right|=k$ = 10, the total number of variables is $n=M+k-1$. The Pareto-optimal solution corresponds to $x_i^{\ast }=0.5\left(x_i^{\ast }\in x_M\right)$ and the objective function values must satisfy the ${\sum }_{m=1}^M(f_m^{\ast }{)}^2=1$.

The DTLZ3 problem (Deb et al., 2002) is described as follows:

(23) $$\mathrm{D}\mathrm{T}\mathrm{L}\mathrm{Z}3:\{\begin{array}{l}\mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_1\left(x\right)=\left(1+g\left(x_M\right)\right)\cos \left({\displaystyle \frac{x_1\pi }{2}}\right)\dots \cos \left({\displaystyle \frac{x_{M-1}\pi }{2}}\right)\\ \mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_2\left(x\right)=\left(1+g\left(x_M\right)\right)\cos \left({\displaystyle \frac{x_1\pi }{2}}\right)\dots \sin \left({\displaystyle \frac{x_{M-1}\pi }{2}}\right)\\ \hfill \begin{array}{cc}⋮& ⋮\end{array}\hfill \\ \mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_M\left(x\right)=\left(1+g\left(x_M\right)\right)\sin \left({\displaystyle \frac{x_1\pi }{2}}\right)\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}0\le x_i\le 1,i=1,2,\dots ,n\\ \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}g\left(x_M\right)=100\left(\left|x_M\right|+\sum _{x_i\in x_M}{\left(x_i-0.5\right)}^2\right)-\cos \left(20\pi \left(x_i-0.5\right)\right)\end{array}$$here, $M$ = 3, $\left|x_M\right|=k$ = 10, the total number of variables is n = M + k − 1. The Pareto-optimal solution corresponds to $x_i^{\ast }=0.5\left(x_i^{\ast }\in x_M\right)$.

The DTLZ6 problem (Deb et al., 2002) is described as follows:

(24) $$\mathrm{D}\mathrm{T}\mathrm{L}\mathrm{Z}6:\{\begin{array}{l}\mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_1\left(x\right)=x_1\\ \mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_2\left(x\right)=x_2\\ \begin{array}{l}\hfill \begin{array}{ll}⋮& ⋮\end{array}\hfill \\ \mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_{M-1}\left(x\right)=x_{M-1}\end{array}\\ \begin{array}{l}\mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{f}_M\left(x\right)=\left(1+g\left(x_M\right)\right)h\left(f_1,f_2,\dots ,f_{M-1},g\right)\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}0\le x_i\le 1,i=1,2,\dots ,n\\ \begin{array}{l}\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}g\left(x_M\right)=1+{\displaystyle \frac{g}{\left|x_M\right|}\sum _{x_i\in x_M}{x}_i}\\ h=M-\sum _{i=1}^{M-1}\left[{\displaystyle \frac{f_i}{1+g}\left(1+\sin \left(3\pi f_i\right)\right)}\right]\end{array}\end{array}\end{array}$$here, M = 3, |x_M| = k = 20, the total number of variables is n = M + k − 1.

Experiment sets

To evaluate the performance of MOICSO for multi-objective problems, we compared it with MODA (Mirjalili, 2016) and MOWCA (Sadollah et al., 2015). In order to do meaningful statistical analysis, each algorithm runs for 10 times. The maximum evaluation count is 20,000. The population size of MOICSO was 100. The experimental results include the best value, the worst value, the mean value, the median value and the variance of the convergence measure and the diversity measure. Parameters for MOICSO are the same as ICSO in the section “Experiment and results of ICSO”. For the other specific parameters for involved algorithms, we can follow parameters settings of the original literatures of MODA (Mirjalili, 2016) and MOWCA (Sadollah et al., 2015).

Experiment results and analysis

Figures 10–12 show the optimal front obtained by three algorithms on two objective test problems. The continuous solid line in the figure represents the Pareto optimal front of the theory, and the circle represents the non-dominated solution set obtained by each algorithm. Figures 13–20 show theoretical Pareto front and the optimal solution set obtained by three algorithms on three objective test problems.

For ZDT1 problem, it can be seen from Table 13 that MOICSO has better performance than MODA and MOWCA in diversity. However, MOWCA got the best performance in convergence. It can be seen from Fig. 10 that MODA cannot converge to the Pareto front of ZDT1.

For ZDT3 problem, it can be seen from Table 14 that MODA has better performance than MOICSO and MOWCA in diversity. However, it can be seen from Fig. 11 that MODA cannot converge to the Pareto front of ZDT3. MOWCA got the best performance in convergence. MOICSO outperforms MODA in terms of convergence and MOWCA in terms of diversity.

For ZDT6 problem, MOICSO outperforms the other two comparison algorithms in terms of convergence and diversity, as seen from Fig. 12 and Table 15. In particular, the performance of MOICSO in convergence metric is three order of magnitude better than that of MODA and MOWCA, as shown in Table 15. The performance of MOICSO in diversity metric is two order of magnitude better than that of MOWCA and one order of magnitude better than that of MODA, as shown in Table 15.

For DTLZ1 problem, it can be seen from Fig. 14 and Table 16 that MOICSO has better performance than MODA and MOWCA in both convergence and diversity.

For DTLZ2 problem, it can be seen from Fig. 16 and Table 17 that MOICSO outperforms the other two comparison algorithms in convergence. However, MODA obtains the best result among all the algorithms in diversity.

For DTLZ3 problem, it can be seen from Fig. 18 and Table 18 that MOICSO outperforms the other two comparison algorithms in both convergence and diversity. In particular, the performance of MOICSO in convergence metric is one order of magnitude better than that of MODA and MOWCA.