STATE ESTIMATION METHOD FOR MULTI-STAGE VOLTAGE SAG

Abstract

A state estimation method for a multi-stage voltage sag comprises analyzing characteristics of multi-stage voltage sags due to different causes. A method for determining the cause of sudden change time of an amplitude of the multi-stage voltage sag is based on a jump characteristic of the sag amplitude. Calculation methods of a relay protection action matrix and a faulty line set are utilized to preliminarily obtain a faulty line set based on the fault clearing time of a relay protection device and other characteristics to effectively reduce a calculation amount of sag state estimation. Based on a substantive characteristic that different events cause a change of system impedance, a method for inferring a cause of an event in each stage of the multi-stage voltage sag is used to improve the accuracy of the state estimation. The state estimation method reduces the difficulty in applying to the multi-stage voltage sag.

Description

CROSS REFERENCE TO THE RELATED APPLICATION

This application claims priority of Chinese application No. 202210201077.X, filed on Mar. 03, 2022, the entire content of which is incorporated herein by reference.

TECHNICAL FIELD

The present disclosure relates to the technical field of voltage sag state estimation, specifically, to a state estimation method for a multi-stage voltage sag.

BACKGROUND

Existing voltage sag state estimation methods are mainly used for a single-stage rectangular sag caused by a short circuit, and there is no state estimation method suitable for a multi-stage voltage sag. Unfortunately, existing researches show that the multi-stage voltage sag occurs more often, and the multi-stage voltage sag cannot be ignored. For example, with the development of a new power system, the proportion of distributed generation (DG) devices connected to the system is getting higher. During the short circuit fault, a sag may cause a trip of a single DG device or cascading trips of a plurality of DG devices, thus forming a complex multi-stage voltage sag. In addition, staged protection commonly used in the power system, such as distance protection, may also cause the multi-stage voltage sag due to its action characteristics.

In voltage sag state estimation, a sag amplitude of an unmonitored bus needs to be estimated by using a short-circuit calculation formula. A short-circuit calculation method is a well-known method. Short-circuit calculation methods corresponding to different fault types are as follows.

It is assumed that ^m is a to-be-solved bus, a line ^l_ij between buses i and j is faulty, f is a fault point, P is a normalized distance from the point f to the bus i . A value range of p is [0, 1], representing f from i to the j according to the following expression:

$\begin{array}{cc}p=\frac{L_{if}}{L_{if}},p\in \left[0,1\right]& \text{­­­(1)}\end{array}$

Where L_if represents the distance from the bus i to the fault point f, and L_ij represents the length of the line l_ij .

Three-Phase Short Circuit Fault

A three-phase voltage amplitude caused by a three-phase fault of the bus m is as follows:

$\begin{array}{cc}f\left(Z_{mf}^1,Z_{ff}^1\right)=V_m^{\text{pref}}-\frac{Z_{mf}^1}{Z_{ff}^1+R_f}{V}_f^{\text{pref}}& \text{­­­(2)}\end{array}$

In the above formula, R_f represents fault resistance, and

$V_m^{\text{pref}}$

and

$V_f^{\text{pref}}$

respectively represent voltages of the bus ^m and the fault point f before the fault. The voltage of the bus before a sag can be obtained through power flow calculation by using a power system simulation tool, and the voltage of the fault point before the sag can be expressed as follows:

$\begin{array}{cc}{V}_f^{\text{pref}}=\left(1-k\right)V_i^{\text{pref}}+k{V}_j^{\text{pref}}& \text{­­­(3)}\end{array}$

Single-Phase-to-Ground Fault

It is assumed that the single-phase-to-ground fault occurs in phase A. A voltage amplitude of each phase of the bus ^m is as follows:

$\begin{array}{cc}\left\{\begin{array}{l}{V}_{m,A}=V_{m,A}^{\text{pref}}-\frac{Z_{mf}^0+Z_{mf}^1+Z_{mf}^2}{Z_{ff}^0+Z_{ff}^1+Z_{ff}^2+3R_f}{V}_{f,A}^{\text{pref}}\hfill \\ V_{m,B}={\alpha }^2{V}_{m,A}^{\text{pref}}-\frac{Z_{mf}^0+{\alpha }^2{Z}_{mf}^1+\alpha Z_{mf}^2}{Z_{ff}^0+Z_{ff}^1+Z_{ff}^2+3R_f}{V}_{f,A}^{\text{pref}}\hfill \\ V_{m,C}=\alpha V_{m,A}^{\text{pref}}-\frac{Z_{mf}^0+\alpha Z_{mf}^1+{\alpha }^2{Z}_{mf}^2}{Z_{ff}^0+Z_{ff}^1+Z_{ff}^2+3R_f}{V}_{f,A}^{\text{pref}}\hfill \end{array}\right)& \text{­­­(4)}\end{array}$

In the above formulas, α = e^j120° represents a rotation factor.

(3) Two-Phase Fault

It is assumed that the two-phase fault occurs in phases B and C. The voltage amplitude of each phase of the bus ^m is as follows:

$\begin{array}{cc}\left\{\begin{array}{l}{V}_{m,A}=V_{m,A}^{\text{pref}}-\frac{Z_{mf}^0+Z_{mf}^1+Z_{mf}^2}{Z_{ff}^0+Z_{ff}^1+Z_{ff}^2+3R_f}{V}_{f,A}^{\text{pref}}\hfill \\ V_{m,B}={\alpha }^2{V}_{m,A}^{\text{pref}}-\frac{Z_{mf}^0+{\alpha }^2{Z}_{mf}^1+\alpha Z_{mf}^2}{Z_{ff}^0+Z_{ff}^1+Z_{ff}^2+3R_f}{V}_{f,A}^{\text{pref}}\hfill \\ V_{m,C}=\alpha V_{m,A}^{\text{pref}}-\frac{Z_{mf}^0+\alpha Z_{mf}^1+{\alpha }^2{Z}_{mf}^2}{Z_{ff}^0+Z_{ff}^1+Z_{ff}^2+3R_f}{V}_{f,A}^{\text{pref}}\hfill \end{array}\right)& \text{­­­(5)}\end{array}$

Two-Phase-to-Ground Fault

It is assumed that the two-phase-to-ground fault occurs in phases B and C. The voltage amplitude of each phase of the bus ^m is as follows:

$\begin{array}{cc}\left\{\begin{array}{l}{V}_{m,A}=V_{m,A}^{\text{pref}}-\frac{Z_{mf}^1\left(Z_{ff}^0+Z_{ff}^2+2R_f\right)-Z_{mf}^2\left(Z_{ff}^0+R_f\right)-Z_{mf}^0\left(Z_{ff}^2+R_f\right)}{\left(Z_{ff}^0{Z}_{ff}^1+Z_{ff}^1{Z}_{ff}^2+Z_{ff}^2{Z}_{ff}^0\right)+2R_f\left(Z_{ff}^0+Z_{ff}^1+Z_{ff}^2\right)+3R_f^2}{V}_{f,A}^{\text{pref}}\hfill \\ V_{m,B}={\alpha }^2{V}_{m,A}^{\text{pref}}-\frac{{\alpha }^2{Z}_{mf}^1\left(Z_{ff}^0+Z_{ff}^2+2R_f\right)-\alpha Z_{mf}^2\left(Z_{ff}^0+R_f\right)-Z_{mf}^0\left(Z_{ff}^2+R_f\right)}{\left(Z_{ff}^0{Z}_{ff}^1+Z_{ff}^1{Z}_{ff}^2+Z_{ff}^2{Z}_{ff}^0\right)+2R_f\left(Z_{ff}^0+Z_{ff}^1+Z_{ff}^2\right)+3R_f^2}{V}_{f,A}^{\text{pref}}\hfill \\ V_{m,C}=\alpha V_{m,A}^{\text{pref}}-\frac{\alpha Z_{mf}^1\left(Z_{ff}^0+Z_{ff}^2+2R_f\right)-{\alpha }^2{Z}_{mf}^2\left(Z_{ff}^0+R_f\right)-Z_{mf}^0\left(Z_{ff}^2+R_f\right)}{\left(Z_{ff}^0{Z}_{ff}^1+Z_{ff}^1{Z}_{ff}^2+Z_{ff}^2{Z}_{ff}^0\right)+2R_f\left(Z_{ff}^0+Z_{ff}^1+Z_{ff}^2\right)+3R_f^2}{V}_{f,A}^{\text{pref}}\hfill \end{array}\right)& \text{­­­(6)}\end{array}$

The sag is a most serious disturbance to a sensitive device and cannot be avoided. Voltage sag levels at different positions in the system should be accurately reflected. Therefore, a monitoring device must be installed to monitor the voltage sag. The power grid is huge and complex. Therefore, the monitoring device can be installed only on some buses in the system due to the monitoring cost and information processing capability. How to use monitoring data from limited monitoring devices to estimate a sag level (mainly including the amplitude and duration of the voltage sag) of the unmonitored bus is an important research task, which is also referred to as the voltage sag state estimation. The existing voltage sag state estimation methods are mainly used for the single-stage sag caused by the short circuit fault and cannot be used for state estimation of the multi-stage voltage sag. As the multi-stage voltage sag occurs more often, it is important to propose a state estimation method suitable for the multi-stage voltage sag.

Difficulties of multi-stage state estimation are as follows. 1. The multi-stage voltage sag is mainly caused by the disconnection of the DG device from the power grid or trips of protection devices on two sides of a line at different times during the short circuit fault. Therefore, figuring out how to determine the causes of sags at different stages is the first challenge. 2. The position of the short circuit fault is an important factor affecting the voltage amplitude. Therefore, how to preliminarily obtain a possible faulty line set and reduce a calculation amount of MVSSE is the second challenge. 3. During the short circuit fault, the voltage amplitude at the same fault position is also affected by the quantity and sequence of DG devices disconnected from the power grid, the tripping condition of a protection device on a faulty line, and other events. Therefore, how to infer a specific cause of each stage of a sag event and improve the accuracy of the MVSSE is the third challenge.

Explanation of Terms Used in the Present Disclosure

Voltage sag: The International Institute of Electrical and Electronics Engineers (IEEE) defines voltage sag as a power quality phenomenon in which an effective value of a supply voltage drops rapidly to a value within the range of 0.1 to 0.9 p.u. and the drop lasts 0.5 cycles to 1 min.

Voltage sag state estimation: It is a method for estimating a voltage sag state (including a sag frequency, a voltage sag amplitude, or a duration of an unmonitored bus in a certain period of time) of the unmonitored bus based on the monitoring data of the limited monitoring devices. The voltage sag state in the present disclosure is a voltage sag amplitude and voltage sag duration in a single sag event.

Sudden voltage sag change point: It is a point at which a voltage amplitude in an effective value waveform of a voltage sag suddenly changes significantly (beyond a normal voltage fluctuation range). The sudden voltage sag change point is formed in the case of the short circuit fault, disconnection of the DG device from the power grid after the short circuit fault, and action of a relay protection device on one side of the faulty line after the short circuit fault. The sudden voltage sag change point can be determined through waveform point detection.

SUMMARY

