# Multi-Swarm Multi-Objective Optimizer Based onp-Optimality Criteria for Multi-Objective Portfolio Management

Authors: Yabao Hu; Hanning Chen; Maowei He; Liling Sun; Rui Liu; Hai Shen
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Portfolio management is an important technology for reasonable investment, fund management, optimal asset allocation, and effective investment. Portfolio optimization problem (POP) has been recognized as an NP-hard problem involving numerous objectives as well as constraints. Applications of evolutionary algorithms and swarm intelligence optimizers for resolving multi-objective POP (MOPOP) have attracted considerable attention of researchers, yet their solutions usually convert MOPOP to POP by means of weighted coefficient method. In this paper, a multi-swarm multi-objective optimizer based onp-optimality criteria calledp-MSMOEAs is proposed that tries to find all the Pareto optimal solutions by optimizing all objectives at the same time, rather than through the above transforming method. The proposedp-MSMOEAs extended original multiple objective evolutionary algorithms (MOEAs) to cooperative mode through combiningp-optimality criteria and multi-swarm strategy. Comparative experiments ofp-MSMOEAs and several MOEAs have been performed on six mathematical benchmark functions and two portfolio instances. Simulation results indicate thatp-MSMOEAs are superior for portfolio optimization problem to MOEAs when it comes to optimization accuracy as well as computation robustness.


---

## Body

## 1. Introduction

In real life, most optimization problems are multi-objective optimization problem (MOP). In MOP, numerous contradictory objectives which are subject to several certain constraints must be optimized at the same time. In general, the common approach to solve MOP is to find Pareto optimal set. And Pareto optimality theory indicates that if no vectorx∈C makes F(x)<F(x∗), vector x∗∈C is Pareto optimality. Pareto optimal set could improve at least one objective without deteriorating other objectives. Currently, MOP has made progress in both theory and application [1, 2], yet there are still many challenges because of its complexity. As an effective method to solve MOP, MOEAs could obtain solutions with good diversity [3]. Consequently, they are widely found in various practical applications, such as power dispatching [4], financial management [5], electric machine designing [6], and spectrum allocation [7].As a problem-solving technique, the main strengths of MOEAs are their fast convergence and efficient search ability. However, it is still a challenge to overcome the local convergence or balance diversity and convergence of swarms in MOP for the researchers. Many researchers have attempted to settle this issue with the ideas of optimality criterion and multi-swarm strategy. Optimality criterion could find the best solution according to the different needs of the problem. Yet the conventional Pareto method does not have the advantage of such flexibility. Meanwhile, multi-swarm strategy is prominent in improving diversity.Ink-optimality employed in literature [8], the dominating solution could have an inferior performance on some particular objectives, which is acceptable to decision maker. The number of objectives is prescribed previously by means of the adjustment ofk value. This is difficult for other methods. In literature [9], Optimality Criteria method is employed to seek the optimum through finding a solution that satisfies some prespecified criteria, which are postulated to the corresponding optimal result for the problem. However, only applying optimality criteria, multi-objective algorithms may suffer from excessive loss of diversity. Indeed, multi-swarm strategy could restrict such rapid convergence and increase diversity effectively due to cooperation and exchange between swarms. In literatures [10–13], several MOEAs introduce a multi-swarm strategy. These proposed algorithms contain multiple slave swarms, and the quantity of slave swarms is equal to that of objective functions. During evolution, every slave swarm is dedicated to optimizing a certain objective to discover its non-dominated solutions. An improved particle swarm optimization (PSO) algorithm is proposed in literature [14]. It is shown that whole swarm is stochastically separated into some small-scale sub-swarms. The swarm is reorganized stochastically everyR generations. In that case, the good information obtained by each sub-swarm could be exchanged. In hybrid multi-swarm PSO algorithm shown in literature [15], PSO and differential evolution approach are employed during evolution. It is worth noting that all above algorithms utilize only one strategy to improve the algorithm and suffer from the balance of diversity and convergence, which have the potential for improvement.In this paper, two classes of strategies combining MOEAs are adopted to solve MOP,p-optimality criteria [3], and multi-swarm strategy.P-optimality criteria are employed in selection operator during evolution. These criteria have the ability of determining the most feasible solutions among ones located in the same non-domination rank. Therefore, they have the potential to better the convergence of feasible solutions in late stages of evolution. Meanwhile, competition and cooperation techniques among sub-swarms are designed in multi-swarm strategy. Distribution list and replacement list are devised, guaranteeing the interaction between swarms. Multiple swarms are utilized for optimizing objectives, and each separated sub-swarm employs the MOEAs combiningp-optimality criteria to find out all the non-dominated solutions. Multi-swarm strategy enables improving the diversity of feasible solutions and preventing local convergence. A multi-swarm multi-objective optimizer combining the above strategies, called multi-swarm multi-objective optimization evolutionary algorithms based onp-optimality criteria (p-MSMOEAs), is proposed. Several groups of experiment are conducted to evaluate the performance ofp-MSMOEAs and MOEAs.Simultaneously,p-MSMOEAs are employed to solve MOPOP. Markowitz puts forward the classic mean-variance (MV) portfolio model, which became the theoretical basis of modern portfolios [16, 17]. In Markowitz’s approach, the expected return of portfolio and the risk of portfolio (represented by mean and variance of assets, respectively) are described as two criteria of portfolio model. In the presence of the above two criteria, there is a tradeoff between risk and return. Since Markowitz proposed the MV theory, a large number of studies have been done to extend or modify the basic model in different directions [18–20]. However, the study of multi-objective model considering expected cost of assets is a little rare. In this paper, a three-objective portfolio model considering expected return, expected cost, and risk is designed to solve MOPOP. In the proposed model, expected return is measured by mean of assets and should be maximized, risk is weighed by semi-variance of assets and should be minimized, and expected cost is represented by Euclidean distance of weight vectors and should be minimized. Afterwards, it is optimized byp-MSMOEAs.The remainder of this paper is structured as follows. Section2 presents three MOEAs. In Section 3, after an introduction top-optimality criteria and multi-swarm strategy,p-MSMOEAs in this study are pictured. In Section 4, performance ofp-MSMOEAs and MOEAs on six commonly used multi-objective benchmark functions will be shown. In Section 5, fulfillment ofp-MSMOEAs for MOPOP is presented. Finally, Section 6 summarizes the paper.


## 2. Multi-Objective Optimization Evolutionary Algorithms

### 2.1. NSGA-II Algorithm

NSGA-II algorithm named by Deb et al. has been identified as a computationally fast MOEA. And because of its simplicity, excellent diversity, and convergence of feasible solutions, this algorithm is proved to be one of the most efficient MOEAs. It is prominent in two aspects: fast non-dominated sorting for individuals and elitist selection [21]. The fast non-dominated sorting depends on the indexes of non-dominated rank as well as crowding distance. Further information about process of NSGA-II algorithms could be seen in [22].


### 2.2. MODE Algorithm

MODE algorithm is a simple and powerful MOEA for MOP over a continuous domain. The outstanding advantages of MODE are its speed and robustness. MODE is mainly composed of three components: mutation, Pareto-based evaluation, and selection. Among them, Pareto-based evaluation is the same as NSGA-II. In addition, each vector of the individual undergoes a mutation process with certain mutation probabilitypm. At last, there is a parameter σcrowd in selection operator that could indicate the distance between a solution and its surrounding solution in objective space. Further information about the process of MODE algorithms could be seen in [23].


### 2.3. MOEA/D Algorithm

The main idea of MOEA/D algorithm lies in decomposing MOP into multiple scalar sub-problems and optimizing them at the same time [24]. Every sub-problem utilizes the information of its neighboring sub-problems, which reduces its computational complexity of each generation. Scalarizing functions which are provided with uniformly distributed weight vectors is the fitness evaluation condition of MOEA/D [25]. In the computational experiments of this paper, the Tchebycheff function is employed to decompose. Further information about the process of MOEA/D algorithms could be seen in [24].


## 2.1. NSGA-II Algorithm

NSGA-II algorithm named by Deb et al. has been identified as a computationally fast MOEA. And because of its simplicity, excellent diversity, and convergence of feasible solutions, this algorithm is proved to be one of the most efficient MOEAs. It is prominent in two aspects: fast non-dominated sorting for individuals and elitist selection [21]. The fast non-dominated sorting depends on the indexes of non-dominated rank as well as crowding distance. Further information about process of NSGA-II algorithms could be seen in [22].


## 2.2. MODE Algorithm

MODE algorithm is a simple and powerful MOEA for MOP over a continuous domain. The outstanding advantages of MODE are its speed and robustness. MODE is mainly composed of three components: mutation, Pareto-based evaluation, and selection. Among them, Pareto-based evaluation is the same as NSGA-II. In addition, each vector of the individual undergoes a mutation process with certain mutation probabilitypm. At last, there is a parameter σcrowd in selection operator that could indicate the distance between a solution and its surrounding solution in objective space. Further information about the process of MODE algorithms could be seen in [23].


## 2.3. MOEA/D Algorithm

The main idea of MOEA/D algorithm lies in decomposing MOP into multiple scalar sub-problems and optimizing them at the same time [24]. Every sub-problem utilizes the information of its neighboring sub-problems, which reduces its computational complexity of each generation. Scalarizing functions which are provided with uniformly distributed weight vectors is the fitness evaluation condition of MOEA/D [25]. In the computational experiments of this paper, the Tchebycheff function is employed to decompose. Further information about the process of MOEA/D algorithms could be seen in [24].


## 3. Multi-Swarm Multi-Objective Optimizer Based onp-Optimality Criteria

### 3.1. Introduction ofP-Optimality Criteria

P-optimality criteria named by Emiliano (2014) are a new kind of optimality criteria to solve MOP. These criteria have the ability of determining the most feasible solutions among the ones located in the same non-domination rank.As described in [3], the criteria are defined as follows. Let K′≠∞ be finite set of feasible solutions. A vector a∗∈K′ is an optimal solution if a∈K′ and Pi(a∗)≥Pi(a), where Pi(a) is the probability that a is better than other solutions from K′ in terms of an objective function (fi, i=1,2,…,k).Pi(a) is calculated as follows:(1)Pia=1-Sia+1K′whereSi(a) is the rank ofa according to the objective functionfi through the quick-sort algorithm.Thence, the aim is searching a vectora∗⊆K′ that maximizes fa=k-1∙∑i=1kPi(a) or equivalentlyfa=∑i=1kPi(a).Ifa∗ maximizes f(a), the following formula should be maximized as well.(2)fa=k-∑i=1kPia=∑i=1k1-Pia=∑i=1k1-PiaInspired by thep-norm, the following function is considered:(3)p-functiona=∑i=1k1-Piap1/pwhere p>0, andpi(a) is calculated as (1). Then,p-optimality criteria could be defined as follows.LetK′≠∞ be the finite set of feasible solutions. A vector a∗∈K′ is optimal if a∈K′ andp-function(a∗) ≤p-function(a).P-optimality criteria have the ability of determining the most feasible solutions among ones located in the same non-domination rank. Section 3.3 will illustrate the evolution of MOEAs combiningp-optimality criteria.


### 3.2. Introduction of Multi-Swarm Strategy

