A method for identifying equivalent inertia constants of power systems based on data regression

Abstract

The present invention discloses a method for identifying equivalent inertia constants of a power system based on data regression, namely, the least square method (LSM), which belongs to the technical field of power systems. The present invention includes identification of equivalent inertia constants of synchronous machines and identification of virtual inertia; the moment of inertia is an inherent property of the synchronous generator SG, which stores kinetic energy through the rotor, and the SG naturally uses this kinetic energy to compensate for power imbalance, as shown in the formula: #imgabs0# where ω SG is the angular frequency of the synchronous motor, p m,SG and p e,SG is the mechanical power and electromagnetic power; p D,SG The active power change is caused by the damping effect. The present invention can achieve good accuracy in identifying the moment of inertia and virtual inertia under different operating conditions through the cooperation of various parts, and can be used to quantify the inertia support effect of time-varying adaptive inertia.

Description

Technical Field

[0001] The invention belongs to the technical field of power systems, and in particular relates to a method for identifying equivalent inertia constants of power systems based on data regression. Background Art

[0002] On the one hand, in power systems, a method will be needed to accurately detect and identify the equivalent inertia constant of devices connected to the grid. In recent years, various online inertia identification methods have been proposed. These methods are based on the assumption of the availability of phasor measurement unit (PMU) data, which are currently widely used by transmission system operators. However, these methods require events to occur at the node that trigger significant frequency changes. Another feature of these methods is that they focus on obtaining fast estimates, but they may not be accurate depending on the transient estimates.

[0003] On the other hand, the inertia market requires more accurate inertia information under different operating conditions rather than instantaneous estimates of its inertial response. In this regard, inertia identification methods based on PMU measurement data have been proposed. However, the accuracy of these methods is highly dependent on the accuracy of the model parameters of the estimation equipment, which need to be adjusted case by case. To this end, a method for identifying the equivalent inertia constant of a power system based on data regression is designed. Summary of the invention

[0004] The purpose of the present invention is to provide a method for identifying equivalent inertia constants of a power system based on data regression to solve the problems raised in the above background technology.

[0005] In order to solve the above technical problems, the present invention is achieved through the following technical solutions:

[0006] The invention discloses a method for identifying equivalent inertia constant of a power system based on data regression, which is used for identifying equivalent inertia constant and virtual inertia.

[0007] Furthermore, the moment of inertia is a fundamental property of the synchronous motor SG. By storing kinetic energy in the rotor, the SG naturally uses this kinetic energy to compensate for power imbalance, as shown in the formula:

[0008]

[0009] in: is the angular frequency of the synchronous motor, p m,SG and p e,SG is the mechanical power and electromagnetic power; p D,SG is the change in active power due to the damping effect, and:

[0010] p D,SG =D SG (ω Grid -ωSG )

[0011] Where: D SG is the damping coefficient; M SG is the mechanical starting time of the synchronous motor SG in seconds;

[0012] D SG <<M SG Always applies to synchronous machines; start-up time M SG is related to the machine's inertia constant as shown in the formula:

[0013] M SG =2H RI

[0014] Where: H RI is the equivalent inertia constant.

[0015] Furthermore, for p m and p D,SG The influence of p m Split into three components:

[0016] p m =p UC +p PFC +p SFC

[0017] Where: p UC The power set point obtained by solving the unit commitment problem; p PFC and p SFC Active power regulated by primary frequency control (PFC) and secondary frequency control (SFC), respectively;

[0018] Then And substitute p m =p UC +p PFC +p SFC And differentiate it, we get:

[0019]

[0020] in: is the rate of change of frequency, i.e. RoCoF, obtained by PMU measurement; the time derivative of RoCoF, i.e. Obtained through PMU measurement using a PI filter;

[0021] In the formula, and is always true, i.e. in the time window that characterizes the inertial response of the generator; therefore, for an appropriate time window, we obtain:

[0022]

[0023] Further, as shown in the following formula:

[0024] P e,IBR =P UC +P PFC +P SFC +P VI -P D,IBR

[0025] Where: P VI is the power increased by the virtual inertia:

[0026]

[0027] The formula shows that H VI pass Since the virtual inertial support is the result of the control action, H can be adaptively adjusted by appropriately measuring the system state. VI , in order to better deal with frequency deviation problems after sudden unexpected events.

[0028] Furthermore, based on the LMS least squares method, the inertia estimation formula needs to be modified to cope with the measurement error: For the i The measurement data obtained are:

[0029]

[0030] Because of its fractional structure, it can cause serious numerical problems; when the numerator When the denominator is very small, but not 0, the effect will be significantly amplified due to the infinitesimal denominator, resulting in a peak, causing H * There are serious errors in the estimation of Rewritten as:

[0031]

[0032] For a set of measurements t1, t2, ..., t N-1 ,t N Define the following vector:

[0033]

[0034]

[0035] Where: T represents the matrix transpose; A and b are collected and The matrix of data, time window [t1,t N ], the data size N of each matrix is ​​determined by the length of the time window T = t N-t1 and PMU sampling rate.

[0036] Furthermore, the least square method is used to find the equivalent inertia constant H, where H is the inertia constant of the synchronous machine:

[0037] Minimize(AH * -b) 2

[0038] Transform the above formula:

[0039] H * =(A T A) -1 A T b

[0040] Formula H * =(A T A) -1 A T b is the expression used to determine the equivalent inertia constant, taking care to avoid pathological A T A, it is convenient to adjust its diagonal elements with small diagonal elements. In the tests solved in preparation for this work, ε = 10 -6 To prevent A T A is singular.

[0041] The present invention has the following beneficial effects:

[0042] 1. The present invention, through the coordination of a series of parts, can identify the moment of inertia and virtual inertia with good accuracy under different operating conditions, and can be used to quantify the inertia support effect of time-varying adaptive inertia.

[0043] 2. The present invention can be applied to various operating conditions, including stable operating conditions, and is independent of the device and can be used to quantify the virtual inertial response of asynchronous devices. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings required for describing the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other accompanying drawings can be obtained based on these accompanying drawings without paying creative work.

[0045] Figure 1 This is a schematic diagram of the frequency response of the power system with high renewable energy penetration;

[0046] Figure 2 This is a schematic diagram of the WSCC9 bus system of the present invention;

[0047] Figure 3It is a schematic diagram of the SG dynamic trajectory of the WSCC9 bus system of the present invention under normal operating conditions;

[0048] Figure 4 It is a schematic diagram of the identification result of the equivalent inertia constant of the SG under normal operating conditions of the present invention;

[0049] Figure 5 It is a schematic diagram of SG trajectory after the load of the WSCC9 bus system of the present invention is suddenly increased;

[0050] Figure 6 It is a schematic diagram of the identification result of the equivalent inertia constant of SG after a sudden load increase of the present invention;

[0051] Figure 7 Schematic diagram of the trajectory of the WSCC9 bus system for different virtual inertial providers;

[0052] Figure 8 Schematic diagram of virtual inertia technology based on VSG inertia identification results of the present invention based on the least squares method. DETAILED DESCRIPTION

[0053] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0054] See also Figure 1-8 As shown, the present invention is a method for identifying equivalent inertia constants of power systems based on data regression. It includes:

[0055] 1. Calculation of equivalent inertia constant.

[0056] 1.1 Moment of inertia

[0057] Moment of inertia is a fundamental property of the SG, which stores kinetic energy through the rotor. The SG naturally uses this kinetic energy to compensate for power imbalance, as shown in the formula:

[0058]

[0059] in: is the angular frequency of the synchronous motor, p m,SG and p e,SG is the mechanical power and electromagnetic power; p D,SG is the change in active power due to the damping effect, and:

[0060] p D,SG =D SG (ωGrid -ω SG )(2)

[0061] Where: D SG is the damping coefficient; M SG is the mechanical start-up time of the SG, in seconds;

[0062] D SG <<M SG Always applies to synchronous motors; starting time M SG is related to the machine's inertia constant as shown in the formula:

[0063] M SG =2H RI (3)

[0064] Where: H RI is the equivalent inertia constant.

[0065] In formula (1), p e,SG Can be measured directly; ω SG Can be obtained through PMU measurement; m It cannot be measured directly. Therefore, we need to solve H by formula (1): RI , for p m and p D,SG The impact of; provides a good solution, first of all, p m Split into three parts:

[0066] p m =p UC +p PFC +p SFC (4)

[0067] Where: p UC The power reference point given by the current optimization scheduling result; p PFC and p SFC Specifically, it is the active power regulated by primary frequency control (PFC) and secondary frequency control (SFC).

[0068] Then, substitute equation (4) into (1) and differentiate to obtain:

[0069]

