CN115173410A - Power grid rotational inertia threshold calculation method for improving GWO algorithm - Google Patents
Power grid rotational inertia threshold calculation method for improving GWO algorithm Download PDFInfo
- Publication number
- CN115173410A CN115173410A CN202210904002.8A CN202210904002A CN115173410A CN 115173410 A CN115173410 A CN 115173410A CN 202210904002 A CN202210904002 A CN 202210904002A CN 115173410 A CN115173410 A CN 115173410A
- Authority
- CN
- China
- Prior art keywords
- inertia
- representing
- model
- wolf
- frequency
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 80
- 238000004364 calculation method Methods 0.000 title claims abstract description 49
- 238000000034 method Methods 0.000 claims abstract description 140
- 238000005457 optimization Methods 0.000 claims abstract description 62
- 230000004044 response Effects 0.000 claims abstract description 45
- 238000012544 monitoring process Methods 0.000 claims abstract description 7
- 241000282461 Canis lupus Species 0.000 claims description 166
- 230000006870 function Effects 0.000 claims description 100
- 230000008569 process Effects 0.000 claims description 51
- 238000012546 transfer Methods 0.000 claims description 46
- 230000008859 change Effects 0.000 claims description 22
- 230000033228 biological regulation Effects 0.000 claims description 18
- 238000013178 mathematical model Methods 0.000 claims description 17
- 230000010354 integration Effects 0.000 claims description 16
- 239000011159 matrix material Substances 0.000 claims description 15
- 230000001360 synchronised effect Effects 0.000 claims description 9
- 238000013016 damping Methods 0.000 claims description 8
- 230000007423 decrease Effects 0.000 claims description 6
- 230000005619 thermoelectricity Effects 0.000 claims description 6
- 238000004458 analytical method Methods 0.000 claims description 5
- 230000006872 improvement Effects 0.000 claims description 5
- 230000001105 regulatory effect Effects 0.000 claims description 5
- 238000006467 substitution reaction Methods 0.000 claims description 5
- 238000009826 distribution Methods 0.000 claims description 4
- 229910052739 hydrogen Inorganic materials 0.000 claims description 4
- 238000010248 power generation Methods 0.000 claims description 4
- 230000005611 electricity Effects 0.000 claims description 3
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 claims description 2
- 230000003449 preventive effect Effects 0.000 abstract description 2
- 238000004088 simulation Methods 0.000 description 13
- 230000008901 benefit Effects 0.000 description 4
- 238000009825 accumulation Methods 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 3
- 238000010276 construction Methods 0.000 description 2
- 230000009467 reduction Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000009466 transformation Effects 0.000 description 2
- UFHFLCQGNIYNRP-UHFFFAOYSA-N Hydrogen Chemical compound [H][H] UFHFLCQGNIYNRP-UHFFFAOYSA-N 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 239000001257 hydrogen Substances 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000009897 systematic effect Effects 0.000 description 1
Images
Classifications
-
- H—ELECTRICITY
- H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
- H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
- H02J3/00—Circuit arrangements for ac mains or ac distribution networks
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/004—Artificial life, i.e. computing arrangements simulating life
- G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]
-
- H—ELECTRICITY
- H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
- H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
- H02J3/00—Circuit arrangements for ac mains or ac distribution networks
- H02J3/24—Arrangements for preventing or reducing oscillations of power in networks
- H02J3/241—The oscillation concerning frequency
-
- H—ELECTRICITY
- H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
- H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
- H02J3/00—Circuit arrangements for ac mains or ac distribution networks
- H02J3/38—Arrangements for parallely feeding a single network by two or more generators, converters or transformers
- H02J3/381—Dispersed generators
-
- H—ELECTRICITY
- H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
- H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
- H02J2203/00—Indexing scheme relating to details of circuit arrangements for AC mains or AC distribution networks
- H02J2203/20—Simulating, e g planning, reliability check, modelling or computer assisted design [CAD]
-
- H—ELECTRICITY
- H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
- H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
- H02J2300/00—Systems for supplying or distributing electric power characterised by decentralized, dispersed, or local generation
- H02J2300/20—The dispersed energy generation being of renewable origin
Landscapes
- Engineering & Computer Science (AREA)
- Power Engineering (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Evolutionary Computation (AREA)
- General Health & Medical Sciences (AREA)
- Biomedical Technology (AREA)
- Biophysics (AREA)
- Computational Linguistics (AREA)
- Data Mining & Analysis (AREA)
- Life Sciences & Earth Sciences (AREA)
- Artificial Intelligence (AREA)
- Molecular Biology (AREA)
- Computing Systems (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Software Systems (AREA)
- Health & Medical Sciences (AREA)
- Supply And Distribution Of Alternating Current (AREA)
Abstract
The invention discloses a power grid rotational inertia threshold calculation method for improving a GWO algorithm, which comprises the following steps of: monitoring power grid equipment and operating parameters, and constructing a system frequency response model; converting the model into a standard state space form, and calculating the frequency of the disturbed system; constructing an inertia threshold optimization model of an inertia power grid; solving the optimization model by adopting an improved GWOO algorithm, and outputting a calculation result of the inertia threshold of the power grid; the method for calculating the power grid rotational inertia threshold value of the improved GWO algorithm can calculate and simulate the disturbed system frequency fluctuation condition of a power grid more accurately, and can greatly improve the convergence and the solving precision of a solving model. And further, a more accurate inertia threshold calculation result is provided for the power grid, and a basis is provided for a follow-up power grid to adopt preventive control aiming at the frequency safety problem.
Description
Technical Field
The invention relates to the technical field of power system operation and control, in particular to a power grid rotational inertia threshold calculation method for improving a GWO algorithm.
Background
The rotational inertia of the power system is an important factor affecting the frequency safety of the power system. With the increasing of the new energy power generation ratio in recent years, on one hand, the synchronous unit with the same capacity is replaced, so that the system inertia is reduced; on the other hand, part of the new energy source unit with inertia support capability can provide virtual inertia support for the power grid, and the change of the operation mode of the new energy source unit can also cause the provided inertia support level to change. The safe and stable operation of the power system is affected by too high or too low frequency of the system, for example, the too high frequency may cause a high frequency tripping, and the too low frequency may cause a low frequency load shedding action. It is therefore necessary to evaluate the required threshold value of the moment of inertia of the power system according to the frequency safety requirements of the power system and to implement further control according to the threshold calculation results.
At present, regarding the calculation of an inertia threshold, on one hand, the corresponding minimum system inertia is calculated according to the safety constraint of the system frequency change rate in the initial stage after disturbance, a calculation model is over simplified, and the precision is lower; on the other hand, safety constraints such as the maximum value and the minimum value of the disturbed system frequency are considered, a corresponding optimization model is constructed, and the inertia threshold is optimized. Because when the system frequency change rate and the frequency maximum/minimum value constraint are considered at the same time, simulation and calculation of a disturbed system frequency curve need to be considered. In the literature on inertia threshold calculation, the method for calculating the system frequency mainly includes two categories: time domain simulation method and equivalent model method. The time domain simulation method is convenient for small systems, but is difficult to be applied to rapid simulation calculation of complex large systems. Regarding the equivalent model method, part of the method is to use a complex system equivalent as a single-machine load model, and an approximate analytic solution of frequency can be obtained, but the accuracy is relatively poor because the influence of non-linear factors such as dead zone and saturation in a power grid and the influence of power regulation and control dynamics of various control resources are not considered. From the above, the solution of the inertia threshold optimization model will involve a complex system frequency calculation process, and belongs to a complex nonlinear optimization problem, so that a proper solution algorithm is necessary.
Aiming at the analysis, the method researches how to consider the complex nonlinear factors in the power system and accurately calculates the frequency fluctuation curve of the system; and based on the frequency calculation result, the inertia threshold value required by the power system is calculated, which is an urgent problem to be solved in the research of maintaining the frequency safety of a novel power system.
Disclosure of Invention
This section is for the purpose of summarizing some aspects of embodiments of the invention and to briefly introduce some preferred embodiments. In this section, as well as in the abstract and title of the application, simplifications or omissions may be made to avoid obscuring the purpose of the section, the abstract and the title, and such simplifications or omissions are not intended to limit the scope of the invention.
The present invention has been made in view of the above-mentioned problems.
Therefore, the technical problem solved by the invention is as follows: and the frequency safety in the operation process of the power grid is guaranteed.
In order to solve the technical problems, the invention provides the following technical scheme: a power grid rotational inertia threshold calculation method for improving a GWO algorithm comprises the following steps:
monitoring power grid equipment and operating parameters, and constructing a system frequency response model;
converting the model into a standard state space form, and calculating the frequency of the disturbed system;
constructing an inertia threshold optimization model of an inertia power grid;
and solving the optimization model by adopting an improved GWO algorithm, and outputting a calculation result of the inertia threshold of the power grid.
As an optimal solution of a method for calculating a grid moment of inertia threshold value of a GWO algorithm, the method comprises the following steps: the power grid system to which the inertia threshold calculation method is applied generally includes: the system comprises an alternating current power grid, a conventional unit, a new energy station and a corresponding frequency modulation controller.
As a preferred scheme of the method for calculating the grid moment of inertia threshold value by improving the GWO algorithm, the method comprises the following steps: the conventional unit is generally provided with a speed regulator, the new energy station part is provided with a frequency modulation controller and comprises a virtual inertia support system, and the conventional unit and the new energy station belong to a power generation unit in an alternating current power grid; and calculating the reference value of the inertia parameter in the virtual inertia support system of the new energy station by the inertia threshold optimization method.
As an optimal solution of a method for calculating a grid moment of inertia threshold value of a GWO algorithm, the method comprises the following steps: the constructing of the system frequency response model comprises:
constructing a frequency feedback frequency modulation model related to the conventional unit according to the monitored information of the capacity, the reserved spare, the frequency modulation controller structure and the like of the station of the conventional unit; the difference mainly exists on water and electricity and thermoelectricity in concrete structure, wherein to thermoelectricity unit, the mathematical model that its structure corresponds includes frequency modulation blind spot link, speed regulator and prime mover link, the power regulation saturation link, specifically as follows:
wherein Δ f represents a system frequency deviation; Δ f act Representing the amount of frequency deviation across the frequency modulation dead zone; delta P g.ref Representing an active regulation reference value based on frequency feedback; delta P g Representing an actual power adjustment; delta P g.up And Δ P g.down Respectively representing the upper limit and the lower limit of active power adjustment; t is G Represents the governor time constant; t is CH Representing the turbine time constant; t is RH Represents a reheat time constant; f HP Represents a reheat coefficient; k is a radical of p Expressing a unit power adjustment coefficient;
the method comprises the following steps of constructing a mathematical model corresponding to a frequency modulation controller of the hydroelectric generating set, similarly comprising a frequency modulation dead zone link, a speed regulator and prime motor link and a power regulation saturation link, wherein the mathematical model is mainly different from a thermal power generating set in that a transfer function corresponding to the speed regulator and the prime motor link is as follows:
wherein R is T Is a temporary rate of decline; r is P Is the permanent rate of decline;
according to the monitored information of the capacity, the reserved spare and the frequency modulation controller structure of the new energy station and the like, a frequency feedback frequency modulation model related to the new energy station is constructed, and the method comprises the following steps: the frequency modulation dead zone, the wave filter, primary frequency modulation and inertia support link, several parts of power saturation link, wherein it is primary frequency modulation and inertia support link to be great to distinguish with conventional unit, specifically as follows:
wherein, Δ P re.ref Indicating a power regulation reference value of the new energy station participating in frequency modulation; k is re And H re Respectively representing a frequency droop coefficient and a virtual inertia coefficient of a frequency modulation controller of the new energy station, namely the frequency modulation controller consists of a primary frequency modulation part and a virtual inertia control part;
the rotor equation of motion is embodied as follows:
wherein H g Representing the inertia constant of the generator; delta P gm Representing the variation of the mechanical power of the generator; delta P ge Representing the variation of electromagnetic power at the terminal; d g Representing the damping coefficient of the generator; Δ f g Representing generator rotor frequency deviation;
equating the whole power grid into a synchronous generator based on a rotor motion equation, wherein the frequency dynamics of the equivalent generator is similar to the rotor motion equation, and the equivalence method comprises the following steps:
wherein, Δ P sm Representing the mechanical power of the equivalent generator; delta P se Representing the generator terminal electromagnetic power of the equivalent generator; d s Representing the damping coefficient of the equivalent generator; h s Representing the inertia constant of the equivalent generator;
based on the mathematical models of the frequency modulation controllers of the conventional unit and the new energy station and the rotor motion equation of the system equivalent generator, a corresponding system frequency response model is constructed, and the method specifically comprises the following steps:
calculating the inertia constant of the equivalent generator of the system through the inertia constant and the capacity of each conventional unit:
wherein,
and with
Respectively representing the inertia constant and the capacity of the ith conventional unit; n is a radical of sg Representing the total number of conventional units;
based on the frequency modulation link based on frequency feedback of conventional units such as hydropower, thermal power and the like, the actual power regulating quantity delta P of each station is adjusted g Summing and corresponding system equivalent generator rotor motion equation mechanical power variation, namely:
adjusting the actual power of each station by delta P based on the frequency modulation link of the new energy station based on the frequency feedback controller re Summing the predicted power disturbance quantity of the N-1 fault, and corresponding to the variation quantity of the electromagnetic power at the equivalent generator end of the system, namely:
wherein, Δ P dis Representing the active disturbance quantity corresponding to the expected N-1 fault; n is a radical of re Representing the total number of new energy stations with frequency modulation control function;
and calculating various variables in the rotor motion equation based on the rotor motion equation of the system equivalent generator, according to the inertia and other information of the conventional unit, the power adjustment quantity of all the units based on frequency feedback and the expected fault disturbance information, and finally completing the construction of a system frequency response model in a transfer function form.
As a preferred scheme of the method for calculating the grid moment of inertia threshold value by improving the GWO algorithm, the method comprises the following steps: converting the system frequency response model into a standard state space form comprises:
and (3) formally arranging the transfer function model according to the relation between input and output:
wherein F(s) represents an nth order transfer function; y and U respectively represent output and input variables corresponding to the transfer function model; a is a 0 ～a n And b 0 ～b n Coefficients representing terms of respective orders in the transfer function; s represents the laplacian operator;
and converting the relation of the s domain about the input variable and the output variable into a high-order differential equation form:
a n Y (n) +a n-1 Y (n-1) +L+a 1 Y′+a 0 Y＝b n U (n) +b n-1 U (n-1) +L+b 1 U′+b 0 U
wherein Y is (n) Represents the nth derivative of Y; in the same way, U (n) Represents the nth derivative of U; y 'and U' respectively represent corresponding first-order derivatives;
converting the differential equation into a state space model form by a variable substitution method, and enabling x to be 0 ＝Y；x 1 = Y', analogy to x n ＝Y (n) (ii) a In the same way, u 0 ＝U；u 1 = U'; analogize to u n ＝U (n) . The following relationship may be finally formed
By the above transformation, a composition comprisingA state space model of 2 × n state variables, the model being composed of 2 × n differential equations, and an algebraic equation; when considering a given model input u 0 ＝u set The state space model is a system of equations comprising 2 × (n + 1) variables, 2 × (n + 1) equations;
known based on model input, i.e. u 0 ＝u set The algebraic equations in the state-space model, in combination, may ultimately be transformed into the form of a standard state-space model, as follows:
wherein, A sub Representing a state matrix, x, corresponding to a transfer function F(s) sub Representing state variables therein, corresponding to x in the state space model 0 :x n-1 And u 0 :u n-1 ；
The system frequency response model according to the constructed transfer function form comprises a plurality of transfer functions, and the transfer functions are written into a state space model only comprising state variables; performing row writing on nonlinear links such as dead zones, power saturation and the like in the established system frequency response model through an algebraic equation; finally, collecting state variables of all transfer functions and algebraic variables in all algebraic equations, and integrally converting the system frequency response model into a state space model form as follows:
the system frequency model not only comprises linear links such as a transfer function and the like, but also comprises a nonlinear link, so that f represents a function corresponding to the state equation; g represents a function corresponding to an algebraic equation; the x table state variable matrix specifically comprises a system frequency difference value delta f and the like in a system frequency response model; y represents an algebraic variable matrix; u represents the input control variable matrix, e.g. the system active disturbance Δ P in the system frequency response model dist And so on.
Based on the initial operation state of the power grid and the expected fault disturbance, assigning initial values to state variables and algebraic variables in the state space model; the specific principle is as follows: giving an initial value of delta P to an input variable u (representing active disturbance) in a state space model dist And because the system initially operates in a stable state, the initial values of the other state variables are assigned to be 0, then the initial values of the corresponding algebraic variables are calculated according to the algebraic equation based on the initial values of the state variables, and finally the initial values of all the variables in the converted state space model can be obtained.
As an optimal solution of a method for calculating a grid moment of inertia threshold value of a GWO algorithm, the method comprises the following steps: the calculating the disturbed system frequency comprises:
the total fluctuation time length of the given frequency to be calculated is recorded as t sum And the numerical integration step is denoted as t bc ；
Based on the state space model and the initial variable value thereof, the specific process of calculating the system frequency by adopting the trapezoidal integral and Newton iteration method is as follows:
taking the value x according to the state variable at the current moment k (ii) a At the initial time, x k ＝x 0 (ii) a Output variable y k ，y k From x k Calculating to obtain the state variable x at the next moment by using the thought of trapezoidal integration k+1 The following are:
wherein: in the equation x k+1 Is a waiting amount, y k+1 Can be composed of x k+1 And (4) calculating. The trapezoidal integration method can reduce error accumulation in the successive iteration process, so that the method for calculating the system frequency by using the trapezoidal integration mode has certain advantages.
Solving variable x by adopting Newton iteration thought according to the solving principle of the nonlinear equation set k+1 The method comprises the following steps:
wherein: i =0,1,2, \8230;
representing the result of solving the equation at time k using Newton iterations after the ith iteration, and
representing the iteration step length obtained by the (i + 1) th iteration calculation;
represents the error after the (i + 1) th iteration; err (r) th Indicating a set error threshold; when the iteration error is less than the error threshold, the iteration is considered to be converged, i.e., by x k Iterative computation yields x k+1 ；
Continuously carrying out iterative computation on the system frequency at each moment until the total time t of frequency stability analysis is reached sum The requirements of (1).
As a preferred scheme of the method for calculating the grid moment of inertia threshold value by improving the GWO algorithm, the method comprises the following steps: the method for constructing the inertia threshold optimization model of the inertia power grid comprises the following steps:
the objective function with the minimum integral inertia provided by the conventional unit and the new energy unit in the optimization model is specifically as follows:
it should be noted that, since the rotor inertia of the conventional unit itself remains unchanged under the same unit combination, the sum of the virtual inertias provided by the new energy station only needs to be guaranteed to be minimum in the objective function.
The variables to be decided of the optimization model are virtual inertia coefficients provided by each new energy station:
the lowest point of the system frequency and the change rate of the power grid after the generator N-1 fault disturbance occurs in the optimization model meet the safety requirement, and the related constraint of the virtual inertia coefficient of the new energy station in the specified range is as follows:
wherein,
representing a virtual inertia coefficient provided by the ith new energy station;
the capacity of the ith new energy station is represented;
representing the moment of inertia provided by the jth conventional station; n is a radical of re Representing the number of new energy stations; Δ f max Expressing the maximum value of the system frequency deviation and taking a positive number; Δ f rocof Representing the rate of change of the system frequency;
representing a system frequency maximum deviation limit;
to representA system frequency rate of change limit;
representing the maximum limit value of the inertia of the new energy station;
and uniformly combining the objective function, the decision variables and the constraint conditions to construct the inertia threshold optimization model of the inertia power grid.
As a preferred scheme of the method for calculating the grid moment of inertia threshold value by improving the GWO algorithm, the method comprises the following steps: and converting the inertia threshold optimization model into a structure suitable for solving by a GWOO algorithm, wherein the definition of a fitness function is as follows:
wherein F (-) represents a fitness function; m represents a normal number 1 to 2 orders of magnitude greater than x; g (-) is an introduced piecewise function and represents a penalty term; when Δ f max And Δ f rocof After the limit is out of limit, the punishment item is a negative number, and the absolute value of the punishment item is rapidly increased due to the value of M; and when the requirement is met, the penalty term is 0;
in an optimization model structure suitable for a GWOO optimization algorithm, constraint conditions are as follows:
wherein,
representing the virtual inertia coefficient provided by the ith new energy station,
and representing the maximum limit value of the inertia of the new energy station.
As an optimal solution of a method for calculating a grid moment of inertia threshold value of a GWO algorithm, the method comprises the following steps: the location updating strategy based on the wolf α in the Levy flight improvement GWO algorithm specifically comprises the following steps: according to the original GWO algorithm: selecting a wolf with the maximum fitness in a wolf group at the beginning of each iteration of the algorithm as a wolf alpha in the current iteration process, and correspondingly updating the position and the fitness of the alpha; on the basis, according to a Levy flight principle and in combination with the optimal position in the historical search process of the wolf pack, the position of the wolf alpha is randomly updated once, and the method specifically comprises the following steps:
wherein,
represents the alpha wolf position after Levy flight; levy (λ) represents a lewy random number, which can be generally simulated by using a Mantegna algorithm, and can be specifically expressed as follows in combination with GWO:
wherein,
σ v ＝1；
wherein N (-) represents a normal distribution; Γ (·) represents a gamma function; the lambda is generally in the range of 1 to 3;
initializing a wolf group position according to a GWOO algorithm execution flow, calculating the wolf group individual fitness, selecting a leading wolf according to the fitness calculation result and updating a position coordinate, randomly updating an alpha wolf position according to Levy flight, updating the position of each individual in the wolf group according to the distance relation between the leading wolf and the wolf group individual, checking whether the wolf group individual coordinate exceeds a parameter value range, correcting the coordinate and entering next iteration until the maximum iteration number is reached, wherein the specific process comprises the following steps:
initializing the position of a wolf group: first the number of wolf clusters and the maximum number of iterations are given. Variable to be optimized
The number of the wolf clusters corresponds to the dimensionality of the position coordinates of the wolf clusters, and the initial coordinates of individuals in the wolf clusters are randomly or self-defined according to the value range of the variables to be optimized;
calculating the fitness of individuals in the wolf group: adding a penalty function corresponding to the constraint condition into the target function as the fitness of each wolf head in each iteration process;
updating the positions of the wolf α, β, δ: selecting three wolfs with higher fitness from all individuals of the wolf group, sequentially using the three wolfs as alpha, beta and delta wolfs in the iteration process, and correspondingly updating the position information;
randomly updating the position of the wolf alpha at one time according to Levy flight;
updating the position of the wolf group individual: depending on the distance between the wolf group individual and the wolf alpha, beta and delta, the position of each wolf group individual is updated as follows:
where, represents a vector dot product; a and C represent random coefficients, A =2a × r 1 -a，C＝2r 2 (ii) a a represents the step length, and can be generally set to be 0.01 times of the parameter value range;
and
representing the position coordinates of the k-th head wolf in the wolf group in the l and l +1 iteration processes;
represents the position of the alpha wolf in the first iteration process; d α Representing a distance variable between the alpha wolf and the wolf group individual;
is obtained by
And then, judging whether the position meets the range requirement of the parameter to be decided, if not, pulling the wolf position back to the boundary of the range, and then entering next iteration until the maximum iteration number is reached.
As a preferred scheme of the method for calculating the grid moment of inertia threshold value by improving the GWO algorithm, the method comprises the following steps: the grid inertia threshold calculation process comprises the following steps:
wherein,
representing a required inertia threshold of the power grid;
representing the inertia of the conventional unit;
representing the inertia of the new energy source unit obtained by optimization;
and with
And respectively representing the capacities of the conventional unit and the new energy station.
The invention has the beneficial effects that: according to the power grid rotational inertia threshold value calculation method for improving the GWO algorithm, provided by the invention, a corresponding system frequency response model is constructed by considering complex nonlinear factors such as generator power adjustment dynamic, frequency modulation dead zone, power saturation and the like, and then a numerical integration mode is adopted to calculate and simulate the system frequency fluctuation condition after the power grid is disturbed more accurately; by considering the participation condition of the new energy station, a corresponding system inertia threshold value optimization model is constructed, and a GWO algorithm improved based on the Levy flight principle is provided to solve the optimization model, so that the convergence and the solving precision of the solved model can be greatly improved. And a more accurate inertia threshold calculation result is provided for the power grid, and a basis is provided for the follow-up power grid to adopt preventive control aiming at the frequency safety problem.
Drawings
In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required to be used in the description of the embodiments are briefly introduced below, and it is obvious that the drawings in the description below are only some embodiments of the present invention, and it is obvious for those skilled in the art that other drawings can be obtained according to the drawings without inventive labor. Wherein:
fig. 1 is an overall flowchart of a method for calculating a grid inertia threshold value of a GWO algorithm according to a first embodiment of the present invention.
Fig. 2 is an exemplary system illustrating an inertia threshold calculation method in a simulation experiment of a grid rotational inertia threshold calculation method for improving a GWO algorithm according to a second embodiment of the present invention;
fig. 3 is a mathematical model corresponding to a frequency modulation controller of a hydroelectric generating set and a thermodynamic generating set in a simulation experiment of a method for calculating a rotational inertia threshold of a power grid for improving a GWO algorithm according to a second embodiment of the present invention;
fig. 4 is a mathematical model corresponding to a new energy station frequency modulation controller in a simulation experiment of a power grid rotational inertia threshold calculation method for improving a GWO algorithm according to a second embodiment of the present invention;
fig. 5 is a system frequency response model constructed in a simulation experiment of a method for calculating a grid rotational inertia threshold value by improving a GWO algorithm according to a second embodiment of the present invention;
fig. 6 is a flow of implementing the improved GWO algorithm in a simulation experiment of a method for calculating a grid moment of inertia threshold value by improving the GWO algorithm according to a second embodiment of the present invention;
fig. 7 shows a convergence condition of a wolf colony of a GWO algorithm in a simulation experiment of a method for calculating a grid moment of inertia threshold value by improving the GWO algorithm according to a second embodiment of the present invention.
Detailed Description
In order to make the aforementioned objects, features and advantages of the present invention more comprehensible, embodiments accompanying figures of the present invention are described in detail below, and it is apparent that the described embodiments are a part, not all or all of the embodiments of the present invention. All other embodiments, which can be obtained by a person skilled in the art without making any creative effort based on the embodiments in the present invention, shall fall within the protection scope of the present invention.
In the following description, numerous specific details are set forth in order to provide a thorough understanding of the present invention, but the present invention may be practiced in other ways than those specifically described and will be readily apparent to those of ordinary skill in the art without departing from the spirit of the present invention, and therefore the present invention is not limited to the specific embodiments disclosed below.
Furthermore, the references herein to "one embodiment" or "an embodiment" refer to a particular feature, structure, or characteristic that may be included in at least one implementation of the present invention. The appearances of the phrase "in one embodiment" in various places in the specification are not necessarily all referring to the same embodiment, nor are separate or alternative embodiments mutually exclusive of other embodiments.
The present invention will be described in detail with reference to the drawings, wherein the cross-sectional views illustrating the structure of the device are not enlarged partially in general scale for convenience of illustration, and the drawings are only exemplary and should not be construed as limiting the scope of the present invention. In addition, the three-dimensional dimensions of length, width and depth should be included in the actual fabrication.
Also in the description of the present invention, it should be noted that the terms "upper, lower, inner and outer" and the like indicate orientations or positional relationships based on the orientations or positional relationships shown in the drawings, which are only for convenience of description and simplification of description, but do not indicate or imply that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and thus, should not be construed as limiting the present invention. Furthermore, the terms first, second, or third are used for descriptive purposes only and are not to be construed as indicating or implying relative importance.
The terms "mounted, connected and connected" in the present invention are to be understood broadly, unless otherwise explicitly specified or limited, for example: can be fixedly connected, detachably connected or integrally connected; they may be mechanically, electrically, or directly connected, or indirectly connected through intervening media, or may be interconnected between two elements. The specific meanings of the above terms in the present invention can be understood in specific cases to those skilled in the art.
Example 1
Referring to fig. 1, a first embodiment of the present invention provides a method for calculating a grid rotational inertia threshold value of an improved GWO algorithm, including:
s1: monitoring power grid equipment and operating parameters, and constructing a system frequency response model;
further, the power grid system to which the inertia threshold calculation method is applied generally includes: the system comprises an alternating current power grid, a conventional unit, a new energy station and a corresponding frequency modulation controller.
The conventional unit is generally provided with a speed regulator, the new energy station part is provided with a frequency modulation controller and comprises a virtual inertia support system, and the conventional unit and the new energy station belong to a power generation unit in an alternating current power grid; and calculating the reference value of the inertia parameter in the virtual inertia support system of the new energy station by the inertia threshold optimization method.
Specifically, constructing the system frequency response model includes:
constructing a frequency feedback frequency modulation model related to the conventional unit according to the monitored information of the capacity, the reserved spare, the frequency modulation controller structure and the like of the station of the conventional unit; the difference mainly exists on water and electricity and thermoelectricity in concrete structure, wherein to thermoelectricity unit, the mathematical model that its structure corresponds includes frequency modulation blind spot link, speed regulator and prime mover link, the power regulation saturation link, specifically as follows:
wherein Δ f represents a system frequency deviation; Δ f act Representing the amount of frequency deviation across the frequency modulation dead zone; delta P g.ref Representing an active regulation reference value based on frequency feedback; delta P g Representing an actual power adjustment; delta P g.up And Δ P g.down Respectively representing the upper limit and the lower limit of active power adjustment; t is G Represents the governor time constant; t is CH Representing the turbine time constant; t is RH Represents a reheat time constant; f HP Represents a reheat coefficient; k is a radical of p Expressing a unit power adjustment coefficient;
the method comprises the following steps of constructing a mathematical model corresponding to a frequency modulation controller of the hydroelectric generating set, similarly comprising a frequency modulation dead zone link, a speed regulator and prime motor link and a power regulation saturation link, wherein the mathematical model is mainly different from a thermal power generating set in that a transfer function corresponding to the speed regulator and the prime motor link is as follows:
wherein R is T Is a temporary rate of decline; r P Is the permanent rate of decline;
according to the monitored information such as the capacity, the reserved spare and the frequency modulation controller structure of the new energy station, a frequency feedback frequency modulation model related to the new energy station is constructed, and the method comprises the following steps: the frequency modulation dead zone, the wave filter, primary frequency modulation and inertia support link, several parts of power saturation link, wherein it is primary frequency modulation and inertia support link to be great to distinguish with conventional unit, specifically as follows:
wherein, Δ P re.ref Indicating a power regulation reference value of the new energy station participating in frequency modulation; k is re And H re Respectively representing a frequency droop coefficient and a virtual inertia coefficient of a frequency modulation controller of the new energy station;
the rotor equation of motion is in the following specific form:
wherein H g Representing the inertia constant of the generator; delta P gm Representing the variation of the mechanical power of the generator; delta P ge Representing the variation of the electromagnetic power at the machine end; d g Representing a damping coefficient of the generator; Δ f g Representing generator rotor frequency deviation;
equating the whole power grid into a synchronous generator based on a rotor motion equation, wherein the frequency dynamics of the equivalent generator is similar to the rotor motion equation, and the equivalence method comprises the following steps:
wherein, Δ P sm Representing the mechanical power of the equivalent generator; delta P se Representing the generator terminal electromagnetic power of the equivalent generator; d s Representing the damping coefficient of the constant value generator; h s Representing the inertia constant of the equivalent generator;
and constructing a corresponding system frequency response model based on the mathematical models of the frequency modulation controllers of the conventional unit and the new energy station and the rotor motion equation of the system equivalent generator, wherein the corresponding system frequency response model is specifically as follows:
calculating the inertia constant of the equivalent generator of the system through the inertia constant and the capacity of each conventional unit:
wherein,
and
respectively representing the inertia constant and the capacity of the ith conventional unit; n is a radical of sg Representing the total number of conventional units;
based on the frequency modulation link based on frequency feedback of conventional units such as hydropower, thermal power and the like, the actual power regulating quantity delta P of each station is adjusted g Summing and corresponding system equivalent generator rotor motion equation mechanical power variation, namely:
adjusting the actual power of each station by the delta P based on the frequency modulation link of the new energy station based on the frequency feedback controller re Summing the predicted power disturbance quantity of the N-1 fault, and corresponding to the variation quantity of the electromagnetic power at the equivalent generator end of the system, namely:
wherein, Δ P dis Representing the active disturbance quantity corresponding to the expected N-1 fault; n is a radical of re Representing the total number of new energy stations with frequency modulation control functionality.
And calculating various variables in the rotor motion equation based on the rotor motion equation of the system equivalent generator, according to the inertia and other information of the conventional unit, the power adjustment quantity of all the units based on frequency feedback and the expected fault disturbance information, and finally completing the construction of a system frequency response model in a transfer function form.
It should be noted that through the above steps, complex characteristics such as power adjustment dynamics, frequency modulation dead zone, power saturation and the like of the conventional unit and the new energy unit are all considered in the constructed system frequency response model, and a more accurate basic model can be provided for subsequent frequency calculation.
S2: converting the model into a standard state space form, and calculating the frequency of the disturbed system;
further, converting the system frequency response model into a standard state space form comprises:
and (3) formally arranging the transfer function model according to the relation between input and output:
wherein F(s) represents an nth order transfer function; y and U respectively represent output and input variables corresponding to the transfer function model; a is 0 ～a n And b 0 ～b n Coefficients representing terms of various orders in the transfer function; s represents the laplacian operator;
and converting the relation of the s domain about the input variable and the output variable into a high-order differential equation form:
a n Y (n) +a n-1 Y (n-1) +L+a 1 Y′+a 0 Y＝b n U (n) +b n-1 U (n-1) +L+b 1 U′+b 0 U
wherein Y is (n) Represents the nth derivative of Y; in the same way, U (n) Represents the nth derivative of U; y 'and U' respectively represent corresponding first-order derivatives;
converting the differential equation into a state space model form by a variable substitution method, and enabling x to be 0 ＝Y；x 1 = Y', analogy to x n ＝Y (n) (ii) a In the same way, u 0 ＝U；u 1 = U'; analogize to u n ＝U (n) . The following relationship may be finally formed
Obtaining a state space model containing 2 x n state variables through the conversion, wherein the model consists of 2 x n differential equations and an algebraic equation; when considering a given model input u 0 ＝u set The state space model is a system of equations comprising 2 × (n + 1) variables, 2 × (n + 1) equations;
known based on model input, i.e. u 0 ＝u set In combination with algebraic equations in the state space model, the final transformation into a standard state space model can be as follows:
wherein A is sub Representing a state matrix, x, corresponding to a transfer function F(s) sub Representing state variables therein, corresponding to x in the state space model 0 :x n-1 And u 0 :u n-1 ；
The system frequency response model according to the constructed transfer function form comprises a plurality of transfer functions, and the transfer functions are written into a state space model only comprising state variables; performing row writing on nonlinear links such as dead zones, power saturation and the like in the established system frequency response model through an algebraic equation; finally, collecting state variables of all transfer functions and algebraic variables in all algebraic equations, and integrally converting the system frequency response model into a state space model form as follows:
wherein, due to beingThe system frequency model not only comprises linear links such as a transfer function and the like, but also comprises a nonlinear link, so that f represents a function corresponding to the state equation; g represents a function corresponding to an algebraic equation; the x table state variable matrix specifically comprises a system frequency difference value delta f and the like in a system frequency response model; y represents an algebraic variable matrix; u represents the input control variable matrix, e.g. the system active disturbance Δ P in the system frequency response model dist And the like.
Assigning initial values to state variables and algebraic variables in the state space model based on the initial operation state of the power grid and expected fault disturbance; the specific principle is as follows: giving an initial value of delta P to an input variable u (representing active disturbance) in a state space model dist Because the system initially runs in a stable state, the initial values of the other state variables are assigned to be 0, then the initial values of the corresponding algebraic variables are calculated according to the algebraic equation based on the initial values of the state variables, and finally the initial values of all the variables in the converted state space model can be obtained.
Further, calculating the disturbed system frequency comprises:
the total fluctuation time length of the given frequency to be calculated is recorded as t sum The numerical integration step is denoted as t bc ；
Based on the state space model and the initial variable value thereof, the specific process of calculating the system frequency by adopting the trapezoidal integral and Newton iteration method is as follows:
taking the value x according to the state variable at the current moment k (ii) a At the initial time, x k ＝x 0 (ii) a Output variable y k ，y k From x k Calculating to obtain the state variable x at the next moment by using the thought of trapezoidal integral k+1 The following are:
wherein: x in the equation k+1 Is a waiting amount, y k+1 Can be composed of x k+1 And (4) calculating. The trapezoidal integration method can reduce the error accumulation in the successive iteration process, so the trapezoidal integration is utilizedThe fractional mode has certain advantages for calculating the system frequency.
Solving variable x by adopting Newton iteration thought according to nonlinear equation system solving principle k+1 The method comprises the following steps:
wherein: i =0,1,2, \8230;
represents the result of solving the equation at time k using Newton's iteration after the ith iteration, an
Representing the iteration step length obtained by the i +1 th iteration calculation;
represents the error after the (i + 1) th iteration; err (r) th Indicating a set error threshold; when the iteration error is less than the error threshold, the iteration is considered to be converged, i.e., by x k Iterative computation to get x k+1 ；
Continuously carrying out iterative computation on the system frequency at each moment until the total time t of frequency stability analysis is reached sum The requirements of (2).
S3: constructing an inertia threshold value optimization model of the inertia power grid;
specifically, the method for constructing the inertia threshold optimization model of the inertia power grid comprises the following steps:
the objective function with the minimum integral inertia provided by the conventional unit and the new energy unit in the optimization model is specifically as follows:
it should be noted that, since the rotor inertia of the conventional unit itself remains unchanged under the same unit combination, the sum of the virtual inertias provided by the new energy station only needs to be guaranteed to be minimum in the objective function.
The variables to be decided of the optimization model are virtual inertia coefficients provided by each new energy station:
the lowest point of the system frequency and the change rate of the power grid after the generator N-1 fault disturbance occurs in the optimization model meet the safety requirement, and the related constraint of the virtual inertia coefficient of the new energy station in the specified range is as follows:
wherein,
representing a virtual inertia coefficient provided by the ith new energy station;
the capacity of the ith new energy station is represented;
represents the moment of inertia provided by the jth conventional station; n is a radical of re Representing the number of new energy stations; Δ f max Representing the maximum value of the system frequency deviation, and taking a positive number; Δ f rocof Representing the rate of change of the system frequency;
representing a system frequency maximum deviation limit;
representing a system frequency change rate limit;
representing the maximum limit value of the inertia of the new energy station;
and uniformly combining the objective function, the decision variable and the constraint condition to construct the inertia threshold optimization model of the inertia power grid.
It should be noted that, on the premise of keeping the inertia constant of the conventional unit unchanged, the inertia support capacity of the new energy station is considered in the frequency modulation model to optimize the system inertia threshold, so that the post-fault system frequency fluctuation characteristic under the condition of high occupancy of new energy in the novel power system can be better reflected.
S4: and solving the optimization model by adopting an improved GWO algorithm, and outputting a calculation result of the power grid inertia threshold.
Specifically, the inertia threshold optimization model is converted into a structure suitable for solving by a GWO algorithm, where a fitness function is specifically defined as follows:
wherein F (·) represents a fitness function; m represents a normal number 1 to 2 orders of magnitude greater than x; g (-) is an introduced piecewise function and represents a penalty term; when Δ f max And Δ f rocof After the limit is out of limit, the punishment item is a negative number, and the absolute value of the punishment item is rapidly increased due to the value of M; and when the requirement is met, the penalty term is 0;
in an optimization model structure suitable for a GWOO optimization algorithm, constraint conditions are as follows:
wherein,
representing the virtual inertia coefficient provided by the ith new energy station,
and representing the maximum limit value of the inertia of the new energy station.
Furthermore, the location updating strategy based on the wolf α in the Levy flight improvement GWO algorithm is as follows: according to the original GWO algorithm: selecting a wolf with the maximum fitness in a wolf group at the beginning of each iteration of the algorithm as a wolf alpha in the current iteration process, and correspondingly updating the position and the fitness of the alpha; on the basis, according to a Levy flight principle and in combination with the optimal position in the historical search process of the wolf pack, the position of the wolf alpha is randomly updated once, and the method specifically comprises the following steps:
wherein,
represents the alpha wolf position after Levy flight; levy (λ) represents a lavian random number, which can be generally simulated by using a Mantegna algorithm, and can be specifically expressed as follows in combination with GWO:
wherein,
σ v ＝1；
wherein N (-) represents a normal distribution; Γ (·) represents a gamma function; the lambda is generally in the range of 1 to 3;
initializing a wolf pack position according to a GWO algorithm execution flow, calculating the individual fitness of the wolf pack, selecting a leading wolf according to the fitness calculation result, updating a position coordinate, randomly updating an alpha wolf position according to Levy flight, updating the position of each individual in the wolf pack according to the distance relation between the leading wolf and the wolf pack individual, checking whether the wolf pack individual coordinate exceeds a parameter value range, correcting the coordinate, and entering next iteration until the maximum iteration number is reached, wherein the specific process comprises the following steps:
initializing a wolf pack position: first the number of wolf clusters and the maximum number of iterations are given. Variable to be optimized
The number of the wolf clusters corresponds to the dimensionality of the position coordinates of the wolf clusters, and the initial coordinates of individuals in the wolf clusters are randomly or self-defined according to the value range of the variables to be optimized;
calculating the fitness of individuals in the wolf group: adding a penalty function corresponding to the constraint condition into the target function as the fitness of each wolf head in each iteration process;
updating the positions of the wolf α, β, δ: selecting three wolfs with higher fitness from all individuals of the wolf group, sequentially using the three wolfs as alpha, beta and delta wolfs in the iteration process, and correspondingly updating the position information;
randomly updating the position of the wolf alpha at one time according to Levy flight;
updating the position of the wolf group individual: depending on the distance between the wolf group individual and the wolf alpha, beta and delta, the position of each wolf group individual is updated as follows:
where, represents a vector dot product; a and C represent random coefficients, A =2a × r 1 -a，C＝2r 2 (ii) a a represents the step length, and can be generally set to be 0.01 times of the parameter value range;
and
representing the position coordinates of the k-th wolf in the wolf group in the processes of the l-th iteration and the l +1 iteration;
represents the position of the alpha wolf in the first iteration process; d α Representing a distance variable between the alpha wolf and a wolf group individual;
in obtaining
And then, judging whether the position meets the range requirement of the parameter to be decided, if not, pulling the wolf position back to the boundary of the range, and then entering next iteration until the maximum iteration number is reached.
Specifically, the power grid inertia threshold calculation process includes:
wherein,
representing a required inertia threshold of the power grid;
representing the inertia of the conventional unit;
representing the inertia of the new energy source unit obtained by optimization;
and
respectively representing the capacities of the conventional unit and the new energy station.
It should be noted that by introducing a Levy flight process in a traditional GWO algorithm, random position update is added to the Alpha wolf position based on the historical optimal position in each iteration optimization process, so that the search range of the wolf colony in the optimization process can be improved, the probability of convergence of the wolf colony at the local optimal point is reduced, and the quality of model solution is improved.
Example 2
Referring to fig. 2 to 7, a simulation example of a power grid rotational inertia threshold calculation method for improving a GWO algorithm is provided for an embodiment of the present invention, and in order to verify the beneficial effects of the present invention, scientific demonstration is performed through simulation experiments.
The embodiment comprises an alternating current power grid, a conventional unit, a new energy station and a corresponding controller; the conventional unit and the new energy station are both connected to an alternating current power grid, virtual inertia support parameters in a frequency modulation controller of the new energy station can be regulated, and a regulation reference value is a calculated result of the method.
Specifically, a specific implementation process of the power grid inertia threshold calculation method based on the improved GWO algorithm is explained by taking fig. 2 as an example. Fig. 2 is a modified IEEE 10 machine 39 node system for explaining the implementation process of the inertia threshold calculation method provided by the present invention, the system includes 7 synchronous machine sets G1, G2, G4, G5, G6, G7, G8, wherein the unit power adjustment coefficients of all the generator primary frequency modulation controllers are 25;3 new energy stations PV1, PV2 and PV3, wherein the photovoltaic unit is provided with a virtual inertia control and primary frequency modulation control system. The method mainly aims at the N-1 fault of a generator G5 in a figure, minimizes the virtual inertia coefficients of frequency modulation controllers PV1, PV2 and PV3 of 3 new energy stations on the premise of ensuring that the lowest point of system frequency and the frequency change rate meet the safety requirements after the fault, and finally calculates the inertia threshold value required by the system by combining the inertia constants of the conventional unit under the current unit combination.
The method adopts the following technical scheme:
monitoring the capacities of G1-G10 conventional units and new energy stations, reserved spare and frequency modulation controller structures and other information, equating the whole system to be a synchronous unit, establishing a system frequency response model according to the monitoring information and expected active disturbance, and calculating a disturbed system frequency fluctuation curve by a numerical integration method; according to the provided system frequency calculation method, a system inertia threshold value optimization model which takes the minimum virtual inertia of the whole PV1, PV2 and PV3 of the new energy station as a target, meets the requirements with the lowest point of the system frequency and the frequency change rate after the N-1 fault and takes the virtual inertia coefficient of each station in the required range as a constraint condition is established; and then, superposing the constraint condition with an objective function in a penalty function mode, converting the optimization model into a structure which is suitable for being solved by a GWOO algorithm, finally solving the optimization model according to the GWOO algorithm based on Levy flight improvement to obtain the minimum inertia required by the new energy stations PV1, PV2 and PV3, and combining the inertia threshold required by a conventional unit inertia calculation system.
The method specifically comprises the following steps:
s1: respectively establishing frequency modulation controller numerical models based on frequency feedback of all stations by monitoring station capacity, reserved standby, frequency modulation controller structure and other information of conventional units G1, G2, G4, G5, G6, G7 and G8 and new energy stations PV1, PV2 and PV 3; then, combining a rotor motion equation of a synchronous generator, equating the whole power grid to be a synchronous machine, and further establishing a system frequency response model considering frequency modulation dead zones, power regulation limits and generator power regulation dynamics;
specifically, a frequency feedback frequency modulation model related to a conventional unit is constructed according to monitored information of capacities, reserved reserves, structures of frequency modulation controllers and the like of conventional unit fields G1, G2, G4, G5, G6, G7 and G8, the specific structure mainly has differences between hydroelectric power and thermal power, mathematical models corresponding to structures of the hydroelectric power and hydroelectric power units comprise a frequency modulation dead zone link, a speed regulator and prime motor link and a power regulation saturation link, and the specific structure is shown in FIG. 3. Wherein in the power gridIn the middle, the speed governor parameter of the same type unit is basically similar, wherein: t is a unit of G =0.5s is the time constant of the speed regulator, T CH =7s is the turbine time constant, T RH =7s is a reheat constant, F HP =0.3 is reheat coefficient, T R =3.3s is the reset time constant, R T Temporary reduction rate, R, =0.16 T =0.025 is the permanent reduction rate, unit power regulation factor
30, 40 respectively.
According to the monitored information of the capacity, the reserved spare and the frequency modulation controller structure of the new energy station and the like, the frequency feedback frequency modulation model related to the new energy station is constructed, and the frequency feedback frequency modulation model also comprises the following steps: the method comprises the following steps of a frequency modulation dead zone, a filter, a primary frequency modulation and inertia supporting link and a power saturation link, wherein the primary frequency modulation and inertia supporting link is greatly different from a conventional unit, and specifically shown in figure 4, wherein a virtual inertia coefficient H in a frequency modulation controller of a new energy station re And optimizing parameters to be optimized of the model for the subsequent inertia threshold value.
Furthermore, the specific form of the rotor motion equation of the conventional units G1, G2, G4, G5, G6, G7, G8 is as follows:
wherein H g Represents an inertia constant of the generator; delta P gm Representing the variation of the mechanical power of the generator; delta P ge Representing the variation of electromagnetic power at the terminal; d g Representing the damping coefficient of the generator; Δ f g Representing the generator rotor frequency deviation.
The whole power grid (mainly comprising the conventional generators G1, G2, G4, G5, G6, G7 and G8) is equivalent to a synchronous generator based on a rotor motion equation, the frequency dynamics of the equivalent generator is similar to that of the rotor motion equation, and the method comprises the following steps:
wherein, Δ P sm Representing the mechanical power of the equivalent generator; delta P se Representing the generator terminal electromagnetic power of the equivalent generator; d s Representing the damping coefficient of the equivalent generator; h s Representing the inertia constant of the equivalent generator.
The method comprises the following steps of constructing a corresponding system frequency response model based on a mathematical model of a frequency modulation controller of the conventional unit and the new energy station and a rotor motion equation of a system equivalent generator, and specifically comprising the following steps:
(1) Calculating the inertia constant of the equivalent generator of the system through the inertia constant and the capacity of each conventional unit G1, G2, G4, G5, G6, G7 and G8:
wherein,
and
respectively representing the inertia constant and the capacity of the ith conventional unit; n is a radical of sg Representing the total number of conventional units, N sg ＝7。
(2) Based on the frequency modulation link based on frequency feedback of the conventional units such as hydropower, thermal power and the like, the actual power of each station is adjusted by the quantity delta P g Summing and corresponding system equivalent generator rotor motion equation mechanical power variation, namely:
(3) Based on the frequency modulation link of the new energy station based on the frequency feedback controller, the actual power regulating quantity delta P of each station is adjusted re And in this case N-1 failure of generator G5Summing the power disturbance quantity, and corresponding to the variation of the electromagnetic power at the equivalent generator end of the system, namely:
wherein, Δ P G5 Representing the active disturbance quantity corresponding to the N-1 fault of the generator G5; n is a radical of re Indicates the total number of new energy stations with frequency modulation control function, in this example N re ＝3。
Finally, combining the above transfer functions and algebraic equations, a system frequency response model in the form of a transfer function can be obtained, as shown in fig. 5.
It should be noted that through the above steps, complex characteristics such as power adjustment dynamics, frequency modulation dead zone, power saturation and the like of the conventional unit and the new energy unit are considered in the constructed system frequency response model, and a more accurate basic model can be provided for subsequent frequency calculation.
S2: converting the established system frequency response model in the form of the transfer function into a standard state space model form; then, based on the initial operation state of the power grid and the expected fault disturbance (corresponding to the N-1 fault of the generator G5 in the embodiment), initial values are given to the state variable and the algebraic variable in the obtained state space model; then based on the initial variable value, setting the required calculated frequency fluctuation duration and the numerical integration step length, and calculating the system frequency fluctuation curve after the power grid is disturbed by adopting a trapezoidal integration and Newton iteration method;
the transfer function form system frequency response model established in S1 is converted into a standard state space model form, where a rotor motion equation is taken as an example (corresponding to the rotor motion equation transfer function in fig. 5), specifically as follows
(1) And (3) performing formal arrangement on the transfer function model according to the relation between input and output:
wherein F(s) represents a transfer function of order 1; Δ f and Δ P represent output and input variables corresponding to the transfer function model, respectively; s represents the laplacian operator.
(2) And converting the relation of the s domain about the input variable and the output variable into a high-order differential equation form:
2H s Δf′+D s Δf＝ΔP
where Δ f' represents the first derivative to Δ f.
(3) Converting the differential equation into a state space model form by a variable substitution method, and enabling x 0 ＝Δf；x 1 ＝Δf′，u 0 = Δ P; the following relationship may be finally formed
(4) Based on model input known, i.e. u 0 For known variables, combining algebraic equations in the state-space model above may eliminate algebraic variables in the model above, leaving only state variables, which may eventually be transformed into the form of a standard state-space model, as follows:
(5) The system frequency response model according to the constructed transfer function form comprises a plurality of transfer functions, is converted according to the processes of (1) to (4) in the same way, and is written into a state space model only comprising state variables; writing the nonlinear links such as dead zones, power saturation and the like in the established system frequency response model through an algebraic equation; finally, the state variables of all transfer functions and the algebraic variables in all algebraic equations are collected, and the system frequency response model can be integrally converted into a state space model form as follows:
the system frequency model not only comprises linear links such as a transfer function and the like, but also comprises a nonlinear link, so that f represents a function corresponding to the state equation; g represents a function corresponding to an algebraic equation; the x table state variable matrix specifically comprises a system frequency difference value delta f and the like in a system frequency response model; y represents an algebraic variable matrix; u represents the input control variable matrix, e.g. the system active disturbance Δ P in the system frequency response model dist And the like.
Based on the initial operation state of the power grid and the expected fault disturbance, the specific principle of assigning initial values to the state variable and the algebraic variable in the obtained state space model is as follows: giving an initial value of delta P to an input variable u (representing active disturbance) in a state space model dist Because the system initially runs in a stable state, the initial values of the other state variables are assigned to be 0, then the initial values of the corresponding algebraic variables are calculated according to the algebraic equation based on the initial values of the state variables, and finally the initial values of all variables of the state space model can be obtained.
The total time duration of the fluctuation of the given frequency to be calculated is denoted t sum And the numerical integration step is denoted as t bc For this example, t sum ＝10s，t bc ＝0.01s。
Based on the state space model and the initial variable value thereof, the specific process of calculating the system frequency by adopting the trapezoidal integral and Newton iteration method is as follows:
(1) Taking the value x according to the state variable at the current moment k (at the initial time, x k ＝x 0 ) And an output variable y k (by x) k Calculated), the state variable x at the next moment is calculated by utilizing the thought of trapezoidal integral k+1 The following are:
wherein: in the equation x k+1 As a quantity to be sought, y k+1 Can be composed of x k+1 And (4) calculating. The trapezoidal integration method can reduce the gradual iteration processThe error accumulation in the method is carried out, so that the calculation of the system frequency by using the trapezoidal integration mode has certain advantages.
(2) Solving the variable x in the equation by adopting the Newton iteration thought according to the solving principle of the nonlinear equation system k+1 The method comprises the following steps:
wherein: i =0,1,2, \ 8230;
representing the result of solving the equation at time k using Newton iterations after the ith iteration, and
representing the iteration step length obtained by the i +1 th iteration calculation;
representing the error after the (i + 1) th iteration; err (r) th Indicating a set error threshold; when the iteration error is less than the error threshold, the iteration is considered to be converged, i.e., by x k Iterative computation to get x k+1 。
(3) According to (1) and (2), continuously carrying out iterative calculation on the system frequency at each moment until the total time length t of frequency stability analysis is reached sum Requirement of =10 s.
S3: based on a system frequency calculation result, the minimum integral inertia provided by a conventional unit and a new energy unit is taken as a target, the lowest point of system frequency fluctuation and the frequency change rate after the power grid has a generator N-1 fault disturbance meet the safety requirement of the power grid, and the value range of the virtual inertia coefficient of the new energy station is taken as a constraint condition, a system inertia threshold value optimization model is constructed, and the minimum virtual inertia parameter required by each new energy station frequency modulation controller is calculated.
The objective function with the minimum integral inertia provided by the conventional unit and the new energy unit in the optimization model is specifically as follows;
wherein N is re And =3. It should be noted that, since the rotor inertia of the conventional unit itself remains unchanged under the same unit combination, the sum of the virtual inertias provided by the new energy station only needs to be guaranteed to be minimum in the objective function.
And optimizing the virtual inertia coefficient provided by each new energy station by the to-be-decided variable of the model:
and
the lowest point of the system frequency and the change rate of the generator G5 subjected to the N-1 fault disturbance in the optimization model meet the safety requirements, and the related constraints of the virtual inertia coefficient of the new energy station in the specified range are as follows:
wherein,
representing a virtual inertia coefficient provided by the ith new energy station;
representing the capacity of the ith new energy station;
represents the moment of inertia provided by the jth conventional station; n is a radical of hydrogen re Representing the number of new energy stations; Δ f max Represents the maximum value (here, a positive number) of the systematic frequency deviation; Δ f rocof Representing the rate of change of the system frequency;
representing a system frequency maximum deviation limit;
representing a system frequency change rate limit;
and representing the maximum limit value of the inertia of the new energy station.
Is 1 (per unit value), the corresponding virtual inertia coefficient reference value is 50; wherein the system frequency minimum and the constraint of the change rate can be calculated through S1 and S2.
And (3) constructing a system inertia threshold optimization model, namely combining the objective function, the decision variables and the constraint conditions. On the premise of keeping the inertia constant of a conventional unit unchanged, the inertia supporting capacity of the new energy station is considered in the frequency modulation model to optimize the system inertia threshold, so that the system frequency fluctuation characteristic after the fault under the condition of high new energy occupancy in the novel power system can be better reflected.
S4: according to the constructed power grid inertia threshold calculation model and a gray wolf optimization (GWOO) algorithm principle, defining a fitness function as a penalty function of adding safety constraints on a frequency lowest point and a change rate to an objective function in the optimization model, only reserving a value range of a variable to be decided as a constraint condition, and converting the inertia threshold calculation model in S3 into a structure suitable for gray wolf optimization; and then, based on a Levy flight principle, improving a wolf pack position updating strategy in the GWOO model, and improving the convergence and solving accuracy of the model. Finally, according to GWOO execution flow, initializing a wolf group position, calculating the fitness of wolf group individuals, selecting a leading wolf in the wolf group according to a fitness function, updating position coordinates, randomly updating the position of part of the leading wolf according to levy flight, updating the position coordinates of the wolf group individuals according to the distance between the wolf group individuals and the leading wolf, checking whether the wolf group position coordinates meet the parameter range requirement or not, entering next iteration, and finally calculating the inertia threshold value required by the whole power grid by combining inertia parameters of a conventional unit under the current unit combination, wherein the convergence point position coordinates correspond to the lowest virtual inertia parameters required by each new energy field station after the wolf group reaches the iteration times.
Specifically, a fitness function is defined as a penalty function of adding a frequency minimum point and frequency change rate safety constraint to an objective function of an S3 constructed inertia threshold optimization model by combining with a GWO algorithm principle, only a value range of a variable to be decided is reserved as a constraint condition, and the inertia threshold optimization model of the S3 is converted into a structure suitable for being solved by the GWO algorithm, wherein the definition of the fitness function is specifically as follows:
wherein F (-) represents a fitness function; m represents a sufficiently large normal number, which in this example takes the value 100; g (-) is an introduced piecewise function, representing a penalty term, i.e. when Δ f max And Δ f rocof After the limit is out of limit, the punishment item is a negative number, and the absolute value of the punishment item is rapidly increased due to the value of M; and when the requirement is met, the penalty term is 0.
In an optimization model structure suitable for a GWOO optimization algorithm, constraint conditions are as follows:
the location updating strategy based on the wolf α in the Levy flight improvement GWO algorithm in S4 specifically includes: according to the original GWOO algorithm, the wolf with the maximum fitness in a wolf group is selected at the beginning of each iteration of the algorithm and is used as the wolf head alpha in the current iteration process, and the position and the fitness of the alpha are correspondingly updated. On the basis, according to a Levy flight principle and in combination with the optimal position in the historical search process of the wolf pack, the position of the wolf head alpha is randomly updated once, and the method specifically comprises the following steps:
wherein a is a =0.01;
represents the alpha wolf position after Levy flight; in this example, λ is 1.5, levy (1.5) represents a lewy random number, which can be generally simulated by using Mantegna algorithm, and can be specifically expressed as follows in combination with GWO:
wherein,
σ v ＝1；
wherein N (-) represents a normal distribution; Γ (·) represents a gamma function.
It should be noted that random position updating is added to the Alpha wolf position once based on the historical optimal position in each iteration optimization process by introducing the Levy flight process in the traditional GWOO algorithm, so that the search range of the wolf colony in the optimization process can be improved, the probability of convergence of the wolf colony at the local optimal point is reduced, and the quality of model solution is improved.
The method comprises the following steps of S4, executing a flow according to a GWO algorithm, initializing a wolf group position, calculating individual fitness of the wolf group, selecting a leading wolf according to a fitness calculation result, updating a position coordinate, randomly updating an alpha wolf position according to Levy flight, updating the position of each individual in the wolf group according to the distance relationship between the leading wolf and the wolf group individual, checking whether the wolf group individual coordinate exceeds a parameter value range, correcting the coordinate, entering next iteration until the maximum iteration frequency is reached, and the specific process is as follows:
(1) Initializing the position of a wolf group: firstly, the number of wolf groups is given to be 50; the maximum number of iterations is 100. Variables to be optimized: (
And
) The number of the binary coordinates is corresponding to the dimensionality of the position coordinates of the wolf pack, and the initial coordinates of the individuals in the wolf pack are randomly or self-defined generated according to the value range (0-1) of the variable to be optimized.
(2) Calculating the fitness of individuals in the wolf group: adding a penalty function corresponding to the constraint condition into the target function as the fitness of each wolf in each iteration process;
(3) Updating the positions of the wolf α, β, δ: selecting three wolfs with higher fitness from all individuals of the wolf group, sequentially using the three wolfs as alpha, beta and delta wolfs in the iteration process, and correspondingly updating the position information;
(4) Randomly updating the position of the wolf alpha at one time according to Levy flight;
(5) Updating the position of the wolf group individual: depending on the distance between the wolf group individual and the wolf alpha, beta and delta, the position of each wolf group individual is updated as follows:
where, represents a vector dot product; a and C represent random coefficients, A =2a × r 1 -a，C＝2r 2 (ii) a a represents the step length, and can be generally set to be 0.01 times of the parameter value range;
and
representing the position coordinates of the k-th wolf in the wolf group in the processes of the l-th iteration and the l +1 iteration;
represents the position of the alpha wolf in the first iteration process; d α Represents the distance variable between the alpha wolf and the wolf group individual.
(5) In obtaining
And then, judging whether the position is between 0 and 1, if not, pulling the wolf position back to the boundary of the range, and then entering next iteration according to the flow until the maximum iteration number is reached.
The flow of the model solution based in part on the improved GWO algorithm is shown in fig. 6.
Based on this example, in the optimization process, different initial position ranges of the wolf pack are given, the final result is displayed under different initial positions, and the wolf packs are finally converged to the vicinity of the same position, which illustrates that by introducing the Levy flight process in the conventional GWO algorithm, the model solving quality can be improved, and the probability of convergence of the wolf pack at the local optimum point in the initial stage is reduced, and the specific result is shown in fig. 7.
The calculation process of the inertia threshold value required by the whole power grid is as follows:
wherein,
representing a required inertia threshold of the power grid;
representing the inertia of the conventional unit;
representing the inertia of the new energy unit obtained by optimization;
and
respectively representing the capacities of the conventional unit and the new energy station.
The simulation experiment fully illustrates the feasibility and the beneficial effects of the method, the minimum virtual inertia coefficient of each new energy station is obtained, the inertia threshold value required by the power grid in the current operation scene is calculated by combining the constant inertia constant of the conventional unit in the current unit combination mode, and a reference basis is provided for the safe and stable control of the power grid frequency.
It should be noted that the above-mentioned embodiments are only for illustrating the technical solutions of the present invention and not for limiting, and although the present invention has been described in detail with reference to the preferred embodiments, it should be understood by those skilled in the art that modifications or equivalent substitutions may be made on the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, which should be covered by the claims of the present invention.
Claims (10)
1. A power grid rotational inertia threshold calculation method for improving a GWO algorithm is characterized by comprising the following steps:
monitoring power grid equipment and operating parameters, and constructing a system frequency response model;
converting the model into a standard state space form, and calculating the frequency of the disturbed system;
constructing an inertia threshold value optimization model of the inertia power grid;
and solving the optimization model by adopting an improved GWO algorithm, and outputting a calculation result of the power grid inertia threshold.
2. The method for calculating the grid moment of inertia threshold value for improving the GWO algorithm as claimed in claim 1, wherein the grid system to which the method for calculating the inertia threshold value is applied generally comprises: the system comprises an alternating current power grid, a conventional unit, a new energy station and a corresponding frequency modulation controller.
3. The method for calculating the grid moment of inertia threshold value for improving the GHO algorithm as claimed in claim 2, comprising: the conventional unit is generally provided with a speed regulator, the new energy station part is provided with a frequency modulation controller and comprises a virtual inertia support system, and the conventional unit and the new energy station belong to a power generation unit in an alternating current power grid; and calculating a reference value of an inertia parameter in the virtual inertia supporting system of the new energy station by the inertia threshold optimization method.
4. The method for calculating the grid moment of inertia threshold for improving the GWO algorithm as claimed in claim 3, wherein the constructing the system frequency response model comprises:
constructing a frequency feedback frequency modulation model related to the conventional unit according to the monitored information such as the capacity, the reserved spare, the frequency modulation controller structure and the like of the conventional unit station; the difference mainly exists on water and electricity and thermoelectricity in concrete structure, wherein to thermoelectricity unit, the mathematical model that its structure corresponds includes frequency modulation blind spot link, speed regulator and prime mover link, the power regulation saturation link, specifically as follows:
wherein Δ f represents a system frequency deviation; Δ f act Representing the amount of frequency deviation across the frequency modulation dead zone; delta P g.ref Representing an active regulation reference value based on frequency feedback; delta P g Representing an actual power adjustment; delta P g.up And Δ P g.down Respectively representing the upper limit and the lower limit of active power adjustment; t is a unit of G Represents the governor time constant; t is CH Representing the turbine time constant; t is RH Represents a reheat time constant; f HP Represents a reheat coefficient; k is a radical of p Expressing a unit power adjustment coefficient;
a mathematical model corresponding to a frequency modulation controller of the hydroelectric generating set is constructed, the mathematical model also comprises a frequency modulation dead zone link, a speed regulator and prime motor link and a power regulation saturation link, and the mathematical model is mainly different from a thermal power generating set in that a transfer function corresponding to the speed regulator and the prime motor link is as follows:
wherein R is T Is a temporary rate of decline; r P Is the permanent rate of decline;
according to the monitored information of the capacity, the reserved spare and the frequency modulation controller structure of the new energy station and the like, a frequency feedback frequency modulation model related to the new energy station is constructed, and the method comprises the following steps: the system comprises a frequency modulation dead zone, a filter, a primary frequency modulation and inertia supporting link and a power saturation link, wherein the primary frequency modulation and inertia supporting link is greatly different from a conventional unit, and the system specifically comprises the following parts:
wherein, Δ P re.ref Indicating a power regulation reference value of the new energy station participating in frequency modulation; k re And H re Respectively representing a frequency droop coefficient and a virtual inertia coefficient of a frequency modulation controller of the new energy station;
the rotor equation of motion is embodied as follows:
wherein H g Representing the inertia constant of the generator; delta P gm Representing the variation of the mechanical power of the generator; delta P ge Representing the variation of electromagnetic power at the terminal; d g Representing the damping coefficient of the generator; Δ f g Representing generator rotor frequency deviation;
equating the whole power grid into a synchronous generator based on a rotor motion equation, wherein the frequency dynamic of the equivalent generator is similar to the rotor motion equation, and the method comprises the following steps:
wherein, Δ P sm Representing the mechanical power of the equivalent generator; delta P se Representing the generator terminal electromagnetic power of the equivalent generator; d s Representing the damping coefficient of the equivalent generator; h s Represents the inertia constant of the constant value generator;
and constructing a corresponding system frequency response model based on the mathematical models of the frequency modulation controllers of the conventional unit and the new energy station and the rotor motion equation of the system equivalent generator, wherein the corresponding system frequency response model is specifically as follows:
calculating the inertia constant of the equivalent generator of the system through the inertia constant and the capacity of each conventional unit:
wherein,
and
respectively representing the inertia constant and the capacity of the ith conventional unit; n is a radical of sg Representing the total number of conventional units;
based on the frequency modulation link based on frequency feedback of conventional units such as hydropower, thermal power and the like, the actual power regulating quantity delta P of each station is adjusted g Summing and corresponding system equivalent generator rotor motion equation mechanical power variation, namely:
adjusting the actual power of each station by delta P based on the frequency modulation link of the new energy station based on the frequency feedback controller re Summing the predicted power disturbance quantity of the N-1 fault, and corresponding to the variation quantity of the electromagnetic power at the equivalent generator end of the system, namely:
wherein, Δ P dis Representing the active disturbance quantity corresponding to the expected N-1 fault; n is a radical of re Representing the total number of new energy stations with frequency modulation control functionality.
5. The method for calculating the grid moment of inertia threshold for improving the GWO algorithm as claimed in claim 4, wherein converting the system frequency response model into a standard state space form comprises:
and (3) performing formal arrangement on the transfer function model according to the relation between input and output:
wherein F(s) represents an nth order transfer function; y and U respectively represent output and input variables corresponding to the transfer function model; a is 0 ～a n And b 0 ～b n Coefficients representing terms of respective orders in the transfer function; s represents the laplacian operator;
and (3) converting the relation of the s domain about the input variable and the output variable into a high-order differential equation form:
a n Y (n) +a n-1 Y (n-1) +L+a 1 Y′+a 0 Y＝b n U (n) +b n-1 U (n-1) +L+b 1 U′+b 0 U
wherein Y is (n) Represents the nth derivative of Y; in the same way, U (n) Represents the nth derivative of U; y 'and U' respectively represent corresponding first-order derivatives;
converting the differential equation into a state space model form by a variable substitution method, and enabling x 0 ＝Y；x 1 = Y', analogy to x n ＝Y (n) (ii) a In the same way, u 0 ＝U；u 1 = U'; analogize to u n ＝U (n) (ii) a The following relationship may be finally formed
Obtaining a state space model containing 2 x n state variables after conversion, wherein the model consists of 2 x n differential equations and an algebraic equation; when considering a given model input u 0 ＝u set The state space model is a system of equations comprising 2 × (n + 1) variables, 2 × (n + 1) equations;
based on model input known, i.e. u 0 ＝u set The algebraic equations in the state-space model, in combination, may ultimately be transformed into the form of a standard state-space model, as follows:
wherein, A sub Representing a state matrix, x, corresponding to a transfer function F(s) sub Representing state variables therein, corresponding to x in the state space model 0 :x n-1 And u 0 :u n-1 ；
The system frequency response model in the form of the constructed transfer function comprises a plurality of transfer functions, and the transfer functions are written into a state space model only comprising state variables; performing row writing on nonlinear links such as dead zones, power saturation and the like in the established system frequency response model through an algebraic equation; finally, collecting state variables of all transfer functions and algebraic variables in all algebraic equations, and integrally converting the system frequency response model into a state space model form as follows:
wherein f represents a function corresponding to the state equation; g represents a function corresponding to an algebraic equation; x is a state variable matrix; y represents an algebraic variable matrix; u represents an input control variable matrix;
assigning initial values to state variables and algebraic variables in the state space model based on the initial operation state of the power grid and expected fault disturbance; the specific principle is as follows: giving an initial value of Δ P to an input variable u in a state space model dist And because the system initially operates in a stable state, the initial values of the other state variables are assigned to be 0, then the initial values of the corresponding algebraic variables are calculated according to the algebraic equation based on the initial values of the state variables, and finally the initial values of all the variables in the converted state space model can be obtained.
6. The method for calculating grid moment of inertia threshold for improving GHO algorithm of claim 5, wherein the calculating the disturbed system frequency comprises:
the total fluctuation time length of the given frequency to be calculated is recorded as t sum And the numerical integration step is denoted as t bc ；
Based on the state space model and the initial variable value thereof, the specific process of calculating the system frequency by adopting the trapezoidal integral and Newton iteration method is as follows:
taking value x according to the state variable at the current moment k (ii) a At the initial time, x k ＝x 0 (ii) a Output variable y k ，y k From x k Calculating to obtain the state variable x at the next moment by using the thought of trapezoidal integration k+1 The following are:
wherein: in the equation x k+1 Is a waiting amount, y k+1 Can be composed of x k+1 Calculating to obtain;
solving variable x by adopting Newton iteration thought according to the solving principle of the nonlinear equation set k+1 The method comprises the following steps:
wherein: i =0,1,2, \ 8230;
representing the result of solving the equation at time k using Newton iterations after the ith iteration, and
representing the iteration step length obtained by the i +1 th iteration calculation;
representing the error after the (i + 1) th iteration; err (r) th Indicating a set error threshold; when the iteration error is less than the error threshold, the iteration is considered to be converged, i.e., by x k Iterative computation yields x k+1 ；
Continuously carrying out iterative computation on the system frequency at each moment until the total time t of frequency stability analysis is reached sum The requirements of (1).
7. The method for calculating the grid moment of inertia threshold value for improving the GHO algorithm of claim 6, wherein the constructing the model for optimizing the inertia grid inertia threshold value comprises:
the objective function with the minimum overall inertia provided by the conventional unit and the new energy unit in the optimization model is specifically as follows:
the variables to be decided of the optimization model are virtual inertia coefficients provided by each new energy station:
the lowest point of the system frequency and the change rate of the power grid after the generator N-1 fault disturbance occurs in the optimization model meet the safety requirement, and the related constraint of the virtual inertia coefficient of the new energy station in the specified range is as follows:
wherein,
representing a virtual inertia coefficient provided by the ith new energy station;
representing the capacity of the ith new energy station;
represents the moment of inertia provided by the jth conventional station; n is a radical of re Representing the number of new energy stations; Δ f max Representing the maximum value of the system frequency deviation, and taking a positive number; Δ f rocof Representing the rate of change of the system frequency;
representing a system frequency maximum deviation limit;
representing a system frequency change rate limit;
representing the maximum limit value of the inertia of the new energy station;
and uniformly combining the objective function, the decision variables and the constraint conditions to construct the inertia threshold optimization model of the inertia power grid.
8. The method for calculating the grid moment of inertia threshold for improving the GWO algorithm as claimed in claim 7, wherein the method comprises:
and converting the inertia threshold optimization model into a structure suitable for solving by a GWOO algorithm, wherein the definition of a fitness function is specifically as follows:
wherein F (·) represents a fitness function; m represents a normal number 1 to 2 orders of magnitude greater than x; g (-) is an introduced piecewise function and represents a penalty term; when Δ f max And Δ f rocof After the limit is out of limit, the punishment item is a negative number, and the absolute value of the punishment item is rapidly increased due to the value of M; and when the requirement is met, the penalty term is 0;
in an optimization model structure suitable for a GWOO optimization algorithm, constraint conditions are as follows:
wherein,
representing the virtual inertia coefficient provided by the ith new energy station,
and representing the maximum limit value of the inertia of the new energy station.
9. The method for calculating the grid moment of inertia threshold for improving the GWO algorithm as claimed in claim 8, comprising: the location updating strategy based on the wolf α in the Levy flight improvement GWO algorithm specifically comprises the following steps:
original GWO algorithm: selecting a wolf with the maximum fitness in a wolf group at the beginning of each iteration of the algorithm as a wolf alpha in the current iteration process, and correspondingly updating the position and the fitness of the alpha;
on the basis, according to a Levy flight principle and in combination with the optimal position in the historical search process of the wolf pack, the position of the wolf head alpha is randomly updated once, and the method specifically comprises the following steps:
wherein,
represents the alpha wolf position after Levy flight; levy (λ) represents a lavian random number, which can be generally simulated by using a Mantegna algorithm, and can be specifically expressed as follows in combination with GWO:
wherein,
σ v ＝1；
wherein N (·) represents a normal distribution; Γ () represents a gamma function; the lambda generally ranges from 1 to 3;
the specific process of improving the GWO algorithm execution flow is as follows:
initializing the position of a wolf group: firstly, the number of wolf groups and the maximum iteration number are given; variable H to be optimized r i e The number of the wolf clusters corresponds to the dimensionality of the position coordinates of the wolf clusters, and the initial coordinates of individuals in the wolf clusters are generated randomly or in a self-defined mode according to the value range of the variables to be optimized;
calculating the fitness of individuals in the wolf group: adding a penalty function corresponding to the constraint condition into the target function as the fitness of each wolf head in each iteration process;
updating the positions of the wolf alpha, beta and delta: selecting three wolfs with higher fitness from all individuals of the wolf group, sequentially using the three wolfs as alpha, beta and delta wolfs in the iteration process, and correspondingly updating the position information;
randomly updating the position of the wolf alpha once according to Levy flight;
updating the position of the wolf group individual: depending on the distance between the wolf pack individual and the wolf α, β, δ, the position of each wolf pack individual is updated as follows:
where, represents a vector dot product; a and C represent random coefficients, A =2a × r 1 -a，C＝2r 2 (ii) a a represents the step length, and can be generally set to be 0.01 times of the parameter value range;
and
representing the position coordinates of the k-th head wolf in the wolf group in the l and l +1 iteration processes;
represents the position of the alpha wolf in the first iteration process; d α Representing a distance variable between the alpha wolf and the wolf group individual;
is obtained by
And then, judging whether the position meets the range requirement of the parameter to be decided, if not, pulling the wolf position back to the boundary of the range, and then entering next iteration until the maximum iteration number is reached.
10. The method for calculating grid inertia threshold for improving GWO algorithm of claim 9, wherein the grid inertia threshold calculation process comprises:
wherein,
representing a required inertia threshold of the power grid;
representing the inertia of the conventional unit;
representing the inertia of the new energy unit obtained by optimization;
and
and respectively representing the capacities of the conventional unit and the new energy station.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210904002.8A CN115173410A (en) | 2022-07-29 | 2022-07-29 | Power grid rotational inertia threshold calculation method for improving GWO algorithm |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210904002.8A CN115173410A (en) | 2022-07-29 | 2022-07-29 | Power grid rotational inertia threshold calculation method for improving GWO algorithm |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN115173410A true CN115173410A (en) | 2022-10-11 |
Family
ID=83476968
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210904002.8A Pending CN115173410A (en) | 2022-07-29 | 2022-07-29 | Power grid rotational inertia threshold calculation method for improving GWO algorithm |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115173410A (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN118052153A (en) * | 2024-04-16 | 2024-05-17 | 上海叁零肆零科技有限公司 | Natural gas pipe network data optimization method, storage medium and electronic equipment |
-
2022
- 2022-07-29 CN CN202210904002.8A patent/CN115173410A/en active Pending
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN118052153A (en) * | 2024-04-16 | 2024-05-17 | 上海叁零肆零科技有限公司 | Natural gas pipe network data optimization method, storage medium and electronic equipment |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Oshnoei et al. | Novel load frequency control scheme for an interconnected two-area power system including wind turbine generation and redox flow battery | |
| Saez-de-Ibarra et al. | Co-optimization of storage system sizing and control strategy for intelligent photovoltaic power plants market integration | |
| Pan et al. | Fractional order fuzzy control of hybrid power system with renewable generation using chaotic PSO | |
| CN106786807B (en) | A kind of wind power station active power control method based on Model Predictive Control | |
| Pahasa et al. | Coordinated control of wind turbine blade pitch angle and PHEVs using MPCs for load frequency control of microgrid | |
| CN108512258B (en) | Wind power plant active scheduling method based on improved multi-agent consistency algorithm | |
| CN111697578B (en) | Multi-target energy storage-containing regional power grid operation control method | |
| CN107065556A (en) | A kind of automatic search method of reactor core unit Variable power optimization of operation strategy scheme | |
| CN112632774A (en) | Data-driven wind power plant frequency control method based on dynamic modal decomposition | |
| CN114037191A (en) | Virtual power plant optimal scheduling method, device, equipment and medium based on big data | |
| CN115912466A (en) | Active power distribution network island division method and system based on information gap decision theory | |
| CN115173410A (en) | Power grid rotational inertia threshold calculation method for improving GWO algorithm | |
| CN116169698A (en) | Distributed energy storage optimal configuration method and system for stable new energy consumption | |
| CN115313380A (en) | New energy hydrogen production system coordination control method adaptive to hydrogen load fluctuation | |
| Kheshti et al. | Improving frequency regulation of wind‐integrated multi‐area systems using LFA‐fuzzy PID control | |
| CN114336592B (en) | Wind power plant AGC control method based on model predictive control | |
| CN115764931A (en) | Automatic power generation control method, system, equipment and medium for power system | |
| CN117039882A (en) | Resource aggregation regulation and control method and system based on binary consistency algorithm | |
| CN115114854A (en) | Two-stage self-organizing optimization aggregation method and system for distributed resources of virtual power plant | |
| CN114884136A (en) | Active power distribution network robust optimization scheduling method considering wind power correlation | |
| CN115099590A (en) | Active power distribution network economic optimization scheduling method and system considering light load uncertainty | |
| Liu et al. | A comprehensive control strategy for photovoltaic virtual synchronous generator considering frequency regulation capability | |
| CN112636366A (en) | Wind power plant dynamic frequency control method based on control process data fitting | |
| Xu et al. | Design of three-stage start-up optimal strategy based on fuzzy fractional-order proportion integration differentiation controller of pumped storage unit under low water head | |
| CN114865649B (en) | Wind-solar-storage integrated station reactive power regulation method and device and electronic equipment |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |