A system area inertia identification method based on power grid area division

Abstract

The application discloses a system area inertia identification method based on power grid area division, comprising the following steps: monitoring bus frequency data of each key position of a power grid, decomposing each bus frequency oscillation component and extracting key components based on a Prony algorithm, and dividing the power grid into areas; monitoring and processing regional boundary tie lines, regional load power and regional bus frequency data, and equivalent dividing areas into a synchronous generator based on monitoring information; constructing a dynamic relationship model between electromagnetic power variation and frequency deviation of the equivalent generator, and deducing the relationship between the model and the inertia of the equivalent generator; identifying the dynamic model parameters by using a subspace method, and estimating the regional inertia based on the identified model; the method provided by the application has high accuracy, can accurately estimate system inertia, can further improve the accuracy of power grid regional inertia identification, and is beneficial to providing more reliable reference basis for power grid system frequency stability control.

Description

Technical Field

[0001] The technical field of this invention is the field of power system operation and control technology, and in particular, it relates to a system area inertia identification method based on power grid area division. Background Technology

[0002] Rotational inertia is a crucial factor affecting the frequency security and stability of power systems, playing a vital role in suppressing rapid frequency changes, especially in the initial stage after a disturbance. However, with the integration of large-scale renewable energy units, the (equivalent) inertia of the power system may vary with changes in operating scenarios and power plant control modes (some renewable energy plants have frequency regulation capabilities). Therefore, online tracking of system inertia is necessary to ensure the safe operation of the power system.

[0003] Existing inertia identification methods can be broadly categorized into two types based on whether or not disturbances are actively applied: inertia identification based on applied disturbances and inertia identification based on conventional power grid disturbance data. The former has the advantage of known disturbance excitation and high signal-to-noise ratio of the monitoring data; however, it requires disturbance application for each identification, potentially impacting power grid operation, and some large-disturbance methods may not be suitable for frequent online identification. Identification methods based on conventional disturbance data, on the other hand, do not require active disturbance application, better meeting the power grid's requirements for online identification; however, the disturbance signal cannot be directly obtained, and conventional power grid disturbance signals generally have small amplitudes and relatively low signal-to-noise ratios, thus requiring more refined identification methods. In terms of identification scope, methods can be divided into two forms: network-wide and regional identification. Power grid frequencies exhibit certain spatial distribution characteristics; bus frequencies at different locations not only share the same common fluctuation trend (COI frequency) but also have their own mutually oscillating fluctuation components. From a stability perspective, not only is the overall power grid inertia not too low, but the inertia distribution across different regions also needs to be reasonable. Therefore, it is necessary to use conventional disturbance data for online regional inertia identification. Compared with the active excitation method, the signal-to-noise ratio of monitoring data under conventional disturbance is relatively low, so maximizing the signal-to-noise ratio of the monitoring data will help improve the identification accuracy. In the process of regional inertia identification, the power monitoring data of boundary tie lines is an important factor affecting the inertia identification results. Different zoning methods may correspond to different sets of boundary tie lines, so dividing the region appropriately to maximize the power fluctuation amplitude of the regional tie lines used for identification will help improve the identification accuracy. Frequency distribution oscillation has a significant impact on the power fluctuation of regional tie lines, so grouping buses and generators with similar oscillation characteristics and close locations into the same region, and generally having relatively large tie line power fluctuations between two regions with different main oscillation components, is crucial.

[0004] Based on the above analysis, researching strategies for dividing the power grid region according to its operating status and further identifying regional inertia based on the division results is an urgent problem to be solved in the study of power system frequency stability control. Summary of the Invention

[0005] The purpose of this section is to outline some aspects of embodiments of the present invention and to briefly describe some preferred embodiments. Simplifications or omissions may be made in this section, as well as in the abstract and title of this application, to avoid obscuring the purpose of these documents; however, such simplifications or omissions should not be construed as limiting the scope of the invention.

[0006] In view of the above-mentioned problems, the present invention is proposed.

[0007] Therefore, the technical problem solved by this invention is to combine power grid area division with inertia identification, thereby improving the accuracy of area rotational inertia identification, which is of great significance for the subsequent frequency stability control of the power system.

[0008] To solve the above technical problems, the present invention provides the following technical solution: a system area inertia identification method based on power grid area division, comprising: monitoring the bus frequency data at each key location of the power grid, decomposing and extracting key components of the frequency oscillation of each bus based on the Prony algorithm, and dividing the power grid into areas by comparing the main oscillation components of the bus frequency at each location with the electrical distance; The boundary tie line of the divided area, the load power within the area, and the bus frequency data within the area are monitored and processed. Based on the monitoring information, the divided area is equivalent to a synchronous generator. A dynamic relationship model between the electromagnetic power change and frequency deviation of the equivalent generator is constructed, and the relationship between the model and the inertia of the equivalent generator is derived. The parameters of the dynamic model are identified using the subspace method, and the regional inertia is estimated based on the identified model.

[0009] As a preferred scheme for system area inertia identification based on power grid area division, the monitoring of bus frequency data at key locations includes: monitoring as many bus frequency data points as possible with significant location differences within the power grid using PMU equipment, with data sampling intervals of [missing information]. The total length of the data collection is .

