A method for power system regional division considering the temporal and spatial distribution characteristics of frequency

Abstract

The present invention discloses a method for power system regional division that takes into account the frequency and time-space distribution characteristics. The method models the equipment side and the network side of the new energy power system respectively to obtain a closed-loop frequency response model of the system, decomposes the frequency response into multiple components based on the spectral decomposition of the system hysteresis matrix; analyzes the response characteristics such as the lowest frequency point and the average change rate of each frequency component to quantify the frequency and time-space distribution characteristics; measures the frequency of the PMU configuration node when disturbances occur at different locations of the system through simulation, and uses the measured frequency response curve after modal decomposition and an improved spectral clustering algorithm to construct a similarity matrix and a Laplace matrix to extract the eigenvector of the modal frequency; finally, uses the DBSCAN clustering algorithm to cluster frequency curves with high similarity to achieve system partitioning; the present invention can quickly and accurately provide real-time system partitioning results, ensuring the consistency of regional frequency responses under various disturbances.

Description

Technical Field

[0001] The present invention relates to the field of novel power system frequency and space distribution characteristics, and in particular to a power system area division method taking frequency and space distribution characteristics into account. Background Art

[0002] The frequency response of power systems exhibits temporal and spatial distribution. Renewable energy sources such as wind turbines and photovoltaics cannot actively provide inertia support to the system. The integration of these new energy sources into the system further skews the inertia distribution, making the spatial and temporal frequency distribution characteristics of the system more pronounced after disturbances. Existing inertia assessments are shifting from system-level to regional-level approaches. Accurately zoning the system is a prerequisite for regional inertia assessment. Therefore, a fully programmable, real-time zoning method is needed to quickly and accurately generate real-time system zoning results and ensure consistent regional frequency responses under various disturbances.

[0003] Most existing zoning methods divide regions based on the correlation between the frequency response characteristics of individual generators, grouping nodes with similar frequency responses into a single region. However, most of these methods fail to consider the system zoning after the integration of renewable energy sources. With the increasing number of renewable energy units in the system, this zoning method based on generator coherence is becoming increasingly limited. To ensure the accuracy of zoning results, it is necessary to consider the impact of system operating conditions and the location of disturbances on zoning, thereby implementing real-time zoning. By analyzing the spatiotemporal distribution of power system frequency responses, it is possible to identify the impact of disturbance type and location on the frequency response of each unit. Through simulation, the frequencies of PMU configuration nodes can be measured when disturbances occur at different system locations. Using the measured frequency response curves after modal decomposition, an improved spectral clustering algorithm is used to construct a similarity matrix and Laplace matrix, extracting the eigenvectors of the modal frequencies. Subsequently, the DBSCAN clustering algorithm is used to cluster and partition the system, grouping frequency curves with high similarity together to achieve system zoning. Summary of the Invention

[0004] The purpose of the present invention is to overcome the defects existing in the above-mentioned background technology. The present invention proposes a method for power system regional division taking into account the spatiotemporal distribution characteristics of frequency. The method takes into account the influence of system operation conditions and the location of disturbance on the partition. The PMU configuration node frequency when disturbance occurs at different positions of the system is measured through simulation. The measured frequency response curve after modal decomposition is used to construct a similarity matrix and Laplace matrix through an improved spectral clustering algorithm to extract the eigenvector of the modal frequency; then, the DBSCAN clustering algorithm is used to cluster and partition the system, and frequency curves with high similarity are clustered together to achieve system partitioning.

[0005] The present invention adopts the following technical solutions to solve the above technical problems:

[0006] A method for dividing power system regions taking into account the temporal and spatial distribution characteristics of frequency is provided, wherein the steps are as follows:

[0007] Step S1: First, the equipment side and the network side of the new energy power system are modeled separately to obtain the system closed-loop frequency response model. Secondly, the frequency response is decomposed into multiple components based on the spectral decomposition of the system hysteresis matrix:

[0008] It is the common mode frequency, and each machine remains consistent;

[0009] is the differential mode frequency, each machine oscillates relative to each other;

[0010] Step S2: for the modal frequencies decomposed in step S1, further analyzing the response characteristics of the lowest frequency point and the average rate of change of each frequency component, and quantifying the spatiotemporal distribution characteristics of the frequency;

[0011] Step S3: The PMU configuration node frequency when disturbances occur at different locations of the system is measured through simulation. The PMU is a synchronized phasor measurement device. The measured frequency response curve after modal decomposition is used to construct a similarity matrix and a Laplace matrix through an improved spectral clustering algorithm to extract the eigenvectors of the modal frequencies. Subsequently, the system is clustered and partitioned using the DBSCAN clustering algorithm, and frequency curves with high similarity are clustered together to achieve system partitioning. The DBSCAN is a density-based clustering algorithm.

[0012] Furthermore, the step S1 is specifically as follows:

[0013] Step S1-1: First, the equipment side and the network side of the new energy power system are modeled separately to obtain a closed-loop frequency response model of the system. The dynamics of the power generation equipment can be represented by the frequency input-active power output transfer function:

[0014]

[0015] in, represents the perturbation of a variable; and are active power vector and frequency vector respectively, which are and For elements, , is the number of devices; is the frequency-active power transfer function matrix of all devices. Its negative sign is because the active power output of the power generation equipment is positive when the frequency drops;

[0016] Step S1-2: The relationship between the network side frequency and active power can be described by the power flow equation, which can be expressed as:

[0017]

[0018] Where, is the system nominal frequency, is a diagonal matrix with the Laplace operator as the diagonal element, The network Laplace matrix to eliminate passive nodes;

[0019] Connect the device to the input and output of the network and superimpose the disturbance on the active power of the device , the closed-loop system frequency response can be obtained:

[0020] .

[0021] Furthermore, the frequency response of the closed-loop system in step S1-2 is decomposed as follows:

[0022] To facilitate the decomposition of the closed-loop system frequency response in equation (3), it is rewritten as:

[0023]

[0024] Where, For the unit matrix, For about The hysteresis matrix will be The spectral decomposition of is recorded as:

[0025]

[0026] in, They are The eigenvalues ​​and left and right eigenvectors of , usually with subscript The ordinal number representing the characteristic value is different from the ordinal number representing the device Differentiate, they are Function; 0, and 1 for The eigenvalues ​​and left and right eigenvectors of The eigenvalues ​​and eigenvectors are as follows:

[0027]

[0028] Substituting equations (5)-(6) into equation (4), the node frequency response is It can also be decomposed into n components:

[0029]

[0030] in, is the change in the global frequency response, which represents the consistent part of the frequency response of all nodes. For the The change in the frequency response of each node, is the total number of nodes, is the weight coefficient, Indicates the The weight coefficient of each device, is the power change, and respectively The left and right eigenvectors of For the The characteristic value of each node, Indicates the The transfer function of a device, represents the inverse transfer function of the system;

[0031] From the above formula, we can see that All elements in are the same, so it is the consistent part of the frequency response of each node, representing the global frequency, and the remaining components represent the difference in frequency distribution. Therefore, the first component is called the frequency common mode component, referred to as the common mode frequency, and the remaining components are collectively referred to as the frequency differential mode component, referred to as the differential mode frequency.

[0032] Furthermore, the specific process of quantifying the spatiotemporal distribution characteristics of the frequency in step S2 is as follows:

[0033] Step 2-1: During the period from the initial disturbance to the frequency minimum, the high-order common-mode frequency trajectory is very close to the second-order damped sine curve. Based on the second-order damped sine curve, a more accurate frequency characteristic of the common-mode frequency minimum and average rate of change can be obtained. Therefore, the device transfer function structure corresponding to the damped sine curve can be selected to unify and simplify the models of various devices. The transfer function structure can be expressed as:

[0034]

[0035] in, 、 as well as are the effective inertia, effective damping and effective frequency modulation coefficient, is the system common mode inertia;

[0036] Step S2-2, solve the unified structural parameters through the following parameter optimization problem, which is described as:

[0037]

[0038] in, is the initial moment of disturbance, The final moment, which can usually be selected as , is the moment of lowest frequency, For the i Devices at time t The change in active power at the moment, u Represents wind power and photovoltaic new energy equipment and nodes; For the i Devices at time t The change in active power at the moment, cm represents common mode; and Respectively represent i Device and i Frequency-active power transfer function of wind power and photovoltaic new energy equipment; It represents the variation of the common-mode frequency response, that is, the part of the frequency response that is consistent across all nodes;

[0039] Depend on , as well as The above optimization problem can be described as approximating the linear combination of these three items. , the optimal coefficients are , as well as ;

[0040] Step S2-3: After solving the optimization problem, the equipment is simplified into a unified structure as follows:

[0041]

[0042] Where, , as well as They are the system common mode inertia, common mode damping and common mode frequency modulation coefficient respectively; considering the power step disturbance:

[0043]

[0044] in, is the equivalent perturbation size, then Equation (10) can be rewritten as:

[0045]

[0046] In the formula The time domain analytical expression can be directly obtained and the lowest frequency point and the average rate of change can be calculated.

[0047] Furthermore, the lowest point of equation (12) in step S2-3 is solved as follows:

[0048] According to the relationship between the system common mode inertia, common mode damping and common mode frequency modulation coefficient in formula (12), the damping ratio, undamped oscillation frequency, damped oscillation frequency and attenuation coefficient of the second-order system are obtained as follows: , , , ;

[0049] but The moment of reaching the lowest point and the lowest point They are:

[0050]

[0051] Where, represents the initial power of the system, represents the gain coefficient of the system, represents the equivalent moment of inertia of the system, represents the attenuation coefficient of the system, Represents the damped oscillation frequency of the system.

[0052] Furthermore, the average rate of change of frequency over a certain period of time is calculated as follows:

[0053] Pick , is a positive integer, then hour, Average rate of change of frequency over time for:

[0054]

[0055] Where, represents the initial power of the system, represents the common mode inertia of the system, It represents the damping ratio of the system, and e is the base of the natural logarithm.

[0056] Furthermore, the specific process of clustering and partitioning the system in step S3 is as follows:

[0057] Step S3-1, first establish a Laplace matrix using a spectral clustering algorithm, then standardize the Laplace matrix, and obtain the eigenvalues ​​of the Laplace matrix to establish a characteristic matrix;

[0058] Step S3-2, optimizing it using Harris Eagle optimization algorithm;

[0059] Step S3-3: Apply the EDR algorithm to the frequency curves of all nodes in the region, perform deletion, insertion, and replacement operations on the frequency time series points of two nodes in the same region within the sampling time window, perform the EDR algorithm comparison on all nodes in the region one by one, and finally find the node with the least number of operations in the region; the frequency curve can best reflect the dynamic response characteristics of the frequency of each node in the region. The EDR algorithm is a redistribution edit distance algorithm, which is an algorithm used to measure the similarity between two time series;

[0060] Step S3-4: Finally, the DBSCAN clustering algorithm is used to cluster the frequency curves with small distance and high similarity together, thereby realizing system partitioning.

[0061] Furthermore, the characteristic matrix in step S3-1 is specifically shown as follows:

[0062]

[0063] Where: is the tie line admittance; node With node Conductivity of the lines between them; It is the imaginary unit, used to represent the imaginary part of a complex number; is an element of an adjacency matrix, used to represent the nodes in the power system a and b Is there a connection between them? represents the set of edges between nodes in the power system, and Represents nodes respectively a and b voltage; Representation node and nodes The line conductance between For nodes With node The line susceptance between them; is the normalized Laplace matrix; is the identity matrix; is the sum of the admittance matrices of the node and all its connected edges; is the degree matrix; is the number of nodes in the area.

[0064] Furthermore, the Harris Hawk optimization algorithm optimization process includes: first setting the initial population number and the maximum number of iterations, then traversing all clustering results, and saving the global maximum silhouette coefficient partition result.

[0065] Compared with the prior art, the present invention adopts the above technical solution and has the following beneficial effects:

[0066] (1) The present invention proposes a method for dividing power system regions by taking into account the frequency-time-space distribution characteristics. Taking into account the modal characteristics of the frequency response of the new energy power system, the frequency-time-space distribution characteristics are quantified by modal frequency decomposition. The response characteristics such as the lowest point and average change rate of each frequency component are further analyzed, and the time distribution characteristics of the modal frequency components in the power system are quantified.

[0067] (2) The present invention proposes a method for dividing power system regions taking into account the spatiotemporal distribution characteristics of frequency. The frequencies of PMU configuration nodes when disturbances occur at different locations of the system are measured through simulation. The measured frequency response curve after modal decomposition is used to construct a similarity matrix and a Laplace matrix through an improved spectral clustering algorithm to extract the eigenvectors of the modal frequencies.

[0068] (3) The present invention proposes a method for dividing power system regions taking into account the spatiotemporal distribution characteristics of frequency. The DBSCAN clustering algorithm is used to cluster frequency curves with high similarity to achieve system partitioning. The present invention provides a real-time partitioning method that can be fully programmed and can quickly and accurately provide real-time partitioning results of the system, ensuring the consistency of regional frequency responses under various disturbance conditions. BRIEF DESCRIPTION OF THE DRAWINGS

[0069] Figure 1 is a flow chart of the steps of the present invention;

[0070] Figure 2 This is a frequency closed-loop dynamic response diagram of the new energy power system of the present invention;

[0071] Figure 3 is a common mode frequency schematic diagram of the present invention;

[0072] Figure 4 is a differential mode frequency schematic diagram of the present invention;

[0073] Figure 5 This is the Harris Eagle algorithm optimization flow chart of the present invention. DETAILED DESCRIPTION

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

[0075] A method for dividing power system regions taking into account the temporal and spatial distribution characteristics of frequency is provided, wherein the steps are as follows:

[0076] Step S1: First, the equipment side and the network side of the new energy power system are modeled separately to obtain the system closed-loop frequency response model. Secondly, the frequency response is decomposed into multiple components based on the spectral decomposition of the system hysteresis matrix:

[0077] It is the common mode frequency, and each machine remains consistent;

[0078] is the differential mode frequency, each machine oscillates relative to each other;

[0079] Step S2: for the modal frequencies decomposed in step S1, further analyzing the response characteristics of the lowest frequency point and the average rate of change of each frequency component, and quantifying the spatiotemporal distribution characteristics of the frequency;

[0080] Step S3: The PMU configuration node frequency when disturbances occur at different locations of the system is measured through simulation. The PMU is a synchronized phasor measurement device. The measured frequency response curve after modal decomposition is used to construct a similarity matrix and a Laplace matrix through an improved spectral clustering algorithm to extract the eigenvectors of the modal frequencies. Subsequently, the system is clustered and partitioned using the DBSCAN clustering algorithm, and frequency curves with high similarity are clustered together to achieve system partitioning. The DBSCAN is a density-based clustering algorithm.

[0081] Furthermore, the step S1 is specifically as follows:

[0082] Step S1-1: First, the equipment side and the network side of the new energy power system are modeled separately to obtain a closed-loop frequency response model of the system. The dynamics of the power generation equipment can be represented by the frequency input-active power output transfer function:

[0083] (1)

[0084] in, represents the perturbation of a variable; and are active power vector and frequency vector respectively, which are and For elements, , is the number of devices; is the frequency-active power transfer function matrix of all devices. Its negative sign is because the active power output of the power generation equipment is positive when the frequency drops;

[0085] Step S1-2: The relationship between the network side frequency and active power can be described by the power flow equation, which can be expressed as:

[0086] (2)

[0087] Where, is the system nominal frequency, is a diagonal matrix with the Laplace operator as the diagonal element, The network Laplace matrix to eliminate passive nodes;

[0088] Connect the device to the input and output of the network and superimpose the disturbance on the active power of the device , the closed-loop system frequency response can be obtained, such as Figure 2 As shown:

[0089] (3).

[0090] Furthermore, the frequency response of the closed-loop system in step S1-2 is decomposed as follows:

[0091] To facilitate the decomposition of the closed-loop system frequency response in equation (3), it is rewritten as:

[0092] (4)

[0093] Where, For the unit matrix, For about The hysteresis matrix will be The spectral decomposition of is recorded as:

[0094] (5)

[0095] in, They are The eigenvalues ​​and left and right eigenvectors of , usually with subscript The ordinal number representing the characteristic value is different from the ordinal number representing the device Differentiate, they are Function; 0, and 1 for The eigenvalues ​​and left and right eigenvectors of The eigenvalues ​​and eigenvectors are as follows:

[0096] (6)

[0097] Substituting equations (5)-(6) into equation (4), the node frequency response is It can also be decomposed into n components:

[0098] (7)

[0099] in, is the change in the global frequency response, which represents the consistent part of the frequency response of all nodes. For the The change in the frequency response of each node, is the total number of nodes, is the weight coefficient, Indicates the The weight coefficient of each device, is the power change, and respectively The left and right eigenvectors of For the The characteristic value of each node, Indicates the The transfer function of a device, represents the inverse transfer function of the system;

[0100] From the above formula, we can see that All elements in are the same, so it is the consistent part of the frequency response of each node, representing the global frequency, and the remaining components represent the difference in frequency distribution, so the first component is called the frequency common mode component, referred to as the common mode frequency. Figure 3 As shown, the remaining components are collectively referred to as frequency differential mode components, referred to as differential mode frequencies. Figure 4 shown.

[0101] Furthermore, the specific process of quantifying the spatiotemporal distribution characteristics of the frequency in step S2 is as follows:

[0102] Step 2-1: During the period from the initial disturbance to the frequency minimum, the high-order common-mode frequency trajectory is very close to the second-order damped sine curve. Based on the second-order damped sine curve, a more accurate frequency characteristic of the common-mode frequency minimum and average rate of change can be obtained. Therefore, the device transfer function structure corresponding to the damped sine curve can be selected to unify and simplify the models of various devices. The transfer function structure can be expressed as:

[0103] (8)

[0104] in, 、 as well as are effective inertia, effective damping and effective frequency modulation coefficient respectively;

[0105] Step S2-2, solve the unified structural parameters through the following parameter optimization problem, which is described as:

[0106] (9)

[0107] in, is the initial moment of disturbance, Generally, you can choose , is the moment of lowest frequency;

[0108] Depend on , as well as The above optimization problem can be described as approximating the linear combination of these three items. , the optimal coefficients are , as well as ;

[0109] Step S2-3: After solving the optimization problem, the equipment is simplified into a unified structure as follows:

[0110] (10)

[0111] Where, , as well as They are the system common mode inertia, common mode damping and common mode frequency modulation coefficient respectively; considering the power step disturbance:

[0112] (11)

[0113] in, is the equivalent perturbation size, then Equation (10) can be rewritten as:

[0114] (12)

[0115] In the formula The time domain analytical expression can be directly obtained and the lowest frequency point and the average rate of change can be calculated.

[0116] Furthermore, the lowest point of equation (12) in step S2-3 is solved as follows:

[0117] According to the relationship between the system common mode inertia, common mode damping and common mode frequency modulation coefficient in formula (12), the damping ratio, undamped oscillation frequency, damped oscillation frequency and attenuation coefficient of the second-order system are obtained as follows: , , , ;

[0118] but The moment of reaching the lowest point and the lowest point They are:

[0119] (13)

[0120] Where, represents the initial power of the system, represents the gain coefficient of the system, represents the equivalent moment of inertia of the system, represents the attenuation coefficient of the system, Represents the damped oscillation frequency of the system.

[0121] Furthermore, the average rate of change of frequency over a certain period of time is calculated as follows:

[0122] Pick , is a positive integer, then hour, Average rate of change of frequency over time for:

[0123] (14)

[0124] Where, represents the initial power of the system, represents the common mode inertia of the system, It represents the damping ratio of the system, and e is the base of the natural logarithm.

[0125] Furthermore, the specific process of clustering and partitioning the system in step S3 is as follows:

[0126] Step S3-1, first establish a Laplace matrix using a spectral clustering algorithm, then standardize the Laplace matrix, and obtain the eigenvalues ​​of the Laplace matrix to establish a characteristic matrix;

[0127] Step S3-2, optimizing it using Harris Eagle optimization algorithm;

[0128] Step S3-3: Apply the EDR algorithm to the frequency curves of all nodes in the region, perform deletion, insertion, and replacement operations on the frequency time series points of two nodes in the same region within the sampling time window, perform the EDR algorithm comparison on all nodes in the region one by one, and finally find the node with the least number of operations in the region; the frequency curve can best reflect the dynamic response characteristics of the frequency of each node in the region. The EDR algorithm is a redistribution edit distance algorithm, which is an algorithm used to measure the similarity between two time series;

[0129] Step S3-4: Finally, the DBSCAN clustering algorithm is used to cluster the frequency curves with small distance and high similarity together, thereby realizing system partitioning.

[0130] Furthermore, the characteristic matrix in step S3-1 is specifically shown as follows:

[0131] (15)

[0132] Where: is the tie line admittance; node With node Conductivity of the lines between them; It is the imaginary unit, used to represent the imaginary part of a complex number; represents the set of edges between nodes in the power system; Representation node and nodes The line conductance between For nodes With node The line susceptance between them; is the normalized Laplace matrix; is the identity matrix; is the sum of the admittance matrices of the node and all its connected edges; is the degree matrix; is the number of nodes in the area.

[0133] Furthermore, the Harris Hawk optimization algorithm optimization process includes: first setting the initial population number and the maximum number of iterations, then traversing all clustering results, and saving the global maximum silhouette coefficient partition result. The specific process is as follows: Figure 5 shown.

[0134] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the above embodiments, or replace some or all of the technical features therein with equivalents. However, these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for dividing power system regions taking into account the temporal and spatial distribution characteristics of frequency, characterized in that: The method comprises the following steps: Step S1: First, the equipment side and the network side of the new energy power system are modeled separately to obtain the system closed-loop frequency response model. Secondly, the frequency response is decomposed into multiple components based on the spectral decomposition of the system hysteresis matrix: Δω1(s) is the common mode frequency, which remains consistent for each machine; Δω k (s)(k=2,...,n) is the differential mode frequency, each machine oscillates relative to each other; Step S2: for the modal frequencies decomposed in step S1, further analyzing the response characteristics of the lowest frequency point and the average rate of change of each frequency component, and quantifying the spatiotemporal distribution characteristics of the frequency; Step S3: Measure the PMU configuration node frequencies when disturbances occur at different locations of the system through simulation. The PMU is a synchronized phasor measurement device. Using the measured node frequency response curves after modal decomposition, a similarity matrix and a Laplace matrix are constructed using an improved spectral clustering algorithm to extract the eigenvectors of the modal frequencies. Subsequently, the system is clustered and partitioned using the DBSCAN clustering algorithm, grouping modal frequency curves with high similarity together to achieve system partitioning. The DBSCAN algorithm is a density-based clustering algorithm. The specific process of clustering and partitioning the system in step S3 is as follows: Step S3-1, first establish a Laplace matrix using a spectral clustering algorithm, then standardize the Laplace matrix, and obtain the eigenvalues ​​of the Laplace matrix to establish a characteristic matrix; Step S3-2, optimizing it using Harris Eagle optimization algorithm; Step S3-3: Apply the EDR algorithm to the frequency curves of all nodes in the region, perform deletion, insertion, and replacement operations on the frequency time series points of two nodes in the same region within the sampling time window, perform the EDR algorithm comparison on all nodes in the region one by one, and finally find the node with the least number of operations in the region; the frequency curve can best reflect the dynamic response characteristics of the frequency of each node in the region. The EDR algorithm is a redistribution edit distance algorithm, which is an algorithm used to measure the similarity between two time series; Step S3-4: Finally, the DBSCAN clustering algorithm is used to cluster the frequency curves with small distance and high similarity together, thereby realizing system partitioning.

2. The method for dividing power system regions taking into account the temporal and spatial distribution characteristics of frequency according to claim 1, characterized in that: The step S1 is specifically as follows: Step S1-1: First, the equipment side and the network side of the new energy power system are modeled separately to obtain a closed-loop frequency response model of the system. The dynamics of the power generation equipment can be represented by the frequency input-active power output transfer function: ΔP(s)=-G(s)Δω(s) (1) Among them, Δ represents the perturbation of the variable; ΔP and Δω are the active power vector and frequency vector respectively, which are expressed as ΔP i and Δω i is an element, i = 1, 2, ..., n, where n is the number of devices; G(s) is the frequency-active power transfer function matrix of all devices, and its negative sign is because the active power output of the power generation equipment is positive when the frequency drops; Step S1-2: The relationship between the network side frequency and active power can be described by the power flow equation, which can be expressed as: ΔP(s)=ω0L r s -1 Give (s) (2) Where ω0 is the nominal frequency of the system, s is a diagonal matrix with the Laplace operator as the diagonal element, and L r The network Laplace matrix to eliminate passive nodes; Connect the device to the input and output of the network and superimpose the disturbance ΔP on the active power of the device D =[ΔP D,1 ,...,ΔP D,n ] T , the closed-loop system frequency response can be obtained: Dω(s)=(G(s)+ω0L r s -1 ) -1 ΔP D (s) (3).

3. The method for dividing power system regions taking into account the temporal and spatial distribution characteristics of frequency according to claim 2, characterized in that: The frequency response of the closed-loop system in step S1-2 is decomposed as follows: To facilitate the decomposition of the closed-loop system frequency response in equation (3), it is rewritten as: Where I is the unit matrix, H(s) is the hysteresis matrix about Δω(s), and the spectral decomposition of H(s) is recorded as: Among them, λ k (s),v k (s),u k (s) are the eigenvalues ​​and left and right eigenvectors of H(s), respectively. k = 1, ..., n. The subscript k represents the ordinal number of the eigenvalue, which is distinguished from the ordinal number of the device. They are all functions of s. 0, ω / n, and 1 are L r The eigenvalues ​​and left and right eigenvectors of ; accordingly, it is easy to know that H(s) has the following eigenvalues ​​and eigenvectors: Substituting equations (5)-(6) into equation (4), the node frequency response Δω(s) can also be decomposed into n components: Among them, Δω1(s) is the change of the global frequency response, which represents the consistent part of the frequency response of all nodes, and Δω k (s) is the change in the frequency response of the kth node, n is the total number of nodes, w is the weight coefficient, and w i Represents the weight coefficient of the i-th device, ΔP D (s) is the power change, u k (s) and v k (s) are the left and right eigenvectors of H(s), λ k (s) is the eigenvalue of the kth node, G i (s) represents the transfer function of the i-th device, G -1 (s) represents the inverse transfer function of the system; It can be seen from the above formula that all elements in Δω1 are the same, so it is the consistent part of the frequency response of each node, representing the global frequency, and the remaining components represent the frequency distribution differences. Therefore, the first component is called the frequency common mode component, referred to as the common mode frequency, and the remaining components are collectively referred to as the frequency differential mode component, referred to as the differential mode frequency.

4. The method for dividing power system regions taking into account the temporal and spatial distribution characteristics of frequency according to claim 1, characterized in that: The specific process of quantizing the spatiotemporal distribution characteristics of the frequency in step S2 is as follows: Step 2-1: During the period from the initial disturbance to the frequency minimum, the high-order common-mode frequency trajectory is very close to a second-order damped sine curve. Based on the second-order damped sine curve, a more accurate frequency characteristic of the common-mode frequency minimum and average rate of change can be obtained. Therefore, the device transfer function structure corresponding to the damped sine curve can be selected to unify and simplify the models of various devices. The transfer function structure can be expressed as: Among them, G u 、D u and K u are effective inertia, effective damping and effective frequency modulation coefficient, J u is the system common mode inertia; Step S2-2, solve the unified structural parameters through the following parameter optimization problem, which is described as: Among them, t0 is the initial time of disturbance, t f For the final moment, select 1.5t nadir , t nadir is the moment of lowest frequency, ΔP u,cm,i (t) is the active power change of the i-th device at time t, u represents the wind power and photovoltaic new energy equipment and nodes; ΔP cm,i (t) is the active power change of the i-th device at time t, cm represents the common mode; G i (s) and G u,i (s) represents the frequency-active power transfer function of the i-th device and the i-th wind power and photovoltaic new energy devices respectively; Δω cm (s) represents the variation of the common-mode frequency response, i.e., the consistent part of the frequency response of all nodes; ΔP u，cm，i (s) is given by sΔω cm , Δω cm and Δω cm / s, the above optimization problem can be described as using the linear combination of these three terms to approximate ΔP cm，i (s), the optimal coefficients are J u,i , D u,i as well as Step S2-3: After solving the optimization problem, the equipment is simplified into a unified structure as follows: Where, J cm , D cm and K cm They are the system common mode inertia, common mode damping and common mode frequency modulation coefficient respectively; considering the power step disturbance: Where P0 is the equivalent perturbation size, then Equation (10) can be rewritten as: Where Δω′ cm (s) The time domain analytical expression can be directly obtained and the lowest frequency point and the average rate of change can be calculated.

5. The method for dividing power system regions taking into account the temporal and spatial distribution characteristics of frequency according to claim 4, characterized in that: The solution to the lowest point of equation (12) in step S2-3 is as follows: According to the relationship between the system common mode inertia, common mode damping and common mode frequency modulation coefficient in formula (12), the damping ratio, undamped oscillation frequency, damped oscillation frequency and attenuation coefficient of the second-order system are obtained as follows: σ=-ζω n ; Then Δω′ cm (t) The time when the lowest point is reached Δt nadir and the lowest point Δω nadir They are: Where P0 represents the initial power of the system, K us Represents the gain coefficient of the system, J us represents the equivalent moment of inertia of the system, σ represents the attenuation coefficient of the system, ω d Represents the damped oscillation frequency of the system.

6. The method for dividing power system regions taking into account the temporal and spatial distribution characteristics of frequency according to claim 5, characterized in that: The average rate of change of frequency over a period of time is calculated as follows: Take t1 = t nadir / n t , n t is a positive integer, then when 0<ζ≤3, the average rate of change of frequency from 0 to t1 is for: Where P0 represents the initial power of the system, J cm represents the common mode inertia of the system, ζ represents the damping ratio of the system, and e is the base of the natural logarithm.

7. The method for dividing power system regions taking into account the temporal and spatial distribution characteristics of frequency according to claim 1, characterized in that: The characteristic matrix in step S3-1 is specifically shown as follows: Y ab =G ab +jB ab Where: Y ab is the tie line admittance; G ab The line conductance between nodes a and b; j is an imaginary unit used to represent the imaginary part of a complex number; A ab is an element of the adjacency matrix, which is used to indicate whether there is a connection between nodes a and b in the power system; E represents the edge set between nodes in the power system, V a and V b Represent the voltages of nodes a and b respectively; W ab represents the line conductance between node a and node b; B ab is the line susceptance between node a and node b; L sym is the normalized Laplace matrix; I is the identity matrix; d ab is the sum of the admittance matrices of the node and all its connected edges; D is the degree matrix; i is the number of nodes in the region.

8. The method for dividing power system regions taking into account the temporal and spatial distribution characteristics of frequency according to claim 1, characterized in that: The Harris Hawk optimization algorithm optimization process includes: first setting the initial population number and the maximum number of iterations, then traversing all clustering results, and saving the global maximum silhouette coefficient partition result.

Images