The bus voltage margin refers to the buffer between the actual operating voltage of the bus and its upper voltage limit. This metric is used to assess the bus’s ability to remain within limits during grid disturbances and can be expressed as Equation 1. F 11 = U max − U U max − U min × 100 % (1) where U is the operating voltage of the evaluated bus and U _max and U _min are the maximum and minimum voltage limits, respectively.

The comprehensive voltage qualification rate refers to the proportion of time within an operational period during which the user side voltage remains within the acceptable range. This metric is used to measure the availability of power quality as experienced by the users and can be expressed as Equation 2. F 12 = ∑ i = 1 M t s t , i M T × 100 % (2) where M is the number of users; t _st,i is the number of hours that the voltage for user i remains within the acceptable range during an operational period; and T is the total assessment period.

The total harmonic distortion (THD) of the current refers to the ratio of the root mean square (RMS) value of the harmonic content to the RMS value of the fundamental component in the current flowing through electrical loads. When the total harmonic current distortion in the distribution network exceeds acceptable limits, it can cause excessive heating in transformers and motors, leading to damage to electrical equipment (Zuo et al., 2015). Therefore, this paper considers the harmonic current margin of the distribution network following the integration of distributed PV and ESS, which can be expressed as follows: T H D = ∑ h = 2 ∞ I h 2 I 1 − I P V − I E S S × 100 % F 13 = 0 , δ U > 5 % 5 % − T H D 5 % − 3 % , 3 % < δ U ≤ 5 % 1 , δ U ≤ 3 % (3) where I _h is the hth harmonic current; I ₁ is the fundamental current; I _PV is the fundamental current injected by the PV generation system; I _ESS is the fundamental current injected by the ESS; and 5% and 3% represent the limit and warning values, respectively (IEEE Recommended Practices and Requirements, 1992).

Voltage deviation refers to the relative difference between the actual operating voltage at a node and the system’s nominal voltage. When the voltage deviation exceeds acceptable limits, it can reduce the performance or even cause damage to electrical equipment, such as motors (Zuo et al., 2015). Therefore, this paper considers the voltage deviation margin of the distribution network, which can be expressed as follows: δ U = U i − U N U N × 100 % F 14 = 0 , δ U > 7 % 7 % − δ U 7 % − 4.2 % , 4.2 % < δ U ≤ 7 % 1 , δ U ≤ 4.2 % (4) where U _i is the actual operating voltage at node i; U _N is the nominal system voltage; and 7% and 4.2% represent the limit and warning values, respectively (General Administration of Quality Supervision, 2008b).

Voltage fluctuation refers to the ratio of the difference between the maximum and minimum voltages at a node during an operational period to the system’s nominal volt-age. When the voltage fluctuation exceeds acceptable limits, it can cause malfunction or even damage to electronic instruments and equipment (Zuo et al., 2015). Therefore, this paper considers the voltage fluctuation margin of the distribution network, which is expressed as follows: d = Δ U i U N × 100 % F 15 = 0 , δ U > 3 % 3 % − d 3 % − 1.8 % , 1.8 % < δ U ≤ 3 % 1 , δ U ≤ 1.8 % (5) where ΔU _i is the difference between the maximum and minimum voltages at node i during the operational period; 3% and 1.8% represent the limit and warning values, respectively (General Administration of Quality Supervision, 2008a).

Power supply reliability refers to the ability of a power system to provide a continuous electricity supply. The increasing prevalence of sensitive equipment has heightened the impact of voltage sags, which can disrupt the normal operation of user equipment and even cause damage, leading to significant losses and hazards for users. Therefore, this study modifies the power supply reliability to account for the impact of voltage sag events, which is expressed as follows: F 21 = 1 − t t d + t t d − s a g T × 100 % (6) where t _td represents the average number of outage hours for users in the region and t _td-sag is the equivalent average number of outage hours due to voltage sag events for users in the region.

The load loss rate refers to the proportion of the interrupted load to the total load during extreme events. This metric is used to assess the ability of the distribution network to maintain a power supply to critical loads during such events and is expressed as Equation 7. F 22 = P z d / P l o a d × 100 % (7) where P _zd represents all the interrupted loads in the system and P _load represents all the loads in the system.

