A method, system, device and medium for node inertia identification of a power system

Abstract

The application discloses a kind of node inertia identification method, system, equipment and medium of power system, it is related to power system inertia identification technical field, including after detecting load step disturbance, the active power and frequency data of target node are preprocessed, and fixed disturbance window before and adjustable disturbance window after are respectively set, power and frequency sequence are intercepted, dynamic parameterization model is established and discrete domain transfer relationship is obtained by nonlinear parameter estimation iterative optimization model parameter, whether model is reduced to first order is judged using fitting degree index and is re-identified, and discrete domain transfer relationship is converted into continuous domain transfer relationship, the order reduction processing is executed to second-order model, inertia time constant and damping coefficient are calculated by applying step excitation, re-optimization process, until meeting the condition output node inertia.The present application cooperates nonlinear parameter solution and adaptive window searching, while considering identification speed, significantly improves the accuracy and robustness of node inertia identification.

Description

Technical Field

[0001] This invention relates to the field of power system inertia identification technology, specifically to a method, system, device, and medium for identifying the node inertia of a power system. Background Technology

[0002] As new energy power generation technologies mature, the penetration rate of new energy power sources such as wind power and photovoltaic power connected to the grid through power electronic converters continues to rise in the power system, gradually replacing traditional synchronous generator sets. This structural change has led to an overall decrease in the rotational inertia reserve of the power system, and the distribution of inertia among different nodes of the power grid shows a significant unevenness. The traditional method of using the overall inertia of the system as a single evaluation indicator is no longer able to accurately reflect the differences in the actual inertia support capacity of each node during frequency disturbance events.

[0003] Existing methods for identifying nodal inertia based on autoregressive moving average models with exogenous inputs essentially extract inertia parameters using the transfer function relationship derived from the swing equations of synchronous generators. This requires extracting a finite-length time window of data before and after a system disturbance event for model parameter estimation. Current methods face two significant technical bottlenecks: First, the selection of the post-disturbance data window length lacks effective theoretical guidance. A window that is too short provides insufficient data information to support accurate modeling, while a window that is too long may introduce post-disturbance dynamics that deviate from the ideal transfer function assumption. Existing solutions often indirectly avoid this problem by combining random windows or batch comparing residuals, failing to provide an explicit criterion between window length and the validity of the identification results based on the physical information inherent in the transfer function. Second, solving for model parameters is essentially a nonlinear optimization problem, and current research mostly uses a linear least squares framework for approximate solutions, neglecting the non-convex nature of the residual function in the parameter space. This results in the accuracy of inertia identification being limited by the approximation error of the solution method itself. Summary of the Invention

[0004] In view of the above-mentioned existing problems, the present invention provides a method, system, device and medium for identifying the node inertia of a power system, in order to solve the problems of lack of effective selection criteria for the length of the data window after disturbance and insufficient accuracy of nonlinear solution of model parameters in the prior art.

[0005] To address the aforementioned technical problems, a method for identifying the nodal inertia of a power system is proposed, including: The system detects load step disturbance events, preprocesses the active power and frequency data of the target node, and sets the window length before and after the disturbance. Based on the current window length, it extracts the active power and frequency sequences, establishes a dynamic parameterized model, obtains initial values ​​for the model parameters, and uses a nonlinear parameter identification method to iteratively optimize the model parameters with the frequency sequence as the fitting target, obtaining the discrete domain transfer relationship of the dynamic parameterized model. It calculates the goodness-of-fit index between the dynamic parameterized model output and the frequency sequence. When the goodness-of-fit index exceeds a preset threshold, the model order is changed from second-order to first-order, and the linear parameter identification method is executed again to obtain the first-order discrete domain transfer relationship. If the current model is second-order, the discrete domain transfer relation is maintained; otherwise, the second-order discrete domain transfer relation is maintained. The discrete domain transfer relation is converted into a continuous domain transfer relation. If the current model is second-order, the continuous domain transfer relation is reduced to a first-order continuous domain transfer relation. A step excitation is applied to the first-order continuous domain transfer relation, and the inertial time constant and damping coefficient are calculated. If the absolute value of the damping coefficient is lower than the preset threshold, the inertial time constant is used as the nodal inertia identification result; otherwise, the window length adaptive adjustment mechanism is started according to the damping coefficient, the window length after the disturbance is adjusted, and the dynamic parameterization model is established again until the absolute value of the damping coefficient is lower than the preset threshold, and then the nodal inertia identification result is output.

[0006] As a preferred embodiment of the node inertia identification method for a power system according to the present invention, the preprocessing includes downsampling the collected raw data and applying low-pass filtering to the downsampled data, wherein the passband frequency of the low-pass filtering is determined based on the lowest oscillation mode in the power system. The filtered data is detrended by eliminating the trend term through differential sampling of adjacent sampling points, and the rated capacity parameters of the corresponding power generation unit of the target node and the rated frequency parameters of the power grid are obtained. The active power data and frequency data after detrending are normalized to obtain the normalized active power sequence and frequency sequence.

[0007] As a preferred embodiment of the node inertia identification method for a power system according to the present invention, the method of obtaining the initial values ​​of the model parameters includes initializing the dynamic parameterized model as a second-order autoregressive exogenous input structure that ignores random disturbance terms. Using the active power and frequency sequences extracted from the preprocessed data, the initial values ​​of the model parameters are solved by linear regression parameter estimation.

[0008] As a preferred embodiment of the nodal inertia identification method for a power system according to the present invention, the nonlinear parameter identification method includes: constructing an optimization objective function characterized by the sum of squared residuals of prediction errors, and using the initial values ​​of the parameters of the dynamic parameterized model as the starting point of the iteration. In each iteration step, the partial derivative matrix of the residual vector with respect to the parameter vector is calculated, an adjustment factor is introduced, and the parameter update amount is solved to update the parameter vector. The adjustment factor is dynamically adjusted according to the decreasing trend of the objective function until the convergence criterion is met, at which point the iteration stops, the parameter estimation vector is output, and the discrete domain transitivity is obtained based on the current parameter estimation vector.

[0009] As a preferred embodiment of the nodal inertia identification method for a power system according to the present invention, the calculation of the fitting degree index between the output of the dynamic parameterized model and the frequency sequence includes calculating the fitting degree index based on the degree of deviation between the output sequence of the dynamic parameterized model and the frequency sequence. When the fit index is higher than the preset threshold, the model order is changed from second order to first order, and the obtained parameter estimation vector is used as the initial parameter to execute the nonlinear parameter identification method again to obtain the first-order discrete domain transitivity.

