Nodal inertia identification method using data preprocessing and disturbance initial point determination

Abstract

The invention belongs to the field of power system data processing and relates to a node inertia identification method using data preprocessing and disturbance initial point determination, and the steps are: first, a quadratic penalty function is set, a data preprocessing method is performed based on the second-order penalty function, the range of the optimal solution is set, and the penalty function is used to remove the frequency oscillation data at the initial moment; secondly, the required identification data is retained as an optimal data set, and based on the optimal data set, an algorithm based on the disturbance initial point determination is added and a variable-order power function is fitted to improve the calculation accuracy; finally, a variable-order power function is fitted based on the disturbance initial point determination algorithm, so as to improve the identification accuracy under small and large disturbances, and the method is suitable for the inertia test scenario of multiple RoCoFs, and the accurate evaluation of the node inertia is realized.

Description

Technical Field

[0001] The invention belongs to the field of power system data processing, and relates to an equivalent inertia time constant identification method, in particular to a node inertia identification method using data preprocessing and disturbance initial point determination. Background Art

[0002] The large-scale access of new energy sources such as wind power and photovoltaic power has become an inevitable trend in the development of future power systems, but the large-scale access of new energy sources with weak / no inertia support capabilities has brought severe challenges to the safe operation of power systems. When a large disturbance event occurs, the system's ability to suppress the rate of change of frequency (RoCoF) becomes weaker, which is easy to trigger underfrequency load shedding, resulting in an increased risk of major power outages. In order to meet the requirements of the Paris Agreement to reduce greenhouse gas emissions by 40% by 2030, the global power system is in the midst of a dramatic transformation of decarbonization and decentralization, and it is imperative to build a new power system with new energy as the main body. A high proportion of new energy replaces traditional synchronous generators, resulting in a reduction in the overall inertia of the power system, which is prone to produce large frequency deviations and frequency change rates when the system fails, and is more likely to trigger the action of the Underfrequency Load Shedding (UFLS) relay or the RoCoF relay, thereby exacerbating the risk of system power outages. In order to ensure that the system as a whole or partially meets the inertia level for safe operation, it is of great significance to perform spatiotemporal identification of the power grid inertia and achieve overall optimal configuration of the system inertia.

[0003] Monitoring systems based on multiple PMUs have been widely studied and applied in the field of power system node inertia identification. Some scholars have proposed a synchronous machine inertia time constant estimation method based on an extended Kalman filter, but it requires immediate start-up at the moment of system disturbance, otherwise it will be affected by the system frequency oscillation and introduce large errors. In addition, some scholars have proposed fitting the system frequency curve based on the linear least squares method, and using differential equations based on the finite difference method to fit the data, reducing the adverse effects of transient disturbances and noise on RoCoF measurement, but not considering the influence of the time window length on inertia estimation. In order to avoid this influence, some scholars have used the system inertia dynamic evaluation method of the sliding data window, but it is only applied to the inertia evaluation of a single synchronous machine. When facing the synchronous measurement and identification of the equivalent inertia of multiple synchronous motor nodes, it is greatly affected by the selected sliding data window length. In order to improve the versatility of inertia identification in different power systems or power generation units, some scholars have adopted a first-order feedback control filtering method based on the limiting link to reduce the error of the RoCoF value at the initial moment by preprocessing the frequency data, but the preprocessing method also leads to the loss of effective data and reduces the accuracy of identification. Therefore, the existing system inertia identification method is greatly affected by the quality of PMU monitoring data and RoCoF calculation errors, and the data preprocessing effect that potentially causes inertia assessment errors needs to be improved.

[0004] In terms of the inertia identification method of system nodes, some scholars proposed a calculation method for the inertia time constant of system nodes based on the synchronous machine swing equation, and then obtained the spatiotemporal distribution of inertia of each node of the large system; some scholars further derived the inverse mathematical relationship between the propagation speed of electromechanical waves and inertia based on the electromechanical wave theory, and indirectly obtained the spatiotemporal distribution of inertia of the system. However, the accuracy of these two types of evaluation methods for inertia evaluation under different working conditions remains to be verified. In order to improve the accuracy of inertia identification, a distributed state estimation method can be used to achieve iterative optimization of the initial values ​​of the identification model parameters, but this method has the problems of high identification order and complex calculation process. If two sliding data windows are added to filter out useless high-frequency components in the data before and after the disturbance, the order of the identification model can be reduced and the inertia identification accuracy can be improved. Based on this, an inertia identification method based on fixed-order polynomial frequency fitting can be used to further reduce the oscillation component of RoCoF, but the identification accuracy is still subject to the reasonable setting of the disturbance occurrence time and the polynomial order. Therefore, the accuracy of the existing power system inertia identification method still has a lot of room for improvement.

