# A Multistrategy-Based Multiobjective Differential Evolution for Optimal Control in Chemical Processes

Authors: Bin Xu; Xu Chen; Xiuhui Huang; Lili Tao
Journal: Complexity
 (2018)

---

## Abstract

Optimal control problems with multiple conflicting objectives in chemical processes are quite challenging. To solve such problems, we put forward a multistrategy-based multiobjective differential evolution, in which (1) a hybrid selection strategy is incorporated from the motivation of no single strategy outperforming all other ones in every stage; (2) a multipopulation strategy is applied to represent the main population and current optimum, and a cyclic crowding estimation is developed to maintain these optimum; and (3) a multimutation strategy is constructed to improve both exploration and exploitation ability. The effectiveness and efficiency of the proposed algorithm are validated by comparisons with some representative multiobjective evolutionary algorithms over 12 test instances. Moreover, the proposed algorithm is applied to solve 3 multiobjective optimal control problems in chemical processes. The obtained results indicate the efficiency and effectiveness of the proposed algorithm for solving multiobjective optimal control problems.


---

## Body

## 1. Introduction

Almost all optimal control problems (OCPs) in engineering have been considered to be single-objective optimization problems. In practical, for most of OCPs in chemical process, they involve multiple conflicting objectives [1]. For instance, some chemical processes are characterized by improving both product quality and yield with low energy consumption [2]. For such multiobjective optimal control problems (MOOCPs), there does not exist a trajectory which optimizes all objectives simultaneously. On the contrary, there exist many optimal trajectories, called Pareto optimal solutions in a multiobjective optimization community, which represent some trade-off between the conflicting objectives.One common and straightforward manner to solve a given MOOCP is transforming the original problem into a single-objective optimal control problem (SOOCP) using some weighted sum methods and then optimizing it [3–5]. One intrinsic characteristic common in the single-objective-based methods is that they aim at locating just one optimal solution rather than multiple optimal solutions in a single run. Thus, many runs should be conducted to obtain the multiple optimal solutions. Furthermore, the relation between a uniform grid of weights and the spread on the Pareto front is unclear.With the development of evolutionary algorithms (EAs), multiobjective evolutionary algorithms (MOEAs) have been widely used in plenty of fields [6]. One advantage of using MOEAs to solve such problems is that they can find a set of representative optimal solutions in a single run [7]. Recently, some researchers have demonstrated that MOEAs not only are effective for common multiobjective problems (MOPs) but also can be extended to solve MOOCPs. Sarkar and Modak used the enhanced version of a nondominated sorting genetic algorithm (NSGA-II) to solve two fed-batch bioreactors [8]. Logist et al. used the freely available toolkit named ACADO, to solve different types of MOOCPs appearing in chemical engineering with different objective numbers [1]. Chen et al. improved the performance of MODE by using a ranking-based mutation operator and then applied the new MODE to solve 4 MOOCPs in chemical process [9]. Fan et al. proposed a multiobjective differential evolution with performance metric-based self-adaptive mutation operator for typical MOOCPs in chemical and biochemical processes [10]. Results showed that the proposed new approach is suitable for solving actual problems and can obtain some Pareto optimum for decision makers.When solving the MOOCP by MOEAs, the principle of it, as shown in Figure1, can be summarized as follows: firstly, transforming the MOOCP into common MOPs and then MOEAs can be applied to solve the transformed problems. Therefore, questions arising naturally in solving MOOCP are twofold: how to perform transforming and how to design effective MOEAs. In this paper, we focus on the second issue and take the most used control vector parameterization (CVP) method to meet the first one.Figure 1
The schematic graph to illustrate the principle of solving an OCP by EAs.The last two decades have witnessed major developments in EA-based algorithms for MOPs. Differential evolution (DE) [11], a simple yet effective nature-inspired stochastic algorithm, has been widely and successfully applied for solving various complex problems [9, 12–14]. However, it must be emphasized that DE was originally proposed for single-objective problems, while the main goal of a multiobjective optimization algorithm is to find a set of trade-off solutions. Thus, to solve MOPs effectively via DE, the following two issues should be considered:
(1)
How to select and/or retain the best solutions? That is, how to perform elitism.
(2)
How to promote diversity into the population.In fact, the basic DE can be directly applied to solve MOPs by replacing the selection component [6]. Thus, one key issue in designing a MODE is designing a proper selection component. In general, the selection component includes assigning a fitness value for each candidate solution and selecting the promising ones for the new population. As pointed out in [15], different selection components in DE makes an algorithm differ from others. During the last three decades, plenty of selection components with different strategies have been developed for solving MOPs. One notable difference among these strategies is how the Pareto dominance concept is used. From this aspect, we classify the current selection strategies into three categories:


### 1.1. No Pareto Dominance Used

It combines all the multiple objectives into a single one, such as weighted sum approach, and then the original EAs can be applied, such as [16–19] and so on. Among all these algorithms, MOEA/D is a recent multiobjective evolutionary algorithmic framework which is based on conventional aggregation approaches [18]. MOEA/D has a fast execution speed and some improved MOEA/D-based algorithms are proposed [7, 20]. However, they need to provide a set of uniformly distributed weights and a decomposition method. This extra information highly affects MOEA/D’s performance.


### 1.2. Use Pareto Dominance as a Ranking Method

In this group, Pareto ranking, also called the nondominated sorting method [21, 22], is used for fitness assignment for the conflicting objectives first and then the selection based on the ranking numbers can be performed. Much research uses this method as a selection strategy [21–29].


### 1.3. Use Pareto Dominance as a Filter

Here, a filter means finding a set of better solutions from a given set. In this situation, Pareto dominance is often used in tournament selection to find a better solution among several candidates or to identify a set of nondominated solutions [30–37]. Since this group of strategies puts more efforts on the elitists, it has a very fast convergence speed. However, it easily gets the search stuck at local optima for multimodal problems.According to the no free lunch theorem [38], it is impossible for a single-selection strategy to outperform all other strategies on all problems at different search stages. In other words, depending on several selection strategies can be effective for different problems during different search stages. Motivated by these observations, we propose a new MOEA named DEHC using DE with a hybrid selection strategy (HSS) and cyclic crowding estimation (CCE) for MOOCPs. For the hybrid selection strategy, the decomposition-based strategy [18], nondominated sorting strategy [22], and nondominated neighbor-based strategy [33] are adopted as the selection strategy. As for the diversity estimation, in this paper, we introduce an external population to store all the found nondominated solutions and then use a cyclic crowding estimation to maintain this population. Furthermore, a multiple-mutation strategy-based mutation component is used as the main reproduction operator.The remainder of this paper is arranged as follows. Section2 briefly presents a basic concept of MOOCP and conventional DE. Section 3 presents our DEHC in detail. Section 4 presents simulation and comparative results of DEHC with other competing algorithms. In Section 5, DEHC is applied to solve 3 MOOCPs in chemical processes. Section 6 draws conclusions.


## 1.1. No Pareto Dominance Used

It combines all the multiple objectives into a single one, such as weighted sum approach, and then the original EAs can be applied, such as [16–19] and so on. Among all these algorithms, MOEA/D is a recent multiobjective evolutionary algorithmic framework which is based on conventional aggregation approaches [18]. MOEA/D has a fast execution speed and some improved MOEA/D-based algorithms are proposed [7, 20]. However, they need to provide a set of uniformly distributed weights and a decomposition method. This extra information highly affects MOEA/D’s performance.


## 1.2. Use Pareto Dominance as a Ranking Method

In this group, Pareto ranking, also called the nondominated sorting method [21, 22], is used for fitness assignment for the conflicting objectives first and then the selection based on the ranking numbers can be performed. Much research uses this method as a selection strategy [21–29].


## 1.3. Use Pareto Dominance as a Filter

Here, a filter means finding a set of better solutions from a given set. In this situation, Pareto dominance is often used in tournament selection to find a better solution among several candidates or to identify a set of nondominated solutions [30–37]. Since this group of strategies puts more efforts on the elitists, it has a very fast convergence speed. However, it easily gets the search stuck at local optima for multimodal problems.According to the no free lunch theorem [38], it is impossible for a single-selection strategy to outperform all other strategies on all problems at different search stages. In other words, depending on several selection strategies can be effective for different problems during different search stages. Motivated by these observations, we propose a new MOEA named DEHC using DE with a hybrid selection strategy (HSS) and cyclic crowding estimation (CCE) for MOOCPs. For the hybrid selection strategy, the decomposition-based strategy [18], nondominated sorting strategy [22], and nondominated neighbor-based strategy [33] are adopted as the selection strategy. As for the diversity estimation, in this paper, we introduce an external population to store all the found nondominated solutions and then use a cyclic crowding estimation to maintain this population. Furthermore, a multiple-mutation strategy-based mutation component is used as the main reproduction operator.The remainder of this paper is arranged as follows. Section2 briefly presents a basic concept of MOOCP and conventional DE. Section 3 presents our DEHC in detail. Section 4 presents simulation and comparative results of DEHC with other competing algorithms. In Section 5, DEHC is applied to solve 3 MOOCPs in chemical processes. Section 6 draws conclusions.


## 2. Basic Concepts

### 2.1. MOOCP: Multiobjective Optimal Control Problem

Consider the systemẋ=fxt,ut,t with state variables xt∈ℝD, initial condition xt0=x0, and control variables ut∈ℝK. Then the main goal of solving the problem is to find the optimal control variables u∗t,t∈t0,tf that drive the plant along the trajectory x∗t,t∈t0,tf such that all the cost functions are minimized where the final time tf is fixed. The structure of a typical MOOCP for a continuous process can be generally described as follows:
(1)minJ1,J2,…JM,subject to
(2)dxdt=fxt,ut,t,t∈t0,tf,L≤u≤U.Each of the objective functionsJii=1,2,…,M in (1) can be formulated as follows:
(3)Jiut=φixtf,tf+∫t0tfψixt,ut,tdt,where M is the objective number, Ji⋅ is the total cost function, ϕi⋅ is the final time performance index, ψi⋅ is the integrated performance during the operation, t0 is the initial time, and tf is the final time; L∈ℝK and U∈ℝK are the lower and upper boundary of control variables, respectively.


### 2.2. Multiobjective Optimization

Multiobjective optimization algorithms are aimed at optimizing conflicting objectives simultaneously. To aid descriptions, some definitions regarding multiobjective optimization are introduced.Definition 1.
Solutionu1 is said to dominate solution u2, denoted as u1≻u2, if and only if ∀i∈1,…,M, Jiu1≤Jiu2, and ∃j∈1,…,M, Jju1<Jju2. Reciprocally, solution u2 is said to be dominated by solution u1, denoted as u2≺u1.Definition 2.
Solutionu∗ is said to be a Pareto optimum only if ¬∃u∈S, u≻u∗. All Pareto optima constitute a Pareto optimal set, denoted as X∗.Definition 3.
Pareto optimal front is defined asPF=Ju∗∣u∗∈X∗, that is, the mapping of the Pareto optimal set in the objective space.


### 2.3. Differential Evolution Algorithm

