Parameter planning method for power system including synchronous generators and virtual synchronous machines

Abstract

The present invention discloses a method for parameter planning of a power system including a synchronous generator and a virtual synchronous machine, comprising the following steps: S1: constructing a differential equation mathematical model of a synchronous generator; S2: constructing a differential equation mathematical model of a virtual synchronous machine; S3: connecting the synchronous generator and the virtual synchronous machine to the system, and constructing a differential equation mathematical model of a power system including the synchronous generator and the virtual synchronous machine; S4: designing the optimal inertia and optimal damping of the synchronous generator and the virtual synchronous machine by using an arithmetic algorithm based on the minimum system inertia and damping support requirements. Under the premise of ensuring the minimum system inertia and damping support requirements, the present invention obtains the optimal inertia and optimal damping configuration of each unit through an arithmetic algorithm, reduces the negative impact caused by parameter coupling between units, and improves the dynamic performance of the system.

Description

Technical Field

[0001] The invention relates to the field of novel power systems, and in particular to a parameter planning method for a power system containing synchronous generators and virtual synchronous machines. Background Art

[0002] With the development of new energy power generation technology at home and abroad, the proportion of virtual synchronous machines in the power system continues to increase, resulting in the system having low inertia and low damping characteristics, which makes its robustness worse, especially the frequency support capability is significantly reduced, and also poses an important challenge to the power supply planning of the entire power system. In order to rationally plan the parameter configuration of synchronous generators and virtual synchronous machines in the power system and improve the dynamic performance of the power system, the present invention proposes a power system parameter planning method containing synchronous generators and virtual synchronous machines. The method is based on the minimum system inertia and damping support requirements, and plans the optimal inertia and optimal damping configuration of each unit through an arithmetic algorithm, reduces the negative impact of parameter coupling between units, ensures the frequency stability of the power system, and improves the dynamic performance of the system. Summary of the invention

[0003] To achieve the above purpose, the technical solution provided by the present invention is:

[0004] S1: Construct the differential equation mathematical model of synchronous generator;

[0005] S2: Construct the differential equation mathematical model of the virtual synchronous machine;

[0006] S3: Connect the synchronous generator and the virtual synchronous machine to the system and construct a mathematical model of differential equations of the power system including the synchronous generator and the virtual synchronous machine;

[0007] S4: Based on the minimum system inertia and damping support requirements, the optimal inertia and optimal damping of the synchronous generator and virtual synchronous machine are designed using an arithmetic algorithm;

[0008] The specific steps of step S1 are as follows:

[0009] S1-1: Get the initial information of the synchronous generator, including: the transient reactance X of the synchronous generator on the d-axis and q-axis d ′、X q ′, the subtransient reactance X of the synchronous generator on the d-axis d ″, the transient resistance R of the synchronous generator i ′, the rated angular frequency of the synchronous generator ω tr , the initial admittance of the synchronous generator Y i The real part G i and the imaginary part B i ;

[0010] S1-2: Connect the synchronous generator to the AC grid and update the parameters, including the corrected admittance Y at the node where the synchronous generator is located i ′, virtual injection current I at the node where the synchronous generator is located i ′, where the modified admittance Y at the node where the synchronous generator is located is i 'for:

[0011]

[0012] In formula (1), δ is the power angle of the synchronous generator;

[0013] Virtual injection current I at the node where the synchronous generator is located i ′ can be expressed as the real part I xi ′ and the imaginary part I yi ′, that is, I i ′=I xi ′+jI yi ′, the specific expression is:

[0014]

[0015] In formula (2), E qi ′ is the transient electromotive force of the synchronous generator on the q axis;

[0016] S1-3: Based on the modified admittance Y at the node where the synchronous generator is located i ′ and virtual injection current I i ′, use Gaussian elimination method to directly solve the above linear equations (2) to obtain the actual voltage V at the node where the synchronous generator is located t The real part V x and the imaginary part V y , and then according to the virtual injection current I at the node where the synchronous generator is located i ′、initial admittance Y i and the actual voltage V t , and the injected current I at the node where the synchronous generator is located is obtained t The real part I x and the imaginary part I y , and then get the actual active power P of the node where the synchronous generator is located t =V x I x ;

