Fatigue minimising power reference control of a de-rated wind farm

Modern wind farms (cluster of wind turbines) can be required to control the total power output to meet a set-point, and would then profit by minimising the structural loads and thereby the cost of energy. In this paper, we propose a new control strategy for a derated wind farm with the objective of maintaining a desired reference power production for the wind farm, while minimising the sum of fatigues on the wind turbines in steady-state. The controller outputs a vector of power references for the individual turbines. It exploits the positive correlation between fatigue and added turbulence to minimise fatigue indirectly by minimising the added turbulence. Simulated results for a wind farm with three turbines demonstrate the efficacy of the proposed solution by assessing the damage equivalent loads.


Introduction
As more wind generation is installed, there is growing interest in actively controlling wind turbines or wind farms (clusters of wind turbines) to meet power set-points provided by the operator while minimizing the cost of energy. This can be done by de-rating the turbines to extract the requested energy from the wind while reducing structural loads. Most research on wind farm (WF) control has focused on maximising total WF active power [1]. Following a reference for total WF active power is also very relevant for ancillary services [2]. In the latter case where the WF is de-rated there is a possibility to distribute the single wind turbine (WT) power references in a way that 'minimises' fatigue loading for the WTs in the farm. In a number of papers it is assumed that reducing average loads, e.g. tower thrust, also reduces fatigue loads [3], [4]. As explained in [1] this is not a safe assumption. Turbulence is a better proxy for fatigue as a reduction of turbulence at a WT will reduce the fatigue if the average and frequency contents for the wind is unchanged. In [5] this idea is pursued for WF control and it is shown that the added turbulence can be reduced. However, in [5] the relation to fatigue is not included in the model, so the fatigue effect of changing power references is hidden. Furthermore, the model used in [5] assumes instant propagation of wakes. In this work, we also explore fatigue minimisation indirectly by minimising the added turbulence generated by wakes. The contributions of the paper are as follows. We investigate two control strategies using numerical optimisation for minimising the added turbulence. The first is based on standard optimisation libraries for static constrained non-linear optimisation. However, such a static optimisation may be cumbersome if the number of variables become very large and is generally not fit for distributed calculations. Therefore, we also investigate a novel dynamic controller based on penalty functions and gradient search. The dynamic approach to solving the optimisation problem which we suggest here is an extremum seeking controller similar to the control suggested in [6]. The controller has the benefit that calculations can be distributed across the WTs in the WF, thus eliminating the single-pointof-failure a higher level supervisory controller constitutes. Extremum seeking control has been investigated in wind power production applications in [7] (single WT) and [8], [9] (multiple WTs), among others. In these works the objective is obtaining maximum power extraction. However, as mentioned, in this work we consider a power reference control for a de-rated wind farm with the objective of maintaining the desired power reference while minimising the added turbulence on the WTs in the wind farm. An additional contribution of this paper is the evaluation of the suggested controllers by using SimWindFarm 1 which is a engineering type WF model including a dynamic wind flow model with wake propagation and WT models with a 2 inertia drive train model and a one degree of freedom (DOF) tower model. This allows evaluation of the effect of power reference changes from the controller on the drive train and tower fatigue. To the authors knowledge this evaluation of the effects of added turbulence minimisation on fatigue in wind farm control is new.
The outline of the paper is as follows. In Section 2, we start by describing the model of added turbulence which has been used throughout this work. This section also introduces the optimisation problem which we wish to solve along with the proposed solutions. Section 3 gives the results of numerical simulations performed in the SimWindFarm simulation environment for MATLAB Simulink. In Section 4 we discuss the obtained results and in Section 5 we conclude and give future perspectives.
Nomenclature For a function f : R n × R m → R we denote by ∇ x f (x, y) the partial derivative of f (·, ·) with respect to x, that is, ∇ x f (x, y) = (∂f (x, y)/∂x 1 , . . . , ∂f (x, y)/∂x n ) T . By 1 we denote a column vector consisting of ones and of the appropriate dimension.

