Analysis of the Spatiotemporal Difference Mechanism of Power System Node Frequency and Frequency Protection Correction Methods

Abstract

This invention relates to the field of power system frequency protection technology, specifically to the analysis of the spatiotemporal frequency differences at power system nodes and a method for frequency protection correction. It addresses the shortcomings of existing single-machine theories, which fail to capture the spatiotemporal characteristics of frequency, and multi-machine theories, which lack time-domain analysis of dynamic frequency response, thus hindering the mechanistic explanation of the system's spatiotemporal frequency characteristics. First, a two-region equivalent frequency response model is established, equating the higher-order model to a combination of typical lower-order systems. Mode shape analysis is used to decouple the model, leading to analytical solutions. Second, the expressions for the frequency offset and the maximum value of the rate of change of frequency (ROCOF) in each region are analytically obtained. Finally, based on the timing and numerical differences of the frequency offset and the maximum ROCOF value, an action value correction strategy is proposed for frequency protection control. This invention achieves a refined time-domain solution for all three frequency components, including different frequencies or oscillation modes, improving analytical accuracy and solving the problem of accurate time-domain analysis of higher-order models in two regions.

Description

Technical Field

[0001] This invention relates to the field of power system frequency protection technology, specifically to a method for analyzing the spatiotemporal differences in power system node frequencies and for frequency protection correction. Background Technology

[0002] In traditional power systems, frequency regulation resources are relatively evenly distributed, with smooth frequency changes and small differences between nodes, and the system frequency can be characterized by the inertial center frequency. However, with the increasing penetration rate of new energy sources, the spatial distribution differences of frequency regulation resources such as system inertia and primary frequency regulation reserves are increasing, and the output of new energy units such as wind and solar power is volatile. This leads to a significant spatiotemporal distribution characteristic in the dynamic response of the system frequency during large-scale disturbances. In order to more accurately describe the dynamic changes of the system frequency, it is necessary to perform targeted modeling and calculation of the system to reflect its spatiotemporal distribution characteristics, thereby achieving better control and protection of the system.

[0003] Currently, the commonly used methods for calculating the frequency response of a system are numerical simulation and analytical methods. Since the calculation of the spatiotemporal frequency distribution characteristics of different nodes requires consideration of the influence of the network structure, the larger the power grid and the more complex the network structure, the higher the system coupling and the greater the difficulty of analysis. Therefore, most mainstream solutions employ numerical simulation. Numerical simulation can adapt to the calculation of the frequency response characteristics of various complex network nodes and has high accuracy. However, it cannot explain the mechanism of the generation of frequency spatiotemporal distribution characteristics and is not conducive to summarizing the general laws of the system's frequency spatiotemporal distribution characteristics, thus possessing certain limitations.

[0004] Analytical methods can be categorized into two types: those based on single-machine models and those based on multi-machine models. Currently, the most researched approach is based on single-machine models. While single-machine models and analytical methods can quickly calculate the system's average frequency response characteristics, their high aggregation and neglect of tie lines prevent them from accurately representing the spatial distribution of frequencies. For multi-machine model-based analytical methods, the difficulty in solving for the eigenvalues ​​of higher-order model parameters arises.

[0005] In summary, existing single-machine theories cannot reflect the spatiotemporal characteristics of frequency, while multi-machine theories, although some can intuitively obtain frequency curves with spatiotemporal distribution characteristics, lack time-domain analysis of frequency dynamic response and cannot achieve a mechanistic explanation of the system's spatiotemporal frequency characteristics. Summary of the Invention

[0006] The purpose of this invention is to provide a method for analyzing the spatiotemporal differences in node frequencies in power systems and for frequency protection correction. This method addresses the problems that existing single-machine theories cannot reflect the spatiotemporal characteristics of frequencies, and multi-machine theories lack time-domain analysis of dynamic frequency response, thus failing to provide a mechanistic explanation of the spatiotemporal characteristics of system frequencies.

[0007] The technical solution adopted by this invention to solve its technical problem is: a method for analyzing the mechanism of spatiotemporal frequency differences at power system nodes and correcting frequency protection, comprising the following steps:

[0008] S1. Construct a two-region equivalent frequency response model.

[0009] Based on the typical aggregated single-machine SFR model, a frequency response model of a two-region interconnected system is established by analyzing the power flow characteristics between regions.

[0010] S2. Order Reduction and Decoupling of Two-Region Equivalent Frequency Response Model

[0011] By breaking the closed-loop feedback and using equivalent substitution, the high-order frequency response model is equivalent to a combination of typical low-order systems. The modal analysis method is used to decouple the model, and then the time-domain analytical solution is performed.

[0012] S3. Solve for the expression of frequency shift and maximum ROCOF value in each region.

[0013] By analyzing the analytical formula of the regional frequency and combining the value range of typical parameters of the power system, the maximum expression of the frequency offset and the rate of change of frequency (ROCOF) of each region is obtained analytically.

[0014] S4. Based on the frequency offset and the time and value difference of the ROCOF maximum value, propose an action value correction strategy in frequency protection control.

[0015] Furthermore, in step S1, the matrix equation of the equivalent frequency response model for the two regions is: In the formula, H is the system inertial time constant matrix, D is the system damping coefficient matrix, K is the system nodal admittance parameter matrix, and ΔP d Let ΔP be the system load power change matrix. M Let be the system unbalanced acceleration power matrix, Δδ be the system phase angle change matrix, Δf be the system frequency change matrix, s be the complex variable in the Laplace transform, and ω0 be the rated angular velocity; In the formula, H1 is the inertial time constant of the equivalent generator in region 1, H2 is the inertial time constant of the equivalent generator in region 2, D1 is the equivalent damping coefficient of region 1, D2 is the equivalent damping coefficient of region 2, and F... H1 For the power ratio of the high-pressure cylinder of the prime mover in region 1, F H2 T represents the power ratio of the high-pressure cylinder of the prime mover in region 2. R1 Let T be the reheat time constant of the prime mover in region 1. R2 Let R1 be the reheat time constant of the prime mover in region 2, R2 be the droop coefficient in region 1, and R2 be the droop coefficient in region 2. ΔP d1 Let ΔP be the load jump power in region 1. d2Let Δf1 be the load change power of region 2, Δf2 be the frequency change of region 1, and Δf2 be the frequency change of region 2.

