Four-Objective Optimization of Irreversible Atkinson Cycle Based on NSGA-II

Variation trends of dimensionless power density (PD) with a compression ratio and thermal efficiency (TE) are discussed according to the irreversible Atkinson cycle (AC) model established in previous literature. Then, for the fixed cycle temperature ratio, the maximum specific volume ratios, the maximum pressure ratios, and the TEs corresponding to the maximum power output (PO) and the maximum PD are compared. Finally, multi-objective optimization (MOO) of cycle performance with dimensionless PO, TE, dimensionless PD, and dimensionless ecological function (EF) as the optimization objectives and compression ratio as the optimization variable are performed by applying the non-dominated sorting genetic algorithm-II (NSGA-II). The results show that there is an optimal compression ratio which will maximize the dimensionless PD. The relation curve of the dimensionless PD and compression ratio is a parabolic-like one, and the dimensionless PD and TE is a loop-shaped one. The AC engine has smaller size and higher TE under the maximum PD condition than those of under the maximum PO condition. With the increase of TE, the dimensionless PO will decrease, the dimensionless PD will increase, and the dimensionless EF will first increase and then decrease. There is no positive ideal point in Pareto frontier. The optimal solutions by using three decision-making methods are compared. This paper analyzes the performance of the PD of the AC with three losses, and performs MOO of dimensionless PO, TE, dimensionless PD, and dimensionless EF. The new conclusions obtained have theoretical guideline value for the optimal design of actual Atkinson heat engine.

The WF's heat absorption rate in the cycle is The WF's heat release rate in the cycle is where v C ( p C ) is the specific heat of WF at constant volume (pressure) and m  is the molar flow rate of the WF. According to References [28,46,47], the compression and expansion efficiencies in two adiabatic processes 1 2 → and 3 4 → are defined to represent the IIL of the cycle The compression ratio γ and the maximum temperature ratio τ of the AC are defined as: / T T τ = (6) In additional, the entropy change of the working fluid equals zero after a cycle, one has: According to the property of isentropic process, one has: For the actual AC, the HTL between cylinder wall and WF cannot be negligible. According to Reference [48], the heat absorption rate in process 2 3 → is where A and B are the heat released rate by fuel combustion and the HTL coefficient, respectively, and 0 T is ambient temperature.
Equation (11) shows that the cycle heat absorption rate includes two parts. The total heat absorption rate by WF is equal to the difference between the heat release rate by fuel combustion and the HTL rate. Therefore, the HTL rate is: The WF's heat release rate in the cycle is where C v (C p ) is the specific heat of WF at constant volume (pressure) and . m is the molar flow rate of the WF.
According to References [28,46,47], the compression and expansion efficiencies in two adiabatic processes 1 → 2 and 3 → 4 are defined to represent the IIL of the cycle The compression ratio γ and the maximum temperature ratio τ of the AC are defined as: In additional, the entropy change of the working fluid equals zero after a cycle, one has: According to the property of isentropic process, one has: From Equations (3)-(8), one has: For the actual AC, the HTL between cylinder wall and WF cannot be negligible. According to Reference [48], the heat absorption rate in process 2 → 3 is where A and B are the heat released rate by fuel combustion and the HTL coefficient, respectively, and T 0 is ambient temperature. Equation (11) shows that the cycle heat absorption rate includes two parts. The total heat absorption rate by WF is equal to the difference between the heat release rate by fuel combustion and the HTL rate. Therefore, the HTL rate is: where T 0 is the ambient temperature and B 1 = B/2. In the actual AC, the FL between the piston and cylinder wall should also be considered. According to the treatment method of Otto cycle in Reference [49], the friction force is where µ and x are the friction coefficient and the displacement of piston, respectively. The FL power is obtained The piston average speed v is used to replace the piston movement speed v, where x 1 is the position of the piston at the maximum volume, x 2 is the position of the piston at the minimum volume, and ∆t 12 is the time consumed by the power stroke. The cycle PO is obtained: where b = µx 2 2 /(∆t 12 ) 2 . The TE is written as: According to Reference [31], the PD is defined as: According to the v 1 /v 4 = T 1 /T 4 , one has: There are HTL, FL, and IIL in the actual irreversible AC. The entropy generation rate caused by HTL and FL are, respectively: The entropy generation rate due to the IIL is calculated by the entropy increase rates in processes 2s → 2 and 4s → 4 mC v ln(T 2 /T 2s ) mC p ln(T 4 /T 4s ) (23) After the power stroke, the WF is discharged to the environment in the exhaust stroke. The entropy generation rate caused in this process is: The total entropy generation rate of the AC is: According to the Reference [50], the EF is defined as: According to the treatment method of reversible Atkinson cycle by Chen et al. [31], the dimensionless PO, dimensionless PD, and dimensionless EF are defined as, respectively: when γ, T 1 , and τ are given, the temperatures at each state point can be solved, and then the numerical solutions of P, η, P d , and E can be obtained.

