*4.1. Improved PSO Algorithm*

The particle swarm optimization (PSO) algorithm is an iterative optimization algorithm. By converting the target into a certain number of particles, the position is updated in each iteration, and the optimal solution is searched through continuous iteration. The updating method is as follows:

$$w\_{id}^{k+1} = wv\_{id}^k + c\_1 r\_1 \left(p\_{id}^{(k)} - \mathbf{x}\_{id}^{(k)}\right) + c\_2 r\_2 \left(\mathbf{g}\_d^{(k)} - \mathbf{x}\_{id}^{(k)}\right) \tag{13}$$

$$
\pi\_{id}^{(k+1)} = \pi\_{id}^{(k)} + v\_{id}^{(k+1)} \tag{14}
$$

where *w* represents inertia weight; *c*1, *c*<sup>2</sup> are acceleration coefficients of particle motion respectively. *<sup>r</sup>*1,*r*<sup>2</sup> are randomly selected in the range of 0 to 1; *<sup>P</sup>*(*k*) *id* represents the *d*dimensional component of the particle numbered *i* in the optimal position vector at time *k*; and *g* (*k*) *<sup>d</sup>* represents the *d*-dimensional component of the optimal position of all particles at time *k*.

In the process of iterative optimization, the traditional muti-objective particle swarm optimization algorithm is prone to fall into local optimal and appear 'Premature convergence'. Therefore, this paper cross-mutates the prescribed bit vectors of particles to prevent them from falling into the local optimum.

$$X\_{i,d} = X\_{min} + (X\_{max} - X\_{min})r\tag{15}$$

where *Xmin* represents the minimum position variable of the particle; *r* represents any value between 0 and 1.

Random cross variation is carried out when the probability distribution *P* < *Pm* is satisfied. When the fitness of particles tends to be the same or locally optimal, the *w* will increase. Otherwise, the *w* will decrease as the particles tend to disperse [36].

$$w = \begin{cases} w\_{\max} - \frac{(w\_{\max} - w\_{\min})\left(f\_i - f\_{\text{avg}}\right)}{f\_{\max} - f\_{\text{avg}}} & f\_i \ge f\_{\text{avg}}\\ & w\_{\max} \, f\_i \le f\_{\text{avg}} \end{cases} \tag{16}$$

where *fi* is the fitness value of the particle *i*, *wmax*,*wmin* are the maximum and minimum of *w*, and *fmax*, *fmin*, *favg* are the average and maximum, minimum, and average fitness of all particles at present, respectively.

In order to make the optimal solution distributed evenly in a certain range, the Pareto solution set should be optimized step by step. In this paper, a dynamic image Pareto solution set updating method is adopted, as is shown in Figure 6. In the initial distribution stage of particles, referring to all particles of the original comparison rule, an image *ViD* is established with the optimal particle as an optimal Pareto solution set. With the progress of iteration, the average value of particles is obtained between every two images, and the optimal solution is used to make the speed and position of particles updated at this time so that the particle population moves towards the target direction. In the whole process, new particles constantly exchange information with image particles and update data in image particles continuously until the end of iteration:

$$W\_{iD}^{t+1} = wv\_{iD}^{t+1} + c\_1 rand\_1 \left(p\_{iD} - X\_{iD}^{t+1}\right) + c\_2 rand\_2 \left(Average\left(pbest\_p\right) - x\_{iD}^t\right) \tag{17}$$

$$
\mathfrak{x}\_{iD}^{k+1} = \mathfrak{x}\_{iD}^{k} + \mathfrak{v}\_{iD}^{k+1} \tag{18}
$$

where *t* is the number of iterations; *D* is the dimension of the decision variable; *piD* is the best historical value of single particle *i*; and *pbestp* is the optimal particle of the current Pareto solution set in the image. *Average pbestp* is the superior particle of the optimal Pareto solution set in the two images.

**Figure 6.** Dynamic image Pareto solution updates schematic.

There are differences in the selection of the optimal solutions between the multiobjective particle swarm optimization and single-objective particle swarm optimization, and the results obtained are complementary dominated Pareto solutions, which cannot be obtained by direct comparison of the particle fitness function. Therefore, in this section, the first 20% Pareto solutions with lower crowding distances and higher priority orders are randomly selected to guide the iterative updating of the particle population.

