Photovoltaic penetration rate is defined as the ratio of the maximum photovoltaic output power to the maximum load output power: D P V = F MAX F L , MAX × 100 % (1)

In the formula 1: D P V represents the photovoltaic penetration rate; F MAX represents the maximum photovoltaic output power; F L , MAX represents the maximum load output power.

People have different criteria for judging the level of photovoltaic penetration. Generally, when it is below 20%, it is considered a low-penetration stage, where the scale of photovoltaic integration into the grid is small, and its impact on the power grid is minimal. As the penetration of photovoltaics increases, when the photovoltaic penetration rate reaches 20%–80%, it becomes necessary to address the issue of enhancing the carrying capacity of photovoltaics in concentrated areas. At this stage, the significant integration of photovoltaics may lead to power reverse flow in the power system. When the photovoltaic penetration rate exceeds 100%, it can be considered a high-penetration stage. At this point, there is a substantial reverse flow of photovoltaic power into the power system, and the role of photovoltaic energy shifts to the supply side. The difficulty of on-site absorption of photovoltaics increases, requiring interventions such as energy storage to enhance the capacity for photovoltaic integration.

The on-site absorption of photovoltaic power is primarily influenced by the load and energy storage. In this paper, we define the on-site absorption rate of photovoltaic power as: φ = ∑ t = 6 19 P P V , f a c t , t − P L O A D , t − P E S S , t P P V , f a c t , t × 100 % (2)

In the formula 2: Where φ represents the on-site absorption rate of photovoltaic power; P P V , f a c t , t represents the actual generation of photovoltaics at time t. P L O A D , T represents the distribution network load during the photovoltaic generation period at time t, and P E S S , t represents the load for energy storage charging during the photovoltaic generation period at time t.

Monte Carlo sampling is used to analyze the electric vehicle charging model.

(1)Probability Distribution of Daily Driving Distance

Processing behavioral data on electric vehicle usage, it is determined that the daily driving distance follows a log-normal distribution, i.e.,: f s x = 1 x 1 2 π exp − ln ⁡ x − μ s 2 2 σ s 2 (11)

In the formula 11: Where μ s taking the mean as 3.2 and σ s the variance as 0.88.

(2)Probability Distribution of Initial Charging Time

The charging time of electric vehicles f t x follows a normal distribution, i.e.,: f t x = 1 σ s 1 2 π exp − x − μ s 2 2 σ s 2 , μ s − 12 < x < 24 1 σ s 1 2 π exp − x + 24 − μ s 2 2 σ s 2 , 0 < x < μ s − 12 (12)

In the formula 12: Where μ s taking the mean as 17.6 and σ s the variance as 3.4.

The upper layer involves multi-energy storage optimization configuration, with the objective function being the minimization of equipment investment costs, equipment operating costs, and grid power purchase costs. f u p = w 1 × f 1 + w 2 × f 2 + w 3 × f 3 + w 4 × f 4 (13)

In the formula 13: Where f 1 represents equipment investment costs; f 2 represents equipment operating costs; f 3 is the grid power purchase cost; f 4 represents the cost of mobile energy storage migration; w 1 、 w 2 、 w 3 、 w 4 is a random number between 0 and 1, and w 1 + w 2 + w 3 + w 4 = 1 .

(1) Minimize equipment investment costs

The equipment investment cost includes one-time investment costs for both energy storage devices and photovoltaic equipment. f 1 = Q ∑ z = 1 Z N z C z Q = q 1 + q y 1 + q y − 1 (14)

In the formula 14: Q represents the capital recovery factor; q represents the annual interest rate; Z represents equipment type; C z represents the investment cost of equipment type Z; N z represents the number of equipment type Z.

(2) Minimize equipment operating costs

The operating cost of equipment refers to the costs associated with regular maintenance and repair of equipment damage. f 2 = ∑ t = 1 T ∑ z = 1 Z C z l o s s t + A (15)

In the formula 15: C z l o s s represents the daily operational maintenance cost of equipment z; A represents the cost of repairing equipment damage, where t is the operating time.

(3) Minimize grid power purchase costs

In the formula 16: P t l i n e represents the power purchased from the grid at time t, and w t l i n e represents the electricity price at time t.