[0010] As a preferred scheme for system area inertia identification based on power grid area division, the decomposition and key component extraction of each bus frequency oscillation component based on the Prony algorithm includes: First, a parameterized model for signal decomposition is established; Specifically, the frequency monitoring data of a bus is decomposed into a sum of multiple damped oscillation components at each sampling point: in, Indicates the discrete (sampling) time step; This represents the signal value at the nth moment, specifically the nth value in the bus frequency monitoring data; Indicates the order of the oscillatory components contained therein; and These represent the initial amplitude and phase of the k-th oscillation component, respectively. and This represents the attenuation coefficient and angular frequency; and These are variables introduced to describe the characteristics of each oscillation component for the convenience of subsequent analysis; This represents a noise signal; the purpose of subsequent Prony algorithm steps is to decompose the bus frequency monitoring signal into the above form, i.e., to extract parameters. , , and (k=1,2… ); A sample function matrix is ​​constructed using monitoring data, as follows: in, The total amount obtained from monitoring Bus system frequency data at each moment; Represents the sample function matrix; This represents the result of multiplying the original data, corresponding to the matrix. Elements in; The coefficients corresponding to each oscillation mode are determined using the least squares method, specifically including... ,in The solution process for the remaining coefficients is as follows: Based on the above sample function matrix, the equation can be written as follows: in, In this formula It is known that it contains A linear equation and For each unknown, the coefficients corresponding to the oscillation modes are determined using the least squares method. Solve the problem and record the result as follows: ; The attenuation coefficient and oscillation frequency parameters of each oscillation mode are calculated based on the coefficients corresponding to each oscillation mode, as follows: Based on the coefficients corresponding to the above oscillation modes and its decaying oscillation characteristic parameters for each oscillation mode The relationship between them: Can be The parameters are calculated, and the result is recorded as follows: ; Based on oscillation characteristic parameters The initial amplitude and phase parameters of each signal, including each oscillation mode, are calculated as follows: Based on the recursive difference equation and The estimated values ​​of the original data sequence are calculated as follows: Given Continuous recursion ; Based on the relationship between the original data sequence and the oscillation characteristic parameters constructed above: as well as Can be The estimation is performed, and the result is recorded as follows: ; According to parameters and The calculation results for the parameters , , and (k=1,2… The calculation is performed as follows: Here, Im and Re represent taking the real part and imaginary part, respectively.

[0011] As a preferred embodiment of a system region inertia identification method based on power grid region division, wherein: the region division of the power grid includes: For a single bus frequency oscillation component and Parameters (representing the initial amplitude of the oscillation component), extract them. Larger Group (pre-set parameters) oscillation components; The results obtained by pairwise comparison of each busbar Group oscillation components, if In the group of oscillating components, there exists (Preset parameters) oscillation components above group and If the parameters are similar, then the two bus lines are grouped into the same group; where the determination is... and The principle for ensuring that the parameters are similar is to calculate the relative deviation between the two parameters mentioned above, and if it is less than the set deviation threshold... If so, then the damped oscillation characteristics of the two sets of oscillation components are considered to be similar; Calculate the electrical distance between every two buses within each group. If the electrical distance between one bus and two or more other buses in a group exceeds a set threshold, the group is considered closed. When this happens, the bus is removed from the group, and its frequency oscillation component is compared with the frequency oscillation components of buses in other groups. The preset parameters are then appropriately reduced. Find a group with similar frequency oscillation characteristics and add the bus to that group to complete the partitioning; if no group can be found that meets the above parameters... For groupings with similar group requirements, the busbar and its nearby unmonitored busbars, generators, loads, etc., are divided into a separate area.

[0012] As a preferred embodiment of a system area inertia identification method based on power grid area division, the method involves: monitoring the power fluctuation data of tie lines between areas, the frequency data of key buses within each area, and the key load data based on the area division results, and calculating the frequency and terminal electromagnetic power information of the equivalent motor, as detailed below: Calculate the equivalent electromagnetic power of the motor based on the power of the inter-regional tie line and the load within the region: in, This indicates the total number of connecting lines between this area and other areas; This indicates the total load in the area; Indicates the power of the tie line; This indicates the load power. In actual power grids, since it is difficult to monitor all loads in a region, it is possible to selectively monitor major loads and loads with large fluctuations in the region.

[0013] As a preferred embodiment of a system area inertia identification method based on power grid area division, the area COI frequency is calculated based on the key bus frequency data within the divided area. in, This indicates the frequency deviation of the monitored bus. Indicates the number of buses with frequency monitoring capabilities; This represents the weighting coefficient, which can be assigned values ​​based on the relative fluctuation range of the frequency and the quality of the monitoring data, such as... var indicates that the variance is calculated; The change in mechanical power of the equivalent generator is analogous to the structure of a synchronous generator governor, which is a primary frequency regulation stage based on frequency feedback.

[0014] As a preferred scheme for system area inertia identification based on power grid area division, the construction of the dynamic relationship model includes: adding a primary frequency regulation stage to the equivalent generator and converting it into a transfer function form, as follows: in, Represents the Laplace operator; This represents the unit power regulation coefficient (takes a negative value); This represents the transfer function that describes the dynamics of primary frequency modulation power regulation; This represents the transfer function.

[0015] As a preferred embodiment of a system region inertia identification method based on power grid region division, the relationship between the derived model and the equivalent generator inertia includes: According to the final value theorem, calculate the initial slope of the unit step response curve corresponding to the model: in, The unit step response representing the derivative of the equivalent motor frequency deviation; This represents the unit step response of the equivalent motor frequency deviation; it is worth noting that the Laplace transform corresponding to the differential element is s. Numerically satisfies (This indicates the generator's ability to track power commands.) Based on the above calculations, the functional relationship between the initial slope of the unit step response curve of the model and the regional inertia can be obtained. According to the final value theorem, calculate the steady-state value of the unit step response curve corresponding to the model: Based on the above calculations, the functional relationship between the model's unit step response curve and the regional primary frequency modulation coefficient and damping coefficient can be obtained.

[0016] As a preferred embodiment of a system region inertia identification method based on power grid region division, the identification of dynamic model parameters using the subspace method includes: The standard discrete form structure of the parameterized model to be identified is as follows: in, , , and This represents the coefficient matrix corresponding to the discrete system. This indicates the error caused by noise; This represents the coefficient matrix corresponding to the error term; the purpose of subspace identification is to... and Identify the corresponding coefficient matrix based on measurement results over a period of time; During the process of identifying the frequency deviation and electromagnetic power change model using the subspace algorithm, the calculated frequency deviation of the equivalent generator in the region is used as the output of the model to be identified, denoted as y; and the electromagnetic power change data of the equivalent generator in the region is used as the input of the model to be identified, denoted as u.

[0017] The subspace algorithm is used to construct the Hankel matrices corresponding to the forward and backward historical data, as follows: in, and Corresponding to the input respectively The Hankel matrix composed of data from different time periods; and Corresponding to the input respectively The Hankel matrix composed of data from different time periods; and Indicates the number of rows and columns of a matrix; During the execution of the subspace algorithm, the numerical relationship between the Hankel matrix and the parameter matrix of the model to be identified is established, as follows: in, and It is similar in form and The error term Hankel matrix; , , All are formed by transforming and combining the model coefficient matrix; , The forms are respectively with , Similar (the difference lies in the coefficients) and ), which will not be described in detail here; During the execution of the subspace identification algorithm, orthogonal projection is used to construct a correlation matrix to reduce the impact of noise, specifically as follows: based on the noise With input data The weak correlation between them can be reduced by the following transformation to decrease the noise effect: in, ; Representation and matrix Orthogonal complement; Indicates about The orthogonal projection operator of a matrix can be calculated using RQ decomposition; During the execution of the subspace identification algorithm, the matrix constructed above is... The transformation and singular value decomposition are performed as follows: A weighting coefficient matrix is ​​introduced using the N4SID method. and Further construct the matrix Then to Perform singular value decomposition, and based on the resulting singular value matrix, approximate it into a block matrix form: During the execution of the subspace identification algorithm, the model parameter matrix is ​​estimated based on the above singular value decomposition results. The specific process is as follows: Estimate the correlation matrix of the model parameters defined above based on the singular value decomposition results. , recorded as ,as follows: according to For model parameters Calculations are performed: Since the actual system is a single-input single-output system, therefore The first line in ,Pick The first line is the identification result. ; according to For model parameters Perform calculations: remove The first line is recorded as Remove the last line and denote it as ,according to ,according to Using least squares estimation, we can obtain ; according to For model parameters and Calculation: First, transform the numerical relationships between the above matrices to obtain: in, express The orthogonal complement satisfies ; for The generalized inverse matrix; only The quantity to be determined. In the above equation... middle and The elements are linearly related to each other, according to and The estimation results were finally estimated using the least squares method. and ; The final system model obtained by executing the subspace identification algorithm is as follows: in: Represents state variables; This represents the input variable, specifically the change in the terminal electromagnetic power of the equivalent generator in the corresponding region. This represents the output variable, specifically the frequency deviation of the equivalent generator in the corresponding region. The unit step response of the system model obtained from the above identification is calculated as follows: At the initial moment, i.e., when k=0, let ;and For the output variable at k=0 and state variables Perform calculations; At time k=1, let , combined calculate and ; In subsequent moments Always keep it at 1, iteratively calculate at each time step. and ; Finally obtained The values ​​taken at each time point are the unit step response curves of the model.

[0018] As a preferred embodiment of a system region inertia identification method based on power grid region division, wherein: the estimation of region inertia based on the identified model includes: Let the steady-state value of the unit step response curve of the obtained model be... Based on the functional relationship between the model's step response and the equivalent motor's primary frequency regulation coefficient and damping coefficient, the following can be derived: ; Let the unit step response curve of the obtained model be... The mean slope within the phase is Based on the functional relationship between the model's step response and the equivalent motor inertia, the following can be derived: .

[0019] The beneficial effects of this invention are as follows: The system region inertia identification method based on power grid region division provided by this invention first uses the Prony algorithm to extract the oscillation components of the bus frequency in the power grid, and then, based on the power grid region division, it has high accuracy; it proposes a method for partitioned inertia identification based on conventional random disturbance data of the power grid, and introduces a subspace method with high accuracy and robustness and low algorithm complexity to identify the system state space model, thereby enabling more accurate estimation of system inertia; it proposes a power grid region inertia identification method that considers the power grid region division results, which can further improve the accuracy of power grid region inertia identification and is beneficial to providing a more reliable reference for the frequency stability control of the power grid system. Attached Figure Description

[0020] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. Wherein: Figure 1 This is an overall flowchart of a system area inertia identification method based on power grid area division according to the first embodiment of the present invention; Figure 2 This is a specific implementation flow of the Prony algorithm in a system region inertia identification method based on power grid region division as described in the first embodiment of the present invention; Figure 3 This is a diagram illustrating the regional inertia equivalence method in a system regional inertia identification method based on power grid regional division as described in the first embodiment of the present invention. Figure 4 This is a specific implementation flow of the subspace identification algorithm in the system region inertia identification method based on power grid region division described in the first embodiment of the present invention; Figure 5 : This refers to the topology of the detailed power grid model used in the simulation example of a system region inertia identification method based on power grid region division as described in the second embodiment of the present invention. Detailed Implementation

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

[0022] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0023] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in different places in this specification does not necessarily refer to the same embodiment, nor is it a single or selective embodiment that is mutually exclusive with other embodiments.

