Method for predicting wind speed in the rotor plane for a wind turbine equipped with a LiDAR sensor

Abstract

The present invention is a method for predicting the wind speed in the rotor plane (PR) of a wind turbine (1), by accounting for an induction factor used in a wind evolution model implemented by a Kalman filter. The invention also is a method for controlling a wind turbine (1), a computer program product, a LiDAR sensor (2) and a wind turbine (1), which uses the wind prediction determined with the method according to the invention.

Claims

1. A method for predicting wind speed in a rotor plane of a wind turbine equipped with a LiDAR sensor, comprising: a) measuring the wind speed in at least one measurement plane distant from the wind turbine by use of the LiDAR sensor; b) determining a wind induction factor, representing a wind deceleration coefficient between the measurement plane and the plane of the rotor; c) determining a delay index k.sub.d0 between the measurement plane and the rotor plane of the wind turbine by use of the induction factor and the wind speed measurement in the measurement plane by use of an equation $k_{d0}=\frac{2x_1}{\left(U_{x1}+U_0\right)T_s};$ with U.sub.0=a.sub.0,x.sub.1U.sub.x.sub.1, with x.sub.1 being the distance between the measurement plane and the rotor plane, T.sub.s being the measurement sampling period, U.sub.x1 being an average wind speed measured in the measurement plane, U.sub.0 being an average wind speed in the rotor plane and a.sub.0,x1 being the induction factor between the measurement plane and the rotor plane; d) constructing a wind evolution model between the measurement plane and the plane of the rotor, the wind evolution model connecting the wind speed in the plane of the rotor at a predetermined time to the measured wind speed in the measurement plane at times prior to the predetermined time, by use of the induction factor with the prior times being determined by use of the delay index; and e) determining the wind speed prediction in the rotor plane by use of the wind evolution model and a Kalman filter.

2. A prediction method as claimed in claim 1, wherein a wind induction factor is determined by carrying out the following steps: i) measuring wind speed in at least three measurement planes distant from the wind turbine by use of the LiDAR sensor; ii) determining at least two wind induction factors between two of the measurement planes using the wind speed measurements in the measurement planes and a linear Kalman filter; and iii) determining the wind induction factor between a measurement plane and the rotor plane of the wind turbine by use of the determined induction factors between two measurement planes and using a linear Kalman filter.

3. A prediction method as claimed in claim 2, wherein the wind speed measurement step comprises reconstructing a wind field in the measurement plane which is used in other steps of the method as for wind speed measurement in the measurement plane.

4. A prediction method as claimed in claim 2, wherein the wind evolution model is expressed as follows:
u.sub.0(k+p)=Ũ.sub.x.sub.1(k−k.sub.d0+p).sup.Tr(k|k), with
Ũ.sub.x.sub.1(k−k.sub.d0+p)=[ũ.sub.x.sub.1(k−k.sub.d0+p)ũ.sub.x.sub.1(k−k.sub.d0+p−1)ũ.sub.x.sub.1(k−k.sub.d0+p+1) . . . ũ.sub.x.sub.1(k−k.sub.d0+p+N.sub.d)].sup.T and ũ.sub.x.sub.1(k)=a.sub.0,x.sub.1u.sub.x.sub.1(k), with u.sub.0 being the wind in rotor plane, k being discretized time, p being a future time step, k.sub.d0 being the delay index, r being a state vector determined by the Kalman filter, x.sub.1 being a measurement plane, N.sub.d being an order of the wind evolution model, u.sub.x1 being wind speed measured in the measurement plane and a.sub.0,x1 being the induction factor between the measurement plane and the rotor plane.

5. A prediction method as claimed in claim 1, wherein the wind speed measurement step comprises reconstructing a wind field in the measurement plane which is used in other steps of the method as for wind speed measurement in the measurement plane.

6. A prediction method as claimed in claim 5, wherein the wind evolution model is expressed as follows:
u.sub.0(k+p)=Ũ.sub.x.sub.1(k−k.sub.d0+p).sup.Tr(k|k), with
Ũ.sub.x.sub.1(k−k.sub.d0+p)=[ũ.sub.x.sub.1(k−k.sub.d0+p)ũ.sub.x.sub.1(k−k.sub.d0+p−1)ũ.sub.x.sub.1(k−k.sub.d0+p+1) . . . ũ.sub.x.sub.1(k−k.sub.d0+p+N.sub.d)].sup.T and ũ.sub.x.sub.1(k)=a.sub.0,x.sub.1u.sub.x.sub.1(k), with u.sub.0 being the wind in rotor plane, k being discretized time, p being a future time step, k.sub.d0 being the delay index, r being a state vector determined by the Kalman filter, x.sub.1 being a measurement plane, Nd being an order of the wind evolution model, u.sub.x1 being wind speed measured in the measurement plane and a.sub.0,x1 being the induction factor between the measurement plane and the rotor plane.

7. A prediction method as claimed in claim 1, wherein the wind evolution model is expressed as follows:
u.sub.0(k+p)=Ũ.sub.x.sub.1(k−k.sub.d0+p).sup.Tr(k|k), with
Ũ.sub.x.sub.1(k−k.sub.d0+p)=[ũ.sub.x.sub.1(k−k.sub.d0+p)ũ.sub.x.sub.1(k−k.sub.d0+p−1)ũ.sub.x.sub.1(k−k.sub.d0+p+1) . . . ũ.sub.x.sub.1(k−k.sub.d0+p+N.sub.d)].sup.T and ũ.sub.x.sub.1(k)=a.sub.0,x.sub.1u.sub.x.sub.1(k), with u.sub.0 being the wind in rotor plane, k being discretized time, p being a future time step, k.sub.d0 being the delay index, r being a state vector determined by the Kalman filter, x.sub.1 being a measurement plane, N.sub.d being an order of the wind evolution model, u.sub.x1 being wind speed measured in the measurement plane and a.sub.0,x1 being the induction factor between the measurement plane and the rotor plane.