To resolve the above problems, the present disclosure provides a state estimation method for a multi-stage voltage sag to analyze characteristics of multi-stage voltage sags due to different causes, consider the sudden change time of each stage and fault clearing time of a relay protection device in the monitoring data of the multi-stage voltage sag, infer a specific event that may occur in each stage of the multi-stage voltage sag, and finally obtain a multi-stage state estimation method that meets project requirements. The technical solutions are as follows:

A state estimation method for a multi-stage voltage sag includes the following steps:

• Step 1: analyzing the causes of a multi-stage voltage sag, determining characteristics of multi-stage voltage sags due to different causes, and identifying causes of voltage sags at different sudden change points accordingly;
• Step 2: building an action matrix of a relay protection device based on specified fault clearing time at different positions of a power grid;
• Step 3: obtaining one or more possible faulty lines by solving to form a faulty line set to reduce the calculation amount of state estimation of the multi-stage voltage sag;
• Step 4: calculating system impedance matrices before and after a fault, an impedance matrix in the case of a relay protection action on one side of the faulty line, and an impedance matrix when a DG device is disconnected from the power grid;
• Step 5: using an optimization model to infer a specific faulty line from the faulty line set obtained in step 3 and inferring a specific event of each stage of the multi-stage voltage sag based on a voltage sag type; and
• Step 6: estimating the duration of each sag stage based on the voltage sag type and the time length between adjacent sudden voltage sag change points and estimating a sag amplitude based on an event inference result in step 5.

Further, in step 1, the causes and corresponding characteristics of the multi-stage voltage sag are as follows:

Cause I: During a short circuit fault, one or more DG devices are in a power system trip, resulting in a loss of a power supply of the power grid. The corresponding characteristic is that a voltage amplitude on the right side of a corresponding sudden voltage sag change point when the voltage sag occurs is smaller than that on the left side.

Cause II: During the short circuit fault, relay protection devices on two sides of the faulty line trip at different times, resulting in a change in the topology of the power grid. The corresponding characteristic is that the voltage amplitude on the right side of the corresponding sudden voltage sag change point when the voltage sag occurs is greater than that on the left side.

Further, the identifying causes of voltage sags at different sudden change points specifically comprises:

Assuming that for a waveform segment containing the voltage sag, a total of s+1 sudden voltage sag change points are detected by using a waveform point detection method, which is respectively expressed as MP₀, MP₁,...,MP_x,...,MP_s and correspond to s + 1 sudden change time points t₀, t₁,..., t_x,..., ts, these sudden voltage sag change points divide a waveform of the voltage sag into ^S segments. Voltages of the different waveform segments are respectively expressed as u₁, u₂,...,u_x,...,u_s, and an ^mth voltage u_m is calculated as follows:

$u_m={\displaystyle {\sum }_{j=ff*\left(t_{m-1}+1/f\right)}^{j=ff*t_m}{x}_j/\left(\left(t_m-t_{m-1}+1/f\right)*ff\right);j\in \left[1,k\right],}$

Where ff represents a sampling rate in units of piece/second, f represents a power-frequency current frequency in units of Hz, x_j represents a j ^th value in an effective value sequence, m ∈ [1, s], and k represents a total quantity of sampling points.

A binary vector W is formed:

$\left\{\begin{array}{l}W=\left[w_1,w_2,\text{L}\text{,}{w}_x,\text{L},w_s\right],x\in \left[1,s\right]\hfill \\ w_x=\left\{\begin{array}{l}0,\text{if}{u}_x<u_{x-1}\text{or}x=1\hfill \\ 1,\text{if}{u}_x>u_{x-1}\hfill \end{array}\right)\hfill \end{array}\right).$

The value of the element W_x in the vector W is 0 or 1. When the value of the W_x is 0, it indicates that the voltage on a right side of MP_x-1 corresponding to a time point t_x-1 is less than that on the left side, which means that the multi-stage voltage sag is formed due to the cause I at this time point, where x > 1. On the contrary, when the value of the w_x is 1, it means that the multi-stage voltage sag is formed due to the cause II at the time point t_x-1, where x > 1.

Further, the building of an action matrix of a relay protection device in step 2 specifically includes:

Step 2.1: building a basic action matrix, in other words, a universal matrix describing the fault removal behavior of the relay protection device is built:

$\text{Γ=}\left(\begin{array}{ccccc}0& {\gamma }_{12}& {\gamma }_{13}& \text{L}& {\gamma }_{1n}\\ {\gamma }_{21}& 0& {\gamma }_{23}& \text{L}& {\gamma }_{2n}\\ {\gamma }_{31}& {\gamma }_{32}& 0& \text{O}& \text{M}\\ \text{M}& \text{M}& \text{O}& 0& {\gamma }_{\left(n-1\right)n}\\ {\gamma }_{n1}& {\gamma }_{n2}& \text{L}& {\gamma }_{n\left(n-1\right)}& 0\end{array}\right),$

Where n represents the quantity of buses in the system; ^γ_ij and ^γ_ji respectively represent the fault clearing time of protection devices close to i -side and j -side buses on a line l_ij in which i, j ∈ [1, n]; i ≠ j ; ^γ_ij, ^γ_ji = 0 represents that there is no physical connection between the i -side bus and the j -side bus.

Step 2.2: improving the universal matrix and decoupling an improved universal matrix into two matrices:

$\text{Γ}\widetilde{\text{Ν}}\text{=}\left(\begin{array}{ccccc}0& 0& 0& \text{0}& 0\\ {\gamma }_{21}& 0& 0& \text{0}& 0\\ {\gamma }_{31}& {\gamma }_{32}& 0& \text{0}& \text{0}\\ \text{M}& \text{M}& \text{O}& 0& 0\\ {\gamma }_{n1}& {\gamma }_{n2}& \text{L}& {\gamma }_{n\left(n-1\right)}& 0\end{array}\right)\text{and}$

$\text{ΓΔ=}\left(\begin{array}{ccccc}0& {\gamma }_{12}& {\gamma }_{13}& \text{L}& {\gamma }_{n1}\\ 0& 0& {\gamma }_{23}& \text{L}& {\gamma }_{n2}\\ 0& 0& 0& \text{O}& \text{M}\\ \text{0}& \text{0}& \text{0}& 0& {\gamma }_{n\left(n-1\right)}\\ 0& 0& \text{0}& 0& 0\end{array}\right).$

The upper and lower triangular matrices respectively represent fault clearing time of protective devices of the same type on two terminal buses of the line. The lower triangular element in the lower triangular matrix represents the fault clearing time of the protection device of the i -side bus of the line l_ij, and ^γ_ij = 0 and i > j represent that there is no physical connection between the i -side bus and the j -side bus. The upper triangular element in the upper triangular matrix represents a parameter on the other side of the line.

To represent two-stage protection, two other similar matrices are constructed to represent the parameters of backup protection:

$\text{Λ}\widetilde{\text{Ν}}\text{=}\left(\begin{array}{ccccc}0& 0& 0& \text{0}& 0\\ {\lambda }_{21}& 0& 0& \text{0}& 0\\ {\lambda }_{31}& {\lambda }_{32}& 0& \text{0}& \text{0}\\ \text{M}& \text{M}& \text{O}& 0& 0\\ {\lambda }_{n1}& {\lambda }_{n2}& \text{L}& {\lambda }_{n\left(n-1\right)}& 0\end{array}\right)\text{and}$

$\text{ΛΔ=}\left(\begin{array}{ccccc}0& {\lambda }_{12}& {\lambda }_{13}& \text{L}& {\lambda }_{1n}\\ 0& 0& {\lambda }_{23}& \text{L}& {\lambda }_{2n}\\ 0& 0& 0& \text{O}& \text{M}\\ \text{0}& \text{0}& \text{0}& 0& {\lambda }_{n\left(n-1\right)}\\ 0& 0& \text{0}& 0& 0\end{array}\right).$

Step 2.3: determining a cooperation relationship of the protection devices:

For the two-stage protection, there are four cooperation relationships between the main protection and the backup protection in the line, and the improved protection action matrix is expressed as the following cooperation relationships:

$\text{Γ}\tilde{\text{N}}+\text{ΓΔ=}{\text{Θ}}_{\text{I}},$

$\text{Γ}\tilde{\text{N}}+\text{ΛΔ=}{\text{Θ}}_{\text{II}},$

$\text{Λ}\tilde{\text{N}}+\text{ΛΔ=}{\text{Θ}}_{\text{III}},\text{and}$

$\text{Λ}\tilde{\text{N}}+\text{ΓΔ=}{\text{Θ}}_{\text{IV}}.$

Where Θ_I represents that the main protection on the two sides of the faulty line cooperates to remove the fault; Θ_II represents that the backup protection on the two sides of the faulty line cooperates to remove the fault; Θ_III and Θ_IV represent that the main protection on one side of the faulty line cooperates with the backup protection on the other side of the faulty line to remove the fault.

Step 2.4: correcting fault clearing time in the action matrix:

The fault clearing time is corrected as follows:

$\left\{\begin{array}{c}{\gamma }_{ij}={\gamma }_{ij,set}+{\delta }_{1,ij}\\ {\lambda }_{ij}={\lambda }_{ij,set}+{\delta }_{2,ij}\end{array};i,j\in \left[1,n\right];i\ne j.\right)$

Where ^γ_ij, _set and ^λ_ij,_set respectively represent specified values of the main protection and the backup protection; ^δ₁,_ij and ^δ₂,_ij respectively represent deviations between actual fault clearing time and the specified values; and the deviations are random numbers within [0, δ], where δ represents a maximum error value during the testing or historical operation of relay protection devices of the same model.

Further, step 3 specifically includes steps 3.1 and 3.2.

Step 3.1: defining the following four voltage sag types based on the sag cause identification in step 1:

• Type I: single-stage rectangular sag,
• Type II: multi-stage voltage sag due to cause I,
• Type III: multi-stage voltage sag due to cause II, and
• Type IV: multi-stage voltage sag due to cause I and cause II.

Step 3.2: obtaining the faulty line set for different types of sags.

For type I and type II:

Assuming that time points at which first and last sudden voltage sag change points of the waveform of the voltage sag are detected in step 1 are t₀ and t_s respectively, a time length from the occurrence of the fault to removal of the fault in the system is t_s -t₀. If two elements of symmetric positions of a main diagonal line of a matrix in the four matrices in step 2.3 are equal to t_s -t₀ within an error threshold, lines corresponding to these two elements are possible faulty lines, and therefore, a solution model of the faulty line set LF is as follows:

$\left\{\begin{array}{l}LF=\left\{\text{L},l_{ij},\text{L}\right\}\\ \text{s}.\text{t}.LF\subseteq LN\\ \sqrt{\left|{\theta }_{ij}-\left(t_s-t_0\right)\right|+\left|{\theta }_{ij}-\left(t_s-t_0\right)\right|}\le \\ \sqrt{\max \left({\delta }_{1,ij},{\delta }_{1,ji}\right)+\max \left({\delta }_{2,ij},{\delta }_{2,ji}\right)}\le \sqrt{2\delta }\\ i,j\in \left[1,n\right];i\ne j\end{array}\right),$