Like other EAs, DE begins with a randomly initialized population in the search space and then adopts trial vector generation and selection operators sequentially at each generation to move the main population toward the global optimum.The trial vector generation comprises two operators: mutation and crossover operators. The following are the five frequently used mutation strategies in the literature:(1)
DE/rand/1
(4)vig=xr1g+F⋅xr2g−xr3g
(2)
DE/best/1
(5)vig=xbestg+F⋅xr1g−xr2g
(3)
DE/current-to-best/1
(6)vig=xig+F⋅xbestg−xig+F⋅xr1g−xr2g
(4)
DE/best/2
(7)vig=xbestg+F⋅xr1g−xr2g+F⋅xr3g−xr4g
(5)
DE/rand/2
(8)vig=xr1g+F⋅xr2g−xr3g+F⋅xr4g−xr5gwhere indices ri, i=1,…,5 are mutually exclusive integers randomly generated within 1,N, respectively, that are also different from index i; xbestg is the best vector in the population at generation g; F is the scale factor and is usually chosen between 0 and 1.The most frequently used binomial crossover operation is performed as follows:(9)ui,jg=vi.jg,ifrand≤CRorj=sn,xi,jg,otherwise,j=1,…,n,where sn is an arbitrary number in 1,2,…,n and ui,j is the jth element of the ith new trial vector.After generating the trial vectors, a greedy selection operator is performed as follows:(10)xig+1=uig,iffuig≤fxig,xig,otherwise,where xig+1 is the ith solution in the main population of the next generation.


## 2.1. MOOCP: Multiobjective Optimal Control Problem

Consider the systemẋ=fxt,ut,t with state variables xt∈ℝD, initial condition xt0=x0, and control variables ut∈ℝK. Then the main goal of solving the problem is to find the optimal control variables u∗t,t∈t0,tf that drive the plant along the trajectory x∗t,t∈t0,tf such that all the cost functions are minimized where the final time tf is fixed. The structure of a typical MOOCP for a continuous process can be generally described as follows:
(1)minJ1,J2,…JM,subject to
(2)dxdt=fxt,ut,t,t∈t0,tf,L≤u≤U.Each of the objective functionsJii=1,2,…,M in (1) can be formulated as follows:
(3)Jiut=φixtf,tf+∫t0tfψixt,ut,tdt,where M is the objective number, Ji⋅ is the total cost function, ϕi⋅ is the final time performance index, ψi⋅ is the integrated performance during the operation, t0 is the initial time, and tf is the final time; L∈ℝK and U∈ℝK are the lower and upper boundary of control variables, respectively.


## 2.2. Multiobjective Optimization

Multiobjective optimization algorithms are aimed at optimizing conflicting objectives simultaneously. To aid descriptions, some definitions regarding multiobjective optimization are introduced.Definition 1.
Solutionu1 is said to dominate solution u2, denoted as u1≻u2, if and only if ∀i∈1,…,M, Jiu1≤Jiu2, and ∃j∈1,…,M, Jju1<Jju2. Reciprocally, solution u2 is said to be dominated by solution u1, denoted as u2≺u1.Definition 2.
Solutionu∗ is said to be a Pareto optimum only if ¬∃u∈S, u≻u∗. All Pareto optima constitute a Pareto optimal set, denoted as X∗.Definition 3.
Pareto optimal front is defined asPF=Ju∗∣u∗∈X∗, that is, the mapping of the Pareto optimal set in the objective space.


## 2.3. Differential Evolution Algorithm

Like other EAs, DE begins with a randomly initialized population in the search space and then adopts trial vector generation and selection operators sequentially at each generation to move the main population toward the global optimum.The trial vector generation comprises two operators: mutation and crossover operators. The following are the five frequently used mutation strategies in the literature:(1)
DE/rand/1
(4)vig=xr1g+F⋅xr2g−xr3g
(2)
DE/best/1
(5)vig=xbestg+F⋅xr1g−xr2g
(3)
DE/current-to-best/1
(6)vig=xig+F⋅xbestg−xig+F⋅xr1g−xr2g
(4)
DE/best/2
(7)vig=xbestg+F⋅xr1g−xr2g+F⋅xr3g−xr4g
(5)
DE/rand/2
(8)vig=xr1g+F⋅xr2g−xr3g+F⋅xr4g−xr5gwhere indices ri, i=1,…,5 are mutually exclusive integers randomly generated within 1,N, respectively, that are also different from index i; xbestg is the best vector in the population at generation g; F is the scale factor and is usually chosen between 0 and 1.The most frequently used binomial crossover operation is performed as follows:(9)ui,jg=vi.jg,ifrand≤CRorj=sn,xi,jg,otherwise,j=1,…,n,where sn is an arbitrary number in 1,2,…,n and ui,j is the jth element of the ith new trial vector.After generating the trial vectors, a greedy selection operator is performed as follows:(10)xig+1=uig,iffuig≤fxig,xig,otherwise,where xig+1 is the ith solution in the main population of the next generation.


## 3. Our Proposed DEHC

Because of the conflicting objectives, the DE can not be applied to solve MOPs directly. To overcome this limitation, the proposed DEHC has the following features:(1)
In the mutation step, (4) and (5) constitute a strategy candidate population pool and a random strategy is determined from this candidate pool for each solution.
(2)
In the selection step, the greedy strategy-based selection component is replaced by a new selection component based on hybrid selection strategy.
(3)
An external population is introduced to store the best solutions found so far and to provide best solutions when strategy (5) is adopted.The overall flowchart of DEHC is shown in Figure2. Below, we elaborate on some of the main steps in detail.Figure 2
The flowchart of DEHC.


### 3.1. Mutation Component with Two Strategies

For any MOPs, there is always more than one best solution, which are also called nondominated solutions. In our approach, an external population, denoted asE, is introduced and the best solutions found are stored in it. This has two purposes:
(1)
The solutions in the external population are the final results reported to the performer(s).
(2)
The external population is the best solution candidate pool for the best solution based on mutation strategy.In DEHC, (4) and (5) are selected to constitute a strategy candidate pool. In the mutation step, a strategy is randomly selected from the pool for each solution. To perform the best solution-based strategy (5), a best solution must be determined. Inspired by Coello et al.’s method [39], the best solution is determined as follows: first, the crowding distance of every solution in E is calculated [22], and then two solutions are randomly selected from the population. Their crowding distances are compared and the greater will be the final best solution for participating in the mutation operator. Note that a new best solution needs to be selected for each solution if the mutation strategy is (5). The purposes for this extension are twofold:
(1)
Strategy (4) has advantage of putting more effort into improving the search engine’s exploration ability.
(2)
Strategy (5) has advantage of putting more effort into improving the search engine’s exploitation ability. Furthermore, it has advantage of putting more effort into exploiting the less crowded area.


### 3.2. Selection Component Based on Hybrid Selection Strategy

Selection component is the most essential step for any MOEAs because it is responsible for guiding the selection process at the various stages of the algorithm toward a Pareto optimal front. Many selection components with different selection strategies have been proposed for solving MOPs. Here, in DEHC, we use a new selection component called hybrid selection strategy (HSS) from the view of multipopulation. The HSS consist of three selection strategies and each strategy has its own subpopulation. These three strategies are called decomposition-based strategy [18], nondominated sorting strategy [22], and nondominated neighbor-based strategy [33]. These three strategies fall into the above three groups in Sections 1.1, 1.2, and 1.3, respectively. The reasons for hybridizing these three selection strategies are as follows:
(1)
Decomposition-based strategy has high execution speed; however, it needs providing a set of uniformly distributed weight. Although determining a set of weight is an easy job for problems with 2 objectives and 3 objectives, it becomes a tough thing as the number of objectives increases. Furthermore, it needs a decomposition approach. Most important of all, these weight values and decomposition approach highly affect the final results.
(2)
Nondominated sorting strategy can maintain the diversity of the main population. But it has aOmN2 computational complexity (where N is the population size). So, less population size means lower computational time.
(3)
Nondominated neighbor-based strategy has a very fast convergence speed because it puts more efforts on the best solutions. This high selection pressure may mislead the search into local optima.The proposed HSS works as follows: generating the offspring for each subpopulation and then assigning the fitness value of each solution in the combination of its subpopulation and offspring subpopulation according to the corresponding selection strategy and finally selecting a proper number of best solutions to form its subpopulation of next generation. This process is described in Figure3.Figure 3
The main selection process using HSS.The main characteristic of HSS is that it utilizes the advantage of every selection strategy such as the high execution speed of the decomposition based-strategy, high population diversity of the nondominated sorting strategy, and fast convergence speed of the nondominated neighbor-based strategy. So, different selection strategies have different trade-offs between the conflicting objectives and the hybrid strategy can take advantage of all these trade-off. Although each strategy has its own subpopulation, these subpopulations are not independent. In the mutation step, all 3 subpopulations are combined together to generate new offspring individuals. Furthermore, all the best solutions found are stored in an external population and this population provides guiding information in performing the best solution-based mutation strategy in Section3.1.


### 3.3. Maintenance of External Population Based on Cyclic Crowding Estimation

Maintaining an external population is beneficial [29, 30], but two issues arise in the maintenance. One is how to add the currently found best solutions into the population and the other is how to remove the inferior ones. Denote the external population as E and the currently found best solutions in the offspring as Z. Then we can add Z into E by Algorithm 1.Algorithm 1: Pseudocode of CCE.
1:for each zi∈Z do
2:if zi is dominated by any member of E then
3:   discardzi
4:end if
5:if zi dominates a set of members Dzi from E then
6:E=E\Dzi;
7:E=E∪zi
8:end if
9:  if AHC andzi are non-dominated with each other then
10:E=E∪zi
11:end if
12:end forAfter addingZ into E, we need to check whether its size exceeds the size limit since the number of nondominated solutions can be huge and the computational complexity of maintaining all of them is high. Hence, we truncate the population to a predefined number. To do this, the crowding distance is the most commonly applied mechanism as for its simplification and free of parameters [22, 29, 33]. This method removes all redundant solutions at once. However, the removing of a solution will lead to the change of its neighborhood’s crowding distance value. From this perspective, the one-time operator does not provide good result in all cases.To overcome this drawback, CCE is presented. The main characteristic of CCE is to remove the most crowded solution one by one or cyclically instead of retaining the least crowed solutions at a time. This process can be described by Algorithm2.Algorithm 2: The main truncating process by the CCE process.
Input: The external population E with Q solutions
Output: The truncated external population E with N solutions
1: LetK=Q−N;
2:for each xi∈Edo
3:  Calculate the Euclidean distancedi,j between xi and xj, where j=i,…,Q;
4:  Sortdi,j, j=i,…,Q, in ascending order;
5:  Add the firstK+1 values and the corresponding solutions as K+1
records intoD;
6:end for
7:while K>0 do
8:  Find the record with minimal distance inD and denote this value as
Vab and the corresponding two solutions as Pa and Pb, respectively;
9:  Remove this record fromD;
10:  Find the record with minimal distance values associated withPa and
Pb in D, respectively; denote the values as Va and Vb, respectively;
11:if Va≤Vb then
12:   RemovePa from E
13:   Remove all records associated withPa in D
14:else
15:   RemovePb from E
16:   Remove all records associated withPb in D
17:end if
18:K=K−1
19:end while


### 3.4. Computational Complexity of DEHC

As pointed above, the DEHC differs from other MOEAs in three main features, that is, multiple-mutation operators, hybrid selection strategy, and cyclic crowding estimation. To use multiple-mutation strategies, we need to find the current best vector in each generation loop. This procedure depends solely onM-independent sortings at N solutions, which has a OMNlogN computational complexity.For the hybrid selection strategy, each selection strategy has a subpopulation with sizeNs=1/3N. The decomposition-based strategy has a time complexity of OMNsT, where T is the neighborhood size. For the nondominated sorting strategy, the computational complexity is OMNs2. As for the nondominated neighbor-based strategy, the worst computational complexity is ONslgNs. So, the worst total computational complexity of HSS is
(11)OMNsT+OMNs2+ONslogNs.As for the CCE, in the worst case, the time complexity for calculating the distances isNQ2; the time complexity for adding the K+1 values is OKQ; the time complexity for removing K solutions from E is OKQ2.Based on the above analysis, the procedure of HSS does not increase any burden on the runtime complexity than using the single-selection strategy. As for the multiple-mutation strategies, it needs higher computational complexity, but it does not dramatically impose any serious burden on the runtime complexity than the single-mutation strategy-based operator. As for the CCE procedure, it needs much higher computational complexity than the complexity of crowding distance, which is onlyOMQlogQ.


