CN112564132B - Wind power primary frequency modulation potential uncertainty modeling method - Google Patents

Wind power primary frequency modulation potential uncertainty modeling method Download PDF

Info

Publication number
CN112564132B
CN112564132B CN202011469856.5A CN202011469856A CN112564132B CN 112564132 B CN112564132 B CN 112564132B CN 202011469856 A CN202011469856 A CN 202011469856A CN 112564132 B CN112564132 B CN 112564132B
Authority
CN
China
Prior art keywords
frequency modulation
power
wind
primary frequency
average injection
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.)
Active
Application number
CN202011469856.5A
Other languages
Chinese (zh)
Other versions
CN112564132A (en
Inventor
林子钊
阎诚
刘雪飞
田启东
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shenzhen Power Supply Bureau Co Ltd
Original Assignee
Shenzhen Power Supply Bureau Co Ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Shenzhen Power Supply Bureau Co Ltd filed Critical Shenzhen Power Supply Bureau Co Ltd
Priority to CN202011469856.5A priority Critical patent/CN112564132B/en
Publication of CN112564132A publication Critical patent/CN112564132A/en
Application granted granted Critical
Publication of CN112564132B publication Critical patent/CN112564132B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • 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
    • 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]
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A30/00—Adapting or protecting infrastructure or their operation

Landscapes

  • Engineering & Computer Science (AREA)
  • Power Engineering (AREA)
  • Control Of Eletrric Generators (AREA)
  • Wind Motors (AREA)

Abstract

The invention discloses a wind power primary frequency modulation potential uncertainty modeling method which includes the steps of obtaining actual measurement operation data of historical rotor rotating speed, operation wind speed, pitch angle and power under a fan load shedding operation state, identifying wind energy of a fan by using coefficient model parameters, calculating prediction errors of wind power parameter frequency modulation release maximum rotor kinetic energy average injection power and primary frequency modulation average injection power by combining a historical wind speed prediction sequence and a wind energy utilization coefficient model, carrying out kernel density estimation to obtain edge distribution of the historical wind speed prediction sequence and the wind energy utilization coefficient model, building a combined distribution model between the prediction errors by using a dynamic Copula function, generating a probability density function releasing the sum of the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power by a discrete convolution method, and finally calculating uncertainty intervals under different confidence levels. The method can be used for carrying out refined modeling on the uncertain interval of the wind power frequency modulation potential, and has important significance for reasonably optimizing system frequency modulation resources and reducing redundant standby.

Description

