A power system inertia identification method and device, a terminal and a storage medium

Abstract

The application discloses a power system inertia identification method and device, a terminal and a storage medium, and relates to the technical field of power system inertia identification. According to the scheme, the three-phase electrical quantity data at the generator node and the grid connection point are subjected to sliding window processing, and the preliminary frequency estimation value is calculated in combination with the frequency domain characteristic state quantity. The preliminary frequency estimation values of adjacent time windows before and after the target time window are used to construct a frequency change rate characteristic vector. Then, the power disturbance characteristic quantity and the frequency change rate characteristic vector are uniformly included in a target function containing an inertia coupling term, and joint optimization solving is performed on the dynamic phasor parameter and the frequency estimation result, thereby realizing bidirectional closed-loop coupling of frequency identification and inertia estimation, effectively improving the overall consistency of the frequency identification result and the inertia calculation result, and enabling the method to provide more stable and reliable frequency change information and inertia characteristic result, thereby effectively supporting stable control, fault analysis and operation mode adjustment of the system.

Description

Technical Field

[0001] This application relates to the field of power system inertia identification technology, and in particular to a power system inertia identification method, device, terminal and storage medium. Background Technology

[0002] With the high proportion of renewable energy connected to the grid, the rapid increase in power electronic interface equipment, and the continuous decline in the proportion of traditional synchronous power sources, the power system's operating status exhibits stronger sensitivity to frequency fluctuations, making the online sensing of the system's equivalent inertia increasingly urgent. Especially in scenarios involving large disturbances, black starts, islanded operation, and high-proportion renewable energy transmission, frequency changes occur more rapidly and adjustment times are shorter. Operators and control systems need timely and reliable frequency change information and inertia characterization results to support stability control, fault analysis, and operational mode adjustments. Therefore, developing more stable inertia identification technologies that better conform to physical processes during dynamic disturbances has become a practical requirement in the field of power system monitoring and control.

[0003] Traditional inertia identification typically employs a sequential process: first, frequency identification; then, the rate of frequency change; and finally, inertia estimation in conjunction with disturbance quantities. While this method is relatively intuitive and effective under steady-state or slowly varying conditions, it often has significant limitations in actual power grid disturbance processes. For example, the independent stages of frequency identification and inertia calculation prevent inertia information from influencing the front-end frequency identification process, resulting in a lack of physical consistency and ultimately unstable and unreliable identification results. Summary of the Invention

[0004] This application provides a method, device, terminal, and storage medium for power system inertia identification, which aims to improve the accuracy and reliability of power system inertia identification under dynamic disturbance conditions.

[0005] To achieve the above-mentioned objectives, the first aspect of this application provides a method for power system inertia identification, comprising:

[0006] Collect three-phase electrical quantity data at the generator node grid connection point;

[0007] For each phase of the three-phase electrical quantity data, the frequency domain response at each joint observation frequency point is calculated by using a preset sliding time window and combining several joint observation frequency points set at equal intervals within the current window. Based on the frequency domain response, the frequency domain characteristic state quantity of the current window is determined. Then, based on the frequency domain response and the frequency domain characteristic state quantity, the preliminary frequency estimate of the current window is calculated. By moving the sliding time window and performing cyclic calculations, the preliminary frequency estimate of each window is obtained sequentially.

[0008] Based on the center sampling time of each window, and combined with the three-phase electrical quantity data corresponding to the center sampling time, the power disturbance characterization quantity at the center sampling time is calculated;

[0009] For a selected target time window, two adjacent time windows are determined, and a frequency change rate representation vector is constructed based on the preliminary frequency estimates of the adjacent time windows.

[0010] Based on the power disturbance characterization quantity and the frequency change rate characterization vector of the target time window, and combined with the preset extended frequency domain model, a unified objective function containing inertia coupling terms is constructed. By solving the unified objective function, the corresponding dynamic phasor parameter vector and frequency estimation result for each corresponding dynamic phasor parameter within the target time window are obtained. The extended frequency domain model is configured with the correspondence between the dynamic phasor parameters to be identified and the Taylor coefficients.

[0011] Preferably, the frequency domain characteristic state quantities include: frequency domain bending state quantities and frequency domain asymmetric state quantities;

[0012] The formula for calculating the frequency domain bending state quantity is as follows:

[0013] ;

[0014] The formula for calculating the frequency domain asymmetric state quantity is:

[0015] ;

[0016] In the formula, This refers to the frequency domain bending state quantity. The frequency domain asymmetric state quantity, This indicates the operation of taking the real part. To prevent positive numbers with a denominator of zero, As the reference frequency, This is the frequency offset. This represents the frequency domain response.

[0017] Preferably, before calculating the frequency domain response at each joint observation frequency point by combining several joint observation frequency points set at equal intervals within the current window, the method further includes:

[0018] For each phase of the three-phase electrical quantity data, a sampling sequence within the current window is extracted through a preset sliding time window, and a normalized position variable is determined based on the sampling sequence.

[0019] Obtain the prior estimate value corresponding to the current window. When the current window is not the first window, the prior estimate value is the offset estimate value and frequency domain characteristic state quantity of the previous window. The offset estimate value is the relative deviation between the preliminary frequency estimate value and the nominal frequency of the corresponding window. When the current window is the first window, the prior estimate value is the obtained preset value.

[0020] Substituting the normalized location variable and the prior estimate into a preset bending adaptive window weighting function, the window weighting coefficients corresponding to each time window are obtained.

[0021] The sampling sequence is weighted point by point according to the window weighting coefficients to obtain a weighted sampling sequence.

[0022] Preferably, the bending adaptive window weighting function is specifically:

[0023] ;

[0024] In the formula, n and q both represent the index of the sampling point, and N is the total number of sampling points within the window. For the normalized position variable corresponding to the nth sampling point, Let q be the normalized position variable corresponding to the sampling point. , and This is the window adjustment coefficient. The offset prediction value in the prior estimate, The frequency domain bending state quantity in the prior estimate. The frequency domain asymmetric state quantity in the prior estimate. This is an adjustable window adaptive coefficient.

[0025] Preferably, calculating the preliminary frequency estimate for the current window based on the frequency domain response and the frequency domain characteristic state variables includes:

[0026] Based on the frequency domain response and the frequency domain characteristic state quantity, the frequency weight of each joint observation frequency point is calculated using a preset frequency weight formula, and the weighted frequency domain response is calculated based on the frequency weight and the frequency domain response.

[0027] The real and imaginary parts of the weighted frequency domain response corresponding to each joint observation frequency point are stacked in a preset order to form a real joint observation vector;

[0028] For the sampling sequence within the sliding time window, a Taylor series expansion is performed at the center of the window to construct a dynamic phasor model with local time variables as independent variables. Then, based on the centered position variable of the sampling sequence and the bending adaptive window weighting function, the frequency domain weighting transformation coefficients of each order are obtained through a preset frequency domain weighting transformation coefficient calculation formula. The frequency domain weighting transformation coefficients are used to characterize the contribution of the temporal changes of different Taylor coefficient orders at each joint observation frequency point.

[0029] Frequency point offset compensation is performed on the frequency domain weighted transform coefficients based on the frequency domain characteristic state variables, and then an extended frequency domain model is constructed based on the compensated frequency domain weighted transform coefficients.

[0030] The preliminary dynamic phasor parameter vector and the corresponding preliminary dynamic phasor Taylor coefficients of the current window are obtained by solving the extended frequency domain model using the least squares method. Based on the preliminary dynamic phasor parameter vector and the preliminary dynamic phasor Taylor coefficients, the preliminary frequency estimate of the current window is calculated.

[0031] Preferably, constructing the extended frequency domain model based on the compensated frequency domain weighted transform coefficients includes:

[0032] Based on the compensated frequency domain weighted transform coefficients, the main observation row vector, the mirror observation row vector, and the real observation submatrix corresponding to a single frequency point are constructed, and the real observation submatrixes of each frequency point are combined to form an extended observation matrix;

[0033] Obtain the correspondence between the dynamic phasor parameters to be identified and the Taylor coefficients of the dynamic phasor;

[0034] The real-valued joint observation vector, the extended observation matrix, and the dynamic phasor parameters to be identified are used to construct the corresponding observation equation for joint optimization, thereby obtaining the extended frequency domain model.

[0035] Preferably, after obtaining the corresponding dynamic phasor parameter vector and frequency estimation results within the target time window, the method further includes:

[0036] For the dynamic phasor parameter vector of each phase within any time window, the dynamic phasor parameter vector is restored to the Taylor coefficients of the corresponding dynamic phasors of each order.

[0037] Based on the Taylor coefficients of each order of dynamic phasors, a positive order transformation is performed to obtain a unified positive order set of Taylor coefficients of dynamic phasors.

[0038] Based on the unified positive-sequence dynamic phasor Taylor coefficient set, the zero-order positive-sequence dynamic phasor Taylor coefficient and the first-order positive-sequence dynamic phasor Taylor coefficient are extracted, and combined with the preset unified main frequency calculation formula, the unified main frequency is calculated. Then, through iterative processing, the unified main frequency corresponding to each time window is calculated in sequence to obtain the unified main frequency sequence.

[0039] The unified main frequency sequence is used as the frequency sequence to be cleaned. The median reference frequency is extracted from the frequency sequence to be cleaned. Based on the median reference frequency and the peak decision threshold, the peak points in the frequency sequence to be cleaned are corrected to obtain the frequency estimate after peak filtering.