## 3.1. Mutation Component with Two Strategies

For any MOPs, there is always more than one best solution, which are also called nondominated solutions. In our approach, an external population, denoted asE, is introduced and the best solutions found are stored in it. This has two purposes:
(1)
The solutions in the external population are the final results reported to the performer(s).
(2)
The external population is the best solution candidate pool for the best solution based on mutation strategy.In DEHC, (4) and (5) are selected to constitute a strategy candidate pool. In the mutation step, a strategy is randomly selected from the pool for each solution. To perform the best solution-based strategy (5), a best solution must be determined. Inspired by Coello et al.’s method [39], the best solution is determined as follows: first, the crowding distance of every solution in E is calculated [22], and then two solutions are randomly selected from the population. Their crowding distances are compared and the greater will be the final best solution for participating in the mutation operator. Note that a new best solution needs to be selected for each solution if the mutation strategy is (5). The purposes for this extension are twofold:
(1)
Strategy (4) has advantage of putting more effort into improving the search engine’s exploration ability.
(2)
Strategy (5) has advantage of putting more effort into improving the search engine’s exploitation ability. Furthermore, it has advantage of putting more effort into exploiting the less crowded area.


## 3.2. Selection Component Based on Hybrid Selection Strategy

Selection component is the most essential step for any MOEAs because it is responsible for guiding the selection process at the various stages of the algorithm toward a Pareto optimal front. Many selection components with different selection strategies have been proposed for solving MOPs. Here, in DEHC, we use a new selection component called hybrid selection strategy (HSS) from the view of multipopulation. The HSS consist of three selection strategies and each strategy has its own subpopulation. These three strategies are called decomposition-based strategy [18], nondominated sorting strategy [22], and nondominated neighbor-based strategy [33]. These three strategies fall into the above three groups in Sections 1.1, 1.2, and 1.3, respectively. The reasons for hybridizing these three selection strategies are as follows:
(1)
Decomposition-based strategy has high execution speed; however, it needs providing a set of uniformly distributed weight. Although determining a set of weight is an easy job for problems with 2 objectives and 3 objectives, it becomes a tough thing as the number of objectives increases. Furthermore, it needs a decomposition approach. Most important of all, these weight values and decomposition approach highly affect the final results.
(2)
Nondominated sorting strategy can maintain the diversity of the main population. But it has aOmN2 computational complexity (where N is the population size). So, less population size means lower computational time.
(3)
Nondominated neighbor-based strategy has a very fast convergence speed because it puts more efforts on the best solutions. This high selection pressure may mislead the search into local optima.The proposed HSS works as follows: generating the offspring for each subpopulation and then assigning the fitness value of each solution in the combination of its subpopulation and offspring subpopulation according to the corresponding selection strategy and finally selecting a proper number of best solutions to form its subpopulation of next generation. This process is described in Figure3.Figure 3
The main selection process using HSS.The main characteristic of HSS is that it utilizes the advantage of every selection strategy such as the high execution speed of the decomposition based-strategy, high population diversity of the nondominated sorting strategy, and fast convergence speed of the nondominated neighbor-based strategy. So, different selection strategies have different trade-offs between the conflicting objectives and the hybrid strategy can take advantage of all these trade-off. Although each strategy has its own subpopulation, these subpopulations are not independent. In the mutation step, all 3 subpopulations are combined together to generate new offspring individuals. Furthermore, all the best solutions found are stored in an external population and this population provides guiding information in performing the best solution-based mutation strategy in Section3.1.


## 3.3. Maintenance of External Population Based on Cyclic Crowding Estimation

Maintaining an external population is beneficial [29, 30], but two issues arise in the maintenance. One is how to add the currently found best solutions into the population and the other is how to remove the inferior ones. Denote the external population as E and the currently found best solutions in the offspring as Z. Then we can add Z into E by Algorithm 1.Algorithm 1: Pseudocode of CCE.
1:for each zi∈Z do
2:if zi is dominated by any member of E then
3:   discardzi
4:end if
5:if zi dominates a set of members Dzi from E then
6:E=E\Dzi;
7:E=E∪zi
8:end if
9:  if AHC andzi are non-dominated with each other then
10:E=E∪zi
11:end if
12:end forAfter addingZ into E, we need to check whether its size exceeds the size limit since the number of nondominated solutions can be huge and the computational complexity of maintaining all of them is high. Hence, we truncate the population to a predefined number. To do this, the crowding distance is the most commonly applied mechanism as for its simplification and free of parameters [22, 29, 33]. This method removes all redundant solutions at once. However, the removing of a solution will lead to the change of its neighborhood’s crowding distance value. From this perspective, the one-time operator does not provide good result in all cases.To overcome this drawback, CCE is presented. The main characteristic of CCE is to remove the most crowded solution one by one or cyclically instead of retaining the least crowed solutions at a time. This process can be described by Algorithm2.Algorithm 2: The main truncating process by the CCE process.
Input: The external population E with Q solutions
Output: The truncated external population E with N solutions
1: LetK=Q−N;
2:for each xi∈Edo
3:  Calculate the Euclidean distancedi,j between xi and xj, where j=i,…,Q;
4:  Sortdi,j, j=i,…,Q, in ascending order;
5:  Add the firstK+1 values and the corresponding solutions as K+1
records intoD;
6:end for
7:while K>0 do
8:  Find the record with minimal distance inD and denote this value as
Vab and the corresponding two solutions as Pa and Pb, respectively;
9:  Remove this record fromD;
10:  Find the record with minimal distance values associated withPa and
Pb in D, respectively; denote the values as Va and Vb, respectively;
11:if Va≤Vb then
12:   RemovePa from E
13:   Remove all records associated withPa in D
14:else
15:   RemovePb from E
16:   Remove all records associated withPb in D
17:end if
18:K=K−1
19:end while


## 3.4. Computational Complexity of DEHC

As pointed above, the DEHC differs from other MOEAs in three main features, that is, multiple-mutation operators, hybrid selection strategy, and cyclic crowding estimation. To use multiple-mutation strategies, we need to find the current best vector in each generation loop. This procedure depends solely onM-independent sortings at N solutions, which has a OMNlogN computational complexity.For the hybrid selection strategy, each selection strategy has a subpopulation with sizeNs=1/3N. The decomposition-based strategy has a time complexity of OMNsT, where T is the neighborhood size. For the nondominated sorting strategy, the computational complexity is OMNs2. As for the nondominated neighbor-based strategy, the worst computational complexity is ONslgNs. So, the worst total computational complexity of HSS is
(11)OMNsT+OMNs2+ONslogNs.As for the CCE, in the worst case, the time complexity for calculating the distances isNQ2; the time complexity for adding the K+1 values is OKQ; the time complexity for removing K solutions from E is OKQ2.Based on the above analysis, the procedure of HSS does not increase any burden on the runtime complexity than using the single-selection strategy. As for the multiple-mutation strategies, it needs higher computational complexity, but it does not dramatically impose any serious burden on the runtime complexity than the single-mutation strategy-based operator. As for the CCE procedure, it needs much higher computational complexity than the complexity of crowding distance, which is onlyOMQlogQ.


## 4. Experimental Study on Test Functions

This section presents the numerical results that we have conducted. First, we describe a set of MOPs and state the quality indicators for evaluating DEHC’s performance. Then we evaluate DEHC by comparison with some other algorithms.


### 4.1. Test Problems and Performance Indicators

Twelve frequently used test MOPs from the literature are used here to evaluate the performance of DEHC. The first 5 are ZDT problems with two objectives and the next 7 are DTLZ problems with three objectives [22, 40]. Among all the test problems, ZDT1, ZDT2, ZDT3, and ZDT6 have 30 decision variables and the others have 10.In order to make a quantitative assessment on the performance of DEHC, two performance indicators are employed to evaluate the obtained approximation sets. The first is an inverted generational distance indicator (IGD) [41]. Given an approximation set A and a reference set R, IGD is defined as follows:
(12)IGDA,R=∑v∈Rdv,AR,where dv,A is the minimum Euclidean distance between v and the points in R. If R is large enough to represent the PF very well, IGDA,R could measure both the diversity and convergence of in a sense. To have a low value of IGDA,R, A must be very close to the PF and cannot miss any part of the whole PF.The second is a spread indicator (IS), which measures the distribution of an approximation set [22]. However, this indicator works well only for bi-objective problems and can not be applied directly to problems with more than two objectives. Inspired by the work in [34, 35, 42], an extended IS is used here
(13)IS=∑i=1Mdie+∑i=1N∑j=1M−1dij−d¯∑i=1Mdie+M−1Nd¯,where M is the objective number and die is the Euclidean distance between the ith extreme solutions in A and R and di,j is the Euclidean distance between the ith solution and its jth nearest solution in A. This one is a harmonic distance-based indicator because it takes all M−1 nearest neighbors around one solution into consideration. So, this works especially effectively for problems whose Pareto fronts consist of curved surfaces. A smaller value of this indicator means a better distribution. Particularly, a value of zero indicates that all the Pareto optima are equidistantly spaced.


### 4.2. Comparative Study

