Power system frequency dynamic characteristic analysis method considering photovoltaic frequency modulation effect

Abstract

The invention relates to the field of power system frequency stability analysis and control, and discloses a power system frequency dynamic characteristic analysis method considering a photovoltaic frequency modulation effect, and the method comprises the steps: building a power system linearization model considering a photovoltaic primary frequency modulation effect; substituting the power system linearization model into an input linear system equation, performing approximate decoupling on the linear system equation to obtain a cubic equation, and deriving an approximate solution of the power system linearization model about a dominant characteristic value of a frequency time domain by solving the cubic equation; and setting a frequency dynamic index, deducing an analytical expression of the frequency dynamic index based on the approximate solution of the dominant characteristic value, and analyzing the influence of the parameters of the power system linearization model on the power system frequency dynamic response characteristics. The method is beneficial to frequency stability analysis and quantitative evaluation of frequency control measures after the power system is connected to photovoltaic power.

Description

Technical Field

[0001] The present invention relates to the field of power system frequency stability analysis and control, and in particular to a method for analyzing the dynamic characteristics of power system frequency taking into account the effect of photovoltaic frequency modulation. Background Art

[0002] Against the backdrop of the "dual carbon" goal, the construction of my country's new power system is progressing steadily. The penetration rate of renewable energy in the new power system continues to increase, and the system's rotational inertia is gradually decreasing. In addition, the output of renewable energy units is volatile, intermittent, and random, which brings safety and stability problems to the power system. Among them, frequency instability is one of the main consequences of the reduction in inertia of the power system. The uncertainty and low inertia characteristics of a high proportion of renewable energy and a high proportion of power electronic equipment make the dynamic mechanism of power system frequency more complex, and new characteristics are reflected in the frequency response. Frequency safety accidents that have occurred at home and abroad in recent years are closely related to this.

[0003] At present, the methods for dynamic analysis of power system frequency are mainly divided into full-time domain simulation method, artificial intelligence-based method and equivalent model method. The full-time domain simulation method takes into account the dynamic characteristics of units, loads, networks and controllers, and uses numerical integration methods to solve the differential algebraic equations of the system, reflecting the frequency dynamic behavior of the system after disturbance as realistically as possible, so as to judge the frequency stability of the system. Although the full-time domain simulation method can describe the power system model in detail and has high calculation accuracy, it is difficult to explain the stability mechanism problem. The frequency prediction method based on artificial intelligence has the characteristics of strong adaptability, small calculation amount and fast calculation speed, but because it is a black box model, it cannot explain the internal mechanism of frequency stability. The equivalent model method simplifies the full-state time domain simulation and assumes that the frequency response of each unit in the system is consistent. The average system frequency model (Average System Frequency, ASF) and the system frequency response model (low order System Frequency Response, SFR) both belong to the equivalent model method.

[0004] Current frequency stability analysis methods often have difficulty in balancing accuracy, computational speed, and interpretability. In addition, the frequency control strategy of new energy units is very different from that of synchronous machines. Traditional simplified models for studying frequency dynamics require further in-depth analysis. Therefore, it is urgent to conduct research on the frequency safety and stability of power systems. Summary of the invention

[0005] The present invention provides a method for analyzing the dynamic characteristics of power system frequency taking into account the photovoltaic frequency modulation effect to solve the above problems.

[0006] The present invention is achieved through the following technical solutions:

[0007] A method for analyzing the dynamic characteristics of power system frequency considering the effect of photovoltaic frequency modulation, comprising:

[0008] Establish a linear model of the power system considering the primary frequency regulation of photovoltaic power;

[0009] Bringing the power system linearization model into a linear system equation with input, and approximately decoupling the linear system equation into a cubic equation, and deriving an approximate solution of the dominant eigenvalue of the power system linearization model in the frequency and time domain by solving the cubic equation;

[0010] A frequency dynamic index is set, and an analytical expression of the frequency dynamic index is derived based on an approximate solution of the dominant eigenvalue, and the influence of the parameters of the power system linearization model on the frequency dynamic response characteristics of the power system is analyzed.

[0011] As an optimization, the specific process of establishing a linearized model of the power system considering the primary frequency regulation of photovoltaic power generation is as follows:

[0012] A1. Obtain the synchronous machine model, prime mover-speed regulator model and photovoltaic power generation system model respectively;

[0013] A2. Linearizing the integrated model consisting of the synchronous machine model, the prime mover-speed regulator model and the photovoltaic power generation system model at the equilibrium point to obtain a sixth-order differential equation;

[0014] A3. Obtain the network equations of the synchronous machine and the photovoltaic power generation system, and based on the network equations, obtain the active power output by the synchronous machine and the photovoltaic power generation system, and then respectively bring the active power of the synchronous machine and the active power of the photovoltaic power generation system into the sixth-order differential equation to obtain the Jacobi matrix of the linearization of the power system, which is the linearization model of the power system.

[0015] As an optimization, in A1, the synchronous machine model is specifically:

[0016]

[0017] Where δ is the generator rotor angle, is the time derivative of δ, ω 0 is the synchronous speed of the generator, ω is the generator speed, is the time derivative of ω, P m is the mechanical power of the prime mover, P e is the electromagnetic power, M is the generator, and D is the damping coefficient;

[0018] The prime mover-speed governor model is specifically:

[0019]

[0020] Where, T CH is the main steam inlet volume and steam chamber time constant, T RH is the reheater time constant, R is the governor droop coefficient, y and z are the prime mover-governor state variables, Respectively represent the time derivatives of y and z;

[0021] The photovoltaic power generation system model is specifically as follows:

[0022]

[0023] Where P ref is the reference value of the active output of the photovoltaic unit, f ref is the system frequency reference value, P pV is the actual active power output of the photovoltaic unit, f is the average frequency of the system, k f is the frequency control gain coefficient of the photovoltaic unit, T f is the frequency measurement time constant, k p , k i is the PI controller parameter, x p is the active state variable of the photovoltaic power generation system, x f is the frequency state variable, V pv is the photovoltaic node voltage, |V pv | is the voltage amplitude of the photovoltaic node, ΔP ref It is the change of the reference value of the active output of the photovoltaic unit.