[0005] In summary, there is an urgent need for an accurate identification method of inertia at power system nodes. Summary of the invention

[0006] The purpose of the present invention is to provide a node inertia identification method using data preprocessing and disturbance initial point determination to overcome the problems existing in the above-mentioned background technology.

[0007] The present invention solves the technical problem by the following technical solutions:

[0008] The node inertia identification method using data preprocessing and disturbance initial point determination is characterized by comprising the following steps:

[0009] Step 1: Collect discrete frequency data of each node in the power system in real time based on PMU;

[0010] Step 2, setting a second-order penalty function, performing data calculation based on the second-order penalty function on the frequency data collected by the PMU in step 1, and obtaining penalty function data corresponding to each frequency data;

[0011] Step 3, setting the inequality constraint condition M, and using the inequality penalty function and the penalty function data obtained in step 2 to eliminate the frequency oscillation data in the frequency data at the initial moment;

[0012] Step 4, retaining the data obtained after removing the frequency oscillation data in step 3 as the optimal data set Q;

[0013] Step 5, adding variable order power function fitting to the optimal data set in step 4;

[0014] Step 6: Based on step 5, an algorithm based on the initial point determination of disturbance is added to improve the high-order fitting algorithm at time 0+;

[0015] Step 7: Based on the above steps, calculate the inertia time constant identification value at each order on the basis of initialization, and output the optimal value.

[0016] Moreover, the second-order penalty function is:

[0017]

[0018] Where: σ is the penalty factor of the second-order function; P I (t,σ) represents the optimal solution corresponding to each frequency data; f(t) represents the value of the corresponding data under the initial fitting curve; c 2 i (t) is the square distance from the point represented by the monitoring data to the initial fitting curve.

[0019] Moreover, the method of performing data preprocessing based on a second-order penalty function on the frequency data collected by the PMU is:

[0020] Step 2.1, the frequency data collected by the PMU is fitted with a third-order power function through a third-order fitting function, and the third-order fitting function is:

[0021] f(t)=at 3 +bt2 +dt+c

[0022] Where: t is the data monitoring point extracted under the corresponding time window length, a, b, c, d are identification parameters;

[0023] Step 2.2, for any frequency data point at any time (t i ,f i ), t i is time, f i is the frequency, and the square distance from this point to the fitting curve f(t) is calculated by the following formula:

[0024]

[0025] Derivative of the above formula:

[0026]

[0027] The six real number solutions obtained by the derivation are opt Substitute the third-order fitting function to obtain six c 2 i (t opt ), where c 2 i The minimum value of (t) is the square distance that needs to be calculated;

[0028] Step 2.3: For all data within a time window, label them in chronological order {x1, x2…, x k}, the second-order penalty function corresponding to the point within the data window length is:

[0029]

[0030] Where σ is the penalty factor of the second-order function; P I (t,σ) represents the optimal solution corresponding to each frequency data; f(t) represents the value of the corresponding data under the initial fitting curve; c 2 i (t) is the square distance from the point represented by the monitoring data to the initial fitting curve.

[0031] Moreover, the method of setting the inequality constraint condition M and using the inequality penalty function and the calculated data obtained in step 2 to remove the frequency oscillation data in the frequency data at the initial moment is:

[0032] Step 3.1, select the point x corresponding to the frequency oscillation peak k , as the reference point of the constraint, and x k The square distance to the calculation is used as the inequality constraint for selecting the fitting data points, and then there is the inequality constraint condition M:

[0033]

[0034] Where: k Represents the penalty factor corresponding to the kth data; x k is the point corresponding to the peak value of frequency oscillation;

[0035] Step 3.2, define the inequality penalty function as:

[0036]

[0037] σ k+1 = rσ k

[0038] stc 2 i (x i )≤M, i=1, 2, ..., n

[0039] Where: min f(x) is the best linear fitting line identified in the first place, R represents a real number, and r is the expansion constant.