In order to assess DEHC, the final results are compared with those obtained by MOEA/D [18], NSGA-II [22], and NNIA [33] since they can be treated as algorithms with a single-selection strategy. Moreover, a recently proposed multiobjective DE (MODE) for dynamic optimization in chemical engineering is also selected as competing algorithm [43]. For all the 5 algorithms, the population sizes are set to 100 for the 2-objective and 300 for 3-objective problems, respectively. The maximal FFEs (fitness function evaluations) are set to 250.00 and 750.00 for the 2-objective and 3-objective problems, respectively. In DEHC, F=0.5 and CR=0.1. As for the other algorithms, their other parameters are kept the same as those in the corresponding references.For each test problem, 25 independent runs are conducted and the statistical values are analyzed in detail. Table1 shows the mean and standard deviation (Std) of IGD by 5 algorithms. It can be observed in Table 1 that DEHC is able to achieve the lowest values on most problems, which means that DEHC converges better than any other algorithms.Table 1
The IGD comparison results of DEHC with 4 MOEAs on 12 problems.
MOEA/D
NSGA-II
NNIA
MODE
DEHC
Mean ± Std
Mean ± Std
Mean ± Std
Mean ± Std
Mean ± Std
ZDT1
6.15E − 3 ± 4.40E − 3
5.19E − 3 ± 2.10E − 4
4.57E − 3 ± 2.68E − 4
3.56E − 2 ± 3.48 E − 3
3.87E − 3 ± 5.75E − 5
ZDT2
4.37E − 3 ± 1.39E − 3
5.16E − 3 ± 2.25E − 4
4.66E − 3 ± 2.42E − 4
6.78E − 2 ± 7.49E − 3
3.92E − 3 ± 4.55E − 5
ZDT3
1.71E − 2 ± 1.14E − 2
1.30E − 2 ± 1.29E − 2
5.19E − 3 ± 2.05E − 4
5.02E − 2 ± 5.22E − 3
4.61E − 3 ± 9.56E − 5
ZDT4
9.74E − 3 ± 3.71E − 3
6.36E − 3 ± 1.15E − 3
5.61E − 3 ± 1.44E − 3
3.88E − 1 ± 3.12E − 1
3.83E − 3 ± 5.70E − 5
ZDT6
5.52E − 3 ± 8.22E − 4
7.94E − 2 ± 6.84E − 3
3.64E − 3 ± 8.87E − 5
2.30E − 1 ± 2.02E − 2
3.054E − 3 ± 2.17 E − 5
DTLZ1
1.53E − 2 ± 1.58E − 4
1.79E − 2 ± 1.70E − 3
2.10E − 2 ± 2.22E − 2
6.99E − 1 ± 5.62E − 1
1.13E − 2 ± 7.71E − 5
DTLZ2
3.84E − 2 ± 8.70E − 5
4.00E − 2 ± 9.06E − 4
4.03E − 2 ± 8.82E − 4
6.34E − 2 ± 2.08E − 3
3.04E − 2 ± 1.39E − 4
DTLZ3
3.85E − 2 ± 2.88E − 4
4.31E − 2 ± 4.10E − 3
1.04E + 0 ± 4.99E + 0
3.45E + 0 ± 1.43E + 0
3.51E − 2 ± 1.51E − 3
DTLZ4
2.51E − 1 ± 3.05E − 1
3.89E − 2 ± 6.71E − 4
3.94E − 2 ± 8.05E − 4
6.48E − 2 ± 1.87E − 3
3.03E − 2 ± 1.06E − 4
DTLZ5
4.61E − 3 ± 1.21E − 5
1.78E − 3 ± 4.75E − 5
1.86E − 3 ± 8.18E − 5
5.48E − 3 ± 2.94E − 4
1.48E − 3 ± 2.12E − 5
DTLZ6
5.51E − 3 ± 3.68E − 4
4.50E − 2 ± 2.23E − 2
6.07E − 2 ± 2.37E − 2
3.44E + 0 ± 1.90E − 1
1.36E − 3 ± 6.58E − 6
DTLZ7
1.65E − 1 ± 1.49E − 1
5.36E − 2 ± 5.94E − 2
4.24E − 2 ± 1.56E − 3
7.72E − 2 ± 3.32E − 3
3.14E − 2 ± 2.96E − 4Table2 shows the mean and Std of IS by 5 algorithms. It can be observed in Table 2 that DEHC is able to obtain a better distributed Pareto optima set than any other algorithms on most problems. Although MOEA/D achieves better Std values on DTLZ1, DTLZ3 and NNIA achieve better Std values on DTLZ3 and DEHC achieves better mean values than MOEA/D and NNIA.Table 2
The IS comparison results of DEHC with 4 MOEAs on 12 problems.
MOEA/D
NSGA-II
NNIA
MODE
DEHC
Mean ± Std
Mean ± Std
Mean ± Std
Mean ± Std
Mean ± Std
ZDT1
3.78E − 1 ± 1.23E − 1
4.90E − 1 ± 3.20E − 2
3.52E − 1 ± 2.89E − 2
4.23E − 1 ± 1.34E − 1
1.43E − 1 ± 1.39E − 2
ZDT2
2.43E − 1 ± 1.68E − 1
4.62E − 1 ± 3.33E − 2
3.42E − 1 ± 3.21E − 2
5.81E − 1 ± 3.71E − 2
1.49E − 1 ± 1.13E − 2
ZDT3
8.86E − 1 ± 1.48E − 2
6.08E − 1 ± 3.97E − 2
5.30E − 1 ± 1.94E − 2
7.47E − 1 ± 6.26E − 2
4.50E − 1 ± 1.34E − 2
ZDT4
6.28E − 1 ± 3.27E − 1
3.60E − 1 ± 2.70E − 2
3.28E − 1 ± 2.56E − 2
1.33E + 0 ± 3.33E − 1
1.23E − 1 ± 1.18E − 2
ZDT6
1.98E − 1 ± 3.10E − 2
5.87E − 1 ± 2.92E − 2
3.23E − 1 ± 2.51E − 2
7.17E − 1 ± 2.65E − 1
1.37E − 1 ± 1.27E − 2
DTLZ1
6.14E − 1 ± 6.95E − 3
4.63E − 1 ± 3.54E − 2
4.92E − 1 ± 8.08E − 2
9.80E − 1 ± 4.17E − 1
8.65E − 2 ± 5.43E − 3
DTLZ2
6.47E − 1 ± 4.90E − 3
4.42E − 1 ± 2.54E − 2
4.55E − 1 ± 2.15E − 2
5.85E − 1 ± 2.41E − 2
9.50E − 2 ± 3.22E − 3
DTLZ3
6.57E − 1 ± 1.01E − 2
4.53E − 1 ± 3.24E − 2
4.81E − 1 ± 2.95E − 2
1.23E + 0 ± 3.24E − 1
1.94E − 1 ± 2.94E − 2
DTLZ4
7.82E − 1 ± 1.60E − 1
4.16E − 1 ± 2.03E − 2
4.34E − 1 ± 1.82E − 2
5.93E − 1 ± 2.53E − 2
9.54E − 2 ± 3.85E − 3
DTLZ5
1.15E + 0 ± 3.99E − 3
3.54E − 1 ± 1.24E − 2
3.68E − 1 ± 1.82E − 2
8.59E − 1 ± 1.01E − 2
1.66E − 1 ± 7.53E − 3
DTLZ6
8.43E − 1 ± 1.00E − 1
4.86E − 1 ± 3.79E − 2
5.35E − 1 ± 4.98E − 2
3.73E − 1 ± 2.63E − 2
1.59E − 1 ± 7.99E − 3
DTLZ7
1.04E+ 0 ± 4.34E − 2
5.15E − 1 ± 3.58E − 2
4.87E − 1 ± 2.47E − 2
4.65E − 1 ± 4.13E − 2
2.39E − 1 ± 8.57E-3In order to determine whether DEHC is significantly better than other algorithms, statistical significance tests ofIGD and IS are conducted between DEHC and MOEA/D, NSGA-II, NNIA, and MODE via Mann–Whitney rank sum test, respectively. Table 3 gives the comparison results. h=1 indicates a significant difference between the compared algorithms with a level of significance β=0.05, whereas h=0 indicates that the performance is not significantly different. It can be observed from the table that DEHC performs significantly better on almost all of the problems except for ZDT1 in terms of IGD with a level of significance β. As for ZDT1, the results do not show a significant difference between DEHC and NSGA-II in terms of IGD indicator. As for IS, the results show that DEHC performs better than other algorithms with a level of significance β.Table 3
The distribution ofIGD and IS using Mann–Whitney rank sum test.
(DEHC, MOEA/D)
(DEHC, NSGA-II)
(DEHC, NNIA)
(DEHC, MODE)
IGD
IS
IGD
IS
IGD
IS
IGD
IS
ZDT1
1
1
0
1
1
1
1
1
ZDT2
1
1
1
1
1
1
1
1
ZDT3
1
1
1
1
1
1
1
1
ZDT4
1
1
1
1
1
1
1
1
ZDT6
1
1
1
1
1
1
1
1
DTLZ1
1
1
1
1
1
1
1
1
DTLZ2
1
1
1
1
1
1
1
1
DTLZ3
1
1
1
1
1
1
1
1
DTLZ4
1
1
1
1
1
1
1
1
DTLZ5
1
1
1
1
1
1
1
1
DTLZ6
1
1
1
1
1
1
1
1
DTLZ7
1
1
1
1
1
1
1
1To provide more information about the convergence performance of the 5 MOEAs, we calculateIGD and IS values of found optimal solutions every 25 generations. Because in our study, all the 5 algorithms employ the same initialized population to avoid difference at the initialization stage, we do not calculate the indicators for the initial population but simply assume that they are at 1. Figures 4–7 illustrate the 2 indicators of the above 5 algorithms versus the evolution time over 25 runs in all the 12 instances.Figure 4
MedianIGD values versus the evolution time. (a) ZDT1, (b) ZDT2, (c) ZDT3, (d) ZDT4, (e) ZDT6, and (f) DTLZ1.
(a)
(b)
(c)
(d)
(e)
(f)Figure 5
MedianIGD values versus the evolution time. (a) DTLZ2, (b) DTLZ3, (c) DTLZ4, (d) DTLZ5, (e) DTLZ6, and (f) DTLZ7.
(a)
(b)
(c)
(d)
(e)
(f)Figure 6
MedianIS values versus the evolution time. (a) DTLZ2, (b) DTLZ3, (c) DTLZ4, (d) DTLZ5, (e) DTLZ6, and (f) DTLZ7.
(a)
(b)
(c)
(d)
(e)
(f)Figure 7
MedianIS values versus the evolution time. (a) DTLZ2, (b) DTLZ3, (c) DTLZ4, (d) DTLZ5, (e) DTLZ6, and (f) DTLZ7.
(a)
(b)
(c)
(d)
(e)
(f)For illustration, we also show some typical fronts obtained by the 5 algorithms. Figure8 presents the fronts on DTLZ1. Please note that these fronts correspond to the run with a median value with respect to IS. This figure clearly demonstrates that for DTLZ1, the 5 algorithms can visually converge to Pareto optimal fronts. However, the distributions have an obvious difference. DEHC provides the most uniformly distributed fronts among the four.Figure 8
The typical Pareto front obtained by 5 algorithms on DTLZ1. (a) MOEA/D, (b) NSGA-II, (c) NNIA, (d) MODE, and (e) DEHC.
(a)
(b)
(c)
(d)
(e)Figure9 presents the fronts on DTLZ6 and these fronts correspond to the run with a median value with respect to IGD. DTLZ6’s Pareto front is a line in three-dimensional space. Figure 9 demonstrates that NSGA-II and NNIA have some difficulties in converging to the optimal front and MOEA/D gets a point far from the real Pareto optimal front. Among all the 5 algorithms, DEHC seems to perform best on DTLZ6 because it converges to a curve exactly.Figure 9
The typical Pareto front obtained by 5 algorithms on DTLZ6. (a) MOEA/D, (b) NSGA-II, (c) NNIA, (d) MODE, and (e) DEHC.
(a)
(b)
(c)
(d)
(e)From the above comparisons and analyses, we can draw the conclusion that in terms ofIGD and IS, DEHC generates approximation sets with competitive diversity, uniformity, and proximity to the Pareto front. Our proposed DEHC can deal with MOPs with different types of Pareto front, that is, continuous, and disjoint, having 2 and 3 objectives.


### 4.3. Some Studies on Scalability of DEHC

To study how the performance of DEHC varies as the number of decision variable increases, we have tried on ZDT1 with different numbers of decision variables, that is, 30, 40, 50, 60, 70, 80, 90, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1500, 2000, 2500, and 3000. In these experiments, MOEA/D, NSGA-II, NNIA, and MODE are also conducted with the above numbers of decision variable. All the parameters are kept the same as in Section4.2. Figure 10 gives the mean IGD values versus the numbers of decision variables. We found that among the 25 independent runs for each number of decision variables, DEHC achieves lower IGD values on ZDT1 with different numbers of decision variable, especially from 60 to 1000. Furthermore, the mean IGD values almost linearly scale up as the number of decision variables increases. This experiment indicates that DEHC is able to handle large-scale MOPs efficiently and this ability may profit from HSS because a different selection strategy is effective for different problems during different stages of the search.Figure 10
MeanIGD on ZDT1 as the number of decision variable increases.


### 4.4. Discussion about CCE