[0010] As a preferred embodiment of the node inertia identification method for a power system according to the present invention, the step of converting the discrete domain transitivity into a continuous domain transitivity includes using a bilinear transformation relationship to convert the current discrete domain transitivity into a continuous domain transitivity. When the current model is second-order, the state space realization of the continuous domain transitivity is constructed, the controllable Glam matrix and the observable Glam matrix are calculated, the system state importance metric is determined using the controllable Glam matrix and the observable Glam matrix, and minor states are eliminated according to the state importance metric to obtain the first-order continuous domain transitivity. Before order reduction, the pole location characteristics of the continuous domain transitivity are examined to determine system stability. If the system is determined to be unstable, the order reduction is terminated and the model is directly modified to first order.

[0011] As a preferred embodiment of the nodal inertia identification method for a power system according to the present invention, the calculation of the inertia time constant and damping coefficient includes applying a step excitation with an amplitude equal to the change in disturbance power to the first-order continuous domain transmission relationship to obtain a frequency step response. The initial frequency change rate is extracted from the frequency step response, the inertial time constant is calculated based on the initial frequency change rate and the power change, and the damping coefficient is extracted from the continuous domain transmission relationship. When the absolute value of the damping coefficient is lower than the preset threshold, the inertial time constant is output as the nodal inertia identification result. When the absolute value of the damping coefficient is greater than or equal to the preset threshold, the sign of the damping coefficient is determined. If the sign is positive, the length of the perturbation window is increased by a preset step size. If the sign is negative, the length of the perturbation window is decreased by a preset step size. The adjusted length of the perturbation window is then returned to the model building and parameter identification process.

[0012] The beneficial effects of this preferred technical solution are as follows: by combining the nonlinear parameter estimation with damping factor adaptive adjustment and the automatic window search guided by damping coefficient feedback, a clear physical criterion for window length selection is provided under the premise of reducing the computational overhead of the low-order model, so that the node inertia identification can achieve higher accuracy and robustness while ensuring speed.

[0013] As a preferred embodiment of the node inertia identification system for a power system according to the present invention, it is characterized by including a data acquisition and preprocessing module, a dynamic parameterized model construction and initialization module, a nonlinear parameter identification module, a model order adaptive adjustment and continuous processing module, and an inertia constant extraction and window adaptive adjustment module.

[0014] The data acquisition and preprocessing module is responsible for detecting load step disturbance events. When the output power step amplitude of the power generation unit connected to the target node exceeds 2% of the steady-state power before the disturbance, the disturbance time is recorded and the active power and frequency data of the target node are acquired. After the acquisition is completed, the data is downsampled to reduce the sampling frequency to 100Hz. The high-frequency electromechanical oscillation components and noise are filtered out by a Butterworth filter with the passband frequency set according to the lowest oscillation mode of the system. The trend term is removed by the first-order difference between adjacent sampling points. The rated capacity of the power generation unit and the rated frequency of the power grid are obtained and the per-unit conversion is completed to provide a standardized active power sequence and frequency sequence.

[0015] The dynamic parameterized model construction and initialization module is used to extract the active power sequence as input and the frequency sequence as the fitting target from the preprocessed data according to the set fixed window before the disturbance and the adjustable window after the disturbance, and establish an ARMAX model as the dynamic parameterized model. The initial order of the model is set to second order, and after ignoring the moving average term, it degenerates into an ARX structure. The initial values ​​of the autoregressive parameters and exogenous input parameters are solved by the least squares method using the extracted data sequence. The obtained initial parameter values ​​are used as the iterative starting point for nonlinear parameter identification.

[0016] The nonlinear parameter identification module is used to iteratively optimize the ARMAX model parameters using the Levenberg-Marquardt algorithm, starting from the initial values ​​of the model parameters. In each iteration, the prediction error residual vector and its Jacobian matrix are calculated. A damping factor is introduced to solve for the parameter update amount and update the parameter vector. The damping factor is dynamically adjusted based on the decrease in the sum of squared residuals until the convergence condition is met. The module outputs the parameter estimate vector and the corresponding discrete transfer function, thus accurately extracting the dynamic relationship between the model input and output.

[0017] The model order adaptive adjustment and continuousization module is used to calculate the bestfit of the model output sequence and the frequency sequence. When the bestfit is higher than 98%, the model is reduced to the first order and nonlinear parameter identification is re-executed to obtain the first-order discrete transfer function. The discrete transfer function is then converted into a continuous transfer function through the Tustin bilinear transform. A controllable canonical form is constructed for the second-order model and the controllable and observable Gram matrices are calculated. The equilibrium truncation order reduction is completed by eliminating minor states based on Hankel singular values. Before the order reduction, the real parts of the poles are checked to confirm the system stability.

[0018] The inertial constant extraction and window adaptive adjustment module is used to apply a step excitation with an amplitude equal to the change in perturbation power to a first-order continuous transfer function, obtain the time-domain frequency response through inverse Laplace transform, extract the frequency change rate at time zero, and calculate the inertial time constant in combination with the power change. At the same time, the damping coefficient is extracted from the transfer function. When the absolute value of the damping coefficient is lower than a preset threshold, the identification result is directly output. Otherwise, the direction of window increase or decrease is determined according to the sign of the damping coefficient, the length of the window after perturbation is adjusted according to a preset step size, and the feedback is sent to the model building module to re-execute the identification process until the convergence condition is met.

[0019] A computer device includes a memory and a processor, the memory storing a computer program, the processor executing the computer program to implement the steps of a method for identifying the node inertia of a power system.

[0020] A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of a method for identifying node inertia in a power system.

[0021] The beneficial effects of this invention are as follows: This invention achieves automated capture of disturbance moments by setting a power step amplitude threshold as the disturbance detection criterion; it sequentially performs downsampling, Butterworth low-pass filtering with a passband frequency set according to the system's lowest oscillation mode, first-order differential detrending, and per-unit processing on the collected data, thus suppressing high-frequency interference and signal drift while retaining the low-frequency information required for inertia identification, making the data from nodes with different capacities comparable; by fixing the window length before the disturbance and only adaptively adjusting the window length after the disturbance, window optimization focuses on the data segment that has the greatest impact on the identification results, reducing redundant degrees of freedom; in the model initialization stage, the moving average term is ignored and the least squares method is used to obtain reliable initial parameters, providing a good starting point for nonlinear iteration; and the Levenberg-Marquardt algorithm is used to non-linearly optimize the model parameters. Linear estimation, by leveraging the dynamic balance of damping factor gradient descent and the search path along the Gauss-Newton direction, brings the parameter estimation results closer to the true optimal value, overcoming the accuracy limitations of linear solution methods on non-convex residual functions. A goodness-of-fit index is used to determine whether to reduce the model to first order, reducing the number of parameters to be identified and compressing computation time while maintaining accuracy. Furthermore, a bilinear transformation preserves the amplitude and phase characteristics of the mapping from the discrete domain to the continuous domain, and a balance truncation reduction method based on Hankel singular values ​​is employed, with pole sign checks added to prevent unstable systems from entering inertia calculations. By constructing an explicit mapping relationship between the damping coefficient and the effectiveness of the window length, the direction of window increase or decrease is determined by the sign of the damping coefficient, and the optimal window is gradually approximated with a preset step size, transforming window selection from empirical trial-and-error to a closed-loop automatic search guided by clear criteria. Attached Figure Description