Methodology
The fatigue of the WF is here defined as the sum of fatigue over all the WTs in the WF. A central hypothesis here is that the positive correlation between fatigue of a WT and turbulence experienced by the WT can be used for fatigue reduction. The developed WF controller exploits this correlation to minimise the fatigue of the WF indirectly by minimising the added turbulence.

Model of Added Turbulence Intensity
The model of added turbulence intensity which will be used in this work is the one suggested in the IEC 61400-1 [10] which comes from [11] and is given by where σ add,j is the added standard deviation in the wind field at the jth WT which is in wake; U j is the effective wind speed at the jth WT; s ij is the spacing in rotor diameters between the wake generating WT and the WT in wake; C T,i is the thrust coefficient of the wake generating WT. With a model in place which captures the effect of interest, we will now describe the optimisation problem to be solved by the farm controller.

Optimisation problem
In this work we seek to solve the following constrained optimisation problem min F f arm (P ) , s.t. : where F f arm (·) is the total fatigue of the WF; P = (P 1 , . . . , P n ) with P i being the power production of the ith WT; P dem is the demanded reference power production of the WF; P i,min is the minimum power production of the ith WT before shutdown; P i,max is the maximal power production of the ith WT and 0 < P i,min < P i,max . We first simplify the problem (2) by assuming that the total fatigue of the entire WF is given as the sum of fatigues on the individual WTs.
That is where F i is the fatigue of the ith WT. The next approximation is that fatigue is proportional to turbulence. That is where I ef f,i (·) is the effective turbulence intensity experienced by the ith WT; I amb is the ambient turbulence intensity; I add,i (·) ≥ 0 is the added turbulence intensity experienced by the ith WT. Note that frequent control updates will likely add to the fatigue and penalty on control updates could be included in the fatigue expression. However, here we investigate which effect minimisation of added turbulence has on the fatigue as our main hypothesis is that minimising added turbulence will minimise total fatigue loading. Since I amb is unaffected by P , we approximate the optimisation problem (2) to the problem of minimising the added turbulence induced by wakes. The new optimisation problem thus becomes min J(P ) = min i I add,i (P ) , s.t. : We immediately see that the feasibility set of the problem (5) is convex. So (5) is a nonlinear optimisation problem over a convex constraint set. We are now ready to introduce the investigated control strategy.