The above three algorithms introduced in Section2 use the analogy of single swarm. However, only with combiningp-optimality criteria, the above multi-objective algorithms may suffer from excessive loss of diversity. Indeed, multi-swarm strategy could restrict such rapid convergence and increase diversity effectively due to cooperation and exchange between swarms.The linchpin of multi-swarm strategy lies in two aspects: on the one hand,K individuals with superior performance are selected from sub-swarms to compose a distribution list. Meanwhile, each sub-swarm prepares a replacement list comprised ofK individuals with inferior performance.K is predefined, and the individuals in distribution list and replacement list are selected by non-domination rank and crowding distance.The whole swarm is evolved in the form of predetermined number sub-swarms, each of them performs the same multi-objective algorithm at early stage. After some predefined generations of evolution, the individuals in distribution list from one sub-swarm are sent to the adjacent sub-swarm to replace individuals in its replacement list. Taking this order, the last sub-swarm performs information transformation with first sub-swarm, as shown in Figure1.Figure 1
Hierarchical interaction topology.Figure1 illustrates hierarchical interaction topologies among individuals. Ring is adopted in swarm level, and ring as well as star topologies are simultaneously employed in individual level. Different sub-swarms are arranged in a unidirectional ring. In other words, every sub-swarm could accept individuals in its distribution list from adjacent sub-swarm to replace individuals of its replacement list. The simplicity of the structure makes this interaction topology very fast in interaction [26].


### 3.3. Description of the Proposedp-MSMOEAs

Multi-swarm non-dominated sorting genetic algorithm II based onp-optimality criteria, namedp-MSNSGA-II, is demonstrated in this part.


#### 3.3.1.p-MSNSGA-II Algorithm

At the beginning ofp-MSNSGA-II algorithm, the maximum number of cycles (MCN), the number of swarm maximum cycles (SMCN), and the number of exchange of swarms (EN) are set. The algorithm will loop a predefined number of times (MCN). As an example, the followingtth loop will be described.In initialization,m sub-swarms (S1,S2,…,Sm) ofN solutions are stochastically generated. Each solution hasn real-valued variables. The control parameter ofp-optimality criteriapj ∈ [p-lb,p-ub] (1 ≤j ≤m) is determined for each sub-swarm Sjt (p-lb andp-ub are lower and upper limits ofpj, respectively). Valuesp1,p2, …,pm are isometric, andp1=p-lb andpm =p-ub. Solutions in sub-swarm are sorted and assigned a certain rank based on the non-domination sorting.For sub-swarmSjt, three operators, including selection, recombination, and mutation, are involved in the update process. The current sub-swarm Sjt is represented as parent swarm. Firstly, a group ofq solutions is stochastically selected from Sjt, and then the solution with the least value ofp-function is determined. After repetitions,N solutions are selected. Secondly, the recombination operator and mutation operator in NSGA-II are employed. Finally, the new sub-swarm Sjt′ of sizeN is regarded as the offspring swarm. The combination of parent and offspring swarms Rjt=Sjt∪Sjt′ is carried out. Consequently, non-domination sorting and crowding distance are utilized to create new parent swarm Sjt+1.After predefined generationsSMCN,K individuals with superior performance are selected from sub-swarms into distribution list, as introduced in Section 3.2. AndK individuals with inferior performance are selected from sub-swarms into replacement list. AfterEN times of information transformation among sub-swarms, all sub-swarms merge into one whole swarm. Thereafter, the whole swarm is no longer partitioned.


#### 3.3.2. Flowchart

The general steps ofp-MSNSGA-II algorithm are shown in Algorithm 1. And the flowchart is illustrated in Figure 2.Algorithm 1:Pseudocode for p-MSNSGA-II.
Main steps ofp-MSNSGA-II ( 1 ) Set the number of swarm (m), the size of each swarm (N), the lower and upper bound of p (p-lb, p-ub), the maximum number of cycles (MCN), the number of swarm maximum cycles (SMCN) and the number of exchange of swarms (EN) ( 2 )   repeat ( 3 ) Initialize the sub-swarms S1,S2,…,Sm, assign p1,p2,…,pm for each swarm ( 4 )  repeat ( 5 ) The initialized solutions are sorted based on non-domination ( 6 ) For each sub-swarm Sjt Selection operator based on p-optimality criteria are adopted to select superior individuals Recombination operator and mutation operator are used Offspring swarm is Sjt′ created Non-domination and crowding distance are used by Sjt and Sjt′ to create Sjt+1 If mod (I,SMCN)=0 Exchange among m sub-swarms End if ( 7 ) End for ( 8 )  Until cycle = SMCN∗EN ( 9 ) Sub-swarms are integrated together and sorted based on non-domination ( 10 )  Until cycle = MCNFigure 2
The flowchart of thep-MSNSGA-II algorithm.


#### 3.3.3. Otherp-MSMOEAs

MODE and MOEA/D based on multi-swarm strategy andp-optimality criteria (namedp-MSMODE algorithm andp-MSMOEA/D algorithm, respectively) are similar to the abovep-MSNSGA-II algorithm.Thep-MSMODE algorithm andp-MSMOEA/D algorithm also start by stochastically generatingm initial sub-swarms. Forp-MSMODE algorithm (orp-MSMOEA/D algorithm), before evolution operation, the superior individuals are chosen from each sub-swarm based onp-optimality criteria and then are placed in the intermediate mating pool (or temporary elite population) for subsequent evolution operation. After some predefined generations of evolution,K individuals with superior (inferior) performance are selected for information transformation same asp-MSNSGA-II algorithm. AfterEN times of information transformation among sub-swarms, all sub-swarms merge into one whole swarm.


## 3.1. Introduction ofP-Optimality Criteria

P-optimality criteria named by Emiliano (2014) are a new kind of optimality criteria to solve MOP. These criteria have the ability of determining the most feasible solutions among the ones located in the same non-domination rank.As described in [3], the criteria are defined as follows. Let K′≠∞ be finite set of feasible solutions. A vector a∗∈K′ is an optimal solution if a∈K′ and Pi(a∗)≥Pi(a), where Pi(a) is the probability that a is better than other solutions from K′ in terms of an objective function (fi, i=1,2,…,k).Pi(a) is calculated as follows:(1)Pia=1-Sia+1K′whereSi(a) is the rank ofa according to the objective functionfi through the quick-sort algorithm.Thence, the aim is searching a vectora∗⊆K′ that maximizes fa=k-1∙∑i=1kPi(a) or equivalentlyfa=∑i=1kPi(a).Ifa∗ maximizes f(a), the following formula should be maximized as well.(2)fa=k-∑i=1kPia=∑i=1k1-Pia=∑i=1k1-PiaInspired by thep-norm, the following function is considered:(3)p-functiona=∑i=1k1-Piap1/pwhere p>0, andpi(a) is calculated as (1). Then,p-optimality criteria could be defined as follows.LetK′≠∞ be the finite set of feasible solutions. A vector a∗∈K′ is optimal if a∈K′ andp-function(a∗) ≤p-function(a).P-optimality criteria have the ability of determining the most feasible solutions among ones located in the same non-domination rank. Section 3.3 will illustrate the evolution of MOEAs combiningp-optimality criteria.


## 3.2. Introduction of Multi-Swarm Strategy

The above three algorithms introduced in Section2 use the analogy of single swarm. However, only with combiningp-optimality criteria, the above multi-objective algorithms may suffer from excessive loss of diversity. Indeed, multi-swarm strategy could restrict such rapid convergence and increase diversity effectively due to cooperation and exchange between swarms.The linchpin of multi-swarm strategy lies in two aspects: on the one hand,K individuals with superior performance are selected from sub-swarms to compose a distribution list. Meanwhile, each sub-swarm prepares a replacement list comprised ofK individuals with inferior performance.K is predefined, and the individuals in distribution list and replacement list are selected by non-domination rank and crowding distance.The whole swarm is evolved in the form of predetermined number sub-swarms, each of them performs the same multi-objective algorithm at early stage. After some predefined generations of evolution, the individuals in distribution list from one sub-swarm are sent to the adjacent sub-swarm to replace individuals in its replacement list. Taking this order, the last sub-swarm performs information transformation with first sub-swarm, as shown in Figure1.Figure 1
Hierarchical interaction topology.Figure1 illustrates hierarchical interaction topologies among individuals. Ring is adopted in swarm level, and ring as well as star topologies are simultaneously employed in individual level. Different sub-swarms are arranged in a unidirectional ring. In other words, every sub-swarm could accept individuals in its distribution list from adjacent sub-swarm to replace individuals of its replacement list. The simplicity of the structure makes this interaction topology very fast in interaction [26].


## 3.3. Description of the Proposedp-MSMOEAs

Multi-swarm non-dominated sorting genetic algorithm II based onp-optimality criteria, namedp-MSNSGA-II, is demonstrated in this part.


### 3.3.1.p-MSNSGA-II Algorithm

At the beginning ofp-MSNSGA-II algorithm, the maximum number of cycles (MCN), the number of swarm maximum cycles (SMCN), and the number of exchange of swarms (EN) are set. The algorithm will loop a predefined number of times (MCN). As an example, the followingtth loop will be described.In initialization,m sub-swarms (S1,S2,…,Sm) ofN solutions are stochastically generated. Each solution hasn real-valued variables. The control parameter ofp-optimality criteriapj ∈ [p-lb,p-ub] (1 ≤j ≤m) is determined for each sub-swarm Sjt (p-lb andp-ub are lower and upper limits ofpj, respectively). Valuesp1,p2, …,pm are isometric, andp1=p-lb andpm =p-ub. Solutions in sub-swarm are sorted and assigned a certain rank based on the non-domination sorting.For sub-swarmSjt, three operators, including selection, recombination, and mutation, are involved in the update process. The current sub-swarm Sjt is represented as parent swarm. Firstly, a group ofq solutions is stochastically selected from Sjt, and then the solution with the least value ofp-function is determined. After repetitions,N solutions are selected. Secondly, the recombination operator and mutation operator in NSGA-II are employed. Finally, the new sub-swarm Sjt′ of sizeN is regarded as the offspring swarm. The combination of parent and offspring swarms Rjt=Sjt∪Sjt′ is carried out. Consequently, non-domination sorting and crowding distance are utilized to create new parent swarm Sjt+1.After predefined generationsSMCN,K individuals with superior performance are selected from sub-swarms into distribution list, as introduced in Section 3.2. AndK individuals with inferior performance are selected from sub-swarms into replacement list. AfterEN times of information transformation among sub-swarms, all sub-swarms merge into one whole swarm. Thereafter, the whole swarm is no longer partitioned.


### 3.3.2. Flowchart

The general steps ofp-MSNSGA-II algorithm are shown in Algorithm 1. And the flowchart is illustrated in Figure 2.Algorithm 1:Pseudocode for p-MSNSGA-II.
Main steps ofp-MSNSGA-II ( 1 ) Set the number of swarm (m), the size of each swarm (N), the lower and upper bound of p (p-lb, p-ub), the maximum number of cycles (MCN), the number of swarm maximum cycles (SMCN) and the number of exchange of swarms (EN) ( 2 )   repeat ( 3 ) Initialize the sub-swarms S1,S2,…,Sm, assign p1,p2,…,pm for each swarm ( 4 )  repeat ( 5 ) The initialized solutions are sorted based on non-domination ( 6 ) For each sub-swarm Sjt Selection operator based on p-optimality criteria are adopted to select superior individuals Recombination operator and mutation operator are used Offspring swarm is Sjt′ created Non-domination and crowding distance are used by Sjt and Sjt′ to create Sjt+1 If mod (I,SMCN)=0 Exchange among m sub-swarms End if ( 7 ) End for ( 8 )  Until cycle = SMCN∗EN ( 9 ) Sub-swarms are integrated together and sorted based on non-domination ( 10 )  Until cycle = MCNFigure 2
The flowchart of thep-MSNSGA-II algorithm.


