Reliability evaluation method for power distribution network containing network construction type energy storage and flexible resources

Abstract

The invention discloses a reliability evaluation method for a power distribution network containing network construction type energy storage and flexible resources, and the method comprises the steps: firstly constructing a wind power and photovoltaic combined output probability distribution model considering the output uncertainty and spatial-temporal correlation of new energy, and a typical load scene; then establishing a'operation-fault 'two-state shutdown model, and generating a system operation state through sequential Monte Carlo simulation; simulating the power distribution network according to the fault information, and analyzing the power output characteristics of the constructed network type energy storage and the voltage support effect of the constructed network type energy storage on the near-area distributed new energy; and finally, constructing an optimal load shedding model considering the network-forming type energy storage power support effect, calculating the fault consequence of the system, and calculating the reliability index of the system according to the fault consequence. According to the method, the improvement effect of network construction type energy storage on the reliability of the power distribution network is analyzed; and meanwhile, the accurate reliability evaluation capability is still kept in a high-permeability scene of the distributed power supply, and a theoretical support is provided for network construction type energy storage planning investment for improving the reliability of the power distribution network.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of reliability evaluation of power systems with high proportion of new energy access, and particularly relates to a distribution network reliability evaluation method with grid-forming energy storage and flexible resources. BACKGROUND

[0002] Massive distributed new energy and flexible load access to distribution networks increase the randomness of operation and voltage out-of-limit risk while improving system flexibility. Existing distribution network reliability evaluation methods mainly include analytical methods and Monte Carlo simulation methods, both of which have significant limitations. Analytical methods (such as the minimal path method and fault mode and consequence analysis) have high computational efficiency, but they cannot quantify the dynamic support characteristics of grid-forming energy storage: on the one hand, their simplified models cannot represent the key role of active-frequency response and reactive-voltage support of energy storage virtual synchronous machine control in system transient stability; on the other hand, they ignore the correlation of new energy output (such as wind-solar complementary characteristics) and the inhibitory effect of energy storage active voltage regulation on new energy anti-disconnection, resulting in optimistic evaluation results. Actual cases show that a coastal distribution network with high photovoltaic penetration rate underestimated the amount of load shedding by 32% due to the neglect of energy storage support.

[0003] Current distribution network reliability evaluation models can handle the randomness of wind and solar output, but still have two major bottlenecks: first, the new energy disconnection model relies on fixed thresholds and does not consider the ability of grid-forming energy storage to maintain grid-connected point voltage through dynamic reactive power increase during faults, resulting in a high disconnection misjudgment rate of up to 30%; second, there is a lack of explicit correlation mechanism between energy storage state and system reliability, and a coupling model between energy storage over-discharge risk and load shedding has not been established. Although recent research has introduced Copula functions to characterize wind-solar correlation and used mixed integer programming (MILP) to optimize load shedding, it has not solved the core contradiction: how to accurately represent the new energy disconnection problem and the reactive power and voltage support capacity of grid-forming energy storage.

[0004] In addition, although intelligent optimization algorithms (such as deep reinforcement learning) attempt to break through, there is still a problem of missing models: existing architectures do not model the "power-voltage-energy" three-dimensional constraints of grid-forming energy storage, the spatio-temporal correlation of wind and solar, and the topology reconstruction dynamics. Especially in scenarios with insufficient actual measurement data, the model generalization ability is severely limited. It is worth noting that although the transient support value of grid-forming energy storage has been proven, its collaborative mechanism in distribution network reliability evaluation is still a blank: grid-connected control energy storage lacks active support capability and easily causes new energy cascading disconnection during faults. Therefore, it is urgent to build an evaluation framework that integrates multi-source random scenario generation, energy storage dynamic support modeling, and efficient solution engine to provide theoretical support for the planning and investment of high-proportion new energy distribution networks. SUMMARY

[0005] Invention purposes: The purpose of the present application is to provide a power distribution network reliability evaluation method containing network-type energy storage and flexibility resources. A wind power and photovoltaic output modeling method considering the volatility and correlation of distributed new energy is proposed, and a power distribution network reliability index system containing power supply fault self-healing rate and flexible load availability coefficient is established. The evaluation method of power distribution network reliability containing network-type energy storage and distributed new energy is studied, and the influence of network-type energy storage access location and capacity on system reliability improvement is analyzed, which lays a foundation for the planning and investment of new power distribution network reliability improvement. The problem of reliability misjudgment caused by ignoring the active support ability of energy storage, the dynamic characteristics of new energy off-grid and the insufficient calculation efficiency of massive scenarios in traditional evaluation is solved.

[0006] Technical scheme: The power distribution network reliability evaluation method containing network-type energy storage and flexibility resources of the present application comprises the following steps:

[0007] Step 1, based on the historical output data of wind power and photovoltaic in the power distribution network, the wind power and photovoltaic joint output scene and the load typical scene are constructed, and the wind power and photovoltaic output data and the system typical load curve are obtained;

[0008] Step 2, for the passive elements in the power distribution network system, a "running-failure" two-state outage model of the element is established, and the running state of the power distribution network system is generated through sequential Monte Carlo simulation;

[0009] Step 3, according to the fault information in the running state of the power distribution network system, the power distribution network is simulated, the power support and voltage support of the network-type energy storage are analyzed, and the system state transition data are obtained;

[0010] Step 4, based on the wind power and photovoltaic output data and the system typical load curve of step 1 and the system state transition data of step 3, an optimal load shedding model considering network-type energy storage is constructed, and the fault result of the power distribution network system is calculated;

[0011] Step 5, the reliability index of the power distribution network system is calculated according to the fault result.

[0012] Further, step 1 specifically comprises the following steps:

[0013] Step 1.1, collect the historical output data of wind power and photovoltaic in the power distribution network;

[0014] Step 1.2, estimate the probability density of wind and light output based on the non-parametric kernel density function of Gaussian kernel. The non-parametric kernel density function of Gaussian kernel is expressed as:

[0015]

[0016] In the formula, P wind and P pv are the historical sample values of wind power and photovoltaic output respectively; Xwind and X pv f(X) represents the output variables of wind power and solar power, respectively; wind ) and f(X pv ) represent the data samples in variable X respectively. wind and X pv The probability density function estimate at point n; n is the number of samples; h is the bandwidth; and K(·) is the Gaussian kernel function.

[0017] According to probability theory, f(X) wind ) and f(X pv The cumulative distribution function F of wind power output was calculated. x (X wind ) and cumulative distribution function of photovoltaic output F y (X PV );

[0018] Step 1.3: Establish the joint distribution function of wind and solar power output for each hour (24 hours per day) using the Frank-Copula function, as shown in equation (2):

[0019] F(X wind ,X pv ,λ)=C[F x (X wind ),F y (X pv ),λ] (2)

[0020] Among them, F x (X wind ) and F y (X pv ) are the cumulative probability distribution functions of wind power and solar power output, respectively, and C is the Frank-Copula function; λ represents the correlation coefficient between wind power output and solar power output: λ>0 indicates that F x (X wind ) and F y (X pv A positive correlation indicates F; conversely, a negative correlation indicates F. x (X wind ) and F y (X pv The negative correlation is given by λ∈(-1,1) and λ≠0.

[0021] The Frank-Copula function is represented as:

[0022]

[0023] In the formula, u and v represent the two independent variables of the binary Frank-Copula function, i.e., u = F x (X wind ), v = Fy (X pv )。

[0024] The cumulative probability distribution function time period of wind power and photovoltaic power output is divided into 24, the cumulative probability distribution of wind power and photovoltaic power output in the first hour period is brought into the Frank-Copula function, and the correlation coefficient λ1 of the cumulative probability distribution of wind power and photovoltaic power output in the first hour period is solved by using the maximum likelihood estimation method, and the correlation coefficients λ2-λ 24 .

[0025] Step 1.4, generating wind power and photovoltaic power output data considering wind-solar correlation, and establishing a wind-solar joint output matrix H. The specific steps are as follows:

[0026] Step 1.4.1, in the first period, two random numbers s1 and s2 are generated in the interval [0-1];

[0027] Step 1.4.2, let u in the Frank-Copula function be s1, and the value v is solved according to formula (4), that is:

[0028]

[0029] Step 1.4.3, since u=F x (X wind ), v=F y (X pv ), the output data of wind power and photovoltaic power in the first period is calculated by using formula (5):

[0030]

[0031] Step 1.4.4, repeat steps 1.4.1-1.4.2 for m times to obtain the joint output data of wind power and photovoltaic power in the first period under m scenarios;

[0032] Step 1.4.5, repeat steps 1.4.1-1.4.4 to calculate the second period (formula (4) uses λ2), the third period (formula (4) uses λ3), …, and the 24th period (formula (4) uses λ 24 ), each period has m groups of wind-solar joint output scenarios;

[0033] Step 1.4.6, establish the wind-solar joint output matrix H as shown in formula (6).

[0034]

[0035] Step 1.5, the wind and light combined output scene is classified by using K-means++ clustering algorithm, and the typical scenes of wind power and photovoltaic output are obtained, specifically:

[0036] Randomly select the ith row h i of the wind and light combined output matrix H k , i = [1, 2, …, N]; k = [1, 2, …, l], wherein, l is the number of cluster centers;

[0037] The shortest distance D(H) between each output scene and the existing cluster center is calculated:

[0038]

[0039] In the formula, j = [1, 2, …, 48]; h ij represents the ith row and jth column element in the output matrix H, h kj represents the kth row and jth column element in the output matrix H;

[0040] The probability P(H) of each output scene being selected as the next cluster center is calculated, and the node corresponding to the maximum probability value is selected as the next cluster center;

[0041]

[0042] Repeat formula (7) and formula (8) until l cluster centers required by the wind and light combined output scene are selected, and find the minimum distance of each node to each cluster center based on D(H), divide the node and the cluster center with the shortest distance into a region, and finally obtain the typical wind and light combined output scene.

[0043] Further, step 2 is: the passive elements include transformers and lines, and the “operation-failure” two-state outage model of the transformers and lines is established, the continuous working time M TTR and the average repair time M TTF are defined as:

[0044] M TTR = 1 / δ (9-1)

[0045] M TTF = 1 / μ (9-2)

[0046] Wherein, δ is the failure rate of the element, μ is the repair rate of the element;

[0047] Repeat sampling the duration of each element state within the given simulation time to obtain the state transition process of all elements, combine the state transition processes of all elements in the system to establish the state transition process of the system, and the probability of the system being in a state X j ​ frequency and duration The specific calculations are as follows:

[0048]

[0049] In the formula: T is the given simulation duration, and D is the sum of the durations of the system in all states; dk For system state X j The duration of the kth occurrence; N d Let X be the system state within a given simulation duration T. j Number of times it appears.

[0050] Furthermore, step 3 specifically involves: analyzing the results of each system state in chronological order; when the system state includes a fault, simulating the distribution network based on the specific fault information to obtain the voltage drop at each node of the distribution network during the transient process; and determining whether any new energy units have disconnected from the grid by judging whether the voltage at the distributed new energy port is lower than the low voltage threshold; wherein, the voltage threshold is 0.8U. N U N This is the rated voltage.

[0051] Furthermore, step 4 specifically includes the following steps:

[0052] Step 4.1: Obtain basic data of the distribution network based on the example system, including network structure and parameters, rated capacity and parameters of distributed new energy sources, capacity and parameters of grid-type energy storage, load weight coefficients of each node, and load power; Based on the wind power and photovoltaic output data obtained in Step 1, the typical load curve of the system, the system state transition data obtained in Step 3, and the new energy disconnection information in each system state, construct an optimal load shedding model considering grid-type energy storage.

[0053] Step 4.2, the optimal load shedding objective function is:

[0054]

[0055] In the formula, and This represents the load level weighting coefficient and load power of the i-th node at time t; This represents the load reduction variable at node i at time t. 1 indicates that the load is completely cut off, and 0 indicates that the load is not cut off;

[0056] The constraints include: power system network parameter constraints, system power flow constraints, distribution network topology radial constraints, distributed renewable energy power and voltage constraints, and grid-based energy storage power and voltage constraints.

[0057] Step 4.3. The optimal cut load model is a mixed integer linear programming model MILP, which is solved by calling the solver CPLEX in Matlab to obtain the optimal cut load results.

[0058] Further, in step 4.2, the constraint conditions are specifically as follows:

[0059] Power system network parameter constraints:

[0060]

[0061] In the formula: variable represents the active power flowing into node j of line i-j at time t; and represent the active power and the reactive power of the substation at time t; represents the reactive power flowing into node j of line i-j at time t; V i t represents the voltage of node j at time t; represents the 0 / 1 variable of whether line i-j is faulty at time t; M is a relaxation large integer; V i min and V i max respectively represent the lower limit and the upper limit of the node voltage; parameter r ij and x ij represent the line impedance; represents the line capacity; represents the upper limit of the substation output; Ω L , Ω s and Ω B respectively represent the branch set, the substation set and the node set;

[0062] System power flow constraints

[0063]

[0064] In the formula: variable represents the active power of the distributed power supply at time t; and represent the active power flowing out of node j of line j-k, the discharging and charging power of the energy storage system at time t; and respectively represent the reactive power flowing out of node j of line j-k, the reactive power injection of the DG and the reactive power release of the energy storage at time t; represents the current on line i-j at time t; α j and β j respectively represent the branch set with node j as the terminal node and the line set with node j as the head node; Ω DG , Ω ESSDG and ES represent the set of distributed generators and energy storage respectively;

[0065] Distribution network topology radial constraint

[0066]

[0067] wherein: and is the virtual power of line i-j flowing into node j and the virtual power of line j-k flowing out of node j; is a 0 / 1 variable representing the type of subgraph that node i belongs to, if node i belongs to the subgraph powered only by distributed generators and passive island subgraph is 1, otherwise 0; in the virtual network, the nodes with 1 are regarded as the source nodes of subgraph, set the virtual power output without restriction, N node is the number of system nodes; the set of potential root nodes Ω Sr is the set of substation, distributed generator nodes and the first and last nodes of fault lines, the rest of the nodes are virtual loads with value 1;

[0068] Distributed new energy power and voltage constraint:

[0069]

[0070] wherein: is the upper limit of DG capacity; is the under-voltage protection threshold of DG, λ pf is the power factor of distributed generator, is the operating state of DG at node i at time t, if the DG at node i has grid-connected operation capability, is 1;

[0071] Grid-forming energy storage power and voltage constraint

[0072]

[0073]

[0074] wherein: E represents the effective value of grid-forming energy storage output voltage; is the effective value of node voltage; δ i,t is the power angle; X f represents the filter reactance; is the state of charge of energy storage at time t; is the capacity value of energy storage system; η d and η c respectively represent the discharge and charging efficiency of energy storage at time t; and respectively represent the state of energy storage discharge and charging at node i, assuming that the energy storage is in the discharge state, is in the state of 1, is in the state of 0; Δt is the energy storage charging and discharging time step.

[0075] Further, step 5 is specifically: calculating the reliability index of the system according to the fault consequences, and the index includes:

[0076] System average interruption frequency index SAIFI:

[0077]

[0078] In the formula: f i and N c,i respectively represent the annual average outage frequency of the load node i and the total number of users; Ω B represents the set of load nodes of the distribution network;

[0079] System average interruption duration index SAIDI:

[0080]

[0081] In the formula: K i represents the annual average outage hours of the load node i;

[0082] System average service availability index ASAI:

[0083]

[0084] System energy not supplied expectation EENS:

[0085]

[0086] In the formula: P i,cut represents the average load power shortage of node i; T c,i represents the average duration of power shortage of node i.

[0087] The application also discloses a distribution network reliability evaluation system containing network-type energy storage and distributed new energy, which comprises:

[0088] The design unit: based on Frank-Copula, wind and light are jointly modeled to obtain typical output scenarios considering wind and light randomness and correlation;

[0089] The simulation unit: the transient simulation of the distribution network containing network-type energy storage is carried out to obtain the transient voltage distribution results of each node of the distribution network;

[0090] The modeling unit: the minimum loss of load objective function and the constraint condition are constructed;

[0091] Optimization solving unit: the optimal load shedding model is solved by calling the commercial solver CPLEX in Matlab, and the optimal load shedding result is obtained.

[0092] Calculation unit: according to the optimal load shedding result, the system reliability index is calculated.

[0093] The application also discloses a computer device, including a memory, a processor and a computer program stored in the memory, and the processor executes the computer program to realize the steps of the method.

[0094] The application also discloses a computer readable storage medium, which stores a computer program / instruction, and the computer program / instruction is executed by a processor to realize the steps of the method.

[0095] The application also discloses a computer program product, which includes a computer program / instruction, and the computer program / instruction is executed by a processor to realize the steps of the method.

[0096] Beneficial effects: compared with the prior art, the application has the following remarkable advantages:

[0097] (1) the application breaks through the limitation of traditional methods ignoring the space-time correlation of new energy by constructing a wind-solar joint probability model based on a Frank-Copula function, and realizes high-precision generation of output scenarios (typical scenario clustering error ≤ 5%).

[0098] (2) the application establishes a network-type energy storage "power-voltage-energy" three-dimensional constraint model, and for the first time, the active frequency regulation and reactive voltage regulation characteristics of virtual synchronous machine control are integrated into reliability evaluation. A collaborative optimization load shedding strategy is proposed: embedding the energy storage transient voltage support mechanism in the mixed integer programming (MILP), solving the 30% new energy off-grid misjudgment problem caused by the fixed threshold. Combined with K-means scenario clustering, 90% of redundant calculation is reduced, and the time consumption of million-level scenario evaluation is compressed from more than 8 hours to 1.5 hours, and the calculation efficiency is improved by 5.3 times.

[0099] (3) by quantifying the mapping relationship between network-type energy storage capacity and reliability, it is found that when the energy storage capacity configuration exceeds 80% of the new energy installed capacity, the system reliability can be overall superior to the field with extensive configuration of low voltage ride through control. The energy storage configuration criterion based on reliability-economy collaborative optimization is proposed, and it is pointed out that deploying network-type energy storage at key nodes can significantly improve the island continuous power supply capability without modifying existing new energy equipment, and avoid the additional investment brought by low voltage ride through function. The empirical results show that the reliability evaluation error of the method is reduced from more than 15% of the traditional method to 4.8%, which provides accurate and efficient decision basis for scientific planning and economic deployment of energy storage in high-proportion new energy distribution networks. BRIEF DESCRIPTION OF DRAWINGS

[0100] Figure 1 For the method flowchart of the application;

[0101] Figure 2 For the wind and light combined typical scene; wherein, (a) is the typical output diagram of wind power; (b) is the typical output diagram of photovoltaic;

[0102] Figure 3 For the system operation sample generation diagram;

[0103] Figure 4 For the grid-connected control model of network-type energy storage;

[0104] Figure 5 For the grid-connected control model of network-type energy storage;

[0105] Figure 6 For the transient voltage support characteristics of network-type energy storage; (a) is the dynamic change of network-type energy storage reactive power output under different reactive power integral coefficients; (b) is the dynamic change of network-type energy storage grid-connected point voltage under different reactive power integral coefficients;

[0106] Figure 7 For the improved IEEE33 node distribution network system diagram. DETAILED DESCRIPTION

[0107] The technical solutions of the application will be further described below with reference to the drawings.

[0108] As Figure 1 shown, a power distribution network reliability evaluation method containing network-type energy storage and distributed new energy of the application comprises the following steps:

[0109] Step 1, generating wind power, photovoltaic combined output scene and load typical scene;

[0110] Step 2, for line, transformer and other passive elements, establishing element “operation-failure” two-state outage model, generating system operation state through sequential Monte Carlo simulation;

[0111] Step 3, simulating the distribution network for fault information, analyzing the power support and voltage support of network-type energy storage;

[0112] Step 4, constructing optimal load shedding model considering network-type energy storage, calculating the fault consequences of the system;

[0113] Step 5, calculating the reliability index of the system according to the fault consequences.

[0114] Further, step 1, generating wind power, photovoltaic combined output scene and load typical scene;

[0115] Based on the historical output data of wind turbines and photovoltaics in a certain region, the output data of wind turbines and photovoltaics in different scenarios in the same region are obtained, and an excel file of wind-solar combined output is obtained.

[0116] According to the historical data in the wind-solar combined output data excel file, the wind-solar output probability density is estimated based on the non-parametric kernel density function of Gaussian kernel. The expression formula of the non-parametric kernel density function of Gaussian kernel is:

[0117]

[0118] In the formula, P wind and P pv are the historical sample values of wind power and photovoltaic output respectively; X wind and X pv represent the wind power and photovoltaic output variables respectively; f(X wind ) and f(X pv ) represent the probability density function estimates of the data samples at the variables X wind and X pv respectively; n is the sample number; h is the bandwidth, and K(·) is the Gaussian kernel function.

[0119] According to the probability density of wind power and photovoltaic output, the cumulative distribution functions of wind power and photovoltaic in each period are solved, and the wind-solar combined distribution function is established through Frank-Copula function:

[0120] F(X wind ,X pv ,λ)=C[F x (X wind ),F y (X pv ),λ] (2)

[0121] Wherein, F x (X wind ) and F y (X pv ) are the cumulative probability distribution functions of wind power and photovoltaic output respectively, and C is the Frank-Copula function; λ represents the correlation coefficient of wind power output and photovoltaic output: λ>0 represents that F x (X wind ) and F y (X pv ) are positively correlated; otherwise, F x (X wind ) and F y (X pv ) are negatively correlated, λ∈(-1,1) and λ≠0.

[0122] The Frank-Copula function is expressed as:

[0123]

[0124] where u and v represent two independent variables of the binary Frank-Copula function, i.e. u=F x (X wind ), v=F y (X pv ).

[0125] The cumulative probability distribution function time period of wind power and photovoltaic power output is divided into 24, the cumulative probability distribution of wind power and photovoltaic power output in the first hour period is brought into the Frank-Copula function, and the correlation coefficient λ1 of the cumulative probability distribution of wind power and photovoltaic power output in the first hour period is solved by using the maximum likelihood estimation method, and the correlation coefficients λ2-λ 24 of the cumulative probability distribution of wind power and photovoltaic power output in 24 time periods are obtained in turn.

[0126] The wind power and photovoltaic power output data considering the correlation of wind and light are generated, and the wind and light joint output matrix H is established:

[0127] (1) two random numbers s1 and s2 are generated in the interval [0-1] in the first period;

[0128] (2) let u=s1 in the Frank-Copula function, and the value v is solved according to formula (4):

[0129]

[0130] (3) since u=F x (X wind ), v=F y (X pv ), the output data of wind power and photovoltaic power in the first period is calculated by using formula (5):

[0131]

[0132] (4) repeat steps (1) and (2) for m times to obtain the joint output data of wind power and photovoltaic power in the first period under m scenarios;

[0133] (5) repeat steps (1), (2), (3) and (4) to calculate the second period (formula (4) uses λ2), the third period (formula (4) uses λ3), …, until the 24th period (formula (4) uses λ 24 ), each period has m groups of wind and light joint output scenarios;

[0134] (6) the wind and light joint output matrix H shown in formula (6) is established.

[0135]

[0136] Further, the K-means++ clustering algorithm is used to cluster similar output scenarios, and typical wind-solar combined output scenarios are obtained according to the clustering results. The steps of classifying wind-solar combined output scenarios by the K-Means++ algorithm are as follows:

[0137] Randomly select the ith row h i of the wind-solar combined output data matrix H as the kth initial cluster center h k , i = [1, 2, …, N]; k = [1, 2, …, 4], wherein 4 is the number of cluster centers;

[0138] Calculate the shortest distance D(H) between each node and the existing cluster center:

[0139]

[0140] In the formula, j = [1, 2, …, 48]; h ij denotes the jth column element in the ith row of the output data matrix H, h kj denotes the jth column element in the kth row of the output data matrix H, and here m = 1000.

[0141] Calculate the probability P(H) of each node being selected as the next cluster center, and select the node corresponding to the maximum probability value as the next cluster center;

[0142]

[0143] Repeat formula (7) and formula (8) until 4 cluster centers required for wind-solar combined output scenarios are selected, and find the minimum distance of each node to each cluster center based on D(H), divide the node and the cluster center with the shortest distance into a region, as shown in formula (9), and finally obtain the results of 4 typical wind-solar combined output scenarios. Figure 2

[0144] Step 2: For passive elements such as lines and transformers, establish an element “operation-failure” two-state outage model, and generate system operation states by sequential Monte Carlo simulation; the system operation state generation process is shown in formula (10). Figure 3

[0145] Establish the “operation-failure” two-state outage model of transformers and lines, and define the continuous working time M TTR and the average repair time M TTF as:

[0146] M TTR = 1 / δ (9-1)

[0147] M TTF = 1 / μ (9-2)​​

[0148] Where, δ is the failure rate of the element, μ is the repair rate of the element. The line failure rate is 0.0015 times / (year·km), and the failure repair time is 5 hours / km.

[0149] The duration of each element state is repeatedly sampled within a given simulation time limit to obtain the state transition process of all elements. The state transition processes of all elements in the system are combined to establish the state transition process of the system. The probability of the system being in a certain state X j Frequency and duration The specific calculation is as follows:

[0150]

[0151] In the formula: T is the given simulation time, which is the sum of the durations of the system in all states; D dk is the duration of the kth occurrence of the system state X j ; N d is the number of occurrences of the system state X j within the given simulation time limit T.

[0152] Step 3: Simulate the distribution network according to the fault information, analyze the power support and voltage support of the grid-type energy storage, and the control topology of the grid-type energy storage is as Figure 4 shown. Wherein: U and U ref are the output voltages of the converter; θ and E are the phase angle and voltage of the active control and reactive control outputs, respectively; i sd and i sq are the d-axis and q-axis components of the PCC point current; C dc is the DC side bus capacitor; U dc is the DC side bus voltage; R f and L f are the filter resistance and filter inductance, respectively; C f is the filter capacitance; R g and L g are the equivalent inductance and resistance of the AC power grid, respectively; P ref and Q ref are the given values of active power and reactive power; U t is the voltage of the new energy PCC point; u td and u tq are the d-axis and q-axis components of the PCC point voltage; θ PLL is the PLL output phase angle; ω is the PLL output angular velocity; i dref and i qref are the d-axis and q-axis components of the given value of the converter output current; e​dref and e qref are the d-axis and q-axis components of the converter output voltage reference, respectively.

[0153] The grid-forming energy storage realizes the active support of the distribution network through active-frequency control and reactive-voltage control, including:

[0154] The active-frequency control is realized by simulating the dynamic characteristics of the synchronous generator rotor, and the formula is:

[0155]

[0156] The reactive-voltage control is realized by simulating the characteristics of the synchronous generator excitation system, and the formula is:

[0157]

[0158] where D p and D q are the droop coefficients of active power and reactive power, respectively; T GFM is the virtual inertia time constant; ω n and ω are the rated value and given value of angular velocity, respectively; K q is the reactive power integral coefficient.

[0159] From the external characteristics, the grid-forming converter is equivalent to a controlled voltage source with low series impedance, as shown in Figure 5 , where E ref and θ ref are the given values of the output voltage and phase angle of the grid-forming converter, respectively, and e ref is the reference value of the output voltage of the grid-forming converter. It actively builds the voltage-frequency support framework of the power grid by autonomously adjusting the amplitude and phase angle of the output voltage.

[0160] Taking the grid-forming energy storage integrated into an infinite power source as an example, suppose a voltage drop disturbance occurs on the infinite power source side, and the reactive power output of the grid-forming new energy and the voltage dynamics of the grid-connected point under different voltage adjustment coefficients k q are shown in Figure 6 .

[0161] As Figure 3 , the consequences of each system state of the above system are analyzed in chronological order. When the system state contains a fault, the distribution network is simulated according to the specific fault information to obtain the voltage drop of each node in the distribution network during the transient process. By judging whether the voltage at the port of the distributed new energy is lower than the low voltage threshold (generally 0.8U N , U N is the rated voltage), it is determined whether the new energy unit is off the grid.

[0162] Step 4: Obtain the basic data of the power distribution network according to the example system, including the network structure and parameters, the rated capacity and parameters of distributed new energy, the capacity and parameters of grid-forming energy storage, the load weight coefficient of each node, and the load power (where the load weight coefficient of each node is divided according to the user level); obtain the wind power and photovoltaic output data according to step 1, the typical load curve of the system, the system state transition process obtained in step 3, and the new energy off-grid information in each system state obtained in step 4.

[0163] The input system state sample data (including wind and light output scene, fault element, and energy storage output) is shown in Table 1, the load weight coefficient and load power are set, the load weight coefficient of each node is divided according to the user level, the input power system network parameters are converted into mathematical constraints, and the formula is as follows: Figure 3

[0164]

[0165]

[0166] In the formula, the variable represents the active power flowing into node j of line i-j at time t; and represent the active power and reactive power of the substation at time t; represents the reactive power flowing into node j of line i-j at time t; represents the voltage of node j at time t; represents the 0 / 1 variable of whether line i-j is faulty at time t; M is a large integer; V i min and V i max represent the lower limit and upper limit of the node voltage respectively; the parameters r ij and x ij represent the line impedance; represents the line capacity; represents the upper limit of the substation output. Ω L , Ω s and Ω B represent the branch set, the substation set and the node set respectively.

[0167] The physical characteristics of the system under fault are converted into a strict mathematical optimization model, the core of which is to build a multi-period mixed integer programming problem with the objective of minimizing the weighted load shedding, while embedding the power support capacity of grid-forming energy storage, the new energy off-grid mechanism, the network radial topology constraint and the linearized power flow equation, which is specifically expressed as:

[0168] (1) The objective function is:

[0169] ​

[0170] where, and denotes the level weight coefficient and the load power of the i-th node at time t; denotes the load shedding variable of the i-th node at time t, 1 denotes full load shedding, and 0 denotes no load shedding.

[0171] (2) System power flow constraints:

[0172]

[0173] where: variable denotes the active power of the distributed generator at time t; and denotes the active power of the line j-k flowing out of the node j, the discharging and charging power of the energy storage system at time t; and denotes the reactive power of the line j-k flowing out of the node j, the reactive power injected by the DG, and the reactive power released by the energy storage at time t; denotes the current on the line i-j at time t; α j and β j denote the branch set with the node j as the terminal node and the line set with the node j as the head node; Ω DG , Ω ESS denote the distributed generator set and the energy storage set, respectively.

[0174] (3) Distribution network topology radial constraint:

[0175]

[0176]

[0177] where: and are the virtual power of the line i-j flowing into the node j and the virtual power of the line j-k flowing out of the node j; denotes the 0 / 1 variable of the type of subgraph in which the node i is located, and if the node i belongs to a subgraph powered only by a distributed generator and a passive island subgraph, then is 1, otherwise it is 0. In the virtual network, the nodes with are considered as subgraph source nodes, and the virtual power output is set to be unrestricted. N node denotes the number of system nodes. The set of potential root nodes Ω Sr denotes the set of substations, distributed generator nodes, and fault line head and tail nodes, and the remaining nodes are virtual loads with a value of 1.

[0178] (4) Distributed generation power and voltage constraints:

[0179]

[0180] wherein: represents the upper limit of DG capacity; represents the DG under-voltage protection threshold, λ pf represents the power factor of the distributed power supply, represents the DG operating state at node i at time t, if the DG at node i has grid-connected operation capability, then is 1.

[0181] (5) Network-constructed energy storage power and voltage constraints:

[0182]

[0183]

[0184] wherein: E represents the effective value of the output voltage of the network-constructed energy storage; is the effective value of the node voltage; δ i,t is the power angle; X f represents the filter reactance; represents the state of charge of the energy storage at time t; is the capacity value of the energy storage system; η d and η c respectively represent the discharging and charging efficiencies of the energy storage at time t; and are 0 / 1 variables representing the discharging and charging states of the energy storage at node i, assuming that the energy storage is in a discharging state, is 1 state, is 0 state; Δt is the charging and discharging time step of the energy storage.

[0185] The mixed integer model is used to solve the reliability model of the distribution network containing network-constructed energy storage and distributed flexible resources to analyze the calculation of the distribution network system reliability index. The accuracy of the reliability evaluation of the distribution network system under different linear solution models of the distribution network reliability is analyzed.

[0186] Step 5, calculating the reliability index of the system according to the fault consequences.

[0187] The system average interruption frequency index (SAIFI), the load average interruption duration index (SAIDI), the average service availability index (ASAI), and the expected energy not supplied (EENS) are respectively defined as:

[0188] SAFIF represents the average number of power outages per user in the distribution network per unit time, and the formula is:

[0189]

[0190] f i and N c,i respectively represent the annual average outage frequency and the total number of users of load node i; Ω B represents the set of load nodes of the distribution network.

[0191] SADID represents the average outage frequency per user per unit of time in the distribution network, and the formula is:

[0192]

[0193] T c,i represents the annual average outage hours of users at load node i.

[0194] ASAI represents the ratio of the actual total power supply time per unit of time to the power supply time required by users, and the formula is:

[0195]

[0196] EENS represents the expected value of the total system power loss per unit of time due to faults in the distribution network, and the formula is:

[0197]

[0198] P i,cut represents the average load power shortage of node i; T c,i represents the average duration of power shortage of node i.

[0199] The following three scenarios are set for the reliability evaluation scheme of the distribution network: Scenario 1: The system is not configured with network-type energy storage, and the distributed new energy is not configured with low-voltage ride-through control; Scenario 2: The system is not configured with network-type energy storage, and the distributed new energy is configured with low-voltage ride-through control; Scenario 3: The distribution network is configured with network-type energy storage, and the energy storage capacity is 10% of the installed capacity of new energy, and the distributed new energy is not configured with low-voltage ride-through control. The reliability results of the system in the three scenarios are shown in Table 1.

[0200] Table 1 Reliability results of the system in three scenarios

[0201]

[0202] From Table 1, it can be seen that:

[0203] 1) The system's reliability level in Scenario 2 is significantly higher than in Scenario 1. This indicates that after configuring low-voltage ride-through control for distributed renewable energy sources, when a system fault causes the voltage level at the renewable energy source's port to fall below the threshold value, the renewable energy source can ride through the fault and continue to supply power to the load, thereby improving system reliability. The indicator with the largest improvement is Expected Energy Shortage (EENS). This is because distributed power sources are grid-connected devices. In Scenario 1, without low-voltage ride-through control, a system fault would cause a large number of distributed renewable energy sources to disconnect from the grid, potentially leading to islanded power supply and a large energy shortage.

[0204] 2) The system's reliability level in scenario 3 is higher than in scenario 1 but lower than in scenario 2, indicating that the configuration of grid-based energy storage has a significant effect on improving system reliability. Its main functions are: 1. Grid-based energy storage can provide transient voltage support during faults, maintaining its own voltage level and that of nearby renewable energy grid connection points, ensuring that renewable energy sources without low-voltage ride-through configuration do not disconnect from the grid, thus continuously supplying power to the load; 2. Grid-based energy storage can also provide active power support during and after faults, further enhancing the system's power recovery capability.

[0205] While in terms of pure technical indicators, configuring grid-based energy storage (Scenario 3) is slightly less effective than fully configuring low-voltage ride-through control (Scenario 2), the former typically has a significant advantage from an economic perspective. Configuring low-voltage ride-through control requires retrofitting all renewable energy units, resulting in high investment costs; whereas grid-based energy storage only requires retrofitting the storage equipment with grid-based control. Furthermore, in addition to improving reliability, grid-based energy storage can participate in various ancillary services such as peak shaving and frequency regulation, offering multiple benefits to the power grid.

[0206] The reliability results of the distribution network with different grid-type energy storage capacity ratios (energy storage ratio = (grid-type energy storage capacity / new energy installed capacity) × 100%) are further analyzed below, as shown in Table 2.

[0207] Table 2 Impact of different energy storage grid configurations on system reliability

[0208]

[0209] As can be seen from Table 2:

[0210] 1) The capacity of grid-type energy storage increases, and the system reliability level is higher.