[0016] Furthermore, in step S1, the output active power table for each region is as follows: P G1 For the active power output of generator 1 in area 1, P G2 For the active power output of the generator in area 2, E q1 E represents the node voltage within the equivalent generator in region 1. q2 Y is the node voltage within the equivalent generator in region 2. 12 G represents the mutual admittance between equivalent nodes in the two regions. 11 For the equivalent node self-conductance of region 1, G 22 Let φ be the self-conductance of the equivalent node in region 2, and φ be the mutual admittance angle between the equivalent nodes in the two regions; Δδ = δ1 - δ2, where δ1 is the phase angle inside the equivalent generator in region 1, and δ2 is the phase angle inside the equivalent generator in region 2; simplifying equation (3) to ΔP G Let Δδ be the system active power output matrix, ΔP be the phase angle change matrix, and ΔP be the phase angle change matrix. G1 Let ΔP be the change in active power output in region 1. G2 Δδ1 represents the change in active power output in region two, Δδ2 represents the change in phase angle in region one, and Δδ2 represents the change in phase angle in region two; ΔP G =[ΔP G1 ,ΔP G2 ] T Δδ=[Δδ1,Δδ2] T In the formula, The initial phase angle of the system is δ 10 and δ 20 Δδ0=δ 10 -δ 20 In high-voltage power grid lines, the resistance of the line is much smaller than the reactance, and the conductance to ground is usually negligible. Therefore, φ = π. Equation (4) can be further simplified to... In the formula, k = E q1 E q2 |B 12 | (7), B 12 Let be the mutual susceptance between equivalent nodes in the two regions; for a single machine equivalent to region i, the rotor motion equation is: In the formula, ΔP Mi Let ΔP be the unbalanced acceleration power in region i. Gi Let ΔP be the external input power of region i; in equation (6) Gi Substituting into equation (8) and using matrix form, we obtain equation (1).

[0017] Furthermore, in step S2, ΔPM It can be broken down into a proportional element and a first-order inertia. The proportional element part is combined with the damping matrix D to form the equivalent damping matrix D', i.e. Using a closed-loop substitution method, the center frequency ΔF in the ASF model after the two regions are aggregated is used to replace Δf. i , Further reorganization of the two regions In the formula, X represents the aggregated H, D, and R. -1 Three parameters, l i S is the conversion factor for region i. i S represents the total rated capacity of the generators in region i. B Let i be the power base value for region i; after order reduction, equation (1) becomes In the formula, Further, in step S2, the phase angle change matrix Δδ in natural coordinates is transformed into the modal matrix y in modal coordinates using the modal matrix, and the motion equations are decoupled; specifically: first, assuming the system is viscous damped, according to equation (12), its mode shape equation is (K-λ 2 H)Φ=0(14); where Φ=[J1,J2] is the eigenvector matrix of equation (12), J1 and J2 are eigenvectors, and λ is the eigenvalue; the eigenvalue is obtained from equation (14). The corresponding feature vector is Secondly, compare Δδ in natural coordinates with the principal coordinates y = [y1, y2]. T The connection is established through Φ, i.e., Δδ=Φy(17); where y1 and y2 are two vibration modes; finally, substitute equation (17) into equation (12) and multiply by the transpose of Φ. T Get Φ T HΦys 2 +Φ T DΦys+Φ T KΦy=ω0Φ T ΔP F +ω0Φ T ΔP d (18), after simplification, its time-domain equation is: In the formula, M is a diagonal matrix containing the inertial constants of each region, C is a matrix of the combined coefficients of damping coefficients and inertial constants, Z is a diagonal matrix of the combined coefficients of inertial constants and nodal admittances, and F(t) is a column vector of unbalanced acceleration power and disturbance power. Expanding equation (19), we get... At this point, the equations change from two-degree-of-freedom vibration equations to two single-degree-of-freedom vibration equations, thus achieving decoupling.

[0018] Furthermore, in step S2, the time-domain analytical solution process involves using the concept of a two-machine converged SFR model for solution. Using the superposition theorem, in ΔP F and ΔP d The results are obtained by solving the problem separately when the two excitation forces act individually. and Then add them together to get The time-domain analytical expression of the equivalent frequency response model for the two regions is as follows:

[0019] Furthermore, solve The steps are as follows: Aggregate the two-region equivalent frequency response models into a low-order simplified SFR model, with the transfer function as: In the formula, ω n1 Let y1 be the natural angular frequency of mode y1. The damping coefficient of mode y1 is given by equation (23). The inverse Laplace transform of equation (23) yields the following results. ω r1 Let be the damped oscillation angular frequency of mode y1, and α be the coefficient of the damped oscillation term. In the formula, The initial phase angle is the derivative of mode y1.

[0020] Furthermore, solve The steps are as follows: First, ΔP d When used alone, it performs an inverse Laplace transform. In the formula, Φ ij For the elements in the transformation matrix Φ, ωn2 is the natural angular frequency of mode y2, ζ is the damping coefficient of mode y2, and ω r2 The damped oscillation angular frequency is obtained by performing an inverse Laplace transform on equation (27). Then, in equation (21), ΔP F Solving for the case when it acts alone. Simplify equation (30) to

[0021] In the formula, M is the common factor derived from factorization, and A is the factor containing ω after factorization. n1 The coefficient of the s term in the numerator, B, is the factor containing ω after factorization. n1 The coefficient of the constant term in the numerator, C, is the factor containing ω after factorization. n2 The coefficient of the s term in the numerator, D, is the factor containing ω after factorization. n2 Let the coefficient of the constant term in the numerator be... Performing an inverse Laplace transform on equation (31) yields: In the formula, ρ1 is the amplitude of the oscillation of factor Y1(t), and ρ2 is the amplitude of the oscillation of factor Y2(t). The initial angles for the synthesis of trigonometric functions from factor Y1(t) Assuming the initial angles of the trigonometric functions composed of factor Y1(t), we can solve the above equations simultaneously to obtain...

[0022] In the formula, γ1 is the modal component y 2F In the derivative of ω r1 γ is the amplitude of the factor of the damping angular frequency, and γ2 is the modal component y. 2F In the derivative of ω r2 The amplitude of the factor of the damping angular frequency.

[0023] Furthermore, in step S3, the maximum value of ROCOF is calculated as follows: the moment when ROCOF is at its maximum in the region opposite to the disturbance is... To determine the offset angle after differentiation; for the disturbance occurrence region i, the maximum ROCOF value occurs at the initial moment of the disturbance occurrence; for the opposite region 3-i, the maximum ROCOF value occurs at the moment shown in equation (38); the maximum values ​​of the two are respectively ROCOF imax It is the maximum ROCOF value in region i, ROCOF (3-i)max It is the maximum ROCOF value in region 3-i; The offset angle is determined by differentiation.

[0024] Furthermore, in step S3, the maximum frequency offset is calculated as follows: Δf COI extreme point t COI With Δf di extreme point t di for Let k be a non-negative integer, and let t COI =t di have to k * To satisfy t COI =t di The nonnegative constant of the condition, k * After rounding down, we have [k] * At this point, we have: The maximum frequency offsets for each region are respectively

[0025] Δf i (t1) is the frequency change of region i at time t1, Δf i (t2) is the frequency change of region i at time t2, Δf i(t3) is the frequency change of region i at time t3, Δf 3-i (t1) is the frequency change of region 3-i at time t1, Δf (3-i) (t2) is the frequency change of region 3-i at time t2, Δf (3-i) (t3) is the frequency change of region 3-i at time t3.

