A power system inertia estimation method and device

Abstract

The application discloses a power system inertia estimation method and device, and solves the technical problem that the existing power system inertia estimation method based on TKEO fixed threshold triggering results in low accuracy of power system inertia estimation. The method comprises the following steps: collecting synchronous phasor measurement unit, PMU active power and steady-state active power data before disturbance; obtaining two normalized TKEO energies through double-channel TKEO energy feature extraction; generating double-channel coherent energy through causal time delay alignment and coherent fusion; completing IS divergence sequential detection and double-criterion decision of shared LTA relying on the energy; and outputting a disturbance event time, an event confidence weight and a sample-by-sample positive CUSUM increment, so as to complete weighted ARMAX identification and extract a power system inertia constant estimation value.

Description

Technical Field

[0001] This invention relates to the field of power system parameter evaluation technology, and in particular to a method and apparatus for estimating the inertia of a power system. Background Technology

[0002] With the large-scale replacement of traditional synchronous generators by inverter-type power sources such as wind power and photovoltaics, the equivalent inertia level of the system continues to decline. Inertia is the first line of defense for the system against frequency abrupt changes—the inertia response stage after a disturbance occurs (<10s). The system RoCoF (Rate of Change of Frequency) and the inertia constant H satisfy the quantitative relationship described by the oscillation equation.

[0003] When inertia is insufficient, RoCoF exceeds the protection relay setting threshold, triggering a cascading trip. Traditional inertia measurement relies on large disturbances of known magnitude (such as generator tripping), but these events occur infrequently in modern power grids and cannot support continuous online monitoring. Therefore, how to automatically detect disturbance event windows suitable for inertia estimation using electrical quantity measurements under normal operating conditions without relying on large disturbances has become a key upstream technology for ensuring frequency security in low-inertia power grids.

[0004] Existing power system inertia estimation methods based on fixed TKEO threshold triggering rely on the core principle of offline calibration of a fixed TKEO energy threshold. Disturbance event identification is achieved by comparing the instantaneous TKEO energy at a single sampling point with the threshold. Once a disturbance is identified, the corresponding data is extracted and input into an ARMAX model for parameter identification, ultimately deriving the system inertia constant. However, this fixed threshold only applies to a single grid operating condition. In reality, frequency baselines shift significantly during grid condition transitions, with frequency fluctuations reaching a standard deviation of 0.02Hz during peak wind power generation periods and only 0.002Hz during light-load nighttime conditions—a difference of an order of magnitude. Under high-fluctuation conditions, the peak TKEO energy generated by normal grid operation noise easily approaches the preset threshold, triggering false triggers of non-genuine active-frequency disturbances. Under low-fluctuation conditions, effective disturbances such as megawatt-level load switching, which have identifiable value, have TKEO energy peaks far below the fixed threshold, resulting in missed detections of effective disturbance events. Meanwhile, this method relies solely on single-point sampling values ​​to complete binary decisions, failing to accumulate weak temporal disturbance characteristics from continuous sampling points. A large number of small, effective disturbance signals suitable for inertia identification are directly filtered out and discarded. False triggers send pure noise data lacking electromechanical dynamic information into the ARMAX identification stage, while missed detections and weak disturbance filtering significantly reduce high signal-to-noise ratio effective samples. Both situations lead to significant deviations in the ARMAX model parameter identification results. After order reduction of the discrete transfer function and mapping calculation of the synchronous generator oscillation equation coefficients, the solved inertia constant will have a large deviation, ultimately resulting in low accuracy in power system inertia estimation. Summary of the Invention

[0005] This invention provides a method and apparatus for estimating the inertia of a power system, which solves the technical problem that the existing power system inertia estimation method based on TKEO fixed threshold triggering has low accuracy.

[0006] The first aspect of this invention provides a method for estimating the inertia of a power system, comprising:

[0007] Acquire data from the synchronous phasor measurement unit, PMU active power data, and steady-state active power before disturbance;

[0008] Based on the synchronous phasor measurement unit data, the PMU active power data and the steady-state active power before the disturbance, dual-channel TKEO energy feature extraction is performed to obtain normalized frequency channel TKEO energy and normalized active channel TKEO energy.

[0009] The normalized frequency channel TKEO energy and the normalized active channel TKEO energy are fused by causal delay alignment to obtain dual-channel coherent energy.

[0010] The IS divergence sequential detection and dual-criteria joint decision are performed on the coherent energy of the dual channels using shared LTA, and the perturbation event time, event confidence weight and sample-by-sample positive part CUSUM increment are output.

[0011] Based on the event time, the event confidence weight, and the per-sample positive CUSUM increment, a weighted ARMAX model identification and inertia constant extraction are performed to obtain an estimated value of the power system inertia constant.

[0012] Optionally, the step of extracting dual-channel TKEO energy features based on the synchronous phasor measurement unit data, the PMU active power data, and the pre-disturbance steady-state active power to obtain normalized frequency channel TKEO energy and normalized active channel TKEO energy includes:

[0013] Based on the frequency data in the synchronous phasor measurement unit data and the rated frequency, calculate the frequency deviation;

[0014] The active power deviation is calculated based on the PMU active power data and the steady-state active power before the disturbance.

[0015] The frequency deviation and the active power deviation are smoothed and denoised respectively to obtain the smoothed frequency deviation and the smoothed active power deviation.

[0016] Based on the smoothed frequency deviation and the smoothed active power deviation, calculate the frequency channel TKEO energy and the active channel TKEO energy.

[0017] The frequency channel TKEO energy and the active channel TKEO energy are normalized to obtain normalized frequency channel TKEO energy and normalized active channel TKEO energy.

[0018] Optionally, the coherent energy fusion of the normalized frequency channel TKEO energy and the normalized active channel TKEO energy, causally time-delay aligned, to obtain dual-channel coherent energy includes:

[0019] The normalized active channel TKEO energy is delayed and aligned according to causal time delay to obtain the time-delay aligned active channel TKEO energy.

[0020] Calculate the geometric mean of the TKEO energy of the active channel after time delay alignment and the TKEO energy of the normalized frequency channel;

[0021] Based on the geometric mean, the coherent energy of the two channels is determined.

[0022] Optionally, the step of performing sequential IS divergence detection and dual-criteria joint decision on the coherent energy of the dual channels using shared LTA, and outputting the perturbation event time, event confidence weight, and per-sample positive part CUSUM increment, includes:

[0023] Based on the dual-channel coherent energy, the short-term window mean and the long-term window mean are calculated;

[0024] Calculate the STA / LTA ratio based on the short window mean and the long window mean;

[0025] Calculate the adaptive threshold based on the long-term window mean;

[0026] The IS divergence is calculated using the dual-channel coherent energy and the long-time window mean.

[0027] Based on the IS divergence, determine the IS-CUSUM cumulative amount;

[0028] Calculate the corrected CUSUM decision threshold;

[0029] Based on the STA / LTA ratio, the adaptive threshold, the IS-CUSUM cumulative amount, and the modified CUSUM decision threshold, a dual-criteria joint decision is performed to determine the timing of the disturbance event.

[0030] The event confidence weight is calculated using the IS-CUSUM cumulative amount and the modified CUSUM decision threshold.

[0031] Calculate the per-sample positive CUSUM increment based on the time of the disturbance event.

[0032] Optionally, the step of performing weighted ARMAX model identification and inertia constant extraction based on the event time, the event confidence weight, and the per-sample positive CUSUM increment to obtain an estimated value of the power system inertia constant includes:

[0033] Based on the event time, extract the window frequency deviation and window active power deviation within the effective period of the disturbance.

[0034] The event confidence weight and the per-sample positive part CUSUM increment are calculated by performing a Gaussian time nearest neighbor and information density weighted calculation to determine the sample weight;

[0035] Using the frequency deviation within the window and the active power deviation within the window as inputs, the ARMAX model parameters are identified by substituting the sample weights into the weighted least squares method.

[0036] A second-order discrete transfer function is constructed using the ARMAX model parameters, and then the second-order discrete transfer function is converted into a first-order continuous transfer function.

[0037] The coefficients of the first-order continuous transfer function are compared with the first-order theoretical transfer function corresponding to the swing equation of the synchronous generator to determine the coefficient mapping relationship.

