A Methodology for Robust Load Reduction in Wind Turbine Blades Using Flow Control Devices

Decades of wind turbine research, development and installation have demonstrated reductions in levelized cost of energy (LCOE) resulting from turbines with larger rotor diameters and increased hub heights. Further reductions in LCOE by up-scaling turbine size can be challenged by practical limitations such as the square-cube law: where the power scales with the square of the blade length and the added mass scales with the volume (the cube). Active blade load control can disrupt this trend, allowing longer blades with less mass. This paper presents the details of the development of a robust load control system to reduce blade fatigue loads. The control system, which we coined sectional lift control or SLC, uses a lift actuator model to emulate an active flow control device. The main contributions of this paper are: (1) Methodology for SLC design to reduce dynamic blade root moments in a neighborhood of the rotor angular frequency (1P). (2) Analysis and numerical evidence supporting the use of a single robust SLC for all wind speeds, without the need for scheduling on wind speed or readily available measurements such as collective pitch or generator angular speed. (3) Intuition and numerical evidence to demonstrate that the SLC and the turbine controller do not interact. (4) Evaluation of the SLC using a full suite of fatigue and turbine performance metrics.


Introduction
With advancements in wind turbine technologies, turbine sizes have increased significantly in the last few decades, and large rotors are becoming increasingly important for utility-scale wind power. This is because larger turbines exhibit relatively lower electricity generation costs due to the higher wind speeds and the larger swept areas. However, larger turbines suffer from more fatigue damage due to unsteady aerodynamic loads [1]. Mitigating these loads using conventional turbine design approaches makes them heavy and expensive.
A well-known method to reduce the fatigue loads on the turbine blades is individual pitch control (IPC) [2][3][4]. IPC generates pitch angle commands for the pitch actuator on each blade, based on the feedback of dynamic loads. IPC can be implemented using different axes transformations like Clarke transformation [5] and multi-blade coordinate transformation [6][7][8] (also known as Coleman transformation). These transformations perform projections onto a set of orthogonal axes for both the measured loads and pitch commands. IPC can also be designed such that individual blades have separate control loops that utilize the local load on the particular blade to generate independent pitch commands for each blade without communicating with each other [9,10]. This approach is known as single-blade control or individual (i.e., independent) blade control. Reference [11] showed that the two transformation-based approaches and individual blade control are fundamentally similar and result in similar load reduction. IPC can be implemented in a feedback-feedforward framework where the feedback control loop utilizes sensed load while feedforward control utilizes wind estimates to improve disturbance rejection [12,13]. The wind estimate can be provided by a Kalman filter or LiDAR, with the latter providing true feedforward control. Load reduction using IPC has also been explored with additional local flow sensors [14], using advance control design techniques like model predictive control [15,16] and to reduce both fatigue and ultimate loads [17].
Although IPC results in a reduction in loads, its use may lead to excessive wear of the pitch actuators [18]. This concern grows with turbine size due to longer and potentially heavier blades. Thus, on-blade, active load control devices offer an attractive alternative to increase rotor size while disrupting the growth of turbine capital costs. An active onblade load control system could mitigate loads before they develop, which would allow turbines to be built with less material and at a lower cost. Alternatively, blade lengths can be extended for constant load envelope, leading to increased annual energy production (AEP) and reduced levelized cost of energy (LCOE). Active load control solutions have been extensively researched, including trailing edge flaps, microtabs, shape changing blades, and others [18,19]. Among these actuators, trailing edge flaps have been studied the most [20,21] and compared with IPC [22]. They have also been tested on a full-scale wind turbine [23]. However, some of the challenges with these approaches lie in their mechanical complexity, costs, and maintainability in real-world conditions.
In this research, we focus on the development of a new load control system based on a Dielectric Barrier Discharge (DBD) plasma actuator ( Figure 1) due to Arctura [24], which is mechanically simple with no moving parts [24]. We are part of a team [25] developing a controllable Gurney flap (GF) for wind turbine blades using DBD plasma actuators. Our team seeks to integrate actuator development, control systems and design methods to determine the greatest potential impact of the technology on reducing LCOE through a combination of modeling, design optimization and testing (wind tunnel and field rotor demonstration). In the current paper, we describe a specific methodology for control system design to support the development of wind turbines with controllable GFs. This methodology is comprised of the following building blocks:

1.
Modeling the effect of the DBD plasma actuator as a dynamic change in local lift coefficient at the location where the actuator is placed. We refer to this modeling assumption as sectional lift actuation (SLA). 2.
A control system architecture based on the well-known multi-blade coordinate transformation [6] to reduce dynamic loads of the blade-root flapwise bending moments. This control system, which drives the sectional lift actuators, is coined sectional lift control (SLC).

