A method, system and device for environment excited subsynchronous oscillation mode identification

Abstract

This invention discloses a method, system, and device for identifying subsynchronous oscillation modes under environmental excitation, belonging to the field of power system oscillation detection. Kalman filtering (KF) is used to estimate the coefficients of an autoregressive (AR) model to fit measured environmental data. The damping factor and frequency of the subsynchronous oscillation (SSO) mode can be directly obtained by solving for the roots of the characteristic polynomial corresponding to the AR model. Furthermore, multi-model segmentation filtering (MMPF) is applied to KF to select the correct AR model order. This invention combines the KF algorithm with MMPF, enabling the algorithm to directly identify weak SSO mode parameters from environmental data. The method has simple operation steps and short response time, thus possessing good application prospects.

Description

Technical Field

[0001] This invention relates to the field of power system oscillation detection, and specifically to a method, system, and device for identifying subsynchronous oscillation modes under environmental excitation. Background Technology

[0002] With the development of power systems, the proportion of renewable energy generation has surged, leading to increasingly frequent subsynchronous oscillations (SSOs). Type 3 wind turbines (WTGs) are among the most widely used renewable energy generation devices, and wind farms composed of these turbines have experienced numerous SSO events globally. SSOs pose a significant threat to the safety of electrical equipment and the stability of the power system, potentially causing severe damage to wind turbines, leading to large-scale tripping, and ultimately resulting in major power outages. Identifying SSO modes is crucial for early warning of SSOs and the effective development of protection and control strategies; therefore, accurate and timely identification of SSO modes is of great importance.

[0003] Existing measured SSO response signals can be divided into two types: disturbance excitation response and operating environment noise excitation response. Under disturbance excitation, the SSO component content is large, and it can be easily detected using methods such as Prony and Discrete Fourier Transform. Operating environment noise excitation is a continuous noise disturbance, such as random fluctuations in energy output and load switching. Under this condition, the SSO periodic amplitude is small and the signal-to-noise ratio is extremely low, making it difficult to directly identify the SSO modal parameters using the aforementioned common identification methods. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention proposes a method, system, and device for identifying subsynchronous oscillation modes under environmental excitation, solving the problem that SSO mode parameters are difficult to identify under environmental noise excitation.

[0005] The objective of this invention can be achieved through the following technical solutions:

[0006] Firstly, this application proposes a method for identifying subsynchronous oscillation modes under environmental excitation, comprising the following steps:

[0007] An AR (Autoregressive) model is constructed based on power system data under SSO mode.

[0008] Calculate the parameters of AR models at different orders using Kalman filtering (KF).

[0009] The order of the AR model is determined using the MMPF (Multi-Model Partition Filter).

[0010] The parameters and order of the AR model are used to calculate the eigenvalues ​​through the characteristic equation;

[0011] Damping and frequency in the SSO mode are determined by eigenvalues.

[0012] In some embodiments, calculating the parameters of AR models at different orders using KF includes the following steps:

[0013] Convert the AR model into a state-space equation:

[0014] x p (k)=x p (k-1) (2)

[0015] y(k)=H p (k)x p (k)+v(k) (3)

[0016] Where y(k) is the measurement signal at time k; x p (k) is the state vector, x p (k)=[a1(k) a2(k) … a p (k)] T ; a1(k), a2(k), ..., a p (k) represents the time-varying coefficients of the AR model; H p (k) is the observation matrix, H p (k) = [y(k-1) y(k-2) … y(kp)]; v(k) is white noise;

[0017] Initialization parameters:

[0018]

[0019] In the formula, P p (1) is P p The initial value of (k) represents the covariance matrix of the state vector estimation error; I p×p It is a p×p identity matrix;

[0020] Based on the state vector and its error covariance matrix P p (k–1), calculate the predicted state vector and its error covariance matrix at time k using the following formula:

[0021]

[0022] P p (k|k-1)=P p (k-1) (6)