[0040] Moreover, the method of retaining the data obtained after removing the frequency oscillation data based on the inequality constraint M in step 3 as the optimal data set Q is:

[0041] Based on the quadratic penalty function, the N frequency data are preprocessed, and the optimal data set Q is selected as

[0042]

[0043] The optimal data set Q = {x1, x2, ..., x i} as the frequency fitting point and proceed to the next step of data identification.

[0044] Moreover, the method of adding variable-order power function fitting to the optimal data set in step 4 is:

[0045] After frequency data preprocessing, the points corresponding to the optimal data set Q of each time window length are taken as sample points {(x1, y1), (x2, y2), …, (x i ,y i )}, these sample points are fitted with a high-order polynomial:

[0046] y=a n x n +a n-1 x n-1 +a n-2 x n-2 +...+a1x+a0

[0047] Where: {a0, a1, ...a n} as a parameter.

[0048] Moreover, the method of adding an algorithm based on the initial point determination of disturbance to the optimal data set to improve the high-order fitting algorithm at the 0+ moment is:

[0049] In the process of high-order fitting, the perturbation initial point is added to determine the fitting strategy so that the fitting curve first satisfies the relationship between the sample points (x1, y1), that is:

[0050] y1=a n x1 n +a n-1 x1 n-1 +a n-2 x1 n-2 +...+a1x1+a0

[0051] y=a n x n +a n-1 x n-1 +a n-2 x n-2 +...+a1x+a0

[0052] On the basis of satisfying the above formula, a high-order fit is performed again, and differentiation is performed to obtain the RoCoF value, and then the inertia value is calculated:

[0053]

[0054] Where: is the rate of change of the polynomial at the initial moment.

[0055] Moreover, based on the above steps, the method of calculating the inertia time constant identification value of each order in sequence under initialization and outputting the optimal value is:

[0056] The measurement error η is defined as the inertia time constant identification value H identify With the setting value H ref Comparison of differences, that is:

[0057]

[0058] The inertia corresponding to the minimum error is the optimal output inertia.

[0059] The advantages and beneficial effects of the present invention are:

[0060] The present invention effectively reduces the influence of frequency oscillation on inertia identification accuracy by preprocessing frequency data, more accurately retains valid data in the inertia identification process, and optimizes data quality;

[0061] The present invention improves the influence of the high-order fitting algorithm at the 0+ moment on the system identification by adding a variable-order fitting algorithm under the initial point judgment of the disturbance, and improves the identification accuracy under small and large disturbances;

[0062] The present invention combines the typical power system topology structure and presents the spatiotemporal characteristics of power grid inertia in the form of a heat map, thereby realizing the visualization of inertia. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 is the node inertia response system diagram;

[0064] Figure 2 It is a data preprocessing flow chart;

[0065] Figure 3 It is a flow chart of node inertia identification based on variable order fitting under second order penalty function and disturbance initial point determination;

[0066] Figure 4 It is a 10-machine 39-node system that deploys PMU. DETAILED DESCRIPTION

[0067] The present invention is further described in detail below through specific examples. The following examples are only illustrative and not restrictive, and the protection scope of the present invention cannot be limited thereto.

[0068] See also Figures 1 to 3 The node inertia identification method (temporal and spatial Power Grid Inertia Perception, PGIP) using data preprocessing and disturbance initial point determination mainly includes the following three parts: (1) real-time acquisition of discrete power frequency data points of each node in the power system; (2) variable-order fitting based on data preprocessing and disturbance initial point determination; (3) identification of equivalent inertia time constant. First, based on the data preprocessing method of the second-order penalty function, inequality constraints are set, and the penalty function is used to eliminate the frequency oscillation data at the initial moment; secondly, the required identification data is retained as the optimal data set, and based on the optimal data set, an algorithm based on disturbance initial point determination is added to improve the influence of the high-order fitting algorithm at the 0+ moment on the system identification, and improve the calculation accuracy; finally, based on the algorithm for disturbance initial point determination, variable-order power function fitting is performed, which improves the identification accuracy under small and large disturbances, and is suitable for the inertia test scenario of multiple RoCoFs.

[0069] The node inertia identification method using data preprocessing and disturbance initial point determination is innovative in that it includes the following steps:

[0070] Step 1: Collect discrete frequency data of each node in the power system in real time based on PMU;

[0071] Step 2, setting a second-order penalty function, and performing data preprocessing based on the second-order penalty function on the frequency data collected by the PMU;

[0072] The second-order penalty function is:

[0073]

[0074] Where: σ is the penalty factor of the second-order function; P I (t,σ) represents the optimal solution corresponding to each data; f(t) represents the value of the corresponding data under the initial fitting curve; c 2 i (t) is the square distance from the point represented by the monitoring data to the initial fitting curve;

[0075] Step 3, set the inequality constraint condition M, and use the data in step 2 and the inequality penalty function to eliminate the frequency oscillation data in the frequency data at the initial moment;

[0076] Step 4, retaining the data after removing the frequency oscillation data as the optimal data set;

[0077] Step 5: Add an algorithm based on the initial point of disturbance to the optimal data set to improve the high-order fitting algorithm at the 0+ moment and obtain the RoCoF value at the initial moment of disturbance;

[0078] Step 6: Based on step 5, a variable-order power function is added to fit the RoCoF value at each node under each order power function;

[0079] Step 7, based on the above steps, initialize the required fitting data, the initial fitting order n is 1, the initial error η0 is 100%, the initial inertia H0 is 1, calculate the inertia time constant identification value under each order, and output the optimal value.

[0080] 1. Power system node inertia

[0081] Typically, the inertia of a synchronous machine is defined as the resistance of its rotating mass to speed changes. In contrast, for a grid node, its equivalent inertia can be contributed by the inertia of each generator set. After the power system suffers a power loss (e.g., a synchronous machine goes offline after a fault), the frequency response with different time scales is as follows: Figure 1 As shown in Figure 1, in the initial stage (usually within 10 seconds) when the system generates a high RoCoF, the inertia response plays a dominant role in stabilizing the system frequency. In this stage, the rotational kinetic energy stored in the synchronous machine rotor is released, which reduces the transient unbalanced power of the system and prevents the system frequency from falling to a certain extent.

[0082] The traditional power system inertia identification method fits the inertia time constant based on the frequency-power data measured by the PMU, and is mainly aimed at the inertia identification of synchronous machines or power grid areas. The main steps are as follows:

[0083] Obtain the electromagnetic active power {P e,1 ,P e,2 …,P e,i}, and measure the corresponding frequency change rate f through PMU i , thus obtaining a single-input single-output model between power and frequency. Therefore, the rotor motion equation of the stepper i node is written as:

[0084]

[0085] Where: f i is the frequency of node i of the synchronous machine; ΔP u,i is the mechanical power change of synchronous machine i; ΔP e,i is the change in electromagnetic power of synchronous machine i; H i is the equivalent inertia of synchronous machine i; D i is the equivalent damping of synchronous machine i; Δf i is the speed offset of synchronous machine i.

[0086] After Laplace transformation of equation (1), we get:

[0087]

[0088] By fitting the first-order model of the RoCoF, the inertia time constant can be identified. If the power change of synchronous machine i is set to a specific value k i , then the denominator is a linear function. By performing an inverse Laplace transform on equation (2) and taking the slope change rate at time t = 0, the inertia time constant of synchronous machine i can be identified:

[0089]

[0090] The method of setting a second-order penalty function and performing data preprocessing based on the second-order penalty function on the frequency data collected by the PMU is as follows:

[0091] First, define a second-order penalty function:

[0092]

[0093] Where: t is the data monitoring point extracted under the corresponding time window length; σ is the penalty factor of the second-order function; P I (t,σ) represents the optimal solution corresponding to each data; f(t) represents the value of the corresponding data under the initial fitting curve; c 2i (t) is the square distance from the point represented by the monitoring data to the initial fitting curve. The square term of the second-order penalty function can reduce the amount of square root calculation of the algorithm and improve the speed of the algorithm.

[0094] Within a fixed data window length, the c obtained at the non-oscillation point 2 i (t) is slightly different from the distance to the actual frequency curve, and its error distance is within the error tolerance range of the penalty function to remove data. Therefore, in the initial calculation of c 2 i (t) The third-order fitting function is used in the process:

[0095] f(t)=at 3 +bt 2 +dt+c (5)

[0096] Where: a, b, c, d are identification parameters.

[0097] For any data point (t i ,f i ),t i is time, f i is the frequency. From formula (5), we can know that the square distance from this point to the fitting curve f(t) is:

[0098]

[0099] Derivative of the formula:

[0100]

[0101] Solve the six real numbers obtained from equation (7) t opt Substitute the fitting curve f(t) into the equation and get six c2 i(t opt ), where the minimum value of c2 i(t) is the square distance to be calculated.

[0102] For all data within a time window, they are labeled {x1, x2…, x k}, from the second-order penalty function, we can know that the second-order penalty function corresponding to the point within the data window length is:

[0103]

[0104] Define the inequality penalty function:

[0105]

[0106] σ k+1 = rσ k (10)

[0107] stc 2 i (x i )≤M,i=1,2,...,n (11)

[0108] Where: minf(x) is the best linear fitting line identified in the first place. k is a given initial point, selected as the point corresponding to the peak value of the frequency oscillation. k Represents the penalty factor corresponding to the kth data, and r is the expansion constant.

[0109] Due to the suppression of inertia after a fault, the frequency fluctuation will gradually decrease. Therefore, for each cycle of fluctuation, the corresponding penalty function standard is different. The fluctuation amplitude can be increased by expanding the constant, and the standard of the penalty function can be reduced, thereby increasing the accuracy of data preprocessing.

[0110] x k The square of the distance to the linear function is used as the inequality constraint for selecting better data points, and then there is the inequality constraint condition M:

[0111]

[0112] The data set Q selected by data preprocessing based on formula (12) is

[0113]

[0114] The screened data set Q = {x1, x2, ..., x i} as the frequency fitting point, and then proceed to the next step of data identification. The data preprocessing flow chart is shown in Figure 2.

[0115] 3. RoCoF identification based on high-order fitting under disturbance initial point determination

[0116] After data preprocessing, the points corresponding to the data set Q of each time window length are taken as sample points {(x1, y1), (x2, y2), …, (x i ,y i )}. Fit these sample points with a high-order polynomial:

[0117] y=a n x n +a n-1 x n-1 +a n-2 x n-2 +...+a1x+a0 (14)

[0118] By differentiating equation (14), we can get the RoCoF value, and then calculate the inertia value:

[0119]

[0120] It can be seen from the above formula that the determination of the value of a1 is of great significance to the accuracy of inertia identification, and the data of the sample point (x1, y1) at time t = 0 has a greater impact on a1. Therefore, in the process of high-order fitting of formula (15), the perturbation initial point determination fitting strategy is added to make the fitting curve first satisfy the relationship between the sample points (x1, y1), that is:

[0121] y1=a n x1 n +a n-1 x1 n-1 +a n-2 x1 n-2 +...+a1x1+a0 (16)

[0122] Within a time window, the perturbation initial point determination fitting strategy adapts the a0 value to the initial state of the frequency change. In order to obtain the optimal parameters for high-order fitting, the square of the residual is used as an indicator for evaluating the fitting accuracy:

[0123]

[0124] In order to minimize the fitting coefficient ε, the following n+1 partial derivatives of ε are calculated and each partial derivative is set to 0:

[0125]

[0126] When a0 is determined, equation (18) can be expanded into:

[0127]

[0128] The remaining n fitting coefficients obtained from the n equations in formula (19) are the optimal solution that minimizes ε, which can fit the measured frequency curve under disturbance to the greatest extent. Based on the fitting curve, RoCoF can be calculated by formula (15).

[0129] The measurement error η is defined as the inertia time constant identification value H identify With the setting value H ref Comparison of differences, that is:

[0130]

[0131] Where: η is the measurement error, H identify is the inertia time constant identification value, H ref is the set value.

[0132] 4. Method implementation process

[0133] The overall measurement process of the PGIP inertia evaluation method is as follows: Figure 3As shown, it mainly includes the following three steps:

[0134] 1) Use the PMU measurement device to measure the power and frequency data of each node when a disturbance or failure occurs in the system;

[0135] 2) The frequency data obtained by PMU monitoring is preprocessed by equations (5) to (13), and high-order power functions are fitted by equations (16) to (19). Based on equations (14) to (15), the RoCoF values ​​obtained under various order power functions at the initial moment of the disturbance are obtained.

[0136] 3) Initialize the data (initial fitting order n is 1, initial error η0 is 100%, initial inertia H0 is 1), and then use formula (3) to calculate the inertia time constant identification value at each node of the power system under PGIP. Compare this value with the fitting error of the previous high-order power function: if the relative error is small, select the inertia value obtained by fitting under the new first-order power function; if the relative error is large, abandon it and enter the next cycle (verified by the example, when the fitting order n is higher than 10th order in the case of disturbance, the fitting curve of the oscillation frequency remains almost unchanged, so the present invention sets the 10th order fitting as the termination condition).

