In this work, the total investment operating cost of the BESSs, the total load fluctuation of ADN, the total voltage fluctuation of ADN, and the BESSs configuration optimization model based on Pareto can be established as follows: min ⁡ F x = min F 1 ,   F 2 , F 3 (1) where F x consists of objective functions F 1 ,   F 2 , F 3 , which are the total investment operating cost of the BESSs, total load fluctuation of ADN, and total voltage fluctuation of ADN; x represents the decision space consisting of the optimal installation node location, capacity, and hourly power optimization variables of BESSs and HESSs.

U i min ≤ U i ≤ U i max (13) where U i max and U i min represent the maximum and minimum of the node voltage respectively. E B E S S s min ≤ E B E S S s ≤ E B E S S s max (14) where, E B E S S s max and E B E S S s min represent the maximum and minimum of the installed capacity of the BESSs, respectively.

N B E S S s min ≤ N B E S S s ≤ N B E S S s max , ⋯ N B E S S , 1 ≠ N B E S S , 2 (15) where N B E S S s max and N B E S S s min represent the maximum and minimum of the location of the BESSs installation node, respectively.

R a = ∑ m = 1 M d D m ∙ ∑ t = 1 T P w i n d t + P P V t − P c / d t − P l o a d t − P l o s s t ∑ m = 1 M d D m ∙ ∑ t = 1 T P w i n d t + P P V t (16) R a ≤ R a ,   max (17) where R a represents the rate of wind and solar abandonment, which can be calculated by Eq. 16; R a ,   max represents the maximum rate of wind and solar abandonment.

NSGA-II has diverse populations and high-accuracy optimization results. Firstly, the algorithm uses a fast non-dominated sorting method to classify the population into different ranking levels. Secondly, the crowding distance between individuals with the same ranking level is calculated, and individuals with a larger crowding distance are selected. The results can be evenly distributed in the target space, thus maintaining the diversity of the population. Finally, the optimal solution in the Pareto solution set is effectively saved by the elite strategy (Wu et al., 2014).

(1) Non-dominated solution set ranking

The non-dominated solution in the population can be determined by calculating the value of the objective function. According to the non-dominated level of each individual in the population, the whole population is sorted hierarchically.

(2) Calculate the crowding distance

To further select individuals at the same non-dominated level, the crowding distance between individuals needs to be calculated. L d = L d + f n i + 1 − f n i − 1 f n max − f n min (18) where L d represents the crowding distance; f n i + 1 and f n i − 1 represent the nth fitness value of the i + 1 th and i − 1 th individual at the same level, respectively; f n max and f n min represent the maximum and minimum value of the nth fitness function.

Then through the binary tournament selection method, two individuals are randomly selected, and the individuals with higher non-dominated levels can enter the next population. If the level is the same, choose the more crowded population.

(3) Elite strategy selection method

The simulated binary method is used to cross and polynomially mutate the population. Finally, the parent population and offspring population are combined to form a new population. The above steps should be repeated until the maximum number of iterations is reached, and output the Pareto optimal solution set.

In the process of NSGA-II algorithm optimization, the Pareto solution set will be continuously updated and put into a limited scale storage pool. When the algorithm iteratively obtains a new non-dominated solution, it needs to compare it with the non-dominated solution set in the storage pool, and then determine whether the new non-dominated solution is liberated into the storage pool. The judgment process is as follows (Faramarzi et al., 2020):

(1) If the new solution dominates the storage pool, it is replaced with the new solution;

To increase the diversity of solutions, when the number of non-dominated solutions in the archive set is higher than the upper limit, the redundant solutions must be removed. In this work, we use the crowded distance ranking method to select a densely distributed set of non-dominated solutions and eliminate redundant solutions by the roulette wheel method.

A multi-objective gray target decision based on the entropy weight method (EWM) is designed. Firstly, the evaluation indexes of each solution are set to build the sample matrix, and the decision matrix can be built by normalizing the sample matrix. Then, the weights of evaluation indexes of all Pareto non-dominated solutions and their distances from the target are obtained based on EWM and their distances from the target, and the solution with the closest distance to the target is determined as the best compromise solution for the BESS optimal allocation scheme (He et al., 2021).