[0023] In the formula, (k|k–1) represents the prediction of the value at time k based on the measurement at time k–1;

[0024] Calculate the Kalman gain:

[0025] Kalman gain K p The formula for calculating (k) is as follows:

[0026]

[0027] In the formula, R(k) is the covariance of the measurement noise v(k), which can be defined as a constant 0.001 based on test observations;

[0028] Using updated measurement data y(k) and Kalman gain K p (k) is used to correct the predicted state vector;

[0029]

[0030] The error covariance matrix of the best state estimate is updated as follows:

[0031] P p (k)=P p (k|k-1)-K p (k)H p (k)P p (k|k-1) (9).

[0032] In some embodiments, determining the order of the AR model using MMPF includes the following steps:

[0033] x determined by the minimum mean square error p Best estimate of (k):

[0034]

[0035] in It is the estimated parameter vector of the j-th order AR model at time k using equation (8) by KF; p j (k) is the probability of selecting the j-th order AR model as the correct model;

[0036] The probability of the j-th order AR model at time k is calculated using the following formula:

[0037]

[0038] In the formula,

[0039]

[0040]

[0041]

[0042] Among them, equation (13) P can be calculated using equation (5) and equation (14) j (k) can be calculated using equation (9);

[0043] At each time step, the model with the highest probability is the one that MMPF selects as the correct model.

[0044] In some embodiments, the parameters and order of the AR model are used to calculate eigenvalues ​​through characteristic equations, including the following steps:

[0045] The characteristic equation of the AR model is shown below:

[0046] z p -a1(k)z p-1 -...-a p-1 (k)za p (k)=0 (15)

[0047] The pole z of the above polynomial i The eigenvalues ​​are the same as those of the power system in the discrete time domain;

[0048] Then, the eigenvalues ​​in the continuous domain are calculated using the following formula:

[0049] s i =f s lnz i =α i +jω i (16)

[0050] In the formula, f s The sampling rate is i = 1, ..., n, where n represents the total number of oscillation modes in the system; the eigenvalue s i The real part a i Let s represent the damping of the i-th oscillation mode. i The imaginary part ω i This represents the angular frequency of the i-th oscillation mode.

[0051] In some embodiments, determining the damping and frequency of the SSO mode through eigenvalues ​​includes the following steps:

[0052] The oscillation frequency f of the i-th SSO mode SSOi (k) can be calculated using the following formula:

[0053] f SSOi (k)=imag[s i (k)] / 2π (17)

[0054] The damping a of the i-th SSO mode SSOi (k) can be calculated using the following formula:

[0055] α SSOi(k)=re[s i (k)] (18).

[0056] In some embodiments, the construction of the AR model includes the following steps:

[0057] The p-order AR model AR(p) is shown in equation (1):

[0058]

[0059] In the formula, y(k) is the measurement signal at time k, which can be the active power of a node in the power system; a1(k), a2(k), ..., a p (k) represents the time-varying coefficients of the AR model; v(k) represents white noise.

[0060] Secondly, this application proposes a subsynchronous oscillation mode identification system applying environmental excitation as described in any one of the first aspects, comprising:

[0061] Model building module: Constructs AR models based on power system data under SSO modes;

[0062] Parameter calculation module: Calculates the parameters of AR models at different orders using KF;

[0063] Order calculation module: Determines the order of the AR model using MMPF;

[0064] Eigenvalue calculation module: Calculates eigenvalues ​​from the parameters and order of the AR model using characteristic equations;

[0065] Output module: Determines the damping and frequency under SSO mode through eigenvalues.

[0066] Thirdly, this application proposes a terminal device, including a memory, a processor, and a computer program stored in the memory and capable of running on the processor. The memory stores the computer program capable of running on the processor. When the processor loads and executes the computer program, it employs a subsynchronous oscillation mode identification method based on environmental excitation as described in any of the first aspects.