[0017] S1-4: According to the rotor motion equation of the synchronous generator, the differential equation of the synchronous generator is obtained:

[0018]

[0019] In formula (3), ω is the angular frequency of the synchronous generator, P tris the output active power reference value of the synchronous generator, D t is the damping coefficient of the synchronous generator, T t is the inertia coefficient of the synchronous generator, E sv is the stator voltage of the excitation system in the synchronous generator, E ev is the excitation voltage of the excitation system in the synchronous generator, I d is the current of the synchronous generator on the d-axis, T stp is the subtransient time constant of the excitation system in the synchronous generator, E rev is the rated excitation voltage of the excitation system in the synchronous generator, G t is the proportional gain parameter of the excitation system in the synchronous generator, V tr is the reference voltage at the node where the synchronous generator is located, Amp is the amplitude symbol, T tc is the time constant of the excitation system in the synchronous generator;

[0020] The specific steps of step S2 are as follows:

[0021] S2-1: Get the initial information of the virtual synchronous machine, including: the line equivalent reactance X of the virtual synchronous machine f , the output reactive power reference value Q at the node where the virtual synchronous machine is located xr , the initial admittance of the virtual synchronous machine Y j The real part G j and the imaginary part B j ;

[0022] S2-2: Connect the virtual synchronous machine to the AC grid and update the parameters, including: the modified admittance Y at the node where the virtual synchronous machine is located j ′, the virtual injection current I at the node where the virtual synchronous machine is located j ′, where the modified admittance Y at the node where the virtual synchronous machine is located is j 'for:

[0023]

[0024] Virtual injection current I at the node where the virtual synchronous machine is located j ′=I xj ′+jI yj ′, where I xj ′ and I yj ′ are virtual injection currents I j The real and imaginary parts of ′ are expressed as follows:

[0025]

[0026] In formula (5), ρ is the power angle of the virtual synchronous machine, E d is the electromotive force of the virtual synchronous machine on the d-axis;

[0027] S2-3: Based on the modified admittance Y at the node where the virtual synchronous machine is located j ′ and virtual injection current I j ′, use Gaussian elimination method to directly solve the above linear equations (5) and get the actual voltage V at the node where the virtual synchronous machine is located t The real part of ′ V′ x and the imaginary part V′ y , and then according to the virtual injection current I at the node where the virtual synchronous machine is located j ′、initial admittance Y j and the actual voltage V t ′, and the injected current I at the node where the virtual synchronous machine is located is obtained t The real part of ′I′ x and the imaginary part I′ y , and then get the actual active power P′ of the node where the virtual synchronous machine is located t =V′ x I′ x and actual reactive power Q′ t =V′ y I′ y ;

[0028] S2-4: According to the rotor motion equation of the virtual synchronous machine, the differential equation of the virtual synchronous machine is obtained:

[0029]

[0030] In formula (6), σ is the power angle of the virtual synchronous machine, ω x is the angular frequency of the virtual synchronous machine, ω xr is the rated angular frequency of the virtual synchronous machine, P xr is the output active power reference value at the node where the virtual synchronous machine is located, D x is the damping coefficient of the virtual synchronous machine, J x is the inertia coefficient of the virtual synchronous machine;

[0031] The specific steps of step S3 are as follows:

[0032] S3-1: according to the initial information of the synchronous generator and equation (3) in step S1-4, the differential equation is solved by using the improved Euler method to obtain a mathematical model of the differential equation of the synchronous generator connected to the power system;

[0033] S3-2: according to the initial information of the virtual synchronous machine and equation (6) in step S2-4, the differential equation is solved by using the improved Euler method to obtain a differential equation mathematical model of the virtual synchronous machine connected to the power system;

[0034] S3-3: Combining the differential equation mathematical models of step S3-1 and step S3-2, constructing a differential equation mathematical model of a power system including a synchronous generator and a virtual synchronous machine;

[0035] The specific steps of step S4 are as follows:

[0036] S4-1: Considering the minimum system inertia and damping support requirements, the objective function Obj is set as:

[0037]

[0038] In formula (7), M is the number of synchronous generators, and N is the number of virtual synchronous machines;

[0039] S4-2: Set the penalty function Pen to:

[0040] Pen=p·(|J total -N.J. x -M.T t |+|D total -N·D x -M.D. t |) (8)

[0041] In formula (8), p is the penalty coefficient, J total is the total inertia of the system, D total is the total damping of the system;

[0042] S4-3: Combine steps S4-1 and S4-2 to design the fitness function of the arithmetic algorithm as:

[0043] fitness=Obj+Pen (9)

[0044] S4-4: Set the damping D of the synchronous generator t , Inertia of synchronous generator T t , the damping D of the virtual synchronous machine x , inertia J of virtual synchronous machine x As the quantity to be optimized and generate the initial population;

[0045] S4-5: Set the mathematical optimization accelerator. The specific expression is:

[0046]

[0047] In formula (10), n is the current iteration number, n max is the maximum number of iterations, M min and M max Define a random number r for the minimum / maximum value of the corresponding function value of the mathematical optimization accelerator 1 , when r 1When it is <MOA>, enter the exploration stage in step S4-6; otherwise, enter the development stage in step S4-7.

[0048] S4-6: The specific expression for the exploration stage is:

[0049]

[0050] In Equation (11), r 2 is a random number uniformly distributed over [0, 1], x ij (n + 1) is the j-th position of the i-th solution in the next iteration, best(x j ) is the j-th position of the optimal solution in the iteration process so far, λ is a small constant 0.01 to prevent the denominator from being zero, OM j and ON j respectively represent the upper bound value and the lower bound value of the optimal value at the j-th position to prevent the solution from going out of bounds after updating the individual, ξ is the control parameter for adjusting the exploration process, MOC is the mathematical optimization coefficient, and its calculation formula is:

[0051]

[0052] In Equation (12), θ is the sensitive parameter;

[0053] S4-7: The specific expression for the development stage is:

[0054]

[0055] In Equation (13), r 3 is a random number uniformly distributed over [0, 1];

[0056] S4-8: Substitute the individual into Equation (9), calculate the fitness value, find the set of values g best with the smallest fitness value at the current iteration number, compare the individual best values in all iterations, select the set of values with the smallest fitness value as the global best value f best , and determine whether the maximum iteration number is reached at this time. If it is reached, jump to step S4-9; if not, jump to step S4-5;

[0057] S4-9: End the arithmetic algorithm, consider the minimum system inertia and damping support requirements, use the arithmetic algorithm to design the optimal inertia and optimal damping of the synchronous generator and the virtual synchronous machine, and substitute them into the differential equation mathematical model of the power system to optimize the dynamic characteristics of the power system.

[0058] Compared with the prior art, the principle and advantages of this solution are as follows:

[0059] The present invention discloses a parameter planning method for a power system including a synchronous generator and a virtual synchronous machine. First, a differential equation mathematical model of a synchronous generator is constructed; second, a differential equation mathematical model of a virtual synchronous machine is constructed; then, the synchronous generator and the virtual synchronous machine are connected to the system to construct a differential equation mathematical model of a power system including the synchronous generator and the virtual synchronous machine; finally, based on the minimum system inertia and damping support requirements, an arithmetic algorithm is used to design the optimal inertia and optimal damping of the synchronous generator and the virtual synchronous machine. Under the premise of ensuring the minimum system inertia and damping support requirements, the present invention obtains the optimal inertia and optimal damping configuration of each unit through an arithmetic algorithm, reduces the negative impact caused by parameter coupling between units, and improves the dynamic performance of the system. BRIEF DESCRIPTION OF THE DRAWINGS

[0060] Figure 1 It is a flow chart of a method for planning parameters of a power system including synchronous generators and virtual synchronous machines in an embodiment of the present invention;

[0061] Figure 2 It is a schematic diagram of two synchronous generators and one virtual synchronous machine connected to a power system in an embodiment of the present invention;

[0062] Figure 3 is a curve showing the change of the fitness value of the system after the arithmetic algorithm is adopted in the embodiment of the present invention;