[0022] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying 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.

[0023] Figure 1 This is a flowchart illustrating a method for identifying the node inertia of a power system, as provided in one embodiment of the present invention.

[0024] Figure 2 The present invention provides a system scheme flowchart for a nodal inertia identification system for a power system according to an embodiment of the present invention. Detailed Implementation

[0025] To make the above-mentioned objects, features, and advantages of the present invention more apparent and understandable, specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the protection scope of the present invention.

[0026] Example 1, referring to Figure 1 As an embodiment of the present invention, a method for identifying the nodal inertia of a power system is provided, comprising: S100: Detects load step disturbance events, performs preprocessing on the active power and frequency data of the target node, and sets the window length before and after the disturbance.

[0027] S200: Based on the current window length, extract the active power sequence and frequency sequence, establish a dynamic parameterized model, obtain the initial values ​​of the model parameters, and use a nonlinear parameter identification method to iteratively optimize the model parameters with the frequency sequence as the fitting target, thereby obtaining the discrete domain transfer relationship of the dynamic parameterized model.

[0028] S300: Calculate the goodness of fit index between the output of the dynamic parameterized model and the frequency sequence. When the goodness of fit index is higher than the preset threshold, change the model order from second order to first order and execute the linear parameter identification method again to obtain the first-order discrete domain transfer relationship. Otherwise, maintain the second-order discrete domain transfer relationship.

[0029] S400: Convert the discrete domain transfer relation into a continuous domain transfer relation. When the current model is second-order, perform order reduction processing on the continuous domain transfer relation to obtain a first-order continuous domain transfer relation. Apply a step excitation to the first-order continuous domain transfer relation and calculate the inertial time constant and damping coefficient. When the absolute value of the damping coefficient is lower than the preset threshold, the inertial time constant is used as the nodal inertia identification result; otherwise, start the window length adaptive adjustment mechanism according to the damping coefficient, adjust the window length after the disturbance, and return to the establishment of the dynamic parameterized model until the absolute value of the damping coefficient is lower than the preset threshold, and then output the nodal inertia identification result.

[0030] It should be noted that this invention combines nonlinear parameter estimation based on dynamic adjustment of damping factor with adaptive window search based on the sign of damping coefficient. This reduces the model order to compress the computational load, while transforming the selection of data window length from empirical trial and error to a closed-loop automatic search process guided by clear physical criteria. This balances identification speed and significantly improves the accuracy and robustness of nodal inertia identification.

[0031] Example 2, refer to Figure 1This is a second embodiment of the present invention, which provides a method for identifying the nodal inertia of a power system, including: In step S100, the detection of the load step disturbance event includes steps S101 to S102: S101: During the operation of the power system, the output power of the power generation unit connected to the target node is continuously monitored. When a step change in output power is detected and the step amplitude exceeds 2% of the steady-state power before the disturbance occurs, the current moment is defined as the disturbance moment.

[0032] S102: Based on the current disturbance time, collect the active power data and frequency data of the target node. The collected data includes data from the period before and after the disturbance, and the data source is real-time measurement data from the phasor measurement unit (PMU).

[0033] The output active power data and frequency data serve as the raw input for the step preprocessing.

[0034] Furthermore, in step S100, the preprocessing includes steps S111 to S114: S111: Downsample the acquired raw data to reduce the data sampling frequency, that is, downsample the data acquired by the phasor measurement unit (PMU) to reduce the data sampling frequency to 100Hz.

[0035] S112: Apply low-pass filtering to the downsampled data. The passband frequency of the low-pass filtering is determined based on the lowest oscillation mode in the power system.

[0036] A Butterworth filter is used for low-pass filtering to remove high-frequency dynamic components, including high-frequency electromechanical oscillations and noise. The order of the Butterworth filter is calculated using the following formula: in, Let the filter order be . For stopband gain, For stopband frequency, For passband gain, This is the passband frequency.

[0037] It should be noted that the specific transfer function of the Butterworth filter is expressed as follows: in, For the filter transfer function, For the Laplace operator, This is the cutoff angular frequency.

[0038] S113: Perform detrending processing on the filtered data, eliminate the trend term by differential sampling between adjacent sampling points, and obtain the rated capacity parameters of the power generation unit corresponding to the target node and the rated frequency parameters of the power grid.

[0039] The trend term is removed by performing a first-order difference on the data from two adjacent sampling points, specifically in the following form: in, Let be the first-order difference value of the active power at the t-th sampling time, representing the change in power between the current time and the previous time. This is the raw active power data at the t-th sampling time. This is the raw active power data at the (t-1)th sampling time. Let be the first-order difference value of the frequency at the t-th sampling time, representing the change in frequency between the current time and the previous time. For the target node frequency data at the t-th sampling time, The target node frequency data is at the (t-1)th sampling time.

[0040] S114: Perform a per-unit transformation on the detrended active power data and frequency data to obtain a per-unit active power sequence and frequency sequence.

[0041] Obtain the rated capacity of the power generation unit corresponding to the measured node and the rated frequency of the power grid, and perform per-unit calculations on the power and frequency. The formula is as follows: in, The frequency data is standardized and dimensionless. The rated frequency of the power grid. To sample the active power data of the target node, The rated capacity of the power generation unit connected to the target node, This is the normalized active power data.

[0042] Furthermore, in step S100, setting the window length before and after the disturbance includes setting the data window length before the disturbance time and the data window length after the disturbance time. The window length before the disturbance remains unchanged throughout the entire identification process, and the selected range is 0.1s to 0.5s. The window length after the disturbance will be continuously adjusted during the identification process, and the initial value is selected to be about 3s. The initial value of the window length after the disturbance will be used as the basis for truncating the data sequence and will be iteratively corrected in the window length adaptive adjustment mechanism.

[0043] In step S200, obtaining the initial values ​​of the model parameters includes steps S201-S202: S201: Based on the set current post-disturbance window length and fixed pre-disturbance window length, extract the corresponding active power sequence from the preprocessed data as the model input, and extract the corresponding frequency sequence as the fitting target for the model output data.