[0024] As an optimization, in A2, the sixth-order differential equation is expressed as:

[0025]

[0026] Where Δω is the column vector formed by the change in generator speed, ΔP e is the column vector formed by the change of electromagnetic power, Δy and Δz are the column vectors formed by the change of prime mover-speed governor state variables, Δx p , Δx f are column vectors formed by the changes in the active state variables and the frequency state variables of the photovoltaic power generation system, Δ|V pv | is the column vector formed by the change in the voltage amplitude of the photovoltaic node, D is the diagonal matrix formed by the damping coefficient, F HP 、T CH 、T RH is the diagonal matrix formed by the prime mover coefficients, K f 、T f are the diagonal matrices formed by the frequency control gain coefficient and time constant of the photovoltaic unit, R is the diagonal matrix formed by the droop coefficient of the speed regulator, M, M roware the diagonal matrix and row vector formed by the generator inertia time constant, K i , K p are the diagonal matrices formed by the PI controller parameters, |V pv | is the column vector formed by the voltage amplitude of the photovoltaic node, x p is the column vector formed by the active state variables of the photovoltaic power generation system, M i represents the i-th diagonal element formed by the generator inertia time constant, Δδ is the column vector formed by the change in the generator rotor angle, are Δδ, Δω, Δx respectively p , Δx f , Δy, Δz derivatives with respect to time.

[0027] As an optimization, the specific process of A3 is:

[0028] A3.1. Obtain the network equation, which is expressed as:

[0029]

[0030] Where Y i =G i +jB i , is the node admittance matrix, B i is the conductivity matrix, G i is the susceptance matrix, i∈[1,4]; S e , S pv are vectors of complex power of synchronous machine and photovoltaic power generation system respectively, . / represents the division of corresponding elements of the matrix; E is the potential inside the generator, E=E q ′*e jδ , e jδ =cosδ+jsinδ. E q ′ is the vector formed by the transient potential of the synchronous machine q axis, S e , S pv 、V pv , the conjugate of E.

[0031] A3.2. The specific expression of the active power output by the synchronous machine and photovoltaic power generation system is:

[0032]

[0033] Among them, N 1 、N 2 , J 1 , J 2 Both are 0, L 1 The matrix ΔQ formed by the change of photovoltaic reactive power output pv For Δ|V pvThe partial derivative matrix of |.

[0034] Where: E i and E j Respectively represent the internal potential amplitude of the i-th and j-th synchronous machines, V i and V j represents the voltage amplitude of the i-th and j-th photovoltaic cells; δ i and δ j Respectively represent the rotor angles of the i-th and j-th synchronous machines; θ i and θ j Respectively represent the voltage phase angle of the i-th and j-th photovoltaic stations; δ ij ,θ ij Respectively represent δ i -δ j ,θ i -θ j ; The subscript 0 indicates taking the equilibrium point value.

[0035] A3.3. Substitute the active power output by the synchronous machine and the photovoltaic power generation system obtained in step A3.2 into the sixth-order differential equation, and the obtained Jacobian matrix is ​​expressed as:

[0036]

[0037]

[0038] Where I is the identity matrix, ω N is the rated speed of the generator.

[0039] As an optimization, the power system linearization model is brought into the linear system equation with input, and the specific process of approximately decoupling the linear system equation into a cubic equation is as follows:

[0040] B1. Establish a linear system equation with input, and bring the Jacobi matrix into the linear system equation. When the mechanical power of the synchronous machine in the power system undergoes a step change, obtain the speed equation of the synchronous machine in the linear system equation, where the speed equation is formed by multiplying the frequency change of the synchronous machine by the first parameter group;

[0041] B2. Simplify the first parameter group that can be 0 in the speed equation, and finally approximately decouple the first parameter group to obtain a cubic equation.

[0042] As an optimization, the specific process in B1 is:

[0043] B1.1. Substitute the Jacobian matrix into the linear system equation to obtain the following formula:

[0044]

[0045] Among them, Δx is the increment of the state variable, the Jacobian matrix J is the state matrix, B is the input matrix, and r is the input vector;

[0046] The zero-state response solution of the linear system equations is expressed as:

[0047]

[0048] Where t is time, λ i is the eigenvalue of the state matrix J, v i ,q i is i The corresponding right and left eigenvectors; e λ,t λ is the natural constant e i t power;

[0049] And Δx=[Δδ Δω Δx p Δx f Δy Δz] T ;

[0050] B1.2. When the mechanical power of the synchronous machine in the power system undergoes a step change, the speed equation of the synchronous machine is:

[0051] (J 21 J 12 +λJ 22 -λ 2 I+f 1 (λ)J 23 J 42 -d 2 Iλ+f 2 (λ)I)Δω=0

[0052] The first parameter group is (J 21 J 12 +λJ 22 -λ 2 I+f 1 (λ)J 23 J 42 -d 2 Iλ+f 2 (λ)I);

[0053]

[0054] Where λ is the eigenvalue of the state matrix J.

[0055] As an optimization, the specific process of B2 is:

[0056] B2.1, let J 21 J 12 =-M -1 D=-d1 I, J 21 J 12 =-UΛU -1 , where υ, is a column vector, Λ is a matrix J 21 J 12 Eigenvalue γ i The diagonal matrix formed, U is J 21 J 12 The first parameter group is equal to 0 and is simplified to the first equation:

[0057]

[0058] Where d = d 1 +d 2 , J 21 J 12 =-M -1 D=-d 1 I;

[0059] B2.2. Simplify the first equation into the second equation:

[0060]

[0061] in,

[0062] Where n g is the number of generators.

[0063] B2.3. In the time domain solution of frequency, the frequency response is mainly composed of several modes corresponding to eigenvalues ​​(J 21 J 12 The eigenvalues ​​of the state matrix J solved when the eigenvalue is zero are determined by the mode of the eigenvalues), these eigenvalues ​​are called dominant eigenvalues, and the frequency response determined by the dominant eigenvalues ​​is called the dominant mode of the frequency. In the frequency dominant mode, the second equation is approximately converted to the third equation:

[0064]

[0065] B2.4. When the eigenvalue is 0, the third equation is approximately simplified to the fourth equation:

[0066] (a 2 λ 2 +b 2 λ+c 2 )λ[a 4 λ 3 +(da 4 +b 4 )2 +(1+db 4 )λ+b 3 +d]=0

[0067] B2.5. Due to (a 2 λ 2 +b 2 λ+c 2 )≠0, therefore, the fourth equation is simplified to a cubic equation:

[0068]

[0069] As an optimization, the frequency dynamic indicators include the maximum frequency change rate, the frequency steady-state deviation and the maximum frequency deviation.

[0070] As an optimization, the maximum frequency change rate is specifically expressed as:

[0071]

[0072] in, Δω i (0) represents the change in the speed of the i-th generator at time 0, ∑M i represents the sum of the inertia time constants of all generators, M row is the row vector of the diagonal matrix of the generator inertia time constant, ΔP m represents the active power disturbance of the synchronous machine;

[0073] The frequency steady-state deviation is specifically expressed as:

[0074]

[0075] H′ 1 , H′ 3 H 1 , H 3 To eliminate the result of the first row;

[0076] The maximum frequency deviation is specifically expressed as:

[0077]

[0078] At time t NF When the frequency of the power system reaches the maximum deviation,

[0079] α=λ i ,λ i =α+jβ,p r 、p i represents the participation factor p mn The real and imaginary parts of mnis the contribution of the nth eigenvalue to the mth state variable of the system.

[0080] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0081] The present invention discloses a method for analyzing the dynamic characteristics of the frequency of a power system considering the effect of photovoltaic frequency modulation. Based on the analysis method of the full-state model, a linear model of the power system considering the effect of photovoltaic primary frequency modulation is established, and the dynamic response of active power control and primary frequency modulation of the photovoltaic system is accurately described, which can effectively analyze the frequency response characteristics of the power system containing photovoltaics and the power system with obvious frequency time and space distribution. The frequency characteristics of the photovoltaic power system are quantitatively analyzed. Compared with the SFR and ASF model methods, the linear model proposed by the present invention has higher accuracy, which is conducive to the frequency stability analysis and quantitative evaluation of frequency control measures after the power system is connected to photovoltaics. BRIEF DESCRIPTION OF THE DRAWINGS

[0082] The drawings described herein are used to provide a further understanding of the embodiments of the present invention, constitute a part of this application, and do not constitute a limitation of the embodiments of the present invention. In the drawings:

[0083] Figure 1 The accompanying drawing is a schematic flow chart of a method for analyzing the dynamic characteristics of power system frequency taking into account the photovoltaic frequency modulation effect provided by the present invention;

[0084] Figure 2 The accompanying drawing is a wiring diagram of a 10-machine 39-node system provided by an embodiment of the present invention;

[0085] Figure 3 The accompanying drawings are schematic diagrams of time domain simulation frequency and frequency responses of various models provided by an embodiment of the present invention;

[0086] Figure 4 The accompanying drawing is a schematic diagram of frequency error comparison of various models provided by an embodiment of the present invention;

[0087] Figure 5 The accompanying drawing is a schematic diagram of system frequency response under different R provided by an embodiment of the present invention;

[0088] Figure 6 The accompanying drawing is a schematic diagram of system frequency response under different photovoltaic penetration rates provided by an embodiment of the present invention;

[0089] Figure 7 The accompanying drawing is a schematic diagram of system frequency response when photovoltaic power generation is provided with and without frequency control according to an embodiment of the present invention;

[0090] Figure 8 The accompanying drawing is a schematic diagram of system frequency response when the virtual inertia ΔM is increased according to an embodiment of the present invention. DETAILED DESCRIPTION

[0091] In order to make the objectives, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with embodiments and drawings. The exemplary embodiments of the present invention and their description are only used to explain the present invention and are not intended to limit the present invention.

[0092] It should be noted that for the same parameter, bold represents a matrix or vector, and non-bold represents a single element in the matrix or vector.

[0093] This embodiment 1 provides a method for analyzing the dynamic characteristics of power system frequency considering the photovoltaic frequency modulation effect, such as Figure 1-Figure 3 As shown, including:

[0094] S1. Establish a linearized model of the power system considering the primary frequency regulation of photovoltaic power;

[0095] In some embodiments, the specific process of S1 is:

[0096] A1. Obtain the synchronous machine model, prime mover-speed regulator model and photovoltaic power generation system model respectively;

[0097] The synchronous machine model adopts a classic second-order model, specifically:

[0098]

[0099] Where δ is the generator rotor angle, is the time derivative of δ, ω 0 is the synchronous speed of the generator, ω is the generator speed, is the time derivative of ω, P m is the mechanical power of the prime mover, P e is the electromagnetic power, M is the generator, and D is the damping coefficient;

[0100] The prime mover-speed governor model also adopts a second-order model, specifically:

[0101]

[0102] Where, T CH is the main steam inlet volume and steam chamber time constant, T RH is the reheater time constant, R is the governor droop coefficient, y and z are the prime mover-governor state variables, Respectively represent the time derivatives of y and z;

[0103] The photovoltaic power generation system is under normal operating conditions. Therefore, the photovoltaic power generation system model is specifically as follows:

[0104]

[0105] Where P ref is the reference value of the active output of the photovoltaic unit, f ref is the system frequency reference value, P pV is the actual active power output of the photovoltaic unit, f is the average frequency of the system, k f is the frequency control gain coefficient of the photovoltaic unit, T f is the frequency measurement time constant, k p , k i is the PI controller parameter, x p is the active state variable of the photovoltaic power generation system, x f is the frequency state variable, V pv is the photovoltaic node voltage, |V pv | is the voltage amplitude of the photovoltaic node, ΔP ref It is the change of the reference value of the active output of the photovoltaic unit.

[0106] A2. Linearizing the integrated model consisting of the synchronous machine model, the prime mover-speed regulator model and the photovoltaic power generation system model at the equilibrium point to obtain a sixth-order differential equation;

