Bayesian optimization setting method for virtual inertia control parameters of wind turbine generator

Abstract

The invention relates to the technical field of power system automation, and particularly discloses a Bayesian optimization setting method for virtual inertia control parameters of a wind turbine generator, and the method comprises the following steps: constructing a disturbed performance index after a power system suffers from an expected fault; constructing a virtual inertia configuration optimization model of the new energy power system by taking the disturbed performance index as an optimization objective function; introducing a constraint condition into a target function in a penalty function form, and converting the virtual inertia configuration optimization model of the new energy power system into an unconstrained optimization model; and iteratively solving an optimal solution of the unconstrained optimization model by adopting a Bayesian optimization algorithm based on a Gaussian process. The method has the advantages that by optimizing the virtual inertia control parameters of the new energy unit, the power angle stability and the frequency stability of the new energy power system under anticipated faults such as an alternating current fault and generator off-grid are improved, and safe and stable operation of the novel power system with new energy as a main body is powerfully guaranteed.

Description

Technical Field

[0001] The invention relates to the technical field of power system automation, in particular to a Bayesian optimization setting method for virtual inertia control parameters of a wind turbine generator set. Background Art

[0002] As wind and solar energy gradually replace traditional thermal power units, new power systems with a high proportion of new energy are also facing increasingly severe challenges in safe and stable operation. Among them, doubly-fed asynchronous wind turbines and direct-drive wind turbines have introduced power electronic converter interfaces and achieved active-reactive decoupling control, which have better operation control capabilities than traditional constant-speed asynchronous wind turbines, and have gradually become the current mainstream wind power generation form. However, it is precisely because of the introduction of the converter interface that wind turbines are decoupled from the power grid in the traditional maximum power tracking control mode and cannot provide inertial support for the system. The equivalent rotational inertia of the power system has decreased, and it has gradually evolved into a "low inertia power system", resulting in a decrease in the system's anti-disturbance ability and frequency stability.

[0003] With the construction and development of new power systems with new energy as the main body, wind and solar new energy units still bear the obligation of active support and stable control of the power grid corresponding to their power generation capacity, so as to avoid the risk of stable operation caused by the decline of system inertia and frequency regulation capacity. In order to improve the active support capability of new energy units for the power grid, domestic and foreign researchers have carried out a lot of research. However, most of the technologies proposed in the existing literature focus on the optimization of the control strategy of the new energy units themselves, but rarely consider the overall perspective of the power system. Although various virtual inertia control strategies can effectively improve the overall equivalent inertia level of the system, different inertia distribution conditions may have different effects on the power angle stability of the system. Summary of the invention

[0004] The purpose of the present invention is to overcome the shortcomings of the prior art and provide a Bayesian optimization setting method for virtual inertia control parameters of a wind turbine set.

[0005] The object of the present invention is achieved by the following technical solution: a Bayesian optimization setting method for virtual inertia control parameters of a wind turbine generator set, the method comprising the following steps:

[0006] Step S1, considering the safety of power angle stability, frequency deviation and frequency change rate in the transient process of the power system, constructing the disturbance performance index of the power system after suffering the expected fault;

[0007] Step S2, taking the disturbance performance index as the optimization objective function, taking the virtual inertia and damping control parameters of the wind turbine as decision variables, and considering the system dynamic model constraints, power angle stability constraints, frequency deviation safety constraints, frequency change rate safety constraints, and feasible domain constraints of the decision variables, to construct a virtual inertia configuration optimization model for the new energy power system;

[0008] Step S3, introducing the constraint condition into the objective function in the form of a penalty function, and transforming the virtual inertia configuration optimization model of the new energy power system into an unconstrained optimization model;

[0009] Step S4: adopting the Bayesian optimization algorithm based on Gaussian process to iteratively obtain the optimal solution of the unconstrained optimization model, and taking it as the optimal solution of the virtual inertia control parameter optimization model of the new energy unit.

[0010] Specifically, the interference performance index is:

[0011]

[0012] Among them, J c represents the performance index of the disturbed system under the cth expected fault scenario; t0 represents the simulation start time, which is generally the 0th second; t s Indicates the simulation end time; Ω SG represents the set of synchronous generators such as hydropower and thermal power in the power system; θ i (t) and θ j (t) respectively represent the phase angles of synchronous generator i and synchronous generator j at time t; ω i (t) represents the deviation of the rotor angular velocity of synchronous generator i relative to the system synchronous speed, reflecting the system frequency offset; It represents the rate of change of the rotor angular velocity of synchronous generator i, reflecting the rate of change of system frequency.

[0013] Specifically, in step S2, minimizing the sum of system performance indicators under multiple expected failure scenarios is used as the optimization objective function, which is as follows:

[0014]

[0015] Among them, J c is the system performance index value under the cth expected failure scenario, Ω c This is the expected fault set.

[0016] Specifically, the system dynamic model constraints are:

[0017]

[0018] In the formula, Represents the power system dynamic model consisting of a set of differential algebraic equations;

[0019]

[0020] θ i =[θ i (t0),…,θ i (t),…,θ i (t s )], t∈[t0, t s ]

[0021]

[0022] ω i =[ω i (t0),…,ω i (t),…,ω i (t s )], t∈[t0, t s ]

[0023] In addition, x0 represents the initial value vector of the system differential variables, which include the state variables of each dynamic link such as the generator power angle and the rotor angular velocity; y0 represents the initial value vector of the system algebraic variables, including the grid node voltage amplitude\phase angle, branch active\reactive power flow and other grid flow variables; z represents the system input variable vector, including the active\reactive demand of each node load, the active\reactive output of each generator and other input variables reflecting the boundary conditions of the grid flow calculation; f c Indicates the cth expected failure scenario.

[0024] Specifically, the power angle stability constraint is:

[0025] -360°<θ i (t)-θ j (t)<360°

[0026]

[0027] In the formula, θ i (t) and θ j (t) represent the phase angles of synchronous generator i and synchronous generator j at time t respectively.

[0028] Specifically, the upper and lower limits of the frequency deviation are constrained as follows:

[0029]

[0030] In the formula, ω i (t) represents the deviation of the rotor angular velocity of synchronous generator i relative to the system synchronous speed; ωIndicates the lower limit of the rotor angular velocity, which is the frequency threshold of the first round of low-frequency load shedding in the power system; Indicates the upper limit of the rotor angular velocity, which is the frequency threshold of the first round of high-frequency cutting; Ω SG represents the set of synchronous generators in the power system; t0 represents the simulation start time, which is usually the 0th second; t s Indicates the simulation end time.

[0031] Specifically, the upper and lower limits of the frequency change rate are constrained as follows:

[0032]

[0033] In the formula, Synchronous generator i The rate of change of the rotor angular velocity; and Respectively represent the lower and upper limits of the rate of change;

[0034]

[0035] Where Δt is the time step of transient simulation.

[0036] Specifically, the upper and lower limit constraints of the decision variables are:

[0037]

[0038] w∈Ω DFIG

[0039] Where M w and D w are the virtual inertia and virtual damping control parameters of the wind power cluster w, are the decision variables of the optimization model, M and are the upper and lower limits of the virtual inertia parameters, D and are the upper and lower limits of the virtual damping parameter, Ω DFIG Represents the collection of doubly-fed induction wind farms.

[0040] Specifically, in step S3, the constraint condition is introduced into the objective function in the form of a penalty function to construct an extended objective function of the virtual inertia configuration optimization model of the new energy power system. Its mathematical expression is as follows:

[0041]

[0042] Among them, p c represents the penalty function associated with the constraints. When all constraints are satisfied, p c = 0, otherwise p c>0, its mathematical expression is as follows:

[0043]

[0044] Penalty term associated with power angle stability constraint The mathematical expression is as follows:

[0045]