[0070] in: is the power change rate, i.e. RoCoF, which can be obtained by PMU measurement; the time derivative of RoCoF, i.e. is the time derivative of RoCoF and can be obtained by PMU measurement using a PI filter.

[0071] In formula (5), under the transient time scale, it can be considered as In practical applications, the power system is constantly subject to small disturbances caused by loads and high-penetration renewable energy power systems with randomness. In this case, the suppression of the power system is mainly through inertial response. In this case, since the frequency deviation is within the PFC dead zone, the frequency control does not work, so For synchronous generators, there is always is always satisfied, so for an appropriate time window, one can rewrite the above as:

[0072]

[0073] With the continuous increase in the penetration rate of renewable energy, SG has been greatly reduced, and the power system has insufficient rotational inertia. In the face of this problem, virtual inertia control and virtual synchronous machines use power electronic control and power-type energy storage to simulate rotational inertia to provide emergency frequency assistance for the system. The inertia response characteristics of the simulated synchronous machine enable new energy and energy storage to provide inertia support for the system. For the above formula, both rotational inertia and virtual inertia are applicable.

[0074] For this purpose, the above formula can be rewritten in a more general form, removing the subscript index:

[0075]

[0076] Where: the superscript * indicates an estimated value, and ξ indicates the measurement error. Although equation (7) is completely correct in theory, it cannot be reasonably applied in practice because its fractional structure will cause serious numerical problems. When the denominator is very small, but not 0, the effect will be significantly amplified due to the infinitesimal denominator, resulting in a peak, causing H * The estimation of has serious error problems. Therefore, based on the existing theory, the application of the least square method (LSM) is proposed for inertial identification based on the above formula and measurement of appropriate time windows, aiming to improve its error problem.

[0077] According to formula (7), for i The measured time series data obtained at the moment can be rewritten as:

[0078]

[0079]

[0080] For a set of measured values ​​t1, t2, ..., t N-1 , t N , we can define the following vectors:

[0081]

[0082] Where: T represents the transpose of the matrix; A and b are in the time window [t1, t N ] collected in the time series data and The matrix vector is composed of the following: where the data length N of each matrix is ​​determined by the length of the time window τ = t N -t1 and the sampling rate of PMU.

[0083] Therefore, the following least square method (LSM) can be proposed to effectively estimate the equivalent inertia constant:

[0084] Minimize(AH * -b) 2 (11)

[0085] H * =(A T A) -1 A T b(12)

[0086] This method can effectively analyze the equivalent inertia constant over a period of time, and assist in the pricing of inertia / fast frequency services.

[0087] 2. Least Squares Inertia Constant Identification

[0088] 2.1 Inertial identification time window

[0089] Figure 1 The frequency response of the power system with high renewable energy penetration is shown in .

[0090] exist Figure 1 In, T inertia Represents the inertial response period after the frequency event; T PFC and T SFC are the time scales of primary and secondary FM; the red rectangles are suitable time windows for identifying inertia. The rationale behind these time windows is as follows.

[0091] 2.1.1. Normal operating conditions with continuous small disturbances: In practice, the power system is continuously subject to small disturbances caused by random changes in load and the randomness of renewable energy. Under normal operating conditions, the power system mainly suppresses these random dynamic synchronous motors through the filtering effect of inertial response. In this case, since the frequency deviation is within the dead band of frequency control, In addition, the time window used in the inertial estimation is long enough to make the proposed estimation robust to environmental noise and other fast fluctuations.

[0092] 2.1.2. Inertial response period after PFC dead-zone frequency event: The characteristics of this time window are and is larger, so the effect of measurement noise on inertial identification can be mitigated compared to normal situations. Figure 1 It shows that the inertial response and the conventional PFC period partially overlap. However, the time scale of PFC, especially by non-synchronous devices, cannot be accurately defined because there are slow, fast and very fast PFC. As a rule of thumb, the estimation time window must be long enough for accurate inertial identification, while short enough to avoid the influence of conventional PFC.

[0093] Fast PFC may overlap with the estimation and therefore be considered as "virtual inertia". Finally, very fast PFC will be filtered out by the technique.

[0094] Note also that we assume that the estimation time window is shorter than the dispatch period of the generators. In any case, the transmission system operators know when the dispatch changes occur, so they can discard the estimates at that time.

[0095] 2.2 Problems with existing formulas