### 3.3.3. Otherp-MSMOEAs

MODE and MOEA/D based on multi-swarm strategy andp-optimality criteria (namedp-MSMODE algorithm andp-MSMOEA/D algorithm, respectively) are similar to the abovep-MSNSGA-II algorithm.Thep-MSMODE algorithm andp-MSMOEA/D algorithm also start by stochastically generatingm initial sub-swarms. Forp-MSMODE algorithm (orp-MSMOEA/D algorithm), before evolution operation, the superior individuals are chosen from each sub-swarm based onp-optimality criteria and then are placed in the intermediate mating pool (or temporary elite population) for subsequent evolution operation. After some predefined generations of evolution,K individuals with superior (inferior) performance are selected for information transformation same asp-MSNSGA-II algorithm. AfterEN times of information transformation among sub-swarms, all sub-swarms merge into one whole swarm.


## 3.3.1.p-MSNSGA-II Algorithm

At the beginning ofp-MSNSGA-II algorithm, the maximum number of cycles (MCN), the number of swarm maximum cycles (SMCN), and the number of exchange of swarms (EN) are set. The algorithm will loop a predefined number of times (MCN). As an example, the followingtth loop will be described.In initialization,m sub-swarms (S1,S2,…,Sm) ofN solutions are stochastically generated. Each solution hasn real-valued variables. The control parameter ofp-optimality criteriapj ∈ [p-lb,p-ub] (1 ≤j ≤m) is determined for each sub-swarm Sjt (p-lb andp-ub are lower and upper limits ofpj, respectively). Valuesp1,p2, …,pm are isometric, andp1=p-lb andpm =p-ub. Solutions in sub-swarm are sorted and assigned a certain rank based on the non-domination sorting.For sub-swarmSjt, three operators, including selection, recombination, and mutation, are involved in the update process. The current sub-swarm Sjt is represented as parent swarm. Firstly, a group ofq solutions is stochastically selected from Sjt, and then the solution with the least value ofp-function is determined. After repetitions,N solutions are selected. Secondly, the recombination operator and mutation operator in NSGA-II are employed. Finally, the new sub-swarm Sjt′ of sizeN is regarded as the offspring swarm. The combination of parent and offspring swarms Rjt=Sjt∪Sjt′ is carried out. Consequently, non-domination sorting and crowding distance are utilized to create new parent swarm Sjt+1.After predefined generationsSMCN,K individuals with superior performance are selected from sub-swarms into distribution list, as introduced in Section 3.2. AndK individuals with inferior performance are selected from sub-swarms into replacement list. AfterEN times of information transformation among sub-swarms, all sub-swarms merge into one whole swarm. Thereafter, the whole swarm is no longer partitioned.


## 3.3.2. Flowchart

The general steps ofp-MSNSGA-II algorithm are shown in Algorithm 1. And the flowchart is illustrated in Figure 2.Algorithm 1:Pseudocode for p-MSNSGA-II.
Main steps ofp-MSNSGA-II ( 1 ) Set the number of swarm (m), the size of each swarm (N), the lower and upper bound of p (p-lb, p-ub), the maximum number of cycles (MCN), the number of swarm maximum cycles (SMCN) and the number of exchange of swarms (EN) ( 2 )   repeat ( 3 ) Initialize the sub-swarms S1,S2,…,Sm, assign p1,p2,…,pm for each swarm ( 4 )  repeat ( 5 ) The initialized solutions are sorted based on non-domination ( 6 ) For each sub-swarm Sjt Selection operator based on p-optimality criteria are adopted to select superior individuals Recombination operator and mutation operator are used Offspring swarm is Sjt′ created Non-domination and crowding distance are used by Sjt and Sjt′ to create Sjt+1 If mod (I,SMCN)=0 Exchange among m sub-swarms End if ( 7 ) End for ( 8 )  Until cycle = SMCN∗EN ( 9 ) Sub-swarms are integrated together and sorted based on non-domination ( 10 )  Until cycle = MCNFigure 2
The flowchart of thep-MSNSGA-II algorithm.


## 3.3.3. Otherp-MSMOEAs

MODE and MOEA/D based on multi-swarm strategy andp-optimality criteria (namedp-MSMODE algorithm andp-MSMOEA/D algorithm, respectively) are similar to the abovep-MSNSGA-II algorithm.Thep-MSMODE algorithm andp-MSMOEA/D algorithm also start by stochastically generatingm initial sub-swarms. Forp-MSMODE algorithm (orp-MSMOEA/D algorithm), before evolution operation, the superior individuals are chosen from each sub-swarm based onp-optimality criteria and then are placed in the intermediate mating pool (or temporary elite population) for subsequent evolution operation. After some predefined generations of evolution,K individuals with superior (inferior) performance are selected for information transformation same asp-MSNSGA-II algorithm. AfterEN times of information transformation among sub-swarms, all sub-swarms merge into one whole swarm.


## 4. Test and Results

To completely evaluate the performance of abovep-MSMOEAs without bias against some certain selected problems, four two-objective as well as two three-objective benchmark functions are utilized. Formulas of those functions could be seen in Appendix.


### 4.1. Evaluation Method

For the purpose of facilitating the quantitative assessment of the performance of proposed algorithms, two performance metrics should be considered: convergence metricγ and diversity metric Δ. More detailed information about them could be seen in [27, 28].


### 4.2. Experimental Setting

Experiments have been executed withp-MSNSGA-II,p-MSMODE,p-MSMOEA/D, NSGA-II, MODE, and MOEA/D. To contrast different algorithms with a fair time metric during experiment, population size is set as 200, and the number of function evaluations (FEs) is set as 40000. Each algorithm runs 10 times, using a PC Intel Core i5-7200U, 2.50 GHz CPU with 8GB of RAM.For NSGA-II, MODE, and MOEA/D, parameter settings are the same as original algorithms described in [22–24], respectively.Forp-MSNSGA-II,p-MSMODE, andp-MSMOEA/D, the whole swarm with 200 individuals is equally partitioned into 4 sub-swarms. And the lower and upper bounds of thep-optimality criteria’s interval are set as 0.5 and 2.0, respectively (according to [3]). The number of swarm maximum cycles (SMCN) and the number of exchange of swarms (EN) are set as 20 and 10, respectively. The algorithm will loop 5 times (MCN). And the other setting is the same as original algorithms.


### 4.3. Results and Discussion

Test results ofp-MSMOEAs and MOEAs on six benchmark functions, including maximum, minimum, average, and standard deviation of the convergence metric (γ) and the diversity metric (Δ) values, are listed in this part. Besides, CPU time is employed to measure the time complexity of the algorithms. To further demonstrate performance ofp-MSMOEAs, results of another multi-swarm algorithm, multi-hive multi-objective artificial bee colony (M2OABC) algorithm proposed in [26], on several of above benchmark functions, are also listed in the following. M2OABC algorithm has been proved to be effective and robust with combining multi-hive strategy.


#### 4.3.1. Two-Objective Functions

Tables1–4 and Figures 3–6 show the optimization results of these algorithms for two-objective functions. In Figures 3–6, the solid lines represent true Pareto front (PF), while the star spots stand for found non-dominated solutions.Table 1
Comparison of performance on SCH2.
SCH2 p-MSNSGA-II p-MSMODE p-MSMOEA/D NSGA-II MODE MOEA/D Converge Metric Max 2.04e-05 1.72e-04 8.03e-04 7.80e-04 3.83e-03 9.96e-03 Min 3.93e-06 1.20e-06 7.05e-04 3.78e-04 3.35e-03 7.29e-03 Avg 5.01e-06 8.73e-05 7.59e-04 4.01e-04 3.54e-03 8.19e-03 Std 2.10e-06 1.61e-05 3.99e-05 1.10e-05 1.98e-04 5.29e-04  Diversity Metric Max 8.62e-02 6.48e-02 1.13e-01 1.16e-01 1.19e-01 1.31e-01 Min 3.06e-02 4.93e-02 7.01e-02 4.09e-02 6.62e-02 7.25e-02 Avg 6.81e-02 5.19e-02 7.73e-02 6.48e-02 7.21e-02 8.59e-02 Std 5.59e-03 4.95e-03 2.57e-03 1.27e-03 5.50e-03 7.60e-03Table 2
Comparison of performance on ZDT3.
ZDT3 p-MSNSGA-II p-MSMODE p-MSMOEA/D NSGA-II MODE MOEA/D M 2 OABC Converge Metric Max 1.28e-04 2.86e-04 4.34e-04 2.26e-02 2.05e-02 3.05e-02 4.64e-03 Min 5.37e-05 2.04e-05 3.96e-04 1.09e-02 5.51e-03 7.40e-03 3.55e-03 Avg 8.59e-05 2.51e-05 4.15e-04 2.03e-02 1.76e-02 8.62e-03 4.07e-03 Std 3.00e-05 4.92e-04 1.59e-05 5.71e-02 7.85e-03 6.30e-02 3.23e-04  Diversity Metric Max 6.40e-02 6.52e-02 1.28e-01 6.59e-02 9.82e-02 1.32e-01 6.71e-02 Min 3.79e-02 5.65e-02 5.65e-02 4.63e-02 5.85e-02 6.55e-02 5.57e-02 Avg 4.62e-02 6.24e-02 7.69e-02 5.21e-02 7.73e-02 8.45e-02 6.27e-02 Std 2.49e-03 2.31e-03 1.42e-02 3.85e-03 7.42e-03 1.80e-02 1.39e-02Table 3
Comparison of performance on ZDT4.
ZDT4 p-MSNSGA-II p-MSMODE p-MSMOEA/D NSGA-II MODE MOEA/D Converge Metric Max 1.14e-02 7.51e-03 9.41e-03 1.21e-01 7.68e-02 9.93e-02 Min 1.01e-02 1.12e-03 1.62e-03 2.40e-02 1.97e-03 3.01e-02 Avg 1.04e-02 2.61e-03 4.77e-03 3.69e-02 1.64e-02 4.78e-02 Std 5.78e-04 1.14e-02 1.56e-02 1.44e-02 2.05e-02 2.41e-02  Diversity Metric Max 7.91e-02 9.87e-02 1.19e-01 1.09e-01 1.68e-01 1.18e-01 Min 6.13e-02 4.90e-02 7.05e-02 6.32e-02 9.05e-02 8.24e-02 Avg 6.63e-02 7.02e-02 9.81e-02 8.93e-02 9.08e-02 1.00e-01 Std 7.36e-03 8.87e-03 6.78e-03 1.12e-02 1.64e-01 6.89e-03Table 4
Comparison of performance on ZDT6.
ZDT6 p-MSNSGA-II p-MSMODE p-MSMOEA/D NSGA-II MODE MOEA/D M 2 OABC Converge Metric Max 3.87e-05 6.45e-04 6.74e-04 7.72e-02 7.37e-02 7.01e-02 1.76e-03 Min 2.39e-05 2.71e-05 3.21e-04 5.31e-02 2.72e-03 5.09e-02 9.64e-05 Avg 2.95e-05 2.89e-04 4.30e-04 6.98e-02 3.71e-02 5.48e-02 5.03e-04 Std 6.08e-06 1.10e-05 7.63e-05 6.52e-03 2.17e-02 6.45e-03 7.40e-01  Diversity Metric Max 7.51e-02 7.32e-02 1.01e-01 1.12e+00 1.24e+00 1.16e+00 1.18e-01 Min 5.87e-02 4.08e-02 6.76e-02 6.12e-01 9.31e-01 6.19e-01 6.08e-02 Avg 6.37e-02 5.68e-02 8.14e-02 8.35e-01 9.99e-01 8.96e-01 7.96e-02 Std 1.26e-03 4.03e-03 2.06e-03 1.08e-01 1.91e-02 1.87e-01 1.93e-02Figure 3
PF obtained byp-MSMOEAs and MOEAs on SCH2.
(a)
PF obtained byp-MSNSGA-II on SCH2 (b)
PF obtained by NSGA-II on SCH2 (c)
PF obtained byp-MSMODE on SCH2 (d)
PF obtained by MODE on SCH2 (e)
PF obtained byp- MSMOEA/D on SCH2 (f)
PF obtained by MOEA/D on SCH2Figure 4
PF obtained byp-MSMOEAs, MOEAs, and M2OABC on ZDT3.
(a)
PF obtained byp-MSNSGA-II on ZDT3 (b)
PF obtained by NSGA-II on ZDT3 (c)
PF obtained byp-MSMODE on ZDT3 (d)
PF obtained by MODE on ZDT3 (e)
PF obtained byp-MSMOEA/D on ZDT3 (f)
PF obtained by MOEA/D on ZDT3  (g)
PF obtained by M2OABC on ZDT3Figure 5
PF obtained byp-MSMOEAs and MOEAs on ZDT4.
(a)
PF obtained byp-MSNSGA-II on ZDT4 (b)
PF obtained by NSGA-II on ZDT4 (c)
PF obtained byp-MSMODE on ZDT4 (d)
PF obtained by MODE on ZDT4 (e)
PF obtained byp-MSMOEA/D on ZDT4 (f)
PF obtained by MOEA/D on ZDT4Figure 6
PF obtained byp-MSMOEAs, MOEAs, and M2OABC on ZDT6.
(a)
PF obtained byp-MSNSGA-II on ZDT6 (b)
PF obtained by NSGA-II on ZDT6 (c)
PF obtained byp-MSMODE on ZDT6 (d)
PF obtained by MODE on ZDT6 (e)
PF obtained byp-MSMOEA/D on ZDT6 (f)
PF obtained by MOEA/D on ZDT6  (g)
PF obtained by M2OABC on ZDT6On SCH2 function, it can be noticed thatp-MSMOEAs have superior performance to MOEAs in terms of γ metric and Δ metric after 40000 FEs from Table 1. Figure 3 shows thatp-MSMOEAs have great potential to discover a well-distributed as well as diverse solution set for SCH2 function. Yet MOEA/D only finds a sparse distribution, although it can basically archive true PF for SCH2.On ZDT3, ZDT4, and ZDT6 functions, Tables2–4 show that performance ofp-MSMOEAs in both γ metric and Δ metric is better than that of MOEAs and lightly better than that of M2OABC. Figures 4–6 show that MOEAs produce poor results on these test functions and they are almost impossible to achieve true PF, whilep-MSMOEAs have great potential to approach true PF.