[0046] The penalty term associated with the frequency deviation constraint The mathematical expression is as follows:

[0047]

[0048] Penalty term associated with frequency change rate constraint The mathematical expression is as follows:

[0049]

[0050] Specifically, the specific steps of step S4 are:

[0051] Step S41: Set the number of iterations N iter , initialize the iteration variable i = 0;

[0052] Step S42: Initialize the sampling point set Constructing and initializing a Gaussian process-based surrogate model Initialize the optimal solution s OPT is a null value, and the extended objective function corresponding to the optimal solution

[0053] Step S43: Considering the feasibility constraints of the decision variables, a set of 20 sampling points is generated using the Latin hypercube sampling method. where s j Represents the value vector of the jth sampling point, reflecting the value of the decision variable;

[0054] Step S44: for the sampling point set All samples in j , respectively, through transient simulation to obtain the trajectory data of the generator power angle θ and the rotor angular velocity ω under the expected fault, and calculate the corresponding extended objective function

[0055] Step S45: for the sampling point set Take the sampling point with the largest extended objective function as the optimal point of the current iteration step And the corresponding extended objective function is

[0056] Step S46: Compare With s OPT The size of make as well as Otherwise, the assignment operation is not performed;

[0057] Step S47: Update sampling point set For the sampling point set With s as the independent variable and α as the dependent variable, the surrogate model is updated by numerical fitting Model parameters of

[0058] Step S48: construct an acquisition function α based on the expected improvement (as shown in the following formula), combine the current Gaussian process proxy model, and use the pseudo-Newton method to calculate the optimal sampling point corresponding to the maximum value of the acquisition function.

[0059]

[0060] In the formula, represents mathematical expectation, max[a, b] represents the maximum value between a and b, then:

[0061]

[0062] In the formula, μ(s) and σ(s) are Gaussian process proxy models, respectively. The extended objective function predicted at the sampling point s The mean and standard deviation of the posterior distribution, Φ and φ are the cumulative density function and probability density function of the posterior distribution respectively; the parameter ξ is a comprehensive balance between the application item and the exploration item;

[0063] Step S49: for sampling points Through transient simulation, the trajectory data of the generator power angle θ and the rotor angular velocity ω under the expected fault are obtained, and the corresponding extended objective function is calculated. At the same time

[0064] Step S410: Let i=i+1, and determine whether i is greater than the maximum number of iterations N. iter If not, return to step S46; if so, end the algorithm process and use the current optimal solution s OPT as the optimal solution to the original problem.

[0065] The present invention has the following advantages:

[0066] The optimization method of the present invention gives full play to the active support capability of virtual inertia control of new energy units, and improves the power angle stability and frequency stability of the new energy power system under anticipated faults such as AC faults and generator disconnection by optimizing the virtual inertia control parameters of the new energy units, thereby effectively ensuring the safe and stable operation of the new power system with new energy as the main body. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] Figure 1 It is a schematic diagram of the optimization method flow of the present invention;

[0068] Figure 2 The flowchart of the Bayesian optimization algorithm of the present invention;

[0069] Figure 3 It is a schematic diagram of the composition structure of the double-fed induction fan of the present invention;

[0070] Figure 4 This is a block diagram of the virtual inertia control of the doubly-fed induction fan of the present invention. DETAILED DESCRIPTION

[0071] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention, that is, the embodiments described are only part of the embodiments of the present invention, rather than all of the embodiments. The components of the embodiments of the present invention described and shown in the drawings herein can be arranged and designed in various different configurations.

[0072] Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the invention claimed for protection, but merely represents selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without making creative work are within the scope of protection of the present invention.

[0073] It should be noted that relational terms such as "first" and "second" are used only to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Moreover, the terms "include", "comprise" or any other variants thereof are intended to cover non-exclusive inclusion, so that a process, method, article or device including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, article or device. In the absence of further restrictions, the elements defined by the sentence "including a..." do not exclude the existence of other identical elements in the process, method, article or device including the elements.

[0074] The present invention is further described below in conjunction with the accompanying drawings, but the protection scope of the present invention is not limited to the following description. Figures 1 to 4 As shown, a Bayesian optimization method for setting virtual inertia control parameters of a wind turbine generator system includes the following steps:

[0075] Step S1, considering the safety of power angle stability, frequency deviation and frequency change rate in the transient process of the power system, constructing the disturbance performance index of the power system after suffering the expected fault;

[0076] Step S2, taking the disturbance performance index as the optimization objective function, taking the virtual inertia and damping control parameters of the wind turbine as decision variables, and considering the system dynamic model constraints, power angle stability constraints, frequency deviation safety constraints, frequency change rate safety constraints, and feasible domain constraints of the decision variables, to construct a virtual inertia configuration optimization model for the new energy power system;

[0077] Step S3, introducing the constraint condition into the objective function in the form of a penalty function, and transforming the virtual inertia configuration optimization model of the new energy power system into an unconstrained optimization model;

[0078] Step S4: adopting the Bayesian optimization algorithm based on Gaussian process to iteratively obtain the optimal solution of the unconstrained optimization model, and taking it as the optimal solution of the virtual inertia control parameter optimization model of the new energy unit.

[0079] Furthermore, the interference performance index is:

[0080]

[0081] Among them, J c represents the performance index of the disturbed system under the cth expected fault scenario; t0 represents the simulation start time, which is generally the 0th second; t s Indicates the simulation end time, which can be 60 seconds; Ω SG represents the set of synchronous generators such as hydropower and thermal power in the power system; θ i (t) and θ j (t) respectively represent the phase angles of synchronous generator i and synchronous generator j at time t; ω i (t) represents the deviation of the rotor angular velocity of synchronous generator i relative to the system synchronous speed, reflecting the system frequency offset; It represents the rate of change of the rotor angular velocity of synchronous generator i, reflecting the rate of change of system frequency. The above performance indicator formula is explained as follows:

[0082] 1. The first term reflects the power angle difference of any pair of generators. When the disturbed system experiences transient power angle instability, the power angle difference between the leading unit and the lagging unit increases over time. Therefore, the worse the system power angle stability, the larger the first term.

[0083] 2. The second term reflects the deviation of the rotor angular velocity of any generator relative to the system synchronous speed. On the one hand, when the disturbed system has transient power angle instability, the angular velocity deviation of the unstable unit will increase over time. Therefore, the worse the system power angle stability, the larger the second term. On the other hand, when the disturbed system has frequency instability, the angular velocity deviation of each generator rotor directly reflects the change in system frequency. Therefore, it will increase (high frequency scenario caused by active power generation greater than load) or decrease (low frequency scenario caused by active power generation less than load) over time. Therefore, the worse the system frequency stability, the larger the second term.

[0084] 3. The third item reflects the frequency change rate. When the system inertia is low, the frequency change rate is relatively large in scenarios where active power imbalance is caused by generator disconnection, DC interlocking, etc. The larger the frequency change rate, the easier it is to cause high-frequency machine cutting or low-frequency load cutting. At the same time, for new energy units such as wind power and photovoltaic power that use power electronic converters as grid-connected interfaces, when the frequency change rate is large, their own grid-related protection will be activated and the unit will be disconnected from the grid to protect the safety of their own components. However, the chain disconnection of new energy sources will cause disordered changes in the system's active power imbalance, which is more likely to cause frequency instability. Therefore, when the system inertia is lower and the frequency stability is worse, the third item is larger.

[0085] Furthermore, in step S2, minimizing the sum of system performance indicators under multiple expected failure scenarios is used as the optimization objective function, which is specifically as follows:

[0086]

[0087] Among them, J c is the system performance index value under the cth expected failure scenario, Ω c This is the expected fault set.

[0088] Furthermore, the system dynamic model constraints are:

[0089]

