Multiobjective model predictive control design for wind turbines and farms

Analysis of wind turbines has traditionally focused on either single wind turbines with varying wind dynamics or large wind farms with simple wind dynamics.¹ Extending these results to more complex models quickly becomes computationally infeasible or leads to solutions that are suboptimal. While power extraction can be accurately maximized for the single turbine case, power losses arise when one or more turbines lie in an upstream turbine's wake. In short, increasing the axial induction factor of the upstream turbine reduces the amount of energy that can be extracted downstream. As such, many of the simple models are inadequate for large-scale wind farms.² However, large-scale wind farm analysis demands simple wind models in order to remain tractable.

In this paper we improve upon prior methods on using lexicographic minimizers,³ wherein the authors apply minimization techniques over a set of heuristics as a method of finding control inputs for wind turbine arrays. These techniques are limited to applications involving fixed-pitch wind turbines. We aim to expand these results with two novel improvements: yaw misalignment and far-field wake estimation. Furthermore, we use a higher fidelity wind velocity profile for the wind control horizon prediction and use estimations derived from the wake meandering model. We apply a previously developed coordinated ascent algorithm to obtain coupling coefficients that capture the effect of far-field wake effects,⁴ and extend these same techniques for offset yaw control to improve power generation across the wind farm as a whole.

Individual horizontal axis wind turbines operate by using the aerodynamic force of wind to spin the rotor blades. This rotational energy is then converted in a generator to create electricity. Turbines are often categorized by their rated power, or the rate at which they will produce power when the wind is in an ideal operating region. This rated power differs for different wind turbine designs based on blade length, materials, and several other factors. Individual power output crucially depends on the upstream wind velocity.

As wind passes through a turbine, it creates a wake that decreases the downstream average wind velocity. The faster the spin of the turbine blades relative to the wind speed, the greater the impact on the downstream wake profile. Thus, for wind farms it is important to control upstream turbines in an efficient manner so that downstream turbines are not adversely effected by upstream wake effects.

This coupling behavior between turbines emphasizes the need to study optimization of overall wind farm power output, rather than greedy local turbine control. Increasing the maximum power available from a wind farm allows for a wider range of power demands to be satisfied. The complexity of the problem increases as the number of turbines increases. As such, attempting to solve for an optimal solution often is impractical.⁵ 

Greedy algorithms which are locally optimal and maximize the power of any individual turbine are suboptimal over the entire system. The complex interdependencies between turbines are due to interactions of the wakes among large numbers of turbines. In addition, the unpredictability and uncertainty of the wind adds further complexity. Creating a model that accurately captures the important factors of these arrays without oversimplifying is vital to solving this problem.

Recent work explores reducing the impact of upstream turbines on downstream and overall power output through wake steering, intentionally misaligning the yaw angle of an upstream wind turbine to deflect its wake away from the downwind turbines.⁶ While locally this leads to reduced power output, exposing the turbines downstream to higher wind velocities increases the power output from the farm overall. Balancing wake steering with other objectives, such as responding to varying power demands and spikes in power generation, is important when evaluating the use of and comparing theoretical algorithms to those currently in practice.

Prior research has been limited by the complexity of wind models. While true turbine dynamics are dictated by fluid dynamics, often simpler models are used to make meaningful approximation results. Furthermore, wind prediction models are often vastly complex and dependent on a variety of factors, such as seasonal temperature changes, past weather, and differing surrounding air pressures. The work presented herein uses a simplified model derived from wind tunnel experiments. Applying these techniques directly on higher complexity models still warrants further research.

The remainder of the paper is organized as follows. We begin by describing the aerodynamic models for both wind turbines and wind dynamics in Sec. II. In Sec. III, we then describe the methodology for control design under a set of constraints: minimizing tracking error, peaks on power generation, and wake deficits. Our proposed coordinated ascent algorithms for calculating the multiobjective weights, with or without yaw misalignment, are further discussed. Section IV includes the model predictive control (MPC) formulation of the wind turbine problem and discusses further assumptions necessary for the proof of stability of this control framework. In Sec. V, a variety of simulations of several farm configurations and wind dynamics validate the results.

