A probabilistic encoding and power system distributed estimation method under spoofing attack

Abstract

The application discloses a kind of probability coding and power system distributed estimation method under spoofing attack, comprising the following steps: step 1: establish power system dynamic model;Step 2: complete the design of flexible distributed state estimator;Step 3: calculate the upper bound of one-step prediction error covariance matrix of the i th node at s time ∑ i,s+1|s ;Step 4: calculate the gain matrix K i,s+1 of the i th node of flexible distributed estimator;Step 5: calculate the state estimation matrix of the i th node at s+1 time step 6: according to the estimator gain matrix K i,s+1 obtained in step 4, calculate the upper bound of estimation error covariance matrix of the i th node at s+1 time ∑ i,s+1|s+1 The estimation method provided by the application solves the problem that the existing state estimation method is difficult to handle the probability coding and spoofing attack of power system.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of power network system control, in particular to a probability coding and power system distributed estimation method under deception attack. BACKGROUND

[0002] With the development of economic society, the demand for electricity by customers continues to grow and shows a diversification trend. Therefore, smart grid enhances system flexibility and efficiency by integrating renewable energy, energy storage technology and advanced electronic components with intelligent technology. However, the high integration results in high complexity of the system, and the system often operates under high load, making it extremely sensitive to external interference. To ensure system safety, real-time monitoring and control of synchronous motor state is crucial. Wide-area monitoring systems collect power data in real time through phasor measurement units, and then use dynamic state estimation technology to accurately estimate the generator state, improving system monitoring capability. However, due to the nonlinear characteristics of the power system, more advanced estimation algorithms need to be developed to address the challenges of complex systems and ensure the safe, stable and economic operation of the power system.

[0003] In sensor networks, wireless communication is often subject to inherent physical limitations, including limited bandwidth and energy scarcity, so transmitting more data within a certain time is crucial for stable system operation. As an efficient data compression technology, the encoding-decoding mechanism has always been an important tool for relieving communication pressure and promoting efficient data transmission. In many encoding schemes, probability coding stands out with its unique feature of ensuring unbiasedness throughout the encoding-decoding process. The overall framework is completed in a probabilistic way, which is more efficient and practical than traditional schemes. At the same time, data security in network communication cannot be ignored. Among numerous network attacks, deception attacks, which inject malicious data to disrupt the system, are a serious security threat. Therefore, while improving data transmission efficiency, we cannot ignore the prevention of such network attacks.

[0004] Currently, there is a lack of research on effectively integrating distributed schemes into probability coding mechanisms and dealing with random deception attacks in the field of state estimation of smart power systems, especially without fully considering the resilience design of estimators under potential disturbances. This means that when faced with complex and variable power system environments, traditional state estimation and filtering methods may not be up to the task, making it difficult to ensure the sustained stability of system performance, and innovative solutions are needed to fill this technical gap. SUMMARY

[0005] The technical problem to be solved by the present application is to provide a probability coding and power system distributed estimation method under deception attack, which can effectively complete the state estimation in the power system under complex scenarios.

[0006] To solve the above technical problems, the specific technical solutions of the present application are as follows:

[0007] Step 1: Establishment of the power system dynamic model

[0008] Based on the smart grid model, select the rotor angle, rotor angular velocity, d-axis transient voltage, and q-axis transient voltage of the power system as state variables to establish a power system dynamic model;

[0009] Step 2: Design of the resilient distributed state estimator

[0010] Under the probability coding mechanism and stochastic deception attacks, complete the design of the resilient distributed state estimator according to the power system dynamic model in Step 1;

[0011] Step 3: Calculate ∑ i,s+1|s

[0012] According to the resilient distributed state estimator in Step 2, calculate the one-step prediction matrix of the i-th sensor node at time s Thereby calculate the upper bound ∑ of the one-step prediction error covariance matrix of the i-th sensor node at time s i,s+1|s ;

[0013] Step 4: Calculate K i,s+1