[0096] Formula (6) shows that the equivalent inertia constant can be determined by PMU measurement. In practical applications and when considering real data, serious numerical problems will arise due to its fractional structure. When the denominator is very small, but not 0, the effect will be significantly amplified due to the infinitesimal denominator, resulting in a peak, causing H * The estimation of has serious error problems. Formula (6) can be rewritten as:

[0097]

[0098] The superscript * indicates the measured / estimated value and ξ indicates the measurement error.

[0099] Although mathematically correct, formula (13) cannot be applied as is.

[0100] 2.3. Least Squares Method Based on Data

[0101] The proposed least square method (LSM) for inertia identification based on equation (13) and an appropriate measurement time window is introduced.

[0102] According to formula (13), for i The measurement data obtained are:

[0103]

[0104]

[0105] For a set of measurements t1, t2, ..., t N-1 ,tN The following vectors can be defined:

[0106]

[0107]

[0108] Where: T represents matrix transpose; A and b are the collected data and The matrix composed of time windows [t1,t N ], the data size N of each matrix is ​​determined by the length of the time window T = t N -t1 and PMU sampling rate.

[0109] We propose to use the following least squares method (LSM) to find H * :

[0110] Minimize(AH * -b) 2 (18)

[0111] You can rewrite it as follows:

[0112] H * =(A T A) -1 A T b(19)

[0113] Formula H * =(A T A) -1 A T b is the expression used to determine the equivalent inertia constant, taking care to avoid pathological A T A, it is convenient to adjust its diagonal elements with small diagonal elements. In the tests solved in preparation for this work, ε = 10 -6 To prevent A T A is singular.

[0114] 3. Case

[0115] This section verifies the proposed least square method (LSM) technique on a modified version of the WSCC9 bus system. Figure 2 As shown, the SGs connected to buses 1 and 3 are replaced by virtual inertia providers, namely virtual synchronous generators (VSGs) and wind power plants (WPPs) that provide small inertia to the system. The modified system has 75.6% non-synchronous generation.

[0116] The bus frequency and output active power of each inertial provider are measured using a PMU with a sampling rate of 50 Hz. and is obtained using a PI filter based technique. All simulations presented in this case study were obtained using the Python based software tool Dome.

[0117] 3.1 Identification of moment of inertia

[0118] SG's H SG =23.64MW S / MVA,D SG = 1.0 MW / MVA. The noise is modeled as an Ornstein-Uhlenbeck random process and added to each measurement. The dynamic wind speed of the WPP is modeled as a Weibull distribution with an exponentially decaying autocorrelation.

[0119] Under normal operating conditions, the power system is subject to continuous wind and load fluctuations. For example, for a 10-second time window, Figure 3 The dynamic trajectory of SG for the modified WSCC9 bus system under normal operating conditions is shown. The assumptions made for are well met in this case.

[0120] The identification results of the equivalent inertia constants of SG with different time window lengths (denoted as T) are as follows: Figure 4 According to the figure, the proposed inertia identification technology is relatively accurate under normal operating conditions with a time window greater than 3s, and the relative error distribution is in the range of [0.7, 2.51]%.

[0121] Next, we consider the transient response of the modified WSCC9 bus system after a large disturbance, i.e., a sudden increase in the load on bus 5 occurring at t = 0.1 s. The transient behavior of SG after the disturbance is shown in Figure 5 shown.

[0122] Figure 6 The results of SG inertia identification for different time windows after the disturbance occurs are shown. The time window of length [0.3, 15] s is long enough to identify the equivalent inertia constant of SG with satisfactory accuracy by the proposed least square method (LSM), i.e., the relative error is less than 1%. Figure 5 It also shows that Figure 3 Compared to the normal operating case shown, the influence of measurement noise is less pronounced and therefore the accuracy of inertial identification is improved.

[0123] 3.2. Virtual Inertia Identification

[0124] We consider different virtual inertia techniques for inertia identification of the VSG on bus 1. For simplicity, noise is not considered in the simulations discussed below.

[0125] For a constant virtual inertia, H VI =20MW S / MVA, the adaptive virtual inertia has the same H under normal operating conditions VI , and varies with the RoCoF of the frequency event. Figure 7 shows the trajectory of the VSG after a sudden load increase.

[0126] in, Figure 7 The upper middle picture is H * The figure below shows the trajectory of COI frequency, that is, the center of inertia frequency. The relative error is 0.1%.

[0127] Adaptive inertia allows additional support through constant inertia control and slightly improves the overall frequency response of the system.

[0128] Figure 8 The identification results of the constant and adaptive virtual inertia in a 15-second time window after a sudden load increase are shown.

[0129] IV. Conclusion

[0130] The present invention proposes a least square method (LSM) problem that can solve the numerical problems affecting such techniques. The proposed method is sufficient to estimate the equivalent inertia constant and virtual inertia. Simulation results show that the proposed inertia identification method has high accuracy.

[0131] The proposed technique can also be used to quantify the inertial support from time-varying inertial control. Therefore, the proposed method is a useful tool to provide information about the availability of inertia to the inertia market. Future work will be devoted to improving the estimation based on the least squares method (LSM), for example using the methods of Tikhonov or Lassos.

[0132] The preferred embodiments of the present invention disclosed above are only used to help explain the present invention. The preferred embodiments do not describe all the details in detail, nor do they limit the present invention to the specific implementation methods described. Obviously, many modifications and changes can be made according to the content of this specification. This specification selects and specifically describes these embodiments in order to better explain the principles and practical applications of the present invention, so that those skilled in the art can understand and use the present invention well. The present invention is limited only by the claims and their full scope and equivalents.

Claims

1. A method for identifying equivalent inertia constants of power systems based on data regression, characterized in that: Used to identify equivalent inertia constants and virtual inertia; The moment of inertia is a fundamental property of a synchronous motor SG. By storing kinetic energy in the rotor, the synchronous motor SG naturally uses this kinetic energy to compensate for power imbalance, as shown in the formula: in: is the first-order time derivative of the angular frequency of the synchronous motor, p m,SG and p e,SG is the mechanical power and electromagnetic power; p D,SG is the active power change due to the damping effect, M SG is the mechanical starting time of the synchronous motor SG in seconds; and: p D,SG =D SG (oh Grid -oh SG ) Where: D SG is the damping coefficient; D SG <<M SG Always applies to synchronous machines; start-up time M SG is related to the machine's inertia constant as shown in the formula: M SG =2H RI Where: H RI is the equivalent inertia constant; P m and p D,SG The influence of p m Split into three components: p m =p UC +p PFC +p SFC Where: p UC The power set point obtained by solving the unit commitment problem; p PFC and p SFC Active power regulated by primary frequency control PFC and secondary frequency control SFC respectively; Then And substitute p m =p UC +p PFC +p SFC And differentiate it, we get: in: is the rate of change of frequency, i.e. RoCoF, obtained by PMU measurement; the time derivative of RoCoF, i.e. Obtained through PMU measurement using a PI filter; In the formula, and is always true, i.e. in the time window that characterizes the inertial response of the generator; therefore, for an appropriate time window, we obtain: Based on the LMS least squares method, the inertia estimation formula needs to be modified to cope with the measurement error: For the i The measurement data obtained are: Formula H * is an expression used to determine the equivalent inertia constant because its fractional structure can cause serious numerical problems; when the numerator When the denominator is very small, but not 0, the effect will be significantly amplified due to the infinitesimal denominator, resulting in a peak, causing H * There are serious errors in the estimation of Rewritten as: For a set of measurements t1, t2, ..., t N-1 ,t N Define the following vector: Where: T represents the matrix transpose; A and b are collected and The matrix of data, time window [t1,t N ], the data size N of each matrix is ​​determined by the length of the time window T = t N -t1 and PMU sampling rate.

2. According to the method for identifying equivalent inertia constant of a power system based on data regression described in claim 1, virtual inertia is also applicable here, as shown in the following formula: P e,IBR =P UC +P PFC +P SFC +P VI -P D,IBR in: P VI is the power increased by the virtual inertia: The formula shows that H VI pass Since the virtual inertial support is the result of the control action, H can be adaptively adjusted by appropriately measuring the system state. VI , in order to better deal with frequency deviation problems after sudden unexpected events.

3. The method for identifying equivalent inertia constants of a power system based on data regression according to claim 1, characterized in that: Use the least squares method to find the equivalent inertia constant H, where H is the inertia constant of the synchronous machine: Minimize(AH * -b) 2 Transform the above formula: H * =(A T A) -1 A T b Be careful to avoid morbid A T A, it is convenient to adjust its diagonal elements with small diagonal elements. In the tests solved in preparation for this work, ε = 10 -6 To prevent A T A is singular.

Images