The line average load rate refers to the average ratio of the maximum load to the capacity of each line in the distribution network. This metric is used to assess the ability of the distribution network to handle future load uncertainties and fluctuations and can be expressed as Equation 8: F 31 = R t o t a l / N b (8) where R _total is the sum of the load rates of all lines in the network; N _b is the number of lines in the system.

The distributed generation penetration rate refers to the proportion of the total capacity of flexible resources connected to the system. This metric is used to assess the ability of the distribution network to regulate various flexible resources in response to disturbances and can be expressed as Equation 9. F 32 = P W + P P V P G + P W + P P V × 100 % (9) where P _W and P _PV are the total capacities of the wind and photovoltaic systems connected to the grid, respectively, and P _G is the total capacity of the thermal power generation units.

The load peak-to-valley difference rate refers to the ratio of the difference between the maximum and minimum loads to the maximum load within a specific time period in the power system. This metric is used to assess the power regulation capacity required by the distribution network to accommodate daily load variations, which can be expressed as Equation 10. F 33 = P H − P L P H × 100 % (10) where P _H is the peak load; P _L is the valley load in the region on a typical day.

It should be noted that the design of the resilience assessment framework in this study, particularly in selecting indicators for resistance and recovery, fully considers the typical disturbances associated with high distributed PV penetration, such as voltage deviations and harmonic currents. These factors are especially prominent given the fluctuations in PV output, which significantly impact the power quality and stability of distribution networks. Consequently, the adjustments and enhancements made to resilience indicators in this context allow the model to more accurately reflect the resilience characteristics of systems under high PV penetration conditions.

Considering the varying operating environments and service users of different types of distributed power sources, it is crucial to consider the perspectives of users from different regions when conducting comprehensive resilience assessments. Therefore, this article adopts the improved AHP. The specific implementation steps are as follows.

Step 1: Determine the hierarchical structure of the comprehensive resilience indicators based on suggestions from all stakeholders, and rank them according to their importance.

The specific steps of entropy weight method are as follows.

Step 1: Indices dimensionless. Since the indicators of resistance, resilience, and adaptability in this article are all positive indicators (i.e., the larger the indicator value, the better), the calculated formula is as Equation 13.

Based on the abovementioned index system, the comprehensive resilience index measurement vector F is obtained, the comprehensive assessment result F _re is expressed as Equations 16, 17: q j = s i ω j ∑ j = j = 1 n s i ω j (16) F = F 11 , F 12 , … , F 33 Q = q 1 , q 2 , … , q n F r e = F × Q t (17) where n represents the number of indices to be evaluated.

The total cost of the model mainly includes the investment conversion cost C _invest , operation and maintenance cost C _maint , and network loss cost C _loss of the distribution network ES device within one operating cycle. This paper further considers the economic loss cost C _sag caused by user voltage sag events. The objective function is as follows: min ⁡ C t o t a l C t o t a l = C i n v e s t + C m a int + C l o s s + C s a g C i n v e s t = T 8760 ∑ i ∈ Ω e s s r e s s c i n v e s t E e s s , i 1 − 1 + r e s s − y e s s C m a int = T 8760 ∑ i ∈ Ω e s s r e s s c m a int E e s s , i 1 − 1 + r e s s − y e s s C l o s s = ∑ i , j ∈ Ω b ∑ t = 1 T c l o s s I i j , t 2 r i j C s a g = ∑ k = 1 K c s a g (22) where Ω _ess is the set of all nodes configured with ESS; R _ess is the discount rate of ESS; Y _es is the service life of ESS; C _invest is the investment cost per unit capacity; E _ess,i is the ESS capacity configured for node i; C _maint is the annual operating and maintenance cost per unit capacity; Ω _b is the set of all branches in the distribution network; C _loss is the unit network loss cost; I _ij,t is the current of branch ij at time t; R _ij is the resistance of branch ij; C _SAG is the average economic loss value of a single voltage sag for the user, and K is the number of voltage sag events that occur within one operating cycle.

Based on the conventional grid elasticity index evaluation system, this paper further considers the impact of distributed PV access on power quality indicators, and proposes a comprehensive resilience index that takes into account both traditional elasticity indicators and power quality. The objective function is to achieve the strongest comprehensive resilience of the distribution network within one operating cycle: max ⁡ F r e (23)

