Power system adaptive robust dynamic state estimation method based on hybrid measurement

Abstract

The invention provides an adaptive robust dynamic state estimation method for a power system based on hybrid measurement, which comprises the following steps of: firstly, realizing fusion of hybrid measurement through a measurement transformation technology, and constructing pseudo measurement based on space-time characteristics of the system to form a hybrid measurement set; furthermore, an adaptive forgetting factor is introduced, and a noise estimation algorithm is improved, so that the time-varying system noise can be estimated more quickly and accurately, and the dynamic tracking capability of the algorithm is improved; secondly, estimating measurement noise by analyzing characteristics of the measurement device, and constructing a robust factor based on standard innovation to suppress bad data; and finally, dynamic state estimation is carried out based on an extended Kalman filtering criterion, and the accurate operation state of the system can be obtained under the influence of unknown model parameters, system state change and bad data.

Description

Technical Field

[0001] The invention relates to the technical field of dispatching automation of power systems containing renewable energy, in particular to a hybrid measurement-based adaptive robust dynamic state estimation method for power systems. Background Art

[0002] State estimation is one of the core functions of the energy management system of the power grid dispatching and control center. By analyzing the collected measurement data, the operation status of the entire power network is calculated and completed, and real-time monitoring of the power system is realized. The development of synchronized phasor measurement units (PMUs) and the construction of wide area measurement systems (WAMSs) provide highly synchronized, high-precision, and high-refresh-frequency measurement data, providing a new way to quickly and accurately obtain the real-time status of the power system. However, due to economic costs and configuration environment, PMUs cannot be widely configured in power grids at present. In order to meet the observability of the system, it is necessary to combine the original data acquisition and monitoring systems (SCADA) and smart electric energy meters in the power grid to form hybrid measurement to improve measurement redundancy. Therefore, it is more practical to study the state estimation of the power system under hybrid measurement. At present, there are large differences between WAMS and SCADA and other measurement systems in terms of data composition, sampling time, refresh frequency, and data accuracy. In order to improve the consistency of measurement information and the accuracy of state estimation, multi-source measurement data fusion is required.

[0003] Relying on the high refresh rate of PMU can realize the real-time state perception of the power system, but pseudo-measurements need to be added only at the PMU sampling time to ensure the rapid global estimation of the system. Linear interpolation is a common pseudo-measurement construction method, but the accuracy of pseudo-measurements is subject to the accuracy of the predicted value, and when the system state fluctuates greatly, pseudo-measurements based on interpolation cannot reflect the dynamic changes of the system. Therefore, the construction of pseudo-measurements needs to make full use of the correlation of measurement data in space and time.

[0004] The state estimation of the power system establishes the relationship between the measurement and the state to be determined through the measurement equation, and solves the operating state of the current time section. It can be divided into two methods: static state estimation and dynamic state estimation. The current application of static state estimation is relatively mature, with the least square method as the main solution method. However, this method only uses the measurement data of a single time section, ignoring the dynamic changes and correlation of the system state at continuous moments. With the development of smart grids, the power system presents more obvious dynamic characteristics. Dynamic state estimation generates state prediction values ​​based on the measurement data of multiple time sections, and estimates the system operating state in combination with the latest measurement values. The results obtained can better reflect the dynamic changes of the system. At present, the dynamic state estimation algorithm is mainly based on the Kalman Filter (KF) design. For nonlinear systems such as power systems, extended Kalman Filter (EKF), Unscented Kalman Filter (UKF) and Cubature Kalman Filter (CKF) algorithms are extended. EKF linearizes nonlinear functions by performing Taylor expansion on nonlinear functions and omitting high-order terms, but it also introduces linearization errors. When the system nonlinearity is strong, it will produce large estimation errors. KF and CKF obtain the statistical characteristics after nonlinear transformation by point set sampling, avoiding linearization errors, but they are not as efficient as EKF in terms of computation. At the same time, the above three algorithms are based on the KF framework and have high requirements on the accuracy of the model. Only when the statistical characteristics of system noise are mastered can good estimation effects be achieved. For power systems with model uncertainty and unknown statistical characteristics of system noise, some literature introduces a forgetting factor based on the Sage-Husa noise estimator, and adjusts the process noise covariance in real time according to the statistical characteristics of the new information to adapt to the dynamic changes of the system. However, the above methods are affected by the value of the forgetting factor. If the degree of forgetting is too large, the estimation will be inaccurate due to system disturbances. If the degree of forgetting is too small, the estimation error will accumulate due to the influence of stale data in the steady state.

[0005] In view of this, it becomes more important to quickly and accurately obtain the real-time state of the entire system in a new power system with greatly enhanced randomness and volatility. In order to reduce the impact of unknown measurement noise and accurately obtain the system state, the present invention innovatively proposes a hybrid measurement state estimation method based on an adaptive robust extended Kalman filter algorithm. The present invention proposes an adaptive robust dynamic state estimation method for power systems based on hybrid measurements to solve the above problems. Summary of the invention

[0006] The purpose of the present invention is to propose an adaptive and robust dynamic state estimation method for power systems based on hybrid measurements; first, the measurement transformation technology is used to linearize the measurement equation, and pseudo-measurements are constructed based on the spatiotemporal correlation between node states to achieve the fusion of hybrid measurements; further, a forgetting factor adaptive method is designed to improve the Sage-Husa noise estimator, and the filtering effect is determined according to the new information and the degree of forgetting is adjusted, which can improve the accuracy of the adaptive estimation of process noise, thereby improving the algorithm's ability to dynamically track the system state. In addition, the measurement noise is estimated by analyzing the statistical characteristics of the measurement device, and an equivalent weight function is constructed according to the standard new information to correct the measurement noise covariance to improve the algorithm's robustness. Finally, dynamic state estimation is performed based on the extended Kalman filter criterion to accurately obtain the operating state of the system.

[0007] To achieve the above object, the technical solution of the present invention is: a method for adaptive robust dynamic state estimation of a power system based on hybrid measurement, comprising the following steps:

[0008] Step S1, using measurement transformation technology to linearize the measurement equation, and constructing pseudo-measurement based on the spatiotemporal correlation between node states to achieve fusion of hybrid measurements;

[0009] Step S2: Determine the filtering effect according to the statistical characteristics of the new information and assign weights to the new data, and adjust the process noise covariance in real time to adapt to the dynamic changes of the system;

[0010] Step S3, estimating the measurement noise by analyzing the statistical characteristics of the measurement device, constructing an equivalent weight function according to the standard innovation, and correcting the measurement noise covariance to improve the algorithm robustness;

[0011] Step S4: Based on the mixed measurement set of step S1, an extended Carr filter algorithm is used to adjust the noise covariance according to steps S2 and S3, perform dynamic state estimation, and obtain the system operation state.

[0012] Preferably, the fusion of hybrid measurements in step S1 is specifically implemented as follows: when SCADA measurements are updated, hybrid measurements are used for state estimation; when non-SCADA measurements are updated, PMU measurements are used in combination with pseudo measurements for state estimation;

[0013] The hybrid measurement is to integrate PMU and SCADA to form a hybrid measurement set, specifically:

[0014] The node voltage phasor and branch current phasor are directly measured by PMU, and the node voltage phasor is used as the state quantity to obtain the linear measurement equation based on PMU:

[0015]

[0016] Where U is the node voltage state quantity, Up and I p are the vectors constructed from the PMU node voltage phasor and the branch current phasor respectively; the superscript p represents the PMU measurement; E p is the identity matrix, Y p is the node-branch admittance matrix, ν U and ν I are the measurement errors of voltage and current respectively;

[0017] The SCADA data is transformed into branch power measurement and node injection power measurement into equivalent branch current phasor measurement and equivalent node injection current phasor measurement. The conversion formula is as follows:

[0018]

[0019] In the formula, and are the equivalent branch current phasor and the equivalent node injection current phasor respectively; P ft , Q ft , P f , Q f The active power and reactive power of branch ft and the injected active power and injected reactive power of node f respectively; U re,f and U im,f are the real and imaginary parts of the voltage at node f, respectively;

[0020] After measurement transformation, the linear measurement equation based on SCADA is obtained:

[0021]

[0022] In the formula, and They are the vectors constructed from the equivalent branch current phasor measurement and the equivalent node injection current phasor measurement obtained by measurement transformation of SCADA data; the superscript s represents SCADA measurement; and are the node-branch admittance matrix and the node admittance matrix, respectively. and They are the measurement errors of node-branch equivalent current and node equivalent current respectively;

[0023] By integrating the measurement equations of PMU and SCADA, the measurement equation based on hybrid measurement is finally obtained:

[0024]

[0025] Where x is the node voltage vector; z p and z s are the measurement vectors of PMU and SCADA respectively; H pand H s Respectively represent the coefficient matrices between PMU and SCADA and the node voltage vector; ν p and ν s represent the measurement noise vectors of PMU and SCADA respectively;

[0026] The pseudo measurement is derived from the PMU measurement of the PMU unobservable node at the non-SCADA update time according to the data space correlation and time correlation, specifically:

[0027] The pseudo-measurement constructed based on the power transfer distribution factor and the pseudo-measurement constructed based on the time correlation are subtracted from the state estimate, and the weight is set using an exponential function. The weighted sum is used to form the pseudo-measurement based on time-space correlation:

[0028]

[0029] In the formula, It is a pseudo-measurement constructed based on spatiotemporal correlation; and They are respectively a pseudo measurement constructed based on the power transfer distribution factor and a pseudo measurement constructed based on the time correlation at time k; and are the corresponding weights respectively; Represents the estimated value of the state at time k-1.

[0030] Preferably, the pseudo measurement constructed based on the power transfer distribution factor is to calculate the indirect measurement of the PMU directly observed node through the PMU measurement of the PMU directly observed node, and construct the power transfer distribution factor based on the indirect measurement to generate the power pseudo measurement of the unobservable branch, specifically as follows:

[0031] Assume that node i is equipped with a PMU and node j is not equipped with a PMU. Then the voltage and current of node j are indirectly calculated through the measurement of node i:

[0032]

[0033] In the formula, Represent the voltages of nodes i and j respectively, They represent the current flowing from node i to node j and the current flowing from node j to node i, respectively. ii , Y ij , Y ji , Y jj Represents the elements corresponding to the node-branch admittance matrix:

[0034]

[0035] Where U, θ are the amplitude and phase angle of the voltage, P, Q are the injected active and reactive powers of the node or branch, G, B are the conductance and susceptance of the node or branch, and the subscripts i, j, ij represent the corresponding node i, node j, and branch ij;

[0036] The change of node voltage is derived when the node injection power changes:

[0037]

[0038] In the formula, ΔU j , Δθ j , ΔP j , ΔQ j They represent the voltage amplitude change, voltage phase angle change, injected active power change and injected reactive power change of node j respectively, and E j 、F j represents the intermediate variable corresponding to node j, E j =(G j +B j G j -1 B j ) -1 , F j =-(G j +B j G j -1 B j ) -1 B j G j -1 ;

[0039] The power change of the PMU unobservable node injected is derived from the power change of the branch indirectly measured by the PMU:

[0040]

[0041] In the formula, ΔP ji , ΔQ ji They represent the active power change and reactive power change of the indirect measurement branch ji respectively. The variable M represents the relationship between the branch power measurement change and the node injection power change. It is composed of branch admittance. The subscript represents the corresponding node and branch, that is, the variable M j-ji It represents the relationship between the measured power change of branch ji and the injected power change of node j;

[0042] By calculating the injected power variation of the PMU unobservable node, the power variation of the PMU unobservable branch is derived:

[0043]

[0044] In the formula, ΔP jl , ΔQ jl They represent the active power change and reactive power change of the unobservable branch jl, respectively. The variable M j-jl It represents the relationship between the power measurement change of branch jl and the power injection change of node j;

[0045] Finally, taking the power measurement of branch jl measured by SCADA at time k as the starting point, the power change of branch ji from time k to time k+1 is calculated based on the PMU measurement. Combining formula (9), the power change of branch jl that cannot be observed by PMU is obtained, and the power pseudo-measurement of branch jl at time k+1 is generated:

[0046]

[0047] In the formula, They represent the active power spatial pseudo-measurement and reactive power spatial pseudo-measurement of branch jl at time k+1 respectively; if time k is the SCADA update time, then is the active power measurement and reactive power measurement of branch jl measured by SCADA at time k; if time k is not a SCADA update time, then It refers to the spatial pseudo-measurement of active power and reactive power of branch jl at time k.

[0048] Preferably, the pseudo measurement constructed based on time correlation is designed based on the change relationship of each injection power of the same node within a period of time, and is specifically as follows:

[0049]

[0050] In the formula, represents the pseudo-measurement of branch jl power at time k+1; if time k is the SCADA update time, then is the power measurement of branch jl measured by SCADA at time k; if time k is not a SCADA update time, then It refers to the power time pseudo measurement of branch jl at time k. represents the indirect power measurement of branch ji by PMU at time k+1 and time k; b ji is the trend coefficient, is the power transformation ratio of branch ji and branch jl at the last two SCADA update times; S ji is the correlation coefficient, which is the correlation coefficient between branch ji and branch jl within a time window L; t s is the SCADA measurement cycle; T s is the time of the most recent SCADA update; n is the number of PMU direct observation nodes connected to node j, Indicates SCADA in T s Moment and T s -The power of branch ji at time ts, Indicates SCADA in T s Moment and T s -The power of branch jl at time ts; and They represent the power data of branch ji and branch jl in the time window L of SCADA respectively.

[0051] Preferably, step S2 is based on adaptively adjusting the forgetting factor and introducing a process noise adaptive estimation of a gradually fading memory exponential weighting method:

[0052] The process noise is the difference between the true state value and the predicted state value. The estimated value of the state is used to approximate the process noise:

[0053]

[0054] Where, the subscript k is used to represent the kth moment, ε k is the estimation error, ω k is the system process noise, x k is the system state value, f(x k ) is the state transfer equation, which reflects the dynamic changes of the distribution network operation state. is the estimated value of the state, is the Jacobian matrix of the state transfer equation, K k+1 is the Kalman gain, v k+1 For new information