(4) Minimize the cost of relocating mobile energy storage

In the formula 17: Where f 4 represents the cost of relocating mobile energy storage; C F U E L represents the unit distance cost, and O D I S T A N C E represents the distance traveled by the energy storage vehicle.

(1) Energy Storage Output Constraint.

In the formula 18: Where P E S , j min represents the lower limit of the charging or discharging power for the jth energy storage unit, P E S , j t represents the charging or discharging power of the jth energy storage unit at time t, P E S , j max represents the upper limit of the charging or discharging power for the jth energy storage unit, and S O C j t represents the state of charge of the energy storage at time t, with a range from 0.2 to 0.9. When t is zero, the initial state of charge S O C j 0 is set to 0.5.

(2) Electricity Purchasing Power Constraint

In the formula 19: Where P t , b u y represents the electricity purchasing power.

The lower level involves 24-h economic dispatch of the distribution network, with the minimization of voltage deviation as the objective function.

In power systems, voltage difference reflects the operational safety of the distribution network, as excessive voltage deviation can impact the operation of electrical equipment and the quality of electrical energy. f l o w = ∑ t = 1 24 ∑ i = 1 N Δ U i , t U i , ⁡ max − U i , ⁡ min 2 (20) Δ U i , t = U i , ⁡ min − U i , t , U i , t < U i , ⁡ min 0 , U i , ⁡ min ≤ U i , t ≤ U i , ⁡ max U i , t − U i , ⁡ max , U i , t ≥ U i , ⁡ max (21)

In the formulas 20, 21: Where Δ U i , t represents the voltage deviation at time t, U i , t represents the voltage at node i at time t, and U i , ⁡ min , U i , ⁡ max represent the upper and lower limits of node i, respectively. In medium and low voltage distribution networks, the permissible range for voltage deviation is −5% to +5%.

(1) Voltage magnitude constraints

In the formula 22: Where U i , ⁡ min represents the voltage lower limit at node i, U i represents the voltage at node i, and U i , ⁡ max represents the voltage upper limit at node i.

(2) Distributed New Energy Output Constraint

In the formula 23: Where P p v , ⁡ min represents the minimum output of photovoltaic power, P p v t represents the photovoltaic power output at time t, and P p v , ⁡ max represents the maximum output of photovoltaic or wind power.

(3) Power Balance Constraint

In the equation 24: P s , i represents the active power output at node i due to the power source; P L , i represents the active power output of the load at node i; U i represents the voltage at node i; U j represents the voltage at node j; Q s , i represents the reactive power output at node i due to the power source; Q L , i represents the reactive power output of the load at node i.

The particle swarm algorithm (Anantathanavit and Munlin, 2013) models each particle considering the current velocity, current position, and a distance-modifying function to pbest and gbest as follows. v i t + 1 = w v i t + c 1 × r a n d × p b e s t i t − x i t + c 2 × r a n d × g b e s t − x i t (25)

In the equation 25: w represents the weighted function; v i t represents the velocity of the i th particle at generation t; c 1 represents the weighting factor; r a n d is a random number between 0 and 1; p b e s t i t represents the best position of the particle at generation t; x i t represents the position of the particle at generation t; c 2 represents the weighting factor; g b e s t represents the best solution. w v i t represents the exploration ability of the particle; c 1 × r a n d × p b e s t i t − x i t represents the personal ability of the particle; c 2 × r a n d × g b e s t − x i t represents the cooperative ability of the particle x i t + 1 = x i t + v i t + 1 (26)

In the formula 26: Where x i t + 1 represents the position of particle i at generation t+1.

GSA (Doraghinejad et al., 2012) originates from Newton’s fundamental theory: the interaction force among particles in the Universe, a force proportional to the particle mass and inversely proportional to the distance between them, is modeled as follows. F i j d t = G t M p i t × M a j t R i j t + ξ x j d t − x i d t (27)

In the formula 27: Where F i j d t represents the gravitational force of particle i on particle j in the d-dimensional space at the tth iteration; G t represents the value of universal gravitational force at the tth iteration; M p i t is the active gravitational mass; M a j t is the passive gravitational mass; R i j t represents the Euclidean distance between i and j; ξ is a constant; x j d t represents the position of particle j in the d-dimensional space at the tth generation; x i d t represents the position of particle i in the d-dimensional space at the tth generation. G t = G 0 × e − ∂ × i t e r max ⁡ i t e r (28)

