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

School of Mechanical Engineering, Shanghai University of Engineering Science, Shanghai 201620, China School of Electrical and Information Engineering, Jiangsu University, Zhenjiang, Jiangsu 212013, China School of Energy and Power Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China College of Engineering, Shanghai Second Polytechnic University, Shanghai 201209, China


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][4][5].One intrinsic characteristic common in the single-objectivebased 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 Figure 1, 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.
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][13][14].However, it must be emphasized that DE was originally proposed for singleobjective 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][17][18][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][22][23][24][25][26][27][28][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][31][32][33][34][35][36][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.Section 2 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.

Basic Concepts
2.1.MOOCP: Multiobjective Optimal Control Problem.Consider the system x = f x t , u t , t with state variables x t ∈ ℝ D , initial condition x t 0 = x 0 , and control variables u t ∈ ℝ K .Then the main goal of solving the problem is to find the optimal control variables u * t , t ∈ t 0 , t f that drive the plant along the trajectory x * t , t ∈ t 0 , t f such that all 2 Complexity the cost functions are minimized where the final time t f is fixed.The structure of a typical MOOCP for a continuous process can be generally described as follows: Each of the objective functions 1) can be formulated as follows: where M is the objective number, J i ⋅ is the total cost function, ϕ i ⋅ is the final time performance index, ψ i ⋅ is the integrated performance during the operation, t 0 is the initial time, and t f is the final time; L ∈ ℝ K and U ∈ ℝ K are the lower and upper boundary of control variables, respectively.Definition 3. Pareto optimal front is defined as PF = J u * | 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 (2) DE/best/1 (3) DE/current-to-best/1 (4) DE/best/2 (5) DE/rand/2 where indices r i , i = 1, … , 5 are mutually exclusive integers randomly generated within 1, N , respectively, that are also different from index i; x g best 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: where sn is an arbitrary number in 1, 2, … , n and u i,j is the jth element of the ith new trial vector.After generating the trial vectors, a greedy selection operator is performed as follows: where is the ith solution in the main population of the next generation.

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.(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 Figure 2. Below, we elaborate on some of the main steps in detail.

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 as E, 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 4 Complexity 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.

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 a O mN 2 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 Figure 3.
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 Section 3.1.

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.
After adding Z 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 ① The 1st selection strategy; ② The 2nd selection strategy; ③ The 3rd selection strategy.
Figure 3: The main selection process using HSS.
5 Complexity solution one by one or cyclically instead of retaining the least crowed solutions at a time.This process can be described by Algorithm 2.

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 on M-independent sortings at N solutions, which has a O MN log N computational complexity.
For the hybrid selection strategy, each selection strategy has a subpopulation with size N s = 1/3 N. The decompositionbased strategy has a time complexity of O MN s T , where T is the neighborhood size.For the nondominated sorting strategy, the computational complexity is O M N s 2 .As for the nondominated neighbor-based strategy, the worst computational complexity is O N s lg N s .So, the worst total computational complexity of HSS is As for the CCE, in the worst case, the time complexity for calculating the distances is N Q 2 ; the time complexity for adding the K + 1 values is O KQ ; the time complexity for removing K solutions from E is O KQ 2 .
Based on the above analysis, the procedure of HSS does not increase any burden on the runtime complexity than The main truncating process by the CCE process.6 Complexity using the single-selection strategy.As for the multiplemutation 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 only O MQ log Q .

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: where d v, A is the minimum Euclidean distance between v and the points in R. If R is large enough to represent the PF very well, IGD A, R could measure both the diversity and convergence of in a sense.To have a low value of IGD A, 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 where M is the objective number and d e i is the Euclidean distance between the ith extreme solutions in A and R and d i,j is the Euclidean distance between the ith solution and its jth nearest solution in A. This one is a harmonic distancebased 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.

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.Table 1 shows the mean and standard deviation (Std) of IG D by 5 algorithms.It can be observed in Table 1 that DEHC is able to achieve the lowest values on most problems, 7 Complexity which means that DEHC converges better than any other algorithms.
Table 2 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.
In order to determine whether DEHC is significantly better than other algorithms, statistical significance tests of IGD 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 β.
To provide more information about the convergence performance of the 5 MOEAs, we calculate IGD 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.
For illustration, we also show some typical fronts obtained by the 5 algorithms.Figure 8 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 9 presents the fronts on DTLZ6 and these fronts correspond to the run with a median value with respect to I GD.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.
From the above comparisons and analyses, we can draw the conclusion that in terms of IGD 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.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 12 Complexity 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 Table 4 that for the five 2objective 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 2objective 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.

Some Studies on Scalability of DEHC. To study how the performance of DEHC varies as the number of decision
To confirm that CCE is really suited for MOPs with many objectives, here, we expand DTLZ2 to 5, 7, and 9 objectives and then study IS. 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.

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 range t 0 , t f is divided into T equal 13 Complexity segments, that is, t 0 , t f = t 0 , t 1 × t 1 , t 2 × ⋯× t T−1 , t T , t T = t f , 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 u t as piecewise constant function.As a result, T × K parameters determine u t over t 0 , t f .In terms of the jth segments; the control variables over the entire time span are as follows: So, the total control variables over the entire time span can be described as the summation below: 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: 14 Complexity Step 1. Divide the time range t 0 , t f into T equal segments and denote the jth segment as t j−1 , t j , where j = 1, 2, … , K and t T = t f .
Step 2. Assume that the value of the ith input variable within the jth segment is u j i , where i = 1, 2, … , m and j = 1, 2, … , T.
Step 3. Once the value of u t is determined within the time range t 0 , t f , a numerical integral method, such as the standard Runge-Kutta method, can be employed to obtain the objective value, that is, J u v .
Step 4. Apply DEHC to obtain the optimal u j i and J u t as an approximation of the continuous problem, where i = 1, 2, … , K and j = 1, 2, … , T.

Case I: Catalyst Mixing Problem in a Tubular Reactor.
This problem considers a steady-state plug flow reactor of fixed length t f .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]: where x 1 and x 2 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. max For this problem, we set T = 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.
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.