Power Density Analysis and Optimization
Based on References [26,27,32], the following parameters are determined: B = 2.2 W/K, T 0 = 300 K, T 1 = 350 K, . m = 1 mol/s, C v = 20.78 J/(mol · K), k = 1.4, τ = 4.28 − 6.28, and b = 20 W. Figures 2 and 3 show the influence of the temperature ratio (τ) on dimensionless PD and compression ratio (P d − γ) and dimensionless PD and TE (P d − η) characteristics, respectively. The shape of the P d − γ curve is parabolic-like one, and there is an optimal value γ P d which makes the P d reach the maximum value (P d ) max . The shape of the P d − η curve is a loop-shape one back to the origin, and there is another optimal value γ η which makes the η reach the maximum value η max . With the increase of τ, the γ P d and the η P d increase. The numerical calculations show that when τ increases from 5.78 to 6.78, the η P d increases from 0.4402 to 0.4684 and increase by 6.41%.          Figure 5 shows the influences of η c , η e , B, and b on P d − η characteristics. The curve 1 reflects the P d − η characteristic when the cycle is completely reversible, and the shape of curve is a parabolic-like one (η P d 0 but (P d ) η = 0). The other curves reflect the P d − η characteristics when one or more irreversibility is considered, the shape of curve is a loop-shaped one back to the origin (both the η P d and (P d ) η are not zero). Comparing curves 1 and 1 , 2 and 2 , 3 and 3 , as well as 4 and 4 in Figure 5, one can see that η P d increases with the decrease of IIL (η c and η e are increased). The numerical calculations show that when B = 2.2 W/K and b = 20 W, and η c and η e increases from 0.94 to 1, η P d increases from 0.4551 to 0.5456, and increases by 19.89%. Comparing curves 1 and 2, 3 and 4, 1 and 2 , as well as 3 and 4 in Figure 5, one can see that η P d decreases with the increase of FL. The numerical calculations show that when η c = η e = 0.94 and B = 2.2 W/K, and b increases from 0 W to 20 W, and η P d decreases from 0.5021 to 0.4551, and decreases by 9.36%. Comparing curves 1 and 3, 2 and 4, 1 and 3 , as well as 2 and 4 in Figure 5, one can be seen that η P d decreases with the increase of HTL. The numerical calculations show that when η c = η e = 0.94 and b = 20 W, and B increases from 0 W/K to 2.    Under the conditions of maximum PO (P max ) and maximum PD ((P d ) max ), Figure 6 shows the relations of maximum specific volume ratio v 4 /v 1 and τ. The (v 4 /v 1 ) P d is smaller than the (v 4 /v 1 ) P when τ is a constant. The numerical calculations show that when τ = 6.28, (v 4 /v 1 ) P is 2.61 and (v 4 /v 1 ) P d is 2.32. Compared with (v 4 /v 1 ) P , (v 4 /v 1 ) P d decreases by 11.1%. The size of the AC engine is smaller under the condition of (P d ) max .
Under conditions of P max and (P d ) max , Figure 7 shows the relations of maximum pressure ratio p 3 /p 1 and τ. It can be seen that (p 3 /p 1 ) P is always smaller than (p 3 /p 1 ) P d when τ is a constant. It means that the decrease of the AC engine size is accompanied by the increase of maximum pressure ratio in the cycle. Figure 8 shows the relations of the η and τ. One can see that when there are three losses, η P d is larger than η P . When τ = 6.28, η P is 0.4389 and η P d is 0.4549. Compared with η P , η P d increases by 3.65%.
relations of maximum specific volume ratio 4 1 / v v and τ . The  It means that the decrease of the AC engine size is accompanied by the increase of maximum pressure ratio in the cycle.        8 show that when there are three losses, compared with the maximum PO condition, η P d at the maximum PD condition increases by 3.65%, while (v 4 /v 1 ) P d at the maximum PD condition decreases by 11.1%. The results show that the TE is larger and the size of the heat engine is smaller when the (P d ) max is taken as the objective.