If we consider a single wind-turbine (see Fig. 1), the power that can be extracted is proportional to the cube of the incoming wind velocity, V and corresponding power coefficient $Cp(λ,β)$⁠; this coefficient directly depends on the pitch angle of the blade β and the $tip−speedratio$ λ. The $tip−speedratio$ is defined as $λ=RΩrV$ where R is the rotor radius and Ω_r is the rotor speed. A visualization of how this power coefficient varies with respect to both the $tip−speedratio$ and β can be seen in Fig. 2. The power of turbine i can be computed as

where $Cp,i(λ,β)$ is derived from wind turbine generator specifications, such as those found in Ref. 8, A_rotor is the cross-sectional area swept by the rotor blades, and ρ is the corresponding air density. Taking physical constraints into account, the power generated by wind turbine i will be determined by

where $Pav,i$ represents the maximum available power and is defined as

Here, $Cp,max$ represents the maximum possible value taken by the power coefficient and P_rated is the rated power for a turbine. The rated output is the peak production at a specific (and usually high) wind speed and varies based on turbine specifications. A plot illustrating the free-flow wind speed to power output relationship is given in Fig. 3. In Region 1 (low wind speeds), wind turbines seek to maximize the energy capture, indicating the first term of the minimization; in Region 2 (high wind speeds), the power is limited by the rated power. To reduce the order of our system dynamics, the remainder of this paper will assume we have a fixed-pitch controller and will presume β is 0. Modern commercial turbines often operate at a variable pitch during rated operation; the results herein can primarily be applied during below-rated operation when turbines traditionally employ fixed-pitch control. Often, turbines are additionally characterized using their $torque coefficient$ defined as

Frictional losses, blade size, and other turbine design selections lead to nonlinearities in $Cp(λ,β)$⁠. Prior work¹⁰ using sampled wind turbine data provides an approximation of $Cp(λ,β)$ as

where coefficients $c1,…,c9$ are computed via linear regression using data given power coefficient curves.⁸ The expression $Cp(λ,β)$ given in (5) is based on grid energy measurements and uses power losses to estimate the nonlinearities.¹¹ In practice, $Cp(λ,β)$ is often reparameterized to lessen complexity. The resulting value is known as the axial induction factor of turbine i, represented as $ui,$ and acts as the control input to the turbine. The axial induction factor is the ratio of the velocity of the air passing through the rotor when compared to the free-stream velocity. While exact behavior of u_i is dependent on turbine characteristics, individual turbines are constrained by the Lanchester–Betz–Joukowsky limit, a conservative upper bound derived from actuator disk theory that is based on the available amount of energy that can be extracted from a turbine, regardless of turbine construction.¹² After reparameterizing both $Cp,i$ and $Cp,max$ in terms of the axial induction factor, the resulting relationships become

and

From analysis of (6), it can be seen that optimizing over all u_i is unnecessary, and rather constraining between 0 and the local maximum at $ui=13$ is sufficient.²⁹ For the remainder of the work, we impose the constraint $0≤ui≤13$⁠.

The value of $1627$ is a loose upper bound for $Cp,max$ derived directly from (6). In practice, this value is conventionally defined as the maximum designed power coefficient for a given turbine. This value often ranges between 0.4 and 0.5, with the NREL 5 MW benchmark turbine having a maximum value of $Cp,max≈0.487$ for β = 0 and $λ≈7.69$ as seen in Ref. 9. Furthermore, there is an additional power loss based on the yaw error (otherwise referred to in this work as the yaw misalignment angle) γ to the incoming wind. The wake adjustment given γ becomes

where η is a parameter between 1 and 2. For the remainder of this work, we approximate η as 1.88 according to Ref. 13. This value captures the nonlinear loss in power for operating a turbine offset to the direction of incoming wind. Traditional turbine operation aims to minimize the yaw misalignment angle γ, thereby lessening the loss in power of a single turbine. It has been shown, however, that intentionally misaligning an upstream turbine can lead to power improvements over the entire farm.^13,14 While the misaligned turbine itself generates suboptimal power, the wake may be steered so as to not directly impinge upon a downstream turbine.

We model the interactions between turbines using the Park wake model. For large-scale wind farms, the Park wake model is one of the most prevalent in the existing literature because it provides a verified simplification of the problem. This model oversimplifies the wind dynamics, but there exist extensions to improve model quality.

Given uniform wind with magnitude, $U∞$ the Park wake model states that each turbine affects the wind velocity downstream, V_i, by a diminishing multiplicative factor, given by

with