WF Reference Control 2.3.1. Static Optimising Control
The first design we consider here is a control which solves the original constrained optimisation problem (5) statically. We use sequential quadratic programming (see e.g. [12]) as implemented in the optimisation library nlopt [13] to solve the optimisation problem (5), and update the power references infrequently (every 24 hours) as the simulated ambient wind field is static (constant mean and variance). To solve the optimisation problem, nlopt needs to be able to evaluate the objective function and its gradient as well as the equality constraint and its gradient, whereas the inequality constraints can be defined as upper/lower bounds of the optimisation problem since they are fixed. The Hessian matrix of the problem is estimated by the algorithm. To be able to evaluate I add,j (P i ) using (1) we need first to find an expression for C T,i (P i ). To this end, we first note that the power P i produced by the ith WT is given by where ρ is the air density; A i is the area swept by the rotor; U i is the local mean wind speed at the WT; C P,i is the power coefficient; i is the power in the wind. Furthermore, from actuator disc theory [14] we have the identities where a i ∈ (0; 1) is the axial induction factor. The maximal efficiency of the WT is attained when C P,i is maximal. This maximum (16/27) is attained when a i = 1/3. Since we in this work consider a de-rated WF, we will here consider WTs operating below this maximum value. As a consequence, we will assume that a i ∈ (0; 1/3). For a i ∈ (0; 1/3) the functions in (7) are both one-to-one which means in this interval there is a one-to-one relation between C P,i and C T,i . Next, we set P i,max in (5) to the value where P i,rated is the rated power of the ith WT and K i is evaluated using wind measurements. The assignment (8) ensures that the feasibility set of (5) only contains points where C P,i ∈ (0, 16/27). In the optimisation algorithm, we can now evaluate C T,i (P i ) by first calculating C P,i as C P,i = P i /K i . Then the one-to-one correspondence between C P,i and C T,i is used to find C T,i . In practice this is done using look-ups in a table generated using (7). Using this C T,i the expression (1) and consequently the objective function in (5) can be evaluated. To be able to evaluate the gradient of the objective function, we first express the power P i as a function of C T,i and calculate the gradient of I add,j (P i ) with respect to P i using the chain rule. That is where by the inverse function theorem [15], we have given that (7) are both one-to-one, we can hence isolate a i in the right equation, which gives where C T,i ∈ (0; 8/9). Substituting the above into the top expression in (7) and using the new expression for C P,i in (6) gives Now, we are ready to calculate the partial derivatives of I add,j and P i with respect to C T,i . These are given by where β ij = 0.8s ij . Using (13) in (9) we obtain It is worth noting that for C T,i ∈ (0; 8/9) it follows that ∂ ∂P i I add,j > 0 which implies that increasing the power production on the wake generating WT increases the added turbulence intensity on WTs in its wake. Evaluation of the equality constraint and its gradient in (5) is straightforward.  (5) as a static optimisation problem at fixed time intervals has the benefit that standard optimisation libraries can be used. However, such an optimisation can be cumbersome if the number of variables become very large and it is generally not fit for distributed calculation which means that a single-point-of-failure is introduced. Here we will introduce a dynamic optimisation algorithm to solve (5) which can be distributed between the WTs. For the purpose of handling the inequality constraints in (5) we here introduce penalty functions into the objective function (see [16] and [17]). Thereby, we can rewrite the constrained problem (5) to an unconstrained one and utilise gradient search to minimise J(·). ByJ(P ) we will denote the objective function obtained by the introduction of penalty functions. To handle the constraints in (5) we introduce the functions s i : R → R and r : R → R (see [16] and [17]) given by where κ > 0 is a gain, which for simplicity is the same for every s i (·) and r(·). Using s i (·) to handle the ith inequality constraint and r(·) to handle the equality constraint the problem (5) is rewritten to where P j,min ≤ x j and P j,max ≥ x j . To solve the unconstrained optimisation problem, the following gradient search is usedṖ where L = diag(L i ) with L i > 0 is a gain matrix. Again, the evaluation of the gradients ∇ P ref s j (P ref ) and ∇ P ref r(P ref ) is straightforward and the evaluation of the gradient ∇ P ref I add,i (P ref ) is performed using (14). In particular, if we partition the vector P as P = (P T ,P T ) T whereP is the vector of power outputs from the WTs which has at least one WT in their wake andP is the vector of power outputs from WTs which has no WTs in their wake then we havė where e = k P k −P dem is the reference tracking error for the WF. The last term in (18) and (19) comes from the fact that ∂ ∂P ref e = 1 if we assume that the internal power reference controller on each WT is able to fulfil P k = P ref,k almost all time. Note that the only global information needed at each WT is the tracking error e. Existence of minima ofJ(·) can be proven by [18,Theorem 9.1.3]. Furthermore, from [17,Theorem 1.25] it follows that if P * is a minimiser ofJ(·) then ∇ P * J (P * ) = 0. Thus all local minimisers ofJ(·) are stable equilibria of (17). However, there are no analytic guarantees that constraints are met during transients.