[0040] A band-constrained low-pass filter is performed on the frequency estimate after peak filtering and the corresponding power perturbation sequence to obtain the filtered power perturbation sequence and the frequency change rate sequence corresponding to different window paths, wherein the window paths include: Hann window path and rectangular window path;

[0041] Based on the filtered power disturbance sequence and the frequency change rate sequence, the inertia identification result used to characterize the power system inertia level is calculated using a preset power system inertia identification result calculation formula. The inertia identification result includes: power disturbance amplitude, extreme values ​​of frequency change rate corresponding to different window paths, and power system equivalent inertia identification result.

[0042] A second aspect of this application provides a power system inertia identification device, comprising:

[0043] The electrical quantity acquisition unit is used to collect three-phase electrical quantity data at the generator node grid connection point;

[0044] The preliminary frequency estimation unit is used to calculate the frequency domain response at each joint observation frequency point for each phase of the three-phase electrical quantity data through a preset sliding time window and in combination with several joint observation frequency points set at equal intervals within the current window. Based on the frequency domain response, the unit determines the frequency domain characteristic state quantity of the current window. Based on the frequency domain response and the frequency domain characteristic state quantity, the unit calculates the preliminary frequency estimate of the current window. By moving the sliding time window and performing cyclic calculations, the unit obtains the preliminary frequency estimate of each window in turn.

[0045] The power disturbance characterization calculation unit is used to calculate the power disturbance characterization quantity at the center sampling time based on the center sampling time of each window and the three-phase electrical quantity data corresponding to the center sampling time.

[0046] The frequency change rate representation vector calculation unit is used to determine two adjacent time windows of the selected target time window, and construct a frequency change rate representation vector based on the preliminary frequency estimates of the adjacent time windows.

[0047] The power system inertia identification unit is used to construct a unified objective function containing inertia coupling terms based on the power disturbance characterization quantity and the frequency change rate characterization vector of the target time window, combined with a preset extended frequency domain model. By solving the unified objective function, the unit obtains the dynamic phasor parameter vector and frequency estimation result for each corresponding dynamic phasor parameter within the target time window. The extended frequency domain model is configured with the correspondence between the dynamic phasor parameters to be identified and Taylor coefficients.

[0048] A third aspect of this application provides a power system inertia identification terminal, comprising: a memory and a processor;

[0049] The memory is used to store program code, which corresponds to the power system inertia identification method provided in the first aspect of this application.

[0050] The processor is used to read and execute the program code to implement the power system inertia identification method.

[0051] The fourth aspect of this application provides a computer-readable storage medium storing program code, which is read and executed by a processor to implement the power system inertia identification method provided in the first aspect of this application.

[0052] As can be seen from the above technical solutions, this application has the following advantages:

[0053] The proposed solution first processes the three-phase electrical quantity data at the generator node grid connection point using a sliding window method and calculates a preliminary frequency estimate by combining it with frequency domain characteristic state quantities. Then, it constructs a frequency change rate representation vector using the preliminary frequency estimates of adjacent time windows before and after the target time window. Next, it integrates the power disturbance representation quantity and the frequency change rate representation vector into an objective function containing an inertia coupling term, and performs joint optimization on the dynamic phasor parameters and frequency estimation results. This achieves a two-way closed-loop coupling between frequency identification and inertia estimation, effectively improving the overall consistency between the frequency identification results and the inertia calculation results. This allows the method to provide more stable and reliable frequency change information and inertia representation results when facing dynamic scenarios such as large disturbances in the power system, black start, islanded operation, and high proportion of renewable energy transmission, effectively supporting the system's stable control, fault analysis, and operation mode adjustment. Attached Figure Description

[0054] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0055] Figure 1 This is a flowchart illustrating an embodiment of a power system inertia identification method provided in this application.

[0056] Figure 2 This is a schematic diagram of the architecture of an embodiment of a power system inertia identification device provided in this application.

[0057] Figure 3 This is a schematic diagram of the architecture of an embodiment of a power system inertia identification terminal provided in this application. Detailed Implementation

[0058] Traditional existing power system inertia identification methods suffer from unstable identification results and low reliability during actual power grid disturbances because frequency identification and inertia calculation are carried out independently in stages. Inertia information cannot constrain the front-end frequency identification process in reverse, resulting in a lack of physical consistency.

[0059] In view of this, embodiments of this application provide a power system inertia identification method, device, terminal, and storage medium to achieve the invention objective of improving the accuracy and reliability of power system inertia identification under dynamic disturbance conditions.

[0060] To make the inventive objectives, features, and advantages of this application more apparent and understandable, the technical solutions in the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the embodiments described below are only some embodiments of this application, and not all embodiments. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0061] First, a detailed description of an embodiment of a power system inertia identification method provided in this application is as follows:

[0062] Please see Figure 1 This application provides an embodiment of a power system inertia identification method, the steps of which include:

[0063] Step 101: Collect three-phase electrical quantity data at the generator node grid connection point;

[0064] Step 102: For each phase of the three-phase electrical quantity data, calculate the frequency domain response at each joint observation frequency point by using a preset sliding time window and combining several joint observation frequency points set at equal intervals within the current window. Based on the frequency domain response, determine the frequency domain characteristic state quantity of the current window. Then, based on the frequency domain response and the frequency domain characteristic state quantity, calculate the preliminary frequency estimate of the current window. By moving the sliding time window and performing cyclic calculations, obtain the preliminary frequency estimate of each window in sequence.

[0065] Step 103: Based on the center sampling time of each window, and combined with the three-phase electrical quantity data corresponding to the center sampling time, calculate the power disturbance characterization quantity at the center sampling time;

[0066] Step 104: For the selected target time window, determine two adjacent time windows of the target time window, and construct a frequency change rate representation vector based on the preliminary frequency estimates of the adjacent time windows.

[0067] Step 105: Based on the power disturbance characterization quantity and the frequency change rate characterization vector of the target time window, and combined with the preset extended frequency domain model, construct a unified objective function containing inertia coupling terms, and obtain the dynamic phasor parameter vector and frequency estimation result for each corresponding term within the target time window by solving the unified objective function.

[0068] For ease of understanding, the following explains some key terms in this embodiment:

[0069] Three-phase electrical quantity data refers to the measured values ​​of common three-phase AC electrical signals such as voltage and current in a power system, typically including instantaneous values ​​or sampled values. These data are fundamental information about the operating status of the power system.

[0070] Sliding time window: refers to a time interval in a continuous sampled data stream where a segment of data of a fixed length is extracted for processing. This time interval slides forward as time progresses, thereby enabling continuous analysis of dynamic processes.

[0071] Joint observation frequency points: In frequency domain analysis, this refers to a set of discrete frequency points pre-defined to capture the frequency characteristics of a signal. By calculating the frequency domain response of the signal at these frequency points, the energy distribution or phase information of the signal at different frequencies can be obtained.

[0072] Frequency domain response refers to the amplitude and phase information of a signal at a specific frequency point, usually obtained through frequency domain analysis methods such as Fourier transform. It reflects the characteristics of the signal in the frequency domain.

[0073] Frequency domain characteristic state quantities: These are quantities extracted from the frequency domain response of a signal that characterize the dynamic frequency properties of the signal. These quantities can be used to describe complex changes in signal frequency, such as curvature and asymmetry.

[0074] Preliminary frequency estimate: This refers to an initial or rough estimate of the signal frequency based on local frequency domain information within a sliding time window. This estimate serves as the basis for subsequent, more accurate frequency identification.

[0075] Power disturbance characterization: refers to an indicator used to quantify changes in system power when disturbances occur in a power system. It can be a change in instantaneous power, active power, or reactive power, reflecting the input-output imbalance of the system's energy.

[0076] Frequency change rate representation vector: This vector, constructed from preliminary frequency estimates across multiple adjacent time windows, describes the rate and trend of frequency change over time. It provides quantitative information about the dynamic changes in frequency.

[0077] Extended frequency domain model: This refers to a mathematical model established in the frequency domain to describe the frequency response and inertia characteristics of a power system under dynamic disturbances. This model typically includes dynamic phasor parameters to be identified and Taylor coefficients, and can couple inertia information with the frequency dynamics.

[0078] Inertia coupling term: In a mathematical model, this term characterizes the influence of power system inertia on frequency dynamics. By introducing this term, the system's inertia characteristics can be directly correlated with the rate of frequency change, enabling joint identification of inertia and frequency.

[0079] Unified objective function: This refers to a comprehensive optimization function that contains multiple variables to be optimized and constraints. By solving this function, multiple related identification results can be obtained simultaneously, such as dynamic phasor parameters and frequency estimation results.

[0080] Dynamic phasor parameter vector: In a dynamic phasor model, this refers to the set of parameters used to describe the dynamic characteristics of a signal, such as amplitude, phase, and frequency. These parameters can accurately characterize the transient behavior of a signal in both the time and frequency domains.

[0081] Frequency estimation result: This refers to the precise estimate of the power system frequency obtained through this method. The result reflects the actual frequency variation of the system under dynamic disturbances.

[0082] The correspondence between the dynamic phasor parameters to be identified and the Taylor coefficients refers to the mathematical mapping relationship between the dynamic phasor parameters and the Taylor coefficients of each order obtained after expanding the signal into Taylor series at the center of the time window in the extended frequency domain model. This relationship allows the dynamic phasor parameters to be characterized and identified through the Taylor coefficients.