[0024] This invention is described in detail with reference to the schematic diagrams. When detailing the embodiments of this invention, for ease of explanation, the cross-sectional views illustrating the device structure may be partially enlarged, not adhering to the usual scale. Furthermore, the schematic diagrams are merely examples and should not be construed as limiting the scope of protection of this invention. In actual fabrication, the three-dimensional spatial dimensions of length, width, and depth should be included.

[0025] Furthermore, in the description of this invention, it should be noted that the terms "upper," "lower," "inner," and "outer," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. These terms are used solely for the convenience of describing the invention and for simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on the invention. In addition, the terms "first," "second," or "third" are used for descriptive purposes only and should not be construed as indicating or implying relative importance.

[0026] Unless otherwise explicitly specified and limited, the terms "installation," "connection," and "joining" in this invention should be interpreted broadly. For example, they can refer to fixed connections, detachable connections, or integral connections; similarly, they can refer to mechanical connections, electrical connections, or direct connections, or indirect connections through an intermediate medium, or internal connections between two components. Those skilled in the art can understand the specific meaning of the above terms in this invention based on the specific circumstances.

[0027] Example 1 Reference Figure 1-4 This is the first embodiment of the present invention, which provides a system region inertia identification method based on power grid region division, including: S1: Monitor the frequency data of the bus at each key location in the power grid, decompose the oscillation components of each bus frequency and extract the key components based on the Prony algorithm, and divide the power grid into regions by comparing the main oscillation components of the bus frequency at each location with the electrical distance. Furthermore, the monitoring of bus frequency data at key locations includes: monitoring as many bus frequency data points as possible within the power grid as possible, using PMU equipment, with data sampling intervals of [missing information]. The total length of the data collection is .