[0211] 2) When the proportion of grid-based energy storage is 90% or more, the system reliability level is higher than that of the scenario where low-voltage ride-through control is fully configured for distributed renewable energy (i.e., scenario 2). This indicates that at this capacity level, the voltage and power support capabilities of grid-based energy storage have been fully utilized, which can effectively replace the low-voltage ride-through function, suppress renewable energy disconnection to a greater extent, and enhance the continuity of power supply during faults.

[0212] In a second aspect, the application further discloses a power distribution network reliability evaluation system containing network-structured energy storage and distributed flexible resources, comprising:

[0213] A modeling unit: a multi-source collaborative reliability evaluation model considering the correlation between wind and light output and the active support characteristics of energy storage is constructed;

[0214] A design unit: wind and light are jointly modeled based on Frank-Copula, then, for passive elements such as lines and transformers, a "running-failure" two-state outage model of the elements is established, and system operating states are generated through sequential Monte Carlo simulation;

[0215] A calculation unit: system state data containing typical wind and light output scenarios are input into the mixed linear reliability evaluation model of the power distribution network containing network-structured energy storage and distributed flexible resources, and system reliability evaluation indexes are calculated and analyzed.

Claims

1. A method for reliability evaluation of power distribution network with meshed energy storage and flexibility resources, characterized in that, The method comprises the following steps: Step 1, based on the historical output data of wind power and photovoltaic in the power distribution network, the wind power and photovoltaic joint output scene and the load typical scene are constructed, and the wind power and photovoltaic output data and the system typical load curve are generated; Step 2, for the passive elements in the power distribution network system, an element "operation-failure" two-state outage model is established, and the operation state of the power distribution network system is generated through sequential Monte Carlo simulation; Step 3, according to the fault information in the operation state of the power distribution network system, the power distribution network is simulated, the power support and voltage support of the network-type energy storage are analyzed, and the system state transition data is obtained; Step 4, based on the wind power and photovoltaic output data, the system typical load curve and the system state transition data in step 1 and step 3, an optimal load shedding model considering the network-type energy storage is constructed, and the fault result of the power distribution network system is calculated; Step 5, according to the fault result, the reliability index of the power distribution network system is calculated.

