**1. Introduction**

In recent years, distributed generation systems using renewable energy technologies such as wind power generation and photovoltaic power generation have become one of the hotspots of research at home and abroad. After the distributed generation (DG) units are connected to the distribution network (DN), the structure, operation mode, and control strategy of the DN will undergo tremendous changes [1,2]. The research shows that DGs provide more flexibility and expansibility for distribution network. When DGs are connected to the distribution network (DN), the performances of DGs are most important for the power quality, reliability, and security of DN. In addition, DGs with renewable energy technology can effectively reduce system line loss and transmission congestion [3–5]. Despite the above advantages, if the placement and sizing of the DGs are improperly selected, it may cause a series of power system safety hazards such as power flow reverse, insufficient voltage stability, and malfunction of the protection device [6,7]. Therefore, the placement and sizing of the DGs has become one of the important subject.

The optimal planning of distributed generation sizing and siting is critical to ensure the operational performance of a distribution network in terms of power quality, voltage stability, reliability and profitability. DG planning problems are often defined as multi-objective and multi-constrained optimization problems. A number of review papers [8–10] have surveyed the optimization techniques for optimal DG planning in power distribution networks. The aforementioned review papers mainly focused on the discussion of various computational methods and metaheuristic algorithms.

Nowadays, industrial systems, such as papermaking, steelmaking, petrochemical, and power generation, are becoming more and more complex. In such cases, data-driven models based on novel nonlinear signal processing and data analysis techniques may provide an attractive alternative [11]. Therefore, it is paramount but challenging to develop effective techniques in modelling, monitoring, and control for complex industrial systems [12,13].

Among many optimization methods, multi-objective heuristic algorithms are the main means to solve multi-objective optimization problems today because they can effectively balance multiple objectives for optimal search [14]. At present, the main objectives of DG planning include the lowest investment cost of DG, the lowest environmental pollution, the optimal voltage quality of the power grid, and the minimum power loss. Partha Kayal et al. aimed to minimize the reactive power loss of the DN, to maximize the system voltage stability, and to establish a multi-objective optimization model for power systems [15]. Although it can effectively reduce the network loss and improve the voltage distribution, planning with only two objectives does not reflect the real condition. R. Li et al. established a multi-objective optimization model considering economy, voltage quality, and network loss and used the distance entropy multi-objective particle swarm optimization algorithm to solve the optimization problem [16]. However, there is a lack of environmental energy-saving factors in economic costs; M. Esmaili and T. Wang et al. also introduced voltage safety margins [17] and pollutant emissions [18] as optimization models, respectively. However, the objective function of optimization varies greatly in the timescale. For example, the voltage and reactive power change rapidly, and the pollutant emission often has a long timescale, which makes it difficult to use a unified model to optimize. Therefore, this paper establishes a multi-objective optimization model for economy, environmental protection, safety, and reliability.

Constraint processing is a key part of engineering optimization problems, and the treatment of constraints is the key to solving constrained multi-objective optimization problems. Commonly used constraint handling methods include rejection of infeasible solutions, penalty functions, and various correction algorithms. In fact, in the iterative process, infeasible solution is always difficult to be rejected, especially in the case of small feasible region, and repeated trial and error will affect the speed of solution; other modified algorithms are always designed for specific problems. The penalty function method [19] is the most classical method to deal with constraints, but the penalty factor has the same weight problem and is not easy to grasp, resulting in large errors in planning results.

In this paper, a new operator is proposed to address the degree of default. This method does not need to set parameters in advance but deals with feasible and infeasible solutions by considering the number and degree of violation of constraints. This method can effectively combine with multi-objective optimization algorithm to improve the diversity of populations and to avoid the algorithm falling into local optima.

