# Multi-Objective Particle Swarm Optimization with Multi-Archiving Strategy

Authors: Qian Zhang; Yanmin Liu; Huayao Han; Meilan Yang; Xiaoli Shu
Journal: Scientific Programming
 (2022)

---

## Abstract

Although the principle of multi-objective particle swarm optimization is simple and the operability is strong, it is still prone to local convergence and the convergence accuracy is not high. In order to solve the above problems, we propose a multi-objective particle swarm optimization algorithm based on multi strategies and archives. This algorithm is mainly divided into three important parts. Firstly, in the phase of sorting the optimal solutions, the solution set is stored in two different archives according to different conditions; secondly, in order to increase the diversity of the optimal solutions, several strategies are adopted in updating archives and maintaining archives’ scale. Finally, Gaussian perturbation strategy is applied to increase the distribution of particles and improve the quality of the optimal solution set. We compare the proposed algorithm with other algorithms and test it with different test indexes, Pareto graphs, and convergence graphs. The results show that this proposed algorithm has remarkable performance and the proposed method has advantages.


---

## Body

## 1. Introduction

With the continuous development of simulating computing, intelligent algorithms are advancing in a variety of directions. As the principle of bionic algorithm comes from natural biological behaviour and is easy to understand, a variety of algorithms emerge in the category of this kind of algorithm. It falls into two main types: one is an intelligent algorithm derived from the evolutionary behaviour of natural organisms. For example, Essam H. Houssein et al. [1] presented a multi-objective optimization algorithm based on the Slime mould algorithm (SMA) called multi-objective SMA (MOSMA); Yutao Yang et al. [2] proposed a general-purpose population-based optimization technique called Hunger Games Search (HGS) and so on. The other is an intelligent algorithm derived from the combination of the two algorithms according to the different characteristics of different populations. For example, Zar Bakht Imtiaz et al. [3] proposed the Multi-Layer Ant Colony Optimization (MLACO), Hongliang Zhang et al. [4] presented a chaotic slap swarm algorithm (SSA) with differential evolution (CDESSA), and so on. In addition, this kind of algorithm is mainly used to solve complex optimization problems in engineering field due to its strong searching ability, such as the electric field [5], the smart homes field [6], and so on.Therefore, swarm intelligence algorithm, as a branch of bionic algorithm, is also developing rapidly; among them, the development of particle swarm algorithm (PSO) [7] is particularly obvious. It is widely used in electric power [8], medical [9], satellite constellation design [10], and other industries, due to its characteristics of easy to understand and easy to operate. It was first proposed by Eberhart and Kennedy in 1995. Its principle is to find the optimal solutions of the first-order linear optimization model by simulating the behaviour of birds looking for food. However, with the dependent of the research, a novel PSO was proposed by Coello in 2004 [11], and applied to resolve multi-objective optimization problems. Therefore, the multi-objective particle swarm optimization (MOPSO) was produced. Up to now, the algorithm presents the characteristics of more and more diverse improvement methods and more and more significant improvement effects. The improvement methods of these algorithms are mainly divided into the following aspects: (1) Since the algorithm is prone to fall into local convergence, researchers help algorithms to get rid of this dilemma through various ways, such as improving mutation strategy, enhancing search strategy, and so on. In order to increase the quality of the optimal solution set in the decision space, Li Guoqing et al. [12] added a grid search strategy to search the particle population based on the particle position on the grid. They used K-means clustering strategy to locate the Pareto optimal solution set in the population, and then searched for high-quality solutions in the decision space by grid search. This method improves the exploration ability of the algorithm, but due to the limitation of K-means strategy, the parameter setting has an impact on the algorithm result. Qu Boyang et al. [13] proposed a grid-guided particle swarm optimizer for solving multimodal multi-objective optimization problems. The researchers divided the mesh adaptively, calculated the mesh density, and randomly deleted particles that exceeded the grid density. In addition, the inertia weight and acceleration coefficients are adjusted adaptively. These methods effectively alleviate the problem of local convergence, but greatly increase the computation of the algorithm. (2) To implement the particle swarm strategy, the population is divided into multiple subgroups, and the parallel operation of multiple subgroups is carried out to enhance the computing capability of the algorithm. Yunfeng Zhang et al. [14] divided the swarm adaptively into multiple subgroups and updated it with different strategies to increase the population diversity. PSO is adaptively divided into several subgroups, and the local optimal particles are found in the subgroups to guide the flight of ordinary particles. These methods can increase the diversity of particles to a certain extent. Zheng Jinhua et al. [15] proposed a dynamic MOPSO based on adversarial decomposition and neighbourhood evolution. The algorithm uses adversarial search direction to alternately update and coevolve the two populations of dynamic PSO, which increases the probability of finding the optimal solution set. (3) PSO is combined with other algorithms through the union of various different algorithms to improve the accuracy of the algorithm. A. Francis Saviour Devaraj et al. [16] combined the firefly algorithm (FA) with PSO to improve the speed and accuracy of the algorithm, so they select the best solution of the objective space by the distance between a point and a specified straight line. Kallol Biswas et al. [17] combined grey wolf optimization (GWO), cellular automata (CA) technique, and PSO to solve the optimization problems. The above two new algorithms combine the PSO with FA and GWO, respectively, to enhance their search ability and solving ability. Finally, the new algorithms are applied to engineering field to reduce resource loss. Based on the above theoretical researches, in order to further improve the searchability of MOPSO, alleviate the problem that the algorithm is easy to fall into local convergence, and improve its convergence accuracy. In this paper, we propose a MOPSO with multi-archiving strategy (MSMOPSO) which has two archives. The purpose is to improve the performance of MOPSO by improving the exploration ability of particles in the algorithm. Firstly, compared with MOPSO, the archive is expanded to increase the number of candidate solutions. Subsequently, the particles stored in archives are selected in different ways. For the first archive, the particles are selected by elite selection strategy, and this strategy is also applied to update this archive in the iterative process. For the second archive, we select particles based on self-organizing strategy, and then update this archive with species migration model. Finally, Gaussian mutation strategy is utilized in perturbing the positions of particles in archives, aiming to improve the solving accuracy of the algorithm by enhancing the diversity of particle population. For the latter, a physical migration model is used to update this archive.The rest of the article consists of the following sections. In the second part, the definition of multi-objective algorithms is introduced, and the prior work for this paper is described. The third part describes the improved MOPSO in detail. The fourth part describes the experimental setup. The fifth part reports the experimental conclusion and carries on the correlation analysis. The sixth part summarizes the article.


## 2. Related Work

### 2.1. Multi-Objective Particles Swarm Optimization

In MOPSO, the velocity and position of each particle are initialized in the solution space in a random way, then other particles search in the direction of global and individual optimality, so the whole search and update process follows the current optimal solutions to obtain the target. Assume that the total number of particles in the population isN, and the search space is D-dimensional. Particles’ positions are expressed as Xi=X1,X2,X3......XN; their velocity is Vi=V1,V2,V3......VN; the whole search process is as follows:(1)Vijt+1=wVijt+c1r1Pijt−Xijt+c2r2Pgjt−Xijt,Xijt+1=Xijt+Vijt+1.Among these,Vijt+1 and Xijt represent the velocity and position components of the ith dimension of the particle in the jth generation, respectively. w represents inertial weight, which enables particles to expand the search space, its value is set to 0.729. c1 and c2 represent social learning coefficients, with a value of 0.7029. r1 and r2 are random numbers within 0,1; Pijt and Pgjt represent the local and global optimal positions experienced by the particle in the jth dimensional component, respectively.Related definitions of multi-objective problems (MOPs) and Pareto font have been proposed in Ref. [18]. The typical feature of MOPSO is that it has multiple optimal value solutions. Therefore, it is necessary to add an external archive to store high-quality particles to form a candidate solution set. In addition, there are various methods to select the optimal solutions, but the core basis of selection is Pareto ranking value for these. If the particle sorting value is one, it proves that this particle has a relatively higher mass in the swarm. The overall process of the algorithm is as follows:Firstly, initialize particle swarm information; then, calculate the fitness values of the particles in the population, in order to determine the optimal position of the individual particles. Thirdly, fitness values of particles are sorted by Pareto ranking strategy, Subsequently, the noninferior solutions are stored in the external archive. At the same time, the particle position is determined according to the initial mesh partition of the objective function values. The fourth step is to use the roulette method to select the optimal particle from the external archive and take it as the global optimal position particle (gbest). Fifth, update the velocity and position of the particle, and the optimal position of the individual particle (pbest) is updated by using the dominance relation. Finally, repeat steps 2–5 to update the external archive, and retain the optimal noninferior solutions. This process ends when MOPSO reaches the iteration stop condition.


## 2.1. Multi-Objective Particles Swarm Optimization

In MOPSO, the velocity and position of each particle are initialized in the solution space in a random way, then other particles search in the direction of global and individual optimality, so the whole search and update process follows the current optimal solutions to obtain the target. Assume that the total number of particles in the population isN, and the search space is D-dimensional. Particles’ positions are expressed as Xi=X1,X2,X3......XN; their velocity is Vi=V1,V2,V3......VN; the whole search process is as follows:(1)Vijt+1=wVijt+c1r1Pijt−Xijt+c2r2Pgjt−Xijt,Xijt+1=Xijt+Vijt+1.Among these,Vijt+1 and Xijt represent the velocity and position components of the ith dimension of the particle in the jth generation, respectively. w represents inertial weight, which enables particles to expand the search space, its value is set to 0.729. c1 and c2 represent social learning coefficients, with a value of 0.7029. r1 and r2 are random numbers within 0,1; Pijt and Pgjt represent the local and global optimal positions experienced by the particle in the jth dimensional component, respectively.Related definitions of multi-objective problems (MOPs) and Pareto font have been proposed in Ref. [18]. The typical feature of MOPSO is that it has multiple optimal value solutions. Therefore, it is necessary to add an external archive to store high-quality particles to form a candidate solution set. In addition, there are various methods to select the optimal solutions, but the core basis of selection is Pareto ranking value for these. If the particle sorting value is one, it proves that this particle has a relatively higher mass in the swarm. The overall process of the algorithm is as follows:Firstly, initialize particle swarm information; then, calculate the fitness values of the particles in the population, in order to determine the optimal position of the individual particles. Thirdly, fitness values of particles are sorted by Pareto ranking strategy, Subsequently, the noninferior solutions are stored in the external archive. At the same time, the particle position is determined according to the initial mesh partition of the objective function values. The fourth step is to use the roulette method to select the optimal particle from the external archive and take it as the global optimal position particle (gbest). Fifth, update the velocity and position of the particle, and the optimal position of the individual particle (pbest) is updated by using the dominance relation. Finally, repeat steps 2–5 to update the external archive, and retain the optimal noninferior solutions. This process ends when MOPSO reaches the iteration stop condition.


## 3. Improved Multi-Objective Particle Swarm Optimization

### 3.1. Multiple Archives

Since the external archive is used to store relatively high-quality particles in MOPSO, the quality of the external archive will directly affect the quality of the optimal solution set. Furthermore, the external archive is stored as follows: particles with sorting value of one obtained by each iteration are stored in it, and then the archive is extended by using elite selection strategy. However, if these particles are too concentrated, the algorithm will inevitably fall into “premature” phenomenon. Therefore, in this part, we use different archiving methods to store potential solutions of particles according to different adaptive conditions.


### 3.2. The First Archive

In the initialization stage, Pareto ranking is carried out on particles first. Meanwhile, if there is only one particle whose sorting value is one in swarm, it will be stored in the first archive. In the iterative process, we compare the fitness values of the particles in the population and the fitness values of the particles in the archive and then adopt the elite selection strategy to select the superior particles. If there are particles like this, they will be stored in archive, and the inferior particles will be deleted from the archive. Through this way, we can achieve the purpose of archiving an update.


### 3.3. The Second Archive

During initialization, Pareto ranking is carried out on particles in the population first. Meanwhile, if there are multiple particles with sorting value of one, they will be stored in the second external archive. Then, in the iterative process, in order to increase the diversity of particles, these particles are regarded as the leaders. The specific operation of the second archive is to select particles around the leaders through the self-organizing grouping strategy. In other words, select particles with excellent performance as potential solutions and put them into the archive to achieve the purpose of archiving expansion. In addition, the physical migration model is used to update the archive to avoid excessive agglomeration of particles. The self-organizing grouping strategy and physical migration model are presented in section2 and section 2.


#### 3.3.1. Archive Expansion Based on Self-Organizing Strategy

A large number of studies show that during the process of multi-objective calculations, the decision variables can be divided into multiple subgroups according to certain rules, which can improve the efficiency of algorithms and the quality of solutions. Therefore, in this stage, we select high-quality particles in swarm through the form of subgroups. The commonly used clustering algorithm isK-means clustering algorithm [19], but the K-means algorithm has three obvious shortcomings: (1) Solutions of the K-means algorithm depends on the selection of the initial point; if the quality of these points are poor, its results will not reach the optimal; (2) The value of K in K-means algorithm needs to be set in advance, and the selection of K value has a great influence on final results; (3) During the K-means clustering, this clustering algorithm is prone to the appearance of empty groups, which makes it fall into local convergence. Thus, we utilize another strategy, namely, self-organizing strategy. This strategy refers to the formation of a certain structure of the objective variables spontaneously in accordance with certain rules, which are widely used in the multimodal MOPs [19]. In MSMOPSO, this strategy is applied to select particles in the population. In addition, when applying it, there is no need to set the number of groups in advance. What is more, it can group the swarm adaptively, which can effectively avoid the occurrence of empty group phenomenon.The clustering mechanism of self-organization strategy applied to the algorithm is as follows: firstly, the entire population is dynamically grouped according to leaders. If there are multiple leaders, the group will be divided according to the Euclidean distance between leaders and other particles in the population. At the same time, if a large number of duplicate particles are generated in the archive due to the proximity distance of two leaders, we will delete the second archive’s duplicate particles. The results show that when the group radius is 1/20, its effect is the best [20]. The specific situation is shown in Figure 1.Figure 1
Self-organizing clustering strategy.As shown in Figure1, take the two-dimensional test function as an example; the particles in the population are marked on the coordinate axis and the vertical axis coordinates represent particles’ fitness values. It is assumed that there are three particles whose sorting vale is one (marked by triangles in Figure 1) and then the Euclidean distance d between the particle in population and the triangle is calculated. If d is less than the previously set value, this particle and triangle belong to the same subgroup (the circular region in the diagram represents the dynamic group established around the triangle and the Euclidean distance). If the value of d is too large, the particle is redundant, and no grouping operation is performed on this particle (such as rhomboid particle in Figure 1). So, this particle was not assigned because the value of d was greater than the previously set value. In addition, if there are overlapping particles between subgroups, such as the quadrangle symbol in Figure 1, this symbol is stored directly in the archive without special grouping. The stored procedure is shown in Figure 2.Figure 2
Archive expansion based on self-organizing strategy.The rectangular table represents the second archive, and particles in it indicate that they are saved at this time. It is obvious from Figure2 that according to the self-organizing clustering strategy, all particles around leaders are stored.


#### 3.3.2. External Storage Update Based on Species Migration Model

Species migration model is used in Biogeography-Based Optimization, which is mainly to simulate the migration of species between multiple habitats. The experimental results show that this model has good convergence and is effective to solve the algorithm [21]. This model is mainly divided into two stages: species migration and species mutation. The migration rules followed are shown in Figure 3:Figure 3
Physical migration model.S stands for the number of species, P represents the migration rate of this habitat, λ represents the immigration rate of this habitat, µ represents the emigration rate of this habitat, E represents the maximum of emigration rate function, I represents the maximum of immigration rate function, S0 represents the equilibrium point, and Smax represents the maximum capacity of habitat, that is, the size of this species can be accommodated. According to Figure 2, when species quantity S is 0, the habitat emigration rate µ is 0, while the immigration rate λ reaches the maximum. When the species number reaches Smax, the species immigration rate λ is 0, and the species emigration rate µ reaches the maximum. When the species number is S0, the immigration rate equals to the emigration rate, indicating that the habitat is in equilibrium at this time, that is, the migration effect for species is the best, and the habitat after migration is of the best quality. λS represents the immigration rate under the equilibrium state, and µS represents the emigration rate under the equilibrium state. Therefore, λS and µS are adopted in the application of MSMOPSO. They are calculated as follows:(2)λS=I1−SSmax,μS=ESSmax.The meanings of related parameters are the same as those in Figure3.In MSMOPSO, the second archive is regarded as the habitat of organisms, particles stored in this archive are regarded as the species living in the habitat, and the initial update of this archive is completed by migrating the species in the habitat. In addition, the maximum value of immigration rate and emigration rate in the above process is 1. Therefore, based on the above formula, MSMOPSO adopts the following formula to calculate emigration rate and immigration rate based on species, so as to complete the preliminary update of the second archive. As shown in Figure3, the probability of habitat mutation is inversely proportional to the probability of species number. The specific form is shown in the following formula: Rk represents the migration rate of particles in the archive, Pk represents the immigration rate of particles in the archive, and the calculation formula is as follows:(3)Rk=1−kSmax,(4)Pk=kSmax.Smax represents the number of species in the archive; K represents the number of species traversed.Besides, species in the migration process cannot improve the diversity of species in the habitat if it only considers the immigration rate and emigration rate. Therefore, it is necessary to select a species in the migrated habitat for mutation disturbance, and then move the disturbed species into the habitat. Referring to this idea, MSMOPSO introduces random disturbance factor “a” at this stage. That is, some random particles in this archive are mutated and disturbed to produce a new incoming particle. This method can increase the diversity of particle population and increase the possibility of obtaining the optimal solution set. Refer to article [22], for the value of random disturbance factor. Therefore, the migration disturbance strategy is shown in the following formula:(5)Xij=aXij.a is the random number between (0, 1.2), Xij represents the particles disturbed by migration. This process is shown in Figure 4.Figure 4
External storage update based on species migration model.In Figure4, the solid line represents the species migration process and the dotted line represents the species disturbance process. The archive on the left side of Figure 4 shows that after the archive expansion (the process shown in Figure 3) is completed, the associated particles are removed based on the emigration rate of individuals. At the same time, according to the cumulative probability of particle emigration rate in the archive, some particles that need to be disturbed are determined; then the disturbed particles are randomly stored in the archive again. The right part of Figure 4 represents the archive after this process has been completed, from which two advantages of applying the species migration model can be clearly seen: first, the storage size within the archive is reduced; second, it increases the diversity of particles in this archive.