Wind power primary frequency modulation potential uncertainty modeling method
Technical Field
The invention relates to the technical field of power system frequency stability control, in particular to a wind power primary frequency modulation potential uncertainty modeling method.
Background
The method aims at the problems that the equivalent rotational inertia of a system is reduced and the risk of system frequency stability is increased due to large-scale wind power access, and challenges are brought to the frequency stability control of the power system. With the diversified development of frequency modulation requirements, the trend of reasonably and efficiently controlling the wind turbine generator to participate in frequency modulation is the future development trend. Therefore, the research on the wind power frequency modulation capability is not limited to the control strategy improvement, and the potential of each wind turbine in the wind power plant actually contributing to the system frequency is considered more, so that a real-time reference is provided for system scheduling.
The traditional control strategy only considers that the energy contained in the rotor is quickly released in a rotor kinetic energy control mode to respond to the frequency change of the system, however, the kinetic energy of the rotor stored by the fan is limited, and the secondary falling of the system frequency can be caused after the frequency modulation is exited, so that the impact is caused on the stability of the power grid. The mainstream development trend of wind power participating in system frequency modulation is that a wind machine can participate in primary frequency modulation through load shedding operation, but the operation state of the wind machine is changed in the load shedding state, and in the frequency response process, both the energy interaction of rotor kinetic energy and the increase of mechanical power exist. How to realize accurate perception of wind power primary frequency modulation potential in the real-time scheduling process of a large-scale wind power access system is a difficult problem to be solved urgently.
The existing research does not consider the influence of uncertainty in the wind power operation process on the frequency modulation potential, and the wind power primary frequency modulation potential indexes have different action characteristics on the system frequency, so that the existing research method cannot simultaneously have different wind power primary frequency modulation potential index difference characteristics to quantify the wind power participation system frequency modulation capacity, and also does not relate to the dynamic correlation research among the wind power primary frequency modulation potential indexes, and therefore, the wind power primary frequency modulation potential cannot be accurately depicted.
Disclosure of Invention
The technical problem to be solved by the invention is to provide a modeling method for uncertainty of primary wind power frequency modulation potential, which can effectively reflect uncertainty of the primary wind power frequency modulation potential.
In order to solve the technical problem, the invention provides a wind power primary frequency modulation potential uncertainty modeling method, which comprises the following steps:
the method comprises the following steps of S1, acquiring actually measured operation data of a historical rotor rotating speed, an operation wind speed, a pitch angle and power in a fan load shedding operation state and preprocessing the actually measured operation data;
s2, identifying wind energy utilization coefficient model parameters of the fan according to the operation data of the fan after load shedding;
s3, calculating a prediction error of the wind power parameter frequency modulation release maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power by combining a historical wind speed prediction sequence and a wind energy utilization coefficient model;
s4, performing kernel density estimation on the prediction error sequence which releases the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power to obtain edge distribution, and establishing a joint distribution model between prediction errors by using a dynamic Copula function;
and S5, generating a probability density function for releasing the sum of the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power by a discrete convolution method, and calculating uncertainty intervals under different confidence levels.
Further, the preprocessing in step S1 specifically includes:
and judging whether data loss and abnormal values exist according to the acquired actually measured operating data, discarding all data at the corresponding moment when the data loss and abnormal values exist in any one of the actually measured data of the rotor rotating speed, the operating wind speed, the pitch angle and the power, and otherwise, reserving the data at the moment.
Further, the step S2 specifically includes:
s21, calculating the tip speed ratio and the wind energy utilization coefficient of the corresponding time point according to the actually measured operation data;
s22, screening the operation data set by judging whether the pitch angle of each moment point is zero or not, and recording the maximum operation wind speed of the fan at the moment point when the pitch angle is zero as a critical wind speed;
s23, calculating a fitting value of sampling points at all times according to the wind energy utilization coefficient model;
and S24, identifying parameters in the wind energy utilization coefficient model.
Further, the tip speed ratio λ is calculated in the following manner:
Figure BDA0002835875900000021
wherein v and omega respectively represent the wind speed and the rotor speed of the fan at the moment, and R is the radius of the wind wheel of the fan.
Actual value C of wind energy utilization coefficient p The calculation method of (A) is as follows:
Figure BDA0002835875900000022
wherein P is the output power of the fan, ρ is the air density, and v and R respectively represent the input wind speed of the fan and the radius of the wind wheel;
the wind energy utilization coefficient model is as follows:
Figure BDA0002835875900000023
where λ and β are the tip speed ratio and pitch angle, respectively, i and j are the order of λ and β, respectively, a ij Is each coefficient;
all the coefficients in the wind energy utilization coefficient model form a matrix A to be solved as shown in the specification C ：
Figure BDA0002835875900000031
Coefficient of wind energy utilization C p The relationship between the pitch angle and the tip speed ratio in the maximum power tracking state is expressed as:
Figure BDA0002835875900000032
when the pitch angle satisfies beta =0, solving the above formula to obtain the optimal tip speed ratio lambda under the maximum power tracking state opt Further simplifying the above formula to obtain:
β＝f b (λ)
calculating an ideal deloading point under a given reserved deloading level according to the running wind speed
Figure BDA0002835875900000033
The method comprises the following steps:
Figure BDA0002835875900000034
/>
wherein, ω is max The maximum rotor speed of the fan;
calculating the actual operating point on the ideal load shedding curve under the given reserved load shedding level by combining the operating wind speed, the actually measured rotor rotating speed and the pitch angle
Figure BDA0002835875900000035
The method comprises the following steps:
Figure BDA0002835875900000036
identifying parameters in the wind energy utilization coefficient model, and determining the matrix A to be solved in the load shedding operation mode C The wind energy utilization coefficient parameter identification problem is expressed as:
Figure BDA0002835875900000037
Figure BDA0002835875900000038
and k represents a sampling point at the kth moment, N is the total length of the sample, and the problem is solved by an interior point method.
Further, in step S3, the method for calculating the predicted value of the average kinetic energy injection power of the maximum rotor released by the wind power parameter modulation includes:
Figure BDA0002835875900000041
the method for calculating the actual value of the average kinetic energy injection power of the maximum rotor released by the wind power parameter modulation comprises the following steps:
Figure BDA0002835875900000042
wherein H W Is the inherent inertia time constant of the fan, T is the time participating in frequency modulation, v and omega are the running wind speed and the actual rotor speed, R is the radius of the wind wheel of the fan, and lambda opt For optimum tip speed ratio, omega max Is the maximum rotor speed, v, of the fan lim Is the critical wind speed;
the method for calculating the prediction error of the average injected power for releasing the maximum rotor kinetic energy comprises the following steps:
ΔP 1 ＝P 1 P -P 1 R ；
the output power of the fan in the maximum power tracking state is represented as:
Figure BDA0002835875900000043
the mode of calculating the predicted value of the primary frequency modulation average injection power is as follows:
Figure BDA0002835875900000044
the way to calculate the actual value of the primary frequency modulation average injection power is:
Figure BDA0002835875900000045
wherein, d% is a reserved load shedding level, and P is an actually sampled fan power value;
the method for calculating the prediction error of the primary frequency modulation average injection power comprises the following steps:
Figure BDA0002835875900000046
further, the step S4 specifically includes:
s41, performing kernel density estimation on the prediction error sequence releasing the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power to obtain the cumulative probability distribution u and v of the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power;
step S42, establishing a joint distribution model between the prediction error cumulative probability distribution u and v by using a dynamic t-Copula function:
Figure BDA0002835875900000051
where ρ is i Is the ith sample dynamic correlation coefficient, upsilon is the degree of freedom,
Figure BDA0002835875900000052
a unitary T-distribution function T as a degree of freedom v υ An inverse function of (·);
step S43, calculating a prediction error sequence joint probability density model for releasing the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power, wherein the method comprises the following steps:
Figure BDA0002835875900000053
wherein, the first and the second end of the pipe are connected with each other,
Figure BDA0002835875900000054
and->
Figure BDA0002835875900000055
Further, obtaining a dynamic t-Copula function maximum likelihood value and an evolution equation through a maximum likelihood estimation method, specifically:
the dynamic t-Copula function time-varying parameter equation is expressed as follows:
Figure BDA0002835875900000056
wherein, theta 1 、θ 2 And theta 3 First, second and third evolution parameters of the time-varying equation,
Figure BDA0002835875900000057
the corrected Logistic function is specifically as follows:
Figure BDA0002835875900000058
the time-varying t-Copula function to-be-estimated parameter expression is integrated into a set
Figure BDA0002835875900000059
In the form of a maximum likelihood estimate, in particularComprises the following steps:
Figure BDA0002835875900000061
wherein, F 1 (. Cndot.) and F 2 The cumulative probability density function of the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power is released.
Further, in step S41, for the known sample, the probability density function and the cumulative probability density function are obtained through kernel density estimation, specifically:
prediction error value X with respect to sample set X = [ X ] 1 ,x 2 ,…,x n ]Is expressed as:
Figure BDA0002835875900000062
where n is the total number of samples, x i For the value of the ith sample, K (-) is the kernel function, and h is the kernel density estimation model bandwidth, expressed as:
Figure BDA0002835875900000063
wherein the content of the first and second substances,
Figure BDA0002835875900000064
is the standard deviation of the sample set;
the cumulative probability density of the variable x is expressed as:
Figure BDA0002835875900000065
wherein G (·) is represented as:
Figure BDA0002835875900000066
further, the step S5 specifically includes:
step S51, generating a probability density function releasing the sum of the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power through a discrete convolution method, wherein for the maximum rotor kinetic energy average injection power release prediction error sample set alpha and the primary frequency modulation average injection power prediction error sample set beta, a convolution expression of alpha + beta based on a Copula function is as follows:
Figure BDA0002835875900000067
wherein, the first and the second end of the pipe are connected with each other,
Figure BDA0002835875900000068
and &>
Figure BDA0002835875900000069
Upper and lower limit values of the sample sets alpha and beta, respectively;
the discretization probability density function of the sum of the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power is as follows:
Figure BDA0002835875900000071
wherein the content of the first and second substances,
Figure BDA0002835875900000072
and &>
Figure BDA0002835875900000073
The upper and lower limit values of the discretized sample set alpha and beta are respectively. P α (. And P) β (. Beta) are respectively the discretized sample alpha and beta probability densities, and the probability density of the discretized sample point for the sample set X is expressed as:
Figure BDA0002835875900000074
where Δ h is the discretization step length, f X (. H) is a continuous probability density function of the sample set X;
and S52, calculating prediction error confidence intervals under different confidence levels according to the sum of the released maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power and the prediction error distribution probability density, calculating a wind power primary frequency modulation potential prediction value through the prediction data of the given wind speed sequence, and overlapping the prediction error intervals to obtain the prediction intervals under different confidence levels.
Further, according to the sum of the released maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power, the prediction error probability distribution probability density is calculated, and prediction error confidence intervals under different confidence levels are specifically as follows:
performing per unit processing on the discrete probability density p (i) with the number of sample points N:
Figure BDA0002835875900000075
calculating the cumulative probability density of each discrete sample point:
Figure BDA0002835875900000076
generating a random variable set U obeying a (0, 1) distribution, and substituting the following formula to obtain a random number set Y obeying a known discrete probability density:
Y＝F -1 (U)
obtaining a cumulative probability distribution function FY (·) of the random number set Y through kernel density estimation, and further calculating an uncertainty interval [ lb, ub ] under the known confidence level alpha:
min ub-lb
Figure RE-GDA0002940562660000081
the embodiment of the invention has the beneficial effects that: the method fully considers the dynamic correlation among the wind power primary frequency modulation potential influence index prediction errors, can realize accurate modeling of the wind power primary frequency modulation potential uncertainty only depending on the operation wind speed sequence, considers the prediction error distribution characteristics of different frequency modulation potential indexes, has simple calculation method and high model solving efficiency, and has important significance for providing reasonable reference for multi-wind motor combined operation coordination control and system scheduling and optimizing frequency modulation spare capacity.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly described below, it is obvious that the drawings in the following description are only some embodiments of the present invention, and for those skilled in the art, other drawings can be obtained according to these drawings without any creative effort.
FIG. 1 is a flow schematic diagram of a wind power primary frequency modulation potential uncertainty modeling method according to an embodiment of the invention.
Fig. 2 is a specific flow diagram of a wind power primary frequency modulation potential uncertainty modeling method according to an embodiment of the invention.
FIG. 3 is a schematic diagram of a dynamic Copula correlation coefficient of a wind power primary frequency modulation potential uncertainty modeling method according to an embodiment of the invention.
FIG. 4 is a schematic diagram of an uncertainty interval of primary wind power frequency modulation potential after the uncertainty modeling method of the primary wind power frequency modulation potential is adopted.
Detailed Description
The following description of the embodiments refers to the accompanying drawings, which are included to illustrate specific embodiments in which the invention may be practiced.
Referring to fig. 1, an embodiment of the present invention provides a wind power primary frequency modulation potential uncertainty modeling method, including:
the method comprises the following steps of S1, acquiring actually measured operation data of a historical rotor rotating speed, an operation wind speed, a pitch angle and power in a fan load shedding operation state and preprocessing the actually measured operation data;
s2, identifying wind energy utilization coefficient model parameters of the fan according to the operation data of the fan after load shedding;
s3, calculating a prediction error of the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power released by wind power parameter frequency modulation through a historical wind speed prediction sequence and a wind energy utilization coefficient model;
s4, performing kernel density estimation on the prediction error sequence releasing the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power to obtain edge distribution, and establishing a joint distribution model between prediction errors by using a dynamic Copula function;
and S5, generating a probability density function for releasing the sum of the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power by a discrete convolution method, and calculating uncertainty intervals under different confidence levels.
Specifically, referring to fig. 2, in step S1, in the embodiment, SCADA collected data of a single fan 2015 in a 1.5MW wind farm in northeast china from 20 days 6 months to 4 days 7 months 4 months is selected as a data sample, and sampling intervals are 1min. The inherent inertia time constant of the fan is 5.04s, the operating wind speed of the fan participating in frequency modulation is 7m/s, the rated wind speed of the fan is 12m/s, and the maximum rotor rotating speed of the fan is 1.2p.u. And judging whether data loss and abnormal values exist according to the acquired data sample at the moment, discarding all data at the corresponding moment when data loss and abnormal values exist in any one of the actually measured data of the rotor rotating speed, the operating wind speed, the pitch angle and the power, and otherwise, reserving the data at the moment.
S2, selecting SCADA acquisition data of 2015, namely, 20 days in 6 months and 3 days in 7 months as a model parameter identification sample, and specifically comprising the following steps:
and S21, calculating the tip speed ratio and the wind energy utilization coefficient of the corresponding time point according to the actually measured operation data. The tip speed ratio may be expressed as:
Figure BDA0002835875900000091
v and omega respectively represent the wind speed and the rotor speed of the fan at the moment, and R is the radius of the fan wheel.
Actual value C of wind energy utilization coefficient p Can be calculated by the following formula:
Figure BDA0002835875900000092
where P is the output power of the fan, ρ is the air density, and v and R represent the fan input wind speed and the rotor radius, respectively.
S22, screening the operation data set by judging whether the pitch angle of each moment point is zero or not, and recording the maximum operation wind speed of the fan at the moment point when the pitch angle is zero as a critical wind speed v lim 。
Step S23, calculating a fitting value of each sampling point at each moment, wherein the wind energy utilization coefficient model can be expressed as:
Figure BDA0002835875900000093
where λ and β are the tip speed ratio and pitch angle, respectively, i and j are the order of λ and β, respectively, a ij For each coefficient, all coefficients in the wind energy utilization coefficient model can form a matrix A to be solved C It can be expressed as:
Figure BDA0002835875900000101
coefficient of wind energy utilization C p The pitch angle and tip speed ratio relationship in the maximum power tracking state may be expressed as:
Figure BDA0002835875900000102
when the pitch angle meets the condition that beta =0, solving the formula can obtain the optimal tip speed ratio lambda under the maximum power tracking state opt Further simplification of the above formula can result in
β＝f b (λ)
Calculating an ideal deloading point at a given reserved deloading level according to the operating wind speed
Figure BDA0002835875900000103
Can be expressed as:
Figure BDA0002835875900000104
wherein, ω is max The maximum rotor speed of the fan.
Calculating the actual operating point on the ideal load shedding curve under the given reserved load shedding level by combining the operating wind speed, the actually measured rotor rotating speed and the pitch angle
Figure BDA0002835875900000105
Can be expressed as:
Figure BDA0002835875900000106
step S24, identifying parameters in the wind energy utilization coefficient model, and determining the coefficient matrix A in the load shedding operation mode C The wind energy utilization coefficient parameter identification problem can be expressed as:
Figure BDA0002835875900000107
Figure BDA0002835875900000108
wherein k represents a sampling point at the kth moment, N is the total length of the sample, and the problem can be solved by an interior point method.
In step S3, a predicted value of the average kinetic energy injection power of the maximum rotor released by the wind power parameter modulation is calculated, which can be expressed as:
Figure BDA0002835875900000111
wherein H W Is the inherent inertia time constant of the fan, v is the running wind speed, R represents the radius of the wind wheel of the fan, omega max Is the maximum rotor speed, v, of the fan lim Representing the critical wind speed, λ opt Is the optimum tip speed ratio.
Calculating the actual value of the average injected power of the maximum rotor kinetic energy released by the wind power parameter frequency modulation, which can be expressed as:
Figure BDA0002835875900000112
wherein H W Is the inherent inertia time constant of the fan, T is the time participating in frequency modulation, v and omega are the running wind speed and the actual rotor speed, R represents the radius of the wind wheel of the fan, and lambda opt For optimal tip speed ratio.
The average injected power prediction error for releasing the maximum rotor kinetic energy can be expressed as:
ΔP 1 ＝P 1 P -P 1 R
the output power of the fan in the maximum power tracking state can be expressed as:
Figure BDA0002835875900000113
calculating the predicted value of the primary frequency modulation average injection power can be expressed as:
Figure BDA0002835875900000114
wherein d% is the reserved load shedding level.
Calculating the actual value of the primary modulation average injection power can be expressed as:
Figure BDA0002835875900000115
and P is the actually sampled fan power value.
The primary modulation average injected power prediction error can be expressed as:
Figure BDA0002835875900000116
in step S4, the predicted value and the actual value of the maximum rotor kinetic energy average injection power released from day 2 of 7 month to day 3 of 7 month in 2015 are selected as samples. Specifically, step S4 includes:
and S41, performing kernel density estimation on the prediction error sequence releasing the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power to obtain the cumulative probability distribution u and v of the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power.
For a known sample, a probability density function and a cumulative probability density function are obtained through kernel density estimation, and the probability density function and the cumulative probability density function are specifically as follows:
prediction error value X with respect to sample set X = [ X ] 1 ,x 2 ,…,x n ]The probability density function of (a) can be expressed as:
Figure BDA0002835875900000121
where n is the total number of samples, x i For the value of the ith sample, K (-) is the kernel function, and h is the kernel density estimation model bandwidth, which can be expressed as:
Figure BDA0002835875900000122
wherein the content of the first and second substances,
Figure BDA0002835875900000123
is the standard deviation of the sample set.
The cumulative probability density of the variable x can be expressed as:
Figure BDA0002835875900000124
wherein G (·) can be represented as:
Figure BDA0002835875900000125
step S42, then, establishing a joint distribution model between the prediction error cumulative probability distribution u and v by using a dynamic t-Copula function, wherein the mode is as follows:
Figure BDA0002835875900000126
where ρ is i Is the ith sample dynamic correlation coefficient, upsilon is the degree of freedom,
Figure BDA0002835875900000127
one-element T distribution function T as degree of freedom upsilon υ Inverse function of (·).
Step S43, calculating a prediction error sequence joint probability density model for releasing the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power, wherein the method comprises the following steps:
Figure BDA0002835875900000131
wherein the content of the first and second substances,
Figure BDA0002835875900000132
and->
Figure BDA0002835875900000133
Obtaining a dynamic t-Copula function maximum likelihood value and an evolution equation by a maximum likelihood estimation method, which comprises the following specific steps:
the dynamic t-Copula function time-varying parametric equation can be expressed as follows:
Figure BDA0002835875900000134
wherein, theta 1 、θ 2 And theta 3 First, second and third evolution parameters of the time-varying equation, respectively.
Figure BDA0002835875900000135
The corrected Logistic function is specifically as follows:
Figure BDA0002835875900000136
the method comprises the steps of combining the expression sets of parameters to be estimated of the time-varying t-Copula function into a set
Figure BDA0002835875900000137
The maximum likelihood estimation is performed, specifically:
Figure BDA0002835875900000138
wherein, F 1 (. And F) 2 (. Cndot.) is the cumulative probability density function of the released maximum rotor kinetic energy average injected power and the primary frequency modulated average injected power, respectively. The correlation coefficients of the time-varying t-Copula function at different time points of the selected sample obtained according to the maximum likelihood estimation result are shown in fig. 3. On the basis, the correlation coefficient of the time-varying t-Copula function can be gradually recurred according to the evolution parameters in the time-varying parameter equation.
Step S5 specifically includes:
step S51, generating a probability density function for releasing the sum of the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power by a discrete convolution method, wherein for the prediction error sample set alpha for releasing the maximum rotor kinetic energy average injection power and the prediction error sample set beta for the primary frequency modulation average injection power, a convolution expression of alpha + beta based on a Copula function is as follows:
Figure BDA0002835875900000139
in the formula (I), the compound is shown in the specification,
Figure BDA0002835875900000141
and &>
Figure BDA0002835875900000142
Upper and lower limits for the sample sets alpha and beta, respectively. The discretization probability density function of the sum of the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power is as follows:
Figure BDA0002835875900000143
/>
wherein the content of the first and second substances,
Figure BDA0002835875900000144
and &>
Figure BDA0002835875900000145
The upper and lower limit values of the discretized sample set alpha and beta are respectively. P α (. And P) β (. Cndot.) are the discretized sample α and β probability densities, respectively, and the probability density for the discretized sample point for sample set X can be expressed as:
Figure BDA0002835875900000146
where Δ h is the discretization step length, f X (. Cndot.) is a continuous probability density function of the sample set X.
And S52, calculating prediction error confidence intervals under different confidence levels according to the sum of the released maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power and the prediction error distribution probability density, calculating a wind power primary frequency modulation potential prediction value through the prediction data of a given wind speed sequence, and overlapping the wind power primary frequency modulation potential prediction value with the prediction error intervals to obtain prediction intervals under different confidence levels.
Calculating prediction error confidence intervals under different confidence levels according to the sum of the released maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power and the prediction error distribution probability density, wherein the method specifically comprises the following steps:
the unitary processing is performed on the discrete probability density p (i) with the number of sample points N, which can be expressed as:
Figure BDA0002835875900000147
the cumulative probability density for each discrete sample point is calculated and can be expressed as:
Figure BDA0002835875900000148
generating a random variable set U obeying a (0, 1) distribution, and substituting the following formula to obtain a random number set Y obeying a known discrete probability density:
Y＝F -1 (U)
the cumulative probability distribution function FY (·) of the random number set Y can be obtained by kernel density estimation, and further the uncertainty interval [ lb, ub ] at the known confidence level α can be calculated:
min ub-lb
Figure BDA0002835875900000151
and taking 144 sampling points which are 7 months and 4 days before 2015 as a verification set. And (3) realizing recursion of the correlation coefficient according to a time-varying parameter equation of a time-varying t-Copula function, namely obtaining a t-Copula function between the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power prediction error released in the next step, obtaining the discrete probability density of the wind power primary frequency modulation potential prediction error through discrete convolution on the basis, solving a corresponding confidence interval and overlapping the confidence interval with the wind power primary frequency modulation potential prediction value, and realizing generation of an uncertainty interval, as shown in fig. 4.
As can be seen from the above description, the embodiments of the present invention have the following beneficial effects: the method fully considers the dynamic correlation among the wind power primary frequency modulation potential influence index prediction errors, can realize accurate modeling of wind power primary frequency modulation potential uncertainty only depending on an operation wind speed sequence, considers the distribution characteristics of prediction errors of different frequency modulation potential indexes, has simple calculation method and high model solving efficiency, and has important significance for providing reasonable reference for multi-wind motor combined operation coordination control and system scheduling and optimizing frequency modulation spare capacity.
While the invention has been described in connection with what is presently considered to be the most practical and preferred embodiment, it is to be understood that the invention is not to be limited to the disclosed embodiment, but on the contrary, is intended to cover various modifications and equivalent arrangements included within the spirit and scope of the appended claims.

