*3.5. INSGA-II: Improved Crowding Distance Operator*

Figure 2 shows that, after solutions a–h pass the truncation strategy, a, d, e, g and h are preserved but the eliminated solution c has a better crowding distance than the solution d, which is due to the traditional crowding distance information being missing, and the truncation strategy eliminates the solution set with small crowding distances at one time, resulting in uneven distribution of the final set. The calculation of the traditional crowding distance only considers the spatial distribution of two adjacent solutions on the axis of the objective function without considering that the effect of deleting the neighborhood solution leads to the elimination of potential high-quality solutions. Finally, a potential crowding distance is proposed. The potential crowding distance in the two-dimensional space is simple to calculate. As shown in Figure 2a, when the solution map is re-target functions 1 and 2, the neighborhoods are the same in the solution set and only the order changes. The truncation process is shown in Figure 2b,c:

(**a**) Individual distribution chart (**b**) Eliminating points chart (**c**) Truncation process chart

**Figure 2.** Pareto-optimal set of truncation strategy.

Therefore, the formula for calculating the potential crowding distance in two-dimensional space is as follows:

$$\begin{cases} \Delta d\_i^1 = \frac{f\_{i-1}^1 - f\_{i-2}^1}{f\_{\max}^1 - f\_{\min}^1} + \frac{f\_{i+2}^2 - f\_{i+1}^2}{f\_{\max}^2 - f\_{\min}^2}, \\ \Delta d\_i^2 = \frac{f\_{i+2}^1 - f\_{i+1}^1}{f\_{\max}^1 - f\_{\min}^1} + \frac{f\_{i-1}^2 - f\_{i-2}^2}{f\_{\max}^2 - f\_{\min}^2}, \\ \end{cases} \tag{36}$$

where Δ*d*<sup>1</sup> *<sup>i</sup>* and <sup>Δ</sup>*d*<sup>2</sup> *<sup>i</sup>* denote the increments of the crowding distance after the adjacent solution of the *k*th solution is deleted, *f* <sup>1</sup> *<sup>i</sup>*−<sup>1</sup> and *<sup>f</sup>* <sup>2</sup> *<sup>i</sup>*−<sup>1</sup> denote the neighborhood solutions mapped to the *k*th solution on the objective function 1, and *pdi* denotes the potential crowding distance of the *k*th solution.

The objective function above three-dimensions does not have the symmetric relationship of the two-dimensional objective function. Each solution no longer has only the same two adjacent solution sets on each objective function axis but has (2–n) different solutions. Therefore, the potential crowding distance of the solution under the 3D objective function is calculated as follows:

$$pd\_i = cd\_i + \max(k\_{\text{\\_}} \Delta d\_i), k \in \Phi \tag{37}$$

$$k\Delta d\_{i} = \sum\_{j=1}^{m} 1^{j} \left[ 1^{l} \left( \frac{f\_{i-1}^{j} - f\_{i-2}^{j}}{f\_{\text{max}}^{j} - f\_{\text{min}}^{j}} \right) + 1^{l} \left( \frac{f\_{i+2}^{j} - f\_{i+1}^{j}}{f\_{\text{max}}^{j} - f\_{\text{min}}^{j}} \right) \right] \tag{38}$$

where Φ denotes the set of adjacent fields mapped to the *i*th solution on each objective function axis; *m* denotes the number of objective functions; *k*Δ*di* denotes the increment generated by the crowding distance of the *i*th solution after the *k*th solution in the neighborhood set is eliminated; and 1*<sup>j</sup>* , 1*<sup>l</sup>* , and 1*<sup>r</sup>* respectively indicate that the *k*th solution belongs to the neighborhood of the *i*th solution or belongs to the left neighborhood solution or the right neighborhood solution on the *j*th objective function.

#### *3.6. INSGA-II: Improved Crowding Distance*

NSGA-II uses the tournament selection operator for population reproduction, and the selection criteria considers the rank and crowding distance. The specific selection process is described in Reference [27]. Since the improved NSGA-II algorithm increases the potential crowding distance, it needs to involve the crowding distance for the selection operation. After some modifications, the selection rules are as follows:

$$new\_k = \begin{cases} \begin{array}{l} \mathbf{x}\_i, if \ cd\_i > cd\_j \text{ and } pd\_i > pd\_j\\ \mathbf{x}\_j, if \ cd\_j > cd\_i \text{ and } pd\_j > pd\_i\\ \mathbf{x}\_i \text{ or } \mathbf{x}\_j, otherwise \end{array} \tag{39}$$