[0083] More specifically, this embodiment provides a method for power system inertia identification, the specific implementation process of which may include the following steps:

[0084] First, three-phase electrical quantity data are collected at the generator node's grid connection point. This data can be real-time sampled values ​​of electrical signals such as voltage and current. For example, by using sensors and data acquisition units installed at the generator's grid connection point, three-phase voltage and current can be continuously sampled, and the sampled data can be stored at a fixed sampling frequency. As one implementation method, analog signals can be digitized by a simple analog-to-digital converter, and the original digital sampled values ​​can be directly recorded.

[0085] Secondly, for each phase of the three-phase electrical quantity data, the frequency domain response at each joint observation frequency point is calculated using a preset sliding time window and several equally spaced joint observation frequency points within the current window. Based on this frequency domain response, the frequency domain characteristic state quantity of the current window is determined. Then, based on this frequency domain response and the frequency domain characteristic state quantity, the preliminary frequency estimate of the current window is calculated. By moving the sliding time window and iterating, the preliminary frequency estimate of each window is obtained sequentially. Specifically, a fixed-length time window can be set, for example, containing N sampling points. This window slides on the data stream with a preset step size. For the sampled data within each window, its frequency domain response can be calculated at the preset joint observation frequency points using methods such as Discrete Fourier Transform or Fast Fourier Transform. For example, the joint observation frequency points can be selected as several frequency points symmetrically distributed to the left and right of the nominal frequency. After obtaining the frequency domain response, based on the amplitude and phase information of these responses, the state quantity characterizing the frequency characteristics of the current window can be extracted using simple weighted averaging or peak detection. For example, the centroid frequency or the frequency corresponding to the maximum amplitude of the frequency domain response can be calculated as a preliminary frequency domain feature. Subsequently, by combining the frequency domain response and these preliminary frequency domain feature state variables, a preliminary frequency estimate can be calculated, for example, based on the frequency corresponding to the maximum amplitude of the frequency domain response. Thus, by continuously moving this sliding time window and repeating the above calculation process, a series of continuous preliminary frequency estimates can be obtained.

[0086] Next, based on the center sampling time of each window and the corresponding three-phase electrical quantity data, the power disturbance characteristic quantity for that center sampling time is calculated. Each sliding time window has a corresponding center sampling time. At that center sampling time, the instantaneous active or reactive power can be calculated using the collected three-phase voltage and current data. For example, the instantaneous voltage and current values ​​at that center sampling time can be directly multiplied to obtain the instantaneous power, which can then be used as the power disturbance characteristic quantity. Alternatively, the average power at that time can be calculated and compared with the system's nominal power to obtain the power disturbance value.

[0087] Then, for the selected target time window, two adjacent time windows are determined. Based on the preliminary frequency estimates of these adjacent time windows, a frequency change rate representation vector is constructed. From a continuous sequence of preliminary frequency estimates, a target time window for which inertia identification is currently required is selected. This target time window typically has one preceding adjacent window and one following adjacent window. The frequency change rate can be simply obtained by calculating the difference between the preliminary frequency estimates of the preceding and following adjacent windows of the target time window and dividing by the corresponding time interval. For example, a vector containing the frequency estimates of these two adjacent windows and their time interval can be constructed.

[0088] Finally, based on the power disturbance characterization and frequency change rate characterization vector within the target time window, and combined with a pre-defined extended frequency domain model, a unified objective function containing an inertia coupling term is constructed. Solving this unified objective function yields the corresponding dynamic phasor parameter vector and frequency estimation results for each element within the target time window. The extended frequency domain model is configured with the correspondence between the dynamic phasor parameters to be identified and the Taylor coefficients. In this step, the previously obtained power disturbance characterization and frequency change rate characterization vector can be used as input. The pre-defined extended frequency domain model can be set as a simplified linear model to describe the relationship between frequency and power disturbance, containing the parameters to be identified. A simple inertia term can be introduced into this model, for example, by directly associating the frequency change rate with the power disturbance. Based on this model, an objective function can be constructed, such as a least-squares objective function, containing the dynamic phasor parameters to be identified and the frequency estimation values, and introducing a coupling term reflecting the influence of inertia. By optimizing the objective function, for example using gradient descent or Newton's method, the corresponding dynamic phasor parameter vectors and the final frequency estimation results within the objective time window can be obtained simultaneously. This extended frequency domain model predefines the mathematical relationship between the dynamic phasor parameters and Taylor coefficients, thus allowing these parameters to be considered simultaneously during the solution process.

[0089] This embodiment integrates power disturbance characterization and frequency change rate characterization into an objective function containing an inertia coupling term, and performs joint optimization on dynamic phasor parameters and frequency estimation results. This achieves bidirectional closed-loop coupling between frequency identification and inertia estimation, effectively improving the overall consistency between frequency identification results and inertia calculation results. This enables the method to provide more stable and reliable frequency change information and inertia characterization results when facing dynamic scenarios such as large disturbances in the power system, black start, islanded operation, and high proportion of renewable energy transmission, effectively supporting the system's stable control, fault analysis, and operation mode adjustment.

[0090] In some of the above-described embodiments, the power system inertia identification method determines the frequency domain characteristic state quantities by calculating the frequency domain response, and calculates the preliminary frequency estimate of the current window based on the frequency domain response and the frequency domain characteristic state quantities. On this basis, this application further proposes that the frequency domain characteristic state quantities include frequency domain bending state quantities and frequency domain asymmetric state quantities.

[0091] Specifically, the frequency-domain curvature state variable is used to characterize the degree of curvature of the frequency response in the frequency domain, i.e., the trend of the rate of frequency change. In power system frequency disturbances, the frequency does not always change linearly; it may exhibit acceleration or deceleration trends. This curvature characteristic is crucial for accurately estimating frequency dynamics. It reflects the second-order characteristics of the frequency-domain response as a function of frequency offset near a reference frequency point, and can capture the nonlinear trend of frequency changes.

[0092] Meanwhile, frequency-domain asymmetric state variables are used to characterize the asymmetry of the frequency response in the frequency domain, i.e., the imbalance of frequency disturbances in the positive and negative directions. In some power system disturbance scenarios, the dynamic characteristics of frequency rise and fall may differ, and this asymmetry is also an important factor affecting the accuracy of frequency estimation. It reflects the first-order characteristics of the frequency domain response as a function of frequency offset near the reference frequency point, and can capture the asymmetric trend of frequency changes.

[0093] By introducing frequency-domain bending state variables and frequency-domain asymmetric state variables as components of the frequency-domain characteristic state variables, this application can more precisely capture the nonlinear and asymmetric dynamic characteristics existing in power system frequency disturbances. Frequency-domain bending state variables can quantify the trend of the frequency change rate, reflecting the degree of frequency acceleration or deceleration; frequency-domain asymmetric state variables can characterize the imbalance of frequency disturbances in the positive and negative directions. The combined use of these two state variables makes the frequency-domain characteristic state variables provide a more comprehensive and accurate description of the actual frequency dynamics. Therefore, the preliminary frequency estimates calculated based on these more accurate frequency-domain characteristic state variables will have higher accuracy and robustness, thus providing more reliable basic data for subsequent power disturbance characterization, frequency change rate construction, and final power system inertia identification, significantly improving the accuracy and reliability of the entire inertia identification method.

[0094] In some of the above embodiments, this application proposes to calculate the frequency domain response at each joint observation frequency point by combining several joint observation frequency points set at equal intervals within the current window, and to determine the frequency domain characteristic state quantity and preliminary frequency estimate based on the frequency domain response. Furthermore, this application proposes to include the following steps before calculating the frequency domain response at each joint observation frequency point by combining several joint observation frequency points set at equal intervals within the current window:

[0095] For each phase of the three-phase electrical quantity data, a sampling sequence within the current window is extracted through a preset sliding time window, and a normalized position variable is determined based on the sampling sequence. A priori estimate corresponding to the current window is obtained. When the current window is not the first window, the priori estimate is the offset estimate and frequency domain characteristic state quantity of the previous window. The offset estimate is the relative deviation between the preliminary frequency estimate and the nominal frequency of the corresponding window. When the current window is the first window, the priori estimate is a preset value obtained. The normalized position variable and the priori estimate are substituted into a preset bending adaptive window weighting function to obtain the window weighting coefficients corresponding to each time window. The sampling sequence is weighted point-by-point according to the window weighting coefficients to obtain a weighted sampling sequence.

[0096] Specifically, for each phase of the three-phase electrical quantity data, a sampling sequence within a preset sliding time window is extracted, and a normalized position variable is determined based on the sampling sequence. This step aims to provide basic data and position references for subsequent weighting processing. The sliding time window is used to dynamically extract data to adapt to the local characteristics of the signal. The sampling sequence is the set of data points within the window. The normalized position variable maps the position of the sampling point within the window to a standardized range (e.g., -1 to 1 or 0 to 1), which helps to improve the generality and computational stability of the subsequent weighting function, making it independent of the specific length of the window or the sampling rate. For example, normalization can be performed based on the ratio of the relative position of the sampling point within the window to the total length of the window.

[0097] For example, by constructing a sliding analysis window of length N, the sampling sequence within the current window is extracted from the voltage sampling data of that phase. ,in And generate normalized position variables. This is used to subsequently construct the window weighting function. Furthermore, the centered position variable can also be obtained from the sampling sequence. This can be used for subsequent frequency domain response calculations at joint observation frequency points; specifically, the centralized location variable... and normalized position variables The calculation formula is as follows:

[0098] ;

[0099] ;

[0100] Subsequently, the prior estimate corresponding to the current window is obtained. When the current window is not the first window, the prior estimate is the offset estimate and frequency domain characteristic state quantity of the previous window. The offset estimate is the relative deviation between the preliminary frequency estimate and the nominal frequency of the corresponding window. When the current window is the first window, the prior estimate is the obtained preset value. The prior estimate is the predictive information used to guide the weighted processing of the current window, aiming to improve the adaptability of the processing. For non-first windows, the preliminary frequency estimate calculated in the previous window is used. With nominal frequency The relative deviation (i.e., the offset estimate) ) and frequency domain characteristic state quantities (such as frequency domain bending state quantities) and frequency domain asymmetric state quantities This serves as prior information. This approach reflects the iterative and adaptive nature of the algorithm, utilizing historical information to optimize the current processing.

[0101] More specifically, for non-first windows, the formula for calculating the prior estimate of the current window can be found in the following example:

[0102] ;

[0103] in, The nominal frequency is indicated by the superscript prev, which represents the quantity that was output from the previous window and mapped to the prior estimate of this window in this step.

[0104] For the first window, since there is no information from the previous window, the prior estimate can be set to a preset empirical value or nominal value to initiate the entire processing flow. These prior values ​​can reflect the local frequency characteristics and distortion degree of the signal, providing a basis for subsequent adaptive weighting.

[0105] Based on this, the normalized position variables and the prior estimates are substituted into a preset curved adaptive window weighting function to obtain the window weighting coefficients corresponding to each time window. The curved adaptive window weighting function is one of the core components of this scheme. It dynamically adjusts the weight of each sampling point based on its normalized position within the window and the prior estimates of the signal (such as frequency shift, curvature, and asymmetry). The design goal of this function is to enable the window shape to adaptively "bend" or adjust to better match the local characteristics of the signal. For example, when there is a frequency shift or harmonic distortion, more relevant data points can be given higher weights, while suppressing noise and irrelevant components. The window weighting coefficients are the weight values ​​corresponding to each sampling point calculated by this function; these weight values ​​will be used to generate subsequent weighted sampling sequences.

[0106] Finally, the sampling sequence is weighted point-by-point according to the window weighting coefficients to obtain a weighted sampling sequence. Point-by-point weighting means multiplying each sampling point in the original sampling sequence by its corresponding window weighting coefficient to obtain a new weighted sampling sequence. This weighting process can effectively highlight the data points within the window that contribute significantly to frequency estimation, while weakening or suppressing the influence of noise, transient components, and data points that do not conform to the current frequency characteristics.

[0107] Specifically, the normalized position variables obtained based on S201 And combined with the estimated offset of the voltage fundamental frequency relative to the nominal frequency of the current window obtained from S202 Prior estimates of frequency domain bending state and frequency domain asymmetric state prior estimates Construct a prototype of a predefined window weighted function. and curved adaptive window weighting function This forms a sequence of window weighting coefficients; specifically, the prototype of the preset window weighting function... and curved adaptive window weighting function The calculation formula is as follows:

[0108] ;

[0109] in, , hour It degenerates into a Hann window. hour Degenerates into a rectangular window. The normalized position variable corresponding to sampling point q is defined in the same way as... same, , and This is the window adjustment coefficient;

[0110] The curved adaptive window weighting function is a mathematical function used to weight sampled sequences. Its core function is to dynamically adjust the contribution of each sampling point to the final frequency domain response calculation based on the signal's local characteristics (reflected by prior estimates) and the relative position of the sampling points within the window. This adaptability allows the weighting process to better adapt to dynamic changes in power system frequency and waveform distortion, thereby improving the accuracy and robustness of frequency estimation. By combining positional information and frequency domain characteristic state variables, this function achieves refined weighting of the sampled data within the window, highlighting effective information and suppressing noise.

[0111] The obtained sampling sequence The obtained bending adaptive window weighting function By performing point-by-point weighting, a weighted sampling sequence is obtained. This provides a foundation for frequency domain response calculation; specifically, the weighted sampling sequence The calculation formula is as follows:

[0112] .

[0113] Based on this, the weighted sampling sequence is obtained This can be combined with the above-mentioned centralized location variables Calculate the frequency domain response at each joint observation frequency point. This is used to extract frequency domain morphological features and construct a joint observation vector; specifically, the frequency domain response... The calculation formula is as follows:

[0114] ;

[0115] in, Represents the m-th joint observation frequency point. The imaginary unit, , ;

[0116] ;

[0117] in, Nominal frequency, Frequency interval, Sampling frequency, As the reference frequency, This represents the frequency offset.

[0118] Based on the frequency domain response obtained from S206, extract the frequency domain bending state quantity of the current window. and frequency domain asymmetric state quantities These two state variables are updated window by window as the sliding analysis window advances, and are calculated using the following formula: The formula for calculating the frequency domain bending state variable is:

[0119] ;

[0120] The formula for calculating the frequency domain asymmetric state quantity is:

[0121] ;

[0122] In the formula, This refers to the frequency domain bending state quantity. The frequency domain asymmetric state quantity, This indicates the operation of taking the real part. To prevent positive numbers with a denominator of zero, This represents the frequency domain response.

[0123] The above technical solution introduces a curved adaptive window weighting process based on prior estimates before calculating the frequency domain response. Specifically, by truncating the sampling sequence and determining the normalized position variables, and combining the offset estimate of the previous window and the frequency domain characteristic state quantities as the prior estimates of the current window (or the preset values ​​of the first window), the weighting function can adaptively adjust the window shape. This adaptive weighting process can effectively suppress the influence of noise and non-stationarity on the frequency domain response calculation, especially when the power system frequency shifts or harmonic distortion exists, it can more accurately capture the local frequency characteristics of the signal. By weighting the sampling sequence point by point, the resulting weighted sampling sequence can provide more accurate and robust input for frequency domain response calculation, thereby significantly improving the accuracy of subsequent frequency domain characteristic state quantities and preliminary frequency estimates. This provides more reliable basic data for the final power system inertia identification, thus improving the accuracy and anti-interference capability of the entire identification method.

[0124] In some of the above embodiments, this application proposes to calculate a preliminary frequency estimate using the frequency domain response and frequency domain characteristic state variables. Based on this, this application further proposes detailed steps for calculating the preliminary frequency estimate of the current window based on the frequency domain response and the frequency domain characteristic state variables, specifically including:

[0125] Based on the frequency domain response and the frequency domain characteristic state quantity, the frequency weight of each joint observation frequency point is calculated using a preset frequency weight formula, and the weighted frequency domain response is calculated based on the frequency weight and the frequency domain response.

[0126] The real and imaginary parts of the weighted frequency domain response corresponding to each joint observation frequency point are stacked in a preset order to form a real joint observation vector;

[0127] For the sampling sequence within the sliding time window, a Taylor series expansion is performed at the center of the window to construct a dynamic phasor model with local time variables as independent variables. Then, based on the centered position variable of the sampling sequence and the bending adaptive window weighting function, the frequency domain weighting transformation coefficients of each order are obtained through a preset frequency domain weighting transformation coefficient calculation formula. The frequency domain weighting transformation coefficients are used to characterize the contribution of the temporal changes of different Taylor coefficient orders at each joint observation frequency point.

[0128] Frequency point offset compensation is performed on the frequency domain weighted transform coefficients based on the frequency domain characteristic state variables, and then an extended frequency domain model is constructed based on the compensated frequency domain weighted transform coefficients.

[0129] The preliminary dynamic phasor parameter vector and the corresponding preliminary dynamic phasor Taylor coefficients of the current window are obtained by solving the extended frequency domain model using the least squares method. Based on the preliminary dynamic phasor parameter vector and the preliminary dynamic phasor Taylor coefficients, the preliminary frequency estimate of the current window is calculated.

[0130] Specifically, when calculating the frequency weights of each jointly observed frequency point, a pre-defined frequency weight formula can be used. This formula aims to quantify the relative importance of different jointly observed frequencies in frequency estimation. For example, weights can be designed based on the distance of the frequency point from the reference frequency, the amplitude of the frequency domain response, or the phase characteristics. Typically, frequencies closer to the reference frequency or with larger amplitude values ​​may be assigned higher weights to highlight their dominant role in frequency estimation. Subsequently, the weighted frequency domain response can be obtained by multiplying the original frequency domain response with the calculated frequency weights or performing other weighting operations. This process helps to enhance the information of frequencies that contribute more to frequency estimation while suppressing the influence of noise or other unimportant frequencies.

[0131] Specifically, based on the frequency domain response obtained in S206 and the frequency domain bending state quantity obtained in S207 and frequency domain asymmetric state quantities Generate frequency weights And organize joint observation vectors This forms the basis for joint observation; specifically, the frequency point weights The calculation formula is as follows:

[0132] ;

[0133] in, , and For frequency point weighting parameters;

[0134] The real and imaginary parts of the weighted frequency domain responses corresponding to each joint observation frequency point are stacked in a preset order to form a real joint observation vector. Its expression is as follows:

[0135] ;

[0136] in, This indicates the operation of taking the imaginary part;