[0063] Figure 4 is an inertia variation curve of two synchronous generators and one virtual synchronous machine after adopting an arithmetic algorithm in an embodiment of the present invention;

[0064] Figure 5 is a damping variation curve of two synchronous generators and one virtual synchronous machine after adopting an arithmetic algorithm in an embodiment of the present invention;

[0065] Figure 6 The angular dynamic frequency curve before and after the planning of the first synchronous generator in the embodiment of the present invention;

[0066] Figure 7 The angular dynamic frequency curve before and after planning of the second synchronous generator in the embodiment of the present invention;

[0067] Figure 8 It is the angular dynamic frequency curve before and after the virtual synchronous machine planning in the embodiment of the present invention. DETAILED DESCRIPTION

[0068] The present invention will be further described below in conjunction with specific embodiments:

[0069] Figure 1 The flowchart of the power system parameter planning method including synchronous generators and virtual synchronous machines is shown. Figure 2The figure shows a schematic diagram of two synchronous generators and one virtual synchronous machine connected to a power system. The parameter planning method includes the following steps:

[0070] S1: Construct the differential equation mathematical model of synchronous generator;

[0071] S2: Construct the differential equation mathematical model of the virtual synchronous machine;

[0072] S3: Connect the synchronous generator and the virtual synchronous machine to the system and construct a differential equation mathematical model of the power system including the synchronous generator and the virtual synchronous machine;

[0073] S4: Based on the minimum system inertia and damping support requirements, the optimal inertia and optimal damping of the synchronous generator and virtual synchronous machine are designed using an arithmetic algorithm;

[0074] The specific steps of step S1 are as follows:

[0075] S1-1: Get the initial information of the synchronous generator, including: the transient reactance X of the synchronous generator on the d-axis and q-axis d ′、X q ′, the subtransient reactance X of the synchronous generator on the d-axis d ″, the transient resistance R of the synchronous generator i ′, the rated angular frequency of the synchronous generator ω tr , the initial admittance of the synchronous generator Y i The real part G i and the imaginary part B i ;

[0076] S1-2: Connect the synchronous generator to the AC grid and update the parameters, including the corrected admittance Y at the node where the synchronous generator is located i ′, virtual injection current I at the node where the synchronous generator is located i ′, where the modified admittance Y at the node where the synchronous generator is located is i 'for:

[0077]

[0078] In formula (14), δ is the power angle of the synchronous generator;

[0079] Virtual injection current I at the node where the synchronous generator is located i ′ can be expressed as the real part I xi ′ and the imaginary part I yi ′, that is, I i ′=I xi ′+jI yi ′, the specific expression is:

[0080]

[0081] In formula (15), E qi ′ is the transient electromotive force of the synchronous generator on the q axis;

[0082] S1-3: Based on the modified admittance Y at the node where the synchronous generator is located i ′ and virtual injection current I i ′, use Gaussian elimination method to directly solve the above linear equations (15) and get the actual voltage V at the node where the synchronous generator is located t The real part V x and the imaginary part V y , and then according to the virtual injection current I at the node where the synchronous generator is located i ′、initial admittance Y i and the actual voltage V t , and the injected current I at the node where the synchronous generator is located is obtained t The real part I x and the imaginary part I y , and then get the actual active power P of the node where the synchronous generator is located t =V x I x ;

[0083] S1-4: According to the rotor motion equation of the synchronous generator, the differential equation of the synchronous generator is obtained:

[0084]

[0085] In formula (16), ω is the angular frequency of the synchronous generator, P tr is the output active power reference value of the synchronous generator, D t is the damping coefficient of the synchronous generator, T t is the inertia coefficient of the synchronous generator, E sv is the stator voltage of the excitation system in the synchronous generator, E ev is the excitation voltage of the excitation system in the synchronous generator, I d is the current of the synchronous generator on the d-axis, T stp is the subtransient time constant of the excitation system in the synchronous generator, E rev is the rated excitation voltage of the excitation system in the synchronous generator, G t is the proportional gain parameter of the excitation system in the synchronous generator, V tr is the reference voltage at the node where the synchronous generator is located, Amp is the amplitude symbol, T tc is the time constant of the excitation system in the synchronous generator;