[0028] Specifically, the decomposition and key component extraction of the frequency oscillation components of each bus based on the Prony algorithm includes: First, a parameterized model for signal decomposition is established; Specifically, the frequency monitoring data of a bus is decomposed into a sum of multiple damped oscillation components at each sampling point: in, Indicates the discrete (sampling) time step; This represents the signal value at the nth moment, specifically the nth value in the bus frequency monitoring data; Indicates the order of the oscillatory components contained therein; and These represent the initial amplitude and phase of the k-th oscillation component, respectively. and This represents the attenuation coefficient and angular frequency; and These are variables introduced to describe the characteristics of each oscillation component for the convenience of subsequent analysis; This represents a noise signal; the purpose of subsequent Prony algorithm steps is to decompose the bus frequency monitoring signal into the above form, i.e., to extract parameters. , , and (k=1,2… ); A sample function matrix is ​​constructed using monitoring data, as follows: in, The total amount obtained from monitoring Bus system frequency data at each moment; Represents the sample function matrix; This represents the result of multiplying the original data, corresponding to the matrix. Elements in; The coefficients corresponding to each oscillation mode are determined using the least squares method, specifically including... ,in The solution process for the remaining coefficients is as follows: Based on the above sample function matrix, the equation can be written as follows: in In this formula It is known that it contains A linear equation and For each unknown, the coefficients corresponding to the oscillation modes are determined using the least squares method. Solve the problem and record the result as follows: ; The attenuation coefficient and oscillation frequency parameters of each oscillation mode are calculated based on the coefficients corresponding to each oscillation mode, as follows: Based on the coefficients corresponding to the above oscillation modes and its decaying oscillation characteristic parameters for each oscillation mode The relationship between them: Can be The parameters are calculated, and the result is recorded as follows: ; Based on oscillation characteristic parameters The initial amplitude and phase parameters of each signal, including each oscillation mode, are calculated as follows: Based on the recursive difference equation and The estimated values ​​of the original data sequence are calculated as follows: Given Continuous recursion ; Based on the relationship between the original data sequence and the oscillation characteristic parameters constructed above: as well as Can be The estimation is performed, and the result is recorded as follows: ; According to parameters and The calculation results for the parameters , , and (k=1,2… The calculation is performed as follows: Here, Im and Re represent taking the real part and imaginary part, respectively.

[0029] Furthermore, the regional division of the power grid includes: For a single bus frequency oscillation component and Parameters (representing the initial amplitude of the oscillation component), extract them. Larger Group (pre-set parameters) oscillation components; The results obtained by pairwise comparison of each busbar Group oscillation components, if In the group of oscillating components, there exists (Preset parameters) oscillation components above group and If the parameters are similar, then the two bus lines are grouped into the same group; where the determination is... and The principle for ensuring that the parameters are similar is to calculate the relative deviation between the two parameters mentioned above, and if it is less than the set deviation threshold... If so, then the damped oscillation characteristics of the two sets of oscillation components are considered to be similar; Calculate the electrical distance between every two buses within each group. If the electrical distance between one bus and two or more other buses in a group exceeds a set threshold, the group is considered closed. When this happens, the bus is removed from the group, and its frequency oscillation component is compared with the frequency oscillation components of buses in other groups. The preset parameters are then appropriately reduced. Find a group with similar frequency oscillation characteristics and add the bus to that group to complete the partitioning; if no group can be found that meets the above parameters... For groupings with similar group requirements, the busbar and its nearby unmonitored busbars, generators, loads, etc., are divided into a separate area.

[0030] It should be noted that this step involves introducing the Prony algorithm to extract the oscillation components of the bus frequencies at various locations in the power grid, and further completing the grid area division. A reasonable grid zoning result will ensure that the power fluctuation amplitude of the interconnecting lines between the monitored areas is relatively large, thereby improving the signal-to-noise ratio of the measurement data used for subsequent area inertia identification, and thus improving the accuracy of the final area inertia identification.

[0031] S2: Monitor and process the boundary tie line of the divided area, the load power within the area, and the bus frequency data within the area, and based on the monitoring information, equate the divided area to a synchronous generator; Specifically, based on the regional division results, the power fluctuation data of the tie lines between each region, the frequency data of the key busbars within each region, and the key load data are monitored, and the frequency and terminal electromagnetic power information of the equivalent motor are calculated, as follows: Calculate the equivalent electromagnetic power of the motor based on the power of the inter-regional tie line and the load within the region: in, This indicates the total number of connecting lines between this area and other areas; This indicates the total load in the area; Indicates the power of the tie line; This indicates the load power. In actual power grids, since it is difficult to monitor all loads in a region, it is possible to selectively monitor major loads and loads with large fluctuations in the region.

[0032] Furthermore, the COI frequency of the region is calculated based on the key bus frequency data within the defined region: in, This indicates the frequency deviation of the monitored bus. Indicates the number of buses with frequency monitoring capabilities; This represents the weighting coefficient, which can be assigned values ​​based on the relative fluctuation range of the frequency and the quality of the monitoring data, such as... var indicates that the variance is calculated; The change in mechanical power of the equivalent generator is analogous to the structure of a synchronous generator governor, which is a primary frequency regulation stage based on frequency feedback.

[0033] It should be noted that this step analyzes the monitoring data and treats the regional power grid as a synchronous generator unit, using the overall inertia of the entire regional power grid as the inertia of the equivalent generator. This facilitates the subsequent derivation of the relationship between the regional inertia and the equivalent generator model of the region, and simplifies the computational complexity of the regional inertia.

[0034] S3: Construct a dynamic relationship model between the electromagnetic power change and frequency deviation of the equivalent generator, and derive the relationship between the model and the inertia of the equivalent generator; Specifically, the construction of the dynamic relationship model includes: adding a frequency regulation stage to the equivalent generator and converting it into a transfer function form, as follows: in, Represents the Laplace operator; This represents the unit power regulation coefficient (takes a negative value); This represents the transfer function that describes the dynamics of primary frequency modulation power regulation; This represents the transfer function.

[0035] Furthermore, the relationship between the derived model and the equivalent generator inertia includes: According to the final value theorem, calculate the initial slope of the unit step response curve corresponding to the model: in, The unit step response representing the derivative of the equivalent motor frequency deviation; This represents the unit step response of the equivalent motor frequency deviation; it is worth noting that the Laplace transform corresponding to the differential element is s. Numerically satisfies (This indicates the generator's ability to track power commands.) Based on the above calculations, the functional relationship between the initial slope of the unit step response curve of the model and the regional inertia can be obtained. According to the final value theorem, calculate the steady-state value of the unit step response curve corresponding to the model: Based on the above calculations, the functional relationship between the model's unit step response curve and the regional primary frequency modulation coefficient and damping coefficient can be obtained.

[0036] It should be noted that this step elucidates the relationship between the overall inertia of the regional power grid and the unit step response of the equivalent model of the regional power grid, greatly simplifying the calculation process of the system's regional inertia.

[0037] S4: Identify the parameters of the dynamic model using the subspace method, and estimate the regional inertia based on the identified model.

[0038] Specifically, the identification of dynamic model parameters using the subspace method includes: The standard discrete form structure of the parameterized model to be identified is as follows: in, , , and This represents the coefficient matrix corresponding to the discrete system. This indicates the error caused by noise; This represents the coefficient matrix corresponding to the error term; the purpose of subspace identification is to... and Identify the corresponding coefficient matrix based on measurement results over a period of time; During the process of identifying the frequency deviation and electromagnetic power change model using the subspace algorithm, the calculated frequency deviation of the equivalent generator in the region is used as the output of the model to be identified, denoted as y; and the electromagnetic power change data of the equivalent generator in the region is used as the input of the model to be identified, denoted as u.

[0039] The subspace algorithm is used to construct the Hankel matrices corresponding to the forward and backward historical data, as follows: in, and Corresponding to the input respectively The Hankel matrix composed of data from different time periods; and Corresponding to the input respectively The Hankel matrix composed of data from different time periods; and Indicates the number of rows and columns of a matrix; During the execution of the subspace algorithm, the numerical relationship between the Hankel matrix and the parameter matrix of the model to be identified is established, as follows: in, and It is similar in form and The error term Hankel matrix; , , All are formed by transforming and combining the model coefficient matrix; , The forms are respectively with , Similar (the difference lies in the coefficients) and ), which will not be described in detail here; During the execution of the subspace identification algorithm, orthogonal projection is used to construct a correlation matrix to reduce the impact of noise, specifically as follows: based on the noise With input data The weak correlation between them can be reduced by the following transformation to decrease the noise effect: in, ; Representation and matrix Orthogonal complement; Indicates about The orthogonal projection operator of a matrix can be calculated using RQ decomposition; During the execution of the subspace identification algorithm, the matrix constructed above is... The transformation and singular value decomposition are performed as follows: A weighting coefficient matrix is ​​introduced using the N4SID method. and Further construct the matrix Then to Perform singular value decomposition, and based on the resulting singular value matrix, approximate it into a block matrix form: During the execution of the subspace identification algorithm, the model parameter matrix is ​​estimated based on the above singular value decomposition results. The specific process is as follows: Estimate the correlation matrix of the model parameters defined above based on the singular value decomposition results. , recorded as ,as follows: according to For model parameters Calculations are performed: Since the actual system is a single-input single-output system, therefore The first line in ,Pick The first line is the identification result. ; according to For model parameters Perform calculations: remove The first line is recorded as Remove the last line and denote it as ,according to ,according to Using least squares estimation, we can obtain ; according to For model parameters and Calculation: First, transform the numerical relationships between the above matrices to obtain: in, express The orthogonal complement satisfies ; for The generalized inverse matrix; only The quantity to be determined. In the above equation... middle and The elements are linearly related to each other, according to and The estimation results were finally estimated using the least squares method. and ; The final system model obtained by executing the subspace identification algorithm is as follows: in: Represents state variables; This represents the input variable, specifically the change in the terminal electromagnetic power of the equivalent generator in the corresponding region. This represents the output variable, specifically the frequency deviation of the equivalent generator in the corresponding region. The unit step response of the system model obtained from the above identification is calculated as follows: At the initial moment, i.e., when k=0, let ;and For the output variable at k=0 and state variables Perform calculations; At time k=1, let , combined calculate and ; In subsequent moments Always keep it at 1, iteratively calculate at each time step. and ; Finally obtained The values ​​taken at each time point are the unit step response curves of the model.

[0040] Furthermore, the estimation of regional inertia based on the identified model includes: Let the steady-state value of the unit step response curve of the obtained model be... Based on the functional relationship between the model's step response and the equivalent motor's primary frequency regulation coefficient and damping coefficient, the following can be derived: ; Let the unit step response curve of the obtained model be... The mean slope within the phase is Based on the functional relationship between the model's step response and the equivalent motor inertia, the following can be derived: .

[0041] Example 2 Reference Figure 5 As an embodiment of the present invention, a system region inertia identification method based on power grid region division is provided. In order to verify the beneficial effects of the present invention, a simulation experiment is conducted for scientific demonstration.

[0042] The technical solution of the system of the present invention is a regional inertia identification system, including: an AC power grid, a PMU data monitoring module, a power grid partitioning module, and a regional inertia identification module; the PMU data monitoring module monitors the bus frequency and transmission line power data of the AC power grid under random disturbances of normal loads, transmits the data to the power grid partitioning module to realize the division of the power grid into regions, and then transmits the region partitioning results and PMU monitoring data to the regional inertia identification module to output the regional inertia identification results of the power grid.

[0043] This simulation example is specifically illustrated in the appendix. Figure 5 The overall structure of the aforementioned regional inertia identification system will be explained using an example. This is achieved through monitoring... Figure 5 Ten frequency data points (39, 2, 25, 28, 21, 23, 20, 19, 12, and 6) were used for the middle bus section. Based on the Prony algorithm, the main oscillation components of each bus frequency were extracted. By comparing the main oscillation components of each bus frequency, buses and related units with similar oscillation components and locations were grouped into the same region (final division results are shown in [link to final division]). Figure 5 The boundary lines in the region); based on the above regional division results, by monitoring the frequency data of key busbars within the region under random load disturbances of the power grid (such as... Figure 5 Busbars 2 and 25 of area 2), power output from the area boundary connecting line (e.g.) Figure 5 Boundary connection lines of Area 2 (26-28, 26-29, 17-16, 18-3, 2-3, 2-1) and key loads within the area (such as...) Figure 5The 25 loads on the busbar in Region 2 are considered as a single synchronous generator unit. Monitoring data is used to calculate the change in electromagnetic power at the generator terminals and the frequency deviation of this equivalent generator. Based on the structure of the equivalent generator in the regional power grid, a model reflecting the dynamic relationship between the frequency deviation and the change in electromagnetic power at the generator terminals is constructed. A method for calculating the theoretical value of the inertia of this equivalent generator is proposed based on the unit step response of this model. The calculated change in electromagnetic power and frequency deviation of the equivalent generator are used as the input and output data of the model to be identified. The parameters of the model are estimated using a subspace identification method. Finally, the unit step response of the obtained model (such as the model reflecting the dynamic relationship between the frequency and electromagnetic power at the generator terminals of the equivalent generator in Region 2) is obtained, thereby estimating the inertia, primary frequency regulation coefficient, and damping coefficient of the region.

[0044] The specific steps include: Step 1: By monitoring the bus frequency data at various key locations in the power grid, the Prony algorithm is used to construct a sample function matrix. The least squares method is used to determine the coefficients corresponding to each oscillation mode. Based on these coefficients, the attenuation coefficient and oscillation frequency parameters of each oscillation mode are calculated. Then, the initial amplitude and phase of each oscillation mode contained in each signal are calculated. This process completes the extraction of the system frequency oscillation components. The main oscillation components of the bus frequency at each location and the electrical location are compared to divide the power grid area. Preferably, the monitoring of bus frequency data at key locations in step 1 specifically involves: monitoring as many bus frequency data points as possible within the power grid as possible, using a PMU device, with data sampling intervals of [missing information]. The total length of the data collection is In this example power grid, the monitored busbars include: busbars 39, 2, 25, 28, 21, 23, 20, 19, 12, and 6. In step 1, the process of extracting the system frequency oscillation components using the Prony algorithm first establishes a parameterized model for signal decomposition, as follows: Taking the system frequency curve measured at bus 39 as an example, it is decomposed into a sum of multiple damped oscillation components at each sampling point: Wherein, the order of the oscillation component Sampling time step ; Indicates the 39th position of the busbar n Frequency monitoring data at each moment; and They represent the first kThe initial amplitude and phase of each oscillation component; and This represents the attenuation coefficient and angular frequency; and These are variables introduced to describe the characteristics of each oscillation component for the convenience of subsequent analysis; This represents a noise signal. The purpose of subsequent steps in the Prony algorithm is to decompose the bus frequency monitoring signal into the above form, i.e., to extract the parameters. , , and (k=1,2…18).

[0045] In step 1, the process of extracting the system frequency oscillation component using the Prony algorithm first requires constructing a sample function matrix using monitoring data, as follows: in, The total amount obtained from monitoring Bus system frequency data at each moment; Represents the sample function matrix; This represents the result of multiplying the original data, corresponding to the matrix. The elements in.

[0046] Step 1 involves using the least squares method to determine the coefficients corresponding to each oscillation mode during the execution of the Prony algorithm. Specifically, this includes... ,in The solution process for the remaining coefficients is as follows: Based on the above sample function matrix, the equation can be written as follows: in In this formula Given that there are 18 linear equations and 18 unknowns, the coefficients corresponding to each oscillation mode are determined using the least squares method. Solve the problem and record the result as follows: Step 1 involves calculating the attenuation coefficient and oscillation frequency parameters of each oscillation mode based on the coefficients corresponding to each oscillation mode during the execution of the Prony algorithm. Specifically, this is done by calculating the coefficients corresponding to each oscillation mode as described above. and its decaying oscillation characteristic parameters for each oscillation mode Relationship Can be The parameters are calculated, and the result is recorded as follows: Step 1 involves using the oscillation characteristic parameters during the execution of the Prony algorithm. The initial amplitude and phase parameters of each signal, including each oscillation mode, are calculated as follows: (1) Based on the recursive difference equation and The estimated values ​​of the original data sequence are calculated as follows: Given Continuous recursion .