[0038] The estimated value of the power system inertia constant is calculated based on the coefficient mapping relationship and the ARMAX model parameters.

[0039] A second aspect of the present invention provides a power system inertia estimation device, comprising:

[0040] The acquisition module is used to acquire data from the synchronous phasor measurement unit, PMU active power data, and steady-state active power before disturbance.

[0041] The feature extraction module is used to perform dual-channel TKEO energy feature extraction based on the synchronous phasor measurement unit data, the PMU active power data and the steady-state active power before the disturbance, to obtain the normalized frequency channel TKEO energy and the normalized active channel TKEO energy.

[0042] The fusion module is used to perform causal delay-aligned coherent energy fusion of the normalized frequency channel TKEO energy and the normalized active channel TKEO energy to obtain dual-channel coherent energy.

[0043] The decision module is used to perform sequential detection of IS divergence of shared LTA and joint decision of dual criteria on the coherent energy of the dual channels, and outputs the perturbation event time, event confidence weight and per-sample positive part CUSUM increment.

[0044] The inertia extraction module is used to perform weighted ARMAX model identification and inertia constant extraction based on the event time, the event confidence weight, and the per-sample positive part CUSUM increment, so as to obtain the estimated value of the power system inertia constant.

[0045] Optionally, the feature extraction module is specifically used for:

[0046] Calculate the frequency deviation based on the frequency data and the rated frequency;

[0047] The active power deviation is calculated based on the PMU active power data and the steady-state active power before the disturbance.

[0048] The frequency deviation and the active power deviation are smoothed and denoised respectively to obtain the smoothed frequency deviation and the smoothed active power deviation.

[0049] Based on the smoothed frequency deviation and the smoothed active power deviation, calculate the frequency channel TKEO energy and the active channel TKEO energy.

[0050] The frequency channel TKEO energy and the active channel TKEO energy are normalized to obtain normalized frequency channel TKEO energy and normalized active channel TKEO energy.

[0051] A third aspect of the present invention provides an electronic device, including a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor performs the steps of the power system inertia estimation method described above.

[0052] The fourth aspect of the present invention provides a computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed, implements the power system inertia estimation method as described above.

[0053] The fifth aspect of the present invention provides a computer program product comprising a computer program stored on a non-transitory computer-readable storage medium, the computer program comprising program instructions, wherein when the program instructions are executed by a computer, the computer performs the steps of the power system inertia estimation method described above.

[0054] As can be seen from the above technical solutions, the present invention has the following advantages:

[0055] The present invention provides a method for estimating the inertia of a power system, which acquires synchronous phasor measurement unit (TPMU) data, PMU active power data, and steady-state active power before disturbance; based on the synchronous phasor measurement unit data, PMU active power data, and steady-state active power before disturbance, dual-channel TKEO energy feature extraction is performed to obtain normalized frequency channel TKEO energy and normalized active channel TKEO energy; coherent energy fusion with causal time delay alignment is performed on the normalized frequency channel TKEO energy and normalized active channel TKEO energy to obtain dual-channel coherent energy; and the dual-channel... The coherent energy is subjected to sequential IS divergence detection and dual-criteria joint decision-making using shared LTA, outputting the disturbance event time, event confidence weight, and per-sample positive part CUSUM increment. Based on the event time, event confidence weight, and per-sample positive part CUSUM increment, weighted ARMAX model identification and inertia constant extraction are performed to obtain the estimated value of the power system inertia constant. Based on the above scheme, after collecting synchronous phasor measurement unit data, PMU active power data, and pre-disturbance steady-state active power, this invention sequentially performs dual-channel TKEO energy feature extraction and causal time delay alignment. The entire process, including dry energy fusion, sequential IS divergence detection and dual-criteria joint decision-making using shared LTA, weighted ARMAX model identification, and inertia constant extraction, relies on dual-channel temporal coherent fusion to replace the traditional single-channel fixed threshold single-point judgment logic. It fuses two TKEO energy streams through time-delay aligned geometric mean fusion, retaining only real disturbance signals conforming to electromechanical timing patterns, filtering out random noise from isolated channels, avoiding noise-triggered invalid identification data, and generating adaptive thresholds based on real-time updated long-term window baselines while accumulating weak temporal disturbance characteristics. This adapts to different load fluctuation conditions in the power grid. The system captures small load disturbances to expand the effective identification samples. The detection process simultaneously outputs the event confidence weight and the per-sample positive CUSUM increment, which are directly used for subsequent identification weighting. The sample weights are distinguished according to the perturbation confidence, and the interference of low signal-to-noise ratio data on ARMAX parameter solving is suppressed, thereby reducing model identification bias. The power system inertia constant is calculated by relying on accurate model parameters. The entire integrated detection and identification process eliminates the estimation errors caused by false triggering, missed detection of effective disturbances, and indiscriminate equal weighting of the existing fixed threshold TKEO scheme, effectively improving the accuracy of power system inertia estimation. Attached Figure Description

[0056] To more clearly illustrate the technical solutions in the embodiments of the present invention 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 the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0057] Figure 1This is a flowchart illustrating the steps of a power system inertia estimation method provided in Embodiment 1 of the present invention.

[0058] Figure 2 This is a structural block diagram of a power system inertia estimation device provided in Embodiment 2 of the present invention. Detailed Implementation

[0059] This invention provides a method and apparatus for estimating the inertia of a power system, which solves the technical problem that existing power system inertia estimation methods based on TKEO fixed threshold triggering have low accuracy in estimating the inertia of the power system.

[0060] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, 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, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention. It should be noted that in the optional embodiments of the present invention, the object information and other related data involved require the permission or consent of the object when the embodiments of the present invention are applied to specific products or technologies, and the collection, use, and processing of related data must comply with relevant laws, regulations, and standards. That is to say, if the embodiments of the present invention involve data related to the object, it needs to be obtained with the authorization and consent of the object, the authorization and consent of relevant departments, and in compliance with relevant laws, regulations, and standards. If personal information is involved in the embodiments, the acquisition of all personal information requires the consent of the individual. If sensitive information is involved, the separate consent of the information subject is required, and the embodiments also need to be implemented with the authorization and consent of the object.

[0061] PMU: Synchronous Phasor Measurement Unit, a device that synchronously acquires voltage, current, and frequency at the generator bus.

[0062] TKEO: Teager-Kaiser Energy Operator, a nonlinear difference operator;

[0063] STA / LTA: Short-Term Average / Long-Term Average, which is the ratio of the short-term average to the long-term average to characterize the degree of transient deviation of the signal;

[0064] CUSUM: Cumulative Sum, a recursive form of the sequential probability ratio test, accumulates the statistical deviation over time, and determines that a change point has occurred when it exceeds a preset threshold;

[0065] IS divergence: Itakura-Saito divergence, used for Bregman divergence of non-negative signals, has scale invariance. );

[0066] ARMAX: Auto-Regressive Moving Average with Xogenous Input, describes a system output as a linear combination of historical output, external input, and noise, and its MA component has built-in noise modeling capability;

[0067] WLS: Weighted Least Squares, a variant of least squares that assigns different weights to samples to reflect their credibility;

[0068] QCD: Quickest Change Detection, a sequential statistical decision framework that minimizes detection latency under the constraint of controlling the false alarm rate;

[0069] RoCoF: Rate of Change of Frequency, the first derivative of the system frequency with respect to time;

[0070] Inertia constant H: The ratio of the rotor kinetic energy to the rated capacity of a synchronous generator;

[0071] HMA: Hybrid Moving Average, a signal smoothing method that combines simple moving average and weighted moving average, has better delay and envelope suppression performance than single smoothing methods;

[0072] Causal delay, Channel TKEO energy peak to The time difference (number of sampling points) of the channel TKEO energy peak is determined by the time constant of the oscillation equation. The value is calculated as 2H / D, where H is the inertia constant of the synchronous generator and D is the damping coefficient of the synchronous generator. This is the electromechanical response time constant, i.e., the oscillation equation time constant;

[0073] Dual-channel coherent energy, the geometric mean of two normalized TKEO energies after causal time delay alignment.