#### 4.3.2. Three-Objective Functions

Figure7 shows the true PF for two three-objective functions. Tables 5 and 6 and Figures 8 and 9 show the optimization results ofp-MSMOEAs, MOEAs, and M2OABC algorithms on DTLZ2 and DTLZ3.Table 5
Comparison of performance on DTLZ2.
DTLZ2 p-MSNSGA-II p-MSMODE p-MSMOEA/D NSGA-II MODE MOEA/D M 2 OABC Converge Metric Max 3.02e-03 7.13e-04 5.67e-03 6.51e-02 7.97e-02 6.27e-02 3.10e-03 Min 2.88e-03 2.30e-05 6.80e-05 3.05e-02 7.50e-04 1.58e-03 2.56e-03 Avg 2.96e-03 1.96e-04 4.06e-04 4.15e-02 1.79e-03 5.51e-03 2.85e-03 Std 6.30e-04 8.48e-05 5.87e-04 8.10e-03 2.11e-02 1.01e-02 1.54e-04  Diversity Metric Max 4.60e-02 4.64e-02 7.61e-02 6.73e-01 4.92e-01 9.98e-01 4.66e-02 Min 3.37e-02 3.75e-02 2.16e-02 3.85e-01 4.23e-01 2.88e-01 3.95e-02 Avg 4.46e-02 4.21e-02 4.76e-02 5.07e-01 4.77e-01 5.39e-01 4.30e-02 Std 1.97e-03 1.88e-03 7.12e-03 5.32e-02 7.81e-02 1.29e-01 2.27e-03Table 6
Comparison of performance on DTLZ3.
DTLZ3 p-MSNSGA-II p-MSMODE p-MSMOEA/D NSGA-II MODE MOEA/D Converge Metric Max 4.10e-02 1.05e-01 8.26e-01 6.56e-01 1.08e+00 1.02e+00 Min 3.48e-02 1.86e-03 6.73e-02 5.08e-02 1.53e-01 1.78e-01 Avg 3.92e-02 4.70e-02 1.84e-01 1.21e-01 3.72e-01 3.60e-01 Std 2.55e-01 1.23e-01 1.24+00 8.65e-01 2.38e+00 2.6+00  Diversity Metric Max 6.97e-02 4.67e-02 1.21e-01 1.24e-01 4.81e-02 1.27e-01 Min 4.21e-02 3.49e-02 2.70e-02 4.54e-02 3.61e-02 3.45e-02 Avg 5.87e-02 3.75e-02 8.00e-02 9.42e-02 4.16e-02 8.02e-02 Std 7.35e-03 2.62e-03 1.52e-02 2.07e-02 2.81e-03 2.06e-02Figure 7
True PF on DTLZ2 and DTLZ3.
(a)
True PF on DTLZ2 (b)
True PF on DTLZ3Figure 8
PF obtained byp-MSMOEAs and MOEAs on DTLZ2.
(a)
PF obtained byp-MSNSGA-II on DTLZ2 (b)
PF obtained by NSGA-II on DTLZ2 (c)
PF obtained byp-MSMODE on DTLZ2 (d)
PF obtained by MODE on DTLZ2 (e)
PF obtained byp-MSMOEA/D on DTLZ2 (f)
PF obtained by MOEA/D on DTLZ2Figure 9
PF obtained byp-MSMOEAs and MOEAs on DTLZ3.
(a)
PF obtained byp-MSNSGA-II on DTLZ3 (b)
PF obtained by NSGA-II on DTLZ3 (c)
PF obtained byp-MSMODE on DTLZ3 (d)
PF obtained by MODE on DTLZ3 (e)
PF obtained byp-MSMOEA/D on DTLZ3 (f)
PF obtained by MOEA/D on DTLZ3On DTLZ2 function, when given 10000 FEs for seven algorithms, performance ofp-MSMOEAs is better than that of MOEAs and is comparable to that of M2OABC as shown in Table 5. Figure 8 shows thatp-MSMOEAs have great potential to obtain a superior PF for DTLZ2, especiallyp-MSNSGA-II andp-MSMODE. Moreover, the performance of MOEAs in Δ metric is a little worse than that ofp-MSMOEAs.On DTLZ3 function, it could be observed from Table6 that performance ofp-MSMOEAs algorithms in both γ metric and Δ metric has considerable competitiveness over this problem. From Figure 9, it can be seen that the fronts obtained fromp-MSMOEAs and MOEAs are found to have a basically outstanding result in terms of convergence, while they do not perform satisfactorily in terms of diversity. However,p-MSMOEAs are better than MOEAs for the problem.


#### 4.3.3. Time Complexity Analysis

In order to demonstrate the difference in the time complexity of six algorithms, Figure10 plots the average CPU time over 10 runs.Figure 10
CPU time of MOEAs andp-MSMOEAs on six benchmark functions.Figure10 shows the CPU time for six algorithms and gives the results on six benchmark functions (SCH2, ZDT3, ZDT4, ZDT6, DTLZ2, and DTLZ3). The results show that the time complexity of the original three MOEAs is basically the same, and NSGA-II and MOEA/D are slightly better than MODE. The time complexity ofp-MSMOEAs is slightly higher than that of MOEAs, but within acceptable limits. In addition, for two-objective functions (SCH2, ZDT3, ZDT4, and ZDT6), the time complexity ofp-MSMOEAs is about 1.32-1.64 times that of the original MOEAs. However, as the number of objectives increases, this difference is more pronounced. For example, for three-objective functions (DTLZ2 and DTLZ3), the time complexity ofp-MSMOEAs is about 1.38-1.92 times that of MOEAs. It is worth mentioning that althoughp-MSMOEAs consume more CPU time, there is a significant improvement in performance compared to MOEAs.


## 4.1. Evaluation Method

For the purpose of facilitating the quantitative assessment of the performance of proposed algorithms, two performance metrics should be considered: convergence metricγ and diversity metric Δ. More detailed information about them could be seen in [27, 28].


## 4.2. Experimental Setting

Experiments have been executed withp-MSNSGA-II,p-MSMODE,p-MSMOEA/D, NSGA-II, MODE, and MOEA/D. To contrast different algorithms with a fair time metric during experiment, population size is set as 200, and the number of function evaluations (FEs) is set as 40000. Each algorithm runs 10 times, using a PC Intel Core i5-7200U, 2.50 GHz CPU with 8GB of RAM.For NSGA-II, MODE, and MOEA/D, parameter settings are the same as original algorithms described in [22–24], respectively.Forp-MSNSGA-II,p-MSMODE, andp-MSMOEA/D, the whole swarm with 200 individuals is equally partitioned into 4 sub-swarms. And the lower and upper bounds of thep-optimality criteria’s interval are set as 0.5 and 2.0, respectively (according to [3]). The number of swarm maximum cycles (SMCN) and the number of exchange of swarms (EN) are set as 20 and 10, respectively. The algorithm will loop 5 times (MCN). And the other setting is the same as original algorithms.


## 4.3. Results and Discussion

Test results ofp-MSMOEAs and MOEAs on six benchmark functions, including maximum, minimum, average, and standard deviation of the convergence metric (γ) and the diversity metric (Δ) values, are listed in this part. Besides, CPU time is employed to measure the time complexity of the algorithms. To further demonstrate performance ofp-MSMOEAs, results of another multi-swarm algorithm, multi-hive multi-objective artificial bee colony (M2OABC) algorithm proposed in [26], on several of above benchmark functions, are also listed in the following. M2OABC algorithm has been proved to be effective and robust with combining multi-hive strategy.


### 4.3.1. Two-Objective Functions