[0137] It should be noted that when the real and imaginary parts of the weighted frequency domain response corresponding to each joint observation frequency point are stacked in a preset order to form a real joint observation vector, the purpose is to convert the complex form of the weighted frequency domain response into a real vector to meet the needs of subsequent real domain mathematical processing such as least squares method.

[0138] Next, for the sampling sequence within the sliding time window, a Taylor series expansion is performed at the center of the window to construct a dynamic phasor model with local time variables as independent variables. Taylor series expansion can represent the instantaneous frequency, amplitude, and other dynamic characteristics of the signal as a series of Taylor coefficients, thereby accurately describing the dynamic behavior of the signal within a local time range. Simultaneously, based on the centered position variable of the sampling sequence and the curved adaptive window weighting function, the frequency domain weighting transformation coefficients of each order are obtained through a preset formula. The centered position variable represents the relative time position of each sampling point in the sampling sequence with respect to the center of the window, and is usually normalized. The curved adaptive window weighting function, as described above, can weight the sampling sequence based on prior estimates (including offset estimates, frequency domain curved state variables, and frequency domain asymmetric state variables), thereby allowing the window shape to adaptively adjust to better adapt to the dynamic characteristics of the signal, effectively reducing spectral leakage and improving frequency estimation accuracy. These frequency-domain weighted transform coefficients are the key bridge connecting the time-domain Taylor coefficients and the frequency-domain observations. They are obtained through pre-defined formulas (e.g., based on the properties of the Fourier transform and the characteristics of the window function) and are used to characterize the contribution of the temporal variations of different Taylor coefficient orders to the frequency-domain response at each joint observation frequency point.

[0139] For example, by performing a Taylor series expansion on the signal within the sliding analysis window at the center of the window, a model can be established using local time variables. Dynamic phasor model with independent variables Take the centralized position variable obtained above. and curved adaptive window weighting function The k-th order frequency domain weighted transform coefficients are calculated using the following formula. This coefficient is used to characterize the contribution of different Taylor coefficient orders in temporal variations at each joint observation frequency point; specifically, the dynamic phasor and the k-th order frequency domain weighted transform coefficients The calculation formula is as follows:

[0140] ;

[0141] in, This refers to the local discrete-time variable relative to the center of the sliding analysis window, i.e., the sampling offset. Let be the Taylor coefficient of the k-th dynamic phasor, where K is the highest order of the Taylor expansion of the dynamic phasor.

[0142] Next, frequency offset compensation is performed on the frequency-domain weighted transform coefficients based on the frequency-domain characteristic state variables. Since the actual signal frequency may deviate from the nominal frequency, the spectral lines of the frequency-domain response will shift. Frequency-domain characteristic state variables (such as frequency-domain bending state variables and frequency-domain asymmetric state variables) provide information about this shift. Compensating the frequency-domain weighted transform coefficients using these state variables can correct the model mismatch caused by the frequency shift, making the constructed extended frequency-domain model more accurately reflect the frequency-domain characteristics of the actual signal. Based on the compensated frequency-domain weighted transform coefficients, the extended frequency-domain model can be constructed. This model is a mathematical model that correlates the weighted frequency-domain response (as observations) with the dynamic phasor parameter vector to be solved (containing Taylor coefficients) through the compensated frequency-domain weighted transform coefficients. It is typically expressed as a system of linear equations, laying the foundation for subsequent least-squares solutions.

[0143] For example, based on the frequency domain weighted transform coefficients, frequency domain bending state variables, and frequency domain asymmetric state variables obtained from the above steps, a normalized frequency offset is generated. Frequency domain bending compensation function and the compensated frequency domain weighted transform coefficients This forms the compensated frequency domain weighted transform coefficients; specifically, the normalized frequency point offset... Frequency domain bending compensation function and the compensated frequency domain weighted transform coefficients The calculation formula is as follows:

[0144] ;

[0145] ;

[0146] ;

[0147] in, and These are the coefficients of the compensation function.

[0148] Based on the obtained compensated frequency domain weighted transform coefficients Construct the main observation row vector Mirror observation row vector and the real observation submatrix corresponding to a single frequency point This forms the submatrix corresponding to the extended observation matrix; specifically, the main observation row vector Mirror observation row vector and real observation submatrix They are as follows:

[0149] ;

[0150] ;

[0151] remember , , and They are respectively of length , and respectively represent real row vectors, and respectively represent and If the real part and imaginary part are given, then:

[0152] ;

[0153] in, It is a 2×(2K+2) real observation submatrix;

[0154] Next, the obtained real-valued observation submatrices at each frequency point are combined to form an extended observation matrix. And establish the dynamic phasor parameter vector to be identified. With dynamic phasor Taylor coefficient The correspondence is used to form a parameterized model corresponding to the unified modeling relationship; specifically, the dynamic phasor parameter vector to be identified With dynamic phasor Taylor coefficient The correspondence is as follows:

[0155] ;

[0156] Next, the obtained real-number joint observation vector The extended observation matrix obtained with S212 and the dynamic phasor parameter vector to be identified The model is unified into an extended frequency domain model, forming a joint optimization solution for the corresponding observation equations; specifically, the extended frequency domain model is as follows:

[0157] ;

[0158] in, and Let be the real and imaginary parameters of the Taylor coefficients of the k-th dynamic phasor, respectively. It is formed by stacking 2×(2K+2) real observation submatrices corresponding to each frequency point along the row direction. The real-valued extended observation matrix, where Z is the number of rows in the extended observation matrix. The real number joint observation vector, where A is a (2K+2)×1 real number parameter vector.

[0159] Finally, the extended frequency domain model is applied using the least squares method. The solution is performed to obtain the initial dynamic phasor parameter vector of the current window and the corresponding initial dynamic phasor Taylor coefficients. and .

[0160] It should be noted that the least squares method is an optimization method aimed at finding a set of parameters that minimizes the sum of squared residuals between the model's predicted values ​​and the actual observed values, thereby robustly obtaining the dynamic phasor parameter vector and Taylor coefficients that best fit the weighted frequency domain response. These solutions contain dynamic information such as the amplitude, phase, frequency, and rate of change of the signal within the current window. Based on the obtained preliminary dynamic phasor parameter vector and preliminary dynamic phasor Taylor coefficients, the preliminary frequency estimate for the current window can be calculated according to the definition of the dynamic phasor model. For example, the zeroth-order Taylor coefficient is usually related to the instantaneous amplitude and phase of the signal, while the first-order Taylor coefficient is related to the instantaneous frequency shift. These coefficients can be used to accurately derive the frequency estimate.

[0161] More specifically, the formula for calculating the preliminary frequency estimate of the current window is as follows:

[0162] ;

[0163] Let the global sampling time index corresponding to the center of the current window be . The preliminary frequency estimate is then used as... Stored as an identifier This serves as the frequency prior for the next window; the sliding analysis window is advanced by one sampling interval along the time axis, shifting the window center to the next sampling time. The process returns to the step of determining the prior estimate for the current window, continuing until all sampling times corresponding to the current phase data have been processed, thereby forming a sequence of preliminary frequency estimates for each sampling time. .

[0164] Through the above technical solution, this application calculates the weighted frequency domain response by introducing a frequency point weighting formula and converts it into a real-number joint observation vector. This allows for more effective utilization of frequency domain information, highlighting key frequencies and suppressing noise. Simultaneously, a Taylor series expansion is performed at the center of the window to construct a dynamic phasor model. Combined with a curved adaptive window weighting function and frequency domain weighted transform coefficients, the dynamic characteristics of the signal can be accurately characterized. Furthermore, frequency point offset compensation is applied to the frequency domain weighted transform coefficients using frequency domain characteristic state variables, effectively correcting the model error caused by frequency offset, enabling the extended frequency domain model to more accurately reflect the actual signal. Finally, solving the extended frequency domain model using the least squares method robustly yields the preliminary dynamic phasor parameter vector and Taylor coefficients, thus calculating a high-precision preliminary frequency estimate. Compared to estimation based solely on the frequency domain response and characteristic state variables, this solution, through refined frequency domain weighting, dynamic modeling, and frequency point offset compensation, significantly improves the accuracy and robustness of the preliminary frequency estimate under dynamic and noisy environments, providing more reliable basic data for subsequent inertia identification.

[0165] In some of the above embodiments, this application obtains the dynamic phasor parameter vector and frequency estimation results for each corresponding time window by solving the extended frequency domain model. Based on this, the application further proposes that after obtaining the dynamic phasor parameter vector and frequency estimation results for each corresponding time window, the following steps are also included:

[0166] For the dynamic phasor parameter vectors of each phase within any time window, the dynamic phasor parameter vectors are restored to the Taylor coefficients of the corresponding dynamic phasors of each order. A positive-order transformation is performed on each of the Taylor coefficients to obtain a unified positive-order dynamic phasor Taylor coefficient set. Based on this unified positive-order dynamic phasor Taylor coefficient set, the zero-order and first-order positive-order dynamic phasor Taylor coefficients are extracted. Combined with a preset unified principal frequency calculation formula, the unified principal frequency is calculated. Then, through iterative processing, the unified principal frequency corresponding to each time window is calculated sequentially to obtain a unified principal frequency sequence. This unified principal frequency sequence is used as the frequency sequence to be cleaned, and the median reference frequency is extracted from this sequence. Based on the median reference frequency and combined with the peak decision threshold, the peak points in the sequence to be cleaned are corrected to obtain the frequency estimate after peak filtering. A band-constrained low-pass filter is performed on the frequency estimate after peak filtering and the corresponding power disturbance sequence to obtain the filtered power disturbance sequence and the frequency change rate sequence corresponding to different window paths, wherein the window paths include: Hann window path and rectangular window path. Based on the filtered power disturbance sequence and the frequency change rate sequence, an inertia identification result characterizing the power system inertia level is calculated using a preset power system inertia identification result calculation formula, wherein the inertia identification result includes: power disturbance amplitude, extreme values ​​of frequency change rate corresponding to different window paths, and power system equivalent inertia identification result.