In the formula 28: Where ∂ represents the descent coefficient; G 0 represents the initial value; i t e r is the current iteration number; max ⁡ i t e r is the maximum iteration number. F i d t = ∑ j ≠ i , j ∈ k b e s t N ran d j F i j d t (29)

In the formula 29: Where F i d t represents the total force experienced by individual i in the d-dimensional space at the tth iteration. a i d t = F i d t M i t (30)

In the formula 30: Where represents a i d t the equation for the acceleration of individual i in the d-dimensional space; M i t represents the mass of individual i at generation t. v i d t + 1 = r a n d × v i d t + a i d t (31)

In the formula 31: Where v i d t + 1 represents the velocity of particle i in the d-dimensional space at generation t+1. x i d t + 1 = x i d t + v i d t + 1 (32)

In the formula 32: Where x i d t + 1 represents the position of particle i in the d-dimensional space at generation t+1.

The PSO-GSA hybrid algorithm combines the individual optimization capability of PSO with the local search ability of GSA such as the formulas 33, 34. The improved convergence of PSO-GSA surpasses that of standalone PSO and GSA. v i t + 1 = w v i t + c 1 × r a n d × a i t + c 2 × r a n d × g b e s t − x i t (33) x i t + 1 = x i t + v i t + 1 (34)

Initially, each particle is considered to have a candidate solution. After initialization, the gravitational force, gravitational constant, and resultant force between particles are calculated. During the iteration process, the algorithm updates to the current best solution, computes the velocity of particles for the (n+1)-th generation, and finally updates the positions of the particles.

The improved convergence of PSO-GSA is superior to that of PSO and GSA, as shown in Figure 2.

The flow chart of the PSO-GSA algorithm is shown in Figure 3.

Simulations were conducted on the IEEE 33-node distribution network using Matlab 2021a software. The system’s base voltage is 12.66 kV, and the maximum load is 3.715 MW. To provide reserve capacity for photovoltaic integration at system nodes, the upper limit of node voltage is set to 1.05, and the lower limit is set to 0.95. Photovoltaic panels are integrated at nodes 9 and 28, while an electric vehicle charging station is added at node 20, as illustrated in Figure 4.

Monte Carlo simulations were employed to model the starting density and charging power of electric vehicles, with a total of 2000 vehicles. The simulation results are presented in Appendix Figure A1 and Figure A2.

This study focuses primarily on the impact of the fixed and mobile energy storage access points and capacities on the integration of photovoltaics. The basic information for both fixed and mobile energy storage is as follows: the energy storage maintenance coefficient is 0.02; the unit capacity investment cost is 1000 CNY/KW; the discount rate is 0.08; the service life is 20 years. For mobile energy storage, the cost per kilometer varies based on the distance traveled each time, and here it is calculated at a monthly cost of 3,000 Yuan. The energy storage electricity prices are 0.31 CNY/kWh from 0:00 to 8:00, 0.84 CNY/kWh from 9:00 to 11:00, 0.31 CNY/kWh from 12:00 to 13:00, 0.84 CNY/kWh from 14:00 to 21:00, and 0.31 CNY/kWh from 22:00 to 24:00.

To achieve coordinated optimization of fixed and mobile energy storage for enhancing the distribution network’s consumption capacity, a PSO-GSA hybrid algorithm is applied to both the upper-layer multi-energy storage optimization configuration and the lower-layer energy storage optimization scheduling. The fixed energy storage locations range from node 2 to 33, with capacities from 0.5 MW to 1 MW. The access nodes for mobile energy storage range from node 2 to 33 (assuming node 1 is the reference node), with capacities from 0.4 MW to 0.9 MW. Fixed energy storage charges during off-peak hours or when photovoltaic energy cannot be accommodated and discharges during peak electricity demand. In contrast, mobile energy storage offers more flexible charge and discharge regulation, responding dynamically to real-time situations in case of emergencies or when fixed energy storage cannot effectively regulate. The PSO-GSA hybrid algorithm is applied to both upper and lower layers with 50 particles each, 50 iterations, an individual learning factor of 0.5, and a global learning factor of 1.5.