[0067] Thirdly, this application proposes a computer-readable storage medium storing a computer program, which, when loaded and executed by a processor, employs a subsynchronous oscillation mode identification method based on environmental excitation as described in any of the first aspects.

[0068] The beneficial effects of this invention are:

[0069] This invention combines the KF algorithm with MMPF, enabling the algorithm to directly identify weak SSO modal parameters from environmental data. The method has simple operation steps and short response time. Attached Figure Description

[0070] The invention will now be further described with reference to the accompanying drawings.

[0071] Figure 1 This is a flowchart illustrating the specific method in a particular embodiment of the present invention;

[0072] Figure 2 The waveform diagram of the second-order system response in the example is shown;

[0073] Figure 3 The diagram illustrates the process of selecting the order of the AR model in the example.

[0074] Figure 4 The diagram shows the estimation of the system damping ratio and frequency by the MMPF-KF algorithm in the example. Detailed Implementation

[0075] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0076] In the description of this specification, references to terms such as "an embodiment," "example," "specific example," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the invention. In this specification, illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0077] This specific implementation discloses a method for identifying subsynchronous oscillation modes under environmental excitation based on MMPF-KF. The detailed flowchart of the method is as follows: Figure 1 As shown, the specific steps include:

[0078] 1. Calculate the parameters of AR models of different orders using KF, including the following steps:

[0079] An AR model of order p, AR(p), is shown in equation (1):

[0080]

[0081] In the formula, y(k) is the measurement signal at time k, which can be the active power of a node in the power system; a1(k), a2(k), ..., a p (k) represents the time-varying coefficients of the AR model; v(k) represents white noise.

[0082] In order for KF to estimate the AR model parameters, AR(p) needs to be transformed into a standard state-space equation:

[0083] x p (k)=x p (k-1) (2)

[0084] y(k)=H p (k)x p (k)+v(k) (3)

[0085] Equation (2) is the process equation, x p (k) is the state vector, x p (k)=[a1(k)a2(k)…a p (k)] T Equation (3) is the measurement equation, H p (k) is the observation matrix, H p (k)=[y(k-1)y(k-2)...y(kp)].

[0086] Using KF estimation, parameters a1(k), a2(k), ..., a p The steps for (k) are as follows:

[0087] (1) Initialize parameters

[0088]

[0089] In the formula, P p (1) is P p The initial value of (k) represents the covariance matrix of the state vector estimation error; I p×p It is a p×p identity matrix; n is a large integer, which can be 100.

[0090] (2) Prediction

[0091] Based on the previously estimated state vector and its error covariance matrix P p (k–1), calculate the predicted state vector and its error covariance matrix at time k using the following formula:

[0092]

[0093]

[0094] In the formula, (k|k–1) represents the prediction of the value of time k based on the measurement at time k–1.

[0095] (3) Calculate Kalman gain

[0096] Kalman gain K p The formula for calculating (k) is as follows:

[0097]

[0098] In the formula, R(k) is the covariance of the measurement noise v(k), which can be defined as a constant 0.001 based on test observations.

[0099] (4)Amendment

[0100] Using the updated measurement data y(k) and Kalman gain K p (k) is used to correct the predicted state vector.

[0101]

[0102] The error covariance matrix of the best state estimate is updated as follows:

[0103] P p (k)=P p (k|k-1)-K p (k)H p (k)P p (k|k-1) (9)

[0104] Then, repeat from step (2) until all values ​​are calculated.

[0105] 2. Accurately determine the order of the AR model using MMPF, including the following steps:

[0106] Assume the order p of the AR model ranges from 1 to M. Then, a set of KF calculations is performed on AR models of different orders. x is determined by the minimum mean square error (MMSE). p The best estimate of (k) can be obtained in the following way:

[0107]