[0044] The dynamically parameterized model is initialized as a second-order autoregressive exogenous input structure that ignores random disturbances. An autoregressive moving average model with exogenous input (ARMAX model) is established as the dynamically parameterized model. The initial order is set to second, and the initial values ​​of the model's parameters are determined using a linear model approximation initialization method. The ARMAX model is represented as follows: in, Let be the output value at the t-th sampling time. For the first The input value at each sampling time point, It is a white noise sequence. Let the order be the autoregressive order. For exogenous input order, The order of the moving average. Given the input delay order, The backward shift factor. , and The parameters to be estimated for the ARMAX model are... It is an autoregressive polynomial. For exogenous input polynomials, It is a moving average polynomial.

[0045] It should be noted that the initial order of the ARMAX model is set to 2, that is... .

[0046] S202: Using the active power sequence and frequency sequence extracted from the preprocessed data, the initial values ​​of the model parameters are solved by the linear regression parameter estimation method.

[0047] The ARMAX model is initialized using a linear model approximation initialization method, ignoring the white noise model of ARMAX, that is, letting The formula is expressed as: in, It is an autoregressive polynomial. The backward shift factor. For exogenous input polynomials, Let be the white noise sequence at the t-th sampling time. After ignoring the moving average term, it degenerates into a regular residual term.

[0048] The regression form is expressed as: in, Let i be the i-th autoregressive parameter. This represents the output value at the ti-th sampling time, i.e., the historical value of the output sequence. For the j-th exogenous input parameter, For the first The input value at each sampling time, i.e., the historical value of the input sequence.

[0049] The matrix form is as follows: in, The vector formed by the output sequence is composed of y(t) at each sampling time arranged in chronological order. This is a linear regression matrix, where each row consists of a regression vector corresponding to the given time step. Each regression vector contains historical output and input values. The parameter vector to be initialized. The residual vector is composed of the residual terms w(t) at each sampling time. These are the initial estimates of the parameter vector. It is the transpose of the linear regression matrix. For matrix The inverse matrix.

[0050] Furthermore, in step S200, the iterative optimization model parameters include steps S211~S214: S211: Construct an optimization objective function characterized by the sum of squared residuals of the prediction error, and use the initial values ​​of the parameters of the dynamically parameterized model as the starting point for iteration.

[0051] The parameters to be identified are represented as follows: in, For the autoregressive parameter vector, For the exogenous input parameter vector, The moving average parameter vector is used to merge the parameters to be identified into a single vector. .

[0052] S212: In each iteration step, calculate the partial derivative matrix of the residual vector with respect to the parameter vector, introduce the adjustment factor, and solve for the parameter update amount to update the parameter vector.

[0053] Based on the data window length L after the disturbance, data after the disturbance is extracted, and nonlinear estimation of each parameter is performed. The process is expressed by the following formula: The expanded form is represented as: in, Let be the prediction error at the t-th sampling time, representing the deviation between the actual output value and the model prediction value. Given a parameter vector θ, this is the model's predicted output value at the t-th sampling time. This is the parameter vector of the ARMAX model, containing autoregressive parameters, exogenous input parameters, and moving average parameters. For sampling time index, Let be the parameter of the k-th moving average. The input delay order.

[0054] S213: Dynamically adjust the value of the adjustment factor according to the downward trend of the objective function until the convergence criterion is met, then stop the iteration and output the parameter estimate vector.

[0055] It should be noted that an objective function is established, and at the r-th iteration point, a first-order Taylor expansion is performed on the residual vector, and the parameters are updated according to the given iteration step size. .

[0056] The objective function is expressed as: in, Let be the optimal estimate of the parameter vector, representing the parameter values ​​that minimize the objective function. To find the values ​​of the parameter vector θ that minimize the objective function, To optimize the objective function, The sum of the squares of the L2 norm of the residual vector w(θ) is the sum of the squares of each element. The total number of sampling points in the data sequence. This is the index of the starting sampling time for parameter estimation.

[0057] The first-order Taylor expansion formula is expressed as: in, To update the parameter vector to The residual vector after, This represents the current value of the parameter vector at the r-th iteration. Let be the parameter update amount, i.e., the adjustment step size of the parameter vector in the r-th iteration. For the r-th iteration, in the parameter vector The residual vector calculated at that point, Let be the partial derivative matrix of the residual vector with respect to the parameter vector, i.e., the Jacobian matrix. Let w(θ) be the partial derivative of the residual vector w(θ) with respect to the parameter vector θ, which is the definition of the Jacobian matrix. The starting sampling time The row vector formed by the partial derivatives of the prediction error with respect to each component of the parameter vector θ. For the first The row vector formed by the partial derivatives of the prediction error with respect to each component of the parameter vector at each sampling time. Let be a row vector composed of the partial derivatives of the prediction error at the Nth sampling time with respect to each component of the parameter vector. It is the transpose of the parameter vector θ.

[0058] Update parameters based on the given iteration step size The formula is expressed as: in, For the parameter point at the r-th iteration The Jacobian matrix calculated at this point is the partial derivative matrix of the residual vector with respect to the parameter vector. Jacobian matrix The transpose of the matrix, This is the damping factor at the r-th iteration, used to adjust the weight distribution of the search direction of the iteration step between the gradient descent direction and the Gaussian-Newton direction. The default initial value is set to 0.001. It is an identity matrix with the same dimensions as the vector of parameters to be identified. Let be the parameter update amount in the r-th iteration, i.e., the adjustment step size of the parameter vector in this iteration. For the parameter point at the r-th iteration The residual vector calculated at that point, This represents the current value of the parameter vector at the r-th iteration. This represents the parameter vector value after the (r+1)th iteration.

[0059] It should be further explained that, This is the damping factor, with a default initial value of 0.001. It needs to be adjusted during iteration. The decrease in damping factor is continuously adjusted using the damping Gauss-Newton method. If the factor decreases, it is reduced to half its original value; if it does not decrease, it is increased to twice its original value, until the convergence condition is met. The convergence condition is expressed as follows: in, The parameter update amount in the r-th iteration The norm of represents the magnitude of change in the parameter vector. This is a preset convergence threshold for the parameter update norm, with a default value of 0.01. Let the objective function value at the parameter point be the value after the (r+1)th iteration update. Let $\mathbf{r}$ be the objective function value at the parameter point in the $r$-th iteration. The preset convergence threshold is set to the change in the objective function. The default value is 0.01. When the decrease in the objective function value is less than or equal to the current threshold, the iteration is considered to have converged.

[0060] S214: Finally, the parameter estimation vector is obtained. Based on the current parameter estimation vector, the discrete domain transitivity is obtained, expressed by the formula: in, The transform domain representation of the output sequence is the transform domain expression of the frequency sequence. The transform domain representation of the input sequence is the transform domain expression of the active power sequence. It is the discrete-domain transfer function, that is, the discrete transfer relationship between the system input and output.