[0086] The specific steps of step S2 are as follows:

[0087] S2-1: Get the initial information of the virtual synchronous machine, including: the line equivalent reactance X of the virtual synchronous machine f , the output reactive power reference value Q at the node where the virtual synchronous machine is located xr , the initial admittance of the virtual synchronous machine Y j The real part G j and the imaginary part B j ;

[0088] S2-2: Connect the virtual synchronous machine to the AC grid and update the parameters, including: the modified admittance Y at the node where the virtual synchronous machine is located j ′, the virtual injection current I at the node where the virtual synchronous machine is located j ′, where the modified admittance Y at the node where the virtual synchronous machine is located is j 'for:

[0089]

[0090] Virtual injection current I at the node where the virtual synchronous machine is located j ′=I xj ′+jI yj ′, where I xj ′ and I yj ′ are virtual injection currents I j The real and imaginary parts of ′ are expressed as follows:

[0091]

[0092] In formula (18), ρ is the power angle of the virtual synchronous machine, E d is the electromotive force of the virtual synchronous machine on the d-axis;

[0093] S2-3: Based on the modified admittance Y at the node where the virtual synchronous machine is located j ′ and virtual injection current I j ′, use Gaussian elimination method to directly solve the above linear equations (18) and get the actual voltage V at the node where the virtual synchronous machine is located t The real part of ′ V′ x and the imaginary part V′ y , and then according to the virtual injection current I at the node where the virtual synchronous machine is located j ′、initial admittance Y j and the actual voltage V t ′, and the injected current I at the node where the virtual synchronous machine is located is obtained t The real part of ′I′ x and the imaginary part I′ y , and then get the actual active power P′ of the node where the virtual synchronous machine is located t =V′ x I′ x and actual reactive power Q′t =V′ y I′ y ;

[0094] S2-4: According to the rotor motion equation of the virtual synchronous machine, the differential equation of the virtual synchronous machine is obtained:

[0095]

[0096] In formula (19), σ is the power angle of the virtual synchronous machine, ω x is the angular frequency of the virtual synchronous machine, ω xr is the rated angular frequency of the virtual synchronous machine, P xr is the output active power reference value at the node where the virtual synchronous machine is located, D x is the damping coefficient of the virtual synchronous machine, J x is the inertia coefficient of the virtual synchronous machine;

[0097] The specific steps of step S3 are as follows:

[0098] S3-1: according to the initial information of the synchronous generator and equation (16) in step S1-4, the differential equation is solved by using the improved Euler method to obtain a mathematical model of the differential equation of the synchronous generator connected to the power system;

[0099] S3-2: according to the initial information of the virtual synchronous machine and equation (19) in step S2-4, the differential equation is solved by using the improved Euler method to obtain a mathematical model of the differential equation of the virtual synchronous machine connected to the power system;

[0100] S3-3: Combining the differential equation mathematical models of step S3-1 and step S3-2, constructing a differential equation mathematical model of a power system including a synchronous generator and a virtual synchronous machine;

[0101] The specific steps of step S4 are as follows:

[0102] S4-1: Considering the minimum system inertia and damping support requirements, the objective function Obj is set as:

[0103]

[0104] In formula (20), M is the number of synchronous generators, and N is the number of virtual synchronous machines;

[0105] S4-2: Set the penalty function Pen to:

[0106] Pen=p·(|J total -N.J. x -M.T t |+|D total -N·D x -M.D t|) (21)

[0107] In Equation (21), p is the penalty coefficient, and J total is the total system inertia, and D total is the total system damping;

[0108] S4-3: Combining Step S4-1 and Step S4-2, design the fitness function fitness of the arithmetic algorithm as follows:

[0109] fitness = Obj + Pen (22)

[0110] S4-4: Take the damping D t of the synchronous generator, the inertia T t of the synchronous generator, the damping D x of the virtual synchronous machine, and the inertia J x of the virtual synchronous machine as the quantities to be optimized and generate the initial population;

