To study the propagation mechanism of transient energy at the fault port in the EES during the occurrence of a fault, this paper establishes an unconstrained propagation model for transient energy: N x x ¨ + D x , x ˙ x ˙ + H x + G x , x ˙ = u (1) where x ∈ R n , x ˙ ∈ R n , a n d x ¨ ∈ R n are the state phasors of each node, the derivative of state phasors, and the second derivative of the corresponding system in the transmission process of the fault port energy under the transient fault state, respectively; N x ∈ R n × n , D x , x ˙ x ˙ , H x ∈ R n , and G x , x ˙ ∈ R n are the inertia influence relation matrix, forward energy propagation term, reverse energy propagation term, and energy loss term when the fault port energy is transmitted in the interconnected system network under the transient fault state, respectively; and u ∈ R n is the control variable of the influence of power supply and load connected to each node on the energy of the fault port in the transient fault state.

The energy action model of node s can be transformed into the following form: { N i x s x ¨ s + D s x s , x ˙ s x ˙ s + H s x ˙ s + G s x s , x ˙ s + A s x s , x ˙ s , x ¨ s = u s A i x s , x ˙ s , x ¨ s = ∑ r = 1 , r ≠ s n N s r x s x ¨ s + N s s x s − N s x s x ¨ s + ∑ r = 1 , r ≠ s n D s r x s , x ˙ r x s , x ˙ r + D s r x s , x ˙ s − D s x s , x ˙ s x ˙ s + H s ¯ x s − H s x s (2) where x s , x ˙ s , x ¨ s , H s ¯ x s , G s x s , x ˙ s , and u s represent the s-th component of vectors x , x ˙ , x ¨ , H ¯ x , G x , x ˙ , and u , respectively; N s r x and D s r x , x ˙ are the s-th and r-th components of the matrix N x and D x , x ˙ , respectively; and A s x , x ˙ , x ¨ ∈ R is the energy transfer subsystem cross-linking term of node s.

Setting x s = x s 1   x s 2 T = x s   x ˙ s T s = 1,2 , ⋯ , n , the port transient energy transitive relation network of the aforementioned formula can be transformed into the following form: x ˙ s = B s x s + C s f s x s , x ˙ s + g s x s u s + h s x s , x ˙ s , x ¨ s y ˙ s = D s x s (3) where x s is the state phasor of node s, y s is the output of node s, and B s = 0   1 0   0 ; C s = 0 1 ;   D s = 1   0 0   1 ; f s x s , x ˙ s = N s − 1 x s − D s x s , x ˙ s x ˙ s − H s x s − G s x s , x ˙ s ; g s x s = N s − 1 x s ; h s x s , x ˙ s , x ¨ s = − N s − 1 x s A s x s , x ˙ s , x ¨ s (4)

During the transient energy balance control process of the EES fault port, when constrained by load fluctuations, wind and photovoltaic output fluctuations, transmission line capacity, and other constraints, each constraint condition is uniformly represented as follows: λ x , t = 0 (5) where x ∈ R n represents the variables of each node in the EES; λ : R n → R m represents the state variable constraint function corresponding to each node; and m is the dimension of the constraint condition acting on the transient energy transmission of the port in the transient energy balance control process.

The derivative of the aforementioned equation is taken to obtain λ ˙ x , t = ∂ λ x , t ∂ x x ˙ + ∂ λ x , t ∂ t (6)

The definition is as follows: K λ x , t = ∂ λ x , t ∂ x = ∂ λ ∂ x Q ˙ t = ∂ λ x , t ∂ t (7) where K λ x , t is the m × n -dimensional Jacobian matrix, Q ˙ t is the change vector of the state variable constraint condition, and its size depends on the change rate of the state variable. Eq. 6 is represented as follows: λ ˙ x , t = K λ x , t x ˙ + Q ˙ t (8)

The relationship between the effect of each state constraint on the transient energy transfer path and mode of the port p and its corresponding x in the spatial coordinate system composed of each node in the EES is as follows: p = I x (9)

Eq. 5 can be rewritten as follows: λ p , t = λ I x , t (10)

At this time, the Jacobian matrix is as follows: K = K λ x , t = ∂ λ ∂ p ∂ I x ∂ x (11)

Assuming δ λ p as the contribution of the constraint condition to the changes in the transient energy transfer joint path and mode of the fault port, the following can be obtained: δ λ = δ λ δ p δ p = 0 (12)

Taking the Lagrange multiplier f ∈ R m into account, we can obtain the following: δ λ δ p δ p T f = 0 (13)

Set δ p T G 2 = 0 (14) where G 2 represents the action of multidimensional constraints on the change of the transient energy transfer path and mode at fault ports; δ p is the variation of various parameters under the action of G 2 after the transient energy propagation of the fault port ends. From Eqs 13, 14, it can be concluded that δ p T G 2 − δ λ δ p δ p T f = 0 (15)

Transforming the state constraint G 2 into the system state space of the transient energy balance control process at the fault port, we obtain the following: G 1 = K T f = K λ T x , t f (16)

In summary, the system model for the transient energy balance control process of fault ports considering n state constraints is as follows: N x x ¨ + D x , x ˙ x ˙ + H x + G x , x ˙ = u + K Φ T x , t f (17)

The variables of each node in the network defined in the aforementioned system are as follows: x = x 1 x 2 , x 1 ∈ R n − m , x 2 ∈ R m (18)

Eq. 18 is substituted into Eq. 5 to obtain Φ x 1 , Ω x , t , t = 0 (19) where x 2 = Ω x 1 , t , and Eq. 18 can be expressed by variable x 1 as x = x 1 Ω x 1 , t (20)

From the derivation of Formula 20, we obtain x ˙ = x ˙ 1 ∂ Ω x 1 , t ∂ x 1 + ∂ Ω x 1 , t ∂ t = J n − m   0 ∂ Ω x 1 , t ∂ x 1   J m x ˙ 1 0 + 0 ∂ Ω x 1 , t ∂ t = U θ ˙ + I (21) where U = J n − m   0 ∂ Ω x 1 , t ∂ x 1   J m ∈ R n × n , θ ˙ = x ˙ 1 0 ∈ R n , and I = 0 ∂ Ω x 1 , t ∂ t ∈ R n .

Eq. 21 can be obtained by calculating the second derivative of x: x ¨ = U θ ¨ + U ˙ θ ˙ + I ˙ . (22)

From Eqs 20–22, it can be concluded that N x U θ ¨ + U ˙ θ ˙ + I · + D x , x ˙ U θ ˙ + I ˙ + H x + G x , x ˙ = u + K Φ T x , t f (23)

By separating the state constraint terms in the network node variable parameters, the node s model can be obtained as follows: N s x s x ¨ s + D s x s , x ˙ s x ˙ s + H s x s + G s x s , x ˙ s + A s x s , x ˙ s , x ¨ s − τ Φ s * = u s (24) where A s x s , x ˙ s , x ¨ s = ∑ r = 1 r ≠ s n N s r x U E x ¨ s r + U ˙ E x ˙ s r + I ˙ r + N s r x U E x ¨ s r + U ˙ E x ˙ s r + I ˙ s − N s x r x ¨ s + ∑ r = 1 r ≠ s n D r s x s , x ˙ r U E x ˙ s r + I r + D s r x , x ˙ s U E x ˙ s r + I r − D s x s , x ˙ s x ˙ s + H s ¯ x − H s x s + τ Φ s * − τ Φ s (25)

We establish an energy correlation model for adjacent nodes in the EES, namely, the degree of topological overlap c s r x s , x r : c s r x s , x r = b s r 1 − b s r (30) b s r = | exp { − 2 × x s − x r 2 2 ρ max r ∈ N s x s − x r 2 2 } | γ (31) where x s and x r are the state variables of two adjacent nodes which provide direct energy exchange to each other, s ≠ r ; b s r represents the energy exchange level; x s − x r 2 represents the Euclidean distance of adjacent nodes; ρ refers to the homogenization element of heterogeneous nodes connected to hybrid energy; and γ is the penalty element for the energy interaction exceeding limit.

Let w s represent a set of nodes within the partition area where node s is located that has direct energy injection or cascading energy interaction, and ϑ s y s is the set of adjoining nodes that provide direct energy exchange with set w s . Then, the energy interaction degree with high-order topological prior for nodes c w s x s , x ϑ s is as follows: c w s x s , x ϑ s = c N 1 x s , x N 1 + c N 2 x s , x N 2 + ⋯ + c N i x s , x N i , (32) where c w s x s , x ϑ s represents all energy correlation values c N i related to node s added together in the partition area where node s is located; c N 1 , c N 2 , … , c N i indicates the energy correlation between node s and all adjacent nodes that have an energy interaction.

In summary, we establish a high-order prior energy model for node energy correlation: E h x w | ϒ = ∑ s ∈ S , r ∈ N s b s r + ∑ u ≠ s , r b s u b r u min ∑ u ≠ s b s u ∑ u ≠ r b r u + 1 − b s r , (33) where E h x w | ϒ is the prior energy of higher-order topological structures in the grid region where node s is located; ϒ = ρ , γ is the higher-order prior parameters for network partitioning.

Based on this, a Gaussian likelihood estimation model for the partition area is constructed to identify the degree of energy imbalance in each partition area: P X | Y , θ = ∏ s = 1 N P x s | y s , θ ∏ r ∈ ϑ s P x r | y r , θ w y r w r w y r = ∑ s ∈ S ∑ r ∈ ϑ s x s − x r = ∑ s ∈ S ∑ r ∈ ϑ s x s + x r − 2 x s x r (34) where x r = x r | r ∈ ϑ s is the neighborhood node state variable set of node s; θ = μ l , σ l 2 l ∈ Λ , where μ l , and σ l 2 are the probability distribution mean and variance of energy imbalance within a partition area; w r is the probability calculation weight of energy imbalance within the partition area, and w r = ∑ r ∈ N s w y r . The smaller the w r , the larger the estimated value Y * of the energy imbalance within its corresponding partition area, and the weaker the energy balance ability of the partition area.

The specific steps for identifying EES energy balance weak grids based on a prior knowledge model are as follows:

Step 1:

Initialize the parameter node set w s , the homogenization element of heterogeneous nodes connected to hybrid energy ρ , and the penalty element for the energy interaction exceeding limit γ .

Step 2:

Calculate the energy correlation between node s and all adjacent nodes with the energy interaction c s r x s , x r in the partition area where node s is located.

Step 3:

Calculate the prior energy of the topological structure γ E h x w s γ of the partition area where node s is located.

Step 4:

Repeat steps 2 to 3 until s = S .

Step 5:

Sort the estimated values Y * of the energy imbalance within the partition area by adopting the Gaussian likelihood estimation model and identify the weaker energy balance ability partition area.

The flow chart for grid identification with weak energy balance capability is shown in Supplementary Figure S1.

Supplementary Figure S4 shows the equivalent power angle instability curve and its corresponding transient stability margin index change curve at node 1 of the sending-end system shown in Supplementary Figure S3, without considering the participation of energy storage in regulation, when a fault occurs at node 1 with a three-phase short circuit, and this fault is not removed within the limit removal time.

As shown in Supplementary Figure S4, when the fault occurs along the output line of channel 1 of the sending-end system and the system becomes unstable, the power angle of the system at node 1 exceeds the power angle stability limit after the first swing and the sending-end system loses synchronization with the receiving-end system. During the instability process of the sending-end system, the equivalent power angle change curve of the system at node 1 is a continuously increasing oscillation process. Meanwhile, as shown in the variation curve of the stability margin index of the system in Supplementary Figure S4, the curve shows that the stability margin index of the sending-end changes significantly during the first swing of the power angle swing of node 1, indicating that the system will lose synchronization with the receiving-end system.

According to the EES dynamic grid partitioning model, as mentioned earlier, the grid partitioning of the sending-end power grid shown in Supplementary Figure S2 is carried out, and the partitioning results are shown in Supplementary Figure S5.

Supplementary Figure S5 is shows that when a transient fault happens in the output channel of the sending-end system, due to the large startup methods of new energy sources, such as photovoltaic and wind power in the system, and the lack of energy storage, when dividing the grid of the sending-end network, the new energy sources are all divided in the same grid as hydroelectric or thermal power units to ensure that the transient energy balance characteristics within the grid meet the system stability requirements as much as possible.

When the sending-end system loses stability under the scenario of the three-phase short circuit fault at node 1 in Supplementary Figure S3, the energy balance weak grid identification method put forward in this article is adopted to calculate the estimated energy imbalance of each grid. This article mainly calculates the estimated value of energy imbalance within the grid and ranks the calculation results of multiple grids in order to identify the grid with the largest energy imbalance and the weakest energy balance ability. The comparison between the calculation results and the maximum frequency deviation of each grid during the fault time period can be seen in Supplementary Table S1.

Supplementary Table S1 shows that except for partition area 1, where the fault point is located, the variation pattern of the estimated value of energy imbalance calculated using the method put forward in this article is mostly consistent with the variation pattern of the maximum grid frequency deviation during the fault time period, which verifies the effectiveness of the energy balance weak partition area identification method mentioned earlier. In addition, when the sending-end system becomes unstable, the frequency deviation of each partition area in the system is relatively large. This also indicates that when the energy storage system is not configured, the sending-end system will not only lose synchronization with the receiving-end system when facing large transient energy injection but also cause significant energy oscillations inside the sending-end system. If the startup mode of thermal and hydroelectric units existing in the power grid is small at this time and the response speed can hardly reach the requirement level of suppressing transient energy propagation within an effective time, the sending-end system is likely to undergo splitting or even collapse.

Considering the involvement of multi-energy storage equipment involved in energy regulation, at the time the fault occurs at node 1 of channel 1 of the sending-end system in Supplementary Figure S2, the equivalent power angle instability curve and its corresponding transient stability margin index change curve at node 1 are shown in Supplementary Figure S6.

Supplementary Figure S6 shows that at the time of fault occurrence to channel 1 in the sending-end system, the system protection and safety control devices do not cut off the fault within the limit cutting time. The power angle of node 1 in the system exceeds the power angle stability limit after the first swing. However, due to the rapid absorption of transient energy by the energy storage device configured in the sending-end system after the system’s transient energy exceeds the limit, the amplitude of the power angle swing in the system decreases rapidly during the second swing and makes the subsequent oscillation process converge quickly.

Meanwhile, the variation curve of the stability margin index in Supplementary Figure S6 shows that the stability margin index in the sending-end system quickly decreases to within the stability threshold after a jump in the first swing of the power angle swing at node 1. It can be seen that after configuring battery energy storage, the ability of the sending-end system to maintain transient stability has been significantly improved, but the system will still enter an unstable state, causing significant energy impacts on the synchronous power supply, new energy power supply, and load of the sending-end system.

The grid division results when considering the participation of energy storage devices in regulation are shown in Supplementary Figure S7.

Supplementary Figure S7 shows that compared to the grid division in Supplementary Figure S5, the transient energy balance grid division results shown in Supplementary Figure S7 not only consider the support role of traditional hydropower and thermal power units for new energy sources but also consider the support role of energy storage systems for new energy sources and loads.

The energy balance weak grid identification method proposed in this article is used to calculate the estimated energy imbalance values of each grid after re-partitioning. The comparison between the calculated results and the maximum frequency deviation of each grid during the fault time period is shown in Supplementary Table S2.

Supplementary Table S2 shows that compared with the estimated energy imbalance value in Supplementary Table S1, after considering the participation of energy storage devices, the estimated values of energy imbalance within each grid have been reduced, and the energy balance ability has been improved. This once again verifies the effectiveness of the energy balance weak grid identification method mentioned previously. Moreover, when the sending-end system is unstable due to the transient energy support effect of the energy storage system, the frequency deviation of each grid area of the sending-end system increases and decreases significantly. On the other hand, the simulation results also indicate that although the energy storage system could improve the transient energy balance ability of the system to a certain extent, relying on the energy absorption or release characteristics of a simple energy storage system still cannot effectively maintain the stability circumstance of the sending-end system during serious transient fault at the outlet of the sending channel.