where $Δd$ is the horizontal displacement from the turbine in question, K is the so-called roughness coefficient, D is the diameter of the wave propagation cone, and u_i is the axial induction factor. At a known location in the intersection of wakes $1,…,i$⁠, the aggregate wind velocity is determined by

where p_ij is the proportion of turbine j lying in the ith wake, measured in cross-sectional area between the velocity cone of the two turbines. A visualization of this cross-sectional overlap can be seen in Fig. 4. As the upstream wind changes, turbine wake effects propagate downstream through a wind farm. A time delay results between upstream controller actions and subsequent downstream wake effects that can lead to suboptimal power extraction. We model this behavior by having the wake changes flow through the farm according to the underlying average wind speed. We assume Taylor's frozen turbulence hypothesis. That is, wake effects propagate downstream as a unified front until they interact with downstream turbines in the wind farm.^15,62

Furthermore, turbines with a nonzero yaw angle γ also skew the wakes when operated at an offset yaw angle. Specifically, the edges of the wakes are deformed by the following model:

where the initial wake deflection angle $ξinit(γ)$ is

and the Betz-optimal thrust coefficient C_T is $89$⁠. The numeric coefficients vary based on exact turbine specifications. For further details, see Refs. 13, 14, and 16. A visualization of this wake skew can be seen in Fig. 5. Equations (13) and (14) capture two effects of offset yaw turbine control; (13) describes how the centerline of the wake drifts based on the distance downstream from the turbine and (14) describes how the wake expands. Both of these effects are estimated from the wake meandering model in Refs. 16 and 17.

We address two different dynamic wind models applicable to the multiobjective optimization wind farm problem. In practice, when locations are being explored as potential sites for wind farms, basic wind data are collected and analyzed. This includes, but is not limited to, average wind velocities, wind directions, and maximum wind speeds. This information can be visually represented as a wind rose. One such example for the Lamma Island wind farm is shown in Fig. 6. Given this statistical information, one wind model that has been proposed is based on a Markov chain model. There has been extensive research on using Markov chains as a stochastic method for extracting statistical wind patterns, with work presented in Ref. 19 offering major improvements over earlier models.^20–22 

The process of developing a Markov chain model can be divided into four steps: state definition, transition matrix estimation, state simulation, and wind speeds simulation.¹⁹ 

The range of the wind speed training dataset is discretized into M finite disjoint intervals. A similar number of intervals can be derived for possible incoming velocity angles. These intervals are then categorized into M states by gridding the space between 0 and V_max. Finer gridding yields higher resolution, but also greater complexity. The wind speed data $(v1,v2,…,vM)$ is then converted into a series of states $(i1,i2,…,iM)$⁠, by assigning intervals. For this work, the wind speed was partitioned in 0.5 m/s intervals and the wind direction in 0.25° bins.

The state time series $(i1,i2,…,iM)$⁠. is treated as a homogeneous Markov chain.^38–41 Maximum likelihood estimates for the transition probabilities are used to estimate the underlying distribution; these are based on empirical frequencies obtained from real-world data. In this work, we form our estimators (transition probabilities) by taking the total number of occurrences of each state transition relative to this large dataset (i.e., we use the frequency plug-in estimator); this forms an estimated transition matrix $P̂$⁠. The same process of constructing estimators is repeated with wind velocity directions.

Using $P̂$⁠, a series of updated states for velocity and direction for the following time instant are generated.

We assume that the wind speeds in each state follow a uniform distribution within that state for simplicity. Wind speed states are converted into wind speeds by $v¯=Vl+κn(Vu−Vl)$⁠, where V_u and V_l represent the upper and lower bounds of the intervals defined in the state characterization step. Assuming that the speeds obey a uniform distribution, we assign κ_n as a random variable selected from a uniform distribution with range $[0,1]$⁠.

Even with wind information that has no noise, it proves to be computationally complex to derive a model from these vast sets of data. These datasets are often large, with the information provided from the United States Eastern/Western Wind Datasets being 10-min sample interval, time-series wind profile measurements across several years. One generalized simplification is to estimate the wind velocity profile using a Weibull distribution.²³ 

The Weibull distribution leads to a probability distribution function over wind velocities V given by

where k > 0 is known as the shape factor and ranges between 1 and 3 according to Ref. 24. The scale factor, Λ, is an additional parameter used to properly scale the normalized distribution to wind velocities. An example of the Weibull distribution is shown in Fig. 7. We approximate wind dynamics with a Weibull wind model in our simulations in Sec. V.

In this work, we aim to design a set of power reference levels that maximize the total available power in a given wind farm. By taking into account far-field wake effects and wake steering techniques, we obtain a methodology that improves on prior MPC analysis of wind turbine arrays. Furthermore, we impose a set of constraints on this formulation. We aim for the assigned power references to satisfy the power demanded by the grid while still minimizing the peaks in power generation.

For the present analysis, it can be assumed that every wind turbine includes an internal power control strategy that follows the power curve in Fig. 2. In this circumstance, the dynamic behavior from the power set-point to the generated power is modeled as a first order system delivering a power from (2) as follows:

with τ the time constant of the dynamic behavior of a first-order model and $Pr,i$ the set point demand for turbine i demanded by the wind farm controller.

The centralized wind farm system has been shown to be appropriately modeled as a linear system^3,25–27 implying the continuous time dynamical model to be controlled can be given by

where

and 1 represents a row vector with 1 in every element, and 0 represents a column of zeros to make the matrix appropriate dimension; $xt=[Pg,t,ξ]$ is a vector of system states where ξ is an error term used to ensure a zero steady-state tracking error; and $ut$ is a vector of control inputs. Error dynamics control under linear state space modeling has been shown to provide promising results, with only small discrepancies likely caused by linearization.²⁸ The system is then discretized based on a sampling time T_s [for example, to implement a model predictive control (MPC) scheme, as discussed later] and is denoted by

where $k∈ℤ+$ is the discrete time step, $xk=[Pg,k,ξ]$ is a vector of system states where ξ captures error to ensure a zero steady-state tracking error, and $uk$ is a vector of control inputs. The centralized wind farm controller has inputs of power demand, P_dem, and has access to the measurements of the power generated $Pg,i$ and available power $Pav,i$ at each wind turbine, while the objective is to provide inputs in the form of desired power reference levels to each turbine, $Pr,i$⁠.

The basic multiobjective problem is formed by identifying three control objectives in a model predictive control framework.

With this objective, we aim to minimize the difference between P_dem and the total power generated, giving us

We want to minimize peaks in power generation to the outgoing system, specifically the difference between the current time and the prior set of $Pr,i$ of the sum of all turbines giving

This objective helps to mitigate load on the grid.

The sequence of $Pr,i$ over the control horizon can be found by properly setting the axial induction factors u_i that minimize the wind deficits. It is straightforward to show that there exists a one-to-one relation between $Pr,i$ and u_i. Additionally, it has been shown^{4,14,29,30,32} that the most downstream turbines should contribute more to total generated power. We choose to model this as a weighted sum

where the weights $ωi∈[0,1]$ are chosen such that they represent the coupling between turbine pairs.

Prior work^3,26,33,34 has shown that the application of a near-field wake model leads to a selection of weights that satisfy $ω1≥ω2⋯≥ωN$⁠, where the turbines are ordered such that the frontmost turbine (e.g., the turbine closest to the free-stream wind direction) is given index 1, and the rearmost turbine (e.g., the turbine most affected by wakes) is given index N. If there exist a set of frontmost turbines (in the case of a grid, for example), we use symmetry to designate identical values of ω and assign each an index less than any downstream turbine. Applying MPC techniques using this ω_i assignment leads to control inputs that track power demands and improve performance.

We aim to improve upon this weight assignment. The authors in Refs. 4 and 29 show that wake deficits depend not only on turbine placement in the farm, but also turbine characteristics such as blade length, turbine configuration, and free-stream wind velocity. They further show that the power extraction of a large wind farm can be improved by selecting a set of u_i that capture far-stream effects. Specifically, one approach originally introduced in Refs. 4 and 32 is an initialized cyclic coordinated-ascent algorithm, referred to as ICyCA-WF. Rather than arbitrary assigning the weights ω_i, we instead propose ICyCA-WF-ω (Algorithm 1), which uses the output of the coordinated ascent algorithm to derive a set of weights. This leads to a solution with $ωi≥ωN$ for all i, but offers no further guarantees. This weight assignment captures turbine specifications, wind speed data, wind turbine placement, and offers an improvement to maximum available total power $Pav,total$⁠. The result is a set of weights calculated assuming $γi=0$ for all i. If the power demanded is larger than what can be provided under this condition, then Algorithm 2 provides a set of weights in which the turbines are operated at incidence angle γ_i. Prior MPC analysis on wind turbine arrays²⁶ arbitrary assigned ω_i such that it captured the number of upstream wake effects; turbines downstream were given lower values of ω_i. This directly corresponds to applying a near-field estimation of wake effects, where the local measurement of turbulence only encompasses the nearest upstream turbine. While these values are useful for demonstrating optimal control using a temporal version of Bellman's equation for simple wake models, the authors in Refs. 33 and 34 demonstrated the near-field assumption is a poor estimator for true wind dynamics. Both Algorithm 1 and Algorithm 2 use far-field estimation of wake effects, where turbine wake information is used from all upstream turbines. It has been shown that these methods lead to noticeable power improvements across the wind farm as a whole.^4,32