3.
A feedback control algorithm designed to achieve load reduction and robustness to both (a) parametric uncertainty in the design model and (b) unmodeled dynamics. This paper has two major contributions. Firstly, the methodology for design of the control system algorithm (Section 4). This methodology uses system identification tools to determine a dynamic model appropriate for controller design, and then a method for robust control synthesis in the frequency domain to obtain the control algorithm. Secondly, the application of this methodology to the design of a sectional lift controller for a 3.4 MW wind turbine. This 3.4 MW turbine is based on the IEA 3.4 MW reference turbine [26] but features a re-designed rotor and a re-tuned baseline turbine controller that accounts for the presence of the Gurney flap [27]. The designed SLC achieves reductions in blade-root flapwise damage equivalent loads ranging from 4% to 8% in below-rated wind speeds, to 8% to 12% in above-rated wind speeds with no adverse effect on other wind turbine loads, power production and pitch actuator activity (Section 5). The resulting SLC also reduces blade-tip deflections. We demonstrate, by example, that a single robust SLC designed for an above-rated wind speed suffices to reduce loads for all wind conditions. We give evidence that this property likely generalizes to other turbines with similar dynamic behavior. This paper is organized as follows. Section 2 provides the specifications of the 3.4 MW turbine on which this study is conducted. The model of the SLA and its implementation on the turbine are discussed in Section 3. Section 4 describes the procedure followed to develop the SLC. This involves establishing the objective and the evaluation criteria of the control system, identification of a linear system for control design, study and evaluation of the identified systems and finally, synthesizing the SLC based on the previous steps. The controller is then evaluated using the industry standard [28] and the impact of the SLC on various loads, blade-tip deflection, as well as on turbine performance is presented in Section 5. Conclusions are given in Section 6.

Turbine Specifications
The IEA-Task 37, 3.4 MW wind turbine [26] is re-designed and used as a test bench for this study. The wind turbine is an IEC class 3A, land-based wind turbine with rated electrical power of 3.37 MW. This size turbine is the "sweet spot" of the land-based wind market, which is focused on 3 to 4 MW turbines. This makes this size turbine very relevant for both new turbine installations and for retrofit applications of today. The methods described also apply to the design of larger turbines, such as those being installed offshore, where load reduction can impact both turbine and balance of station costs.
The OpenFAST model of the turbine has been made publicly available by NREL. The OpenFAST model of the wind turbine is modified and implemented in the simulation software SIMULINK [29].
We now discuss the overall turbine design methodology for the present study, which examines a controllable Gurney flap (GF) based on DBD plasma actuators. A detailed description of the complete design process and implementation for this 3.4 MW turbine is presented in a companion work [27] and summarized here. The GF is placed on the pressure side of the blade near the trailing edge, as shown in Figure 2. As a result, the local lift coefficient will increase due to the presence of the GF. Therefore, the design process begins with an aerodynamic re-design of the 3.4 MW blade based on computed sectional lift and drag accounting for the GF. The effects of the new aerodynamic design are then propagated to the re-design of the blade structure and re-tuning of the baseline turbine control system, where iterations take place until structural requirements for the blades are satisfied while the turbine baseline controller is re-calibrated with the new GF-based aerostructural design. The specifications of the re-designed 3.4 MW turbine with controllable GF are presented in Table 1 and are detailed further in [27].

Sectional Lift Actuators
For this study, we utilize sectional lift actuators (SLAs) as a proxy for plasma-based flow actuators. SLAs have been used before to investigate load mitigation in wind turbines using active flow control devices [7,30]. The SLAs are assumed to be capable of modifying the local lift coefficient at a section of the wind turbine blade. Plasma-based flow actuators have limited authority. In this paper, we assume that the SLAs can change the sectional lift coefficients within the range −0.2 to +0.2 from the baseline value. These limits on the sectional lift coefficient are within reach of future plasma technology [7,25]. The SLAs are assumed to be effective for angles of attack ranging from −15 degrees to +20 degrees. Figure 3 describes the effect of SLAs on the lift coefficient at a particular blade section. The SLAs provide a means to modify the local aerodynamic force at specific locations on the turbine blades, which, in turn, can be used to modify the bending moment across the blade and, particularly, at the blade root.
A single SLA is implemented on each of the three rotor blades. The SLAs are located in the 24% outer span of the blades. Specifically, the SLA is placed at a blade span of 47.9 m and extends to the blade tip at 63 m. Figure 4 shows the actuated section on a turbine blade.
OpenFAST is a turbine simulation tool developed by NREL [31]. The OpenFAST model of the turbine is modified to implement SLAs on the blades. The OpenFAST code is modified to enable it to accept three sectional lift commands as external inputs (one for each blade). The AERODYN module of OpenFAST, which computes the aerodynamic forces on various components of the wind turbine, is modified to change the local lift coefficient at the blade elements subject to the SLA. The modified OpenFAST model of the turbine is implemented in the simulation software SIMULINK along with a feedback controller. The control system, described in Section 4.5, provides the SLA commands u 1 , u 2 , u 3 to the rotor blades.

Sectional Lift Control
The purpose of the load control system is to reduce blade fatigue caused by out-ofplane dynamic loading due to the combined effects of rotor tilt, wind shear, tower shadow, yaw misalignment, and turbulence [32]. Blade fatigue correlates highly with the blade bending moments; thus, it is generally accepted [33] that reducing such moments translates into a reduction of the fatigue damage equivalent loads as defined in the IEC standard [28].
Sectional lift control (SLC) seeks to reduce the dynamic flapwise bending moments of the blades by proper modulation of local lift coefficient along the blade span. These dynamic (unsteady) moments are the primary contributor to blade fatigue due to out-ofplane loading. The moments at the blade root are of particular interest as this location usually experiences larger amplitude variations than any other location on the span of the blade. Therefore, the objective of the SLC is to reduce the fluctuations of blade-root flapwise bending moments in the presence of wind disturbances by proper control of the on-blade SLAs. Figure 5 depicts the feedback configuration used to calculate SLA commands using the measurements of the blade-root flapwise bending moments in real time.

167
The SLC is designed to reduce the amplitude variations of the blade-root flapwise 168 bending moments. As explained in Section 4.5, emphasis is placed on reducing the 169 bending moments at the rotor frequency, the so-called one-per-revolution (1P) and its  [35]. A similar approach has been successfully implemented in prior load and vibration 178 control problems for offshore wind turbines [36,37]. Details of the design methodology 179 are given in Section 4.6.