[0111] S4-5: Set the mathematical optimization accelerator, and the specific expression is:

[0112]

[0113] In Equation (23), n is the current iteration number, and n max is the maximum iteration number, M min and M max are the minimum / maximum values of the corresponding function values of the mathematical optimization accelerator. Define a random number r 1 , when r 1 <MOA, enter the exploration stage in Step S4-6, otherwise enter the exploitation stage in Step S4-7;

[0114] S4-6: The specific expression of the exploration stage is:

[0115]

[0116] In Equation (24), r 2 is a random number uniformly distributed in [0, 1], x ij (n + 1) is the j-th position of the i-th solution in the next iteration, best(x j ) is the j-th position of the optimal solution in the iteration process so far, λ is a small constant 0.01 to prevent the denominator from being zero, OM j and ON j respectively represent the upper and lower bound values of the optimal value at the j-th position to prevent the solution from exceeding the boundary after updating the individual. ξ is the control parameter for adjusting the exploration process, and MOC is the mathematical optimization coefficient, and its calculation formula is:

[0117]

[0118] In formula (25), θ is a sensitive parameter;

[0119] S4-7: The specific expression of the development stage is:

[0120]

[0121] In formula (26), r 3 is a random number that follows a uniform distribution between [0, 1];

[0122] S4-8: Substitute the individual into formula (22), calculate the fitness value, and find a set of values ​​g with the minimum fitness value under the current number of iterations best , compare the best individual values ​​in all iterations, and select the set of values ​​with the smallest fitness value as the global best value f best , and determine whether the maximum number of iterations has been reached at this time. If so, jump to step S4-9; if not, jump to step S4-5;

[0123] S4-9: End the arithmetic algorithm, consider the minimum system inertia and damping support requirements, use the arithmetic algorithm to design the optimal inertia and optimal damping of the synchronous generator and the virtual synchronous machine, and substitute them into the differential equation mathematical model of the power system to optimize the dynamic characteristics of the power system.

[0124] Figure 3 The figure shows the fitness value change curve of the system after the arithmetic algorithm is adopted. After the arithmetic algorithm is adopted, the minimum fitness value is 6.81×10 -4 , indicating that this method effectively reduces the frequency and power fluctuations of the power system and significantly improves the dynamic characteristics of the power system.

[0125] Figure 4 and Figure 5 They are the inertia and damping change curves of two synchronous generators and one virtual synchronous machine after the arithmetic algorithm is adopted. Before the arithmetic algorithm is adopted, the inertias of the two synchronous generators and one virtual synchronous machine are designed to be 150, 150, and 700, and the dampings are 150, 150, and 700, respectively. After the arithmetic algorithm is adopted, the corresponding inertia combinations are 390, 187, and 423, and the corresponding damping combinations are 356, 256, and 388.

[0126] Figure 6 , Figure 7 and Figure 8These are the angular dynamic frequency curves of the first, second synchronous generators and virtual synchronous machine before and after planning, respectively. The maximum angular frequency amplitudes without planning are 1.00068 (pu), 1.00035 (pu) and 1.00015 (pu), and the adjustment times are 12.35s, 11.3s and 10s, respectively; the maximum angular frequency amplitudes after planning are 1.00025 (pu), 1.00024 (pu) and 1.0001 (pu), and the adjustment times are 8.2s, 6.7s and 8s, respectively. The maximum angular frequency amplitudes before and after planning are reduced by 0.00043 (pu), 0.00011 (pu) and 0.00005 (pu), respectively, and the adjustment times are reduced by 4.15s, 4.6s and 2s, respectively.

[0127] From the above analysis, it can be seen that the present invention constructs a model of connecting synchronous generators and virtual synchronous machines to the power system. Under the premise of ensuring the minimum system inertia and damping support requirements, the optimal inertia and optimal damping configuration of each unit are obtained through an arithmetic algorithm, thereby reducing the negative impact of parameter coupling between units and improving the dynamic performance of the system.