In order to demonstrate an improved performance of CCE, we replace CCE in DEHC with other four truncating methods and then compare them with DEHC. These four methods are crowding distance [22], an improved version of crowding distance [44], crowding entropy [35], and harmonic distance [31]. We denote the algorithms with the four methods as DEHCD, DEHCD2, DEHCE, and DEHHD, respectively. For all the five algorithms, 25 independent runs are conducted to see the consistency of the algorithms.It can be observed from Table4 that for the five 2-objective problems, MOECD2 performs very well because MOECD2 achieves the best values for all problems except ZDT4. As for DEHC, it follows DEHCD2 because it has the best IS on ZDT4 and ranks 2nd for the other four 2-objective problems. As for the seven 3-objective problems, DEHC generates the most uniformly distributed Pareto front among all the five algorithms. From this experiment, we can deduce that CCE is able to maintain a uniformly distributed result especially for many-objective problems.Table 4
Comparison results of CCE and other truncating methods onIS.
DEHCD
DEHCD2
DEHCE
DEHHD
DEHC
Mean ± Std
Mean ± Std
Mean ± Std
Mean ± Std
Mean ± Std
ZDT1
2.774E − 1 ± 2.383E − 2
1.406E − 1 ± 1.199E − 2
3.623E − 1 ± 3.416E − 2
4.254E − 1 ± 3.093E − 2
1.425E − 1 ± 1.390E − 2
ZDT2
2.885E − 1 ± 2.474E − 2
1.420E − 1 ± 1.063E − 2
3.702E − 1 ± 3.298E − 2
4.273E − 1 ± 3.389E − 2
1.492E − 1 ± 1.134E − 2
ZDT3
5.021E − 1 ± 1.324E − 2
4.427E − 1 ± 1.064E − 2
5.565E − 1 ± 2.896E − 2
5.803E − 1 ± 3.071E − 2
4.497E − 1 ± 1.341E − 2
ZDT4
2.769E − 1 ± 2.869E − 2
1.241E − 1 ± 1.099E − 2
3.898E − 1 ± 3.178E − 2
4.487E − 1 ± 3.238E − 2
1.229E − 1 ± 1.175E − 2
ZDT6
3.049E − 1 ± 2.767E − 2
1.231E − 1 ± 1.030E − 2
3.841E − 1 ± 3.080E − 2
6.152E − 1 ± 2.915E − 1
1.373E − 1 ± 1.275E − 2
DTLZ1
3.860E − 1 ± 2.015E − 2
3.885E − 1 ± 1.938E − 2
3.986E − 1 ± 1.911E − 2
2.381E − 1 ± 2.997E − 2
8.648E − 2 ± 5.427E − 3
DTLZ2
3.604E − 1 ± 9.635E − 3
3.557E − 1 ± 1.431E − 2
3.682E − 1 ± 1.879E − 2
2.445E − 1 ± 8.816E − 3
9.496E − 2 ± 3.217E − 3
DTLZ3
3.707E − 1 ± 2.048E − 2
3.718E − 1 ± 1.773E − 2
3.732E − 1 ± 1.904E − 2
1.065E + 0 ± 6.659E − 1
1.939E − 1 ± 2.940E − 2
DTLZ4
3.578E − 1 ± 1.537E − 2
3.568E − 1 ± 1.741E − 2
3.612E − 1 ± 1.920E − 2
2.369E − 1 ± 1.005E − 2
9.538E − 2 ± 3.852E − 3
DTLZ5
3.505E − 1 ± 1.374E − 2
1.962E − 1 ± 5.629E − 3
3.929E − 1 ± 1.847E − 2
3.996E − 1 ± 2.107E − 2
1.657E − 1 ± 7.528E − 3
DTLZ6
3.258E − 1 ± 1.225E − 2
1.904E − 1 ± 8.794E − 3
3.936E − 1 ± 1.432E − 2
4.098E − 1 ± 1.460E − 2
1.592E − 1 ± 7.995E − 3
DTLZ7
4.400E − 1 ± 1.966E − 2
4.440E − 1 ± 2.472E − 2
4.468E − 1 ± 1.997E − 2
2.575E − 1 ± 1.166E − 2
2.395E − 1 ± 8.573sE − 3To confirm that CCE is really suited for MOPs with many objectives, here, we expand DTLZ2 to 5, 7, and 9 objectives and then studyIS. Figure 11 gives the mean IS on DTLZ2 with 3, 5, 7, and 9 objectives. It can be observed that for many objectives, CCE still works very well and it provides the most competitive results.Figure 11
MeanIS on DTLZ2 as the number of objectives increases.


## 4.1. Test Problems and Performance Indicators

Twelve frequently used test MOPs from the literature are used here to evaluate the performance of DEHC. The first 5 are ZDT problems with two objectives and the next 7 are DTLZ problems with three objectives [22, 40]. Among all the test problems, ZDT1, ZDT2, ZDT3, and ZDT6 have 30 decision variables and the others have 10.In order to make a quantitative assessment on the performance of DEHC, two performance indicators are employed to evaluate the obtained approximation sets. The first is an inverted generational distance indicator (IGD) [41]. Given an approximation set A and a reference set R, IGD is defined as follows:
(12)IGDA,R=∑v∈Rdv,AR,where dv,A is the minimum Euclidean distance between v and the points in R. If R is large enough to represent the PF very well, IGDA,R could measure both the diversity and convergence of in a sense. To have a low value of IGDA,R, A must be very close to the PF and cannot miss any part of the whole PF.The second is a spread indicator (IS), which measures the distribution of an approximation set [22]. However, this indicator works well only for bi-objective problems and can not be applied directly to problems with more than two objectives. Inspired by the work in [34, 35, 42], an extended IS is used here
(13)IS=∑i=1Mdie+∑i=1N∑j=1M−1dij−d¯∑i=1Mdie+M−1Nd¯,where M is the objective number and die is the Euclidean distance between the ith extreme solutions in A and R and di,j is the Euclidean distance between the ith solution and its jth nearest solution in A. This one is a harmonic distance-based indicator because it takes all M−1 nearest neighbors around one solution into consideration. So, this works especially effectively for problems whose Pareto fronts consist of curved surfaces. A smaller value of this indicator means a better distribution. Particularly, a value of zero indicates that all the Pareto optima are equidistantly spaced.


## 4.2. Comparative Study

In order to assess DEHC, the final results are compared with those obtained by MOEA/D [18], NSGA-II [22], and NNIA [33] since they can be treated as algorithms with a single-selection strategy. Moreover, a recently proposed multiobjective DE (MODE) for dynamic optimization in chemical engineering is also selected as competing algorithm [43]. For all the 5 algorithms, the population sizes are set to 100 for the 2-objective and 300 for 3-objective problems, respectively. The maximal FFEs (fitness function evaluations) are set to 250.00 and 750.00 for the 2-objective and 3-objective problems, respectively. In DEHC, F=0.5 and CR=0.1. As for the other algorithms, their other parameters are kept the same as those in the corresponding references.For each test problem, 25 independent runs are conducted and the statistical values are analyzed in detail. Table1 shows the mean and standard deviation (Std) of IGD by 5 algorithms. It can be observed in Table 1 that DEHC is able to achieve the lowest values on most problems, which means that DEHC converges better than any other algorithms.Table 1
The IGD comparison results of DEHC with 4 MOEAs on 12 problems.
MOEA/D
NSGA-II
NNIA
MODE
DEHC
Mean ± Std
Mean ± Std
Mean ± Std
Mean ± Std
Mean ± Std
ZDT1
6.15E − 3 ± 4.40E − 3
5.19E − 3 ± 2.10E − 4
4.57E − 3 ± 2.68E − 4
3.56E − 2 ± 3.48 E − 3
3.87E − 3 ± 5.75E − 5
ZDT2
4.37E − 3 ± 1.39E − 3
5.16E − 3 ± 2.25E − 4
4.66E − 3 ± 2.42E − 4
6.78E − 2 ± 7.49E − 3
3.92E − 3 ± 4.55E − 5
ZDT3
1.71E − 2 ± 1.14E − 2
1.30E − 2 ± 1.29E − 2
5.19E − 3 ± 2.05E − 4
5.02E − 2 ± 5.22E − 3
4.61E − 3 ± 9.56E − 5
ZDT4
9.74E − 3 ± 3.71E − 3
6.36E − 3 ± 1.15E − 3
5.61E − 3 ± 1.44E − 3
3.88E − 1 ± 3.12E − 1
3.83E − 3 ± 5.70E − 5
ZDT6
5.52E − 3 ± 8.22E − 4
7.94E − 2 ± 6.84E − 3
3.64E − 3 ± 8.87E − 5
2.30E − 1 ± 2.02E − 2
3.054E − 3 ± 2.17 E − 5
DTLZ1
1.53E − 2 ± 1.58E − 4
1.79E − 2 ± 1.70E − 3
2.10E − 2 ± 2.22E − 2
6.99E − 1 ± 5.62E − 1
1.13E − 2 ± 7.71E − 5
DTLZ2
3.84E − 2 ± 8.70E − 5
4.00E − 2 ± 9.06E − 4
4.03E − 2 ± 8.82E − 4
6.34E − 2 ± 2.08E − 3
3.04E − 2 ± 1.39E − 4
DTLZ3
3.85E − 2 ± 2.88E − 4
4.31E − 2 ± 4.10E − 3
1.04E + 0 ± 4.99E + 0
3.45E + 0 ± 1.43E + 0
3.51E − 2 ± 1.51E − 3
DTLZ4
2.51E − 1 ± 3.05E − 1
3.89E − 2 ± 6.71E − 4
3.94E − 2 ± 8.05E − 4
6.48E − 2 ± 1.87E − 3
3.03E − 2 ± 1.06E − 4
DTLZ5
4.61E − 3 ± 1.21E − 5
1.78E − 3 ± 4.75E − 5
1.86E − 3 ± 8.18E − 5
5.48E − 3 ± 2.94E − 4
1.48E − 3 ± 2.12E − 5
DTLZ6
5.51E − 3 ± 3.68E − 4
4.50E − 2 ± 2.23E − 2
6.07E − 2 ± 2.37E − 2
3.44E + 0 ± 1.90E − 1
1.36E − 3 ± 6.58E − 6
DTLZ7
1.65E − 1 ± 1.49E − 1
5.36E − 2 ± 5.94E − 2
4.24E − 2 ± 1.56E − 3
7.72E − 2 ± 3.32E − 3
3.14E − 2 ± 2.96E − 4Table2 shows the mean and Std of IS by 5 algorithms. It can be observed in Table 2 that DEHC is able to obtain a better distributed Pareto optima set than any other algorithms on most problems. Although MOEA/D achieves better Std values on DTLZ1, DTLZ3 and NNIA achieve better Std values on DTLZ3 and DEHC achieves better mean values than MOEA/D and NNIA.Table 2
The IS comparison results of DEHC with 4 MOEAs on 12 problems.
MOEA/D
NSGA-II
NNIA
MODE
DEHC
Mean ± Std
Mean ± Std
Mean ± Std
Mean ± Std
Mean ± Std
ZDT1
3.78E − 1 ± 1.23E − 1
4.90E − 1 ± 3.20E − 2
3.52E − 1 ± 2.89E − 2
4.23E − 1 ± 1.34E − 1
1.43E − 1 ± 1.39E − 2
ZDT2
2.43E − 1 ± 1.68E − 1
4.62E − 1 ± 3.33E − 2
3.42E − 1 ± 3.21E − 2
5.81E − 1 ± 3.71E − 2
1.49E − 1 ± 1.13E − 2
ZDT3
8.86E − 1 ± 1.48E − 2
6.08E − 1 ± 3.97E − 2
5.30E − 1 ± 1.94E − 2
7.47E − 1 ± 6.26E − 2
4.50E − 1 ± 1.34E − 2
ZDT4
6.28E − 1 ± 3.27E − 1
3.60E − 1 ± 2.70E − 2
3.28E − 1 ± 2.56E − 2
1.33E + 0 ± 3.33E − 1
1.23E − 1 ± 1.18E − 2
ZDT6
1.98E − 1 ± 3.10E − 2
5.87E − 1 ± 2.92E − 2
3.23E − 1 ± 2.51E − 2
7.17E − 1 ± 2.65E − 1
1.37E − 1 ± 1.27E − 2
DTLZ1
6.14E − 1 ± 6.95E − 3
4.63E − 1 ± 3.54E − 2
4.92E − 1 ± 8.08E − 2
9.80E − 1 ± 4.17E − 1
8.65E − 2 ± 5.43E − 3
DTLZ2
6.47E − 1 ± 4.90E − 3
4.42E − 1 ± 2.54E − 2
4.55E − 1 ± 2.15E − 2
5.85E − 1 ± 2.41E − 2
9.50E − 2 ± 3.22E − 3
DTLZ3
6.57E − 1 ± 1.01E − 2
4.53E − 1 ± 3.24E − 2
4.81E − 1 ± 2.95E − 2
1.23E + 0 ± 3.24E − 1
1.94E − 1 ± 2.94E − 2
DTLZ4
7.82E − 1 ± 1.60E − 1
4.16E − 1 ± 2.03E − 2
4.34E − 1 ± 1.82E − 2
5.93E − 1 ± 2.53E − 2
9.54E − 2 ± 3.85E − 3
DTLZ5
1.15E + 0 ± 3.99E − 3
3.54E − 1 ± 1.24E − 2
3.68E − 1 ± 1.82E − 2
8.59E − 1 ± 1.01E − 2
1.66E − 1 ± 7.53E − 3
DTLZ6
8.43E − 1 ± 1.00E − 1
4.86E − 1 ± 3.79E − 2
5.35E − 1 ± 4.98E − 2
3.73E − 1 ± 2.63E − 2
1.59E − 1 ± 7.99E − 3
DTLZ7
1.04E+ 0 ± 4.34E − 2
5.15E − 1 ± 3.58E − 2
4.87E − 1 ± 2.47E − 2
4.65E − 1 ± 4.13E − 2
2.39E − 1 ± 8.57E-3In order to determine whether DEHC is significantly better than other algorithms, statistical significance tests ofIGD and IS are conducted between DEHC and MOEA/D, NSGA-II, NNIA, and MODE via Mann–Whitney rank sum test, respectively. Table 3 gives the comparison results. h=1 indicates a significant difference between the compared algorithms with a level of significance β=0.05, whereas h=0 indicates that the performance is not significantly different. It can be observed from the table that DEHC performs significantly better on almost all of the problems except for ZDT1 in terms of IGD with a level of significance β. As for ZDT1, the results do not show a significant difference between DEHC and NSGA-II in terms of IGD indicator. As for IS, the results show that DEHC performs better than other algorithms with a level of significance β.Table 3
The distribution ofIGD and IS using Mann–Whitney rank sum test.
(DEHC, MOEA/D)
(DEHC, NSGA-II)
(DEHC, NNIA)
(DEHC, MODE)
IGD
IS
IGD
IS
IGD
IS
IGD
IS
ZDT1
1
1
0
1
1
1
1
1
ZDT2
1
1
1
1
1
1
1
1
ZDT3
1
1
1
1
1
1
1
1
ZDT4
1
1
1
1
1
1
1
1
ZDT6
1
1
1
1
1
1
1
1
DTLZ1
1
1
1
1
1
1
1
1
DTLZ2
1
1
1
1
1
1
1
1
DTLZ3
1
1
1
1
1
1
1
1
DTLZ4
1
1
1
1
1
1
1
1
DTLZ5
1
1
1
1
1
1
1
1
DTLZ6
1
1
1
1
1
1
1
1
DTLZ7
1
1
1
1
1
1
1
1To provide more information about the convergence performance of the 5 MOEAs, we calculateIGD and IS values of found optimal solutions every 25 generations. Because in our study, all the 5 algorithms employ the same initialized population to avoid difference at the initialization stage, we do not calculate the indicators for the initial population but simply assume that they are at 1. Figures 4–7 illustrate the 2 indicators of the above 5 algorithms versus the evolution time over 25 runs in all the 12 instances.Figure 4
MedianIGD values versus the evolution time. (a) ZDT1, (b) ZDT2, (c) ZDT3, (d) ZDT4, (e) ZDT6, and (f) DTLZ1.
(a)
(b)
(c)
(d)
(e)
(f)Figure 5
MedianIGD values versus the evolution time. (a) DTLZ2, (b) DTLZ3, (c) DTLZ4, (d) DTLZ5, (e) DTLZ6, and (f) DTLZ7.
(a)
(b)
(c)
(d)
(e)
(f)Figure 6
MedianIS values versus the evolution time. (a) DTLZ2, (b) DTLZ3, (c) DTLZ4, (d) DTLZ5, (e) DTLZ6, and (f) DTLZ7.
(a)
(b)
(c)
(d)
(e)
(f)Figure 7
MedianIS values versus the evolution time. (a) DTLZ2, (b) DTLZ3, (c) DTLZ4, (d) DTLZ5, (e) DTLZ6, and (f) DTLZ7.
(a)
(b)
(c)
(d)
(e)
(f)For illustration, we also show some typical fronts obtained by the 5 algorithms. Figure8 presents the fronts on DTLZ1. Please note that these fronts correspond to the run with a median value with respect to IS. This figure clearly demonstrates that for DTLZ1, the 5 algorithms can visually converge to Pareto optimal fronts. However, the distributions have an obvious difference. DEHC provides the most uniformly distributed fronts among the four.Figure 8
The typical Pareto front obtained by 5 algorithms on DTLZ1. (a) MOEA/D, (b) NSGA-II, (c) NNIA, (d) MODE, and (e) DEHC.
(a)
(b)
(c)
(d)
(e)Figure9 presents the fronts on DTLZ6 and these fronts correspond to the run with a median value with respect to IGD. DTLZ6’s Pareto front is a line in three-dimensional space. Figure 9 demonstrates that NSGA-II and NNIA have some difficulties in converging to the optimal front and MOEA/D gets a point far from the real Pareto optimal front. Among all the 5 algorithms, DEHC seems to perform best on DTLZ6 because it converges to a curve exactly.Figure 9
The typical Pareto front obtained by 5 algorithms on DTLZ6. (a) MOEA/D, (b) NSGA-II, (c) NNIA, (d) MODE, and (e) DEHC.
(a)
(b)
(c)
(d)
(e)From the above comparisons and analyses, we can draw the conclusion that in terms ofIGD and IS, DEHC generates approximation sets with competitive diversity, uniformity, and proximity to the Pareto front. Our proposed DEHC can deal with MOPs with different types of Pareto front, that is, continuous, and disjoint, having 2 and 3 objectives.