[0137] 4. Method performance verification

[0138] An IEEE 39-node system was built in the PSCAD simulation platform, in which the power and frequency at each synchronous machine node were measured by the PMU deployed therein. The accuracy of inertia identification and the effectiveness of spatiotemporal perception of the proposed PGIP were verified and demonstrated in the MATLAB platform. The present invention compares the performance of the traditional inertia identification method based on the least squares method, the ARMAX inertia identification algorithm, the inertia identification algorithm containing only second-order penalty function data preprocessing, and the PGIP algorithm proposed in the present invention through simulation examples. Figure 4 It is a 10-machine 39-node system. The base capacity of the system of the present invention is selected as 500MVA, and the rated frequency is 60Hz.

[0139] Error evaluation of inertia identification algorithm under small disturbance

[0140] Except for the PGIP method proposed in the present invention, the other three methods all have significant errors. If only the data preprocessing method based on the second-order penalty function is used to eliminate the data of frequency oscillation, the overall identification error under small disturbance can be controlled at about 11% when the time window length is 4s, while the errors of the least square method and ARMAX are about 25%. On the basis of the data preprocessing of the second-order penalty function, PGIP further determines the frequency change at the precise initial moment by the initial point of the disturbance, improves the influence of the high-order fitting algorithm of the initial point at 0+ on the system identification, and further reduces the identification error of the data preprocessing method (when the time window length is 4s) to 5.1%.

[0141] At the same time, since the value range of the inertia time constant of the synchronous machine is mostly 3.5s to 4s, as the time window length is closer to the inertia time constant of the specific synchronous machine, the error of the results obtained by the four methods is lower. This means that according to the synchronous machine motor of the actual system, it is particularly important to reasonably select a time window length that is closer to the inertia time constant of the synchronous machine to ensure the identification accuracy.

[0142] Error evaluation of inertia identification algorithm under large disturbance

[0143] Compared with the small disturbance example, the identification accuracy of the ARMAX algorithm is significantly improved under large disturbance (with the same synchronous machine number and time window length), but the inertia identification result based on the PGIP algorithm has the lowest error, which is slightly better than the inertia identification method based on the ARMAX algorithm and the data preprocessing, and significantly better than the inertia identification method based on the least squares method.

[0144] Taking the synchronous machine G1 as an example, when the time window length is selected as 4s, the G1 inertia identification error based on the PGIP algorithm is only 0.7%, which is significantly improved compared with the identification error of 14.75% under the traditional least squares method, 2.1% under the ARMAX algorithm, and 2.95% under data preprocessing.

[0145] Although the embodiments and drawings of the present invention are disclosed for illustrative purposes, those skilled in the art will appreciate that various substitutions, changes and modifications are possible without departing from the spirit and scope of the present invention and the appended claims. Therefore, the scope of the present invention is not limited to the contents disclosed in the embodiments and drawings.

Claims

1. A node inertia identification method using data preprocessing and disturbance initial point determination is characterized by: The steps include: Step 1: Collect discrete frequency data of each node in the power system in real time based on PMU; Step 2, setting a second-order penalty function, performing data calculation based on the second-order penalty function on the frequency data collected by the PMU in step 1, and obtaining penalty function data corresponding to each frequency data; Step 3, setting the inequality constraint condition M, and using the inequality penalty function and the penalty function data obtained in step 2 to eliminate the frequency oscillation data in the frequency data at the initial moment; Step 4, retaining the data obtained after removing the frequency oscillation data in step 3 as the optimal data set Q; Step 5, adding variable order power function fitting to the optimal data set in step 4; Step 6: Based on step 5, an algorithm based on the initial point determination of disturbance is added to improve the high-order fitting algorithm at time 0+; Step 7: Based on the above steps, calculate the inertia time constant identification value at each order on the basis of initialization, and output the optimal value.

2. The node inertia identification method using data preprocessing and disturbance initial point determination according to claim 1 is characterized in that: The second-order penalty function is: Where: σ is the penalty factor of the second-order function; P I (t,σ) represents the optimal solution corresponding to each frequency data; f(t) represents the value of the corresponding data under the initial fitting curve; c 2 i (t) is the square distance from the point represented by the monitoring data to the initial fitting curve.