The distribution system satisfies power balance as Equation 24 P G , t + P W T , t + P P V , t + P E S , t = P L , t (24) where P _G,t , P _WT,t , P _PV,t , P _ES,t are the outputs of thermal power units, wind turbines, PV, and ESS in the distribution network at time t, respectively.

The ESS operating satisfies Equations 25, 26: − μ i P i , ⁡ max c h ≤ P i , t ≤ μ i P i , ⁡ max d i s e i , t + Δ t = e i , t − P i , t Δ t S O C i , ⁡ min ≤ S O C i , t ≤ S O C i , ⁡ max S O C i , 0 = S O C i , T (25) S O C i , t = S O C i , 0 + ∫ 0 τ P i . t d t E e s s , i (26) where μ _i is a 0-1 variable, in which μ _i = 1 indicates that ESS is configured at node i, and μi = 0 indicates that ESS is not configured at node i; P ^ch _i,max , P ^dis _i,max are the maximum values of ESS charging and discharging power at node i, respectively; P _i,t is the exchange power between ESS and distribution network at node i at time t; e _i,t , e _i,t+Δt are the capacities of ESS at node i at time t, t+Δt, respectively; SOC _i,t , SOC _i,0, SOC _i,T are the state of charge of the ESS at node i at time t, 0, and T, respectively; SOC _i,min, SOCi _,max are the upper and lower limits of the state of charge of the ESS at node i, respectively; E _ess,i is the capacity of ESS at node i.

Node power constraints as Equations 27, 28: p j = ∑ k ∈ δ j P j k − ∑ i ∈ π j P i j − I i j 2 r i j + g j V j 2 (27) q j = ∑ k ∈ δ j Q j k − ∑ i ∈ π j Q i j − I i j 2 r i j + b j V j 2 (28)

Ohm’s law constraint (Equation 29): V j 2 = V i 2 − 2 P i j r i j + Q i j x i j + I i j 2 r i j 2 + x i j 2 (29)

Node voltage constraint as Equation 30: V i , ⁡ min ≤ V i ≤ V i , ⁡ max (30) where V _i,min and V _i,max are the minimum and maximum voltage amplitudes operating at node i.

Branch current constraint as Equation 31: 0 ≤ I i j 2 ≤ I i j , ⁡ max 2 (31) where I _ij,t is the current flowing through branch ij; I _ij,max is the maximum current allowed by branch ij.

To ensure the safe operation of the distribution network, the comprehensive resilience index should meet the constraint conditions of not exceeding the limit as Equation 32: F r e ≤ F ¯ (32) where F ¯ is the comprehensive resilience indicator limits.

Perform second-order cone relaxation transformation on the ES optimization configuration model with X _ij = I² _ij and Y _j = V² _j as Equations 33–36: p j = ∑ k ∈ δ j P j k − ∑ i ∈ π j P i j − X i j r i j + g j Y j q j = ∑ k ∈ δ j Q j k − ∑ i ∈ π j Q i j − X i j r i j + b j Y j (33) Y j = Y i − 2 P i j r i j + Q i j x i j + X i j r i j 2 + x i j 2 (34) X i j Y j ≥ P i j 2 + Q i j 2 (35) V i , ⁡ min ≤ Y i ≤ V i , ⁡ max 0 ≤ X i j ≤ I i j , ⁡ max 2 (36)

The cost of network loss can be linearized as Equation 37: C l o s s = ∑ i , j ∈ Ω b ∑ t = 1 T C l o s s X i j , t r i j 0 ≤ X i j ≤ I i j , ⁡ max 2 (37)

Nonlinear constraints containing 0-1 variables can be transformed using the Big-M method as Equation 38: − A 1 ≤ P i , t ≤ A 2 − M 1 1 − μ i + P i , ⁡ max c h ≤ A 1 ≤ M 1 1 − μ i + P i , ⁡ max d i s − M 2 1 − μ i + P i , ⁡ max c h ≤ A 2 ≤ M 2 1 − μ i + P i , ⁡ max d i s − M 1 μ i ≤ A 1 ≤ M 1 μ i − M 2 μ i ≤ A 2 ≤ M 2 μ i (38) where M ₁ and M ₂ are relatively large constants; A ₁ and A ₂ are auxiliary variables.