Power system inertia identification method based on ARMAX and adaptive recursive optimization

Abstract

An electric power system inertia identification method based on ARMAX and adaptive recursive optimization comprises the steps that active power and frequency data signals at a bus are collected in real time, the data signals are preprocessed, and a standard data set is formed; establishing an active-frequency dynamic response model of the new energy power system, and deducing frequency response transfer functions of the system in large and small disturbance scenes; an ARMAX model is constructed, and the coupling relation between parameters of the ARMAX model and the inertia is established; and solving parameters to be identified in the ARMAX model by adopting an adaptive recursive least square algorithm. Designing a fuzzy controller to carry out online correction on the adaptive forgetting factor, and calculating to obtain an optimal parameter vector; and substituting the optimal parameter vector into the ARMAX model, and calculating the inertia time constant of each unit. The method solves the problems of data saturation and noise sensitivity of a traditional method; and meanwhile, the inertia identification precision is improved through a time-varying gain matrix and a covariance iteration mechanism, and full-scene coverage from noise-like small disturbance to transient large disturbance is realized.

Description

Technical Field

[0001] The present invention relates to the technical field of power system parameter identification, and particularly relates to a power system inertia identification method based on ARMAX and adaptive recursive optimization. Background Art

[0002] New energy power system inertia identification refers to identifying the changes in inertia in the system through real-time monitoring and data analysis means. With the large-scale grid connection of new energy power generation equipment such as wind energy and solar energy, the number of synchronous machines in the power system has decreased significantly, and the system inertia level has decreased accordingly. This reduction in inertia results in a slower response speed of the system to frequency fluctuations and increases the instability risk of the system in the event of sudden faults.

[0003] Currently, inertia identification technology, as the core means to real-time perceive the dynamic characteristics of the system and quantify the virtual inertia support ability, has become a key technical requirement for ensuring the frequency safety of new power systems. However, the existing inertia identification methods generally have the following technical bottlenecks: First, new energy units are connected to the grid through power electronic interfaces, and their inertia response mechanisms are significantly different from those of traditional synchronous machines. Existing linear models are difficult to accurately represent the non-linear dynamic characteristics of new power systems; Second, random noise interference and data saturation phenomena in grid operation easily lead to insufficient accuracy of existing identification technologies; Third, the existing technologies have insufficient adaptability to multi-disturbance scenarios and are difficult to cover the inertia dynamic tracking requirements of various complex working conditions such as small noise-like disturbances and transient large disturbances. Summary of the Invention

[0004] To solve the problem of insufficient accuracy of traditional identification methods in new energy power systems due to noise interference, non-linear dynamics, and complex disturbance scenarios. The present invention proposes a power system inertia identification method based on ARMAX and adaptive recursive optimization. This method constructs an ARMAX model and combines an adaptive recursive optimization algorithm to real-time track the changes in the dynamic working conditions of the power system, and uses a dynamic forgetting factor and a fuzzy controller to balance the parameter update speed and noise resistance ability, solving the problems of data saturation and noise sensitivity of traditional methods; At the same time, the inertia identification accuracy is improved through a time-varying gain matrix and a covariance iteration mechanism, realizing full-scenario coverage from small noise-like disturbances to transient large disturbances, and providing high-precision and strong-robustness real-time data support for system frequency safety early warning and grid dispatching.

[0005] The technical solution adopted by the present invention is as follows:

[0006] A power system inertia identification method based on ARMAX and adaptive recursive optimization, comprising the following steps:

[0007] Step 1: Real-time collect the active power and frequency data signals at the bus, and preprocess the data signals to form a standard data set;

[0008] Step 2: Establish an active power-frequency dynamic response model for the new energy power system, and derive the frequency response transfer function of the system under large and small disturbance scenarios;

[0009] Step 3: Construct an ARMAX model, and establish the coupling relationship between the parameters of the ARMAX model and the inertia size;

[0010] Step 4: Use the adaptive recursive least squares algorithm to solve the parameters to be identified in the ARMAX model. Design a fuzzy controller to correct the adaptive forgetting factor online, and calculate the optimal parameter vector;

[0011] Step 5: Substitute the optimal parameter vector into the ARMAX model to calculate the inertia time constant of each unit.

[0012] In step 1, the active power and frequency data signals at the bus are collected in real time through PMU. The collected data signals are normalized based on the rated capacity of the unit and the rated frequency of the system respectively, and the signal mean value is removed to form a standard data set. Step 1 includes the following steps:

[0013] First, in order to eliminate the noise components introduced by random interference, the collected original active power and frequency data signals are processed by applying a Butterworth low-pass filter. The transfer function of the Butterworth low-pass filter is represented by T(s):

[0014]

[0015] In formula (1), s represents the Laplace operator; ω c represents the cut-off frequency; n represents the filter order.

[0016] This operation can effectively filter out the high-frequency noise in the data, making the subsequent analysis more accurate and reliable. This filtering process helps to retain the main signals in the data while removing the high-frequency noise that may affect the identification accuracy;

[0017] Secondly, after the data processed by the low-pass filter, per-unit value conversion is performed to convert the data processed by the low-pass filter into a proportional representation form P pu (t), f pu (t), as specifically implemented in formula (2):

[0018]

[0019] In formula (2), P pu (t) represents the per-unit value of active power; f pu (t) represents the per-unit value of frequency; t represents the time variable; P base (t) represents the actual value of active power; f base(t) represents the actual frequency value; S base represents the rated capacity; f base represents the rated frequency. This step helps in the standardization of data, making it independent of specific dimensions and facilitating subsequent analysis and comparison. Finally, based on the processed and transformed data signal and its difference from the corresponding mean value, as the state observation input signal for the identification method, it helps the algorithm to more accurately capture the changes in the system state, thereby improving the accuracy and reliability of identification. The specific formula is shown in Equation (3):

[0020]

[0021] In Equation (3), △P(t) is the power deviation signal; △f(t) is the frequency deviation signal; represents the power mean value, represents the frequency mean value, and M is the time window length.

[0022] In Step 2, during the inertial dynamic response of new energy power electronic devices, when a large power deficit disturbance occurs, the system will initiate a response process similar to that of a traditional power grid: the inertia support mechanism is activated first, and the synchronous units, in combination with new energy devices and energy storage systems, jointly provide inertial support to effectively suppress the rate of change of frequency (RoCoF); after the system frequency deviation exceeds the primary frequency regulation action threshold, the synchronous units, new energy virtual inertia modules, and energy storage droop control units will increase the active power output to prompt the frequency to quickly return to the steady-state operation range;

[0023] For power electronic interface devices, their inertia support ability has a linear relationship with the equivalent inertia time constant. Based on the modeling idea of the synchronous machine rotor motion equation, a first-order differential equation can be constructed to describe its dynamic process. For a single generator set i operating in the dynamic regulation state, its electromechanical transient characteristics can be characterized in incremental form, as shown in Equation (4):