[0014] According to the upper bound ∑ of the one-step prediction error covariance matrix obtained in Step 3 i,s+1|s , calculate the resilient distributed estimator gain matrix K of the i-th sensor node at time s + 1 i,s+1 ;

[0015] Step 5: Calculate

[0016] Substitute the estimator gain matrix K obtained in Step 4 i,s+1 into the resilient distributed state estimator in Step 2 to calculate the state estimation matrix of the i-th sensor node at time s + 1 Judge whether s + 1 reaches the total duration T. If s + 1 < T is satisfied, continue to execute Step 6; otherwise, end the operation;

[0017] Step 6: Calculate ∑ i,s+1|s+1

[0018] According to the estimator gain matrix K obtained in Step 4 i,s+1 , calculate the upper bound Σ of the estimation error covariance matrix of the i-th sensor node at time s + 1 i,s+1|s+1 , meanwhile, let s = s + 1, execute Step 2 until s + 1 reaches the total duration T.

[0019] Preferably, in Step 1, the calculation formula of the power system dynamic model is as follows:

[0020]

[0021] wherein, is a random nonlinear function based on the smart grid model,

[0022] wherein, x s = [δ s ω s E′ q,s E′ d,s ] T is the state vector of the power system at time s, wherein δ s is the rotor angle of the power system at time s, ω s is the rotor angular velocity of the power system at time s, E′ q,s is the q-axis transient voltage of the power system at time s, E′ d,s is the d-axis transient voltage of the power system at time s; μ s = [V s θ s T m,s E fd,s ] T is the known input vector of the power system at time s, wherein V s is the terminal bus voltage amplitude of the power system at time s, θ s is the terminal bus voltage phase angle of the power system at time s, T m,s is the mechanical torque input of the power system at time s, E fd,s is the damping coefficient of the power system at time s;

[0023] x s+1 is the state vector of the power system at time s+1;

[0024] ε s is a process noise with zero mean and variance R s >0 at time s;

[0025] z i,s is the actual measurement value of the i th sensor at time s;

[0026] g i,s = [f i,s P i,s Q i,s ] T is the ideal measurement value of the i th sensor at time s, wherein f i,s is the terminal frequency measured by the i th sensor node of the power system at time s, P i,s is the active power measured by the i th sensor node of the power system at time s, Q i,s is the reactive power measured by the i th sensor node of the power system at time s;

[0027] is the measurement noise of the ith sensor node at time s, which satisfies zero mean and covariance k >0.

[0028] The specific steps of the design of the recursive state estimator in step 2 are as follows:

[0029] Step 2.1: Based on the probability encoding scheme, calculate the probability encoding-decoding measurement value according to the measurement value of the ith sensor node in step 1, and the calculation formula of the probability encoding-decoding measurement value is as follows:

[0030]

[0031] In the formula, is the decoding measurement value of the ith sensor node at time s;

[0032] represents the decoding error of the ith sensor node at time s under the probability encoding mechanism, and the specific calculation formula is as follows:

[0033]

[0034] In the formula, is the encoding probability, is the encoding level, represents the encoding interval of the encoder, and k represents the numerical description corresponding to the encoding level.

[0035] Step 2.2: Considering the random occurrence of spoofing attacks in network communication, calculate the corresponding transmission measurement of the ith sensor node according to the decoding measurement value of the ith sensor node obtained in step 2.1, that is, the final received output measurement value of the estimator, and the calculation formula of the final received output measurement value of the estimator is as follows:

[0036]

[0037] In the formula, represents the malicious signal sent by the attacker to the ith sensor node at time s, where r i,s satisfies is a known scalar;

[0038] η i,s is the probability of spoofing attack of the ith sensor node, and the mathematical expectation is

[0039] Step 2.3. Designing the probabilistic encoding mechanism and the resilient distributed state estimator in the dynamic model of the power system under random deception attacks according to the output measurement values finally received by the estimator obtained in step 2.2, the calculation formula of the resilient distributed state estimator is as follows:

[0040]

[0041] wherein, represents the one-step prediction value of the state vector x s+1 of the i th sensor node at time s;

[0042] respectively represent the estimated value of the state vector x s of the i th sensor node at time s;

[0043] respectively represent the estimated value of the state vector x s+1 of the i th sensor node at time s+1;

[0044] K i,s+1 represents the gain of the distributed state estimator of the i th sensor node at time s+1;

[0045] ΔK i,s+1 =ξ i,s C i represents the gain mutation of the distributed state estimator of the i th sensor node at time s+1, wherein, in the formula, C i is a known matrix, and ξ i,s is a Gaussian white noise sequence with a mean of zero and a variance of ;

[0046] represents the predicted value of the measurement value of the i th sensor node at time s;

[0047] is the mathematical expectation of the probability η i,s of the deception attack of the i th sensor node;

[0048] I represents a known unit matrix with appropriate dimensions;

[0049] a ij represents the connection coefficient of the i th sensor node and the j th sensor node;

[0050] represents the adjacent node set of the i th sensor node.

[0051] Preferably, the calculation formula of the one-step prediction error covariance matrix upper bound Σ i,s+1|s at time s in step 3 is as follows:

[0052]

[0053] In the formula, Σ i,s|s It is the upper bound of the estimation error covariance matrix of the i-th sensor node at time s;

[0054] A i,s yes right Point to ask about x s The matrix obtained by the partial derivatives;

[0055] U s and L s It is a matrix of known appropriate dimensions;

[0056] λ1 is a known positive scalar. It is the reciprocal of λ1;

[0057] R s It is the covariance matrix of the process noise at time s;

[0058] It is Σ i,s|s The inverse matrix;

[0059] They are A i,s U s L s The transpose of .

[0060] Preferably, the gain matrix K in step 4 i,s+1 The specific calculation formula is as follows:

[0061]

[0062] In the formula,

[0063]

[0064] In the formula, ρ1 is the first intermediate parameter in the gain matrix;

[0065] ρ2 is the second intermediate parameter in the gain matrix;

[0066] ρ3 is the third intermediate parameter in the gain matrix;

[0067] λ² is a known positive scalar. It is the reciprocal of λ²;

[0068] B i,s+1 It is g i,s right Point to ask about x s The matrix obtained by the partial derivatives;

[0069] F s+1 and G s+1 are matrices of known appropriate dimensions;

[0070] O i is an intermediate matrix;

[0071] and are intermediate matrices generated in the process of handling the spoofing attack;

[0072] W i,s+1 is the covariance matrix of the process noise of the i-th sensor node at time s+1;

[0073] are the transpose matrices of B i,s+1 , F s+1 , G s+1 , O i , respectively;

[0074] is the inverse matrix of ; s+1|s

[0075] represents the Hadamard product operation.

[0076] Preferably, the estimation error covariance matrix upper bound in step 6 is i,s+1|s+1 The specific calculation formula is as follows:

[0077]

[0078] wherein, is the estimation error covariance upper bound of the i-th sensor node at time s+1; i,s+1|s+1

[0079] are the transpose matrices of K i,s+1 and. i,s+1

[0080] Beneficial effects

[0081] Compared with the prior art, the present application has the following advantages and positive effects:

[0082] The elastic distributed state estimation method comprehensively considers the probability encoding mechanism, random spoofing attack and estimator gain mutation and the like complex conditions, and compared with the traditional encoding scheme, the probability encoding scheme ensures the unbiasedness in the whole encoding and decoding process in a probability form, thereby improving the encoding and decoding precision.

[0083] ​​​The present application is based on the security problem in the data transmission process, considers the random spoofing attack, and ensures the security of the transmitted data when encountering network attacks. In addition, the traditional method usually ignores the disturbance of the estimator, but this often cannot be ignored in actual engineering, so the present application designs a flexible estimator to overcome external disturbance.

[0084] The present application comprehensively and deeply considers the influence of the probability encoding mechanism, random spoofing attack and estimator gain mutation, improves the performance of state estimation, and has the advantages of easy online solution and implementation. BRIEF DESCRIPTION OF DRAWINGS

[0085] The present application will be further described in detail below in combination with the drawings and specific embodiments.

[0086] Figure 1 is a flow chart of the power system distributed estimation method described in the present application;

[0087] Figure 2 is a schematic diagram of the sensor network topology used in the present application;

[0088] Figure 3 is a schematic diagram of the attack sequence of the four sensor nodes in the power system;

[0089] Figure 4 is a schematic diagram of the actual rotor angle and the estimation of the first sensor node and the estimation of the second sensor node in the power system;

[0090] Figure 5 is a schematic diagram of the actual rotor angle and the estimation of the third sensor node and the estimation of the fourth sensor node in the power system;

[0091] Figure 6 is a schematic diagram of the actual rotor angular velocity and the estimation of the first sensor node and the estimation of the second sensor node in the power system;

[0092] Figure 7 is a schematic diagram of the actual rotor angular velocity and the estimation of the third sensor node and the estimation of the fourth sensor node in the power system;

[0093] Figure 8 is a schematic diagram of the actual q-axis transient voltage and the estimation of the first sensor node and the estimation of the second sensor node in the power system;

[0094] Figure 9 is a schematic diagram of the actual q-axis transient voltage and the estimation of the third sensor node and the estimation of the fourth sensor node in the power system;

[0095] Figure 10is a schematic diagram of the actual d-axis transient voltage and the estimation of the first sensor node and the estimation of the second sensor node in the power system;

[0096] Figure 11 is a schematic diagram of the actual d-axis transient voltage and the estimation of the third sensor node and the estimation of the fourth sensor node in the power system;

[0097] Figure 12 is a schematic diagram of the logarithm of the mean square error of the estimated value of the system state in the power system and the logarithm of the upper bound of the estimation error covariance corresponding thereto. DETAILED DESCRIPTION

[0098] The present application will be further described below in conjunction with the accompanying drawings and specific embodiments. It should be understood that these embodiments are only used to illustrate the present application and are not used to limit the scope of the present application. In addition, it should be understood that those skilled in the art can make various modifications or changes to the present application after reading the content taught by the present application, and these equivalent forms also fall within the scope defined by the appended claims of the present application.

[0099] Embodiment 1:

[0100] Please refer to Figure 1 Embodiments of the present application relate to a power system distributed estimation method under probability coding and fraud attack, and a flowchart of the power system distributed estimation method is as shown in Figure 1 , which comprises the following steps:

[0101] Step 1: Based on the smart grid model, the rotor angle, the rotor angular velocity, the d-axis transient voltage and the q-axis transient voltage of the power system are selected as the state variables, and a power system dynamic model is established.

[0102] The calculation formula of the power system dynamic model is as follows:

[0103]

[0104] In the formula, is a random nonlinear function based on the smart grid model, wherein x s =[δ s ω s E′ q,s E′ d,s ] T is the state vector of the power system at time s, wherein δ s is the rotor angle of the power system at time s,

[0105] ω s is the rotor angular velocity of the power system at time s, E′ q,s is the q-axis transient voltage of the power system at time s, and E′ d,sis the d-axis transient voltage of the power system at time s.

[0106] μ s = [v s θ s T m,s E fd,s ] T is the known input vector of the power system at time s, where v s is the terminal bus voltage magnitude of the power system at time s, θ S is the terminal bus voltage phase angle of the power system at time s, T M,s is the mechanical torque input of the power system at time s, E fd,s is the damping coefficient of the power system at time s.

[0107] x s+1 is the state vector of the power system at time s+1.

[0108] ε s is the process noise with zero mean and variance R s > 0 at time s.

[0109] z i,s is the actual measurement of the i-th sensor at time s.

[0110] g i,s = [f i,s P i,s Q i,s ] T is the ideal measurement of the i-th sensor at time s, where f i,s is the terminal frequency measured by the i-th sensor node of the power system at time s, P i,s is the active power measured by the i-th sensor node of the power system at time s, Q i,s is the reactive power measured by the i-th sensor node of the power system at time s.

[0111] is the measurement noise with zero mean and covariance W k > 0 in the measurement of the i-th sensor node at time s.

[0112] Step 2: Under the probability encoding mechanism and random spoofing attack, the design of the resilient distributed state estimator is completed according to the dynamic model of the power system in step 1, and the specific steps are as follows:

[0113] Step 2.1 Based on the probability encoding scheme, the probability encoding-decoding measurement value is calculated according to the measurement value of the i-th sensor node in step 1, and the calculation formula of the probability encoding-decoding measurement value is as follows:

[0114]

[0115] wherein, is the decoded measurement of the ith sensor node at time s;

[0116] represents the decoding error of the ith sensor node at time s under the probabilistic encoding mechanism, and the specific calculation formula is as follows:

[0117]

[0118] wherein, is the encoding probability;

[0119] is the encoding level, wherein k represents the numerical description corresponding to the encoding level;

[0120] represents the encoding interval of the encoder;

[0121] Step 2.2 considers the random occurrence of spoofing attacks in network communication, and calculates the corresponding transmission measurement of the ith sensor node according to the decoded measurement of the ith sensor node obtained in step 2.1, that is, the output measurement value finally received by the estimator, and the calculation formula is as follows:

[0122]

[0123] wherein, represents the malicious signal sent by the attacker to the ith sensor node at time s;

[0124] r i,s satisfies is a known scalar;

[0125] η i,s is the probability of spoofing attack of the ith sensor node, and the mathematical expectation is

[0126] Step 2.3 designs the probabilistic encoding mechanism and the resilient distributed state estimator in the dynamic model of the power system under random spoofing attack according to the output measurement value finally received by the estimator obtained in step 2.2, and the calculation formula of the resilient distributed state estimator is as follows:

[0127]

[0128] wherein, represents the one-step prediction value of the state vector x s+1 of the ith sensor node at time s;

[0129] respectively represent the one-step prediction value of the state vector x sThe estimated value;

[0130] These represent the state vector x of the i-th sensor node at time s+1. s+1 The estimated value;

[0131] K i,s+1 The distributed state estimator gain represents the value of the i-th sensor node at time s+1.

[0132] The probability η of the i-th sensor node being subjected to a spoofing attack. i,s The mathematical expectation;

[0133] ΔK i,s+1 =ξ i,s C i Represents the abrupt change in the distributed state estimator gain of the i-th sensor node at time s+1, where C i Given a matrix ξ i,s It is a symmetric property with zero mean and variance of 0. Gaussian white noise sequence;

[0134] This represents the predicted value of the measurement of the i-th sensor node at time s;

[0135] I represents a known identity matrix of appropriate dimensions;

[0136] a ij This represents the connection coefficient between the i-th sensor node and the j-th sensor node;

[0137] This represents the set of adjacent nodes of the i-th sensor node.

[0138] Step 3: Based on the elastic distributed state estimator in Step 2, calculate the one-step prediction matrix of the i-th sensor node at time s. Therefore, the upper bound of the one-step prediction error covariance matrix ∑ for the i-th sensor node at time s is calculated. i,s+1|s The calculation formula is as follows:

[0139]

[0140] In the formula, ∑ i,s|s It is the upper bound of the estimation error covariance matrix of the i-th sensor node at time s;

[0141] A i,s yes right Point to ask about x s The matrix obtained by the partial derivatives;

[0142] U s and Ls It is a matrix of known appropriate dimensions;

[0143] λ1 is a known positive scalar. It is the reciprocal of λ1;

[0144] R s It is the covariance matrix of the process noise at time s;

[0145] It is Σ i,s|s The inverse matrix;

[0146] They are A i,s U s L s The transpose of .

[0147] Step 4: Based on the upper bound Σ of the one-step prediction error covariance matrix obtained in Step 3. i,s+1|s Calculate the elastic distributed estimator gain matrix K of the i-th sensor node at time s+1. i,s+1 The specific calculation formula is as follows:

[0148]

[0149] In the formula,

[0150]

[0151] In the formula, ρ1 is the first intermediate parameter in the gain matrix;

[0152] ρ2 is the second intermediate parameter in the gain matrix;

[0153] ρ3 is the third intermediate parameter in the gain matrix;

[0154] λ² is a known positive scalar. It is the reciprocal of λ²;

[0155] B i,s+1 It is g i,s right Point to ask about x s The matrix obtained by the partial derivatives;

[0156] F s+1 and G s+1 It is a matrix of known appropriate dimensions;

[0157] O i It is an intermediate matrix;

[0158] ψ and φ are intermediate matrices generated during the process of dealing with deception attacks;

[0159] W i,s+1 is the covariance matrix of the process noise of the i-th sensor node at s+1 moment;

[0160] are the transpose matrices of B i,s+1 , F s+1 , and G s+1 , respectively; O i

[0161] is the inverse matrix of ∑ s+1|s

[0162] represents the Hadamard product operation.

[0163] Step 5: Substitute the estimator gain matrix K i,s+1 obtained in step 4 into the elastic distributed state estimator in step 2 to calculate the state estimation matrix X of the i-th sensor node at s+1 moment.

[0164] Step 6: According to the estimator gain matrix K i,s+1 obtained in step 4, calculate the estimation error covariance matrix upper bound ∑ i,s+1|s+1 of the i-th sensor node at s+1 moment, and set s = s+1, execute step 2 until s+1 reaches the total time T. i,s+1|s+1 The specific calculation formula is as follows:

[0165]

[0166] In the formula, ∑ i,s+1|s+1 is the estimation error covariance upper bound of the i-th sensor node at s+1 moment;

[0167] are the transpose matrices of K i,s+1 and ΔK i,s+1 , respectively.

[0168] Embodiment 2:

[0169] In this embodiment, the estimation method provided by the present application is used to perform simulation experiments on the IEEE 39-bus New England power system.

[0170] Selected parameters: the initial state vector of the power system is x 0|0 = [0.4100] T ​​where the first component, the second component, the third component and the fourth component represent rotor angle, rotor angular velocity, q-axis transient voltage and d-axis transient voltage respectively;

[0171] The mean values of the system noise and the measurement noise are both 0, and the covariance is R k = 10 -4 I and Q k = 10 -4 I;

[0172] The encoding interval in the probabilistic encoding scheme is

[0173] The occurrence probabilities of the four sensor nodes cheating attacks are η 1,s = 0.3, η 2,s = 0.21, η 3,s = 0.2 and η 4,s = 0.18 respectively;

[0174] The upper bound of the Euclidean norm of the attack signal is

[0175] The occurrence probability ξ i,s of the gain disturbance satisfies the mean value of 0 and the variance

[0176] The intermediate matrix C i in the gain disturbance is expressed as follows:

[0177]

[0178] Other related parameters are set as U s = I, L s = 0.01I, F s = 0.01I, G s = 0.01I respectively.

[0179] In addition, the mean square error is used to further measure the performance of the elastic distributed state estimation method provided by the application, and the specific process is as follows:

[0180]

[0181] Please refer to Figures 3-12 , and the experimental results are shown in Figures 3 to 12 , wherein Figure 2 is Figure 3 is a schematic diagram of attack sequence of four sensor nodes cheating attacks in a power system, wherein the vertical coordinates 0, 1, 2 and 3 positions represent the attack sequence of the four sensor nodes respectively;

[0182] Figures 4 to 11is an estimation schematic diagram of actual rotor angle, actual rotor angular velocity, actual q-axis transient voltage, actual d-axis transient voltage and four sensor nodes in the power system, wherein the solid line is the actual value and the dashed line is the corresponding estimated value;

[0183] Figure 12 is a schematic diagram of the logarithm of the estimation value mean square error of the system state in the power system and the logarithm of the corresponding estimation error covariance upper bound, wherein the solid line is the logarithm of the estimation value mean square error and the dashed line is the logarithm of the corresponding estimation error covariance upper bound.

[0184] Based on the above simulation results, it can be verified that the method proposed in the application is an effective solution to the problem of smart power system state estimation under the probability encoding mechanism and random deception attack.

[0185] According to the comprehensive results of embodiment 2, the probability encoding mechanism and the flexible distributed state estimation method applied to the smart power system under the random deception attack adopted by the application solve the problem that the existing state estimation method cannot effectively cope with the smart power system state estimation under the probability encoding mechanism and the random deception attack, and has great practical application value.

Claims

1. A distributed estimation method for power systems under probabilistic coding and spoofing attacks, characterized in that, The method includes the following steps: Step 1: Establishment of the power system dynamic model Based on the smart grid model, select the rotor angle, rotor angular velocity, d-axis transient voltage, and q-axis transient voltage of the power system as state variables, and establish the power system dynamic model; Step 2: Design of the resilient distributed state estimator Under the probability coding mechanism and random deception attacks, complete the design of the resilient distributed state estimator according to the power system dynamic model in Step 1; Step 3: Calculate ∑ i,s+1|s Based on the elastic distributed state estimator in step 2, calculate the one-step prediction matrix of the i-th sensor node at time s. Therefore, the upper bound of the one-step prediction error covariance matrix ∑ for the i-th sensor node at time s is calculated. i,s+1|s ; Step 4: Calculate K i,s+1 Based on the upper bound of the one-step prediction error covariance matrix obtained in step 3, ∑ i,s+1|s Calculate the elastic distributed estimator gain matrix K of the i-th sensor node at time s+1. i,s+1 ; Step 5: Calculation The estimator gain matrix K obtained in step 4 i,s+1 Substitute this into the elastic distributed state estimator in step 2 to calculate the state estimation matrix of the i-th sensor node at time s+1. Judge whether s + 1 reaches the total duration T. If s + 1 < T is satisfied, continue to execute Step 6, otherwise end the operation; Step 6: Calculate Σ i,s+1|s+1 Based on the estimator gain matrix K obtained in step 4 i,s+1 Calculate the upper bound of the estimation error covariance matrix of the i-th sensor node at time s+1, ∑ i,s+1|s+1 Meanwhile, let s = s + 1, and execute step 2 until s + 1 reaches the total duration T.

2. The method according to claim 1, characterized in that, The power system dynamic model in Step 1 is calculated as follows: In the formula, It is a stochastic nonlinear function established based on the smart grid model, where x in the formula is... s =[δ s ω s E′ q,s E′ d,s ] T Let δ be the state vector of the power system at time s, where δ s ω is the rotor angle of the power system at time s. s E′ is the rotor angular velocity of the power system at time s. q,s It is the q-axis transient voltage of the power system at time s, E′ d,s It is the d-axis transient voltage of the power system at time s, where μ is the voltage at time s. s =[V s θ s T m,s E fd,s ] T V is the known input vector of the power system at time s, where V s θ is the amplitude of the terminal bus voltage at time s in the power system. s T is the phase angle of the terminal bus voltage at time s in the power system. m,s It is the mechanical torque input of the power system at time s, E fd,s It is the damping coefficient of the power system at time s; x s+1 It is the state vector of the power system at time s+1; ε s It is zero mean and variance R at time s s Process noise >0; z i,s It is the actual measurement value of the i-th sensor at time s; g i,s =[f i,s P i,s Q i,s ] T It is the ideal measurement value of the i-th sensor at time s, where f i,s P is the terminal frequency measured by the i-th sensor node at time s in the power system. i,s Q is the active power measured by the i-th sensor node at time s in the power system. i,s It is the reactive power measured by the i-th sensor node at time s in the power system; The measurement process of the i-th sensor node at time s satisfies zero mean and covariance W. k Measurement noise >0.