Tables1–4 and Figures 3–6 show the optimization results of these algorithms for two-objective functions. In Figures 3–6, the solid lines represent true Pareto front (PF), while the star spots stand for found non-dominated solutions.Table 1
Comparison of performance on SCH2.
SCH2 p-MSNSGA-II p-MSMODE p-MSMOEA/D NSGA-II MODE MOEA/D Converge Metric Max 2.04e-05 1.72e-04 8.03e-04 7.80e-04 3.83e-03 9.96e-03 Min 3.93e-06 1.20e-06 7.05e-04 3.78e-04 3.35e-03 7.29e-03 Avg 5.01e-06 8.73e-05 7.59e-04 4.01e-04 3.54e-03 8.19e-03 Std 2.10e-06 1.61e-05 3.99e-05 1.10e-05 1.98e-04 5.29e-04  Diversity Metric Max 8.62e-02 6.48e-02 1.13e-01 1.16e-01 1.19e-01 1.31e-01 Min 3.06e-02 4.93e-02 7.01e-02 4.09e-02 6.62e-02 7.25e-02 Avg 6.81e-02 5.19e-02 7.73e-02 6.48e-02 7.21e-02 8.59e-02 Std 5.59e-03 4.95e-03 2.57e-03 1.27e-03 5.50e-03 7.60e-03Table 2
Comparison of performance on ZDT3.
ZDT3 p-MSNSGA-II p-MSMODE p-MSMOEA/D NSGA-II MODE MOEA/D M 2 OABC Converge Metric Max 1.28e-04 2.86e-04 4.34e-04 2.26e-02 2.05e-02 3.05e-02 4.64e-03 Min 5.37e-05 2.04e-05 3.96e-04 1.09e-02 5.51e-03 7.40e-03 3.55e-03 Avg 8.59e-05 2.51e-05 4.15e-04 2.03e-02 1.76e-02 8.62e-03 4.07e-03 Std 3.00e-05 4.92e-04 1.59e-05 5.71e-02 7.85e-03 6.30e-02 3.23e-04  Diversity Metric Max 6.40e-02 6.52e-02 1.28e-01 6.59e-02 9.82e-02 1.32e-01 6.71e-02 Min 3.79e-02 5.65e-02 5.65e-02 4.63e-02 5.85e-02 6.55e-02 5.57e-02 Avg 4.62e-02 6.24e-02 7.69e-02 5.21e-02 7.73e-02 8.45e-02 6.27e-02 Std 2.49e-03 2.31e-03 1.42e-02 3.85e-03 7.42e-03 1.80e-02 1.39e-02Table 3
Comparison of performance on ZDT4.
ZDT4 p-MSNSGA-II p-MSMODE p-MSMOEA/D NSGA-II MODE MOEA/D Converge Metric Max 1.14e-02 7.51e-03 9.41e-03 1.21e-01 7.68e-02 9.93e-02 Min 1.01e-02 1.12e-03 1.62e-03 2.40e-02 1.97e-03 3.01e-02 Avg 1.04e-02 2.61e-03 4.77e-03 3.69e-02 1.64e-02 4.78e-02 Std 5.78e-04 1.14e-02 1.56e-02 1.44e-02 2.05e-02 2.41e-02  Diversity Metric Max 7.91e-02 9.87e-02 1.19e-01 1.09e-01 1.68e-01 1.18e-01 Min 6.13e-02 4.90e-02 7.05e-02 6.32e-02 9.05e-02 8.24e-02 Avg 6.63e-02 7.02e-02 9.81e-02 8.93e-02 9.08e-02 1.00e-01 Std 7.36e-03 8.87e-03 6.78e-03 1.12e-02 1.64e-01 6.89e-03Table 4
Comparison of performance on ZDT6.
ZDT6 p-MSNSGA-II p-MSMODE p-MSMOEA/D NSGA-II MODE MOEA/D M 2 OABC Converge Metric Max 3.87e-05 6.45e-04 6.74e-04 7.72e-02 7.37e-02 7.01e-02 1.76e-03 Min 2.39e-05 2.71e-05 3.21e-04 5.31e-02 2.72e-03 5.09e-02 9.64e-05 Avg 2.95e-05 2.89e-04 4.30e-04 6.98e-02 3.71e-02 5.48e-02 5.03e-04 Std 6.08e-06 1.10e-05 7.63e-05 6.52e-03 2.17e-02 6.45e-03 7.40e-01  Diversity Metric Max 7.51e-02 7.32e-02 1.01e-01 1.12e+00 1.24e+00 1.16e+00 1.18e-01 Min 5.87e-02 4.08e-02 6.76e-02 6.12e-01 9.31e-01 6.19e-01 6.08e-02 Avg 6.37e-02 5.68e-02 8.14e-02 8.35e-01 9.99e-01 8.96e-01 7.96e-02 Std 1.26e-03 4.03e-03 2.06e-03 1.08e-01 1.91e-02 1.87e-01 1.93e-02Figure 3
PF obtained byp-MSMOEAs and MOEAs on SCH2.
(a)
PF obtained byp-MSNSGA-II on SCH2 (b)
PF obtained by NSGA-II on SCH2 (c)
PF obtained byp-MSMODE on SCH2 (d)
PF obtained by MODE on SCH2 (e)
PF obtained byp- MSMOEA/D on SCH2 (f)
PF obtained by MOEA/D on SCH2Figure 4
PF obtained byp-MSMOEAs, MOEAs, and M2OABC on ZDT3.
(a)
PF obtained byp-MSNSGA-II on ZDT3 (b)
PF obtained by NSGA-II on ZDT3 (c)
PF obtained byp-MSMODE on ZDT3 (d)
PF obtained by MODE on ZDT3 (e)
PF obtained byp-MSMOEA/D on ZDT3 (f)
PF obtained by MOEA/D on ZDT3  (g)
PF obtained by M2OABC on ZDT3Figure 5
PF obtained byp-MSMOEAs and MOEAs on ZDT4.
(a)
PF obtained byp-MSNSGA-II on ZDT4 (b)
PF obtained by NSGA-II on ZDT4 (c)
PF obtained byp-MSMODE on ZDT4 (d)
PF obtained by MODE on ZDT4 (e)
PF obtained byp-MSMOEA/D on ZDT4 (f)
PF obtained by MOEA/D on ZDT4Figure 6
PF obtained byp-MSMOEAs, MOEAs, and M2OABC on ZDT6.
(a)
PF obtained byp-MSNSGA-II on ZDT6 (b)
PF obtained by NSGA-II on ZDT6 (c)
PF obtained byp-MSMODE on ZDT6 (d)
PF obtained by MODE on ZDT6 (e)
PF obtained byp-MSMOEA/D on ZDT6 (f)
PF obtained by MOEA/D on ZDT6  (g)
PF obtained by M2OABC on ZDT6On SCH2 function, it can be noticed thatp-MSMOEAs have superior performance to MOEAs in terms of γ metric and Δ metric after 40000 FEs from Table 1. Figure 3 shows thatp-MSMOEAs have great potential to discover a well-distributed as well as diverse solution set for SCH2 function. Yet MOEA/D only finds a sparse distribution, although it can basically archive true PF for SCH2.On ZDT3, ZDT4, and ZDT6 functions, Tables2–4 show that performance ofp-MSMOEAs in both γ metric and Δ metric is better than that of MOEAs and lightly better than that of M2OABC. Figures 4–6 show that MOEAs produce poor results on these test functions and they are almost impossible to achieve true PF, whilep-MSMOEAs have great potential to approach true PF.


### 4.3.2. Three-Objective Functions

Figure7 shows the true PF for two three-objective functions. Tables 5 and 6 and Figures 8 and 9 show the optimization results ofp-MSMOEAs, MOEAs, and M2OABC algorithms on DTLZ2 and DTLZ3.Table 5
Comparison of performance on DTLZ2.
DTLZ2 p-MSNSGA-II p-MSMODE p-MSMOEA/D NSGA-II MODE MOEA/D M 2 OABC Converge Metric Max 3.02e-03 7.13e-04 5.67e-03 6.51e-02 7.97e-02 6.27e-02 3.10e-03 Min 2.88e-03 2.30e-05 6.80e-05 3.05e-02 7.50e-04 1.58e-03 2.56e-03 Avg 2.96e-03 1.96e-04 4.06e-04 4.15e-02 1.79e-03 5.51e-03 2.85e-03 Std 6.30e-04 8.48e-05 5.87e-04 8.10e-03 2.11e-02 1.01e-02 1.54e-04  Diversity Metric Max 4.60e-02 4.64e-02 7.61e-02 6.73e-01 4.92e-01 9.98e-01 4.66e-02 Min 3.37e-02 3.75e-02 2.16e-02 3.85e-01 4.23e-01 2.88e-01 3.95e-02 Avg 4.46e-02 4.21e-02 4.76e-02 5.07e-01 4.77e-01 5.39e-01 4.30e-02 Std 1.97e-03 1.88e-03 7.12e-03 5.32e-02 7.81e-02 1.29e-01 2.27e-03Table 6
Comparison of performance on DTLZ3.
DTLZ3 p-MSNSGA-II p-MSMODE p-MSMOEA/D NSGA-II MODE MOEA/D Converge Metric Max 4.10e-02 1.05e-01 8.26e-01 6.56e-01 1.08e+00 1.02e+00 Min 3.48e-02 1.86e-03 6.73e-02 5.08e-02 1.53e-01 1.78e-01 Avg 3.92e-02 4.70e-02 1.84e-01 1.21e-01 3.72e-01 3.60e-01 Std 2.55e-01 1.23e-01 1.24+00 8.65e-01 2.38e+00 2.6+00  Diversity Metric Max 6.97e-02 4.67e-02 1.21e-01 1.24e-01 4.81e-02 1.27e-01 Min 4.21e-02 3.49e-02 2.70e-02 4.54e-02 3.61e-02 3.45e-02 Avg 5.87e-02 3.75e-02 8.00e-02 9.42e-02 4.16e-02 8.02e-02 Std 7.35e-03 2.62e-03 1.52e-02 2.07e-02 2.81e-03 2.06e-02Figure 7
True PF on DTLZ2 and DTLZ3.
(a)
True PF on DTLZ2 (b)
True PF on DTLZ3Figure 8
PF obtained byp-MSMOEAs and MOEAs on DTLZ2.
(a)
PF obtained byp-MSNSGA-II on DTLZ2 (b)
PF obtained by NSGA-II on DTLZ2 (c)
PF obtained byp-MSMODE on DTLZ2 (d)
PF obtained by MODE on DTLZ2 (e)
PF obtained byp-MSMOEA/D on DTLZ2 (f)
PF obtained by MOEA/D on DTLZ2Figure 9
PF obtained byp-MSMOEAs and MOEAs on DTLZ3.
(a)
PF obtained byp-MSNSGA-II on DTLZ3 (b)
PF obtained by NSGA-II on DTLZ3 (c)
PF obtained byp-MSMODE on DTLZ3 (d)
PF obtained by MODE on DTLZ3 (e)
PF obtained byp-MSMOEA/D on DTLZ3 (f)
PF obtained by MOEA/D on DTLZ3On DTLZ2 function, when given 10000 FEs for seven algorithms, performance ofp-MSMOEAs is better than that of MOEAs and is comparable to that of M2OABC as shown in Table 5. Figure 8 shows thatp-MSMOEAs have great potential to obtain a superior PF for DTLZ2, especiallyp-MSNSGA-II andp-MSMODE. Moreover, the performance of MOEAs in Δ metric is a little worse than that ofp-MSMOEAs.On DTLZ3 function, it could be observed from Table6 that performance ofp-MSMOEAs algorithms in both γ metric and Δ metric has considerable competitiveness over this problem. From Figure 9, it can be seen that the fronts obtained fromp-MSMOEAs and MOEAs are found to have a basically outstanding result in terms of convergence, while they do not perform satisfactorily in terms of diversity. However,p-MSMOEAs are better than MOEAs for the problem.


### 4.3.3. Time Complexity Analysis