[0061] In step S300, the calculation of the fitting index between the output of the dynamic parameterized model and the frequency sequence includes steps S301-S302: S301: Calculate the goodness-of-fit index based on the degree of deviation between the output sequence of the dynamic parameterization model and the frequency sequence; The model output sequence is calculated using the obtained parameter estimate vector, and the bestfit index is used to evaluate the goodness of fit between the model output sequence and the truncated frequency sequence. The bestfit index is used to calculate the goodness of fit between the ARMAX model output data and the fitting target, expressed by the formula: in, This is the actual output value vector. The average of the actual output value vector. Output value vectors for ARMAX models.

[0062] S302: The preset threshold is set to 98%. When the fit index is higher than the preset threshold (i.e., Bestfit value > 98%), it indicates that the second-order model has accurately captured the dynamic characteristics of the frequency sequence. The model order is changed from second-order to first-order, and the obtained parameter estimation vector is used as the initial parameter. The nonlinear parameter identification method is executed again to obtain the first-order discrete domain transfer relationship. When the fit index is lower than or equal to the preset threshold (i.e., Bestfit value ≤ 98%), it indicates that there is still room for optimization in the characterization of the frequency sequence by the second-order model. The order of the ARMAX model is maintained at second-order, and the obtained discrete transfer function is directly transferred.

[0063] Furthermore, in step S400, the conversion of the discrete domain transitive relation into a continuous domain transitive relation includes steps S401 to S406: S401: Use bilinear transformation to convert the current discrete domain transitive relation into a continuous domain transitive relation.

[0064] The discrete transfer function is converted into a continuous transfer function using the Tustin bilinear transform method. During the transformation, the backward shift factor is replaced with a rational fraction relating the Laplace operator and the discrete sampling step size. The sampling step size is taken as the sampling period after downsampling. The formula is expressed as: in, The backward shift factor. This is the discrete sampling step size, i.e., the sampling period after downsampling.

[0065] S402: The continuous transfer function obtained after the bilinear transformation is standardized to normalize the highest-order coefficients of the Laplace operator in the numerator and denominator polynomials, and a controllable canonical form state-space realization of the system is constructed. The transfer function is then transformed into a state equation form described by the state matrix, input matrix, and output matrix, expressed as: in, The time derivative of the state variable vector represents the rate of change of the state variables. For the state variable vector, Y is the input variable, and Y is the output variable. and The coefficients of the denominator polynomial of the continuous transfer function correspond to the coefficients of the constant term and the linear term, respectively. and are the coefficients of the numerator polynomial of the continuous transfer function, corresponding to the coefficients of the constant term and the linear term, respectively.

[0066] S403: Extract the state matrix, input matrix, and output matrix from the controllable canonical form. For continuous and stable systems, construct the controllable Gram matrix and the observable Gram matrix, respectively.

[0067] Among them, the controllable Gram matrix is ​​obtained by solving the Lyapunov equation formed by the state matrix and its transpose, and the input matrix and its transpose, and measures the ease with which each state component is excited by the input signal; the observable Gram matrix is ​​obtained by solving the Lyapunov equation formed by the transpose and original matrix of the state matrix, and the transpose and original matrix of the output matrix, and measures the degree to which each state component is reflected in the output signal.

[0068] Define matrices A, B, and C. For a continuous and stable system, the controllable Gram matrix and the observable Gram matrix can be solved. The process formula is expressed as follows: in, The controllable Gram matrix measures the ease with which each state component is excited by the input signal. To observe the Gram matrix, we need to measure the degree to which each state component is represented in the output signal. The system state matrix, For the input matrix, For the output matrix, This is a transpose operation.

[0069] S404: When the current model is second-order, construct the state space implementation of the continuous domain transitivity, calculate the controllable Girahm matrix and the observable Girahm matrix, and use the controllable Girahm matrix and the observable Girahm matrix to determine the system state importance metric.

[0070] Furthermore, it should be noted that a state transformation matrix is ​​defined, and the controllable Gram matrix and the observable Gram matrix are simultaneously transformed using the current transformation matrix, such that the two transformed Gram matrices are both equal to the same diagonal matrix, which satisfies: in, This is the state transformation matrix, used to perform coordinate transformations in the state space, converting the system to the equilibrium coordinate system. It is the inverse of the state transformation matrix. This is the inverse of the state transformation matrix. It is the transpose of the state transformation matrix. It is a diagonal matrix.

[0071] Furthermore, the diagonal elements of the diagonal matrix are the Hankel singular values. A second-order system contains two singular values, each corresponding to one of the two state components. The Hankel singular value represents the system state importance of each state component; the larger the value, the more dominant the state component is in the system's input-output behavior.

[0072] In the second-order transfer function, the diagonal matrix is ​​represented as: in, This is the first Hankel singular value, corresponding to the system state importance metric for the first state component. This is the second Hankel singular value, corresponding to the system state importance metric for the second state component.

[0073] Before reducing the order, the pole location characteristics of the continuous domain transitivity are examined to determine the system stability. If the system is determined to be unstable, the order reduction is terminated and the model is directly modified to first order.

[0074] Before performing the order reduction operation, the pole locations of the transfer function are checked: if the real parts of all poles are negative, it indicates that the system is stable and the order reduction process can continue; if there are poles with positive real parts, it indicates that the system is unstable and the order reduction operation cannot be performed. The dynamic parameterized model is directly modified to a first-order model, and the following steps are skipped to directly enter the inertia constant extraction stage.

[0075] S405: Eliminate minor states based on state importance metrics to obtain the first-order continuous domain transitivity.

[0076] Perform coordinate transformations on each matrix in the state-space implementation. Define matrices A, B, and C and their inverses using S403. Update the state matrix, input matrix, and output matrix according to the similarity transformation rules to obtain the system equations in the equilibrium coordinate system. The process formula is expressed as follows: in, This is the state matrix after coordinate transformation. The input matrix after coordinate transformation. This is the output matrix after coordinate transformation. The state variable vector in the equilibrium coordinate system. The step change in active power is the input variable to the system. The system output variable is the frequency response.

[0077] S406: Divide the state matrix in the equilibrium coordinate system into blocks according to the first state and the second state, divide the input matrix into blocks according to the corresponding two state components, and divide the output matrix into blocks according to the corresponding two state components.

