Electric power system inertia online tracking method, system and equipment

Abstract

The invention discloses a power system inertia online tracking method, system and device, and relates to the technical field of power system frequency stability control, and the method comprises the steps: screening an inertia effective time period according to a frequency change rate and an energy coupling coefficient; performing synchronous extrusion wavelet transform on the signal in the effective time period, identifying a dominant low-frequency energy ridge line and constructing an adaptive transmission band, and performing mask filtering and inverse transformation reconstruction to obtain a noise reduction signal; constructing a difference equation based on the noise reduction signal, and recursively identifying a pulse transfer function according to the order of the width dynamic constraint equation of the adaptive transmission band; performing continuous domain conversion on the pulse transfer function, identifying a dominant inertia mode through a damping ratio-frequency composite criterion, and calculating an inertia time constant estimation value; and constructing a dynamic Bayesian network, carrying out iterative solution through maximum likelihood estimation, and outputting an inertia time constant correction value at the current moment. According to the invention, inertial precise continuous tracking of the high-proportion new energy power system under the natural micro-disturbance condition is realized.

Description

Technical Field

[0001] This invention relates to the field of power system frequency stability control technology, and more specifically, to a power system inertial online tracking method, system, and device. Background Technology

[0002] With the large-scale grid connection of new energy power generation and storage devices such as wind power and photovoltaics through power electronic converters, traditional synchronous generators are being gradually replaced, leading to a significant decrease in the equivalent inertia level of the power system. Simultaneously, new energy power generation and storage devices can provide inertial support to the system through virtual inertial control, causing the system's equivalent inertial constant to exhibit significant time-varying and uncertainties. Therefore, online and efficient tracking of the power system's inertia level has become a key technical means to ensure the frequency stability of power systems with a high proportion of new energy sources.

[0003] Existing power system inertial online tracking mainly identifies the first-order inertial time constant by fitting the transfer function between node injected power and frequency response using real-time measurement data.

[0004] For example, the invention patent application CN118316071A discloses a method and system for online identification of power system inertia. The method includes: under power system noise disturbance, synchronously collecting active power and frequency data at the connection bus of generating devices during normal grid operation and preprocessing them to obtain active power change and frequency change data; performing preliminary parameter identification based on the active power-frequency change data to obtain a preliminary state-space model describing the active power-frequency response process of the power system; optimizing the preliminary state-space model to obtain an optimized state-space model; obtaining a first-order state-space model based on the optimized state-space model; solving for the first inertial time constant of each generating device based on the first-order state-space model; obtaining the second inertial time constant of each generating device in the time period to be determined; and reconstructing the node frequency change data using an empirical mode decomposition filter to obtain node inertia metrics and equivalent node inertial time constants. This invention can accurately identify the inertial time constants at the generating device end and obtain the equivalent node inertial time constants without understanding the internal structure of the system, thus understanding the spatiotemporal characteristics of inertia distribution within the system.

[0005] For example, the invention patent application CN118427782A discloses a method for inertia extraction and correction based on WRARMAX and LSTM. Specifically, this method includes: acquiring active power and frequency data of a power system at different times; preprocessing the active power and frequency data to obtain active power and frequency offset value sequences; constructing an ARMAX model; identifying the active power and frequency offset value sequences using the recursive least squares method to obtain the discrete transfer function of the ARMAX model; transforming the discrete transfer function of the ARMAX model into a continuous transfer function using the bilinear transformation method; obtaining step response curves at different times based on the continuous transfer function; determining the inertial time constant of the power system at different times based on the step response curves to obtain an inertial time constant sequence; and constructing an LSTM model to correct the inertial time constant sequence. This approach solves the problems of ARMAX models failing to reflect the correlation between time series and the high computational cost and time consumption caused by large amounts of input data and numerous calculations.

[0006] The above-disclosed technical solutions have at least the following technical problems: First, these methods rely on a pre-defined model structure (such as a first-order or higher-order transfer function), while the equivalent dynamic characteristics of high-proportion renewable energy power systems are complex and variable. A fixed model structure cannot accurately describe the real physical processes of the system, leading to model mismatch errors. Second, these methods usually require sufficiently strong perturbation excitation or active injection of test signals to stimulate system dynamics, making it difficult to achieve continuous inertial tracking under normal operating conditions with only natural micro-perturbations. Furthermore, transfer function identification methods treat inertia as a fixed parameter for extraction, ignoring the dynamic evolution characteristics of inertia on a continuous time scale, making it difficult to truly reflect the time-varying nature of the inertia of renewable energy systems.

[0007] To address the above problems, this invention proposes a solution. Summary of the Invention

[0008] To overcome the aforementioned deficiencies of the prior art, embodiments of the present invention provide a method, system, and device for online inertial tracking of power systems, which solves the problem that existing methods suffer from poor adaptability and difficulty in achieving continuous and accurate tracking due to the fixed model structure under natural micro-disturbance conditions.

[0009] To achieve the above objectives, the present invention provides the following technical solution: A method for online inertial tracking in a power system includes the following steps: screening effective inertial time periods based on the rate of change of frequency and energy coupling coefficient; performing synchronous squeezing wavelet transform on the signal within the effective time period to identify the dominant low-frequency energy ridge and construct an adaptive passband; obtaining a denoised signal through mask filtering and inverse transform reconstruction; constructing a difference equation based on the denoised signal; recursively identifying the pulse transfer function after considering the order of the dynamic constraint equation based on the width of the adaptive passband; performing continuous domain transformation on the pulse transfer function; identifying the dominant inertial mode and calculating the estimated value of the inertial time constant through the damping ratio-frequency composite criterion; constructing a dynamic Bayesian network; iteratively solving the network through maximum likelihood estimation; and outputting the corrected value of the inertial time constant at the current moment.