Where ^θ_ij represents an element in an i ^th row and a j ^th column in the matrices Θ_I to Θ_IV, and LN represents a set of lines in an intersection of sag domains of a bus of a monitoring device.

For type III and type IV:

Assuming that the time point at which the first sudden voltage sag change point of the waveform of the voltage sag is detected in step 1 is t₀ and time points of sudden voltage sag change points of two corresponding relay protection actions are t_x-1 and t_s respectively, a time length from the occurrence of the fault to the protection actions of the relay protection devices on the two sides of the faulty line are t_x-1-t₀ and t_s -t₀ respectively. In this case, the solution model of the faulty line set is as follows:

$\left\{\begin{array}{l}LF=\left\{\text{L},l_{ij},\text{L}\right\}\\ \text{s}.\text{t}.LF\subseteq LN\\ \sqrt{\left|{\theta }_{ij}-\left(t_{x-1}-t_0\right)\right|+\left|{\theta }_{ij}-\left(t_s-t_0\right)\right|}\le \\ \sqrt{\max \left({\delta }_{1,ij},{\delta }_{1,ji}\right)+\max \left({\delta }_{2,ij},{\delta }_{2,ji}\right)}\le \sqrt{2\delta }\\ i,j\in \left[1,n\right];i\ne j;s\ge 2;w_x=1\end{array}\right).$

Further, step 4 specifically includes steps 4.1, 4.2, 4.3, and 4.4.

Step 4.1: calculating the system impedance matrix before the fault:

A system admittance _matri_x Y^SE is expressed as a sum of a line admittance matrix

$Y_L^{SE}$

and a generator admittance matrix

$Y_G^{SE}$

, as shown below:

$Y^{SE}=Y_G^{SE}+Y_L^{SE}.$

Assuming there are ⁿ buses in the power system, the line admittance matrix

$Y_L^{SE}$

is calculated as follows based on a line topology relationship and an impedance parameter:

$Y_L^{SE}=\left(\begin{array}{cccc}{\alpha }_{11}^{se}& {\alpha }_{12}^{se}& \text{L}& {\alpha }_{1n}^{se}\\ {\alpha }_{21}^{se}& {\alpha }_{22}^{se}& \text{L}& {\alpha }_{2n}^{se}\\ \text{M}& \text{M}& \text{O}& \text{M}\\ {\alpha }_{n1}^{se}& {\alpha }_{n2}^{se}& \text{L}& {\alpha }_{nn}^{se}\end{array}\right),$

Where se = 1,2,0 represents a positive sequence, a negative sequence, and a zero sequence, ^α_ij represents mutual admittance of nodes i and j, ^α_ii represents self-admittance of the node i, and i ≠ j; i, j ∈ [1, n] .

The matrix

$Y_G^{SE}$

is a diagonal matrix, and an element value on a diagonal line is equal to the self-admittance of a corresponding generator, as shown below:

$Y_G^{SE}=\left(\begin{array}{cccc}{\beta }_{11}^{se}& 0& \text{L}& 0\\ 0& {\beta }_{22}^{se}& \text{L}& 0\\ \text{M}& \text{M}& \text{O}& \text{M}\\ 0& 0& \text{L}& {\beta }_{nn}^{se}\end{array}\right).$

Where β_ii = 0 represents that there is no generator for the bus.

The system impedance matrix is calculated as follows:

$Z_n^{SE}={\left(Y_G^{SE}\right)}^{-1}=\left(\begin{array}{cccc}{Z}_{11}^{se}& Z_{12}^{se}& \text{L}& Z_{1n}^{se}\\ Z_{21}^{se}& Z_{22}^{se}& \text{L}& Z_{2n}^{se}\\ \text{M}& \text{M}& \text{O}& \text{M}\\ Z_{n1}^{se}& Z_{n2}^{se}& \text{L}& Z_{nn}^{se}\end{array}\right).$

Step 4.2: calculating the system impedance matrix after the short circuit fault:

Mutual impedance

$Z_{mf}^{se}$

between a fault position f_l and a to-be-solved node ^m after the circuit short fault, and self-impedance

$Z_{ff}^{se}$

of the fault position f_l are respectively obtained by solving the following two formulas:

$Z_{mf}^{se}=\left(1-p\right)Z_{mi}^{se}+p{Z}_{mj}^{se}=g_{mf1}\left(i,j,p\right),\text{and}$

$Z_{ff}^{se}=\left\{\begin{array}{c}{\left(1-p\right)}^2{Z}_{ii}^{se}+p^2{Z}_{jj}^{se}+\\ 2p\left(1-p\right)Z_{ij}^{se}+p\left(1-p\right)z_{ij}^{se}\end{array}\right\}=g_{ff1}\left(i,j,p\right).$

Where

$Z_{mi}^{se},Z_{mj}^{se},Z_{ij}^{se},Z_{ii}^{se},\text{and}{Z}_{jj}^{se}$

are elements in the matrix

$z_{ij}^{se}$

represents impedance of the line ^l_ij, where the impedance is further expressed as functions g_mf1 (i, j, p) and g _ff1 (i, j, p). The two functions g_mf1 (i, j, p) and g _ff1 (i, j, p) respectively represent mutual impedance between the fault position and the target node ^m and self-impedance of the fault position when the short circuit fault occurs at the position f_l having a distance P away from the i -side of the line ^l_ij .

Step 4.3: calculating the impedance matrix in the case of the relay protection action on one side of the faulty line:

Assuming that after the short circuit fault occurs on the line l_ij, a j -side protection device acts to cut off the line on the corresponding side, and the system impedance is calculated as follows:

First, a branch with impedance of -z_ij is added between the i -side bus and the j -side bus in the original system. In this case, the system impedance matrix

$Z_{AT}^{SE}$

is corrected according to the following formula:

$\begin{array}{c}{Z}_{AT}^{SE}=Z_N^{SE}-\frac{\text{Δ}Z\text{Δ}{Z}^T}{-z_{ij}^{se}+Z_{ii}^{se}+Z_{jj}^{se}-2Z_{ij}^{se}}\\ =\left(\begin{array}{cccc}{Z}_{11,\text{AT}}^{se}& Z_{12,\text{AT}}^{se}& \text{L}& Z_{1n,\text{AT}}^{se}\\ Z_{21,\text{AT}}^{se}& Z_{22,\text{AT}}^{se}& \text{L}& Z_{2n,\text{AT}}^{se}\\ \text{M}& \text{M}& \text{O}& \text{M}\\ Z_{n1,\text{AT}}^{se}& Z_{n2,\text{AT}}^{se}& \text{L}& Z_{nn,\text{AT}}^{se}\end{array}\right).\end{array}$

Where ΔZ represents a process quantity, which is calculated according to the following formula:

$\text{Δ}Z={\left[\begin{array}{cccc}{Z}_{1j}^{se}-Z_{1i}^{se}& Z_{2j}^{se}-Z_{2i}^{se}& \cdots & Z_{nj}^{se}-Z_{ni}^{se}\end{array}\right]}^{\text{T}}.$

Subsequently, a branch with impedance of pz_ij is added to the i -side bus. In this case, the system impedance matrix is further corrected as

$Z_{ARP}^{SE}$

according to the following formula:

$Z_{ARP}^{SE}=\left(\begin{array}{ccccc}{Z}_{11,AT}^{se}& Z_{11,AT}^{se}& L& Z_{11,AT}^{se}& Z_{11,AT}^{se}\\ Z_{11,AT}^{se}& Z_{11,AT}^{se}& L& Z_{11,AT}^{se}& Z_{11,AT}^{se}\\ M& M& O& M& M\\ Z_{11,AT}^{se}& Z_{11,AT}^{se}& L& Z_{11,AT}^{se}& Z_{11,AT}^{se}\\ Z_{11,AT}^{se}& Z_{11,AT}^{se}& L& Z_{11,AT}^{se}& Z_{11,AT}^{se}+p{z}_{ij}\end{array}\right)$

Compared with the

$Z_{AT}^{SE},$

the

$Z_{ARP}^{SE}$

adds one line and one column to represent the mutual impedance

$Z_{mf}^{se}$

between the fault position f_l and each target bus ^m or the self-impedance

$Z_{ff}^{se}$

of the fault position. The impedance is further expressed as functions g_{mf 2} (i, j, p, d) and g_{ff 2} (i, j, p, d), where the two functions respectively represent mutual impedance between the fault position and the target node ^m and self-impedance of the fault position when a d -side protection device on the line acts to cut off a part of the line after the short circuit fault occurs at the position f_l far from the terminal P of the node i on the line l_ij, as shown in the following formulas:

$Z_{mf}^{se}=Z_{im,\text{AT}}^{se}=g_{mf2}\left(i,j,p,d\right);d=i\text{or}j\text{and}$

$Z_{ff}^{se}=Z_{ii,\text{AT}}^{se}+p{z}_{ij}=g_{ff2}\left(i,j,p,d\right);d=i\text{or}j.$

Step 4.4: calculating the impedance matrix when the DG device is disconnected from the power grid:

Assuming that DG devices of all buses in a bus set h are disconnected from the power grid, the diagonal matrix first needs to be corrected as follows:

$Y_G^{SE}=\left(\begin{array}{cccc}{\beta }_{11}^{se}& 0& \text{L}& 0\\ 0& {\beta }_{22}^{se}& \text{L}& 0\\ \text{M}& \text{M}& \text{O}& \text{M}\\ 0& 0& \text{L}& {\beta }_{nn}^{se}\end{array}\right);{\beta }_{ii}^{se}=0;\forall i\in h.$

The system impedance matrix before the fault is calculated according to the above formula. Finally, system impedance when the DG device is disconnected from the power grid in the case of the short circuit fault is calculated. The impedance is further expressed as functions g_{mf 3} (i, j, p, h) and g _ff3 (i, j, p, h), where the two functions respectively represent mutual impedance between the fault position and the target node ^m and self-impedance of the fault position after the short circuit fault occurs at the position f_l far from the terminal P of the node i on the line l_ij, and the DG devices of all the buses in the bus set h are disconnected from the power grid, as shown in the following two formulas:

$Z_{mf}^{se}=g_{mf3}\left(i,j,p,h\right)\text{and}$

$Z_{ff}^{se}=g_{ff3}\left(i,j,p,h\right).$

Further, in step 5, the following four optimization models are available based on the sag types:

(1) For type I, the following optimization model is used to infer the faulty line and its short circuit condition using the variables i, j, and P :

$\left\{\begin{array}{l}\max \text{-}\sqrt{\left|f\left(g_{mf1}\left(i,j,p\right),g_{ff1}\left(i,j,p\right)\right)\right|-u_1}\hfill \\ \text{s}\text{.t}\text{.}i,j\in \left[1,n\right];i\ne j\hfill \end{array}\right),$

where f (·) represents a short-circuit calculation function.

(2) For type II, the following optimization model is used to infer the faulty line and its short circuit condition and a time point and a sequence of disconnecting the DG device from the power grid using the variables i, j, p, and ^h_q:

$\left\{\begin{array}{l}\max \text{-}\left(\begin{array}{c}\sqrt{\left|f\left(g_{mf1}\left(i,j,p\right),g_{ff1}\left(i,j,p\right)\right)\right|-u_1}+\\ {\displaystyle \sum _{x=2}^s{\displaystyle \sum _{q=1}^{s-1}\sqrt{\left|f\left(g_{mf3}\left(i,j,p,h_q\right),g_{ff3}\left(i,j,p,h_q\right)\right)\right|-u_x}}}\end{array}\right)\hfill \\ \text{s}\text{.t}\text{.}i,j\in \left[1,n\right];i\ne j;h_1\subseteq h_2\subseteq \cdots \subseteq h_{s\text{-1}}\hfill \end{array}\right),$

where ^s represents a quantity of stages of the multi-stage voltage sag and ^h_q represents a set of DG devices disconnected from the power grid during a q-stage voltage sag.

(3) For type III, the following optimization model is used to infer the faulty line and its short circuit condition and a tripping time point and sequence of the relay protection device using the variables i, j, p, and d :

$\left\{\begin{array}{l}\max \text{-}\left(\begin{array}{c}\sqrt{\left|f\left(g_{mf1}\left(i,j,p\right),g_{ff1}\left(i,j,p\right)\right)\right|-u_1}+\\ \sqrt{\left|f\left(g_{mf2}\left(i,j,p,d\right),g_{ff2}\left(i,j,p,d\right)\right)\right|-u_2}\end{array}\right)\hfill \\ \text{s}\text{.t}\text{.}i,j\in \left[1,n\right];i\ne j\hfill \end{array}\right).$

(4) For type IV, the following optimization model is used to infer the faulty line and its short circuit condition, the tripping time point and sequence of the relay protection device, and the time point and the sequence of disconnecting the DG device from the power grid using the variables i, j, p, d, and ^h_q :

$\left\{\begin{array}{l}\max \text{-}\left(\begin{array}{c}\begin{array}{c}\sqrt{\left|f\left(g_{mf1}\left(i,j,p\right),g_{ff1}\left(i,j,p\right)\right)\right|-u_1}+\\ {\displaystyle \sum _{x=2}^o{\displaystyle \sum _{q=1}^{o-1}\sqrt{\left|f\left(g_{mf3}\left(i,j,p,h_q\right),g_{ff3}\left(i,j,p,h_q\right)\right)\right|-u_x}}}\end{array}\\ \begin{array}{c}+\left|f\left(g_{mf2}\left(i,j,p,d\right),g_{ff2}\left(i,j,p,d\right)\right)\right|-u_{o+1}+\\ {\displaystyle \sum _{x=o+2}^s{\displaystyle \sum _{h=o}^{s-2}\sqrt{\left|f\left(g_{mf3}\left(i,j,p,h_q\right),g_{ff3}\left(i,j,p,h_q\right)\right)\right|-u_x}}}\end{array}\end{array}\right)\hfill \\ \text{s}\text{.t}\text{.}i,j\in \left[1,n\right];i\ne j;h_1\subseteq h_2\subseteq \cdots \subseteq h_{s\text{-2}}\hfill \end{array}\right).$

Where o represents a time point at which a relay protection device on the side of the faulty line first acts to trip and ^u₀₊₁ represents an amplitude of an o+1-stage voltage sag.

Further, in step 6, the following four voltage sag estimation methods are available based on the sag types:

1) For type I, a sag amplitude of any unmonitored bus ^m is estimated according to the following formula:

$f\left(g_{mf1}\left(i,j,p\right),g_{ff1}\left(i,j,p\right)\right).$

2) For type II, a sag amplitude of any unmonitored bus ^m is estimated according to the following formula:

$\left\{\begin{array}{l}{u}_1=f\left(g_{mf1}\left(i,j,p\right),g_{ff1}\left(i,j,p\right)\right)\hfill \\ u_x=f\left(g_{mf3}\left(i,j,p,h_{x-1}\right),g_{ff3}\left(i,j,p,h_{x-1}\right)\right)\hfill \end{array}\right).$

3) For type III, a sag amplitude of any unmonitored bus ^m is estimated according to the following formula:

$\left\{\begin{array}{l}{u}_1=f\left(g_{mf1}\left(i,j,p\right),g_{ff1}\left(i,j,p\right)\right)\hfill \\ u_2=f\left(g_{mf2}\left(i,j,p,d\right),g_{ff2}\left(i,j,p,d\right)\right)\hfill \end{array}\right).$

4) For type IV, a sag amplitude of any unmonitored bus ^m is estimated according to the following formula:

$\left\{\begin{array}{l}{u}_1=f\left(g_{mf1}\left(i,j,p\right),g_{ff1}\left(i,j,p\right)\right)\hfill \\ u_x=f\left(g_{mf3}\left(i,j,p,h_{x-1}\right),g_{ff3}\left(i,j,p,h_{x-1}\right)\right)\hfill \\ u_y=f\left(g_{mf2}\left(i,j,p,d\right),g_{ff2}\left(i,j,p,d\right)\right)\hfill \\ u_z=f\left(g_{mf3}\left(i,j,p,h_{z-1}\right),g_{ff3}\left(i,j,p,h_{z-1}\right)\right)\hfill \\ x\in \left[2,y-1\right];z\in \left[y+1,s-2\right]\hfill \end{array}\right),$

Where ^y and z respectively represent y^th and z^th stages of the multi-stage voltage sag.

The present disclosure has the following beneficial effects:

1) The present disclosure analyzes characteristics of multi-stage voltage sags due to different causes and proposes a method for determining a cause of sudden change time of an amplitude of the multi-stage voltage sag based on a jump characteristic of the sag amplitude.

2) The present disclosure proposes calculation methods of a relay protection action matrix and a faulty line set to preliminarily obtain a faulty line set based on fault clearing time of a relay protection device and other characteristics to effectively reduce the calculation amount of sag state estimation.

3) Based on a substantive characteristic that different events cause a change of system impedance, the present disclosure proposes a method for inferring a cause of an event in each stage of the multi-stage voltage sag to improve the accuracy of the state estimation.

4) Based on a cause inference result of the event, the present disclosure proposes a state estimation method for the multi-stage voltage sag to resolve the problem that it is difficult to apply an existing state estimation method to the multi-stage voltage sag.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 shows a short-circuit calculation model of a power system;

FIG. 2 shows a short-circuit calculation model of a power system after a faulty line is cut off on one side according to the present disclosure; and

FIG. 3 is a process of building a state estimation model of a multi-stage voltage sag according to the present disclosure.

DETAILED DESCRIPTION OF THE EMBODIMENTS

The present disclosure will be further described below in conjunction with the accompanying drawings and specific embodiments.

The technical solution of the present disclosure is mainly divided into six steps: identifying causes of a multi-stage voltage sag, building an action matrix of a relay protection device, obtaining a faulty line set by solving equations, obtaining an impedance matrix by solving equations, and inferring a specific event at each stage of the multi-stage voltage sag and a state estimation method for the multi-stage voltage sag. Details of each step and its substeps are as follows:

1. The causes of the multi-stage voltage sag are identified.

There Are Two Causes of the Multi-Stage Voltage Sag

Cause 1: During a short circuit fault, one or more DG device are in a power system trip, resulting in a loss of a power supply from a power grid.

Cause 2: During the short circuit fault, relay protection devices on two sides of a faulty line trip at different times, resulting in a change of the topology of the power grid.

Characteristics of Multi-Stage Voltage Sags Due to Different Causes

Cause 1 means that after the fault, the system further loses the power supply, so a monitored voltage sag amplitude is further reduced. The corresponding characteristic is that a voltage amplitude on a right side of the corresponding sudden voltage sag change point when cause 1 occurs is smaller than that on a left side.

On the contrary, cause 2 means that after the fault occurs, a protection device on one side of the line partially isolates the fault from the system, so the monitored voltage sag amplitude increases. The corresponding characteristic is that a voltage amplitude on the right side of the corresponding sudden voltage sag change point when the cause 2 occurs is greater than that on the left side.

Cause Identification of Voltage Sags at Different Sudden Change Points

It is assumed that for a waveform segment containing the voltage sag, a total of s+1 sudden voltage sag change points are detected by using a waveform point detection method, which are respectively expressed as MP₀,MP₁,...,MP_x,...,MP_s and correspond to s +1 sudden change time points t₀, t₁,..., t_x,..., t_s . These sudden voltage sag change points divide a waveform of the voltage sag into ^S segments. Voltages of different waveform segments are respectively expressed as ^u₁, ^u₂, .. .,^u_x,.. ., u_s, and an m (m ∈ [1, s]) ^th voltage u_m can be calculated as follows:

$\begin{array}{cc}{u}_m={\displaystyle {\sum }_{j=ff*\left(t_{m-1}+1/f\right)}^{j=ff*t_m}{x}_j}/\left(\left(t_m-t_{m-1}+1/f\right)*ff\right);j\in \left[1,k\right]& \text{­­­(7)}\end{array}$

In the above formula, ff represents a sampling rate in units of piece/second, f represents a power-frequency current frequency in units of Hz, and ^x_j represents a j ^th value in an effective value sequence. Binary vector W shown in formula (8) can be formed according to the formula (7):

$\begin{array}{cc}\left\{\begin{array}{l}W=\left[w_1,w_2,\text{L}\text{,}{w}_x,\text{L,}{w}_s\right],x\in \left[1,s\right]\hfill \\ w_x=\left\{\begin{array}{l}0,\text{if}{u}_x<u_{x-1}\text{or}x=1\hfill \\ 1,\text{if}{u}_x>u_{x-1}\hfill \end{array}\right)\hfill \end{array}\right)& \text{­­­(8)}\end{array}$

The value of the element w_x in the vector W is 0 or 1. When the value is 0, it indicates that a voltage on a right side of MP_x-1 corresponding to time point t_x-1 (x > 1) is less than that on a left side, which means that the multi-stage voltage sag is formed due to cause 1 at this time point. On the contrary, when the value is 1, it means that the multi-stage voltage sag is formed due to cause 2 at the time point t_x-1 (x > 1).

2. Action matrix of the relay protection device.

Actual tripping time of the relay protection device in a sag event can be roughly obtained based on multi-stage voltage sag causes at time points corresponding to different sudden voltage sag change points. In order to use the information, a method for building the action matrix of the relay protection device based on specified fault clearing time at different positions of the power grid is proposed in step 2.

Basic Action Matrix

There are different types of relay protection, including low-voltage protection, over-current protection, and distance protection. Each type of protection has its own trigger rule and fault clearing time to clear the fault. The present disclosure provides a universal matrix for describing a fault removal behavior of the relay protection device, as shown in the following formula:

$\begin{array}{cc}\text{Γ=}\left(\begin{array}{ccccc}0& {\gamma }_{12}& {\gamma }_{13}& \text{L}& {\gamma }_{1n}\\ {\gamma }_{21}& 0& {\gamma }_{23}& \text{L}& {\gamma }_{2n}\\ {\gamma }_{31}& {\gamma }_{32}& 0& \text{O}& \text{M}\\ \text{M}& \text{M}& \text{O}& 0& {\gamma }_{\left(n-1\right)n}\\ {\gamma }_{n1}& {\gamma }_{n2}& \text{L}& {\gamma }_{n\left(n-1\right)}& 0\end{array}\right)& \text{­­­(9)}\end{array}$