(1) Establishing the sample matrix

To evaluate the similarity and equilibrium of each solution, this method adds two indicators in the sample matrix. One is the Euclidean distance (ED) between each solution and the ideal point, and the other is the Mahalanobis distance (MD) between each solution and the equilibrium point. Therefore, the sample matrix of effects containing five evaluation metrics can be expressed as: Z = v g k n × m + 2 = F , E , M = F 1 1 ⋯ F m 1 E 1 M 1 F 1 2 ⋯ F m 2 E 2 M 2 ⋯ F 1 n ⋱ ⋯ ⋯ F m n ⋯ E n ⋯ M n (19) where E and M denote the ideal point distance matrix and equilibrium point distance matrix of all solutions, respectively. The ED between the k th solution and the ideal point can be calculated by Eq. 20, and the MD between the k th solution and the equilibrium point can be calculated by Eq. 21. E k = ∑ g + 1 m + 2 F g k − P g 2 (20) M k = ∑ g = 1 m F g k − u g T ∑ − 1 F g k − u g (21) where P g denotes the ideal point of the g th objective; u g represents the average of all solutions under the g th objective; ∑ − 1 is the covariance matrix.

(2) Design of decision matrix

According to the principle of rewarding good and punishing bad, the following operator is designed: q g = 1 n ∑ g = 1 n X g k (22) where X g k represents the evaluation index of the k th solution at the g th objective.

Then the decision matrix is established, as follows: V = v g k n × m + 2 = q g − F g k max max 1 ≤ k ≤ n F g k − q g , q g − min 1 ≤ k ≤ n F g k (23)

The target is identified in the gray decision region formed by the decision matrix as follows: v g 0 = max ⁡ { v g k 1 ≤ k ≤ n } (24)

(3) Grey target decision based on EWM

First, the target weight of the k th scenario under the g th objective can be calculated by Eq. 25 w g k = F g k ∑ k = 1 n Z g k ,   i f   F g k ≥ 0 (25)

Next, the entropy values of the g targets can be calculated as follows: E V g = 1 ln ⁡ n ∑ k = 1 n w g k ∙ ln ⁡ w g k ,   i f   E V g > 0 (26)

Among them, the smaller the entropy value of the indicator tends to contain more information, so the entropy weight is larger. The entropy weight of the g th target is calculated as follows. w g = 1 − E V g ∑ g = 1 m + 2 1 − E V g (27)

Then, the distance of each solution to the target can be calculated as follows: M D B k = v k − v 0 = ∑ g = 1 m + 2 w g F g k − u g T ∑ − 1 F g k − u g (28)

Finally, the target distance of each solution is sorted, and the solution closest to the target is selected as the best decision. The specific process of Table 1 NSGA-II solving the optimization problem of BESSs is shown in Figure 1. The basic principles and procedures of NSGA-II are as follows:

The extended IEEE-33 node test system topology is shown in Figure 3. Wind power No. 1, No. 2, and No. 3 are connected at nodes 8, 24, and 17, respectively, and PV is connected at node 31. Table 1 shows the main parameters of ADN.

Table 3 shows the specific configuration scheme of BESSs. Table 4 shows the numerical changes of specific indicators in ADN before and after the access of the BESSs configuration scheme obtained by NSGA-II. By comparing the data in Table 4, after accessing BESSs, specifically, daily voltage fluctuation decreased by 4.8484 p. u. (the improvement rate is 38.54%), and the load fluctuation decreased by 2.4941 MW (the improvement rate is 37.77%). The total cost of BESSs is $-1.92e+03, which means a profit of $ 1.92e+03 a year after BESSs configuration. Moreover, Figure 4 compares the changes of three technical indicators of the ADN before and after the access of BESSs. As shown in Figure 4A, after BESSs is connected to the active distribution network, the node voltage level generally has an upward trend (at 1p.u. nearby). Meanwhile, it can be seen from Figure 4B that the average voltage fluctuation of nodes in the ADN decreases in 1 day. This means that BESSs can effectively improve the voltage quality of ADN. At last Figure 5 shows the charge and discharge power of BESSs in a day.