[0026] Further, in step S4, the correction strategy is: for the maximum frequency offset, substitute the time of equation (42) into equation (22) for comparison, and take its maximum or minimum value according to the different disturbances; when the disturbance is positive, i.e. ΣΔP di When the value is greater than 0, the maximum value of the frequency deviation is taken. Multiply by 1.1-1.2 times the sensitivity, then add the initial frequency f. N This will then serve as the upper frequency limit f for this region. max When the disturbance is negative, i.e., ΣΔP di When <0, take the maximum value of its frequency deviation. Multiply by 1.1-1.2 times the sensitivity, then add the initial frequency f. N f, as the lower frequency limit of this region min ,Right now For frequency protection with frequency change as the threshold, disturbances are classified as follows: ① Disturbances occurring at low inertia H i_l Region ② The disturbance occurs at high inertia H i_h Region; for low inertia H il For the region, take 0 from ① according to formula (39) + The absolute value of the time point, multiplied by 1.1-1.2 times the sensitivity, is used as the maximum frequency change rate limit K for that region. Hi_l For high inertia H ih The region, according to formula (39) for ①, 0 + Time and t of ② max The values ​​at different times are compared, and the value with the larger absolute value in equation (39) is multiplied by 1.1-1.2 times the sensitivity to obtain the maximum frequency change rate limit K for that region. Hi_h ,Right now

[0027] The beneficial effects of this invention are as follows: Based on the generator electromagnetic equation and rotor motion equation, combined with the network parameter model, a frequency response model between regions is established. By employing closed-loop equivalent substitution and mode shape analysis, the order of the power grid frequency control system is reduced and the decoupling between various frequency components is achieved. Furthermore, based on the superposition theorem and the concept of equivalent convergence, a refined time-domain solution for all three frequency components containing different frequencies or oscillation forms is achieved, improving analytical accuracy and solving the problem of accurate time-domain analysis of high-order models in two regions. Based on the solved analytical expressions of the frequency response in the two regions, the frequency response characteristics of the two regions are analyzed, and the spatiotemporal distribution characteristics and related expressions of the maximum value of frequency offset and ROCOF are obtained. Based on these conclusions, a correction strategy for the action value of system frequency protection is proposed. Attached Figure Description

[0028] Figure 1 A system diagram of two interconnected regions;

[0029] Figure 2 The equivalent circuit diagram is for the Y type.

[0030] Figure 3 This is the Δ-type equivalent circuit diagram;

[0031] Figure 4 This is a diagram of a single-machine SFR model;

[0032] Figure 5 A diagram showing the equivalent frequency response model of the two interconnected regions;

[0033] Figure 6 Model order reduction processing Figure 1 ;

[0034] Figure 7 Model order reduction processing Figure 2 ;

[0035] Figure 8 For ΔP F Solution flowchart for when used alone;

[0036] Figure 9 Numerical simulation graphs for typical values;

[0037] Figure 10 This is a flowchart for frequency over-limit detection. Detailed Implementation

[0038] This invention first establishes a two-region equivalent frequency response model. By breaking the closed-loop feedback and using equivalent substitution, the high-order model is equivalent to a combination of typical low-order systems. Modal analysis is then used to decouple the model, allowing for analytical solutions and overcoming the challenge of analytically solving high-order systems. Secondly, through mathematical analysis of the regional frequency expressions, combined with the typical parameter ranges of the power system, the expressions for the maximum values ​​of frequency offset and ROCOF in each region are analytically obtained. The conclusion and corresponding mechanism explanation for the spatiotemporal differences in frequency offset and ROCOF across different regions are derived and discussed. Furthermore, based on the timing and numerical differences of the frequency offset and ROCOF maximum values, an action value correction strategy is proposed for frequency protection control. The specific methods of this invention include the following:

[0039] S1. Construct a two-region equivalent frequency response model.

[0040] To characterize the spatiotemporal frequency characteristics of multiple regions, a frequency response model of a two-region interconnected system is established based on a typical aggregated single-machine SFR model and by analyzing the power flow characteristics between regions. This provides a foundation for solving the analytical expression of the frequency response of the two regions and analyzing the spatiotemporal distribution characteristics.

[0041] (1) Equivalent model of regional interconnection

[0042] Two interconnected regional systems such as Figure 1 As shown in the figure, E qi δ represents the node voltage within the equivalent generator in the region. i U is the phase angle within the equivalent generator in the region. i and θ i These are the phase angles of the voltages at the interconnection ports of the two regions. Figure 1 The system shown can be equivalent to Figure 2 The Y-type equivalent circuit shown is transformed into the following: Figure 3 The diagram shows a delta-type equivalent circuit. The self-admittance and mutual admittance between regions are set as... In the formula, Y ii Y is the self-admittance of the equivalent node in region i. 12 G represents the mutual admittance between equivalent nodes in the two regions. ii For the equivalent node self-conductance of region i, B ii For the equivalent node self-susceptance of region i, G 12 B represents the mutual conductance between the equivalent nodes of the two regions. 12 φ is the mutual susceptance between the equivalent nodes of the two regions. i Let φ be the self-admittance angle of the equivalent node in region i, and φ be the mutual admittance angle between the equivalent nodes in the two regions. i0 For the admittance to ground on the i-th side of the Δ-type equivalent circuit node, y 12Let be the admittance between the two nodes of the delta-type equivalent circuit, j be the imaginary unit in the complex coordinate system, and i take values ​​of 1 and 2. The output power of the two regions satisfies the power flow equation P. Gi +jQ Gi =E qi ∑Y ij E qj i = 1, 2 (2), where P Gi For the active power output of generator i in region i, Q Gi Y represents the active and reactive power output of generator i in region i. ij For self-guided nano-Y ii Or mutual conductance Y 12 .

[0043] Let Δδ = δ1 - δ2, substitute it into equation (2) and retain the real part, then we get the expression for the output active power of each region as follows: The initial phase angles of the system are δ 10 δ 20 Linearization is performed at the system's operating point, δ1 = δ 10 ,δ2=δ 20 Δδ0=δ 10 -δ 20 Simplify the above formula to In the formula, ΔP G Both ΔP and Δδ are two-dimensional column vectors. G =[ΔP G1 ,ΔP G2 ] T Δδ=[Δδ1,Δδ2] T ΔP G Let ΔP be the system's active power output matrix. G1 Let ΔP be the change in active power output in region 1. G2 Let Δδ be the change in active power output in region two, Δδ be the system phase angle change matrix, Δδ1 be the phase angle change in region one, and Δδ2 be the phase angle change in region two; k 11 k 12 k 21 k 22 For ΔP G The elements k in the constant coefficient transformation matrix of Δδ 11 k is the parameter value in the first row and first column of the matrix. 12 k is the parameter value in the first row and second column of the matrix. 21 k is the parameter value in the second row and first column of the matrix. 22 This is the parameter value in the second row and second column of the matrix.