In order to demonstrate the difference in the time complexity of six algorithms, Figure10 plots the average CPU time over 10 runs.Figure 10
CPU time of MOEAs andp-MSMOEAs on six benchmark functions.Figure10 shows the CPU time for six algorithms and gives the results on six benchmark functions (SCH2, ZDT3, ZDT4, ZDT6, DTLZ2, and DTLZ3). The results show that the time complexity of the original three MOEAs is basically the same, and NSGA-II and MOEA/D are slightly better than MODE. The time complexity ofp-MSMOEAs is slightly higher than that of MOEAs, but within acceptable limits. In addition, for two-objective functions (SCH2, ZDT3, ZDT4, and ZDT6), the time complexity ofp-MSMOEAs is about 1.32-1.64 times that of the original MOEAs. However, as the number of objectives increases, this difference is more pronounced. For example, for three-objective functions (DTLZ2 and DTLZ3), the time complexity ofp-MSMOEAs is about 1.38-1.92 times that of MOEAs. It is worth mentioning that althoughp-MSMOEAs consume more CPU time, there is a significant improvement in performance compared to MOEAs.


## 4.3.1. Two-Objective Functions

Tables1–4 and Figures 3–6 show the optimization results of these algorithms for two-objective functions. In Figures 3–6, the solid lines represent true Pareto front (PF), while the star spots stand for found non-dominated solutions.Table 1
Comparison of performance on SCH2.
SCH2 p-MSNSGA-II p-MSMODE p-MSMOEA/D NSGA-II MODE MOEA/D Converge Metric Max 2.04e-05 1.72e-04 8.03e-04 7.80e-04 3.83e-03 9.96e-03 Min 3.93e-06 1.20e-06 7.05e-04 3.78e-04 3.35e-03 7.29e-03 Avg 5.01e-06 8.73e-05 7.59e-04 4.01e-04 3.54e-03 8.19e-03 Std 2.10e-06 1.61e-05 3.99e-05 1.10e-05 1.98e-04 5.29e-04  Diversity Metric Max 8.62e-02 6.48e-02 1.13e-01 1.16e-01 1.19e-01 1.31e-01 Min 3.06e-02 4.93e-02 7.01e-02 4.09e-02 6.62e-02 7.25e-02 Avg 6.81e-02 5.19e-02 7.73e-02 6.48e-02 7.21e-02 8.59e-02 Std 5.59e-03 4.95e-03 2.57e-03 1.27e-03 5.50e-03 7.60e-03Table 2
Comparison of performance on ZDT3.
ZDT3 p-MSNSGA-II p-MSMODE p-MSMOEA/D NSGA-II MODE MOEA/D M 2 OABC Converge Metric Max 1.28e-04 2.86e-04 4.34e-04 2.26e-02 2.05e-02 3.05e-02 4.64e-03 Min 5.37e-05 2.04e-05 3.96e-04 1.09e-02 5.51e-03 7.40e-03 3.55e-03 Avg 8.59e-05 2.51e-05 4.15e-04 2.03e-02 1.76e-02 8.62e-03 4.07e-03 Std 3.00e-05 4.92e-04 1.59e-05 5.71e-02 7.85e-03 6.30e-02 3.23e-04  Diversity Metric Max 6.40e-02 6.52e-02 1.28e-01 6.59e-02 9.82e-02 1.32e-01 6.71e-02 Min 3.79e-02 5.65e-02 5.65e-02 4.63e-02 5.85e-02 6.55e-02 5.57e-02 Avg 4.62e-02 6.24e-02 7.69e-02 5.21e-02 7.73e-02 8.45e-02 6.27e-02 Std 2.49e-03 2.31e-03 1.42e-02 3.85e-03 7.42e-03 1.80e-02 1.39e-02Table 3
Comparison of performance on ZDT4.
ZDT4 p-MSNSGA-II p-MSMODE p-MSMOEA/D NSGA-II MODE MOEA/D Converge Metric Max 1.14e-02 7.51e-03 9.41e-03 1.21e-01 7.68e-02 9.93e-02 Min 1.01e-02 1.12e-03 1.62e-03 2.40e-02 1.97e-03 3.01e-02 Avg 1.04e-02 2.61e-03 4.77e-03 3.69e-02 1.64e-02 4.78e-02 Std 5.78e-04 1.14e-02 1.56e-02 1.44e-02 2.05e-02 2.41e-02  Diversity Metric Max 7.91e-02 9.87e-02 1.19e-01 1.09e-01 1.68e-01 1.18e-01 Min 6.13e-02 4.90e-02 7.05e-02 6.32e-02 9.05e-02 8.24e-02 Avg 6.63e-02 7.02e-02 9.81e-02 8.93e-02 9.08e-02 1.00e-01 Std 7.36e-03 8.87e-03 6.78e-03 1.12e-02 1.64e-01 6.89e-03Table 4
Comparison of performance on ZDT6.
ZDT6 p-MSNSGA-II p-MSMODE p-MSMOEA/D NSGA-II MODE MOEA/D M 2 OABC Converge Metric Max 3.87e-05 6.45e-04 6.74e-04 7.72e-02 7.37e-02 7.01e-02 1.76e-03 Min 2.39e-05 2.71e-05 3.21e-04 5.31e-02 2.72e-03 5.09e-02 9.64e-05 Avg 2.95e-05 2.89e-04 4.30e-04 6.98e-02 3.71e-02 5.48e-02 5.03e-04 Std 6.08e-06 1.10e-05 7.63e-05 6.52e-03 2.17e-02 6.45e-03 7.40e-01  Diversity Metric Max 7.51e-02 7.32e-02 1.01e-01 1.12e+00 1.24e+00 1.16e+00 1.18e-01 Min 5.87e-02 4.08e-02 6.76e-02 6.12e-01 9.31e-01 6.19e-01 6.08e-02 Avg 6.37e-02 5.68e-02 8.14e-02 8.35e-01 9.99e-01 8.96e-01 7.96e-02 Std 1.26e-03 4.03e-03 2.06e-03 1.08e-01 1.91e-02 1.87e-01 1.93e-02Figure 3
PF obtained byp-MSMOEAs and MOEAs on SCH2.
(a)
PF obtained byp-MSNSGA-II on SCH2 (b)
PF obtained by NSGA-II on SCH2 (c)
PF obtained byp-MSMODE on SCH2 (d)
PF obtained by MODE on SCH2 (e)
PF obtained byp- MSMOEA/D on SCH2 (f)
PF obtained by MOEA/D on SCH2Figure 4
PF obtained byp-MSMOEAs, MOEAs, and M2OABC on ZDT3.
(a)
PF obtained byp-MSNSGA-II on ZDT3 (b)
PF obtained by NSGA-II on ZDT3 (c)
PF obtained byp-MSMODE on ZDT3 (d)
PF obtained by MODE on ZDT3 (e)
PF obtained byp-MSMOEA/D on ZDT3 (f)
PF obtained by MOEA/D on ZDT3  (g)
PF obtained by M2OABC on ZDT3Figure 5
PF obtained byp-MSMOEAs and MOEAs on ZDT4.
(a)
PF obtained byp-MSNSGA-II on ZDT4 (b)
PF obtained by NSGA-II on ZDT4 (c)
PF obtained byp-MSMODE on ZDT4 (d)
PF obtained by MODE on ZDT4 (e)
PF obtained byp-MSMOEA/D on ZDT4 (f)
PF obtained by MOEA/D on ZDT4Figure 6
PF obtained byp-MSMOEAs, MOEAs, and M2OABC on ZDT6.
(a)
PF obtained byp-MSNSGA-II on ZDT6 (b)
PF obtained by NSGA-II on ZDT6 (c)
PF obtained byp-MSMODE on ZDT6 (d)
PF obtained by MODE on ZDT6 (e)
PF obtained byp-MSMOEA/D on ZDT6 (f)
PF obtained by MOEA/D on ZDT6  (g)
PF obtained by M2OABC on ZDT6On SCH2 function, it can be noticed thatp-MSMOEAs have superior performance to MOEAs in terms of γ metric and Δ metric after 40000 FEs from Table 1. Figure 3 shows thatp-MSMOEAs have great potential to discover a well-distributed as well as diverse solution set for SCH2 function. Yet MOEA/D only finds a sparse distribution, although it can basically archive true PF for SCH2.On ZDT3, ZDT4, and ZDT6 functions, Tables2–4 show that performance ofp-MSMOEAs in both γ metric and Δ metric is better than that of MOEAs and lightly better than that of M2OABC. Figures 4–6 show that MOEAs produce poor results on these test functions and they are almost impossible to achieve true PF, whilep-MSMOEAs have great potential to approach true PF.


## 4.3.2. Three-Objective Functions

Figure7 shows the true PF for two three-objective functions. Tables 5 and 6 and Figures 8 and 9 show the optimization results ofp-MSMOEAs, MOEAs, and M2OABC algorithms on DTLZ2 and DTLZ3.Table 5
Comparison of performance on DTLZ2.
DTLZ2 p-MSNSGA-II p-MSMODE p-MSMOEA/D NSGA-II MODE MOEA/D M 2 OABC Converge Metric Max 3.02e-03 7.13e-04 5.67e-03 6.51e-02 7.97e-02 6.27e-02 3.10e-03 Min 2.88e-03 2.30e-05 6.80e-05 3.05e-02 7.50e-04 1.58e-03 2.56e-03 Avg 2.96e-03 1.96e-04 4.06e-04 4.15e-02 1.79e-03 5.51e-03 2.85e-03 Std 6.30e-04 8.48e-05 5.87e-04 8.10e-03 2.11e-02 1.01e-02 1.54e-04  Diversity Metric Max 4.60e-02 4.64e-02 7.61e-02 6.73e-01 4.92e-01 9.98e-01 4.66e-02 Min 3.37e-02 3.75e-02 2.16e-02 3.85e-01 4.23e-01 2.88e-01 3.95e-02 Avg 4.46e-02 4.21e-02 4.76e-02 5.07e-01 4.77e-01 5.39e-01 4.30e-02 Std 1.97e-03 1.88e-03 7.12e-03 5.32e-02 7.81e-02 1.29e-01 2.27e-03Table 6
Comparison of performance on DTLZ3.
DTLZ3 p-MSNSGA-II p-MSMODE p-MSMOEA/D NSGA-II MODE MOEA/D Converge Metric Max 4.10e-02 1.05e-01 8.26e-01 6.56e-01 1.08e+00 1.02e+00 Min 3.48e-02 1.86e-03 6.73e-02 5.08e-02 1.53e-01 1.78e-01 Avg 3.92e-02 4.70e-02 1.84e-01 1.21e-01 3.72e-01 3.60e-01 Std 2.55e-01 1.23e-01 1.24+00 8.65e-01 2.38e+00 2.6+00  Diversity Metric Max 6.97e-02 4.67e-02 1.21e-01 1.24e-01 4.81e-02 1.27e-01 Min 4.21e-02 3.49e-02 2.70e-02 4.54e-02 3.61e-02 3.45e-02 Avg 5.87e-02 3.75e-02 8.00e-02 9.42e-02 4.16e-02 8.02e-02 Std 7.35e-03 2.62e-03 1.52e-02 2.07e-02 2.81e-03 2.06e-02Figure 7
True PF on DTLZ2 and DTLZ3.
(a)
True PF on DTLZ2 (b)
True PF on DTLZ3Figure 8
PF obtained byp-MSMOEAs and MOEAs on DTLZ2.
(a)
PF obtained byp-MSNSGA-II on DTLZ2 (b)
PF obtained by NSGA-II on DTLZ2 (c)
PF obtained byp-MSMODE on DTLZ2 (d)
PF obtained by MODE on DTLZ2 (e)
PF obtained byp-MSMOEA/D on DTLZ2 (f)
PF obtained by MOEA/D on DTLZ2Figure 9
PF obtained byp-MSMOEAs and MOEAs on DTLZ3.
(a)
PF obtained byp-MSNSGA-II on DTLZ3 (b)
PF obtained by NSGA-II on DTLZ3 (c)
PF obtained byp-MSMODE on DTLZ3 (d)
PF obtained by MODE on DTLZ3 (e)
PF obtained byp-MSMOEA/D on DTLZ3 (f)
PF obtained by MOEA/D on DTLZ3On DTLZ2 function, when given 10000 FEs for seven algorithms, performance ofp-MSMOEAs is better than that of MOEAs and is comparable to that of M2OABC as shown in Table 5. Figure 8 shows thatp-MSMOEAs have great potential to obtain a superior PF for DTLZ2, especiallyp-MSNSGA-II andp-MSMODE. Moreover, the performance of MOEAs in Δ metric is a little worse than that ofp-MSMOEAs.On DTLZ3 function, it could be observed from Table6 that performance ofp-MSMOEAs algorithms in both γ metric and Δ metric has considerable competitiveness over this problem. From Figure 9, it can be seen that the fronts obtained fromp-MSMOEAs and MOEAs are found to have a basically outstanding result in terms of convergence, while they do not perform satisfactorily in terms of diversity. However,p-MSMOEAs are better than MOEAs for the problem.