At present, the planning methods of multi-objective optimization problems can be divided into weighted single-objective optimization algorithm and multi-objective optimization algorithm. In order to obtain the unique solution, scholars transform the multi-objective into single-objective optimization model by the weight method and solve it by combining with the single-objective optimization algorithm. This method makes the calculation convenient by reducing the dimension of the problem. For example, H. Su introduced improved simulated annealing particle swarm optimization [20] and established a multi-objective model to maximize DG utilization and to minimize system losses and environmental pollution. Finally, the DG injection model is optimized by this algorithm. However, the selection of weights depends on experience. Multi-objective functions are both interrelated and independent. The

subjectivity of weight selection makes the calculation result unsatisfactory. Aiming at the above problems, multi-objective optimization algorithm is introduced to solve the DG planning problem [21]. J. Neale solved the reconfiguration problem of the DN with the integrated DG by using non-dominated sorting genetic algorithm-II (NSGA-II) [22]. This method can effectively improve the operation performance of DN. A. Ameli used the multi-objective particle swarm optimization (MOPSO) [23] with adaptive grid method to optimize the distribution of energy supply system electric vehicle (EV) in DG technology. The results show the effectiveness of MOPSO programming; R. S. Maciel proposed an evolutionary particle swarm optimization algorithm (MEPSO) [24] as a multi-objective optimization algorithm for DG, of which the planning scheme can ultimately improve the objectives. However, in the multi-objective optimization algorithm, there is a common problem; that is, the evaluation information of solution set is insufficient due to crowding distance, which leads to the elimination of potential high-quality solutions in the truncation process. The distribution of the final algorithm planning results is uneven.

In addition, the mutation operator of NSGA-II will gradually lose the ability of local search as the dimension of solution variable increases, so the population updating strategy in the fireworks algorithm [25] is introduced to simultaneously improve the mutation operator and the diversity of solution set. Based on the effective Pareto-optimal set of multi-objective optimization problems, decision makers cannot get a unique solution. Therefore, this paper introduces an unbiased compromise strategy [26] to choose a best compromise from the Pareto-optimal set.

According to the above analysis, a lack of the operation state research planning of DN with the integrated DG considers the relay protection. Actually, when the DG is connected to the DN, the DN that is originally powered by a single system power becomes a network with multiple power structures, causing a change in the magnitude and direction of the short-circuit current during the failure. The action affects the safe and stable operation of the DN. Therefore, this paper increases the short-circuit current constraint condition and compares the result with the optimization result without considering the short-circuit constraint, so as to prove the necessity of considering the relay protection condition. At present, with the development of technology, distributed energy resources are bound to penetrate into the fields of industry, commerce, and urban and rural residents. The problem of DG planning with timing characteristics needs to be solved urgently. Therefore, the paper studies the planning of various examples and plans the static state and dynamic running state of the IEEE 33-bus system separately to provide a reasonable DG planning for decision makers.

The main contributions of this work are as follows:


The proposed algorithm is applied to solve the static and dynamic planning of the DN with the integrated DG units. The rest of this article is organized as follows. Section 2 presents the mathematical programming model. Section 3 briefly introduces the algorithms needed in the planning scheme and proposes INSGA-II. Section 4 provides the comparison of the simulation results of the proposed algorithm in various cases and the performance of the algorithm. Section 5 concludes our work.

#### **2. Problem Formulation**

#### *2.1. Objective Functions*

Three objectives should be taken into account when establishing multi-criteria optimization model of DN with the integrated DG units: (1) Improving the energy-saving benefit of the system; (2) reducing system line losses; and (3) reducing the node voltage deviation.

#### 1. Maximizing annual energy-saving benefit

Maximizing annual energy-saving benefit of DN with the integrated DG units can be expressed as follows:

$$\max f\_1 = Z\_{\text{cost}}^{\text{NODG}} - Z\_{\text{cost}}^{\text{DG}} \tag{1}$$

where *Z*NODG cost is the total annual cost of DN without DGs and *<sup>Z</sup>*NODG cost is the total annual cost of DN with DGs.

Total annual costs of DN excluding DGs can be expressed as follows:

$$\mathbf{Z}\_{\rm cost}^{\rm NODG} = \mathbf{C}\_{\rm loss} + \mathbf{C}\_{\rm b} \tag{2}$$

where *C*loss is the annual loss of DN without DGs and *C*<sup>b</sup> is the annual purchase cost of DN without DGs.

$$\mathcal{C}\_{\text{loss}} = \sum\_{j=1}^{k} \left( \mathbb{C}\_{\text{p}} \cdot \tau\_{\text{max}} \cdot \Delta P\_{j} \right) \tag{3}$$

where *k* is the number of branches in the DN; *C*<sup>p</sup> is the unit price of electricity consumed per unit (\$/kWh); τmax is the annual maximum load loss hour (h) of the DN; and Δ*Pj* is the active power loss (kW) of the *j*th branch.

$$\mathcal{C}\_{\mathbf{b}} = \mathcal{C}\_{\mathbf{p}} \cdot P\_{\text{load}} \cdot T\_{\text{max}} \tag{4}$$

where *P*load is the DN total active load and *T*max is the DN annual maximum load utilization hours.

The total annual cost of DN with DGs can be expressed as follows:

$$\mathbf{Z}\_{\rm cost}^{\rm DG} = \mathbf{C}\_{\rm dgm} + \mathbf{C}\_{\rm ploss} + \mathbf{C}\_{\rm B} - \mathbf{C}\_{\rm sub} \tag{5}$$

where *C*dgm is the annual investment and maintenance cost of DGs; *C*ploss is the annual cost of DN with DGs (calculated in the same way as *C*loss); *C*<sup>B</sup> is the annual purchase cost of DN with DGs; and *C*sub is the financial subsidy of new energy generation.

$$\mathbb{C}\_{\text{dgm}} = \sum\_{i=1}^{n} \left( \frac{r(1+r)^t}{(1+r)^t - 1} \cdot \mathbb{C}\_{d\text{g}i} + \mathbb{C}\_{mi} \right) \cdot P\_{d\text{g}i} \tag{6}$$

where *r* is the annual interest rate; *t* is the planning period; *Cdgi* is the *i*th distributed generation equipment investment cost (\$/kWh); *Cmi* is the annual operation and maintenance cost (\$/kWh) for the *i*th distributed generation; and *Pdgi* is the installed capacity (kW) of the *i*th distributed generation.

$$\mathcal{C}\_{\rm B} = \mathcal{C}\_{\rm P} \cdot \left( P\_{\rm load} - \sum\_{i=1}^{n} \lambda\_i \cdot S\_{d\rm gi} \right) \cdot T\_{\rm max} \tag{7}$$

In order to reflect the benefits of DGs in environmental protection, the government has put forward policy support to DGs, that is, financial subsidies for distributed generation. It can be expressed as follows:

$$\mathcal{C}\_{\text{sub}} = \mathcal{C}\_{\text{sp}} \cdot \sum\_{i=1}^{n} P\_{d\text{ç}i} \cdot T\_{\text{max}} \tag{8}$$

where *C*sp is the subsidy amount (\$) of distributed generation unit.

2. Minimize line losses

The existence of line losses will lead to line heating, which will accelerate the aging of insulated lines, will reduce the insulation of lines, and will ultimately lead to the risk of leakage. Reducing line losses can improve energy efficiency of electrical equipment or processes.

Lowering the line losses can improve the energy utilization efficiency of the energy-using equipment or process, and it is also one of the important measures of energy saving. The objective function can be expressed as follows [27]:

$$\min f\_2 = \sum\_{i,j \in n} g\_{ij} \left( \mathcal{U}\_i^2 + \mathcal{U}\_j^2 - 2\mathcal{U}\_i \mathcal{U}\_j \cos \theta\_{ij} \right) \tag{9}$$

where *n* is the total number of DN nodes; *gij* is the admittance of the branch (*i*, *j*); *Ui* and *Uj* are the voltage amplitudes of branches *i* and *j* respectively; and θ*ij* is the voltage phase angle.