[0024]

[0025] In Equation (4), H i represents the inertia time constant of generator i; △f i represents the rotor frequency deviation of generator i; △P im and △P ie respectively represent the mechanical power increment and electromagnetic power increment of generator i; D represents the damping coefficient.

[0026] During the period when the new energy power system is subjected to small noise-like disturbances, the primary frequency regulation function is in the dead zone and not activated. At this time, △P im = 0. By performing Laplace transform on Equation (4), the equivalent model transfer function of the active power - frequency dynamic characteristics under small disturbance excitation can be derived as shown in Equation (5):

[0027]

[0028] In Equation (5), G m (s) represents the equivalent transfer function of the active power - frequency dynamic characteristics of the generator under small - disturbance excitation; s is the Laplace operator; △f i (s) is the frequency - domain form of the change in the generator frequency; △P ie (s) is the frequency - domain form of the change in the power on the generator bus;

[0029] During a large - disturbance in the new - energy power system, the primary frequency regulation starts to adjust to suppress the frequency deviation. The open - loop transfer - function model of the frequency response under large - disturbance conditions is derived, and its mathematical expression is as shown in Equation (6):

[0030]

[0031] In Equation (6), G g (s) is the dynamic transfer function of the primary frequency - regulation governor.

[0032] In step 3, using the pre - processed active - power and frequency data signals as input / output, an ARMAX model is constructed, and the order of the ARMAX model is determined by the Akaike Information Criterion (AIC). Then, the ARMAX model is rewritten in the form of a transfer function to establish the coupling relationship between the ARMAX model parameters and the inertia size.

[0033] Specifically, it includes the following steps:

[0034] S3.1: When applying the ARMAX model to identify the inertia of the new - energy power system, it is necessary to preliminarily set the order of the model based on prior knowledge and structural - identification methods;

[0035] The present invention uses the Akaike Information Criterion (AIC) as the basis for optimizing the model order. This criterion optimizes the generalization ability of the model by introducing a complexity penalty factor, and its mathematical expression can be represented by AIC(η):

[0036]

[0037] In Equation (7): N represents the total number of data samples; η represents the estimated value of the order of the ARMAX model, represents the prediction - error variance of the η - order model.

[0038] The mathematical essence of the AIC criterion is reflected in the relative entropy approximation measure between the model-estimated probability distribution and the true distribution. This criterion realizes the optimal order selection by balancing the model parameter complexity and generalization ability. When the AIC value reaches the minimum, the corresponding model order η min is the order of the optimal identification model.

[0039] S3.2: Due to the continuous random noise disturbance in power grid operation, the traditional inertia identification method based on physical mechanism has the problem of limited identification accuracy because it ignores the nonlinear characteristics of the dynamic process and noise interference. At the same time, the system frequency deviation is highly correlated with the power disturbance data. Therefore, using the input / output-based system identification method can overcome the accuracy limitation of the traditional method and then construct a higher-precision ARMAX model. The general expression of the ARMAX model is shown in formula (8):

[0040]

[0041] In formula (8): A(q), B(q), C(q), D(q), and F(q) are all polynomial backward shift operators; y(t) represents the system output; u(t - n k ) represents the time-delay input; t represents time; n k represents the time-delay parameter; e(t) represents the zero-mean white noise disturbance;

[0042] The general mathematical form of A(q), B(q), C(q), D(q), and F(q) can be represented by R(q) as follows:

[0043]

[0044] In formula (9): r1, r2,..., r nr are polynomial coefficients; q represents the backward shift operator; n r represents the order of the polynomial R(q); R ∈ {A, B, C, D, F} corresponds to the coefficient r ∈ {a, b, c, d, f}, and the polynomial order is defined as n R , such as n A is the order of A(q). The system input u(t) and output y(t) are related by the time-delay parameter n k ; the typical value n k ∈ {0, 1}.

[0045] As a typical extended structure in the field of system identification, the core mechanism of the ARMAX model lies in that the system dynamic response and the noise transfer process share common poles. Specifically, when the parameters D(q) and F(q) in formula (8) are set to 1, the poles of both the system dynamic characteristics and the noise transfer function are constrained by the polynomial A(q). The discrete transfer function expression of the ARMAX model is shown in formula (10):

[0046]

[0047] In Equation (10): z is a complex variable in the discrete domain; A(z) is an autoregressive polynomial; B(z) is an exogenous input polynomial; C(z) is a moving average polynomial; u(t) represents the system input; e(t) represents a white noise disturbance with zero mean.

[0048] The determining link B(z) / A(z) characterizes the discrete transfer characteristics of the dynamic response of the generator set's power-frequency; the stochastic link C(z) / A(z) characterizes the dynamic action mechanism of the noise disturbance on the deterministic state components of the system. B(z) / A(z) is expressed as shown in Equation (11):

[0049]

[0050] In Equation (11): G(z) is a discrete transfer function; and are the parameters to be identified in the discrete transfer function model; n a is the order of the polynomial A(z); n b is the order of the polynomial B(z).

[0051] It should be particularly noted that in practical applications, the transfer function representation of the system to be identified is usually constructed based on the deterministic dynamic link B(z) / A(z) of the model.

[0052] Similarly, the stochastic link C(z) / A(z) is specifically expressed as shown in Equation (12):

[0053]

[0054] In Equation (12): are the moving average coefficients, reflecting the influence of the noise historical values on the current output; n c is the order of the polynomial C(z). In the new energy power system, the stochastic link separates the dynamic characteristics of the noise to avoid its mixing into the determining link, thereby improving the accuracy of the model in a noisy environment.

[0055] S3.3: Convert the discrete transfer function into a continuous-domain model through inverse Laplace transform, and establish an operator mapping based on the approximate relationship in Equation (13);

[0056] z = e Ts ≈ 1 + Ts (13);

[0057] In Equation (13): z is a complex variable in the discrete domain; e Ts is the mapping relationship from the discrete domain to the continuous domain; s is the Laplace operator; T is the sampling period.

[0058] The standard expression of the continuous-domain transfer function obtained after conversion is as shown in Formula (14):

[0059]

[0060] In Formula (14): G(s) represents the continuous-domain transfer function; a n , a n-2 , …, a0 and b n-1 , b n-2 , …, b0 are the parameters to be identified in the continuous-domain transfer function and need to be identified through system identification algorithms; s n , s n-1 , s n-2 represent the power terms of the complex frequency variable s in the continuous domain, and the highest order is n.

[0061] In Step 4, the adaptive recursive least squares algorithm is used to solve the parameters to be identified in the ARMAX model. In each correction period, the performance index and its change gradient are calculated using the residuals between the actual output and the predicted output; the adaptive forgetting factor is corrected online by designing a fuzzy controller, and the optimal parameter vector is obtained through iterative calculation.