where *newxk* denotes the selected *k*th descendant solution, and *xi* and *xj* respectively denote two different solutions in the parent. Improved selection strategy comprehensively considers the current crowding distance and potential crowding distance. Therefore, it can select a better solution and can avoid the defects of insufficient diversity of the solution set when the truncation strategy is executed.

#### *3.7. Unbiased Compromise Strategy*

In order to solve the cumbersome scheme of the Pareto-optimal set, this paper uses an unbiased compromise strategy based on fuzzy set, which uses unbiased member parameter ω to evaluate the quality of each individual in the Pareto-optimal set. The calculation formula is as follows:

$$\boldsymbol{\omega}^\* = \max\_{i=1,\ldots,m} \left( \frac{\sum\_{j=1}^m \alpha\_j^i}{\sum\_{i=1}^m \sum\_{j=1}^m \alpha\_j^i} \right) \tag{40}$$

where ω*<sup>i</sup> <sup>j</sup>* denotes the unbiased parameter of the *i*th solution in the Pareto-optimal set on the *j*th objective function; *f* min *<sup>j</sup>* and *<sup>f</sup>* max *<sup>j</sup>* denote the minimum and maximum values of the *j*th objective function corresponding to the Pareto-optimal set, respectively; *f<sup>i</sup> <sup>j</sup>* denotes the value of the *j*th objective function corresponding to the *i*th solution in the Pareto-optimal set; *n* denotes the number of objective functions; *m* denotes the number of solutions in the Pareto-optimal set; and ω<sup>∗</sup> denotes the unbiased member parameters of the comprehensive optimal solution.

#### *3.8. Optimized Configuration Process Based on Improved NSGA-II*

The proposed algorithm optimizes the placement and sizing process configuration of DG, as follows:

**Algorithm 1. The algorithm of the proposed INSGA-II**

Input: Set IEEE-33 distribution network-related parameters, the number of target functions, the number of variables, the size of the population, the maximum number of iterations, etc.

Outputs: Comprehensive optimal solution *x*<sup>∗</sup> .

<sup>1:</sup> Initialization population *P*(*xc*)*t*, set the number of iterations t = 0;

<sup>2:</sup> *P*(*xs*)*<sup>t</sup>* is calculated by forward and backward substitution method, and then, the objective function values *f*1, *f*2, *f*3, and *S*\_*t* are calculated;

<sup>3:</sup> Algorithm iteration start, while *t* < *t*max do;

<sup>4:</sup> The solution set is classified by *S*\_*t*, feasible solution to perform non-dominated sorting strategy, and the infeasible solution to calculate *DCI* and sorted;

<sup>5:</sup> Generation of children *Q*(*xc*)*<sup>t</sup>* from the parent *P*(*xc*)*<sup>t</sup>* by improving selection operations and improving genetic manipulation;

<sup>6:</sup> *Q*(*xc*)*<sup>t</sup>* performs power flow calculation, mixes with *P*(*xc*)*<sup>t</sup>* to form a mixed population *R*(*xc*)*<sup>t</sup>* and *R*(*xc*)*<sup>t</sup>* for non-dominated sorting strategy, and calculates the crowding distance according to Equations (37) and (38);

<sup>7:</sup> Selecting a new parent population *R*(*xc*)*<sup>t</sup>* of size from *P*(*xc*)*t*+<sup>1</sup> by a truncation strategy;

<sup>8:</sup> *t* = *t* + 1;

<sup>9:</sup> End while;

<sup>10:</sup> Using unbiased compromise strategy to select comprehensive optimal solution *x*<sup>∗</sup> ;

<sup>11:</sup> Return *x*<sup>∗</sup> .

#### **4. Experiment and Results**

#### *4.1. Experiment Setup and Description*

In order to better verify the effectiveness of the proposed algorithm, IEEE 33-, 69-, and 118-bus systems are introduced as simulation examples to compare performance with current mainstream multi-objective optimization algorithms, considering the dynamic and static characteristics of the DGs and DN. Hardware parameters used in the experiment are as follows: Intel (R) Core (TM) i5-3337 CPU 1.80 GHz, memory: 4.00 GB, and simulation software: Matlab-2014a (The MathWorks, Inc.3 Apple Hill Drive, Natick, MA, USA, 2014).

#### *4.2. IEEE 33-Bus System Case Simulation Experiment*