Applying discrete MPC techniques leads to the following optimization problem to solve for a set of axial induction factors at time k.

Problem 1.

where w_j represent weights to scale each of the control objectives to the same order of magnitude; their exact values are listed in  Appendix A.

Given the objectives discussed in Sec. III, we now formulate the wind turbine control problem using a multiobjective MPC framework..^42–60 At each time instant k, our general objective can now be stated as

where $J*$ represents the optimal solution at any given time k over all turbines $i=1,…,N$⁠. To solve this optimization problem, we expand and convexify (27) over the three objectives, leading to the following minimization problem:

for some positive definite matrices P, Q, R. These matrices contain w_i and guarantee that no single constraint dominates the other simply by a matter of scale. Furthermore, by convexifying the problem in this manner, the minimum is the same [$Ji(uk,i,γk,i)=Ji*$]. However, the numerical value of nonoptimal results may differ. Imposing this linear quadratic structure has been used in prior work^35–37 and allows for the construction of closed-form Lyapunov functions. For the sake of notation, we will denote the sum in (28) as $l(u,γ)$⁠.

To formulate our main MPC problem, we now consider optimizing over a finite horizon under some additional constraints,

Here, $T$ represents the finite horizon over which the MPC problem is solved, the constraints on $U$ are dictated by the Lanchester–Betz–Joukowsky limit, and Γ is bounded. Further, we assume that there is a bounded rate of change for the yaw angle, labeled $Γ¯$⁠. We impose this final constraint, as physical turbines in general are not operated beyond a critical yaw angle, where air no longer streams smoothly over the blade's upper surface and due to their size and weight are slow to control. The power loss of a single turbine for offset yaw alignment beyond this value cannot be feasibly recovered by power gains from downstream turbines. In addition, large yaw misalignment angles lead to mechanical loads that can damage turbines. The velocity constraint directly correlates with having all turbines operating in the active region. Combining (26) and (29) provides the full MPC problem,

Problem 2.

The additional constraints imposed are the active region constraint, which implies the turbine is operating and providing at least P_min; the bounded demand constraint, wherein the power request by the grid is no greater than the maximum available power of the farm; and a terminal stage stability constraint. The result of the minimization (30) provides a set of control variables that meet the power demands under the finite horizon $T$⁠.

For a stabilizing solution to exist, the underlying model (22) must be stabilizable with respect to both B₁ and B₂. While this can be checked explicitly from the system realization given in (30), it is important to note that individual turbines are often stopped because of safety constraints. The system realization of individual wind turbines in terms of λ and β approximates an Linear Parameter Varying (LPV) model.²⁵ The stability of an individual wind turbine depends on the operating parameters therein. As such, if the incoming wind velocity is too fast or the blades are spinning beyond the safety regime, the turbine may approach an area of unsafe operation. When this occurs, the continuous time first order linear system model (18) no longer accurately approximates turbines operating in the active region, as $Pg,t$ may require an operating torque that will damage the turbine. To ensure the model maintains a stabilizability condition, we assume that the turbines are operating within the active region but implement a safety cutoff if wind speeds are beyond the safe operating threshold or torque of the turbine blades.

We assume that (30) has output equations y_k = Cx_k, with C such that (30) is observable. This is a reasonable assumption as there exist sensors allowing measurement of $λ,Ω,$ and V_i. This means that the axial induction factor can be calculated, which implies P_g is known, which in turn implies our controller scheme can be centralized. We will first assume full state information, i.e., C = I. Application of this assumption proves difficult, as sensor measurements are noisy due to placement of the anemometer within the wake. Furthermore, the wind speed is nonuniform over the cross-sectional area spanned by the turbine. To address this complication, one can either integrate some form of wind estimation or consider time averages across a large horizon.