## 4.3. Some Studies on Scalability of DEHC

To study how the performance of DEHC varies as the number of decision variable increases, we have tried on ZDT1 with different numbers of decision variables, that is, 30, 40, 50, 60, 70, 80, 90, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1500, 2000, 2500, and 3000. In these experiments, MOEA/D, NSGA-II, NNIA, and MODE are also conducted with the above numbers of decision variable. All the parameters are kept the same as in Section4.2. Figure 10 gives the mean IGD values versus the numbers of decision variables. We found that among the 25 independent runs for each number of decision variables, DEHC achieves lower IGD values on ZDT1 with different numbers of decision variable, especially from 60 to 1000. Furthermore, the mean IGD values almost linearly scale up as the number of decision variables increases. This experiment indicates that DEHC is able to handle large-scale MOPs efficiently and this ability may profit from HSS because a different selection strategy is effective for different problems during different stages of the search.Figure 10
MeanIGD on ZDT1 as the number of decision variable increases.


## 4.4. Discussion about CCE

In order to demonstrate an improved performance of CCE, we replace CCE in DEHC with other four truncating methods and then compare them with DEHC. These four methods are crowding distance [22], an improved version of crowding distance [44], crowding entropy [35], and harmonic distance [31]. We denote the algorithms with the four methods as DEHCD, DEHCD2, DEHCE, and DEHHD, respectively. For all the five algorithms, 25 independent runs are conducted to see the consistency of the algorithms.It can be observed from Table4 that for the five 2-objective problems, MOECD2 performs very well because MOECD2 achieves the best values for all problems except ZDT4. As for DEHC, it follows DEHCD2 because it has the best IS on ZDT4 and ranks 2nd for the other four 2-objective problems. As for the seven 3-objective problems, DEHC generates the most uniformly distributed Pareto front among all the five algorithms. From this experiment, we can deduce that CCE is able to maintain a uniformly distributed result especially for many-objective problems.Table 4
Comparison results of CCE and other truncating methods onIS.
DEHCD
DEHCD2
DEHCE
DEHHD
DEHC
Mean ± Std
Mean ± Std
Mean ± Std
Mean ± Std
Mean ± Std
ZDT1
2.774E − 1 ± 2.383E − 2
1.406E − 1 ± 1.199E − 2
3.623E − 1 ± 3.416E − 2
4.254E − 1 ± 3.093E − 2
1.425E − 1 ± 1.390E − 2
ZDT2
2.885E − 1 ± 2.474E − 2
1.420E − 1 ± 1.063E − 2
3.702E − 1 ± 3.298E − 2
4.273E − 1 ± 3.389E − 2
1.492E − 1 ± 1.134E − 2
ZDT3
5.021E − 1 ± 1.324E − 2
4.427E − 1 ± 1.064E − 2
5.565E − 1 ± 2.896E − 2
5.803E − 1 ± 3.071E − 2
4.497E − 1 ± 1.341E − 2
ZDT4
2.769E − 1 ± 2.869E − 2
1.241E − 1 ± 1.099E − 2
3.898E − 1 ± 3.178E − 2
4.487E − 1 ± 3.238E − 2
1.229E − 1 ± 1.175E − 2
ZDT6
3.049E − 1 ± 2.767E − 2
1.231E − 1 ± 1.030E − 2
3.841E − 1 ± 3.080E − 2
6.152E − 1 ± 2.915E − 1
1.373E − 1 ± 1.275E − 2
DTLZ1
3.860E − 1 ± 2.015E − 2
3.885E − 1 ± 1.938E − 2
3.986E − 1 ± 1.911E − 2
2.381E − 1 ± 2.997E − 2
8.648E − 2 ± 5.427E − 3
DTLZ2
3.604E − 1 ± 9.635E − 3
3.557E − 1 ± 1.431E − 2
3.682E − 1 ± 1.879E − 2
2.445E − 1 ± 8.816E − 3
9.496E − 2 ± 3.217E − 3
DTLZ3
3.707E − 1 ± 2.048E − 2
3.718E − 1 ± 1.773E − 2
3.732E − 1 ± 1.904E − 2
1.065E + 0 ± 6.659E − 1
1.939E − 1 ± 2.940E − 2
DTLZ4
3.578E − 1 ± 1.537E − 2
3.568E − 1 ± 1.741E − 2
3.612E − 1 ± 1.920E − 2
2.369E − 1 ± 1.005E − 2
9.538E − 2 ± 3.852E − 3
DTLZ5
3.505E − 1 ± 1.374E − 2
1.962E − 1 ± 5.629E − 3
3.929E − 1 ± 1.847E − 2
3.996E − 1 ± 2.107E − 2
1.657E − 1 ± 7.528E − 3
DTLZ6
3.258E − 1 ± 1.225E − 2
1.904E − 1 ± 8.794E − 3
3.936E − 1 ± 1.432E − 2
4.098E − 1 ± 1.460E − 2
1.592E − 1 ± 7.995E − 3
DTLZ7
4.400E − 1 ± 1.966E − 2
4.440E − 1 ± 2.472E − 2
4.468E − 1 ± 1.997E − 2
2.575E − 1 ± 1.166E − 2
2.395E − 1 ± 8.573sE − 3To confirm that CCE is really suited for MOPs with many objectives, here, we expand DTLZ2 to 5, 7, and 9 objectives and then studyIS. Figure 11 gives the mean IS on DTLZ2 with 3, 5, 7, and 9 objectives. It can be observed that for many objectives, CCE still works very well and it provides the most competitive results.Figure 11
MeanIS on DTLZ2 as the number of objectives increases.


## 5. Applications to Solve MOOCPs