## 3.1. Multiple Archives

Since the external archive is used to store relatively high-quality particles in MOPSO, the quality of the external archive will directly affect the quality of the optimal solution set. Furthermore, the external archive is stored as follows: particles with sorting value of one obtained by each iteration are stored in it, and then the archive is extended by using elite selection strategy. However, if these particles are too concentrated, the algorithm will inevitably fall into “premature” phenomenon. Therefore, in this part, we use different archiving methods to store potential solutions of particles according to different adaptive conditions.


## 3.2. The First Archive

In the initialization stage, Pareto ranking is carried out on particles first. Meanwhile, if there is only one particle whose sorting value is one in swarm, it will be stored in the first archive. In the iterative process, we compare the fitness values of the particles in the population and the fitness values of the particles in the archive and then adopt the elite selection strategy to select the superior particles. If there are particles like this, they will be stored in archive, and the inferior particles will be deleted from the archive. Through this way, we can achieve the purpose of archiving an update.


## 3.3. The Second Archive

During initialization, Pareto ranking is carried out on particles in the population first. Meanwhile, if there are multiple particles with sorting value of one, they will be stored in the second external archive. Then, in the iterative process, in order to increase the diversity of particles, these particles are regarded as the leaders. The specific operation of the second archive is to select particles around the leaders through the self-organizing grouping strategy. In other words, select particles with excellent performance as potential solutions and put them into the archive to achieve the purpose of archiving expansion. In addition, the physical migration model is used to update the archive to avoid excessive agglomeration of particles. The self-organizing grouping strategy and physical migration model are presented in section2 and section 2.


### 3.3.1. Archive Expansion Based on Self-Organizing Strategy

A large number of studies show that during the process of multi-objective calculations, the decision variables can be divided into multiple subgroups according to certain rules, which can improve the efficiency of algorithms and the quality of solutions. Therefore, in this stage, we select high-quality particles in swarm through the form of subgroups. The commonly used clustering algorithm isK-means clustering algorithm [19], but the K-means algorithm has three obvious shortcomings: (1) Solutions of the K-means algorithm depends on the selection of the initial point; if the quality of these points are poor, its results will not reach the optimal; (2) The value of K in K-means algorithm needs to be set in advance, and the selection of K value has a great influence on final results; (3) During the K-means clustering, this clustering algorithm is prone to the appearance of empty groups, which makes it fall into local convergence. Thus, we utilize another strategy, namely, self-organizing strategy. This strategy refers to the formation of a certain structure of the objective variables spontaneously in accordance with certain rules, which are widely used in the multimodal MOPs [19]. In MSMOPSO, this strategy is applied to select particles in the population. In addition, when applying it, there is no need to set the number of groups in advance. What is more, it can group the swarm adaptively, which can effectively avoid the occurrence of empty group phenomenon.The clustering mechanism of self-organization strategy applied to the algorithm is as follows: firstly, the entire population is dynamically grouped according to leaders. If there are multiple leaders, the group will be divided according to the Euclidean distance between leaders and other particles in the population. At the same time, if a large number of duplicate particles are generated in the archive due to the proximity distance of two leaders, we will delete the second archive’s duplicate particles. The results show that when the group radius is 1/20, its effect is the best [20]. The specific situation is shown in Figure 1.Figure 1
Self-organizing clustering strategy.As shown in Figure1, take the two-dimensional test function as an example; the particles in the population are marked on the coordinate axis and the vertical axis coordinates represent particles’ fitness values. It is assumed that there are three particles whose sorting vale is one (marked by triangles in Figure 1) and then the Euclidean distance d between the particle in population and the triangle is calculated. If d is less than the previously set value, this particle and triangle belong to the same subgroup (the circular region in the diagram represents the dynamic group established around the triangle and the Euclidean distance). If the value of d is too large, the particle is redundant, and no grouping operation is performed on this particle (such as rhomboid particle in Figure 1). So, this particle was not assigned because the value of d was greater than the previously set value. In addition, if there are overlapping particles between subgroups, such as the quadrangle symbol in Figure 1, this symbol is stored directly in the archive without special grouping. The stored procedure is shown in Figure 2.Figure 2
Archive expansion based on self-organizing strategy.The rectangular table represents the second archive, and particles in it indicate that they are saved at this time. It is obvious from Figure2 that according to the self-organizing clustering strategy, all particles around leaders are stored.


### 3.3.2. External Storage Update Based on Species Migration Model

Species migration model is used in Biogeography-Based Optimization, which is mainly to simulate the migration of species between multiple habitats. The experimental results show that this model has good convergence and is effective to solve the algorithm [21]. This model is mainly divided into two stages: species migration and species mutation. The migration rules followed are shown in Figure 3:Figure 3
Physical migration model.S stands for the number of species, P represents the migration rate of this habitat, λ represents the immigration rate of this habitat, µ represents the emigration rate of this habitat, E represents the maximum of emigration rate function, I represents the maximum of immigration rate function, S0 represents the equilibrium point, and Smax represents the maximum capacity of habitat, that is, the size of this species can be accommodated. According to Figure 2, when species quantity S is 0, the habitat emigration rate µ is 0, while the immigration rate λ reaches the maximum. When the species number reaches Smax, the species immigration rate λ is 0, and the species emigration rate µ reaches the maximum. When the species number is S0, the immigration rate equals to the emigration rate, indicating that the habitat is in equilibrium at this time, that is, the migration effect for species is the best, and the habitat after migration is of the best quality. λS represents the immigration rate under the equilibrium state, and µS represents the emigration rate under the equilibrium state. Therefore, λS and µS are adopted in the application of MSMOPSO. They are calculated as follows:(2)λS=I1−SSmax,μS=ESSmax.The meanings of related parameters are the same as those in Figure3.In MSMOPSO, the second archive is regarded as the habitat of organisms, particles stored in this archive are regarded as the species living in the habitat, and the initial update of this archive is completed by migrating the species in the habitat. In addition, the maximum value of immigration rate and emigration rate in the above process is 1. Therefore, based on the above formula, MSMOPSO adopts the following formula to calculate emigration rate and immigration rate based on species, so as to complete the preliminary update of the second archive. As shown in Figure3, the probability of habitat mutation is inversely proportional to the probability of species number. The specific form is shown in the following formula: Rk represents the migration rate of particles in the archive, Pk represents the immigration rate of particles in the archive, and the calculation formula is as follows:(3)Rk=1−kSmax,(4)Pk=kSmax.Smax represents the number of species in the archive; K represents the number of species traversed.Besides, species in the migration process cannot improve the diversity of species in the habitat if it only considers the immigration rate and emigration rate. Therefore, it is necessary to select a species in the migrated habitat for mutation disturbance, and then move the disturbed species into the habitat. Referring to this idea, MSMOPSO introduces random disturbance factor “a” at this stage. That is, some random particles in this archive are mutated and disturbed to produce a new incoming particle. This method can increase the diversity of particle population and increase the possibility of obtaining the optimal solution set. Refer to article [22], for the value of random disturbance factor. Therefore, the migration disturbance strategy is shown in the following formula:(5)Xij=aXij.a is the random number between (0, 1.2), Xij represents the particles disturbed by migration. This process is shown in Figure 4.Figure 4
External storage update based on species migration model.In Figure4, the solid line represents the species migration process and the dotted line represents the species disturbance process. The archive on the left side of Figure 4 shows that after the archive expansion (the process shown in Figure 3) is completed, the associated particles are removed based on the emigration rate of individuals. At the same time, according to the cumulative probability of particle emigration rate in the archive, some particles that need to be disturbed are determined; then the disturbed particles are randomly stored in the archive again. The right part of Figure 4 represents the archive after this process has been completed, from which two advantages of applying the species migration model can be clearly seen: first, the storage size within the archive is reduced; second, it increases the diversity of particles in this archive.


## 3.3.1. Archive Expansion Based on Self-Organizing Strategy

A large number of studies show that during the process of multi-objective calculations, the decision variables can be divided into multiple subgroups according to certain rules, which can improve the efficiency of algorithms and the quality of solutions. Therefore, in this stage, we select high-quality particles in swarm through the form of subgroups. The commonly used clustering algorithm isK-means clustering algorithm [19], but the K-means algorithm has three obvious shortcomings: (1) Solutions of the K-means algorithm depends on the selection of the initial point; if the quality of these points are poor, its results will not reach the optimal; (2) The value of K in K-means algorithm needs to be set in advance, and the selection of K value has a great influence on final results; (3) During the K-means clustering, this clustering algorithm is prone to the appearance of empty groups, which makes it fall into local convergence. Thus, we utilize another strategy, namely, self-organizing strategy. This strategy refers to the formation of a certain structure of the objective variables spontaneously in accordance with certain rules, which are widely used in the multimodal MOPs [19]. In MSMOPSO, this strategy is applied to select particles in the population. In addition, when applying it, there is no need to set the number of groups in advance. What is more, it can group the swarm adaptively, which can effectively avoid the occurrence of empty group phenomenon.The clustering mechanism of self-organization strategy applied to the algorithm is as follows: firstly, the entire population is dynamically grouped according to leaders. If there are multiple leaders, the group will be divided according to the Euclidean distance between leaders and other particles in the population. At the same time, if a large number of duplicate particles are generated in the archive due to the proximity distance of two leaders, we will delete the second archive’s duplicate particles. The results show that when the group radius is 1/20, its effect is the best [20]. The specific situation is shown in Figure 1.Figure 1
Self-organizing clustering strategy.As shown in Figure1, take the two-dimensional test function as an example; the particles in the population are marked on the coordinate axis and the vertical axis coordinates represent particles’ fitness values. It is assumed that there are three particles whose sorting vale is one (marked by triangles in Figure 1) and then the Euclidean distance d between the particle in population and the triangle is calculated. If d is less than the previously set value, this particle and triangle belong to the same subgroup (the circular region in the diagram represents the dynamic group established around the triangle and the Euclidean distance). If the value of d is too large, the particle is redundant, and no grouping operation is performed on this particle (such as rhomboid particle in Figure 1). So, this particle was not assigned because the value of d was greater than the previously set value. In addition, if there are overlapping particles between subgroups, such as the quadrangle symbol in Figure 1, this symbol is stored directly in the archive without special grouping. The stored procedure is shown in Figure 2.Figure 2
Archive expansion based on self-organizing strategy.The rectangular table represents the second archive, and particles in it indicate that they are saved at this time. It is obvious from Figure2 that according to the self-organizing clustering strategy, all particles around leaders are stored.


## 3.3.2. External Storage Update Based on Species Migration Model

Species migration model is used in Biogeography-Based Optimization, which is mainly to simulate the migration of species between multiple habitats. The experimental results show that this model has good convergence and is effective to solve the algorithm [21]. This model is mainly divided into two stages: species migration and species mutation. The migration rules followed are shown in Figure 3:Figure 3
Physical migration model.S stands for the number of species, P represents the migration rate of this habitat, λ represents the immigration rate of this habitat, µ represents the emigration rate of this habitat, E represents the maximum of emigration rate function, I represents the maximum of immigration rate function, S0 represents the equilibrium point, and Smax represents the maximum capacity of habitat, that is, the size of this species can be accommodated. According to Figure 2, when species quantity S is 0, the habitat emigration rate µ is 0, while the immigration rate λ reaches the maximum. When the species number reaches Smax, the species immigration rate λ is 0, and the species emigration rate µ reaches the maximum. When the species number is S0, the immigration rate equals to the emigration rate, indicating that the habitat is in equilibrium at this time, that is, the migration effect for species is the best, and the habitat after migration is of the best quality. λS represents the immigration rate under the equilibrium state, and µS represents the emigration rate under the equilibrium state. Therefore, λS and µS are adopted in the application of MSMOPSO. They are calculated as follows:(2)λS=I1−SSmax,μS=ESSmax.The meanings of related parameters are the same as those in Figure3.In MSMOPSO, the second archive is regarded as the habitat of organisms, particles stored in this archive are regarded as the species living in the habitat, and the initial update of this archive is completed by migrating the species in the habitat. In addition, the maximum value of immigration rate and emigration rate in the above process is 1. Therefore, based on the above formula, MSMOPSO adopts the following formula to calculate emigration rate and immigration rate based on species, so as to complete the preliminary update of the second archive. As shown in Figure3, the probability of habitat mutation is inversely proportional to the probability of species number. The specific form is shown in the following formula: Rk represents the migration rate of particles in the archive, Pk represents the immigration rate of particles in the archive, and the calculation formula is as follows:(3)Rk=1−kSmax,(4)Pk=kSmax.Smax represents the number of species in the archive; K represents the number of species traversed.Besides, species in the migration process cannot improve the diversity of species in the habitat if it only considers the immigration rate and emigration rate. Therefore, it is necessary to select a species in the migrated habitat for mutation disturbance, and then move the disturbed species into the habitat. Referring to this idea, MSMOPSO introduces random disturbance factor “a” at this stage. That is, some random particles in this archive are mutated and disturbed to produce a new incoming particle. This method can increase the diversity of particle population and increase the possibility of obtaining the optimal solution set. Refer to article [22], for the value of random disturbance factor. Therefore, the migration disturbance strategy is shown in the following formula:(5)Xij=aXij.a is the random number between (0, 1.2), Xij represents the particles disturbed by migration. This process is shown in Figure 4.Figure 4
External storage update based on species migration model.In Figure4, the solid line represents the species migration process and the dotted line represents the species disturbance process. The archive on the left side of Figure 4 shows that after the archive expansion (the process shown in Figure 3) is completed, the associated particles are removed based on the emigration rate of individuals. At the same time, according to the cumulative probability of particle emigration rate in the archive, some particles that need to be disturbed are determined; then the disturbed particles are randomly stored in the archive again. The right part of Figure 4 represents the archive after this process has been completed, from which two advantages of applying the species migration model can be clearly seen: first, the storage size within the archive is reduced; second, it increases the diversity of particles in this archive.


## 4. Archive Scale Maintenance

### 4.1. First Archive Scale Maintenance

In the first archive, if the archive size exceeds the limit, its particles are placed in an adaptive grid, and then the archive size is maintained according to the density of the particles in the adaptive grid. That is, if the particles in the grid reach a preset value of maximum density (the maximum density set in MSMOPSO is 10), firstly, judge whether there are duplicate lines, if so, then perform the deduplication eliminatiny, otherwise, conduct deletion operation for the particles in it.


### 4.2. Second Archive Scale Maintenance

In the second archive phase, if the archive size exceeds the limit, the archive size is reduced in two ways. One way is to move some particles out of this archive through the emigration rate in the physical migration model. Specifically, a random number is generated to compare with the corresponding emigration rate of particles within the habitat. If the former is less than the latter, it means that this particle needs to be emigrated. The second method is to randomly delete some particles in this archive, if the size of habitat still exceeds the limit after species migration, and this process continues until the size of this archive does not exceed the set value. The specific process is shown as the pseudo code in Table1:Table 1
The pseudocode of the second archive.
Algorithm: The pseudocode of the second archive1. %gleaders represent the array of storing particles whose sort value is 12. %k represents the traversed particles in the second archive3. %Sbest represents an array generated by calculating the fitness values of particles in the population4. %Sbest1 represents the updated second archive5. %Max represents the population size of particles in Sbest6. %Max1 represents the population size of particles in Sbest17. %repsize represents the preset maximum size of the archive8. %Rk represents the particle migration rate9. %Pk represents the particle immigration rate10. %rand1 and rand2 represent the random number arrays11.  if length (gleaders (:,1)) > 112.   Apply the self-organizing policy to extend the second archive according to Figure1.13.   fork = 1:Max14.   Calculate the particle migration rateRk in the second archive according to formula (3).15.    Calculate the particle immigration ratePk in the second archive according to formula (4)16.    Generate random number arraysrand1 and rand217.    Ifrand1 (k) < Rk18.     Delete the particlek in the archive19.    end20.    Ifrand2 (k) < cumsum (Pk)21.     Perturbed the particlek and produced the new particle according to formula (5).22.     In the second archive, replace the particlek with this new particle23.    end24.   end25.   Generate the updated second archive,Sbest126.   Ifmax1 > repsize27.    Delete the extra parts28.   endThe significance of bold values is presented in the form of “ %+X” in Table 1.


### 4.3. Mutation Strategy