Results
Advanced farm simulators like SOWFA [19] includes WT models with many DOFs and CFD type flow simulation. However, simulations with these models is very slow even when using a few WTs, short simulation time and super computers [20]. Consequently, in this work we use the SimWindFarm toolbox for Matlab Simulink from the Aeolus project [21]. SimWindFarm is a suitable compromise between accuracy and computing speed for this investigation. To be able to evaluate the efficacy of the developed control, SimWindFarm is extended to include the capability of simulating added turbulence using the model (1). Multiple wakes at a turbine are merged using a weighted average based on the area of the wake. The WF used in the simulations consists of three NREL 5MW WTs (see [22]) in a row with the wind direction along the row. Three farm layouts with three different distances between the WTs have been simulated, corresponding to approximately 3, 5 and 10 rotor diameters. Four scenarios in total are simulated and the results are then compared. In the first scenario, the power reference for the WF is evenly distributed between the individual WTs throughout the simulation and added turbulence is not simulated. In the second scenario, the power reference for the WF is also evenly distributed and added turbulence is simulated. Comparison of these two scenarios allows assessment of the implementation of the model (1) ) are determined from the shaft torsion and tower bending moment time series using the rainflow algorithm in Matlab provided by Adam Niesłony [23], where n i is the number of times the load with amplitude S i occurs. The number of cycles for equivalent load N eq = 631·10 6 corresponding approximately 20 years for a 1 Hz signal and a Wöhler exponent of m = 3 are used in all DEL calculations.

Simulation of Added Turbulence
The comparison between the simulations with and without added turbulence are summarised in Figure 1. The top row of the figure shows the result of damage equivalent load (DEL) calculations on the shaft torsion moment (STM) and the bottom row is for the tower bending moment (TBM). The blue columns in Figure 1 are the DEL with no added turbulence in the simulation and the red columns are the DEL with the added turbulence simulated. The leftmost column of plots give the results for the most upwind WT , the middle column is for the middle WT and the rightmost column is for the most downwind WT. Finally, 3D, 5D and 10D indicates that the results are for simulations with a spacing of 3, 5 or 10 rotor diameters between the WTs.

Turbulence Minimising Control
The simulations are performed with the following parameters of the dynamic turbulence minimising controller L = 5 · 10 11 · I 3 , κ = 5 · 10 −14 , x i = 5 · 10 6 and x i = 1 · 10 6 for every i = 1, 2, 3. The difference in magnitude between L i and κ is due to the fact that the power is in the order of 10 6 [W], while turbulence intensity is in the order of 10 −1 [·] and the gradient of the added turbulence with respect to power is in the order of 10 −8 [W −1 ]. Due to the latter, the gradient of the objective function in the static optimisation is also multiplied with a factor of 8 · 10 13 . The WF demand P dem is 12 [MW] throughout the simulations. Figure 2 compares the mean added turbulence on the two downwind WTs between the open loop even distribution and the turbulence minimising controllers. Table 1 gives the relative reduction in the mean added turbulence when using the turbulence minimising controllers compared to the open loop even distribution. Figure 3 compares the DEL between the turbulence minimising controllers and the open loop even distribution in the layout with 3D inter-turbine distance. Table 2 gives the relative reduction in DEL when using the turbulence minimising controllers compared to the open loop even distribution. In Figure 4, the last 1000 seconds of the time series representing  3D 5D 10D Static Opt mean( i σ add,i ) reduction 5.6% 5.6% 5.68% Dynamic Opt mean( i σ add,i ) reduction 7.9% 8.5% 8.2% Table 1 Table 3 summarises the root mean square (over time) of the reference tracking error of the WF in the simulations, where the values are given both as absolute values and relative to the reference demand.

Simulation Environment
From the model for added turbulence intensity given in (1), we see that the added turbulence decreases with the distance between the WTs. As can be seen in Figure 1 the results fit well with our assumption that fatigue and added turbulence is positively correlated. In particular, as the distance between WTs increases, the turbulence -and also the calculated DEL -decreases. Lastly, it is seen that the most downwind WT has the largest DEL compared to the middle WT. This is also expected since the most downwind WT experiences added turbulence from both upwind WTs. This difference diminishes with added distance between the WTs.