[0074] Please see Figure 1 , Figure 1 The flowchart illustrates the steps of a power system inertia estimation method provided in Embodiment 1 of the present invention.

[0075] This invention provides a method for estimating the inertia of a power system, comprising:

[0076] Step 101: Obtain the synchronous phasor measurement unit data, PMU active power data, and steady-state active power before the disturbance.

[0077] The synchronous phasor measurement unit (TPMU) collects three-phase voltage, three-phase current, and frequency timing measurement data synchronously at a sampling rate of 30 times per second. All data have a unified timestamp, enabling synchronous recording of the electrical operating status of the power grid at the same moment.

[0078] The PMU active power data is the generator active power time-series data obtained by real-time calculation of the real-time three-phase voltage and three-phase current values ​​collected by the synchronous phasor measurement unit, which can reflect the small fluctuations in active power in real time.

[0079] The steady-state active power before the disturbance is the constant active power benchmark value during the period when the power grid does not generate active power disturbance and the frequency and active power remain stable. It serves as the reference standard for calculating the active power deviation.

[0080] It should be noted that the synchronous phasor measurement unit deployed at the generator bus position synchronously collects three-phase voltage, three-phase current, and frequency timing data at a sampling rate of 30 times / second. This data constitutes the synchronous phasor measurement unit data. The active power data of the PMU is obtained by real-time calculation of the voltage and current output by the synchronous phasor measurement unit. The constant active power value during the stable operation phase of the power grid before the disturbance is extracted as the steady-state active power before the disturbance.

[0081] Step 102: Based on the synchronous phasor measurement unit data, PMU active power data and steady-state active power before disturbance, perform dual-channel TKEO energy feature extraction to obtain the normalized frequency channel TKEO energy and the normalized active channel TKEO energy.

[0082] It should be noted that the frequency deviation and active power deviation are first solved by combining the rated frequency and the steady-state active power before the disturbance. The two sets of deviation signals are smoothed and denoised by using a hybrid moving average. Then, the two original instantaneous energies are calculated by using the TKEO operator. The numerical scaling is completed by using the long time window sliding peak normalization method. Finally, two sets of normalized energies are output.

[0083] Furthermore, step 102 may include the following sub-steps:

[0084] S21. Calculate the frequency deviation based on the frequency data in the synchronous phasor measurement unit data and the rated frequency;

[0085] S22. Calculate the active power deviation based on the PMU active power data and the steady-state active power before the disturbance;

[0086] S23. Perform smoothing and noise reduction on the frequency deviation and active power deviation respectively to obtain the smoothed frequency deviation and smoothed active power deviation.

[0087] S24. Based on the smoothed frequency deviation and smoothed active power deviation, calculate the TKEO energy of the frequency channel and the TKEO energy of the active channel.

[0088] S25. Normalize the frequency channel TKEO energy and the active channel TKEO energy to obtain the normalized frequency channel TKEO energy and the normalized active channel TKEO energy.

[0089] It should be noted that the generator bus three-phase voltage, three-phase current, and frequency data synchronously acquired by the PMU are processed by calculating the frequency deviation and active power deviation, and then HMA smoothing is performed to obtain the smoothed frequency deviation signal. and active power deviation signal Specifically, and The calculation formula is: = - ( =50Hz), = - ,in This represents the steady-state active power before the disturbance. To eliminate high-frequency measurement noise, [the following is used]: and Smoothing was performed using the HMA method.

[0090] Furthermore, the smoothed frequency deviation and active power deviation The instantaneous energy is calculated using the discrete TKEO operator and then normalized by sliding peak normalization to obtain the normalized dual-channel TKEO energy. and Specifically, the formula for calculating TKEO energy is:

[0091] (1)

[0092] (2)

[0093] The two TKEO outputs have different dimensions. To unify them to a comparable scale, sliding peak normalization is used to map them to [0,1]:

[0094] , (3)

[0095] in, The real-time frequency timing value of the power grid is collected by the synchronous phasor measurement unit at the nth sampling time. The rated frequency of the power system is fixed at 50Hz and used as the benchmark for frequency deviation calculation. The real-time active power time-series value is calculated and output by the synchronous phasor measurement unit at the nth sampling time. For frequency channel TKEO energy, TKEO energy is the active channel. To normalize the frequency channel TKEO energy, the original frequency channel TKEO energy is mapped to a dimensionless energy value in the [0,1] interval through long-time window sliding peak normalization. To normalize the active channel TKEO energy, the original active channel TKEO energy is mapped to a dimensionless energy value in the [0,1] interval by long-time window sliding peak normalization, N. l The length of the long time window is consistent with the LTA window length in step 104, so that the normalized baseline and the noise baseline are tracked synchronously.

[0096] It is worth mentioning that the present invention can use RMS normalization (Root Mean Square) instead of moving peak normalization (Equation 3): R-value quantiles are used instead. (Equation 9); MOESP+PEM replaces ARMAX (i.e., weights are retained in the PEM stage); Multi-PMU spatial voting expansion - each bus runs independently in steps 101-103, and the fusion center triggers confirmation of the entire network event based on K / N PMUs simultaneously.

[0097] Step 103: Perform causal delay-aligned coherent energy fusion on the normalized frequency channel TKEO energy and the normalized active channel TKEO energy to obtain dual-channel coherent energy.

[0098] It should be noted that the cross-correlation function is first used to solve the electromechanical response causal delay corresponding to the two normalized energies. Based on this delay, the normalized active channel TKEO energy is time-shifted and aligned. Then, the geometric mean of the time-aligned active channel energy and the normalized frequency channel TKEO energy is calculated. The calculation result is used as the coherent energy of the two channels and sent to the subsequent sequential detection stage of IS divergence (Itakura-Saito Divergence) of shared LTA (Long-Term Average).

[0099] Furthermore, step 103 may include the following sub-steps:

[0100] S31. The normalized active channel TKEO energy is delayed and aligned according to the causal time delay to obtain the time-delay aligned active channel TKEO energy.

[0101] S32. Calculate the geometric mean of the active channel TKEO energy and the normalized frequency channel TKEO energy after time delay alignment.

[0102] S33. Determine the coherent energy of the dual-channel system based on the geometric mean.

[0103] It should be noted that the a priori lower bound H of the inertia constant of the synchronous generator is... min Upper bound of damping coefficient D max and PMU sampling rate f s By converting the electromechanical response time constant of the oscillation equation, we obtain... Channel TKEO Energy Peak Leading Causal delay of the channel Specifically, The calculation formula is:

[0104] (4)

[0105] in The electromechanical response time constant derived from the oscillation equation. ∈[0.1,0.3] represents the empirical coefficients for relative position. This is the a priori lower bound for inertia (in seconds). This represents the upper limit of damping (unit: pu / Hz). PMU sampling rate (in Hz). This indicates rounding up to the nearest integer.

[0106] when , When prior knowledge is unavailable, the cross-correlation function of the two-channel energy during the initial steady-state phase can be used for calibration. The search interval is given by ±50% of the value calculated by equation (4).

[0107] Furthermore, the normalized TKEO energy after causal time delay alignment and Through geometric mean fusion, dual-channel coherent energy is obtained. (n). Specifically, The formula for calculating (n) is:

[0108] (5)

[0109] Equation (5) will Channel TKEO energy delay Later and Channel TKEO energy alignment multiplication square root.

[0110] It is worth mentioning that the present invention can use a weighted arithmetic mean instead of the geometric mean: replace equation (5) with the time-delay aligned arithmetic mean. It has slightly higher sensitivity to real events but loses the single-channel noise suppression properties, making it suitable for scenarios with high coherence between two-channel noise. The weighted allocation coefficient for the normalized active channel TKEO energy is limited to the range (0,1).