3. Minimize the total voltage deviation

With the increase of load, the system voltage stability will deteriorate gradually, even voltage collapse, resulting in a system in danger. Therefore, the voltage deviation is one of the important indexes to evaluate the operation safety and power quality of the system. In addition, the increase of network node voltage can effectively reduce the reactive power loss of the system. The objective function can be expressed as follows:

$$\min f\_3 = \sum\_{i=0}^{N} \frac{\left| \mathcal{U}\_i - \mathcal{U}\_i^{specific} \right|}{\mathcal{U}\_i^{\max} - \mathcal{U}\_i^{\min}} \tag{10}$$

where *Ui* is the node *<sup>i</sup>*th voltage real value of the DN when DG is incorporated and *<sup>U</sup>speci fied <sup>i</sup>* is the voltage rated value. This paper assumes that the rated voltage is 1 (p.u).

#### *2.2. Influence of DG on Relay Protection of DN*

DN is generally equipped with three-stage current protection. After the DG is connected to the DN, the system changes from single power supply to multiple power supply. When the transmission line is short-circuited, the magnitude and direction of the short-circuit current will change.

Figure 1 shows a typical simple DN. L1 and L2 represent feeders; T1 represents a transformer; PD1, PD2, PD3, and PD4 are protective devices of the corresponding feeders; and DG is connected to the DN from node C.

**Figure 1.** Distribution network system diagram of distributed generation (DG).

Assuming that the DG is connected, the short-circuit fault can be divided into the following three cases:

1. The adjacent feeder connected to the DG is short-circuited

When the short-circuit fault F1 occurs at the end of the adjacent feeder, the power source *S* and DG provide the overlapping short-circuit current to PD1 and PD2. Compared without DG connection, the short-circuit current increases and its value may be larger than the original I segment setting value of PD1, which makes PD1 malfunction, resulting in the loss of coordination between PD1 and PD2. Simultaneously, PD3 will also flow through the reverse current provided by DG. When the capacity of DG is larger, the value of reverse current may be larger than the original I segment setting value of PD3, resulting in PD3 malfunction.

2. The upstream of DG access is short-circuited

When the DG upstream feeder AC is short-circuit fault F2, PD4 can operate normally. However, when PD3 is activated, the island downstream of the DG will form an island operation. At the same time, the DG will always provide short-circuit current to the fault point, which will affect the reliability of the system.

3. The downstream of DG access is short-circuited

When DG downstream feeder CD short-circuit fault F3 occurs, protections PD1 and PD2 can operate normally. Although the short-circuit current through PD3 is only supplied by the system power *S*, the fault current will be smaller than before due to DG connection, which may cause the backup protection of PD3 to refuse to operate. The short-circuit current through PD4 is provided by the system power *S* and DG together. The increase of short-circuit current may increase the protection distance of PD4, thus losing the selectivity.

When there are multiple DGs in the system, the analysis method is similar. From the above analysis, it can be seen that the access of DG will have adverse effects on the reliability of the relay protection and system of DN, so it is necessary to consider the restraint of relay protection when optimizing the configuration of DG. In the traditional three-stage current protection, the I-stage instantaneous current quick-break protection is the most important, so this paper focuses on the analysis of the impact of DG on the I-stage of the original current protection of the PD in the DN.

#### *2.3. Constraint Condition*

DG connected to a distribution network has great influence on power flow, voltage distribution, branch current, and power quality of the system. In order to design a reasonable programming model, the following equality constraints and inequality constraints are included in this paper.

1. Equality constraint

Below are the power flow equations constraints.