Claims (8)

1. A wind power primary frequency modulation potential uncertainty modeling method is characterized by comprising the following steps:
s1, acquiring actually measured operation data of a historical rotor rotating speed, an operation wind speed, a pitch angle and power in a fan load shedding operation state and preprocessing the actually measured operation data;
s2, identifying wind energy utilization coefficient model parameters of the fan according to the operation data of the fan after load shedding;
s3, calculating a prediction error of the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power released by wind power parameter frequency modulation through a historical wind speed prediction sequence and a wind energy utilization coefficient model;
s4, performing kernel density estimation on the prediction error sequence for releasing the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power to obtain edge distribution, and establishing a joint distribution model between prediction errors by using a dynamic Copula function;
s5, generating a probability density function for releasing the sum of the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power by a discrete convolution method, and calculating uncertainty intervals under different confidence levels;
the step S2 specifically includes:
s21, calculating the tip speed ratio and the wind energy utilization coefficient of the corresponding time point according to the actually measured operation data;
s22, screening the running data set by judging whether the pitch angle of each moment point is zero or not, and recording the running maximum wind speed of the fan at the moment point when the pitch angle is zero as a critical wind speed;
s23, calculating a fitting value of each sampling point at each moment according to the wind energy utilization coefficient model;
s24, identifying parameters in the wind energy utilization coefficient model;
the calculation mode of the tip speed ratio lambda is as follows:
Figure FDA0004031789570000011
v and omega respectively represent the wind speed and the rotor speed of the fan at the moment, and R is the radius of a wind wheel of the fan;
actual value C of wind energy utilization coefficient p The calculation method is as follows:
Figure FDA0004031789570000012
wherein P is the output power of the fan, ρ is the air density, and v and R respectively represent the input wind speed of the fan and the radius of the wind wheel;
the wind energy utilization coefficient model is as follows:
Figure FDA0004031789570000013
where λ and β are the tip speed ratio and pitch angle, respectively, i and j are the order of λ and β, respectively, a ij Is each coefficient;
all the coefficients in the wind energy utilization coefficient model form a matrix A to be solved as shown in the specification C ：
Figure FDA0004031789570000021
Coefficient of wind energy utilization C p The relationship between the pitch angle and the tip speed ratio in the maximum power tracking state is expressed as:
Figure FDA0004031789570000022
when the pitch angle satisfies beta =0, solving the above formula to obtain the optimal tip speed ratio lambda under the maximum power tracking state opt And further simplifying the formula to obtain:
β＝f b (λ)
calculating an ideal deloading point at a given reserved deloading level according to the operating wind speed
Figure FDA0004031789570000023
The method comprises the following steps:
Figure FDA0004031789570000024
wherein, ω is max The maximum rotor speed of the fan;
calculating the actual operating point on the ideal load shedding curve under the given reserved load shedding level by combining the operating wind speed, the actually measured rotor rotating speed and the pitch angle
Figure FDA0004031789570000025
The method comprises the following steps:
Figure FDA0004031789570000026
identifying parameters in the wind energy utilization coefficient model, and determining the matrix A to be solved in the load-shedding operation mode C The wind energy utilization coefficient parameter identification problem is expressed as:
Figure FDA0004031789570000027
Figure FDA0004031789570000028
and k represents a sampling point at the kth moment, N is the total length of the sample, and the problem is solved by an interior point method.
2. The wind power primary frequency modulation potential uncertainty modeling method according to claim 1, characterized in that the preprocessing in the step S1 specifically comprises:
and judging whether data loss and abnormal values exist according to the acquired actually measured operation data, discarding all data at the corresponding moment when the phenomena of data loss and abnormal values exist in any one of the actually measured data of the rotor rotating speed, the operation wind speed, the pitch angle and the power, and otherwise, reserving the data at the moment.
3. The wind power primary frequency modulation potential uncertainty modeling method according to claim 1, wherein in the step S3, the mode of calculating the predicted value of the maximum rotor kinetic energy average injection power released by wind power parameter modulation is as follows:
Figure FDA0004031789570000031
the method for calculating the actual value of the average kinetic energy injection power of the maximum rotor released by the wind power parameter modulation comprises the following steps:
Figure FDA0004031789570000032
wherein H W The method is characterized in that the method is a fan inherent inertia time constant, T is the time participating in frequency modulation, v and omega are the running wind speed and the actual rotor speed, R is the fan wind wheel radius, and lambda is opt For optimum tip speed ratio, omega max Is the maximum rotor speed, v, of the fan lim Is the critical wind speed;
the method for calculating the prediction error of the average injected power for releasing the maximum rotor kinetic energy comprises the following steps:
ΔP 1 ＝P 1 P -P 1 R ；
the output power of the fan in the maximum power tracking state is represented as:
Figure FDA0004031789570000033
the mode of calculating the predicted value of the primary frequency modulation average injection power is as follows:
Figure FDA0004031789570000034
the way to calculate the actual value of the primary frequency modulation average injection power is:
Figure FDA0004031789570000041
wherein, d% is a reserved load shedding level, and P is an actually sampled fan power value;
the method for calculating the prediction error of the primary frequency modulation average injection power comprises the following steps:
ΔP 2 ＝P 2 P -P 2 R 。
4. the wind power primary frequency modulation potential uncertainty modeling method according to claim 3, characterized in that the step S4 specifically comprises:
s41, performing kernel density estimation on the prediction error sequence of the released maximum rotor kinetic energy average injection power and the released primary frequency modulation average injection power to obtain the cumulative probability distribution u and v of the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power;
step S42, establishing a joint distribution model between the prediction error cumulative probability distribution u and v by using a dynamic t-Copula function:
Figure FDA0004031789570000045
wherein ρ i Is the ith sample dynamic correlation coefficient, upsilon is the degree of freedom,
Figure FDA0004031789570000046
a unitary T-distribution function T as a degree of freedom v υ An inverse function of (·);
step S43, calculating a prediction error sequence joint probability density model for releasing the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power, wherein the method comprises the following steps:
Figure FDA0004031789570000042
wherein, the first and the second end of the pipe are connected with each other,
Figure FDA0004031789570000043
and->
Figure FDA0004031789570000044
5. The wind power primary frequency modulation potential uncertainty modeling method according to claim 4, characterized in that a dynamic t-Copula function maximum likelihood value and an evolution equation are obtained by a maximum likelihood estimation method, specifically:
the dynamic t-Copula function time-varying parameter equation is expressed as:
Figure FDA0004031789570000051
wherein, theta 1 、θ 2 And theta 3 First, second and third evolution parameters of the time-varying equation,
Figure FDA0004031789570000052
is a corrected Logistic functionThe body is as follows:
Figure FDA0004031789570000053
the method comprises the steps of combining the expression sets of parameters to be estimated of the time-varying t-Copula function into a set
Figure FDA0004031789570000054
The maximum likelihood estimation is performed, specifically:
Figure FDA0004031789570000055
wherein, F 1 (. And F) 2 The cumulative probability density function of the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power is released.
6. The wind power primary frequency modulation potential uncertainty modeling method according to claim 4, characterized in that in step S41, a probability density function and an accumulated probability density function are obtained through kernel density estimation for a known sample, specifically:
prediction error value X with respect to sample set X = [ X = 1 ,x 2 ,…,x n ]Is expressed as:
Figure FDA0004031789570000056
where n is the total number of samples, x i For the value of the ith sample, K (-) is the kernel function, and h is the kernel density estimation model bandwidth, expressed as:
Figure FDA0004031789570000057
wherein the content of the first and second substances,
Figure FDA0004031789570000058
is the standard deviation of the sample set;
the cumulative probability density of the variable x is expressed as:
Figure FDA0004031789570000059
wherein G (·) is represented as:
Figure FDA00040317895700000510
7. the wind power primary frequency modulation potential uncertainty modeling method according to claim 4, wherein the step S5 specifically comprises:
step S51, generating a probability density function for releasing the sum of the maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power by a discrete convolution method, wherein for the prediction error sample set alpha for releasing the maximum rotor kinetic energy average injection power and the prediction error sample set beta for the primary frequency modulation average injection power, a convolution expression of alpha + beta based on a Copula function is as follows:
Figure FDA0004031789570000061
wherein, the first and the second end of the pipe are connected with each other,
Figure FDA0004031789570000062
and &>
Figure FDA0004031789570000063
Upper and lower limit values of sample sets α and β, respectively;
the discretization probability density function of the sum of the released maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power is as follows:
Figure FDA0004031789570000064
wherein the content of the first and second substances,
Figure FDA0004031789570000065
and &>
Figure FDA0004031789570000066
Respectively representing upper and lower limit values of alpha and beta of the discretized sample set; p α (. And P) β (. Cndot.) are the discretized sample α and β probability densities, respectively, and the probability density for the discretized sample point for sample set X is represented as:
Figure FDA0004031789570000067
where Δ h is the discretization step length, f X () is a continuous probability density function of the sample set X;
and S52, calculating prediction error confidence intervals under different confidence levels according to the sum of the released maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power and the prediction error distribution probability density, calculating a wind power primary frequency modulation potential prediction value through the prediction data of the given wind speed sequence, and overlapping the prediction error intervals to obtain the prediction intervals under different confidence levels.
8. The wind power primary frequency modulation potential uncertainty modeling method according to claim 7, characterized by calculating prediction error confidence intervals under different confidence levels according to the sum of the released maximum rotor kinetic energy average injection power and the primary frequency modulation average injection power, specifically:
performing per-unit processing on the discrete probability density p (i) with the number of sample points being N:
Figure FDA0004031789570000068
calculating the cumulative probability density of each discrete sample point:
Figure FDA0004031789570000071
generating a random variable set U obeying a (0, 1) distribution, and substituting the following formula to obtain a random number set Y obeying a known discrete probability density:
Y＝F -1 (U)
obtaining a cumulative probability distribution function FY (·) of the random number set Y through kernel density estimation, and further calculating an uncertainty interval [ lb, ub ] under the known confidence level alpha:
min ub-lb
Figure FDA0004031789570000072
/>
CN202011469856.5A 2020-12-15 2020-12-15 Wind power primary frequency modulation potential uncertainty modeling method Active CN112564132B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011469856.5A CN112564132B (en) 2020-12-15 2020-12-15 Wind power primary frequency modulation potential uncertainty modeling method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011469856.5A CN112564132B (en) 2020-12-15 2020-12-15 Wind power primary frequency modulation potential uncertainty modeling method