Where n represents a quantity of buses in the system, ^γ_ij and ^γ_ji (i, j ∈ [1, n]; i ≠ j) respectively represent fault clearing time of protection devices close to i -side and j -side buses on line l_ij in units of ms, and ^γ_ij, ^γ_ji= 0 represents that there is no physical connection between the i -side bus and the j -side bus. The fault clearing time is a specified value of the relay protection device recorded in the software of a power distribution automation system and other power enterprises.

The Action Matrix is Improved

In practical application, staged protection is widely used in a power distribution system, such as two-stage protection. The two-stage protection means that two types of protection are provided on one line, namely, main protection and backup protection. The two types of protection have different fault clearing time, and the main protection has shorter fault clearing time. When the short circuit fault occurs, the main protection is triggered first, and the backup protection works when the main protection cannot be triggered. In general, the fault can be removed by the main protection independently, and sometimes may be removed by the backup protection independently or by both the main protection and the backup protection. Considering that it is difficult for the basic action matrix defined in the formula (9) to present a cooperation relationship of the staged protection, the present disclosure improves the basic action matrix and decouples an improved basic action matrix into two matrices, as shown in equations (10) and (11):

$\begin{array}{cc}\text{Γ}\tilde{\text{N}}\text{=}\left(\begin{array}{ccccc}0& 0& 0& 0& 0\\ {\gamma }_{21}& 0& 0& 0& 0\\ {\gamma }_{31}& {\gamma }_{32}& 0& 0& 0\\ \text{M}& \text{M}& \text{O}& 0& 0\\ {\gamma }_{n1}& {\gamma }_{n2}& \text{L}& {\gamma }_{n\left(n-1\right)}& 0\end{array}\right)& \text{­­­(10)}\end{array}$

$\begin{array}{cc}\text{ΓΔ=}\left(\begin{array}{ccccc}0& {\gamma }_{12}& {\gamma }_{13}& \text{L}& {\gamma }_{1n}\\ 0& 0& {\gamma }_{23}& \text{L}& {\gamma }_{2n}\\ 0& 0& 0& \text{O}& \text{M}\\ 0& 0& 0& 0& {\gamma }_{\left(n-1\right)n}\\ 0& 0& 0& 0& 0\end{array}\right)& \text{­­­(11)}\end{array}$

In the above formulas, the upper and lower triangular matrices respectively represent fault clearing time of protection devices of the same type on two terminal buses of the line. The lower triangular element in the lower triangular matrix represents the fault clearing time of the protection device of the i -side bus of the line l_ij, and ^γ_ij = 0(i > j) represents that there is no physical connection between the i -side bus and the j -side bus. The upper triangular element in the upper triangular matrix represents a parameter on the other side of the line, which has a similar meaning. The formulas (10) and (11) represent parameters of the main protection. To represent the two-stage protection, two other similar matrices are constructed to represent parameters of the backup protection, as shown in formulas (12) and (13):

$\begin{array}{cc}\text{Λ}\tilde{\text{N}}=\left(\begin{array}{ccccc}0& 0& 0& 0& 0\\ {\lambda }_{21}& 0& 0& 0& 0\\ {\lambda }_{31}& {\lambda }_{32}& 0& 0& 0\\ \text{M}& \text{M}& \text{O}& 0& 0\\ {\lambda }_{n1}& {\lambda }_{n2}& \text{L}& {\lambda }_{n\left(n-1\right)}& 0\end{array}\right)& \text{­­­(12)}\end{array}$

$\begin{array}{cc}\text{ΛΔ}=\left(\begin{array}{ccccc}0& {\lambda }_{12}& {\lambda }_{13}& \text{L}& {\lambda }_{1n}\\ 0& 0& {\lambda }_{23}& \text{L}& {\lambda }_{2n}\\ 0& 0& 0& \text{O}& \text{M}\\ 0& 0& 0& 0& {\lambda }_{\left(n-1\right)n}\\ 0& 0& 0& 0& 0\end{array}\right)& \text{­­­(13)}\end{array}$

Cooperation Relationship of the Protection Devices

For the widely used two-stage protection, there are four cooperation relationships between the main protection and the backup protection in the line. The improved protection action matrix can conveniently express these cooperation relationships as follows:

$\begin{array}{cc}\text{Γ}\tilde{\text{N}}\text{+}\text{ΓΔ=}{\text{Θ}}_{\text{I}}& \text{­­­(14)}\end{array}$

$\begin{array}{cc}\text{Γ}\tilde{\text{N}}\text{+}\text{ΛΔ=}{\text{Θ}}_{\text{II}}& \text{­­­(15)}\end{array}$

$\begin{array}{cc}\text{Λ}\tilde{\text{N}}\text{+}\text{ΛΔ=}{\text{Θ}}_{\text{III}}& \text{­­­(16)}\end{array}$

$\begin{array}{cc}\text{Λ}\tilde{\text{N}}\text{+}\text{ΓΔ=}{\text{Θ}}_{\text{IV}}& \text{­­­(17)}\end{array}$

Relationship 1: As shown in the formula (14), Θ_I represents that the main protection on two sides of the faulty line cooperates with each other to remove the fault.

Relationship 2: As shown in the formula (15), Θ_II represents that the backup protection on the two sides of the faulty line cooperates with each other to remove the fault.

Relationship ¾: As shown in the formulas (16) and (17), Θ_III and Θ_IV represent that the main protection on one side of the faulty line cooperates with the backup protection on the other side of the faulty line to remove the fault.

Fault Clearing Time in the Action Matrix is Corrected

The fault clearing time used in the formulas (9) to (17) is a specified value obtained from a power company. Theoretically, the fault clearing time of the main protection or the backup protection on the two sides of the line should be the same and equal to the specified value. However, in actual operation, due to error, the elements in these matrices may deviate from specified values, so necessary corrections need to be made, as shown in formula (18):

$\begin{array}{cc}\left\{\begin{array}{l}{\gamma }_{ij}={\gamma }_{ij,set}+{\delta }_{1,ij}\hfill \\ {\lambda }_{ij}={\lambda }_{ij,set}+{\delta }_{2,ij}\hfill \end{array}\right);i,j\in \left[1,n\right];i\ne j& \text{­­­(18)}\end{array}$

In the above formula, ^γ_ij,_set and ^λ_ij,_set respectively represent specified values of the main protection and the backup protection. ^δ_1,ij and ^δ₂,_ij respectively represent deviations between actual fault clearing time and the specified values. The deviations are random numbers within [0, δ ], where δ represents a maximum error value during the testing or historical operation of relay protection devices of the same model.

3. The faulty line set is obtained by solving equations.

One or more possible faulty lines can be obtained by solving based on the sag cause identification in step 1 and the content of the protection action matrix in step 2 to reduce the calculation amount of the state estimation of the multi-stage voltage sag.

Sag Type

The following four voltage sag types are defined based on the sag cause identification in step 1:

• Type 1: single-stage rectangular sag,
• Type 2: multi-stage voltage sag due to the cause 1,
• Type 3: multi-stage voltage sag due to the cause 2, and
• Type 4: multi-stage voltage sag due to the cause 1 and the cause 2.

The Faulty Line Set is Obtained for Different Types of Sags by Solving Equations

Type 1 and 2: In these two types of voltage sag events, the relay protection devices on the two sides of the faulty line act simultaneously at a time point of the last sudden voltage sag change to remove the fault. Assuming that time points at which first and last sudden voltage sag change points of the waveform of the voltage sag are detected in step 1 are t₀ and t_s respectively, the time length from the occurrence of the fault to removal of the fault in the system is t_s-t₀ . If two elements of symmetric positions of the main diagonal line of a matrix in formulas (14) to (17) are equal to t_s-t₀ within an error threshold, lines corresponding to these two elements are possible faulty lines. Therefore, a solution model of the faulty line set LF is as follows:

$\begin{array}{cc}\left\{\begin{array}{l}LF=\left\{\text{L},l_{ij},\text{L}\right\}\\ \text{s}\text{.t}\text{.}LF\subseteq LN\\ \sqrt{\left|{\theta }_{ij}-\left(t_s-t_0\right)\right|+\left|{\theta }_{ji}-\left(t_s-t_0\right)\right|}\le \\ \sqrt{\max \left({\delta }_{1,ij},{\delta }_{1,ji}\right)+\max \left({\delta }_{2,ij},{\delta }_{2,ji}\right)}\le \sqrt{2\delta }\\ i,j\in \left[1,n\right];i\ne j\end{array}\right)& \text{­­­(19)}\end{array}$

Where ^θ_ij represents an element in an i ^th row and a j ^th column in the matrices Θ_I to Θ_IV and LN represents a set of lines in the intersection of sag domains of a bus of a monitoring device.

Type 3 and 4: In these two types of voltage sag events, the relay protection devices on the two sides of the faulty line cannot act simultaneously to remove the fault. It is assumed that the time point at which the first sudden voltage sag change point of the waveform of the voltage sag is detected in step 1 and time points of sudden voltage sag change points of two corresponding relay protection actions are t₀, t_x-1, and t_s respectively. The time length from the occurrence of the fault to the protection actions of the relay protection devices on the two sides of the faulty line are t_x-1-t₀ and t_s -t₀, respectively. In this case, the solution model of the faulty line set is as follows:

$\begin{array}{cc}\left\{\begin{array}{l}LF=\left\{\text{L},l_{ij},\text{L}\right\}\\ \text{s}\text{.t}\text{.}LF\subseteq LN\\ \sqrt{\left|{\theta }_{ij}-\left(t_s-t_0\right)\right|+\left|{\theta }_{ji}-\left(t_s-t_0\right)\right|}\le \\ \sqrt{\max \left({\delta }_{1,ij},{\delta }_{1,ji}\right)+\max \left({\delta }_{2,ij},{\delta }_{2,ji}\right)}\le \sqrt{2\delta }\\ i,j\in \left[1,n\right];i\ne j\end{array}\right)& \text{­­­(20)}\end{array}$

4. The impedance matrix is obtained by solving equations.

Whether caused by cause 1 or cause 2, the multi-stage voltage sag is essentially a change of a system impedance matrix due to an event during the fault.

System Impedance Matrix Before the Fault

System admittance matrix Y^SE may be expressed as the sum of line admittance matrix

$Y_L^{SE}$

and generator admittance matrix

$Y_G^{SE}$

, as shown below:

$\begin{array}{cc}{Y}^{SE}=Y_G^{SE}+Y_L^{SE}& \text{­­­(21)}\end{array}$

Assuming there are ⁿ buses in the power system, the line admittance matrix

$Y_L^{SE}$

can be calculated as follows based on a line topology relationship and an impedance parameter:

$\begin{array}{cc}{Y}_L^{SE}=\left(\begin{array}{cccc}{\alpha }_{11}^{se}& {\alpha }_{12}^{se}& \text{L}& {\alpha }_{1n}^{se}\\ {\alpha }_{21}^{se}& {\alpha }_{22}^{se}& \text{L}& {\alpha }_{2n}^{se}\\ \text{M}& \text{M}& \text{O}& \text{M}\\ {\alpha }_{n1}^{se}& {\alpha }_{n2}^{se}& \text{L}& {\alpha }_{nn}^{se}\end{array}\right)& \text{­­­(22)}\end{array}$

In the above formula, se=1,2,0 represents a positive sequence, a negative sequence, and a zero sequence, α _ij (i ≠ j; i, j ∈ [1,n]) represents mutual admittance of nodes i and j, and α_ii (i ∈ [1, n]) represents self-admittance of the node i.

The matrix

$Y_G^{SE}$

is a diagonal matrix, and an element value on a diagonal line is equal to self-admittance of a corresponding generator, as shown below:

$\begin{array}{cc}{Y}_G^{SE}=\left(\begin{array}{cccc}{\beta }_{11}^{se}& 0& \text{L}& 0\\ 0& {\beta }_{22}^{se}& \text{L}& 0\\ \text{M}& \text{M}& \text{O}& \text{M}\\ 0& 0& \text{L}& {\beta }_{nn}^{se}\end{array}\right)& \text{­­­(23)}\end{array}$

where β_ii = 0(i ∈ [1, n]) represents that there is no generator for the bus.

The system impedance matrix can be calculated as follows according to the formula (21):

$\begin{array}{cc}{Z}_N^{SE}={\left(Y_G^{SE}\right)}^{-1}=\left(\begin{array}{cccc}{Z}_{11}^{se}& Z_{12}^{se}& \text{L}& Z_{1n}^{se}\\ Z_{21}^{se}& Z_{22}^{se}& \text{L}& Z_{2n}^{se}\\ \text{M}& \text{M}& \text{O}& \text{M}\\ Z_{n1}^{se}& Z_{n2}^{se}& \text{L}& Z_{nn}^{se}\end{array}\right)& \text{­­­(24)}\end{array}$

System Impedance Matrix After the Short Circuit Fault

When the short circuit fault occurs in the system, system impedance can be calculated using a short-circuit calculation model of a complex power system shown in FIG. 1. Mutual impedance

$Z_{mf}^{se}$

between fault position f_l and to-be-solved node ^m after the circuit short fault and self-impedance

$Z_{ff}^{se}$

of the fault position f_l can be respectively obtained by solving the following two formulas:

$\begin{array}{cc}{Z}_{mf}^{se}=\left(1-p\right)Z_{mi}^{se}+p{Z}_{mj}^{se}=g_{mf1}\left(i,j,p\right)& \text{­­­(25)}\end{array}$

$\begin{array}{cc}{Z}_{ff}^{se}=\left\{\begin{array}{c}{\left(1-p\right)}^2{Z}_{ii}^{se}+p^2{Z}_{jj}^{se}+\\ 2p\left(1-p\right)Z_{ij}^{se}+p\left(1-p\right)z_{ij}^{se}\end{array}\right\}=g_{ff1}\left(i,j,p\right)& \text{­­­(26)}\end{array}$

Where

$Z_{mi}^{se},Z_{mj}^{se},Z_{ij}^{se},Z_{ii}^{se},\text{and}{Z}_{jj}^{se}$

are elements in the matrix

$Z_N^{SE};\text{and}{z}_{ij}^{se}$

represents impedance of the line ^l_ij, where the impedance is further expressed as functions g_mf1 (i, j,p) and g _ff1 (i, j,p), and the two functions respectively represent mutual impedance between the fault position and the target node ^m and self-impedance of the fault position when the short circuit fault occurs at the position ^f_l far from terminal P of the node i on the line ^l_ij .

Impedance Matrix in the Case of a Relay Protection Action on One Side of the Faulty Line

Assuming that after the short circuit fault occurs on the line l_ij, a j -side protection device acts to cut off the line on the corresponding side, and the system impedance can be calculated by using the model shown in FIG. 2.

First, a branch with impedance of -z_ij is added between the i -side bus and the j -side bus in the original system. In this case, the system impedance matrix

$Z_{AT}^{SE}$

can be corrected according to the following formula:

$\begin{array}{cc}\begin{array}{c}{Z}_{AT}^{SE}=Z_N^{SE}-\frac{\text{Δ}Z\text{Δ}{Z}^T}{-z_{ij}^{se}+Z_{ii}^{se}+Z_{jj}^{se}-2Z_{ij}^{se}}\\ =\left(\begin{array}{cccc}{Z}_{11,\text{AT}}^{se}& Z_{12,\text{AT}}^{se}& \text{L}& Z_{1n,\text{AT}}^{se}\\ Z_{21,\text{AT}}^{se}& Z_{22,\text{AT}}^{se}& \text{L}& Z_{2n,\text{AT}}^{se}\\ \text{M}& \text{M}& \text{O}& \text{M}\\ Z_{n1,\text{AT}}^{se}& Z_{n2,\text{AT}}^{se}& \text{L}& Z_{nn,\text{AT}}^{se}\end{array}\right)\end{array}& \text{­­­(27)}\end{array}$

Where ΔZ can be calculated according to the following formula:

$\begin{array}{cc}\text{Δ}Z={\left[\begin{array}{cccc}{Z}_{1j}^{se}-Z_{1i}^{se}& Z_{2j}^{se}-Z_{2i}^{se}& \cdots & Z_{nj}^{se}-Z_{ni}^{se}\end{array}\right]}^T& \text{­­­(28)}\end{array}$

Then, a branch with impedance of ^pz_ij is added to the i -side bus. In this case, the system impedance matrix is further corrected as

$Z_{ARP}^{SE}$

according to the following formula:

$\begin{array}{cc}{Z}_{ARP}^{SE}=\left(\begin{array}{ccccc}{Z}_{11,AT}^{se}& Z_{12,AT}^{se}& \text{L}& Z_{1n,AT}^{se}& Z_{1i,AT}^{se}\\ Z_{21,AT}^{se}& Z_{22,AT}^{se}& \text{L}& Z_{2n,AT}^{se}& Z_{2i,AT}^{se}\\ \text{M}& \text{M}& \text{O}& \text{M}& \text{M}\\ Z_{n1,AT}^{se}& Z_{n2,AT}^{se}& \text{L}& Z_{nn,AT}^{se}& Z_{ni,AT}^{se}\\ Z_{i1,AT}^{se}& Z_{i2,AT}^{se}& \text{L}& Z_{in,AT}^{se}& Z_{ii,AT}^{se}+p{z}_{ij}\end{array}\right)& \text{­­­(29)}\end{array}$

Compared with the

$Z_{AT}^{SE},Z_{ARP}^{SE}$

adds one line and one column to represent the mutual impedance

$Z_{mf}^{se}$

between the fault position f and each target bus ^m or the self-impedance

$Z_{ff}^{se}$

of the fault position. The impedance is further expressed as functions g _{mf 2} (i, j, p, d) and g_{ff 2} (i, j, p, d), where the two functions respectively represent mutual impedance between the fault position and the target node ^m and self-impedance of the fault position when a d -side protection device on the line acts to cut off a part of the line after the short circuit fault occurs at the position f far from the terminal P of the node i on the line l_ij, as shown in the following formulas:

$\begin{array}{cc}{Z}_{mf}^{se}=Z_{im,\text{AT}}^{se}=g_{mf2}\left(i,j,p,d\right);d=i\text{or}j& \text{­­­(30)}\end{array}$

and

$\begin{array}{cc}{Z}_{ff}^{se}=Z_{ii,\text{AT}}^{se}+p{z}_{ij}=g_{ff2}\left(i,j,p,d\right);d=i\text{or}j& \text{­­­(31)}\end{array}$

Impedance Matrix When the DG Device is Disconnected From the Power Grid

When the short circuit fault occurs, DG devices on some nodes may be disconnected from the power grid due to insufficient low-voltage ride through capabilities. Assuming that DG devices of all buses in bus set h are disconnected from the power grid, the matrix shown in the formula (23) first needs to be corrected as follows:

$\begin{array}{cc}{Y}_G^{SE}=\left(\begin{array}{cccc}{\beta }_{11}^{se}& 0& \text{L}& 0\\ 0& {\beta }_{22}^{se}& \text{L}& 0\\ \text{M}& \text{M}& \text{O}& \text{M}\\ 0& 0& \text{L}& {\beta }_{nn}^{se}\end{array}\right);{\beta }_{ii}^{se}=0;\forall i\in h& \text{­­­(32)}\end{array}$

Then, the system impedance matrix before the fault is calculated by substituting the above formula into the formula (24). Finally, system impedance when the DG device is disconnected from the power grid in the case of the short circuit fault is calculated according to the formulas (1) to (26). The impedance is further expressed as functions g _{mf 3} (i, j, p, h) and g _ff3 (i, j, p,h), where the two functions respectively represent mutual impedance between the fault position and the target node ^m and self-impedance of the fault position after the short circuit fault occurs at the position f far from the terminal P of the node i on the line l_ij, and the DG devices of all the buses in the bus set h are disconnected from the power grid, as shown in the following formulas:

$\begin{array}{cc}{Z}_{mf}^{se}=g_{mf3}\left(i,j,p,h\right)& \text{­­­(33)}\end{array}$

and

$\begin{array}{cc}{Z}_{ff}^{se}=g_{ff3}\left(i,j,p,h\right)& \text{­­­(34)}\end{array}$

5. The specific event at each stage of the multi-stage voltage sag is inferred.

After the system impedance matrices under different situations are obtained by solving equations in step 4, an optimization model in step 5 may be used to infer a specific faulty line from the faulty line set obtained in step 3, and the specific event of each stage of the multi-stage voltage sag is inferred. There are also the following four inference models based on the sag types:

Type 1

For this type of sag, the following optimization model may be used to infer the faulty line and its short circuit condition (in other words, i, j, and p ):

$\begin{array}{cc}\left\{\begin{array}{l}\max \text{-}\sqrt{\left|f\left(g_{mf1}\left(i,j,p\right),g_{ff1}\left(i,j,p\right)\right)\right|-u_1}\\ \text{s}\text{.t}\text{.}i,j\in \left[1,n\right];i\ne j\end{array}\right)& \text{­­­(35)}\end{array}$

Where f ( ) represents a short-circuit calculation function, and calculation formulas of different types of short circuit faults are shown in the formulas (2) to (6).

Type 2

For this type of sag, the following optimization model may be used to infer the faulty line and its short circuit condition, and a time point and a sequence of disconnecting the DG device from the power grid (in other words, i, j, p, and h_q)_:

$\begin{array}{cc}\left\{\begin{array}{l}\max \text{-}\left(\begin{array}{c}\sqrt{\left|f\left(g_{mf1}\left(i,j,p\right),g_{ff1}\left(i,j,p\right)\right)\right|-u_1}+\\ {\displaystyle \sum _{x=2}^s{\displaystyle \sum _{q=1}^{s-1}\sqrt{\left|f\left(g_{mf3}\left(i,j,p,h_q\right),g_{ff3}\left(i,j,p,h_q\right)\right)\right|-u_x}}}\end{array}\right)\\ \text{s}\text{.t}\text{.}i,j\in \left[1,n\right];i\ne j;h_1\subseteq h_2\subseteq \cdots \subseteq h_{s-1}\end{array}\right)& \text{­­­(36)}\end{array}$