[0055] When the number of sampling moments in the sample space is N, the unbiased estimate of the process noise covariance is:

[0056]

[0057] In the formula, Q k is the process noise ω k The covariance matrix of is the mean of N ω samples, P is the estimation error covariance,

[0058] In these N samplings, the weight factor d is introduced to give a larger weight to the latest sampled data. The weight factor d at time k is k for:

[0059]

[0060] Where b k is the forgetting factor. The greater the system process noise changes, the greater the k The smaller the value of ;

[0061] Modify the process noise covariance to:

[0062]

[0063] The adaptive adjustment of the forgetting factor is specifically as follows:

[0064] The intermediate variable η is designed by the ratio of the innovation variance and the innovation covariance matrix trace. k Measure the overall filtering effect of the current state estimate:

[0065]

[0066] Adaptively adjust the forgetting factor according to the filtering effect:

[0067]

[0068] In the formula, R k is the covariance matrix of the measurement noise, H k is the Jacobian matrix of the measurement equation, is the prediction error covariance matrix, b1 and b2 are the lower limit of the forgetting factor and the adaptive dynamic range respectively, and c1 and c2 are parameters to measure the filtering effect.

[0069] Preferably, the lower limit of the forgetting factor b1 is 0.6, b1+b2 is 0.95-0.995; c1=2, c2=3; when η k ≤2, that is, when the sum of the new standard deviation is less than 2 times the theoretical standard deviation, the filtering effect is good and the forgetting factor b is reduced. k To improve the degree of forgetting, the process noise can quickly track the accurate value; when 2<η k ≤3, the filtering effect is judged to be average, and the forgetting factor b is maintained k Stable, making the process noise stable estimation; when η k >3, indicating that the filtering process is disturbed, the forgetting factor b k If it is 1, the process noise is not updated.

[0070] Preferably, estimating the measurement noise by analyzing the statistical characteristics of the measurement device specifically includes: adaptively estimating the measurement noise according to the latest measurement before each state estimation:

[0071]

[0072] In the formula, R(i,i) represents the i-th element of the diagonal of R, z p 、z s Respectively represent the corresponding PMU measurement and SCADA measurement; σ m and σ a They represent the amplitude error e measured by PMU respectively. mand phase error e a The standard deviation of , t is the intermediate variable, and σ P and σ Q are the standard deviations of SCADA active and reactive power measurement errors, U re and U im are the real and imaginary parts of the node voltage.

[0073] Preferably, the equivalent weight function is constructed according to the standard innovation information to correct the measurement noise covariance as follows:

[0074] Standard new information is used to identify and process bad data. The expression of the standard new information is:

[0075]

[0076] In the formula, v k For the new news, Indicates standard news.

[0077] Based on the standard innovation, a robustness factor is introduced to correct the measurement noise. The robustness factor is generated using the IGGⅢ weight function:

[0078]

[0079] In the formula, Represents the robustness factor The i-th element of express The i-th element of , k0 and k1 are the upper and lower thresholds, which are parameters to measure the quality of the measurement value;

[0080] Get the corrected measurement noise covariance matrix for:

[0081]

[0082] In the formula, R k represents the measurement noise covariance matrix at the kth moment, Represents the robustness factor.

[0083] Preferably, the parameters for measuring the quality of the measurement values ​​are k0=3, k1=4. When the measured data is greater than 3 times the standard deviation, it is judged as bad data and subjected to anti-error processing to reduce its weight in state estimation; when the measured data is greater than 4 times the standard deviation, the bad data is eliminated.

[0084] Preferably, the Kalman filter in step S4 is divided into two processes: state prediction and state filtering. The specific process is shown as follows:

[0085] State prediction process:

[0086]

[0087] State filtering process:

[0088]

[0089] Where, the subscript k is used to represent the kth moment, is the state prediction value; is the estimated value of the state; It represents the state transfer equation, reflecting the dynamic changes of the distribution network operation state. Represents the system measurement equation, reflecting the relationship between state quantity and quantity measurement, is the prediction error covariance matrix; is the estimated error covariance matrix; I represents the identity matrix, K k is the Kalman gain; Q k and are the process noise ω adjusted and corrected in step S2 and step S3 respectively k and measurement noise ν k The covariance matrix of is the Jacobian matrix of the measurement equation, which is obtained according to formula (4) under mixed measurement; is the Jacobian matrix of the state transfer equation, which is obtained according to the Holt exponential smoothing method:

[0090]

[0091] In the formula, is the state prediction value, is the state estimation value; a and b are the horizontal component and vertical component of Holt exponential smoothing respectively; α and β are the smoothing coefficients.

[0092] Compared with the prior art, the present invention has the following beneficial effects:

[0093] The invention is based on the hybrid measurement of the power system adaptive robust dynamic state estimation method. The proposed pseudo measurement construction method and the adaptive robust state estimation method can effectively perform global state estimation when the PMU configuration does not meet the system observability, and can accurately obtain the system operation state under the influence of unknown model parameters, system state changes and bad data. It has high application value. BRIEF DESCRIPTION OF THE DRAWINGS

[0094] Figure 1 This is a schematic diagram of the inter-node power flow in step S12 of the present invention;

[0095] Figure 2 The hybrid measurement state estimation architecture of step S12 of the present invention;

[0096] Figure 3 It is the overall flow chart of the present invention;

[0097] Figure 4 The pseudo-measurement result of the voltage of the unobservable node 31 in the simulation verification of the present invention is compared;

[0098] Figure 5 The comparison of the state estimation results of node 3 during the transient period in the simulation verification of the present invention;

[0099] Figure 6 The figure is a comparison of the state estimation results of node 3 under bad data in the simulation verification of the present invention. DETAILED DESCRIPTION

[0100] The following is combined with Figure 1-6 , the technical solution of the present invention is specifically described.

[0101] The present invention proposes a method for estimating the dynamic state of an electric power system based on hybrid measurement. Figure 3 , including the following steps:

[0102] Step S1: Use measurement transformation technology to linearize the measurement equation, and construct pseudo-measurement based on the spatiotemporal correlation between node states to achieve the fusion of hybrid measurements;

[0103] Step S2: Determine the filtering effect according to the statistical characteristics of the new information and assign weights to the new data, and adjust the process noise covariance in real time to adapt to the dynamic changes of the system;

[0104] Step S3: Estimate the measurement noise by analyzing the statistical characteristics of the measurement device, construct an equivalent weight function based on the standard innovation, and correct the measurement noise covariance to improve the algorithm's robustness;

[0105] Step S4: Based on the mixed measurement set of step S1, the extended Carr filter algorithm is used to adjust the noise covariance according to steps S2 and S3, perform dynamic state estimation, and obtain the system operation state.

[0106] Furthermore, step S1 is specifically as follows:

[0107] Although PMU measurement data has the characteristics of high precision and high refresh rate, it is currently difficult to ensure system observability by relying solely on PMU measurement due to the configuration cost and technical requirements of the device. In order to meet actual application needs, it is necessary to integrate PMU measurement and traditional SCADA measurement data to form hybrid measurement.