[0107] Considering the classical model of synchronous machine, the second-order model of prime mover-speed regulator and the unified model of photovoltaic system, the following sixth-order differential equation can be obtained by linearization at the equilibrium point:

[0108]

[0109] Where Δω is the column vector formed by the change in generator speed, ΔP e is the column vector formed by the change of electromagnetic power, Δy and Δz are the column vectors formed by the change of prime mover-speed governor state variables, Δx p , Δx f are column vectors formed by the changes in the active state variables and the frequency state variables of the photovoltaic power generation system, Δ|V pv | is the column vector formed by the change in the voltage amplitude of the photovoltaic node, D is the diagonal matrix formed by the damping coefficient, F HP 、T CH 、T RH is the diagonal matrix formed by the prime mover coefficients, K f 、T f are the diagonal matrices formed by the frequency control gain coefficient and time constant of the photovoltaic unit, R is the diagonal matrix formed by the droop coefficient of the speed regulator, M, M row are the diagonal matrix and row vector formed by the generator inertia time constant, K i , K p are the diagonal matrices formed by the PI controller parameters, |V pv | is the diagonal matrix corresponding to the voltage amplitude of the photovoltaic node, xp is the diagonal matrix formed by the active state variables of the photovoltaic power generation system, M i represents the i-th diagonal matrix formed by the generator inertia time constant, Δδ is the column vector formed by the change in the generator rotor angle, are Δδ, Δω, Δx respectively p , Δx f , Δy, Δz derivatives with respect to time.

[0110] A3. Obtain the network equations of the synchronous machine and the photovoltaic power generation system, and based on the network equations, obtain the active power output by the synchronous machine and the photovoltaic power generation system, and then respectively bring the active power of the synchronous machine and the active power of the photovoltaic power generation system into the sixth-order differential equation to obtain the Jacobi matrix of the linearization of the power system, which is the linearization model of the power system.

[0111] The specific process of A3 is as follows:

[0112] A3.1. Obtain the network equation. The system network equation including the photovoltaic unit is:

[0113]

[0114] Where Y i =G i +jB i , is the node admittance matrix, B i is the conductivity matrix, G i is the susceptance matrix, i∈[1,4]; S e , S pv are vectors of complex power of synchronous machine and photovoltaic power generation system respectively, . / represents the division of corresponding elements of the matrix; E is the potential inside the generator, E=E q ′*e jδ , e jδ =cosδ+jsinδ. E q ′ is the vector formed by the transient potential of the synchronous machine q axis, S e , S pv 、V pv , the conjugate of E.

[0115] A3.2. If only the synchronous machine electromagnetic power Pe and the photovoltaic unit active output Ppv are considered, the expression of the synchronous machine and photovoltaic output active power is further obtained as follows: the synchronous machine output active power ΔP e And the active power output of the photovoltaic power generation system ΔP pv The specific expression is:

[0116]

[0117] Since in practice N 1 、N 2 , J 1 , J 2 The value is very small, so the relevant terms are ignored in the formula. Therefore, let N 1 、N 2 , J 1 , J 2 Both are 0, L 1 The matrix ΔQ formed by the change of photovoltaic reactive power output pv For Δ|V pv The partial derivative matrix of |.

[0118] Where: E i and E j Respectively represent the internal potential amplitude of the i-th and j-th synchronous machines, V i and V j represents the voltage amplitude of the i-th and j-th photovoltaic cells; δ i and δ j Respectively represent the rotor angles of the i-th and j-th synchronous machines; θ i and θ j Respectively represent the voltage phase angle of the i-th and j-th photovoltaic stations; δ ij ,θ ij Respectively represent δ i -δ j ,θ i -θ j ; The subscript 0 indicates taking the equilibrium point value.