An IEEE 33-bus system topology diagram is shown in Figure 3, where the system includes 33 nodes and 32 branches [30]. The 0 node is assumed to be a balanced node with a voltage reference of 12.66 kV, total power load of 3.17 MW, and reactive power 2.30 MVAr.

**Figure 3.** IEEE 33-bus system topology diagram.

Multi-objective function parameter setting: The maximum installation capacity of each node of DG is 1 MW, the annual maximum load utilization time and the annual maximum line loss hour of DGs are 4500 h, the investment cost of micro-turbine (MT) unit equipment is 0.07 million \$/kW, annual operation and maintenance cost is 0.009 \$/kWh, wind turbine generator (WTG) group unit equipment investment cost is 0.114 million \$/kW, annual operation and maintenance costs are 0.004 \$/kWh, MT and WTG unified as a PQ model for processing, power factor is 0.9, the unit price of the loss is 0.071 \$/kWh, the financial subsidy per unit of electricity is 0.019 \$/kWh, and the conversion factor of the annual equipment investment cost of the distributed power source is obtained by 3% of the annual interest rate. We use the algorithm of the proposed INSGA-II for the following series of analysis.

#### 4.2.1. Analysis of Optimization Results with Different Parameters

As shown in Table 1, in order to improve the operation speed of the algorithm and the accuracy of the solution, experimental comparisons of the experimental hyperparameters are made. Compared with experiment 1 and experiment 2, the optimization results are improved with the increase of iteration times.


**Table 1.** Simulation results of INSGA-II with different parameters.

Comparing experiments 1, 3, and 4, the optimization effects are improved with the increase of population size but the time increases. The population selection of 100 and the number of iterations of 100 are taken as the hyperparameters of the algorithm.

Through large numbers of experiments, subsequent experiments decided to set improved NSGA-II parameters: population size is 100, number of iterations is 100, and crossover probability is 0.7.

#### 4.2.2. Analysis of DG Installation Capacity Optimization Results

The DG optimization configuration results of the improved NSGA-II are shown in Table 2.


**Table 2.** INSGA-II optimization results.

As shown in Table 2, Case-1 has no DG, Case-2 is installation without short-circuit current constraints (other constraints are considered), and Case-3 is installation with short-circuit current constraints (other constraints are considered). Compared with Case-1, the planning of Case-2 and -3 obtained by INSGA-II can effectively improve the objective function. Case-2 can produce energy-saving benefits of about \$0.197 million, the voltage deviation is improved by 73.24%, and the line loss is reduced by 69.60%. Case-3 can produce energy-saving benefits of about \$0.193 million, the voltage deviation is improved by 74.83%, and the line loss is reduced by 60.49%. Comparison of Case-2 and -3 show that, under the influence of short-circuit current protection constraints, the capacity of DGs decrease and voltage deviation and line losses are improved. In addition, when the DGs are close, the short-circuit current of some branches will be very high, so the Case-3 optimization results do not show that the DG access points are close. The above results confirm the rationality of considering short-circuit current constraints.

As shown in Figure 4a,b, Case-2 and -3 can improve the voltage amplitude of each node and line losses by installing DG appropriately. Because most of the installation nodes of Case-3 are terminal nodes of the system, the supporting effect of the node voltage is stronger [31]. Active network loss is usually determined by node voltage and branch resistance. Distribution network voltage level is lower and R/X value is larger, so it will lead to larger network losses. When DG is reasonably connected, the voltage fading will be improved [31]. In the case of the conditions, Case-3 has a greater effect on node voltage support than Case-2, so the line losses are better.

(**a**) The voltage amplitude of each node (**b**) The line loss of each node

**Figure 4.** Results of nodal voltage optimization.

### *4.3. Performance Analysis of INSGA-II Algorithm*

In order to verify the optimization performance of the proposed algorithm, INSGA-II, NSGA-II, and MOPSO are used to analyze and compare the planning results of the IEEE 33-bus system

independently. The parameters of NSGA-II are set as INSGA-II, the inertia factor is 0.8, the local speed factor is 0.1, and the global speed factor is 0.1; other parameters are the same as INSGA-II.

As shown in Figure 5, NSGA-II and MOPSO have uneven distribution of Pareto solutions due to the defects of crowding distance operators. MPSO optimizes by choosing leaders and by archiving mechanism. Similar to NSGA-II, its selection process also depends on a dominant relationship, so crowding distance is also needed to participate in it. When solving the high-dimensional nonconvex problem, the solution set is easy to fall into the local optimal solution and the population diversity is insufficient. Obviously, the uniformity of Pareto-optimal set obtained by INSGA-II is the best of the three algorithms, which verifies the validity of the potential crowding distance and improves the mutation operator to improve the optimization ability of the algorithm.