[0108] 1. Integrate PMU and SCADA to form a hybrid measurement set

[0109] Since PMU can directly measure the node voltage phasor and branch current phasor, the linear measurement equation based on PMU can be obtained by taking the node voltage phasor as the state quantity:

[0110]

[0111] Where U is the node voltage state quantity, U p and I p are the vectors constructed from the PMU node voltage phasor and the branch current phasor respectively; the superscript p represents the PMU measurement; E p is the identity matrix, Y p is the node-branch admittance matrix, ν U and ν I are the measurement errors of voltage and current respectively.

[0112] After integrating SCADA measurement data, power measurement will be introduced, resulting in nonlinear measurement equations. In order to improve calculation efficiency, it is necessary to perform measurement transformation on SCADA data, converting branch power measurement and node injection power measurement into equivalent branch current phasor measurement and equivalent node injection current phasor measurement. The conversion formula is as follows:

[0113]

[0114] In the formula, and are the equivalent branch current phasor and the equivalent node injection current phasor respectively; P ft , Q ft , P f , Q f The active power and reactive power of branch ft and the injected active power and injected reactive power of node f respectively; U re,f and U im,f are the real and imaginary parts of the voltage at node f, respectively.

[0115] After measurement transformation, the linear measurement equation based on SCADA is obtained:

[0116]

[0117] In the formula, and They are the vectors constructed from the equivalent branch current phasor measurement and the equivalent node injection current phasor measurement obtained by measurement transformation of SCADA data; the superscript s represents SCADA measurement; and are the node-branch admittance matrix and the node admittance matrix, respectively. and They are the measurement errors of node-branch equivalent current and node equivalent current respectively;

[0118] By integrating the measurement equations of PMU and SCADA, the measurement equation based on hybrid measurement is finally obtained:

[0119]

[0120] Where x is the node voltage vector; z p and z s are the measurement vectors of PMU and SCADA respectively; H p and H s Respectively represent the coefficient matrices between PMU and SCADA and the node voltage vector; ν p and ν s Represent the measurement noise vectors of PMU and SCADA respectively.

[0121] 2. Constructing pseudo-measurements based on spatiotemporal correlations:

[0122] When the PMU configuration does not meet the system observability, dynamic state estimation based only on PMU measurements cannot reflect the global state of the system. Although the hybrid measurement formed by combining SCADA data can meet the system observability, due to the slow update frequency of SCADA, in order to quickly track the dynamic changes of the power system, it is still necessary to construct pseudo measurements for PMU unobservable nodes at non-SCADA update times to support the global state estimation of the system, where observable nodes include direct observation nodes and indirect observation nodes.

[0123] Figure 1 A schematic diagram of data flow between nodes is given. According to the spatial and temporal correlation of data, pseudo-measurements of PMU unobservable nodes are derived from PMU measurements.

[0124] Assuming that node i is equipped with a PMU and node j is not equipped with a PMU, the voltage and current of node j can be indirectly calculated from the measurement of node i:

[0125]

[0126] In the formula, Represent the voltages of nodes i and j respectively, They represent the current flowing from node i to node j and the current flowing from node j to node i, respectively. ii , Y ij , Y ji , Y jj Represents the element corresponding to the node-branch admittance matrix.

[0127] Based on indirect measurement, a power transfer distribution factor is constructed to generate pseudo-measurement of power of unobservable branches. According to the DC power flow calculation formula, the relationship between branch power and node injection power and node voltage can be obtained:

[0128]

[0129] Where U, θ are the amplitude and phase angle of the voltage, P, Q are the injected active and reactive powers of the node or branch, G, B are the conductance and susceptance of the node or branch, and the subscripts i, j, ij represent the corresponding node i, node j, and branch ij.

[0130] When the node injection power changes, the change of node voltage can be derived:

[0131]

[0132] In the formula, ΔU j , Δθ j , ΔP j , ΔQ j They represent the voltage amplitude change, voltage phase angle change, injected active power change and injected reactive power change of node j respectively, and E j 、F j Represents the intermediate variable corresponding to node j. The relationship between power and voltage is established through variables E and F. j =(G j +B j G j -1 B j ) -1 , F j =-(G j +B j G j -1 B j ) -1 B j G j -1 .

[0133] According to equations (6) and (7), the power change of the branch power indirectly measured by the PMU can be deduced as follows:

[0134]

[0135] In the formula, ΔP ji , ΔQ ji They represent the active power change and reactive power change of the indirect measurement branch ji respectively. The variable M represents the relationship between the branch power measurement change and the node injection power change. It is composed of branch admittance. The subscript represents the corresponding node and branch, that is, the variable M j-ji It represents the relationship between the power measurement change of branch ji and the power injection change of node j.

[0136] Similarly, by calculating the injected power change of the PMU unobservable node, the power change of the PMU unobservable branch can be derived:

[0137]

[0138] Where ΔP jl , ΔQ jl They represent the active power change and reactive power change of the unobservable branch jl, respectively. The variable M j-jl It represents the relationship between the power measurement change of branch jl and the power injection change of node j.

[0139] Finally, taking the power measurement of branch jl measured by SCADA at time k as the starting point, the power change of branch ji from time k to time k+1 is calculated according to the PMU measurement. Combining formula (9), the power change of branch jl that cannot be observed by PMU can be obtained, thereby generating the power pseudo-measurement of branch jl at time k+1.

[0140]

[0141] In the formula, They represent the active power spatial pseudo-measurement and reactive power spatial pseudo-measurement of branch jl at time k+1 respectively; if time k is the SCADA update time, then is the active power measurement and reactive power measurement of branch jl measured by SCADA at time k; if time k is not a SCADA update time, then It refers to the spatial pseudo-measurement of active power and reactive power of branch jl at time k.

[0142] The obtained PMU unobservable branch power pseudo-measurement is integrated with PMU measurement for state estimation, which can achieve full observability of the power system.

[0143] Since DC power flow is a simplified power network model that ignores line resistance and parallel branches, it cannot provide the same accuracy and detail as the AC power flow model. Therefore, as the number of estimates increases, there may be a problem of estimation divergence. Therefore, a pseudo-measurement construction method based on time correlation is introduced. Many components in the power system, represented by generators, are dynamic components. These components respond according to inputs and historical states, making the current state of the power system closely related to the previous state, and the system state has a strong time correlation. Similarly, there is an obvious time correlation between the various injected powers of the same node. According to the changing relationship between the various injected powers of the same node over a period of time, a pseudo-measurement construction method based on time correlation is designed, which is expressed as follows:

[0144]

[0145]