The purpose of mutation strategy is to perturb the positions of a certain number of particles, so as to increase the possibility of these particles finding unexplored regions, and search for a relatively high-quality solution set in these regions. Gaussian mutation strategy has strong global searchability, so it is helpful to alleviate the phenomenon of local convergence. Therefore, MSMOPSO uses Gaussian perturbation mutation strategy to mutate the position information of particles in order to increase the probability of obtaining optimal solutions. In this stage, firstly, a nonlinear function is used to determine the number of variant particles in the archive, and its number varies with the number of iterations. Subsequently, the number of mutation of each particle was determined adaptively according to the number of iterations, and the number of mutation was guaranteed not to exceed the preset value (i.e., the maximum and minimum limit of the number of mutation). Assuming that the number of iterations at this time ist, the number of mutant particles and the number of mutations required by the mutant particles are shown in the following formula, which is similar to Ref. [23].(6)nump=1−ttmax5/mutrate∗lengthposition,1,npt=unp−lnp∗1−2∗ttmax2+lnpunp−lnp0≤t≤tmax2,unp−lnp∗2∗t−maxtmax2+lnpunp−lnptmax2≤t≤tmax.In the above formulae,t represents the current number of iterations; tmax represents the maximum number of iterations in MSMOPSO; mutrate is the rate of variation with respect to the number of particles; nump represents the number of mutated particles; length (position,1) is the number of particles in the population; unp represents the maximum number of disturbances of the particle; lnp represents the minimum perturbation of the particle; np refers to the number of disturbances per variant particle. In addition to increase the randomness of the mutation, the variation dimension of the particle is determined randomly. Similar to Ref. [24], in dimensions to be changed, the position of components of the decision variable is changed by combining the upper and lower proportions of perturbation and Gaussian variation. It is nothing that in MSMOPSO, ∆d is regarded as the perturbation range of component position, that is, the variation range of particles is limited, so as to avoid producing too many particles that are too far away from the variation particles. The calculation formulae are as follows:(7)r2=Gaussian0,19,ld=rld∗Xjmin,ud=rud∗Xjmax,Δdr2=ud−ld∗2∗r22+ldud−ld0≤r2≤0.5,ud−ld∗1−2∗r2−12+ldud−ld0.5≤r2≤1,Xij=Xij+Δdr2.r2 represents the random number generated through Gaussian distribution, but the value is limited here, that is, the value cannot exceed 1; rld refers to the lower limit of variation; rud refers to the upper limit of variation; Xjmin is the lower bound of the jth dimension of particle X; Xjmax represents the upper bound of the jth dimension of particle X; ld refers to the disturbance lower bound of the decision component in the disturbance process; ud refers to the disturbance upper bound of the decision component in the disturbance process; ∆d represents the perturbation step term of the decision variable. In addition, in order to prevent the position of the disturbed component beyond the set boundary, it is necessary to reset the component beyond the boundary.According to the above part, MSMOPSO uses two archiving methods to store the better particles. Then, in the iteration process, the archived particles are updated, and these archived particles are applied to guide the flight direction and step size of other particles in the population. Therefore, the speed and position of particles are updated through the following formulae in MSMOPSO:(8)Vijt+1=wVijt+c1r1pbestijt−Xijt+c2r2repijt−Xijt,Xijt+1=Xijt+Vijt+1.In the iteration process, if the Pareto ranking value of the archived particle increases, it indicates that the particle performs better than other particles after updating (archived updating) and mutation (mutation strategy). Such particles are regarded aspbest particles, which guide other particles in flight, rep stands for particles in the archive.


### 4.4. Basic Framework of the Algorithm

MOPSO puts potential solutions in a single archive, so the diversity of solutions is limited. This algorithm improves the diversity of candidate solution set by adding archiving strategy. Its flow steps are as follows:Step 1.
Initialize the particle population, and set related parameters;Step 2.
Calculate the position and velocity of particles;Step 3.
Calculate the fitness value of particles in the population, and perform Pareto ranking on the fitness value;Step 4.
Mark the particles whose sorting value is one, and select pbest according to the fitness value of the particles;Step 5.
Determine whether there are multiple particles with a sorting value of one. If there are more than one particle, they will be stored in the second archive; otherwise, the particle is stored in the first archive;Step 6.
Update the first and second archives, respectively. The first archive uses the elite selection strategy to select related particles in the particle population. In the second archive, the physical migration model is used to carry out the in-migration, out-migration, and mutation operations for the particles;Step 7.
Select gbest particles in the corresponding archiving strategy, and guide other particles to carry out flight operations;Step 8.
Judge whether the algorithm reaches the iteration stop condition at this time. If the iteration stop condition is reached, the operation of the whole algorithm will be terminated. Otherwise, the mutation strategy is applied to the particles in the archive and step2 is returned.
The flow framework of the algorithm is shown in Figure5:Figure 5
The flow chart of MSMOPSO.


## 4.1. First Archive Scale Maintenance

In the first archive, if the archive size exceeds the limit, its particles are placed in an adaptive grid, and then the archive size is maintained according to the density of the particles in the adaptive grid. That is, if the particles in the grid reach a preset value of maximum density (the maximum density set in MSMOPSO is 10), firstly, judge whether there are duplicate lines, if so, then perform the deduplication eliminatiny, otherwise, conduct deletion operation for the particles in it.


## 4.2. Second Archive Scale Maintenance

In the second archive phase, if the archive size exceeds the limit, the archive size is reduced in two ways. One way is to move some particles out of this archive through the emigration rate in the physical migration model. Specifically, a random number is generated to compare with the corresponding emigration rate of particles within the habitat. If the former is less than the latter, it means that this particle needs to be emigrated. The second method is to randomly delete some particles in this archive, if the size of habitat still exceeds the limit after species migration, and this process continues until the size of this archive does not exceed the set value. The specific process is shown as the pseudo code in Table1:Table 1
The pseudocode of the second archive.
Algorithm: The pseudocode of the second archive1. %gleaders represent the array of storing particles whose sort value is 12. %k represents the traversed particles in the second archive3. %Sbest represents an array generated by calculating the fitness values of particles in the population4. %Sbest1 represents the updated second archive5. %Max represents the population size of particles in Sbest6. %Max1 represents the population size of particles in Sbest17. %repsize represents the preset maximum size of the archive8. %Rk represents the particle migration rate9. %Pk represents the particle immigration rate10. %rand1 and rand2 represent the random number arrays11.  if length (gleaders (:,1)) > 112.   Apply the self-organizing policy to extend the second archive according to Figure1.13.   fork = 1:Max14.   Calculate the particle migration rateRk in the second archive according to formula (3).15.    Calculate the particle immigration ratePk in the second archive according to formula (4)16.    Generate random number arraysrand1 and rand217.    Ifrand1 (k) < Rk18.     Delete the particlek in the archive19.    end20.    Ifrand2 (k) < cumsum (Pk)21.     Perturbed the particlek and produced the new particle according to formula (5).22.     In the second archive, replace the particlek with this new particle23.    end24.   end25.   Generate the updated second archive,Sbest126.   Ifmax1 > repsize27.    Delete the extra parts28.   endThe significance of bold values is presented in the form of “ %+X” in Table 1.


## 4.3. Mutation Strategy

The purpose of mutation strategy is to perturb the positions of a certain number of particles, so as to increase the possibility of these particles finding unexplored regions, and search for a relatively high-quality solution set in these regions. Gaussian mutation strategy has strong global searchability, so it is helpful to alleviate the phenomenon of local convergence. Therefore, MSMOPSO uses Gaussian perturbation mutation strategy to mutate the position information of particles in order to increase the probability of obtaining optimal solutions. In this stage, firstly, a nonlinear function is used to determine the number of variant particles in the archive, and its number varies with the number of iterations. Subsequently, the number of mutation of each particle was determined adaptively according to the number of iterations, and the number of mutation was guaranteed not to exceed the preset value (i.e., the maximum and minimum limit of the number of mutation). Assuming that the number of iterations at this time ist, the number of mutant particles and the number of mutations required by the mutant particles are shown in the following formula, which is similar to Ref. [23].(6)nump=1−ttmax5/mutrate∗lengthposition,1,npt=unp−lnp∗1−2∗ttmax2+lnpunp−lnp0≤t≤tmax2,unp−lnp∗2∗t−maxtmax2+lnpunp−lnptmax2≤t≤tmax.In the above formulae,t represents the current number of iterations; tmax represents the maximum number of iterations in MSMOPSO; mutrate is the rate of variation with respect to the number of particles; nump represents the number of mutated particles; length (position,1) is the number of particles in the population; unp represents the maximum number of disturbances of the particle; lnp represents the minimum perturbation of the particle; np refers to the number of disturbances per variant particle. In addition to increase the randomness of the mutation, the variation dimension of the particle is determined randomly. Similar to Ref. [24], in dimensions to be changed, the position of components of the decision variable is changed by combining the upper and lower proportions of perturbation and Gaussian variation. It is nothing that in MSMOPSO, ∆d is regarded as the perturbation range of component position, that is, the variation range of particles is limited, so as to avoid producing too many particles that are too far away from the variation particles. The calculation formulae are as follows:(7)r2=Gaussian0,19,ld=rld∗Xjmin,ud=rud∗Xjmax,Δdr2=ud−ld∗2∗r22+ldud−ld0≤r2≤0.5,ud−ld∗1−2∗r2−12+ldud−ld0.5≤r2≤1,Xij=Xij+Δdr2.r2 represents the random number generated through Gaussian distribution, but the value is limited here, that is, the value cannot exceed 1; rld refers to the lower limit of variation; rud refers to the upper limit of variation; Xjmin is the lower bound of the jth dimension of particle X; Xjmax represents the upper bound of the jth dimension of particle X; ld refers to the disturbance lower bound of the decision component in the disturbance process; ud refers to the disturbance upper bound of the decision component in the disturbance process; ∆d represents the perturbation step term of the decision variable. In addition, in order to prevent the position of the disturbed component beyond the set boundary, it is necessary to reset the component beyond the boundary.According to the above part, MSMOPSO uses two archiving methods to store the better particles. Then, in the iteration process, the archived particles are updated, and these archived particles are applied to guide the flight direction and step size of other particles in the population. Therefore, the speed and position of particles are updated through the following formulae in MSMOPSO:(8)Vijt+1=wVijt+c1r1pbestijt−Xijt+c2r2repijt−Xijt,Xijt+1=Xijt+Vijt+1.In the iteration process, if the Pareto ranking value of the archived particle increases, it indicates that the particle performs better than other particles after updating (archived updating) and mutation (mutation strategy). Such particles are regarded aspbest particles, which guide other particles in flight, rep stands for particles in the archive.


## 4.4. Basic Framework of the Algorithm

MOPSO puts potential solutions in a single archive, so the diversity of solutions is limited. This algorithm improves the diversity of candidate solution set by adding archiving strategy. Its flow steps are as follows:Step 1.
Initialize the particle population, and set related parameters;Step 2.
Calculate the position and velocity of particles;Step 3.
Calculate the fitness value of particles in the population, and perform Pareto ranking on the fitness value;Step 4.
Mark the particles whose sorting value is one, and select pbest according to the fitness value of the particles;Step 5.
Determine whether there are multiple particles with a sorting value of one. If there are more than one particle, they will be stored in the second archive; otherwise, the particle is stored in the first archive;Step 6.
Update the first and second archives, respectively. The first archive uses the elite selection strategy to select related particles in the particle population. In the second archive, the physical migration model is used to carry out the in-migration, out-migration, and mutation operations for the particles;Step 7.
Select gbest particles in the corresponding archiving strategy, and guide other particles to carry out flight operations;Step 8.
Judge whether the algorithm reaches the iteration stop condition at this time. If the iteration stop condition is reached, the operation of the whole algorithm will be terminated. Otherwise, the mutation strategy is applied to the particles in the archive and step2 is returned.
The flow framework of the algorithm is shown in Figure5:Figure 5
The flow chart of MSMOPSO.


## 5. Experimental Settings

In this experiment, MATLAB (R2020b) is used to test MSMOPSO; other simulations were run on Windows 10 platform using Intel i7-7700, 3.60 GHz processor with 16 GB RAM. Other test indicators are set as follows.


### 5.1. Experimental Comparison Algorithm and Comparison Function

We compared the experimental results with MOPSO and MOEA series algorithms. MOPSO series (MOPSOs) include: MOPSO based on clustering strategy: MPSOD [23], SMPSO [25], dMOPSO [26], and NMPSO [27]; classical MOPSO :  MOPSOCD [28] and CMOSPO [29]. MOEA series algorithms (MOEAs): decomposition-based MOEA : MOEAD [30]; classical MOEAs : ARMOEA [31], NSGAIISDR [32], RVEA [33], GrEA [34], and gNSGAII [35]. These comparison algorithms were iterated on ZDT1-4, ZDT6 [36], UF1-10 [37], and DTLZ1-7 [38] test functions for 10000 times, respectively. Meanwhile, these algorithms were run for 30 times, respectively. The code of the comparison algorithm was provided by PlatEMO [39]. Then, the performance of these algorithms is compared and analyzed by comprehensive evaluation indexes. See Section 5.2 for evaluation indicators. Then, the mean and standard deviation of the evaluation indexes in 30 runs are obtained to comprehensively demonstrate the superiority of the improved algorithm from the statistical aspect. The related parameter settings are shown in Table 2, where D represents the dimension of the given problems, N means the group size, FEs signifies fitness evolution number, and Problems means test problems.Table 2
Experimental parameter settings.
ProblemsDNFEsZDT1-ZDT43020010000ZDT63020010000DTLZ1-DTLZ71210010000UF1-UF73020010000UF8-UF103020010000To fairly test the characteristics of algorithms, each algorithm uses the same parameter settings when testing the benchmark function. As shown in Table3, the parameters for each algorithm are set to their default values. As Reported by Arcuri and Fraser (2013), using default algorithm parameter values effectively reduces the risk of better parameterization bias [12].Table 3
Algorithms’ parameter settings for MSMOPSO and other algorithms.
AlgorithmsParameters settingCommon settingsN:See Table.2Maxiteration=200Number of independent runs=30SMOPSOw∈0.1,0.5,c1,c2∈1.5,2.5,pm=1/n,ηm=20MPSODw∈0.1,0.5,c1,c2∈1.5,2.5,pc=0.9,F=0.5,CR=0.5,pm=1/n,ηm=20MOPSOCDw∈0.1,0.5,c1,c2∈1.5,2.5NMPSOw∈0.1,0.5,c1,c2,c3∈1.5,2.5,pm=1/n,ηc=ηm=20CMOPSOr=10dMOPSOw∈0.1,0.5,c1,c2∈1.5,2.5ARMOEApc=1.0,pm=1/n,ηc=20,ηm=20MOEADpc=1.0,pm=1/n,ηc=20,ηm=20NSGAIISDRpc=1.0,pm=1/n,ηc=20,ηm=20gNSGAIIpc=1.0,pm=1/n,ηc=20,ηm=20RVEAα=2,fr=0.1GrEApc=1.0,pm=1/n,ηc=20,ηm=20,div=9MSMOPSOc1=2,c2=2,w=0.4,unp=5;lnp=1,rud=0.9,rld=0.02,α∈0,1.2


### 5.2. Comprehensive Evaluation Indexes

In this algorithm, Inverted Generation Distance (IGD) [40] evaluation index and Hyper Volume (HV) [41] evaluation index are adopted to comprehensively evaluate the performance of the algorithm.


### 5.3. IGD Evaluation Index

This algorithm adopts IGD evaluation index, which evaluates the convergence and distribution performance of the algorithm mainly by calculating the minimum distance sum between each individual on the real Pareto front surface and the individual set obtained by the algorithm. The smaller the value is, the better the convergence and distribution performance of the algorithm is. IGD evaluation index is calculated as follows:(9)IGDP,Q=∑v∈Pdv,QP.Among them, theP is for uniform distribution on the true Pareto point set, P is for distribution in the true Pareto set of points on the surface of the individual number. Q is the Pareto optimal solution set obtained by the algorithm and dv,Q is the minimum Euclidean distance between individual V and population Q in P.


### 5.4. HV Evaluation Index

HV evaluation index refers to the hyper volume index, that is, the volume of the region in the target space surrounded by the nondominated solution set and the reference points obtained by the algorithm. This index can simultaneously evaluate the convergence and distribution of the solution set in space. Convergence refers to the approximation degree between the solution set and the real Pareto front, distribution refers to the wide distribution of the whole solution set in the target space. A larger HV value indicates that the algorithm can dominate a larger area in the target space and has better diversity. The formula of HV index is as follows:(10)HV=δUi=1Svi,where δ represents the Lebesgue measure, which is used to measure volume. S is the number of nondominated solution set, vi represents the supervolume consisting of the reference point and the ith solution in the nondominant solution set.


## 5.1. Experimental Comparison Algorithm and Comparison Function