[0010] In a preferred embodiment, the selection of effective inertial time periods specifically involves: acquiring power and frequency signals from the power generation device; calculating the frequency change rate based on the frequency signal, and marking the disturbance trigger moment when its absolute value exceeds a preset threshold; forming a candidate response interval by extending forward and backward from the disturbance trigger moment; calculating the energy coupling coefficient between the change in active power and the frequency change rate within the candidate response interval, and selecting effective inertial time periods accordingly; and determining the usage duration of the time period in subsequent steps based on the disturbance amplitude and duration within the effective inertial time period.

[0011] In a preferred embodiment, constructing the adaptive passband includes: performing synchronous squeezing wavelet transform on the power and frequency signals to obtain the time-frequency energy distribution; identifying the dominant low-frequency energy ridge based on the energy spectral density; and constructing the adaptive passband with the energy ridge as the center.

[0012] In a preferred embodiment, the step of obtaining the denoised signal through mask filtering and inverse transform reconstruction specifically involves: using the adaptive passband to perform mask filtering on the time-frequency energy distribution, retaining the time-frequency coefficients within the passband, and eliminating out-of-band noise and non-inertial mode components; performing synchronous squeezing wavelet inverse transform on the filtered time-frequency coefficients to reconstruct the denoised power fluctuation component and frequency fluctuation component.

[0013] In a preferred embodiment, the recursive identification of the pulse transfer function specifically involves: constructing a difference equation based on the denoised power fluctuation component and frequency fluctuation component; limiting the search range of the order of the difference equation according to the width of the adaptive passband; within the search range, determining the order of the equation at the current time step using the information content criterion; and obtaining the pulse transfer function at the current time step by updating the equation parameters through online learning based on the determined equation order.

[0014] In a preferred embodiment, the continuous domain transformation includes: obtaining the adaptive passband boundary frequency of the current time step and the center frequency corresponding to the dominant low-frequency energy ridge; performing pre-distortion correction on the boundary frequency and the center frequency respectively according to a preset sampling period; and performing continuous domain transformation on the pulse transfer function model based on the corrected frequency parameters to obtain the distortion-corrected continuous domain transfer function.

[0015] In a preferred embodiment, identifying the dominant inertial mode and calculating the estimated inertial time constant specifically involves: decomposing the continuous domain transfer function into independent dynamic modal components characterized by different poles; identifying the dominant inertial mode that conforms to electromechanical oscillation characteristics from the independent dynamic modal components based on a composite criterion of damping ratio and natural oscillation frequency, and extracting the decay time constant corresponding to the mode; establishing an analytical mapping relationship between the decay time constant and the equivalent rotational inertia, and calculating the estimated inertial time constant for the current time step.

[0016] In a preferred embodiment, the construction of the dynamic Bayesian network, which is solved iteratively by maximum likelihood estimation, and outputs the correction value of the inertial time constant at the current moment, specifically involves: determining whether the current system is in a steady state or a transient state based on the instantaneous offset of the center frequency of the energy ridge between adjacent time steps, and constructing Bayesian networks with different prior distributions accordingly; obtaining the dominant mode energy residual and signal quality index at the current time step, and dynamically adjusting the likelihood probability distribution parameters of the Bayesian network observation nodes; using the usage duration and adaptive bandwidth as constraints, and solving iteratively by maximum likelihood estimation to output the correction value of the inertial time constant at the current moment.

[0017] The system of the online inertial tracking method for power systems includes: a data filtering module for filtering effective inertial time periods based on the frequency change rate and energy coupling coefficient; a time-frequency processing module for performing synchronous squeezing wavelet transform on the signal within the effective time period, identifying the dominant low-frequency energy ridge and constructing an adaptive passband, and obtaining a denoised signal through mask filtering and inverse transform reconstruction; a model identification module for constructing a difference equation based on the denoised signal, and recursively identifying the pulse transfer function after considering the order of the dynamic constraint equation based on the width of the adaptive passband; a modal analysis module for performing continuous domain transformation on the pulse transfer function, identifying the dominant inertial mode through the damping ratio-frequency composite criterion, and calculating the estimated value of the inertial time constant; and a result correction module for constructing a dynamic Bayesian network, iteratively solving it through maximum likelihood estimation, and outputting the corrected value of the inertial time constant at the current moment.

[0018] An electronic device, comprising a memory and a processor: the memory for storing a program; the processor for executing the program to implement the various steps of the power system inertial online tracking method.

[0019] The technical effects and advantages of the present invention regarding the online inertial tracking method, system, and equipment for power systems are as follows: 1. This invention extracts the dominant low-frequency energy ridge line by synchronous squeezing wavelet transform and constructs an adaptive passband. The time-frequency coefficients are then masked and filtered using this passband to remove out-of-band noise and non-inertial mode components. The denoised signal is then reconstructed through inverse transform. This solves the problems of low signal-to-noise ratio and difficulty in extracting inertial response signals under natural micro-perturbation conditions in existing methods, and provides a high-quality data foundation for subsequent model identification.

[0020] 2. This invention determines the system's operating state based on the instantaneous offset of the center frequency of the energy ridge between adjacent time steps, constructs a dynamic Bayesian network with different prior distributions, and introduces the dominant mode energy residual and signal quality index to dynamically adjust the observation noise. It combines the use duration and bandwidth as constraints to perform maximum likelihood estimation iterative solution, which solves the problems of existing methods ignoring the time-varying characteristics of inertia and lacking a dynamic correction mechanism, and achieves high-precision smooth tracking of the inertial time constant. Attached Figure Description