[0119] A3.3, the active power output by the synchronous machine and the photovoltaic power generation system obtained in step A3.2 is introduced into the sixth-order differential equation (that is, the non-state variable P e and P PV Bringing back the original system state equation, we get the linearized Jacobian matrix J of the system. The Jacobian matrix J is expressed as:

[0120]

[0121] Where I is the identity matrix, ω N is the rated speed of the generator.

[0122] S2, bringing the power system linearization model into a linear system equation with input, and approximately decoupling the linear system equation into a cubic equation, and deriving an approximate solution of the dominant eigenvalue of the power system linearization model in the frequency and time domain by solving the cubic equation;

[0123] In some embodiments, the specific process of S2 is:

[0124] B1. Establish a linear system equation with input, and bring the Jacobi matrix into the linear system equation. When the mechanical power of the synchronous machine in the power system undergoes a step change, obtain the speed equation of the synchronous machine in the linear system equation, where the speed equation is formed by multiplying the frequency change of the synchronous machine by the first parameter group;

[0125] The specific process is:

[0126] B1.1. First, the linear system equation can be expressed as:

[0127]

[0128] Where: x is the state variable, A is the state matrix, B is the input matrix, and r is the input vector.

[0129] Since the state matrix of the present invention is a Jacobian matrix, the Jacobian matrix is ​​introduced into the linear system equation to obtain the following formula:

[0130]

[0131] Among them, Δx is the increment of the state variable, the Jacobian matrix J is the state matrix, B is the input matrix, and r is the input vector;

[0132] The zero-state response solution of the linear system equations is expressed as:

[0133]

[0134] Where t is time, λ i is the eigenvalue of the state matrix J, v i ,q i is i The corresponding right and left eigenvectors; e λ,t λ is the natural constant e i t power;

[0135] And the state variable Δx = [Δδ Δω Δx p Δx f Δy Δz] T ;

[0136] B1.2. When the mechanical power of the synchronous machine in the power system undergoes a step change, the speed equation of the synchronous machine is:

[0137] (J 21 J 12 +λJ 22 -λ 2 I+f 1 (λ)J 23 J 42 -d2 Iλ+f 2 (λ)I)Δω=0

[0138] In order to solve the zero-state response of the power system, the eigenvalue problem of the Jacobi matrix J is analyzed. The eigenvalue and eigenvalue vector of the Jacobi matrix J satisfy the following relationship:

[0139] JΔx=λΔx

[0140] The part about speed (speed equation of synchronous machine) is obtained by sorting out:

[0141] (J 21 J 12 +λJ 22 -λ 2 I+f 1 (λ)J 23 J 42 -d 2 Iλ+f 2 (λ)I)Δω=0

[0142] Where f 1 (λ) and f 2 (λ) is a scalar function, and its expression is:

[0143]

[0144] Where λ is the eigenvalue of the state matrix J.

[0145] Since Δω≠0 in frequency dynamics, the solution of the eigenvalue can be transformed into the solution of the first parameter group (J 21 J 12 +λJ22-λ2I+f1λJ23[42-d2Iλ+f2λI, that is, solve the following column problem:

[0146] det[J 21 J 12 +λJ 22 -λ 2 I+f 1 (λ)J 23 J 42 +f 2 (λ)I]=0.

[0147] That is, B2, simplifying the first parameter group that can be 0 in the speed equation, and finally approximately decoupling the first parameter group to obtain a cubic equation.

[0148] The specific process is:

[0149] B2.1, let J 21 J 12 =-M -1 D=-d1 I, J 21 J 12 =-UΛU -1 , where υ, is a column vector, Λ is the eigenvalue γ i The diagonal matrix formed is then the first parameter group is equal to 0, which simplifies to the first equation:

[0150]

[0151] Where d = d 1 +d 2 ,υ′=U -1 υ, J 21 J 12 =-M -1 D=-d 1 I;

[0152] When the power system damping is uniform, J 22 Let J be a diagonal matrix. 21 J 12 =-M -1 D=-d 1 I; J 23 J 42 is a rank-one matrix, let in is a column vector; J 21 J 12 Non-diagonal matrix, can be diagonalized: J 21 J 12 =-UΛU -1 , where Λ is the eigenvalue γ i The diagonal matrix formed.

[0153] B2.2. Simplify the first equation into the second equation

[0154] Since the above formula in B2.1 can be simplified to:

[0155]

[0156] Where d = d 1 +d 2 ,υ′=U -1 υ,

[0157] Regarding the rank-one corrected determinant, we have the following conclusions:

[0158] det(I+uv T )=1+u T v

[0159] So when det[Λ+λ 2 I+λdI-f 2 When (λ)I]≠0, the above equation can be transformed into the second equation:

[0160]

[0161] in,

[0162] Where n g is the number of generators.

[0163] B2.3. In the time domain solution of frequency, the frequency response is mainly determined by J 21 J 12 The mode of the eigenvalue of the state matrix J solved when the eigenvalue is zero is determined. These eigenvalues ​​are called dominant eigenvalues. The frequency response determined by the dominant eigenvalue is called the dominant mode of the frequency. Under the frequency dominant mode, the second equation is approximately converted into a third equation.

[0164] If there is a dominant mode, an The value of is much larger than the other terms, then the formula can be approximated as a scalar equation:

[0165]

[0166] B2.4. When the eigenvalue is 0, the third equation is approximately simplified to the fourth equation.

[0167] Because J 21 J 12 There must be a zero eigenvalue, so there is γ i = 0. And because is small, and the relevant terms can be ignored. i =0, the above equation can be approximately simplified to the fourth equation:

[0168] (a 2 λ 2 +b 2 λ+c 2 )λ[a 4 λ 3 +(da 4 +b 4 ) 2 +(1+db 4 )λ+b 3 +d]=0

[0169] B2.5. Due to (a 2 λ 2 +b 2 λ+c 2 )≠0, therefore, the fourth equation is simplified to a cubic equation.

[0170] It can be seen that the eigenvalue formula is approximately decoupled into a quadratic equation and a cubic equation problem, where (a 2 λ 2 +b 2 λ+c 2 )≠0. So we only need to focus on the solution of the cubic equation. The solution of the cubic equation can be obtained using the Shengjin formula, and its discriminant is:

[0171]

[0172] When the discriminant Δ>0, the equation has one real root and a pair of conjugate complex roots; when Δ<0, the equation has three unequal real roots. When Δ=0, the equation has three real roots, one of which is a double root. The solution to this cubic equation is the characteristic root of the frequency-dominant mode, that is, the dominant eigenvalue. In actual examples, the solution to the dominant eigenvalue is generally a real number and a pair of conjugate complex roots.

[0173] S3, set the frequency dynamic index, and derive the analytical expression of the frequency dynamic index based on the approximate solution of the dominant eigenvalue, and analyze the influence of the parameters of the power system linearization model on the dynamic response characteristics of the power system frequency. Based on the S2 linearization model and its time domain solution, analyze the system frequency response characteristics under the action of photovoltaic frequency modulation.

[0174] In some embodiments, the frequency dynamic indicators include a maximum frequency change rate, a frequency steady-state deviation, and a maximum frequency deviation.

[0175] Maximum frequency change rate:

[0176] When the synchronous machine mechanical power step change occurs in the power system, the power system inertia center frequency (COI) deviation expression is:

[0177]

[0178] Substitute Δω i (t) time domain expression, we can get:

[0179]

[0180] In the formula, E λ for The diagonal matrix formed by Δω i (t) represents the change in the speed of the i-th generator at time t, ∑M i represents the sum of the inertia time constants of all generators, M row is the row vector of the diagonal matrix of the generator inertia time constant, ΔP m represents the active power disturbance of the synchronous machine;

[0181] Pω The right eigenvalue vector v i The part of the matrix P related to the angular frequency; Q is the left eigenvector q i The matrix formed.

[0182] When t=0, the rate of change of the system inertia center frequency at time 0, that is, the maximum frequency change rate, can be obtained:

[0183]

[0184] It can be seen from the expression that the maximum change rate of the system inertia center frequency is only related to the total inertia of the synchronous machine and the size of the active disturbance, and has nothing to do with the speed regulator and photovoltaic frequency control parameters. It is inversely proportional to the total inertia of the synchronous machine and directly proportional to the active disturbance of the synchronous machine.

[0185] Frequency steady-state deviation:

[0186] The steady-state deviation expression of the system state quantity can be expressed as:

[0187] Δx(∞)=-PΛ′ -1 QBr=-J′ -1 Br

[0188] Where: J′ is the reduced-order matrix of J with zero eigenvalues ​​eliminated, Λ′ -1 is a diagonal matrix composed of eigenvalues. Since there must be a zero eigenvalue in J, the zero eigenvalue can be eliminated by setting the power angle of a synchronous machine as the reference. At this time, the Jacobian matrix J becomes J′. Compared with the Jacobian matrix J before reduction, the sub-block J 12 , J 21 becomes J′ 12 , J′ 21 , whose expression is:

[0189]

[0190] Where: H′ 1 , H′ 3 H 1 , H 3 To eliminate the result of the first row, let J -1 ′ corresponds to J in J′ 22 The sub-block of position is K 22 , take the part related to the synchronous machine speed ω:

[0191] Δω(∞)=-K 22 M -1 ΔP m

[0192] Pair Block K 22 Perform analysis. Block inversion:

[0193]

[0194] Get Π 1 The expression is:

[0195]

[0196] The sub-blocks Corresponding to ΔP m The part related to the speed ω, i.e. Π 11 ,Π 12 , Π 13 The lower right corner of the block is K 22 .

[0197] right Under small disturbance, it is approximately considered that the photovoltaic node voltage V pv The amplitude is constant, Substituting in, we get

[0198]

[0199] Among them J 1 ' 2 and J 2 ' 1 The expression of does not contain the damping term D, is a fully negative rank-one matrix. Therefore, it can be seen that reducing R and increasing K f Both are equivalent to increasing the system damping D, thereby reducing the system frequency steady-state error after disturbance.

[0200] Further calculations can give the analytical expression for the frequency steady-state deviation:

[0201]

[0202] Where:

[0203]

[0204] It can be seen that the steady-state frequency deviation has nothing to do with the inertia M, that is, the inertia M of the synchronous machine does not affect the steady-state frequency deviation.

[0205] Maximum frequency deviation:

[0206] The dominant eigenvalue is used to calculate the system frequency. The system frequency in the dominant mode is approximately the system inertia center frequency. The time domain analytical expression for Δω is:

[0207]

[0208] Where Δω jiis the i-th component of the speed deviation of the j-th synchronous machine, p r and p i are the participation factors p ji The real and imaginary parts, p mn Quantifies the contribution of the nth eigenvalue to the mth state variable of the system. mn The expression is p mn =v mn q n Br, q n represents the nth left eigenvector, where v mn is an element in the right eigenvector. By taking the derivative of the Δω expression with respect to time and setting it equal to 0, we can obtain the time t at which the frequency reaches the maximum deviation. NF :

[0209]

[0210] We can further get the maximum frequency deviation, that is, the expression of the lowest frequency point of the system:

[0211]

[0212] In the embodiment, the Figure 2 The 10-machine 39-node classic example is shown in Table 1. Some parameters of the synchronous machine are shown in Table 1. The remaining parameters are shown in Table 1. CH =0.3s, T RH =7s, F HP =0.3. Photovoltaic control parameter T f =0.02s, k p =0.1, k i =2.

[0213] Table 1 Partial parameters of synchronous machine

[0214]

[0215] (1) Model accuracy verification

[0216] Replace the synchronous machines at nodes 32 and 33 with photovoltaics. Reduce the active output of the synchronous machine by 100MW at t=0. Draw the system inertia center frequency under time domain simulation and compare it with the system frequency response under the linear model proposed in this invention, the ASF model, and the SFR model, as shown in Figure 2. Figure 3 The error of the linear model frequency analytical solution, the error of the linear model dominant mode frequency, and the errors of the ASF model and SFR model are shown in Figure 4 shown.

[0217] The goodness of fit of each model is shown in Table 2.

[0218] Table 2 Fit of each model

[0219]

[0220] On the one hand, the SFR model of equal value aggregation cannot consider the dynamics of all units separately because it aggregates all synchronous machines and speed regulators into one machine. The ASF model aggregates all synchronous machines into one machine, but retains the dynamics of all speed regulators, so its model accuracy is slightly better than SFR. The linear model proposed in this invention retains the dynamics of all synchronous machines and speed regulators, so its accuracy is improved compared with the SFR model and the ASF model, and it is suitable for systems with obvious time and space distribution. On the other hand, both the SFR model and the ASF model ignore the network and assume that the electromagnetic power P e The linear model proposed in the present invention takes into account the influence of the network, so the accuracy is improved. In addition, the processing of new energy in the SFR model and the ASF model (such as power closed loop and reference frequency closed loop) is simpler than that in the present invention. The frequency modulation of the photovoltaic power generation system of the present invention adopts a second-order model. In processing systems containing new energy, the present invention has better accuracy.

[0221] (2) Influence of synchronous machine parameters on power system frequency response

[0222] When R=∞, the speed regulator does not participate in frequency modulation. Changing the speed regulator coefficient R, the system frequency response is compared. Figure 5 As shown. Figure 5 It can be seen that the change of R does not change the frequency change rate at time 0. Reducing R is equivalent to increasing the system damping, thereby reducing the steady-state frequency deviation of the system, which is consistent with the results of theoretical analysis. At the same time, reducing R will lead to an increase in the real and imaginary parts of the dominant eigenvalues, thereby reducing the time domain solution of the frequency deviation. Therefore, when R decreases, the lowest point of the system frequency and the steady-state frequency value will improve.

[0223] (3) Impact of PV penetration and frequency regulation capability on power system frequency response

[0224] The system loses 10% active power, and the PV unit does not participate in the system frequency regulation, that is, k f = 0, gradually change the system photovoltaic penetration rate, the system frequency response is as follows Figure 6 It is found that when the penetration rate of photovoltaic power generation that does not participate in frequency regulation increases, the frequency stability problem of the power system will worsen.

[0225] Normally, when the system active power changes by 10%, the system frequency should not change by more than 1Hz. When the PV penetration rate is 30%, the system minimum frequency is lower than 49Hz, which does not meet the frequency safety and stability requirements of the power system. Therefore, in this example, when the PV unit has no frequency control, the PV penetration rate should not exceed 30%.

[0226] When the photovoltaic penetration rate is 30%, the system response curve when the photovoltaic unit frequency control is added is compared with that when the photovoltaic unit frequency control is not added. Figure 7 Obviously, after the photovoltaic unit participates in frequency regulation, the system frequency response characteristics are improved.

[0227] Further analysis shows that k f The change does not affect the frequency change rate at time 0, increasing k f This is equivalent to increasing the system damping D, thereby reducing the steady-state frequency deviation of the system. At the same time, increasing k f The frequency drop depth can be improved, and the system frequency will stabilize faster.

[0228] (4) Impact of virtual inertia on power system frequency response

[0229] Research shows that energy storage can be used as virtual inertia to reduce grid frequency fluctuations and improve the dynamic characteristics of frequency in the case of active power shortage. Assuming that there is enough energy storage in the grid to provide virtual inertia support ΔM, the system frequency response is compared by changing the virtual inertia ΔM. Figure 8 As shown in the figure. The larger the total inertia of the system, the smaller the rate of change of the system inertia center frequency at time 0. After adding the virtual inertia ΔM in the figure, the frequency response curve decreases more slowly. From the theoretical analysis part, it can be seen that the frequency steady-state error has nothing to do with the inertia M. Therefore, when the virtual inertia ΔM is increased, the system frequency will eventually stabilize at the same value.

[0230] The specific implementation methods described above further illustrate the objectives, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above description is only a specific implementation method of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A method for analyzing the dynamic characteristics of power system frequency considering the effect of photovoltaic frequency modulation, characterized in that: include: Establish a linear model of the power system considering the primary frequency regulation of photovoltaic power; Bringing the power system linearization model into a linear system equation with input, and approximately decoupling the linear system equation into a cubic equation, and deriving an approximate solution of the dominant eigenvalue of the power system linearization model in the frequency and time domain by solving the cubic equation; A frequency dynamic index is set, and an analytical expression of the frequency dynamic index is derived based on an approximate solution of the dominant eigenvalue, and the influence of the parameters of the power system linearization model on the frequency dynamic response characteristics of the power system is analyzed.

2. A method for analyzing the dynamic characteristics of power system frequency considering photovoltaic frequency modulation according to claim 1, characterized in that: The specific process of establishing a linearized model of the power system considering the primary frequency regulation of photovoltaic power generation is as follows: A1. Obtain the synchronous machine model, prime mover-speed regulator model and photovoltaic power generation system model respectively; A2. Linearizing the integrated model consisting of the synchronous machine model, the prime mover-speed regulator model and the photovoltaic power generation system model at the equilibrium point to obtain a sixth-order differential equation; A3. Obtain the network equations of the synchronous machine and the photovoltaic power generation system, and based on the network equations, obtain the active power output by the synchronous machine and the photovoltaic power generation system, and then respectively bring the active power of the synchronous machine and the active power of the photovoltaic power generation system into the sixth-order differential equation to obtain the Jacobi matrix of the linearization of the power system, which is the linearization model of the power system.

3. A method for analyzing the dynamic characteristics of power system frequency considering photovoltaic frequency modulation according to claim 2, characterized in that: In A1, the synchronous machine model is specifically: Where δ is the generator rotor angle, is the time derivative of δ, ω0 is the synchronous speed of the generator, ω is the generator speed, is the time derivative of ω, P m is the mechanical power of the prime mover, P e is the electromagnetic power, M is the generator, and D is the damping coefficient; The prime mover-speed governor model is specifically: Where, T CH is the main steam inlet volume and steam chamber time constant, T RH is the reheater time constant, R is the governor droop coefficient, y and z are the prime mover-governor state variables, Respectively represent the time derivatives of y and z; The photovoltaic power generation system model is specifically as follows: Where P ref is the reference value of the active output of the photovoltaic unit, f ref is the system frequency reference value, P pV is the actual active power output of the photovoltaic unit, f is the average frequency of the system, k f is the frequency control gain coefficient of the photovoltaic unit, T f is the frequency measurement time constant, k p , k i is the PI controller parameter, x p is the active state variable of the photovoltaic power generation system, x f is the frequency state variable, V pv is the photovoltaic node voltage, |V pv | is the voltage amplitude of the photovoltaic node, △P ref It is the change of the reference value of the active output of the photovoltaic unit.

4. A method for analyzing the dynamic characteristics of power system frequency considering photovoltaic frequency modulation according to claim 3, characterized in that: In A2, the sixth-order differential equation is expressed as: Where Δω is the column vector formed by the change in generator speed, ΔP e is the column vector formed by the change of electromagnetic power, Δy and Δz are the column vectors formed by the change of prime mover-speed governor state variables, Δx p , Δx f are column vectors formed by the changes in the active state variables and the frequency state variables of the photovoltaic power generation system, Δ|V pv | is the column vector formed by the change in the voltage amplitude of the photovoltaic node, D is the diagonal matrix formed by the damping coefficient, F HP , T CH , T RH is the diagonal matrix formed by the prime mover coefficients, K f , T f are the diagonal matrices formed by the frequency control gain coefficient and time constant of the photovoltaic unit, R is the diagonal matrix formed by the droop coefficient of the speed regulator, M, M row are the diagonal matrix and row vector formed by the generator inertia time constant, K i , K p are the diagonal matrices formed by the PI controller parameters, |V pv | is the column vector formed by the voltage amplitude of the photovoltaic node, x p is the column vector formed by the active state variables of the photovoltaic power generation system, M i represents the i-th diagonal element formed by the generator inertia time constant, Δδ is the column vector formed by the change in the generator rotor angle, are Δδ, Δω, Δx respectively p , Δx f , Δy, Δz derivatives with respect to time.

5. A method for analyzing the dynamic characteristics of power system frequency considering photovoltaic frequency modulation according to claim 4, characterized in that: The specific process of A3 is as follows: A3.

1. Obtain the network equation, which is expressed as: Where Y i =G i +jB i , is the node admittance matrix, B i is the conductivity matrix, G i is the susceptance matrix, i∈[1,4]; S e , S pv are vectors of complex power of synchronous machine and photovoltaic power generation system respectively, . / represents the division of corresponding elements of matrix; E is the potential inside the generator, E=E q ′*e jδ , e jδ =cosδ+jsinδ,E q ′ is the vector formed by the transient potential of the synchronous machine q axis, S e , S pv 、V pv , the conjugate of E. A3.

2. The specific expression of the active power output by the synchronous machine and photovoltaic power generation system is: H 2,ij =E i0 V j0 (G 2,ij sin(δ i0 -θ j0 )-B 2,ij cos(δ i0 -θ j0 )) H 3,ij =V i0 E j0 (G 3,ij sin(θ i0 -d j0 )-B 3,ij cos(θ i0 -d j0 )) Among them, N1, N2, J1, and J2 are all 0, and L1 is the change in photovoltaic reactive power output ΔQ pv For Δ|V pv |The partial derivative matrix formed. Where: E i and E j Respectively represent the internal potential amplitude of the i-th and j-th synchronous machines, V i and V j represents the voltage amplitude of the i-th and j-th photovoltaic cells; δ i and δ j Respectively represent the rotor angles of the i-th and j-th synchronous machines; θ i and θ j Respectively represent the voltage phase angle of the i-th and j-th photovoltaic stations; δ ij ,θ ij Respectively represent δ i -δ j ,θ i -θ j ; The subscript 0 indicates taking the equilibrium point value. A3.

3. Substitute the active power output by the synchronous machine and the photovoltaic power generation system obtained in step A3.2 into the sixth-order differential equation, and the obtained Jacobian matrix is ​​expressed as: Where I is the identity matrix, ω N is the rated speed of the generator.

6. A method for analyzing the dynamic characteristics of power system frequency considering photovoltaic frequency modulation according to claim 5, characterized in that: The specific process of bringing the power system linearization model into the linear system equation with input and approximating and decoupling the linear system equation into a cubic equation is as follows: B1. Establish a linear system equation with input, and bring the Jacobi matrix into the linear system equation. When the mechanical power of the synchronous machine in the power system undergoes a step change, obtain the speed equation of the synchronous machine in the linear system equation, where the speed equation is formed by multiplying the frequency change of the synchronous machine by the first parameter group. B2. Simplify the first parameter group that can be 0 in the speed equation, and finally approximately decouple the first parameter group to obtain a cubic equation.

7. A method for analyzing the dynamic characteristics of power system frequency considering photovoltaic frequency modulation according to claim 6, characterized in that: The specific process in B1 is: B1.

1. Substitute the Jacobian matrix into the linear system equation to obtain the following formula: Among them, Δx is the increment of the state variable, the Jacobian matrix J is the state matrix, B is the input matrix, and r is the input vector; The zero-state response solution of the linear system equations is expressed as: Where t is time, λ i is the eigenvalue of the state matrix J, v i ,q i is i The corresponding right and left eigenvectors, e λ,t λ is the natural constant e i t power; And △x=[△δ △ω △x p △x f △y △z] T ; B1.

2. When the mechanical power of the synchronous machine in the power system undergoes a step change, the speed equation of the synchronous machine is: (J 21 J 12 +λJ 22 -l 2 I+f1(λ)J 23 J 42 -d2Iλ+f2(λ)I)△ω=0 The first parameter group is (J 21 J 12 +λJ 22 -λ 2 I+f1(λ)J 23 J 42 -d2Iλ+f2(λ)I); Where λ is the eigenvalue of the state matrix J.

8. A method for analyzing the dynamic characteristics of power system frequency considering photovoltaic frequency modulation according to claim 7, characterized in that: The specific process of B2 is: B2.1, let J 21 J 12 =-M -1 D=-d1I, J 21 J 12 =-UΛU -1 , where υ, is a column vector, Λ is a matrix J 21 J 12 Eigenvalue γ i The diagonal matrix formed, U is J 21 J 12 The non-singular transformation matrix of the diagonal matrix Λ is transformed, and the first parameter group is equal to 0, which is simplified to the first equation: Where, d = d1 + d2, J 21 J 12 =-M -1 D = -d1I; B2.

2. Simplify the first equation into the second equation: in, Where n g is the number of generators; B2.

3. In the time domain solution of frequency, the frequency response is mainly determined by J 21 J 12 The mode of the eigenvalue of the state matrix J solved when the eigenvalue is zero is determined. These eigenvalues ​​are called dominant eigenvalues. The frequency response determined by the dominant eigenvalue is called the dominant mode of the frequency. In the frequency dominant mode, the second equation is approximately converted to the third equation: B2.4, when the J 21 J 12 When the eigenvalue is 0, the third equation is approximately simplified to the fourth equation: (a2λ 2 +b2λ+c2)λ[a4λ 3 +(da4+b4)λ 2 +(1+db4)λ+b3+d]=0 B2.5, due to (a2λ 2 +b2λ+c2)≠0, therefore, the fourth equation is simplified to a cubic equation:

9. A method for analyzing the dynamic characteristics of power system frequency considering photovoltaic frequency modulation according to claim 8, characterized in that: The frequency dynamic indicators include the maximum frequency change rate, the frequency steady-state deviation and the maximum frequency deviation.

10. A method for analyzing the dynamic characteristics of power system frequency considering photovoltaic frequency modulation according to claim 9, characterized in that: The maximum frequency change rate is specifically expressed as: in, Δω i (0) represents the change in the speed of the i-th generator at time 0, ∑M i represents the sum of the inertia time constants of all generators, M row is the row vector of the diagonal matrix of the generator inertia time constant, ΔP m represents the active power disturbance of the synchronous machine; The frequency steady-state deviation is specifically expressed as: H′1 and H′3 are the results of eliminating the first row of H1 and H3 respectively; The maximum frequency deviation is specifically expressed as: At time t NF When the frequency of the power system reaches the maximum deviation, α=λ i ,λ i =α+jβ,p r 、p i represents the participation factor p mn The real and imaginary parts of mn is the contribution of the nth eigenvalue to the mth state variable of the system.

Images