[0108] in p is the estimated parameter vector of the j-th order AR model at time k using equation (8) by KF. j (k) represents the probability of selecting the j-th order AR model as the correct model. Typically, there is no prior knowledge of the probability, and the initial value is p. j (1) It can be set to 1 / M.

[0109] The probability of the j-th order AR model at time k is calculated using the following formula:

[0110]

[0111] In the formula,

[0112]

[0113]

[0114]

[0115] Among them, L j (k), ε j (k) and Q j (k) refers to formulas, as shown in equations (12), (13), and (14). Equation (13) P can be calculated using equation (5) and equation (14) j (k) can be calculated using equation (9).

[0116] At each time step, the model with the highest probability is the one selected as the correct model by MMPF. The highest probability is typically close to 1, while other probabilities are typically close to 0. All the key functions (KFs) required in MMPF can be run in parallel, which is a key feature of MMPF; therefore, the computation time for identifying SSO modes does not increase.

[0117] 3. After determining the AR model parameters and order, the method for calculating eigenvalues ​​using the characteristic equation to determine the damping and frequency of the SSO mode is as follows:

[0118] The characteristic equation of the AR model is shown below:

[0119] z p -a1(k)z p-1 -...-a p-1 (k)za p (k)=0 (15)

[0120] The pole z of the above polynomial i The eigenvalues ​​are the same as those of the power system in the discrete-time domain. Then, the eigenvalues ​​in the continuous domain are calculated using the following formula:

[0121] s i =f s lnz i =α i +jω i (16)

[0122] In the formula, f s The sampling rate is i = 1, ..., n, where n represents the total number of oscillation modes in the system; the eigenvalue s i The real part a iLet s represent the damping of the i-th oscillation mode. i The imaginary part ω i This represents the angular frequency of the i-th oscillation mode.

[0123] The oscillation frequency f of the i-th SSO mode SSOi (k) can be calculated using the following formula:

[0124] f SSOi (k)=imag[s i (k)] / 2π (17)

[0125] The damping a of the i-th SSO mode SSOi (k) can be calculated using the following formula:

[0126] α SSOi (k)=re[s i (k)] (18)

[0127] Figures 2 to 4 A computational example is presented to demonstrate the effectiveness and accuracy of the proposed process. This example is tested using a classical second-order system:

[0128]

[0129] In the formula, F(t) is the random external excitation signal. y(t) is the random response to the excitation. ξ and f are the damping ratio and natural frequency of the second-order system, respectively. The damping ratio can be calculated by the following formula:

[0130]

[0131] In the formula, α is the damping of the second-order system.

[0132] In this example, the driving signal F(t) is Gaussian white noise with a variance of 1; the damping ratio ξ is set to 2% and the natural frequency f is set to 25Hz; the maximum model order M of MMPF is set to 4; and the sampling rate is 200Hz. Figure 2 The waveform diagram of the second-order system response is shown. Figure 3 The diagram illustrates the process of selecting the order of the AR model. Figure 4 The graphs showing the estimation of the system damping ratio and frequency by the MMPF-KF algorithm are presented.

[0133] like Figure 2 As shown, the response excited by Gaussian white noise also exhibits noise-like properties. Figure 3 As shown, MMPF quickly selects 2 as the correct model order. Figure 4 As shown in (a), the estimated average damping ratio is 2.02%, which is very close to the true damping ratio of the second-order system, and the estimation process converges within 3 seconds. Figure 4As shown in (b), the average estimated frequency is 24.92 Hz, which is very close to the true natural frequency of the second-order system. Furthermore, the frequency estimation process converges within 1 second, which is much faster than the damping ratio estimation. These results demonstrate that the MMPF-KF method can accurately and rapidly identify modal parameters from the environmental response of a second-order system, validating the effectiveness of the proposed method.

[0134] This application discloses a subsynchronous oscillation mode identification system based on environmental excitation, comprising:

[0135] Model building module: Constructs AR models based on power system data under SSO modes;

