Method and System for Dynamic State Estimation of Distribution Network Based on Adaptive Set-Membership Filtering

Abstract

The present invention discloses a dynamic state estimation method for a distribution network based on adaptive set membership filtering, including: 1. Collecting distribution parameters and measurement information; 2. Establishing a dynamic state estimation model considering unknown but bounded noise, transforming the hybrid measurement data, linearizing the measurement function, and establishing a linear dynamic state estimation model in combination with a linear state prediction equation; 3. Time update, determining a prior estimate ellipsoid containing the current state ellipsoid based on the previous state ellipsoid and the system state equation; 4. Adaptive bad data detection, defining an adaptive detection threshold to identify and eliminate bad data in the measurements; 5. Measurement correction, using the measurement function and the measured value to correct the state prediction value, and determining a posterior estimate ellipsoid containing the current prior estimate ellipsoid and the measurement update ellipsoid. The present invention can achieve high-precision dynamic state estimation in an environment where the noise is unknown but bounded, can effectively identify and eliminate bad data when the measurement set contains bad data, and is of great significance for formulating the operation plan of the power system and preventing power system accidents.

Description

Technical Field

[0001] The present invention relates to the field of distribution network state estimation, and particularly to a method and system for dynamic state estimation of a distribution network based on adaptive set membership filtering. Background Art

[0002] The increasing scale of distributed power grid connection has exacerbated the randomness and volatility of the operation of the distribution network. The unknown but bounded characteristics of system noise have become increasingly prominent, and its accurate modeling has become more and more difficult. The distribution management system faces the challenge of comprehensively and accurately monitoring and perceiving the operation state of the system. Implementing accurate dynamic state estimation based on the existing measurement system is of great significance for formulating the operation plan of the power system and preventing power system accidents. The development of advanced measurement technologies has made the measurement set of the distribution network composed of data from different measurement systems in the long term. However, the existence of differences in the sampling periods of multi-source measurement data makes the probability of collecting all types of real-time measurements at the same moment small. In most cases, only one or two types of measurements can be obtained. In this measurement environment, the estimation results obtained by the static state estimation method based only on the measurement data at the current moment are difficult to reflect the actual state of the distribution network. The dynamic state estimation method combines prior state prediction and posterior measurement correction, and can achieve system observability in the environment of insufficient measurements, improve the state estimation frequency and the utilization rate of high-frequency measurements. Therefore, studying a dynamic state estimation method that integrates multi-source measurement data is crucial for tracking the operation state of the system.

[0003] The current dynamic state estimation algorithms are mainly the Kalman filter based on the Gaussian noise assumption and its derivative algorithms, such as the extended Kalman filter, the adaptive Kalman filter, etc. These methods have good estimation accuracy in the preset noise environment. However, in the actual power grid where only the boundary information of the noise can be obtained, the state estimation performance will be greatly reduced. For this reason, the set membership filtering algorithm based on the unknown but bounded noise assumption has been favored by researchers and has been gradually applied to the power system, such as the extended set membership filtering, the set membership filtering algorithm based on double-layer optimization, and the set membership filtering method based on the event trigger mechanism, etc. However, most of the current power system state estimation algorithms based on set membership filtering ignore the situation that the measurement set contains bad data, which will lead to the estimation result deviating from the true value. Summary of the Invention

[0004] In order to overcome the deficiencies existing in the above-mentioned prior art, the present invention proposes a method for dynamic state estimation of a distribution network based on adaptive set membership filtering. It is expected to accurately track the state of the distribution network system in the environment where the noise is unknown but bounded, effectively filter out bad measurements, and provide an accurate data basis for the situation awareness and operation planning of the distribution network.

[0005] The present invention solves the above technical problems through the following technical means:

[0006] A dynamic state estimation method for distribution networks based on adaptive set membership filtering, comprising the following steps:

[0007] Step 1. Collect distribution network parameter information and measurement data;

[0008] Step 2. Transform the hybrid measurement data, linearize the measurement function, and establish a linear dynamic state estimation model in combination with the linear state prediction equation;

[0009] Step 3. Update the set membership filtering time; determine a prior estimation ellipsoid containing the current state ellipsoid according to the previous state ellipsoid and the system state equation;

[0010] Step 4. Adaptive bad data detection, define an adaptive detection threshold to identify and eliminate bad data in the measurement; specifically:

[0011] Step 4.1. Calculate the adaptive detection threshold for bad data at time k, and the calculation formula is as follows:

[0012]

[0013] Where: σ ε,k is the equivalent measurement noise standard deviation, L w is the sliding window length, is the mean value of the measurement values within the sliding window, is the variance of the absolute value sequence of the differences between measurements within the sliding window, and the calculation formula is as follows:

[0014]

[0015]

[0016] Where: is the measurement value at time i, is the absolute value of the difference between the measurement at time i - 1 and the measurement at time i within the sliding window, is within the sliding window mean;

[0017] Step 4.2. Bad data identification, compare the absolute value of the measurement difference with the adaptive detection threshold size, if the absolute value of the measurement difference is less than then it is determined that the measurement data is good; if the absolute value of the measurement difference exceeds then it is identified as bad data;

[0018] Step 4.3. Bad data replacement, if the measurement z k+1 at time k + 1 is identified as bad data, then assign the mean value of the remaining measurements within the sliding window to z k+1 , as follows:

[0019]

[0020] Step 5. Measurement calibration: Use the measurement function and the measured value to calibrate the state prediction value, determine a posteriori estimation ellipsoid that includes the prior estimation ellipsoid and the measurement update ellipsoid at the current moment, and estimate the dynamic state of the distribution network using the posteriori estimation ellipsoid.

[0021] Furthermore, the specific execution process of the said Step 1 is as follows:

[0022] Step 1.1: Collect information such as the distribution network topology and line parameters;

[0023] Step 1.2: Collect system measurement data and establish an original measurement set as follows:

[0024] z 0 =[U i ,θ i ,I ij ,θ ij ,P ij ,Q ij ,P i ,Q i T

[0025] where U i and θ i are respectively the amplitude and phase angle measurements of the node i voltage phasor measurement, I ij and θ ij are respectively the amplitude and phase angle measurements of the branch ij current phasor measurement, P ij and Q ij are respectively the active and reactive power measurements of the branch ij, P i and Q i are respectively the active and reactive power injection power measurements of the node i.

[0026] Furthermore, the specific execution process of the said Step 2 is as follows:

[0027] Step 2.1: Multisource measurement linear transformation;

[0028] Step 2.1.1: Set the state variables as the real and imaginary parts of the node voltage phasor as follows:

[0029] x=[e i ,f i T

[0030] where: e i is the real part phasor of the node voltage, f i is the imaginary part phasor of the node voltage; ​​

[0031] Step 2.1.2. Convert the bus voltage phasor into the real and imaginary parts of the equivalent voltage as follows:

[0032]

[0033] where: e i , f i are the real and imaginary part measurements of the equivalent node i voltage respectively;

[0034] Step 2.1.3. Convert the branch current phasor into the real and imaginary parts of the equivalent branch current as follows:

[0035]

[0036] where: are the real and imaginary part measurements of the equivalent branch ij current respectively;

[0037] Step 2.1.4. Convert the branch active power and reactive power into the real and imaginary parts of the equivalent branch current as follows:

[0038]

[0039] where: are the real and imaginary part measurements of the equivalent current of branch ik respectively;

[0040] Step 2.1.5. Convert the node injection power into the real and imaginary parts of the equivalent node injection current as follows:

[0041]

[0042] where: are the real and imaginary part measurements of the equivalent injection current of node i respectively;

[0043] Step 2.2. Build a dynamic state estimation model based on the ellipsoid set;

[0044] Step 2.2.1. Establish a prediction equation as follows:

[0045] x k+1 = F k x k + w k

[0046] where: x k+1 , x k are the state variables at times k + 1 and k respectively; F k is the state transition matrix at time k; w k is the process noise;

[0047] Step 2.2.2. Establish a measurement equation as follows:

[0048] z k = H k x k + v k

[0049] where: z k is the measurement variable at time k; H k is the measurement matrix at time k; v k is the measurement noise;

[0050] Step 2.2.3, establish the ellipsoidal set models of process noise and measurement noise, as follows:

[0051]

[0052]

[0053] where: W k and V k are the ellipsoidal sets of process noise and measurement noise respectively, is a given real positive definite matrix, representing the shape matrix of W k and V k respectively;

[0054] Step 2.3, set the current time k = 0, the stop estimation time k max , use the ellipsoidal set as the shape of the feasible set of state variables, and set the initial set of system state variables, as follows:

[0055]

[0056] where: x0 is the initial state, is the estimated value of x0, representing the center of the ellipsoid, is a given real positive definite matrix, representing the ellipsoidal shape matrix.

[0057] Furthermore, the specific execution process of the said Step 3 is as follows:

[0058] Step 3.1, construct the ellipsoidal set of state variables at time k, as follows:

[0059]

[0060] where: is the posteriori estimated value of x k , representing the center of the ellipsoid, is a given real positive definite matrix, representing the ellipsoidal shape matrix;

[0061] Step 3.2, construct the prior estimated ellipsoidal set of state variables at time k + 1, as follows:

[0062]

[0063] Wherein: is the prior estimated value of x, representing the center of the ellipsoid, k+1 ; is a given real positive definite matrix, representing the ellipsoid shape matrix;

[0064] Step 3.3: Construct the state transition matrix ellipsoid set at time k, and calculate its center and shape matrix as follows:

[0065]

[0066] Where: F k ={F k x k : x k ∈EL k} is the state transition matrix ellipsoid set at time k, is the center of the ellipsoid, and Q fk is the ellipsoid shape matrix;

[0067] Step 3.4: Calculate the center and shape matrix of the prior estimated ellipsoid set of the state variable at time k+1 as follows:

[0068]

[0069] Where: tr(Q) is the sum of the squares of the semi-axis lengths of the ellipsoid Q, representing the size of the ellipsoid shape matrix.

[0070] Furthermore, the specific execution process of step 5 is as follows:

[0071] Step 5.1: Construct the posterior estimated ellipsoid set of the state variable at time k+1 as follows:

[0072]

[0073] Wherein: is the estimated value of x k+1 , representing the center of the ellipsoid, is a given real positive definite matrix, representing the ellipsoid shape matrix;

[0074] Step 5.2: Construct the measurement update ellipsoid set according to the measurement equation as follows:

[0075]

[0076] Where: Z k+1 is the measurement update ellipsoid at time k+1;

[0077] Step 5.3: Construct the measurement update ellipsoid set according to the measurement equation as follows:

[0078]

[0079] where: Z k+1 is the measurement update ellipsoid at time k + 1;

[0080] Step 5.4, construct a minimum ellipsoid that satisfies as follows:

[0081]

[0082] where: for any 0 ≤ ρ k+1 ≤ 1, the above equation holds;

[0083] Step 5.5, rewrite the above equation:

[0084]

[0085] where: 0 ≤ δ k+1 ≤ 1, is an intermediate variable. If we let then we can obtain the measurement correction ellipsoid EL k+1 ;

[0086] Step 5.6, solve for the center and shape matrix of the posterior estimate ellipsoid set of the state variable at time k + 1 according to the measurement correction iteration formula, as follows:

[0087]

[0088]

[0089]

[0090]

[0091] where: the optimal value of the parameter ρ k+1 is obtained by solving Equation (29):

[0092]

[0093] Step 5.7, determine whether the state estimation calculation stops. If k < k max , then let k = k + 1 and jump to Step 3 to continue the calculation; otherwise, end the calculation.

[0094] The present invention also discloses a distribution network dynamic state estimation system based on adaptive set membership filtering, which is applied to the above method and includes:

[0095] A data acquisition module for collecting distribution network parameter information and measurement data;

[0096] A data transformation module, which is used to transform the mixed measurement data, linearize the measurement function, and establish a linear dynamic state estimation model by combining a linear state prediction equation;

[0097] A set membership filtering module, which is used for set membership filtering time update; determining a prior estimation ellipsoid containing the current moment state ellipsoid according to the previous moment state ellipsoid and the system state equation;

[0098] An outlier data adaptive detection module, which is used for outlier data adaptive detection, defining an adaptive detection threshold to identify and eliminate outlier data in the measurement; specifically:

[0099] Step 4.1: Calculate the outlier data adaptive detection threshold at time k, and the calculation formula is as follows:

[0100]

[0101] Where: σ ε,k is the equivalent measurement noise standard deviation, L w is the sliding window length, is the mean value of the measurement values within the sliding window, is the variance of the absolute value sequence of the differences between measurements within the sliding window, and the calculation formula is as follows:

[0102]

[0103]

[0104] Where: is the measurement value at time i, is the absolute value of the difference between the measurement at time i - 1 and the measurement at time i within the sliding window, is within the sliding window mean value of;

[0105] Step 4.2: Outlier data identification, comparing the absolute value of the measurement difference with the adaptive detection threshold of the size, if the absolute value of the measurement difference is less than then it is determined that the measurement data is good; if the absolute value of the measurement difference exceeds then it is identified as outlier data;

[0106] Step 4.3: Outlier data replacement, if the measurement z at time k + 1 k+1 is identified as outlier data, then assign the mean value of the remaining measurements within the sliding window to z k+1 , as follows:

[0107]

[0108] The measurement correction module is used for measurement correction. It corrects the state prediction value by using the measurement function and the measurement value, determines a posteriori estimation ellipsoid that includes the prior estimation ellipsoid and the measurement update ellipsoid at the current moment, and estimates the dynamic state of the distribution network by using the posteriori estimation ellipsoid.

[0109] Furthermore, the specific execution process of the data acquisition module is as follows:

[0110] Step 1.1: Collect information such as the distribution network topology and line parameters;

[0111] Step 1.2: Collect system measurement data and establish an original measurement set as follows:

[0112] z 0 =[U i ,θ i ,I ij ,θ ij ,P ij ,Q ij ,P i ,Q i T

[0113] Where U i 、θ i are respectively the amplitude and phase angle measurements of the node i voltage phasor measurement, I ij 、θ ij are respectively the amplitude and phase angle measurements of the branch ij current phasor measurement, P ij 、Q ij are respectively the active and reactive power measurements of the branch ij, P i 、Q i are respectively the active and reactive power injection power measurements of the node i.

[0114] Furthermore, the specific execution process of the data transformation module is as follows:

[0115] Step 2.1: Multisource measurement linear transformation;

[0116] Step 2.1.1: Set the state variables as the real and imaginary parts of the node voltage phasor as follows:

[0117] x=[e i ,f i T

[0118] Where: e i is the real part phasor of the node voltage, f i is the imaginary part phasor of the node voltage;

[0119] Step 2.1.2: Convert the bus voltage phasor to the real and imaginary parts of the equivalent voltage as follows: ​​

[0120]

[0121] where: e i and f i are the measured real and imaginary parts of the voltage of equivalent node i, respectively;

[0122] Step 2.1.3. Convert the branch current phasor into the real and imaginary parts of the equivalent branch current as follows:

[0123]

[0124] where: are the measured real and imaginary parts of the current of equivalent branch ij, respectively;

[0125] Step 2.1.4. Convert the active and reactive powers of the branch into the real and imaginary parts of the equivalent branch current as follows:

[0126]

[0127] where: are the measured real and imaginary parts of the equivalent current of branch ik, respectively;

[0128] Step 2.1.5. Convert the node injection power into the real and imaginary parts of the equivalent node injection current as follows:

[0129]

[0130] where: are the measured real and imaginary parts of the equivalent injection current of node i, respectively;

[0131] Step 2.2. Construct a dynamic state estimation model based on the ellipsoidal set;

[0132] Step 2.2.1. Establish a prediction equation as follows:

[0133] x k+1 = F k x k + w k

[0134] where: x k+1 and x k are the state variables at times k + 1 and k, respectively; F k is the state transition matrix at time k; w k is the process noise;

[0135] Step 2.2.2. Establish a measurement equation as follows:

[0136] z k = H k xk +v k

[0137] where: z k is the measurement variable at time k; H k is the measurement matrix at time k; v k is the measurement noise;

[0138] Step 2.2.3, establish the ellipsoidal set models of process noise and measurement noise, as follows:

[0139]

[0140]

[0141] where: W k and V k are the ellipsoidal sets of process noise and measurement noise respectively, is a given real positive definite matrix, representing the shape matrices of W k and V k respectively;

[0142] Step 2.3, set the current time k = 0, the stop estimation time k max , use the ellipsoidal set as the shape of the feasible set of state variables, and set the initial set of system state variables, as follows:

[0143]

[0144] where: x0 is the initial state, is the estimated value of x0, representing the ellipsoid center, is a given real positive definite matrix, representing the ellipsoid shape matrix.

[0145] Furthermore, the specific execution process of the set membership filtering module is as follows:

[0146] Step 3.1, construct the ellipsoidal set of state variables at time k, as follows:

[0147]

[0148] where: is the posteriori estimated value of x k , representing the ellipsoid center, is a given real positive definite matrix, representing the ellipsoid shape matrix;

[0149] Step 3.2, construct the prior estimated ellipsoidal set of state variables at time k+1, as follows:

[0150]

[0151] where: is xk+1 The prior estimated value of, representing the center of the ellipsoid, is a given real positive definite matrix, representing the ellipsoid shape matrix;

[0152] Step 3.3: Construct the state transition matrix ellipsoid set at time k, and calculate its center and shape matrix, as shown in the following formula:

[0153]

[0154] where: F k ={F k x k : x k ∈EL k} is the state transition matrix ellipsoid set at time k, is the center of the ellipsoid, Q fk is the ellipsoid shape matrix;

[0155] Step 3.4: Calculate the center and shape matrix of the prior estimated ellipsoid set of the state variable at time k + 1, as shown in the following formula:

[0156]

[0157] where: tr(Q) is the sum of the squares of the semi-axis lengths of the ellipsoid Q, representing the size of the ellipsoid shape matrix.

[0158] Furthermore, the specific execution process of the measurement correction module is as follows:

[0159] Step 5.1: Construct the posterior estimated ellipsoid set of the state variable at time k + 1, as shown in the following formula:

[0160]

[0161] where: is the estimated value of x k+1 representing the center of the ellipsoid, is a given real positive definite matrix, representing the ellipsoid shape matrix;

[0162] Step 5.2: Construct the measurement update ellipsoid set according to the measurement equation, as shown in the following formula:

[0163]

[0164] where: Z k+1 is the measurement update ellipsoid at time k + 1;

[0165] Step 5.3: Construct the measurement update ellipsoid set according to the measurement equation, as shown in the following formula:

[0166]

[0167] where: Z k+1is the measurement update ellipsoid at time k+1;

[0168] Step 5.4, construct a minimum ellipsoid that satisfies as follows:

[0169]

[0170] where: for any 0 ≤ ρ k+1 ≤ 1, the above equation holds;

[0171] Step 5.5, rewrite the above equation:

[0172]

[0173] where: 0 ≤ δ k+1 ≤ 1, is an intermediate variable. If we let then we can obtain the measurement correction ellipsoid EL k+1 ;

[0174] Step 5.6, solve for the center and shape matrix of the posterior estimate ellipsoid set of the state variable at time k+1 according to the measurement correction iteration formula, as follows:

[0175]

[0176]

[0177]

[0178]

[0179] where: the optimal value of the parameter ρ k+1 is obtained by solving equation (29):

[0180]

[0181] Step 5.7, determine whether the state estimation calculation stops. If k < k max , then let k = k + 1 and jump to step three to continue the calculation; otherwise, end the calculation.

[0182] The advantages of the present invention are as follows:

[0183] 1. The present invention addresses the problems of incompatible multi-source measurement data forms and high computational complexity of nonlinear state estimation models. By linearly transforming the multi-source measurement set into a measurement set that only contains voltage and current phasors, linearizing the measurement function, and combining the state prediction equation, a linear dynamic state estimation model is established, reducing the computational complexity of state estimation for nonlinear power systems.

[0184] 2. In view of the problem that the performance of traditional dynamic state estimation algorithms based on the Gaussian noise assumption is greatly reduced in an unknown but bounded noise environment, the present invention uses an adaptive set membership filtering algorithm to solve the dynamic state estimation model, which can achieve high-precision dynamic state estimation, solve the problem of system unobservability caused by insufficient measurements during the sampling interval of low-frequency measurements, and improve the utilization rate of high-frequency measurements.

[0185] 3. In view of the problem that it is difficult for the set membership filtering method to monitor bad measurement data, the present invention introduces an adaptive detection threshold for bad data to pre-screen the measurement data adaptively and filter out bad measurements. When the measurement set contains bad data, the proposed algorithm can adaptively adjust the detection threshold according to the measurement sequence, effectively filter out bad data, and improve the robustness of state estimation. BRIEF DESCRIPTION OF THE DRAWINGS

[0186] Figure 1 It is a flowchart for solving the dynamic state estimation model based on adaptive set membership filtering of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0187] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the embodiments of the present invention. Obviously, the described embodiments are some, rather than all, of the embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0188] Embodiment 1

[0189] In this embodiment, a method for dynamic state estimation of a distribution network based on adaptive set membership filtering is provided. First, in view of the differences in various forms of measurement data, a multi-source measurement linear transformation strategy is proposed to simplify the multi-type measurement set into a measurement set containing only voltage and current phasors. Secondly, the measurement function is linearized, and a dynamic state estimation model based on an ellipsoidal set is established in combination with the linear state prediction equation. Then, an adaptive set membership filtering algorithm considering bad data detection is used to solve the proposed dynamic state estimation model, and the state quantities at different times are iteratively solved through time update, adaptive bad data detection, and measurement correction. Specifically, it is carried out according to the following steps:

[0190] Step 1: Collect the parameter information and measurement data of the distribution network;

[0191] Step 1.1: Collect information such as the distribution network topology and line parameters;

[0192] Step 1.2: Collect the system measurement data and establish an original measurement set, as shown in Equation (1):

[0193] z 0 =[Ui , θ i , I ij , θ ij , P ij , Q ij , P i , Q i T

[0194] In Equation (1), U i , θ i are respectively the magnitude and phase angle measurements of the node i voltage phasor measurement, and I ij , θ ij are respectively the magnitude and phase angle measurements of the branch ij current phasor measurement, and P ij , Q ij are respectively the active and reactive power measurements of branch ij, and P i , Q i are respectively the active and reactive power injection power measurements of node i;

[0195] Step 2: Establish a dynamic state estimation model considering unknown but bounded noise;

[0196] Step 2.1: Linear transformation of multi-source measurements;

[0197] Step 2.1.1: Set the state variables as the real and imaginary parts of the node voltage phasor, as shown in Equation (2):

[0198] x = [e i , f i T

[0199] In Equation (2): e i is the real part phasor of the node voltage, and f i is the imaginary part phasor of the node voltage;

[0200] Step 2.1.2: Convert the bus voltage phasor to the real and imaginary parts of the equivalent voltage, as shown in Equation (3):

[0201]

[0202] In Equation (3): e i , f i are respectively the real and imaginary part measurements of the equivalent node i voltage;

[0203] Step 2.1.3: Convert the branch current phasor to the real and imaginary parts of the equivalent branch current, as shown in Equation (4):

[0204]

[0205] In Equation (4): ​​They are the real and imaginary part measurements of the equivalent branch ij current respectively.

[0206] Step 2.1.4: Convert the active and reactive powers of the branch into the real and imaginary parts of the equivalent branch current, as shown in Equation (5):

[0207]

[0208] In Equation (5): They are the real and imaginary part measurements of the equivalent current of branch ik respectively.

[0209] Step 2.1.5: Convert the nodal injection power into the real and imaginary parts of the equivalent nodal injection current, as shown in Equation (6):

[0210]

[0211] In Equation (6): They are the real and imaginary part measurements of the equivalent injection current of node i respectively.

[0212] Step 2.2: Build a dynamic state estimation model based on the ellipsoidal set.

[0213] Step 2.2.1: Establish a prediction equation to predict the state quantity at the current time through the posterior state quantity at the previous time, as shown in Equation (7):

[0214] x k+1 = F k x k + w k

[0215] In Equation (7): x k+1 and x k are the state variables at times k + 1 and k respectively; F k is the state transition matrix at time k; w k is the process noise.

[0216] Step 2.2.2: Establish a measurement equation to calculate the state quantity at the current time through the measurement quantity at the current time, as shown in Equation (8):

[0217] z k = H k x k + v k

[0218] In Equation (8): z k is the measurement variable at time k; H k is the measurement matrix at time k; v k is the measurement noise.

[0219] Step 2.2.3, establish the ellipsoidal set models of process noise and measurement noise, as shown in Equations (9) and (10):

[0220]

[0221]

[0222] In Equations (9) and (10): W k and V k are the ellipsoidal sets of process noise and measurement noise respectively, are given real positive definite matrices, representing the shape matrices of W k and V k respectively;

[0223] Step 2.3, set the current time k = 0, the stop estimation time k max , use the ellipsoidal set as the shape of the feasible set of state variables, and set the initial set of system state variables, as shown in Equation (11):

[0224]

[0225] In Equation (11): x0 is the initial state, is the estimated value of x0, representing the ellipsoid center, is a given real positive definite matrix, representing the ellipsoid shape matrix;

[0226] Step Three, set membership filtering time update;

[0227] Determine a prior estimated ellipsoid that includes the current time state ellipsoid based on the previous time state ellipsoid and the prediction equation.

[0228] Step 3.1, construct the ellipsoidal set of state variables at time k, as shown in Equation (12):

[0229]

[0230] In Equation (12): is the posterior estimated value of x k , representing the ellipsoid center, is a given real positive definite matrix, representing the ellipsoid shape matrix;

[0231] Step 3.2, construct the prior estimated ellipsoidal set of state variables at time k+1, as shown in Equation (13):

[0232]

[0233] In Equation (13): is the prior estimated value of x k+1 , representing the ellipsoid center, For a given real positive definite matrix, representing the ellipsoid shape matrix;

[0234] Step 3.3: Construct the ellipsoidal set of the state transition matrix at time k, and calculate its center and shape matrix, as shown in Equation (14):

[0235]

[0236] In Equation (14): F k ={F k x k :x k ∈EL k} is the ellipsoidal set of the state transition matrix at time k, is the center of the ellipsoid, Q fk is the shape matrix of the ellipsoid;

[0237] Step 3.4: Calculate the center and shape matrix of the prior estimate ellipsoidal set of the state variable at time k+1, as shown in Equation (15):

[0238]

[0239] In Equation (15): tr(Q) is the sum of the squares of the semi-axis lengths of the ellipsoid Q, representing the size of the shape matrix of the ellipsoid;

[0240] Step Four: Set membership filtering for adaptive bad data detection;

[0241] According to the characteristic that bad data and normal data have weak correlation in time, an adaptive detection mechanism is introduced to identify and replace bad data in the measurement.

[0242] The basic idea for establishing the adaptive detection threshold of bad data is as follows: When the power system is operating normally, since the measurement values at adjacent sampling times change very little, it can be considered that the absolute value of the difference between the measurements at time k and time k-1, |z k -z k-1 | only contains noise signals; when there is bad data at time k, |z k -z k-1 | will also contain mutation signals in addition to noise signals. Based on this, bad data can be identified by comparing the size of the adaptive detection threshold with the measurement difference.

[0243] Step 4.1: Calculate the adaptive detection threshold of bad data at time k, and the calculation formula is as shown in Equation (16):

[0244]

[0245] In Equation (16): σ ε,k is the equivalent measurement noise standard deviation, L w is the sliding window length, is the mean value of the measurement values within the sliding window, is the variance of the sequence of absolute values of the differences between measurements within the sliding window, and the calculation formulas are as shown in Equations (17) and (18):

[0246]

[0247]

[0248] In Equations (17) and (18): is the measurement value at time i, is the absolute value of the difference between the measurement at time i - 1 and the measurement at time i within the sliding window, is within the sliding window mean value;

[0249] In Equation (18): and The calculation formulas are as follows:

[0250]

[0251] Step 4.2, Bad data identification, compare the absolute value of the measurement difference with the adaptive detection threshold If the absolute value of the measurement difference is less than then it is determined that the measurement data is good; if the absolute value of the measurement difference exceeds then it is identified as bad data;

[0252] Analyzing Equations (16) to (18), it can be obtained that when there is no bad data, the absolute value of the measurement difference between adjacent times only contains noise signals, the value is small, and the detection threshold the value is high, which can ensure that the absolute value of the measurement difference is within range, thus determining that the measurement data is good; when there is bad data, there will be an absolute value of the measurement difference containing mutation signals and noise signals, the value rises, and the detection threshold the value decreases, making the absolute value of the measurement difference exceed range, thus identifying it as bad data.

[0253] Step 4.3, Bad data replacement, if the measurement z k+1 at time k + 1 is identified as bad data, then assign the mean value of the remaining measurements within the sliding window to z k+1 , as shown in Equation (19):

[0254]

[0255] Step Five, Set - membership filtering measurement correction;

[0256] Determine a minimum ellipsoid that includes the prior estimate ellipsoid and the measurement update ellipsoid at the current time based on the measurement value and the measurement function at the current time.

[0257] Step 5.1. Construct the posterior estimate ellipsoid set of the state variable at time k + 1, as shown in Equation (20):

[0258]

[0259] In Equation (20): is the estimated value of x k+1 indicating the center of the ellipsoid, is a given real positive definite matrix representing the ellipsoid shape matrix;

[0260] Step 5.2. Construct the measurement update ellipsoid set according to the measurement equation, as shown in Equation (21):

[0261]

[0262] In Equation (21): Z k+1 is the measurement update ellipsoid at time k + 1;

[0263] Step 5.3. Construct the measurement update ellipsoid set according to the measurement equation, as shown in Equation (22):

[0264]

[0265] In Equation (22): Z k+1 is the measurement update ellipsoid at time k + 1;

[0266] Step 5.4. Construct a minimum ellipsoid that satisfies as shown in Equation (23):

[0267]

[0268] In Equation (23): For any 0 ≤ ρ k+1 ≤ 1, the above equation holds.

[0269] Step 5.5. Rewrite Equation (23) as Equation (24):

[0270]

[0271] In Equation (24): 0 ≤ δ k+1 ≤ 1, is an intermediate variable. If we let then we can obtain the measurement correction ellipsoid EL k+1 ;

[0272] Step 5.6: Solve for the center and shape matrix of the posterior estimate ellipsoid set of the state variables at time k+1 according to the measurement correction iteration formula, as shown in Equations (25) to (28):

[0273]

[0274]

[0275]

[0276]

[0277] In Equation (27): The optimal value of the parameter ρ k+1 is obtained by solving Equation (29):

[0278]

[0279] Step 5.7: Determine whether the state estimation calculation stops. If k < k max , then set k = k + 1 and jump to Step 3 to continue the calculation; otherwise, end the calculation.

[0280] In summary, the present invention first performs a linear transformation on multi-source measurements to simplify the multi-type measurement set into a measurement set that only includes voltage and current phasors; then linearizes the measurement function and establishes a dynamic state estimation model in combination with the linear state prediction equation; finally, uses the adaptive set membership filtering algorithm considering bad data detection to solve the proposed dynamic state estimation model, and obtains the state quantities at different times through time update, bad data adaptive detection, and measurement correction iteration.

[0281] Embodiment 2

[0282] This embodiment discloses a distribution network dynamic state estimation system based on adaptive set membership filtering, which is applied to the above method and includes:[[]]

[0283] A data acquisition module for acquiring distribution network parameter information and measurement data;

[0284] A data transformation module for transforming the hybrid measurement data, linearizing the measurement function, and establishing a linear dynamic state estimation model in combination with the linear state prediction equation;

[0285] A set membership filtering module for set membership filtering time update; determining a prior estimate ellipsoid that includes the current time state ellipsoid according to the previous time state ellipsoid and the system state equation;

[0286] A bad data adaptive detection module for bad data adaptive detection, defining an adaptive detection threshold to identify and eliminate bad data in the measurements;

[0287] A measurement correction module is used for measurement correction. It corrects the state prediction value by using the measurement function and the measurement value, determines a posteriori estimation ellipsoid that includes the prior estimation ellipsoid and the measurement update ellipsoid at the current moment, and estimates the dynamic state of the distribution network by using the posteriori estimation ellipsoid.

[0288] Each of the above modules executes each step in Embodiment 1, which will not be elaborated here.

[0289] The above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that: they can still modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements on some of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of each embodiment of the present invention.

Claims

1. A distribution network dynamic state estimation method based on adaptive set membership filtering, characterized in that: The following steps are involved: Step 1. Collect distribution network parameter information and measurement data; Step 2. Transform the mixed measurement data, linearize the measurement function, and establish a linear dynamic state estimation model based on the linear state prediction equation; Step 3. Update the set membership filter time; determine a priori estimated ellipsoid that includes the current state ellipsoid based on the state ellipsoid at the previous moment and the system state equation; Step 4. Adaptive detection of bad data: define adaptive detection thresholds to identify and eliminate bad data in the measurement; specifically: Step 4.1: Calculate the adaptive detection threshold of bad data at time k. The calculation formula is as follows: Where: ε,k is the equivalent measurement noise standard deviation, L w is the sliding window length, is the mean value of the measurement in the sliding window, is the variance of the absolute value sequence of the difference between measurements within the sliding window, and is calculated as follows: in: is the measured value at time i, is the absolute value of the difference between the measurement at time i-1 and the measurement at time i within the sliding window, In the sliding window The mean of Step 4.2: Identify bad data and compare the absolute value of the measured difference with the adaptive detection threshold If the absolute value of the measured difference is less than The measurement data is judged to be good; if the absolute value of the measurement difference exceeds It is identified as bad data; Step 4.3: Replace bad data. If the measurement z at time k+1 is k+1 If it is identified as bad data, the mean of the remaining measurements in the sliding window is assigned to z k+1 , as follows: Step 5. Measurement correction: Use the measurement function and the measurement value to correct the state prediction value, determine a posterior estimation ellipsoid that includes the current moment prior estimation ellipsoid and the measurement update ellipsoid, and use the posterior estimation ellipsoid to estimate the dynamic state of the distribution network.

2. The method for dynamic state estimation of a distribution network based on adaptive set membership filtering according to claim 1, characterized in that: The specific execution process of step 1 is as follows: Step 1.1, collect information such as distribution network topology and line parameters; Step 1.2: Collect system measurement data and create an original measurement set as follows: z 0 =[U i ,θ i ,I ij ,θ ij ,P ij ,Q ij ,P i ,Q i ] T Among them, U i ,θ i are the amplitude and phase angle measurements of the voltage phasor of node i, I ij ,θ ij They are the amplitude and phase angle of the branch ij current phasor measurement, P ij , Q ij They are respectively the active and reactive power measurements of branch ij, P i , Q i are the active and reactive injected power measurements of node i respectively.

3. The method for dynamic state estimation of distribution network based on adaptive set membership filtering according to claim 1, characterized in that: The specific execution process of step 2 is: Step 2.1, multi-source measurement linear transformation; Step 2.1.

1. Set the state variables to the real and imaginary phasors of the node voltages as follows: x=[e i ,f i ] T Among them: e i is the real part phasor of the node voltage, f i is the imaginary part phasor of node voltage; Step 2.1.2: Convert the bus voltage phasor into the real and imaginary parts of the equivalent voltage as follows: Among them: e i 、f i They are the real and imaginary part measurements of the voltage at the equivalent node i respectively; Step 2.1.

3. Convert the branch current phasor into the real and imaginary parts of the equivalent branch current as follows: in: They are respectively the real and imaginary part measurements of the equivalent branch ij current; Step 2.1.4: Convert the branch active power and reactive power into the real and imaginary parts of the equivalent branch current as follows: in: They are respectively the real and imaginary part measurements of the equivalent current of branch ik; Step 2.1.

5. Convert the node injection power into the real and imaginary parts of the equivalent node injection current as follows: in: are the real and imaginary part measurements of the equivalent injected current at node i respectively; Step 2.2, construct a dynamic state estimation model based on the ellipsoid set; Step 2.2.1: Establish the prediction equation as follows: x k+1 =F k x k +w k Where: x k+1 、x k are the state variables at time k+1 and k respectively; F k is the state transfer matrix at time k; w k is the process noise; Step 2.2.2: Establish the measurement equation as follows: z k =H k x k +v k Where: z k is the measured variable at time k; H k is the measurement matrix at time k; v k To measure noise; Step 2.2.3: Establish the ellipsoid set model of process noise and measurement noise as follows: Where: W k and V k are the ellipsoid sets of process noise and measurement noise respectively, is a given real positive definite matrix, representing W k and V k The shape matrix of Step 2.3: Set the current time k = 0 and stop estimating the time k max , using the ellipsoid set as the shape of the state variable feasible set, setting the initial set of system state variables as follows: Among them: x0 is the initial state, is the estimated value of x0, representing the center of the ellipsoid, Given a real positive definite matrix, represents the ellipsoidal shape matrix.

4. A method for dynamic state estimation of a distribution network based on adaptive set membership filtering according to any one of claims 1 to 3, characterized in that: The specific execution process of step 3 is: Step 3.1: Construct the k-time state variable ellipsoid set as follows: in: is x k The posterior estimate of , which represents the center of the ellipsoid, For a given real positive definite matrix, represents the ellipsoid shape matrix; Step 3.2: Construct the k+1 time state variable prior estimation ellipsoid set as follows: in: is x k+1 The a priori estimate of , representing the center of the ellipsoid, For a given real positive definite matrix, represents the ellipsoid shape matrix; Step 3.3: Construct the k-time state transition matrix ellipsoid set and calculate its center and shape matrix as follows: Among them: F k ={F k x k :x k ∈EL k } is the ellipsoid set of the state transfer matrix at time k, is the center of the ellipsoid, Q fk is the ellipsoid shape matrix; Step 3.4, calculate the center and shape matrix of the prior estimated ellipsoid set of state variables at time k+1, as follows: Where: tr(Q) is the sum of the squares of the lengths of the semi-axes of the ellipsoid Q, representing the size of the ellipsoid shape matrix.

5. A method for dynamic state estimation of a distribution network based on adaptive set membership filtering according to any one of claims 1 to 3, characterized in that: The specific execution process of step 5 is as follows: Step 5.1: Construct the posterior estimation ellipsoid set of the state variables at time k+1, as follows: in: is x k+1 An estimate of , representing the ellipsoid center, For a given real positive definite matrix, represents the ellipsoid shape matrix; Step 5.2: Construct the measurement update ellipsoid set according to the measurement equation, as follows: Where: Z k+1 Update the ellipsoid for the measurement at time k+1; Step 5.3: Construct the measurement update ellipsoid set according to the measurement equation, as follows: Where: Z k+1 Update the ellipsoid for the measurement at time k+1; Step 5.4: Build a satisfying The minimum ellipsoid is as follows: Where: for any 0≤ρ k+1 ≤1, the above formula is valid; Step 5.5: Rewrite the above formula: Where: 0≤δ k+1 ≤1, is an intermediate variable, if The measurement correction ellipsoid EL can be obtained k+1 ; Step 5.6: Solve the center and shape matrix of the ellipsoid set of the state variables at time k+1 according to the measurement correction iterative formula, as follows: Where: parameter ρ k+1 The optimal value of is obtained by solving formula (29): Step 5.7: Determine whether the state estimation calculation is stopped. If k <k max , then set k = k + 1 and jump to step 3 to continue the calculation; otherwise, end the calculation.

6. A distribution network dynamic state estimation system based on adaptive set membership filtering, characterized in that: include: Data acquisition module, used to collect distribution network parameter information and measurement data; The data transformation module is used to transform the mixed measurement data, linearize the measurement function, and establish a linear dynamic state estimation model by combining the linear state prediction equation; The ensemble membership filter module is used to update the ensemble membership filter time; a priori estimated ellipsoid containing the current state ellipsoid is determined based on the state ellipsoid at the previous moment and the system state equation; The bad data adaptive detection module is used for the adaptive detection of bad data and defines the adaptive detection threshold to identify and eliminate bad data in the measurement; specifically: Step 4.1: Calculate the adaptive detection threshold of bad data at time k. The calculation formula is as follows: Where: ε,k is the equivalent measurement noise standard deviation, L w is the sliding window length, is the mean value of the measurement in the sliding window, is the variance of the absolute value sequence of the difference between measurements within the sliding window, and is calculated as follows: in: is the measured value at time i, is the absolute value of the difference between the measurement at time i-1 and the measurement at time i within the sliding window, In the sliding window The mean of Step 4.2: Identify bad data and compare the absolute value of the measured difference with the adaptive detection threshold If the absolute value of the measured difference is less than The measurement data is judged to be good; if the absolute value of the measurement difference exceeds It is identified as bad data; Step 4.3: Replace bad data. If the measurement z at time k+1 is k+1 If it is identified as bad data, the mean of the remaining measurements in the sliding window is assigned to z k+1 , as follows: The measurement correction module is used for measurement correction. It uses the measurement function and the measurement value to correct the state prediction value, determines a posterior estimation ellipsoid including the current moment prior estimation ellipsoid and the measurement update ellipsoid, and uses the posterior estimation ellipsoid to estimate the dynamic state of the distribution network.

7. A distribution network dynamic state estimation system based on adaptive set membership filtering according to claim 6, characterized in that: The specific execution process of the data acquisition module is as follows: Step 1.1, collect information such as distribution network topology and line parameters; Step 1.2: Collect system measurement data and create an original measurement set as follows: z 0 =[U i ,θ i ,I ij ,θ ij ,P ij ,Q ij ,P i ,Q i ] T Among them, U i ,θ i are the amplitude and phase angle measurements of the voltage phasor of node i, I ij ,θ ij They are the amplitude and phase angle of the branch ij current phasor measurement, P ij , Q ij They are respectively the active and reactive power measurements of branch ij, P i , Q i are the active and reactive injected power measurements of node i respectively.

8. The distribution network dynamic state estimation system based on adaptive set membership filtering according to claim 6, characterized in that: The specific execution process of the data conversion module is as follows: Step 2.1, multi-source measurement linear transformation; Step 2.1.

1. Set the state variables to the real and imaginary phasors of the node voltages as follows: x=[e i ,f i ] T Among them: e i is the real part phasor of the node voltage, f i is the imaginary part phasor of node voltage; Step 2.1.2: Convert the bus voltage phasor into the real and imaginary parts of the equivalent voltage as follows: Among them: e i 、f i They are the real and imaginary part measurements of the voltage at the equivalent node i respectively; Step 2.1.

3. Convert the branch current phasor into the real and imaginary parts of the equivalent branch current as follows: in: They are respectively the real and imaginary part measurements of the equivalent branch ij current; Step 2.1.4: Convert the branch active power and reactive power into the real and imaginary parts of the equivalent branch current as follows: in: They are respectively the real and imaginary part measurements of the equivalent current of branch ik; Step 2.1.

5. Convert the node injection power into the real and imaginary parts of the equivalent node injection current as follows: in: are the real and imaginary part measurements of the equivalent injected current at node i respectively; Step 2.2, construct a dynamic state estimation model based on the ellipsoid set; Step 2.2.

1. Establish the prediction equation as follows: x k+1 =F k x k +w k Where: x k+1 、x k are the state variables at time k+1 and k respectively; F k is the state transfer matrix at time k; w k is the process noise; Step 2.2.2: Establish the measurement equation as follows: z k =H k x k +v k Where: z k is the measured variable at time k; H k is the measurement matrix at time k; v k To measure noise; Step 2.2.3: Establish the ellipsoid set model of process noise and measurement noise as follows: Where: W k and V k are the ellipsoid sets of process noise and measurement noise respectively, is a given real positive definite matrix, representing W k and V k The shape matrix of Step 2.3: Set the current time k = 0 and stop estimating the time k max , using the ellipsoid set as the shape of the state variable feasible set, setting the initial set of system state variables as follows: Among them: x0 is the initial state, is the estimated value of x0, representing the center of the ellipsoid, Given a real positive definite matrix, represents the ellipsoidal shape matrix.

9. A distribution network dynamic state estimation system based on adaptive set membership filtering according to any one of claims 6 to 8, characterized in that: The specific execution process of the set membership filtering module is as follows: Step 3.1: Construct the k-time state variable ellipsoid set as follows: in: is x k The posterior estimate of , which represents the center of the ellipsoid, For a given real positive definite matrix, represents the ellipsoid shape matrix; Step 3.2: Construct the k+1 time state variable prior estimation ellipsoid set as follows: in: is x k+1 The a priori estimate of , representing the center of the ellipsoid, For a given real positive definite matrix, represents the ellipsoid shape matrix; Step 3.3: Construct the k-time state transition matrix ellipsoid set and calculate its center and shape matrix as follows: Among them: F k ={F k x k :x k ∈EL k } is the ellipsoid set of the state transfer matrix at time k, is the center of the ellipsoid, Q fk is the ellipsoid shape matrix; Step 3.4, calculate the center and shape matrix of the prior estimated ellipsoid set of state variables at time k+1, as follows: Where: tr(Q) is the sum of the squares of the lengths of the semi-axes of the ellipsoid Q, representing the size of the ellipsoid shape matrix.

10. A distribution network dynamic state estimation system based on adaptive set membership filtering according to any one of claims 6 to 8, characterized in that: The specific execution process of the measurement and correction module is as follows: Step 5.1: Construct the posterior estimation ellipsoid set of the state variables at time k+1, as follows: in: is x k+1 An estimate of , representing the ellipsoid center, For a given real positive definite matrix, represents the ellipsoid shape matrix; Step 5.2: Construct the measurement update ellipsoid set according to the measurement equation, as follows: Where: Z k+1 Update the ellipsoid for the measurement at time k+1; Step 5.3: Construct the measurement update ellipsoid set according to the measurement equation, as follows: Where: Z k+1 Update the ellipsoid for the measurement at time k+1; Step 5.4: Build a satisfying The minimum ellipsoid is as follows: Where: for any 0≤ρ k+1 ≤1, the above formula is valid; Step 5.5, rewrite the above formula: Where: 0≤δ k+1 ≤1, is an intermediate variable, if The measurement correction ellipsoid EL can be obtained k+1 ; Step 5.6: Solve the center and shape matrix of the ellipsoid set of the state variables at time k+1 according to the measurement correction iterative formula, as follows: Where: parameter ρ k+1 The optimal value of is obtained by solving formula (29): Step 5.7: Determine whether the state estimation calculation is stopped. If k <k max , then set k = k + 1 and jump to step 3 to continue the calculation; otherwise, end the calculation.

Images