[0021] Figure 1 This is a flowchart illustrating an online inertial tracking method for a power system according to the present invention. Figure 2 This is a schematic diagram of the system structure of an online inertial tracking method for power systems according to the present invention; Figure 3 A structural block diagram of an exemplary electronic device provided for implementing embodiments of the present disclosure; Figure 4 This is the IEEE 39-node power distribution system used in the embodiments of the present invention; Figure 5 This is a schematic diagram of the dominant inertial mode recognition in this invention; Figure 6 This is a comparison curve before and after the correction of the inertial time constant of the present invention. Detailed Implementation

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

[0023] Example 1, Figure 1 This invention provides a method for online inertial tracking in a power system, comprising the following steps: S1, the effective time period of inertia is selected based on the rate of change of frequency and the energy coupling coefficient; In this embodiment, the step of selecting the effective inertial time period based on the frequency change rate and energy coupling coefficient specifically involves: S1.1, the specific method for collecting the power and frequency signals of the power generation device is as follows: collecting the active power signal at the connection bus of the power generation device. and frequency signals In this embodiment, the data sampling frequency is set to... This system can meet the analysis needs of electromechanical dynamic processes in power systems. The collected data is cached in the form of a sliding data window, with the window length set to... The overlap rate between adjacent data windows is 50%.

[0024] S1.2, the step of calculating the frequency change rate based on the frequency signal, and marking the disturbance trigger moment when its absolute value exceeds a preset threshold, specifically involves: The acquired frequency signal Numerical differentiation is performed to calculate the rate of change of frequency. This embodiment uses the central difference method to improve calculation accuracy.

[0025] in, At the current sampling time, the th calculation of the rate of change of frequency is performed. A discrete time point, The frequency of the next time step represents the frequency of the current time step. The frequency measurement value at the next sampling interval, The frequency of the previous moment represents the frequency of the current moment. The frequency measurement value at the previous sampling interval, For discrete-time index, it represents the first index in the sampled sequence. One sampling point, The sampling interval is denoted as .

[0026] The preset threshold is a frequency change rate threshold; a frequency change rate threshold is set. At a certain moment satisfy When the time is reached, this moment is marked as the disturbance trigger moment. The selection of this threshold is based on the statistical characteristics of frequency fluctuations during normal operation of the power system: the frequency change rate under normal fluctuations is usually less than 0.005 Hz / s, so 0.01 Hz / s can effectively identify disturbance events with potential inertial response characteristics, while avoiding false triggering caused by measurement noise.

[0027] S1.3, the process of forming a candidate response interval by extending forward and backward from the disturbance trigger moment is specifically as follows: The disturbance trigger time marked by step S1.2 Extend forward from the center Seconds, extending backwards Seconds, forming a candidate response range .

[0028] In this embodiment, based on the typical timescale of electromechanical oscillations in a power system (typically within the range of 0.5 to 5 seconds), the following parameters are set: Second, Extending the time forward by 2 seconds can capture the steady-state operating state before the disturbance occurs, serving as a benchmark for response analysis; extending the time backward by 5 seconds can fully cover the main dynamic process of system frequency recovery after the disturbance, ensuring that inertial response information is not truncated.

[0029] S1.4, the calculation of the energy coupling coefficient between the change in active power and the rate of change of frequency within the candidate response interval, and the selection of the effective inertial time period based on this coefficient, specifically involves: Within the candidate response interval, it is further determined whether that time period truly contains valid inertial response information. To this end, an energy coupling coefficient is introduced. This is used to measure the degree of energy coupling between changes in active power and frequency. The energy coupling coefficient is... It can be represented as:

[0030]

[0031] in, For active power fluctuation components, For frequency fluctuation components, This is the active power signal. It is a frequency signal. and These are the average values ​​of active power and frequency within 1 second before the start of the candidate response interval, representing the steady-state baseline value before the disturbance.

[0032] Furthermore, an energy coupling coefficient threshold is set. If the energy coupling coefficient is greater than the preset energy coupling coefficient threshold, the current candidate response interval is determined to be an inertial effective time period; otherwise, the data of that time period is discarded and the system continues to wait for the next trigger moment.

[0033] Furthermore, based on the disturbance amplitude and duration within the effective inertial time period, the usage duration of this time period in subsequent steps is determined. The disturbance amplitude can be characterized by the maximum offset of the frequency signal or active power change, and the duration is the length of time the frequency change rate remains above a threshold. When the disturbance amplitude is large or the duration is long, the usage duration of the effective inertial time period can be appropriately extended to ensure that the subsequent time-frequency analysis and model identification processes can fully utilize the inertial response information; when the disturbance amplitude is small or the duration is short, the usage duration of the effective inertial time period can be shortened accordingly to reduce the impact of noise signals on subsequent calculation results.

[0034] It should be noted that the energy coupling coefficient is designed based on the physical essence of the generator rotor motion equation, and can quantify the degree of energy coupling between the change in active power and the rate of frequency change. Only when the energy coupling coefficient is greater than the preset energy coupling coefficient threshold is it determined to be an effective inertial period, thereby eliminating non-inertial disturbances and strong noise periods at the source. Compared with the unfiltered full-time processing method, this step allows subsequent identification to focus only on data segments that truly contain inertial response, avoiding invalid calculations, and reducing the inertial estimation error by about 15% (as shown in Table 1).

[0035] Table 1 Comparison of estimation errors with and without time period filtering

[0036] S2 performs synchronous squeezing wavelet transform on the signal within the effective time period, identifies the dominant low-frequency energy ridge and constructs an adaptive passband, and obtains the denoised signal through mask filtering and inverse transform reconstruction. In this embodiment, the step of performing synchronous squeezing wavelet transform on the signal within the effective time period, identifying the dominant low-frequency energy ridge and constructing an adaptive passband, and obtaining the denoised signal through mask filtering and inverse transform reconstruction is specifically as follows: S2.1, the construction of the adaptive passband includes: S2.1.1, the synchronous squeezing wavelet transform of the power and frequency signals to obtain the time-frequency energy distribution specifically involves: For signal (can represent active power signal) or frequency signal Perform continuous wavelet transform:

[0037] in, These are wavelet transform coefficients, representing the signal at different scales. Translation The wavelet transform result at the location, This is the scale factor, corresponding to the reciprocal of the frequency. The translation factor corresponds to time. For the mother wavelet function, Indicates complex conjugation.

[0038] In this embodiment, the Morlet wavelet is selected as the mother wavelet, and the parameter settings are shown in Table 2. Its expression is:

[0039] in, For the center frequency, Represents a time variable.

[0040] Table 2. Parameter settings for synchronous squeezing wavelet transform

[0041] Furthermore, after continuous wavelet transform, the instantaneous frequency is calculated. :

[0042] It should be noted that the core idea of ​​the synchronous compression is to use wavelet coefficients... The energy concentration of the time-frequency representation is improved by redistributing the coefficients near the same frequency according to the instantaneous frequency.

[0043] Results of synchronous squeeze wavelet transform for:

[0044] in, Angular frequency represents the discretized frequency point. For frequency resolution, As a discretization scale, The scale interval is used. Through synchronous compression, the energy that was originally ambiguous on the time-frequency plane is redistributed to more precise frequency positions, resulting in the time-frequency energy distribution. , where ω is the angular frequency and t is time.

[0045] S2.1.2, the step of identifying the dominant low-frequency energy ridge based on the energy spectral density and constructing an adaptive passband centered on the energy ridge specifically includes: Define energy spectral density :

[0046] The dominant low-frequency energy ridge Ridge extraction refers to the frequency trajectory at each time t that results in a local maximum in the energy spectral density. This can be achieved by solving the following optimization problem:

[0047] The first term maximizes the total energy along the ridge line, ensuring that the ridge line is located in an energy-concentrated region; the second term penalizes violent fluctuations in the ridge line, ensuring the smoothness of the ridge line. The regularization parameter balances the objectives of maximizing energy and ensuring smoothness. In this embodiment, the regularization parameter is set as follows: This value is determined based on the statistical characteristics of typical electromechanical oscillation signals, and can maintain the smoothness of the ridge line without losing important frequency variation details.

[0048] In the actual calculations, a dynamic programming algorithm is used to solve the problem. First, local maxima are found at each time step along the frequency direction. Then, dynamic programming connects these points to form a continuous ridge trajectory. For the active power signal and the frequency signal, the corresponding energy ridges are extracted respectively. and Typically, due to the strong correlation between power and frequency response, the two ridge lines largely coincide. To improve robustness, this embodiment uses a weighted average of the two as the final dominant low-frequency energy ridge line. .

[0049] Furthermore, based on the identified dominant low-frequency energy ridges Construct an adaptive passband centered on [the target]. :

[0050] in, The ridge half-width is defined as the frequency offset at which the energy around the ridge drops to half of the peak energy, i.e., the half-power bandwidth. This is the bandwidth expansion factor, used to appropriately widen the passband range to ensure that complete inertial response information is preserved.

[0051] S2.2, the noise-reduced signal obtained through mask filtering and inverse transform reconstruction specifically includes: S2.2.1, using the adaptive passband constructed in step S2.1.2 The time-frequency coefficients obtained by synchronous squeezing wavelet transform Perform mask filtering. The filtering rule is: retain the time-frequency coefficients within the passband, and set the time-frequency coefficients outside the passband to zero.

[0052] For the active power signal and the frequency signal, perform the above mask filtering operation respectively to obtain the filtered time-frequency coefficients. and .

[0053] S2.2.2, the step of performing synchronous squeezing wavelet inverse transform on the filtered time-frequency coefficients to reconstruct the denoised power fluctuation component and frequency fluctuation component is as follows: For the filtered time-frequency coefficients Perform a synchronous squeezing wavelet inverse transform to reconstruct the time-domain signal. The formula for the synchronous squeezing wavelet inverse transform is:

[0054] in, To reconstruct the wavelet, the dual wavelet of the mother wavelet is usually taken. This indicates taking the real part.

[0055] In practical calculations, the inverse transform described above can be achieved using numerical integration. First, each frequency component is integrated in the time direction, and then the results are summed in the frequency direction to obtain the reconstructed time-domain signal.

[0056] Furthermore, the filtered power time-frequency coefficients Perform an inverse transform to obtain the power fluctuation component after noise reduction. For the filtered frequency time coefficients Perform an inverse transform to obtain the noise-reduced frequency fluctuation components. .

[0057] It should be noted that, compared to traditional wavelet transform or Fourier transform, synchronous squeezing wavelet transform can "squeeze" fuzzy energy on the time-frequency plane onto a more precise frequency ridge, greatly improving the energy concentration of the time-frequency representation. This makes the extraction of the dominant low-frequency energy ridge more accurate, laying a solid foundation for the subsequent construction of an adaptive passband. An adaptive passband is constructed centered on the energy ridge, and mask filtering is used to remove out-of-band components, achieving precise filtering driven by the physical characteristics of the signal. Compared to low-pass filtering with a fixed cutoff frequency, the passband of this invention dynamically adjusts with the frequency of the signal itself, rather than being a fixed band, fully preserving the dominant frequency components of the inertial response, avoiding signal distortion, and specifically processing inertial modes to effectively separate non-inertial components.

[0058] S3, construct a difference equation based on the noise-reduced signal, and recursively identify the pulse transfer function according to the order of the adaptive passband width dynamic constraint equation; In this embodiment, the step of constructing a difference equation based on the noise-reduced signal and recursively identifying the pulse transfer function according to the order of the adaptive passband width dynamic constraint equation is as follows: S3.1, the difference equation constructed based on the denoised power fluctuation component and frequency fluctuation component is as follows: The power fluctuation component after noise reduction As input signal Frequency fluctuation component As output signal Construct a difference equation describing the active-frequency dynamic characteristics of the system:

[0059] in, For discrete-time indexing, The output at the current moment, i.e. , The input at the current moment, i.e. , This outputs the order, specifically the order of the autoregressive component. The input order is the order of the exogenous input component. and The model parameters to be identified. The model residuals represent unmodeled dynamics and noise.

[0060] Furthermore, a shift operator is introduced. (satisfy The above difference equation can be expressed in transfer function form:

[0061] in:

[0062]

[0063] The corresponding pulse transfer function is:

[0064] S3.2, the step of limiting the search range of the order of the difference equation based on the width of the adaptive passband specifically involves: Constructed adaptive passband Its bandwidth is:

[0065] Determine the model order search range based on the bandwidth:

[0066]

[0067] in, For reference bandwidth, and This is an empirical coefficient. For the maximum allowed order, This is for rounding down.

[0068] S3.3, within the search range, the information content criterion is used to determine the order of the equation at the current time step, specifically as follows: Within the order range The inner loop iterates through different order combinations and calculates the information content criterion. This embodiment adopts the Akaike information content criterion:

[0069]

[0070] Where N is the number of data points used in the current time step. For the model residual variance, The number of model parameters serves as a penalty for model complexity. For the summation index, For the first Each residual value.

[0071] Choose the order combination with the smallest AIC value as the optimal order at the current time step:

[0072] in, For the optimal output order, This is the optimal input order.

[0073] In this embodiment, in addition to the AIC criterion, the Bayesian Information Criterion (BIC) can also be used as an alternative:

[0074] This embodiment uses AIC by default, while retaining BIC as an optional configuration.

[0075] S3.4 Based on the determined order of the equations, the equation parameters are updated through online learning to obtain the pulse transfer function at the current time step.

[0076] After determining the optimal order, the model parameters are updated online using the recursive least squares method.

[0077] Rewrite the difference model in linear regression form:

[0078]

[0079]

[0080] in, For the regression vector, For the parameter vector to be identified, the superscript It is the transpose operator. This is a discrete-time index.

[0081] It should be noted that the update formula for the recursive least squares method is common knowledge and will not be elaborated here. Through recursive updates, real-time online identification of model parameters can be achieved.

[0082] Through the above recursive update, the updated parameter vector can be obtained at each sampling time k. Then, the pulse transfer function for the current time step is constructed:

[0083] in, This is the estimated value of the numerator polynomial at the current time step. This is the estimated value of the denominator polynomial at the current time step.

[0084] In this embodiment, to reduce computational load, it is unnecessary to perform an order search at every sampling time. Each order search... The process is executed once per second (i.e., every 100 sampling points), with intermediate time steps using the most recently determined order. This ensures that the model can adapt to the slow time-varying characteristics of the system while avoiding the computational burden caused by frequent order searches.

[0085] It should be noted that this invention introduces the physical features (bandwidth) extracted in step S2 into the model identification stage, establishing an intrinsic link between the physical characteristics of the signal and the complexity of the mathematical model. A wider bandwidth indicates a more complex system dynamic and requires a higher model order; conversely, a narrower bandwidth allows for description using a lower order. Compared to methods that only use the AIC criterion to search within a fixed range [0,10], this step dynamically limits the order range through bandwidth, shrinking the order range from the unconstrained 1~20 to a finite interval with clear physical meaning. This reduces the order search space by approximately 60% and the computational load by approximately 40%.

[0086] S4. Perform continuous domain transformation on the pulse transfer function, identify the dominant inertial mode and calculate the estimated value of the inertial time constant through the damping ratio-frequency composite criterion; In this embodiment, the continuous-domain transformation of the pulse transfer function, the identification of the dominant inertial mode using the damping ratio-frequency composite criterion, and the calculation of the estimated inertial time constant specifically include: S4.1, the continuous domain transformation includes: S4.1.1, obtain the adaptive passband boundary frequency of the current time step, and the center frequency corresponding to the dominant low-frequency energy ridge, specifically: The frequency range of the adaptive passband is:

[0087] in, This is the lower boundary frequency of the passband. The upper boundary frequency of the passband. It is a time variable.

[0088] The dominant low-frequency energy ridge is extracted based on the time-frequency energy distribution, and its corresponding center frequency is denoted as:

[0089] To reduce the impact of instantaneous fluctuations, this embodiment averages the frequency parameters within the previous short time window, obtaining the following results: , and .

[0090] S4.1.2, due to the identified pulse transfer function The discrete transfer function is described in the discrete domain (z-domain), while the physical properties of the system (such as inertia) are usually more intuitively expressed in the continuous domain (s-domain). Therefore, it is necessary to convert the discrete transfer function into a continuous transfer function.

[0091] This embodiment uses the bilinear transformation method to realize the mapping from the discrete domain to the continuous domain. The basic relationship is as follows:

[0092] in, For the Laplace transform complex frequency domain variables, The sampling period is This is the shift operator.

[0093] Since the bilinear transformation introduces frequency nonlinear distortion, this embodiment uses the center frequency obtained in step S4.1.1 for pre-distortion correction. First, the center frequency is converted into a continuous angular frequency:

[0094] The corrected angular frequency is obtained based on the frequency mapping relationship of the bilinear transform:

[0095] Therefore, the frequency correction coefficient is obtained:

[0096] S4.2, Based on the corrected frequency parameters, a continuous-domain transformation is performed on the pulse transfer function model to obtain the distortion-corrected continuous-domain transfer function, specifically: After obtaining the corrected frequency parameters, the pulse transfer function obtained in step S3 is... Perform continuous domain transformation.

[0097] pulse transfer function variables in Replace it with the following according to the bilinear transformation relationship:

[0098] After algebraic simplification, the transfer function in the continuous field is obtained:

[0099] in, Let be the order of the denominator polynomial. Let the order of the numerator polynomial be denoted by . These are the coefficients of the denominator polynomial and the numerator polynomial, respectively. This represents the complex frequency domain variable of the Laplace transform.

[0100] S4.3, the identification of the dominant inertial mode and the calculation of the estimated inertial time constant specifically includes: S4.3.1 decomposes the continuous domain transfer function into independent dynamic modal components characterized by different poles, specifically as follows: Transfer function in the continuous domain The transfer function can be decomposed into a superposition of several independent dynamic modal components. For a linear time-invariant system, the transfer function can be expressed by partial fractional decomposition as follows:

[0101] in, For the current time step index, For the summation index, Indicates the first The system poles can be real numbers or conjugate complex numbers. For the first The residues at each pole represent the degree of participation of that mode. For directly transitive terms, they exist when the order of the numerator equals the order of the denominator, indicated by the subscript. The index identifier identifies different poles and their corresponding residues.

[0102] In this embodiment, a numerical method is used to solve the problem. The poles are found. For higher-order systems (n≤10), the roots of the denominator polynomial can be directly calculated; for even higher-order systems, numerically stable algorithms such as QR decomposition can be used. Solving for the pole set yields the solution. Each pole represents an independent dynamic mode component.

[0103] For each pair of conjugate complex poles ,in, The attenuation coefficient determines the decay rate of the mode. Let be the damped oscillation frequency. Calculate the natural oscillation frequency. Damping ratio :

[0104]

[0105] S4.3.2, based on a composite criterion of damping ratio and natural oscillation frequency, the dominant inertial mode conforming to electromechanical oscillation characteristics is identified from the independent dynamic modal components, and the decay time constant corresponding to the mode is extracted, specifically: In this embodiment, the preset composite criterion based on damping ratio and natural oscillation frequency is specifically as follows:

[0106] in, and This represents the damping ratio threshold. Based on the typical characteristics of electromechanical oscillations in power systems, the inertial response modes typically exhibit weakly or moderately damped oscillations, with damping ratios generally between 0.05 and 0.3. If multiple modes satisfy the condition, the mode with the highest oscillation energy is selected as the dominant inertial mode.

[0107] After identifying the dominant inertial mode, the attenuation coefficient corresponding to that mode is extracted. And calculate the corresponding decay time constant:

[0108] S4.3.3, establish the analytical mapping relationship between the decay time constant and the equivalent moment of inertia, and calculate the estimated value of the current time step inertial time constant, specifically as follows: According to the small-signal analysis theory of power systems, there is an approximate mapping relationship between the system's equivalent inertial time constant and the decay time constant of the dominant dynamic mode. Therefore, this embodiment calculates the estimated value of the inertial time constant at the current time step using the following relationship. :

[0109] in, The mapping coefficients can be obtained through historical system operating data or offline simulation calibration. In this embodiment, the values ​​are taken based on the calibration results of a typical system. , This is the frequency correction coefficient, used to eliminate the distortion caused by the bilinear transformation.

[0110] It should be noted that a single frequency criterion can easily misclassify control modes or noise modes as inertial modes, and a single damping ratio criterion is also difficult to distinguish between different modes with similar frequencies. The dual constraint mechanism of this invention simultaneously filters modes from both frequency and damping dimensions, significantly improving identification reliability. The frequency filtering range directly adopts the adaptive passband from step S2, rather than a fixed passband, allowing the criterion to adapt to changes in the inertial response frequency under different operating conditions. When a change in system operation causes a shift in the oscillation frequency, the passband is adjusted accordingly, ensuring that the identified mode is always the dominant inertial mode under the current operating condition.

[0111] S5. Construct a dynamic Bayesian network and solve iteratively through maximum likelihood estimation to output the correction value of the inertial time constant at the current moment.

[0112] In this embodiment, the construction of the dynamic Bayesian network, through iterative solution using maximum likelihood estimation, and the output of the current inertial time constant correction value, specifically includes: S5.1 Calculate the instantaneous offset of the center frequency of the energy ridge between adjacent time steps, set a state discrimination threshold, if the instantaneous offset is less than or equal to the state discrimination threshold, the system is determined to be in steady state; otherwise, if the instantaneous offset is greater than the state discrimination threshold, the system is determined to be in transient state.

[0113] S5.2, the construction of Bayesian networks with different prior distributions specifically includes: When the system is in steady state, the inertial time constant changes slowly and exhibits strong temporal correlation. In this case, a Gaussian prior distribution centered on the previous inertial correction value is constructed. When the system is in transient state, the inertial time constant may change rapidly, and it is not advisable to rely too much on historical information. In this case, a uniform prior distribution or a wide-variance Gaussian distribution is constructed to allow for a larger range of variation. The observation nodes of the dynamic Bayesian network correspond to the estimated current time-step inertial time constant obtained in step S4. .

[0114] S5.3, Obtain the dominant mode energy residual and signal quality index at the current time step, and dynamically adjust the likelihood probability distribution parameters of the Bayesian network observation nodes, specifically: This embodiment introduces two metrics to quantify the estimation quality of the current time step: the dominant mode energy residual. and signal quality indicators The dominant mode energy residual is defined as the ratio of the energy removed after mask filtering in step S2-4 within the frequency band of the dominant inertial mode to the original energy. The smaller the ratio, the better the signal quality and the purer the inertial response components; the larger the ratio, the more noise or interference components there are, and the lower the reliability of the estimation results. The signal quality index comprehensively evaluates the signal-to-noise ratio and data integrity, and its value ranges from [0,1]. The larger the value, the better the quality.

[0115] Furthermore, based on the aforementioned dominant mode energy residuals and signal quality indices, the likelihood probability distribution parameters of the Bayesian network observation nodes are dynamically adjusted. It is assumed that the observation error follows a Gaussian distribution:

[0116]

[0117] in, This represents the true value of the current time step inertial time constant. For reference energy residuals, this embodiment uses 0.2. To observe the noise variance, The standard deviation of the basic observation noise is taken as 0.1s in this embodiment. and For adjustment coefficients, This is the index for the current time step.

[0118] S5.4, using the usage duration and adaptive bandwidth as constraints, construct the confidence radius:

[0119] in, The confidence radius at the current time step. Based on the confidence radius, For usage duration, For reference duration, The current time step bandwidth, This is the reference bandwidth.

[0120] Furthermore, based on the above parameters, the posterior mean (i.e., the maximum likelihood estimate) of the dynamic Bayesian network can be expressed as:

[0121] in, This is the inertial correction value from the previous moment. The prior variance is 0.01 in steady state or 0.09 in transient state. To observe the noise variance.

[0122] Furthermore, considering physical characteristic constraints, the results are truncated:

[0123] Arrange the inertial time constant correction values ​​output at each time step in chronological order to obtain the real-time evolution trajectory of the equivalent inertia of the power system.

[0124] Example 2, Figure 2 The present invention provides a system for an online inertial tracking method for a power system, comprising: The data filtering module is used to filter effective inertial time periods based on the rate of change of frequency and the energy coupling coefficient. The time-frequency processing module is used to perform synchronous squeezing wavelet transform on the signal within the effective time period, identify the dominant low-frequency energy ridge and construct an adaptive passband, and obtain the noise-reduced signal through mask filtering and inverse transform reconstruction. The model identification module is used to construct a difference equation based on the noise-reduced signal and recursively identify the pulse transfer function according to the order of the adaptive passband width dynamic constraint equation. The modal analysis module is used to perform continuous domain transformation on the pulse transfer function, identify the dominant inertial mode through the damping ratio-frequency composite criterion, and calculate the estimated value of the inertial time constant. The result correction module is used to construct a dynamic Bayesian network, solve iteratively through maximum likelihood estimation, and output the correction value of the inertial time constant at the current moment.

[0125] Example 3, an electronic device, such as Figure 3 As shown, the electronic device includes a memory and a processor: the memory is used to store a program; the processor is used to execute the program to implement any of the embodiments in Example 1.

[0126] Example 4: To verify the effectiveness of the method of the present invention, a simulation model of the IEEE 10-machine 39-node standard test system was built in PSCAD / EMTDC. The topology is as follows: Figure 4 As shown, the system comprises 10 synchronous generators, 39 buses, and 12 transformers. All generators are equipped with IEEE standard excitation systems, speed governors, and power system stabilizers, enabling them to accurately reflect the dynamic characteristics of the actual power system.

[0127] Figure 4 In the diagram, ① to ⑩ represent 10 generator nodes within the system, which are the electrical connection points for the power generation units; numbers 30 to 39 are dedicated access buses for generators, corresponding one-to-one with generators ① to ⑩; nodes 1 to 15 and 17 to 29 are the tie buses and load buses in the system, connecting various transformers, loads, and transmission lines; number 16 is the test bus for applying load step disturbance in this embodiment.

[0128] This embodiment simulates a high-proportion renewable energy scenario. Three synchronous generators (buses 30, 32, and 37) in the system are replaced with doubly-fed wind farms of equal capacity, and virtual inertial control is configured. The remaining seven power generation nodes are retained as synchronous generators, achieving a renewable energy penetration rate of 30%. The system's base capacity is 100 MVA, and the rated frequency is 60 Hz.

[0129] A load step disturbance (increasing the load by 10% for 0.1s) was set at test bus 16. Active power and frequency data were collected at each generator bus, with the sampling frequency set to 100Hz. Inertial estimation was performed using the following three methods: Method A (the method of this invention): adopts synchronous squeezing wavelet transform for noise reduction + adaptive order constraint recursive identification + damping ratio-frequency composite criterion + dynamic Bayesian network correction.

[0130] Method B (State-space model method): A state-space model is established based on noise-like data. A first-order state-space model is obtained through model optimization and order reduction, and then the inertial time constant is solved.

[0131] Method C (ARMAX-LSTM method): Based on the recursive least squares method, the discrete transfer function of the ARMAX model is identified, and after being converted into a continuous transfer function by bilinear transformation, the inertial time constant is extracted. Then, the residual is corrected by using an LSTM neural network to evaluate the estimation results.

[0132] The comparison results are shown in Table 3: Table 3 Comparison of inertial estimation errors of different methods in the IEEE 39-node system

[0133] As shown in Table 3, under the IEEE 39-node standard test system, the error indicators of the method of this invention are significantly better than those of existing methods. The mean absolute error is reduced by 65% ​​compared to the state-space model method and by 55% compared to the ARMAX-LSTM method.

[0134] Since the electronic device described in this embodiment is the device used to implement the method in Embodiment 1 of the present invention, those skilled in the art can understand the specific implementation method and various variations of the electronic device in this embodiment based on the method described in Embodiment 1 of this application. Therefore, how the electronic device implements the method in the embodiments of this application will not be described in detail here. Any device used by those skilled in the art to implement the method in the embodiments of this application falls within the scope of protection of this application.

[0135] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0136] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented, in whole or in part, in the form of a computer program product.