[0167] Specifically, restoring the dynamic phasor parameter vector to the Taylor coefficients of the corresponding phase dynamic phasors means converting the frequency domain parameters obtained by solving the extended frequency domain model into zero-order, first-order, second-order, etc. Taylor coefficients in the time domain through a preset mapping relationship or inverse transformation. These Taylor coefficients characterize the local dynamic characteristics of the signal within the time window and are the basis for subsequent forward-order transformation and frequency calculation.

[0168] Performing positive-sequence transformations on the Taylor coefficients of each dynamic phasor to obtain a unified set of positive-sequence dynamic phasor Taylor coefficients is intended to eliminate potential imbalances in three-phase power systems. Positive-sequence transformation is a method of decomposing three-phase imbalances into symmetrical components, where the positive-sequence components represent the system's primary operating modes. By performing positive-sequence transformations on the Taylor coefficients of each phase, positive-sequence Taylor coefficients reflecting the overall characteristics of the system can be obtained, thus providing accurate input for calculating the unified system dominant frequency.

[0169] Based on the unified positive-sequence dynamic phasor Taylor coefficient set, the zero-order and first-order positive-sequence dynamic phasor Taylor coefficients are extracted. Combined with a pre-defined unified dominant frequency calculation formula, the unified dominant frequency is calculated. This process is then iteratively repeated to calculate the unified dominant frequency for each time window, resulting in a unified dominant frequency sequence. The zero-order positive-sequence dynamic phasor Taylor coefficient typically corresponds to the steady-state component or average value of the frequency, while the first-order positive-sequence dynamic phasor Taylor coefficient is related to the rate of change of the frequency. The pre-defined unified dominant frequency calculation formula comprehensively utilizes these two coefficients to calculate the unified dominant frequency for the current time window. For example, the zero-order coefficient can be used as a reference frequency, combined with the first-order coefficient to correct or predict the frequency. By performing this process on all time windows, a continuous unified dominant frequency sequence reflecting the overall frequency change of the system can be obtained.

[0170] A unified master frequency sequence is used as the frequency sequence to be cleaned. A median reference frequency is extracted from this sequence, and peaks in the sequence are corrected based on this median reference frequency and a peak decision threshold, resulting in peak-filtered frequency estimates. The unified master frequency sequence may still contain outliers caused by measurement noise or transient events. The median reference frequency is an effective method for smoothing a sequence and is insensitive to outliers. By calculating the median of the sequence as a reference and setting a peak decision threshold (e.g., a deviation from the median reference frequency exceeding a preset percentage), peaks exceeding the threshold can be identified and corrected. Correction methods can include replacing the median reference frequency, the previous valid value, or a value obtained through interpolation, thereby obtaining a smoother and more reliable peak-filtered frequency estimate sequence.

[0171] Band-constrained low-pass filtering is performed on the frequency estimates and corresponding power disturbance sequences after peak filtering to obtain the filtered power disturbance sequence and frequency change rate sequences corresponding to different window paths, including Hann window paths and rectangular window paths. The frequency estimates and power disturbance sequences after peak filtering may still contain high-frequency noise, which can affect the accuracy of inertia identification. Band-constrained low-pass filtering aims to remove these high-frequency components while retaining low-frequency information related to inertia dynamics. Low-pass filters can be implemented in various ways, such as Butterworth filters and Chebyshev filters, and their cutoff frequencies should be set according to the typical frequency range of power system inertia dynamics. Different window paths (such as Hann windows and rectangular windows) are used to weight the data in signal processing to reduce spectral leakage. Hann windows have a smooth transition band and can effectively suppress sidelobes, while rectangular windows have a narrower main lobe but higher sidelobes. By comparing the frequency change rate sequences under different window paths, the robustness of the identification results and the sensitivity to window function selection can be evaluated.

[0172] Based on the filtered power disturbance sequence and frequency change rate sequence, an inertia identification result characterizing the power system's inertia level is calculated using a pre-defined formula for power system inertia identification. This inertia identification result includes: power disturbance amplitude, extreme values ​​of frequency change rates corresponding to different window paths, and the equivalent inertia identification result of the power system. The core of inertia identification lies in establishing the relationship between power disturbance and frequency change rate. The pre-defined formula for calculating the power system inertia identification result is typically based on a simplified form of the power angle equation of the power system. By analyzing the filtered power disturbance sequence and frequency change rate sequence, the power disturbance amplitude (e.g., the maximum power change after the disturbance occurs) and extreme values ​​of the frequency change rate (e.g., the maximum rate of frequency decrease or increase) can be extracted. Substituting these parameters into the inertia identification formula yields the equivalent inertia identification result of the power system, which is a key indicator for measuring the power system's ability to resist frequency changes.

[0173] For example, the dynamic phasor parameter vectors of phases A, B, and C at sampling time i obtained from the preceding steps. , and The dynamic phasor parameter vectors of each phase are restored to the Taylor coefficients of the corresponding dynamic phasors of each phase. , and This forms the Taylor coefficient set of dynamic phasors for phases A, B, and C.

[0174] Next, based on the obtained Taylor coefficient sets of dynamic phasors for phases A, B, and C, a positive-sequence transformation is performed on the Taylor coefficients of each order of dynamic phasors to form a unified positive-sequence Taylor coefficient set; specifically, the Taylor coefficients of the k-th order positive-sequence dynamic phasors at sampling time i... The calculation formula is as follows:

[0175] ;

[0176] in, It is the orthogonal rotation factor.

[0177] Next, based on the obtained set of Taylor coefficients for the positive-order dynamic phasors, the Taylor coefficients for the zeroth-order positive-order dynamic phasors are extracted. Taylor coefficients of the first-order positive-sequence dynamical phasor And calculate the unified main frequency result. Specifically, the unified main frequency result The calculation formula is as follows:

[0178] ;

[0179] Next, the above steps are repeated for all sampling times to obtain a unified main frequency sequence composed of the unified main frequency results from each sampling time. This serves as the input required for constructing the frequency sequence to be cleaned in subsequent steps.

[0180] Next, the obtained unified main frequency sequence As a frequency sequence to be cleaned, a median reference frequency is generated based on the frequency sequence to be cleaned. The frequency reference value required to form the peak decision; specifically, the median reference frequency. The calculation formula is as follows:

[0181] ;

[0182] in, This is the half-window width for median filtering. This indicates taking the midpoint of the set;

[0183] Based on the obtained median reference frequency and the frequency sequence to be cleaned, a consistency ratio coefficient is generated. and peak decision threshold This forms the decision threshold required for peak filtering; specifically, the consistency ratio coefficient... and peak decision threshold The calculation formula is as follows:

[0184] ;

[0185] in, The cumulative distribution function of the standard normal distribution. The quantile function represents the standard normal distribution. These are the upper quartiles of the standard normal distribution. ≈0.6745; This is the lower limit of the system's preset peak decision threshold;

[0186] Next, based on the obtained median reference frequency and peak decision threshold, the peaks in the frequency sequence to be cleaned are corrected to obtain the frequency estimate after peak removal. To perform frequency band constrained filtering; specifically, the frequency estimate after peak filtering... The calculation formula is as follows:

[0187] ;

[0188] Based on the frequency estimate after peak filtering obtained from S603 and its corresponding power perturbation sequence, a band-constrained low-pass filter is performed to obtain the filtered power perturbation sequence. and the frequency change rate sequences corresponding to the Hann window path and the rectangular window path. , In order to perform inertia identification; specifically, the band-constrained low-pass filter The transfer function is as follows:

[0189] , , ;

[0190] Where s is the Laplace variable, This is the cutoff angular frequency of the low-pass filter. This is the cutoff frequency of the low-pass filter;

[0191] Finally, based on the obtained filtered power perturbation sequence And the frequency change rate sequences corresponding to the Hann window path and the rectangular window path are used to calculate the power perturbation amplitude. Extreme values ​​of frequency change rate and and the results of power system equivalent inertia identification This forms an inertia identification result characterizing the inertia level of the power system; specifically, the power disturbance amplitude... Extreme values ​​of frequency change rate and and the results of power system equivalent inertia identification The calculation formula is as follows:

[0192] ;

[0193] ;

[0194] , ;

[0195] in, and The start and end sampling indices are extracted for the time interval corresponding to the power plateau before the disturbance. and This provides the start and end sampling indices for the time interval corresponding to the power platform extraction after the disturbance. and This refers to the start and end sampling indices corresponding to the search time interval for the extreme value of the rate of change. The inertia identification results are obtained based on the Hann window path. The inertia identification results are obtained based on the rectangular window path.