The extended IEEE-69 node test system topology is shown in Figure 6. Wind power No. 1, No. 2, and No. 3 are connected at nodes 24, 34, and 44, respectively, and PV is connected at node 64.

Table 5 shows the specific configuration scheme of BESSs. By comparing the data in Table 6, after accessing BESSs, specifically, the daily voltage fluctuation decreased by 8.095 p. u. (the improvement rate is 22.058%), and the load fluctuation decreased by 2.4481 MW (the improvement rate is 37.07%). Moreover, Figure 7 compares the changes in three technical indicators of the ADN before and after the access of BESSs. The total cost of BESSs is $ -1.2e+03. Figure 7 shows the three technical indicators before and after the access of BESSs. In Figure 7A, the node voltage level of node 57-node 65 has been greatly improved, and the curve in Figure 7B also shows that the voltage fluctuation of ADN has been effectively reduced after BESSs access. It can be concluded that by installing BESSs in the ADN, the power quality in the ADN has been improved, and the revenue can be $1.2e+03. And Figure 8 shows the charge and discharge power in a day under the extended IEEE-69 node system.

In this section, we use multi-objective particle swarm optimization (MOPSO) as a comparison algorithm to verify the superiority of NSGA-II. Table 7 shows the result obtained by two algorithms.

By comparing the results in Table 7, it can be seen that NSGA-II has the smallest daily voltage fluctuation and daily load fluctuation except for the total cost of BESSs. Meanwhile, the cost results obtained by the two algorithms are not much different. Therefore, it can be concluded that NSGA-II can obtain a more reasonable Pareto non-dominated solution set and provide the best candidate for decision makers.

Figure 9 shows the Pareto front distribution obtained by NSGA-II in the simulation experiments based on the extended IEEE-33 node test system and the extended IEEE-69 node test system. It can be seen from Figure 9 that the NSGA-II can obtain a widely distributed Pareto front. This proves the effectiveness and powerful optimizing ability of NSGA-II algorithm.

To verify the rationality of the EWM-based grey target decision, this work uses the method in Reference (Hu et al., 2014) to map the a th objective function corresponding to the b th non-dominated solution to the two-dimensional plane of ( a × i ). The equation is as follows: L i , a = K f b , a − f a min f a max − f a min (32) where K is a random parameter; f b , a represents the a th objective function of the b th non-dominated solution.

Figure 10 shows a two-dimensional plane ( a × i ) mapped the Pareto front of the extend IEEE-33 node system and the extended IEEE-69 node system through the high-dimensional Cartesian coordinates. The abscissa is the optimized objective function, and the ordinate is the mapped objective function value. The objective functions are connected by a virtual line. The optimal compromise solution of the grey target decision based on EWM is represented by the red solid line. It can be seen that the compromise solution is in a reasonable range, which can avoid a greater preference for a certain goal due to subjective weight. And the distribution of the non-dominated solution set is more extensive, which indicates that NSGA-II algorithm has excellent optimization ability.

It should be noted that at 1 p.m. of Figure 4C and 1 and 2 p.m. of Figure 7C, the load level of the ADN with BESS drops to a very low level. This is because when configuring BESSs, the cost of configuring BESSs is considered an important factor. At 1 and 2 p.m. every day, it is in the peak load period. At this time, the load of the ADN increases, and the BESSs will deliver more electricity to the ADN to obtain more benefits while ensuring the stability of the ADN. The net load P L of the ADN can be calculated by Eq. 33. As shown in Eq. 33, the increase of electric energy transmitted by BESSs to ADN will lead to a decrease in the net load level of the ADN. P L = P l o a d − P w i n d − P P V + P B E S S s (33) where P l o a d , P w i n d and P P V represent the load demand, wind poaxwer output, and PV output; P B E S S s represents the power of BESSs, and its value is negative when it outputs electricity to the ADN.