2. The method of claim 1, wherein, Step 1 specifically comprises the following steps: Step 1.1, collect the historical output data of wind power and photovoltaic in the power distribution network; Step 1.2, estimate the probability density of wind and light output based on the non-parametric kernel density function of Gaussian kernel, which is expressed as: where P wind and P pv are the historical sample values of wind and PV power output, respectively; X wind and X pv represent the wind and PV power output variables, respectively; f(X wind ) and f(X pv ) represent the probability density function estimates of the data samples at variables X wind and X pv , respectively; n is the sample number; h is the bandwidth, and K(·) is the Gaussian kernel function. According to the probability theory method, the wind power output cumulative distribution function F x (X wind ) and the photovoltaic output cumulative distribution function F y (X PV ) are calculated from f(X wind ) and f(X pv ). Step 1.3, the joint distribution function of wind and light output per hour is established through Frank-Copula function, see formula (2): F(X wind ,X pv ,λ) = C[F x (X wind ), F y (X pv ), λ] (2) where F x (X wind ) and F y (X pv ) are the cumulative probability distribution functions of wind power and photovoltaic power, respectively, C is the Frank-Copula function; λ represents the correlation coefficient of wind power and photovoltaic power: λ>0 indicates that F x (X wind ) and F y (X pv ) are positively correlated; otherwise, F x (X wind ) and F y (X pv ) are negatively correlated, and λ∈(-1,1) and λ≠0; The Frank-Copula function is expressed as: wherein u and v represent two arguments of a bivariate Frank-Copula function, i.e. u = F x (X wind ), v = F y (X pv ). The cumulative probability distribution function time period of wind power and photovoltaic output is divided into 24, the cumulative probability distribution of wind power and photovoltaic output in the first hour period is brought into a Frank-Copula function, and the correlation coefficient λ1 of the cumulative probability distribution of wind power and photovoltaic output in the first hour period is solved by using a maximum likelihood estimation method, and the correlation coefficients λ2-λ 24 of the cumulative probability distribution of wind power and photovoltaic output in 24 time periods are sequentially obtained. Step 1.4, generate wind power and photovoltaic output data considering wind and light correlation, and establish wind and light joint output matrix H; the specific steps are as follows: Step 1.4.1, in the first time period, generate two random numbers s1 and s2 in the interval [0-1]; Step 1.4.2, let u=s1 in the Frank-Copula function, and solve the numerical value v according to formula (4), that is: Step 1.4.3, since u = F x (X wind ), v = F y (X pv ), the output data of wind power and photovoltaic in the first period is calculated by using formula (5): Step 1.4.4, repeat steps 1.4.1-1.4.2 for m times to obtain the joint output data of wind power and photovoltaic in the first time period under m scenarios; Step 1.4.5, repeat steps 1.4.1-1.4.4 to calculate the wind and light joint output scene of the second time period, the third time period, …, and the 24th time period, each time period has m groups of wind and light joint output scenes; Step 1.4.6, establish the wind and light joint output matrix H as shown in formula (6); Step 1.5, the K-means++ clustering algorithm is used to classify the wind and light joint output scene, and the typical scene of wind power and photovoltaic output is obtained, which is specifically: Randomly select the ith row h of the wind and light combined output matrix H i As the kth initial clustering center h k , i = [1, 2, …, N]; k = [1, 2, …, l], wherein, l is the number of clustering centers; Calculate the shortest distance D(H) between each output scene and the existing clustering center: where j = [1, 2,..., 48]; h ij denotes the element in the i-th row and j-th column of the output matrix H, h kj denotes the element in the k-th row and j-th column of the output matrix H; Calculate the probability P(H) of each output scene being selected as the next clustering center, and select the node corresponding to the maximum probability value as the next clustering center; Repeat formula (7) and formula (8) until the required l clustering centers of the wind and light joint output scene are selected, and find the minimum distance of each node to each clustering center based on D(H), divide the node and the clustering center with the shortest distance into a region, and finally obtain the typical wind and light joint output scene.