Four Objective Optimization and Decision-Making Based on NSGA-II Algorithm
The NSGA-II algorithm [51] is a MOO algorithm based on genetic algorithm, which is based on Pareto optimal solution. When the γ is used as the optimization variable, the dimensionless PO, TE, dimensionless PD, and dimensionless EF cannot be optimized at the same time. Pareto put forward the concept of Pareto domination in 1986. It is impossible to optimize the solution for any objective without making other objectives worse. Since there is no optimal solution to make multiple objectives reach the optimal at the same time, the MOO algorithm gives a series of non-inferior solutions. Compared with other solutions, these non-inferior solutions have the least conflict of objectives, which can provide a better choice space for decision makers. These solution sets are called the optimal Pareto solution sets, and the corresponding objective functions are called the Pareto frontier. The specific algorithm flow chart is shown in Figure 9. There are multiple feasible optimal solutions in Pareto frontier. The decision-making methods such as linear programming technique for multidimensional analysis of preference (LINMAP), technique for order preferences by similarity to ideal solution (TOPSIS) and Shannon entropy are used to select the suitable solution from Pareto frontier. According to Reference [51], the deviation index D is introduced to select the most suitable method. without making other objectives worse. Since there is no optimal solution to make multiple objectives reach the optimal at the same time, the MOO algorithm gives a series of non-inferior solutions. Compared with other solutions, these non-inferior solutions have the least conflict of objectives, which can provide a better choice space for decision makers. These solution sets are called the optimal Pareto solution sets, and the corresponding objective functions are called the Pareto frontier. The specific algorithm flow chart is shown in Figure 9. There are multiple feasible optimal solutions in Pareto frontier. The decision-making methods such as linear programming technique for multidimensional analysis of preference (LINMAP), technique for order preferences by similarity to ideal solution (TOPSIS) and Shannon entropy are used to select the suitable solution from Pareto frontier. According to Reference [51], the deviation index D is introduced to select the most suitable method. In order to obtain the optimization design variable of the cycle, the program is conducted by using the "gamultiobj" function of MATLAB. Setting the population "populationsize" as 500 and the algebra "generations" as 1000, the Pareto frontier corresponding to the MOO and the optimal solutions by using three decision-making methods are obtained, as shown in Figure 10. The color on the Pareto frontier edge indicates the size of d P . The positive triangle represents the positive ideal point, the inverted triangle represents the negative ideal point, and the square represents the point corresponding to the LINMAP and TOPSIS decision-making method (the optimal points are the same); the diamond represents the point corresponding to the Shannon entropy decision-making method. According to Figure 10, one can see that, with the increase of TE, the dimensionless PO decreases, the dimensionless PD increases, and the dimensionless EF first increases and then In order to obtain the optimization design variable of the cycle, the program is conducted by using the "gamultiobj" function of MATLAB. Setting the population "populationsize" as 500 and the algebra "generations" as 1000, the Pareto frontier corresponding to the MOO and the optimal solutions by using three decision-making methods are obtained, as shown in Figure 10. The color on the Pareto frontier edge indicates the size of P d . The positive triangle represents the positive ideal point, the inverted triangle represents the negative ideal point, and the square represents the point corresponding to the LINMAP and TOPSIS decision-making method (the optimal points are the same); the diamond represents the point corresponding to the Shannon entropy decision-making method. According to Figure 10, one can see that, with the increase of TE, the dimensionless PO decreases, the dimensionless PD increases, and the dimensionless EF first increases and then decreases. There is no point on the Pareto frontier which will make the dimensionless PO, TE, dimensionless PD, and dimensionless EF reach the maximum values at the same time, i.e., the Pareto frontier does not include positive ideal point.  Table 1 shows the comparisons of the optimal solutions gained by using the MOO to optimize the performance of the irreversible AC model with P , η , d P , and E as the optimization objectives. It can be seen, from Table 1, that the results gained by using LINMAP and TOPSIS decision-making methods are the same. Compared with results gained by using Shannon entropy decision-making method, the optimal compression ratios gained by using LINMAP and TOPSIS decision-making methods are smaller. The D by using Shannon entropy decision-making method is the largest. In the actual decision-making process, the optimal decision-making method should be selected according to different design requirements.