8. A prediction method as claimed in claim 7, wherein an order of the wind evolution model Nd is 2.

9. A method of controlling a wind turbine equipped with a LiDAR sensor, comprising steps of: a) determining the wind speed prediction in the rotor plane of the wind turbine by use of the method as claimed in claim 1; and b) controlling the wind turbine as a function of the wind speed in the rotor plane.

10. A non-transiently recorded computer program product, comprising code instructions for carrying out the method of claim 1, when the program is executed on a processor which performs processing for the LiDAR sensor.

11. A LiDAR sensor for a wind turbine, comprising a processing unit which implements a non transiently recorded method which performs a method as claimed in claim 1.

12. A wind turbine comprising a LiDAR sensor as claimed in claim 11, wherein the LiDAR sensor disposed on a nacelle of the wind turbine or in a nose of the wind turbine.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) Other features and advantages of the method according to the invention will be clear from reading the description hereafter of embodiments given by way of non-limitative example, with reference to the accompanying drawings wherein:

(2) FIG. 1 illustrates a wind turbine equipped with a LiDAR sensor according to an embodiment of the invention;

(3) FIG. 2 illustrates the steps of the method for predicting the wind speed in the rotor plane according to an embodiment of the invention;

(4) FIG. 3 illustrates the steps of the method for predicting the wind speed in the rotor plane according to a second embodiment of the invention;

(5) FIG. 4 illustrates a comparative curve in relation to a reference wind for an example of a prediction at 0s of the wind in the rotor plane for a method according to the prior art and for the prediction method according to an embodiment of the invention;

(6) FIG. 5 illustrates a comparative curve in relation to a reference wind for an example of the prediction at 1s of the wind in the rotor plane for a method according to an embodiment of the invention; and

(7) FIG. 6 illustrates a comparative curve in relation to a reference wind for an example of the prediction at 2s of the wind in the rotor plane for a method according to an embodiment of the invention.

DETAILED DESCRIPTION OF THE INVENTION

(8) The present invention relates to a method for predicting the wind speed in the rotor plane of a wind turbine equipped with a LiDAR sensor. The method enables prediction of the wind speed because it allows determination for at least one future time the wind speed in the rotor plane of the turbine.

(9) According to the invention, the LiDAR sensor allows measurement of the wind speed in at least one measurement plane upstream from the wind turbine. There are several types of LiDAR sensors, for example scanning LiDAR, continuous wave LiDAR or pulsed LiDAR sensors. Within the context of the invention, a pulsed LiDAR is preferably used. However, the other LiDAR technologies may also be used while remaining within the scope of the invention.

(10) LiDAR sensors allow fast measurement. Therefore, using such a sensor enables fast and continuous determination of the wind speed prediction. For example, the sampling rate of the LiDAR sensor can range between 1 and 5 Hz (or even more in the future), and it can be 4 Hz. Furthermore, the LiDAR sensor allows obtaining data relative to the wind upstream from the wind turbine. This data is related to the wind that is going to hit the turbine. The LiDAR sensor can therefore be used for predicting the wind speed in the rotor plane of the wind turbine.

(11) FIG. 1 schematically shows, by way of non-limitative example, a horizontal-axis wind turbine 1 equipped with a LiDAR sensor 2 for the method according to an embodiment of the invention. LiDAR sensor 2 is used to measure the wind speed at a given distance in multiple measurement planes PM (only two measurement planes are shown). Knowing the wind measurement in advance a priori allows providing much data. This figure also shows axes x, y and z. The reference point of this coordinate system is the center of the rotor. Direction x is the longitudinal direction corresponding to the direction of the rotor axis, upstream from the wind turbine. This direction also corresponds to the measurement direction of LiDAR sensor 2. Direction y, perpendicular to direction x, is the lateral direction located in a horizontal plane (directions x, y form a horizontal plane). Direction z is the vertical direction (substantially corresponding to the direction of tower 4) pointing up, axis z is perpendicular to axes x and y. The rotor plane is indicated by the rectangle in dotted line PR. It is defined by directions y, z for a zero value of x. Measurement planes PM are planes formed by directions y, z at a distance from rotor plane PR (for a non-zero value of x). Measurement planes PM are parallel to rotor plane PR.

(12) Conventionally, a wind turbine 1 allows converting the kinetic energy of the wind into electrical or mechanical energy. To convert the wind energy to electrical energy, it is made up of the following elements:

(13) a tower 4 allowing a rotor (not shown) to be positioned at a sufficient height to enable motion thereof (necessary for horizontal-axis wind turbines) or this rotor to be positioned at a height enabling it to be driven by a stronger and more regular wind than at ground level 6. Tower 4 generally houses part of the electrical and electronic components (modulator, control, multiplier, generator, etc.),

(14) a nacelle 3 mounted at the top of tower 4, housing mechanical, pneumatic and some electric and electronic components (not shown) necessary for operating the machine. Nacelle 3 can rotate to orient the machine in the right direction,

(15) the rotor, fastened to the nacelle, comprises several blades 7 (generally three) and the nose of the wind turbine. The rotor is driven by the energy from the wind and it is connected by a mechanical shaft, directly or indirectly (via a gearbox and mechanical shaft system), to an electrical machine (electric generator) (not shown) that converts the energy recovered to electrical energy. The rotor is potentially provided with control systems such as variable-angle blades or aerodynamic brakes,

(16) a transmission made up of two shafts (mechanical shaft of the rotor and mechanical shaft of the electrical machine) connected by a transmission (gearbox) (not shown).

(17) As can be seen in FIG. 1, which is an example embodiment of a pulsed LiDAR sensor, the LiDAR sensor 2 used comprises four beams or measurement axes (b1, b2, b3, b4). By way of non-limitative example, the method according to the invention also operates with a LiDAR sensor comprising any number of beams. The LiDAR sensor performs a punctual measurement at each intersection point of a measurement plane PM and a beam (b1, b2, b3, b4). These measurement points are represented by black circles in FIG. 1. Processing the measurements at these measurement points allows determination of the wind speed in measurement planes PM. The wind modelling method described in French patent application FR-3,068,139 (WO-2018/234,409) can therefore notably be applied.