[0137] Those skilled in the art will recognize that the modules and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0138] In addition, the functional modules in the various embodiments of this application can be integrated into one processing module, or each module can exist physically separately, or two or more modules can be integrated into one module.

[0139] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

[0140] In conclusion, the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for online inertial tracking in a power system, characterized in that, Includes the following steps: The effective time period of inertia is selected based on the rate of change of frequency and the energy coupling coefficient; Synchronous squeezing wavelet transform is performed on the signal within the effective time period to identify the dominant low-frequency energy ridge and construct an adaptive passband. The denoised signal is then obtained by mask filtering and inverse transform reconstruction. A difference equation is constructed based on the noise-reduced signal, and the pulse transfer function is recursively identified according to the order of the adaptive passband width dynamic constraint equation. The pulse transfer function is transformed into the continuous domain, and the dominant inertial mode is identified by the damping ratio-frequency composite criterion, and the estimated value of the inertial time constant is calculated. A dynamic Bayesian network is constructed, and the problem is solved iteratively through maximum likelihood estimation to output the correction value of the inertial time constant at the current moment.

2. The power system inertial online tracking method according to claim 1, characterized in that, The specific effective time period for screening inertia is as follows: Collect power and frequency signals from the power generation device; The frequency change rate is calculated based on the frequency signal, and when its absolute value exceeds a preset threshold, it is marked as the disturbance trigger moment. Candidate response intervals are formed by extending forward and backward from the moment the disturbance is triggered. Calculate the energy coupling coefficient between the change in active power and the rate of change of frequency within the candidate response interval, and select the effective inertial time period accordingly; Based on the disturbance amplitude and duration within the effective inertial time period, determine the duration of this time period to be used in subsequent steps.