[0090] In the formula, It represents the dynamic model of the power system composed of a group of differential algebraic equations; its calculation process is the transient simulation process of the power system, which is usually achieved by building a model of the target system in the relevant simulation software, setting the initial flow conditions and fault scenario conditions, and then performing transient simulation; in addition:

[0091]

[0092] θ i =[θ i (t0),…,θ i (t),…,θ i (t s )], t∈[t0, t s ]

[0093]

[0094] ω i =[ω i (t0),…,ω i (t),…,ω i (t s )], t∈[t0, t s ]

[0095] In addition, x0 represents the initial value vector of the system differential variables, which include the state variables of each dynamic link such as the generator power angle and the rotor angular velocity; y0 represents the initial value vector of the system algebraic variables, including the grid node voltage amplitude\phase angle, branch active\reactive power flow and other grid flow variables; z represents the system input variable vector, including the active\reactive demand of each node load, the active\reactive output of each generator and other input variables reflecting the boundary conditions of the grid flow calculation; f c Indicates the cth expected failure scenario.

[0096] Furthermore, the power angle stability constraint is:

[0097] -360°<θ i (t)-θ j (t)<360°

[0098]

[0099] In the formula, θ i (t) and θ j (t) represents the phase angle of synchronous generator i and synchronous generator j at time t respectively; the power angle stability constraint considers that the power angle difference between any two synchronous generators does not exceed 360° at any time of simulation;

[0100] Furthermore, the upper and lower limits of the frequency deviation are constrained as follows:

[0101]

[0102] In the formula, ω i (t) represents the deviation of the rotor angular velocity of synchronous generator i relative to the system synchronous speed; ω Indicates the lower limit of the rotor angular velocity, which is the frequency threshold of the first round of low-frequency load shedding in the power system; Indicates the upper limit of the rotor angular velocity, and takes the frequency threshold of the first round of high-frequency cutting action; the upper and lower limits of the frequency deviation constraints consider that at any time of the simulation, the rotor angular velocity of any synchronous generator is within the safe range; Ω SG represents the set of synchronous generators in the power system; t0 represents the simulation start time, which is usually the 0th second; t s Indicates the simulation end time.

[0103] Furthermore, the upper and lower limits of the frequency change rate are constrained as follows:

[0104]

[0105] In the formula, represents the rate of change of the rotor angular velocity of synchronous generator i; and Respectively represent the lower and upper limits of the rate of change; the upper and lower limits of the frequency deviation constraints consider that at any time of the simulation, the rate of change of the rotor angular velocity of any synchronous generator is within a safe range;

[0106]

[0107] Where Δt is the time step of transient simulation, which is generally 50 microseconds for electromagnetic transient simulation.

[0108] Furthermore, the upper and lower limits of the decision variables are:

[0109]

[0110] w∈Ω DFIG

[0111] Where M w and D w are the virtual inertia and virtual damping control parameters of the wind power cluster w, are the decision variables of the optimization model, M and are the upper and lower limits of the virtual inertia parameters, D and are the upper and lower limits of the virtual damping parameter, Ω DFIG Represents the collection of doubly-fed induction wind farms.

[0112] Furthermore, in step S3, the constraint condition is introduced into the objective function in the form of a penalty function to construct an extended objective function of the virtual inertia configuration optimization model of the new energy power system. Its mathematical expression is as follows:

[0113]

[0114] Among them, p crepresents the penalty function associated with the constraints. When all constraints are satisfied, p c = 0, otherwise p c >0, its mathematical expression is as follows:

[0115]

[0116] Penalty term associated with power angle stability constraint The mathematical expression is as follows:

[0117]

[0118] The penalty term associated with the frequency deviation constraint The mathematical expression is as follows:

[0119]

[0120] Penalty term associated with frequency change rate constraint The mathematical expression is as follows:

[0121]

[0122] Furthermore, the specific steps of step S4 are:

[0123] Step S41: Set the number of iterations N iter , initialize the iteration variable i = 0;