Turbulence Minimising Control
From the results listed in the previous section, we have the following observations. First, from Figure 2 and Figure 4 the reduction in i σ add,i as listed in Table 1 seems solely to be due to the fact that WT 1 is producing significantly less when the controllers are used. This reduces the added turbulence at WT 2 significantly. On the other hand, the reduced output from WT 1 has the consequence that WT 2 increases production which in turn increases the added turbulence at WT 3. Secondly, from Figure 4 and Table 3 we see that the required increase in power production from WT 3 compared to the open loop scenario means that it becomes more difficult to track the farm power demand with the turbulence minimising controls since the wind deficit is the largest at WT 3. This has the consequence that WT 3 more often than in open loop cannot produce the required power and the dynamic controller uses WT 1 (and to some extend WT 2) to compensate frequently. Thirdly, the reduction in DEL on the tower bending moment (TBM) as listed in Table 2 is mostly due to the reduced DEL on WT 1 (see Figure 3) which in turn is due to the reduced output from WT 1 (see Figure 4). However, as can also be seen from Figure 3 (and Table 2) the DEL on the shaft torsion moment (STM) is increased on WT 3 for both controllers. This is most likely due to the fact that WT 3 more often than in the open loop case cannot meet the power reference required by the controller which is likely to produce more actuation from the pitch controller on WT 3. Furthermore, for the dynamic controller the STM DEL is significantly increased on WT 1 due to the fact that it needs to compensate for these reduced output of WT 3. The latter explains the large difference between the additional STM DEL between the two controllers as evident from Table 2. The trends discussed above from Figures 2, 3 and 4 are also present in the two other layouts (5D and 10D) so for simplicity we have left them out of the presentation. In Table 3 it can be seen that the open loop WF is best to meet the farm power demand with the closed loop static optimisation coming second. This is also due to the fact that for the closed loop scenarios WT 3 more often than in the open loop scenario is unable to meet the required power reference. In conclusion, it is clear that the added turbulence can be reduced by a static optimization and even more by a dynamic. This in general reduces DEL for the tower (TBM) but never for the drive train (STM). This calls for an interpretation. One possibility is that tower and drive train have some fundamental differences.
The tower trust will always be varying with the wind speed. In partial load, the pitch and trust coefficient C T is fairly constant and the trust increases, and varies, with the wind speed squared. In full load, the trust decreases with the pitch which varies with the wind speed. On top of this, the trust will also have changes due to changes in power reference. In contrast, there is a major difference between partial and full load for the drive train. In full load, the generator torque is constant or almost constant (depending on use of constant torque or constant power) and the rotor torque only changes a little with the pitch control by the speed controller to keep speed within limits. Consequently, the drive train torque is fairly constant. However, in partial load maximum C P is used, the pitch is constant, and the drive train torque varies with wind speed squared. In summary, tower fatigue is accumulating in both partial and full load whereas drive train fatigue stems mostly from drive torque, and then generator power, changes in partial load or from power reference changes. Both fatigue loads must then be reduced by reduced turbulence but the drive train load suffers even more from increased power changes due to being in partial load or power reference changes when the power is unevenly distributed between the WTs compared to an even distribution. This case is clearly demonstrated in Figure 4.

Conclusion
The starting point of this paper is the thesis that turbulence minimising control actions can be used to reduce the fatigue on WTs in a de-rated WF. The central premise of this thesis is the fact that the added turbulence from the wake of upwind WTs adds to the fatigue of downwind WTs. This was shown to be true by comparing two simulations: one without the added turbulence model and one including the added turbulence model. Next, two turbulence minimising controllers were designed and tested in numerical simulations. The results showed that even though the controls were able to reduce the added turbulence significantly and the fatigue on the tower bending moment slightly, this came at a high cost in the form of additional fatigue on the rotor shaft. This additional fatigue is more pronounced when the updates in power reference for the WTs are frequent as seen when using the dynamic optimising control. We conclude that reductions in the tower fatigue can be obtained with a turbulence minimising controller. However, it comes at a cost of additional fatigue on the rotor shaft and this additional fatigue becomes high if power reference updates are frequent. Future work on the topic could include a more detailed WT simulation model with additional DOFs to the ones available in SimWindFarm, to reveal a more full picture of how the reduced turbulence affects the fatigue of the WTs. Furthermore, to avoid a large variance in the fatigue between turbines, it could be considered to include minimisation on the variance in the optimisation problem. Lastly, it should be considered to introduce a penalty on control updates as frequent control updates seems to subtract from the benefit of minimising added turbulence.