By construction $Ĵ*<∞$⁠. Furthermore, we assume that there exists a bounded constant $κ>0$ such that $Ĵ*≤κ||x||22$ for all x built from $uk,i,γk,i∈U,Γ$⁠. This form of $Ĵ$ is chosen specifically since this assumption is always satisfied if the constraints in (30) are polytopic.

If (30) is feasible at time t + T with optimal values $u*,γ*$⁠, then a feasible solution at time $t+T+1$ can be obtained by a timeshift of the prior optimal solution. For example, if $∑τ=tt+Tl(u(τ),γ(τ))$ has an optimal solution, then $∑τ=t+1t+T+1l(u(τ),γ(τ))=∑τ=tt+T+1l(u(τ),γ(τ))−l(u*(t),γ*(t))$⁠. This is equivalent to saying that the feasibility is maintained, regardless of the starting time t and the final cost is dependent on the interval length T.

The proof of stability of (22) follows a traditional Lyapunov argument, wherein a general form of a Lyapunov function is written and shown to satisfy the input-to-state stability criterion. The complete details are given in  Appendix B.

Three different wind farm layouts are considered to test the proposed control strategy. First, we consider a simple 3-turbine string, then a 7-turbine string (to accentuate the effects of wake steering). The 5 MW NREL benchmark turbines are used and spaced 630 m apart in the x direction. For the Park wake model, the coefficient of air density ρ is set to 1.225 (kg/m³) and roughness coefficient K = 0.04. The AEOLUS SimWindFarm (SWF) Simulink toolbox was used for simulating the wind speed at wind farm grid points in two dimensions.³¹ Wake effects within the wind farm are modeled according to the dynamic wake meandering model for given ambient turbulent intensity and wind speed direction.

The dynamic model (18) was discretized using a sampling time T_s of 0.01 s to form the discrete model (22) and a prediction horizon $T$ was chosen as 4 s. Each wind turbine was modeled with the time constant (τ) set at 0.08 s. Yaw dynamics were kept slower than the torque dynamics.

A simple turbine configuration is that of the 3-turbine string, shown in Fig. 8.

We test this configuration of turbines with two different wind configurations. Since there are a few turbines in the string, the effects of wake steering are negligible and γ is fixed to 0.

We assume a constant wind profile of 10 m/s traveling from left to right and set a simulated operation run time of 200 s. For this simulation, the goal is to show consistent power tracking of the demanded power of the farm as a whole. The simulation results can be seen in Fig. 9. $Pav,tot$ represents the maximum total available power. This term varies as wakes travel through the farm, decreasing as the wakes negatively affect downstream turbine power extraction and increasing when it is beneficial for upstream turbines to reduce their axial induction factors to improve the farm as a whole. The power demanded P_dem was artificially constructed, so that it is less than the total available power of the farm and does not contain large peaks of power demand. P_g represents the power generated by the farm from the controls resulting from (30). As the wind front passes through the farm, the wind at downstream turbines is impacted due to turbulence and wake effects.

The average difference across the simulation between P_dem and P_g is 0.6%.

We assume a Weibull wind profile with average wind speed 10 m/s, but which allows variations over the time horizon. For this simulation, the goal is to show consistent power tracking of the demanded power of the farm as a whole despite variations in wind speed. The wakes, modeled by the Park wake model (9)–(10), travel through the farm at 10 m/s. Here, P_dem is fixed at a constant value (3 MW), and has a simulated run time of 160 s, so that the wake effects from all upstream turbines impact the downstream turbines. The average error between P_dem and $Pg,total$ across the simulation was 4.12%, with most losses occurring due to dynamic wakes caused by the variable incoming wind.

We now investigate the 7-turbine offset string, shown in Figs. 10 and 11. Using traditional approaches wherein γ = 0, there will always exist wake interference from upstream turbines. This leads to a velocity degradation for downstream turbines, which means the amount of available power decreases.