#### *4.2. Multi-Attribute Decision-Making Based on TOPSIS Method*

After solving the lower multi-objective optimization problem, the optimal solution obtained by the IMOPSO algorithm is a set of Pareto solutions, and the selection of the optimal solution is essentially a multi-attribute decision problem. TOPSIS is a method for ordering by similarity to an ideal solution, it selects the optimal solution set and the worst solution set through the established initial decision matrix, then compares the distance between the two solution sets and the evaluation index with the optimal solution set and finally sorts them to evaluate the pros and cons of the scheme. The TOPSIS method has high strict requirements in selecting weights. In this paper, the information entropy method is used to determine the weight of each target value. The information entropy method determines *w* by the difference of the target value in the Pareto solution, improves the accuracy of the final decision, reduces the difference, and ensures the objectivity of the decision. By using the TOPSIS method, we can determine a set of optimal Pareto solutions to guide us to choose an energy storage configuration scheme.

The optimal solution of the Pareto solution set obtained is selected from *X*<sup>1</sup> ∼ *XN* and combined into *N* alternative schemes. The scheme *Xi* is selected from *N* records. It is the composition of some optimal solutions in the Pareto solution set. *gm*(*Xi*) represents the value of the *m*th attribute of the scheme *Xi*. Since each attribute is different, it should be unified and changed into the same type. The new attribute value is *Gm*(*Xi*) , which can be expressed as:

$$G\_m(X\_i) = \frac{g\_m(X\_i)}{\sqrt{\frac{1}{N} \sum\_{i=1}^{N} g\_m^2(X\_i)}}\tag{19}$$

$$d(\mathbf{x}\_i) = \frac{d\_+(\mathbf{x}\_i)}{d\_+(\mathbf{x}\_i) + d\_-(\mathbf{x}\_i)} \tag{20}$$

$$d\_{+}(\mathbf{x}\_{i}) = \sqrt{\sum\_{m=1}^{n} \left[\lambda\_{m}\mathbf{g}\_{m}^{\prime}(\mathbf{x}\_{i}) - \lambda\_{m}\mathbf{g}\_{m+}^{\prime}\right]^{2}}\tag{21}$$

$$d\_{-}(\mathbf{x}\_{i}) = \sqrt{\sum\_{m=1}^{n} \left[\lambda\_{m}\mathbf{g}\_{m}^{\prime}(\mathbf{x}\_{i}) - \lambda\_{m}\mathbf{g}\_{m-}^{\prime}\right]^{2}}\tag{22}$$

where *d*(*xi*) is the relative distance of scheme *xi*; *d*+(*xi*) represents the distance between scheme *xi* and the optimal solution. *d*−(*xi*) represents the distance from solution *xi*, the negative worst solution. *λ<sup>m</sup>* indicates the weight value of *gm*(*Xi*), which is randomly set between 0 and 1. *g <sup>m</sup>*<sup>+</sup> and *g <sup>m</sup>*<sup>−</sup> indicate the optimal and worst values of all schemes *gm*.

#### **5. Analysis and Discussion**

### *5.1. Case Description*

In this paper, the proposed scheme is tested on the IEEE-33 node distribution network [37]. In addition, the structure of the system is shown in Figure 7.

In this paper, the rated voltage of the selected distribution network is 12.66 kV, and the total load is 3715 kW + j2300 kvar. The upper and lower limits of the node voltage are specified as not exceeding ±5% of the rated voltage. Node 1 is a balance node, which is connected to the upper-level distribution network for power transmission. Taking into account the actual work and construction of photovoltaic power generation connected to the distribution network, photovoltaic power generation is connected to Node 9, and its installed power generation capacity is 1.077 MW (29% penetration rate).

**Figure 7.** The topology diagram of IEEE-33 bus system.

The typical PV output curve selected by the method in Section 2 is shown in Figure 8. The typical daily load curve in this area is shown in Figure 9.

**Figure 8.** Photovoltaic typical sunrise force curve.

**Figure 9.** Typical daily characteristic curve of load.

In this paper, the battery is used as the energy storage system for research and introduces the time-of-use pricing strategy proposed in [38]. The specific time-of-use price is shown in Table 2. The energy storage control parameters are shown in Table 3. The specific setting parameters of the energy storage configuration optimization simulation are shown in Table 4.