180
The effect of SLC on the turbine is evaluated using damage equivalent loads (DEL)

SLC Design and Evaluation Criteria
The SLC is designed to reduce the amplitude variations of the blade-root flapwise bending moments. As explained in Section 4.5, emphasis is placed on reducing the bending moments at the rotor frequency, the so-called one-per-revolution (1P) and its neighboring frequencies. This control problem is posed such that the objective is to reduce the closedloop sensitivity transfer matrix [34] at the plant output (i.e., bending moments) in a given bandwidth of interest. The design methodology is based on a frequency-domain loop shaping technique, where performance is achieved by "shaping" the magnitude of the appropriate frequency response function. The design is then tested for robust stability and, if necessary, the robustness of the control system is augmented using the H ∞ loopshaping method developed by Glover and McFarlane [35]. A similar approach has been successfully implemented in prior load and vibration control problems for offshore wind turbines [36,37]. Details of the design methodology are given in Section 4.6.
The effect of SLC on the turbine is evaluated using damage equivalent loads (DEL) [38] of various turbine loads. The DELs are calculated using simulation data obtained from OpenFAST under various wind conditions such as uniform wind, wind shear and turbulence. The software MLife [39] provided by NREL is utilized for the calculation of DELs. The loads are classified into primary loads and secondary loads. The criteria also include turbine performance metrics such as power, rotor speed and pitch activity. The evaluation criteria are summarized in Table 2.

Loads and Deflections
Blade-root flapwise bending moments Primary load

Multi-Blade Coordinate Transformation
It is known that loads on the turbine blades have the largest component at the 1P harmonic frequency. Thus, a load control system is designed focusing on the harmonic components of the load at 1P and neighboring frequencies. The multi-blade coordinate transformation (MBC) [6] is used to extract the one-per-revolution (1P) harmonic load components from the blade root moment sensors [8]. The MBC assumes that the three turbine blades are identical to each other and located 120 degrees apart. The MBC transformation requires the azimuth angle of the rotor (θ) for the calculation of harmonic load components. Figure 6 shows the wind turbine augmented with the MBC transformations. components. Figure 6 shows the wind turbine augmented with the MBC transformations. The forward MBC transformation is used to map the flapwise root bending moments of the three blades, M 1 , M 2 and M 3 , from the rotating frame to M cos and M sin , the corresponding loads in fixed frame. The forward MBC transformation can be written as: where θ is the azimuth angle (in radians) of the turbine rotor. Under steady wind 198 conditions, the forward MBC transforms the 1P harmonic loads in the rotating frame 199 into constant loads in fixed frame.

200
The right inverse of the forward MBC transformation, denoted as "Inverse MBC" block in the Fig. 6, is utilized to calculate the sectional-lift actuator commands for the three blades, i.e. the rotating frame signals u 1 , u 2 and u 3 , by transforming the control commands in the fixed frame u cos and u sin . The fixed frame control commands are calculated by a feedback controller described in Section 4.6. The inverse MBC transformation can be written as follows: The inverse MBC transformation commands the sectional lift actuator of each blade with 201 1P cosine and sine carrier signals modulated by the control signals u cos and u sin . When 202 these latter signals are steady (constant), and the rotor frequency is constant, the rotating 203 frame control signals u 1 , u 2 , u 3 are pure harmonics at 1P. Otherwise, the rotating frame 204 control signals will have spectral content at frequencies other than 1P as well.

205
The configuration shown in Fig. 6 is useful for developing an SLC system that The forward MBC transformation is used to map the flapwise root bending moments of the three blades, M 1 , M 2 and M 3 , from the rotating frame to M cos and M sin , the corresponding loads in a fixed frame. The forward MBC transformation can be written as: where θ is the azimuth angle (in radians) of the turbine rotor. Under steady wind conditions, the forward MBC transforms the 1P harmonic loads in the rotating frame into constant loads in a fixed frame. The right inverse of the forward MBC transformation, denoted as "Inverse MBC" block in Figure 6, is utilized to calculate the sectional-lift actuator commands for the three blades, i.e., the rotating frame signals u 1 , u 2 and u 3 , by transforming the control commands in the fixed frame u cos and u sin . The fixed frame control commands are calculated by a feedback controller described in Section 4.6. The inverse MBC transformation can be written as follows: The inverse MBC transformation commands the sectional lift actuator of each blade with 1P cosine and sine carrier signals modulated by the control signals u cos and u sin . When these latter signals are steady (constant), and the rotor frequency is constant, the rotating frame control signals u 1 , u 2 , u 3 are pure harmonics at 1P. Otherwise, the rotating frame control signals will have spectral content at frequencies other than 1P as well.
The configuration shown in Figure 6 is useful for developing an SLC system that reduces blade-root loads in the rotating frame at 1P and its neighboring frequencies; that is, reduce the blade-root loads in the frequency range 1P ± ω o , with ω o < 1P. In order to reduce the 1P harmonic loads using the MBC transformations, one would need to design a control system to lower the fixed-frame components of the load in the frequency interval [0, ω o ]. Recall that the MBC transforms a 1P rotating-frame frequency component to a 0 rad/s fixed-frame frequency component [8,11].