[0124] Step S42: Initialize the sampling point set Constructing and initializing a Gaussian process-based surrogate model Initialize the optimal solution s OPT is a null value, and the extended objective function corresponding to the optimal solution

[0125] Step S43: Considering the feasibility constraints of the decision variables, a set of 20 sampling points is generated using the Latin hypercube sampling method. where s j Represents the value vector of the jth sampling point, reflecting the value of the decision variable;

[0126] Step S44: for the sampling point set All samples in j , respectively, through transient simulation to obtain the trajectory data of the generator power angle θ and the rotor angular velocity ω under the expected fault, and calculate the corresponding extended objective function

[0127] Step S45: for the sampling point set Take the sampling point with the largest extended objective function as the optimal point of the current iteration step And the corresponding extended objective function is

[0128] Step S46: Compare With s OPT The size of make as well as Otherwise, the assignment operation is not performed;

[0129] Step S47: Update sampling point set For the sampling point set With s as the independent variable and α as the dependent variable, the surrogate model is updated by numerical fitting Model parameters of

[0130] Step S48: construct an acquisition function α based on the expected improvement (as shown in the following formula), combine the current Gaussian process proxy model, and use the pseudo-Newton method to calculate the optimal sampling point corresponding to the maximum value of the acquisition function.

[0131]

[0132] In the formula, represents mathematical expectation, max[a, b] represents the maximum value between a and b, then:

[0133]

[0134] In the formula, μ(s) and σ(s) are Gaussian process proxy models, respectively. The extended objective function predicted at the sampling point s The mean and standard deviation of the posterior distribution of , Φ and φ are the cumulative density function (CDF) and probability density function (PDF) of the posterior distribution, respectively. The first term in the above formula can measure the expected improvement of the sampling point s on the optimization problem, which is called the application term (Exploitation); the second term measures the uncertainty of the Gaussian process surrogate model at the sampling point s, which is called the exploration term (Exploration); the parameter ξ is a comprehensive balance between the application term and the exploration term;

[0135] Step S49: for sampling points Through transient simulation, the trajectory data of the generator power angle θ and the rotor angular velocity ω under the expected fault are obtained, and the corresponding extended objective function is calculated. At the same time

[0136] Step S410: Let i=i+1, and determine whether i is greater than the maximum number of iterations N. iter If not, return to step S46; if so, end the algorithm process and use the current optimal solution s OPT as the optimal solution to the original problem.

[0137] The present invention takes the doubly-fed induction wind turbine and the classic virtual inertia control strategy as examples, but it is adaptable and extensible to various new energy generator sets such as permanent magnet direct-drive wind turbines that use power electronic converters as grid-connected interfaces, as well as various improved virtual inertia controllers in the future;

[0138] The structure of the doubly fed induction fan is as follows: Figure 3 As shown;

[0139] The block diagram of virtual inertia control of double-fed induction fan is as follows Figure 4 As shown;

[0140] In the figure, P MPPT It represents the active power signal generated by the pitch angle control link of the wind turbine according to the wind conditions, f m represents the voltage frequency of the wind power grid connection point measured by the phase-locked loop (PLL), f ref Indicates the frequency reference value (i.e. the rated frequency of the system), P VI Indicates the active power control amount generated by virtual inertia control, P VD represents the active power control quantity generated by virtual damping control, P ref Indicates comprehensive P MPPT , P VI and P VD The active power control signal of the three, Q ref Represents the reactive power control signal generated by the reactive power control link of the wind turbine, P ref and Q ref The voltage reference signal V of the rotor-side converter is generated after pulse width modulation (PWM). ref .

[0141] The above is only a preferred embodiment of the present invention and does not limit the present invention in any form. Any technician familiar with the art can make many possible changes and modifications to the technical solution of the present invention by using the above-mentioned technical content without departing from the scope of the technical solution of the present invention, or modify it into an equivalent embodiment of equivalent changes. Therefore, any changes, modifications, equivalent changes and modifications made to the above embodiments based on the technology of the present invention without departing from the content of the technical solution of the present invention belong to the protection scope of the present technical solution.