(18) Preferably, LiDAR sensor 2 can be mounted on nacelle 3 of wind turbine 1.

(19) According to the invention, the wind prediction method comprises the following steps:

(20) 1) Wind speed measurement

(21) 2) Induction factor determination

(22) 3) Delay index determination

(23) 4) Construction of a wind evolution model

(24) 5) Determination of the wind speed prediction.

(25) These steps are carried out in real time.

(26) FIG. 2 schematically illustrates, by way of non-limitative example, the steps of the prediction method according to an embodiment of the invention. The first step is a step (MES) of measuring the wind speed u.sub.x in at least one measurement plane by use of the LiDAR sensor. The second step (IND) uses measurements u.sub.x to determine an induction factor a.sub.PR between the measurement plane and the rotor plane. The third step (RET) uses induction factor a.sub.PR and measurements u.sub.x to determine a delay index k.sub.d. The fourth step is a step of constructing a wind evolution model (MOD) from measurements u.sub.x, induction factor a.sub.PR and delay index k.sub.d. A Kalman filter (KAL), notably an adaptive Kalman filter, is then applied to the wind evolution model (MOD) to obtain a prediction of the wind speed in the rotor plane u.sub.0.

(27) FIG. 3 schematically illustrates, by way of non-limitative example, the steps of the prediction method according to a second embodiment of the invention. The first step is a step (MES) of measuring the wind speed in at least one measurement plane by means of the LiDAR sensor. In relation to the embodiment of FIG. 1, a step (REC) of reconstructing the wind from the measurements prior to determining wind speed u.sub.x in the measurement plane is added. The third step (IND) uses the reconstructed measurements u.sub.x to determine an induction factor a.sub.PR between the measurement plane and the rotor plane. The fourth step (RET) uses induction factor a.sub.PR and reconstructed measurements u.sub.x to determine a delay index k.sub.d. The fifth step is a step of constructing a wind evolution model (MOD) from reconstructed measurements u.sub.x, induction factor a.sub.PR and delay index k.sub.d. A Kalman filter (KAL), notably an adaptive Kalman filter, is then applied to the wind evolution model to obtain a prediction of the wind speed in the rotor plane u.sub.0.

(28) 1. Wind Speed Measurement

(29) In this step, the wind speed is continuously measured in at least one measurement plane distant from the wind turbine, by using the LiDAR sensor. Thus, the wind speed can be known upstream from the wind turbine in at least one measurement plane. For this step, the three wind components (longitudinal, lateral and vertical) are considered, as well as the wind speed variations in the measurement plane (the wind speed increases with height for example).

(30) According to an embodiment of the invention, the wind speed can be measured in at least three measurement planes to facilitate induction factor determination (see step 2). Preferably, the wind speed can be measured in four measurement planes to improve the accuracy of knowledge of the wind upstream from the wind turbine and, therefore, the accuracy of estimation of the wind induction factor in the rotor plane.

(31) According to an implementation of the invention, the measurement planes can be located within a longitudinal distance (along axis x in FIG. 1) preferably ranging between 50 and 400 m from the rotor plane. It is thus possible to determine the evolution of the wind speed over a long distance upstream from the wind turbine, which increases the accuracy of estimation of the wind speed prediction in the rotor plane.

(32) Alternatively, the measurement planes may be located at a shorter or a longer distance in relation to the aforementioned preferred range.

(33) According to a non-limitative example embodiment, the LiDAR sensor can perform measurements for ten measurement planes, which can notably be located at distances of 50, 70, 90, 100, 110, 120, 140, 160, 180 and 200 m from the rotor plane respectively.

(34) According to an implementation of the invention corresponding to the embodiment of FIG. 3, the wind measurement step can comprise a wind reconstruction step allowing the raw measurements of the LiDAR sensor to be interpreted. Various wind reconstruction methods can be implemented. By way of non-limitative example, the wind reconstruction method described in patent application FR-3,068,139 (WO-2018/234,409), whose main steps are noted hereafter, may be implemented: gridding the space located upstream from the LiDAR sensor, the grid pattern comprising estimation points and measurement points; measuring the wind amplitude and direction at the different measurement points; estimating the wind amplitude and direction at any time for all the estimation points, using a recursive least-squares method of a cost function; and reconstructing the incident wind field in three dimensions and in real time on all the discretized points.

(35) 2. Induction Factor Determination

(36) This step determines a wind induction factor between a measurement plane and the rotor plane. The induction factor is a wind deceleration coefficient in the induction zone of the wind turbine (that is upstream from the turbine). Indeed, the wind is slowed down upstream from the turbine by the presence of the turbine and its operation: in other words, the power drawn from the wind by the turbine causes wind deceleration upstream from the turbine. Thus, the induction factor is representative of a physical phenomenon and it provides an indication relative to the resource of the wind turbine. The induction factor is calculated between two planes upstream from the turbine, by definition, it corresponds to the speed ratio between these planes. By denoting the induction factor, u the wind speed, and d1 and d2, the respective distances of the two planes considered in relation to the rotor plane, the induction factor between the planes located at distances d1 and d2 from the rotor plane can be written:

(37) $a_{d1,d2}=\frac{u_{d1}}{u_{d2}}$

(38) When one of the planes is considered in the rotor plane, d1=0.

(39) For this step, the induction factor can be determined by any method. By way of non-limitative example, the method of determining an induction factor is disclosed in the French patent application FR-1,871,455, whose steps are described below, can be implemented.

(40) The method of determining a wind induction factor between a measurement plane and the rotor plane of a wind turbine comprises the following steps:

(41) A—Wind Speed Measurement

(42) This step is identical to step 1) of the prediction method according to the invention for at least three measurement planes. Preferably, a single series of measurements can be performed for steps 1) and 2)A.

(43) B—Determining Induction Factors Between Measurement Planes