[0047] (2) Based on the relationship between the original data sequence constructed above and the oscillation characteristic parameters as well as Can be The estimation is performed, and the result is recorded as follows: Step 1: During the execution of the Prony algorithm, based on the parameters... and The calculation results for the parameters , , and (k=1,2… The calculation is performed as follows: Here, Im and Re represent taking the real part and imaginary part, respectively.

[0048] The process of dividing the power grid into regions based on the extraction results of the main oscillation components of the signal and the electrical location in step 1 is as follows: (1) For a single bus frequency oscillation component and Parameters (representing the initial amplitude of the oscillation component), extract them. Larger Group (pre-set parameters) oscillation components, in this example =5; (2) Compare each busbar pairwise. Based on the 5 sets of oscillation components extracted in (1), if there is a pairwise comparison of the 5 sets of oscillation components, then the pairwise comparison is performed. (In this example) Oscillation components with values ​​of 3 or more and If the parameters are similar, then the two bus lines are grouped into the same group; where the determination is... and The principle for ensuring that the parameters are similar is to calculate the relative deviation between the two parameters mentioned above, and if it is less than the set deviation threshold... If so, it is considered that the damped oscillation characteristics of the two sets of oscillation components are similar.

[0049] (3) Calculate the electrical distance between every two buses in each group obtained in (2). When there is an electrical distance between one bus and two or more other buses in a certain group that is greater than the set threshold, When this happens, the busbar is removed from the group, and its frequency oscillation component is compared with the frequency oscillation components of other buses in the group according to the process in (2), and the preset parameters are appropriately reduced. Find a group with similar frequency oscillation characteristics and add the bus to that group to complete the partitioning; if no group can be found that meets the above parameters... For groupings with similar group requirements, the busbar and its nearby unmonitored busbars, generators, loads, etc., are divided into a separate area.

[0050] In this example, the final partitioning result is (see...) Figure 5 (The boundary lines in the diagram are described using generators as individuals): Region 1 includes generator G1; Region 2 includes generators G8 and G10; Region 3 includes generators G2 and G3; Region 4 includes generators G4, G5, G6, G7, and G9.

[0051] Step 1 above introduces the Prony algorithm to extract the oscillation components of the bus frequency at each location in the power grid and further completes the division of the power grid area. The reasonable power grid zoning result will ensure that the power fluctuation amplitude of the interconnection lines between the monitored areas is relatively large, thereby improving the signal-to-noise ratio of the measurement data used for subsequent area inertia identification, and thus improving the accuracy of the final area inertia identification.

[0052] Step 2: Based on the rotor motion equation of the synchronous generator, this region is equivalent to a parameterized synchronous generator with a fixed structure in terms of frequency stability analysis. Based on the regional division results, the frequency and terminal electromagnetic power information of the equivalent generator are calculated by monitoring the power fluctuation data of the tie lines between regions, the frequency data of the key bus in the region, and the fluctuation data of the key load. Preferably, in step 2, regarding the frequency stability analysis problem, the region is equivalent to a synchronous generator unit based on the rotor motion equation of the synchronous generator. The structure of this equivalent generator unit is as follows: in, This indicates the deviation of the angular velocity from the rated value; , These represent the changes in mechanical power and electromagnetic power, respectively. Indicates the damping coefficient; It represents the inertial constant.

[0053] Step 2 describes processing the power fluctuation data of the interconnection lines between regions based on the regional division results (specifically including each boundary interconnection line, such as...). Figure 5 Lines 2-3, 1-2, 8-9, 4-5, 4-14, 18-3, 25-26, 17-16, 16-15, 16-21, 16-24, 28-26, and 29-26 in the region, and key bus frequency data within the region (such as...). Figure 5 Busbars 39, 2, 25, 28, 21, 23, 20, 19, 12, and 6) and critical loads (such as...) Figure 5 The load data at busbars 3, 4, 8, 15, and 20 are monitored, and the frequency and terminal electromagnetic power information of the equivalent motor are calculated to... Figure 5 The entire central area 2 is equivalent to one generator set, as detailed below: (1) The equivalent electromagnetic power of the motor is calculated based on the power of the inter-regional tie lines and the load within the region. The inter-regional tie lines include: lines 26-28, 26-29, 17-16, 18-3, 2-3, and 2-1, totaling 6 lines. The load power fluctuation at bus 25 of the key load data in this region is statistically analyzed. Therefore, the equivalent electromagnetic power of the motor in this region is as follows: In actual power grids, since it is difficult to monitor all loads in a region, it is possible to selectively monitor major loads and loads with large fluctuations in the region.

[0054] (2) The COI frequency of the region is calculated based on the frequency data of the key bus in the region. The key bus monitored in region 2 includes bus 2 and bus 15. Therefore, the frequency of the equivalent motor in this region can be approximately calculated as follows: in, and This indicates the frequency deviation between monitored busbars 2 and 25; and This represents the weighting coefficient, which can be assigned values ​​based on the relative fluctuation range of the frequency and the quality of the monitoring data, such as... var indicates that the variance is calculated.

[0055] (3) The mechanical power change of the equivalent generator is compared with the structure of the synchronous generator speed governor as a primary frequency regulation link based on frequency feedback.

[0056] Step 2 above analyzes the monitoring data and treats the regional power grid as a synchronous generator unit. The overall inertia of the entire regional power grid is used as the inertia of the equivalent generator, which facilitates the subsequent derivation of the relationship between the regional inertia and the equivalent generator model of the region, and simplifies the calculation complexity of the regional inertia.

[0057] Step 3: Based on the obtained equivalent generator model, establish a parameterized transfer function model that reflects the dynamic relationship between the generator frequency deviation and the change in terminal electromagnetic power, and obtain the unit step response of the model. Further calculate the functional relationship between parameters such as the equivalent generator inertia and the unit step response of the model. Preferably, in step 3, a parameterized transfer function model reflecting the dynamic relationship between the equivalent generator frequency deviation and the change in terminal electromagnetic power is established based on the obtained equivalent generator model. The specific process is as follows: A primary frequency regulation stage is added to the regional equivalent generator model established in step 2, and it is then converted into a transfer function form, as follows: in, Represents the Laplace operator; This represents the unit power regulation coefficient (takes a negative value); This represents the transfer function that describes the dynamics of primary frequency modulation power regulation; This represents the transfer function.

[0058] Step 3 involves obtaining the unit step response of the model and further calculating the functional relationship between parameters such as the equivalent generator inertia and the model's unit step response, as follows: (1) According to the final value theorem, calculate the initial slope of the unit step response curve corresponding to the model: in, The unit step response representing the derivative of the equivalent motor frequency deviation; This represents the unit step response of the equivalent motor frequency deviation; it is worth noting that the Laplace transform corresponding to the differential element is s. Numerically satisfies (This indicates the generator's ability to track power commands.) Based on the above calculations, the functional relationship between the initial slope of the unit step response curve of the model and the regional inertia can be obtained.

[0059] (2) According to the final value theorem, calculate the steady-state value of the unit step response curve corresponding to the model: Based on the above calculations, the functional relationship between the model's unit step response curve and the regional primary frequency modulation coefficient and damping coefficient can be obtained.

[0060] The above steps elucidate the relationship between the overall inertia of the regional power grid and the unit step response of the equivalent model of the regional power grid, greatly simplifying the calculation process of the system's regional inertia.

[0061] Step 4: Based on the frequency deviation data and terminal electromagnetic power variation data of the equivalent generator in the region obtained in Step 2, the parameters of the frequency deviation and electromagnetic power deviation model constructed in Step 3 are identified through processes such as constructing forward and backward historical data Hankel matrices using the subspace method, establishing numerical relationships between matrices, constructing correlation matrices using orthogonal projection to reduce noise influence, singular value decomposition, and system parameter matrix estimation. The unit step response of the model is then obtained, and the identification result of the regional inertia is finally obtained based on the functional relationship between the inertia and the unit step response obtained in Step 3.

[0062] As a preferred embodiment, the standard discrete form structure of the parameterized model to be identified in step 4 is as follows: in, , , and This represents the coefficient matrix corresponding to the discrete system. This indicates the error caused by noise; This represents the coefficient matrix corresponding to the error term; the purpose of subspace identification is to... and Identify the corresponding coefficient matrix based on the measurement results over a period of time.