3. The method of claim 1, wherein, Step 2 is specifically: passive elements include transformers and lines, to establish the "operation-failure" two-state outage model of passive elements, the continuous working time M of the element TTR and the average repair time M TTF is defined as: M TTR = 1 / δ (9-1) M TTF = 1 / μ (9-2) Wherein, δ is the failure rate of the element, and μ is the repair rate of the element; The state transition process of all elements in the system is combined to establish the state transition process of the system; the probability of the system being in a state X j frequency and duration are defined as follows:​ where: T is a given simulation time duration, which is the sum of the durations of the system in all states; D dk is the duration of the system in state X j is the duration of the kth occurrence; N d is the number of occurrences of the system in state X j during the given simulation time duration T.

4. The method of claim 1, wherein, The step 3 is specifically: analyzing each system state of the system in time sequence; when the system state contains a fault, simulating the power distribution network according to the specific fault information, obtaining the voltage drop of each node of the power distribution network in the transient process, and determining whether the new energy unit is off-grid by judging whether the voltage of the distributed new energy port is lower than the low voltage ride through threshold; wherein the voltage threshold is 0.8U N , and U N is the rated voltage.

5. The method of claim 1, wherein, Step 4 specifically comprises the following steps: Step 4.1, obtain the basic data of the power distribution network according to the example system, including the network structure and parameters, the rated capacity and parameters of distributed new energy, the capacity and parameters of grid-forming energy storage, the load weight coefficient of each node and the load power; obtain the wind power and photovoltaic output data, the typical load curve of the system, the system state transition data obtained in step 3 and the new energy off-grid information in each system state according to step 1, and construct an optimal load shedding model considering grid-forming energy storage; Step 4.2, the optimal load shedding objective function is: In the formula, and represents the level weight coefficient and the load power of the i-th node at time t; represents the load reduction variable of the i-th node at time t, where 1 represents full load shedding, and 0 represents no load shedding. The constraint conditions include: power system network parameter constraint, system power flow constraint, distribution network topology radial constraint, distributed new energy power and voltage constraint and grid-forming energy storage power and voltage constraint; Step 4.3, the optimal load shedding model is a mixed integer linear programming model MILP, which is solved by calling the solver CPLEX in Matlab to obtain the optimal load shedding result.