[0062] Specifically, it includes the following steps:

[0063] S4.1: To accurately evaluate the equivalent inertia parameter of the system, a suitable identification algorithm needs to be selected to determine the parameters to be estimated in the ARMAX model. The traditional recursive least squares (RLS) method uses a fixed weight allocation mechanism in the parameter update process. As the identification process continues, it is prone to cause the over-accumulation effect of historical data, resulting in a decrease in the algorithm's convergence speed and inaccurate parameter estimation. Therefore, the present invention proposes to introduce a dynamic forgetting factor optimization strategy to enhance the algorithm's tracking ability for time-varying parameters by adjusting the historical data decay rate in real time, thereby improving the model identification accuracy under non-stationary working conditions. The construction of its optimization objective function is as shown in Formula (15):

[0064]

[0065] In Formula (15): J(θ) is the objective function of parameter estimation; y(i) is the actual output of the system at the i-th moment;

[0066] represents the information vector composed of historical input and output data, which is defined as:

[0067]

[0068] y(i - 1), y(i - 2), …, y(i - n a ) represent the system outputs at various past moments; na Represents the dynamic order of the autoregressive part; u(i - 1), u(i - 2), …, u(i - n b ) represents the system inputs at various past moments; n b Represents the dynamic order of the exogenous input part; Represents the vector space composed of parameters to be identified; The dynamic forgetting factor λ i Characterizes the exponential decay rate of the weight of historical data over time, and is a positive real number less than 1;

[0069] To minimize J(θ), take the first derivative of formula (15) and set it equal to 0, and verify that its second derivative is greater than 0. At this time, the estimated value of θ can be obtained As shown in formula (17):

[0070]

[0071] In formula (17): Represents the corrected objective function; Represents the information vector at time i; Represents the transpose of the information vector at time i; y(i) represents the system output at time i; t represents time.

[0072] S4.2: Introduce an adaptive forgetting factor improvement mechanism. The online recurrence formula of the ARMAX model parameters is as shown in formula (18):

[0073]

[0074] In formula (18): Represents the parameter estimation value at the previous moment; L(t) is a time-varying gain matrix, responsible for dynamically adjusting the parameter update rate; P(t) represents the covariance iteration matrix, used to quantify the statistical characteristics of the estimation error; P(t - 1) represents the covariance iteration matrix; Represents the information vector at time t; Represents the transpose of the information vector at time t; λ t Represents the forgetting factor, which is dynamically optimized by comprehensively evaluating the algorithm tracking performance and anti-noise robustness indicators;

[0075] S4.3: To improve the comprehensive efficiency of the algorithm, construct a fuzzy controller based on the performance index J and its change gradient △J to achieve adaptive optimization of the forgetting factor. Among them, the mathematical expressions of J(k) and △J(k) at the kth sampling moment are as shown in formula (19):

[0076]

[0077] In formula (19): J(k) represents the residual performance index at the kth sampling moment; Is the model prediction output;

[0078] is the residual between the actual output and the predicted output; △J(k) represents the change gradient of the residual performance index; J(k - m) represents the performance index m sampling periods ago; n characterizes the scale of residual sample screening; m corresponds to the discretization parameter of the time constant;

[0079] During the parameter update period, the preprocessing module parses the performance index J and its change gradient △J in real time based on Equation (19), and generates a dynamic compensation amount △λ through a fuzzy inference mechanism.

[0080] The fuzzy rule table is shown in Table 1. The fuzzy sets corresponding to the performance index and its change gradient are: negative large (PB), negative small (PS), zero (ZO), positive small (NS), positive large (NB); the fuzzy sets of the corresponding output dynamic compensation amount are: negative large (PB), negative medium (PM), negative small (PS), zero (ZO), positive small (NS), positive medium (NM), positive large (NB). The calculation of the fuzzy controller output △λ is based on Mamadani fuzzy inference. When both J and △J are NB, it indicates that the estimated error of the system equivalent inertia is extremely large and the error is still increasing rapidly. At this time, the output △λ is PB, greatly reducing the forgetting factor; when both J and △J are PB, it indicates that the estimated result of the system equivalent inertia is very close to the true value and the error remains small. At this time, the output △λ is NB, and the adaptive forgetting factor is set to the maximum value, increasing the correction effect of new data on the inertia estimation result, thereby realizing the online adaptive adjustment of the forgetting factor.

[0081] Table 1 Fuzzy rule table

[0082]

[0083] S4.4: Further, the efficient online estimation of the parameter vector θ can be realized based on the adaptive recursive least squares algorithm through Equation (18). Formulas (15) - (19) are the implementation process of the adaptive recursive least squares algorithm.

[0084] In step 5, the optimal parameter vector is substituted into the ARMAX model, and the model is transformed into the form of a continuous transfer function based on the forward Euler transform, and then a unit step excitation signal is input, and the inertia time constant of each unit is calculated by analyzing the initial segment slope.

[0085] Specifically: After substituting the optimal parameters obtained by the identification algorithm in step 3 into the ARMAX model, the conversion of the transfer function from discrete to continuous is realized based on the Euler approximation method of formula (8). Subsequently, a unit step excitation signal is input, and the derivative value at the initial moment of the step response curve is extracted, and then the inertia of each unit is calculated. The calculation formula of the inertia of each unit is shown in formula (20):

[0086]

[0087] In Equation (20), H is the inertia of each unit; △P = 1 p.u. is the input unit step power disturbance signal; is the rate of change of frequency deviation at the initial moment of the step response curve, that is, the slope of the initial segment; t is time.

[0088] The mathematical representation of the equivalent inertia of the power system can be defined as the algebraic sum of the inertia contributions of synchronous units, and its physical essence is reflected as the vector sum superposition of the inertia support capabilities of each unit. The specific mathematical expression is shown in Equation (21):

[0089]

[0090] In Equation (21), E sys represents the equivalent inertia of the power system; N i represents the number of synchronous units participating in inertia response in the system. The inertia time constant and its rated capacity of the i-th device are denoted as H G,i and S G,i ; M j represents the number of renewable new energy units with virtual inertia support capabilities in the system. The equivalent virtual inertia constant and its rated capacity of the j-th device are corresponding to H V,j and S V,j .

[0091] Based on the theoretical framework of the equivalent inertia time constant, the present invention uses H sys to represent the overall inertia level of the system, and its mathematical representation is shown in Equation (22):

[0092]

[0093] After obtaining the inertia of each unit, finally, the global inertia parameter integration is carried out according to Equation (22) to complete the dynamic identification of the equivalent inertia of the primary system.

[0094] A power system inertia identification method based on ARMAX and adaptive recursive optimization of the present invention has the following technical effects:

[0095] 1) By identifying the equivalent inertia parameters of the new energy power system in real time and accurately, the present invention provides key data support for frequency safety analysis, virtual inertia control strategy optimization, and emergency frequency regulation, and effectively improves the frequency stability and anti-disturbance ability of the high-proportion new energy power grid.

[0096] 2) By constructing an ARMAX model, the present invention solves the problem that the traditional mechanism model ignores non-linear dynamics and noise interference; combined with the Akaike information criterion (AIC) to optimize the model order, the generality and robustness of the model in the scenarios of quasi-steady-state noise disturbance and transient large disturbance are significantly improved.

[0097] 3) Aiming at the historical data saturation and noise sensitivity problems of traditional recursive least squares (RLS) method in time-varying parameter tracking, the present invention proposes a dynamic forgetting factor optimization strategy. Through the fuzzy logic regulator, the residual performance index and its changing gradient are analyzed in real time, and the forgetting factor weight is adaptively adjusted, thereby achieving the coordinated optimization of the algorithm convergence speed and noise robustness.

[0098] 4) The present invention achieves full-scenario coverage from quasi-steady-state noise disturbances to transient large disturbance conditions, solving the problem of insufficient identification stability of traditional methods when the sampling window length and load mutation intensity change, and provides reliable technical support for real-time inertia perception and frequency safety warning of new energy power systems. BRIEF DESCRIPTION OF THE DRAWINGS

[0099] The present invention will be further described below in conjunction with the accompanying drawings and embodiments:

[0100] Figure 1 This is a flowchart of equivalent inertia identification of the new power system described in the present invention.

[0101] Figure 2 This is the topology diagram of the IEEE10 machine 39-node simulation example system used in the present invention.

[0102] Figure 3 This is a comparison chart of the measured and estimated frequency deviation of generator G10 in the simulation example.

[0103] Figure 4 It is a comparison chart of the inertia response fitting results of the generator G10 in the simulation example under various comparison models.

[0104] Figure 5 This is a comparison chart of the identification results errors of each comparison model under different disturbance levels.

[0105] Figure 6 This is a comparison chart of the identification result errors of each comparison model under different time window lengths. DETAILED DESCRIPTION

[0106] A new power system inertia identification method based on ARMAX and adaptive recursive optimization includes the following steps:

[0107] Step 1: First, in order to eliminate the noise components introduced by random interference, the collected original data is first processed by applying a Butterworth low-pass filter. This operation can effectively filter out the high-frequency noise in the data, making the subsequent analysis more accurate and reliable. This filtering process helps to retain the main signals in the data while removing the high-frequency noise that may affect the identification accuracy. Secondly, the data processed by the low-pass filter then needs to be subjected to per-unit conversion, converting the data into a proportional representation relative to a certain reference value. This step helps with the standardization of the data, making it independent of the specific dimension and facilitating subsequent analysis and comparison. Finally, based on the processed and converted data information and its difference from the corresponding mean value, it is input as the state observation information of the identification algorithm. This helps the algorithm to more accurately capture the changes in the system state, thereby improving the accuracy and reliability of the identification.

[0108] Step 2: In the inertial dynamic response of new energy power electronic devices, when a large power deficit disturbance occurs, the system will initiate a response process similar to that of a traditional power grid: the inertia support mechanism is activated first, and the synchronous units, in conjunction with new energy devices and energy storage systems, jointly provide inertial support to effectively suppress the rate of change of frequency (RoCoF); after the system frequency deviation exceeds the primary frequency regulation action threshold, the synchronous units, new energy virtual inertia modules, and energy storage droop control units will increase the active power output to prompt the frequency to quickly return to the steady-state operation range.

[0109] For a power electronic interface device, its inertia support ability has a linear relationship with the equivalent inertia time constant. Based on the modeling idea of the synchronous machine rotor motion equation, a first-order differential equation can be constructed to describe its dynamic process. For a single generator set i operating in the dynamic regulation state, its electromechanical transient characteristics can be characterized in incremental form, as shown in Equation (4):

[0110]

[0111] In the formula, H i represents the inertia time constant of generator i, △f i represents the rotor frequency deviation of generator i, △P im , △P ie represent the mechanical power increment and electromagnetic power increment of generator i respectively, and D represents the damping coefficient.

[0112] During the period when the system is subjected to small disturbances similar to noise, the primary frequency regulation function is in the dead zone and not activated. At this time, △P im = 0. By performing a Laplace transform on Equation (4), the equivalent model transfer function of the active power-frequency dynamic characteristics under small disturbance excitation can be derived, as shown in Equation (5):

[0113]

[0114] where s is the Laplace operator; △f i (s) is the frequency-domain form of the change in generator frequency; △P ie (s) is the frequency-domain form of the change in power on the generator busbar.

[0115] During a large disturbance in the system, primary frequency regulation starts to adjust to suppress frequency deviation. The open-loop transfer function model of the frequency response under large disturbance conditions is derived, and its mathematical expression is shown in Equation (6):

[0116]

[0117] where G g (s) is the dynamic transfer function of the primary frequency regulator. In addition, the inertia response mechanism of new energy power electronic devices is mathematically equivalent to that of synchronous units at the power-frequency dynamic characteristic level. Based on this characteristic, during the inertia support process, the dynamic characteristics of power electronic devices can be simulated by the equivalent modeling method of synchronous machines.

[0118] Step 3: In the power grid operation, there is continuous random noise disturbance. The traditional inertia identification method based on physical mechanism has the problem of limited identification accuracy due to ignoring the nonlinear characteristics of the dynamic process and noise interference. At the same time, the system frequency deviation is highly correlated with the power disturbance data. Therefore, using the system identification method based on input / output can overcome the accuracy limitation of the traditional method and then construct a higher-precision frequency response model. The model expression is shown in Equation (8):

[0119]

[0120] where the general mathematical forms of A(q), B(q), C(q), D(q), and F(q) are: where R ∈ {A, B, C, D, F}, the corresponding coefficient r ∈ {a, b, c, d, f}, and q is the backward shift operator. The polynomial order is defined as n R (such as n A is the order of A(q)). The system input u(t) and output y(t) are related by the time-delay parameter n k (typical value n k ∈ {0, 1}), and e(t) represents white noise disturbance with zero mean.

[0121] When applying the ARMAX model to identify the inertia of a new energy power system, it is necessary to initially set the model order based on prior knowledge and the structure identification method. The present invention uses the Akaike Information Criterion (AIC) as the basis for optimizing the model order. This criterion optimizes the model generalization ability by introducing a complexity penalty factor, and its mathematical expression is shown in Equation (7):

[0122]

[0123] where η represents the estimated order of the ARMAX model, and represents the prediction error variance of the η-order model. The mathematical essence of the AIC criterion is reflected in the approximate measure of the relative entropy between the estimated probability distribution of the model and the true distribution. This criterion realizes the optimal order selection by balancing the model parameter complexity and the generalization ability. When the AIC value reaches the minimum, the corresponding model order η min is the order of the optimal identification model.