Design Model via System Identification
A model for the dynamic system with fixed frame input [u cos u sin ] and fixed frame outputs [M cos M sin ] is obtained for the design of the SLC. This system is shown in Figure 6.
To determine a mathematical model for this system, it is assumed that the wind speed is constant and uniform over the rotor disk. Under this assumption, the open-loop plant model may be assumed (approximately) time invariant, albeit depending on the constant wind speed V impinging on the rotor.
The primary focus of the work described in this article is on designing and utilizing SLC in above-rated wind speeds, when the loads are relatively large. The rated and cut-out wind speeds for the 3.4 MW are 9.8 and 25 m/s, respectively [26]. The design model for developing a robust SLC solution is obtained at V = 18 m/s, which is close to the midpoint between rated and cut-out wind speeds. A linear state-space model of the open-loop plant at V = 18 m/s is obtained using the nonlinear simulation data from the modified OpenFAST model of the turbine. The simulations assume uniform wind without any shear or turbulence.
To identify the plant model, it is necessary to excite the plant with appropriate inputs and record the response. Filtered pseudo-random binary signals (PRBS) are chosen as the input to the open-loop plant in Figure 6. A second-order low-pass Butterworth filter with a 3-dB bandwidth of 20 rad/s given by is used. The frequency response of the low-pass filter is shown in Figure 7. Note that the largest 1P frequency of the rotor is 1.22 rad/s, which corresponds to the rated rotor speed. Thus, a filter with 20 rad/s bandwidth is suitable to identify the relevant openloop dynamics. Two independent PRBS signals are filtered by G LP and the resulting signals applied at the input [u cos u sin ] of the open-loop plant in Figure 6. The OpenFAST simulations are run for a total of 3000 s, and the first 500 s (98 revolutions) are discarded to remove any initialization effects, resulting in 2500 s (490 revolutions) of usable data. A time step of 0.01 s is chosen for the simulation, which allows estimation of frequency responses up to 50 Hz or 314 rad/s, theoretically. Figure 8 shows the filtered PRBS inputs signals. The left plots show a small interval (15 s) of the time series signals at the two plant inputs. The corresponding plots on the right show the magnitude of the Fourier transform of the signals. The magnitude of the transfer function of the low pass filter G LP is shown with a red line. This magnitude is calibrated to match the signal spectrum at low frequency. It can be seen that the frequency content of the signals is focused in the range of 0-10 rad/s and slowly tapers off at higher frequencies as desired, following the magnitude of the low-pass filter.  Table 3 shows the structural degrees of freedom enabled in the OpenFAST model for system identification. Other degrees of freedom and aerodynamic behavior options like wind shear and tower-blade interaction are also relevant to the control problem but are not enabled in system identification. These latter effects are considered unmodeled disturbances/dynamics and are handled with the robust control method selected for the design of the feedback controller. Table 3. OpenFAST degrees of freedom enabled for system identification. First side-to-side tower bending-mode DOF

TwSSDOF2
Second side-to-side tower bending-mode DOF The outputs [M cos M sin ] are logged and the prediction error method (pem function in MATLAB [40]) is used to identify a state-space model from the simulation data. The 2500-s window of usable data is divided into two sets, a 2000-s long "training" dataset and a 500-s long "validation" dataset. The training dataset is used to obtain models of varying orders, which are evaluated for best fit to the validation dataset. The model mismatch is measured with the normalized root mean square error (NRMSE), defined by where y v is the time series validation data, y m is the time series predicted by the model, and the 2 norm of the time series x is defined by The NRMSE for models of different orders, for the case of V = 18 m/s, is shown in Table 4. The model with n = 6 states is selected as it has the smallest NRMSE with the validation data. Further insight into the model's ability to capture essential dynamics at low frequency (up to 3-5 rad/s, approximately) can be obtained from Figure 9. This figure shows the frequency response of the transfer function from the input u cos to the output M cos (see block diagram in Figure 6). The frequency response calculated using the sixth order model shows good fit of the empirical transfer function obtained from the time-series data. The empirical transfer function is obtained using the etfe command provided by the System Identification Toolbox [40] in MATLAB. We define the four entries of this matrix as follows: The four entries of the transfer matrix G(s) are shown in Appendix A. Figure 10 shows the Bode plots of the entries of G(jω). It follows from the figure that G 11 ≈ G 22 and G 21 ≈ −G 12 , for frequencies ω ∈ [0, 10] rad/s. In addition, in this frequency range, we have the imaginary part of G 11 (jω)G * 21 (jω) as approximately zero. As a result of this structure, the frequency response of G(jω) is essentially a rotation matrix modulo of a non-unitary magnitude. This structure implies that the two singular values of G(jω) approximately satisfy σ 1 (G(jω)) ≈ σ 2 (G(jω)) ≈ |G 11 (jω)| 2 + |G 21 (jω)| 2 (7) in the interval ω ∈ [0, 10] rad/s. This, in turn, implies that in this frequency range, inverting the frequency response matrix is a fairly robust operation as its condition number is nearly one. Figure 11 depicts the two singular values of the frequency response matrix G(jω). The Euclidean norm of the first column of the frequency response matrix G(jω), defined in the right-hand side of Equation (7), is also shown with a red dashed line. The condition number (ratio of maximum to minimum singular value) is shown in Figure 12. These data suggest that inverting the frequency response matrix does not pose a challenge. The inverse of the frequency response matrix is used in the development of the SLC algorithm presented in Section 4.6. More specifically, this transfer matrix is inverted at ω = 0 to decouple the control loops at low frequencies. Magnitude (kN.m) Figure 11. Singular values σ 1 (G(jω)) and σ 2 (G(jω)) (solid blue) and the norm of the first column of frequency response |G 11 (jω)| 2 + |G 21 (jω)| 2 (dashed red).

Variation of Design Model with Wind Speed
Design models for wind speeds other than V = 18 m/s may be obtained using the system identification procedure described earlier. This section discusses the main properties of the models obtained at various wind speeds and the implications of any model variations on the design of the SLC algorithm for load reduction in above-rated wind conditions (Region-3).
The Bode plots of the four individual transfer matrices are shown in Figure 13 for a range of wind speeds in Region-3.
The singular values of the frequency response matrix G(jω) are shown in Figure 14. With the chosen scales, the frequency responses are clustered together. Note that the frequency responses change slope at the same corner frequencies despite wind speed variations. Thus, for practical purposes, the linear model does not appear sensitive to wind speed variations in Region-3.
A plausible explanation for this behavior is that in above-rated wind conditions, the turbine's pitch controller is engaged to regulate the rotor angular speed to its rated value. The presence of such feedback loop will reduce the sensitivity to parameter variations of turbine responses, such as the response of blade-root moments to changes in the SLA inputs. As a result, it is natural to expect little variation of the transfer matrix G(s) with the wind speed V. The practical implication of this observation is that scheduling the control law on wind speed, or other more practical measurement indicative of wind speed, might not have major benefits over a fixed controller. Additional data supporting this observation is provided in Section 5.4.

Sectional Lift Control Architecture
The architecture of the SLC is shown in Figure 15. The control algorithm utilizes the outputs [M cos M sin ] of the augmented plant to calculate the sectional lift commands in fixed frame [u cos u sin ] . The control algorithm has four basic building blocks. First, the two loops are decoupled at DC, and at low frequencies, by the Inverse of DC Matrix block, which is implemented inverting the DC-gain matrix G(0) of the identified augmented plant. The Robust Controller block is optional and nominally set to the identity matrix. This block can be designed using H ∞ loop-shaping [35] to increase the robustness of the feedback system if desired. The third block is a MIMO diagonal integrator used to achieve disturbance rejection at low frequency; from DC to ω 0 , approximately. The control signals calculated by the integrators are then modified to take the limited actuator authority of the SLAs into account. This operation is performed by the Saturation Management block in the figure. The block contains a nonlinear scaling, which ensures that the controller commands do not exceed the actuator limits. Back calculation-based integral anti-windup is also implemented as shown by the blocks with gain K aw in Figure 15. The saturation management block outputs the fixed frame actuator signals [u cos u sin ] , which are then converted to the SLA actuator commands u 1 , u 2 and u 3 . Details on the design of each block are given in Section 4.6.