[0146] In the formula, represents the pseudo-measurement of branch jl power at time k+1; if time k is the SCADA update time, then is the power measurement of branch jl measured by SCADA at time k; if time k is not a SCADA update time, then It refers to the power time pseudo measurement of branch jl at time k. represents the indirect power measurement of branch ji by PMU at time k+1 and time k; b ji is the trend coefficient, is the power transformation ratio of branch ji and branch jl at the last two SCADA update times; S ji is the correlation coefficient, which is the correlation coefficient between branch ji and branch jl within a time window L; t s is the SCADA measurement cycle; T s is the time of the most recent SCADA update; n is the number of PMU direct observation nodes connected to node j, Indicates SCADA in T s Moment and T s -The power of branch ji at time ts, Indicates SCADA in T s Moment and T s -The power of branch jl at time ts; and They represent the power data of branch ji and branch jl in the time window L of SCADA respectively.

[0147] The pseudo-measurement constructed based on the power transfer distribution factor and the pseudo-measurement constructed based on the time correlation are subtracted from the state estimate, and the weight is set using an exponential function. The weighted sum is used to form the pseudo-measurement based on time-space correlation:

[0148]

[0149] In the formula, It is a pseudo-measurement constructed based on spatiotemporal correlation; and They are respectively a pseudo measurement constructed based on the power transfer distribution factor and a pseudo measurement constructed based on the time correlation at time k; and are the corresponding weights respectively; Represents the estimated value of the state at time k-1.

[0150] Combining the two pseudo-measurements, based on time correlation and spatial structure relationship respectively, and combining them with the dynamic state estimation results can suppress the problem of pseudo-measurement estimation divergence when the number of state estimations increases, and further improve the estimation accuracy of pseudo-measurements.

[0151] State estimation architecture based on hybrid measurements Figure 2 As shown, t0 is the first SCADA sampling time, t p is the PMU measurement period, t s is the SCADA measurement cycle. When SCADA measurements are updated, hybrid measurements are used for state estimation. When non-SCADA measurements are updated, PMU measurements combined with pseudo measurements are used for state estimation.

[0152] Furthermore, step S2 is specifically as follows:

[0153] 1. Adaptive estimation of process noise using the exponential weighting method with fading memory

[0154] The Sage-Husa noise estimation uses the latest value of the new information to adjust the process noise covariance sequence. By introducing the fading memory exponential weighting method, the new data is weighted according to the degree of system change, which can better reflect the real-time changes of the system process noise.

[0155] The process noise is the difference between the true state value and the predicted state value. Since the true state of the system is unknown, the estimated value of the state is used to approximate the process noise:

[0156]

[0157] Where, the subscript k is used to represent the kth moment, ε k is the estimation error, ω k is the system process noise, x k is the system state value, f(x k ) is the state transfer equation, which reflects the dynamic changes of the distribution network operation state. is the estimated value of the state, is the Jacobian matrix of the state transfer equation, K k+1 is the Kalman gain, v k+1 For new breath;

[0158] When the number of sampling moments in the sample space is N, the unbiased estimate of the process noise covariance is:

[0159]

[0160] In the formula, Q k is the process noise ω k The covariance matrix of is the mean of N ω samples, and P is the estimation error covariance.

[0161] In these N samplings, the latest sampling data can better reflect the time-varying characteristics of the current process noise, so the weight factor d is introduced to give a larger weight to the latest sampling data. The weight factor d at time k is k for:

[0162]

[0163] Where b k is the forgetting factor. The greater the system process noise changes, the greater the k The smaller the value of .

[0164] Therefore, the process noise covariance can be modified as:

[0165]

[0166] 2. Adaptive adjustment of forgetting factor

[0167] In order to better adapt to the dynamic changes of time-varying noise systems, a forgetting factor adaptive adjustment method is designed to improve the algorithm's ability to estimate process noise under various operating conditions.

[0168] The intermediate variable η is designed by the ratio of the innovation variance and the innovation covariance matrix trace. k Measure the overall filtering effect of the current state estimate:

[0169]

[0170] Adaptively adjust the forgetting factor according to the filtering effect:

[0171]

[0172] In the formula, R k is the covariance matrix of the measurement noise, H k is the Jacobian matrix of the measurement equation, is the prediction error covariance matrix, b1 and b2 are the lower limit of the forgetting factor and the adaptive dynamic range, respectively. In order to quickly adapt to the changes in the statistical characteristics of process noise, the lower limit of the forgetting factor b1 can be 0.6, and b1+b2 can be 0.95-0.995; c1 and c2 are parameters to measure the quality of the filtering effect. In this paper, c1=2 and c2=3 are taken. When η k ≤2, that is, when the sum of the new standard deviation is less than 2 times the theoretical standard deviation, the filtering effect is good and the forgetting factor b is reduced. k To improve the degree of forgetting, the process noise can quickly track the accurate value; when 2<η k ≤3, the filtering effect is judged to be average, and the forgetting factor b is maintained k Stable, making the process noise stable estimation; when η k>3, indicating that the filtering process is disturbed, the forgetting factor b k If it is 1, the process noise is not updated. At the same time, if the system changes suddenly, the prediction value generated based on the stable system will produce a large error, and the process noise covariance should be increased to reduce the weight of the prediction value in the Kalman filter.

[0173] Furthermore, step S3 is specifically as follows:

[0174] 1. Estimating measurement noise by analyzing the statistical characteristics of the measurement device

[0175] The amplitude and phase of PMU measurements are usually assumed to follow independent zero-mean Gaussian distributions. is an element measured by the PMU, and its amplitude and phase measurement errors are e m and e a , according to the Euclidean distance, The measurement error It can be expressed as:

[0176]

[0177] In the formula, express The true value of z T express The amplitude of, j represents the imaginary part, θT represents The phase angle.

[0178] Assume that the amplitude error e measured by the PMU is m and phase error e a The mean is μ m and μ a , with standard deviation σ m and σ a The phasor error is The mean value of can be expressed as:

[0179]

[0180] In the formula, About cose a Performing Taylor expansion and ignoring higher-order terms yields So the phasor error The variance of is expressed as:

[0181]

[0182] Because the amplitude and phase of PMU measurements are usually assumed to follow independent zero-mean Gaussian distributions, μ m and μ a are all equal to 0, and e m and e aIndependent of each other, phasor error The variance can be further simplified as:

[0183]

[0184] Since the PMU amplitude error is generally a few thousandths of the true amplitude value, after identifying and processing the bad PMU data, the phasor amplitude measurement value z can be used to replace the unknown phasor amplitude true value z T In the known m and e a In the case of statistical properties, the value of t can be approximately estimated by mathematical statistics methods.

[0185] After the power measurement of SCADA is transformed, the measurement variance also changes. According to the error propagation theory, the measurement variance of the equivalent current phasor can be obtained: for:

[0186]

[0187] In the formula, σ I is the standard deviation of the equivalent current measurement error, σ P and σ Q are the standard deviations of SCADA active and reactive power measurement errors, respectively.

[0188] The measurement noise covariance matrix R is set by the diagonal matrix of the measurement error variance, and the measurement noise can be adaptively estimated based on the latest measurement before each state estimation:

[0189]

[0190] In the formula, R(i,i) represents the i-th element of the diagonal of R, z p 、z s They represent the corresponding PMU measurements and SCADA measurements respectively;

[0191] 2. Construct an equivalent weight function based on standard innovation and correct the measurement noise covariance