Publications (2)

Publication Number Publication Date
CN112564132A CN112564132A (en) 2021-03-26
CN112564132B true CN112564132B (en) 2023-04-14

Family

ID=75064658

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011469856.5A Active CN112564132B (en) 2020-12-15 2020-12-15 Wind power primary frequency modulation potential uncertainty modeling method

Country Status (1)

Country Link
CN (1) CN112564132B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113193572B (en) * 2021-04-30 2022-11-11 国网江苏省电力有限公司 Method for realizing fan control by modeling wind power frequency modulation parameter probability distribution
CN113708389B (en) * 2021-09-10 2023-08-04 国网湖南省电力有限公司 Wind farm primary frequency modulation model parameter identification method and system based on actual power response
CN116599163B (en) * 2023-04-27 2024-01-23 华能烟台风力发电有限公司 High-reliability wind farm power control system based on frequency modulation control

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106528912A (en) * 2016-09-19 2017-03-22 国网浙江省电力公司经济技术研究院 Method for estimating frequency regulation capacity of wind power plant

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102063575B (en) * 2011-01-01 2013-01-02 国网电力科学研究院 Method for analyzing influence of output power fluctuation of wind farm on power grid
CN103066620B (en) * 2012-12-24 2014-10-22 中国电力科学研究院 Design method of automatic generation control model under intermittent energy grid-connection
CN103902837B (en) * 2014-04-16 2017-02-15 广西大学 Method for wind speed prediction based on experience Copula function
CN103927695B (en) * 2014-04-22 2017-11-24 国家电网公司 Ultrashort-term wind power prediction method based on self study complex data source
CN107425520B (en) * 2017-06-12 2020-04-21 东南大学 Active power distribution network three-phase interval state estimation method containing node injection power uncertainty
CN111092444B (en) * 2019-12-25 2021-11-09 浙江运达风电股份有限公司 Method and device for controlling inertial energy support of variable-speed constant-frequency wind turbine generator
CN111900743B (en) * 2020-07-28 2021-11-16 南京东博智慧能源研究院有限公司 Wind power frequency modulation potential prediction error distribution estimation method

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106528912A (en) * 2016-09-19 2017-03-22 国网浙江省电力公司经济技术研究院 Method for estimating frequency regulation capacity of wind power plant