[0128] The embodiments described above are only preferred embodiments of the present invention and are not intended to limit the scope of implementation of the present invention. Therefore, all changes made according to the shape and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for parameter planning of a power system including synchronous generators and virtual synchronous machines, characterized in that: The following steps are involved: S1: Construct the differential equation mathematical model of synchronous generator; S1-1: Get the initial information of the synchronous generator, including: the transient reactance X of the synchronous generator on the d-axis and q-axis d ′、X q ′, the subtransient reactance X of the synchronous generator on the d-axis d ″, the transient resistance R of the synchronous generator i ′, the rated angular frequency of the synchronous generator ω tr , the initial admittance of the synchronous generator Y i The real part G i and the imaginary part B i ; S1-2: Connect the synchronous generator to the AC grid and update the parameters, including the corrected admittance Y at the node where the synchronous generator is located i ′, virtual injection current I at the node where the synchronous generator is located i ′, where the modified admittance Y at the node where the synchronous generator is located is i 'for: In formula (1), δ is the power angle of the synchronous generator; Virtual injection current I at the node where the synchronous generator is located i ′ can be expressed as the real part I xi ′ and the imaginary part I yi ′, that is, I i ′=I xi ′+jI yi ′, the specific expression is: In formula (2), E qi ′ is the transient electromotive force of the synchronous generator on the q axis; S1-3: Based on the modified admittance Y at the node where the synchronous generator is located i ′ and virtual injection current I i ′, use Gaussian elimination method to directly solve the above linear equations (2) to obtain the actual voltage V at the node where the synchronous generator is located t The real part V x and the imaginary part V y , and then according to the virtual injection current I at the node where the synchronous generator is located i ′、initial admittance Y i and the actual voltage V t , and the injected current I at the node where the synchronous generator is located is obtained t The real part I x and the imaginary part I y , and then get the actual active power P of the node where the synchronous generator is located t =V x I x ; S1-4: According to the rotor motion equation of the synchronous generator, the differential equation of the synchronous generator is obtained: In formula (3), ω is the angular frequency of the synchronous generator, P tr is the output active power reference value of the synchronous generator, D t is the damping coefficient of the synchronous generator, T t is the inertia coefficient of the synchronous generator, E sv is the stator voltage of the excitation system in the synchronous generator, E ev is the excitation voltage of the excitation system in the synchronous generator, I d is the current of the synchronous generator on the d-axis, T stp is the subtransient time constant of the excitation system in the synchronous generator, E rev is the rated excitation voltage of the excitation system in the synchronous generator, G t is the proportional gain parameter of the excitation system in the synchronous generator, V tr is the reference voltage at the node where the synchronous generator is located, Amp is the amplitude symbol, T tc is the time constant of the excitation system in the synchronous generator; S2: Construct the differential equation mathematical model of the virtual synchronous machine; S2-1: Get the initial information of the virtual synchronous machine, including: the line equivalent reactance X of the virtual synchronous machine f , the output reactive power reference value Q at the node where the virtual synchronous machine is located xr , the initial admittance of the virtual synchronous machine Y j The real part G j and the imaginary part B j ; S2-2: Connect the virtual synchronous machine to the AC grid and update the parameters, including: the modified admittance Y at the node where the virtual synchronous machine is located j ′, the virtual injection current I at the node where the virtual synchronous machine is located j ′, where the modified admittance Y at the node where the virtual synchronous machine is located is j 'for: Virtual injection current I at the node where the virtual synchronous machine is located j ′=I xj ′+jI yj ′, where I xj ′ and I yj ′ are virtual injection currents I j The real and imaginary parts of ′ are expressed as follows: In formula (5), ρ is the power angle of the virtual synchronous machine, E d is the electromotive force of the virtual synchronous machine on the d-axis; S2-3: Based on the modified admittance Y at the node where the virtual synchronous machine is located j ′ and virtual injection current I j ′, use Gaussian elimination method to directly solve the above linear equations (5) and get the actual voltage V at the node where the virtual synchronous machine is located t The real part of ′ V′ x and the imaginary part V′ y , and then according to the virtual injection current I at the node where the virtual synchronous machine is located j ′、initial admittance Y j and the actual voltage V t ′, and the injected current I at the node where the virtual synchronous machine is located is obtained t The real part of ′I′ x and the imaginary part I′ y , and then get the actual active power P′ of the node where the virtual synchronous machine is located t =V′ x I′ x and actual reactive power Q′ t =V′ y I′ y ; S2-4: According to the rotor motion equation of the virtual synchronous machine, the differential equation of the virtual synchronous machine is obtained: In formula (6), σ is the power angle of the virtual synchronous machine, ω x is the angular frequency of the virtual synchronous machine, ω xr is the rated angular frequency of the virtual synchronous machine, P xr is the output active power reference value at the node where the virtual synchronous machine is located, D x is the damping coefficient of the virtual synchronous machine, J x is the inertia coefficient of the virtual synchronous machine; S3: Connect the synchronous generator and the virtual synchronous machine to the system and construct a mathematical model of differential equations of the power system including the synchronous generator and the virtual synchronous machine; S3-1: according to the initial information of the synchronous generator and equation (3) in step S1-4, the differential equation is solved by using the improved Euler method to obtain a mathematical model of the differential equation of the synchronous generator connected to the power system; S3-2: according to the initial information of the virtual synchronous machine and equation (6) in step S2-4, the differential equation is solved by using the improved Euler method to obtain a differential equation mathematical model of the virtual synchronous machine connected to the power system; S3-3: Combining the differential equation mathematical models of step S3-1 and step S3-2, constructing a differential equation mathematical model of a power system including a synchronous generator and a virtual synchronous machine; S4: Based on the minimum system inertia and damping support requirements, the optimal inertia and optimal damping of the synchronous generator and virtual synchronous machine are designed using an arithmetic algorithm; S4-1: Considering the minimum system inertia and damping support requirements, the objective function Obj is set as: In formula (7), M is the number of synchronous generators, and N is the number of virtual synchronous machines; S4-2: Set the penalty function Pen to: In formula (8), p is the penalty coefficient, J total is the total inertia of the system, D total is the total damping of the system; S4-3: Combine steps S4-1 and S4-2 to design the fitness function of the arithmetic algorithm as: fitness=Obj+Pen (9) S4-4: Set the damping D of the synchronous generator t , Inertia of synchronous generator T t , the damping D of the virtual synchronous machine x , inertia of virtual synchronous machine J x As the quantity to be optimized and generate the initial population; S4-5: Set the mathematical optimization accelerator. The specific expression is: In formula (10), n is the current iteration number, and n max is the maximum iteration number, M min and M max are the minimum / maximum values of the corresponding function values of the mathematical optimization accelerator. Define a random number r1. When r1 < MOA, enter the exploration stage in step S4-6; otherwise, enter the development stage in step S4-7; S4-6: The specific expression of the exploration stage is: In formula (11), r2 is a random number that follows a uniform distribution between [0, 1], and x ij (n+1) is the jth position of the ith solution in the next iteration, best(x j ) is the jth position of the best solution in the iteration so far, λ is a small constant of 0.01 to prevent the denominator from being equal to zero, OM j and ON j They represent the upper and lower bounds of the optimal value at the jth position, respectively, to prevent the solution from crossing the boundary after the individual is updated. ξ is the control parameter for adjusting the exploration process. MOC is the mathematical optimization coefficient, and its calculation formula is: In formula (12), θ is a sensitive parameter; S4-7: The specific expression of the development stage is: In formula (13), r3 is a random number that obeys the uniform distribution of [0, 1]; S4-8: Substitute the individual into formula (9), calculate the fitness value, and find a set of values ​​g with the minimum fitness value under the current number of iterations best , compare the best individual values ​​in all iterations, and select the set of values ​​with the smallest fitness value as the global best value f best , and determine whether the maximum number of iterations has been reached at this time. If so, jump to step S4-9; if not, jump to step S4-5; S4-9: End the arithmetic algorithm, consider the minimum system inertia and damping support requirements, use the arithmetic algorithm to design the optimal inertia and optimal damping of the synchronous generator and the virtual synchronous machine, and substitute them into the differential equation mathematical model of the power system to optimize the dynamic characteristics of the power system.

Images