[0192] Considering that bad data is inevitable in practical applications, in order to avoid the adverse impact of bad data on state estimation, standard new information is used to identify and process bad data. The new information is the difference between the measured value and the predicted value:

[0193]

[0194] In the formula, v k For new information, z k is the measured value, is the state prediction value;

[0195] The theoretical value of the innovation covariance can be expressed as:

[0196]

[0197] Normalize the news to get the standard news:

[0198]

[0199] In the formula, Indicates standard news.

[0200] Based on the standard new information, the robustness factor is introduced to correct the measurement noise. The robustness factor is generated using the IGGⅢ weight function:

[0201]

[0202] In the formula, Represents the robustness factor The i-th element of express The i-th element of , k0 and k1 are the upper and lower limit thresholds, which are parameters to measure the quality of the measurement value. In order to ensure the suppression of bad data under measurement redundancy, this paper takes k0=3, k1=4, that is, when the measurement data is greater than 3 times the standard deviation, it is judged as bad data, and anti-error processing is performed to reduce its weight in state estimation; when the measurement data is greater than 4 times the standard deviation, the bad data is eliminated.

[0203] Compared with the standard residual, the standard innovation takes into account the predicted state and avoids residual contamination. Finally, the corrected measurement noise covariance matrix is ​​obtained for:

[0204]

[0205] In the formula, R k represents the measurement noise covariance matrix at the kth moment.

[0206] Furthermore, step S4 is specifically as follows:

[0207] Dynamic state estimation based on Kalman filtering

[0208] When performing dynamic state estimation on nonlinear power systems, the extended Kalman filter algorithm is usually used. The extended Kalman filter is based on the dynamic state estimation model. According to the statistical characteristics of noise, it follows the linear unbiased minimum mean square error criterion to filter out random interference in the measurement information and obtain the best estimate of the state variable. Kalman filtering can be divided into two processes: state prediction and state filtering. The specific process is shown as follows:

[0209] State prediction process:

[0210]

[0211] State filtering process:

[0212]

[0213] Where, the subscript k is used to represent the kth moment, is the state prediction value; is the state estimation value; f(.) represents the state transfer equation, which reflects the dynamic changes of the distribution network operation state. Represents the system measurement equation, reflecting the relationship between state quantity and quantity measurement, is the prediction error covariance matrix; is the estimated error covariance matrix; I represents the identity matrix, K k is the Kalman gain; Q k and are the process noise ω adjusted and corrected in step S2 and step S3 respectively k and measurement noise ν k The covariance matrix of is the Jacobian matrix of the measurement equation, which can be obtained according to formula (4) under mixed measurement; is the Jacobian matrix of the state transfer equation, which can be obtained according to the Holt exponential smoothing method:

[0214]

[0215] In the formula, is the state prediction value, is the state estimation value; a and b are the horizontal component and vertical component of Holt exponential smoothing respectively; α and β are the smoothing coefficients.

[0216] Based on the above model and process design, the present invention uses the IEEE39 node system for simulation verification.

[0217] The proposed algorithm is implemented and run in MATLAB 2020a. The true value is obtained by performing power flow calculation on the simulation case, and the measured value is obtained by adding the measurement error of zero-expectation Gaussian distribution to the true value. The standard deviation of the measurement error simulated in this paper is shown in Table 1. PMU samples at intervals of 0.01s, SCADA samples at intervals of 1s, and the interval of each dynamic state estimation is 0.01s.

[0218] Table 1 Standard deviation of measurement error

[0219]

[0220] In order to measure the accuracy of the estimation results, the mean absolute error (MAE), root mean square error (RMSE) and mean vector error (MVE) are introduced as indicators for evaluating the estimation accuracy of the algorithm. The definitions of the three evaluation indicators are as follows:

[0221]

[0222]

[0223] In the formula, x i is the i-th component of the true value of the state quantity; is the i-th component of the estimated value of the state quantity; is the i-th component of the state quantity measurement value; n is the dimension of the state quantity.

[0224] (1) Pseudo-measurement structure simulation verification

[0225] PMUs are configured at nodes 3, 8, 11, 19, 23, and 26. At this time, the system is not completely observable. According to step S12, pseudo measurements of nodes where PMUs are not observable are constructed, and compared with method 1 (using the estimation result of the previous moment as the pseudo measurement of the current moment) and method 2 (pseudo measurement construction method based on power transfer distribution factor). Figure 4 shown.

[0226] from Figure 4 It can be seen that method one has large fluctuations near the true value; method two considers the measurement of adjacent nodes and has a better estimation effect at first, but as the number of estimates increases, there will be a gradual deviation; the method of the present invention combines the spatiotemporal correlation between nodes, which to a certain extent suppresses the deviation of the estimation results of unobservable nodes and can be closer to the true value.

[0227] Table 2 shows the state estimation results after introducing pseudo measurements.

[0228] Table 2 Comparison of MAE indicators of state estimation under three pseudo-measurement construction methods

[0229]

[0230] It can be seen from Table 2 that the state estimation using the pseudo-measurement construction method proposed in the present invention has higher estimation accuracy for both observable nodes and unobservable nodes than Method 1 and Method 2, indicating that the proposed method not only has higher accuracy for the state estimation of unobservable nodes, but also has minimal impact on the state estimation of observable nodes.

[0231] (2) State estimation simulation verification

[0232] Based on the mixed measurement of step S1, steps S2 and S3 are used to adaptively update the noise covariance matrix, step S4 is used for state estimation, and the generalized maximum likelihood iterative extended Kalman filter (GMIEKF) algorithm, adaptive H∞Extended Kalman filter (AHEKF) algorithm and EKF algorithm are used for comparative verification.

[0233] First, we test under model uncertainty, assuming that the model noise parameters are unknown, and set the initial variances of process noise Q and measurement noise R to 10 -5 , 10 -5 ,Table 3 gives the state estimation results under model uncertainty.

[0234] Table 3 Comparison of evaluation indicators of four state estimation algorithms under unknown model parameters

[0235]

[0236] It can be seen from Table 2 that under model uncertainty, the IAREKF proposed in the present invention has the smallest estimation error, which shows that the method of the present invention can estimate the model noise quickly and accurately.

[0237] Next, we test the system in transient state. We add a short circuit fault to node 4 of the IEEE 39-node system and remove the fault after 0.1s. During this transient period, we perform state estimation. Figure 5 The state estimation results of the voltage amplitude and phase angle of node 3 are given. Table 4 gives the evaluation indicators of the four state estimation algorithms.

[0238] Table 6 Comparison of evaluation indicators of four state estimation algorithms during transient period

[0239]

[0240] Combination Figure 5 As can be seen from Table 4, the method of the present invention can quickly forget previous statistical information through adaptive adjustment of the forgetting factor during the transient period, accurately adjust the process noise covariance to adapt to the new operating conditions, and still maintain a high estimation effect during the transient period.

[0241] Finally, we tested with bad data, randomly adding bad data to the voltage measurement value (bad data is set to have an amplitude greater than 3% and a phase angle greater than 3°). Figure 6 The state estimation results of the voltage amplitude and phase angle of node 3 are given. Table 5 gives the evaluation indicators of four state estimation algorithms under bad data.