Also Published As

Publication number Publication date
CN112564132A (en) 2021-03-26

Similar Documents

Publication Publication Date Title
CN112564132B (en) Wind power primary frequency modulation potential uncertainty modeling method
CN111900743B (en) Wind power frequency modulation potential prediction error distribution estimation method
CN109038686B (en) Rolling optimization scheduling method based on wind power output prediction error
CN108317040B (en) Method, device, medium, equipment and wind generating set for correcting yaw to wind
CN107045574B (en) SVR-based effective wind speed estimation method for low wind speed section of wind generating set
CN114865701A (en) Wind storage combined frequency modulation method based on adaptive model predictive control
CN110991701A (en) Wind power plant fan wind speed prediction method and system based on data fusion
CN115510677A (en) Wind power plant generating capacity evaluation method and system
CN115660412A (en) Scheduling risk quantitative evaluation method and device considering source load fluctuation and rescheduling
CN115189401A (en) Day-ahead-day coordinated optimization scheduling method considering source load uncertainty
CN115030866A (en) Wind power plant group control system
CN110649638B (en) Optimization method of energy storage system for compensating wind power prediction error
CN116562422A (en) Combined clearing method for electric energy and climbing reserve and terminal
CN111160653B (en) Distributed energy storage system wind power consumption capability monitoring method based on cloud computing
CN113468767B (en) Method and system for evaluating generating capacity of offshore wind turbine
CN116742622B (en) Photovoltaic power generation-based power generation amount prediction method and system
CN116599163B (en) High-reliability wind farm power control system based on frequency modulation control
CN112653132B (en) Method, system, device and medium for judging stability of offshore wind power-containing power system
CN113193572B (en) Method for realizing fan control by modeling wind power frequency modulation parameter probability distribution
CN117313415A (en) Wind power plant fan blade extension transformation site selection evaluation method, device and storage medium
CN110502807B (en) Robust optimization-based transmission equipment limit passing power and utilization rate evaluation method
Wei et al. Calibration power curve of wind generator and forecasting model of wind power unit
CN115660249A (en) Comprehensive evaluation method for wind power prediction error
CN117543617A (en) Combined clearing method and system for frequency modulation auxiliary service market and energy market
CN115907406A (en) Medium-and-long-term power and electric quantity balancing method and device, electronic equipment and storage medium

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
GR01 Patent grant
GR01 Patent grant