**Figure 5.** Pareto-optimal set of three algorithms.

In order to verify the performance of the improved algorithm more intuitively, the convergence metric [29] and the spacing metric [32] are used to compare the three optimization algorithms.

(1) C\_mertic

Convergence metric can be expressed as follows:

$$I\_C(X', X'') = \frac{|\{a' \in X''; \exists a' \in X' : a' \prec a''\}|}{|X''|} \tag{41}$$

where *X* and *X* are respectively two solution sets, *a* and *a* are solutions belonging to two solution sets, *p* is the dominant symbol, and |*X* | is the number of solutions in the solution set *X* . If *IC*(*X* , *X*) is 1, all solutions in the solution set *X* dominate the solution in the solution set *X* , and if *IC*(*X* , *X*) is equal to 0, all solutions in the solution set *X* are not dominated by the solution set *X* .

As shown in Table 3, about 32.39% of NSGA-II and 33.04% of MOPSO are dominated by INSGA-II. In contrast, only about 6.41% and 7.92% of INSGA-II are dominated by NSGA-II and MOPSO, respectively. In addition, the computational complexity O (mn3) is similar, so the running time of the three algorithms is almost the same, which proves that the improved crowding distance operator can effectively improve the quality of the solution set and can ensure operation efficiency.


**Table 3.** Convergence metric of each algorithm.

(2) S\_mertic

Spacing metric [27] can be expressed as follows:

$$I\_S = \sqrt{\frac{1}{N-1} \sum\_{i=1}^{N} \left(\overline{d} - d\_i\right)^2} \tag{42}$$

where *N* represents the number of solutions, *di* represents the shortest distance from the *i*th individual to the rest of the solution, and *d* represents the mean of all individuals *di*. A smaller spacing metric means that the solution distribution in the Pareto solution is more uniform. The zero value of the interval metric means that all solutions in the Pareto-optimal set are equally spaced.

As shown in Figure 6, in order to visually describe the uniformity of the distribution of the solution sets of different algorithms, the diversity indices obtained by each algorithm running 30 times independently are represented by box diagrams, in which each box represents the distribution of the measurement values of the distance between the solution sets of different algorithms. The upper quartile line at the top of each box diagram and the lower quartile line at the bottom represent the boundary values except outliers. If there is an abnormal value, use "+" to identify it. The median value is the red line in the rectangle box. Compared with the other two algorithms, INSGA-II has smaller median values (see red line in the middle) and minimum values (see quartile line below). Therefore, it is further explained that the distribution of INSGA-II solutions is more uniform and has better diversity.

**Figure 6.** Three algorithms' spacing metric.

#### (3) Analysis of the relationship between objective functions

As shown in Figure 7a, the energy-saving benefit has a nonlinear relationship with line loss, and the degree of line loss improvement will have an extreme point, which will not decrease with the increase of DG capacity but will be damaged. In view of Figure 7b, the voltage increases with the increase of energy-saving benefits. The results of DG optimization show that the capacity of DG-mounted nodes is too large to reduce the total voltage deviation correction effect. Therefore, the energy-saving benefits have a linear relationship with the voltage deviation when the permeability of DG is satisfied. Line loss and voltage deviation also show a nonlinear relationship in Figure 7c, and the two restrict each other.

(**a**) The energy-saving benefit relationship with line loss

(**b**) The energy-saving benefit relationship with voltage deviation

(**c**) The line loss relationship with voltage deviation

#### *4.4. Case Analyses of IEEE 33-, 69-, and 118-Bus Systems*

In order to prove the improvement degree of the NSGA-II algorithm on network optimization results of different node systems, the IEEE 33-, 69- and 118-bus systems are introduced. The parameters of the IEEE 69- and 118-bus systems are described in References [33,34].

Table 4 shows the optimization results of the INSGA-II algorithm on three examples. For the objective function, with the increase of network complexity, the improvement degree of the objective function of the three DNs increases. For example, IEEE 33-bus, 69-bus, and 118-bus are improved by 74.38%, 77.39%, and 77.32%, respectively, in terms of voltage deviation. On the hand, it also proves the importance of reasonable installation of DG. On the other hand, most of the installation nodes are the end of the system, which also proves that DGs can improve the voltage deviation better than other nodes.


**Table 4.** Results of different network optimization.

#### *4.5. Experiments on Load and DG Output Considering Annual Time Series Changes*

With the increase of installations and capacity of DG in the future, considering only static load and DG operation status, it may not be applicable to other load levels of permeability constraints and may even lead to reverse power flow, voltage limit, original system protection measure failures, and other phenomena, so that the system is threatened.

To this end, this section optimizes the configuration scheme to consider the DG annual output change and load change trend and the revised IEEE 33-bus system-related data reference [34]. The data was recorded every hour, so it is assumed that the DG output is equal every hour and that the load demand is constant. The four-season output level and load demand of wind turbines are shown in Figure 8a,b.

(**b**) The four-season load demand

WK

**Figure 8.** Time series data of wind turbine output and load.

As seen in Figure 8a,b, the wind speed is higher in spring and winter and the output of wind turbine is also improved. On the contrary, the output of wind turbine in summer and autumn is small and the seasonal data of the load has changed but the overall change is not large. On the one hand, if the optimization is based on the summer and autumn season time data, it may lead to waste of DG capacity installation. On the other hand, optimizing the configuration in spring and winter will destroy the network penetration rate and cause economic and security crisis. Here, we choose to optimize the four seasons separately and consider them comprehensively. The four seasonal DG optimization results using INSGA-II are shown in Table 5.


**Table 5.** Four quarterly independent optimization results.

As shown in Table 5, due to the relatively stable output of the wind turbines in summer and autumn, it is possible to seek a relatively high-quality configuration result. The spring and winter seasons are relatively unsatisfactory because of the large fluctuations in output. The combined installation capacity of the dynamic optimization scheme should meet the DG permeability constraint of the remaining time period. In each season, it is quite difficult for a DG configuration scheme to improve the optimization function while satisfying the different constraints of 24 h. Therefore, the installation capacity of each season is almost close to the maximum installation capacity, which is a compromise that is comprehensively considered to optimize the time-series data of each objective function. In order to meet the annual operating constraints, the combination of the smallest installation capacity on each quarter node is selected as the planning scheme.

### **5. Conclusions**

Considering the influence of DG access on the economy, security, and reliability of DN, this paper establishes a multi-objective optimization model to minimize the line losses and voltage deviation and to maximize the annual energy benefit. For DN with integrated DG units, equality constraints and inequality constraints, which include power-flow equality, nodal voltage, branch current, DG capacity, and short-circuit current, are considered in the optimization model. The principle of installation of DG is summarized through experiments and power system theory, and the violation constraint indicators are proposed to solve the infeasible solution.

The mutation operator, crowding distance operator, and selection operator of traditional NSGA II are improved to increase the population diversity and consistency in the optimal allocation. The introduction of an unbiased, compromised strategy can quickly give decision makers a comprehensive plan to weigh each goal.

The different examples of the IEEE 33-, 69-, and 118-bus systems are analyzed and optimized, and different super parameters are compared experimentally. The performance of the proposed INSGA-II is verified by comparing NSGA-II and MOPSO algorithms. The IEEE 69- and 118-bus systems are selected as the verification case. The proposed algorithm is applied to solve the static and dynamic characteristics of the DN with the integrated DGs. Finally, the IEEE 33-bus system is modified and optimized by using the four season data of wind turbine and load. The results indicate that the proposed method can achieve better precision and diversity and can provide a good configuration plan for decision makers under the premise of meeting the annual penetration rate.

In practice, the choice of the best site may not always be feasible due to many reality constraints. However, the optimization and analysis here suggest that considering multi-objectives helps to decide siting and sizing of DG units for the decision maker. The optimization planning of a multi-distributed generation access distribution system combined with geographic information system (GIS) technology will be further investigated.

**Author Contributions:** Conceptualization, W.L. and F.L.; methodology, W.L.; software, F.L.; validation, Y.L., F.L., and W.D.; formal analysis, W.L.; investigation, F.L.; resources, F.L.; data curation, W.D.; writing—original draft preparation, F.L.; writing—review and editing, W.L. and Y.L.; visualization, W.D.; supervision, W.L.; project administration, W.L.; funding acquisition, W.L.

**Funding:** This research was funded by Natural Science Foundation of Heilongjiang Province, China, grant number E201332.

**Acknowledgments:** The authors would like to thank the editors and the reviewers for their constructive comments and suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **Abbreviations**