3. The method according to claim 2, characterized in that, The specific steps for the design of the resilient distributed state estimator in Step 2 are as follows: Step 2.1: Based on the probability coding scheme, calculate the probability coding-decoding measurement value according to the measurement value of the i-th sensor node in Step 1. The calculation formula for the probability coding-decoding measurement value is as follows: In the formula, It is the decoded measurement value of the i-th sensor node at time s; This represents the decoding error corresponding to the i-th sensor node at time s under the probabilistic coding mechanism. The specific calculation formula is as follows: In the formula, It is the encoded probability; It represents the coding level, where k represents the numerical description corresponding to the coding level; Represents the encoder's encoding interval; Step 2.2 Considering the randomly occurring deception attacks in network communication, calculate the corresponding transmission measurement of the i-th sensor node according to the decoded measurement value of the i-th sensor node obtained in Step 2.1, that is, the output measurement value finally received by the estimator. The calculation formula for the output measurement value finally received by the estimator is as follows: In the formula, Let r represent the malicious signal emitted by the attacker at time s to the i-th sensor node, where r i,s satisfy It is a known scalar; η i,s Let be the probability that the i-th sensor node is subjected to a spoofing attack, and its expected value is . Step 2.3 According to the output measurement value finally received by the estimator obtained in Step 2.2, design the resilient distributed state estimator in the power system dynamic model under the probability coding mechanism and random deception attacks. The calculation formula for the resilient distributed state estimator is as follows: In the formula, x represents the state vector of the i-th sensor node at time s. s+1 One-step prediction value; These represent the state vector x of the i-th sensor node at time s, respectively. s The estimated value; These represent the state vector x of the i-th sensor node at time s+1. s+1 The estimated value; K i,s+1 The distributed state estimator gain represents the value of the i-th sensor node at time s+1. ΔK i,s+1 =ξ i,s C i This represents the sudden change in the distributed state estimator gain of the i-th sensor node at time s+1. Where, C i Given a matrix ξ i,s It is a symmetric property with zero mean and variance of 0. Gaussian white noise sequence; The probability η of the i-th sensor node being subjected to a spoofing attack. i,s The mathematical expectation; This represents the predicted value of the measurement of the i-th sensor node at time s; I represents the identity matrix of a known appropriate dimension; a ij This represents the connection coefficient between the i-th sensor node and the j-th sensor node; This represents the set of adjacent nodes of the i-th sensor node.

4. The method according to claim 3, characterized in that, In step 3, the upper bound of the one-step prediction error covariance matrix at time s is ∑ i,s+1|s The calculation formula is as follows: In the formula, Σ i,s|s It is the upper bound of the estimation error covariance matrix of the i-th sensor node at time s; It is Σ i,s|s The inverse matrix; A i,s yes right Point to ask about x s The matrix obtained by the partial derivatives; U s and L s It is a matrix of known appropriate dimensions; λ1 is a known positive scalar. It is the reciprocal of λ1; R s It is the covariance matrix of the process noise at time s; They are A i,s U s L s The transpose of .

5. The method according to claim 4, characterized in that, The gain matrix K in step 4 i,s+1 The specific calculation formula is as follows: In the formula, In the formula, ρ1 is the first intermediate parameter in the gain matrix; ρ2 is the second intermediate parameter in the gain matrix; ​ λ² is a known positive scalar. It is the reciprocal of λ²; B i,s+1 It is g i,s right Point to ask about x s The matrix obtained by the partial derivatives; F s+1 and G s+1 It is a matrix of known appropriate dimensions; O i It is an intermediate matrix; ψ and It is an intermediate matrix generated during the process of dealing with deception attacks; W i,s+1 It is the covariance matrix of the process noise of the i-th sensor node at time s+1; They are B i,s+1 F s+1 G s+1 , O i The transpose of the matrix; It is ∑ s+1|s The inverse matrix; ​ 6. The method according to claim 5, characterized in that, The upper bound of the estimation error covariance matrix in step 6, ∑ i,s+1|s+1 The specific calculation formula is as follows: In the formula, ∑ i,s+1|s+1 It is the upper bound of the estimation error covariance of the i-th sensor node at time s+1; They are K i,s+1 ΔK i,s+1 The transpose of .

Images