## 4.3.3. Time Complexity Analysis

In order to demonstrate the difference in the time complexity of six algorithms, Figure10 plots the average CPU time over 10 runs.Figure 10
CPU time of MOEAs andp-MSMOEAs on six benchmark functions.Figure10 shows the CPU time for six algorithms and gives the results on six benchmark functions (SCH2, ZDT3, ZDT4, ZDT6, DTLZ2, and DTLZ3). The results show that the time complexity of the original three MOEAs is basically the same, and NSGA-II and MOEA/D are slightly better than MODE. The time complexity ofp-MSMOEAs is slightly higher than that of MOEAs, but within acceptable limits. In addition, for two-objective functions (SCH2, ZDT3, ZDT4, and ZDT6), the time complexity ofp-MSMOEAs is about 1.32-1.64 times that of the original MOEAs. However, as the number of objectives increases, this difference is more pronounced. For example, for three-objective functions (DTLZ2 and DTLZ3), the time complexity ofp-MSMOEAs is about 1.38-1.92 times that of MOEAs. It is worth mentioning that althoughp-MSMOEAs consume more CPU time, there is a significant improvement in performance compared to MOEAs.


## 5. Application for Multi-Objective Portfolio Management

### 5.1. Introduction of MOPOP

MOPOP has always had an indispensable place in modern risk management. Its ultimate goal is to find an optimal way of distributing a set of available assets, and the budget is scheduled. This is why MOPOP is favored by investors in terms of determining portfolio strategies. However, there is not a single portfolio that could fully satisfy the needs of all investors. In fact, investors' preference for risk-return ultimately determines a portfolio is an optimal one or not.


#### 5.1.1. Two-Objective Portfolio Model

Markowitz [16] proposed a formal two-objective portfolio model, which is MV portfolio model. He used the mean returns of assets and covariance of them to describe the investment return and risk, respectively. His model considers two conflicting aspects: maximizing expected return of a portfolio while minimizing its risk [28]. In the existence of risk and return, the optimal solution is not a value yet a set of optimal portfolios, which weighs between above two aspects.The MV model considered in this paper could be shown as follows [16]:(4)minσ2=∑i=1N∑j=1Nwiwjσijminimizingriskmaxr=∑i=1NwirimaximizingexpectedreturnSubjectto∑i=1Nwi=10≤wi≤1,i=1,2,…,Nbudgetconstraintwhere σ2 is variance of portfolios andwi andwj are the weights of assetsi andj in all assets, respectively.σij is covariance between above two assets.N denotes quantity of available assets, andr expresses expected return whileri represents expected return of the asseti.Although being popular in the past, the above MV model has an assumption that the return of assets is normally distributed. Unfortunately, the conditions are a bit harsh in real life and are rarely satisfied. There is now a widespread recognition of the fact that portfolios in reality do not follow a multivariate normal distribution. Skewness started to be considered in POP [29]. This implies the risk criteria of MV model could replace the variance with semi-variance. Markowitz indeed recommended that models considering semi-variance are preferable [30].


#### 5.1.2. Three-Objective Portfolio Model

In this section, a three-objective portfolio model, return-risk-cost model, is proposed. The two main innovations are the addition of transaction costs and replacement of the risk criterion. There exist three criteria in the proposed model: expected return (measured by mean of assets) that should be maximized, risk (semi-variance of assets) that should be minimized, and expected cost (Euclidean distance of weight vectors) that should be minimized.The MV model has a fact-based function limitation; that is, the weight of each asset in the portfolio should be a non-negative real number and the sum of them should be 1. In the modified Markowitz model by Fernandez and Gomez, the upper and lower limit constraints are added. The constraints of the modified model could be shown as follows:(5)∑i=1Nwi=1aizi≤wi≤bizi,i=1,2,…,N0≤ai≤1,0≤bi≤1zi=1,forwi>00,otherwise∑i=1Nzi=KThe semi-variance is described in [29] as follows:(6)ΣijB=Emin⁡Ri-B,0∙min⁡Rj-B,0=1T∙∑t=1Tmin⁡Rit-B,0∙min⁡Rjt-B,0whereB indicates the return after comparison;Rit expresses the return of asseti in periodt.An indirect method is employed in this paper to represent transaction costs. There are two assumptions that the transaction cost is related to the quantity of an asset and the distance between current portfolio as well as the expected portfolio is considered in transaction cost. According to the above assumptions, the distance which is the number of different weights in two different portfolios could be quantified as the Euclidean distance:(7)dWi,Wj=∑n=1Nwni-wnj2with higher values meaning higher transaction costs.Hence, the return–risk-cost portfolio model proposed in this paper could be formalized as follows:(8)maxE=∑i=1NwirimaximizingexpectedreturnminR=∑i=1N∑j=1NwiwjΣijminimizingriskminC=∑n=1Nwni-wnj2minimizingexpectedcostAmong them, budget constraints are given as follows:(9)Subjectto∑i=1Nwi=1aizi≤wi≤bizi,i=1,2,…,N0≤ai≤1,0≤bi≤1zi=1,forwi>00,otherwise∑i=1Nzi=KΣijB=1T∙∑t=1Tmin⁡Rit-B,0∙min⁡Rjt-B,0It is worth mentioning that the number of model objectives can also be increased if there are suitable indicators.


### 5.2. Applications for Portfolio Problem

A numerical example is provided to demonstrate the effectiveness ofp-MSMOEAs for solving MOPOP in this part. A Cesarone’s study shows that limiting the size of portfolio is capable of improving performance [31]. So historical daily data of 12 assets from Shanghai Stock Exchange are employed in the experiment, which are collected by every stock's month rates from January 2010 to December 2016.Runningp-MSMOEAs in MATLAB software, under the mentioned above models, the efficient frontiers of thep-MSMOEAs are calculated, as shown in Figures 11 and 12. In addition, since the true PF of the portfolio problems is not clear, all the PF ofp-MSMOEAs algorithms are integrated to obtain a PF as the reference PF (true PF) [32]. The hypervolume indicator has been utilized to evaluate the performance ofp-MSMOEAs. The results of comparison are listed in Tables 7 and 8.Table 7
Comparison of performance on MV portfolio model.
two-objective portfolio model p-MSNSGA-II p-MSMODE p-MSMOEA/D Hypervolume indicator Max 1.07e-02 3.65e-02 2.98e-02 Min 4.20e-03 1.97e-02 9.25e-03 Avg 8.74e-03 2.17e-02 1.35e-02 Std 5.46e-02 6.77e-03 7.50e-02Table 8
Comparison of performance on return-risk-cost portfolio model.
three-objective portfolio model p-MSNSGA-II p-MSMODE p-MSMOEA/D Hypervolume indicator Max 1.49e-01 5.41e-01 6.37e-01 Min 3.63e-02 9.86e-02 7.02e-02 Avg 8.17e-02 1.07e-01 2.83e-01 Std 3.95e-02 2.59e-01 3.91e-01Figure 11
PF obtained byp-MSMOEAs on MV portfolio model.
(a)
PF obtained byp-MSNSGA-II on MV portfolio model (b)
PF obtained byp-MSMODE on MV portfolio model (c)
PF obtained byp-MSMOEA/D on MV portfolio modelFigure 12
PF obtained byp-MSMOEAs on return-risk-cost portfolio model.
(a)
PF obtained byp-MSNSGA-II on return-risk-cost portfolio model (b)
PF obtained byp-MSMODEI on return-risk-cost portfolio model (c)
PF obtained byp-MSMOEA/D on return-risk-cost portfolio modelFigure11 gives the generated non-dominated solutions of MV model usingp-MSMOEAs. Risk and expected return are almost positively correlated, which means the greater risk the investors can accept, the greater expected return they can obtain. Investors could select a portfolio approach based on their preference for risk to get a corresponding return.Obviously, in Figure11, efficient curves ofp-MSNSGA-II andp-MSMOEA/D are smooth and continuous, whilep-MSMODE's is converse. From Table 7, it shows that performance ofp-MSNSGA-II,p-MSMOEA/D, andp-MSMODE in hypervolume indicator is decreasing in order.Figure12 gives the generated non-dominated solutions of return-risk-cost model usingp-MSMOEAs. The distributions of solutions are located in three-dimensional coordinate graphs corresponding to return, risk, and cost. As can be seen, risk and expected return are also almost positively correlated when cost is not considered. Similarly, this is also the case with cost and risk, cost and expected benefit without considering another variable, which is consistent with the reality. When considering these three variables, investors could select a portfolio approach based on their preference for risk and cost to receive corresponding benefits.Apparently, in Figure12, the distribution of non-dominated solutions in Figure 12(a) as well as Figure 12(c), which resolved byp-MSNSGA-II andp-MSMODE, is more uniform and more diverse in the whole searching space. From Table 8, it can be perceived that the performance ofp-MSNSGA-II is better thanp-MSMODE andp-MSMOEA/D in terms of hypervolume indicator.


## 5.1. Introduction of MOPOP

MOPOP has always had an indispensable place in modern risk management. Its ultimate goal is to find an optimal way of distributing a set of available assets, and the budget is scheduled. This is why MOPOP is favored by investors in terms of determining portfolio strategies. However, there is not a single portfolio that could fully satisfy the needs of all investors. In fact, investors' preference for risk-return ultimately determines a portfolio is an optimal one or not.


### 5.1.1. Two-Objective Portfolio Model