[0044] Because the resistance of a high-voltage power grid line is much smaller than its reactance, and its conductance to ground is usually negligible, therefore G 12=0, φ=π, further simplifying equation (4) to In the formula, k = E q1 E q2 |B 12 | (7).

[0045] (2) Interconnected area frequency response model

[0046] All generator sets within a region can be approximated as a single-unit model reflecting the average frequency response characteristics, i.e., a single-unit SFR model, whose frequency control block diagram is as follows: Figure 4 As shown. H i Let D be the inertial time constant of the equivalent generator in region i. i Let F be the equivalent damping coefficient of region i. Hi T represents the work ratio of the high-pressure cylinder of the prime mover in region i. Ri Let R be the reheat time constant of the prime mover in region i. i Let X be the adjustment coefficient for region i. The aggregation equation for all parameters is: X i =∑ m∈Vi k im X im (8); In the formula, X i The frequency regulation parameters of the equivalent generator portion after aggregation in region i are given. X im Let m be a partial frequency regulation parameter of the m-th generator within region i, i.e. F Him T represents the work ratio of the high-pressure cylinder of the prime mover of the m-th generator in region i. Rim Let H be the prime mover reheat time constant of the m-th generator in region i. im Let D be the inertial time constant of the m-th generator in region i. im R is the damping coefficient of the m-th generator in region i. im Let be the droop coefficient of the m-th generator in region i, and s be the complex variable in the Laplace transform. im and l im , where is the conversion factor for different parameters of region i. S im Let m be the rated capacity of the m-th generator in region i.

[0047] For a single machine with region i equivalent, the following is adopted: Figure 4 The low-order SFR model shown represents the frequency response in this region, and its rotor motion equation is: In the formula, ΔP Mi Let ΔP be the unbalanced acceleration power in region i. di Let ΔP be the load jump power in region i. Gi Let Δf be the external input power of region i.i Let ω0 be the frequency change in region i, and let ω0 be the rated angular velocity.

[0048] ΔP in equation (6) Gi Substituting into equation (11) and expressing it in matrix form, we obtain the matrix equation of the two-region equivalent frequency response model as follows: The corresponding equivalent frequency response control block diagram for the interconnection of the two regions is as follows: Figure 5 As shown in the formula, H is the system inertial time constant matrix, D is the system damping coefficient matrix, K is the system nodal admittance parameter matrix, and ΔP d Let ΔP be the system load power change matrix. M This represents the unbalanced acceleration power matrix of the system.

[0049] Thus, a two-region equivalent frequency response model has been established. Due to the presence of multiple closed loops, the system is a high-order complex system, making it difficult to solve the time-domain expression using traditional factorization combined with inverse Laplace transform. Figure 5 As shown by the red line, the interaction of factors such as tidal currents between regions leads to mutual coupling between Δf1 and Δf2, making it difficult to solve for Δf = f(t). To solve the problem of solving the above-mentioned high-order highly coupled model, the problem of system order reduction and decoupling is solved by equivalent closed-loop decomposition and mode shape analysis. The time-domain solution of the model parameter expression is obtained by separating variables and equivalent aggregation.

[0050] S2. Order Reduction and Decoupling of Two-Region Equivalent Frequency Response Model

[0051] S2.1, Model Order Reduction

[0052] S2.1.1 First, let ΔP M It is decomposed into the sum of a proportional element and a first-order inertia. The proportional element part is combined with the damping matrix D to form the equivalent damping matrix D', i.e. Its transformed form is as follows Figure 6 As shown.

[0053] S2.1.2, due to ΔP fi It contains a first-order inertial element and has low-pass filter properties, although the multi-region frequency response Δf i It contains higher-order harmonics, but the feedback power value ΔP of the synchronous speed controller is affected. fi Since the impact of the response in this region is very small, a closed-loop substitution method is adopted to break Δf. i With ΔP i The control closed loop uses the center frequency ΔF in the ASF model after the two regions are aggregated to replace Δf. i Its simplified equivalent diagram is as follows: Figure 7 As shown, the expression for ΔF is: In the formula, the parameters are still aggregated according to formulas (8) and (9), and the two regions are further sorted out to obtain In the formula, X represents the aggregated H, D, and R. -1 Three parameters, l i Then S is the conversion factor for region i. i S represents the total rated capacity of the generators in region i. B Let be the power base value for region i.

[0054] After processing, equation (12) becomes In the formula, S2.2 Transformation and Decoupling The above processing reduces the order of equation (17) to 6th order and transforms it into a combination of typical first and second order elements, reducing the complexity in form. However, the coupling problem between Δf1 and Δf2 is still not solved, and it is still difficult to solve it directly. From the form, equation (17) is a two-degree-of-freedom vibration equation with damping and excitation forces. Therefore, modal analysis is used to process it. Its essence is to use the modal matrix to transform the phase angle change matrix Δδ in natural coordinates into the modal matrix y in modal coordinates and achieve decoupling of the motion equation.

[0055] S2.2.1 First, assuming the system is viscous damped, we ignore damping and force to solve for eigenvalues ​​and eigenvectors. According to equation (17), its mode shape equation is (K-λ 2 H)Φ=0(19). Where Φ=[J1,J2] is the eigenvector matrix of equation (17), J1 and J2 are eigenvectors, and λ is the eigenvalue. The eigenvalue is obtained from equation (19). The corresponding feature vector is

[0056] S2.2.2 According to the transformation relationship, Δδ in natural coordinates and principal coordinates y = [y1, y2] T The connection is established through Φ, i.e., Δδ=Φy(22). In the formula, y1 and y2 are two vibration modes.

[0057] S2.2.3 Substitute transformation formula (22) into formula (17) and multiply by the transpose matrix Φ. T Get Φ T HΦys 2 +Φ T DΦys+Φ T KΦy=ω0Φ T ΔP F +ω0Φ T ΔP d (23). After simplification, its time-domain equation is: In the formula, M is a diagonal matrix containing the inertial constants of each region, C is a matrix of the combined coefficients of damping coefficients and inertial constants, Z is a diagonal matrix of the combined coefficients of inertial constants and nodal admittances, and F(t) is a column vector of unbalanced acceleration power and disturbance power.

[0058] From equation (25), it can be seen that C is a non-diagonal matrix, and there is still coupling between y1 and y2. This is inconsistent with the initial assumption that the system is viscous damped, i.e., C = εM + ηZ. In this case, it is impossible to continue solving for y1 and y2. However, if the system satisfies D1' / 2H1 = D2' / 2H2, then C in equation (25) is a diagonal matrix. At this time, the assumption holds, and equation (24) can also achieve decoupling. In actual large-scale power systems, damping is generally uniformly distributed, and this condition naturally holds in most cases. Therefore, the decoupling problem between y1 and y2 is solved. Under this condition, expanding (24) yields equation (26).

[0059] As can be seen, the equation has changed from a two-degree-of-freedom vibration equation to two single-degree-of-freedom vibration equations. y1 and y2 can be solved separately, and the coupling problem is solved.

[0060] S2.3 Time-domain analytical solution

[0061] In equation (26), variables y1 and y2 are decoupled, allowing for individual solutions. Furthermore, analytical expressions for Δδ1 and Δδ2 are obtained based on equation (22). To solve for the expressions for Δf1 and Δf2, equations (11) and (22) show that... Therefore Δf i The expression is solved and To obtain.

[0062] S2.3.1, Solve

[0063] The first equation in equation (26) is formally the same as the ASF model, so the idea of ​​a two-machine converged SFR model is used to solve it. In this way, the two-region equivalent frequency response models are converged into a low-order simplified SFR model, and the transfer function is: ω n1 Let y1 be the natural angular frequency of mode y1. Let be the damping coefficient of mode y1. Performing an inverse Laplace transform on equation (27), we obtain... ω r1 Let be the damped oscillation angular frequency of mode y1, and α be the coefficient of the damped oscillation term. In the formula, The initial phase angle is the derivative of mode y1.

[0064] Solve S2.3.2

[0065] The second equation in equation (26) is more complex in form. Using the superposition theorem, in ΔP... F and ΔP d The results are obtained by solving the problem separately when the two excitation forces act individually. and Then add them together to get

[0066] S2.3.2.1, First ΔP d When acting alone, perform the inverse Laplace transform, let Φ ij For the elements in the transformation matrix Φ, In the formula, ω n2 Let ω be the natural angular frequency of mode y2, ζ be the damping coefficient of mode y2, and ω be the frequency of mode y2. r2 Let be the damped oscillation angular frequency. Taking the inverse Laplace transform of equation (31) yields...

[0067] S2.3.2.2, then regarding ΔP in equation (26) F Solving for the case when it acts alone.

[0068] S2.3.2.3, Combining equations (15) and (18), it can be seen that the order of equation (34) is 6, therefore there are 6 factor terms, making direct solution quite difficult. Observing that the equation contains a summation part, and that the equation contains the central frequency expression of the ASF model, its form is similar to solving the SFR model, therefore the idea of ​​SFR model aggregation is adopted, such as... Figure 5 As shown, this makes equation (34) contain only two pairs of conjugate factors, reducing the order of the equation and thus greatly simplifying the calculation. In the equation, Therefore, equation (34) simplifies to In the formula, M is the common factor derived from factorization, and A is the factor containing ω after factorization. n1 The coefficient of the s term in the numerator, B, is the factor containing ω after factorization. n1 The coefficient of the constant term in the numerator, C, is the factor containing ω after factorization. n2 The coefficient of the s term in the numerator, D, is the factor containing ω after factorization. n2 The coefficient of the constant term in the numerator of the factor. Equation (36) simplifies the fourth-order expression into the form of a sum of two second-order factors. Now, we solve for the two second-order factors separately, letting... Performing an inverse Laplace transform on equation (36) yields: In the formula, ρ1 is the amplitude of the oscillation of factor Y1(t), and ρ2 is the amplitude of the oscillation of factor Y2(t). The initial angles for the synthesis of trigonometric functions from factor Y1(t) Let Y1(t) be the initial angle for the trigonometric function composition. Combining the above equations, we get... In the formula, γ1 is the modal component y 2F In the derivative of ω r1 γ is the amplitude of the factor of the damping angular frequency, and γ2 is the modal component y. 2F In the derivative of ω r2 The amplitude of the factor of the damping angular frequency.

[0069] At this point, Two components and All solutions have been completed. In summary, The time-domain analytical expression of the equivalent frequency response model for the two regions is obtained as follows: Equation (43) completes the analytical expression of the parameters of Δf1 and Δf2 in the time domain of the two regions, and characterizes the frequency dynamic characteristics of the interconnected system of the two regions from the overall level of the two regions. It is of great significance for the analysis of the generation mechanism of the frequency spatiotemporal distribution characteristics of the two regions and their engineering application.

[0070] S3. Solve for the expression of frequency shift and maximum ROCOF value in each region.

[0071] Thus, the time-domain analytical expression of the equivalent frequency response model for the two regions is obtained. The following analysis will be conducted at the time-domain level. Through mathematical derivation and the combination of typical value ranges, the spatiotemporal distribution differences of the maximum frequency offset and the maximum ROCOF value are obtained. Relevant expressions for the occurrence time and magnitude of the frequency offset and the maximum ROCOF value in different regions are given, and a correction strategy for frequency-related protection is proposed based on this.

[0072] S3.1, ROCOF Time Domain Analysis

[0073] As shown in equation (43), the frequency responses of the two regions have a symmetrical form, and their frequency responses are both based on the system center frequency superimposed with a certain oscillation frequency. Although the oscillation frequency is relatively complex in form, being a superposition of multiple decaying sine functions, the oscillation frequencies of the two regions have the same period, opposite oscillation directions, and the oscillation amplitude is inversely proportional to the equivalent moment of inertia of each region. Now let the three frequency components in equation (43) be respectively In the formula, Δf COI Let Δf be the system's inertial center frequency. di For region i and mode y 2d Correlated frequency components, ΔfFi For region i and mode y 2F Relevant frequency components.

[0074] When a disturbance occurs in region i, for that disturbed region i, Δf COI With Δf di and Δf Fi The signs are the same, but for the opposite region 3-i, Δf COI With Δf d(3-i) and Δf F(3-i) The signs before are reversed, which makes the difference in the early frequency response between the two regions more obvious. Therefore, taking the disturbance in region 2 as an example, we analyze ROCOF and the maximum frequency shift.

[0075] From equation (44) we get In the formula, The offset angle after differentiation is Δf. Fi The expression is complex and difficult to analyze directly using mathematics. Therefore, a numerical simulation approach is adopted for analysis. The numerical simulation of its typical values ​​is as follows: Figure 8 As shown in the image, it can be seen from the graph that Δf... Fi derivative The value is much smaller than Δf di derivative Its impact on ROCOF is minimal and can therefore be ignored. For the inertial center frequency Δf... COI derivative According to existing research, the maximum value of ROCOF at the inertial center of the system occurs at the initial moment of the disturbance.

[0076] for and According to equation (45), we know and It is an increasing function of ξ. According to equation (32), ξ = 0.5σλ² -1 , where σ=0.5H -1 (D+F H R -1 The value is less than 2 within the typical parameter range of the generator, while ω0 = 2πf N f N This is the system's rated frequency, therefore within the typical range of parameter values. ξ>>1, Approximating -π / 2, we get Equation (47) all reach their maximum values ​​at the initial time.

[0077] Based on the above analysis, when the disturbance occurs in region two, the initial time... and Since both have the same sign and reach their maximum modulus at the initial time, the ROCOF maximum value in region two occurs at the initial time. In the region opposite to where the disturbance occurs, at the initial time... and The signs are opposite, and according to formula (45) we get rate of change at initial time and rate of change at the initial time for

[0078] Therefore, the initial ROCOF of region one is 0. Clearly, the maximum ROCOF value in region one does not occur at the initial time. Since... and Both contain exponential decay terms, and according to equations (29) and (31), ω n2 >>ω n1 Regarding the index portion, Within the typical range, this ratio is close to 1, so the two exponential decay rates are not significantly different. However, for the trigonometric function part, the period of change T1 >> T2, therefore... The rate of change is much smaller than In summary, the changes in ROCOF in the early stages were mainly influenced by... The influence of the attenuation factor. And due to the existence of the attenuation factor, the maximum ROCOF value in region one must occur in... Within the first oscillation cycle, we approximate the first peak moment, i.e. The time when the ROCOF of this region reaches its maximum value is given by the time when the ROCOF of the region opposite the disturbance reaches its maximum value.

[0079] Therefore, for the disturbance region i, the maximum ROCOF value occurs at the initial moment of the disturbance, while for the opposite region 3-i, the maximum ROCOF value occurs at the moment shown in equation (49). The maximum values ​​for both are respectively...

[0080] ROCOF imax It is the maximum ROCOF value in region i, ROCOF (3-i)max It is the maximum ROCOF value in region 3-i.

[0081] S3.2, Time Domain Analysis of Maximum Frequency Offset

[0082] Traditional methods obtain the maximum frequency shift and its occurrence time by solving the extreme points of equation (44). However, as shown in equation (45), its differential equation is a transcendental equation, making it difficult to obtain the parameter solution for the maximum shift using this method. Based on the above analysis, Δf Fi The impact on the overall derivative is relatively small, so we will still only consider Δf. COIWith Δf di The decay rates of the exponential terms of the two are very close, while Δf COI Comparison with Δf di The trigonometric function terms change much more slowly, therefore the system's maximum offset occurs at Δf. di In Δf COI Among the 2-3 extreme points near the first extreme point. Therefore, according to equation (45), Δf COI extreme point t COI With Δf di extreme point t di for k is a non-negative integer. Let t COI =t di We can obtain, k * To satisfy t COI =t di The nonnegative constant of the condition. k * After rounding down, we have [k] * At this point, we have:

[0083] Assuming the disturbance occurs in region i, and β < 0 (i.e., active disturbance), when [k * When t1 is an odd number, it is a local minimum point, and the maximum deviation of the disturbance region i at this point is Δf. i (t1), the maximum deviation in region 3-i may be Δf (3-i) (t2) or Δf (3-i) (t3); when [k * When t1 is an even number, it is a local minimum point. At this point, the maximum frequency deviation of the perturbation region i may be Δf. i (t2) or Δf i (t3), the maximum deviation in region 3-i is Δf 3-i (t1). When β>0, i.e., there is active power disturbance, the situation is reversed. Therefore, the maximum frequency offset in each region is respectively Δf i (t1) is the frequency change of region i at time t1, Δf i (t2) is the frequency change of region i at time t2, Δf i (t3) is the frequency change of region i at time t3, Δf 3-i (t1) is the frequency change of region 3-i at time t1, Δf (3-i) (t2) is the frequency change of region 3-i at time t2, Δf (3-i) (t3) is the frequency change of region 3-i at time t3.

[0084] In summary, the frequency shift and ROCOF after the disturbance exhibit distinct spatiotemporal distribution characteristics. Regarding ROCOF, the ROCOF in the disturbance region reaches its maximum value at the initial moment, while the ROCOF at more distant spatial nodes reaches its maximum after half a cycle of oscillation. The maximum ROCOF value for each region can be calculated using equation (50). As for the maximum frequency shift, the spatiotemporal distribution difference of the maximum frequency shift over a wide area is relatively small, but the differences in the maximum shift time between regions are more significant. The maximum shift in each region occurs at the extreme point near the center frequency extreme point of the additional oscillation, and can be calculated using equations (54) and (55) according to the disturbance type.

[0085] S4. Frequency protection action value correction strategy

[0086] The above analysis shows that when a disturbance occurs in the power system, frequency offset and ROCOF exhibit spatiotemporal distribution characteristics in different regions, with significant differences in the timing and magnitude of the maximum frequency offset and ROCOF values. This difference can affect the system's frequency protection. If the entire network continues to use a uniform inertial center frequency for protection action settings, it may lead to maloperation or failure of protection devices.

[0087] The frequency-dependent protection detection process proposed in this invention is as follows: Figure 10 As shown, the system determines whether the frequency index exceeds the limit by detecting the offset of the tracking frequency on the converter output side before and after a fault, as well as the frequency change rate. First, the system frequency is acquired to obtain the system frequency change curve, and the ROCOF curve is obtained based on the frequency curve. These curves are then compared with the protection setpoints to determine whether the system frequency is operating within the set range. In the formula, f is the measured value of the system frequency. min f is the lower limit of frequency protection operation. max R is the upper limit of frequency protection action. f Here, f is the calculated rate of change of system frequency, and K is the setpoint for frequency change protection. When an out-of-limit condition is detected, i.e., f > f... max or f <f min When |R|>K, the relevant protection action will be activated.

[0088] Therefore, for frequency protection with the maximum frequency offset as the threshold, the frequency protection action value of the present invention is modified as follows:

[0089] For the maximum frequency offset, the time in equation (53) is substituted into equation (43) and compared. The maximum or minimum value is taken according to the different disturbances, as shown in equations (54) and (55). When the disturbance is positive, i.e., ΣΔP di When the value is greater than 0, the maximum value of the frequency deviation is taken. Multiply by 1.1-1.2 times the sensitivity, then add the initial frequency f. NThis will then serve as the upper frequency limit f for this region. max When the disturbance is negative, i.e., ΣΔP di When <0, take the maximum value of its frequency deviation. Multiply by 1.1-1.2 times the sensitivity, then add the initial frequency f. N f, as the lower frequency limit of this region min ,Right now

[0090] For frequency protection with frequency change as the threshold, under the condition that the overall system disturbance remains unchanged and there is a difference in inertia H between the two regions, there are two disturbance scenarios: ① The disturbance occurs at the lower inertia H. il Region; ② The disturbance occurs at high inertia H ih Region. According to equations (45) and (50), the low inertia H il The maximum value of ROCOF in the region must occur in case ① (0). + At that moment, and the high inertia H ih The maximum value of ROCOF in the region may occur at the initial time in case ① or at time t in case ②. max At any given moment, this is mainly related to the magnitude of the difference in inertia H between the two regions and the inertia H of the local region.

[0091] Therefore, the correction method for the frequency protection action value with frequency change as the threshold is as follows:

[0092] For low inertia H il For the region, take 0 from ① according to formula (50). + The absolute value of the time point, multiplied by 1.1-1.2 times the sensitivity, is used as the maximum frequency change rate limit K for that region. Hi_l For high inertia H ih For the region, apply equation (50) to ① for 0. + Time and t of ② max The values ​​at different times are compared, and the value with the larger absolute value is taken. This value is then multiplied by 1.1-1.2 times the sensitivity to obtain the maximum frequency change rate limit K for that region. Hi_h .

[0093] This invention establishes a frequency response model between regions based on the generator electromagnetic equation and rotor motion equation, combined with a network parameter model. By employing closed-loop equivalent substitution and mode shape analysis, it achieves order reduction and decoupling between frequency components in the power grid frequency control system. Furthermore, based on the superposition theorem and the concept of equivalent convergence, it realizes a refined time-domain solution for all three frequency components containing different frequencies or oscillation modes, improving analytical accuracy and solving the problem of accurate time-domain analysis of high-order models in two regions. Based on the solved analytical expressions for the frequency response of the two regions, the frequency response characteristics of the two regions are analyzed, and the spatiotemporal distribution characteristics and related expressions of the maximum value of frequency offset and ROCOF are derived. Based on these conclusions, a correction strategy for the operating value of system frequency protection is proposed.

Claims

1. A method for analyzing the spatiotemporal frequency differences at power system nodes and for frequency protection correction, characterized in that, Includes the following steps: S1. Construct a two-region equivalent frequency response model. Based on the typical aggregated single-machine SFR model, a frequency response model of a two-region interconnected system is established by analyzing the power flow characteristics between regions. S2. Order Reduction and Decoupling of Two-Region Equivalent Frequency Response Model By breaking the closed-loop feedback and using equivalent substitution, the high-order frequency response model is equivalent to a combination of typical low-order systems. The modal analysis method is used to decouple the model, and then the time-domain analytical solution is performed. S3. Solve for the expression of frequency shift and maximum ROCOF value in each region. By analyzing the analytical formula of the regional frequency and combining the value range of typical parameters of the power system, the maximum expression of the frequency offset and the rate of change of frequency (ROCOF) of each region is obtained analytically. S4. Based on the frequency offset and the time and value difference of the ROCOF maximum value, propose an action value correction strategy in frequency protection control.

2. The method for analyzing the spatiotemporal differences in power system node frequencies and correcting frequency protection according to claim 1, characterized in that, In step S1, the matrix equation of the two-region equivalent frequency response model is: In the formula, H is the system inertial time constant matrix, D is the system damping coefficient matrix, K is the system nodal admittance parameter matrix, and ΔP d Let ΔP be the system load power change matrix. M Let be the system unbalanced acceleration power matrix, Δδ be the system phase angle change matrix, Δf be the system frequency change matrix, s be the complex variable in the Laplace transform, and ω0 be the rated angular velocity; In the formula, H1 is the inertial time constant of the equivalent generator in region 1, H2 is the inertial time constant of the equivalent generator in region 2, D1 is the equivalent damping coefficient of region 1, D2 is the equivalent damping coefficient of region 2, and F... H1 For the power ratio of the high-pressure cylinder of the prime mover in region 1, F H2 T represents the power ratio of the high-pressure cylinder of the prime mover in region 2. R1 Let T be the reheat time constant of the prime mover in region 1. R2 Let R1 be the reheat time constant of the prime mover in region 2, R2 be the droop coefficient in region 1, and R2 be the droop coefficient in region 2. ΔP d1 Let ΔP be the load jump power in region 1. d2 Let Δf1 be the load change power of region 2, Δf2 be the frequency change of region 1, and Δf2 be the frequency change of region 2.

3. The method for analyzing the spatiotemporal differences in power system node frequencies and correcting frequency protection according to claim 2, characterized in that, In step S1, the output active power meters for each region are as follows: P G1 For the active power output of generator 1 in area 1, P G2 For the active power output of the generator in area 2, E q1 E represents the node voltage within the equivalent generator in region 1. q2 Y is the node voltage within the equivalent generator in region 2. 12 G represents the mutual admittance between equivalent nodes in the two regions. 11 For the equivalent node self-conductance of region 1, G 22 Let φ be the self-conductance of the equivalent node in region 2, and φ be the mutual admittance angle between the equivalent nodes in the two regions; Δδ = δ1 - δ2, where δ1 is the phase angle inside the equivalent generator in region 1, and δ2 is the phase angle inside the equivalent generator in region 2; simplifying equation (3) to ΔP G Let Δδ be the system active power output matrix, ΔP be the phase angle change matrix, and ΔP be the phase angle change matrix. G1 Let ΔP be the change in active power output in region 1. G2 Δδ1 represents the change in active power output in region two, Δδ2 represents the change in phase angle in region one, and Δδ2 represents the change in phase angle in region two; ΔP G =[ΔP G1 ,ΔP G2 ] T Δδ=[Δδ1,Δδ2] T In the formula, The initial phase angle of the system is δ 10 and δ 20 Δδ0=δ 10 -δ 20 In high-voltage power grid lines, the resistance of the line is much smaller than the reactance, and the conductance to ground is usually negligible. Therefore, φ = π. Equation (4) can be further simplified to... In the formula, k = E q1 E q2 |B 12 |(7), B 12 Let be the mutual susceptance between equivalent nodes in the two regions; for a single machine equivalent to region i, the rotor motion equation is: In the formula, ΔP Mi Let ΔP be the unbalanced acceleration power in region i. Gi Let be the external input power of region i; ΔP in equation (6) Gi Substituting into equation (8) and using matrix form, we obtain equation (1).

4. The method for analyzing the spatiotemporal differences in power system node frequencies and correcting frequency protection according to claim 3, characterized in that, In step S2, ΔP M It can be broken down into a proportional element and a first-order inertia. The proportional element part is combined with the damping matrix D to form the equivalent damping matrix D', i.e. Using a closed-loop substitution method, the center frequency ΔF in the ASF model after the two regions are aggregated is used to replace Δf. i , Further reorganization of the two regions In the formula, X represents the aggregated H, D, and R. -1 Three parameters, l i S is the conversion factor for region i. i S represents the total rated capacity of the generators in region i. B Let i be the power base value for region i; after order reduction, equation (1) becomes In the formula, 5. The method for analyzing the spatiotemporal differences in power system node frequencies and correcting frequency protection according to claim 4, characterized in that, In step S2, the phase angle transformation matrix Δδ in natural coordinates is transformed into the modal matrix y in modal coordinates using the modal matrix, and the motion equations are decoupled; specifically: first, assuming the system is viscous damped, according to equation (12), its mode shape equation is (K-λ 2 H)Φ=0(14); where Φ=[J1,J2] is the eigenvector matrix of equation (12), J1 and J2 are eigenvectors, and λ is the eigenvalue; The eigenvalues ​​obtained from equation (14) are: The corresponding feature vector is Secondly, compare Δδ in natural coordinates with the principal coordinates y = [y1, y2]. T The connection is established through Φ, i.e., Δδ=Φy(17); where y1 and y2 are two vibration modes; finally, substitute equation (17) into equation (12) and multiply by the transpose of Φ. T Get Φ T HΦys 2 +Φ T DΦys+Φ T KΦy=ω0Φ T ΔP F +ω0Φ T ΔP d (18), after simplification, its time-domain equation is: In the formula, M is a diagonal matrix containing the inertial constants of each region, C is a matrix of the combined coefficients of damping coefficients and inertial constants, Z is a diagonal matrix of the combined coefficients of inertial constants and nodal admittances, and F(t) is a column vector of unbalanced acceleration power and disturbance power. Expanding equation (19), we get... At this point, the equations change from two-degree-of-freedom vibration equations to two single-degree-of-freedom vibration equations, thus achieving decoupling.

6. The method for analyzing the spatiotemporal differences in power system node frequencies and correcting frequency protection according to claim 5, characterized in that, In step S2, the time-domain analytical solution is performed by using the idea of ​​a two-machine converged SFR model. Using the superposition theorem, in ΔP F and ΔP d The results are obtained by solving the problem separately when the two excitation forces act individually. and Then add them together to get The time-domain analytical expression of the equivalent frequency response model for the two regions is as follows:

7. The method for analyzing the spatiotemporal frequency differences of power system nodes and correcting frequency protection according to claim 6, characterized in that, Solve The steps are as follows: Aggregate the two-region equivalent frequency response models into a low-order simplified SFR model, with the transfer function as: In the formula, ω n1 Let y1 be the natural angular frequency of mode y1. The damping coefficient of mode y1 is given by equation (23). The inverse Laplace transform of equation (23) yields the following results. ω r1 Let be the damped oscillation angular frequency of mode y1, and α be the coefficient of the damped oscillation term. In the formula, The initial phase angle is the derivative of mode y1.

8. The method for analyzing the spatiotemporal frequency differences of power system nodes and correcting frequency protection according to claim 7, characterized in that, Solve The steps are as follows: First, ΔP d When used alone, it performs an inverse Laplace transform. In the formula, Φ ij For the elements in the transformation matrix Φ, ω n2 Let ω be the natural angular frequency of mode y2, ζ be the damping coefficient of mode y2, and ω be the frequency of mode y2. r2 Given the damped oscillation angular frequency, we obtain the inverse Laplace transform of equation (27). (29); then, regarding ΔP in equation (21) F Solving for the case when it acts alone. Simplify equation (30) to In the formula, M is the common factor derived from factorization, and A is the factor containing ω after factorization. n1 The coefficient of the s term in the numerator, B, is the factor containing ω after factorization. n1 The coefficient of the constant term in the numerator, C, is the factor containing ω after factorization. n2 The coefficient of the s term in the numerator, D, is the factor containing ω after factorization. n2 Let the coefficient of the constant term in the numerator be... Performing an inverse Laplace transform on equation (31) yields: In the formula, ρ1 is the amplitude of the oscillation of factor Y1(t), and ρ2 is the amplitude of the oscillation of factor Y2(t). The initial angles for the synthesis of trigonometric functions from factor Y1(t) Assuming the initial angles of the trigonometric functions synthesized from factor Y1(t), we can obtain the following equations: In the formula, γ1 is the modal component y 2F In the derivative of ω r1 γ is the amplitude of the factor of the damping angular frequency, and γ2 is the modal component y. 2F In the derivative of ω r2 The amplitude of the factor of the damping angular frequency.

9. The method for analyzing the spatiotemporal differences in power system node frequencies and correcting frequency protection according to claim 8, characterized in that, In step S3, the maximum value of ROCOF is calculated as follows: the moment when ROCOF is at its maximum in the region opposite to the disturbance is... To calculate the offset angle after differentiation; for the disturbance region i, where i takes the values ​​1 and 2, the maximum ROCOF value occurs at the initial moment of the disturbance; for the opposite region 3-i, the maximum ROCOF value occurs at the moment shown in equation (38); the maximum values ​​of the two are respectively ROCOF imax It is the maximum ROCOF value in region i, ROCOF (3-i)max It is the maximum ROCOF value in region 3-i; To determine the offset angle after differentiation; In step S3, the maximum frequency offset is calculated as follows: Δf COI extreme point t COI With Δf di extreme point t di for Let k be a non-negative integer, and let t COI =t di have to k * To satisfy t COI =t di The nonnegative constant of the condition, k * After rounding down, we have [k] * At this point, we have: The maximum frequency offsets for each region are as follows: Δf i (t1) is the frequency change of region i at time t1, Δf i (t2) is the frequency change of region i at time t2, Δf i (t3) is the frequency change of region i at time t3, Δf 3-i (t1) is the frequency change of region 3-i at time t1, Δf (3-i) (t2) is the frequency change of region 3-i at time t2, Δf (3-i) (t3) is the frequency change of region 3-i at time t3.

10. The method for analyzing the spatiotemporal differences in power system node frequencies and correcting frequency protection according to claim 9, characterized in that, In step S4, the correction strategy is as follows: for the maximum frequency offset, substitute the time in equation (42) into equation (22) for comparison, and take the maximum or minimum value according to the different disturbances; when the disturbance is positive, i.e. ΣΔP di When the value is greater than 0, the maximum value of the frequency deviation is taken. Multiply by 1.1-1.2 times the sensitivity, then add the initial frequency f. N This will then serve as the upper frequency limit f for this region. max When the disturbance is negative, i.e., ΣΔP di When <0, take the maximum value of its frequency deviation. Multiply by 1.1-1.2 times the sensitivity, then add the initial frequency f. N f, as the lower frequency limit of this region min ,Right now For frequency protection with frequency change as the threshold, disturbances are classified as follows: ① Disturbances occurring at low inertia H i_l Region ② The disturbance occurs at high inertia H i_h Region; for low inertia H i_l For the region, take 0 from ① according to formula (39) + The absolute value of the time point, multiplied by 1.1-1.2 times the sensitivity, is used as the maximum frequency change rate limit K for that region. Hi_l For high inertia H i_h The region, according to formula (39) for ①, 0 + Time and t of ② max The values ​​at different times are compared, and the value with the larger absolute value in equation (39) is multiplied by 1.1-1.2 times the sensitivity to obtain the maximum frequency change rate limit K for that region. Hi_h ,Right now

Images