**Table 2.** Time-of-use electricity price table.

**Table 3.** Energy storage control parameter table.


**Table 4.** Simulation parameter settings.


In order to study the actual effect of energy storage configuration, we first analyzed the specific benefits of a photovoltaic distribution network connecting to energy storage configuration and demonstrated that energy storage still has good benefits in the high-light volt distribution network. Then, we compared the photovoltaic distribution network scenarios under different permeability and analyzed and compared the change of photovoltaic permeability with the corresponding change of optimal energy storage configuration scheme. The specific analysis content is introduced in the following section.

### *5.2. Energy Storage Optimization Scenario Division*

Analyze the effectiveness of the method proposed in this paper, set different conditions, divide it into four scenarios, and compare them one by one to verify the feasibility of the method:


storage system is charged and discharged at a constant power regardless of the high or low electricity price.

Scenario 4: The optimal configuration result of energy storage in Scenario 2 is used as the constraint condition of this scenario, and the traditional multi-objective PSO algorithm is used to simulate and analyze the lower model in the optimal configuration model of the energy storage double-level. Node voltage curves and load curves in different scenarios are shown in Figures 10 and 11 below, and Table 5 shows the optimization results of different scenarios.

**Figure 10.** System node voltage curve in different scenarios.



**Figure 11.** Load curve in different scenarios.

By comparing Scenario 1 and Scenario 2, it can be found that the voltage amplitude curve of the photovoltaic distribution network is smoother after the energy storage is connected, and the voltage fluctuation and load fluctuation are reduced to a large extent, which indicates that the BESS plays a good role in suppressing the node voltage fluctuation and load fluctuation when it is connected to the distribution network. Compared with Scenario 3, the load fluctuation range of Scenario 2 is smaller, and the load smoothing capacity is better. At the same time, the total cost of Scenario 2 is 803.27 yuan lower than that of Scenario 3 economically, which verifies the good characteristics of the model proposed in this paper.

The optimal Pareto solution set distribution of Scenario 2 and Scenario 4 is shown in Figure 12. Scenario 2 adopts the improved IMOPSO algorithm in this paper to solve the inner model, and Scenario 4 adopts the unimproved MOPSO algorithm to solve it. The Pareto solution set in Scene 2 is more evenly distributed than that in Scene 4 due to the introduction of particle cross mutation, adaptive inertia weight, and the Pareto solution set update method of the dynamic image. Moreover, Scene 2 adopts multi-attribute decisions based on the TOPSIS method, resulting in a more diverse solution set.

**Figure 12.** Pareto solution set distribution.

Through the comparison of scenes, it is obvious that the optimization results of the IMOPSO algorithm are obviously better than the MOPSO algorithm, and the search accuracy is higher. In order to compare the performance of the two algorithms, the external solution set and the spacing S are used in this paper to measure the optimization performance of the two algorithms. The S index refers to whether the particles in the Pareto solution set are evenly distributed in space. The mean variance of the particle density distance is used in this paper to characterize the uniformity and global nature of the population particles, as shown in Equation (23).

$$S = \sqrt{\frac{1}{N} \sum\_{i=1}^{N} [I(\mathbf{x}\_i) - \overline{I}]^2} \tag{23}$$

where *I* represents the average of all particles *I*(*xi*) in the Pareto solution set.

According to the different internal environments of the two algorithms, after 20 cycles, the node voltage and load fluctuations in the optimization target are taken as the research object, as shown in Figure 13 and Table 6.

**Figure 13.** Convergence curves of external solutions for different objectives.



By combining and comparing the charts, it was concluded that the IMOPSO algorithm proposed in this paper reduces the number of iterations in the node voltage fluctuation and load fluctuation, and the convergence performance is obviously better than the MOPSO algorithm. In addition, the improved algorithm and Pareto solution set update strategy make the solution set distribution more uniform and the type of solution set more diverse, and the improved MOPSO has better robustness and convergence.

#### *5.3. Energy Storage Benefit Analysis under Different Photovoltaic Permeability*

In order to verify the effectiveness of the dual-layer multi-objective optimal configuration model of the energy storage system proposed in this paper in the high-light volt permeability distribution network, the upper limit of photovoltaic power generation permeability was set at 60%, and the verification started from 30% permeability. Using the optimal configuration strategy of the BESS, the curve as shown in Figure 14 was obtained.