3. The power system inertial online tracking method according to claim 2, characterized in that, The construction of the adaptive passband includes: The power and frequency signals are synchronously squeezed wavelet transformed to obtain the time-frequency energy distribution. The dominant low-frequency energy ridge is identified based on the energy spectral density, and an adaptive passband is constructed centered on the energy ridge.

4. The power system inertial online tracking method according to claim 3, characterized in that, The process of obtaining the denoised signal through mask filtering and inverse transform reconstruction is as follows: The time-frequency energy distribution is masked using the adaptive passband to retain the time-frequency coefficients within the passband while eliminating out-of-band noise and non-inertial mode components. The filtered time-frequency coefficients are subjected to synchronous squeezing wavelet inverse transform to reconstruct the denoised power fluctuation component and frequency fluctuation component.

5. The power system inertial online tracking method according to claim 4, characterized in that, The recursive identification pulse transfer function is specifically as follows: Based on the power fluctuation component and frequency fluctuation component after noise reduction, a difference equation is constructed; The search range for the order of the difference equation is limited based on the width of the adaptive passband. Within the search range, the information content criterion is used to determine the order of the equation at the current time step; Based on the determined order of the equations, the pulse transfer function at the current time step is obtained by updating the equation parameters through online learning.

6. The power system inertial online tracking method according to claim 5, characterized in that, The continuous domain transformation includes: Obtain the adaptive passband boundary frequency of the current time step, and the center frequency corresponding to the dominant low-frequency energy ridge. According to the preset sampling period, pre-distortion correction is performed on the boundary frequency and the center frequency respectively; Based on the corrected frequency parameters, a continuous domain transformation is performed on the pulse transfer function model to obtain the distortion-corrected continuous domain transfer function.