We compared the experimental results with MOPSO and MOEA series algorithms. MOPSO series (MOPSOs) include: MOPSO based on clustering strategy: MPSOD [23], SMPSO [25], dMOPSO [26], and NMPSO [27]; classical MOPSO :  MOPSOCD [28] and CMOSPO [29]. MOEA series algorithms (MOEAs): decomposition-based MOEA : MOEAD [30]; classical MOEAs : ARMOEA [31], NSGAIISDR [32], RVEA [33], GrEA [34], and gNSGAII [35]. These comparison algorithms were iterated on ZDT1-4, ZDT6 [36], UF1-10 [37], and DTLZ1-7 [38] test functions for 10000 times, respectively. Meanwhile, these algorithms were run for 30 times, respectively. The code of the comparison algorithm was provided by PlatEMO [39]. Then, the performance of these algorithms is compared and analyzed by comprehensive evaluation indexes. See Section 5.2 for evaluation indicators. Then, the mean and standard deviation of the evaluation indexes in 30 runs are obtained to comprehensively demonstrate the superiority of the improved algorithm from the statistical aspect. The related parameter settings are shown in Table 2, where D represents the dimension of the given problems, N means the group size, FEs signifies fitness evolution number, and Problems means test problems.Table 2
Experimental parameter settings.
ProblemsDNFEsZDT1-ZDT43020010000ZDT63020010000DTLZ1-DTLZ71210010000UF1-UF73020010000UF8-UF103020010000To fairly test the characteristics of algorithms, each algorithm uses the same parameter settings when testing the benchmark function. As shown in Table3, the parameters for each algorithm are set to their default values. As Reported by Arcuri and Fraser (2013), using default algorithm parameter values effectively reduces the risk of better parameterization bias [12].Table 3
Algorithms’ parameter settings for MSMOPSO and other algorithms.
AlgorithmsParameters settingCommon settingsN:See Table.2Maxiteration=200Number of independent runs=30SMOPSOw∈0.1,0.5,c1,c2∈1.5,2.5,pm=1/n,ηm=20MPSODw∈0.1,0.5,c1,c2∈1.5,2.5,pc=0.9,F=0.5,CR=0.5,pm=1/n,ηm=20MOPSOCDw∈0.1,0.5,c1,c2∈1.5,2.5NMPSOw∈0.1,0.5,c1,c2,c3∈1.5,2.5,pm=1/n,ηc=ηm=20CMOPSOr=10dMOPSOw∈0.1,0.5,c1,c2∈1.5,2.5ARMOEApc=1.0,pm=1/n,ηc=20,ηm=20MOEADpc=1.0,pm=1/n,ηc=20,ηm=20NSGAIISDRpc=1.0,pm=1/n,ηc=20,ηm=20gNSGAIIpc=1.0,pm=1/n,ηc=20,ηm=20RVEAα=2,fr=0.1GrEApc=1.0,pm=1/n,ηc=20,ηm=20,div=9MSMOPSOc1=2,c2=2,w=0.4,unp=5;lnp=1,rud=0.9,rld=0.02,α∈0,1.2


## 5.2. Comprehensive Evaluation Indexes

In this algorithm, Inverted Generation Distance (IGD) [40] evaluation index and Hyper Volume (HV) [41] evaluation index are adopted to comprehensively evaluate the performance of the algorithm.


## 5.3. IGD Evaluation Index

This algorithm adopts IGD evaluation index, which evaluates the convergence and distribution performance of the algorithm mainly by calculating the minimum distance sum between each individual on the real Pareto front surface and the individual set obtained by the algorithm. The smaller the value is, the better the convergence and distribution performance of the algorithm is. IGD evaluation index is calculated as follows:(9)IGDP,Q=∑v∈Pdv,QP.Among them, theP is for uniform distribution on the true Pareto point set, P is for distribution in the true Pareto set of points on the surface of the individual number. Q is the Pareto optimal solution set obtained by the algorithm and dv,Q is the minimum Euclidean distance between individual V and population Q in P.


## 5.4. HV Evaluation Index

HV evaluation index refers to the hyper volume index, that is, the volume of the region in the target space surrounded by the nondominated solution set and the reference points obtained by the algorithm. This index can simultaneously evaluate the convergence and distribution of the solution set in space. Convergence refers to the approximation degree between the solution set and the real Pareto front, distribution refers to the wide distribution of the whole solution set in the target space. A larger HV value indicates that the algorithm can dominate a larger area in the target space and has better diversity. The formula of HV index is as follows:(10)HV=δUi=1Svi,where δ represents the Lebesgue measure, which is used to measure volume. S is the number of nondominated solution set, vi represents the supervolume consisting of the reference point and the ith solution in the nondominant solution set.


## 6. Comparison of Experimental Results

MSMOPSO amplifies the archiving strategy, selects different archiving strategies according to different adaptability conditions; then, in the second archive, self-organizing dynamic clustering strategy and species migration model are applied to update it; finally, Gaussian perturbation strategy is applied to change the positions of particles in the population, so the optimal solution set is obtained. In order to further detect the superiority of the performance of this algorithm, it is compared with other evolutionary algorithms, and experimental results are described below.


### 6.1. Comparison of IGD Data

Table4 shows the results of IGD evaluation indexes obtained by this algorithm and other MOPSOs. The standard deviation part is marked by parentheses. It can be obtained from Table 4: among the 22 test functions, 12 of the results of MSMOPSO are better than those of other algorithms, SMPSO has 2 better results, MPSOD has 2 better results, MOPSOCD has one better result, NMPSO has 2 better results, dMOPSO has 0 better result, and CMOPSO has 3 better results. It shows that MSMOPSO has better convergence performance and particle distribution than MOPSOs. In addition, considering the standard deviation of IGD index, the performance of MSMOPSO is more stable than other MOPSOs in the test functions DTLZ1, UF1, UF2, UF6-8, UF10, ZDT2, and ZDT4.Table 4
Compared with IGD index of MOPSOs.
FunctionsSMPSOMPSODMOPSOCDNMPSOCMOPSOdMOPSOMSMOPSODTLZ11.67e + 01 (1.92e + 01)3.87e + 01 (6.23e + 00)3.93e + 01 (1.13e + 01)3.17e + 01 (5.86e + 00)6.39e + 01 (1.61e + 01)9.90e + 00 (9.18e + 00)3.27e − 02 (2.08e − 03)DTLZ29.29e − 02 (1.02e − 02)5.67e − 02 (6.22e − 04)1.00e − 01 (9.27e − 03)7.69e − 02 (3.18e − 03)5.82e − 02 (1.17e − 03)1.45e − 01 (9.05e − 03)2.11e − 01 (2.27e − 02)DTLZ33.11e + 01 (2.98e + 01)1.15e + 02 (1.41e + 01)1.10e + 02 (2.84e + 01)9.31e + 01 (2.72e + 01)1.46e + 02 (3.56e + 01)7.74e + 01 (6.43e + 01)1.40e + 02 (2.45e + 01)DTLZ44.17e − 01 (2.98e + 01)6.26e − 02 (6.85e − 03)3.00e − 01 (5.48e − 02)1.55e − 01 (1.76e − 01)1.22e − 01 (2.24e − 01)3.12e − 01 (2.20e − 02)3.30e − 01 (4.23e − 02)DTLZ55.71e − 03 (3.98e − 04)5.35e − 02 (5.34e − 03)2.28e − 02 (6.99e − 03)1.13e − 02 (1.37e − 03)7.32e − 03 (6.70e − 04)4.42e − 02 (6.78e − 03)6.38e − 02 (1.78e − 02)DTLZ61.17e + 00 (9.91e − 01)2.15e − 01 (2.38e − 01)7.85e − 02 (2.24e − 01)1.45e − 02 (3.06e − 03)6.99e − 02 (2.50e − 01)3.28e − 02 (2.39e − 04)1.15e − 02 (2.62e − 03))DTLZ71.86e − 01 (1.52e − 01)1.65e − 01 (2.10e − 02)8.72e − 02 (6.07e − 03)7.04e − 02 (4.99e − 03)1.68e − 01 (2.06e − 01)1.66e − 01 (1.33e − 01)1.47e − 01 (2.15e − 02)UF13.61e − 01 (7.22e − 02)2.67e − 01 (4.33e − 02)6.52e − 01 (1.56e − 01)1.39e − 01 (8.31e − 02)9.59e − 02 (1.36e − 02)6.52e − 01 (1.01e − 01)8.79e − 02 (1.45e − 02)UF29.97e − 02 (8.61e − 03)1.17e − 01 (1.11e − 02)1.40e − 01 (1.38e − 02)8.44e − 02 (6.99e − 03)6.14e − 02 (6.58e − 03)9.58e − 02 (5.86e − 03)4.04e − 02 (4.69e − 03)UF34.56e − 01 (5.18e − 02)5.00e − 01 (1.43e − 02)3.82e − 01 (7.66e − 2)3.66e − 1 (4.91e − 2)3.96e − 01 (3.27e − 02)3.33e − 01 (5.62e − 03)2.07e − 01 (2.31e − 02)UF41.10e − 01 (3.28e − 03)9.91e − 02 (3.92e − 03)8.06e − 02 (6.81e − 03)6.55e − 02 (9.54e − 03)1.11e − 01 (7.75e − 03)1.35e − 01 (5.92e − 03)7.35e − 02 (1.49e − 02)UF52.92e + 00 (4.91e − 01)2.72e + 00 (2.77e − 01)3.85e + 00 (4.64e − 01)1.64e + 00 (4.96e − 01)8.93e − 01 (1.95e − 01)3.19e + 00 (2.93e − 01)8.66e − 01 (3.20e − 01)UF61.27e + 00 (4.19e − 01)1.33e + 00 (2.01e − 01)3.13e + 00 (8.08e − 01)6.90e − 01 (1.13e − 01)3.98e − 01 (5.98e − 02)2.30e + 00 (5.43e − 01)4.23e − 01 (7.62e − 02)UF73.29e − 01 (1.20e − 01)2.38e − 01 (5.67e − 02)5.91e − 01 (1.19e − 01)2.13e − 01 (1.62e − 01)1.44e − 01 (1.25e − 01)3.59e − 01 (8.49e − 02)4.38e − 02 (4.96e − 03)UF83.90e − 01 (4.54e − 02)5.54e − 01 (5.12e − 02)7.66e − 01 (1.64e − 01)4.51e − 01 (5.65e − 02)5.93e − 01 (9.13e − 02)3.73e − 01 (3.90e − 02)2.88e − 01 (1.06e − 02)UF95.61e − 01 (3.92e − 02)6.68e − 01 (3.63e − 02)8.75e − 01 (1.39e − 01)4.63e − 01 (6.11e − 02)8.90e − 01 (1.33e − 01)5.72e − 01 (3.80e − 02)2.94e − 01 (7.12e − 02)UF102.71e + 00 (5.13e − 01)4.12e + 00 (3.75e − 01)5.17e + 00 (9.18e − 01)1.65e + 00 (3.36e − 01)4.50e + 00 (5.24e − 01)9.40e − 01 (2.83e − 02)3.31e − 01 (5.78e − 02)ZDT15.39e − 02 (6.45e − 02)9.43e − 02 (3.56e − 02)4.89e − 03 (1.04e − 02)3.03e − 02 (1.30e − 02)3.20e − 03 (5.76e − 04)5.26e − 02 (1.53e − 02)5.52e − 03 (7.89e − 04)ZDT21.30e − 01 (2.00e − 01)2.00e − 01 (1.46e − 01)1.68e − 01 (2.40e − 01)4.84e − 02 (7.60e − 02)2.81e − 03 (3.91e − 04)3.98e − 02 (1.07e − 02)6.22e − 03 (1.05e − 03)ZDT31.84e − 01 (9.98e − 02)2.11e − 01 (6.69e − 02)6.59e − 02 (6.10e − 02)8.81e − 02 (2.57e − 02)3.64e − 03 (5.89e − 04)3.50e − 02 (7.06e − 03)2.03e − 01 (5.17e − 03)ZDT49.02e + 01 (3.23e + 01)2.28e + 02 (2.12e + 01)1.57e + 02 (3.25e + 01)1.28e + 02 (3.21e + 01)1.72e + 02 (3.04e + 01)3.07e + 01 (4.78e + 01)5.57e − 03 (8.40e − 04)ZDT61.50e + 00 (1.82e + 00)2.04e + 00 (1.45e + 00)1.62e + 00 (1.91e + 00)2.35e − 03 (2.30e − 04)2.75e − 01 (2.71e − 01)3.18e − 01 (1.02e − 01)2.17e − 03 (2.59e − 04)The bold values in Table 4 represent: under the premise of taking the mean as the main comparison object, the optimal values in the comparison results are compared for different detection functions.Table5 shows the results of IGD evaluation indexes obtained by MSMOPSO and other MOEAs. The standard deviation part is marked by parentheses. As shown in Table 5, 14 of the 22 test functions of MSMOPSO have better results compared with the solving mean of other algorithms, ARMOEA has 3 better results, MOEAD has 1 better result, NSGAIISDR has 2 better results, gNSGAII has 0 better result, RVEA has 1 better result, and GrEA has 1 better result. It shows that MSMOPSO has better convergence performance and particle distribution than MOEAs. In addition, in the test functions DTLZ1, UF1, UF2, UF6-8, UF10, ZDT1, ZDT2, ZDT4, and ZDT6, the standard deviation of IGD index of MSMOPSO is smaller than that of other algorithms, indicating that the convergence performance and particle distribution performance of MSMOPSO are more stable than other MOEAs.Table 5
Compared with IGD index of MOEAs.
FunctionsARMOEAMOEADNSGAIISDRgNSGAIIRVEAGrEAMSMOPSODTLZ13.16e + 00 (1.80e + 00)4.81e + 00 (2.66e + 00)2.60e + 00 (1.40e + 00)2.14e + 01 (5.60e + 00)5.43e + 00 (2.42e + 00)2.84e + 00 (1.78e + 00)3.27e − 02 (2.08e − 03)DTLZ25.50e − 02 (1.70e − 04)5.49e − 02 (2.37e − 04)4.66e − 01 (2.61e − 02)3.86e − 01 (5.48e − 03)5.58e − 02 (7.40e − 04)6.61e − 02 (1.82e − 03)2.11e − 01 (2.27e − 02)DTLZ39.43e + 00 (3.92e + 00)1.54e + 01 (9.61e + 00)8.08e + 00 (3.39e + 00)5.02e + 01 (1.41e + 01)1.66e + 01 (5.95e + 00)8.70e + 00 (4.40e + 00)1.40e + 02 (2.45e + 01)DTLZ43.84e − 01 (3.21e − 01)5.86e − 01 (3.13e − 01)5.21e − 01 (9.07e − 02)3.95e − 01 (1.07e − 02)1.05e − 01 (1.48e − 01)1.57e − 01 (2.46e − 01)3.30e − 01 (4.23e − 02)DTLZ55.73e − 03 (1.39e − 04)3.25e − 02 (6.41e − 04)3.33e − 02 (5.91e − 03)3.11e − 01 (9.11e − 03)7.77e − 02 (9.51e − 03)2.07e − 02 (1.03e − 03)6.38e − 02 (1.78e − 02)DTLZ68.75e − 03 (1.98e − 02)1.99e − 01 (4.08e − 01)7.08e − 02 (6.31e − 02)7.62e − 01 (3.37e − 01)1.22e − 01 (1.07e − 01)2.21e − 02 (5.76e − 04)1.15e − 02 (2.62e − 03)DTLZ72.32e − 01 (2.55e − 01)1.74e − 01 (1.19e − 01)8.71e − 02 (5.40e − 03)4.02e − 01 (3.05e − 03)1.23e − 01 (1.11e − 02)1.03e − 01 (7.95e − 02)1.47e − 01 (2.15e − 02)UF11.22e − 01 (2.88e − 02)3.21e − 01 (7.76e − 02)1.42e − 01 (3.06e − 02)3.35e − 01 (1.99e − 02)1.75e − 01 (4.00e − 02)1.23e − 01 (3.12e − 02)8.79e − 02 (1.45e − 02)UF27.60e − 02 (1.43e − 02)2.09e − 01 (7.56e − 02)8.44e − 02 (5.55e − 03)2.95e − 01 (2.30e − 02)1.20e − 01 (8.96e − 03)6.91e − 02 (5.77e − 03)4.04e − 02 (4.69e − 03)UF34.33e − 01 (2.85e − 02)3.36e − 01 (3.32e − 02)4.81e − 01 (1.55e − 02)4.46e − 01 (1.63e − 02)4.94e − 01 (1.61e − 02)4.62e − 01 (1.62e − 02)2.07e − 01 (2.31e − 02)UF48.29e − 02 (3.18e − 03)1.16e − 01 (5.37e − 03)9.93e − 02 (4.02e − 03)2.20e − 01 (1.43e − 03)1.523e − 01 (4.80e − 03)8.55e − 02 (3.18e − 03)7.35e − 02 (1.49e − 02)UF58.19e − 01 (2.17e − 01)1.47e + 00 (3.21e − 01)1.41e + 00 (3.28e − 01)8.71e − 01 (2.66e − 01)1.16e + 00 (3.73e − 01)9.04e − 01 (2.68e − 01)8.66e − 01 (3.20e − 01)UF65.32e − 01 (1.00e − 01)5.52e − 01 (1.53e − 01)6.85e − 01 (1.28e − 01)5.46e − 01 (1.13e − 01)6.98e − 01 (1.02e − 01)6.55e − 01 (1.52e − 01)4.23e − 01 (7.62e − 02)UF72.29e − 01 (1.32e − 01)4.18e − 01 (1.35e − 01)1.73e − 01 (7.94e − 02)3.12e − 01 (4.69e − 02)2.45e − 01 (6.35e − 02)1.42e − 01 (8.44e − 02)4.38e − 02 (4.96e − 03)UF83.67e − 01 (5.34e − 02)5.87e − 01 (2.74e − 01)3.88e − 01 (3.59e − 02)5.11e − 01 (1.98e − 02)4.16e − 01 (3.58e − 02)3.14e − 01 (3.13e − 02)2.88e − 01 (1.06e − 02)UF94.94e − 01 (7.59e − 02)5.25e − 01 (8.58e − 02)5.33e − 01 (1.02e − 01)5.87e − 01 (4.80e − 02)4.78e − 01 (6.58e − 02)4.66e − 01 (8.41e − 02)2.94e − 01 (7.12e − 02)UF101.13e + 00 (2.93e − 01)7.22e − 01 (7.97e − 02)2.37e + 00 (5.73e − 01)1.25e + 00 (3.83e − 01)1.72e + 00 (5.17e − 01)2.06e + 00 (5.34e − 01)3.31e − 01 (5.78e − 02)ZDT16.41e − 02 (4.55e − 02)1.85e − 01 (7.29e − 02)1.39e − 01 (2.37e − 02)3.28e − 01 (2.38e − 02)2.01e − 01 (2.99e − 02)5.34e − 02 (1.05e − 02)5.52e − 03 (7.89e − 04)ZDT26.97e − 01 (3.20e − 02)5.73e − 01 (5.82e − 02)2.05e − 01 (4.00e − 02)2.40e − 01 (2.59e − 02)2.87e − 01 (4.18e − 02)7.94e − 02 (1.45e − 02)6.22e − 03 (1.05e − 03)ZDT35.40e − 02 (1.75e − 02)1.58e − 01 (6.72e − 02)1.06e − 01 (1.45e − 02)3.38e − 01 (3.42e − 02)1.89e − 01 (1.42e − 02)4.84e − 02 (1.42e − 02)2.03e − 01 (5.17e − 03)ZDT44.79e + 01 (6.75e + 00)1.73e + 01 (5.14e + 00)4.28e + 01 (6.28e + 00)4.79e + 01 (7.66e + 00)5.57e + 01 (1.08e + 1)3.13e + 01 (5.09e + 00)5.57e − 03 (8.40e − 04)ZDT63.46e + 00 (1.57e − 01)2.21e + 00 (6.45e − 01)3.60e + 00 (2.04e − 01)3.32e + 00 (2.12e − 01)3.60e + 00 (1.74e − 01)2.99e + 00 (1.89e − 01)2.17e − 03 (2.59e − 04)The bold values in Table 5 represent under the premise of taking the mean as the main comparison object, the optimal values in the comparison results are compared for different detection functions.Based on the comparison of IGD indexes in Table4 and Table 5, it is found that MSMOPSO performs relatively well both for MOPSOs and MOEAs, especially when detection functions are DTLZ1, ZDT4, ZDT6, and UF series (UFs). It is shown that the Gaussian mutation strategy increases the possibility of finding better solutions, and more optimal solutions are stored by the multi-archiving mechanism, which can further improve the exploration ability of the algorithm. This can enable MSMOPSO to find as many noninferior solutions as possible. In addition, according to the statistical results, the variance of IGD index value of the proposed algorithm for MOPs of DTLZ1, UF2, UF7, UF8, ZDT4, and ZDT6 is smaller, indicating that the algorithm has certain stability in solving specific optimization problems, and has some advantages in convergence to the optimal solution set.