Synthesis of the Control Algorithm
To specify the control algorithm in Figure 15, it is necessary to choose the frequency ω 0 , the robust controller, the saturation management scheme and the anti-windup gain K aw . The selection of ω 0 and, if necessary, the robust controller is done using a loop shaping design procedure [35] whereby the augment plant G(s) is pre-and post-compensated with weights W 1 (s) and W 2 (s), respectively, so that the frequency shaped plant G s (s) = W 2 (s)G(s)W 1 (s) (8) has desired closed-loop properties as measured by the frequency response of its (two) singular values. In our method, these weights are taken as where I denotes the 2 × 2 identity matrix. The weight W 1 can be seen in Figure 15 following the robust controller block, while the weight W 2 is shown following the output of the augmented plant. After a desirable loop shape for G s (s) is achieved by selecting ω 0 , the resulting loop is tested for robust stability. If robust stability is satisfactory, then the robust controller block is taken to be identity, and the design of the linear elements of the control algorithm is complete. Otherwise, a robust controller with transfer matrix K s (s) is obtained by solving the following H ∞ optimization problem where the minimization is done over all controllers that stabilize G s . The controller K s is generally of high-order, and its order can and should be reduced to obtain the final robust controller in Figure 15. Let us elaborate on the selection of the weights W 1 and W 2 . The post-compensator W 2 renders the transfer matrix W 2 G(s) non-dimensional and equal to the identity matrix at s = 0. That is, post-compensation is used to decouple and normalize the two control channels at DC. The calculation of W 2 = G(0) −1 , and its use as a post-compensator, is numerically stable and remarkably robust given that the condition number of G(jω) is close to one at low frequencies (cf. Figure 12). The singular values of W 2 G(jω) are shown in Figure 16. Note that they are exactly one at ω = 0, and remain close to one for low frequencies. Thus, in order to reduce the fixed-frame moments M cos and M sin for frequencies ranging from DC to a low-frequency ω 0 , it suffices to select an integral precompensator W 1 , as shown in Equation (9). For all practical purposes, ω 0 is the cross-over frequency of magnitude for G s . Hence, if robust stability can be achieved with the robust controller K s = I, ω 0 will be the closed-loop bandwidth. That is, ω 0 defines the frequency range where load reduction of fixed-frame moments will occur. Given this argument, we recommend taking 0.2P ≤ ω 0 < 1P, where 1P is the rated turbine rotor frequency. To complete the design method, we need to define the criteria for robust stabilization of the shaped plant G s . We use the normalized coprime factor stability margin [34] defined by Based on experience, we take the minimum acceptable stability margin to be 0.4. That is, if s > 0.4, the robust controller block in Figure 15 is set to be the identity matrix. Otherwise, we design the robust controller from Equation (11) and possibly include a controller order reduction step. Note that our threshold for the stability margin exceeds the recommended value s > 0.33 [41] to increase robustness to modeling uncertainty.
The significance of the normalized coprime factor stability margin in Equation (12) may be understood with the aid of the block diagram in Figure 17. In this figure, ∆ 1 and ∆ 2 are dynamic perturbations that represent inverse multiplicative uncertainty at plant output and direct multiplicative uncertainty at plant input. It is well-known that the block ∆ 1 can capture modeling uncertainty and/or parametric variations, e.g., pole migrations due to changes in operating conditions such as the wind speed. On the other hand, the block ∆ 2 may capture unmodeled dynamics, such as unsteady aerodynamic effects due to sectional lift actuation using plasma-based devices. See [41] for further details on uncertainty modeling. It turns out that the feedback loop in Figure 17 is (internally stable) if and only if Loosely speaking, if s > 0.4, then the feedback system in Figure 17 can tolerate modeling uncertainty as large as 57% at all frequencies without losing stability. Further insight on the interpretation of s follows from the bounds it provides for classical single-input single-output (SISO) stability margins. While this is not our case, these SISO bounds are informative. More specifically, a given normalized coprime stability margin s provides the bounds on the gain and phase margin shown in Equation (14) [34].
For example, a coprime stability margin s = 0.4 guarantees a gain margin of 2.33 and a phase margin of 47 degrees, which are acceptable values for SISO loops. Now, we provide the details of our SLC design for the 3.4 MW turbine. The design model is given by the augmented plan G(s) as identified at 18 m/s wind speed. The choice of the bandwidth ω 0 is based on the reduction in the damage equivalent load of the flapwise root bending moments on the turbine blades. In order to choose the bandwidth, controllers were designed and tested with turbulent wind. The ω 0 value that provides the most load reduction is selected. Based on this analysis, the bandwidth is chosen to be ω 0 = 0.3 × 1P = 0.37 rad/s. Once a numerical value for ω 0 is known, the shaped plant G s in Equation (8) is known and we may compute its stability margin using Equation (12). For our plant, we obtained s = 0.62. This implies that the control system is robust and H ∞ loop-shaping based modifications are not needed. Thus, the final controller is a MIMO integral controller of the form Figure 18 shows the singular value plots of the open-loop transfer matrix L(s) = G(s)K(s) for the control loop broken at plant output, the closed-loop sensitivity S(s) = (I + L(s)) −1 and complementary sensitivity T(s) = L(s)S(s). Load reduction is expected to occur for frequencies up to 0.4-0.5 rad/s because the sensitivity is less than one in this range of frequencies. The peak value of the sensitivity is S ∞ = 1.11, and the peak value of the complementary sensitivity is T ∞ = 1, well below recommended values [34].
The linear component of the control algorithm in Figure 15 was designed using the plant model identified at a wind speed of 18 m/sec. The same design process can be used to design linear controllers for plants at other wind speeds. The controllers thus obtained can be utilized to obtain an SLC scheduled on wind speed. Section 5.4 demonstrates the results of this approach. No significant benefits are observed with the scheduled SLC in terms of load reduction or stability margin or performance to justify the added complexity of scheduling control parameters on wind speed. Therefore, the results presented in the rest of the paper are based on the controller designed at 18 m/s, as shown in Equation (15).

Saturation Management
In this paper, we assume that the actuator limits are ±0.2. These limits on the sectional lift coefficient (i.e., actuator) are within reach of future plasma technology [7,25]. Following Reference [42], nonlinear scaling is applied to the control signals calculated by the controller, as shown on Figure 15. The purpose of the scaling is to modify the output of the controllers so that they remain within the actuator limits while providing suitable load reduction. Denote the inputs to the nonlinear scaling block as v cos (t) and v sin (t), while the outputs are u cos (t) and u sin (t). The scaling can be described as follows It should be noted that the scaled inputs always have the following property.
The constraint shown in Equation (17) on the fixed frame control signal implies that the corresponding control signals in the rotating frame are always within the range of ±0.2, which satisfies the actuator constraint. This scaling method is non-conservative (i.e., it uses the full actuator authority) because the magnitude of each actuator command satisfies max |u i | = u 2 cos + u 2 sin . Furthermore, scaling the control signals in the fixed frame, prior to applying the MBC transformation to obtain the actual actuator commands, ensures that at every time instant, the sum of all the actuator commands is zero, which is necessary to avoid undesirable changes on the collective lift on the blades.
When the output of the scaling function saturates, the system in Figure 15 becomes open-loop. The two integrators in the controller would windup unless this saturation is detected and acted upon. A simple back-calculation anti-windup scheme [43] is implemented to ensure the integrator outputs v cos (t) and v sin (t) track the scaled signals u cos (t) and u sin (t) to avoid integrator windup. The anti-windup gain K aw is non-dimensional and determines the bandwidth of the anti-windup loops. A unit value for K aw implies a bandwidth comparable to the SLC. Typically, this gain would need to be much larger than one. After trial and error, we take K aw = 100. The parameters of the final SLC controller are shown in Table 5. Robust Controller I (Identity Matrix)

Loads and Performance Evaluation
The performance of the control system shown in Figure 15, with the parameters given in Table 5, is analyzed by comparing the evaluation criteria in Table 2 with and without SLC. Damage equivalent load (DEL) [38] provides a measure of the fatigue due to various dynamic loads on a wind turbine. The DELs of the moments are analyzed based on the IEC standard design load case 1.2 [28]. Simulations are run for the full operational range of the turbine. The software TURBSIM [44] is used to generate eight different turbulent wind files for each mean wind speed. The DELs for various loads are calculated using MLife [39]. Turbine performance is also evaluated to determine what impact, if any, the SLC has on the power produced, the rotor angular speed and the collective pitch activity.