Firstly, without the addition of any energy storage, gradually increasing photovoltaics until the voltage exceeds the limit at nodes 9 or 28, marks the maximum photovoltaic capacity that the distribution network can bear. The maximum capacity is determined to be 2.9 MW.

Considering the future large-scale integration of photovoltaics and the transition of photovoltaic energy from the demand side to the supply side, there may be reverse power flows. In such scenarios, energy storage can be flexibly adjusted to enhance photovoltaic energy integration, reduce the risk of voltage exceeding limits, and improve the stability of the power system. When there is a sudden increase in photovoltaics and fixed energy storage devices cannot regulate effectively, flexible adjustments can be made using mobile energy storage. The following case considers an extreme photovoltaic output scenario of 2.9 MW and a charging station output of 0.6 MW.

To validate the effectiveness of the proposed model and method, a comparison is made across four different scenarios. Scenario One: integration of photovoltaics without energy storage; Scenario Two: integration of photovoltaics with optimized configuration of fixed energy storage; Scenario Three: integration of photovoltaics with coordinated optimization of fixed and mobile energy storage; Scenario Four: integration of photovoltaics, electric vehicle charging station, and coordinated optimization of fixed and mobile energy storage. The analysis includes voltage offset, multi-energy storage operating costs, and on-site photovoltaic integration rate, as shown in Table 4-1.

According to Table 1, compared to Scenario One, Scenario Two, which adds fixed energy storage, reduces the voltage offset by 0.0010 and increases the on-site photovoltaic integration rate by 3.01%. Scenario Three, with the addition of both fixed and mobile energy storage, reduces the voltage offset by 0.0018 and increases the on-site photovoltaic integration rate by 5.77% compared to Scenario One. In Scenario Four, with the addition of an electric vehicle charging station load, the voltage offset is 0.0033, the on-site photovoltaic integration rate is 71.39%. Compared to Scenario One, the voltage offset increases by 0.1638, and the curtailment rate decreases by 0.8081%. These results indicate that fixed energy storage adjustment has limitations, and through coordinated optimization of fixed and mobile energy storage, the on-site photovoltaic integration can be increased, and voltage offset can be reduced.

To further illustrate the improvement in power grid stability through the coordinated optimization of fixed and mobile energy storage, a comparative analysis is conducted among Scenario One, Scenario Two, and Scenario Three.

As shown in Figure 5, Figure 6, and Figure 7, in extreme photovoltaic Scenario One, voltage exceeds the limit at nodes 9 and 28 when photovoltaics are added. From the node voltage diagram in Scenario Two, it is evident that by adding fixed energy storage, only node 9 experiences a voltage limit exceedance, and the degree of voltage offset is smaller compared to Scenario One. The node voltage diagram in Scenario Three indicates that through the coordinated action of fixed and mobile energy storage, all nodes are within the range of 0.95 p. u. to 1.05 p. u., demonstrating the efficient synergy between fixed and mobile energy storage. This synergy can significantly enhance the capacity for photovoltaic integration.

As shown in Figure 8, for the charge and discharge strategy of fixed energy storage, during 3:00–7:00 when the grid load is relatively low, the energy storage system remains in the charging state. During 10:00–14:00 when the load is high and there is sufficient photovoltaic output, fixed energy storage discharges at a lower power. From 19:00 to 22:00, when the load is high, fixed energy storage discharges, and when fixed energy storage cannot meet the load requirements, coordinated operation with mobile energy storage is employed to jointly provide power support to the grid.

As illustrated in Figure 9, due to the uncertainty of photovoltaic output, there are two charging methods for the charge and discharge strategy of mobile energy storage: one is during 3:00–7:00 when the electricity price is lower, mobile energy storage utilizes grid electricity for charging; the other is during 14:00–16:00 when the load is low and photovoltaics cannot fully integrate, mobile energy storage is charged at the access node. During peak electricity demand periods at 10:00–14:00 and 19:00–22:00, if fixed energy storage cannot effectively regulate the grid voltage, coordinated discharge of mobile and fixed energy storage is implemented to maintain the stable operation of the power system.