[0078] The block-based formula for the state matrix is ​​as follows: in, This is the system state matrix in the equilibrium coordinate system after coordinate transformation. This is the top-left submatrix of the state matrix, corresponding to the dynamic characteristics of the first state component itself. The upper right submatrix of the state matrix represents the dynamic coupling effect of the second state component on the first state component. The lower left submatrix of the state matrix represents the dynamic coupling effect of the first state component on the second state component. This is the lower right submatrix of the state matrix, corresponding to the dynamic characteristics of the second state component itself. This is the input matrix in the equilibrium coordinate system after coordinate transformation. The block submatrix above the input matrix represents the weights of the first state components affected by the input signal. This is the block submatrix below the input matrix, corresponding to the weights of the second-state components affected by the input signal. This is the output matrix in the equilibrium coordinate system after coordinate transformation. The left-hand submatrix of the output matrix represents the contribution weights of the first state components to the output signal. This is the block submatrix on the right side of the output matrix, corresponding to the contribution weight of the second state component to the output signal.

[0079] Compare the two Hankel singular values ​​obtained from S404. Assuming the second singular value is much smaller than the first, the second state is determined to be a minor state relative to the first state and is discarded. The first state component corresponding to the larger singular value and its associated submatrix are retained. The reduced-order dynamical system equations are constructed, and the process formula is expressed as follows: in, Let be the time derivative of the reduced-order state variable vector, representing the rate of change of the system state after the reduction. This is the reduced state variable vector. This is the reduced-order state matrix, which is the submatrix corresponding to the first state component in the original block, describing the dynamic characteristics of the reduced-order system. This is the reduced-order input matrix, i.e., the submatrix of the original block's input matrix corresponding to the first state component, describing the effect of the reduced-order input signal on the state. For the output variables of the reduced-order system, This is the output matrix after order reduction, which is the submatrix of the output matrix in the original block corresponding to the first state component. It describes the mapping relationship from the state variables to the output after order reduction. This is a reduced-order first-order continuous transfer function that describes the continuous-domain transfer relationship from input to output.

[0080] The corresponding first-order continuous transfer function is derived from the reduced state matrix, input matrix, and output matrix and transformed into standard form.

[0081] Furthermore, in step S400, the calculation of the inertial time constant and damping coefficient includes steps S411~S415: S411: Extract the continuous transfer function and decompose the inertial process.

[0082] The first-order continuous transfer function obtained after order reduction can be decomposed into the sum of the through-gain term and the first-order inertial element. Since the order reduction process has described the system response as an inertial-dominated process, the coefficient of the first-order term in the numerator approaches zero, and the through-gain term can be ignored. The transfer function can be further simplified to the standard form containing only the first-order inertial element, where the numerator is the constant gain and the denominator is the sum of the Laplace operator and the poles.

[0083] The continuous transfer function is expressed as: in, The first-order continuous transfer function obtained after order reduction describes the continuous-domain transfer relationship between the system input and output. The coefficients of the first-order term in the numerator of the transfer function are... The constant term coefficients of the transfer function numerator are... This is the constant term in the denominator of the transfer function, i.e., the pole values ​​of the transfer function.

[0084] The transfer function can be decomposed into the formula for the through gain and the first-order inertial process as follows: in, To decompose the molecular gain of the first-order inertial link.

[0085] The simplified transfer function is expressed as the ratio of frequency change to power change, and its coefficients are compared with those of the standard first-order inertial transfer function, which includes the inertial time constant and damping coefficient. This establishes the correspondence between gain and inertial time constant, and between poles and damping coefficients and inertial time constant.

[0086] The final transfer function that needs to be identified is expressed as: in, To simplify the transfer function to contain only a first-order inertial element, This represents the frequency change in the complex frequency domain. This represents the power change in the complex frequency domain. The nodal inertial time constant characterizes the node's ability to resist frequency changes. is the nodal damping coefficient, which characterizes the damping effect of the system itself when the frequency changes.

[0087] S412: Apply a step excitation with an amplitude equal to the change in perturbation power to the first-order continuous domain transfer relation to obtain a frequency step response.

[0088] Applying a step input with an amplitude equal to the change in disturbance power to the simplified first-order inertial transfer function, multiplying the Laplace transform of the step input with the transfer function in the complex frequency domain yields the complex frequency domain expression of the step response. Performing an inverse Laplace transform on the complex frequency domain expression yields the time domain expression of the frequency response, which is a function that increases exponentially with time and approaches its steady-state value. The process formula is as follows: in, The magnitude of the change in disturbance power. is the base of the natural logarithm. This is an exponentially decaying term that describes the dynamic process by which the frequency response gradually approaches its steady-state value over time.

[0089] Taking the time derivative of the time-domain expression reveals the rate of frequency change over time. The rate of change decays exponentially. Letting time approach zero, we take the initial rate of frequency change, which is equal to the product of the transfer function gain and the step amplitude. The formula is as follows: in, It is a function of the rate of change of frequency with time, and is the derivative of the frequency response with respect to time. At the initial time (t=0) + The rate of change of frequency is the maximum rate of change of the frequency response.

[0090] S413: Extract the initial frequency change rate from the frequency step response, and calculate the inertial time constant based on the initial frequency change rate and power change.

[0091] According to the physical definition of inertia, the inertial time constant is determined by the ratio of the change in disturbance power to the rate of change of frequency at the initial moment. Substituting the obtained rate of change of frequency at the initial moment into the definition, the calculated value of the inertial time constant is obtained, expressed by the formula: in, The calculated inertial time constant, This represents the change in disturbance power.

[0092] S414: Extract the damping coefficient from the continuous domain transfer relationship. When the absolute value of the damping coefficient is lower than the preset threshold, output the inertial time constant as the nodal inertia identification result.

[0093] Based on the established correspondence between poles, damping coefficients, and inertial time constants, the corresponding damping coefficients are calculated using the calculated inertial time constants and transfer function pole values. The formula is as follows: in, is the damping coefficient.

[0094] The absolute value of the damping coefficient is compared with the preset reference threshold. If the absolute value of the damping coefficient is less than the reference threshold, it indicates that the current perturbation window length can enable the identification model to accurately reflect the pure inertial response characteristics of the node, and the window length is reasonable. The inertial time constant calculated in the step is used as the final node inertia identification result output, and the identification process ends.

[0095] S415: When the absolute value of the damping coefficient is greater than or equal to a preset threshold, determine the sign of the damping coefficient. If the sign is positive, increase the length of the perturbation window by a preset step size; if the sign is negative, decrease the length of the perturbation window by a preset step size. Return the adjusted length of the perturbation window to the model building and parameter identification process. The formula is expressed as: in, To identify the window length, To obtain the new window length, The iteration step size for window updates.

[0096] Specifically, when the absolute value of the damping coefficient is greater than or equal to the preset threshold, the direction of window length adjustment is determined based on the sign of the damping coefficient; if the damping coefficient is negative, it indicates that the current window length after the disturbance is too short and insufficient to fully cover the dominant period of the inertial response, so the window length is increased according to the preset iteration step size; if the damping coefficient is positive, it indicates that the current window length after the disturbance is too long and has introduced dynamic components that deviate from the ideal inertial response in the later stage of the disturbance, so the window length is decreased according to the preset iteration step size.