$$\begin{cases} \begin{aligned} P\_{D\bar{G}i} - P\_{Li} - \mathcal{U}\_i \sum\_{\substack{k \in I\\k \neq i}}^{N\_1} \mathcal{U}\_k (G\_{\bar{k}} \cos \theta\_{\bar{i}k} + B\_{\bar{i}k} \sin \theta\_{\bar{i}k}) &= 0\\ Q\_{D\bar{G}i} - Q\_{Li} - \mathcal{U}\_i \sum\_{k \in I} \mathcal{U}\_k (G\_{\bar{k}} \sin \theta\_{\bar{i}k} + B\_{\bar{i}k} \cos \theta\_{\bar{i}k}) &= 0 \end{aligned} \tag{11}$$

where *PDGi* and *QDGi* are the *i*th node active and reactive power of the DGs injected; *PLi* and *QLi* are the *i*th node active and reactive powers of load output; *Uk* is the voltage amplitude of all *k*th nodes connected to *i*th nodes; *Gik* is branch conductance; *Bik* is the branch susceptance; θ*ik* is the difference between the voltage angles of *i*th node and *k*th node; and *N*<sup>1</sup> is the number of branches associated with the *i*th node.

2. Inequality constraints

In order to make the planning scheme meet the power system operation standards, it is usually required that the node voltage, current, power, and permeability satisfy the constraints after DGs are connected to the grid.

*Processes* **2019**, *7*, 955

a. Node voltage constraints

$$
\mathcal{U}\_i^{\text{min}} \le \mathcal{U}\_i \le \mathcal{U}\_i^{\text{max}} \tag{12}
$$

where *Ui* represents the voltage amplitude of the *i*th node, and *U*min *<sup>i</sup>* and *<sup>U</sup>*max *<sup>i</sup>* represent the upper and lower limits of the *i*th nodes voltage amplitude, respectively.

b. Branch current constraints

$$I\_{ij} \le I\_{ij}^{\max} \tag{13}$$

where *Iij* is the current of branch (*i*, *j*) and *I* max *ij* is the maximum current allowed of branch (*i*, *j*). c. Single DG access capacity constraints

$$\begin{cases} \begin{array}{c} P\_{DGi}^{\min} \le P\_{DGi} \le P\_{DGi}^{\max} \\ Q\_{DGi}^{\min} \le Q\_{DGi} \le Q\_{DGi}^{\max} \end{array} \end{cases} \tag{14}$$

where *PDGi* and *QDGi* are the active power and reactive power of the DG connected to the *i*th node, respectively.

d. Distributed generation access total capacity constraints

$$\sum\_{i}^{N} S\_{DGi} \le 0.3 S\_L \tag{15}$$

where *SDGi* is the total capacity of distributed generation and *SL* is the maximum load value of DN.

⎧

⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩

e. Short-circuit current constraints

$$\begin{array}{lcl} K\_{\text{set},i}^{\text{III}} = \frac{I\_{k,i+1,\text{min}}^{(3)}}{I\_{\text{set},i}^{\text{III}}} \ge 1.2\\ I\_{\text{set},i}^{1} > I\_{k,i.\text{max}}^{(3)}\\ I\_{\text{set},i}^{1} > I\_{r k,i-1.\text{max}}^{(3)}\\ I\_{\text{set},i}^{\text{II}} > I\_{r k,i-1.\text{max}}^{(3)}\\ I\_{\text{set},i}^{\text{III}} > I\_{r k,i-1.\text{max}}^{(3)}\\ I\_{\text{set},i}^{1} > I\_{r k,i-1.\text{max}}^{(3)}\\ \end{array} \tag{16}$$

where *K*III *sen*,*<sup>i</sup>* is the sensitivity coefficient of current III segment protection of branch *i* as backup protection; *I* I *set*,*i* , *I* II *set*,*i* , and *I* III *set*,*<sup>i</sup>* are the setting values of the current I, II, and III protections of branch *i*, respectively; *I* I *set*,j is the setting value of the current I segment protection of the branch *j* of the feeder line of the DG; *I* (3) *<sup>k</sup>*,*i*−1.max is the maximum reverse current through branch *i* when three-phase short circuit occurs at the end of the upstream of branch *i*; *I* (3) *<sup>k</sup>*,*i*.max and *I* (3) *<sup>k</sup>*,*j*.max are the maximum short-circuit currents through branch *i* of the terminal three-phase short circuit of branch *i* and branch *j*, respectively; and *I* (2) *<sup>k</sup>*,*i*+1.max is the minimum short-circuit current through branch *i* of the two-phase short circuit at the end of the downstream of branch *i*.

