Method for Identifying Distribution of Grid Inertia Coefficient Based on Disturbance Measurement Data

Abstract

The present invention discloses a method for identifying the distribution of grid inertia coefficients based on disturbance measurement data. By analyzing the frequency time-domain data measured by PMUs at each node of the system, disturbance events are extracted, and the frequency response curves are fitted with variable orders to calculate the inertia coefficients corresponding to each node. Then, the above steps are repeated according to the next disturbance event, and the results obtained are weighted and averaged with the results of the previous step to make the identified values of the inertia coefficients at each node approach the true values. This method avoids the influence of the disturbance position and size on the results. At the same time, this method only requires the frequency data of the measurement nodes, has strong operability, can effectively identify the distribution of inertia coefficients of a power system with a high proportion of new energy access, and effectively improves the ability of system operators to evaluate the most vulnerable nodes in disturbance events, ensuring the safe operation of the power system.

Description

Technical Field

[0001] The present invention relates to the technical field of power system operation state perception, and particularly relates to a method for identifying the distribution of grid inertia coefficients based on disturbance measurement data. Background Art

[0002] Grid frequency stability is the basis for high-quality power supply. When there is unbalanced power in the power system, the synchronous machine rotor spontaneously injects or absorbs rotational inertia into the system to offset frequency fluctuations and provide frequency support for the system. In modern complex power grids, new energy power generation represented by wind power and photovoltaic power penetrates into the grid through power electronic devices, directly reducing the grid inertia reserve and resulting in a decrease in the frequency support ability of the grid during the moment of being disturbed.

[0003] The inertia coefficient is an index characterizing the ability of the grid to resist frequency changes and represents the rotational inertia that a generator of unit capacity can provide. Existing inertia research mostly focuses on the identification of the overall inertia coefficient at the system level. However, the internal inertia coefficient of the grid has spatial distribution characteristics, that is, the inertia coefficients of different nodes and regions may be different. Therefore, there may be significant differences in the frequency responses of different nodes in the disturbed system. The traditional central inertia coefficient and central frequency cannot accurately describe this difference, and nodes with a lower inertia coefficient level may fall below the minimum value allowed by the system. This feature is particularly obvious in power grids with weak grid connections and high power electronic penetration rates. Summary of the Invention

[0004] In order to solve the problem that the existing method for identifying the central inertia coefficient of the power system is difficult to effectively identify the ability of each node of the grid to resist frequency disturbances online, the present invention provides a method for identifying the distribution of grid inertia coefficients based on disturbance measurement data. The inertia coefficients of each node in the grid are identified through the node PMU measurement data, and the traditional overall inertia of the grid is further refined to the distribution of inertia at each node, quantifying the frequency stability degree of the node, so as to facilitate the grid dispatcher to accurately grasp the fine distribution of inertia caused by the change of the grid topology structure, and help the dispatching and operation personnel better evaluate the frequency vulnerable nodes of the grid and ensure the safe operation of the power system.

[0005] In order to achieve the above technical objectives, the technical solution of the present invention is as follows:

[0006] A method for identifying the distribution of grid inertia coefficients based on disturbance measurement data, comprising the following steps:

[0007] S1, obtaining the real-time frequency data sequence of each node of the power grid system.

[0008] S2. According to the acquired real-time frequency data sequence, calculate the real-time rate of change of frequency RoCoF for each node. If the RoCoF of a certain node exceeds the preset threshold, then define the moment when the RoCoF exceeds the preset threshold as the disturbance occurrence moment t0 of this node.

[0009] S3. Calculate the equivalent inertia coefficient without power source nodes: First, take the synchronous machine node with the largest capacity in the power grid system as the calibration node, and then calculate the equivalent inertia coefficient of the nodes without power source according to the following formula,

[0010]

[0011] where H i and f i are the equivalent inertia coefficient and node frequency of the i-th node without power source. H J and f J are the inertia coefficient and node frequency of the calibration node J. Where the power source refers to high-power power generation equipment directly connected to the main grid.

[0012] S4. Calculate the equivalent inertia coefficient with power source nodes: First, perform multi-order polynomial fitting on the frequency curve and unbalanced active power curve of the nodes with power source, and then calculate the equivalent inertia coefficient of the nodes with power source according to the following formula,

[0013]

[0014] where B 0-j is the constant term of the unbalanced power fitting polynomial at the moment t0, and A 1-j is the first-order constant term coefficient of the frequency fitting polynomial at the moment t0.

[0015] S5. Record the equivalent inertia coefficients of all nodes at the current disturbance moment as a set of system inertia coefficient values, and compare the current number of disturbances with the preset number K, and respectively execute the following processes according to the comparison results:

[0016] If it does not reach K, return to step S1 and execute in a loop.

[0017] If it is equal to K, aggregate the set of K system inertia coefficient values recorded as the identification result of the system inertia coefficient, and return to step S1 and execute in a loop.

[0018] If it exceeds K, starting from the current number of disturbances, count backwards K system inertia coefficient values and aggregate them all as the identification result of the system inertia coefficient, and return to step S1 and execute in a loop.

[0019] For the method described above, in step S1, after acquiring the frequency data sequence, it further includes a process of using a low-pass filter to filter out high-frequency noise.

[0020] In the method described above, in step S3, for the selection of the calibration node, the synchronous motor node with the largest capacity in the power grid system is selected. Or after removing all power electronic power supplies in the power grid system, the system center frequency f is calculated based on the frequencies, capacities, and inertia coefficients of the remaining generators, i.e., synchronous motors. coi , and then compare f coi with the Pearson correlation coefficient (R(f coi , f j ) of each synchronous motor frequency response curve, and use the synchronous motor node with the highest correlation coefficient as the calibration node, where the inertia coefficient of the calibration node is:

[0021] H J = H(Max[R(f coi , f j )])

[0022] where the power electronic power supply refers to the power supply connected to the grid through a power electronic converter.

[0023] In the method described above, in step S4, the multi-order polynomial fitting of the frequency curve and the unbalanced active power curve of the power supply node includes:

[0024] Perform an n-order polynomial fitting on the frequency f(t) and unbalanced active power ΔP(t) of the power supply node through the following formula

[0025]

[0026] where A0, A1,..., A n are the fitting polynomial coefficients of the frequency data, and B0, B1,..., B n are the fitting polynomial coefficients of the power deviation data. t is the time t0 when the disturbance occurs.

[0027] Then, through the fitting of the frequency curve and the unbalanced active power curve, the constant term B of the unbalanced power fitting polynomial at time t0 is obtained 0-j , and the first-order constant term coefficient A of the frequency fitting polynomial 1-j .

[0028] In the method described above, when performing polynomial fitting on the frequency time series, the determination process of n is as follows: Set the tolerance ε, start increasing the order from the fifth order, i.e., n = 5, until the absolute value of the difference between the calculation result and the previous order result is less than ε.

[0029] In the method described above, in step S5, if it is equal to K, then aggregate the set of K system inertia coefficient values recorded as the identification result of the system inertia coefficient, which includes the following steps:

[0030] Under the k-th perturbation, the set of system inertia coefficient values is: where P is the total number of nodes in the system, k is the number of current perturbations, and k = K.

[0031] Then, the identification result H of the system inertia coefficient after a total of K perturbation events Node is aggregated based on the following formula:

[0032]

[0033] For the method described above, in step S5, if it exceeds K, starting from the current number of perturbations, the sets of the previous K system inertia coefficient values are counted backwards and all aggregated as the identification result of the system inertia coefficient, which includes the following steps:

[0034] Under the (k + n)-th perturbation, the identification result H of the system inertia coefficient is updated through the following formula Node to obtain the updated identification result of the system inertia coefficient

[0035]

[0036] The technical effects of the present invention are as follows

[0037] 1. The present invention defines the node inertia coefficient. Based on the frequency time-domain data measured by the PMUs of each node in the system, curve fitting with variable orders is performed on the frequency response curve under perturbations, and the node inertia coefficient corresponding to each node is calculated. Compared with the traditional method for identifying the central inertia coefficient of the system, the node inertia coefficient distribution proposed by the present invention can give the frequency anti-perturbation ability at different spatial node positions of the system, providing powerful auxiliary decision-making information for power grid operators.

[0038] 2. By calculating a new node inertia index through new perturbation events and weighted averaging it with the old node inertia, the node inertia coefficient is updated accordingly. As the number of perturbation events increases, the node inertia coefficient of the power grid can be recursively updated based on new data, effectively eliminating the deviation of the single identification result.

[0039] 3. The present invention utilizes the frequency response time-series data under perturbation events such as power grid load switching measured by PMUs. Due to the frequent occurrence of such perturbations, it has stronger operability and engineering practicability compared with the perturbation-based method for identifying the inertia coefficient.

[0040] 4. The present invention only requires frequency data, avoiding the influence of specific parameters such as the magnitude and location of perturbations and the power grid network structure on the identification.

[0041] The following further illustrates the present invention with reference to the accompanying drawings. Description of the Drawings

[0042] Figure 1 This is the topology and generator parameters of the IEEE 3-machine 9-bus system in the embodiment of the present invention.

[0043] Figure 2 This is the schematic flow chart of the embodiment of the present invention.

[0044] Figure 3 This is the frequency response curve diagram under disturbance in the embodiment of the present invention.

[0045] Figure 4 This is the schematic diagram of the relationship between the fitting order and the error in the embodiment of the present invention.

[0046] Figure 5 This is the diagram of the identification result and the weighted average result of the node inertia coefficient in the embodiment of the present invention.

[0047] Figure 6 This is the frequency response and the central inertia frequency diagram under disturbance in the embodiment of the present invention. Detailed Embodiment

[0048] The method disclosed in this embodiment includes the following steps:

[0049] S1. Obtain the real-time frequency data sequence of each node in the power grid system. In this embodiment, the real-time frequency data sequence is obtained by using a PMU, that is, a phasor measurement unit. At the same time, in order to improve the accuracy, after obtaining the frequency data sequence, this embodiment also uses a low-pass filter to filter out the high-frequency noise of the frequency data sequence.

[0050] S2. According to the obtained real-time frequency data sequence, calculate the real-time rate of change of frequency RoCoF of each node. If the RoCoF of a certain node exceeds the preset threshold, the moment when the RoCoF exceeds the preset threshold is defined as the disturbance occurrence moment t0 of this node.

[0051] S3. Calculate the equivalent inertia coefficient of the nodes without power sources: First, take the synchronous motor node with the largest capacity in the power grid system as the calibration node, and then calculate the equivalent inertia coefficient of the nodes without power sources according to the following formula,

[0052]

[0053] where H i and f i are the equivalent inertia coefficient and the node frequency of the i-th node without power source. H J and f JTo calibrate the inertia coefficient and node frequency of node J, and at the same time, the power supply nodes mentioned here refer to nodes without large - power power supplies directly connected to the main grid, such as synchronous motors used in thermal power and hydropower, or energy storage, photovoltaic, permanent - magnet synchronous wind turbines, etc.

[0054] Meanwhile, in this embodiment, in addition to using the synchronous - motor node with the largest capacity in the power - grid system as the calibration node, a system model can also be built in the simulation software. Then, power - electronic power supplies in the power - grid system, such as those connected to the grid through power - electronic converters like energy storage, photovoltaic, and permanent - magnet synchronous wind turbines, are removed. Then, based on the frequencies, capacities, and inertia coefficients of the remaining generators, i.e., synchronous motors, the system center frequency f is calculated. coi , and then compare f coi with the Pearson correlation coefficient (R(f coi, f j ) of each synchronous - motor frequency - response curve, and use the node with the highest correlation coefficient as the calibration node, where the inertia coefficient of the calibration node is H J = H(Max[R(f coi , f j )])

[0055] S4. Calculate the equivalent inertia coefficient of the node with a power supply:

[0056] First, perform an n - order polynomial fitting on the frequency f(t) and unbalanced active power ΔP(t) of the node with a power supply through the following formula

[0057]

[0058] where A0, A1, …, A n are the fitting polynomial coefficients of the frequency data, and B0, B1, …, B n are the fitting polynomial coefficients of the power - deviation data. t is the time t0 when the disturbance occurs.

[0059] The determination process of n - order is as follows: Set the tolerance ε, start increasing the order from five - order until the absolute value of the difference between the calculation result and the result of the previous order is less than ε.

[0060] Then, through the fitting of the frequency curve and the unbalanced active - power curve, the constant term B 0-j of the unbalanced - power fitting polynomial at time t0 and the first - order constant - term coefficient A 1-j of the frequency fitting polynomial are obtained.

[0061] Then, the equivalent inertia coefficient H j ,

[0062]

[0063] Among them, B 0-j is the constant term of the unbalanced power fitting polynomial at time t0, and A 1-j is the first-order constant term coefficient of the frequency fitting polynomial at time t0.

[0064] S5. Record the equivalent inertia coefficients of all nodes at the current disturbance moment as a set of system inertia coefficient values, compare the current number of disturbances with the preset number K, and perform the following processes respectively according to the comparison results:

[0065] If it does not reach K, return to step S1 and execute in a loop.

[0066] If it is equal to K, aggregate the set of K system inertia coefficient values recorded as the identification result of the system inertia coefficient, and return to step S1 to execute in a loop. That is, at the k-th disturbance, the set of system inertia coefficient values is: where P is the total number of nodes in the system, k is the current number of disturbances, and k = K.

[0067] Then the identification result H of the system inertia coefficient after a total of K disturbance events Node is aggregated based on the following formula:

[0068]

[0069] If it exceeds K, starting from the current number of disturbances, count backwards K system inertia coefficient values and aggregate them all as the identification result of the system inertia coefficient, and return to step S1 to execute in a loop. That is, at the (k + n)-th disturbance, the identification result H of the system inertia coefficient is updated through the following formula Node to obtain the updated identification result of the system inertia coefficient

[0070]

[0071] For example, if the current is the (K + 1)-th new disturbance event, calculate the new node inertia coefficient and perform iterative update on H Node as follows:

[0072]

[0073] Thus, as the number of disturbance events increases, the grid node inertia coefficient can be updated according to new data, effectively eliminating the deviation of the single identification result.

[0074] The following further illustrates the embodiments of the present invention in combination with simulations:

[0075] The topological structure built in the PowerFactory simulation software in this embodiment is asFigure 1 The IEEE-3 machine 9-node system shown is used as the application scenario of the method. In the IEEE-3 machine 9-node system, 20 different load shedding disturbances are generated through Latin hypercube sampling as the data source of this embodiment. The load shedding amount each time is 2% - 7% of the load at that position. The generated data is used to verify the feasibility and effectiveness of a method for identifying the node inertia distribution based on PMU measurement data proposed by the present invention.

[0076] The method steps of this embodiment are as Figure 2 shown. Specifically, in this example:

[0077] First, after each disturbance, the frequency response time series of each node after the disturbance event is collected. The frequencies of each node in the system under a typical disturbance are as Figure 3 shown.

[0078] The node frequency data collected from each simulation is preprocessed and the rate of change of frequency RoCoF is calculated. In this example, the judgment criterion for the disturbance occurrence moment is that the absolute value of RoCoF of f COI at this moment is greater than 0.2 Hz / s. At t = 0.01 s, RoCoF fCOI = -0.224 Hz / s, and thus the power grid disturbance occurrence moment is determined.

[0079] According to the method described in step S2, the Pearson correlation coefficients between the frequencies of the three generator nodes and the system center and the system center frequency f COI are shown in the following table.

[0080]

[0081] The Pearson correlation coefficients between the frequencies f i of each node and f COI within 1000 ms from the disturbance occurrence moment are shown in the following table. By comparison, the frequency of generator node 2 is close to the system center frequency. Therefore, node 2 is selected as the calibration node.

[0082] Through the fitting of the frequency data, the system inertia coefficient is calculated. Taking the Figure 3 frequency data as an example, its corresponding fitting order error is as Figure 4 shown. The tolerance ε = 1% is set. When the order reaches the 10th order, the relative errors of each node are within the tolerance. Therefore, a 10th-order polynomial is used for fitting.

[0083] The same method is used to calculate for the remaining 19 disturbance events, and each event is iteratively updated according to the method described in step S5. Finally, the inertia coefficient distribution of each node in the system under 20 disturbance data is as Figure 5 shown. The inertia coefficients of each node are shown in the following table.

[0084]

[0085] To further illustrate the reliability of the identification results, a three-phase ground short-circuit disturbance event with a duration of 0.1 seconds is set between nodes 8 and 9 of the simulation system. The frequency response curves of each node after the disturbance are as Figure 6 shown. The identification results of the inertia coefficient magnitudes of each node are consistent with the degree of frequency fluctuation. Through analysis, it can be seen that the generator unit connected to node 3 has the highest single-machine inertia, the generator unit connected to node 1 has the lowest single-machine inertia, and the inertia of other nodes decreases in a gradient between the two. The identification results meet the theoretical expectations, verifying the practicability of the method.

[0086] In summary, the present invention further refines the parameter of the overall inertia of the traditional power grid system to the distribution of inertia at each node, enhancing the perception of power grid operators on the inertia parameters of the power grid and providing favorable reference information for their accurate judgment of weak nodes in the power grid and for dispatching operations.

[0087] The above embodiments are only used to illustrate the technical concept and features of the present invention, and their purpose is to enable those skilled in the art to understand the content of the present invention and implement it accordingly, and cannot be used to limit the protection scope of the present invention. Any equivalent changes or modifications made according to the spirit and essence of the present invention should be covered within the protection scope of the present invention.

Claims

1. A method for identifying the distribution of grid inertia coefficients based on perturbed measurement data, characterized in that It includes the following steps: S1. Obtain the real-time frequency data sequence of each node in the power grid system; S2. According to the obtained real-time frequency data sequence, calculate the real-time rate of change of frequency (RoCoF) of each node. If the RoCoF of a certain node exceeds the preset threshold, then define the moment when the RoCoF exceeds the preset threshold as the disturbance occurrence moment t0 of this node; S3. Calculate the equivalent inertia coefficient of the nodes without power sources: First, take the synchronous motor node with the largest capacity in the power grid system as the calibration node, and then calculate the equivalent inertia coefficient of the nodes without power sources according to the following formula, Among them, H i and f i are the equivalent inertia coefficient and node frequency of the i-th power-source-free node; H J and f J are the inertia coefficient and node frequency of the calibration node J; where the power source refers to a high-power power generation device directly connected to the main grid; S4. Calculate the equivalent inertia coefficient of the nodes with power sources: First, perform multi-order polynomial fitting on the frequency curve and unbalanced active power curve of the nodes with power sources, and then calculate the equivalent inertia coefficient of the nodes with power sources according to the following formula, Among which B 0-j is the constant term of the unbalanced power fitting polynomial at time t0, and A 1-j is the first-order constant term coefficient of the frequency fitting polynomial at time t0; S5. Record the equivalent inertia coefficients of all nodes at the current disturbance moment as a set of system inertia coefficient values, compare the current number of disturbances with the preset number K, and respectively execute the following processes according to the comparison results: If it does not reach K, return to step S1 and execute in a loop; If it is equal to K, aggregate the set of K system inertia coefficient values recorded as the identification result of the system inertia coefficient, and return to step S1 and execute in a loop; If it exceeds K, starting from the current number of disturbances, count backwards K system inertia coefficient values and aggregate them all as the identification result of the system inertia coefficient, and return to step S1 and execute in a loop.

2. The method according to claim 1, wherein In step S1, after obtaining the frequency data sequence, it further includes the process of using a low-pass filter to filter out high-frequency noise.

3. The method according to claim 1, characterized in that, In the said step S3, for the selection of the calibration node, the synchronous motor node with the largest capacity in the power grid system is selected; or after removing all the power electronic power supplies in the power grid system, the system center frequency f is calculated based on the frequencies, capacities and inertia coefficients of the remaining generators, i.e., synchronous motors. coi , and then f coi is compared with the Pearson correlation coefficients (R(f coi , f j ) of the frequency response curves of each synchronous motor, and the synchronous motor node with the highest correlation coefficient is used as the calibration node, where the inertia coefficient of the calibration node is: H J = H(Max[R(f coi , f j )]) Wherein the power electronic power source refers to the power source connected to the grid through a power electronic converter.

4. The method according to claim 1, wherein In step S4, the multi-order polynomial fitting of the frequency curve and unbalanced active power curve of the nodes with power sources includes: Perform n-order polynomial fitting on the frequency f(t) and unbalanced active power ΔP(t) of the nodes with power sources through the following formula wherein, A0, A1, …, A n are the fitting polynomial coefficients of the frequency data, and B0, B1, …, B n are the fitting polynomial coefficients of the power deviation data; t is the time t0 when the disturbance occurs; Then, by fitting the frequency curve and the unbalanced active power curve, the constant term B of the unbalanced power fitting polynomial at time t0 is obtained 0-j , and the first-order constant term coefficient A of the frequency fitting polynomial 1-j .

5. The method according to claim 4, characterized in that, When performing polynomial fitting on the frequency time series, the determination process of n is: Set the tolerance ε, start from the fifth order (i.e., n = 5) and increase the order until the absolute value of the difference between the calculation result and the result of the previous order is less than ε.

6. The method according to claim 1, characterized in that, In step S5, if it is equal to K, aggregating the set of K system inertia coefficient values recorded as the identification result of the system inertia coefficient includes the following steps: The set of system inertia coefficient values under the k-th disturbance is as follows: where P is the total number of nodes in the system, k is the number of the current disturbance, and k = K; The identification result H of the system inertia coefficient after a total of K disturbance events Node is aggregated based on the following formula:

7. The method according to claim 6, characterized in that, In step S5, if it exceeds K, starting from the current number of disturbances, counting backwards K system inertia coefficient values and aggregating them all as the identification result of the system inertia coefficient includes the following steps: Under the (k + n)-th perturbation, the identification result H of the system inertia coefficient is updated by the following formula to obtain the updated identification result of the system inertia coefficient Node ​

Images