Primary Load
The primary objective of the SLC is to reduce the flapwise root bending moments on the turbine blades. The SLC impact on these loads is evaluated with steady and turbulent wind. Figure 19 shows the effect of the load control system on the blade flapwise root bending moment for 18 m/s steady wind with vertical shear exponent α = 0.2. The top plot shows the time series of the actuator command where the SLC is switched on at 400 s. The middle plot shows the corresponding time series of the bending moment. A 37% reduction in the peak-to-peak value of the moment signal is observed. The frequency domain plot at the bottom of Figure 19 shows a 43% reduction of the 1P harmonic component. There are very small changes in the 2P, 3P and other higher components.
The change in the bending moment signals is also observed for the case of turbulent wind and vertical wind shear. This is shown in Figure 20. A 32% reduction is observed in the 1P harmonic component along with reductions in a frequency range within 1P ± ω 0 = 1.22 ± 0.37 rad/s. This is expected because the bandwidth of the SLC is chosen as ω 0 = 0.37 rad/s. Figures 19 and 20 also show the theoretical lower bounds for load reduction. These bounds are obtained from the following result, which is proven in Appendix B. Lemma 1. Let G(0) denote the DC-gain of the MBC augmented plant defined in Figure 6 at wind speed V. Assume that the singular values of G(0) are σ 1 ≥ σ 2 . Assume also that the feedback system in Figure 15 Figure 15. Then, the 1P closed-loop response M cl 1P is bounded below as follows: where u max = 0.2 is the authority of the SLA.
Thus, if M ol 1P > σ 1 u max , then d 2 > σ 1 u max , which in turn implies that G(0) −1 d 2 > u max . This last inequality shows that M cos = M sin = 0 is not achievable regardless of the controller used. On the other hand, when G(0) −1 d 2 ≤ u max , the controller in Figure 15 does not saturate and due to the integral control component, the closed-loop response vanishes in steady state, achieving full load reduction, i.e., M cl 1P = M 2 cos + M 2 sin = 0.

57%
Theoretical bound from Equation (19) 32% Figure 20. SLC impact on flapwise root bending moment for 18 m/s turbulent wind with vertical shear exponent α = 0.2 and turbulence intensity TI = 16%. Achieved load reduction at 1P is 32% compared with the 57% theoretical upper bound. Note that the achieved load reduction is smaller than the theoretical upper bound from Equation (19) (which assumes steady wind) due to the presence of turbulence.
The DEL for the blade flapwise root bending moment is also calculated for turbulent wind conditions with vertical shear, according to the IEC standard. Figure 21 shows this DEL with and without SLC.
It can be seen that SLC is able to reduce the DEL for the blade flapwise root bending moment, i.e., the primary load targeted by the control system. To quantify the impact of the SLC on the DEL, the following % change used: The % change in DEL for a range of mean wind speeds is shown in the bottom plot of Figure 21. It can be seen that the load control system achieves a reduction up to around 12% depending on the wind speed. Higher reductions in DEL are observed for above-rated wind speeds as compared to below-rated wind speeds.

Secondary Loads
The effect of SLC on the secondary loads on the wind turbine is evaluated as well.
Ideally, the secondary loads should either decrease or not be affected. DEL and % change in DEL are used to measure the effect of SLC on these secondary loads. Figure 22 shows the DEL for the blade root edgewise bending moment with and without SLC. The SLC has a very small effect on the edgewise moments. This is also evident by the fact that % change in DEL is almost zero for all wind speeds. A similar conclusion can be drawn about the tower base moments. Figure 23 shows the DEL and % change in DEL for the tower base fore-aft moment, while Figure 24 shows the same criteria for the tower base side-side moment. Note that the effect of the SLC on these loads is not significant.
Next, we compare the blade tip deflection with and without SLC. Figure 25 shows the variation of the mean, standard deviation and the 95th percentile of blade tip deflection with wind speed. These metrics are calculated using the data from the eight random simulations; thus, the 95th percentile is a more reliable indicator of large deflections than the actual maximum deflection. As expected, tip defection grows up to the rated wind speed (9.8 m/s) and then decays because the turbine "pitches-to-feather" to regulate rotor speed. As expected, the SLC does not affect the mean value of the tip displacement. However, there is reduction in the spread (standard deviation) and the 95th percentile of the displacement values under the SLC. This implies that reducing dynamic fluctuating loads leads to a reduction in dispersion and peak values of tip deflection.