[0063] In step 4, during the subspace algorithm's identification of the frequency deviation and electromagnetic power change model, the calculated regional equivalent generator frequency deviation needs to be used as the output of the model to be identified, corresponding to the value obtained in step 2. , recorded as y The electromagnetic power variation data of the equivalent generators in the monitored area are used as input to the model to be identified, corresponding to the values ​​obtained in step 2. , recorded as u .

[0064] In step 4, the subspace algorithm is executed to construct the Hankel matrix corresponding to the forward and backward historical data, as follows: in, and Corresponding to the input respectively The Hankel matrix composed of data from different time periods; and Corresponding to the input respectively The Hankel matrix composed of data from different time periods; and This represents the number of rows and columns of a matrix, in this example. The value is 500; The value is 1000.

[0065] In step 4, during the execution of the subspace algorithm, the numerical relationship between the Hankel matrix and the parameter matrix of the model to be identified is established, as follows: in, and It is similar in form and The error term Hankel matrix; , , All are formed by transforming and combining the model coefficient matrix; , The forms are respectively with , Similar (the difference lies in the coefficients) and ).

[0066] Step 4 involves constructing a correlation matrix using orthogonal projection during the subspace identification algorithm to reduce the impact of noise, specifically as follows: based on the noise... With input data The weak correlation between them can be reduced by the following transformation to decrease the noise effect: in, ; Representation and matrix Orthogonal complement; Indicates about The orthogonal projection operator of a matrix can be calculated using RQ decomposition.

[0067] Step 4 involves processing the constructed matrix during the execution of the subspace identification algorithm. The transformation and singular value decomposition are performed as follows: A weighting coefficient matrix is ​​introduced using the N4SID method. and Further construct the matrix Then to Perform singular value decomposition, and based on the resulting singular value matrix, approximate it into a block matrix form: Step 4 involves estimating the model parameter matrix based on the singular value decomposition results during the execution of the subspace identification algorithm. The specific process is as follows: (1) Estimate the correlation matrix of the model parameters defined above based on the singular value decomposition results. , recorded as ,as follows: (2) According to For model parameters Calculations are performed: Since the actual system is a single-input single-output system, therefore The first line in ,Pick The first line is the identification result. .

[0068] (3) According to For model parameters Perform calculations: remove The first line is recorded as Remove the last line and denote it as ,according to ,according to Using least squares estimation, we can obtain .

[0069] (4) According to For model parameters and Calculation: First, transform the numerical relationships between the above matrices to obtain: in, express The orthogonal complement satisfies ; for The generalized inverse matrix; only The quantity to be determined. In the above equation... middle and The elements are linearly related to each other, according to and The estimation results were finally estimated using the least squares method. and .

[0070] Step 4, executing the subspace identification algorithm, finally yields the following system model: in: Represents state variables; This represents the input variable, specifically the change in the terminal electromagnetic power of the equivalent generator in the corresponding region. This represents the output variable, specifically the frequency deviation of the equivalent generator in the corresponding region.

[0071] Step 4 involves calculating the unit step response of the system model obtained from the above identification, as follows: (1) At the initial moment, i.e. k When =0, let ;and ,right k Output variable at time =0 and state variables Perform calculations; (2) In k At time =1, let , combined calculate and ; (3) In subsequent moments Always keep it at 1, iteratively calculate at each time step. and ; (4) Finally obtained The values ​​taken at each time point are the unit step response curves of the model.

[0072] As described in step 4, the system inertia and other parameters are calculated based on the unit step response of the identified model, as follows: (1) Let the steady-state value of the unit step response curve of the obtained model be denoted as _____. Based on the functional relationship between the model step response obtained in step 3 and the equivalent motor primary frequency regulation coefficient and damping coefficient, the following can be derived: ; (2) Let the unit step response curve of the obtained model be in The mean slope within the phase is Based on the functional relationship between the model step response and the equivalent motor inertia obtained in step 3, we can derive: .

[0073] This simulation experiment fully demonstrates the feasibility and beneficial effects of the method described in this invention. It combines power grid area division with inertia identification, improves the accuracy of area rotational inertia identification, and is of great significance for the subsequent frequency stability control of the power system.

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

Claims

1. A method for identifying system region inertia based on power grid region division, characterized in that, include: Monitor the frequency data of busbars at key locations in the power grid, decompose the oscillation components of each busbar frequency and extract key components based on the Prony algorithm, and divide the power grid into regions by comparing the main oscillation components of the busbar frequency at each location with the electrical distance. The boundary connection lines of the divided area, the load power within the area, and the bus frequency data within the area are monitored and processed. Based on the monitoring information, the divided area is equivalent to a synchronous generator. A dynamic relationship model between the electromagnetic power change and frequency deviation of the equivalent generator is constructed, and the relationship between the model and the inertia of the equivalent generator is derived. Constructing the dynamic relationship model involves adding a primary frequency regulation stage to the equivalent generator and transforming it into a transfer function form, as follows: in, Represents the Laplace operator; This represents the unit power regulation coefficient, and is taken as a negative value. This represents the transfer function that describes the dynamics of primary frequency modulation power regulation; Represents the transfer function; The relationship between the derived model and the equivalent generator inertia includes: According to the final value theorem, calculate the initial slope of the unit step response curve corresponding to the model: in, The unit step response representing the derivative of the equivalent motor frequency deviation; The unit step response represents the equivalent motor frequency deviation; s is the Laplace transform corresponding to the differential element. ; Based on the above calculations, the functional relationship between the initial slope of the unit step response curve of the model and the regional inertia can be obtained. According to the final value theorem, calculate the steady-state value of the unit step response curve corresponding to the model: Based on the above calculations, the functional relationship between the model's unit step response curve and the regional primary frequency modulation coefficient and damping coefficient can be obtained. The parameters of the dynamic relationship model are identified using the subspace method, and the regional inertia is estimated based on the identified model.

2. The system region inertia identification method based on power grid region division as described in claim 1, characterized in that, Monitoring bus frequency data at key locations includes: monitoring as many bus frequency data points as possible within the power grid as possible, using PMU equipment, with data sampling intervals of [missing information]. The total length of the data collection is .