[0136] Parameter calculation module: Calculates the parameters of AR models at different orders using KF;

[0137] Order calculation module: Determines the order of the AR model using MMPF;

[0138] Eigenvalue calculation module: Calculates eigenvalues ​​from the parameters and order of the AR model using characteristic equations;

[0139] Output module: Determines the damping and frequency under SSO mode through eigenvalues.

[0140] This application also discloses a terminal device, including a memory, a processor, and a computer program stored in the memory and capable of running on the processor. When the processor executes the computer program, it employs any of the environmental excitation subsynchronous oscillation mode identification methods described in the above embodiments.

[0141] The terminal device can be a computer device such as a desktop computer, a laptop computer, or a cloud server. The terminal device includes, but is not limited to, a processor and a memory. For example, the terminal device may also include input / output devices, network access devices, and buses.

[0142] The processor can be a central processing unit (CPU). Of course, depending on the actual use, it can also be other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), off-the-shelf programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The general-purpose processor can be a microprocessor or any conventional processor, etc., and this application does not limit it in this regard.

[0143] The memory can be an internal storage unit of the terminal device, such as a hard disk or RAM of the terminal device, or an external storage device of the terminal device, such as a plug-in hard disk, smart memory card (SMC), secure digital card (SD), or flash memory card (FC) equipped on the terminal device. Furthermore, the memory can be a combination of internal storage units and external storage devices of the terminal device. The memory is used to store computer programs and other programs and data required by the terminal device. The memory can also be used to temporarily store data that has been output or will be output. This application does not limit this.

[0144] In this terminal device, any of the environmental excitation subsynchronous oscillation mode identification methods in the above embodiments can be stored in the memory of the terminal device and loaded and executed on the processor of the terminal device for convenient use.

[0145] This application also discloses a computer-readable storage medium, which stores a computer program, wherein when the computer program is executed by a processor, it employs any of the environmental excitation subsynchronous oscillation mode identification methods described in the above embodiments.

[0146] The computer program can be stored in a computer-readable medium. The computer program includes computer program code, which can be in the form of source code, object code, executable file, or certain middleware. The computer-readable medium includes any entity or device capable of carrying computer program code, recording media, USB flash drive, portable hard drive, magnetic disk, optical disk, computer memory, read-only memory (ROM), random access memory (RAM), electrical carrier signals, telecommunication signals, and software distribution media, etc. It should be noted that the computer-readable medium includes, but is not limited to, the above-mentioned components.

[0147] The subsynchronous oscillation mode identification method of any of the above embodiments can be stored in the computer-readable storage medium and loaded and executed on the processor to facilitate the storage and application of the above method.

[0148] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the claimed invention.

Claims

1. A method for identifying subsynchronous oscillation modes under environmental excitation, characterized in that, Includes the following steps: Constructing an AR model based on power system data under SSO mode; Parameters of AR models at different orders are calculated using KF. Determine the order of the AR model using MMPF; The parameters and order of the AR model are used to calculate the eigenvalues ​​through the characteristic equation; Damping and frequency in the SSO mode are determined by eigenvalues.