In this section, DEHC is applied to solve three MOOCPs in chemical processes taken from the literature. To apply DEHC for MOOCPs, the MOOCP is first transformed into common MOPs by the commonly used CVP approach. For the CVP method, the time ranget0,tf is divided into T equal segments, that is, t0,tf=t0,t1×t1,t2×⋯×tT−1,tT, tT=tf, and the control variables can be approximated by basic functions, such as piece-wise constant function and piece-wise linear functions [43, 45]. Here, we assume ut as piecewise constant function. As a result, T×K parameters determine ut over t0,tf. In terms of the jth segments; the control variables over the entire time span are as follows:
(14)ckt=uk,t∈tj−1,tj,0,otherwise,j=1,2,…,T.So, the total control variables over the entire time span can be described as the summation below:(15)ut=∑j=1Tcjt,t∈tj−1,tj.Ones the control variables are discretized, DEHC can be applied to find the best control profile of the discrete time system as an approximation of the continuous problem. So, the major steps to optimize the MOOCP are outlined as the following steps:Step 1.
Divide the time ranget0,tf into T equal segments and denote the jth segment as tj−1,tj, where j=1,2,…,K and tT=tf.Step 2.
Assume that the value of theith input variable within the jth segment is uij, where i=1,2,…,m and j=1,2,…,T.Step 3.
Once the value ofut is determined within the time range t0,tf, a numerical integral method, such as the standard Runge-Kutta method, can be employed to obtain the objective value, that is, Juv.Step 4.
Apply DEHC to obtain the optimaluij and Jut as an approximation of the continuous problem, where i=1,2,…,K and j=1,2,…,T.


### 5.1. Case I: Catalyst Mixing Problem in a Tubular Reactor

This problem considers a steady-state plug flow reactor of fixed lengthtf. The reactor is packed with two catalysts which are required to stimulate a series of reactions (one reversible and one irreversible S1↔S2→S3). These assumptions give rise to the following model [46]:
(16)dx1dt=u10x2−x1,dx2dt=ux1−10x2−1−ux2,xt0=1,0T,0≤u≤1,tf=12,where x1 and x2 are the concentrations of S1 and S2, u is the fraction of catalyst A, and t is the spatial coordinate. The objective is to determine the optimal mixing policy of the two catalysts in order to maximize the production S3 at the reactor outlet and minimize the amount of the most expensive catalyst A.
(17)maxJ1=1−x1tf−x2tf,minJ2=∫0tfudt.For this problem, we setT=10, the population size is N=50, and the maximum generation is 100. The obtained Pareto optimal fronts by DEHC and MODE are illustrated in Figure 12. Figure 12 shows that both the two algorithms achieve a wide variety of solutions which also uniformly spread along the Pareto optimal front, but DEHC provides a more uniformly and widely distributed result along the front.Figure 12
The Pareto optimal front for case I.When considering maximizing the production of S3 only, Vassiliadis et al. obtained a best value of 0.0480238 by a single-objective optimization method [46] and Chen et al. obtained a best value of 0.04798 by a multiobjective DE with a ranking-based mutation operator [9]. The result of our obtained Pareto front is 0.04800, which is very close to 0.04798. To assist decision making, five Pareto optimal reaction temperatures, that is, u1, u2, u3, u4, and u5 in Figure 12, with equal distribution along the front are shown in Table 5 and the corresponding Pareto optimal trajectories of these are illustrated in Figure 13. It can be observed that the trajectories of mixing policy of the catalysts in the entire interval gradually decrease from 1 to 0, that is, from u1 to u5, leading to the reduction of the use of catalyst A.Table 5
Five solutions selected from the front for case I.
Segment
u1
u2
u3
u4
u5
0.0–1.2
1.0000
1.0000
0.7555
0.2663
0.0000
1.2–2.4
0.5477
0.3310
0.0000
0.0392
0.0000
2.4–3.6
0.1534
0.1307
0.0447
0.0143
0.0000
3.6–4.8
0.2784
0.0953
0.0096
0.0191
0.0000
4.8–6.0
0.2354
0.0000
0.0229
0.0214
0.0000
6.0–7.2
0.2482
0.0840
0.0000
0.0000
0.0000
7.2–8.4
0.1683
0.0000
0.0000
0.0038
0.0000
8.4–9.6
0.0657
0.0043
0.0000
0.0000
0.0000
9.6–10.8
0.0011
0.0000
0.0003
0.0000
0.0000
10.8–12.0
0.0000
0.0000
0.0000
0.0000
0.0000Figure 13
The Pareto optimal trajectories of five selected solutions for case I. (a) u1, (b) u2, (c) u3, (d) u4, and (e) u5.
(a)
(b)
(c)
(d)
(e)


### 5.2. Case II: Optimal Operation of a Fed Batch Reactor

The third case is a batch reactor based on the two reaction systems [47]:
(18)A+B→k1C,B+B→k2D.This process can be described by the following mechanistic model:(19)dAdt=−k1AB−AVu,dBdt=−k1AB−2k2B2+bfeed−BVu,dCdt=k1AB−CVu,dDdt=2k2B−DVu,dVdt=u,0≤u≤0.01,tf=120min,where A, B, C, and D are the concentrations of A, B, C, and D, respectively; V is the current reaction volume, u is the reactant feed rate, and k1=0.5 and k2=0.5 are the reaction rate constants. At the start of reaction, the reactor contains A0=0.2 mole/L of A; no BB0=0 and is fed to 50% (V0=0.5 m3). The objective is to maximize the amount of product C and minimize the amount of by-product D.For this problem, we setT=10, the population size is N=50, and the maximum generation is 100. The obtained Pareto fronts by DEHC and MODE are showed in Figure 14. It can be observed form Figure 14 that both DEHC and MODE achieve very high-quality results in approximating the true Pareto front. However, DEHC achieves a higher distribution behavior.Figure 14
The Pareto optimal front for case II.The obtained Pareto optimal front is illustrated in Figure14. When the complete focus is put on maximizing the amount of product C, Chen et al. obtained a best value of 0.07936 and we obtained 0.08079 [9]. To assist decision making, five Pareto optimal reaction temperatures along the front in Figure 14 with equal distribution are shown in Table 6 and the corresponding Pareto optimal trajectories of these are illustrated in Figure 15. It can be observed from the multiple trajectories that as the reduction of the amount of by-product D, the feed rate in the entire time interval gradually reduces to the low value.
(20)maxJ1=CtfVtf,minJ2=DtfVtf.Table 6
Five solutions selected from the front for case II.
Segment
u1
u2
u3
u4
u5
0.0–12.0
0.0100
0.0100
0.0096
0.0031
0.0000
12.0–24.0
0.0100
0.0100
0.0029
0.0017
0.0000
24.0–36.0
0.0100
0.0075
0.0044
0.0001
0.0000
36.0–48.0
0.0100
0.0073
0.0029
0.0027
0.0000
48.0–60.0
0.0100
0.0099
0.0045
0.0006
0.0000
60.0–72.0
0.0100
0.0029
0.0041
0.0062
0.0000
72.0–84.0
0.0100
0.0094
0.0031
0.0000
0.0000
84.0–96.0
0.0100
0.0070
0.0036
0.0020
0.0000
96.0–108.0
0.0100
0.0071
0.0072
0.0011
0.0000
108.0–120.0
0.0100
0.0002
0.0009
0.0000
0.0000Figure 15
The Pareto optimal trajectories of five selected solutions for case II. (a) u1, (b) u2, (c) u3, (d) u4, and (e) u5.
(a)
(b)
(c)
(d)
(e)


### 5.3. Case III: Foreign Protein Production Using Recombinant Bacteria

The fourth problem is a model of foreign protein production using recombinant bacteria [48]. In this case, there are 7 state variables which are reaction volume (x1, L), cell density (x2, g/L), nutrient concentration (x3, g/L), foreign protein concentration (x4, g/L), inducer concentration (x5, g/L), inducer shock factor on the cell growth rate (x6), and the inducer recovery factor on the cell growth rate (x7). The two control variables are glucose feed rate (u1, L/h) and inducer feed rate (u2, L/h). The model was described as follows [49]:
(21)dx1dt=u1+u2,dx2dt=ux2+λx2,dx3dt=Cnfu1x1+λx3−μYx2,dx4dt=Rfpx2+λx4,dx5dt=Cifu2x1+λx5,dx6dt=−k1x6,dx7dt=k21−x7,xt0=1,0.1,40,0,0,1,0T,0≤u1≤10−2,0≤u2≤10−2,t0=0,tf=10h,where
(22)μ=0.407x30.108+x3+x32/14814.8×x6+0.22x70.22+x5,Rfp=0.095x30.108+x3+x32/14814.8×0.0005+x50.022+x5,λ=−u1+u2x1,k1=k2=0.09x50.034+x5,Cnf=100.0,Cif=4.0,Y=0.51.As the price of the inducer, that is,Pr, is expensive, here, the two objectives are set as maximizing foreign protein production and minimizing the consumption of the inducer by controlling the glucose feed rate and the inducer feed rate, which are modeled as follows:
(23)maxJ1=x1tfx4tf,minJ1=∫0tfu2dt.For this problem, we setT=10, the population size is N=50, and the maximum generation is 100. The obtained Pareto fronts by DEHC and MODE are showed in Figure 16. It can be observed form Figure 16 that both DEHC and MODE generate very close approximations to the Pareto front. But the attained Pareto front with MODE cannot spread along the front as much as DEHC.Figure 16
The Pareto optimal front for case III.To assist decision making, five Pareto optimal reaction temperatures along the front in Figure16 with equal distribution are shown in Table 7 and the corresponding Pareto optimal trajectories of these are illustrated in Figure 17.
(24)maxJ=x1tfx4tf−Pr∫t0tfu2tdt.Table 7
Five solutions selected from the front for case III.
u1
u2
u3
u4
u5
u1(1)
u1(2)
u2(1)
u2(2)
u3(1)
u3(2)
u4(1)
u4(2)
u5(1)
u5(2)
0.0-1.0
0.0000
0.0000
0.0000
0.0000
0.0072
0.0000
0.0079
0.0002
0.0040
0.0000
1.0–2.0
0.0001
0.0000
0.0053
0.0000
0.0003
0.0000
0.0002
0.0017
0.0057
0.0000
2.0–3.0
0.0039
0.0006
0.0050
0.0045
0.0002
0.0006
0.0100
0.0000
0.0098
0.0000
3.0–4.0
0.0000
0.0100
0.0011
0.0095
0.0000
0.0100
0.0000
0.0008
0.0000
0.0000
4.0–5.0
0.0012
0.0100
0.0001
0.0100
0.0063
0.0000
0.0025
0.0000
0.0094
0.0000
5.0–6.0
0.0000
0.0100
0.0020
0.0100
0.0021
0.0000
0.0042
0.0000
0.0044
0.0000
6.0–7.0
0.0004
0.0100
0.0000
0.0000
0.0045
0.0002
0.0068
0.0000
0.0055
0.0000
7.0–8.0
0.0002
0.0100
0.0061
0.0001
0.0002
0.0000
0.0060
0.0001
0.0095
0.0000
8.0–9.0
0.0003
0.0100
0.0058
0.0000
0.0002
0.0000
0.0056
0.0000
0.0094
0.0000
9.0–10.0
0.0013
0.0096
0.0055
0.0000
0.0055
0.0000
0.0096
0.0000
0.0064
0.0000Figure 17
The Pareto optimal trajectories of five selected solutions for case III. (a) u1, (b) u2, (c) u3, (d) u4, and (e) u5.
(a)
(b)
(c)
(d)
(e)Roubos et al. setPr=5 and solve (24) using GA with a population size of 40 and generation of 25000 and obtained an optimum of 0.8149 [50]. We calculate the combined objective values with Pr=5 and obtain an optimal value of 0.81345. Figure 18 also shows the combined objective values along the Pareto front by setting Pr=5.Figure 18
The objective value of the Pareto optimal front for case III.