7. The power system inertial online tracking method according to claim 6, characterized in that, The identification of the dominant inertial mode and the calculation of the estimated inertial time constant specifically involve: The continuous domain transfer function is decomposed into independent dynamic modal components characterized by different poles; A composite criterion based on damping ratio and natural oscillation frequency is used to identify the dominant inertial mode that conforms to electromechanical oscillation characteristics from the independent dynamic modal components, and the decay time constant corresponding to the mode is extracted. An analytical mapping relationship between the decay time constant and the equivalent moment of inertia is established, and the estimated value of the inertial time constant at the current time step is calculated.

8. The power system inertial online tracking method according to claim 7, characterized in that, The construction of the dynamic Bayesian network, through iterative solution using maximum likelihood estimation, outputs the correction value of the inertial time constant at the current moment, specifically as follows: Based on the instantaneous shift of the center frequency of the energy ridge between adjacent time steps, determine whether the current system is in a steady state or a transient state, and construct Bayesian networks with different prior distributions accordingly. Obtain the dominant mode energy residual and signal quality index at the current time step, and dynamically adjust the likelihood probability distribution parameters of the Bayesian network observation nodes; Using the usage duration and adaptive bandwidth as constraints, the system iteratively solves the problem through maximum likelihood estimation and outputs the corrected value of the inertial time constant at the current moment.

9. A system using the power system inertial online tracking method as described in any one of claims 1-8, characterized in that, include: The data filtering module is used to filter effective inertial time periods based on the rate of change of frequency and the energy coupling coefficient. The time-frequency processing module is used to perform synchronous squeezing wavelet transform on the signal within the effective time period, identify the dominant low-frequency energy ridge and construct an adaptive passband, and obtain the noise-reduced signal through mask filtering and inverse transform reconstruction. The model identification module is used to construct a difference equation based on the noise-reduced signal and recursively identify the pulse transfer function according to the order of the adaptive passband width dynamic constraint equation. The modal analysis module is used to perform continuous domain transformation on the pulse transfer function, identify the dominant inertial mode through the damping ratio-frequency composite criterion, and calculate the estimated value of the inertial time constant. The result correction module is used to construct a dynamic Bayesian network, solve iteratively through maximum likelihood estimation, and output the correction value of the inertial time constant at the current moment.

10. An electronic device, characterized in that, The electronic device includes a memory and a processor: The memory is used to store programs; The processor is used to execute the program to implement the various steps of the power system inertial online tracking method as described in any one of claims 1-8.

Images