Conclusions
Based on the irreversible AC model with constant specific heat established in Reference [28], the effects of cycle temperature ratio, HTL, FL, and IIL on PD were analyzed, and the optimization results under maximum PO and maximum PD were compared. By using the NSGA-II algorithm and taking γ as the optimization variable, the corresponding Pareto frontiers with the dimensionless PO, TE, dimensionless PD, and dimensionless EF as the optimization objectives were obtained. The results show that:  Table 1 shows the comparisons of the optimal solutions gained by using the MOO to optimize the performance of the irreversible AC model with P, η, P d , and E as the optimization objectives. It can be seen, from Table 1, that the results gained by using LINMAP and TOPSIS decision-making methods are the same. Compared with results gained by using Shannon entropy decision-making method, the optimal compression ratios gained by using LINMAP and TOPSIS decision-making methods are smaller. The D by using Shannon entropy decision-making method is the largest. In the actual decision-making process, the optimal decision-making method should be selected according to different design requirements.

Conclusions
Based on the irreversible AC model with constant specific heat established in Reference [28], the effects of cycle temperature ratio, HTL, FL, and IIL on PD were analyzed, and the optimization results under maximum PO and maximum PD were compared. By using the NSGA-II algorithm and taking γ as the optimization variable, the corresponding Pareto frontiers with the dimensionless PO, TE, dimensionless PD, and dimensionless EF as the optimization objectives were obtained. The results show that: (1) The relationship curve of cycle P d − γ is parabolic-like one. There is an optimal γ which can maximize the PD. With the decrease of τ and the increases of FL and IIL, the PD of cycle decreases. (2) The relationship curve of cycle P d − η is loop-shaped one. With the decrease of γ and the increases of three losses, the corresponding TE at the maximum PD decreases. (3) The efficiency η P d under the condition of (P d ) max is larger than the efficiency η P under the condition of P max , and the corresponding (v 4 /v 1 ) P d is smaller than (v 4 /v 1 ) P . The AC engine designed under the condition of (P d ) max has smaller size and higher TE. (4) For the results by using MOO, with the increase of TE, the dimensionless PO decreases, the dimensionless PD increases, and the dimensionless EF first increases and then decreases.
There is no point on the Pareto frontier which will maximize P, η, P d and E, i.e., the positive ideal point is not on Pareto frontier.

Acknowledgments:
The authors wish to thank the editor reviewers for their careful, unbiased and constructive suggestions, which led to this revised manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
A Heat released rate by fuel (W) B Heat transfer loss coefficient (W/K) C p Specific heat at constant pressure (J/(mol · K)) C v Specific heat at constant volume (J/(mol · K)) E Ecological function (W) k Specific heat ratio .