[0097] The updated window length replaces the original window length. The data sequence is re-trunculated, a dynamic parameterized model is established, nonlinear parameter identification is performed, the fit is evaluated and the model order is reduced, the continuous transfer function is converted and the inertial constant and damping coefficient are extracted. The above process is repeated until the absolute value of the damping coefficient meets the threshold condition, and the final nodal inertia identification result is output.

[0098] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

[0099] Example 3, referring to Figure 2This is the third embodiment of the present invention, which provides a nodal inertia identification system for a power system, including a data acquisition and preprocessing module, a dynamic parameterized model construction and initialization module, a nonlinear parameter identification module, a model order adaptive adjustment and continuous processing module, and an inertia constant extraction and window adaptive adjustment module.

[0100] The data acquisition and preprocessing module is responsible for detecting load step disturbance events. When the output power step amplitude of the power generation unit connected to the target node exceeds 2% of the steady-state power before the disturbance, the disturbance time is recorded and the active power and frequency data of the target node are acquired. After the acquisition is completed, the data is downsampled to reduce the sampling frequency to 100Hz. The high-frequency electromechanical oscillation components and noise are filtered out by a Butterworth filter with the passband frequency set according to the lowest oscillation mode of the system. The trend term is removed by the first-order difference between adjacent sampling points. The rated capacity of the power generation unit and the rated frequency of the power grid are obtained and the per-unit conversion is completed to provide a standardized active power sequence and frequency sequence.

[0101] The dynamic parameterized model construction and initialization module is used to extract the active power sequence as input and the frequency sequence as the fitting target from the preprocessed data according to the set fixed window before the disturbance and the adjustable window after the disturbance, and establish an ARMAX model as the dynamic parameterized model. The initial order of the model is set to second order, and after ignoring the moving average term, it degenerates into an ARX structure. The initial values ​​of the autoregressive parameters and exogenous input parameters are solved by the least squares method using the extracted data sequence. The obtained initial parameter values ​​are used as the iterative starting point for nonlinear parameter identification.

[0102] The nonlinear parameter identification module is used to iteratively optimize the ARMAX model parameters using the Levenberg-Marquardt algorithm, starting from the initial values ​​of the model parameters. In each iteration, the prediction error residual vector and its Jacobian matrix are calculated. A damping factor is introduced to solve for the parameter update amount and update the parameter vector. The damping factor is dynamically adjusted based on the decrease in the sum of squared residuals until the convergence condition is met. The module outputs the parameter estimate vector and the corresponding discrete transfer function, thus accurately extracting the dynamic relationship between the model input and output.

[0103] The model order adaptive adjustment and continuousization module is used to calculate the bestfit of the model output sequence and the frequency sequence. When the bestfit is higher than 98%, the model is reduced to the first order and nonlinear parameter identification is re-executed to obtain the first-order discrete transfer function. The discrete transfer function is then converted into a continuous transfer function through the Tustin bilinear transform. A controllable canonical form is constructed for the second-order model and the controllable and observable Gram matrices are calculated. The equilibrium truncation order reduction is completed by eliminating minor states based on Hankel singular values. Before the order reduction, the real parts of the poles are checked to confirm the system stability.

[0104] The inertial constant extraction and window adaptive adjustment module is used to apply a step excitation with an amplitude equal to the change in perturbation power to a first-order continuous transfer function, obtain the time-domain frequency response through inverse Laplace transform, extract the frequency change rate at time zero, and calculate the inertial time constant in combination with the power change. At the same time, the damping coefficient is extracted from the transfer function. When the absolute value of the damping coefficient is lower than a preset threshold, the identification result is directly output. Otherwise, the direction of window increase or decrease is determined according to the sign of the damping coefficient, the length of the window after perturbation is adjusted according to a preset step size, and the feedback is sent to the model building module to re-execute the identification process until the convergence condition is met.

[0105] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

[0106] Example 4, the fourth embodiment of the present invention, differs from the previous three embodiments in that: If the aforementioned functions are implemented as software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this invention, essentially, or the part that contributes to the prior art, or a portion 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 described in the various embodiments of this 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.

[0107] The logic and / or steps represented in the flowchart or otherwise described herein, for example, can be considered as a sequenced list of executable instructions for implementing logical functions, and can be embodied in any computer-readable medium for use by, or in conjunction with, an instruction execution system, apparatus, or device (such as a computer-based system, a processor-including system, or other system that can fetch and execute instructions from, an instruction execution system, apparatus, or device). For the purposes of this specification, "computer-readable medium" can be any means that can contain, store, communicate, propagate, or transmit programs for use by, or in conjunction with, an instruction execution system, apparatus, or device.

[0108] More specific examples of computer-readable media (a non-exhaustive list) include: electrical connections (electronic devices) having one or more wires, portable computer disk drives (magnetic devices), random access memory (RAM), read-only memory (ROM), erasable and editable read-only memory (EPROM or flash memory), fiber optic devices, and portable optical disc read-only memory (CDROM). Furthermore, computer-readable media can even be paper or other suitable media on which the program can be printed, because the program can be obtained electronically, for example, by optically scanning the paper or other medium, followed by editing, interpreting, or otherwise processing as necessary, and then stored in computer memory.

[0109] It should be understood that various parts of the present invention can be implemented in hardware, software, firmware, or a combination thereof. In the above embodiments, multiple steps or methods can be implemented in software or firmware stored in memory and executed by a suitable instruction execution system. For example, if implemented in hardware, as in another embodiment, it can be implemented using any one or a combination of the following techniques known in the art: discrete logic circuits having logic gates for implementing logical functions on data signals, application-specific integrated circuits (ASICs) having suitable combinational logic gates, programmable gate arrays (PGAs), field-programmable gate arrays (FPGAs), etc.

Claims

1. A method for identifying nodal inertia in a power system, characterized in that: include, Detect load step disturbance events, perform preprocessing on the active power and frequency data of the target node, and set the window length before and after the disturbance. Based on the current window length, active power and frequency sequences are extracted to establish a dynamic parameterized model, obtain initial values ​​of model parameters, and use a nonlinear parameter identification method to iteratively optimize model parameters with the frequency sequence as the fitting target, thereby obtaining the discrete domain transfer relationship of the dynamic parameterized model. Calculate the goodness of fit index between the output of the dynamic parameterized model and the frequency sequence. When the goodness of fit index is higher than the preset threshold, change the model order from second order to first order and execute the linear parameter identification method again to obtain the first-order discrete domain transitivity. Otherwise, maintain the second-order discrete domain transitivity. The discrete-domain transfer relation is converted into a continuous-domain transfer relation. When the current model is second-order, the continuous-domain transfer relation is reduced to a first-order continuous-domain transfer relation. A step excitation is applied to the first-order continuous-domain transfer relation, and the inertial time constant and damping coefficient are calculated. When the absolute value of the damping coefficient is lower than a preset threshold, the inertial time constant is used as the nodal inertia identification result; otherwise, the window length adaptive adjustment mechanism is activated according to the damping coefficient, the window length after the disturbance is adjusted, and the dynamic parameterization model is established again until the absolute value of the damping coefficient is lower than the preset threshold, and then the nodal inertia identification result is output.