6. The power distribution network reliability evaluation method with meshed network type energy storage and flexibility resources according to claim 5, characterized in that, In step 4.2, the constraint conditions are specifically: Power system network parameter constraint: where variables Pij(t) denotes the active power flowing from line i-j into node j at time t; and Psub(t) and Qsub(t) denote the active power and the reactive power of the substation at time t; Qij(t) denotes the reactive power flowing from line i-j into node j at time t;V i t Vj(t) denotes the voltage of node j at time t; δij(t) is a 0 / 1 variable indicating whether line i-j is faulty at time t;M is a large integer for relaxation;V i min andV i max andV denote the lower and upper bounds of the voltage of node j, respectively; Parameter r ij and x ij denotes line impedance; denotes line capacity; denotes substation output upper limit; Ω L , Ω s , and Ω B denote the set of branches, the set of substations and the set of nodes, respectively; System power flow constraint In the formula: variable Pj(t) represents the active power of the distributed power supply at time t; and Pj(t) represents the active power of the distributed power supply at time t; and Qj(t) represents the reactive power of the distributed power supply at time t; Iij(t) represents the current on the line i-j at time t; α j and β j The first and second parts represent the branch set with node j as the end node and the line set with node j as the start node; Ω DG , Ω ESS represent the distributed power supply set and the energy storage set, respectively; Distribution network topology radial constraint where F ij t and is the virtual power of line i-j flowing into node j and the virtual power of line j-k flowing out of node j; is a 0 / 1 variable indicating the type of subgraph that node i belongs to, if node i belongs to a subgraph powered only by distributed generation and passive island subgraph is 1, otherwise 0; in the virtual network, the nodes with 1 are considered as subgraph source nodes, set the unlimited virtual power output, N node is the number of system nodes; the set of potential root nodes Ω Sr is the set of substation, distributed generation nodes and the first and last nodes of fault lines, the rest of the nodes are virtual loads with a value of 1; Distributed new energy power and voltage constraint: wherein: represents the DG capacity upper limit; represents the DG under-voltage protection threshold, λ pf represents the power factor of the distributed generator, represents the DG operating status at node i at time t, if the DG at node i has grid-connected operation capability, then is 1; Grid-forming energy storage power and voltage constraint E represents the effective value of the network-forming energy storage output voltage; is the node voltage effective value; δ i,t is the power angle; X f represents the filter reactance; represents the state of charge of the energy storage at time t; is the capacity value of the energy storage system; η d and η c respectively represent the discharging and charging efficiencies of the energy storage at time t; and are 0 / 1 variables representing the discharging and charging states of the energy storage at node i, and the energy storage is in a discharging state, is a 1 state, is a 0 state; Δt is the energy storage charging and discharging time step.

7. The method of claim 1, wherein, Step 5 is specifically: calculate the reliability indexes of the system according to the fault consequences, including: System average interruption frequency index SAIFI: where: f i and N c,i respectively represent the annual average outage frequency and the total number of users of load node i; Ω B represents the set of load nodes of the distribution network; System average interruption duration index SAIDI: where: K i represents the annual average outage hours of the load node i; System average service availability index ASAI: System energy not supplied EENS: where: P i,cut represents the average load power deficit of node i; T c,i represents the average duration of power deficit of node i.

8. A computer apparatus comprising a memory, a processor, and a computer program stored on the memory, wherein the computer program, when executed by the processor, causes the processor to perform the method of any one of claims 1 to 7. The processor executes the computer program to realize the steps of the method of claim 1.

9. A computer readable storage medium having stored thereon computer programs / instructions, characterized in that, The computer program / instructions are executed by the processor to realize the steps of the method of claim 1.

10. A computer program product comprising computer programs / instructions, characterized in that, The computer program / instructions are executed by the processor to realize the steps of the method of claim 1.

Images