Power system distributed estimation method under probability coding and spoofing attack

Abstract

The invention discloses an electric power system distributed estimation method under probability coding and spoofing attack. The method comprises the following steps: step 1, establishing an electric power system dynamic model; 2, completing the design of an elastic distributed state estimator; 3, calculating an upper bound sigma i, s + 1s of a one-step prediction error covariance matrix of the ith node at the time s; 4, calculating an elastic distributed estimator gain matrix Ki, s + 1 of the ith node; step 5, calculating a state estimation matrix # imgabs0 # of the ith node at the s + 1 moment; step 6, calculating an upper bound sigma i, s + 1s + 1 of an estimation error covariance matrix of the ith node at the s + 1 moment according to the estimator gain matrix Ki, s + 1 obtained in the step 4. The estimation method provided by the invention solves the problem that the existing state estimation method is difficult to process the occurrence probability coding and spoofing attack of the power system.

Description

Technical Field

[0001] The invention relates to the technical field of power network system control, and in particular to a distributed estimation method for a power system under probability coding and deception attacks. Background Art

[0002] With the development of the economy and society, customers' demand for electricity continues to grow and diversify. To this end, smart grids improve system flexibility and efficiency by integrating renewable energy, energy storage technology, advanced electronic components and intelligent technology. However, the high degree of integration has led to 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 the state of synchronous motors is essential. The wide-area monitoring system collects power data in real time through phasor measurement units, and then uses dynamic state estimation technology to accurately estimate the state of the generator and improve the system monitoring capability. However, due to the nonlinear characteristics of the power system, more advanced estimation algorithms need to be developed to meet the challenges of complex systems and ensure the safe, stable and economical operation of the power system.

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

[0004] At present, there is still a lack of research on effectively integrating distributed solutions into probability coding mechanisms and coping with random deception attacks in the field of state estimation of smart power systems, especially the resilient design of estimators under potential disturbances. This means that in the face of complex and changing power system environments, traditional state estimation and filtering methods may appear to be inadequate and difficult to ensure the continuous stability of system performance. Innovative solutions are urgently needed to fill this technical gap. Summary of the invention

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

[0006] In order to solve the above technical problems, the specific technical solutions of the present invention 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 attack, 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 Thus, 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 , and at the same time, 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] In the formula, It is a random nonlinear function based on the smart grid model.

[0022] Among them, x s =[δ s ω s E′ q,s E′ d,s ] T is the state vector of the power system at time s, where δ 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, where 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 zero mean at time s and has variance R s Process noise >0;