[0242] Table 5 Comparison of evaluation indicators of four state estimation algorithms under bad data

[0243]

[0244] Combination Figure 6 It can be seen from Table 5 that under the influence of bad data, the method of the present invention can still maintain a lower estimation error, showing good robustness of the algorithm.

[0245] It can be concluded from the simulation cases of the present invention that the proposed pseudo-measurement construction method has high accuracy, can effectively estimate the state of unobservable nodes, and will not cause adverse effects on observable nodes. The proposed adaptive robust extended Kalman filter method can accurately obtain the system operation state, has good adaptive adjustment ability under system disturbance, and has good robustness under bad data.

[0246] The invention is based on the hybrid measurement of the power system adaptive robust dynamic state estimation method. The proposed pseudo measurement construction method and the adaptive robust state estimation method can effectively perform global state estimation when the PMU configuration does not meet the system observability, and can accurately obtain the system operation state under the influence of unknown model parameters, system state changes and bad data. It has high application value.

[0247] This patent is not limited to the above-mentioned optimal implementation mode. Anyone can derive other forms of adaptive robust dynamic state estimation methods for power systems based on hybrid measurements under the inspiration of this patent. All equal changes and modifications made according to the scope of the patent application of the present invention should be covered by this patent.

Claims

1. A hybrid measurement-based adaptive robust dynamic state estimation method for power systems, characterized in that: The following steps are involved: Step S1, using measurement transformation technology to linearize the measurement equation, and constructing pseudo-measurement based on the spatiotemporal correlation between node states to achieve fusion of hybrid measurements; Step S2: Determine the filtering effect according to the statistical characteristics of the new information and assign weights to the new data, and adjust the process noise covariance in real time to adapt to the dynamic changes of the system; Step S3, estimating the measurement noise by analyzing the statistical characteristics of the measurement device, constructing an equivalent weight function according to the standard innovation, and correcting the measurement noise covariance to improve the algorithm robustness; Step S4: Based on the mixed measurement set of step S1, an extended Carr filter algorithm is used to adjust the noise covariance according to steps S2 and S3, perform dynamic state estimation, and obtain the system operation state.

2. The method for adaptive robust dynamic state estimation of a power system based on hybrid measurement according to claim 1, characterized in that: The fusion of hybrid measurements in step S1 is specifically achieved by: when SCADA measurements are updated, hybrid measurements are used for state estimation; when non-SCADA measurements are updated, PMU measurements combined with pseudo measurements are used for state estimation; The hybrid measurement is to integrate PMU and SCADA to form a hybrid measurement set, specifically: The node voltage phasor and branch current phasor are directly measured by PMU, and the node voltage phasor is used as the state quantity to obtain the linear measurement equation based on PMU: Where U is the node voltage state quantity, U p and I p are the vectors constructed from the PMU node voltage phasor and the branch current phasor respectively; the superscript p represents the PMU measurement; E p is the identity matrix, Y p is the node-branch admittance matrix, ν U and ν I are the measurement errors of voltage and current respectively; The SCADA data is transformed into branch power measurement and node injection power measurement into equivalent branch current phasor measurement and equivalent node injection current phasor measurement. The conversion formula is as follows: In the formula, and are the equivalent branch current phasor and the equivalent node injection current phasor respectively; P ft , Q ft , P f , Q f The active power and reactive power of branch ft and the injected active power and injected reactive power of node f respectively; U re,f and U im,f are the real and imaginary parts of the voltage at node f, respectively; After measurement transformation, the linear measurement equation based on SCADA is obtained: In the formula, and They are the vectors constructed from the equivalent branch current phasor measurement and the equivalent node injection current phasor measurement obtained by measurement transformation of SCADA data; the superscript s represents SCADA measurement; and are the node-branch admittance matrix and the node admittance matrix, respectively. and They are the measurement errors of node-branch equivalent current and node equivalent current respectively; By integrating the measurement equations of PMU and SCADA, the measurement equation based on hybrid measurement is finally obtained: Where x is the node voltage vector; z p and z s are the measurement vectors of PMU and SCADA respectively; H p and H s Respectively represent the coefficient matrices between PMU and SCADA and the node voltage vector; ν p and ν s represent the measurement noise vectors of PMU and SCADA respectively; The pseudo measurement is derived from the PMU measurement of the PMU unobservable node at the non-SCADA update time according to the data space correlation and time correlation, specifically: The pseudo-measurement constructed based on the power transfer distribution factor and the pseudo-measurement constructed based on the time correlation are subtracted from the state estimate, and the weight is set using an exponential function. The weighted sum is used to form the pseudo-measurement based on time-space correlation: In the formula, It is a pseudo-measurement constructed based on spatiotemporal correlation; and They are respectively a pseudo measurement constructed based on the power transfer distribution factor and a pseudo measurement constructed based on the time correlation at time k; and are the corresponding weights respectively; Represents the estimated value of the state at time k-1.

3. The method for adaptive robust dynamic state estimation of a power system based on hybrid measurement according to claim 2, characterized in that: The pseudo measurement constructed based on the power transfer distribution factor is to calculate the indirect measurement of the PMU directly observed node through the PMU measurement of the PMU directly observed node, and construct the power transfer distribution factor based on the indirect measurement to generate the power pseudo measurement of the unobservable branch, as follows: Assume that node i is equipped with a PMU and node j is not equipped with a PMU. Then the voltage and current of node j are indirectly calculated through the measurement of node i: In the formula, Represent the voltages of nodes i and j respectively, They represent the current flowing from node i to node j and the current flowing from node j to node i, respectively. ii , Y ij , Y ji , Y jj Represents the elements corresponding to the node-branch admittance matrix: Where U, θ are the amplitude and phase angle of the voltage, P, Q are the injected active and reactive powers of the node or branch, G, B are the conductance and susceptance of the node or branch, and the subscripts i, j, ij represent the corresponding node i, node j, and branch ij; The change of node voltage is derived when the node injection power changes: In the formula, ΔU j , Δθ j , ΔP j , ΔQ j They represent the voltage amplitude change, voltage phase angle change, injected active power change and injected reactive power change of node j respectively, and E j 、F j represents the intermediate variable corresponding to node j, E j =(G j +B j G j -1 B j ) -1 , F j =-(G j +B j G j -1 B j ) -1 B j G j -1 ; The power change of the PMU unobservable node injected is derived from the power change of the branch indirectly measured by the PMU: In the formula, ΔP ji , ΔQ ji They represent the active power change and reactive power change of the indirect measurement branch ji respectively. The variable M represents the relationship between the branch power measurement change and the node injection power change. It is composed of branch admittance. The subscript represents the corresponding node and branch, that is, the variable M j-ji It represents the relationship between the power measurement change of branch ji and the power injection change of node j; By calculating the injected power variation of the PMU unobservable node, the power variation of the PMU unobservable branch is derived: In the formula, ΔP jl , ΔQ jl They represent the active power change and reactive power change of the unobservable branch jl, respectively. The variable M j-jl It represents the relationship between the power measurement change of branch jl and the power injection change of node j; Finally, taking the power measurement of branch jl measured by SCADA at time k as the starting point, the power change of branch ji from time k to time k+1 is calculated based on the PMU measurement. Combining formula (9), the power change of branch jl that cannot be observed by PMU is obtained, and the power pseudo-measurement of branch jl at time k+1 is generated: In the formula, They represent the active power spatial pseudo-measurement and reactive power spatial pseudo-measurement of branch jl at time k+1 respectively; if time k is the SCADA update time, then is the active power measurement and reactive power measurement of branch jl measured by SCADA at time k; if time k is not a SCADA update time, then It refers to the spatial pseudo-measurement of active power and reactive power of branch jl at time k.