Case II:
Optimal Operation of a Fed Batch Reactor.The third case is a batch reactor based on the two reaction systems [47]: This process can be described by the following mechanistic model: 15 Complexity reactant feed rate, and k 1 = 0 5 and k 2 = 0 5 are the reaction rate constants.At the start of reaction, the reactor contains A 0 = 0 2 mole/L of A; no B B 0 = 0 and is fed to 50% (V 0 = 0 5 m 3 ).The objective is to maximize the amount of product C and minimize the amount of by-product D.
For this problem, we set T = 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.
The obtained Pareto optimal front is illustrated in Figure 14.When the complete focus is put on maximizing the amount of product C, Chen et al. obtained a best value   16 Complexity 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.[48].In this case, there are 7 state variables which are reaction volume (x 1 , L), cell density (x 2 , g/L), nutrient concentration (x 3 , g/L), foreign protein concentration (x 4 , g/L), inducer concentration (x 5 , g/L), inducer shock factor on the cell growth rate (x 6 ), and the inducer recovery factor on the cell growth rate (x 7 ).The two control variables are glucose feed rate (u 1 , L/h) and inducer feed rate (u 2 , L/h).The model was described as follows [49]:   As the price of the inducer, that is, P r , 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: To assist decision making, five Pareto optimal reaction temperatures along the front in Figure 16 with equal distribution are shown in Table 7 and the corresponding Pareto optimal trajectories of these are illustrated in Figure 17.
Roubos et al. set P r = 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 P r = 5 and obtain an optimal value of 0.81345.Figure 18 also shows the combined objective values along the Pareto front by setting P r = 5.

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 of IGD 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.

Figure 1 :
Figure 1: The schematic graph to illustrate the principle of solving an OCP by EAs.

3 Complexity( 2 )
In the selection step, the greedy strategy-based selection component is replaced by a new selection component based on hybrid selection strategy.

Figure 10 :
Figure 10: Mean IGD on ZDT1 as the number of decision variable increases.

Figure 11 :Figure 12 :
Figure 11: Mean IS on DTLZ2 as the number of objectives increases.

23 For
this problem, we set T = 10, the population size is N = 50, and the maximum generation is 100.The obtained Pareto fronts by DEHC and MODE are showed in Figure16.It can be observed form Figure16that 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 18 :
Figure 18: The objective value of the Pareto optimal front for case III.

Table 1 :
The IGD comparison results of DEHC with 4 MOEAs on 12 problems.

Table 2 :
The IS comparison results of DEHC with 4 MOEAs on 12 problems.

Table 3 :
The distribution of IGD and IS using Mann-Whitney rank sum test.

Table 4 :
Comparison results of CCE and other truncating methods on IS.

Table 5 :
Five solutions selected from the front for case I.

Table 6 :
Five solutions selected from the front for case II.