3. The node inertia identification method using data preprocessing and disturbance initial point determination according to claim 2 is characterized in that: The method of performing data calculation based on the second-order penalty function on the frequency data collected by the PMU in step 1 to obtain the penalty function data corresponding to each frequency data is: Step 2.1, the frequency data collected by the PMU is fitted with a third-order power function through a third-order fitting function, and the third-order fitting function is: f(t)=at 3 +bt 2 +dt+c Where: t is the data monitoring point extracted under the corresponding time window length, a, b, c, d are identification parameters; Step 2.2, for any frequency data point at any time (t i ,f i ), t i is time, f i is the frequency, and the square distance from this point to the fitting curve f(t) is calculated by the following formula: Derivative of the above formula: The six real number solutions obtained by the derivation are opt Substitute the third-order fitting function to obtain six c 2 i (t opt ), where c 2 i The minimum value of (t) is the square distance that needs to be calculated; Step 2.3: For all data within a time window, label them in chronological order {x1, x2…, x k }, the second-order penalty function corresponding to the point within the data window length is: Where σ is the penalty factor of the second-order function; P I (t,σ) represents the corresponding optimization solution under each frequency data; f(x) represents the value of the corresponding data under the initial fitting curve; c2i(t) is the square distance from the point represented by the monitoring data to the initial fitting curve.

4. The node inertia identification method using data preprocessing and disturbance initial point determination according to claim 1 is characterized in that: The method of setting the inequality constraint condition M and using the inequality penalty function and the calculated data obtained in step 2 to remove the frequency oscillation data in the frequency data at the initial moment is: Step 3.1, select the point x corresponding to the frequency oscillation peak k , as the reference point of the constraint, and x k The square distance to the calculation is used as the inequality constraint for selecting the fitting data points, and then there is the inequality constraint condition M: Where: k Represents the penalty factor corresponding to the kth data; x k is the point corresponding to the peak value of frequency oscillation; Step 3.2, define the inequality penalty function as: s k+1 =rσ k s.t.c 2 i (x i )≤M,i=1,2,...,n Where: minf(x) is the best linear fitting line identified initially, R represents a real number, and r is the expansion constant; The method of retaining the data obtained after removing the frequency oscillation data based on the inequality constraint M in step 3 as the optimal data set Q is: Based on the quadratic penalty function, the N frequency data are preprocessed, and the optimal data set Q is selected as The optimal data set Q = {x1, x2, ..., x i } as the frequency fitting point and proceed to the next step of data identification.

5. The node inertia identification method using data preprocessing and disturbance initial point determination according to claim 1 is characterized in that: The method of adding variable-order power function fitting to the optimal data set in step 4 is: After frequency data preprocessing, the points corresponding to the optimal data set Q of each time window length are taken as sample points {(x 1， y1),(x2,y2),…,(x i， y i )}, these sample points are fitted with a high-order polynomial: y=a n x n +a n-1 x n-1 +a n-2 x n-2 +...+a1x+a0 Where: {a0, a1, ...a n } as a parameter.

6. The node inertia identification method using data preprocessing and disturbance initial point determination according to claim 1 is characterized in that: The method of adding an algorithm based on the initial point determination of disturbance based on step 5 to improve the high-order fitting algorithm at time 0+ is: In the process of high-order fitting, the perturbation initial point is added to determine the fitting strategy so that the fitting curve first satisfies the relationship between the sample points (x1, y1), that is: y1=a n x1 n +a n-1 x1 n-1 +a n-2 x1 n-2 +...+a1x1+a0 y=a n x n +a n-1 x n-1 +a n-2 x n-2 +...+a1x+a0 On the basis of satisfying the above formula, a high-order fit is performed again, and differentiation is performed to obtain the RoCoF value, and then the inertia value is calculated: Where: is the rate of change of the polynomial at the initial moment.

7. The node inertia identification method using data preprocessing and disturbance initial point determination according to claim 1 is characterized in that: Based on the above steps, the method of calculating the inertia time constant identification value of each order in sequence under initialization and outputting the optimal value is: The measurement error η is defined as the inertia time constant identification value H identify With the setting value H ref Comparison of differences, that is: The inertia corresponding to the minimum error is the optimal output inertia.

Images