3. Overview formulation

From the above analysis, the multi-objective function and constraints can be described as follows:

$$\min \left[ f\_1(\mathbf{x}\_{\mathcal{S}}, \mathbf{x}\_{\mathcal{C}}), f\_2(\mathbf{x}\_{\mathcal{S}}, \mathbf{x}\_{\mathcal{C}}), \dots, f\_{\mathcal{N}\_{obj}}(\mathbf{x}\_{\mathcal{S}}, \mathbf{x}\_{\mathcal{C}}) \right] \tag{17}$$

$$\text{s.t.} \quad h(\mathbf{x}, \mathbf{x}) = 0,\\ \mathbf{i} = 1, \dots, p \tag{18}$$

$$g\_i(\mathbf{x}\_{\mathbf{s}}, \mathbf{x}\_{\mathbf{c}}) \le 0, \ i = 1, 2, \cdots, q \tag{19}$$

where *fi* is the *i*th objective function, *Nobj* is the number of objective functions, *xs* and *xc* denote the state vector and the control vector respectively, *hi* is the equation constraint, *p* is the number of equality constraints, *gi* is the inequality constraint, and *q* is the number of inequality constraints.

*xc* is composed of *n* variables, and *n* represents the number of nodes of the optimized network, where each variable represents two components: the DG installation placement and capacity. For example, a variable is assigned a value indicating the placement of the variable to installation DG. The value of the variable indicates the installation capacity. If DG is not installed, the corresponding variable value is 0. *xc* can be illustrated as follows:

$$\mathbf{x}\_{\mathbb{C}} = \left[ (L\_{\text{DG1}}, P\_{\text{DG1}}), (L\_{\text{DG2}}, P\_{\text{DG2}}), \dots, (L\_{\text{DGN}}, P\_{\text{DGN}}) \right] \tag{20}$$

where *LDGi* indicates that the *i*th node is the DG installation location and *PDGi* represents the DG installation capacity of the *i*th node.

*xs* is the network parameter that *xc* is calculated from the power flow. Variables such as line losses, node voltage, and branch current are used to calculate multi-objective functions and to determine the satisfaction of constraints. *xs* can be illustrated as follows:

$$\mathbf{x}\_s = [P', \mathbb{Q}', \mathbb{U}', \mathbb{I}', \theta'] \tag{21}$$

where *P* , *Q* , *U* , θ , and *I* represent the network feeder active power vector, reactive power vector, node voltage vector, voltage argument vector, and branch current vector respectively after DG is connected.

#### *2.4. Establish PQ Mode of DG*

The DG injects its power output into the grid through power electronics. Typically, for a PQ model, which shows the active power (P) versus reactive power (Q) called a PQ model. DG is modeled as a constant power factor and negative load model.

In this case, the DG is modeled as a constant power source. *PDG* is the actual power output specified for the DG model. The load at *i*th node with the DG device is modified as Equations (22)–(23).

$$P'\_{\text{load}i} = P\_{\text{load}i} - P\_{\text{DGi}} \tag{22}$$

$$Q\_{DGi} = P\_{DGi} \cdot \tan\left(\cos^{-1}(\varphi)\right) \tag{23}$$

where *Ploadi* is the original network load power of *i*th node, *PDGi* and *QDGi* are the active and reactive power of DG at the *i*th node, *P loadi* is the active power after installing DG, and cosϕ is the power factor.

In general, constrained problems can be solved using either deterministic or stochastic algorithms. However, deterministic approaches such as feasible direction and generalized gradient descent require strong mathematical properties of the objective function such as continuity and differentiability. In cases where these properties are absent, evolutionary computation, such as NSGA-II, offers reliable alternative methods [13].

#### **3. Improved NSGA-II Algorithm**

#### *3.1. NSGA-II Algorithm Overview*

NSGA-II was proposed by Deb, K. in 2002 [28], which is one of the most popular multi-objective optimization algorithms. It uses simulated binary crossover and polynomial variation to introduce non-dominated sorting and crowding distance operator instead of NSGA. The sharing radius of the algorithm ensures the diversity of the population and the uniformity of the Pareto-optimal set and introduces the elite retention strategy and the elimination strategy, which improves the operation speed and stability of the algorithm.

#### *3.2. Dominant, Non-Dominated, Pareto-Optimal Set and Crowding Distance*

NSGA-II introduces Pareto-optimal set and crowding distance to replace the fitness of traditional intelligent optimization algorithm. The solution set is divided into dominated solution set and non-dominated solution set according to the dominant relationship among the solutions. The rank of the solution is from 1 to n: the better the quality, the smaller the value. Solutions of the same rank are divided into pros and cons by comparing the size of the crowding distance.

Multi-objective optimization problems can be described as follow:

$$\min f\_i(\mathbf{x}), \ i = 1, 2, \dots, N\_{obj\_i}, \mathbf{x} \in \mathbf{I} \tag{24}$$

where I is a feasible solution space.

**Definition 1.** *If solution x*<sup>1</sup> dominates *x*2*, then the following definitions must be met.*

$$\begin{array}{c} \forall i, j \in \{1, 2, \dots, N\_{obj}\} \\ \exists\_{i \neq f} f\_i(\mathbf{x}\_1) < f\_i(\mathbf{x}\_2) \lor f\_j(\mathbf{x}\_1) = f\_j(\mathbf{x}\_2) \end{array} \tag{25}$$

**Definition 2.** *If solution x*<sup>1</sup> does not dominate *x*2*, then the following definitions must be met.*

$$\begin{array}{c} \forall i, j \in \{1, 2, \dots, N\_{obj}\} \\ \exists f\_i(\mathbf{x}\_1) \le f\_i(\mathbf{x}\_2) \land f\_j(\mathbf{x}\_1) > f\_j(\mathbf{x}\_2) \end{array} \tag{26}$$

Equations (25) and (26) in Definitions 1 and 2 are derived from Reference [29].

**Definition 3.** *Pareto-optimal set is composed of mutually non-dominated solution sets, and ranks are composed of mutually dominated solution sets.*

**Definition 4.** *Assume that P* = {*x*1, *x*2, ··· , *xn*} *is a Pareto-optimal set and that the crowding distance represents the distribution density of other solution sets around the solution. The larger the solution distance at the same layer, the better the solution distribution, that is, the better the solution set diversity. Its calculation formula is as follows:*

$$cd\_{i}^{k} = \begin{cases} \frac{f\_{i+1}^{k} - f\_{i-1}^{k}}{f\_{\max}^{k} - f\_{\min}^{k}}, if \text{ } ind \text{ex} \{ \mathbf{x}\_{i}^{k} \} \in [2, n - 1] \\\ \infty, \qquad otherwise \end{cases} \tag{27}$$

$$cd\_i = \sum\_{k=1}^{m} cd\_i^k \tag{28}$$

*where cdk <sup>i</sup> is the crowding distance on the kth objective function of xi, m is the number of objective functions, index (xk i ) is the sort index of xi on the kth objective function, f <sup>k</sup> <sup>i</sup>*+<sup>1</sup> *is the value of objective function corresponding to the last solution on the axis of the kth objective function, and f <sup>k</sup> max and f <sup>k</sup> min represent the maximum and minimum values of the kth objective function, respectively.*

#### *3.3. INSGA-II: Improved Mutation Operator*

Because mutations in genetics often lead to worse results, the probability of mutation operations in genetic algorithm (GA) is usually set to be very small. In GA, crossover and mutation operations have two functions: global search and local search. When dealing with low-dimensional convex problems, only the crossover operator can solve these problems well. However, when dealing with high-dimensional nonconvex problems, it often falls into the local optimal solution. For this reason, the mutation operator is improved by referring to the population-updating method in the fireworks algorithm [25]; the new mutation individual is created by Equations (29)–(31).

$$d\_i = vd \cdot rand \tag{29}$$

$$S = randperm(1, d\_i) \tag{30}$$

$$u'\_j = \begin{cases} u\_j^{\text{num}}, \text{if } j \in S \\ u\_{j\prime} \text{otherwise} \end{cases} \tag{31}$$

where *di* denotes the dimension of the *i*th solution for mutation, *vd* denotes the variable dimension, *rand* denotes the random number between [0, 1], *S* denotes the index for performing the mutation operation, *umu <sup>j</sup>* denotes the *<sup>j</sup>*th gene performing the mutation operation, and *<sup>u</sup> <sup>j</sup>* represents the value of the *j*th gene after the mutation operation.

Crossover and mutation operations often produce solutions beyond the feasible space. The general method is to assign solutions less than the lower limit of feasible region to the lower limit, and the solution that exceeds the upper limit of the feasible region is given the upper limit value. This operation is simple to perform but reduces the diversity of solution sets. According to Equation (32), solutions beyond the boundary can be effectively mapped back to the feasible space.

$$\boldsymbol{\mu}\_{i}^{\wedge} = \boldsymbol{\mu}\_{i}^{\text{min}} + |\boldsymbol{\mu}\_{i}|^{\diamond} \big( \boldsymbol{\mu}\_{i}^{\text{max}} - \boldsymbol{\mu}\_{i}^{\text{min}} \big) \tag{32}$$

where *u*max *<sup>i</sup>* and *<sup>u</sup>*min *<sup>i</sup>* represent the upper limit and lower limit of the gene, respectively, and % is the remainder operation.

#### *3.4. INSGA-II: Violation Constrained Index*

In view of the subjectivity of penalty factor selection, the violation constrained index (VCI) is proposed when dealing with infeasible solutions.

The method divides the solution set into feasible solution and infeasible solution by calculating the error degree (ED) of the constraint condition of the solution.

The VCI calculation formula is as follows:

$$\mathbb{E}\,ED(\mathbf{x}\_{s\prime},\mathbf{x}\_{\varepsilon}) = \left| h(\mathbf{x}\_{s\prime},\mathbf{x}\_{\varepsilon}) \right| + \max(g(\mathbf{x}\_{s\prime},\mathbf{x}\_{\varepsilon}),0) \tag{33}$$

$$VCI\_i = \sum\_{j}^{\mathbb{C}} \frac{ED\_j^i - ED\_j^{\min}}{ED\_j^{\max} - ED\_j^{\min}} \tag{34}$$

where *ED*max *<sup>j</sup>* and *ED*min *<sup>j</sup>* represent the maximum and minimum error degree of infeasibility under the *j*th constraint, respectively, and *C* indicates the number of constraints. The classification process of the solution set is performed according to Equation (35).

$$S\_{-}t = \begin{cases} \begin{array}{c} 0 \ \text{\$ } \/ \text{ if } \text{all } \, VCI\_{i} = 0 \\ 1 \ \text{ } \/ \text{ if } \text{all } \, VCI\_{i} \, ! = \, 0 \\ \begin{array}{c} \frac{1}{2} \ \text{ } \/ \text{ otherwise} \end{array} \end{cases} \tag{35}$$

where *S*\_*t* denotes the type of the solution set. The default index can be effectively used in combination with the non-dominated sorting strategy. For example, when the *S*\_*t* is 0 or 1/2, the feasible solution part is given the true rank and the crowding distance, and the rank of the infeasible solution is sorted in descending order according to *VCIi* and continues to be sorted from the feasible solution.