[0196] This embodiment effectively suppresses the impact of frequency measurement jitter on inertia calculation through a dual-window path fusion and multi-layer filtering mechanism, improving the stability and reliability of inertia identification results in engineering applications. Specifically, firstly, peak filtering based on the median reference frequency is performed on the unified main frequency sequence to eliminate short-term abnormal fluctuations; then, a smoothed power disturbance sequence and frequency change rate sequence are obtained through band-constrained low-pass filtering; finally, the inertia identification results are calculated based on the Hann window path and the rectangular window path respectively, and the average value is fused, thereby reducing the sensitivity of a single window function path to frequency jitter and ensuring the robustness of the equivalent inertia identification results under noisy conditions.

[0197] The above is a detailed description of an embodiment of a power system inertia identification method provided in this application. The following is a detailed description of an embodiment of a power system inertia identification device provided in this application.

[0198] Please see Figure 2 This application provides a power system inertia identification device, comprising:

[0199] Electrical quantity acquisition unit 201 is used to acquire three-phase electrical quantity data at the generator node grid connection point;

[0200] The preliminary frequency estimation unit 202 is used to calculate the frequency domain response at each joint observation frequency point for each phase of the three-phase electrical quantity data through a preset sliding time window and in combination with several joint observation frequency points set at equal intervals within the current window. Based on the frequency domain response, the unit determines the frequency domain characteristic state quantity of the current window. Based on the frequency domain response and the frequency domain characteristic state quantity, the unit calculates the preliminary frequency estimate of the current window. The unit obtains the preliminary frequency estimate of each window by moving the sliding time window and performing cyclic calculations.

[0201] The power disturbance characterization calculation unit 203 is used to calculate the power disturbance characterization quantity at the center sampling time based on the center sampling time of each window and the three-phase electrical quantity data corresponding to the center sampling time.

[0202] The frequency change rate representation vector calculation unit 204 is used to determine two adjacent time windows of the selected target time window, and construct a frequency change rate representation vector based on the preliminary frequency estimates of the adjacent time windows.

[0203] The power system inertia identification unit 205 is used to construct a unified objective function containing inertia coupling terms based on the power disturbance characterization quantity and the frequency change rate characterization vector of the target time window, combined with a preset extended frequency domain model. By solving the unified objective function, the corresponding dynamic phasor parameter vector and frequency estimation result are obtained for each corresponding dynamic phasor parameter in the target time window. The extended frequency domain model is configured with the correspondence between the dynamic phasor parameters to be identified and Taylor coefficients.

[0204] In the specific implementation process, the electrical quantity acquisition unit first acquires the three-phase electrical quantity data at the generator node grid connection point, which serves as the basic input for subsequent processing. The preliminary frequency estimation unit extracts data segments for each phase using a sliding time window, setting equally spaced joint observation frequency points within the window. It calculates the frequency domain response at each frequency point through frequency domain analysis. Based on this frequency domain response, it further determines frequency domain characteristic state quantities such as frequency domain bending state quantities and frequency domain asymmetric state quantities. Then, it combines the frequency domain response and frequency domain characteristic state quantities to calculate the preliminary frequency estimate, thus forming a continuous frequency estimation sequence. The power disturbance characterization calculation unit uses the three-phase electrical quantity data corresponding to the sampling time at the center of each window to calculate the power disturbance characterization quantity at that time. This characterization quantity reflects the instantaneous characteristics of the system power disturbance. The frequency change rate characterization vector calculation unit constructs a frequency change rate characterization vector by analyzing the preliminary frequency estimates of adjacent windows of the target time window to quantify the dynamic trend of frequency change. Finally, the power system inertia identification unit inputs the power disturbance characterization quantity and the frequency change rate characterization vector into the extended frequency domain model to construct a unified objective function containing inertia coupling terms. By solving this objective function, the dynamic phasor parameter vector and frequency estimation results are obtained simultaneously. The extended frequency domain model realizes the inverse constraint of inertia information on the frequency identification process through the correspondence between the dynamic phasor parameters to be identified and Taylor coefficients.

[0205] Through the above technical solution, this device integrates frequency identification and inertia calculation into a unified optimization process, avoiding the limitation of unidirectional information flow transmission in the traditional serial processing method. It effectively solves the problem of lack of physical consistency under dynamic disturbance scenarios, significantly improves the stability and engineering applicability of inertia identification results, and provides reliable technical support for the stable control of power systems under scenarios such as large disturbances, black start and high proportion of new energy access.

[0206] In some embodiments, this application proposes a power system inertia identification terminal, which effectively solves the above-mentioned problems through specific hardware configuration. For example... Figure 3 As shown, the terminal includes a memory 33 and a processor 31, which can be connected via a communication bus 34.

[0207] The memory 33 is a hardware device used to store data and instructions. In the power system inertia identification terminal, the memory is mainly used to store program code, which corresponds to the power system inertia identification method described above. It is also used to store intermediate data and final results generated during the method execution. Specifically, the memory can use random access memory (RAM) as the working memory to provide high-speed data read and write capabilities, supplemented by large-capacity non-volatile memory, such as solid-state drives (SSDs) or flash memory, for long-term storage of program code and historical identification data, ensuring data persistence and reliability.

[0208] Processor 31 is a hardware unit that executes computer instructions and is the core of the computer system's computation and control. In the power system inertia identification terminal, the processor is responsible for reading the program code stored in the memory and executing the various steps of the power system inertia identification method according to the instructions in the program code. This includes preprocessing the collected three-phase electrical quantity data, performing complex frequency domain transformations and response calculations, performing frequency estimation, calculating power disturbance characterization quantities, and constructing and solving a unified objective function containing inertia coupling terms. To meet the computational performance and real-time requirements of the power system inertia identification method, the processor can be a high-performance multi-core central processing unit (CPU), graphics processing unit (GPU), or digital signal processor (DSP) to support parallel computing and floating-point operations. In some scenarios with extremely high real-time requirements, field-programmable gate arrays (FPGAs) or application-specific integrated circuits (ASICs) can also be used to implement hardware acceleration for specific algorithms.

[0209] In some embodiments, this application proposes a technical solution for carrying and executing the inertia identification method through a computer-readable storage medium.

[0210] Specifically, this application provides a computer-readable storage medium, which is a non-transitory medium capable of storing digital data and being read by a computer system. This medium may include, but is not limited to, read-only memory (ROM), random access memory (RAM), flash memory, electrically erasable programmable read-only memory (EEPROM), optical discs (such as CD-ROM, DVD, Blu-ray disc), magnetic storage (such as hard disk drive, floppy disk), or solid-state drives (SSD). The medium stores program code, which is a set of instructions written in a specific programming language to guide the computer in performing specific tasks. For example, the program code may exist in the form of machine code, assembly code, or a high-level language (such as C++, Java, Python) compiled or interpreted. When the program code is read and executed by a processor, the processor (e.g., a central processing unit (CPU), graphics processing unit (GPU), digital signal processor (DSP), microcontroller (MCU), or field-programmable gate array (FPGA) configured as a processor) will implement the power system inertia identification method described above according to the preset logic and algorithms in the program code. This means that the entire complex process, from collecting three-phase electrical quantity data at the generator node grid connection point, to calculating the frequency domain response at each joint observation frequency point, determining the frequency domain characteristic state quantity, calculating the preliminary frequency estimate, then calculating the power disturbance characterization quantity, constructing the frequency change rate characterization vector, solving the unified objective function containing inertia coupling terms, and subsequent dynamic phasor parameter restoration, unified main frequency calculation, frequency sequence cleaning, frequency band constraint low-pass filtering, until finally obtaining the power system inertia identification result, can be automatically and efficiently completed on the processor through the program code.

[0211] The above technical solution stores the power system inertia identification method in the form of program code on a computer-readable storage medium, which is then read and executed by a processor, achieving automated and efficient operation of this complex method. This implementation significantly improves the deployment flexibility and maintainability of the method, eliminating the need for hardware modifications for each method update; only the program code needs to be updated. Simultaneously, the processor's high-speed computing power ensures the real-time performance and accuracy of the inertia identification process, effectively avoiding errors that may be introduced by manual operation, thus providing reliable inertia information for the stable operation and control of the power system. Furthermore, the use of standardized computer hardware and software interfaces reduces the system's development and maintenance costs, enabling this advanced inertia identification technology to be more widely applied in practical power systems.

[0212] The above-described embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit them. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of this application.

Claims