2. The method for identifying nodal inertia in a power system as described in claim 1, characterized in that: The preprocessing includes downsampling the collected raw data and applying a low-pass filter to the downsampled data. The passband frequency of the low-pass filter is determined based on the lowest oscillation mode in the power system. The filtered data is detrended by eliminating the trend term through differential sampling of adjacent sampling points, and the rated capacity parameters of the corresponding power generation unit of the target node and the rated frequency parameters of the power grid are obtained. The active power data and frequency data after detrending are normalized to obtain the normalized active power sequence and frequency sequence.

3. The method for identifying the nodal inertia of a power system as described in claim 2, characterized in that: The process of obtaining the initial values ​​of the model parameters includes initializing the dynamically parameterized model as a second-order autoregressive exogenous input structure that ignores random perturbation terms; Using the active power and frequency sequences extracted from the preprocessed data, the initial values ​​of the model parameters are solved by linear regression parameter estimation.

4. The method for identifying nodal inertia in a power system as described in claim 3, characterized in that: The nonlinear parameter identification method includes constructing an optimization objective function characterized by the sum of squared residuals of the prediction error, and using the initial values ​​of the parameters of the dynamically parameterized model as the starting point for iteration. In each iteration step, the partial derivative matrix of the residual vector with respect to the parameter vector is calculated, an adjustment factor is introduced, and the parameter update amount is solved to update the parameter vector. The adjustment factor is dynamically adjusted according to the decreasing trend of the objective function until the convergence criterion is met, at which point the iteration stops, the parameter estimation vector is output, and the discrete domain transitivity is obtained based on the current parameter estimation vector.

5. The method for identifying nodal inertia in a power system as described in claim 4, characterized in that: The calculation of the fitting index between the output of the dynamic parameterized model and the frequency sequence includes calculating the fitting index based on the degree of deviation between the output sequence of the dynamic parameterized model and the frequency sequence. When the fit index is higher than the preset threshold, the model order is changed from second order to first order, and the obtained parameter estimation vector is used as the initial parameter to execute the nonlinear parameter identification method again to obtain the first-order discrete domain transitivity.

6. The method for identifying nodal inertia in a power system as described in claim 5, characterized in that: The step of converting a discrete-domain transitive relation into a continuous-domain transitive relation includes using a bilinear transformation relation to convert the current discrete-domain transitive relation into a continuous-domain transitive relation; When the current model is second-order, the state space realization of the continuous domain transitivity is constructed, the controllable Glam matrix and the observable Glam matrix are calculated, the system state importance metric is determined using the controllable Glam matrix and the observable Glam matrix, and minor states are eliminated according to the state importance metric to obtain the first-order continuous domain transitivity. Before order reduction, the pole location characteristics of the continuous domain transitivity are examined to determine system stability. If the system is determined to be unstable, the order reduction is terminated and the model is directly modified to first order.

7. The method for identifying nodal inertia in a power system as described in claim 6, characterized in that: The calculation of the inertial time constant and damping coefficient includes applying a step excitation with an amplitude equal to the change in disturbance power to the first-order continuous domain transfer relation to obtain a frequency step response; The initial frequency change rate is extracted from the frequency step response, the inertial time constant is calculated based on the initial frequency change rate and the power change, and the damping coefficient is extracted from the continuous domain transmission relationship. When the absolute value of the damping coefficient is lower than the preset threshold, the inertial time constant is output as the nodal inertia identification result. When the absolute value of the damping coefficient is greater than or equal to the preset threshold, the sign of the damping coefficient is determined. If the sign is positive, the length of the perturbation window is increased by a preset step size. If the sign is negative, the length of the perturbation window is decreased by a preset step size. The adjusted length of the perturbation window is then returned to the model building and parameter identification process.

8. A nodal inertia identification system for a power system, employing the nodal inertia identification method for a power system as described in any one of claims 1 to 7, characterized in that, It includes a data acquisition and preprocessing module, a dynamic parameterized model construction and initialization module, a nonlinear parameter identification module, a model order adaptive adjustment and continuous processing module, and an inertial constant extraction and window adaptive adjustment module. The data acquisition and preprocessing module is used to detect load step disturbance events, record the disturbance time and acquire active power and frequency data of the target node, perform downsampling processing on the data to reduce the sampling frequency to 100Hz, filter out noise through a Butterworth filter set according to the lowest oscillation mode of the system at the passband frequency, remove the trend term by using the first-order difference of adjacent sampling points, and obtain the rated capacity of the power generation unit and the rated frequency of the power grid to complete the per-unit conversion. The dynamic parameterized model construction and initialization module is used to extract the active power sequence as input and the frequency sequence as fitting target from the preprocessed data according to the set fixed window before the disturbance and the adjustable window after the disturbance, and establish an ARMAX model as a dynamic parameterized model. The initial values ​​of the autoregressive parameters and exogenous input parameters are solved by the least squares method using the extracted data sequence, and the obtained initial values ​​of the parameters are used as the iterative starting point for nonlinear parameter identification. The nonlinear parameter identification module is used to iteratively optimize the ARMAX model parameters using the Levenberg-Marquardt algorithm, starting from the initial values ​​of the model parameters, and dynamically adjust the damping factor according to the decrease of the sum of squared residuals until the convergence condition is met, and output the parameter estimate vector and the corresponding discrete transfer function. The model order adaptive adjustment and continuousization module is used to calculate the bestfit of the model output sequence and frequency sequence, and to convert the discrete transfer function into a continuous transfer function through the Tustin bilinear transformation. It constructs a controllable canonical form for the second-order model and calculates the controllable and observable Gram matrices. It completes the equilibrium truncation and order reduction by eliminating minor states based on Hankel singular values. The inertial constant extraction and window adaptive adjustment module is used to apply a step excitation with an amplitude equal to the change in disturbance power to a first-order continuous transfer function, obtain the time-domain frequency response through inverse Laplace transform, extract the frequency change rate at time zero and calculate the inertial time constant in combination with the power change, and extract the damping coefficient from the transfer function.

9. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the processor executes the computer program, it implements the steps of the nodal inertia identification method for a power system according to any one of claims 1 to 7.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the steps of the nodal inertia identification method for a power system according to any one of claims 1 to 7.

Images