**Figure 14.** Changes in energy storage capacity and power under different photovoltaic penetration rates.

When the photovoltaic permeability increases from 30% to 60%, the capacity and power of the energy storage system have an obvious rising trend. When the photovoltaic permeability reaches 50%, the growth slows down and tends to remain unchanged. In other words, it is of little significance to increase the capacity of the energy storage system when the permeability reaches a certain level.

Figure 15 below shows the variation trend of energy storage investment and the total cost of distribution network operation under different photovoltaic permeability. It can be clearly seen that the total cost of the system decreases first and then increases when the photovoltaic permeability increases, and the total cost is the minimum when the permeability is 45%. As the cost of photovoltaic power generation decreases with the continuous increase of the permeability but is limited by the load level, the cost of the energy storage system increases with the increase of the capacity. The interaction between the two makes the total cost of the system decrease to the minimum when the photovoltaic permeability is 45%. When the permeability increases again, the system's total cost will keep rising, and the system operation economy will be seriously affected.

Based on the discussion of the above two legends, it is found that the total capacity of the BESS should be controlled in the optimal range according to the actual situation, and the photovoltaic permeability should also be controlled at a certain value so as to ensure the system operation economy while ensuring the safe and stable operation of the system. In order to improve the overall economy of the system, this paper selected 45% photovoltaic permeability to verify and analyze the two-layer programming model of the energy storage system proposed in this paper.

As shown in Figure 16 below, after optimizing the configuration of the energy storage system with 45% photovoltaic permeability, the load curve of the distribution network presents an obvious smoothing trend, and the peak–valley difference decreases. The sufficiency proves that the two-layer optimal configuration model of energy storage can still effectively improve the off-peak load, reduce the peak load of the distribution network, and increase the scheduling flexibility of the distribution network under the condition of high photovoltaic permeability.

**Figure 15.** Total cost curve under different PV penetration rates.

**Figure 16.** Distribution network load curve before and after energy storage configuration optimization under 45% photovoltaic penetration rate.

#### **6. Conclusions**

In order to ensure the power quality and scheduling flexibility of the photovoltaic distribution network with increasing permeability, this paper proposes a joint optimization operation mode of optical storage. Firstly, the PV output model was analyzed, and the scenario planning method was applied. The K-means clustering algorithm was used to divide the output scenarios, and the typical output scenarios were selected for analysis. A BESS two-layer decision model was established, and the improved IMOPSO algorithm was used to solve the two-layer model. The IEEE-33 node example was adopted, and the simulation verification was carried out based on the current feed-in price, selected energy storage parameters, and other parameters. The simulation analysis results are as follows:


#### **7. Future Work**

In this paper, only batteries are considered in the selection of batteries in the energy storage system. However, with a wider application of energy storage, a single energy storage system may not be able to meet the actual demand in the future. In subsequent research, we will combine other types of energy storage for optimization analysis of hybrid energy storage.

At the same time, because of the variety of renewable energy, more and more distributed power is connected to the distribution network. This paper only analyzes access to photovoltaic power generation. In a follow-up study, we will conduct a further study on scenarios with access to various energy sources.

In addition to the voltage fluctuation and load fluctuation considered in this paper, the power system with energy storage access has more indicators to measure security. In a follow-up study, we will also analyze the improvement and influence of energy storage access on various indicators.

**Author Contributions:** Conceptualization, J.Z. and L.Z.; methodology, S.Z.; software, S.Z. and L.Z.; validation, J.Z., S.Z. and L.Z.; formal analysis, J.Y.; investigation, L.L.; resources, J.Y.; data curation, J.Y.; writing—original draft preparation, L.Z.; writing—review and editing, J.Z.; visualization, L.L.; supervision, J.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is supported in part by the National Key Research and Development Program Project (Grant number: 2019YFE0104800), Natural Science Foundation of Henan Province (Grant number: 202300410271) and the Key Project of Science and Technology in Colleges and Universities of Henan Province (Grant number: 162102110130).

**Data Availability Statement:** The original contributions presented in the study are included in the article; further inquiries can be directed to the corresponding author.

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