[0124] As a typical extended structure in the field of system identification, the core mechanism of the ARMAX model lies in the sharing of common poles in the system dynamic response and the noise transfer process. Specifically, when the parameters D(q) and F(q) in formula (8) are set to 1, the poles of both the system dynamic characteristics and the noise transfer function are constrained by the polynomial A(q). The discrete transfer function expression of the ARMAX model is shown in formula (10):

[0125]

[0126] where z is the complex variable in the discrete domain; u(t) is the excitation signal (power deviation); y(t) is the response observation value (frequency deviation); e(t) is the random disturbance term; the deterministic link B(z) / A(z) represents the discrete transfer characteristics of the power-frequency dynamic response of the generator set; the random link C(z) / A(z) represents the dynamic action mechanism of the noise disturbance on the deterministic state components of the system. The expression of B(z) / A(z) is shown in formula (11):

[0127]

[0128] where and are the parameters to be identified in the discrete transfer function model. It should be particularly noted that in practical applications, the transfer function representation of the system to be identified is usually constructed based on the deterministic dynamic link of the model.

[0129] When converting the discrete transfer function to a continuous-domain model through the inverse Laplace transform, an operator mapping needs to be established based on the approximate relationship in formula (13):

[0130] z = e Ts ≈ 1 + Ts (13) (1

[0131] where T is the sampling period. The standard expression of the continuous-domain transfer function obtained after conversion is shown in formula (14):

[0132]

[0133] where a n and a n-2 , …, a0 and b n-1 and b n-2 , …, b0 are the parameters to be identified in the continuous-domain transfer function and need to be identified through a system identification algorithm.

[0134] Step 4: To accurately evaluate the system equivalent inertia parameter, a suitable identification algorithm needs to be selected to determine the parameters to be estimated in the ARMAX model. The traditional Recursive Least Squares (RLS) uses a fixed weight allocation mechanism in the parameter update process. As the identification process continues, it is prone to cause the over-accumulation effect of historical data, resulting in a decrease in the algorithm convergence rate and inaccurate parameter estimation. Therefore, the present invention proposes to introduce a dynamic forgetting factor optimization strategy to enhance the algorithm's tracking ability for time-varying parameters by adjusting the historical data decay rate in real time, thereby improving the model identification accuracy under non-stationary conditions. The optimized objective function is constructed as shown in formula (15):

[0135]

[0136] where the forgetting factor λ i represents the exponential decay rate of the historical data weight over time and is a positive real number less than 1; is the information vector composed of the system input and output data at time t; represents the vector space composed of the parameters to be identified.

[0137] To minimize J(θ), take the first derivative of formula (15) and set it equal to 0, and check that its second derivative is greater than 0. At this time, the estimated value of θ can be obtained as shown in formula (17):

[0138]

[0139] By introducing the dynamic forgetting factor improvement mechanism, the online recurrence formula of the ARMAX model parameters is as shown in formula (18):

[0140]

[0141] where L(t) is the time-varying gain matrix responsible for dynamically adjusting the parameter update rate; P(t) represents the covariance iteration matrix used to quantify the statistical characteristics of the estimation error. The forgetting factor λ is dynamically optimized by comprehensively evaluating the algorithm tracking performance and anti-noise robustness index.

[0142] To improve the comprehensive efficiency of the algorithm, a fuzzy logic regulator based on the performance index J and its change gradient △J is constructed to realize the adaptive optimization of the forgetting factor. Among them, the mathematical expressions of J(k) and △J(k) at the kth sampling moment are as shown in formula (19):

[0143]

[0144] In the formula, n represents the scale of residual sample screening; m corresponds to the discretization parameter of the time constant.

[0145] During the parameter update period, the preprocessing module parses the performance index J and its change gradient ΔJ in real time based on Equation (19), generates a dynamic compensation amount Δλ through a fuzzy inference mechanism, and realizes the online adaptive adjustment of the forgetting factor. Further, an efficient online estimation of the parameter vector θ is realized based on the adaptive recursive least squares algorithm.

[0146] Step Five: After substituting the optimal parameters obtained by the identification algorithm in Step Three into the identification model, the transfer function is converted from discrete to continuous based on the Euler approximation method of Equation (8). Subsequently, a unit step excitation signal is input, and by extracting the derivative value of the step response curve at the initial moment, the inertia of each unit is calculated.

[0147] The mathematical representation of the equivalent inertia of the power system can be defined as the algebraic sum of the inertia contributions of synchronous units. Its physical essence is the vector sum superposition of the inertia support capabilities of each unit, and the specific mathematical expression is shown in Equation (21):

[0148]

[0149] In the formula, N i represents the number of synchronous units participating in inertial response in the system, where the inertia time constant and its rated capacity of the i-th device are denoted as H G,i and S G,i ; M j represents the number of renewable new energy units with virtual inertia support capabilities in the system, where the equivalent virtual inertia constant and its rated capacity of the j-th device are correspondingly H V,j and S V,j .

[0150] Based on the theoretical framework of the equivalent inertia time constant, the present invention uses H sys to represent the overall inertia level of the system, and its mathematical representation is shown in Equation (22):

[0151]

[0152] After obtaining the inertia of each unit, finally, the global inertia parameter integration is carried out according to Equation (22) to complete the dynamic identification of the equivalent inertia of the primary system. The flow chart of a new power system inertia identification method based on ARMAX and adaptive recursive optimization of the present invention is as Figure 1 shown.

[0153] Step 6: The above-mentioned novel power system inertia identification method is verified for its accuracy through simulation examples.

[0154] To multi-dimensionally test the applicability of the present invention, an IEEE 10-machine 39-bus simulation test system is constructed based on the MATLAB / Simulink platform. The topological structure of the simulation system is as Figure 2 shown. The configuration of the simulation system is as follows: The No. 5, No. 6, and No. 9 synchronous generator sets of the original system are replaced with doubly-fed wind turbine generators (connected to the grid through power electronic converters), and the new energy penetration rate reaches 32.2%. Among them, the wind turbine generators are equipped with grid-following virtual inertia control. The simulation benchmark parameters are set as: the system rated frequency is 50 Hz, and the capacity reference value is 100 MVA.

[0155] To verify the measurement accuracy of the present invention, generator G10 is selected as the test object. The active-frequency dynamic response data of the unit tie line port is monitored through the PMU device, and data preprocessing is implemented. The dynamic characteristics matching analysis is carried out between the measured value of the unit terminal frequency and the frequency identification value under the active excitation of the model. The comparison result of the actual value and the identification value of the frequency dynamic deviation is as Figure 3 shown. As Figure 3 can be seen, the curve of the frequency deviation estimation value obtained by using the ARLS-ARMAX identification method of the present invention has a high degree of coincidence with the measured value curve in terms of waveform, verifying that the model has a high fitting ability for the actual output, indicating that it is applicable to the high-precision identification of the system equivalent inertia.

[0156] Through simulation, the adaptation situations of three main system identification methods (controlled autoregressive (CAR) model, transfer function estimation (TFEST), and state space (SSEST) model) in the current power system inertia parameter identification field are compared under different disturbance intensities and sampling time windows.

[0157] The present invention uses the normalized root mean square error to evaluate the fitting degree between the estimated value and the actual value identified by the model, and its mathematical expression is as shown in formula (23):

[0158]

[0159] In the formula, F represents the model fitting degree; y i is the actual value; is the mean value of the actual value; is the estimated value identified by the model;

[0160]

[0161] By comparing the system equivalent inertia identification result with the actual inertia value, the inertia identification deviation ε H can be calculated, and its mathematical expression is as shown in formula (24):

[0162]

[0163] In the formula, H i is the actual inertia value; is the estimated value of the system equivalent inertia.

[0164] In the disturbance configuration, the load frequency fluctuations of the power system under normal conditions are simulated by small noise-like disturbances, while the fault scenarios are simulated by large disturbance conditions. During the experiment, a load disturbance is applied to the system at t = 10 s, and the simulation sampling time interval is Δt = 5 s. In order to explore the adaptability of the model proposed in the present invention under different disturbance magnitudes, the present invention designs the following three analysis scenarios:

[0165] Scenario 1: Select bus nodes such as No. 3, 8, 12, 16, 20, 24, 25, and 39 to configure random noise-like load small disturbances, and construct a load fluctuation test scenario with noise-like characteristics under the quasi-steady state condition of the power grid.

[0166] In Scenario 1, the estimated results of the equivalent inertia of each unit are shown in Table 2.

[0167] Table 2 Estimated results of the equivalent inertia of each unit under four identification methods

[0168]

[0169] The experimental comparison results in Table 2 show that the ARLS-ARMAX model proposed in the present invention has obvious advantages in the inertia identification accuracy, and its error range remains between 0.22% and 3.79%. Under small disturbance conditions, the model can effectively handle the system input noise problem and obtain more accurate identification results. Compared with other methods: although the overall error of the SSEST model is small, its identification effect decreases under complex disturbance conditions; the CAR model has weak adaptability to non-linear disturbances; due to the limitation of the frequency domain analysis method, the maximum error of the TFEST model reaches 6.50%. In addition, affected by the virtual inertia control link, the identification errors of the equivalent inertia of the wind farm (3.1%, 3.79% and 3.15%) are slightly higher than those of the synchronous units. However, compared with the other three identification methods, the identification results of the method proposed in the present invention are basically consistent with the inertia reference value, verifying the effectiveness of the method in identifying the system inertia under the quasi-steady state operation condition of the new energy power system.

[0170] Scenario 2: Ignore the system noise-like disturbance. When t = 10 s, set a sudden increase of 10% active power load at bus 39, and remove the fault after 1 s.

[0171] In Scenario 2, select the synchronous generator G10 as the research object. For the power system equivalent inertia response model constructed based on four different system inertia identification algorithms, the degree of fitting of its dynamic characteristics can be passed throughFigure 4 Comparative analysis of the shown frequency deviation curves. Analysis Figure 4 It can be seen that the overall fitting effect of the ARLS-ARMAX model is significantly better than the other three methods, and the fitting accuracy reaches 98.97%. This indicates that the model can accurately reflect the dynamic process of system frequency regulation, which is closely related to its model structure characteristics. The model proposed in this invention realizes the dynamic parameter update mechanism by integrating the recursive least squares algorithm, and can track the non-linear time-varying characteristics of the system during load mutation in real time. In contrast, the CAR model is limited by the linear autoregressive structure of the fixed order and has insufficient adaptability at the moment of load mutation; although the SSEST model constructs the state space equation based on the subspace identification principle, it is highly sensitive to the disturbances under the non-stationary conditions of the power system; while the transfer function parameterization method of the TFEST model has a phase tracking mismatch problem when describing non-linear dynamic characteristics. It can be seen from the local enlarged figure that at the fault removal moment (t = 11 s), the identification model proposed in this invention still maintains a high fitting degree, further proving its superior performance in dynamic response and verifying its robustness and accuracy under complex working conditions.

[0172] Scenario 3: When t = 10 s, 6%, 9%, 12%, and 15% of the load are sequentially increased at bus 23 and removed 1 s after each increase to simulate the transient large disturbance in the actual power system.

[0173] In Scenario 3, Figure 5 shows the error of the system equivalent inertia identification results obtained by using 4 identification methods under different degrees of large disturbances. From Figure 5 the analysis, it can be seen that under large disturbances of different levels (6% - 15% load mutation), the maximum error of the four identification models under large disturbance conditions does not exceed 7.5%, and all remain within an acceptable accuracy range, with the credibility of the basic measurements. However, specific analysis shows that the error fluctuation range of the TFEST model is the largest, and the difference between the minimum error and the maximum error reaches 3% under different load fluctuation ratios. Although the error fluctuation of the CAR model is the smallest among the four comparison models, the average error is relatively high, and the error distribution under different load fluctuation conditions is about 5%. In contrast, the measurement error and fluctuation range of this invention are both the smallest, showing better identification stability.

[0174] Therefore, compared with the other three algorithms, the ARLS-ARMAX identification model has a wider applicability. Its core advantage lies in the collaborative optimization of the model structure and parameter update mechanism, realizing full-scenario coverage from quasi-steady state noise disturbances to transient large disturbance conditions.

[0175] Based on the operating conditions set in Scenario 1, by changing the duration of the sampling time window, the model identification errors of four identification algorithms in the inertia response stage under different window durations are quantified. Furthermore, the performance differences of the four different identification methods under sampling time windows of different lengths can be studied. The identification result errors are as Figure 6 shown. From Figure 6 the comparison results, it can be seen that the inertia identification result errors of the proposed ARLS-ARMAX model in the present invention are all within 4% under different sampling window lengths, and the degree of fluctuation is small, showing strong adaptability and stability. In contrast, the identification result errors of the TFFEST model almost reach three times that of the ARLS-ARMAX model, indicating its poor adaptability when dealing with different sampling window lengths.

[0176] A novel power system inertia identification method based on ARMAX and adaptive recursive optimization in the present invention solves the problems of data saturation and noise sensitivity of traditional methods under time-varying operating conditions, and realizes real-time identification of the inertia of a novel power system. The present invention combines the ARMAX model with the adaptive recursive optimization algorithm, overcomes the problem of insufficient accuracy of traditional identification methods in new energy power systems due to noise interference, nonlinear dynamics, and complex disturbance scenarios, significantly improves the robustness and adaptability of inertia identification, and provides reliable technical support for power system frequency safety analysis and real-time warning.

Claims

1. A power system inertia identification method based on ARMAX and adaptive recursive optimization, characterized in that It includes the following steps: Step 1: Real-time collect the active power and frequency data signals at the bus, and preprocess the data signals to form a standard data set; Step 2: Establish an active power-frequency dynamic response model for the new energy power system, and deduce the frequency response transfer function of the system under large and small disturbance scenarios; Step 3: Construct an ARMAX model, and establish the coupling relationship between the ARMAX model parameters and the inertia size; Step 4: Use the adaptive recursive least squares algorithm to solve the parameters to be identified in the ARMAX model; Design a fuzzy controller to online correct the adaptive forgetting factor, and calculate the optimal parameter vector; Step 5: Substitute the optimal parameter vector into the ARMAX model to calculate the inertia time constant of each unit.

2. The method for identifying the inertia of a power system based on ARMAX and adaptive recursive optimization according to claim 1, characterized in that: In the said Step 1, the active power and frequency data signals at the bus are real-time collected by the PMU, and the collected data signals are normalized based on the rated capacity of the unit and the rated frequency of the system respectively, and the signal mean value is removed to form a standard data set.

3. The power system inertia identification method based on ARMAX and adaptive recursive optimization according to claim 2, characterized in that: Step 1 includes the following steps: First, apply a Butterworth low-pass filter to the collected original active power and frequency data signals for processing, and the transfer function of the Butterworth low-pass filter is represented by T(s): In Equation (1), s represents the Laplace operator; ω c represents the cut-off frequency; n represents the filter order; Secondly, after the data is processed by low-pass filtering, per-unit conversion is performed to convert the data processed by low-pass filtering into a proportional representation form P pu (t), f pu (t), and the specific implementation is shown in Equation (2): In Equation (2), P pu (t) represents the per-unit value of active power; f pu (t) represents the per-unit value of frequency; t represents the time variable; P base (t) represents the actual value of active power; f base (t) represents the actual value of frequency; S base represents the rated capacity; f base represents the rated frequency; finally, based on the processed and transformed data signal and its difference from the corresponding mean value, as the state observation input signal of the identification method, the specific formula is shown in Equation (3): In Equation (3), △P(t) is the power deviation signal; △f(t) is the frequency deviation signal; represents the power mean value, represents the frequency mean value, and M is the time window length.

4. The method for identifying the inertia of a power system based on ARMAX and adaptive recursive optimization according to claim 1, characterized in that: In the said Step 2, for the power electronic interface device, its dynamic process is described by constructing a first-order differential equation; For a single generator set i operating in the dynamic regulation state, its electromechanical transient characteristics are characterized in incremental form, as shown in Equation (4): In Equation (4), H i represents the inertia time constant of generator i; △f i represents the rotor frequency deviation of generator i; △P im and △P ie respectively represent the mechanical power increment and electromagnetic power increment of generator i; D represents the damping coefficient; During the period when the new energy power system is subjected to noise-like small disturbances, the primary frequency regulation function is in the dead zone and not activated. At this time, △P im = 0. The Laplace transform is performed on formula (4), and the equivalent model transfer function of the active power-frequency dynamic characteristics under small disturbance excitation is derived as shown in formula (5): In formula (5), G m (s) represents the equivalent transfer function of the active - frequency dynamic characteristics of the generator under small - disturbance excitation; s is the Laplace operator; △f i (s) is the frequency - domain form of the change in the generator frequency; △P ie (s) is the frequency - domain form of the change in the power on the generator busbar; During a large disturbance in the new energy power system, the primary frequency regulation starts to adjust the action to suppress the frequency deviation, and the open-loop transfer function model of the frequency response under the large disturbance condition is deduced, and its mathematical expression is shown in Formula (6): In Equation (6), G g (s) is the dynamic transfer function of the primary frequency regulation governor.

5. The power system inertia identification method based on ARMAX and adaptive recursive optimization according to claim 1, characterized in that: In the said Step 3, using the preprocessed active power and frequency data signals as input / output, construct an ARMAX model, and determine the order of the ARMAX model by the Akaike information criterion AIC, and then rewrite the ARMAX model into the transfer function form to establish the coupling relationship between the ARMAX model parameters and the inertia size.

6. The power system inertia identification method based on ARMAX and adaptive recursive optimization according to claim 1, characterized in that: Step 3 includes the following steps: S3.1: When applying the ARMAX model to carry out inertia identification of the new energy power system, based on prior knowledge and the structure identification method, the preliminary setting of the model order is realized; Use the Akaike information criterion AIC as the basis for optimizing the model order. This criterion realizes the optimization of the model generalization ability by introducing a complexity penalty factor, and its mathematical expression is represented by AIC(η): In Equation (7): N represents the total number of data samples; η represents the estimated order of the ARMAX model, which represents the prediction error variance of the η-order model; When the AIC value reaches the minimum, the corresponding model order η min is the order of the optimal identification model; S3.2: Construct an ARMAX model; the general expression of the ARMAX model is shown in Formula (8): In Equation (8): A(q), B(q), C(q), D(q), and F(q) are all polynomial backward shift operators; y(t) represents the system output; u(t - n k ) represents the time-delay input; t represents time; n k represents the time-delay parameter; e(t) characterizes the zero-mean white noise disturbance; The general mathematical forms of A(q), B(q), C(q), D(q), and F(q) can be represented by R(q) as: In formula (9): is the polynomial coefficient; q represents the backward shift operator; n r represents the order of the polynomial R(q); the polynomial order is defined as n R ; the system input u(t) and the output y(t) are related by the time-delay parameter n k ; the typical value of n k ∈{0, 1}; When the parameter D(q) and F(q) in Equation (8) are set to 1, the poles of the system dynamic characteristics and the noise transfer function are both constrained by the polynomial A(q), and the discrete transfer function expression of the ARMAX model is shown in Formula (10): In Equation (10): z is a complex variable in the discrete domain; A(z) is an autoregressive polynomial; B(z) is an exogenous input polynomial; C(z) is a moving average polynomial; u(t) represents the system input; e(t) represents a white noise disturbance with zero mean. The determination link B(z) / A(z) characterizes the discrete transfer characteristics of the dynamic response of the generator set's power-frequency; the stochastic link C(z) / A(z) characterizes the dynamic action mechanism of the noise disturbance on the deterministic state components of the system; B(z) / A(z) is expressed as shown in Equation (11): In Equation (11): G(z) is a discrete transfer function; and are the parameters to be identified for the discrete transfer function model; n a is the order of the polynomial A(z); n b is the order of the polynomial B(z); Similarly, the stochastic link C(z) / A(z) is specifically expressed as shown in Equation (12): In formula (12): is the moving average coefficient, reflecting the influence of the historical value of noise on the current output; n c is the order of the polynomial C(z); S3.3: Convert the discrete transfer function into a continuous-domain model through inverse Laplace transform, and establish an operator mapping based on the approximate relationship in Equation (13); z = e Ts ≈ 1 + Ts(13); In Equation (13): z is a complex variable in the discrete domain; e Ts is the mapping relationship from the discrete domain to the continuous domain; s is the Laplace operator; T is the sampling period; The standard expression of the continuous-domain transfer function obtained after conversion is as shown in Equation (14): In Equation (14): G(s) represents the continuous-domain transfer function; a n , a n-2 , …, a0 and b n-1 , b n-2 , …, b0 are the parameters to be identified for the continuous-domain transfer function and need to be identified through a system identification algorithm; s n , s n-1 , s n-2 represent the power terms of the complex frequency variable s in the continuous domain, and the highest order is n.

7. The power system inertia identification method based on ARMAX and adaptive recursive optimization according to claim 1, wherein: In Step 4, an adaptive recursive least squares algorithm is used to solve the parameters to be identified in the ARMAX model; within each correction period, the performance index and its change gradient are calculated using the residual between the actual output and the predicted output; the adaptive forgetting factor is corrected online by designing a fuzzy controller, and the optimal parameter vector is obtained through iterative calculation.

8. The method for identifying the inertia of a power system based on ARMAX and adaptive recursive optimization according to claim 6, wherein: Step 4 includes the following steps: S4.1: Introduce a dynamic forgetting factor optimization strategy to improve the model identification accuracy under non-stationary conditions. The optimization objective function is constructed as shown in Equation (15): In Equation (15): J(θ) is the objective function of parameter estimation; y(i) is the actual output of the system at the i-th moment; Denote the information vector composed of historical input and output data, which is defined as: y(i - 1), y(i - 2), …, y(i - n a ) represent the system outputs at past various moments; n a represents the dynamic order of the autoregressive part; u(i - 1), u(i - 2), …, u(i - n b ) represent the system inputs at past various moments; n b represents the dynamic order of the exogenous input part; represents the vector space composed of the parameters to be identified; the dynamic forgetting factor λ i characterizes the exponential decay rate of the historical data weight over time and is a positive real number less than 1; To minimize J(θ), take the first derivative of Equation (15) and set it equal to 0, and verify that its second derivative is greater than 0. At this point, the estimated value of θ can be obtained. As shown in Equation (17): In Equation (17): represents the corrected objective function; represents the information vector at time i; represents the transpose of the information vector at time i; y(i) represents the system output at time i; t represents time; S4.2: Introduce an adaptive forgetting factor improvement mechanism. The online recurrence formula of the parameters of the ARMAX model is as shown in Equation (18): In formula (18): represents the parameter estimation value at the previous moment; L(t) is a time-varying gain matrix responsible for dynamically adjusting the parameter update rate; P(t) represents the covariance iteration matrix used to quantify the statistical characteristics of the estimation error; P(t - 1) represents the covariance iteration matrix; represents the information vector at time t; represents the transpose of the information vector at time t; λ t represents the forgetting factor; S4.3: To improve the comprehensive efficiency of the algorithm, a fuzzy controller based on the performance index J and its change gradient △J is constructed to realize the adaptive optimization of the forgetting factor; among them, the mathematical expressions of J(k) and △J(k) at the k-th sampling moment are as shown in Equation (19): In formula (19): J(k) represents the residual performance index at the k-th sampling moment; is the model prediction output; is the residual between the actual output and the predicted output; △J(k) represents the change gradient of the residual performance index; J(k - m) represents the performance index m sampling periods ago; n characterizes the scale of residual sample screening; m corresponds to the time constant discretization parameter.

9. The power system inertia identification method based on ARMAX and adaptive recursive optimization according to claim 8, characterized in that: During the parameter update period, the preprocessing module analyzes the performance index J and its change gradient △J in real time based on Equation (19), and generates a dynamic compensation amount △λ through a fuzzy inference mechanism; The fuzzy sets corresponding to the performance index and its change gradient are: negative large (PB), negative small (PS), zero (ZO), positive small (NS), positive large (NB); the fuzzy sets of the corresponding output dynamic compensation amount are: negative large (PB), negative medium (PM), negative small (PS), zero (ZO), positive small (NS), positive medium (NM), positive large (NB); the calculation of the fuzzy controller output △λ is based on Mamadani fuzzy inference. If both J and △J are NB, it indicates that the estimated error of the system equivalent inertia is extremely large and the error is still increasing rapidly. At this time, the output △λ is PB, and the forgetting factor is greatly reduced; if both J and △J are PB, it indicates that the estimated result of the system equivalent inertia is very close to the true value and the error remains small. At this time, the output △λ is NB, and the adaptive forgetting factor is set to the maximum value to increase the correction effect of new data on the inertia estimation result, thereby realizing the online adaptive adjustment of the forgetting factor.

10. The power system inertia identification method based on ARMAX and adaptive recursive optimization according to claim 9, characterized in that: In step 5, the optimal parameter vector is substituted into the ARMAX model, and the model is transformed into the form of a continuous transfer function based on the forward Euler transform. Then, a unit step excitation signal is input, and the inertia time constant of each unit is calculated by analyzing the slope of the initial segment as follows: After substituting the optimal parameters obtained by the identification algorithm in step 3 into the ARMAX model, the conversion of the transfer function from discrete to continuous is realized based on the Euler approximation method of formula (8). Subsequently, a unit step excitation signal is input, and the derivative value of the step response curve at the initial moment is extracted, and then the inertia of each unit is calculated. The calculation formula of the inertia of each unit is shown in formula (20): In Equation (20), H is the inertia of each unit; △P = 1 p.u. is the input unit step power disturbance signal; is the rate of change of frequency deviation at the initial moment of the step response curve, that is, the slope of the initial segment; t is time; The mathematical representation of the equivalent inertia of the power system can be defined as the algebraic sum of the inertia contributions of the synchronous units, which is reflected as the vector sum superposition of the inertia support capabilities of each unit. The specific mathematical expression is shown in formula (21): In formula (21), E sys represents the equivalent inertia of the power system; N i represents the number of synchronous units participating in inertial response within the system. The inertia time constant and rated capacity of the i-th device are denoted as H G,i and S G,i ; M j represents the number of renewable energy units with virtual inertia support capacity within the system. The equivalent virtual inertia constant and rated capacity of the j-th device are correspondingly H V,j and S V,j ; Based on the theory of equivalent inertia time constant, use H sys to represent the overall inertia level of the system, and its mathematical representation is shown in formula (22): After obtaining the inertia of each unit, finally, the global inertia parameter integration is carried out according to formula (22) to complete the dynamic identification of the equivalent inertia of the primary system.

Images