An Improved NSGA-II Algorithm Based on Adaptive Weighting and Searching Strategy

Abstract

Non-dominated sorting genetic algorithm II is a classical multi-objective optimization algorithm but it suffers from poor diversity and the tendency to fall into a local optimum. In this paper, we propose an improved non-dominated sorting genetic algorithm, which aims to address the issues of poor global optimization ability and poor convergence ability. The improved NSGA-II algorithm not only uses Levy distribution for global search, which enables the algorithm to search a wider range, but also improves the local search capability by using the relatively concentrated search property of random walk. Moreover, an adaptive balance parameter is designed to adjust the respective contributions of the exploration and exploitation abilities, which lead to a faster search of the algorithm. It helps to expand the search area, which increases the diversity of the population and avoids getting trapped in a local optimum. The superiority of the improved NSGA-II algorithm is demonstrated through benchmark test functions and a practical application. It is shown that the improved strategy provides an effective improvement in the convergence and diversity of the traditional algorithm.

1. Introduction

The principle of multi-objective optimization [1] is that when multiple conflicting objectives need to be achieved in a particular situation, the optimization of one objective is often accompanied by degradation of the performance of other objectives, so it is necessary to find a set of solutions [2] that balance multiple sub-objectives to obtain the best optimization for each sub-objective.

Traditional mathematical methods cannot achieve more satisfactory results because multi-objective optimization problems are computationally difficult. Evolutionary algorithms [3] can be a good solution [4]. The most widely used evolutionary multi-objective optimization algorithm for discovering Pareto optimal solutions is known as the non-dominated sorting genetic algorithm II. In 1989, Goldberg proposed the concept of Pareto-based optimal solutions and the calculation of the individual fitness method, where the population evolves toward the Pareto optimal solution through non-inferior solutions and the corresponding selection operations.

The NSGA algorithm, which was first proposed by Professors Srinivas and Deb in the early 1990s, is based on classifying all individuals in a population at different levels. NSGA uses a non-dominated sorting method that allows good individuals to have a greater chance of being inherited by the next generation, and an adaptation sharing strategy that allows individuals to be evenly distributed and maintains population diversity. Indian scientist Deb proposed a non-dominance ranking genetic algorithm based on the NSGA algorithm in 2002 [5] and used an elite strategy by applying the Pareto dominance relationship and crowding distance mechanism to update the population.

NSGA-II has better convergence and distribution [6] because of the use of fast non-dominated sorting and crowded-distance-sorting mechanisms. However, NSGA-II has a slow convergence speed due to the use of an inefficient simulated binary crossover algorithm. It has poor convergence speed and optimization finding ability when solving complex problems [7]. With the NSGA-II algorithm applied to a wider range, some problems have gradually emerged. In the genetic operation, constant genetic parameters lead to a smaller moving space and poorer search performance during the iteration of the algorithm. Meanwhile, the NSGA-II algorithm shows the problem of falling into local optimal solutions.

The principal contributions of this study are summarized as follows:

1.

An improved non-dominated ranking genetic algorithm is proposed to enhance the search ability by using Levy flight and random walk strategy, which can provide a larger search range and better local search ability.

2.

To improve the accuracy and speed of convergence, an adaptive parameter is designed, which balances the exploration and development periods.

In the remainder of the paper, we first introduce some fundamental conceptions of multi-objective optimization problems and discuss some existing works in Section 2. At the beginning of Section 3, we present some background on NSGA-II, then introduce Levy distribution and random walk-in detail, while designing an adaptive weight to balance global exploration and local exploitation. After that, we propose the improved NSGA-II algorithm. In Section 4, we submit the results on standard test functions and a mathematical model for the parallel chillers. At the end of the article, we outline the conclusions of the study.

2. Preliminaries

2.1. Multi-Objective Optimization

In engineering, we often encounter design and decision problems under multiple criteria or multiple design objectives, which are often contradictory for finding the best design solution to satisfy these objectives, which is the multi-objective optimization problem [8]. As an example of minimizing MOPs, define the relevant concepts of MOPs as follows:

 Definition 1. 

Multi-objective optimization problem. For a MOP with n-dimensional decision variables and m-dimensional objectives, the mathematical model is defined as

$$\left\{\begin{array}{c}minF\left(x\right)={\left(f_1\left(x\right),f_2\left(x\right),\dots ,f_{\mathrm{m}}\left(x\right),\right)}^T\hfill \\ M_i\left(x\right)\ge 0,\mathrm{i}\in \{1,2,\dots ,\mathrm{p}\}\hfill \\ N_j\left(x\right)=0,\mathrm{j}\in \{1,2,\dots ,\mathrm{q}\}\hfill \end{array}\right.$$

(1)

where $\mathbf{x}=\left({\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_n\right)\in \mathsf{\Omega }$ is the decision variable, and Ω is the decision space. F(x) is an objective function consisting of multiple mutually contradictory sub-functions, and M(x) and N(x) are constraint functions with independent variable x.

 Definition 2. 

Pareto optimal solutions [9]. In the context of multi-objective optimization, objective functions conflict with each other or cannot be compared in a simple way. The solution of a problem that is optimal in the first objective function may appear to be the worst solution in the second objective function, so people refer to a series of solutions that cannot be simply compared in the computation of multiple objective optimization problems as non-dominated solutions or Pareto optimal solutions.

 Definition 3. 

Pareto dominate. Supposing that $x_1,x_2ϵX$ are two decision variables in the objective space that satisfy the constraints, the definition of the dominance between the two solutions can be expressed as $x_1$ dominating $x_2$ if

$$\left\{\forall \mathrm{i}\in \{1,2,\dots ,\mathrm{m}\},f_1\left(x_1\right)\le f_1\left(x_2\right)\right\}\cap \left\{\exists \mathrm{j}\in \{1,2,\dots ,\mathrm{m}\},f_j\left(x_1\right)<f_j\left(x_2\right)\right\}.$$

(2)

Therefore, a design x is said to be non-dominated [10] if no other feasible design dominates it.

2.2. NSGA-II Algorithm

Deb et al. proposed the fast non-dominated sorting genetic algorithm with elite strategy (NSGA-II) [5] in 2000, which specifically shows its advantages from the following three points: (1) It provides a hierarchy-based non-dominated sorting method to select the optimal solution by hierarchical sorting. Compared with the original method, the operation speed is significantly improved, unnecessary calculations are reduced, and the speed of problem solving is greatly improved. (2) The parent population is fused with the offspring population to retain the better individuals by selecting the optimal solution through an elite strategy. (3) A crowding distance comparison operator is given to solve the original problem that requires a specified shared radius, and makes the optimal solution evenly distributed [9] in space.

In order to determine the degree of superiority of each solution at the identical ranking level in the results of algorithm, a crowding distance is assigned to all solutions. The basic idea is to spread out the already solved Pareto optimal solutions from the search domain to the greatest extent possible, and this process is also a performance indicator to ensure diversity among the solutions. The crowded distance comparison operator can be used as a solution for each stage of the selection algorithm by ranking individuals in terms of crowded distance comparison and taking the top N individuals to drive the Pareto optimal front in an evenly distributed state.

The specific implementation process of NSGA-II algorithm is as follows:

Step 1. Population initialization. The population size is N, and the Gen = 1 is set as the number of evolutionary generations;

Step 2. Based on the execution of selection, crossover and mutation operators, determine whether the first generation of offspring population has been formed. If it has been generated, make the evolutionary generation Gen = 2. If not, perform the non-dominant sorting operation on the initial population and continue to perform the selection, crossover, and variation operations to produce a population of first-generation offspring with evolutionary generation Gen = 2;

Step 3. Generate a new population by integrating parent and offspring populations;

Step 4. First determine whether a new parent population has been formed, and if not, the objective function of individuals in the new population is calculated, and a new parent population is generated by performing operations such as fast non-dominance sorting, calculating crowding [11], and elite strategy. Otherwise, go to step 5;

Step 5. Generate the offspring population by performing selection, crossover, and variation operators on the generated parent population;

Step 6. First determine whether Gen is equal to the maximum evolutionary generation; if not, then make the evolutionary generation Gen = Gen + 1 and return to step 3. Otherwise, the algorithm runs to the end.

3. Improved NSGA-II Algorithm

3.1. Levy Flight Strategy for Global Search

Levy flight [12] belongs to the Markov process. In optimization problems, Levy flight covers a larger area with fewer distances and steps, which is useful to explore the unknown. The levy flight strategy grows faster than the Brownian motion due to the stochastic nature of the step size; therefore, it achieves better results than the Brownian random motion when searching in a large unknown range.

3.1.1. Levy Distribution

Levy flight belongs to random wandering, which means that its trajectory cannot be accurately predicted. Both continuous Brownian motion and Poisson processes are Levy processes, and the most essential characteristic of Levy processes is that they have smooth independent increments.

The mathematical form of the Levy distribution is as follows:

$$f_X\left(x;\sigma ,\mu \right)={\left(\frac{\sigma }{2\pi }\right)}^{1/2}\frac{1}{{(x-\mu )}^{3/2}}exp\left\{-\frac{\sigma }{2(x-\mu )}\right\}\mu <x<\infty$$

(3)

in which $\mu$ is the position parameter, which affects the left and right translation of the distribution curve so that the distribution interval is $\left(\mu ,\infty \right)$, and c is the scale parameter.

3.1.2. Levy Flight

Levy flight is a random walk between short-term search and occasional long-term walk, which drives Levy flight to have an excellent capacity for overall search.

The position update function for Levy flight is as follows:

$$x_i^{{}^{\prime }}\left(t\right)=x_i\left(t\right)+l\times Levy\left(\lambda \right)$$

(4)

in which x denotes the current position and l denotes the parameter to adjust the step size. The Levy distribution is often simulated by the man algorithm in many studies due to its complexity:

$$s=\frac{\mu }{{\left|\nu \right|}^{1/\gamma }}$$

(5)

in which $\mu$ and $\nu$ follow a normal distribution and

$${\sigma }_{\mu }={\left\{\frac{\mathrm{\Gamma }(1+\gamma )sin(\pi \gamma /2)}{\gamma ·\mathrm{\Gamma }[(\gamma +1)/2]·2^{(\gamma -1)/2}}\right\}}^{1/\gamma }$$

(6)

$${\sigma }_{\nu }=1$$

(7)

where $\gamma$ usually takes the value of 1.5.

The long and short hop features of the Levy flight strategy are used to randomly update each position of the population to enhance the search ability of the global optimum, which can increase the diversity of the population distribution and find the global optimum solution more quickly.

3.2. Random Walk Strategy for Local Search

Random walk uses mixed variation and crossover to generate new solutions, which can enhance the diversity of populations to some extent and prevent the attraction of regional extremes, thus improving the capacity of local search of the algorithm. It can be used to accelerate the speed of the search for optimal solutions.

The position update formula for the random walk is as follows:

$$x_i^{\prime }\left(t\right)=x_i\left(t\right)+\epsilon \left(x_j\left(t\right)-x_k\left(t\right)\right)$$

(8)

where $x_j\left(t\right)$ and $x_k\left(t\right)$ are the two random solutions in the $t_{th}$ generation, and $\epsilon$ is the scaling factor, $\epsilon ϵU\left(0,1\right)$.

3.3. Framework and Details of INSGA-II Algorithm

To overcome the disadvantages of slow convergence and falling into a local optimum [13] in the NSGA-II algorithm, we use an adaptive weight [14] to balance the global search and local search ability to accelerate the convergence speed, while using Levy flight and random walk strategies to improve the global and local search ability of the algorithm, respectively.

The search mechanism of an optimization algorithm usually consists of two steps, which are the global exploration and the local exploitation. Exploration denotes the tendency of an algorithm to behave in a highly randomized manner. Significant changes in the solutions prompt further exploration of the variables space to discover their promising positions. Exploitation is the discovery of promising areas followed by the reduction in random behavior so that the algorithm can search around the promising areas.

Since Levy flight has a characteristic of walking between long and short, it can effectively improve the global search ability of the algorithm and make the population individuals evenly distributed in the space. The random walk can focus the search on the locations in the space, where the optimal solution is more likely to occur and make the algorithm converge to the optimal solution.

In the early stage of optimization, individuals in the population should be distributed extensively throughout the entire search space, with individuals converging to the global optimum during the later stage of optimization. Therefore, we design an adaptive weighting factor C in this study to balance the pre-global exploration and post-local exploitation capability of the algorithm.

The value of C decreases from 1 to 0.02 during iterations. With the increase in iterations, the algorithm gradually switches from global search into local exact search that aims to enhance the overall convergence accuracy, which is calculated as

$$C=1-It\left(\frac{0.98}{\mathrm{MaxIt}}\right)$$

(9)

where It is the current number of iterations and MaxIt is the maximum number of iterations. The maximum number of iterations for the two-objective test function was set to 800, while the maximum number of iterations for the three-objective test function was set to 1500.

The offspring population individuals in the iteration are generated with the formula as follows:

$$X(t+1)=C\times \left(X_t+\alpha \otimes Levy\left(\beta \right)\right)$$

(10)

$$X(t+1)=(1-C)\times \left(X_t+\alpha \otimes Guass\left(\beta \right)\right)$$

(11)

The adaptive weight factor C in this paper is considered a balanced treatment of the speed and accuracy of convergence, which also ensures that the algorithm has some global search capability at the late merit search iterations. The intuitive algorithmic flow is shown in Figure 1.

The population distance D between population G and population F is calculated as follows:

$$D=\sum _{P{G}_g\in {\mathrm{\Omega }}_G^{pos}}sum\left(\sum _{P{F}_f\in Q_F^{pos}}min\left(ED\left(P{F}_g,P{F}_f\right)\right)\right)/NC{H}_G$$

where $\mathrm{\Omega }$ is the set of solutions for a population G or F with a non-dominated ranking of 1. G and F are solutions with non-dominated ranking of 1, respectively. $NC{H}_G$ is the number of populations with a non-dominance ranking of 1. ED is the Euclidean distance.

The algorithm results begin to stabilize if the distance between the populations of neighboring generations is less than the acceptable population distance. The INSGA-II algorithm is converged if it satisfies one of the following conditions:

1.

The population distances of neighboring generations are consecutively less than the threshold.

2.

The number of iterations reaches the maximum set value.

The INSGA-II algorithm is implemented as follows:

1.

Initial population and set the evolutionary generation Gen = 1;

2.

Generate the first generation sub-population by the Levy flight strategy and random wandering strategy and make the evolutionary generation Gen = 2;

3.

Combine the parent population with the offspring population to form a new population.

4.

Determine whether a new parent population has been generated; if not, calculate the objective function of individuals in the new population, and perform operations, such as fast non-dominated sorting, calculating crowding, and elite strategy, to generate a new parent population. Otherwise, go to step 5;

5.

Perform adaptive iteration operations on the generated parent population to generate the child population;

6.

Determine whether the termination condition is reached or whether Gen is equal to the maximum evolutionary generation; if neither is satisfied, the evolutionary generation Gen = Gen + 1 and return to step 3. Otherwise, run to the end.

4. Case Studies

For the purpose of evaluating the performance of the proposed INSGA-II algorithm, different two-objective ZDT and three-objective DTLZ test functions are used to verify the optimization capability and to make a comparison with other common multi-objective optimization algorithms.

4.1. Pareto Performance Evaluation Metrics

In contrast to single-objective optimization algorithms, it is not possible to evaluate the performance of multi-objective optimization results in terms of only one metric for each optimization algorithm. In this subsection, we choose two common evaluation metrics [15], which are the inverted generational distance as well as hypervolume, as a way to quantify the performance of the algorithms and to compare the performance of other different methods:

1.

Inverted generational distance

IGD [16] is an average minimum distance between the solution set obtained from the search for the best and the true solution set. The decrease in the value of IGD means that the algorithm converges more effectively and is more closely related to the real Pareto front. The mathematical expression is as follows:

$$\mathrm{I}\mathrm{G}\mathrm{D}=\frac{{\sum }_{x\in P^{*}}mindis(x,P)}{\left|P^{*}\right|}$$

(12)

in which $P^{*}$ is the ideal Pareto front, P is the collection of optimized solutions acquired from the algorithm, $mindis\left(x,P\right)$ is the minimum Euclidean distance between the solution x of true Pareto and the found Pareto solution P obtained by the algorithm, and $|P^{*}|$ is the number of solutions in the solution set $P^{*}$.

2.

Hypervolume

The hypervolume [17] indicator refers to the volume of the area surrounded by the optimal set of solutions obtained by the algorithm and the selected reference points in the target space. With a larger value of HV, it means that the comprehensive performance of the algorithm is better. The mathematical expression is as follows:

$$\mathrm{HV}(X,P)=\bigcup _{x\in X}^{X}v(x,P)$$

(13)

The HV metric allows assessing the convergence and diversity of the dominant solution set. In addition, the HV metric is obtained without real solution sets, so it can be widely used for engineering problems with unknown results in the real world.

4.2. ZDT Test Functions

In order to prove the feasibility and validity of the improved NSGA-II algorithm proposed, the ZDT series of test functions [18] is selected to study the efficiency of the algorithm in solving multiple-objective optimization problems. These test functions are different, with different characteristics such as convex or concave, and their Pareto optimal frontier shapes are different, the specific descriptions and characteristics of the functions are shown in

Table 1. ZDT test functions and characteristics.

Function Name	Mathematical Formula	Function Characteristic
ZDT1	$f_1\left(x\right)=x_1$ $f_2\left(x\right)=g\left(x\right)\left[1-\sqrt{f_1\left(x\right)/g\left(x\right)}\right]$ where $g\left(x\right)=1+9·\left({\sum }_{i=2}^n{x}_i\right)/(n-1)$ $0\le x_i\le 1,i=1,2,\mathrm{L}n$	Convex, continuous
ZDT2	$f_1\left(x\right)=x_1$ $f_2\left(x\right)=g\left(x\right)\left[1-{\left(f_1\left(x\right)/g\left(x\right)\right)}^2\right]$ where $g\left(x\right)=1+9·\left({\sum }_{i=2}^n{x}_i\right)/(n-1)$ $0\le x_i\le 1,i=1,2,\mathrm{L}n$	Concave, continuous
ZDT3	$f_1\left(x\right)=x_1$ $f_2\left(x\right)=g\left(x\right)\left[1-\left(f_1\left(x\right)/g\left(x\right)\right)-\left(f_1\left(x\right)/g\left(x\right)\right)sin\left(10\pi x_1\right)\right]$              where $g\left(x\right)=1+9·\left({\sum }_{i=2}^n{x}_i\right)/(n-1)$ $0\le x_i\le 1,i=1,2,\mathrm{L}n$	Concave, discontinuous
ZDT4	$f_1\left(x\right)=x_1$ $f_2\left(x\right)=g\left(x\right)\left[\sqrt{1-f_1\left(x\right)/g\left(x\right)}\right]$ where $g\left(x\right)=1+10·(n-1)+{\sum }_{i=2}^n\left[x_i^2-10cos\left(4\pi x_i\right)\right]$ $0\le x_1\le 1,-5\le x_i\le 5,i=2,\mathrm{L}n$	Convex, multi-modal

. Therefore, researchers often use these functions when conducting multi-objective optimization tests experiments. The simulation results are shown in Figure 2.

In this section, to demonstrate the capabilities of the improved NSGA-II, we perform experiments on the test functions with the original NSGA-II and other methods, such as MOEA/D [19], IBEA [20] and SPEA2 [21], as comparison algorithms.

The average values of the experimental algorithms for IGD and HV are presented in

Table 2. IGD metric results for ZDT function.

Test Function		Comparison	Algorithm
	NSGA-II	MOEA/D	IBEA	SPEA2	INSGA-II
ZDT1	0.0052	0.0271	0.0041	0.0042	0.0015
ZDT2	0.0027	0.0132	0.0025	0.0032	0.0016
ZDT3	0.0143	0.0471	0.0083	0.0091	0.0032
ZDT4	0.0072	0.0147	0.0065	0.0062	0.0023

and

Table 3. HV metric results for DTLZ function.

Test Function		Comparison	Algorithm
	NSGA-II	MOEA/D	IBEA	SPEA2	INSGA-II
ZDT1	0.7124	0.6951	0.7199	0.7191	0.9115
ZDT2	0.6074	0.5865	0.6175	0.6327	0.8592
ZDT3	0.7331	0.5753	0.7486	0.7362	0.9269
ZDT4	0.7083	0.6608	0.7168	0.7240	0.8866

, respectively. (Each result is the average of the two metrics calculated for 10 independent runs of the same algorithm over the same problem.)

It is shown in Table 2 and Table 3 that the performance index of the INSGA-II algorithm is superior to four comparison algorithms on test functions. It is obvious from the above experiments that INSGA-II has better convergence accuracy and distribution performance.

4.3. DTLZ Test Functions

In order to verify the feasibility of the algorithm proposed in this chapter in the high-dimensional multi-objective problem, the DTLZ standard test set function is selected for simulation experiments. The simulation results are shown in Figure 3 and the specific expressions are shown in

Table 4. DTLZ Test functions and characteristics.

Function Name	Mathematical Formula
DTLZ1	$\begin{array}{cc}& min{f}_1\left(x\right)=\frac{1}{2}{x}_1{x}_2\cdots x_{M-1}\left(1+g\left(x_M\right)\right)\hfill \\ & min{f}_2\left(x\right)=\frac{1}{2}{x}_1{x}_2\cdots \left(1-x_{M-1}\right)\left(1+g\left(x_M\right)\right)\hfill \\ & ⋮\hfill \\ & min{f}_{M-1}\left(x\right)=\frac{1}{2}{x}_1\left(1-x_2\right)\left(1+g\left(x_M\right)\right)\hfill \\ & min{f}_M\left(x\right)=\frac{1}{2}\left(1-x_1\right)\left(1+g\left(x_M\right)\right)\hfill \\ & g\left(x_M\right)=100\left(\left|x_M\right|+\sum _{x_i\in x_M}{\left(x_i-0.5\right)}^2-cos\left(20\pi \left(x_i-0.5\right)\right)\right)\hfill \end{array}$
DTLZ2	$\begin{array}{cc}& min{f}_1\left(x\right)=\left(1+g\left(x_M\right)\right)cos\left(x_1\pi /2\right)\cdots cos\left(x_{M-2}\pi /2\right)cos\left(x_{M-1}\pi /2\right)\hfill \\ & min{f}_2\left(x\right)=\left(1+g\left(x_M\right)\right)cos\left(x_1\pi /2\right)\cdots cos\left(x_{M-2}\pi /2\right)sin\left(x_{M-1}\pi /2\right)\hfill \\ & ⋮\hfill \\ & min{f}_{M-1}\left(x\right)=\left(1+g\left(x_M\right)\right)cos\left(x_1\pi /2\right)sin\left(x_2\pi /2\right)\hfill \\ & min{f}_M\left(x\right)=\left(1+g\left(x_M\right)\right)sin\left(x_1\pi /2\right)\hfill \\ & g\left(x_M\right)=\sum _{x_1\in x_M}{\left(x_i-0.5\right)}^2\hfill \end{array}$
DTLZ3	$\begin{array}{cc}& min{f}_1\left(x\right)=\left(1+g\left(x_M\right)\right)cos\left(x_1\pi /2\right)\cdots cos\left(x_{M-2}\pi /2\right)cos\left(x_{M-1}\pi /2\right)\hfill \\ & min{f}_2\left(x\right)=\left(1+g\left(x_M\right)\right)cos\left(x_1\pi /2\right)\cdots cos\left(x_{M-2}\pi /2\right)sin\left(x_{M-1}\pi /2\right)\hfill \\ & ⋮\hfill \\ & min{f}_{M-1}\left(x\right)=\left(1+g\left(x_M\right)\right)cos\left(x_1\pi /2\right)sin\left(x_2\pi /2\right)\hfill \\ & min{f}_M\left(x\right)=\left(1+g\left(x_M\right)\right)sin\left(x_1\pi /2\right)\hfill \\ & g\left(x_M\right)=100\left(\left|x_M\right|+\sum _{x_x\in x_M}{\left(x_i-0.5\right)}^2-cos\left(20\pi \left(x_i-0.5\right)\right)\right)\hfill \end{array}$
DTLZ4	$\begin{array}{cc}& min{f}_1\left(x\right)=\left(1+g\left(x_M\right)\right)cos\left(x_1^{\alpha }\pi /2\right)\cdots cos\left(x_{M-2}^{\alpha }\pi /2\right)cos\left(x_{M-1}^{\alpha }\pi /2\right)\hfill \\ & min{f}_2\left(x\right)=\left(1+g\left(x_M\right)\right)cos\left(x_1^{\alpha }\pi /2\right)\cdots cos\left(x_{M-2}^{\alpha }\pi /2\right)sin\left(x_{M-1}^{\alpha }\pi /2\right)\hfill \\ & ⋮\hfill \\ & min{f}_{M-1}\left(x\right)=\left(1+g\left(x_M\right)\right)cos\left(x_1^{\alpha }\pi /2\right)sin\left(x_2^{\alpha }\pi /2\right)\hfill \\ & min{f}_M\left(x\right)=\left(1+g\left(x_M\right)\right)sin\left(x_1^{\alpha }\pi /2\right)\hfill \\ & g\left(x_M\right)=\sum _{x_i\in x_M}{\left(x_i-0.5\right)}^2\hfill \end{array}$
DTLZ5	$\begin{array}{cc}& min{f}_1\left(x\right)=\left(1+g\left(x_M\right)\right)cos\left({\theta }_1\pi /2\right)\cdots cos\left({\theta }_{M-2}\pi /2\right)cos\left({\theta }_{M-1}\pi /2\right)\hfill \\ & min{f}_2\left(x\right)=\left(1+g\left(x_M\right)\right)cos\left({\theta }_1\pi /2\right)\cdots cos\left({\theta }_{M-2}\pi /2\right)sin\left({\theta }_{M-1}\pi /2\right)\hfill \\ & ⋮\hfill \\ & min{f}_{M-1}\left(x\right)=\left(1+g\left(x_M\right)\right)cos\left({\theta }_1\pi /2\right)sin\left({\theta }_2\pi /2\right)\hfill \\ & min{f}_M\left(x\right)=\left(1+g\left(x_M\right)\right)sin\left({\theta }_1\pi /2\right)\hfill \\ & g\left(x_M\right)=\sum _{x_1\in x_M}{x}_i^{0.1}\hfill \\ & {\theta }_i=\frac{\pi }{4\left(1+g\left(x_M\right)\right)}\left(1+2g\left(x_M\right)x_i\right)\hfill \end{array}$
DTLZ6	$\begin{array}{cc}& min{f}_1\left(x_1\right)=x_1\hfill \\ & ⋮\hfill \\ & min{f}_{M-1}\left(x_{M-1}\right)=x_{M-1}\hfill \\ & min{f}_M\left(x\right)=\left(1+g\left(x_M\right)\right)h\left(f_1,f_2,\dots ,f_{M-1},g\right)\hfill \\ & g\left(x_M\right)=1+\frac{9}{\left|x_M\right|}\sum _{x_i\in x_M}{x}_i\hfill \\ & h=M-\sum _{i=1}^{M-1}\left[\frac{f_i}{1+g}\left(1+sin\left(3\pi f_i\right)\right)\right]\hfill \end{array}$

.

To prove the capability of the improved NSGA-II, we conduct experiments on the test function with the original NSGA-II and other methods, such as MOEA/D [19], IBEA [20] and SPEA2 [21], as comparison algorithms.

The average values of the experimental algorithms for IGD are presented in

Table 5. IGD metric results for DTLZ function.

Test Function		Comparison	Algorithm
	NSGA-II	MOEA/D	IBEA	SPEA2	INSGA-II
DTLZ1	0.131	0.236	0.097	0.089	0.085
DTLZ2	0.377	0.576	0.284	0.262	0.205
DTLZ3	1.741	2.464	1.074	1.191	0.807
DTLZ4	0.561	0.896	0.320	0.334	0.231
DTLZ5	0.142	0.264	0.092	0.089	0.053
DTLZ6	0.278	0.529	0.108	0.095	0.052

, respectively. (Each result is the average of the two metrics calculated for 10 independent runs of the same algorithm over the same problem.)

It can be seen from Table 5, the performance index of the INSGA-II algorithm is better than those of the four compared algorithms in terms of test functions. From the above experiments, it can be seen that INSGA-II also has better convergence accuracy and distribution performance in high-dimensional multi-objective optimization problems.

4.4. Parallel Chiller System

The proposed algorithm is applied to the chiller load distribution problem for the sake of verifying the effectiveness of the presented algorithm in practical project problems [22]. The cooling capacity required for the central air conditioning system is provided by the chilled water system, which includes chilled water pumps and a chiller system, as shown in Figure 4. Generally, the chiller system is operated jointly by multiple chillers with different refrigeration performance to manufacture chilled water, and the different refrigeration performance determines that the energy consumption of each chiller is also different. A multi-chiller system is composed of two or more chillers connected by parallel or tandem piping to a common distribution system. The chiller system can adjust its own load size so that each chiller operates at the optimal operating point.

The chiller can adjust its own load size to meet the basic end load demand so that the parallel chiller system can operate at optimal performance. In a refrigeration system, the system performs best when the total power of the chiller system is minimal and the end load demand is met. The problem of the load distribution of parallel chillers is to adjust the load distribution of the parallel chillers to achieve energy savings. At a certain wet bulb temperature, the electrical power of a concentric chiller can be presented as an optimized function of the partial load rate [23] as follows:

$$P_i=a_i+b_i\times PL{R}_i+c_i\times PL{R}_i^2+d_i\times PL{R}_i^3$$

(14)

where $P_i$ is the power of the $i_{th}$ chiller, $PL{R}_i$ is the partial load rate of the $i_{th}$ chiller, and $a_i$, $b_i$, and $c_i$ are the fitting coefficients of the energy consumption curve of the chiller.

An optimal distribution problem for parallel chillers is to achieve the lowest energy consumption and the highest cooling capacity. The abstraction of the actual project into a mathematical problem could be represented as the formula below:

$$minP=\sum _{i=1}^n{P}_i$$

(15)

$$maxQ=\sum _{i=1}^{n}PL{R}_i\times Q_i$$

(16)

where P is the total power of multiple parallel chillers, and Q is the cooling capacity of the chiller.

A chiller system with three chillers of 800 RT cooling capacity was selected for testing in Case 1 so as to confirm that the improved NSGA-II algorithm could be implemented to tackle the load distribution issues. To further verify the algorithm in solving the load distribution problem of parallel chillers, a cooling system which consists of six chillers was selected in Case 2, including four same chillers with a cooling capacity 1280 RT and 2 other chillers of 1250 RT. The purpose was to test whether the improved NSGA-II algorithm can search for optimal values in different multi-chiller systems. Through these two classical test cases, we verified the feasibility of the INSGA-II algorithm for solving practical problems. In the two test cases, the performance parameters of the chillers are shown in

Table 6. Chiller parameters for test cases.

Test Case	Chiller	${\mathit{a}}_{\mathit{i}}$	${\mathit{b}}_{\mathit{i}}$	${\mathit{c}}_{\mathit{i}}$	${\mathit{d}}_{\mathit{i}}$	Capacity
Case 1	1	100.95	818.61	−973.43	788.55	800
	2	66.598	606.34	−380.58	275.95	800
	3	130.09	304.58	14.377	99.80	800
Case 2	1	399.345	−122.12	770.46	-	1200
	2	287.116	80.04	700.48	-	1280
	3	−120.505	1525.99	−502.14	-	1280
	4	−19.121	898.76	−98.15	-	1280
	5	−95.029	1202.39	−352.16	-	1250
	6	191.750	224.86	524.04	-	1250

.

We applied the proposed INSGA-II algorithm proposed in this paper to the parallel chiller’s cooling capacity distribution problem. The improved NSGA-II algorithm was compared with the original NSGA-II, and other methods such as MOEA/D, IBEA and SPEA2 algorithms; the simulation and optimization results are shown as follows.

In Figure 5, each point represents the optimal solution that satisfies both Equations (15) and (16).

In

Table 7. Comparison of the optimization results of different algorithms in Case 1.

Load Requirements	NSGA-II	MOEA/D	IBEA	SPEA2	INSGA-II
2160	1589.4	1591.2	1574.9	1575.3	1565.3
1920	1419.4	1421.1	1413.4	1412.4	1406.5
1680	1250.6	1251.3	1247.9	1244.5	1242.5
1440	1007.4	1017.2	998.3	997.1	994.9
1200	923.1	945.8	862.2	854.6	842.7
960	742.8	781.4	715.3	712.0	697.8

and

Table 8. Comparison of the optimization results of different algorithms in Case 2.

Load Requirements	NSGA-II	MOEA/D	IBEA	SPEA2	INSGA-II
6850	4768.3	4819.5	4758.0	4754.6	4743.4
6470	4449.6	4453.4	4432.9	4431.1	4423.3
6090	4185.8	4189.9	4161.6	4152.4	4142.7
5710	3940.7	3951.0	3922.4	3918.2	3904.6
5330	3656.2	3682.5	3641.4	3638.3	3626.2

, the results of the INSGA-II algorithm are compared with the results of the NSGA-II, MOEA/D, IBEA and SPEA2 algorithms. Under the same load, the INSGA-II algorithm is more effective in saving energy, which shows the effectiveness of the INSGA-II algorithm in solving practical problems [24].

5. Conclusions

By analyzing the classical multi-objective evolutionary algorithm NSGA-II, it is found that the balance of convergence and diversity of the algorithm mainly depends on the strategy and method within the algorithm. We propose an improved NSGA-II algorithm which balances the global search ability of the algorithm in the preliminary period and the local exploitation ability in the later period by introducing adaptive parameters, and INSGA-II is able to improve the convergence speed and accuracy. At the same time, we adopt the Levy flight strategy to improve the global exploitation capability and a random walk strategy to improve the local search capability so that it can solve the problems of insufficient global search capability and the tendency to be trapped in a local optimum. To evaluate the performance of the INSGA-II algorithm proposed in this paper, it is compared with the algorithm before improvement and three other multi-objective evolutionary algorithms in the two-objective ZDT test set and three-objective DTLZ test set functions, as well as in practical engineering. The results show that the proposed algorithm has better convergence and diversity. In the future, we will analyze the complex Pareto front properties, design more efficient algorithms, and concentrate on the study of dynamic multi-objective optimization problems.