Claims

1. A Bayesian optimization method for setting virtual inertia control parameters of wind turbines, characterized in that: The method comprises the following steps: Step S1, considering the safety of power angle stability, frequency deviation and frequency change rate in the transient process of the power system, constructing the disturbance performance index of the power system after suffering the expected fault; Step S2, taking the disturbance performance index as the optimization objective function, taking the virtual inertia and damping control parameters of the wind turbine as decision variables, and considering the system dynamic model constraints, power angle stability constraints, frequency deviation safety constraints, frequency change rate safety constraints, and feasible domain constraints of the decision variables, to construct a virtual inertia configuration optimization model for the new energy power system; Step S3, introducing the constraint condition into the objective function in the form of a penalty function, and transforming the virtual inertia configuration optimization model of the new energy power system into an unconstrained optimization model; Step S4: adopting the Bayesian optimization algorithm based on Gaussian process to iteratively obtain the optimal solution of the unconstrained optimization model, and taking it as the optimal solution of the virtual inertia control parameter optimization model of the new energy unit.

2. The Bayesian optimization setting method for virtual inertia control parameters of a wind turbine generator set according to claim 1 is characterized in that: The interference performance index is: Among them, J c represents the performance index of the disturbed system under the cth expected fault scenario; t0 represents the simulation start time, which is generally the 0th second; t s Indicates the simulation end time; Ω SG represents the set of synchronous generators such as hydropower and thermal power in the power system; θ i (t) and θ j (t) respectively represent the phase angles of synchronous generator i and synchronous generator j at time t; ω i (t) represents the deviation of the rotor angular velocity of synchronous generator i relative to the system synchronous speed, reflecting the system frequency offset; It represents the rate of change of the rotor angular velocity of synchronous generator i, reflecting the rate of change of system frequency.

3. The Bayesian optimization setting method for virtual inertia control parameters of a wind turbine generator set according to claim 2 is characterized in that: In step S2, the optimization objective function is to minimize the sum of system performance indicators under multiple expected failure scenarios, which is as follows: Among them, J c is the system performance index value under the cth expected failure scenario, Ω c This is the expected fault set.

4. The Bayesian optimization setting method for virtual inertia control parameters of a wind turbine generator set according to claim 3 is characterized in that: The system dynamic model constraints are: In the formula, Represents the dynamic model of the power system consisting of a set of differential algebraic equations; θ i =[θ i (t0),…,θ i (t),…,θ i (t s )],t∈[t0,t s ] ω i =[ω i (t0),…,ω i (t),…,ω i (t s )],t∈[t0,t s ] In addition, x0 represents the initial value vector of the system differential variables, which include the state variables of each dynamic link such as the generator power angle and the rotor angular velocity; y0 represents the initial value vector of the system algebraic variables, including the grid node voltage amplitude\phase angle, branch active\reactive power flow and other grid flow variables; z represents the system input variable vector, including the active\reactive demand of each node load, the active\reactive output of each generator and other input variables reflecting the boundary conditions of the grid flow calculation; f c Indicates the cth expected failure scenario.

5. The Bayesian optimization setting method for virtual inertia control parameters of a wind turbine generator set according to claim 4 is characterized in that: The power angle stability constraint is: -360°<θ i (t)-θ j (t)<360° In the formula, θ i (t) and θ j (t) represent the phase angles of synchronous generator i and synchronous generator j at time t respectively.

6. The Bayesian optimization setting method for virtual inertia control parameters of a wind turbine generator set according to claim 5 is characterized in that: The upper and lower limits of the frequency deviation are constrained as follows: In the formula, ω i (t) represents the deviation of the rotor angular velocity of synchronous generator i relative to the system synchronous speed; ω Indicates the lower limit of the rotor angular velocity, which is the frequency threshold of the first round of low-frequency load shedding in the power system; Indicates the upper limit of the rotor angular velocity, which is the frequency threshold of the first round of high-frequency cutting; Ω SG represents the set of synchronous generators in the power system; t0 represents the simulation start time, which is usually the 0th second; t s Indicates the simulation end time.