We assume a constant wind profile of 10 m/s. For this simulation, we aim to show an improvement in the maximum power available in a farm. To do this, we compare the maximum power available in the farm and note the difference when wake steering is implemented. The minimization in (30) is solved in MATLAB. The baseline value is set by fixing γ = 0, implying no wake steering, and running ICyCa-WF-ω to obtain a set of control variables that would maximize power. Applying these values gives a vector of powers, $Pg,NoWakeSteering$⁠. For the same turbine configuration, ICyCa-WF-$ωγ$ is used to calculate a set of control variables and yaw angles, which correspond to $Pg,WakeSteering$⁠. The results can be seen in Fig. 12, in which the ratio between the two powers is plotted, if the resulting ratio is less than 1, the turbine provided less power when wake steering was applied. Conversely, if the value $Pg,WakeSteering/Pg,NoWakeSteering$ is larger than 1, then the maximum power available at that turbine improved. By intentionally misaligning the yaw of the frontmost turbine, it provides less power. However, all other turbines improve significantly, leading to the 7-turbine farm having a power improvement of approximately 14%.

In this paper, we present a model predictive control framework for controlling wind turbine arrays. Under a set of underlying assumptions, a method for assigning turbine controls, more specifically axial induction factors, is shown via simulations to validate the results. A proof which demonstrates guaranteed stability of the MPC problem concludes the work.

The next stages of this work will build on the preceding results but construct a dynamic method to handle nonconverging and rapidly varying wind dynamics.

We also seek to extend upon some of the limitations of the work. For instance, we assumed full-state information. This assumption can be relaxed and an observer-based controller designed. The effect of yaw misalignment on power production does include several tradeoffs, most notably mechanical fatigue. The optimization in this work does not take into account these negative effects, and thus introduces large yaw misalignment to maximize power by wake steering. Future work will explore modeling this fatigue and incorporating these dynamics in the MPC framework. Furthermore, while the work presented herein focuses primarily on a fixed-pitch controller, many turbines operate with variable-pitch controllers. An extension to the turbine model and subsequent MPC controller design is the topic of continuing research.

We wish to acknowledge the financial support provided by NSF Grant No. CNS 15-44953, the Coordinated Science Laboratory, and the Department of Industrial and Systems Engineering, both at the University of Illinois at Urbana-Champaign.

The data that support the findings of this study are available from the corresponding author upon request.

As the scales of $J1,J2,J3$ differ, they are normalized with a set of weights so that the objectives are on the same scale, that is on the order of 10⁵ W. Specifically, we use $w1=10, w2=5$⁠, and $w3=10 000$⁠.

From Assumption 4 and (30), we have

This holds $∀t$ since t is arbitrary; we will choose t = 0 for the remainder of the proof. Then

Under Assumptions 1 and 2, there exists a matrix $P̂≻0$ such that the Lyapunov Function $V(x)=xTP̂x$ satisfies the following inequality:⁶¹ 

Theorem 1. $Ĵ*$ is a stable solution to (30)

Proof. We define $V̂=Ĵ*(x)+βV(x)$ for some $β>0$⁠. Showing $V̂(x)$ is a Lyapunov function is equivalent to showing the stability of $Ĵ*$⁠.

Step (1) $V̂(x)>0 ∀x≠0$⁠.

Rather than show this directly, we demonstrate that $V̂$ is both lower and upper bounded by two different quadratic functions in x that both pass through 0 at x = 0.

We obtain the lower bound by considering the definition of $V̂$ and the fact that $Ĵ≥0$⁠. Specifically, $V̂≥βV(x)≥βλmin(P̂)||x||22.$

We obtain the upper bound by noting that $V̂(x)=Ĵ*(x)+βV(x)≤Ĵ*(x)+βλmax(P̂)||x||22$⁠.

Using Assumption 3 we then have $V̂(x)≤κ||x||22+βλmax(P̂)||x||22≤(κ+βλmax(P̂))||x||22.$

Thus $V̂(x)$ is both upper and lower bounded by a pair of quadratic functions, both of which exclusively are 0 at x = 0 as long as $λmin(P̂)≠0$ and $λmax(P̂)≠0$⁠. This is satisfied by construction.

Step (2) $V̂(xk+1)−V̂(xk)<0$⁠.

By definition, $V̂(xk+1)=Ĵ*(xk+1)+βV(xk+1) and V̂(xk)=Ĵ*(xk)+βV(xk)$⁠. Thus

Using (B2) and (B4)

Without loss of generality we suppose $u*(0)=0,γ*(0)=0$ and $J1*,J2*,J3*=0$⁠. Then,

We then define our parameter β as

By this choice of β the P, Q, R terms dominate. Then

Where $ζ<0$ and represents the difference between the J and c terms, and thus the desired result holds.

Step (3) $V̂(0)=0$⁠.

This holds by construction.