Markowitz [16] proposed a formal two-objective portfolio model, which is MV portfolio model. He used the mean returns of assets and covariance of them to describe the investment return and risk, respectively. His model considers two conflicting aspects: maximizing expected return of a portfolio while minimizing its risk [28]. In the existence of risk and return, the optimal solution is not a value yet a set of optimal portfolios, which weighs between above two aspects.The MV model considered in this paper could be shown as follows [16]:(4)minσ2=∑i=1N∑j=1Nwiwjσijminimizingriskmaxr=∑i=1NwirimaximizingexpectedreturnSubjectto∑i=1Nwi=10≤wi≤1,i=1,2,…,Nbudgetconstraintwhere σ2 is variance of portfolios andwi andwj are the weights of assetsi andj in all assets, respectively.σij is covariance between above two assets.N denotes quantity of available assets, andr expresses expected return whileri represents expected return of the asseti.Although being popular in the past, the above MV model has an assumption that the return of assets is normally distributed. Unfortunately, the conditions are a bit harsh in real life and are rarely satisfied. There is now a widespread recognition of the fact that portfolios in reality do not follow a multivariate normal distribution. Skewness started to be considered in POP [29]. This implies the risk criteria of MV model could replace the variance with semi-variance. Markowitz indeed recommended that models considering semi-variance are preferable [30].


### 5.1.2. Three-Objective Portfolio Model

In this section, a three-objective portfolio model, return-risk-cost model, is proposed. The two main innovations are the addition of transaction costs and replacement of the risk criterion. There exist three criteria in the proposed model: expected return (measured by mean of assets) that should be maximized, risk (semi-variance of assets) that should be minimized, and expected cost (Euclidean distance of weight vectors) that should be minimized.The MV model has a fact-based function limitation; that is, the weight of each asset in the portfolio should be a non-negative real number and the sum of them should be 1. In the modified Markowitz model by Fernandez and Gomez, the upper and lower limit constraints are added. The constraints of the modified model could be shown as follows:(5)∑i=1Nwi=1aizi≤wi≤bizi,i=1,2,…,N0≤ai≤1,0≤bi≤1zi=1,forwi>00,otherwise∑i=1Nzi=KThe semi-variance is described in [29] as follows:(6)ΣijB=Emin⁡Ri-B,0∙min⁡Rj-B,0=1T∙∑t=1Tmin⁡Rit-B,0∙min⁡Rjt-B,0whereB indicates the return after comparison;Rit expresses the return of asseti in periodt.An indirect method is employed in this paper to represent transaction costs. There are two assumptions that the transaction cost is related to the quantity of an asset and the distance between current portfolio as well as the expected portfolio is considered in transaction cost. According to the above assumptions, the distance which is the number of different weights in two different portfolios could be quantified as the Euclidean distance:(7)dWi,Wj=∑n=1Nwni-wnj2with higher values meaning higher transaction costs.Hence, the return–risk-cost portfolio model proposed in this paper could be formalized as follows:(8)maxE=∑i=1NwirimaximizingexpectedreturnminR=∑i=1N∑j=1NwiwjΣijminimizingriskminC=∑n=1Nwni-wnj2minimizingexpectedcostAmong them, budget constraints are given as follows:(9)Subjectto∑i=1Nwi=1aizi≤wi≤bizi,i=1,2,…,N0≤ai≤1,0≤bi≤1zi=1,forwi>00,otherwise∑i=1Nzi=KΣijB=1T∙∑t=1Tmin⁡Rit-B,0∙min⁡Rjt-B,0It is worth mentioning that the number of model objectives can also be increased if there are suitable indicators.


## 5.1.1. Two-Objective Portfolio Model

Markowitz [16] proposed a formal two-objective portfolio model, which is MV portfolio model. He used the mean returns of assets and covariance of them to describe the investment return and risk, respectively. His model considers two conflicting aspects: maximizing expected return of a portfolio while minimizing its risk [28]. In the existence of risk and return, the optimal solution is not a value yet a set of optimal portfolios, which weighs between above two aspects.The MV model considered in this paper could be shown as follows [16]:(4)minσ2=∑i=1N∑j=1Nwiwjσijminimizingriskmaxr=∑i=1NwirimaximizingexpectedreturnSubjectto∑i=1Nwi=10≤wi≤1,i=1,2,…,Nbudgetconstraintwhere σ2 is variance of portfolios andwi andwj are the weights of assetsi andj in all assets, respectively.σij is covariance between above two assets.N denotes quantity of available assets, andr expresses expected return whileri represents expected return of the asseti.Although being popular in the past, the above MV model has an assumption that the return of assets is normally distributed. Unfortunately, the conditions are a bit harsh in real life and are rarely satisfied. There is now a widespread recognition of the fact that portfolios in reality do not follow a multivariate normal distribution. Skewness started to be considered in POP [29]. This implies the risk criteria of MV model could replace the variance with semi-variance. Markowitz indeed recommended that models considering semi-variance are preferable [30].


## 5.1.2. Three-Objective Portfolio Model

In this section, a three-objective portfolio model, return-risk-cost model, is proposed. The two main innovations are the addition of transaction costs and replacement of the risk criterion. There exist three criteria in the proposed model: expected return (measured by mean of assets) that should be maximized, risk (semi-variance of assets) that should be minimized, and expected cost (Euclidean distance of weight vectors) that should be minimized.The MV model has a fact-based function limitation; that is, the weight of each asset in the portfolio should be a non-negative real number and the sum of them should be 1. In the modified Markowitz model by Fernandez and Gomez, the upper and lower limit constraints are added. The constraints of the modified model could be shown as follows:(5)∑i=1Nwi=1aizi≤wi≤bizi,i=1,2,…,N0≤ai≤1,0≤bi≤1zi=1,forwi>00,otherwise∑i=1Nzi=KThe semi-variance is described in [29] as follows:(6)ΣijB=Emin⁡Ri-B,0∙min⁡Rj-B,0=1T∙∑t=1Tmin⁡Rit-B,0∙min⁡Rjt-B,0whereB indicates the return after comparison;Rit expresses the return of asseti in periodt.An indirect method is employed in this paper to represent transaction costs. There are two assumptions that the transaction cost is related to the quantity of an asset and the distance between current portfolio as well as the expected portfolio is considered in transaction cost. According to the above assumptions, the distance which is the number of different weights in two different portfolios could be quantified as the Euclidean distance:(7)dWi,Wj=∑n=1Nwni-wnj2with higher values meaning higher transaction costs.Hence, the return–risk-cost portfolio model proposed in this paper could be formalized as follows:(8)maxE=∑i=1NwirimaximizingexpectedreturnminR=∑i=1N∑j=1NwiwjΣijminimizingriskminC=∑n=1Nwni-wnj2minimizingexpectedcostAmong them, budget constraints are given as follows:(9)Subjectto∑i=1Nwi=1aizi≤wi≤bizi,i=1,2,…,N0≤ai≤1,0≤bi≤1zi=1,forwi>00,otherwise∑i=1Nzi=KΣijB=1T∙∑t=1Tmin⁡Rit-B,0∙min⁡Rjt-B,0It is worth mentioning that the number of model objectives can also be increased if there are suitable indicators.


## 5.2. Applications for Portfolio Problem

A numerical example is provided to demonstrate the effectiveness ofp-MSMOEAs for solving MOPOP in this part. A Cesarone’s study shows that limiting the size of portfolio is capable of improving performance [31]. So historical daily data of 12 assets from Shanghai Stock Exchange are employed in the experiment, which are collected by every stock's month rates from January 2010 to December 2016.Runningp-MSMOEAs in MATLAB software, under the mentioned above models, the efficient frontiers of thep-MSMOEAs are calculated, as shown in Figures 11 and 12. In addition, since the true PF of the portfolio problems is not clear, all the PF ofp-MSMOEAs algorithms are integrated to obtain a PF as the reference PF (true PF) [32]. The hypervolume indicator has been utilized to evaluate the performance ofp-MSMOEAs. The results of comparison are listed in Tables 7 and 8.Table 7
Comparison of performance on MV portfolio model.
two-objective portfolio model p-MSNSGA-II p-MSMODE p-MSMOEA/D Hypervolume indicator Max 1.07e-02 3.65e-02 2.98e-02 Min 4.20e-03 1.97e-02 9.25e-03 Avg 8.74e-03 2.17e-02 1.35e-02 Std 5.46e-02 6.77e-03 7.50e-02Table 8
Comparison of performance on return-risk-cost portfolio model.
three-objective portfolio model p-MSNSGA-II p-MSMODE p-MSMOEA/D Hypervolume indicator Max 1.49e-01 5.41e-01 6.37e-01 Min 3.63e-02 9.86e-02 7.02e-02 Avg 8.17e-02 1.07e-01 2.83e-01 Std 3.95e-02 2.59e-01 3.91e-01Figure 11
PF obtained byp-MSMOEAs on MV portfolio model.
(a)
PF obtained byp-MSNSGA-II on MV portfolio model (b)
PF obtained byp-MSMODE on MV portfolio model (c)
PF obtained byp-MSMOEA/D on MV portfolio modelFigure 12
PF obtained byp-MSMOEAs on return-risk-cost portfolio model.
(a)
PF obtained byp-MSNSGA-II on return-risk-cost portfolio model (b)
PF obtained byp-MSMODEI on return-risk-cost portfolio model (c)
PF obtained byp-MSMOEA/D on return-risk-cost portfolio modelFigure11 gives the generated non-dominated solutions of MV model usingp-MSMOEAs. Risk and expected return are almost positively correlated, which means the greater risk the investors can accept, the greater expected return they can obtain. Investors could select a portfolio approach based on their preference for risk to get a corresponding return.Obviously, in Figure11, efficient curves ofp-MSNSGA-II andp-MSMOEA/D are smooth and continuous, whilep-MSMODE's is converse. From Table 7, it shows that performance ofp-MSNSGA-II,p-MSMOEA/D, andp-MSMODE in hypervolume indicator is decreasing in order.Figure12 gives the generated non-dominated solutions of return-risk-cost model usingp-MSMOEAs. The distributions of solutions are located in three-dimensional coordinate graphs corresponding to return, risk, and cost. As can be seen, risk and expected return are also almost positively correlated when cost is not considered. Similarly, this is also the case with cost and risk, cost and expected benefit without considering another variable, which is consistent with the reality. When considering these three variables, investors could select a portfolio approach based on their preference for risk and cost to receive corresponding benefits.Apparently, in Figure12, the distribution of non-dominated solutions in Figure 12(a) as well as Figure 12(c), which resolved byp-MSNSGA-II andp-MSMODE, is more uniform and more diverse in the whole searching space. From Table 8, it can be perceived that the performance ofp-MSNSGA-II is better thanp-MSMODE andp-MSMOEA/D in terms of hypervolume indicator.


## 6. Conclusion

In this paper, a multi-swarm multi-objective optimizer based onp-optimality criteria calledp-MSMOEAs is proposed. Inp-MSMOEAs,p-optimality criteria are employed to ensure algorithms converge to true PF. In addition, multiple swarm cooperative coevolution is adopted, guaranteeing the diversity of the whole population.P-MSMOEAs are simply constructed and easily achieved and have considerable potential to solve complex MOP. With six mathematical benchmark functions,p-MSMOEAs are proven to obtain good distributed PF with respect to optimization accuracy and convergence robust.Additionally, a constrained MOPOP with three criteria, expected return (mean of return), risk (semi-variance of return), and expected cost (the Euclidean distance of weight vectors), is proposed. Andp-MSMOEAs are utilized to resolve the return-risk-cost portfolio model. The proposedp-MSMOEAs are considered to be capable of getting a good diversity of solutions' distribution.



---

*Source: 8418369-2019-01-21.xml*