[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 ith sensor at time s, where f i,s is the terminal frequency measured by the ith sensor node in the power system at time s, P i,s is the active power measured by the i-th sensor node in the power system at time s, Q i,s is the reactive power measured by the i-th sensor node in the power system at time s;

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

[0028] Preferably, the specific steps of designing the recursive state estimator in step 2 are as follows:

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

[0030]

[0031] In the formula, is the decoded measurement value of the i-th sensor node at time s;

[0032] represents the decoding error corresponding to the i-th sensor node at time s under the probability coding mechanism. The specific calculation formula is as follows:

[0033]

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

[0035] Step 2.2 considers the random deception attack in the network communication. According to the decoded measurement value of the i-th sensor node obtained in step 2.1, the corresponding transmission measurement of the i-th sensor node is calculated, that is, the output measurement value finally received by the estimator. The calculation formula of the output measurement value finally received by the estimator is as follows:

[0036]

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

[0038] η i,s is the probability of a spoofing attack on the i-th sensor node, and its mathematical expectation is

[0039] Step 2.3 designs a probabilistic coding mechanism and a resilient distributed state estimator in a dynamic model of a power system under a random deception attack according to the output measurement value finally received by the estimator obtained in step 2.2. The calculation formula of the resilient distributed state estimator is as follows:

[0040]

[0041] In the formula, Represents the state vector x of the i-th sensor node at time s s+1 One-step forecast value of

[0042] They represent the state vector x of the i-th sensor node at time s. s An estimated value of

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

[0044] K i,s+1 represents the distributed state estimator gain 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, where, i is a known matrix, ξ i,s is a zero-mean and a Gaussian white noise sequence;

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

[0047] is the probability that the ith sensor node will be attacked by a spoofing attack. i,s The mathematical expectation of

[0048] I represents the identity matrix of known appropriate dimensions;

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

[0050] Represents the set of neighboring nodes of the i-th sensor node.

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

[0052]

[0053] In the formula, Σ i,s|s 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 Click to ask about x s The matrix obtained by the partial derivative of

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

[0056] λ 1 is a known positive scalar, is λ 1 The reciprocal of

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

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

[0059] A i,s , U s , L s The transposed matrix of .

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

[0061]

[0062] In the formula,

[0063]

[0064] Among them, ρ 1 It is the first intermediate parameter in the gain matrix;

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

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

[0067] λ 2 is a known positive scalar, is λ 2 The reciprocal of

[0068] B i,s+1 Yes i,s right Click to ask about x s The matrix obtained by the partial derivative of

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

[0070] O i is an intermediate matrix;

[0071] ψ and φ are the intermediate matrices generated during the processing of spoofing attacks;

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

[0073] B i,s+1 、F s+1 , G s+1 , O i The transposed matrix of

[0074] is Σ s+1|s The inverse matrix of

[0075] Stands for Hadamard Product operation.

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

[0077]

[0078] 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;

[0079] K i,s+1 , ΔK i,s+1 The transposed matrix of .

[0080] Beneficial Effects

[0081] Due to the adoption of the above technical solution, the present invention has the following advantages and positive effects compared with the prior art:

[0082] The elastic distributed state estimation method of the present invention comprehensively considers complex situations such as probabilistic coding mechanism, random deception attack and estimator gain mutation. Compared with traditional coding schemes, the probabilistic coding scheme ensures unbiasedness in the entire encoding and decoding process in a probabilistic form, thereby improving the encoding and decoding accuracy.

[0083] Based on the security problem in the data transmission process, the present invention takes into account the random deception attacks, 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 cannot be ignored in practical engineering. Therefore, the present invention designs a flexible estimator to overcome the external disturbance.

[0084] The present invention comprehensively and deeply considers the influence of probability coding mechanism, random deception attack and sudden change of estimator gain, improves the performance of state estimation, and has the advantage of being easy to solve and realize online. BRIEF DESCRIPTION OF THE DRAWINGS

[0085] The present invention will be further described in detail below in conjunction with the accompanying drawings and specific embodiments.

[0086] Figure 1 is a flow chart of the distributed estimation method for a power system according to the present invention;

[0087] Figure 2 It is a schematic diagram of the sensor network topology structure adopted by the present invention;

[0088] Figure 3 It is a schematic diagram of the attack sequence of a deception attack on four sensor nodes in a 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 in the power system and the estimation of the first sensor node and the estimation of the second sensor node;

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

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

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

[0097] Fig.12 It 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 covariance of the corresponding estimation error. DETAILED DESCRIPTION

[0098] The present invention 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 invention and are not intended to limit the scope of the present invention. In addition, it should be understood that after reading the content taught by the present invention, those skilled in the art may make various changes or modifications to the present invention, and these equivalent forms fall within the scope limited by the appended claims of the application equally.

[0099] Embodiment 1:

[0100] Please refer to Figure 1 The embodiment of the present invention relates to a distributed estimation method for a power system under probability coding and deception attack. The flow chart of the distributed estimation method for a power system is as follows: Figure 1 As shown, the following steps are included:

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

[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, where x s =[δ s ω s E′ q,s E′ d,s ] T is the state vector of the power system at time s, where δ sis 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, E′ d,s is 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 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.

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

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

[0109] z i,s is the actual measurement value 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 value of the i-th sensor at time s, where f i,s is the terminal frequency measured by the ith sensor node in the power system at time s, P i,s is the active power measured by the i-th sensor node in the power system at time s, Q i,s is the reactive power measured by the i-th sensor node in the power system at time s.

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

[0112] Step 2: Under the probabilistic coding mechanism and random deception attack, the design of the elastic distributed state estimator is completed according to the dynamic model of the power system in step 1. The specific steps are as follows:

[0113] Step 2.1 is based on the probability coding scheme. The probability coding-decoding measurement value is calculated according to the measurement value of the i-th sensor node in step 1. The calculation formula of the probability coding-decoding measurement value is as follows:

[0114]

[0115] In the formula, is the decoded measurement value of the i-th sensor node at time s;

[0116] represents the decoding error corresponding to the i-th sensor node at time s under the probability coding mechanism. The specific calculation formula is as follows:

[0117]

[0118] In the formula, is the encoding probability;

[0119] is the coding level, where k represents the numerical description corresponding to the coding level;

[0120] Represents the encoding interval of the encoder;

[0121] Step 2.2 considers the random deception attack in the network communication. According to the decoded measurement value of the i-th sensor node obtained in step 2.1, the corresponding transmission measurement of the i-th sensor node is calculated, that is, the output measurement value finally received by the estimator. The calculation formula is as follows:

[0122]

[0123] In the formula, represents the malicious signal sent by the attacker to the i-th sensor node at time s;

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

[0125] η i,s is the probability of a deception attack on the i-th sensor node, and its mathematical expectation is

[0126] Step 2.3: Based on the output measurement value finally received by the estimator obtained in step 2.2, a probability coding mechanism and a resilient distributed state estimator in the dynamic model of the power system under random deception attack are designed. The calculation formula of the resilient distributed state estimator is as follows:

[0127]

[0128] In the formula, Represents the state vector x of the i-th sensor node at time s s+1 One-step forecast value of

[0129] Represents the state vector x of the i-th sensor node at time s s An estimated value of

[0130] They represent the state vector x of the i-th sensor node at time s+1 respectively. s+1 An estimated value of

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

[0132] is the probability that the ith sensor node will be attacked by a spoofing attack. i,s The mathematical expectation of

[0133] Δ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, where C i is a known matrix, ξ i,s is a zero-mean and a Gaussian white noise sequence;

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

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

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

[0137] Represents the set of neighboring 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 Thus, the upper bound of the one-step prediction error covariance matrix of 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 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 Click to ask about x s The matrix obtained by the partial derivative of

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

[0143] λ 1 is a known positive scalar, is λ 1 The reciprocal of

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

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

[0146] A i,s , U s , L s The transposed matrix of .

[0147] Step 4: According to the upper bound of the one-step forecast 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 calculation formula is as follows:

[0148]

[0149] In the formula,

[0150]

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

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

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

[0154] λ 2 is a known positive scalar, is λ 2 The reciprocal of

[0155] B i,s+1 Yes i,s right Click to ask about xs The matrix obtained by taking the partial derivative of

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

[0157] O i is an intermediate matrix;

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

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

[0160] are respectively the transpose matrices of B i,s+1 、F s+1 、G s+1 、 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 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, then continue to execute Step 6; otherwise, end the operation;

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

[0165]

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

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

[0168] Embodiment 2:

[0169] This embodiment takes the IEEE 39 bus New England power system as an example and uses the estimation method provided by the present invention to perform simulation experiments.

[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 the rotor angle, the rotor angular velocity, the q-axis transient voltage and the d-axis transient voltage respectively;

[0171] The mean of system noise and measurement noise are both 0, and their covariances are R k =10 -4 I and Q k =10 -4 I;

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

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

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

[0175] The probability of gain disturbance occurring ξ i,s Satisfies the mean and variance of 0

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

[0177]

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

[0179] In addition, the present invention uses mean square error to further measure the performance of the elastic distributed state estimation method proposed by the present invention, as follows:

[0180]

[0181] Please refer to Figure 3-Figure 12 The experimental results are as follows Figures 3 to 12 As shown, Figure 2 yes Figure 3 It is a schematic diagram of the attack sequence of deception attacks on four sensor nodes in the power system, where the vertical coordinates 0, 1, 2, and 3 respectively represent the attack sequences of the four sensor nodes;

[0182] Figures 4 to 11 is a schematic diagram of the actual rotor angle, actual rotor angular velocity, actual q-axis transient voltage, actual d-axis transient voltage, and estimation of four sensor nodes in the power system, where the solid lines are actual values ​​and the dotted lines are the corresponding estimated values;

[0183] Fig.12 It 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 corresponding estimated error covariance, where the solid line is the logarithm of the mean square error of the estimated value and the dotted line is the logarithm of the upper bound of the corresponding estimated error covariance.

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

[0185] According to the comprehensive results of Example 2, the present invention adopts a probability coding mechanism and a resilient distributed state estimation method for intelligent power systems under random deception attacks, which solves the problem that existing state estimation methods cannot effectively cope with the probability coding mechanism and the state estimation of intelligent power systems under random deception attacks, and has great practical application value.

Claims

1. A distributed estimation method for power system under probability coding and deception attack, 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 to establish the power system dynamic model; Step 2: Design of the resilient distributed state estimator Under the probability encoding 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 According to the elastic distributed state estimator in step 2, calculate the one-step prediction matrix of the i-th sensor node at time s Thus, the upper bound of the one-step prediction error covariance matrix of the i-th sensor node at time s is calculated ∑ i,s+1|s ; Step 4: Calculate K i,s+1 According to the upper bound of the one-step prediction error covariance matrix ∑ 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 it into the elastic distributed state estimator in step 2 and 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 According to the estimator gain matrix K obtained in step 4 i,s+1 , calculate the upper bound of the estimated error covariance matrix of the i-th sensor node at time s+1∑ i,s+1|s+1 At the same time, 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, is a random nonlinear function based on the smart grid model, where x s =[δ s ω s E′ q,s E′ d,s ] T is the state vector of the power system at time s, where δ 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, where μ 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 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; x s+1 is the state vector of the power system at time s+1; ε s is zero mean at time s and has variance R s Process noise >0; z i,s 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 is the ideal measurement value of the ith sensor at time s, where f i,s is the terminal frequency measured by the ith sensor node in the power system at time s, P i,s is the active power measured by the i-th sensor node in the power system at time s, Q i,s is the reactive power measured by the i-th sensor node in the power system at time s; is the value of zero mean and covariance W in the measurement process of the i-th sensor node at time s k >0 measurement noise.

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 encoding scheme, calculate the probability encoding-decoding measurement value according to the measurement value of the i-th sensor node in Step 1. The calculation formula for the probability encoding-decoding measurement value is as follows: In the formula, is the decoded measurement value of the i-th sensor node at time s; represents the decoding error corresponding to the i-th sensor node at time s under the probability coding mechanism. The specific calculation formula is as follows: In the formula, is the encoding probability; is the coding level, and k represents the numerical description corresponding to the coding level; Represents the encoding interval of the encoder; 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, represents the malicious signal sent by the attacker to the i-th sensor node at time s, where r i,s satisfy is a known scalar; η i,s is the probability of a deception attack on the i-th sensor node, and its mathematical expectation is Step 2.3: Design the resilient distributed state estimator in the power system dynamic model under the probability encoding mechanism and random deception attacks according to the output measurement value finally received by the estimator obtained in Step 2.

2. The calculation formula for the resilient distributed state estimator is as follows: In the formula, Represents the state vector x of the i-th sensor node at time s s+1 One-step forecast value of They represent the state vector x of the i-th sensor node at time s. s An estimated value of They represent the state vector x of the i-th sensor node at time s+1 respectively. s+1 An estimated value of K i,s+1 represents the distributed state estimator gain of the i-th sensor node at time s+1; Δ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, Among them, C i is a known matrix, ξ i,s is a zero-mean and a Gaussian white noise sequence; is the probability η that the ith sensor node will be attacked by a spoofing attack i,s The mathematical expectation of represents the predicted value of the measurement value of the i-th sensor node at time s; I represents the identity matrix of a known appropriate dimension; a ij represents the connection coefficient between the i-th sensor node and the j-th sensor node; Represents the set of neighboring nodes of the i-th sensor node.

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

5. The method according to claim 4, characterized in that The gain matrix K in step 4 i,s+1 The 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; ρ3 is the third intermediate parameter in the gain matrix; λ2 is a known positive scalar, is the reciprocal of λ2; B i,s+1 Yes i,s right Click to ask about x s The matrix obtained by the partial derivative of F s+1 and G s+1 is a matrix of known appropriate dimensions; O i is an intermediate matrix; ψ and φ are the intermediate matrices generated during the processing of spoofing attacks; W l,s+1 is the covariance matrix of the process noise of the i-th sensor node at time s+1; B i,s+1 、F s+1 , G s+1 , O i The transposed matrix of Yes s+1|s The inverse matrix of "°” represents the Hadamard product operation.

6. The method according to claim 5, characterized in that The upper bound ∑ of the estimated error covariance matrix in step 6 is 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; K i,s+1 , ΔK i,s+1 The transposed matrix of .

Images