7. The Bayesian optimization setting method for virtual inertia control parameters of a wind turbine generator set according to claim 6 is characterized in that: The upper and lower limits of the frequency change rate are: In the formula, represents the rate of change of the rotor angular velocity of synchronous generator i; and Respectively represent the lower and upper limits of the rate of change; Where Δt is the time step of transient simulation.

8. The Bayesian optimization setting method for virtual inertia control parameters of a wind turbine generator set according to claim 7 is characterized in that: The upper and lower bounds of the decision variables are: w∈Ω DFIG Where M w and D w are the virtual inertia and virtual damping control parameters of the wind power cluster w, are the decision variables of the optimization model, M and are the upper and lower limits of the virtual inertia parameters, D and are the upper and lower limits of the virtual damping parameter, Ω DFIG Represents the collection of doubly-fed induction wind farms.

9. The Bayesian optimization setting method for virtual inertia control parameters of a wind turbine generator set according to claim 3 is characterized in that: In step S3, the constraint condition is introduced into the objective function in the form of a penalty function to construct an extended objective function of the virtual inertia configuration optimization model of the new energy power system. Its mathematical expression is as follows: Among them, p c represents the penalty function associated with the constraints. When all constraints are satisfied, p c = 0, otherwise p c >0, its mathematical expression is as follows: Penalty term associated with power angle stability constraint The mathematical expression is as follows: Penalty term associated with frequency deviation constraint The mathematical expression is as follows: Penalty term associated with frequency change rate constraint The mathematical expression is as follows:

10. The Bayesian optimization setting method for virtual inertia control parameters of a wind turbine generator set according to claim 9, characterized in that: The specific steps of step S4 are: Step S41: Set the number of iterations N iter , initialize the iteration variable i = 0; Step S42: Initialize the sampling point set Constructing and initializing a Gaussian process-based surrogate model Initialize the optimal solution s OPT is a null value, and the extended objective function corresponding to the optimal solution Step S43: Considering the feasibility constraints of the decision variables, a set of 20 sampling points is generated using the Latin hypercube sampling method. where s j Represents the value vector of the jth sampling point, reflecting the value of the decision variable; Step S44: for the sampling point set All samples in j , respectively, through transient simulation to obtain the trajectory data of the generator power angle θ and the rotor angular velocity ω under the expected fault, and calculate the corresponding extended objective function Step S45: for the sampling point set Take the sampling point with the largest extended objective function as the optimal point of the current iteration step And the corresponding extended objective function is Step S46: Compare With s OPT The size of make as well as Otherwise, the assignment operation is not performed; Step S47: Update sampling point set For the sampling point set With s as the independent variable and α as the dependent variable, the surrogate model is updated by numerical fitting Model parameters of Step S48: construct an acquisition function α based on the expected improvement (as shown in the following formula), combine the current Gaussian process proxy model, and use the pseudo-Newton method to calculate the optimal sampling point corresponding to the maximum value of the acquisition function. In the formula, represents mathematical expectation, max[a, b] represents the maximum value between a and b, then: In the formula, μ(s) and σ(s) are Gaussian process proxy models, respectively. The extended objective function predicted at the sampling point s The mean and standard deviation of the posterior distribution, Φ and φ are the cumulative density function and probability density function of the posterior distribution respectively; the parameter ξ is a comprehensive balance between the application item and the exploration item; Step S49: for sampling points Through transient simulation, the trajectory data of the generator power angle θ and the rotor angular velocity ω under the expected fault are obtained, and the corresponding extended objective function is calculated. At the same time Step S410: Let i=i+1, and determine whether i is greater than the maximum number of iterations N. iter If not, return to step S46; if so, end the algorithm process and use the current optimal solution s OPT as the optimal solution to the original problem.

Images