(44) This step determines in real time at least two wind induction factors between two measurement planes. In other words, at a minimum, a first induction factor is determined between a first measurement plane and a second measurement plane, and a second induction factor is determined between a third measurement plane and a fourth measurement plane (at maximum one of the third and fourth measurement planes corresponds to one of the first and second measurement planes). In a simplified manner, these induction factors are referred to as induction factors between measurement planes in the description hereafter.

(45) According to the invention, the induction factors between measurement planes are determined by use of the wind speed measurements and of a Kalman filter, notably a linear Kalman filter.

(46) Advantageously, for this step and for the embodiment where the measurement is performed for at least four measurement planes, at least three induction factors of the wind between measurement planes can be determined. It is thus possible to increase the accuracy of the wind deceleration phenomenon upstream from the wind turbine and, therefore, the accuracy of estimation of the wind induction factor in the rotor plane.

(47) According to an embodiment of the invention, an induction factor between measurement planes can be determined for the measurement plane closest to the rotor plane, to obtain information relative to the wind field as close as possible to the rotor plane.

(48) According to an implementation of the invention, the induction factors between measurement planes can be determined for measurement planes having the same spacing. For example, the spacing can be 50 m and a first induction factor can be determined for measurement planes located at 70 and 120 m, and a second induction factor for measurement planes located at 90 and 140 m.

(49) Preferably, the spacing between the measurement planes used for the wind induction factors between measurement planes can be identical to the distance between the rotor plane and the closest measurement plane. Thus, the induction factor model is simplified, which notably facilitates solution of the Kalman filter.

(50) For example, by combining the variants described above, if the first measurement plane is located at 50 m from the rotor plane, a first induction factor can be determined for measurement planes located at 50 and 100 m, a second induction factor for measurement planes located at 70 and 120 m, and a third induction factor for measurement planes located at 90 and 140 m.

(51) According to an embodiment of the invention, the Kalman filter can be used with the different steps described hereafter. The steps are described for an embodiment where a first induction factor a.sub.50,100 is determined for measurement planes located at 50 and 100 m, a second induction factor a.sub.70,120 for measurement planes located at 70 and 120 m, and a third induction factor a.sub.90,140 for measurement planes located at 90 and 140 m.

(52) Below, it is only shown how to estimate a.sub.50,100 in real time with a.sub.70,120 and a.sub.90,140 being obtained exactly in the same way. Since u.sub.50, u.sub.100 are available in real time, the induction factor definition equation could be directly used to determine a.sub.50,100. However, this method involves drawbacks. Indeed, the information on the standard deviation of the estimated wind speeds u.sub.50, u.sub.100 is not used. The standard deviation of estimation a.sub.50,100 can therefore not be known. On the other hand, a calculation stability problem may occur for low speeds, i.e. when u.sub.100 is close to zero.

(53) The induction factor at the time k is denoted by a.sub.50,100(k). It is clear that the variation a.sub.50,100(k)-a.sub.50,100(k−1) is relatively small, therefore it can be expressed as:
a.sub.50,100(k)=a.sub.50,100(k−1)+η(k−1)
where η(k−1) is used to describe the variation of a.sub.50,100(k) over time.

(54) The first induction factor definition equation is rewritten as:
u.sub.100(k)a.sub.50,100(k)=u.sub.50(k).

(55) Since estimations u.sub.50(k), u.sub.100(k) contain noise, a more realistic model of the above equation is:
(u.sub.100(k)+ϵ.sub.100(k))a.sub.50,100(k)=u.sub.50(k)+ϵ.sub.50(k)
where ε.sub.50(k), ε.sub.100(k) are the noises for u.sub.50(k), u.sub.100(k) respectively. The previous equation can then be rewritten as follows:
u.sub.100(k)a.sub.50,100(k)=u.sub.50(k)+ϵ.sub.50(k)−ϵ.sub.100(k)a.sub.50,100(k)

(56) By combining the previous equations, the following equation of state is obtained:

(57) $\{\begin{array}{c}{a}_{50,100}\left(k\right)=a_{50,100}\left(k-1\right)+\eta \left(k-1\right),\\ u_{50}\left(k\right)=u_{100}\left(k\right)a_{50,100}\left(k\right)+\mu \left(k\right)\end{array}\mu \left(k\right)={ϵ}_{100}\left(k\right)a_{50,100}\left(k\right)-{ϵ}_{50}\left(k\right)$

(58) One way of estimating the unknown state vector a.sub.50,100(k) that can account for information on ε(k) and μ(k) applies a Kalman filtering algorithm referred to as Kalman filter. In practice, this filter provides the solution to the following problem:

(59) $\underset{a_{50,100}\left(k\right)}{\min }J\left(k\right)$$\mathrm{with}$$J\left(k\right)={\left(a_{50,100}\left(0\right)-{\stackrel{\_}{a}}_{50,100}\left(0\right)\right)}^T{P}_0^{-1}\left(a_{50,100}\left(0\right)-{\stackrel{\_}{a}}_{50,100}\left(0\right)\right)+\underset{j=1}{\overset{k}{.Math.}}\left({\eta \left(j-1\right)}^T{Q}^{-1}\eta \left(k-1\right)+{\mu \left(j\right)}^T{R}^{-1}\mu \left(j\right)\right)$
where P.sub.0, Q, R are weighting matrices of suitable dimension, a.sub.50,100(0) with an overbar is the mean value of initial state a.sub.50,100(0).

(60) In order to solve the optimization problem using the Kalman filtering algorithm, the following assumptions are made. These assumptions mainly relate to a mathematical interpretation for P.sub.0, Q, R. a.sub.50,100(0) is a random vector not correlated with noises ε(k) and μ(k), a.sub.50,100(0) has a known mean with P0 as the covariance matrix, that is:
P.sub.0=E[(a.sub.50,100(0)−ā.sub.50,100(0))(a.sub.50,100(0)−ā.sub.50,100(0)).sup.T]
with ā.sub.50100(0) being the mean value of the initial state, ε(k) and μ(k) are white noises with zero mean and not correlated with covariance matrices Q and R respectively, i.e.:

(61) $E\left[\eta \left(k\right){\eta \left(j\right)}^T\right]=\{\begin{array}{cc}Q,& \mathrm{if}k=j,\\ 0,& \mathrm{if}k\ne j\end{array}E\left[\mu \left(k\right){\mu \left(j\right)}^T\right]=\{\begin{array}{cc}R,& \mathrm{if}k=j,\\ 0,& \mathrm{if}k\ne j\end{array}E\left[\eta \left(k\right){\mu \left(j\right)}^T\right]=0,\mathrm{for}\mathrm{all}k,j$ It is noted that this assumption also implies that Q and R are positive semi-definite symmetric matrices.

(62) The following notations are adopted: â.sub.50,100(k|k−1) is the estimation of a.sub.50,100(k) given the time measurements k−1. a.sub.50,100(k|k) is the estimation of a.sub.50,100(k) given the time measurements k. P(k|k−1) is the covariance matrix of a.sub.50,100(k) given the time measurements k−1. P(k|k) is the covariance matrix of a.sub.50,100(k) given the time measurements k.

(63) The Kalman filtering algorithm can then be summarized as follows: Time update equation:

(64) $\{\begin{array}{c}{\hat{a}}_{50,100}\left(k.Math.k-1\right)={\hat{a}}_{50,100}\left(k-1.Math.k-1\right)\\ P\left(k.Math.k-1\right)=P\left(k-1.Math.k-1\right)+Q\end{array}\hspace{1em}$ Measurement update equation:

(65) $\{\begin{array}{c}K\left(k\right)=P\left(k.Math.k-1\right){u_{100}\left(k\right)}^T{\left(P\left(k.Math.k-1\right)+u_{100}\left(k\right){{\mathrm{Ru}}_{100}\left(k\right)}^T\right)}^{-1}\\ {\hat{a}}_{50,100}\left(k.Math.k\right)={\hat{a}}_{50,100}\left(k.Math.k-1\right)+K\left(k\right)\left(u_{50}\left(k\right)-u_{100}\left(k\right)\hat{x}\left(k.Math.k-1\right)\right),\\ P\left(k.Math.k\right)=\left(I-K\left(k\right)\right)P\left(k.Math.k-1\right)\end{array}\hspace{1em}$

(66) By carrying out these steps, induction factor a.sub.50,100 can be determined. These steps are repeated to determine induction factors a.sub.70,120 and a.sub.90,140.

(67) C—Determining the induction factor between a measurement plane and the rotor plane

(68) This step determines in real time the wind induction factor between one of the measurement planes and the rotor plane. Thus, the evolution of the wind at the rotor can be represented by accounting for the physical phenomena, in particular the wind deceleration. According to the invention, the wind induction factor between a measurement plane and the rotor plane is determined by use of the induction factors determined in the previous step and using a Kalman filter, notably a linear Kalman filter. To simplify, this induction factor is referred to as induction factor in the rotor plane in the description hereafter.

(69) Preferably, the wind induction factor can be determined between the measurement plane closest to the rotor and the rotor plane.

(70) For example, in this step, the induction factor can be determined between a measurement plane located 50 m from the rotor and the rotor plane.

(71) According to an embodiment of the invention, the Kalman filter can be used by applying the various steps described below. The steps are described for an embodiment for which a.sub.50,100, a.sub.70,120 and a.sub.90,140 have been determined and for which a.sub.0,50 is determined, that is the induction factor between a measurement plane located at 50 m and the rotor plane.

(72) Using a.sub.50,100, a.sub.20,120 and a.sub.90,140, the main idea of the estimation of a.sub.0,50 is to assume that a.sub.0,50, a.sub.50,100, a.sub.70,120 and a.sub.90,140 are a function of distance. The following relation is therefore assumed:

(73) $\{\begin{array}{c}{a}_{0,50}=0x_1+x_2,\\ a_{50,100}=50x_1+x_2,\\ a_{70,120}=70x_1+x_2,\\ a_{90,140}=90x_1+x_2\end{array}\hspace{1em}$
where x.sub.1, x.sub.2 are unknown parameters that need to be determined. Since a.sub.50,100, a.sub.70,120 and a.sub.90,140 change slowly over time, the same applies to x.sub.1, x.sub.2. We can express the relationships:

(74) $\hspace{1em}\{\begin{array}{c}{x}_1\left(k\right)=x_1\left(k-1\right)+{Ϛ}_1\left(k-1\right),\\ x_2\left(k\right)=x_2\left(k-1\right)+{Ϛ}_2\left(k-1\right)\end{array}$
where ζ1(k), ζ2(k) are used to characterize the variation of x.sub.1(k), x.sub.2(k). This relationship is expressed as:

(75) 0$x\left(k\right)=\left[\begin{array}{c}{x}_1\left(k\right)\\ x_2\left(k\right)\end{array}\right],Ϛ\left(k\right)=\left[\begin{array}{c}{Ϛ}_1\left(k\right)\\ {Ϛ}_2\left(k\right)\end{array}\right]$

(76) Then the relationship can be expressed in compact form:
x(k)=x(k−1)+(k−1).

(77) A more realistic model accounting for the noises is:

(78) $\hspace{1em}\{\begin{array}{c}{a}_{50,100}\left(k\right)=\left[\begin{array}{cc}50& 1\end{array}\right]x\left(k\right)+{ϵ}_1\left(k\right),\\ a_{70,120}\left(k\right)=\left[\begin{array}{cc}70& 1\end{array}\right]x\left(k\right)+{ϵ}_2\left(k\right),\\ a_{90,140}\left(k\right)=\left[\begin{array}{cc}90& 1\end{array}\right]x\left(k\right)+{ϵ}_3\left(k\right),\end{array}$
wherein ε1, ε2, ε3 are the noises of estimations a.sub.50,100(k), a.sub.70,120(k) and a.sub.90,140(k) which can then be expressed as:

(79) $y\left(k\right)=\left[\begin{array}{c}{a}_{50,100}\left(k\right)\\ a_{70,120}\left(k\right)\\ a_{90,140}\left(k\right)\end{array}\right],C=\left[\begin{array}{c}\begin{array}{cc}50& 1\end{array}\\ \begin{array}{cc}70& 1\end{array}\\ \begin{array}{cc}90& 1\end{array}\end{array}\right],ϵ\left(k\right)=\left[\begin{array}{c}{ϵ}_1\left(k\right)\\ {ϵ}_2\left(k\right)\\ {ϵ}_3\left(k\right)\end{array}\right]$

(80) Then the relationship can be expressed in compact form as:
y(k)=Cx(k)+ε(k).

(81) By combining the previous equations, the equation of state is as follows:

(82) $\hspace{1em}\{\begin{array}{c}x\left(k\right)=x\left(k-1\right)+Ϛ\left(k-1\right)\\ y\left(k\right)=\mathrm{Cx}\left(k\right)+ϵ\left(k\right)\end{array}$

(83) As for a.sub.50,100(k), a.sub.70,120(k) and a.sub.90,140(k), one way of obtaining x(k) that accounts for noises ζ(k), ε(k) uses the linear Kalman filter technique. The same steps as those described for the previous step are therefore applied.

(84) Once x(k) is estimated, induction factor a.sub.0,50 can be calculated as follows:
a.sub.0,50(k)=[01]x(k).

(85) 3. Delay Index Determination

(86) This step determines a delay index by use of the wind speed measurements (optionally the reconstructed measurements) and by use of the induction factor. The delay index defines the time delay of the wind between the measurement plane and the rotor plane of the wind turbine. The dependence of the delay index on the induction factor allows accounting for physical phenomena (induction) when determining the delay, which makes the wind prediction more representative of real conditions.

(87) According to an embodiment of the invention, the delay index k.sub.d0 can be expressed by an equation:

(88) $k_{d0}=\frac{2x_1}{\left(U_{x1}+U_0\right)T_s}$$U_0=a_{0,x_1}{U}_{x_1}$
with x.sub.1 being the distance between the measurement plane and the rotor plane, T.sub.s being the measurement sampling period, U.sub.x1 being the average wind speed measured in the measurement plane (determined from the measurements), U.sub.0 being the average wind speed in the rotor plane and a.sub.0,x1 being the induction factor between the measurement plane and the rotor plane determined in the previous step.

(89) 4. Construction of a Wind Evolution Model

(90) This step constructs a wind evolution model. A wind evolution model provides the evolution of the wind speed between the measurement plane and the rotor plane of the wind turbine. The wind evolution model connects the wind speed in the rotor plane at a predetermined time t1 (t1 being a future time, the time of the wind speed prediction in the rotor plane) to the wind speed measured in the measurement plane at several times prior to time t1, by use of the induction factor. The times prior to time t1 are determined by use of the delay index. The dependence of the wind evolution model on the induction factor and the delay index allows accounting for physical phenomena (induction) when determining the wind evolution, which makes the wind prediction more representative of real conditions. Furthermore, considering prior times allows accounting for a larger amount of data relative to the wind which improves the prediction.

(91) Preferably, the times prior to time t1 can be times around a time t2 that corresponds to the time t1 less the delay related to the displacement speed of the wind. In other words, time t2 can be determined by use of delay index k.sub.d0 with a formula:
t2=t1−k.sub.d0×T.sub.s.
with T.sub.s being the sampling period of the measurement.

(92) The prior times accounted for in the model can then be times t, such that:
t=t2±n×T.sub.s
with n being a positive or zero integer less than or equal to Nd, which is the order of the wind evolution model.

(93) According to an embodiment of the invention, the wind evolution model can be written as follows:
u.sub.0(k+p)=Ũ.sub.x.sub.1(k−k.sub.d0+p).sup.Tr(k|k)
Ũ.sub.x.sub.1(k−k.sub.d0+p)=[ũ.sub.x.sub.1(k−k.sub.d0+p)ũ.sub.x.sub.1(k−k.sub.d0+p−1)ũ.sub.x.sub.1(k−k.sub.d0+p+1) . . . ũ.sub.x.sub.1(k−k.sub.d0+p+N.sub.d)].sup.T
ũ.sub.x.sub.1(k)=a.sub.0,x.sub.1u.sub.x.sub.1(k)
with u.sub.0 being the wind in the rotor plane, k being the discretized time, p being the future time step, k.sub.d0 being the delay index, r being the state vector determined by the Kalman filter (in step 5), x.sub.1 being the measurement plane, N.sub.d being the order of the wind evolution model, u.sub.x1 being the wind speed measured in the measurement plane (optionally reconstructed) and a.sub.0,x1 being the induction factor between the measurement plane and the rotor plane determined in step 2.

(94) Preferably, in order to limit the computation time, the order of the wind evolution model N.sub.d can be 2.

(95) 5. Determination of the Wind Speed Prediction

(96) In order to apply the wind evolution model constructed in the previous step to determine the wind speed prediction in the rotor plane, a Kalman filter is applied to determine the state vector r of the wind evolution model. Preferably, an adaptive Kalman filter is applied.

(97) For this step, the notations n(k), p(k), J(k), P0, Q and R are specific to this step (the same notations are used independently in step 2).

(98) According to an embodiment of the invention, this step identifies vector r defined by values r.sub.0, r.sub.−1, r.sub.1, . . . , r.sub.−Nd, r.sub.Nd. To characterize the variable nature of the wind evolution, the following model is used: j=−N.sub.d, −N.sub.d+1, . . . , 0, . . . , N.sub.d−1, N.sub.d.

(99) where n.sub.j(k−1) is used to describe the variation of r.sub.j(k).

(100) This equation can be rewritten in compact vector form:

(101) $\begin{array}{cc}& r\left(k\right)=r\left(k-1\right)+\eta \left(k-1\right)\\ \mathrm{où}& \\ & \{\begin{array}{ccccccc}r\left(k\right)=& [r_0\left(k\right)& r_{-1}\left(k\right)& r_1\left(k\right)& \mathrm{.Math.}& r_{-N_d}\left(k\right)& {r_{N_d}\left(k\right)]}^T\\ \eta \left(k\right)=& [{\eta }_0\left(k\right)& {\eta }_{-1}\left(k\right)& {\eta }_1\left(k\right)& \mathrm{.Math.}& \eta -N_d\left(k\right)& {\eta N_d\left(k\right)]}^T\end{array};\end{array}$

(102) This can then bewritten:
Ũ.sub.x.sub.2(k−k.sub.d).sup.Tr(k)=u.sub.x.sub.1(k)
with
U.sub.x.sub.2(k−k.sub.d)=[ũ.sub.x.sub.2(k−k.sub.d)ũ.sub.x.sub.2(k−k.sub.d−1)ũ.sub.x.sub.2(k−k.sub.d+1) . . . ũ.sub.x.sub.2(k−k.sub.d+N.sub.d)].sup.T

(103) It is noted that the values of u in this equation are estimated values containing noises. A more realistic model of this equation can therefore be:
(Ũ.sub.x.sub.2(k−k.sub.d).sup.T+ϵ.sub.x.sub.2(k−k.sub.d).sup.T)r(k)=u.sub.x.sub.1(k)+ϵ.sub.x.sub.1(k)
where ϵ.sub.x1(k) and ϵ.sub.x2(k−k.sub.d) are noises for u.sub.x1(k) and Ũ.sub.x2(k−k.sub.d) respectively. The previous equation can then be rewritten as follows:
Ũ.sub.x.sub.2(k−k.sub.d).sup.Tr(k)=u.sub.x.sub.1(k)+ϵ.sub.x.sub.1(k)−ϵ.sub.x.sub.2(k−k.sub.d).sup.Tr(k)

(104) The following state space equation can then be obtained:

(105) $\begin{array}{cc}\{\begin{array}{cc}r\left(k\right)& =r\left(k-1\right)+\eta \left(k-1\right)\\ {{\tilde{U}}_{x_2}\left(k-k_d\right)}^{T}r\left(k\right)& =u_{x_1}\left(k\right)+\mu \left(k\right)\end{array}& \\ \mathrm{with}& \\ \mu \left(k\right)={ϵ}_{x_1}\left(k\right)-{{ϵ}_{x_2}\left(k-k_d\right)}^{T}r\left(k\right)& \end{array}$
μ(k) being the measurement noise. It needs to be noted that this noise depends on the estimation of state r(k). In other words, the covariance matrix of the measurements is a function of r(k).

(106) One way of estimating the unknown state vector r(k) that accounts for information on η(k) and μ(k) applies a Kalman filter. In practice, the Kalman filter provides the solution to the following problem:

(107) $\begin{array}{cc}\underset{r\left(k\right)}{\min }J\left(k\right)& \\ \mathrm{with}& \\ J\left(k\right)={\left(r\left(0\right)-\tilde{r}\left(0\right)\right)}^T{P}_0^{-1}\left(r\left(0\right)-\tilde{r}\left(0\right)\right)+\underset{j=1}{\overset{k}{.Math.}}\left({\eta \left(j-1\right)}^T{Q}^{-1}\eta \left(j-1\right)+{\mu \left(j\right)}^T{R}^{-1}\mu \left(j\right)\right)& \end{array}$
where P.sub.0, Q, R are weighting matrices of suitable dimension and r(0) is the mean value of r(0).

(108) In order to solve the optimization problem using the Kalman filter, the following assumptions are made. These assumptions mainly relate to a mathematical interpretation for P.sub.0, Q, R. r(0) is a random vector not correlated with noises η(k) and μ(k), r(0) has a known mean r(0) with P0 being the covariance matrix, that is:
P0=E[(r(0)−r(0))(r(0)−r(0)).sup.T] η(k) and μ(k) are white noise processes with zero mean and are not correlated, i.e.:

(109) $E\left[\left(\eta \left(k\right){\eta \left(j\right)}^T\right)\right]=0,\forall k\ne j,E\left[\left(\mu \left(k\right){\mu \left(j\right)}^T\right)\right]=0,\forall k\ne j,E\left[\left(\eta \left(k\right){\mu \left(j\right)}^T\right)\right]=0,\forall k,j$ η(k) has a covariance matrix Q, that is:
E[η(k)η(k).sup.T)]=Q

(110) For the measurement noise, it is obtained:
E[η(k)η(k).sup.T)]=E[(ϵ.sub.x.sub.1(k)−ϵ.sub.x.sub.2(k−k.sub.d).sup.Tr(k))(ϵ.sub.x.sub.1(k)−ϵ.sub.x.sub.2(k−k.sub.d).sup.Tr(k)).sup.T]
E[η(k)η(k).sup.T)]=E(ϵ.sub.x.sub.1(k)ϵ.sub.x.sub.1(k).sup.T)+r(k).sup.TE(ϵ.sub.x.sub.2(k−k.sub.d)ϵ.sub.x.sub.2(k−k.sub.d).sup.T)r(k)
R=E(ϵ.sub.x.sub.1(k)ϵ.sub.x.sub.1(k).sup.T)+r(k).sup.TE(ϵ.sub.x.sub.2(k−k.sub.d)ϵ.sub.x.sub.2(k−k.sub.d).sup.T)r(k)

(111) It can be noted that the covariance matrix of the measurement noise depends on the estimated state r(k).

(112) The following notations are adopted: r(k|k−1) is the estimation of r(k) from the time k−1 r(k|k) is the estimation of r(k) from the measurements of time k P(k|k−1) is the error covariance matrix of r(k) from the measurements of time k−1 P(k|k) is the error covariance matrix of r(k) from the measurements of time k.

(113) The Kalman filter can then be summarized as follows: Prediction

(114) $\hspace{1em}\{\begin{array}{cc}r\left(k|k-1\right)& =r\left(k-1|k-1\right)\\ P\left(k|k-1\right)& =P\left(k-1|k-1\right)+Q\end{array}$ Correction

(115) 0$\hspace{1em}\{\begin{array}{cc}S\left(k\right)& ={\tilde{U}}_{x_2}\left(k-k_d\right)P\left(k|k-1\right){{\tilde{U}}_{x_2}\left(k-k_d\right)}^T+R\\ K\left(k\right)& =P\left(k|k-1\right){{\tilde{U}}_{x_2}\left(k-k_d\right)}^T{S\left(k\right)}^{-1}\\ r\left(k|k\right)& =r\left(k|k-1\right)+K\left(k\right)\left(u_{x_1}\left(k\right)-{\tilde{U}}_{x_2}\left(k-k_d\right)r\left(k|k-1\right)\right)\\ P\left(k|k\right)& =\left(I-K\left(k\right)\right)P\left(k|k-1\right)\end{array}.$

(116) Thus, this Kalman filter allows determination of state vector r, which is then integrated in the wind evolution model determined in the previous step, in order to determine the wind prediction in the rotor plane.

(117) The present invention also relates to a method of controlling a wind turbine equipped with a LiDAR sensor. The following steps are carried out for this method: predicting the wind speed in the rotor plane by use of the wind speed prediction method according to any one of the above variants; and controlling the wind turbine according to the wind speed in the rotor plane.

(118) Accurate real-time prediction of the wind speed in the rotor plane allows suitable wind turbine control in terms of minimization of the effects on the wind turbine structure and maximization of the recovered power. Indeed, by use of this control, the LiDAR allows reduction of the loads on the structure, with the blades and the tower representing 54% of the cost. Therefore, using a LiDAR sensor allows optimizing the wind turbine structure and thus decreasing the costs and maintenance.

(119) According to an implementation of the invention, the inclination angle of at least one of the blades and the electrical recovery torque of the wind turbine generator can be controlled as a function of the wind speed. Other types of regulation devices can also be used.

(120) According to an embodiment of the invention, the inclination angle of at least one of the blades and the electrical recovery torque are determined by use of wind turbine maps as a function of the wind speed at the rotor. For example, the control method described in French patent application FR-2,976,630 A1 corresponding to US published patent application 2012/0,321,463 can be applied.

(121) Furthermore, the invention relates to a computer program product comprising code instructions designed to carry out the steps of one of the methods described above (method of determining the wind speed in the rotor plane, control method). The program is executed on a processor which provides processing for the LiDAR sensor, or on any similar medium connected to the LiDAR sensor or to the wind turbine.

(122) According to an aspect, the present invention also relates to a LiDAR sensor for a wind turbine, comprising a processing unit configured to implement one of the methods described above (method of determining the wind speed in the rotor plane, control method).

(123) According to an implementation of the invention, the LiDAR sensor can be a scanning LiDAR, a continuous wave LiDAR or a pulsed LiDAR sensor. Preferably, the LiDAR sensor is a pulsed LiDAR sensor.

(124) The invention also relates to a wind turbine, notably an offshore or an onshore wind turbine equipped with a LiDAR sensor as described above. According to an embodiment of the invention, the LiDAR sensor can be arranged on the nacelle of the wind turbine or in the nose of the wind turbine. The LiDAR sensor is oriented to perform a measurement of the wind upstream from the wind turbine (that is before the wind turbine and along the longitudinal axis thereof, designated by axis x in FIG. 1). According to an embodiment, the wind turbine can be similar to the turbine illustrated in FIG. 1.

(125) For the embodiment of the control method, the wind turbine can comprise a control, for example for control of the pitch angle of at least one blade of the wind turbine or of the electrical torque, for implementing the control method according to the invention.

Examples

(126) The features and advantages of the prediction method according to the invention will be clear from reading the comparative examples presented hereafter.

(127) For these examples, the wind speed in the rotor plane is predicted with a method according to the prior art wherein Taylor's frozen turbulence hypothesis is applied and the wind speed in the rotor plane is predicted with a method according to the invention, for a wind turbine equipped with a LiDAR sensor.

(128) In order to validate the method according to the invention, it is assumed that the wind turbine is in the measurement plane at 50 m. Thus, by use of the measurements of the LiDAR sensor, a reference is considered which is the averaged measured wind speed in the measurement plane at 50 m, and the method according to the prior art and the method according to the invention is used to estimate the wind speed in the measurement plane at 50 m from a measurement plane distant from the measurement plane at 50 m, for example a measurement plane at 100 m.

(129) FIG. 4 is a curve illustrating the averaged wind speed in the rotor plane RAWS (average wind speed in the rotor plane in the area formed by the turbine blades) in m/s as a function of time T in s. FIG. 4 corresponds to a prediction at the time t=0s (real time). This figure shows the averaged reference wind speed REF, the averaged wind speed predicted with a method according to the prior art AA wherein Taylor's frozen turbulence hypothesis is applied and the averaged wind speed is predicted with a method according to the invention INV. It is noted that the wind speed estimated with the method according to the invention INV is more representative of reference speed REF than the wind speed estimated with the method according to the prior art AA.

(130) FIG. 5 is a curve illustrating the averaged wind speed in the rotor plane RAWS in m/s as a function of time T in s. FIG. 5 corresponds to a prediction at the time t=1 s (prediction at is of the real time). This figure shows the averaged reference wind speed REF and the averaged wind speed predicted with a method according to the invention INV. It is noted that the wind speed estimated with the method according to the invention INV is close to reference speed REF.

(131) FIG. 6 is a curve illustrating the averaged wind speed in the rotor plane RAWS in m/s as a function of time T in s. FIG. 6 corresponds to a prediction at the time t=2s (prediction at 2s of the real time). This figure shows the averaged reference wind speed REF and the averaged wind speed predicted with a method according to the invention INV. It is noted that the wind speed estimated with the method according to the invention INV is close to reference speed REF.

(132) The method according to the invention provides a wind speed prediction in the rotor plane having good accuracy.