3. The system region inertia identification method based on power grid region division as described in claim 2, characterized in that, The decomposition and key component extraction of the frequency oscillation components of each bus based on the Prony algorithm includes: First, a parameterized model for signal decomposition is established; Specifically, the frequency monitoring data of a bus is decomposed into a sum of multiple damped oscillation components at each sampling point: in, This represents the signal value at the nth moment, specifically the nth value in the bus frequency monitoring data; Indicates the order of the oscillatory components contained therein; and These represent the initial amplitude and phase of the k-th oscillation component, respectively. and This represents the attenuation coefficient and angular frequency; and These are variables introduced to describe the characteristics of each oscillation component for the convenience of subsequent analysis; Indicates a noise signal; A sample function matrix is ​​constructed using monitoring data, as follows: in, The total amount obtained from monitoring Bus system frequency data at each moment; Represents the sample function matrix; This represents the result of multiplying the original data, corresponding to the matrix. Elements in; The coefficients corresponding to each oscillation mode are determined using the least squares method, specifically including... ,in The solution process for the remaining coefficients is as follows: Based on the above sample function matrix, the equation can be written as follows: in In the formula It is known that it contains A linear equation and For each unknown, the coefficients corresponding to the oscillation modes are determined using the least squares method. Solve the problem and record the result as follows: ; The attenuation coefficient and oscillation frequency parameters of each oscillation mode are calculated based on the coefficients corresponding to each oscillation mode, as follows: Based on the coefficients corresponding to the above oscillation modes and its decaying oscillation characteristic parameters for each oscillation mode The relationship between them: right The parameters are calculated, and the result is recorded as... ; Based on oscillation characteristic parameters The initial amplitude and phase parameters of each signal, including each oscillation mode, are calculated as follows: Based on the recursive difference equation and The estimated values ​​of the original data sequence are calculated as follows: Given Continuous recursion ; Based on the relationship between the constructed original data sequence and the oscillation characteristic parameters: as well as ,right The estimation is performed, and the result is recorded as follows: ; According to parameters and The calculation results for the parameters , , and (k=1,2… The calculation is performed as follows: Here, Im and Re represent taking the real part and the imaginary part, respectively.

4. The system region inertia identification method based on power grid region division as described in claim 3, characterized in that, The division of the power grid into regions includes: For a single bus, the frequency oscillation component represents the initial amplitude of the oscillation component. and Parameters, extract them Larger Group oscillation components; The results obtained by pairwise comparison of each busbar Group oscillation components, if In the group of oscillating components, there exists More than one group of oscillation components and If the parameters are similar, the two bus lines being compared pairwise will be grouped into the same group; where the determination... and The principle of similar parameters is used for calculation. and The relative deviation of the parameters, when less than the set deviation threshold If so, then the damped oscillation characteristics of the two sets of oscillation components are considered to be similar; Calculate the electrical distance between every two buses in each group. If the electrical distance between one bus A and two or more other buses in a group exceeds a set threshold, the following applies. When this happens, bus A is removed from the group, and the frequency oscillation component of bus A is compared with the frequency oscillation components of buses in other groups. The preset parameters are then appropriately reduced. Find a group with similar frequency oscillation characteristics and add bus A to that group to complete the partitioning; if no group can be found that meets the parameters... For grouping with similar group requirements, bus A and its nearby unmonitored busbars, generators, and loads are divided into a separate area.

5. The system region inertia identification method based on power grid region division as described in claim 4, characterized in that, include: Based on the regional division results, the power fluctuation data of the tie lines between each region, the frequency data of the key busbars within each region, and the key load data are monitored, and the frequency and terminal electromagnetic power information of the equivalent motor are calculated, as follows: Calculate the equivalent electromagnetic power of the motor based on the power of the inter-regional tie line and the load within the region: in, This represents the total number of connection lines between the calculation area and other areas; This indicates the total number of loads within the calculation area; Indicates the power of the tie line; This indicates the load power.

6. The system region inertia identification method based on power grid region division as described in claim 5, characterized in that, include: The COI frequency of the region is calculated based on the key bus frequency data within the divided region: in, This indicates the frequency deviation of the monitored busbar; Indicates the number of buses with frequency monitoring capabilities; This represents the weighting coefficient, which can be assigned values ​​based on the relative fluctuation range of the frequency and the quality of the monitoring data, such as... var indicates that the variance is calculated; The change in mechanical power of the equivalent generator is analogous to the structure of the synchronous generator governor, which is a primary frequency regulation link based on frequency feedback.

7. The system region inertia identification method based on power grid region division as described in claim 6, characterized in that, The identification of parameters of the dynamic relation model using the subspace method includes: The standard discrete form structure of the parameterized model to be identified is as follows: in, , , and This represents the coefficient matrix corresponding to the discrete system. Indicates the error caused by noise; This represents the coefficient matrix corresponding to the error term; In the process of identifying the frequency deviation and electromagnetic power change model by executing the subspace algorithm, the calculated regional equivalent generator frequency deviation is used as the output of the model to be identified, denoted as y; and the monitored regional equivalent generator electromagnetic power change data is used as the input of the model to be identified, denoted as u. The subspace algorithm is used to construct the Hankel matrices corresponding to the forward and backward historical data, as follows: in, and Corresponding to the input respectively The Hankel matrix composed of data from different time periods; and Corresponding to the input respectively The Hankel matrix composed of data from different time periods; and Indicates the number of rows and columns of a matrix; During the execution of the subspace algorithm, the numerical relationship between the Hankel matrix and the parameter matrix of the model to be identified is established, as follows: in, and It is similar in form and The error term Hnakel matrix; , , All are formed by transforming and combining the model coefficient matrix; , The forms are respectively with , similar; During the execution of the subspace identification algorithm, a correlation matrix is ​​constructed using orthogonal projection, as follows: based on noise With input data The weak correlation between them can be mitigated by the following transformation to reduce noise: in, ; Representation and matrix Orthogonal complement; Indicates about The orthogonal projection operator of a matrix can be calculated using RQ decomposition; During the execution of the subspace identification algorithm, the matrix constructed above is... The transformation and singular value decomposition are performed as follows: The weighting coefficient matrix is ​​introduced using the N4SID method. and Further construct the matrix Then to Perform singular value decomposition, and based on the resulting singular value matrix, approximate it into a block matrix form: During the execution of the subspace identification algorithm, the model parameter matrix is ​​estimated based on the above singular value decomposition results. The specific process is as follows: Estimate the correlation matrix of model parameters based on the singular value decomposition results. , recorded as ,as follows: according to For model parameters Calculations are performed: Since the actual system is a single-input single-output system, therefore The first line in ,Pick The first line is the identification result. ; according to For model parameters Perform calculations: remove The first line is recorded as Remove the last line and denote it as ,according to ,according to Using least squares estimation, we can obtain ; according to For model parameters and Calculation: First, transform the numerical relationships between the above matrices to obtain: in, express The orthogonal complement satisfies ; for The generalized inverse matrix; only For the quantity to be determined; in the above formula middle and The elements are linearly related to each other, according to and The estimation results were finally estimated using the least squares method. and ; The final system model obtained by executing the subspace identification algorithm is as follows: in: Represents state variables; This represents the input variable, specifically the change in the terminal electromagnetic power of the equivalent generator in the corresponding region. This represents the output variable, specifically the frequency deviation of the equivalent generator in the corresponding region. The unit step response of the system model obtained from the above identification is calculated as follows: At the initial moment, i.e., when k=0, let ;and For the output variable at k=0 and state variables Perform calculations; At time k=1, let , combined calculate and ; In subsequent moments Always keep it at 1, iteratively calculate at each time step. and ; Finally obtained The values ​​taken at each time point are the unit step response curves of the model.

8. The system region inertia identification method based on power grid region division as described in claim 7, characterized in that, The estimated regional inertia based on the identified model includes: Let the steady-state value of the unit step response curve of the obtained model be... Based on the functional relationship between the model's step response and the equivalent motor's primary frequency regulation coefficient and damping coefficient, we can derive: ; Let the unit step response curve of the obtained model be... The mean slope within the phase is Based on the functional relationship between the model step response and the equivalent motor inertia, we can derive: .

Images