1. A method for identifying the inertia of a power system, characterized in that, include: Collect three-phase electrical quantity data at the generator node grid connection point; For each phase of the three-phase electrical quantity data, the frequency domain response at each joint observation frequency point is calculated by using a preset sliding time window and combining several joint observation frequency points set at equal intervals within the current window. Based on the frequency domain response, the frequency domain characteristic state quantity of the current window is determined. Then, based on the frequency domain response and the frequency domain characteristic state quantity, the preliminary frequency estimate of the current window is calculated. By moving the sliding time window and performing cyclic calculations, the preliminary frequency estimate of each window is obtained sequentially. Based on the center sampling time of each window, and combined with the three-phase electrical quantity data corresponding to the center sampling time, the power disturbance characterization quantity at the center sampling time is calculated; For a selected target time window, two adjacent time windows are determined, and a frequency change rate representation vector is constructed based on the preliminary frequency estimates of the adjacent time windows. Based on the power disturbance characterization quantity and the frequency change rate characterization vector of the target time window, and combined with the preset extended frequency domain model, a unified objective function containing inertia coupling terms is constructed. By solving the unified objective function, the corresponding dynamic phasor parameter vector and frequency estimation result for each corresponding dynamic phasor within the target time window are obtained. The extended frequency domain model is configured with the correspondence between the dynamic phasor parameters to be identified and the dynamic phasor Taylor coefficients. For the dynamic phasor parameter vector of each phase within any time window, the dynamic phasor parameter vector is restored to the Taylor coefficients of the corresponding dynamic phasors of each order. Based on the Taylor coefficients of each order of dynamic phasors, a positive order transformation is performed to obtain a unified positive order set of Taylor coefficients of dynamic phasors. Based on the unified positive-sequence dynamic phasor Taylor coefficient set, the zero-order positive-sequence dynamic phasor Taylor coefficient and the first-order positive-sequence dynamic phasor Taylor coefficient are extracted, and combined with the preset unified main frequency calculation formula, the unified main frequency is calculated. Then, through iterative processing, the unified main frequency corresponding to each time window is calculated in sequence to obtain the unified main frequency sequence. The unified main frequency sequence is used as the frequency sequence to be cleaned. The median reference frequency is extracted from the frequency sequence to be cleaned. Based on the median reference frequency and the peak decision threshold, the peak points in the frequency sequence to be cleaned are corrected to obtain the frequency estimate after peak filtering. A band-constrained low-pass filter is performed on the frequency estimate after peak filtering and the corresponding power perturbation sequence to obtain the filtered power perturbation sequence and the frequency change rate sequence corresponding to different window paths, wherein the window paths include: Hann window path and rectangular window path; Based on the filtered power disturbance sequence and the frequency change rate sequence, the inertia identification result used to characterize the power system inertia level is calculated using a preset power system inertia identification result calculation formula. The inertia identification result includes: power disturbance amplitude, extreme values ​​of frequency change rate corresponding to different window paths, and power system equivalent inertia identification result.

2. The method for power system inertia identification according to claim 1, characterized in that, The frequency domain characteristic state quantities include: frequency domain bending state quantities and frequency domain asymmetric state quantities; The formula for calculating the frequency domain bending state quantity is as follows: ； The formula for calculating the frequency domain asymmetric state quantity is: ； In the formula, This refers to the frequency domain bending state quantity. The frequency domain asymmetric state quantity, This indicates the operation of taking the real part. To prevent positive numbers with a denominator of zero, As the reference frequency, This is the frequency offset. This represents the frequency domain response.

3. The method for power system inertia identification according to claim 2, characterized in that, Before calculating the frequency domain response at each joint observation frequency point, which is set at equal intervals within the current window, the following steps are also included: For each phase of the three-phase electrical quantity data, a sampling sequence within the current window is extracted through a preset sliding time window, and a normalized position variable is determined based on the sampling sequence. Obtain the prior estimate value corresponding to the current window. When the current window is not the first window, the prior estimate value is the offset estimate value and frequency domain characteristic state quantity of the previous window. The offset estimate value is the relative deviation between the preliminary frequency estimate value and the nominal frequency of the corresponding window. When the current window is the first window, the prior estimate value is the obtained preset value. Substituting the normalized location variable and the prior estimate into a preset bending adaptive window weighting function, the window weighting coefficients corresponding to each time window are obtained. The sampling sequence is weighted point by point according to the window weighting coefficients to obtain a weighted sampling sequence.

4. The method for power system inertia identification according to claim 3, characterized in that, The specific bending adaptive window weighting function is as follows: ； In the formula, n and q both represent the index of the sampling point, and N is the total number of sampling points within the window. For the normalized position variable corresponding to the nth sampling point, Let q be the normalized position variable corresponding to the sampling point. , and This is the window adjustment coefficient. The offset prediction value in the prior estimate, The frequency domain bending state quantity in the prior estimate. The frequency domain asymmetric state quantity in the prior estimate. This is an adjustable window adaptive coefficient.

5. The method for power system inertia identification according to claim 3, characterized in that, The calculation of the preliminary frequency estimate for the current window based on the frequency domain response and the frequency domain characteristic state quantities includes: Based on the frequency domain response and the frequency domain characteristic state quantity, the frequency weight of each joint observation frequency point is calculated using a preset frequency weight formula, and the weighted frequency domain response is calculated based on the frequency weight and the frequency domain response. The real and imaginary parts of the weighted frequency domain response corresponding to each joint observation frequency point are stacked in a preset order to form a real joint observation vector; For the sampling sequence within the sliding time window, a Taylor series expansion is performed at the center of the window to construct a dynamic phasor model with local time variables as independent variables. Then, based on the centered position variable of the sampling sequence and the bending adaptive window weighting function, the frequency domain weighting transformation coefficients of each order are obtained through a preset frequency domain weighting transformation coefficient calculation formula. The frequency domain weighting transformation coefficients are used to characterize the contribution of the temporal changes of different dynamic phasor Taylor coefficient orders at each joint observation frequency point. Frequency point offset compensation is performed on the frequency domain weighted transform coefficients based on the frequency domain characteristic state variables, and then an extended frequency domain model is constructed based on the compensated frequency domain weighted transform coefficients. The preliminary dynamic phasor parameter vector and the corresponding preliminary dynamic phasor Taylor coefficients of the current window are obtained by solving the extended frequency domain model using the least squares method. Based on the preliminary dynamic phasor parameter vector and the preliminary dynamic phasor Taylor coefficients, the preliminary frequency estimate of the current window is calculated.

6. The method for power system inertia identification according to claim 5, characterized in that, Based on the compensated frequency domain weighted transform coefficients, the extended frequency domain model is constructed as follows: Based on the compensated frequency domain weighted transform coefficients, the main observation row vector, the mirror observation row vector, and the real observation submatrix corresponding to a single frequency point are constructed, and the real observation submatrixes of each frequency point are combined to form an extended observation matrix; Obtain the correspondence between the dynamic phasor parameters to be identified and the Taylor coefficients of the dynamic phasor; The real-valued joint observation vector, the extended observation matrix, and the dynamic phasor parameters to be identified are used to construct the corresponding observation equation for joint optimization, thereby obtaining the extended frequency domain model.

7. A power system inertia identification device, characterized in that, include: The electrical quantity acquisition unit is used to collect three-phase electrical quantity data at the generator node grid connection point; The preliminary frequency estimation unit is used to calculate the frequency domain response at each joint observation frequency point for each phase of the three-phase electrical quantity data through a preset sliding time window and in combination with several joint observation frequency points set at equal intervals within the current window. Based on the frequency domain response, the unit determines the frequency domain characteristic state quantity of the current window. Based on the frequency domain response and the frequency domain characteristic state quantity, the unit calculates the preliminary frequency estimate of the current window. By moving the sliding time window and performing cyclic calculations, the unit obtains the preliminary frequency estimate of each window in turn. The power disturbance characterization calculation unit is used to calculate the power disturbance characterization quantity at the center sampling time based on the center sampling time of each window and the three-phase electrical quantity data corresponding to the center sampling time. The frequency change rate representation vector calculation unit is used to determine two adjacent time windows of the selected target time window, and construct a frequency change rate representation vector based on the preliminary frequency estimates of the adjacent time windows. The power system inertia identification unit is used to construct a unified objective function containing inertia coupling terms based on the power disturbance characterization quantity and the frequency change rate characterization vector of the target time window, combined with a preset extended frequency domain model. By solving the unified objective function, the unit obtains the dynamic phasor parameter vector and frequency estimation result for each corresponding dynamic phasor within the target time window. The extended frequency domain model is configured with the correspondence between the dynamic phasor parameters to be identified and the dynamic phasor Taylor coefficients. For the dynamic phasor parameter vector of each phase within any time window, the dynamic phasor parameter vector is restored to the Taylor coefficients of the corresponding dynamic phasors of each order. Based on the Taylor coefficients of each order of dynamic phasors, a positive order transformation is performed to obtain a unified positive order set of Taylor coefficients of dynamic phasors. Based on the unified positive-sequence dynamic phasor Taylor coefficient set, the zero-order positive-sequence dynamic phasor Taylor coefficient and the first-order positive-sequence dynamic phasor Taylor coefficient are extracted, and combined with the preset unified main frequency calculation formula, the unified main frequency is calculated. Then, through iterative processing, the unified main frequency corresponding to each time window is calculated in sequence to obtain the unified main frequency sequence. The unified main frequency sequence is used as the frequency sequence to be cleaned. The median reference frequency is extracted from the frequency sequence to be cleaned. Based on the median reference frequency and the peak decision threshold, the peak points in the frequency sequence to be cleaned are corrected to obtain the frequency estimate after peak filtering. A band-constrained low-pass filter is performed on the frequency estimate after peak filtering and the corresponding power perturbation sequence to obtain the filtered power perturbation sequence and the frequency change rate sequence corresponding to different window paths, wherein the window paths include: Hann window path and rectangular window path; Based on the filtered power disturbance sequence and the frequency change rate sequence, the inertia identification result used to characterize the power system inertia level is calculated using a preset power system inertia identification result calculation formula. The inertia identification result includes: power disturbance amplitude, extreme values ​​of frequency change rate corresponding to different window paths, and power system equivalent inertia identification result.

8. A power system inertia identification terminal, characterized in that, include: Memory and processor; The memory is used to store program code, which corresponds to the power system inertia identification method as described in any one of claims 1 to 6; The processor is used to read and execute the program code to implement the power system inertia identification method.

9. A computer-readable storage medium, characterized in that, The computer-readable storage medium contains program code that is read and executed by a processor to implement the power system inertia identification method as described in any one of claims 1 to 6.

Images