A robust MM estimation-based unscented Kalman filter state estimation method and attack method
By using an unscented Kalman filter method based on robust MM estimation, the problems of observational anomalies and network attacks in generator state estimation are solved, achieving efficient and robust state estimation and ensuring power system stability.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- GUIZHOU UNIV
- Filing Date
- 2023-10-12
- Publication Date
- 2026-05-05
AI Technical Summary
Existing technologies cannot effectively handle uncertainties caused by observational anomalies and cyberattacks in generator state estimation, leading to power system stability problems.
An unscented Kalman filter method based on robust MM estimation is adopted. By statistically linearizing the nonlinear measurement function of the power system, a batch regression equation is constructed by combining the unscented Kalman filter algorithm. Robust pre-whitening and data point dispersion processing are performed, and state estimation is carried out by combining iterative weighted least squares method.
When faced with observational anomalies and cyberattacks, it can accurately estimate generator status, improve estimation efficiency and robustness, ensure power system stability, and achieve 95% estimation efficiency and high crash point performance.
Smart Images
Figure CN117272208B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to an unscented Kalman filter state estimation method and attack method based on robust MM estimation, belonging to the field of generator motion state estimation technology. Background Technology
[0002] In power systems, synchronous generators are core components, serving as power generation equipment. The real-time changes in the power grid and the operating state of the power system are significantly reflected in the dynamic variables of synchronous generators. Therefore, monitoring the state variables of synchronous generators to capture the dynamic changes of the power system is crucial. However, accurate monitoring and control of synchronous generators faces numerous challenges, such as generator model defects, system anomalies, and network attacks. The uncertainties brought about by multiple disturbances can lead to power system collapse, fires, and other safety accidents. Therefore, to ensure the stable operation of generators, it is urgent to address the observation and control problems caused by these uncertainties.
[0003] Dynamic state estimation of generators is a common method to eliminate uncertainties in power systems, and its study is of great significance for maintaining the stable operation of the power grid.
[0004] In dynamic state estimation of generators, the traditional Kalman filter algorithm achieves rapid estimation of system state through prior prediction, real-time observation, and dynamic model, but it has limitations in estimating the state of nonlinear systems. Some studies have proposed using the Extended Kalman Filter (EKF) for state estimation. This method uses Taylor series expansion to linearize the nonlinear system model, but the Jacobian determinant has high computational cost and diverges when the model is strongly nonlinear. Some studies have proposed using the Unscented Kalman Filter (UKF) for state estimation. This method approximates the propagation and observation process of the nonlinear function by using sampling points of the unscented transform, thus providing a more accurate state estimate. However, the above methods all have the following drawbacks: (1) they cannot handle observation anomalies; (2) they cannot handle network attacks where measurement data and model parameters are tampered with. Summary of the Invention
[0005] The technical problem to be solved by the present invention is to provide a method for spoofing attacks in the electricity market and a method for assessing the resilience against spoofing attacks, so as to overcome the shortcomings of the prior art.
[0006] The technical solution of this invention is as follows:
[0007] The first aspect is to provide a dynamic state estimation method for unscented Kalman filters based on robust MM estimation, including the following steps:
[0008] S1. Establish a power system state estimation model and discretize the model, initialize the state vector, covariance matrix, process noise covariance matrix and measurement noise covariance matrix, and transmit measurement data in real time through phasor measurement unit;
[0009] S2. Linearize the system measurement function using statistical linearization methods and establish a batch regression equation by combining the data calculated using unscented Kalman filtering;
[0010] S3. Perform robust pre-whitening on the batch regression equation established in step S2;
[0011] S4. Disperse the data points in the batch regression equation obtained in step S3 using a linear transformation method;
[0012] S5. Use the robust MM estimation method on the batch regression equation in step S4 to estimate the state vector and update the covariance matrix.
[0013] Furthermore, the power system state estimation model in step S1 is as follows:
[0014] Differential equations:
[0015]
[0016]
[0017]
[0018]
[0019] Algebraic equations:
[0020]
[0021] Where f0 represents the system frequency, D is the damping coefficient, J is the system inertial constant, and T′ qo Let T′ be the q-axis transient open-circuit time constant. do Let x be the d-axis transient open-circuit time constant. d For d-axis reactance, x q Let x′ be the q-axis reactance. d Let x′ be the d-axis transient reactance. q E is the q-axis transient reactance. fd Excitation voltage, V t For terminal voltage, T m Where δ is the mechanical input torque, ω is the rotor angle, and P is the rotor angular velocity. t e′ represents the active power at the terminal. q e′ represents the transient voltage on the q-axis. d This represents the transient voltage on the d-axis.
[0022] Furthermore, the method for discretizing the model is to use the second-order Runge-Taku method to discretize the model, and the discretized state-space equation is:
[0023]
[0024] Where, x k =[δ,ω,e′ q ,e′ d ] T It is the state vector at time k, z k =[P t ] is the measurement vector, u k =[T m E fd V t ] T It is the control vector, the state vector x k Includes rotor angle δ, rotor angular velocity ω, and transient voltage e′ on the d-axis. d and the transient voltage e′ on the q-axis q P t It is the terminal active power, u k Includes excitation voltage E fd Terminal voltage V t Mechanical input torque T m w k v represents the process noise at time k in the state-space equation. k For the measurement noise at time k in the state-space equation, w k v k The mean is 0 and the data are independent. Further, the method of linearizing the system measurement function using statistical linearization and establishing a batch regression equation based on the data calculated by unscented Kalman filtering includes the following steps:
[0025] S2-1, using the state estimate at time k-1 Covariance Matrix Generate 2n sigma points.
[0026]
[0027] Where n represents the state dimension, Representation matrix In the i-th column, λ = α 2 (n+κ)-n is a composite scaling factor, where α and κ are free parameters controlling the distribution of sigma points. The value of α is related to the dimension of the state vector; the higher the dimension, the smaller the value of α. The range of α is 10. -4 <α<1, κ is 2;
[0028] S2-2, Set weight w i =1 / (2n), and for each sigma point, a set of transformed samples is generated by mapping through the nonlinear function f(·). Through weight w i and Calculate the state prior value at time k. and prior covariance matrix
[0029]
[0030]
[0031] Among them, Q k Let the process noise covariance matrix at time k be denoted as .
[0032]
[0033] S2-3, Based on the prior state value at time k. and prior covariance matrix Generate new sigma points.
[0034]
[0035] S2-4. Map the new sigma using a nonlinear function h(·), generating a set of transformed samples, denoted as . Calculate the prior measurement vector Error covariance matrix and cross covariance matrix
[0036]
[0037]
[0038]
[0039] Among them, R k Let the measurement noise covariance matrix at time k be denoted as .
[0040]
[0041] S2-5. Calculate the Kalman gain.
[0042]
[0043] S2-6, Covariance Matrix
[0044]
[0045] S2-7. Using the sigma point obtained in step S2-3 and the results obtained in step S2-4 For these 2n points Using the least squares method to fit the data (i = 1, ..., n), a linear approximation function of the measurement function at time k is obtained:
[0046]
[0047] in, Let k represent the prior measurement vector at time k.
[0048] S2-8. Based on the prior prediction error at time k of the state vector and Establish batch regression equations.
[0049]
[0050] Where I represents an n-dimensional identity matrix.
[0051] Furthermore, the robust pre-whitening treatment method is as follows:
[0052] Calculate the covariance matrix W of the regression equation. k And the covariance matrix W of the regression equation k Perform Cholesky decomposition:
[0053]
[0054] The batch regression equation is obtained by multiplying both sides by the Cholesky decomposition. get:
[0055] y k =a k x k +b k
[0056] in, Error covariance matrix
[0057] Furthermore, the method for dispersing the data points in the batch regression equation obtained in step S3 through linear transformation in step S4 is as follows:
[0058] Construct the robust pre-whitening regression equation at time k:
[0059]
[0060] Multiply both sides of the regression equation by q i The regression equation obtained after the data points at time k are dispersed is as follows:
[0061]
[0062] Where, q i It is a random number on a uniform distribution U(1,10).
[0063] Furthermore, the method for estimating the state vector is as follows:
[0064] The state estimation equation is iterated using the iterative weighted least squares method. During the (j+1)th iteration... The iteration stops when the state estimation equation is:
[0065]
[0066] Among them, M k For m k,i The matrix form of W j The weighting function ω(u) for the j-th iteration i The matrix form of ω(u) i )=ψ(u i ) / u i C k c at time k k,i =q i *y k,i Matrix form;
[0067] And c = 4.685;
[0068] u i =e i / σ, where σ is the robust scaling estimate obtained by solving the S-estimation. ρ(u i ) is the weight function.
[0069]
[0070] Where c k,i =m k,i x k +d k,i i = 1, ..., 30 The simplified form,
[0071] c k,i =q i *y k,i m k,i =q i *ak,i d k,i =q i *b k,i .
[0072] Secondly: An attack method is provided, the method being:
[0073] From 30c k,i and m k,i We select N1 and N2 elements respectively, where N1 and N2 are both less than 15, and inject the following attack at each time k:
[0074]
[0075]
[0076] Among them, t j Let U(1.5,2.5) represent a random number on a uniform distribution U(1.5,2.5). randn(5,1) represents generating a random number vector of size 5 rows and 1 column, where each element is an independent and identically distributed random number from a standard normal distribution.
[0077] The beneficial effects of this invention are: compared with the prior art,
[0078] 1) This invention linearizes the nonlinear measurement function of the power system using a statistical linearization method, and constructs a batch regression equation based on the data calculated by unscented Kalman filtering. The regression equation is pre-whitened and the data points are dispersed for use in MM estimation. After estimation, a more accurate state vector can be obtained. Thus, the advantages of MM estimation and unscented Kalman filtering are combined. Under the premise of ensuring estimation performance and being able to estimate the state of nonlinear systems, MM estimation can obtain high collapse points, enabling this method to solve the network attack problems of observation anomalies, measurement data and model parameters being tampered with, and also achieve an estimation efficiency of 95%.
[0079] 2) This invention proposes to use a linear transformation method to disperse the data points in the batch regression equation, so that MM estimation can accurately estimate the state vector;
[0080] 3) This invention inherits the high collapse point of MM estimation, and can still estimate accurate results even with 40% erroneous data. Attached Figure Description
[0081] Figure 1 This is a flowchart of the present invention;
[0082] Figure 2 Use MM estimation to fit the dense data points;
[0083] Figure 3 The MM estimation was used to fit the scattered data points;
[0084] Figure 4 This invention provides a comparative estimation method with EKF, UKF, and M-UKF under conditions of no observed anomalies and no network attacks.
[0085] Figure 5 This invention, along with EKF, UKF, and M-UKF, measures rotor angle under conditions of observed anomalies but without network attacks.
[0086] Estimation and comparison;
[0087] Figure 6 This invention, along with EKF, UKF, and M-UKF, measures rotor angle under conditions of observed anomalies but without network attacks.
[0088] Speed estimation comparison;
[0089] Figure 7 This invention, along with EKF, UKF, and M-UKF, temporarily modulates the q-axis in the event of anomalies but without network attacks.
[0090] Comparison of state voltage estimation;
[0091] Figure 8 This invention, along with EKF, UKF, and M-UKF, temporarily modulates the d-axis in the event of anomalies but without network attacks.
[0092] Comparison of state voltage estimation;
[0093] Figure 9 Comparison of the estimation of this invention with M-UKF when N1=5 and N2=4 in network attacks;
[0094] Figure 10 Comparison of the estimation of this invention with M-UKF when N1=0, N2=12 in network attacks;
[0095] Figure 11 The second-order Runge-Takrug method used in this invention and two other discretization methods are compared in terms of the time variation of rotor angular velocity at a time interval of 0.001s.
[0096] Figure 12 The second-order Runge's method used in this invention and two other discretization methods are compared in terms of the time variation of rotor angular velocity at a time interval of 0.03s.
[0097] Figure 13 The second-order Runge-Takrug method used in this invention and two other discretization methods are compared in the time variation diagram of rotor angular velocity at a time interval of 0.05s. Detailed Implementation
[0098] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings.
[0099] Example 1: As shown in the attached document Figure 1 As shown, an unscented Kalman filter algorithm based on MM estimation for dynamic state estimation of power systems is provided, including the following steps:
[0100] S1. Establish a power system state estimation model, and initialize the state vector, covariance matrix, process noise, and measurement noise covariance matrix. Transmit measurement data in real time through the phasor measurement unit.
[0101] Specifically, the establishment of the power system model in step S1 is as follows:
[0102] Differential equations:
[0103]
[0104]
[0105]
[0106]
[0107] Algebraic equations:
[0108]
[0109] Where f0 represents the system frequency, D is the damping coefficient, J is the system inertial constant, and T′ qo Let T′ be the q-axis transient open-circuit time constant. do Let x be the d-axis transient open-circuit time constant. d For d-axis reactance, x q Let x′ be the q-axis reactance. d Let x′ be the d-axis transient reactance. q E is the q-axis transient reactance. fd Excitation voltage, V t For terminal voltage, T m Where δ is the mechanical input torque, ω is the rotor angle, and P is the rotor angular velocity. t e′ represents the active power at the terminal. q e′ represents the transient voltage on the q-axis. d This represents the transient voltage on the d-axis.
[0110] After weighing the accuracy and computational complexity, the state-space equations are discretized using the second-order Runge-Taku method:
[0111]
[0112] Where, x k =[δ,ω,e′ q ,e′ d ]T It is the state vector and the measurement vector z. k =[P t and control vector u k =[T m E fd V t ] T State vector x k Includes rotor angle δ, rotor angular velocity ω, and transient voltage e′ on the d-axis. d and the transient voltage e′ on the q-axis q , z k It is a measurement vector (terminal active power P) t ), control input u k Includes excitation voltage E fd Terminal voltage V t Mechanical input torque T m w k v represents the process noise at time k in the state-space equation. k For the measurement noise at time k in the state-space equation, w k v k The mean is 0 and they are independent.
[0113] The second-order Runge-Taku method considers differential equations:
[0114]
[0115] The specific steps for discretizing differential equations using the second-order Runge-Taku method are as follows:
[0116]
[0117] k1=f(x k-1 ,t k-1 )Δt
[0118] k2=f(x k-1 +k1,t k +Δt)Δt
[0119] This method is equivalent to considering the first and second derivative terms in the Taylor series, and therefore has higher accuracy than the Euler method.
[0120] To compare the Euler method, the second-order Runge-Taku method, and the fourth-order Runge-Taku method, this patent uses three methods to discretize the system model and plots the time variation of the rotor angular velocity at different time intervals, as shown below. Figure 11-13 As shown, Figure 11 The graph shows the time variation of rotor angular velocity for the three methods at a time interval of 0.001s. At a time interval of 0.001s, the accuracy of the three methods is almost the same, which is logical because the smaller the time interval, the smaller the error. Figure 12 The graph shows the time variation of rotor angular velocity for the three methods at a time interval of 0.03s. Compared with the fourth-order Runge-Taco method, the Euler method has the largest error, which is about 0.002; the second-order Runge-Taco method has a very small error, which is less than 0.001. Figure 13 The graphs show the time variation of rotor angular velocity for three methods at a time interval of 0.05s. The errors of the Euler method and the second-order Runge-Taku method are relatively large.
[0121] Under equal conditions, the fourth-order Runge-Taku method has the highest accuracy compared to the Euler method and the second-order Runge-Taku method. However, the fourth-order Runge-Taku method also has the highest computational complexity. Therefore, this patent compares the accuracy of the Euler method and the second-order Runge-Taku method relative to the fourth-order Runge-Taku method under different conditions. Considering both accuracy and computational complexity, the second-order Runge-Taku method with a time interval of 0.03s is chosen to discretize the system model.
[0122] S2. Linearize the system measurement function using statistical linearization methods and establish a batch regression equation by combining the data calculated using unscented Kalman filtering;
[0123] Specifically, in step S2, an unscented Kalman filter algorithm is used to calculate... and
[0124] The prediction phase first uses an unscented transformation, employing the state estimate at time k-1. Covariance Matrix Generate 2n sigma points:
[0125]
[0126] Where n represents the state dimension, Representation matrix In the i-th column, λ = α 2 (n+κ)-n is a composite scaling factor, where α and κ are free parameters controlling the distribution of sigma points. The value of α is related to the dimension of the state vector; the higher the dimension, the smaller the value of α. The typical range of α is 10. -4 <α<1, κ is usually set to 2.
[0127] Set weight w i = 1 / (2n), each sigma point is mapped through a nonlinear function f() to generate a set of transformed samples, denoted as Then the prior state value x at time k k|k-1 and prior covariance matrix It can be calculated using equations (9) and (10):
[0128]
[0129]
[0130] Among them, Q k This represents the process noise covariance matrix at time k.
[0131] During the update phase, the prior state values obtained from equations (9) and (10) are first used. and prior covariance matrix Generate updated sigma points:
[0132]
[0133] Then, each sigma point is mapped through a nonlinear function h(·) to generate a set of transformed samples, denoted as . Next, the prior measurement vector is calculated. Error covariance matrix and cross covariance matrix
[0134]
[0135]
[0136]
[0137] Calculate Kalman gain
[0138]
[0139] Finally, for The final state estimate is obtained by making corrections. And update
[0140]
[0141]
[0142] Specifically, in step S2, a statistical linearization method is used to linearize the system measurement function and establish a batch regression form:
[0143]
[0144] in, Let k represent the prior measurement vector at time k. in, No longer the Jacobian determinant, the prediction error of the state vector is expressed as:
[0145]
[0146] Combining equations (18) and (19), the batch regression form is obtained as follows:
[0147]
[0148] Where I represents an n-dimensional identity matrix, and It is calculated using equations (9) and (12) in the unscented Kalman filter algorithm. The above equation can be simplified as:
[0149]
[0150] S3. Perform robust pre-whitening on the batch regression equation established in step S2;
[0151] Specifically, in step S3: the batch regression equation established in step S2 is subjected to robust pre-whitening processing;
[0152] The batch processing error covariance matrix can be calculated using the following formula.
[0153]
[0154] Among them, S k It is obtained from Cholesky decomposition.
[0155] Pre-whitening transforms the independent and dependent variables into irrelevant whitened variables, making their covariance matrices identity matrices. This simplifies the interpretation and inference process of the regression model, reducing redundant information and computational complexity. Multiplying both sides of equation (21) by... get:
[0156] y k =a k x k +b k (twenty three)
[0157] in, Error covariance matrix That is, b k There is no correlation between any two features.
[0158] S4. Disperse the data points in the batch regression equation obtained in step S3 using a linear transformation method;
[0159] Specifically, in step S4: the data points in the batch regression equation obtained in step S3 are dispersed by a linear transformation method;
[0160] Because the phasor measurement unit (PMU) has a high sampling rate, at least 30 measurement data can be collected at time k. For the batch regression equation y in equation (23)... k =a k x k +b k Each time data is measured, a y-value can be obtained. k a k b k Therefore, for any time k, we can obtain 30 y values. k a k b k Then the regression equation at time k can be expressed as:
[0161]
[0162] This meets the conditions for using MM estimation. However, because the measurement time interval is too short, the measurement data is too dense, and direct use of MM estimation may not yield accurate results.
[0163] Therefore, this invention proposes that before using MM estimation, these 30 data points need to be dispersed through linear transformation to ensure accurate MM estimation results. The method for dispersing the data points is as follows:
[0164] Multiply both sides of each equation by q i get:
[0165] q i *y k,i =q i *a k,i x k +q i *b k,i ,i=1,…,30 (25)
[0166] Where, q i It is a random number on a uniform distribution U(1,10). At this point, the regression equation at time k can be expressed as:
[0167]
[0168] The above formula can be simplified as follows:
[0169] c k,i =m k,i x k +d k,i ,i=1,…,30 (27)
[0170] Where c k,i =qi *y k,i m k,i =q i *a k,i d k,i =q i *b k,i .
[0171] To verify the impact of the data point set on the MM estimation accuracy of the state vector estimation, this application first constructs a function y = 8 - 2.1x, then generates several random numbers near 0 on the x-axis, substitutes these random numbers into the function to obtain the y values, and adds a perturbation term to each y value to obtain... Figure 2 The scatter points generated by the above methods are relatively concentrated, mostly near 0 on the x-axis. If we use MM estimation to obtain the fitted curve of the MM estimation, that is... Figure 2 As shown by the solid line in the graph, the fitted curve estimated by MM differs significantly from the true curve.
[0172] To verify the impact of data point dispersion on the accuracy of MM estimation on state vector estimation, this application first constructs a function y = 8 - 2.1x, then generates several random numbers in the x-axis interval from -3 to 3. Substituting these random numbers into the function yields the y values. Finally, a perturbation term is added to each y value to obtain the result. Figure 3 The scattered points generated by the above methods are quite dispersed. If we use MM estimation to obtain the fitted curve of the MM estimation, that is... Figure 3 As shown by the solid line in the graph, the fitted curve estimated by MM is relatively close to the true curve.
[0173] S5. Use the robust MM estimation method on the batch regression equation in step S4 to estimate the state vector and update the covariance matrix;
[0174] Specifically, in step S5: the robust MM estimation method is used on the batch regression equation in step S4 to estimate the state vector and update the covariance matrix;
[0175] MM combines the advantages of S-estimation and M-estimation, and it provides a robust estimate of the residual scale σ obtained from S-estimation. s Instead of the inrory scaling estimate MAD / 0.6745 (MAD represents the median absolute residual) in the M-estimation, a regression estimation method with high collapse point and high estimation efficiency is obtained. For the regression equation of equation (27), the robust scaling estimate σ is defined as follows:
[0176]
[0177]
[0178] Among them, u i =e i / σ, σ=MAD / 0.6745 is initial value, For the residual, ω(u) i Let represent the weighting function. To ensure a 50% crash point, K = 0.1996, c = 1.5476. Furthermore, the MM estimate is defined as follows:
[0179]
[0180] Where u i =e i / σ, σ is the robust scaling estimate obtained by S estimation, and ρ is the weight function. In order to enable the MM-UKF algorithm to remove the influence of multiple perturbations to the greatest extent, this paper adopts Tukey's double weight function as the weight function and takes c = 4.685.
[0181]
[0182] To estimate x k Find equation (30) with respect to x. k Taking the partial derivative of and setting it equal to 0, we get:
[0183]
[0184] in,
[0185]
[0186] Through a weighting function ω(u) i The solution to equation (32) is given:
[0187]
[0188] Where, ω(u i )=ψ(u i ) / u i From equation (29), in matrix representation, the above equation can be simplified to:
[0189]
[0190] Finally, using iterative weighted least squares, the result of the (j+1)th iteration is:
[0191]
[0192] when Stop iterating when the time comes. This is the estimated state vector at time k. The covariance matrix is then updated using equation (17).
[0193] Figure 4 The method of this application is compared with the estimation of EKF, UKF, and M-UKF under the condition of no observed anomalies and no network attacks. It can be seen that under the condition of no observed anomalies and no network attacks, the estimation of this application is almost the same as that of EKF, UKF, and M-UKF, and they are all very close to the true value.
[0194] and Figure 5 This invention is compared with EKF, UKF, and M-UKF in rotor angle estimation under the condition of observational anomalies but no network attack. Figure 6 This invention is compared with EKF, UKF, and M-UKF in terms of rotor angular velocity estimation under the condition of observational anomalies but no network attack. Figure 7 This invention is compared with EKF, UKF, and M-UKF for q-axis transient voltage estimation under conditions of observational anomalies but no network attacks. Figure 8 This paper compares the d-axis transient voltage estimation of this invention with EKF, UKF, and M-UKF under the condition of observational anomalies but no network attack. It can be seen that the methods of this application and M-UKF are closer to the true values, while M-UKF, UKF, and EKF differ significantly from the true values. This indicates that the method of this application and M-UKF can estimate the state more accurately under the condition of observational anomalies but no network attack.
[0195] Figure 9 A comparison of the estimations of this invention and M-UKF under network attacks N1=5 and N2=4 shows that the method of this patent is closer to the true value, while M-UKF differs more from the true value. This indicates that the method of this application can estimate the state more accurately than M-UKF when there are fewer outliers.
[0196] Figure 10 For the comparison of the estimation of this invention and M-UKF under network attacks N1=0, N2=12, outlier ratio Figure 9 Furthermore, it can be seen that the method of this patent is closer to the true value, while M-UKF differs greatly from the true value. This indicates that when there are many outliers, the method of this application has a lower collapse point. When faced with a large number of outliers, this method can estimate the state more accurately than M-UKF.
[0197] Example 2: This invention provides an attack scheme, and the performance of this invention is tested under this attack scheme, as detailed below:
[0198] As can be seen from equations (9), (12), and (20), tampering with measurement data or model parameters will directly pass to c in equation (36). k,iand m k,i Therefore, for these 30 c k,i and m k,i Select N1 and N2 respectively, where N1 and N2 are both less than 15, and inject the following attack at each time k.
[0199]
[0200] Among them, t j Let U(1.5,2.5) represent a random number on a uniform distribution U(1.5,2.5). randn(5,1) represents generating a random number vector of size 5 rows and 1 column, where each element is an independent and identically distributed random number from a standard normal distribution.
[0201] All aspects not detailed herein are well-known to those skilled in the art. Finally, it should be noted that the above embodiments are merely illustrative of the technical solutions of this invention and not intended to limit it. Although the invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of this invention without departing from the spirit and scope of the invention, and all such modifications and substitutions should be covered within the scope of the claims of this invention.
Claims
1. A dynamic state estimation method for unscented Kalman filters based on robust MM estimation, characterized in that, Includes the following steps: S1. Establish a power system state estimation model and discretize the model, initialize the state vector, covariance matrix, process noise covariance matrix and measurement noise covariance matrix, and transmit measurement data in real time through phasor measurement unit; S2. Linearize the system measurement function using statistical linearization methods and establish a batch regression equation by combining the data calculated using unscented Kalman filtering; S3. Perform robust pre-whitening on the batch regression equation established in step S2; S4. Disperse the data points in the batch regression equation obtained in step S3 using a linear transformation method; S5. Use the robust MM estimation method on the batch regression equation in step S4 to estimate the state vector and update the covariance matrix; The linear transformation method involves multiplying both sides of the regression equation by a random number, where the random number is uniformly distributed. A random number.
2. The unscented Kalman filter dynamic state estimation method based on robust MM estimation according to claim 1, characterized in that, The power system state estimation model in step S1 is as follows: Differential equations: ， ， ， ， Algebraic equations: ， in, Indicates the system frequency. The damping coefficient is... Let be the system's inertial constant. The q-axis transient open-circuit time constant is... The d-axis transient open-circuit time constant is... For d-axis reactance, For q-axis reactance, For d-axis transient reactance, For q-axis transient reactance, For excitation voltage, For terminal voltage, For mechanical input torque, For rotor angle, For rotor angular velocity and For terminal active power, express Transient voltage on the shaft, This represents the transient voltage on the d-axis.
3. The unscented Kalman filter dynamic state estimation method based on robust MM estimation according to claim 2, characterized in that, The method for discretizing the model is to use the second-order Runge-Taku method to discretize the model. The discretized state-space equation is as follows: ， in, It is the state vector at time k. It is a measurement vector. It is the control vector and the state vector. Including rotor angle Rotor angular velocity as well as transient voltage on the shaft and transient voltage on the shaft , It is the terminal active power. Includes excitation voltage Terminal voltage Mechanical input torque , This represents the process noise at time k in the state-space equations. The measurement noise at time k in the state-space equation is... The mean is 0 and they are independent.
4. The unscented Kalman filter dynamic state estimation method based on robust MM estimation according to claim 3, characterized in that, The method of linearizing the system measurement function using statistical linearization and establishing a batch regression equation based on the data calculated by unscented Kalman filtering includes the following steps: S2-1, by using State estimate at time 1 Covariance Matrix generate 1 Sigma point, ， in, Indicates the dimension of the state. Representation matrix The List, It is a composite proportionality factor. and These are the free parameters that control the distribution of sigma points. The value of is related to the dimension of the state vector; the higher the dimension, the better. The smaller the value, the smaller the range of values. , It is 2; S2-2, Setting Weights For each sigma point, a nonlinear function is used. Perform mapping to generate a set of transformed samples, denoted as . Through weight and Calculate the first State prior value at time 1 and prior covariance matrix , ， ， in, Let the process noise covariance matrix at time k be denoted as . ； S2-3, according to the first State prior value at time 1 and prior covariance matrix Generate new sigma points. ， S2-4, Applying a nonlinear function to the new sigma Perform mapping to generate a set of transformed samples, denoted as . Calculate the prior measurement vector Error covariance matrix and cross covariance matrix , ， ， ， in, Let the measurement noise covariance matrix at time k be denoted as . ， S2-5. Calculate the Kalman gain. ， S2-6, Update the covariance matrix ; ； S2-7. Using the sigma point obtained in step S2-3 and the results obtained in step S2-4 For these 2n points ( , A linear approximation function of the measurement function at time k is obtained by fitting using the least squares method: ； in, Let k represent the prior measurement vector at time k. ; S2-8. Prediction error at time k based on the state vector. and Establish batch regression equations. ， in, express 3D identity matrix.
5. The unscented Kalman filter dynamic state estimation method based on robust MM estimation according to claim 4, characterized in that, The robust pre-whitening treatment method is as follows: Calculate the covariance matrix of the regression equation And the covariance matrix of the regression equation Perform Cholesky decomposition: ， The batch regression equation is obtained by multiplying both sides by the Cholesky decomposition. get: ， in, , , Error covariance matrix .
6. The unscented Kalman filter dynamic state estimation method based on robust MM estimation according to claim 5, characterized in that, The method for dispersing the data points in the batch regression equation obtained in step S3 through linear transformation in step S4 is as follows: Construct the robust pre-whitening regression equation at time k: ， Multiply both sides of the regression equation by The regression equation obtained after the data points at time k are dispersed is as follows: ， in, It is a uniform distribution A random number.
7. The unscented Kalman filter dynamic state estimation method based on robust MM estimation according to claim 6, characterized in that, The method for estimating the state vector is as follows: The state estimation equation is iterated using the iterative weighted least squares method. During the (j+1)th iteration... The iteration stops when the state estimation equation is: ， in, for In matrix form, The weighting function for the j-th iteration In matrix form, ; For time k Matrix form; and ; , The robust scaling estimate is obtained by solving the S-estimation. ; For the weight function, ; ， in for The simplified form, ， ， 。 8. An attack method according to claim 7, characterized in that, The method is as follows: From 30 and Select from the options below. Individual and One, and and All less than 15, at each time point All of them are injected with the following attacks: ， in, Indicates uniform distribution A random number on the screen. This indicates the generation of a random number vector of size 5 rows and 1 column, where each element is an independent and identically distributed random number from a standard normal distribution.