[0111] In this embodiment, the normalized active channel TKEO energy is delayed and aligned according to causal delay to obtain the delayed-aligned active channel TKEO energy; the geometric mean of the delayed-aligned active channel TKEO energy and the normalized frequency channel TKEO energy is calculated; based on the geometric mean, the coherent energy of the two channels is determined, and the causal delay is obtained through two calculation paths when performing the delay alignment operation. When the synchronous generator inertia prior lower bound Upper bound of damping coefficient When available, select an empirical coefficient for relative position with a value range of 0.1 to 0.3. Combined with PMU sampling rate Substitute into the formula The calculation is completed and the integer sampling offset points are obtained by rounding up. When the prior parameters of the unit cannot be obtained, the two normalized TKEO energy sequences during the initial steady-state period of the system are extracted to construct a cross-correlation summation function. The offset search interval is defined by ±50% of the parameter estimation delay, and all integer offsets within the interval are traversed. The number of offset points corresponding to the maximum value of the cross-correlation sum is taken as the causal delay. After determining the number of time delay sampling points, the normalized active channel TKEO energy of the complete time series is shifted backward as a whole. After timing matching is completed at each sampling point, the TKEO energy of the active channel after time delay alignment is obtained. Then, the aligned active channel energy at the same sampling time n is extracted. With normalized frequency channel TKEO energy The geometric mean is calculated by multiplying the two sets of values ​​and then taking the second square root of the product. Finally, the geometric mean obtained by taking the square root is directly defined as the dual-channel coherent energy at the current time. (n) completes the dual-channel energy coherent fusion and outputs it to the subsequent IS divergence sequential detection step with shared LTA.

[0112] Step 104: Perform sequential detection of IS divergence and joint decision-making of dual-criteria on the coherent energy of the dual-channel LTA, and output the perturbation event time, event confidence weight and per-sample positive part CUSUM increment.

[0113] The disturbance event time is a sampling point time sequence index where both the STA / LTA ratio exceeds the adaptive threshold and the IS-CUSUM cumulative amount exceeds the modified CUSUM decision threshold are met simultaneously in the dual-criteria joint decision. This time sequence index is bound to the PMU unified synchronization timestamp and fully records the precise time and location when the active power disturbance of the power grid is first identified by the algorithm. It contains two types of information: sampling sequence number and the corresponding actual electrical quantity acquisition time.

[0114] The event confidence weight is a dimensionless value calculated by dividing the accumulated IS-CUSUM in the current identification window by the corrected CUSUM decision threshold. The value is always greater than 0, and the magnitude of the value directly corresponds to the level of confidence of the disturbance signal. The calculation process relies on two sets of parameters: the total amount of accumulated disturbance features and the fixed threshold benchmark. It is used to quantify the degree of significance of the currently identified disturbance deviating from the steady-state energy baseline of the power grid.

[0115] The per-sample positive CUSUM increment is based on the IS divergence single-step CUSUM increment calculated at each independent sampling time in the full time series. The per-point perturbation quantization index is obtained by performing the maximum value operation between the single-step increment and the numerical zero. Only the effective perturbation components that deviate positively from the steady-state baseline are retained. The dataset contains the positive perturbation amplitude data corresponding to all sampling points, and invalid values ​​generated by perturbation-free and negative noise offsets are removed.

[0116] It should be noted that, based on the same long-time window baseline, the short-time window mean and long-time window mean of the dual-channel coherent energy are calculated synchronously. The STA / LTA ratio (Short-Time Average / Long-Time Average) is solved from the two sets of means. An adaptive threshold is generated in real time based on the long-time window mean. The IS divergence is calculated by substituting the dual-channel coherent energy and the long-time window mean into each sampling point. The IS divergence values ​​are continuously accumulated to obtain the IS-CUSUM cumulative amount. A fixed corrected CUSUM is pre-calculated. SUM (cumulative sum) decision threshold: When both the STA / LTA ratio is greater than the adaptive threshold and the IS-CUSUM cumulative amount is greater than the modified CUSUM decision threshold, the current sampling time series is recorded as the disturbance event time. The event confidence weight is calculated using the ratio of the IS-CUSUM cumulative amount to the modified CUSUM decision threshold. Finally, the entire time series is traversed with the disturbance event time as the boundary, and the maximum value of the IS divergence CUSUM increment and the zero value at each time is taken to obtain the sample-by-sample positive part CUSUM increment.

[0117] Furthermore, step 104 may include the following sub-steps:

[0118] S41. Based on the dual-channel coherent energy, calculate the short-term window mean and the long-term window mean;

[0119] S42. Calculate the STA / LTA ratio based on the short-term window mean and the long-term window mean;

[0120] S43. Calculate the adaptive threshold based on the long-term window mean;

[0121] S44. Calculate the IS divergence using dual-channel coherent energy and long-time window averaging;

[0122] S45. Determine the IS-CUSUM cumulative amount based on IS divergence;

[0123] S46. Calculate the corrected CUSUM decision threshold;

[0124] S47. Based on the STA / LTA ratio, adaptive threshold, IS-CUSUM cumulative amount and modified CUSUM decision threshold, perform dual-criteria joint decision to determine the timing of the disturbance event.

[0125] S48. Calculate the event confidence weight using the IS-CUSUM cumulative amount and the modified CUSUM decision threshold;

[0126] S49. Calculate the positive CUSUM increment per sample based on the time of the disturbance event.

[0127] It should be noted that the dual-channel coherent energy (n), by calculating the mean of the short-term window and the long-term window using a dual sliding window architecture and constructing an adaptive threshold, the STA / LTA ratio R(n) and the adaptive threshold are obtained. (n). Specifically, R(n) and The formula for calculating (n) is:

[0128] (6)

[0129] (7)

[0130] Where N s The short time window length (the value satisfies) ), N l The length of the long time window (the value covers the steady-state baseline before the disturbance occurs).

[0131] (8)

[0132] (9)

[0133] In formula (8) To prevent division by zero of small constants, in equation (9) / is the coefficient of variation, and k is the sensitivity coefficient (k∈[3,5]).

[0134] For the IS-CUSUM sequential detection pathway: the dual-channel coherence energy (n) and long-term window mean The single-step CUSUM increment is calculated using IS divergence and recursively accumulated to obtain the IS-CUSUM cumulative amount S(n). Specifically, the formula for calculating S(n) is:

[0135] (10)

[0136] in, The Itakura-Saito (IS) divergence is a statistic used to quantify the coherent energy of a single-time dual-channel circuit. With steady-state baseline mean The degree of deviation between the distribution values; the larger the value, the more significant the current signal deviates from the steady-state baseline of the power grid. To monitor the ratio of energy to steady-state reference energy in real time; The dual-channel coherent energy at the current sampling time (n); The steady-state baseline mean; the reference mean of the IS-CUSUM pathway is directly taken from the steps above. This allows both pathways to share the same statistical benchmark.

[0137] (11)

[0138] In the form of (11) For reference average, from (n) Calculate the single-step CUSUM increment and recursively accumulate it:

[0139] (12)

[0140] (13)

[0141] The drift reference constant c is: .

[0142] CUSUM Decision Threshold F corrected Using a physical parameterization calibration method, a quantitative relationship is established between the threshold and the minimum detectable energy increment, the maximum allowable delay, and the baseline noise level:

[0143] (14)

[0144] in The smallest detectable energy increment (IS divergence unit). Where c is the maximum allowable detection delay (number of sampling points), and c is the drift reference constant (estimated from the steady-state IS divergence during the initialization phase). N is the short time window length. l For the length of the long window, and For the correction factor ( <1, >1). This threshold is used to determine the cumulative amount of IS divergence.

[0145] Furthermore, the STA / LTA ratio R(n) and the IS-CUSUM cumulative amount S(n) are used to jointly determine the disturbance event using a dual-criteria method, and confidence information is extracted to obtain the event time t. event Confidence weight w conf and per-sample information density L + (k). Specifically, the joint decision conditions and confidence quantification formulas are as follows:

[0146] (15)

[0147] (16)

[0148] (17)