Type 3

For this type of sag, the following optimization model may be used to infer the faulty line and its short circuit condition, and a tripping time point and sequence of the relay protection device (in other words, i, j, p, and d ):

$\begin{array}{cc}\left\{\begin{array}{l}\max \text{-}\left(\begin{array}{c}\sqrt{\left|f\left(g_{mf1}\left(i,j,p\right),g_{ff1}\left(i,j,p\right)\right)\right|-u_1}+\\ \sqrt{\left|f\left(g_{mf2}\left(i,j,p,d\right),g_{ff2}\left(i,j,p,d\right)\right)\right|-u_2}\end{array}\right)\\ \text{s}\text{.t}\text{.}i,j\in \left[1,n\right];i\ne j\end{array}\right)& \text{­­­(37)}\end{array}$

Type 4

It is assumed that an o ^th sudden voltage sag change point detected in a waveform segment of the voltage sag corresponds to a single-side protection trip. For this type of sag, the following optimization model may be used to infer the faulty line and its short circuit condition, the tripping time point and sequence of the relay protection device, and the time point and the sequence of the DG device disconnected from the power grid (in other words, i, j, p, d, and ^hq ):

$\begin{array}{cc}\left\{\begin{array}{l}\max \text{-}\left(\begin{array}{c}\sqrt{\left|f\left(g_{mf1}\left(i,j,p\right),g_{ff1}\left(i,j,p\right)\right)\right|-u_1}+\\ {\displaystyle \sum _{x=2}^o{\displaystyle \sum _{q=1}^{o-1}\sqrt{\left|f\left(g_{mf3}\left(i,j,p,h_q\right),g_{ff3}\left(i,j,p,h_q\right)\right)\right|-u_x}}}\\ +\left|f\left(g_{mf2}\left(i,j,p,d\right),g_{ff2}\left(i,j,p,d\right)\right)\right|-u_{o+1}+\\ {\displaystyle \sum _{x=o+2}^s{\displaystyle \sum _{h=0}^{s-2}\sqrt{\left|f\left(g_{mf3}\left(i,j,p,h_q\right),g_{ff3}\left(i,j,p,h_q\right)\right)\right|-u_x}}}\end{array}\right)\\ \text{s}\text{.t}\text{.}i,j\in \left[1,n\right];i\ne j;h_1\subseteq h_2\subseteq \cdots \subseteq h_{s-2}\end{array}\right)& \text{­­­(38)}\end{array}$

It should be noted that the formulas (35) to (38) are optimization models for a single monitoring device, which can be obtained through solving by using an intelligent algorithm such as a quantum genetic algorithm. When there are a plurality of monitoring devices, an objective function in the optimization model is the sum of respective objectives functions of the monitoring devices.

6. The state estimation method for the multi-stage voltage sag is inferred.

In voltage sag state estimation, a sag amplitude and duration of an unmonitored bus must be estimated. Duration of each sag stage can be estimated based on the time length between adjacent sudden voltage sag change points. The sag amplitude is estimated based on an event inference result in step 5. Similarly, there may be the following four voltage sag state estimation methods based on the voltage sag types:

Type 1

For this type of sag, a sag amplitude of any unmonitored bus ^m may be estimated according to the following formula:

$\begin{array}{cc}f\left(g_{mf1}\left(i,j,p\right),g_{ff1}\left(i,j,p\right)\right)& \text{­­­(39)}\end{array}$

Type 2

For this type of sag, a sag amplitude of any unmonitored bus ^m may be estimated according to the following formula:

$\begin{array}{cc}\left\{\begin{array}{l}{u}_1=f\left(g_{mf1}\left(i,j,p\right),g_{ff1}\left(i,j,p\right)\right)\\ u_x=f\left(g_{mf3}\left(i,j,p,h_{x-1}\right),g_{ff3}\left(i,j,p,h_{x-1}\right)\right)\end{array}\right)& \text{­­­(40)}\end{array}$

Type 3

For this type of sag, a sag amplitude of any unmonitored bus ^m may be estimated according to the following formula:

$\begin{array}{cc}\left\{\begin{array}{l}{u}_1=f\left(g_{mf1}\left(i,j,p\right),g_{ff1}\left(i,j,p\right)\right)\\ u_2=f\left(g_{mf2}\left(i,j,p,d\right),g_{ff2}\left(i,j,p,d\right)\right)\end{array}\right)& \text{­­­(41)}\end{array}$

Type 4

For this type of sag, a sag amplitude of any unmonitored bus ^m may be estimated according to the following formula:

$\begin{array}{cc}\left\{\begin{array}{l}{u}_1=f\left(g_{mf1}\left(i,j,p\right),g_{ff1}\left(i,j,p\right)\right)\\ u_x=f\left(g_{mf3}\left(i,j,p,h_{x-1}\right),g_{ff3}\left(i,j,p,h_{x-1}\right)\right)\\ u_y=f\left(g_{mf2}\left(i,j,p,d\right),g_{ff2}\left(i,j,p,d\right)\right)\\ u_z=f\left(g_{mf3}\left(i,j,p,h_{z-1}\right),g_{ff3}\left(i,j,p,h_{z-1}\right)\right)\\ x\in \left[2,y-1\right];z\in \left[y+1,s-2\right]\end{array}\right)& \text{­­­(42)}\end{array}$

The method in the present disclosure analyzes characteristics of multi-stage voltage sags due to different causes, considers sudden change time of each stage and fault clearing time of a relay protection device in the monitoring data of the multi-stage voltage sag, infers a specific event that may occur in each stage of the multi-stage voltage sag, and finally obtains a multi-stage state estimation method that meets project requirements.

Claims

1. A state estimation method for a multi-stage voltage sag, comprising the following steps: step 1: analyzing causes of a multi-stage voltage sag, determining characteristics of multi-stage voltage sags due to different causes, and identifying causes of voltage sags at different sudden change points accordingly;step 2: building an action matrix of a relay protection device based on a specified fault clearing time at different positions of a power grid;step 3: obtaining one or more possible faulty lines through solving to form a faulty line set to reduce a calculation amount of state estimation of the multi-stage voltage sag;step 4: calculating system impedance matrices before and after a fault, an impedance matrix in the case of a relay protection action on one side of the faulty line, and an impedance matrix when a distributed generation (DG) device is disconnected from the power grid;step 5: using an optimization model to infer a specific faulty line from the faulty line set obtained in step 3, and inferring a specific event of each stage of the multi-stage voltage sag based on a voltage sag type; andstep 6: estimating duration of each sag stage based on the voltage sag type and a time length between adjacent sudden voltage sag change points, and estimating a sag amplitude based on an event inference result in step 5.

2. The state estimation method for a multi-stage voltage sag according to claim 1, wherein the causes and corresponding characteristics of the multi-stage voltage sag are as follows: cause I: during a short circuit fault, one or more DG devices in a power system trip, resulting in a loss of a power supply of the power grid; and a corresponding characteristic is that a voltage amplitude on a right side of a corresponding sudden voltage sag change point when the voltage sag occurs is smaller than that on a left side; andcause II: during the short circuit fault, relay protection devices on two sides of the faulty line trip at different time, resulting in a change of a topology of the power grid; and a corresponding characteristic is that the voltage amplitude on the right side of the corresponding sudden voltage sag change point when the voltage sag occurs is greater than that on the left side.

3. The state estimation method for a multi-stage voltage sag according to claim 2, wherein the identifying causes of voltage sags at different sudden change points specifically comprises: assuming that for a waveform segment containing the voltage sag, a total of s + 1 sudden voltage sag change points are detected by using a waveform point detection method, which are respectively expressed as MP0,MP1,...,MPx,...,MPs and correspond to s + 1 sudden change time points t0,t1,...,tx,...,ts, wherein these sudden voltage sag change points divide a waveform of the voltage sag into s segments; and voltages of the different waveform segments are respectively expressed as u1,u2,...,ux,...,us, and an m th voltage um is calculated as follows: um=∑j=ff*tm−1+1/fj=ff*tmxj/tm−tm−1+1/f*ff;j∈1,kwherein ƒƒ represents a sampling rate in units of piece/second; ƒ represents a power-frequency current frequency in units of Hz; xj represents a j th value in an effective value sequence; m ∈ [1,s] ; and k represents a total quantity of sampling points; andforming a binary vector W according to the above formula: W=w1,w2,L ,wx,L  ,ws,x∈1,swx=0,if ux<ux−1 or x=11,if ux<ux−1wherein a value of the element wx in the vector W is 0 or 1; when the value of the wx is 0, it indicates that a voltage on a right side of MPx-1 corresponding to a time point tx-1 is less than that on a left side, which means that the multi-stage voltage sag is formed due to the cause I at this time point, wherein x > 1; on the contrary, when the value of the wx is 1, it means that the multi-stage voltage sag is formed due to the cause II at the time point tx-1, wherein x > 1.

4. The state estimation method for a multi-stage voltage sag according to claim 3, wherein the building an action matrix of a relay protection device in step 2 specifically comprises: step 2.1: building a basic action matrix, in other words, a universal matrix describing a fault removal behavior of the relay protection device: Γ=0γ12γ13Lγ1nγ210γ23Lγ2nγ31γ320OMMMO0γn−1nγn1γn2Lγnn−10wherein n represents a quantity of buses in the system, γij and γji respectively represent fault clearing time of protection devices close to i -side and j -side buses on a line lij, i, j ∈ [1, n];i ≠ j, and γij,γji = 0 represents that there is no physical connection between the i -side bus and the j -side bus;step 2.2: improving the universal matrix and decoupling an improved universal matrix into two matrices: ΓΝ˜=00000γ210000γ31γ32000MMO00γn1γn2Lγnn−10ΓΔ=0γ12γ13Lγ1n00γ23Lγ2n000OM0000γn−1n00000wherein the upper and lower triangular matrices respectively represent fault clearing time of protective devices of a same type on two terminal buses of the line; a lower triangular element in the lower triangular matrix represents the fault clearing time of the protection device of the i -side bus of the line lij, and γij=0 and i>j represent that there is no physical connection between the i -side bus and the j -side bus; and an upper triangular element in the upper triangular matrix represents a parameter on the other side of the line; andto represent two-stage protection, two other similar matrices are constructed to represent parameters of backup protection: ΛΝ˜=00000λ210000λ31λ32000MMO00λn1λn2Lλnn−10ΛΔ=0λ12λ13Lλ1n00λ23Lλ2n000OM0000λn−1n00000step 2.3: determining a cooperation relationship of the protection devices, wherein for the two-stage protection, there are four cooperation relationships between main protection and the backup protection in the line, and the improved protection action matrix is expressed as the following cooperation relationships: ΓΝ˜+ΓΔ=ΘIΓΝ˜+ΛΔ=ΘIIΛΝ˜+ΛΔ=ΘIIIΛΝ˜+ΓΔ=ΘIVwherein ΘI represents that the main protection on the two sides of the faulty line cooperates with each other to remove the fault; ΘII represents that the backup protection on the two sides of the faulty line cooperates with each other to remove the fault; ΘIII and ΘIV represent that the main protection on one side of the faulty line cooperates with the backup protection on the other side of the faulty line to remove the fault; andstep 2.4: correcting fault clearing time in the action matrix, wherein the fault clearing time is corrected as follows: γij=γij,set+δ1,ijλij=λij,set+δ2,ij;i,j∈1,n;i≠jwherein γij,set and ηij,set respectively represent specified values of the main protection and the backup protection; δ1,ij and δ2,ij respectively represent deviations between actual fault clearing time and the specified values; and the deviations are random numbers within [0, δ], wherein δ represents a maximum error value during testing or historical operation of relay protection devices of a same model.

5. The state estimation method for a multi-stage voltage sag according to claim 4, wherein step 3 specifically comprises: step 3.1: defining the following four voltage sag types based on the sag cause identification in step 1: type I: single-stage rectangular sag;type II: multi-stage voltage sag due to the cause I;type III: multi-stage voltage sag due to the cause II; andtype IV: multi-stage voltage sag due to the cause I and the cause II; andstep 3.2: obtaining the faulty line set for different types of sags, wherein for the type I and the type II:assuming that time points at which first and last sudden voltage sag change points of the waveform of the voltage sag are detected in step 1 are t0 and ts respectively, a time length from occurrence of the fault to removal of the fault in the system is ts-t0 ; and if two elements of symmetric positions of a main diagonal line of a matrix in the four matrices in step 2.3 are equal to ts-t0 within an error threshold, lines corresponding to these two elements are possible faulty lines, and therefore, a solution model of the faulty line set LF is as follows: LF=L,lij,Ls.t. LF⊆LNθij−ts−t0+θji−ts−t0≤maxδ1,ij,δ1,ji+maxδ2,ij,δ2,ji≤2δi,j∈1,n;i≠jwherein θij represents an element in an i th row and a j th column in the matrices ΘItoΘIV ,and LN represents a set of lines in an intersection of sag domains of a bus of a monitoring device; andfor the type III and the type IV:assuming that the time point at which the first sudden voltage sag change point of the waveform of the voltage sag detected in step 1 is t0, and time points of sudden voltage sag change points of two corresponding relay protection actions are tx-1 and ts respectively, a time length from the occurrence of the fault to the protection actions of the relay protection devices on the two sides of the faulty line are tx-1-t0 and ts-t0 respectively; and in this case, the solution model of the faulty line set is as follows: LF=L,lij,Ls.t. LF⊆LNθij−ts−t0+θji−ts−t0≤maxδ1,ij,δ1,ji+maxδ2,ij,δ2,ji≤2δi,j∈1,n;i≠j;s≥2;wx=1..

6. The state estimation method for a multi-stage voltage sag according to claim 5, wherein step 4 specifically comprises: step 4.1: calculating the system impedance matrix before the fault, wherein a system admittance matrix YSE is expressed as a sum of a line admittance matrix YLSEand a generator admittance matrix,YGSEas shown below:YSE=YGSE+YLSEassuming there are n buses in the power system, the line admittance matrix YLSEis calculated as follows based on a line topology relationship and an impedance parameter:YLSE=α11seα12seLα1nseα21seα22seLα2nseMMOMαn1seαn2seLαnmsewherein se = 1,2,0 represents a positive sequence, a negative sequence, and a zero sequence, αij represents mutual admittance of nodes i and j, αii represents self-admittance of the node i, and i ≠ j; i, j ∈ [1,n];the matrix YGSEis a diagonal matrix, and an element value on a diagonal line is equal to self-admittance of a corresponding generator, as shown below:YGSE=β11se0L00β22seL0MMOM00Lβnmsewherein βii = 0 represents that there is no generator for the bus; andthe system impedance matrix is calculated as follows: ZNSE=YLSE−1=Z11seZ12seLZ1nseZ21seZ22seLZ2nseMMOMZn1seZn2seLZnmsestep 4.2: calculating the system impedance matrix after the short circuit fault, wherein mutual impedance Zmfsebetween a fault position ƒl and a to-be-solved node m after the circuit short fault, and self-impedance Zffseof the fault position ƒl are respectively obtained through solving according to the following two formulas: Zmfse=1−pZmise+pZmjse=gmf1i,j,pZffse=1−p2Ziise+p2Zjjse+2p1−pZijse+p1−pzijse=gjj1i,j,pwherein Zmise,Zmjse,Zijse,Ziise, andZjjseare elements in the matrixZNSE; andzijserepresents impedance of the linelij, wherein the impedance is further expressed as functions gmf1(i,j,p) and gff1(i,j,p), and the two functions respectively represent mutual impedance between the fault position and the target node m and self-impedance of the fault position when the short circuit fault occurs at the position ƒl having a distance P away from the i -side of the line lij;step 4.3: calculating the impedance matrix in the case of the relay protection action on one side of the faulty line, wherein assuming that after the short circuit fault occurs on the line lij, a j -side protection device acts to cut off the line on the corresponding side, and the system impedance is calculated as follows: first, a branch with impedance of -zij is added between the i -side bus and the j -side bus in the original system, and in this case, the system impedance matrix ZATSEis corrected according to the following formula:ZATSE=ZNSE−ΔZΔZT−zijse+Ziise+Zjjse−2Zijse=Z11,ATseZ12,ATseLZ1n,ATseZ21,ATseZ22,ATseLZ2n,ATseMMOMZn1,ATseZn2,ATseLZnn,ATsewherein ΔZ represents a process quantity, which is calculated according to the following formula: ΔZ=Z1jse−Z1iseZ2jse−Z2iseLZnjse−ZniseTthen a branch with impedance of pzij is added to the i -side bus, and in this case, the system impedance matrix is further corrected as ZARPSEaccording to the following formula:ZARPSE=Z11,ATseZ12,ATseLZ1n,ATseZ1i,ATseZ21,ATseZ22,ATseLZ2n,ATseZ2i,ATseMMOMMZn1,ATseZn2,ATseLZnn,ATseZni,ATseZi1,ATseZi2,ATseLZin,ATseZii,ATse+pzijcompared with the ZATSE, theZARPSEadds one line and one column to represent the mutual impedanceZmfsebetween the fault position ƒl and each target bus m or the self-impedance Zffseof the fault position; and the impedance is further expressed as functions gmf 2 (i,j,p,d) and g ff 2 (i, j, p, d), wherein the two functions respectively represent mutual impedance between the fault position and the target node m and self-impedance of the fault position when a d -side protection device on the line acts to cut off a part of the line after the short circuit fault occurs at the position ƒl far from the terminal P of the node i on the line lij, as shown in the following formulas: Zmfse=Zim,ATse=gmf2i,j,p,d;d=i or jZffse=Zii,ATse+pzij=gff2i,j,p,d;d=i or jstep 4.4: calculating the impedance matrix when the DG device is disconnected from the power grid, wherein assuming that DG devices of all buses in a bus set h are disconnected from the power grid, the diagonal matrix first needs to be corrected as follows: YGSE=β11se0L00β22seL0MMOM00Lβnnse;βiise=0;∀i∈hthe system impedance matrix before the fault is calculated according to the above formula, and finally system impedance when the DG device is disconnected from the power grid in the case of the short circuit fault is calculated; and the impedance is further expressed as functions gmf3 (i, j, p, h) and gff3 (i, j, p, h), wherein the two functions respectively represent mutual impedance between the fault position and the target node m and self-impedance of the fault position after the short circuit fault occurs at the position ƒl far from the terminal P of the node i on the line lij and the DG devices of all the buses in the bus set h are disconnected from the power grid, as shown in the following two formulas: Zmfse=gmf3i,j,p,hZffse=gff3i,j,p,h. .

7. The state estimation method for a multi-stage voltage sag according to claim 6, wherein in step 5, the following four optimization models are available based on the sag types: (1) for the type I, the following optimization model is used to infer the faulty line and its short circuit condition, in other words, i, j, and P ; max -fgmf1i,j,p,gff1i,j,p−u1s.t. i,j∈1,n;i≠j,wherein ƒ (·) represents a short-circuit calculation function;(2) for the type II, the following optimization model is used to infer the faulty line and its short circuit condition, a time point and a sequence of disconnecting the DG device from the power grid, in other words, i, j, p, and hq ; max -fgmf1i,j,p,gff1i,j,p−u1+∑x=2s∑q=1s−1fgmf3i,j,p,hq,gff3i,j,p,hq−uxs.t. i,j∈1,n;i≠j;h1⊆h2⊆⋯⊆hs-1 whereins represents a quantity of stages of the multi-stage voltage sag; and hq represents a set of DG devices disconnected from the power grid during a q-stage voltage sag;(3) for the type III, the following optimization model is used to infer the faulty line and its short circuit condition, and a tripping time point and sequence of the relay protection device, in other words, i, j, P, and d ; max -fgmf1i,j,p,gff1i,j,p−u1+fgmf2i,j,p,d,gff2i,j,p,d−u2s.t. i,j∈1,n;i≠j (4) for the type IV, the following optimization model is used to infer the faulty line and its short circuit condition, the tripping time point and sequence of the relay protection device, and the time point and the sequence of disconnecting the DG device from the power grid, in other words, i, j, p, d, and hq; max -fgmf1i,j,p,gff1i,j,p−u1+∑x=2o∑q=1o−1fgmf3i,j,p,hq,gff3i,j,p,hq−ux+fgmf2i,j,p,d,gff2i,j,p,d−uo+1+∑x=o+2s∑q=os−2fgmf3i,j,p,hq,gff3i,j,p,hq−uxs.t. i,j∈1,n;i≠j;h1⊆h2⊆L⊆hs-2 wherein o represents a time point at which a relay protection device on a side of the faulty line first acts to trip; and uo+1 represents an amplitude of an o+1-stage voltage sag.

8. The state estimation method for a multi-stage voltage sag according to claim 7, wherein in step 6, the following four voltage sag state estimation methods are available based on the sag types: 1) for the type I, a sag amplitude of any unmonitored bus m is estimated according to the following formula: fgmf1i,j,p,gff1i,j,p2) for the type II, a sag amplitude of any unmonitored bus m is estimated according to the following formula: u1=fgmf1i,j,p,gff1i,j,pux=fgmf3i,j,p,hx−1,gff3i,j,p,hx−13) for the type III, a sag amplitude of any unmonitored bus m is estimated according to the following formula: u1=fgmf1i,j,p,gff1i,j,pu2=fgmf2i,j,p,d,gff2i,j,p,d4) for the type IV, a sag amplitude of any unmonitored bus m is estimated according to the following formula: u1=fgmf1i,j,p,gff1i,j,pux=fgmf3i,j,p,hx−1,gff3i,j,p,hx−1uy=fgmf2i,j,p,d,gff2i,j,p,duz=fgmf3i,j,p,hz−1,gff3i,j,p,hz−1x∈2,y−1;z∈y+1,s−2wherein y and z respectively represent yth and zth stages of the multi-stage voltage sag.