4. The method for adaptive robust dynamic state estimation of a power system based on hybrid measurement according to claim 2, characterized in that: The pseudo measurement constructed based on time correlation is designed based on the change relationship of each injection power of the same node within a period of time, and is specifically as follows: In the formula, represents the pseudo-measurement of branch jl power at time k+1; if time k is the SCADA update time, then is the power measurement of branch jl measured by SCADA at time k; if time k is not a SCADA update time, then It refers to the power time pseudo measurement of branch jl at time k. represents the indirect power measurement of branch ji by PMU at time k+1 and time k; b ji is the trend coefficient, is the power transformation ratio of branch ji and branch jl at the last two SCADA update times; S ji is the correlation coefficient, which is the correlation coefficient between branch ji and branch jl within a time window L; t s is the SCADA measurement cycle; T s is the time of the most recent SCADA update; n is the number of PMU direct observation nodes connected to node j, Indicates SCADA in T s Moment and T s -The power of branch ji at time ts, Indicates SCADA in T s Moment and T s -The power of branch jl at time ts; and They represent the power data of branch ji and branch jl in the time window L of SCADA respectively.

5. The method for adaptive robust dynamic state estimation of a power system based on hybrid measurement according to claim 1, characterized in that: The step S2 is based on adaptively adjusting the forgetting factor and introducing the process noise adaptive estimation of the gradually fading memory exponential weighting method: The process noise is the difference between the true state value and the predicted state value. The estimated value of the state is used to approximate the process noise: Where, the subscript k is used to represent the kth moment, ε k is the estimation error, ω k is the system process noise, x k is the system state value, f(x k ) is the state transfer equation, which reflects the dynamic changes of the distribution network operation state. is the estimated value of the state, is the Jacobian matrix of the state transfer equation, K k+1 is the Kalman gain, v k+1 For new breath; When the number of sampling moments in the sample space is N, the unbiased estimate of the process noise covariance is: In the formula, Q k is the process noise ω k The covariance matrix of is the mean of N ω samples, P is the estimation error covariance, In these N samplings, the weight factor d is introduced to give a larger weight to the latest sampled data. The weight factor d at time k is k for: Where b k is the forgetting factor. The greater the system process noise changes, the greater the k The smaller the value of ; Modify the process noise covariance to: The adaptive adjustment of the forgetting factor is specifically as follows: The intermediate variable η is designed by the ratio of the innovation variance and the innovation covariance matrix trace. k Measure the overall filtering effect of the current state estimate: Adaptively adjust the forgetting factor according to the filtering effect: In the formula, R k is the covariance matrix of the measurement noise, H k is the Jacobian matrix of the measurement equation, is the prediction error covariance matrix, b1 and b2 are the lower limit of the forgetting factor and the adaptive dynamic range respectively, and c1 and c2 are parameters to measure the filtering effect.

6. The method for adaptive robust dynamic state estimation of a power system based on hybrid measurement according to claim 5, characterized in that: The lower limit of the forgetting factor b1 is 0.6, b1+b2 is 0.95~0.995; c1=2, c2=3; when η k ≤2, that is, when the sum of the new standard deviation is less than 2 times the theoretical standard deviation, the filtering effect is good and the forgetting factor b is reduced. k To improve the degree of forgetting, the process noise can quickly track the accurate value; When 2<η k ≤3, the filtering effect is judged to be average, and the forgetting factor b is maintained k Stable, making the process noise stable estimation; when η k >3, indicating that the filtering process is disturbed, the forgetting factor b k If it is 1, the process noise is not updated.

7. The method for adaptive robust dynamic state estimation of a power system based on hybrid measurement according to claim 1, characterized in that: The method of estimating the measurement noise by analyzing the statistical characteristics of the measurement device is specifically as follows: before each state estimation, the measurement noise is adaptively estimated according to the latest measurement: In the formula, R(i,i) represents the i-th element of the diagonal of R, z p 、z s Respectively represent the corresponding PMU measurement and SCADA measurement; σ m and σ a They represent the amplitude error e measured by PMU respectively. m and phase error e a The standard deviation of , t is the intermediate variable, and σ P and σ Q are the standard deviations of SCADA active and reactive power measurement errors, U re and U im are the real and imaginary parts of the node voltage.

8. The method for adaptive robust dynamic state estimation of a power system based on hybrid measurement according to claim 1, characterized in that: The equivalent weight function is constructed according to the standard innovation, and the measurement noise covariance is modified as follows: Standard new information is used to identify and process bad data. The expression of the standard new information is: In the formula, v k For the new news, Indicates standard news. Based on the standard innovation, a robustness factor is introduced to correct the measurement noise. The robustness factor is generated using the IGGⅢ weight function: In the formula, Represents the robustness factor The i-th element of express The i-th element of , k0 and k1 are the upper and lower thresholds, which are parameters to measure the quality of the measurement value; Get the corrected measurement noise covariance matrix for: In the formula, R k represents the measurement noise covariance matrix at the kth moment, Represents the robustness factor.

9. The method for adaptive robust dynamic state estimation of a power system based on hybrid measurement according to claim 8, characterized in that: The parameters for measuring the quality of the measurement values ​​are k0=3, k1=4. When the measured data is greater than 3 times the standard deviation, it is judged as bad data and subjected to anti-error processing to reduce its weight in state estimation; when the measured data is greater than 4 times the standard deviation, the bad data is eliminated.

10. The hybrid measurement-based adaptive robust dynamic state estimation method for power systems according to claim 8, characterized in that: The Kalman filter in step S4 is divided into two processes: state prediction and state filtering. The specific process is shown as follows: State prediction process: State filtering process: Where, the subscript k is used to represent the kth moment, is the state prediction value; is the estimated value of the state; It represents the state transfer equation, reflecting the dynamic changes of the distribution network operation state. Represents the system measurement equation, reflecting the relationship between state quantity and quantity measurement, is the prediction error covariance matrix; is the estimated error covariance matrix; I represents the identity matrix, K k is the Kalman gain; Q k and are the process noise ω adjusted and corrected in step S2 and step S3 respectively k and measurement noise ν k The covariance matrix of is the Jacobian matrix of the measurement equation, which is obtained according to formula (4) under mixed measurement; is the Jacobian matrix of the state transfer equation, which is obtained according to the Holt exponential smoothing method: In the formula, is the state prediction value, is the state estimation value; a and b are the horizontal component and vertical component of Holt exponential smoothing respectively; α and β are the smoothing coefficients.

Images