### 6.2. Comparison of HV Data

Table6 shows the solution results of HV evaluation index of MSMOPSO and other MOPSOs. The standard deviation part is marked by parentheses. As shown in Table 6, among the 22 test functions, 10 of the results of MSMOPSO are better than those of other algorithms. SMPSO has 2 better results, MPSOD has 1 better result, MOPSOCD has 0 better result, NMPSO has 5 better results, CMOPSO has 4 better results, and dMOPSO has 0 better result. It indicates that the results of MSMOPSO are more widely distributed in objective space. In addition, in the test functions UF1, UF2, UF7, ZDT2, and ZDT3 of MSMOPSO, the standard deviation of HV index is smaller than that of other algorithms; it illustrates that MSMOPSO has better robustness compared with other algorithms.Table 6
Compared with HV index of MOPSOs.
FunctionsSMPSOMPSODMOPSOCDNMPSOCMOPSOdMOPSOMSMOPSODTLZ11.06e − 01 (2.06e − 01)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)1.20e − 02 (5.58e − 02)8.19e − 01 (4.59e − 03)DTLZ24.55e − 01 (2.06e − 02)5.47e − 01 (2.14e − 03)4.60e − 01 (1.96e − 02)5.59e − 01 (1.08e − 03)5.40e − 01 (2.89e − 03)3.67e − 01 (1.45e − 02)3.07e − 01 (2.74e − 02)DTLZ32.82e − 02 (4.75e − 02)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)DTLZ43.29e − 01 (1.06e − 01)5.38e − 01 (1.11e − 02)3.25e − 01 (6.09e − 02)5.23e − 01 (8.30e − 02)5.04e − 01 (1.12e − 01)3.09e − 01 (4.92e − 02)3.36e − 01 (3.46e − 02)DTLZ51.98e − 01 (4.12e − 04)1.46e − 01 (7.08e − 03)1.84e − 01 (5.62e − 03)1.97e − 01 (5.42e − 04)1.97e − 01 (5.06e − 04)1.50e − 01 (1.05e − 02)1.18e − 01 (1.78e − 02)DTLZ65.68e − 02 (8.97e − 02)1.19e − 01 (5.63e − 02)1.85e − 01 (4.43e − 02)1.96e − 01 (1.07e − 03)1.87e − 01 (5.08e − 02)1.83e − 01 (2.07e − 04)1.90e − 01 (6.17e − 03)DTLZ72.48e − 01 (1.77e − 02)2.16e − 01 (1.57e − 02)2.67e − 01 (3.25e − 03)2.75e − 01 (1.46e − 03)2.63e − 01 (2.13e − 02)2.42e − 01 (1.07e − 02)2.48e − 01 (6.98e − 03)UF12.79e − 01 (6.38e − 02)3.60e − 01 (4.75e − 02)8.04e − 02 (8.04e − 02)5.25e − 01 (4.72e − 02)5.73e − 01 (1.94e − 02)6.90e − 02 (6.61e − 02)5.97e − 01 (2.04e − 02)UF26.10e − 01 (7.96e − 03)5.75e − 01 (1.26e − 02)5.45e − 01 (1.72e − 02)6.17e − 01 (8.01e − 03)6.46e − 01 (7.19e − 03)6.15e − 01 (5.56e − 03)6.73e − 01 (6.53e − 03)UF32.07e − 01 (4.61e − 02)1.74e − 01 (1.31e − 02)2.50e − 01 (5.55e − 2)2.78e − 1 (4.07e − 2)2.62e − 01 (2.03e − 02)3.06e − 01 (7.11e − 03)4.69e − 01 (2.76e − 02)UF42.94e − 01 (8.10e − 03)3.08e − 01 (5.05e − 03)3.31e − 01 (8.75e − 03)3.55e − 01 (1.34e − 02)2.90e − 01 (9.55e − 03)2.55e − 01 (7.03e − 03)3.48e − 01 (8.30e − 03)UF50.00e + 00 (0.00e + 00)0.00e + 00 (0.00e − 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)8.47e − 03 (2.26e − 02)0.00e + 00 (0.00e + 00)3.04e − 03 (1.07e − 02)UF62.03e − 04 (1.02e − 03)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)1.51e − 02 (1.94e − 02)1.35e − 01 (7.59e − 02)0.00e + 00 (0.00e + 00)6.67e − 02 (1.95e − 03)UF72.18e − 01 (9.91e − 02)2.68e − 01 (6.41e − 02)3.06e − 02 (3.55e − 02)3.59e − 01 (1.12e − 01)4.31e − 01 (8.88e − 02)1.92e − 01 (6.16e − 02)5.22e − 01 (7.28e − 03)UF81.65e − 01 (4.59e − 02)5.64e − 02 (1.98e − 02)8.24e − 03 (1.28e − 02)2.86e − 01 (4.64e − 02)1.76e − 02 (1.65e − 02)2.61e − 01 (2.38e − 02)2.94e − 01 (1.52e − 02)UF92.05e − 01 (3.82e − 02)1.06e − 01 (2.27e − 02)3.12e − 02 (3.41e − 02)3.13e − 01 (5.37e − 02)2.34e − 02 (2.83e − 02)2.29e − 01 (1.55e − 02)4.48e − 01 (7.20e − 02)UF100.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)9.10e − 02 (8.14e − 04)2.92e − 01 (2.24e − 02)ZDT16.54e − 01 (8.02e − 02)5.84e − 01 (4.85e − 02)7.18e − 01 (1.51e − 02)6.88e − 01 (1.72e − 02)7.20e − 01 (8.85e − 04)6.58e − 01 (1.70e − 02)7.17e − 01 (1.83e − 03)ZDT23.41e − 01 (1.49e − 01)2.34e − 01 (1.14e − 01)3.20e − 01 (1.59e − 01)4.08e − 01 (7.01e − 02)4.45e − 01 (7.18e − 04)3.85e − 01 (1.69e − 02)4.40e − 01 (2.47e − 03)ZDT35.31e − 01 (7.24e − 02)4.57e − 01 (6.08e − 02)5.59e − 01 (4.33e − 02)5.72e − 01 (9.34e − 03)6.00e − 01 (1.23e − 03)6.01e − 01 (1.12e − 02)6.53e − 01 (2.75e − 03)ZDT40.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)3.44e − 02 (4.87e − 02)7.17e − 01 (1.83e − 03)ZDT61.97e − 01 (1.75e − 01)4.72e − 02 (8.39e − 02)1.69e − 01 (1.62e − 01)3.90e − 01 (1.99e − 04)1.87e − 01 (1.56e − 01)1.23e − 01 (6.93e − 02)3.83e − 01 (1.62e − 03)Table7 shows the solution results of the HV evaluation index of MSMOPSO and other MOEAs, and the standard deviation part is marked by parentheses. As shown in Table 7, 13 of the 22 test functions of this algorithm have better results compared with the solving mean of other algorithms, ARMOEA has 4 better results, MOEAD has 1 better result, NSGAIISDR has 1 better result, GrEA has 1 better result, and RVEA has 1 better result. It shows that the spatial distribution of MSMOPSO is better than other MOEAs. In addition, in the test functions UF1, UF2, UF7, ZDT1, ZDT2, and ZDT3 of this algorithm, the standard deviation of HV index is smaller than that of other algorithms, indicating that this algorithm has good stability compared with other algorithms.Table 7
Compared with HV index of MOEAs.
FunctionsARMOEAMOEADNSGAIISDRgNSGAIIRVEAGrEAMSMOPSODTLZ10.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)8.19e − 01 (4.59e − 03)DTLZ25.56e − 01 (5.76e − 04)5.55e − 01 (9.01e − 04)2.55e − 01 (1.45e − 02)1.59e − 01 (3.54e − 04)5.53e − 01 (1.33e − 03)5.58e − 01 (9.00e − 04)3.07e − 01 (2.74e − 02)DTLZ30.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)DTLZ44.03e − 01 (1.59e − 01)2.92e − 01 (1.63e − 01)2.85e − 01 (5.48e − 02)1.60e − 01 (2.75e − 04)5.31e − 01 (6.62e − 02)5.11e − 01 (1.28e − 01)3.36e − 01 (3.46e − 02)DTLZ51.98e − 01 (2.20e − 04)1.82e − 01 (9.22e − 04)1.80e − 01 (3.08e − 03)6.02e − 02 (3.55e − 05)1.50e − 01 (7.73e − 03)1.88e − 01 (4.74e − 04)1.18e − 01 (1.78e − 02)DTLZ61.96e − 01 (1.63e − 02)1.46e − 01 (6.73e − 02)1.63e − 01 (2.19e − 02)9.79e − 03 (1.74e − 02)1.21e − 01 (3.90e − 02)1.88e − 01 (2.63e − 04)1.90e − 01 (6.17e − 03)DTLZ72.55e − 01 (2.51e − 02)2.45e − 01 (9.61e − 03)2.69e − 01 (1.76e − 03)7.56e − 02 (6.89e − 04)2.48e − 01 (7.73e − 03)2.70e − 01 (1.06e − 02)2.48e − 01 (6.98e − 03)UF15.57e − 01 (2.76e − 02)4.14e − 01 (5.21e − 02)5.09e − 01 (4.17e − 02)3.36e − 01 (5.33e − 02)4.65e − 01 (5.23e − 02)5.31e − 01 (4.03e − 02)5.97e − 01 (2.04e − 02)UF26.32e − 01 (9.92e − 03)5.61e − 01 (3.29e − 02)6.09e − 01 (8.31e − 03)4.30e − 01 (4.61e − 02)5.60e − 01 (1.12e − 02)6.32e − 01 (7.17e − 03)6.73e − 01 (6.53e − 03)UF32.05e − 01 (2.27e − 02)3.06e − 01 (4.61e − 02)1.91e − 01 (1.41e − 02)2.10e − 1 (1.17e − 02)1.78e − 01 (1.12e − 02)2.03e − 01 (1.29e − 02)4.69e − 01 (2.76e − 02)UF43.31e − 01 (4.76e − 03)2.84e − 01 (6.97e − 03)3.08e − 01 (4.66e − 03)2.42e − 01 (1.51e − 03)2.30e − 01 (5.14e − 03)3.27e − 01 (4.17e − 03)3.48e − 01 (8.30e − 03)UF51.42e − 02 (2.84e − 02)0.00e + 00 (0.00e − 00)0.00e + 00 (0.00e + 00)7.42e − 03 (1.82e − 02)5.24e − 03 (1.62e − 02)6.06e − 03 (2.03e − 02)3.04e − 03 (1.07e − 02)UF65.25e − 02 (3.75e − 02)7.61e − 02 (6.40e − 02)1.28e − 02 (1.86e − 02)4.46e − 02 (3.49e − 02)2.03e − 02 (2.10e − 02)1.99e − 02 (2.55e − 02)6.67e − 02 (1.95e − 03)UF73.55e − 01 (8.09e − 02)2.54e − 01 (7.49e − 02)3.60e − 01 (7.25e − 02)2.52e − 01 (3.90e − 02)2.82e − 01 (5.49e − 02)3.97e − 01 (8.12e − 02)5.22e − 01 (7.28e − 03)UF82.96e − 01 (2.05e − 02)1.54e − 01 (6.00e − 02)1.03e − 01 (4.94e − 02)7.58e − 02 (1.93e − 02)7.63e − 02 (4.86e − 02)2.27e − 01 (5.17e − 02)2.94e − 01 (1.52e − 02)UF92.86e − 01 (5.22e − 02)2.88e − 01 (4.81e − 02)2.03e − 01 (8.42e − 02)1.38e − 01 (2.10e − 02)2.61e − 01 (5.10e − 02)2.74e − 01 (7.82e − 02)4.48e − 01 (7.20e − 02)UF100.00e + 00 (0.00e + 00)3.68e − 02 (2.43e − 02)0.00e + 00 (0.00e + 00)1.01e − 07 (5.51e − 07)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)2.92e − 01 (2.24e − 02)ZDT16.47e − 01 (3.05e − 02)5.39e − 01 (6.01e − 02)5.44e − 01 (2.66e − 02)3.51e − 01 (2.92e − 02)4.61e − 01 (3.19e − 02)6.51e − 01 (1.37e − 02)7.17e − 01 (1.83e − 03)ZDT21.09e − 03 (3.71e − 03)8.82e − 02 (1.55e − 02)2.02e − 01 (3.35e − 02)2.20e − 01 (2.45e − 02)1.31e − 01 (2.89e − 02)3.35e − 01 (1.95e − 02)4.40e − 01 (2.47e − 03)ZDT35.88e − 01 (4.01e − 02)5.53e − 01 (5.15e − 02)5.34e − 01 (2.23e − 02)5.84e − 01 (5.63e − 02)4.51e − 01 (1.53e − 02)5.67e − 01 (8.99e − 03)6.53e − 01 (2.75e − 03)ZDT40.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)7.17e − 01 (1.83e − 03)ZDT60.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)3.83e − 01 (1.62e − 03)The bold values in Table 7 represent under the premise of taking the mean as the main comparison object, the optimal values in the comparison results are compared for different detection functions.Based on Table6 and Table 7, it is found that compared with MOPSOs and MOEAs, the HV index value of this algorithm is relatively large in DTLZ1, UF1-3, UF7-10, and ZDT3-4 detection functions. According to the characteristics of HV index, the larger the HV value of this algorithm, the better the spatial distribution of decision variables in it. On this basis, statistical results show that the variance of HV index value of MSMOPSO is smaller for DTLZ1, UF7, and ZDT4 test functions, indicating that it has better stability and significant solving ability when solving these two kinds of MOPs. In conclusion, the application of Gaussian mutation strategy can improve the spatial distribution of particles in this proposed algorithm, and then the multi-archive is used to store the noninferior solutions, further indicating that these strategies can effectively improve the development ability of this algorithm for the target space and effectively alleviate the disadvantage of the algorithm easily falling into local convergence.


### 6.3. Comparison of Box Plots

In this part, some test functions are taken as examples to store the IGD index of each running result of MSMOPSO, which is displayed in the form of box graph and compared with the calculation results of other comparison algorithms. It is found that the distance between the upper and the lower quartile loci of MSMOPSO is closer than that of other algorithms, indicating that the distribution of experimental data is more concentrated, so MSMOPSO is more stable. Meanwhile, for outliers in the data, “+” is used to represent them in the box charts. It can be seen from these plots that there are fewer outliers in the experimental data of MSMOPSO, indicating that the solution results are better. To sum up, the performance of MSMOPSO is more significant than other algorithms, both in terms of solution results and algorithm stability. The box plots are shown in Figure6, where 1 stands for MPSOD, 2 stands for MOPSOCD, 3 stands for SMPSO, 4 stands for NMPSO, 5 stands for ARMOEA, 6 stands for MOEAD, 7 stands for NSGAIISDR, 8 stands for gNSGAII, and 9 stands for MSMOPSO.Figure 6
(a) Boxplot for the statistical results for ZDT1. (b) Boxplot for the statistical results for ZDT2. (c) Boxplot for the statistical results for ZDT4. (d) Boxplot for the statistical results for ZDT6. (e) Boxplot for the statistical results for UF1. (f) Boxplot for the statistical results for UF2. (g) Boxplot for the statistical results for UF3. (h) Boxplot for the statistical results for UF6. (i) Boxplot for the statistical results for UF7. (j) Boxplot for the statistical results for UF8. (k) Boxplot for the statistical results for UF10. (l) Boxplot for the statistical results for DTLZ1.