[0149] In equation (16) w conf ≥1 reflects the cumulative threshold exceedance multiple of CUSUM. Equation (17) defines the identification window [ The positive CUSUM increment L for each sample within the range + (k), where N pre N post These represent the number of samples for the identification window before and after the event time, respectively. This serves as the sampling time sequence index for disturbance events, the sampling point number corresponding to the simultaneous fulfillment of the dual-criteria joint judgment condition, and is bound to the PMU unified synchronization timestamp to mark the precise time sequence location where the power grid active power disturbance is first identified by the algorithm. After the event is confirmed, S(n) is reset to zero, and the algorithm locks T. lock Avoid triggering the time repeatedly.

[0150] It is worth mentioning that the present invention can employ... - Replace IS divergence with other members of the divergence family: Use instead =1 (KL divergence, distribution parameter required) or =2 (Euclidean distance / squared error), IS divergence ( The scale invariance of (=0) is not possessed by other members, making it the preferred feature of this invention.

[0151] In this embodiment, the following conditions are first met during implementation: Short window length Ns Long window length N l With dual-channel coherent energy (n) is the input, which takes the current sampling point forward N times. s N l The short-time window mean is obtained by summing the energy time series and dividing by the corresponding window length. With long window mean The long-term window fully covers the steady-state baseline before the disturbance, and the mean of the short-term window is divided by the superimposed minimum constant. The STA / LTA ratio R(n) is obtained from the long-term mean value, and then the mean of the long-term mean value sequence is calculated. Standard deviation Combining the sensitivity coefficient k with a value range of 3 to 5 and the minimum constant Real-time adaptive threshold generation calculation averaging over a long window The steady-state reference mean, used as the unified standard for IS divergence calculation, is substituted into the Itakura-Saito divergence formula to calculate the difference in distribution between the dual-channel coherent energy and the steady-state baseline. Then, the drift reference constant c, obtained from the long-time-window steady-state IS divergence mean, is subtracted to obtain the single-step CUSUM increment. The IS-CUSUM cumulative amount is obtained by continuously accumulating the single-step increment through a recursive maximum value rule. Substituting the minimum detectable energy increment, maximum allowable detection delay, long and short time window lengths, and satisfying the following conditions... <1、 The two sets of correction coefficients >1 complete the physical parameterization solution, resulting in the corrected CUSUM decision threshold F. corrected Synchronous verification Two judgment conditions are met simultaneously. When both conditions are met, the current sampling time sequence index is recorded as the time t of the disturbance event. event The event confidence weight w is calculated by dividing the accumulated IS-CUSUM at that moment by the corrected CUSUM decision threshold. conf The area centered on the time of the disturbance event includes the preceding N. pre Post-N post The identification window for sampling points is defined. The single-step CUSUM increment L(k) of all sampling points within the window is traversed. The maximum value between each sampling point's increment and zero is taken to obtain the CUSUM increment L for each sample's positive part within the window. + (k) After determining a disturbance event, reset the IS-CUSUM accumulation S(n) and lock it for a fixed duration T. lock This avoids triggering the same decision repeatedly within the same disturbance duration.

[0152] Step 105: Based on the event time, event confidence weight, and per-sample positive CUSUM increment, perform weighted ARMAX model identification and inertia constant extraction to obtain the estimated value of the power system inertia constant.

[0153] It should be noted that, firstly, the frequency deviation and active power deviation within the window of the effective period of the disturbance are extracted based on the event time. Then, the event confidence weight and the per-sample positive part CUSUM increment are integrated. The weight of each sample is calculated by combining Gaussian time nearest neighbor and information density. The ARMAX model parameters are identified by performing weighted least squares operation using the sample weights. Based on the parameters, a second-order discrete transfer function is constructed and transformed into a first-order continuous transfer function. The coefficient mapping relationship is determined by comparing it with the first-order theoretical transfer function of the synchronous generator swing equation. Finally, the estimated value of the power system inertia constant is obtained by combining the mapping relationship and the ARMAX model parameters.

[0154] Furthermore, step 105 may include the following sub-steps:

[0155] S51. Based on the event time, extract the window frequency deviation and window active power deviation within the effective period of the disturbance.

[0156] S52. Perform Gaussian time nearest neighbor and information density combined weighting calculation on the event confidence weight and the per-sample positive part CUSUM increment to determine the sample weight;

[0157] S53. Using the frequency deviation and active power deviation within the window as inputs, the ARMAX model parameters are identified by substituting the sample weights into the weighted least squares method.

[0158] S54. Construct a second-order discrete transfer function using ARMAX model parameters, and convert the second-order discrete transfer function into a first-order continuous transfer function.

[0159] S55. Compare the coefficients of the first-order continuous transfer function with the first-order theoretical transfer function corresponding to the swing equation of the synchronous generator to determine the coefficient mapping relationship.

[0160] S56. Calculate the estimated value of the power system inertia constant based on the coefficient mapping relationship and ARMAX model parameters.

[0161] It should be noted that the event time t event and per-sample information density L + (k), the weight W(k) of each sample within the identification window is obtained by combining Gaussian time nearest neighbors and information density weighting. Specifically, the formula for calculating W(k) is:

[0162] (18)

[0163] in The width parameter of the Gaussian time window. =max k L + (k). To simultaneously utilize the steady-state baseline data before the disturbance and the transient electromechanical response data after the disturbance during the identification phase, with t event Center the identification window [ ]——N pre The pre-event samples provide a steady-state operating benchmark, N post The frequency change process containing inertia information is captured in the sample after each event. Output from the IS-CUSUM pathway in the above steps.

[0164] It is worth mentioning that the confidence level of this invention is only used for window length adjustment: the sample-by-sample weights in equation (18) are replaced with those based on w conf Window length dynamic adjustment — w conf Larger time window length shortened, w conf The window length is extended to near 1. This makes the implementation simpler but doesn't delve into the sample weighting level.

[0165] Furthermore, the input within the recognition window will be... Output The inertia constant is estimated by using the sample weights W(k) and extracting them through weighted least squares ARMAX identification and order reduction. Specifically, The calculation process is as follows:

[0166] Using ARMAX model description arrive The input-output relationship is generally in the form of: , where u(k) = For exogenous input, y(k) = For autoregressive output, For backward shift operators, Let be the model residual noise term at sampling time k, representing random measurement noise and random disturbance components in the active power deviation and frequency deviation time series data that cannot be fitted by the linear dynamic relationship of the model. Let A and B be of order 2 and C be of order 1 to identify the parameter vector. , of which are It is a second-order autoregressive polynomial The first and second order coefficients describe the autoregressive coupling effect of historical frequency deviation on current frequency deviation. For second-order exogenous input polynomials The zeroth, first, and second order coefficients describe the excitation effect of historical active power deviation disturbances on frequency deviation. First-order noise polynomial The first-order coefficients describe the impact of historical residual noise on the current frequency output. To identify the parameter vector, T is the transpose. Let the identification window contain N... w =N pre +N post For each sample, the input and output data within the window are used to form a regression matrix according to the ARMAX regression structure. (Each row corresponds to the historical values ​​of y, u, and e at a given time point), and the target vector is formed by stacking the target values ​​y(k). The diagonal weight matrix W = diag{W(1), ..., W(N) is constructed using the weights in equation (18). w Solve the WLS:

[0167] (19)

[0168] in, The optimal parameter vector of the ARMAX model obtained by weighted least squares solution includes all autoregressive, input, and noise coefficients. The regression matrix has dimension N. w Rows of 6 columns, each row storing the historical output, historical input, and historical noise time series data at a single sampling time, N w The total number of samples in the identification window; T is the transpose; The target vector is a one-dimensional column vector, composed of the in-window frequency deviations at all times within the identification window. Stacked structure; after the second-order discrete transfer function obtained from equation (19) is reduced to the first-order continuous form G'(s)=b' / (s+a'), it is compared with the coefficients of the first-order transfer function of the oscillation equation to extract the estimated power system inertia constant:

[0169] (20)

[0170] Wherein, G'(s) is the first-order complex frequency domain continuous transfer function obtained by order reduction through discrete-continuous transformation, which characterizes the electromechanical dynamic relationship between active power deviation and frequency deviation; b' is the numerator constant term of the first-order continuous transfer function, corresponding to the inertia-related term of the synchronous generator swing equation, and the two can be combined to calculate the system inertia constant; s is the Laplace complex frequency domain operator, used for continuous system dynamics modeling; a' is the coefficient of the first-order linear term in the denominator of the first-order continuous transfer function, corresponding to the damping-related term of the synchronous generator swing equation.

[0171] In this embodiment, the event time t is used as... event Centered on the point, cut off N points forward. pre One preceding steady-state sample, N truncated backwards post A complete identification window is formed by several transient disturbance samples, and the active power deviation within the window corresponding to the sampling point is extracted from the window time sequence. Frequency deviation within the window Set the Gaussian time window width parameter Construct an exponential time weight term centered on the event time, and calculate the maximum value of the positive part CUSUM increment for all samples within the window. Divide the increment of each sampling point by The normalized information density term is obtained by taking the minimum value of the value 1. After adjusting the amplitude by combining the event confidence weight, the time weight term is multiplied by the information density term. The sample weight W(k) corresponding to each sample is calculated according to equation (18). All sample weights are arranged along the diagonal of the matrix to generate a diagonal weight matrix W. The order of the autoregressive and exogenous inputs of the ARMAX model is set to 2, and the noise order is 1. As exogenous input u(k), As the autoregressive output y(k), based on N within the window w =N pre +N post A six-column regression matrix was constructed from the time series data of the sample. With the target vector Substituting the diagonal weight matrix into the weighted least squares formula of equation (19), the six-dimensional ARMAX parameter vector is obtained. A second-order discrete transfer function is generated based on the combination of parameter vectors. Then, it is reduced to a first-order continuous transfer function G'(s)=b' / (s+a') through discrete continuous transformation. The coefficients of this transfer function are compared with those of the first-order theoretical transfer function of the synchronous generator swing equation. The mapping equation between a', b' and the system inertia is established. Finally, the coefficients of the first-order continuous transfer function are substituted into equation (20) to complete the calculation and obtain the estimated value of the power system inertia constant. .

[0172] As a comparison of technical effects, existing technologies can be used as a reference. Under normal operating conditions of the power system, active power disturbance events (including large disturbances and frequent small load fluctuations) can be automatically detected from the continuous frequency and active power signals synchronously measured by the PMU, and inertia estimation can be triggered immediately after the event is confirmed.

[0173] The existing mainstream approaches include: (a) Large disturbance triggering method - using the oscillation equation, H is directly calculated based on the known magnitude of the disturbance and the RoCoF after the disturbance, but large disturbance events are infrequent and cannot work continuously; (b) Continuous environmental data identification method - system identification is performed on random fluctuation data under normal operating conditions for a long window (200~500s), which has a slow response; (c) Event triggering method - the time of occurrence of the disturbance event is detected online first, and then inertia estimation is performed within the event window, which has both continuous online capability and fast response capability, but the sensitivity of existing event detection methods is constrained by a fixed threshold.

[0174] Based on the above, the shortcomings of existing technologies can be roughly divided into three parts:

[0175] Disadvantage 1: Fixed thresholds cannot adapt to changes in operating conditions, and single-point decisions cannot utilize the accumulation of time-series evidence.

[0176] The existing TKEO threshold is 1×10 -7 Offline calibration under specific operating conditions; frequency baseline drift during periods of high wind power generation (σ) when grid operating conditions change. f The TKEO peak value can reach 0.02Hz, while it is only 0.002Hz under light load at night, a difference of an order of magnitude. Under high fluctuation conditions, the peak value of TKEO during normal fluctuations can approach the threshold, leading to false triggering. Under low fluctuation conditions, the peak value of TKEO during MW-level load switching is far below the threshold, leading to missed detection. In addition, existing methods compare the TKEO value of a single sampling point with a fixed threshold to make a binary decision—the accumulation of slight energy deviations (each below the threshold) from multiple consecutive sampling points is sufficient to support the detection decision, but single-point comparison cannot utilize this time-series accumulation information. Another existing method, CUSUM, provides a sequential accumulation framework, but its log-likelihood ratio depends on the specific distribution family assumption (Gaussian / Gamma), resulting in a distribution family mismatch problem when transferred to the TKEO energy space.

[0177] Disadvantage 2: The single signal source lacks redundancy verification and does not explicitly utilize the causal timing differences between the two channels.

[0178] Current TKEO methods only apply to Δf. At the instant of a disturbance event, ΔP is a step-type electrical quantity change (its TKEO instantly generates a spike), while Δf is an integral electromechanical response (requiring an integration process of the rotor motion equations, which takes tens to hundreds of milliseconds to form an energy gradient that can be captured by TKEO). Even with the introduction of a ΔP channel, if it is only fused with the Δf channel in a linear weighted manner, this causal timing difference cannot be utilized. The other two existing methods also only operate on a single signal channel. Furthermore, existing methods use only a single PMU for detection, lacking spatial redundancy and unable to distinguish between local noise spikes and network-wide disturbance events.

[0179] Disadvantage 3: The detection end and the estimation end only have an interface connection, and no data loop is formed at the algorithm layer.

[0180] Current methods initiate ARMAX identification after TKEO triggering, but there is no coupling between the energy of TKEO and the sample weights within the ARMAX identification window—"how confident the detection end is" and "how the estimation end uses the data" are two separate issues. Two other existing methods don't even involve inertia estimation. When the perturbation amplitude is small and the detection confidence is low, ARMAX identification still processes all samples within the window using equal-weighted least squares. Samples in the low signal-to-noise ratio phase drag down the overall identification accuracy, leading to an amplification of the H estimation variance. Furthermore, in another existing method, the STA / LTA and CUSUM pathways are independent at the statistical parameter estimation level—the threshold λ of STA / LTA and the reference mean μ0 of CUSUM are maintained by two independent estimators. Inconsistent parameter updates during the transition period between operating conditions can easily lead to a surge in false alarms or missed detections.

[0181] Therefore, the object of the present invention is:

[0182] The objective of addressing Disadvantage 1 is to construct a sequential detection mechanism that adaptively adjusts the detection threshold according to the power grid operating conditions and can make decisions based on the accumulation of time-series evidence—replacing the fixed threshold with an adaptive baseline and replacing single-point comparison with statistically optimal sequential accumulation.

[0183] To address the second drawback, the objective is to construct a coherent energy fusion method that explicitly utilizes the causal time sequence difference between the two channels ΔP and Δf. This method involves physically calculating the causal time delay between the two channels using the oscillation equation, aligning them, and taking the geometric mean. This ensures that the actual disturbance is only confirmed when both channels are "simultaneously significant" according to the physical causal time sequence.

[0184] To address the third drawback, the objective is to construct a joint identification mechanism that feeds back detection confidence information to the estimation end—making the statistical confidence of the detection end (CUSUM cumulative threshold multiple and point-by-point increment) the sample weight for ARMAX identification, thereby achieving a data closure loop in the detection-estimation algorithm layer; and simultaneously unifying the statistical references of the STA / LTA and CUSUM pathways onto the same baseline estimator.

[0185] Specifically, the present invention provides a method for estimating the inertia of a power system:

[0186] S1: Data Acquisition and Dual-Channel TKEO Energy Feature Extraction – Acquisition from PMU , The HMA method was used for smoothing; the TKEO energy of the two channels was calculated separately and normalized according to the sliding peak value;

[0187] S2: Coherent Energy Fusion with Causal Time Delay Alignment – ​​From the Oscillation Equation Time Constant Calculate the physical causal delay from ΔP to Δf channel TKEO energy peak. ;Will Delay Later and The geometric mean is taken to obtain the coherent energy of the two channels. ;

[0188] S3: Sequential IS divergence detection and dual-criteria joint decision-making with shared LTA – To ensure consistent input, the STA / LTA pathway and the IS-CUSUM pathway are run in parallel, with both pathways sharing the same LTA(n) as a statistical reference; dual criteria (R(n) > 0.05) are used. And S(n)>F corrected When both conditions are met, the disturbance event is confirmed, and the confidence level w is recorded. conf and the point-by-point CUSUM increment within the window {L + (k)};

[0189] S4: Detection Confidence-Weighted ARMAX Joint Identification – Based on {L + (k)} Construct WLS sample weights and use them to identify the ARMAX model; after second-order to first-order reduction and coefficient comparison, extract the inertia constant. .

[0190] As can be seen from the above, the innovation of this invention lies in:

[0191] Innovation Point 1: Introducing a causal time delay alignment and geometric mean fusion mechanism for dual-channel TKEO energy characteristics of ΔP and Δf, so that the confirmation of disturbance is based on the premise that significant energy appears simultaneously in both channels according to physical causal time sequence.

[0192] The electromechanical response time constant from the oscillation equation Calculate the physical delay of the energy peak of channel ΔP leading that of channel Δf. (Equation 4) The geometric mean of the normalized TKEO energy of the two channels after time delay alignment is taken to obtain the coherent energy of the two channels. (Equation 5). The geometric mean structure ensures that the coherent energy only reaches a high value when significant energies appear simultaneously in both channels, thus suppressing noise spikes in the single channel.

[0193] Innovation Point 2: Construct a sequential detection architecture in which the STA / LTA pathway and the IS-CUSUM pathway share the same long-term mean LTA(n) as a statistical reference, enabling the two pathways to synchronously adapt when switching operating conditions.

[0194] LTA(n) simultaneously provides a normalized baseline for the STA / LTA pathway (Equations 8-9) and a divergence reference point for the IS-CUSUM pathway (Equations 11-12). The CUSUM increment is calculated using IS divergence (Equation 10), and its three properties—closed-form expression, scale invariance, and independence from the distribution family assumption—match the structural characteristics of the TKEO energy space.

[0195] Innovation Point 3: Establishing IS-CUSUM single-step incremental L at the detection end+ (k) The joint mechanism of ARMAX identification and backfeeding to the estimation end enables the perturbation information density independently determined by the detection end to directly modulate the contribution weight of each sample in the estimation end.

[0196] The sample weights W(k) are constructed using Gaussian time nearest neighbors and CUSUM incremental normalization two-factor methods (Equation 18), which drive weighted least squares identification (Equation 19). + (k) Output from the IS-CUSUM pathway (independently determined by the detection end), this information is directly entered into the identification weighting without the intervention of the estimation end's post-event residual.

[0197] Compared with the prior art, the advantages of the present invention are as follows:

[0198] Advantage 1 (corresponding to innovation point 1): Compared with the existing method that only utilizes a single-channel TKEO of Δf, this invention utilizes dual-channel information of ΔP and Δf, and calculates the causal delay through the oscillation equation. This method achieves temporal alignment of TKEO energies between two channels. Compared to simple linear weighted fusion, geometric mean fusion (Equation 5) yields a zero result when the energy of either channel is zero, preventing isolated noise spikes in a single channel from triggering detection. This scheme solidifies the physical fact that "ΔP abruptly occurs before Δf" into the algorithm structure, thus suppressing noise events that do not satisfy this causal temporal order.

[0199] Advantage 2 (corresponding to innovation point 2): Compared with the fixed threshold single-point decision of the existing method, the STA / LTA path of the present invention provides an adaptive threshold by tracking the change of operating conditions through LTA(n) (Equation 9), and the IS-CUSUM path accumulates the weak energy deviations of multiple sampling points over time (Equation 13), so that the detection capability is not limited by the fixed threshold under low fluctuation conditions and is not falsely triggered under high fluctuation conditions. Compared with the log-likelihood ratio CUSUM of another existing method, IS divergence (Equation 10) has scale invariance, and the baseline magnitude change caused by the change of operating conditions does not contaminate the accumulated amount. Compared with the independent maintenance of the parameters of the STA / LTA and CUSUM paths of another existing method, sharing LTA(n) (Equation 11) eliminates the false alarms or missed detections caused by the inconsistency of the parameter updates of the two paths during the transition period of the change of operating conditions.

[0200] Advantage 3 (corresponding to innovation point 3): Compared with the existing equal-weight least squares ARMAX identification method, this invention uses IS-CUSUM single-step incremental L + (k) The constructed sample weights (Equation 18) feed back the statistical confidence information from the detection end to the estimation end—in strong perturbation events, high L + (k) Samples are given high weights, focusing on high signal-to-noise ratio sub-windows; in weakly perturbed events, L +The (k) distribution tends to flatten, and the weights degenerate into time nearest neighbor dominance. This mechanism ensures that the identification accuracy is not significantly amplified by the estimation variance due to the drag from low signal-to-noise ratio samples when the perturbation amplitude changes.

[0201] In this embodiment of the invention, a method for estimating the inertia of a power system is provided. This method acquires synchronous phasor measurement unit (TPMU) data, PMU active power data, and steady-state active power before disturbance. Based on the TPMU data, PMU active power data, and steady-state active power before disturbance, dual-channel TKEO energy feature extraction is performed to obtain normalized frequency channel TKEO energy and normalized active channel TKEO energy. The normalized frequency channel TKEO energy and normalized active channel TKEO energy are then fused using causal time delay alignment to obtain dual-channel coherent energy. This invention performs sequential IS divergence detection and dual-criteria joint decision-making on the shared LTA of dual-channel coherent energy, outputting the disturbance event time, event confidence weight, and per-sample positive part CUSUM increment. Based on the event time, event confidence weight, and per-sample positive part CUSUM increment, weighted ARMAX model identification and inertia constant extraction are performed to obtain the estimated value of the power system inertia constant. Based on the above scheme, this invention collects synchronous phasor measurement unit data, PMU active power data, and pre-disturbance steady-state active power, and then sequentially performs dual-channel TKEO energy feature extraction and causal time delay alignment. The entire process, including coherent energy fusion, sequential IS divergence detection and dual-criteria joint decision-making using shared LTA, weighted ARMAX model identification, and inertia constant extraction, relies on dual-channel temporal coherent fusion to replace the traditional single-channel fixed threshold single-point judgment logic. By fusing two TKEO energy streams through time delay alignment geometric mean, it retains only the real disturbance signals that conform to electromechanical timing rules, filters out random noise from isolated channels, avoids noise-induced false triggering of invalid identification data, and generates adaptive thresholds based on real-time updated long-time-window baselines and accumulates weak temporal disturbance characteristics to adapt to different load fluctuation conditions of the power grid. The system fully captures small load disturbances to expand the effective identification samples. The detection process simultaneously outputs the event confidence weight and the per-sample positive CUSUM increment, which are directly used for subsequent identification weighting. The sample weights are distinguished according to the disturbance confidence, and the interference of low signal-to-noise ratio data on ARMAX parameter solving is suppressed, reducing model identification bias. The power system inertia constant is calculated by relying on accurate model parameters. The entire integrated detection and identification process eliminates the estimation errors caused by false triggering, missed detection of effective disturbances, and indiscriminate equal weighting of the existing fixed threshold TKEO scheme, effectively improving the accuracy of power system inertia estimation.

[0202] Please see Figure 2 , Figure 2 This is a structural block diagram of a power system inertia estimation device provided in Embodiment 2 of the present invention.

[0203] The present invention provides a power system inertia estimation device, comprising:

[0204] The acquisition module 201 is used to acquire synchronous phasor measurement unit data, PMU active power data and steady-state active power before disturbance;

[0205] The feature extraction module 202 is used to extract dual-channel TKEO energy features based on the synchronous phasor measurement unit data, PMU active power data and steady-state active power before disturbance, so as to obtain normalized frequency channel TKEO energy and normalized active channel TKEO energy.

[0206] The fusion module 203 is used to perform coherent energy fusion of the normalized frequency channel TKEO energy and the normalized active channel TKEO energy with causal time delay alignment to obtain dual-channel coherent energy.

[0207] Decision module 204 is used to perform sequential detection of IS divergence of shared LTA and joint decision of dual criteria on the coherent energy of dual channels, and outputs the perturbation event time, event confidence weight and per-sample positive part CUSUM increment.

[0208] The inertia extraction module 205 is used to perform weighted ARMAX model identification and inertia constant extraction based on the event time, event confidence weight and sample-by-sample positive CUSUM increment, so as to obtain the estimated value of the power system inertia constant.

[0209] Furthermore, the feature extraction module 202 is specifically used for:

[0210] Calculate the frequency deviation based on the frequency data and the rated frequency;

[0211] The active power deviation is calculated based on the PMU active power data and the steady-state active power before the disturbance.

[0212] The frequency deviation and active power deviation are smoothed and denoised separately to obtain the smoothed frequency deviation and smoothed active power deviation.

[0213] Based on the smoothed frequency deviation and smoothed active power deviation, calculate the TKEO energy of the frequency channel and the TKEO energy of the active channel.

[0214] The frequency channel TKEO energy and the active channel TKEO energy are normalized to obtain the normalized frequency channel TKEO energy and the normalized active channel TKEO energy.

[0215] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the specific working process of the above-described device and module can be referred to the corresponding process in the foregoing method embodiments, and will not be repeated here.

[0216] This invention also provides a computer device, including a memory and a processor, wherein the memory stores a computer program; when the computer program is executed by the processor, the processor performs the steps of the power system inertia estimation method as described in the above embodiments.

[0217] This invention also provides a computer-readable storage medium storing a computer program / instructions thereon, which, when executed by a processor, implements the steps of the power system inertia estimation method as described in the above embodiments.

[0218] This invention also provides a computer program product, including a computer program stored on a non-transitory computer-readable storage medium, the computer program including program instructions, wherein when the program instructions are executed by a computer, the computer performs the steps of the power system inertia estimation method as described in the above embodiments.

[0219] In the several embodiments provided in this application, it should be understood that the disclosed apparatus and methods can be implemented in other ways. For example, the apparatus embodiments described above are merely illustrative; for instance, the division of units is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection shown or discussed may be through some interfaces; the indirect coupling or communication connection between apparatuses or units may be electrical, mechanical, or other forms.

[0220] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

[0221] Furthermore, the functional units in the various embodiments of the present invention can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit. The integrated unit can be implemented in hardware or as a software functional unit.

[0222] If the integrated unit is implemented as a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or all or part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods of the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0223] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention 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 the present invention.

Claims

1. A method for estimating the inertia of a power system, characterized in that, include: Acquire synchronous phasor measurement unit data, PMU active power data, and steady-state active power before disturbance; Based on the synchronous phasor measurement unit data, the PMU active power data and the steady-state active power before the disturbance, dual-channel TKEO energy feature extraction is performed to obtain normalized frequency channel TKEO energy and normalized active channel TKEO energy. The normalized frequency channel TKEO energy and the normalized active channel TKEO energy are fused by causal delay alignment to obtain dual-channel coherent energy. The IS divergence sequential detection and dual-criteria joint decision are performed on the dual-channel coherent energy using shared LTA, and the perturbation event time, event confidence weight and sample-by-sample positive part CUSUM increment are output. Based on the event time, the event confidence weight, and the per-sample positive CUSUM increment, a weighted ARMAX model identification and inertia constant extraction are performed to obtain an estimated value of the power system inertia constant.

2. The power system inertia estimation method according to claim 1, characterized in that, The step of extracting dual-channel TKEO energy features based on the synchronous phasor measurement unit data, the PMU active power data, and the pre-disturbance steady-state active power to obtain normalized frequency channel TKEO energy and normalized active channel TKEO energy includes: Based on the frequency data in the synchronous phasor measurement unit data and the rated frequency, calculate the frequency deviation; The active power deviation is calculated based on the PMU active power data and the steady-state active power before the disturbance. The frequency deviation and the active power deviation are smoothed and denoised respectively to obtain the smoothed frequency deviation and the smoothed active power deviation. Based on the smoothed frequency deviation and the smoothed active power deviation, calculate the frequency channel TKEO energy and the active channel TKEO energy. The frequency channel TKEO energy and the active channel TKEO energy are normalized to obtain normalized frequency channel TKEO energy and normalized active channel TKEO energy.

3. The power system inertia estimation method according to claim 1, characterized in that, The coherent energy fusion of the normalized frequency channel TKEO energy and the normalized active channel TKEO energy, causally time-delay aligned, to obtain dual-channel coherent energy includes: The normalized active channel TKEO energy is delayed and aligned according to causal time delay to obtain the time-delay aligned active channel TKEO energy. Calculate the geometric mean of the TKEO energy of the active channel after time delay alignment and the TKEO energy of the normalized frequency channel; Based on the geometric mean, the coherent energy of the two channels is determined.

4. The power system inertia estimation method according to claim 1, characterized in that, The sequential detection of IS divergence and joint decision-making based on dual criteria for shared LTA of the dual-channel coherent energy outputs the perturbation event time, event confidence weight, and per-sample positive part CUSUM increment, including: Based on the dual-channel coherent energy, the short-term window mean and the long-term window mean are calculated; Calculate the STA / LTA ratio based on the short window mean and the long window mean; Calculate the adaptive threshold based on the long-term window mean; The IS divergence is calculated using the dual-channel coherent energy and the long-time window mean. Based on the IS divergence, determine the IS-CUSUM cumulative amount; Calculate the corrected CUSUM decision threshold; Based on the STA / LTA ratio, the adaptive threshold, the IS-CUSUM cumulative amount, and the modified CUSUM decision threshold, a dual-criteria joint decision is performed to determine the timing of the disturbance event. The event confidence weight is calculated using the IS-CUSUM cumulative amount and the modified CUSUM decision threshold. Calculate the per-sample positive CUSUM increment based on the time of the disturbance event.

5. The power system inertia estimation method according to claim 1, characterized in that, The step of performing weighted ARMAX model identification and inertia constant extraction based on the event time, the event confidence weight, and the per-sample positive CUSUM increment to obtain the power system inertia constant estimate includes: Based on the event time, extract the window frequency deviation and window active power deviation within the effective period of the disturbance. The event confidence weight and the per-sample positive part CUSUM increment are calculated by performing a Gaussian time nearest neighbor and information density weighted calculation to determine the sample weight; Using the frequency deviation within the window and the active power deviation within the window as inputs, the ARMAX model parameters are identified by substituting the sample weights into the weighted least squares method. A second-order discrete transfer function is constructed using the ARMAX model parameters, and then the second-order discrete transfer function is converted into a first-order continuous transfer function. The coefficients of the first-order continuous transfer function are compared with the first-order theoretical transfer function corresponding to the swing equation of the synchronous generator to determine the coefficient mapping relationship. The estimated value of the power system inertia constant is calculated based on the coefficient mapping relationship and the ARMAX model parameters.

6. A power system inertia estimation device, characterized in that, include: The acquisition module is used to acquire data from the synchronous phasor measurement unit, PMU active power data, and steady-state active power before disturbance. The feature extraction module is used to perform dual-channel TKEO energy feature extraction based on the synchronous phasor measurement unit data, the PMU active power data and the steady-state active power before the disturbance, to obtain the normalized frequency channel TKEO energy and the normalized active channel TKEO energy. The fusion module is used to perform causal delay-aligned coherent energy fusion of the normalized frequency channel TKEO energy and the normalized active channel TKEO energy to obtain dual-channel coherent energy. The decision module is used to perform sequential detection of IS divergence of shared LTA and joint decision of dual criteria on the coherent energy of the dual channels, and outputs the perturbation event time, event confidence weight and per-sample positive part CUSUM increment. The inertia extraction module is used to perform weighted ARMAX model identification and inertia constant extraction based on the event time, the event confidence weight, and the per-sample positive part CUSUM increment, so as to obtain the estimated value of the power system inertia constant.

7. The power system inertia estimation device according to claim 6, characterized in that, The feature extraction module is specifically used for: Calculate the frequency deviation based on the frequency data and the rated frequency; The active power deviation is calculated based on the PMU active power data and the steady-state active power before the disturbance. The frequency deviation and the active power deviation are smoothed and denoised respectively to obtain the smoothed frequency deviation and the smoothed active power deviation. Based on the smoothed frequency deviation and the smoothed active power deviation, calculate the frequency channel TKEO energy and the active channel TKEO energy. The frequency channel TKEO energy and the active channel TKEO energy are normalized to obtain normalized frequency channel TKEO energy and normalized active channel TKEO energy.

8. An electronic device, characterized in that, The system includes a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor performs the steps of the power system inertia estimation method as described in any one of claims 1-5.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed, it implements the power system inertia estimation method as described in any one of claims 1-5.

10. A computer program product, characterized in that, The computer program product includes a computer program stored on a non-transitory computer-readable storage medium, the computer program including program instructions, wherein when the program instructions are executed by a computer, the computer performs the steps of the power system inertia estimation method as described in any one of claims 1-5.

Images