2. The method for identifying subsynchronous oscillation modes under environmental excitation according to claim 1, characterized in that, The calculation of parameters of AR models of different orders using KF includes the following steps: Convert the AR model into a state-space equation: (2) (3) in, y ( k ) for at time k The measurement signal; x p ( k ) is the state vector, x p ( k )=[ a 1( k ) a 2( k … a p ( k )] T ; a 1( k ), a 2( k ), ..., a p ( k H represents the time-varying coefficients of the AR model; p ( k H is the observation matrix. p ( k )=[ y ( k -1) y ( k -2) … y ( k - p )]; v ( k () represents white noise; Initialization parameters: (4) In the formula, P p (1) is P p ( k The initial value of ) represents the covariance matrix of the state vector estimation error; I p×p for p × p The identity matrix; Based on the state vector and its error covariance matrix P p ( k –1), calculate at time t using the following formula. k The predicted state vector and its error covariance matrix: (5) (6) In the formula, ( k | k –1) indicates time-based k The measurement at -1 corresponds to time... k Predict the value; Calculate the Kalman gain: Kalman gain K p ( k The calculation formula for ) is as follows: (7) In the formula, R ( k () is a noise measurement v ( k The covariance of () is based on test observations. R ( k The value is set to a constant of 0.001; Utilizing updated measurement data y ( k ) and Kalman gain K p ( k To correct the predicted state vector; (8) The error covariance matrix of the best state estimate is updated as follows: (9)。 3. The method for identifying subsynchronous oscillation modes under environmental excitation according to claim 2, characterized in that, The determination of the order of the AR model using MMPF includes the following steps: x determined by the minimum mean square error p ( k Best estimate: (10) in KF uses equation (8) at time... k place j The estimated parameter vector of the AR model; p j ( k ) is a choice j The probability that the AR model is the correct model; At any moment k Place j The probability of the AR model is calculated using the following formula: (11) In the formula, (12) (13) (14) Among them, equation (13) Using equation (5), P of equation (14) is calculated. j ( k ) Calculate using equation (9); At each time step, the model with the highest probability is the one that MMPF selects as the correct model.

4. The method for identifying subsynchronous oscillation modes under environmental excitation according to claim 3, characterized in that, The process of calculating eigenvalues ​​from the parameters and order of the AR model using characteristic equations includes the following steps: The characteristic equation of the AR model is shown below: (15) The poles of the above polynomial z i The eigenvalues ​​are the same as those of the power system in the discrete time domain; Then, the eigenvalues ​​in the continuous domain are calculated using the following formula: (16) In the formula, f s It is the sampling rate; i =1, ..., n , n Represents the total number of oscillation modes in the system; eigenvalues s i real part a i Indicates the first i Damping of each oscillation mode, s i imaginary part ω i Indicates the first i The angular frequency of each oscillation mode.

5. The method for identifying subsynchronous oscillation modes under environmental excitation according to claim 4, characterized in that, The method of determining the damping and frequency under the SSO mode through eigenvalues ​​includes the following steps: No. i Oscillation frequency of each SSO mode f SSOi ( k Calculated using the following formula: (17) And the first i Damping of each SSO mode a SSOi ( k Calculated using the following formula: (18)。 6. The method for identifying subsynchronous oscillation modes under environmental excitation according to claim 1, characterized in that, The construction of the AR model includes the following steps: p AR model of order AR ( p As shown in equation (1): (1) In the formula, y ( k ) for at time k The measured signal is the active power of the power system node; a 1( k ), a 2( k ), ..., a p ( k ) represents the time-varying coefficients of the AR model; v ( k () represents white noise.

7. A subsynchronous oscillation mode identification system for environmental excitation, characterized in that, The system is used to implement the subsynchronous oscillation mode identification method for environmental excitation as described in any one of claims 1 to 6, and the system comprises: Model building module: Constructs AR models based on power system data under SSO modes; Parameter calculation module: Calculates the parameters of AR models at different orders using KF; Order calculation module: Determines the order of the AR model using MMPF; Eigenvalue calculation module: Calculates eigenvalues ​​from the parameters and order of the AR model using characteristic equations; Output module: Determines the damping and frequency under SSO mode through eigenvalues.

8. A terminal device, comprising a memory, a processor, and a computer program stored in the memory and capable of running on the processor, characterized in that, The memory stores a computer program that can run on the processor. When the processor loads and executes the computer program, it employs a subsynchronous oscillation mode identification method based on environmental excitation as described in any one of claims 1 to 6.

9. A computer-readable storage medium storing a computer program, characterized in that, When the computer program is loaded and executed by the processor, it employs a subsynchronous oscillation mode identification method based on environmental excitation as described in any one of claims 1 to 6.

Images