### 6.4. Comparison of Results

In order to verify the authenticity and validity of the experimental results, the optimal solution set of each algorithm in the objective space is displayed on the front surface of each algorithm, represented by a two-dimensional detection function ZDT4 and a three-dimensional detection function DTLZ1, as shown in Figure7 and Figure 8. The points in the figure represent the optimal solution set of the algorithm on test functions, and the Pareto front is represented by shadow.Figure 7
(a) Results of MPSOD on DTLZ1. (b) Results of MOPSOCD on DTLZ1. (c) Results of SMPSO on DTLZ1. (d) Results of NMPSO on DTLZ1. (e) Results of ARMOEA on DTLZ1. (f) Results of MOEAD on DTLZ1. (g) Results of NSGAIISDR on DTLZ1. (h) Results of gNSGAII on DTLZ1. (i) Results of MSMOPSO on DTLZ1.Figure 8
(a) Results of MPSOD on ZDT4. (b) Results of MOPSOCD on ZDT4. (c) Results of SMPSO on ZDT4. (d) Results of NMPSO on ZDT4. (e) Results of ARMOEA on ZDT4. (f) Results of MOEAD on ZDT4. (g) Results of NSGAIISDR on ZDT4. (h) Results of gNSGAII on ZDT4. (i) Results of MSMOPSO on ZDT4.Comparing Figures7 and 8, it is found that the algorithm is evenly distributed on the ZDT4 and the optimal solution set is completely distributed on the Pareto front, while other test algorithms are not attached to the front surface. In addition, due to the large numerical value of solution results obtained by other test functions, there is a large difference between them and the optimal solution set that forms the front surface, so it is not shown in other comparison diagrams. For the DTLZ1, only MSMOPSO is evenly distributed on the Pareto front surface, while the optimal solution sets of ARMOEA, MOEAD, and NSGAIISDR are too scattered, so they do not converge to the front surface. The solution values of MPSOD, MOPSOCD, SMPSO, NMPSO, and gNSGAII are large and the distribution of the solution set is too scattered, so the Pareto front does not appear.


### 6.5. Complexity of MSMOPSO

In MSMOPSO, the first archive is updated by elite selection strategy, self-organizing clustering strategy and species migration strategy are applied to expand and update the second archive, and finally, the particles of the whole swarm are randomly disturbed by Gaussian mutation strategy. Therefore, the time complexity of this algorithm mainly depends on the above three aspects. For the first archive, the elite strategy is applied to select the particles in the archive and the swarm, so the time complexity of this process isOP∗N∗M (P represents the archive size, N represents the population size, and M represents the dimension of decision variables). For the second archive, self-organizing clustering strategy is used to select and store the particles within swarm according to the grouping center, then physical migration strategy is utilized in updating this archive, so the time complexity of the second archive is divided into two parts, one is a self-organizing clustering strategy, the other is the species migration strategy. The time complexity of the first part is ON∗N∗M, and the second part is ON∗N∗M+P∗N∗M. Because N≥P, the complexity of the second archive is ON∗N∗M. For Gaussian mutation strategy, the computational complexity is ON∗D∗M∗T (T represents the number of iterations of the algorithm). So MSMOPSO’s time complexity is OP∗N∗M+ON∗N∗M+N∗N∗M+ON∗D∗M∗T, that is ON2M+OPNM+ONDMT.In addition, the (tic, toc) order was used to record the running time of MSMOPSO in 1000 evaluations, and the results were compared with MOEAs and MOPSOs. The comparison results are shown in Table8 and Table 9. As can be seen from the table, the calculation time of MSMOPSO for different detection functions is larger than that of other algorithms, but the difference is relatively small. It is worth noting that these results of running time depend on the influence of computer hardware, software, and other factors, so the running result is only for reference.Table 8
Compared with runtime of MOEAs.
FunctionsARMOEAMOEADNSGAIISDRgNSGAIIRVEAGrEAMSMOPSODTLZ17.55e−011.53e + 002.52e−011.51e−012.35e−013.33e−017.13e + 00DTLZ23.81e + 001.50e + 002.72e−011.42e−012.31e−018.10e−017.92e + 00DTLZ34.18e−011.53e + 002.51e−011.51e−012.35e−012.89e−016.44e + 00DTLZ43.07e + 001.55e + 003.00e−011.58e−012.35e−017.60e−011.78e + 00DTLZ53.07e + 001.53e + 002.42e−011.45e−011.98e−016.78e−018.00e + 00DTLZ62.71e + 001.56e + 002.50e−011.55e−012.28e−016.24e−014.78e + 01DTLZ72.97e + 001.50e + 002.43e−011.48e−012.07e−017.36e−015.78e−01UF13.37e−011.74e + 003.04e−011.52e−012.26e−012.73e−016.03e−01UF28.28e−011.86e + 003.12e−011.60e−012.34e−013.13e−011.81e + 00UF33.67e−011.82e + 003.19e−011.567e−012.40e−012.87e−011.34e + 00UF49.05e−011.77e + 003.04e−011.41e−012.34e−013.76e−013.51e + 00UF52.86e−011.75e + 003.10e−011.52e−012.23e−012.60e−011.28e + 00UF62.93e−011.76e + 003.07e−011.52e−012.27e−012.62e−011.37e + 00UF73.79e−011.75e + 003.04e−011.48e−012.27e−012.77e−011.29e + 00UF81.67e + 001.93e + 003.05e−011.54e−012.12e−014.55e−018.47e + 00UF91.17e + 001.93e + 003.08e−011.55e−012.15e−014.16e−018.41e + 00UF108.41e−011.81e + 003.08e−011.56e−012.12e−014.15e−012.06e + 00ZDT15.750e−011.68e + 003.01e−011.53e−012.22e−012.95e−013.90e + 00ZDT23.07e−011.67e + 002.99e−011.44e−012.16e−012.63e−019.56e + 00ZDT37.90e−111.70e + 002.99e−011.49e−012.28e−013.17e−018.05e + 00ZDT42.68e−111.66e + 003.03e−011.54e−012.14e−012.54e−019.29e + 00ZDT62.62e−11.68e + 003.02e−011.51e−112.16e−012.50e−018.53e + 00Table 9
Compared with runtime of MOPSOs.
FunctionsCMOPSOMOPSOCDMPSODNMPSOSMPSOdMOPSOMSMOPSODTLZ12.61e−011.31e−011.30e + 001.91e + 001.42e−013.84e−017.13e + 00DTLZ21.30e + 001.57e−011.30e + 002.90e + 001.61e−013.45e−017.92e + 00DTLZ32.53e−011.27e−011.29e + 009.83e−011.38e−013.78e−016.44e + 00DTLZ49.54e−011.22e−011.07e + 002.96e + 001.49e−013.63e−011.78e + 00DTLZ57.57e−011.28e−018.14e−012.30e + 001.60e−013.61e−018.00e + 00DTLZ67.24e−012.20e−011.18e + 003.68e + 001.63e−013.73e−014.78e + 01DTLZ77.68e−012.31e−016.14e−012.98e + 001.53e−013.24e−015.78e−01UF12.36e−011.28e−011.16e + 005.34e−011.41e−013.39e−016.03e−01UF22.74e−011.34e−011.20e + 008.82e−011.59e−013.70 e−011.81e + 00UF32.62e−012.12e−011.16e + 001.07e + 001.61e−013.48e−011.34e + 00UF42.72e−011.40e−011.27e + 002.77e + 001.53e−13.67e−013.51e + 00UF52.28e−011.32e−011.07e + 004.97e−011.45e−013.61e−011.28e + 00UF62.30e−011.34e−011.10e + 005.11e−011.48e−013.61e−011.37e + 00UF72.52e−011.31e−011.14e + 006.64e−011.45e−013.26e−011.29e + 00UF84.56e−011.34e−019.43e−011.42e + 001.74e−013.24e−018.47e + 00UF94.21e−011.33e−019.57e−011.24e + 001.73e−013.25e−018.41e + 00UF103.53e−011.34e−018.93e−016.05e−011.55e−013.14e−012.06e + 00ZDT18.31e−011.89e−011.00e + 005.94e + 001.41e−012.89e−013.90e + 00ZDT28.31e−011.61e−018.68e−015.60e + 001.43e−012.87e−019.56e + 00ZDT33.42e−011.28e−019.11e−013.64e + 001.39e−013.11e−018.05e + 00ZDT44.60e−011.26e−014.76e−011.19e + 001.38e−013.29e−019.29e + 00ZDT62.34e−011.27e−018.22e−014.87e + 001.43e−012.77e−018.53e + 00


### 6.6. Comparison of IGD Index Convergence

Figure9 shows the change trajectory of IGD index of each algorithm under 10000 evaluation times. Figure 9(a) represents IGD index’s changes of ZDT1 test function between MSMOPSO and evolutionary algorithm ARMOEA, MOEAD, NSGAIISDR, and gNSGAII. It can be seen from Figure 9 that the trajectory of MSMOPSO reaches a stable state before 1000 iterations, indicating that the convergence rate of MSMOPSO is fast. The convergence rates of ARMOEA, MOEAD, and NSGAIISDR are similar. According to the termination status of each algorithm in this process, it can be seen that the IGD index value of MSMOPSO is the minimum, which further illustrates that the MSMOPSO has significant convergence performance. In Figure 9(b), the IGD changes of MSMOPSO compared with MOPSOCD, NMPSO, SMPSO, and MPSOD have been presented. It can be clearly seen from Figure 9(b) that MSMOPSO achieves convergence before the 1000th iteration, and the subsequent iteration process remains stable, while IGD index values of other algorithms tend to be stable slowly after the 7000th iteration. It shows that the convergence speed of MSMOPSO is fast and it has certain stability.Figure 9
(a) Convergence graph of MSMOPSO and MOEAs on ZDT1. (b) Convergence graph of MSMOPSO and MOPSOs on ZDT1.


## 6.1. Comparison of IGD Data