## 5.1. Case I: Catalyst Mixing Problem in a Tubular Reactor

This problem considers a steady-state plug flow reactor of fixed lengthtf. The reactor is packed with two catalysts which are required to stimulate a series of reactions (one reversible and one irreversible S1↔S2→S3). These assumptions give rise to the following model [46]:
(16)dx1dt=u10x2−x1,dx2dt=ux1−10x2−1−ux2,xt0=1,0T,0≤u≤1,tf=12,where x1 and x2 are the concentrations of S1 and S2, u is the fraction of catalyst A, and t is the spatial coordinate. The objective is to determine the optimal mixing policy of the two catalysts in order to maximize the production S3 at the reactor outlet and minimize the amount of the most expensive catalyst A.
(17)maxJ1=1−x1tf−x2tf,minJ2=∫0tfudt.For this problem, we setT=10, the population size is N=50, and the maximum generation is 100. The obtained Pareto optimal fronts by DEHC and MODE are illustrated in Figure 12. Figure 12 shows that both the two algorithms achieve a wide variety of solutions which also uniformly spread along the Pareto optimal front, but DEHC provides a more uniformly and widely distributed result along the front.Figure 12
The Pareto optimal front for case I.When considering maximizing the production of S3 only, Vassiliadis et al. obtained a best value of 0.0480238 by a single-objective optimization method [46] and Chen et al. obtained a best value of 0.04798 by a multiobjective DE with a ranking-based mutation operator [9]. The result of our obtained Pareto front is 0.04800, which is very close to 0.04798. To assist decision making, five Pareto optimal reaction temperatures, that is, u1, u2, u3, u4, and u5 in Figure 12, with equal distribution along the front are shown in Table 5 and the corresponding Pareto optimal trajectories of these are illustrated in Figure 13. It can be observed that the trajectories of mixing policy of the catalysts in the entire interval gradually decrease from 1 to 0, that is, from u1 to u5, leading to the reduction of the use of catalyst A.Table 5
Five solutions selected from the front for case I.
Segment
u1
u2
u3
u4
u5
0.0–1.2
1.0000
1.0000
0.7555
0.2663
0.0000
1.2–2.4
0.5477
0.3310
0.0000
0.0392
0.0000
2.4–3.6
0.1534
0.1307
0.0447
0.0143
0.0000
3.6–4.8
0.2784
0.0953
0.0096
0.0191
0.0000
4.8–6.0
0.2354
0.0000
0.0229
0.0214
0.0000
6.0–7.2
0.2482
0.0840
0.0000
0.0000
0.0000
7.2–8.4
0.1683
0.0000
0.0000
0.0038
0.0000
8.4–9.6
0.0657
0.0043
0.0000
0.0000
0.0000
9.6–10.8
0.0011
0.0000
0.0003
0.0000
0.0000
10.8–12.0
0.0000
0.0000
0.0000
0.0000
0.0000Figure 13
The Pareto optimal trajectories of five selected solutions for case I. (a) u1, (b) u2, (c) u3, (d) u4, and (e) u5.
(a)
(b)
(c)
(d)
(e)


## 5.2. Case II: Optimal Operation of a Fed Batch Reactor

The third case is a batch reactor based on the two reaction systems [47]:
(18)A+B→k1C,B+B→k2D.This process can be described by the following mechanistic model:(19)dAdt=−k1AB−AVu,dBdt=−k1AB−2k2B2+bfeed−BVu,dCdt=k1AB−CVu,dDdt=2k2B−DVu,dVdt=u,0≤u≤0.01,tf=120min,where A, B, C, and D are the concentrations of A, B, C, and D, respectively; V is the current reaction volume, u is the reactant feed rate, and k1=0.5 and k2=0.5 are the reaction rate constants. At the start of reaction, the reactor contains A0=0.2 mole/L of A; no BB0=0 and is fed to 50% (V0=0.5 m3). The objective is to maximize the amount of product C and minimize the amount of by-product D.For this problem, we setT=10, the population size is N=50, and the maximum generation is 100. The obtained Pareto fronts by DEHC and MODE are showed in Figure 14. It can be observed form Figure 14 that both DEHC and MODE achieve very high-quality results in approximating the true Pareto front. However, DEHC achieves a higher distribution behavior.Figure 14
The Pareto optimal front for case II.The obtained Pareto optimal front is illustrated in Figure14. When the complete focus is put on maximizing the amount of product C, Chen et al. obtained a best value of 0.07936 and we obtained 0.08079 [9]. To assist decision making, five Pareto optimal reaction temperatures along the front in Figure 14 with equal distribution are shown in Table 6 and the corresponding Pareto optimal trajectories of these are illustrated in Figure 15. It can be observed from the multiple trajectories that as the reduction of the amount of by-product D, the feed rate in the entire time interval gradually reduces to the low value.
(20)maxJ1=CtfVtf,minJ2=DtfVtf.Table 6
Five solutions selected from the front for case II.
Segment
u1
u2
u3
u4
u5
0.0–12.0
0.0100
0.0100
0.0096
0.0031
0.0000
12.0–24.0
0.0100
0.0100
0.0029
0.0017
0.0000
24.0–36.0
0.0100
0.0075
0.0044
0.0001
0.0000
36.0–48.0
0.0100
0.0073
0.0029
0.0027
0.0000
48.0–60.0
0.0100
0.0099
0.0045
0.0006
0.0000
60.0–72.0
0.0100
0.0029
0.0041
0.0062
0.0000
72.0–84.0
0.0100
0.0094
0.0031
0.0000
0.0000
84.0–96.0
0.0100
0.0070
0.0036
0.0020
0.0000
96.0–108.0
0.0100
0.0071
0.0072
0.0011
0.0000
108.0–120.0
0.0100
0.0002
0.0009
0.0000
0.0000Figure 15
The Pareto optimal trajectories of five selected solutions for case II. (a) u1, (b) u2, (c) u3, (d) u4, and (e) u5.
(a)
(b)
(c)
(d)
(e)


## 5.3. Case III: Foreign Protein Production Using Recombinant Bacteria

The fourth problem is a model of foreign protein production using recombinant bacteria [48]. In this case, there are 7 state variables which are reaction volume (x1, L), cell density (x2, g/L), nutrient concentration (x3, g/L), foreign protein concentration (x4, g/L), inducer concentration (x5, g/L), inducer shock factor on the cell growth rate (x6), and the inducer recovery factor on the cell growth rate (x7). The two control variables are glucose feed rate (u1, L/h) and inducer feed rate (u2, L/h). The model was described as follows [49]:
(21)dx1dt=u1+u2,dx2dt=ux2+λx2,dx3dt=Cnfu1x1+λx3−μYx2,dx4dt=Rfpx2+λx4,dx5dt=Cifu2x1+λx5,dx6dt=−k1x6,dx7dt=k21−x7,xt0=1,0.1,40,0,0,1,0T,0≤u1≤10−2,0≤u2≤10−2,t0=0,tf=10h,where
(22)μ=0.407x30.108+x3+x32/14814.8×x6+0.22x70.22+x5,Rfp=0.095x30.108+x3+x32/14814.8×0.0005+x50.022+x5,λ=−u1+u2x1,k1=k2=0.09x50.034+x5,Cnf=100.0,Cif=4.0,Y=0.51.As the price of the inducer, that is,Pr, is expensive, here, the two objectives are set as maximizing foreign protein production and minimizing the consumption of the inducer by controlling the glucose feed rate and the inducer feed rate, which are modeled as follows:
(23)maxJ1=x1tfx4tf,minJ1=∫0tfu2dt.For this problem, we setT=10, the population size is N=50, and the maximum generation is 100. The obtained Pareto fronts by DEHC and MODE are showed in Figure 16. It can be observed form Figure 16 that both DEHC and MODE generate very close approximations to the Pareto front. But the attained Pareto front with MODE cannot spread along the front as much as DEHC.Figure 16
The Pareto optimal front for case III.To assist decision making, five Pareto optimal reaction temperatures along the front in Figure16 with equal distribution are shown in Table 7 and the corresponding Pareto optimal trajectories of these are illustrated in Figure 17.
(24)maxJ=x1tfx4tf−Pr∫t0tfu2tdt.Table 7
Five solutions selected from the front for case III.
u1
u2
u3
u4
u5
u1(1)
u1(2)
u2(1)
u2(2)
u3(1)
u3(2)
u4(1)
u4(2)
u5(1)
u5(2)
0.0-1.0
0.0000
0.0000
0.0000
0.0000
0.0072
0.0000
0.0079
0.0002
0.0040
0.0000
1.0–2.0
0.0001
0.0000
0.0053
0.0000
0.0003
0.0000
0.0002
0.0017
0.0057
0.0000
2.0–3.0
0.0039
0.0006
0.0050
0.0045
0.0002
0.0006
0.0100
0.0000
0.0098
0.0000
3.0–4.0
0.0000
0.0100
0.0011
0.0095
0.0000
0.0100
0.0000
0.0008
0.0000
0.0000
4.0–5.0
0.0012
0.0100
0.0001
0.0100
0.0063
0.0000
0.0025
0.0000
0.0094
0.0000
5.0–6.0
0.0000
0.0100
0.0020
0.0100
0.0021
0.0000
0.0042
0.0000
0.0044
0.0000
6.0–7.0
0.0004
0.0100
0.0000
0.0000
0.0045
0.0002
0.0068
0.0000
0.0055
0.0000
7.0–8.0
0.0002
0.0100
0.0061
0.0001
0.0002
0.0000
0.0060
0.0001
0.0095
0.0000
8.0–9.0
0.0003
0.0100
0.0058
0.0000
0.0002
0.0000
0.0056
0.0000
0.0094
0.0000
9.0–10.0
0.0013
0.0096
0.0055
0.0000
0.0055
0.0000
0.0096
0.0000
0.0064
0.0000Figure 17
The Pareto optimal trajectories of five selected solutions for case III. (a) u1, (b) u2, (c) u3, (d) u4, and (e) u5.
(a)
(b)
(c)
(d)
(e)Roubos et al. setPr=5 and solve (24) using GA with a population size of 40 and generation of 25000 and obtained an optimum of 0.8149 [50]. We calculate the combined objective values with Pr=5 and obtain an optimal value of 0.81345. Figure 18 also shows the combined objective values along the Pareto front by setting Pr=5.Figure 18
The objective value of the Pareto optimal front for case III.


## 6. Conclusions

In this paper, we put forward an extended multiobjective differential evolution algorithm called DEHC by integrating multistrategies from three aspects, that is, hybrid selection strategy, external archive with cyclic crowding estimation, and a two-strategy based on mutation component. We have systematically studied the performance of DEHC by testing it on a set of MOPs. The numerical results show that DEHC is suitable for solving MOPs with different types of Pareto front and that it is superior to or competitive with MOEA/D, NSGA-II, and NNIA. Our proposed DEHC is able to generate very competitive results in terms ofIGD and IS indicators. Finally, three well-known multiobjective optimal control problems in chemical process have been solved to further demonstrate the effectiveness and applicability of DEHC for real-world complex MOOCPs.The numerical results validate that DEHC is an advantageous approach to solving complex MOPs. This advantage comes from three aspects: (1) the hybrid selection strategy can deal with the conflicting objectives from different aspects which is advantageous to different problems during different search processes, (2) the cyclic crowding estimation is helpful to maintain a uniformly distributed nondominated solutions, and (3) the candidate pool of mutation strategies for DE mutation component is beneficial for improving DEHC’s exploration and exploitation.What remain interesting for future work would be constrained MOOCPs and MOOCPs with uncertainties. Another aspect is developing effective CVP methods to convert the original MOOCPs.



---

*Source: 2317860-2018-09-16.xml*