Turbine Activity
First, we demonstrate that the SLC does not have any impact on turbine performance by evaluating the relevant time series at a particular wind condition. Figure 26 shows the time domain signals from the turbine for a steady wind of 18 m/s with vertical wind shear α = 0.2. The controller is switched on at 400 s. The figure shows that the SLC does not affect the turbine performance. The rotor speed, blade pitch, generator torque and power show similar behavior before and after SLC is activated. This is expected because at any instant, the sum of the lift generated by the actuators on the three blades is zero. Moreover, the sum of the corresponding moments due to these lifts is also zero as the actuators are located at the same spanwise position on the three blades. Thus, the fixed frame quantities, similarly to turbine performance, are not affected by the SLC.  The effect of the SLC on turbine performance is evaluated for turbulent wind as well by analyzing the mean and standard deviation (STD) of the sample distributions resulting from the eight random simulations performed to analyze loads. Figure 27 shows the comparison of the mean and STD of the generator power with and without SLC. This figure also shows the theoretical mean and STD calculated using the turbine's maximum power coefficient (C p,max = 0.47) and the sample statistics of V 3 (mean and STD). The rated wind speed is indicated with a vertical dashed red line at 9.8 m/s. The change in the mean and the STD are plotted at the bottom of the figure. Minor changes (<2%) are observed in the mean and STD of the power with SLC, and it can be concluded that SLC does not affect the turbine power in any significant manner.
Similarly, the mean and STD of the rotor speed are compared, as shown in Figure 28. The SLC results in less than 1% changes in these metrics for the rotor speed. The turbine controller keeps the blade pitch constant for below-rated wind speeds. Thus, to evaluate the impact of the SLC on the pitch system, the mean and the STD of the blade pitch are compared for the above-rated wind speeds, as shown in Figure 29. Again, the SLC shows <1% change in the blade pitch parameters.

Scheduling the SLC with Wind Speed
This section shows results for an SLC controller scheduled with wind speed. The scheduled controller is obtained by applying the system identification procedure in Section 4.3 to determine linear models at 12 wind speeds ranging from 3 to 25 m/s. Then, for each plant model, a MIMO integral controller, as shown in Equation (15), is defined using ω 0 = 0.37 rad/s and the DC-gain G(0), which is a function of the wind speed V. Figure 30 shows the DELs for the flapwise blade root bending moments where we compare the controller designed using the model obtained at 18 m/s (red) with the controller scheduled on wind speed (yellow). Note that the scheduled controller results in nearly identical DELs as compared to the single controller designed at 18 m/s. Thus, the added complexity required to schedule the controller on wind speed is not justified.

Conclusions
This paper described a methodology for the design of a dynamic feedback control system to reduce blade fatigue loads. This approach, which we have coined sectional lift control (SLC), uses the local lift coefficient to model on-blade active flow control devices such as controllable Gurney flaps based on DBD plasma actuators. The SLC architecture is based on the well-known multi-blade coordinate transformation (MBC), which allows the reduction of blade loads within a neighborhood of the fundamental rotor angular frequency (1P).
Salient features of our methodology include: 1. The SLC is model-based and can be designed from simple approximate models (physics-based or data) due to the explicit incorporation of a robust stability margin to both parametric uncertainty and unmodeled dynamics.

2.
The SLC does not interact with the main turbine controller and reduces the controlled loads in below-and above-rated wind conditions, without adverse impacts on turbine power or blade pitch activity. 3.
The methodology is directly applicable to an actual wind turbine equipped with DBD plasma actuators (or other flow control devices) by using system identification tools to obtain a model from the actuator signals (voltage for DBD plasma actuators) to the controlled moments.
The design methodology is demonstrated using an OpenFAST model of a 3.4 MW wind turbine equipped with a controllable Gurney flap, controlled via plasma actuation, to reduce the blade-root flapwise bending moments. The study shows that fatigue loads, as measured by DELs of the flapwise bending moments, can be reduced from 4% to 8% in below-rated wind speeds and from 8% to 12% in above-rated wind conditions. Secondary loads, such as blade-root flapwise bending and tower fore-aft bending moments, are not affected by the SLC. The tip deflection in response to stochastic wind inputs is evaluated as well. The statistical analysis of the data shows that the SLC reduces the 95% percentile of the tip deflection and its dispersion about the mean.
It is also shown that a single SLC designed at a specific above-rated wind speed (18 m/s in this case) suffices to control loads from cut-in to cut-out wind speed. This result suggests that there is no need to gain-schedule the SLC on wind speed or another correlated measurement. For above-rated wind speeds, this hypothesis has a controltheoretic explanation. In this case, the main turbine controller is designed to regulate the power and rotor speed to rated values. This feedback loop reduces the sensitivity of wind turbine input/output transfer functions to parameter (e.g., wind speed) variations. Thus, the transfer function from changes in local lift to blade-root moments is relatively independent of wind speed in the bandwidth of interest.
In the companion work [27], the 3.4 MW turbine re-designed with the controllable GF, controlled with the SLC presented in the current paper, was compared to the baseline reference IEA 3.4 MW turbine from Reference [26]. This comparison shows LCOE reduction of 3.5-5.9%, a range based on analysis of different OPEX assumptions. These LCOE savings are derived from capital costs reductions due to wind turbine load reductions in this re-design scenario as the re-designed rotor with Gurney flap (GF) includes minimal AEP increase. As noted in [27], and in part in Section 5, this GF-based active load control system results not only in blade load reductions but also load reductions in other major components (including the pitch system, drive-train, and yaw system)-significant reduction in both capital costs and operational costs result from these load reductions.

SLC
Sectional Lift Control STD Standard Deviation

Appendix A. Details of the Identified Model
The four transfer functions of the transfer matrix G(s) in Equation (6) are shown in Equations (A1)-(A4) (units are kN·m). The corresponding pole-zero plots of the four transfer functions are shown in Figure A1. Real Axis (rad/sec)

Appendix B. Proof of Lemma 1
Under the assumption of an additive disturbance of 1P fundamental frequency on the flapwise blade-root bending moments (the controlled loads) and a steady-state stable closed-loop system, the 1P amplitudes of the controlled moments satisfy which, after completion of squares, yields By definition M ol 1P = d 2 . Moreover, a standard singular value inequality shows that Substituting Equation (A9) into Equation (A8) and using M ol 1P = d 2 , we get the desired conclusion which completes the proof.