Table4 shows the results of IGD evaluation indexes obtained by this algorithm and other MOPSOs. The standard deviation part is marked by parentheses. It can be obtained from Table 4: among the 22 test functions, 12 of the results of MSMOPSO are better than those of other algorithms, SMPSO has 2 better results, MPSOD has 2 better results, MOPSOCD has one better result, NMPSO has 2 better results, dMOPSO has 0 better result, and CMOPSO has 3 better results. It shows that MSMOPSO has better convergence performance and particle distribution than MOPSOs. In addition, considering the standard deviation of IGD index, the performance of MSMOPSO is more stable than other MOPSOs in the test functions DTLZ1, UF1, UF2, UF6-8, UF10, ZDT2, and ZDT4.Table 4
Compared with IGD index of MOPSOs.
FunctionsSMPSOMPSODMOPSOCDNMPSOCMOPSOdMOPSOMSMOPSODTLZ11.67e + 01 (1.92e + 01)3.87e + 01 (6.23e + 00)3.93e + 01 (1.13e + 01)3.17e + 01 (5.86e + 00)6.39e + 01 (1.61e + 01)9.90e + 00 (9.18e + 00)3.27e − 02 (2.08e − 03)DTLZ29.29e − 02 (1.02e − 02)5.67e − 02 (6.22e − 04)1.00e − 01 (9.27e − 03)7.69e − 02 (3.18e − 03)5.82e − 02 (1.17e − 03)1.45e − 01 (9.05e − 03)2.11e − 01 (2.27e − 02)DTLZ33.11e + 01 (2.98e + 01)1.15e + 02 (1.41e + 01)1.10e + 02 (2.84e + 01)9.31e + 01 (2.72e + 01)1.46e + 02 (3.56e + 01)7.74e + 01 (6.43e + 01)1.40e + 02 (2.45e + 01)DTLZ44.17e − 01 (2.98e + 01)6.26e − 02 (6.85e − 03)3.00e − 01 (5.48e − 02)1.55e − 01 (1.76e − 01)1.22e − 01 (2.24e − 01)3.12e − 01 (2.20e − 02)3.30e − 01 (4.23e − 02)DTLZ55.71e − 03 (3.98e − 04)5.35e − 02 (5.34e − 03)2.28e − 02 (6.99e − 03)1.13e − 02 (1.37e − 03)7.32e − 03 (6.70e − 04)4.42e − 02 (6.78e − 03)6.38e − 02 (1.78e − 02)DTLZ61.17e + 00 (9.91e − 01)2.15e − 01 (2.38e − 01)7.85e − 02 (2.24e − 01)1.45e − 02 (3.06e − 03)6.99e − 02 (2.50e − 01)3.28e − 02 (2.39e − 04)1.15e − 02 (2.62e − 03))DTLZ71.86e − 01 (1.52e − 01)1.65e − 01 (2.10e − 02)8.72e − 02 (6.07e − 03)7.04e − 02 (4.99e − 03)1.68e − 01 (2.06e − 01)1.66e − 01 (1.33e − 01)1.47e − 01 (2.15e − 02)UF13.61e − 01 (7.22e − 02)2.67e − 01 (4.33e − 02)6.52e − 01 (1.56e − 01)1.39e − 01 (8.31e − 02)9.59e − 02 (1.36e − 02)6.52e − 01 (1.01e − 01)8.79e − 02 (1.45e − 02)UF29.97e − 02 (8.61e − 03)1.17e − 01 (1.11e − 02)1.40e − 01 (1.38e − 02)8.44e − 02 (6.99e − 03)6.14e − 02 (6.58e − 03)9.58e − 02 (5.86e − 03)4.04e − 02 (4.69e − 03)UF34.56e − 01 (5.18e − 02)5.00e − 01 (1.43e − 02)3.82e − 01 (7.66e − 2)3.66e − 1 (4.91e − 2)3.96e − 01 (3.27e − 02)3.33e − 01 (5.62e − 03)2.07e − 01 (2.31e − 02)UF41.10e − 01 (3.28e − 03)9.91e − 02 (3.92e − 03)8.06e − 02 (6.81e − 03)6.55e − 02 (9.54e − 03)1.11e − 01 (7.75e − 03)1.35e − 01 (5.92e − 03)7.35e − 02 (1.49e − 02)UF52.92e + 00 (4.91e − 01)2.72e + 00 (2.77e − 01)3.85e + 00 (4.64e − 01)1.64e + 00 (4.96e − 01)8.93e − 01 (1.95e − 01)3.19e + 00 (2.93e − 01)8.66e − 01 (3.20e − 01)UF61.27e + 00 (4.19e − 01)1.33e + 00 (2.01e − 01)3.13e + 00 (8.08e − 01)6.90e − 01 (1.13e − 01)3.98e − 01 (5.98e − 02)2.30e + 00 (5.43e − 01)4.23e − 01 (7.62e − 02)UF73.29e − 01 (1.20e − 01)2.38e − 01 (5.67e − 02)5.91e − 01 (1.19e − 01)2.13e − 01 (1.62e − 01)1.44e − 01 (1.25e − 01)3.59e − 01 (8.49e − 02)4.38e − 02 (4.96e − 03)UF83.90e − 01 (4.54e − 02)5.54e − 01 (5.12e − 02)7.66e − 01 (1.64e − 01)4.51e − 01 (5.65e − 02)5.93e − 01 (9.13e − 02)3.73e − 01 (3.90e − 02)2.88e − 01 (1.06e − 02)UF95.61e − 01 (3.92e − 02)6.68e − 01 (3.63e − 02)8.75e − 01 (1.39e − 01)4.63e − 01 (6.11e − 02)8.90e − 01 (1.33e − 01)5.72e − 01 (3.80e − 02)2.94e − 01 (7.12e − 02)UF102.71e + 00 (5.13e − 01)4.12e + 00 (3.75e − 01)5.17e + 00 (9.18e − 01)1.65e + 00 (3.36e − 01)4.50e + 00 (5.24e − 01)9.40e − 01 (2.83e − 02)3.31e − 01 (5.78e − 02)ZDT15.39e − 02 (6.45e − 02)9.43e − 02 (3.56e − 02)4.89e − 03 (1.04e − 02)3.03e − 02 (1.30e − 02)3.20e − 03 (5.76e − 04)5.26e − 02 (1.53e − 02)5.52e − 03 (7.89e − 04)ZDT21.30e − 01 (2.00e − 01)2.00e − 01 (1.46e − 01)1.68e − 01 (2.40e − 01)4.84e − 02 (7.60e − 02)2.81e − 03 (3.91e − 04)3.98e − 02 (1.07e − 02)6.22e − 03 (1.05e − 03)ZDT31.84e − 01 (9.98e − 02)2.11e − 01 (6.69e − 02)6.59e − 02 (6.10e − 02)8.81e − 02 (2.57e − 02)3.64e − 03 (5.89e − 04)3.50e − 02 (7.06e − 03)2.03e − 01 (5.17e − 03)ZDT49.02e + 01 (3.23e + 01)2.28e + 02 (2.12e + 01)1.57e + 02 (3.25e + 01)1.28e + 02 (3.21e + 01)1.72e + 02 (3.04e + 01)3.07e + 01 (4.78e + 01)5.57e − 03 (8.40e − 04)ZDT61.50e + 00 (1.82e + 00)2.04e + 00 (1.45e + 00)1.62e + 00 (1.91e + 00)2.35e − 03 (2.30e − 04)2.75e − 01 (2.71e − 01)3.18e − 01 (1.02e − 01)2.17e − 03 (2.59e − 04)The bold values in Table 4 represent: under the premise of taking the mean as the main comparison object, the optimal values in the comparison results are compared for different detection functions.Table5 shows the results of IGD evaluation indexes obtained by MSMOPSO and other MOEAs. The standard deviation part is marked by parentheses. As shown in Table 5, 14 of the 22 test functions of MSMOPSO have better results compared with the solving mean of other algorithms, ARMOEA has 3 better results, MOEAD has 1 better result, NSGAIISDR has 2 better results, gNSGAII has 0 better result, RVEA has 1 better result, and GrEA has 1 better result. It shows that MSMOPSO has better convergence performance and particle distribution than MOEAs. In addition, in the test functions DTLZ1, UF1, UF2, UF6-8, UF10, ZDT1, ZDT2, ZDT4, and ZDT6, the standard deviation of IGD index of MSMOPSO is smaller than that of other algorithms, indicating that the convergence performance and particle distribution performance of MSMOPSO are more stable than other MOEAs.Table 5
Compared with IGD index of MOEAs.
FunctionsARMOEAMOEADNSGAIISDRgNSGAIIRVEAGrEAMSMOPSODTLZ13.16e + 00 (1.80e + 00)4.81e + 00 (2.66e + 00)2.60e + 00 (1.40e + 00)2.14e + 01 (5.60e + 00)5.43e + 00 (2.42e + 00)2.84e + 00 (1.78e + 00)3.27e − 02 (2.08e − 03)DTLZ25.50e − 02 (1.70e − 04)5.49e − 02 (2.37e − 04)4.66e − 01 (2.61e − 02)3.86e − 01 (5.48e − 03)5.58e − 02 (7.40e − 04)6.61e − 02 (1.82e − 03)2.11e − 01 (2.27e − 02)DTLZ39.43e + 00 (3.92e + 00)1.54e + 01 (9.61e + 00)8.08e + 00 (3.39e + 00)5.02e + 01 (1.41e + 01)1.66e + 01 (5.95e + 00)8.70e + 00 (4.40e + 00)1.40e + 02 (2.45e + 01)DTLZ43.84e − 01 (3.21e − 01)5.86e − 01 (3.13e − 01)5.21e − 01 (9.07e − 02)3.95e − 01 (1.07e − 02)1.05e − 01 (1.48e − 01)1.57e − 01 (2.46e − 01)3.30e − 01 (4.23e − 02)DTLZ55.73e − 03 (1.39e − 04)3.25e − 02 (6.41e − 04)3.33e − 02 (5.91e − 03)3.11e − 01 (9.11e − 03)7.77e − 02 (9.51e − 03)2.07e − 02 (1.03e − 03)6.38e − 02 (1.78e − 02)DTLZ68.75e − 03 (1.98e − 02)1.99e − 01 (4.08e − 01)7.08e − 02 (6.31e − 02)7.62e − 01 (3.37e − 01)1.22e − 01 (1.07e − 01)2.21e − 02 (5.76e − 04)1.15e − 02 (2.62e − 03)DTLZ72.32e − 01 (2.55e − 01)1.74e − 01 (1.19e − 01)8.71e − 02 (5.40e − 03)4.02e − 01 (3.05e − 03)1.23e − 01 (1.11e − 02)1.03e − 01 (7.95e − 02)1.47e − 01 (2.15e − 02)UF11.22e − 01 (2.88e − 02)3.21e − 01 (7.76e − 02)1.42e − 01 (3.06e − 02)3.35e − 01 (1.99e − 02)1.75e − 01 (4.00e − 02)1.23e − 01 (3.12e − 02)8.79e − 02 (1.45e − 02)UF27.60e − 02 (1.43e − 02)2.09e − 01 (7.56e − 02)8.44e − 02 (5.55e − 03)2.95e − 01 (2.30e − 02)1.20e − 01 (8.96e − 03)6.91e − 02 (5.77e − 03)4.04e − 02 (4.69e − 03)UF34.33e − 01 (2.85e − 02)3.36e − 01 (3.32e − 02)4.81e − 01 (1.55e − 02)4.46e − 01 (1.63e − 02)4.94e − 01 (1.61e − 02)4.62e − 01 (1.62e − 02)2.07e − 01 (2.31e − 02)UF48.29e − 02 (3.18e − 03)1.16e − 01 (5.37e − 03)9.93e − 02 (4.02e − 03)2.20e − 01 (1.43e − 03)1.523e − 01 (4.80e − 03)8.55e − 02 (3.18e − 03)7.35e − 02 (1.49e − 02)UF58.19e − 01 (2.17e − 01)1.47e + 00 (3.21e − 01)1.41e + 00 (3.28e − 01)8.71e − 01 (2.66e − 01)1.16e + 00 (3.73e − 01)9.04e − 01 (2.68e − 01)8.66e − 01 (3.20e − 01)UF65.32e − 01 (1.00e − 01)5.52e − 01 (1.53e − 01)6.85e − 01 (1.28e − 01)5.46e − 01 (1.13e − 01)6.98e − 01 (1.02e − 01)6.55e − 01 (1.52e − 01)4.23e − 01 (7.62e − 02)UF72.29e − 01 (1.32e − 01)4.18e − 01 (1.35e − 01)1.73e − 01 (7.94e − 02)3.12e − 01 (4.69e − 02)2.45e − 01 (6.35e − 02)1.42e − 01 (8.44e − 02)4.38e − 02 (4.96e − 03)UF83.67e − 01 (5.34e − 02)5.87e − 01 (2.74e − 01)3.88e − 01 (3.59e − 02)5.11e − 01 (1.98e − 02)4.16e − 01 (3.58e − 02)3.14e − 01 (3.13e − 02)2.88e − 01 (1.06e − 02)UF94.94e − 01 (7.59e − 02)5.25e − 01 (8.58e − 02)5.33e − 01 (1.02e − 01)5.87e − 01 (4.80e − 02)4.78e − 01 (6.58e − 02)4.66e − 01 (8.41e − 02)2.94e − 01 (7.12e − 02)UF101.13e + 00 (2.93e − 01)7.22e − 01 (7.97e − 02)2.37e + 00 (5.73e − 01)1.25e + 00 (3.83e − 01)1.72e + 00 (5.17e − 01)2.06e + 00 (5.34e − 01)3.31e − 01 (5.78e − 02)ZDT16.41e − 02 (4.55e − 02)1.85e − 01 (7.29e − 02)1.39e − 01 (2.37e − 02)3.28e − 01 (2.38e − 02)2.01e − 01 (2.99e − 02)5.34e − 02 (1.05e − 02)5.52e − 03 (7.89e − 04)ZDT26.97e − 01 (3.20e − 02)5.73e − 01 (5.82e − 02)2.05e − 01 (4.00e − 02)2.40e − 01 (2.59e − 02)2.87e − 01 (4.18e − 02)7.94e − 02 (1.45e − 02)6.22e − 03 (1.05e − 03)ZDT35.40e − 02 (1.75e − 02)1.58e − 01 (6.72e − 02)1.06e − 01 (1.45e − 02)3.38e − 01 (3.42e − 02)1.89e − 01 (1.42e − 02)4.84e − 02 (1.42e − 02)2.03e − 01 (5.17e − 03)ZDT44.79e + 01 (6.75e + 00)1.73e + 01 (5.14e + 00)4.28e + 01 (6.28e + 00)4.79e + 01 (7.66e + 00)5.57e + 01 (1.08e + 1)3.13e + 01 (5.09e + 00)5.57e − 03 (8.40e − 04)ZDT63.46e + 00 (1.57e − 01)2.21e + 00 (6.45e − 01)3.60e + 00 (2.04e − 01)3.32e + 00 (2.12e − 01)3.60e + 00 (1.74e − 01)2.99e + 00 (1.89e − 01)2.17e − 03 (2.59e − 04)The bold values in Table 5 represent under the premise of taking the mean as the main comparison object, the optimal values in the comparison results are compared for different detection functions.Based on the comparison of IGD indexes in Table4 and Table 5, it is found that MSMOPSO performs relatively well both for MOPSOs and MOEAs, especially when detection functions are DTLZ1, ZDT4, ZDT6, and UF series (UFs). It is shown that the Gaussian mutation strategy increases the possibility of finding better solutions, and more optimal solutions are stored by the multi-archiving mechanism, which can further improve the exploration ability of the algorithm. This can enable MSMOPSO to find as many noninferior solutions as possible. In addition, according to the statistical results, the variance of IGD index value of the proposed algorithm for MOPs of DTLZ1, UF2, UF7, UF8, ZDT4, and ZDT6 is smaller, indicating that the algorithm has certain stability in solving specific optimization problems, and has some advantages in convergence to the optimal solution set.


## 6.2. Comparison of HV Data

Table6 shows the solution results of HV evaluation index of MSMOPSO and other MOPSOs. The standard deviation part is marked by parentheses. As shown in Table 6, among the 22 test functions, 10 of the results of MSMOPSO are better than those of other algorithms. SMPSO has 2 better results, MPSOD has 1 better result, MOPSOCD has 0 better result, NMPSO has 5 better results, CMOPSO has 4 better results, and dMOPSO has 0 better result. It indicates that the results of MSMOPSO are more widely distributed in objective space. In addition, in the test functions UF1, UF2, UF7, ZDT2, and ZDT3 of MSMOPSO, the standard deviation of HV index is smaller than that of other algorithms; it illustrates that MSMOPSO has better robustness compared with other algorithms.Table 6
Compared with HV index of MOPSOs.
FunctionsSMPSOMPSODMOPSOCDNMPSOCMOPSOdMOPSOMSMOPSODTLZ11.06e − 01 (2.06e − 01)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)1.20e − 02 (5.58e − 02)8.19e − 01 (4.59e − 03)DTLZ24.55e − 01 (2.06e − 02)5.47e − 01 (2.14e − 03)4.60e − 01 (1.96e − 02)5.59e − 01 (1.08e − 03)5.40e − 01 (2.89e − 03)3.67e − 01 (1.45e − 02)3.07e − 01 (2.74e − 02)DTLZ32.82e − 02 (4.75e − 02)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)DTLZ43.29e − 01 (1.06e − 01)5.38e − 01 (1.11e − 02)3.25e − 01 (6.09e − 02)5.23e − 01 (8.30e − 02)5.04e − 01 (1.12e − 01)3.09e − 01 (4.92e − 02)3.36e − 01 (3.46e − 02)DTLZ51.98e − 01 (4.12e − 04)1.46e − 01 (7.08e − 03)1.84e − 01 (5.62e − 03)1.97e − 01 (5.42e − 04)1.97e − 01 (5.06e − 04)1.50e − 01 (1.05e − 02)1.18e − 01 (1.78e − 02)DTLZ65.68e − 02 (8.97e − 02)1.19e − 01 (5.63e − 02)1.85e − 01 (4.43e − 02)1.96e − 01 (1.07e − 03)1.87e − 01 (5.08e − 02)1.83e − 01 (2.07e − 04)1.90e − 01 (6.17e − 03)DTLZ72.48e − 01 (1.77e − 02)2.16e − 01 (1.57e − 02)2.67e − 01 (3.25e − 03)2.75e − 01 (1.46e − 03)2.63e − 01 (2.13e − 02)2.42e − 01 (1.07e − 02)2.48e − 01 (6.98e − 03)UF12.79e − 01 (6.38e − 02)3.60e − 01 (4.75e − 02)8.04e − 02 (8.04e − 02)5.25e − 01 (4.72e − 02)5.73e − 01 (1.94e − 02)6.90e − 02 (6.61e − 02)5.97e − 01 (2.04e − 02)UF26.10e − 01 (7.96e − 03)5.75e − 01 (1.26e − 02)5.45e − 01 (1.72e − 02)6.17e − 01 (8.01e − 03)6.46e − 01 (7.19e − 03)6.15e − 01 (5.56e − 03)6.73e − 01 (6.53e − 03)UF32.07e − 01 (4.61e − 02)1.74e − 01 (1.31e − 02)2.50e − 01 (5.55e − 2)2.78e − 1 (4.07e − 2)2.62e − 01 (2.03e − 02)3.06e − 01 (7.11e − 03)4.69e − 01 (2.76e − 02)UF42.94e − 01 (8.10e − 03)3.08e − 01 (5.05e − 03)3.31e − 01 (8.75e − 03)3.55e − 01 (1.34e − 02)2.90e − 01 (9.55e − 03)2.55e − 01 (7.03e − 03)3.48e − 01 (8.30e − 03)UF50.00e + 00 (0.00e + 00)0.00e + 00 (0.00e − 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)8.47e − 03 (2.26e − 02)0.00e + 00 (0.00e + 00)3.04e − 03 (1.07e − 02)UF62.03e − 04 (1.02e − 03)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)1.51e − 02 (1.94e − 02)1.35e − 01 (7.59e − 02)0.00e + 00 (0.00e + 00)6.67e − 02 (1.95e − 03)UF72.18e − 01 (9.91e − 02)2.68e − 01 (6.41e − 02)3.06e − 02 (3.55e − 02)3.59e − 01 (1.12e − 01)4.31e − 01 (8.88e − 02)1.92e − 01 (6.16e − 02)5.22e − 01 (7.28e − 03)UF81.65e − 01 (4.59e − 02)5.64e − 02 (1.98e − 02)8.24e − 03 (1.28e − 02)2.86e − 01 (4.64e − 02)1.76e − 02 (1.65e − 02)2.61e − 01 (2.38e − 02)2.94e − 01 (1.52e − 02)UF92.05e − 01 (3.82e − 02)1.06e − 01 (2.27e − 02)3.12e − 02 (3.41e − 02)3.13e − 01 (5.37e − 02)2.34e − 02 (2.83e − 02)2.29e − 01 (1.55e − 02)4.48e − 01 (7.20e − 02)UF100.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)9.10e − 02 (8.14e − 04)2.92e − 01 (2.24e − 02)ZDT16.54e − 01 (8.02e − 02)5.84e − 01 (4.85e − 02)7.18e − 01 (1.51e − 02)6.88e − 01 (1.72e − 02)7.20e − 01 (8.85e − 04)6.58e − 01 (1.70e − 02)7.17e − 01 (1.83e − 03)ZDT23.41e − 01 (1.49e − 01)2.34e − 01 (1.14e − 01)3.20e − 01 (1.59e − 01)4.08e − 01 (7.01e − 02)4.45e − 01 (7.18e − 04)3.85e − 01 (1.69e − 02)4.40e − 01 (2.47e − 03)ZDT35.31e − 01 (7.24e − 02)4.57e − 01 (6.08e − 02)5.59e − 01 (4.33e − 02)5.72e − 01 (9.34e − 03)6.00e − 01 (1.23e − 03)6.01e − 01 (1.12e − 02)6.53e − 01 (2.75e − 03)ZDT40.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)3.44e − 02 (4.87e − 02)7.17e − 01 (1.83e − 03)ZDT61.97e − 01 (1.75e − 01)4.72e − 02 (8.39e − 02)1.69e − 01 (1.62e − 01)3.90e − 01 (1.99e − 04)1.87e − 01 (1.56e − 01)1.23e − 01 (6.93e − 02)3.83e − 01 (1.62e − 03)Table7 shows the solution results of the HV evaluation index of MSMOPSO and other MOEAs, and the standard deviation part is marked by parentheses. As shown in Table 7, 13 of the 22 test functions of this algorithm have better results compared with the solving mean of other algorithms, ARMOEA has 4 better results, MOEAD has 1 better result, NSGAIISDR has 1 better result, GrEA has 1 better result, and RVEA has 1 better result. It shows that the spatial distribution of MSMOPSO is better than other MOEAs. In addition, in the test functions UF1, UF2, UF7, ZDT1, ZDT2, and ZDT3 of this algorithm, the standard deviation of HV index is smaller than that of other algorithms, indicating that this algorithm has good stability compared with other algorithms.Table 7
Compared with HV index of MOEAs.
FunctionsARMOEAMOEADNSGAIISDRgNSGAIIRVEAGrEAMSMOPSODTLZ10.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)8.19e − 01 (4.59e − 03)DTLZ25.56e − 01 (5.76e − 04)5.55e − 01 (9.01e − 04)2.55e − 01 (1.45e − 02)1.59e − 01 (3.54e − 04)5.53e − 01 (1.33e − 03)5.58e − 01 (9.00e − 04)3.07e − 01 (2.74e − 02)DTLZ30.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)DTLZ44.03e − 01 (1.59e − 01)2.92e − 01 (1.63e − 01)2.85e − 01 (5.48e − 02)1.60e − 01 (2.75e − 04)5.31e − 01 (6.62e − 02)5.11e − 01 (1.28e − 01)3.36e − 01 (3.46e − 02)DTLZ51.98e − 01 (2.20e − 04)1.82e − 01 (9.22e − 04)1.80e − 01 (3.08e − 03)6.02e − 02 (3.55e − 05)1.50e − 01 (7.73e − 03)1.88e − 01 (4.74e − 04)1.18e − 01 (1.78e − 02)DTLZ61.96e − 01 (1.63e − 02)1.46e − 01 (6.73e − 02)1.63e − 01 (2.19e − 02)9.79e − 03 (1.74e − 02)1.21e − 01 (3.90e − 02)1.88e − 01 (2.63e − 04)1.90e − 01 (6.17e − 03)DTLZ72.55e − 01 (2.51e − 02)2.45e − 01 (9.61e − 03)2.69e − 01 (1.76e − 03)7.56e − 02 (6.89e − 04)2.48e − 01 (7.73e − 03)2.70e − 01 (1.06e − 02)2.48e − 01 (6.98e − 03)UF15.57e − 01 (2.76e − 02)4.14e − 01 (5.21e − 02)5.09e − 01 (4.17e − 02)3.36e − 01 (5.33e − 02)4.65e − 01 (5.23e − 02)5.31e − 01 (4.03e − 02)5.97e − 01 (2.04e − 02)UF26.32e − 01 (9.92e − 03)5.61e − 01 (3.29e − 02)6.09e − 01 (8.31e − 03)4.30e − 01 (4.61e − 02)5.60e − 01 (1.12e − 02)6.32e − 01 (7.17e − 03)6.73e − 01 (6.53e − 03)UF32.05e − 01 (2.27e − 02)3.06e − 01 (4.61e − 02)1.91e − 01 (1.41e − 02)2.10e − 1 (1.17e − 02)1.78e − 01 (1.12e − 02)2.03e − 01 (1.29e − 02)4.69e − 01 (2.76e − 02)UF43.31e − 01 (4.76e − 03)2.84e − 01 (6.97e − 03)3.08e − 01 (4.66e − 03)2.42e − 01 (1.51e − 03)2.30e − 01 (5.14e − 03)3.27e − 01 (4.17e − 03)3.48e − 01 (8.30e − 03)UF51.42e − 02 (2.84e − 02)0.00e + 00 (0.00e − 00)0.00e + 00 (0.00e + 00)7.42e − 03 (1.82e − 02)5.24e − 03 (1.62e − 02)6.06e − 03 (2.03e − 02)3.04e − 03 (1.07e − 02)UF65.25e − 02 (3.75e − 02)7.61e − 02 (6.40e − 02)1.28e − 02 (1.86e − 02)4.46e − 02 (3.49e − 02)2.03e − 02 (2.10e − 02)1.99e − 02 (2.55e − 02)6.67e − 02 (1.95e − 03)UF73.55e − 01 (8.09e − 02)2.54e − 01 (7.49e − 02)3.60e − 01 (7.25e − 02)2.52e − 01 (3.90e − 02)2.82e − 01 (5.49e − 02)3.97e − 01 (8.12e − 02)5.22e − 01 (7.28e − 03)UF82.96e − 01 (2.05e − 02)1.54e − 01 (6.00e − 02)1.03e − 01 (4.94e − 02)7.58e − 02 (1.93e − 02)7.63e − 02 (4.86e − 02)2.27e − 01 (5.17e − 02)2.94e − 01 (1.52e − 02)UF92.86e − 01 (5.22e − 02)2.88e − 01 (4.81e − 02)2.03e − 01 (8.42e − 02)1.38e − 01 (2.10e − 02)2.61e − 01 (5.10e − 02)2.74e − 01 (7.82e − 02)4.48e − 01 (7.20e − 02)UF100.00e + 00 (0.00e + 00)3.68e − 02 (2.43e − 02)0.00e + 00 (0.00e + 00)1.01e − 07 (5.51e − 07)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)2.92e − 01 (2.24e − 02)ZDT16.47e − 01 (3.05e − 02)5.39e − 01 (6.01e − 02)5.44e − 01 (2.66e − 02)3.51e − 01 (2.92e − 02)4.61e − 01 (3.19e − 02)6.51e − 01 (1.37e − 02)7.17e − 01 (1.83e − 03)ZDT21.09e − 03 (3.71e − 03)8.82e − 02 (1.55e − 02)2.02e − 01 (3.35e − 02)2.20e − 01 (2.45e − 02)1.31e − 01 (2.89e − 02)3.35e − 01 (1.95e − 02)4.40e − 01 (2.47e − 03)ZDT35.88e − 01 (4.01e − 02)5.53e − 01 (5.15e − 02)5.34e − 01 (2.23e − 02)5.84e − 01 (5.63e − 02)4.51e − 01 (1.53e − 02)5.67e − 01 (8.99e − 03)6.53e − 01 (2.75e − 03)ZDT40.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)7.17e − 01 (1.83e − 03)ZDT60.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)0.00e + 00 (0.00e + 00)3.83e − 01 (1.62e − 03)The bold values in Table 7 represent under the premise of taking the mean as the main comparison object, the optimal values in the comparison results are compared for different detection functions.Based on Table6 and Table 7, it is found that compared with MOPSOs and MOEAs, the HV index value of this algorithm is relatively large in DTLZ1, UF1-3, UF7-10, and ZDT3-4 detection functions. According to the characteristics of HV index, the larger the HV value of this algorithm, the better the spatial distribution of decision variables in it. On this basis, statistical results show that the variance of HV index value of MSMOPSO is smaller for DTLZ1, UF7, and ZDT4 test functions, indicating that it has better stability and significant solving ability when solving these two kinds of MOPs. In conclusion, the application of Gaussian mutation strategy can improve the spatial distribution of particles in this proposed algorithm, and then the multi-archive is used to store the noninferior solutions, further indicating that these strategies can effectively improve the development ability of this algorithm for the target space and effectively alleviate the disadvantage of the algorithm easily falling into local convergence.


## 6.3. Comparison of Box Plots

In this part, some test functions are taken as examples to store the IGD index of each running result of MSMOPSO, which is displayed in the form of box graph and compared with the calculation results of other comparison algorithms. It is found that the distance between the upper and the lower quartile loci of MSMOPSO is closer than that of other algorithms, indicating that the distribution of experimental data is more concentrated, so MSMOPSO is more stable. Meanwhile, for outliers in the data, “+” is used to represent them in the box charts. It can be seen from these plots that there are fewer outliers in the experimental data of MSMOPSO, indicating that the solution results are better. To sum up, the performance of MSMOPSO is more significant than other algorithms, both in terms of solution results and algorithm stability. The box plots are shown in Figure6, where 1 stands for MPSOD, 2 stands for MOPSOCD, 3 stands for SMPSO, 4 stands for NMPSO, 5 stands for ARMOEA, 6 stands for MOEAD, 7 stands for NSGAIISDR, 8 stands for gNSGAII, and 9 stands for MSMOPSO.Figure 6
(a) Boxplot for the statistical results for ZDT1. (b) Boxplot for the statistical results for ZDT2. (c) Boxplot for the statistical results for ZDT4. (d) Boxplot for the statistical results for ZDT6. (e) Boxplot for the statistical results for UF1. (f) Boxplot for the statistical results for UF2. (g) Boxplot for the statistical results for UF3. (h) Boxplot for the statistical results for UF6. (i) Boxplot for the statistical results for UF7. (j) Boxplot for the statistical results for UF8. (k) Boxplot for the statistical results for UF10. (l) Boxplot for the statistical results for DTLZ1.


## 6.4. Comparison of Results

In order to verify the authenticity and validity of the experimental results, the optimal solution set of each algorithm in the objective space is displayed on the front surface of each algorithm, represented by a two-dimensional detection function ZDT4 and a three-dimensional detection function DTLZ1, as shown in Figure7 and Figure 8. The points in the figure represent the optimal solution set of the algorithm on test functions, and the Pareto front is represented by shadow.Figure 7
(a) Results of MPSOD on DTLZ1. (b) Results of MOPSOCD on DTLZ1. (c) Results of SMPSO on DTLZ1. (d) Results of NMPSO on DTLZ1. (e) Results of ARMOEA on DTLZ1. (f) Results of MOEAD on DTLZ1. (g) Results of NSGAIISDR on DTLZ1. (h) Results of gNSGAII on DTLZ1. (i) Results of MSMOPSO on DTLZ1.Figure 8
(a) Results of MPSOD on ZDT4. (b) Results of MOPSOCD on ZDT4. (c) Results of SMPSO on ZDT4. (d) Results of NMPSO on ZDT4. (e) Results of ARMOEA on ZDT4. (f) Results of MOEAD on ZDT4. (g) Results of NSGAIISDR on ZDT4. (h) Results of gNSGAII on ZDT4. (i) Results of MSMOPSO on ZDT4.Comparing Figures7 and 8, it is found that the algorithm is evenly distributed on the ZDT4 and the optimal solution set is completely distributed on the Pareto front, while other test algorithms are not attached to the front surface. In addition, due to the large numerical value of solution results obtained by other test functions, there is a large difference between them and the optimal solution set that forms the front surface, so it is not shown in other comparison diagrams. For the DTLZ1, only MSMOPSO is evenly distributed on the Pareto front surface, while the optimal solution sets of ARMOEA, MOEAD, and NSGAIISDR are too scattered, so they do not converge to the front surface. The solution values of MPSOD, MOPSOCD, SMPSO, NMPSO, and gNSGAII are large and the distribution of the solution set is too scattered, so the Pareto front does not appear.


## 6.5. Complexity of MSMOPSO

In MSMOPSO, the first archive is updated by elite selection strategy, self-organizing clustering strategy and species migration strategy are applied to expand and update the second archive, and finally, the particles of the whole swarm are randomly disturbed by Gaussian mutation strategy. Therefore, the time complexity of this algorithm mainly depends on the above three aspects. For the first archive, the elite strategy is applied to select the particles in the archive and the swarm, so the time complexity of this process isOP∗N∗M (P represents the archive size, N represents the population size, and M represents the dimension of decision variables). For the second archive, self-organizing clustering strategy is used to select and store the particles within swarm according to the grouping center, then physical migration strategy is utilized in updating this archive, so the time complexity of the second archive is divided into two parts, one is a self-organizing clustering strategy, the other is the species migration strategy. The time complexity of the first part is ON∗N∗M, and the second part is ON∗N∗M+P∗N∗M. Because N≥P, the complexity of the second archive is ON∗N∗M. For Gaussian mutation strategy, the computational complexity is ON∗D∗M∗T (T represents the number of iterations of the algorithm). So MSMOPSO’s time complexity is OP∗N∗M+ON∗N∗M+N∗N∗M+ON∗D∗M∗T, that is ON2M+OPNM+ONDMT.In addition, the (tic, toc) order was used to record the running time of MSMOPSO in 1000 evaluations, and the results were compared with MOEAs and MOPSOs. The comparison results are shown in Table8 and Table 9. As can be seen from the table, the calculation time of MSMOPSO for different detection functions is larger than that of other algorithms, but the difference is relatively small. It is worth noting that these results of running time depend on the influence of computer hardware, software, and other factors, so the running result is only for reference.Table 8
Compared with runtime of MOEAs.
FunctionsARMOEAMOEADNSGAIISDRgNSGAIIRVEAGrEAMSMOPSODTLZ17.55e−011.53e + 002.52e−011.51e−012.35e−013.33e−017.13e + 00DTLZ23.81e + 001.50e + 002.72e−011.42e−012.31e−018.10e−017.92e + 00DTLZ34.18e−011.53e + 002.51e−011.51e−012.35e−012.89e−016.44e + 00DTLZ43.07e + 001.55e + 003.00e−011.58e−012.35e−017.60e−011.78e + 00DTLZ53.07e + 001.53e + 002.42e−011.45e−011.98e−016.78e−018.00e + 00DTLZ62.71e + 001.56e + 002.50e−011.55e−012.28e−016.24e−014.78e + 01DTLZ72.97e + 001.50e + 002.43e−011.48e−012.07e−017.36e−015.78e−01UF13.37e−011.74e + 003.04e−011.52e−012.26e−012.73e−016.03e−01UF28.28e−011.86e + 003.12e−011.60e−012.34e−013.13e−011.81e + 00UF33.67e−011.82e + 003.19e−011.567e−012.40e−012.87e−011.34e + 00UF49.05e−011.77e + 003.04e−011.41e−012.34e−013.76e−013.51e + 00UF52.86e−011.75e + 003.10e−011.52e−012.23e−012.60e−011.28e + 00UF62.93e−011.76e + 003.07e−011.52e−012.27e−012.62e−011.37e + 00UF73.79e−011.75e + 003.04e−011.48e−012.27e−012.77e−011.29e + 00UF81.67e + 001.93e + 003.05e−011.54e−012.12e−014.55e−018.47e + 00UF91.17e + 001.93e + 003.08e−011.55e−012.15e−014.16e−018.41e + 00UF108.41e−011.81e + 003.08e−011.56e−012.12e−014.15e−012.06e + 00ZDT15.750e−011.68e + 003.01e−011.53e−012.22e−012.95e−013.90e + 00ZDT23.07e−011.67e + 002.99e−011.44e−012.16e−012.63e−019.56e + 00ZDT37.90e−111.70e + 002.99e−011.49e−012.28e−013.17e−018.05e + 00ZDT42.68e−111.66e + 003.03e−011.54e−012.14e−012.54e−019.29e + 00ZDT62.62e−11.68e + 003.02e−011.51e−112.16e−012.50e−018.53e + 00Table 9
Compared with runtime of MOPSOs.
FunctionsCMOPSOMOPSOCDMPSODNMPSOSMPSOdMOPSOMSMOPSODTLZ12.61e−011.31e−011.30e + 001.91e + 001.42e−013.84e−017.13e + 00DTLZ21.30e + 001.57e−011.30e + 002.90e + 001.61e−013.45e−017.92e + 00DTLZ32.53e−011.27e−011.29e + 009.83e−011.38e−013.78e−016.44e + 00DTLZ49.54e−011.22e−011.07e + 002.96e + 001.49e−013.63e−011.78e + 00DTLZ57.57e−011.28e−018.14e−012.30e + 001.60e−013.61e−018.00e + 00DTLZ67.24e−012.20e−011.18e + 003.68e + 001.63e−013.73e−014.78e + 01DTLZ77.68e−012.31e−016.14e−012.98e + 001.53e−013.24e−015.78e−01UF12.36e−011.28e−011.16e + 005.34e−011.41e−013.39e−016.03e−01UF22.74e−011.34e−011.20e + 008.82e−011.59e−013.70 e−011.81e + 00UF32.62e−012.12e−011.16e + 001.07e + 001.61e−013.48e−011.34e + 00UF42.72e−011.40e−011.27e + 002.77e + 001.53e−13.67e−013.51e + 00UF52.28e−011.32e−011.07e + 004.97e−011.45e−013.61e−011.28e + 00UF62.30e−011.34e−011.10e + 005.11e−011.48e−013.61e−011.37e + 00UF72.52e−011.31e−011.14e + 006.64e−011.45e−013.26e−011.29e + 00UF84.56e−011.34e−019.43e−011.42e + 001.74e−013.24e−018.47e + 00UF94.21e−011.33e−019.57e−011.24e + 001.73e−013.25e−018.41e + 00UF103.53e−011.34e−018.93e−016.05e−011.55e−013.14e−012.06e + 00ZDT18.31e−011.89e−011.00e + 005.94e + 001.41e−012.89e−013.90e + 00ZDT28.31e−011.61e−018.68e−015.60e + 001.43e−012.87e−019.56e + 00ZDT33.42e−011.28e−019.11e−013.64e + 001.39e−013.11e−018.05e + 00ZDT44.60e−011.26e−014.76e−011.19e + 001.38e−013.29e−019.29e + 00ZDT62.34e−011.27e−018.22e−014.87e + 001.43e−012.77e−018.53e + 00


## 6.6. Comparison of IGD Index Convergence

Figure9 shows the change trajectory of IGD index of each algorithm under 10000 evaluation times. Figure 9(a) represents IGD index’s changes of ZDT1 test function between MSMOPSO and evolutionary algorithm ARMOEA, MOEAD, NSGAIISDR, and gNSGAII. It can be seen from Figure 9 that the trajectory of MSMOPSO reaches a stable state before 1000 iterations, indicating that the convergence rate of MSMOPSO is fast. The convergence rates of ARMOEA, MOEAD, and NSGAIISDR are similar. According to the termination status of each algorithm in this process, it can be seen that the IGD index value of MSMOPSO is the minimum, which further illustrates that the MSMOPSO has significant convergence performance. In Figure 9(b), the IGD changes of MSMOPSO compared with MOPSOCD, NMPSO, SMPSO, and MPSOD have been presented. It can be clearly seen from Figure 9(b) that MSMOPSO achieves convergence before the 1000th iteration, and the subsequent iteration process remains stable, while IGD index values of other algorithms tend to be stable slowly after the 7000th iteration. It shows that the convergence speed of MSMOPSO is fast and it has certain stability.Figure 9
(a) Convergence graph of MSMOPSO and MOEAs on ZDT1. (b) Convergence graph of MSMOPSO and MOPSOs on ZDT1.


## 7. Conclusions

In this paper, we propose a new MOPSO called MSMOPSO to deal with MOPs. In MSMOPSO, different archiving strategies and Gaussian mutation strategy are applied. Then, it was compared with MOPSOs and MOEAs on different detection functions. In addition, we use box plots, IGD index, HV index, and front diagrams to test its performance. Through IGD and HV indexes, it was found that MSMOPSO had stronger convergence ability and more significant robustness. Through the comparison between box plots and front diagrams, it is found that the solution effect of MSMOPSO is better. What is more, we can see the convergence rate of MSMOPSO is faster by the IGD index trajectory graph. The results show that different update strategies adopted for different archives can improve the shortcoming of MOPSO which is easy to fall into local convergence. And, Gaussian mutation strategy increases the density of the particles in the population.Further research on MSMOPSO can be divided into two parts: algorithm optimization research and practical application research. Most researchers apply MOPSO to solve optimization problems in static environment and then solve MOPs in dynamic environment. Therefore, improving MSMOPSO and using it to solve the optimization problem in a dynamic environment will be one of the focuses of future work. Another is to test the performance of MSMOPSO in real-world areas such as the smart grid [42].



---

*Source: 7372450-2022-05-17.xml*