Real-Time Performance Analyses and Optimal Gain-Scheduling Control of Offshore Wind Turbine Under Ice Creep Loads

At high latitudes, offshore wind turbines often face unfavorable loads by severe ice-induced vibrations during winter season, which may endanger facilities available on the platforms and degrade operating performance of wind turbine. Hence, it is vital to analyze the influences caused by ice loads. In this paper, based on a real-time simulation model, quantifiable analyses of performance losses due to ice creep loads are firstly studied deeply. Under ice creep loads, reliable pitch control is necessary to ensure safety and high-level power-tracking capability of modern offshore wind turbine. However, uncertain influences of ice creep loads are coupled with wind turbine operation and make it a challenge for wind turbine pitch control using traditional Proportional-Integral (PI) controller from the view of industry. As a result, improved pitch control using optimal gain-scheduling strategy is proposed to alleviate impacts of ice loads where the support vector regression algorithm is adopted to represent the strong nonlinear relationship among PI parameters under different operation conditions. For each operation point, PI parameters are optimally tuned by the particle swarm optimization algorithm. Finally, the presented nonlinear optimal gain-scheduling PI (OGS-PI) controller is applied on regulating generation power and reducing tower top displacement caused by ice creep loads based on software of Fatigue, Aerodynamics, Structures, and Turbulence, a high-fidelity wind turbine simulator. Simulation results show that unfavorable influence of ice creep loads to wind turbine operation can be significantly alleviated by the OGS-PI controller, which performs much better than the traditional PI controller.


I. INTRODUCTION
The 21st century witnesses a new era of rapid development of renewable energy technology. As one of the renewable clean energy, wind energy can bring huge social, economic and environmental benefits. It is estimated that more than 20% of the world's electricity demands will be met by wind energy by 2050 [1], [2]. In recent years, as offshore wind power technology has made great progress, all countries regard offshore wind power as an important development direction of renewable energy [3]- [6]. At the end of 2017, the oceanic waters of eleven European countries had approximately 84% (15,780 MW) of all offshore wind farms, and most of the remaining 16% farms were located in China, The associate editor coordinating the review of this manuscript and approving it for publication was Ehsan Asadi .
Vietnam, Japan, South Korea, and the USA [7], [8]. Chinese government proposed an ambitious plan that the capacity of offshore wind power reaches 30GW by 2020 [9]. At present, the technologies of Chinese offshore wind power are in their infancy, which lacks accurately modeling, exhaustive control systems and multi-scenario optimization. Thus, the studies in this paper are forward-looking and will be valuable for development of wind power in Bohai Sea.
When it comes to the development of offshore wind farm, the challenges cannot be neglected. Previous studies have indicated that offshore wind farms may face more economic, operational, and environmental challenges than onshore wind farms [10]- [15]. Among the challenges, the threats posted by ice are our concern in the present study because offshore structures deployed in cold regions have to undergo ice action [16], which could cause performance degradation and even damage to wind turbine (WT). The Bohai Sea undergoes ice period for about three months every year due to the invasion of strong cold wave. Ice loads can be found in the form of creep failure, crushing failure and buckling failure [17]. The creep failure is selected as the research object of this paper due to the long period of time acting on WT. In the cold regions, ice-induced vibrations are normally occurred when an ice sheet against an ocean structure and severe ice-induced vibrations threaten the structural production facility, which may lead to pipeline fracture and flange loosening [18]. Two different case studies were dynamically assessed in the Bohai Sea, involving jacket platforms collapsed due to ice-induced vibrations in 1969 and 1979 [19].
Several cases and theoretical researches demonstrate that it is significant to study how to reduce the ice-induced vibrations against offshore platform. Yue et al. [20] and Wang et al. [21] proposed that the way of adding a tuned mass damper and installing the cones in compliant connection in the shaft to reduce the vibrations. Wang et al. [22] discussed the way of adding ice-breaking cones at the water level, which will decrease the amplitude of ice loads and change the ice breaking frequency to reduce the vibrations. Liu et al. [23] adopted the finite element method analysis to show that the ice-induced vibration can be significantly reduced for JZ20-2MUQ jacket platform with the isolation cone system. In order to achieve an economical and rational design with the consideration of structural and non-structural performances, Karr et al. [24] used an acceleration-oriented design optimization of ice-resistant jacket platforms in the Bohai Gulf.
While most previous studies focused on the physical methods, which were fragile and difficult to maintain under the harsh condition, this paper aims to design and optimize controller to ensure the safety and economy of WT. Compared with the WT without ice loads, it is more significant to study the response of WT to the ice loads. The tower base moment and the tower top displacements (TTD) are usually selected by scholars to represent the tower response to wind and ice force [25]. The TTD is chosen as one of the evaluation indexes in this work because when it exceeds a certain value, the normal operation of the rotor may be affected. In addition to the ice loads, varying pitch angle can also change TTD [26]. Therefore, it is feasible to reduce the displacement by designing an appropriate pitch controller. However, it is well known that pitch controller is used to obtain a steady output power in variable-speed WT when the wind speed is above the rated value. Hence, this study focuses on finding a balance point between steady output power and the displacement reduction.
This paper designs an optimal gain-scheduling Proportional-Integral (OGS-PI) controller for WT pitch angle control. In order to ensure the output power tracking rated power, PI control is adopted for the pitch controller, which has been proved to be a powerful control tool in industrial processes owing to its robust performance and the simplicity for implementation [27]. Nevertheless, the traditional fixed-gain PI control doesn't work well because the partial derivative of output power to pitch angle also varies with the change of operating condition. Thus, the PI control with gain scheduling of pitch angle is applied in this paper to solve the problem of power fluctuation when the wind speed changes sharply. Additionally, adjusting PI control parameters is the key to stabilize processes. However, traditional methods are difficult to deal with multi-objective optimization. Hence, Particle Swarm Optimization (PSO), which can enhance the adaptability of the controller and was made up for the limitations of traditional parameter tuning is adopted in this paper. Under a certain operating condition, the PI parameter is optimized with the cost function of output power and TTD by applying PSO. Even though proposed optimal control strategy has rapid convergence of finding the best PI parameter under a certain operating condition, it is also a challenge under variable operating conditions. To cope with the problem, nonlinear regression of the optimized parameter of different operating conditions based on support vector regression (SVR) is adopted for real-time gain compensation, which can be regarded as another gain scheduling based on wind speed.
Considering the impact of ice loads, the same control strategy is applied to the corresponding controller. More importantly, combined with detection mechanism of wind speed and ice, the WT can operate safely and economically even if the environment changes rapidly. Additionally, a plenty of simulations are based on the Fatigue, Aerodynamics, Structures, and Turbulence (FAST) which uses high-fidelity models and is developed by Nation Renewable Energy Laboratory (NREL) [28], [29]. The bending moments of tower and blades from the simulation models will not have big error because the most of dynamic of the tower have been considered in FAST. The improving performance of WT due to the proposed controller is demonstrated by the simulation results.
The remainder of the article is organized as follows. Section 2 introduces the theoretical model of WT, environmental loads and structural vibration. Section 3 analyzes the specific impacts on WT performance caused by ice loads quantitatively. Section 4 demonstrates the design of gain scheduling pitch controller, including off-line optimization of PI parameters based on PSO, nonlinear fitting of PI parameters based on SVR and the effectiveness of the proposed controller. Conclusions and future work are discussed in Section 5.

II. MODELING
With the purpose of explaining the effect of ice loads on offshore WT theoretically and further optimal controlling, WT and environmental loads are modelled in this section. As the security of WT is our primary concern, the mechanism of tower top vibration under multiple loads is mainly illustrated.

A. WIND TURBINE MODELING
The aerodynamic system of WT is the key for WT to capture wind energy, which generates torque to realize rotation.  According to Betz theory [30], the aerodynamic power of rotor and mechanical torque captured by WT are as follows: where ρ a is the air density; R is the radius of the wind rotor; v w is the wind speed passing through the rotor; C P (λ, β) is the rotor power coefficient; β is the pitch angle; ω rot is the rotation speed of rotor; λ is the tip speed ratio, which can be calculated as ω rot R/v w .
Variable pitch system is a nonlinear servo system, which is adopted to rotate turbine blades around the axis. In order to simplify the system, the pitch actuator is modelled as a first-order inertial link with finite amplitude and speed limit. The dynamic equation is as follows: where β ref is the reference value of pitch angle; τ b is the time constant of the pitch actuator. In the drivetrain system, this paper completely ignores the flexibility characteristics of high-speed shaft and low-speed shaft. Besides, it does not consider the internal friction and other factors, which equates the transmission system to a single mass model, as shown in Fig. 1.
The dynamic equation for the model is as follows: where J rot and J gen are rotational inertia of turbine rotor and generator rotor; T rot and T gen are mechanical torque and generator electromagnetic torque; ω gen is rotation speeds of rotor and generator; B damp is damping coefficient of low speed shaft; N gear is gear ratio of gearbox. Generator is the core mechanism to realize conversion from mechanical energy to electric energy. Due to the ignoration of the converter model which is not the focus of this paper, the generator is regarded as a torque source, which can be modeled as a first-order inertial link: where τ gen is the time constant of the generator system. Additionally, considering the power losses due to the highspeed shaft driving generator during the conversion process, output power of the generator can be given as follows: P gen = ηT gen ω gen (5) where η is the generator efficiency. Besides, as the modeling of WT in this section has strong consistency with the model of NREL WT, the NREL 5MW offshore WT is adopted as the research object, which is supported by the US Department of Energy and has been applied as a reference WT for a large number of experimental projects. The main parameters of the WT are given in Table 1.

B. ENVIRONMENTAL LOADS MODELING
Offshore WT faces more not only chances but also extreme challenges than onshore WT because of the environmental loads. Hence, the theoretical modeling of loads may explain the dynamic performance of offshore WT under multiple loads. Loads of wind, wave and ice on the structure of offshore WT are the main loads considered in this study.

1) WIND LOADS
In order to calculate the wind loads on offshore WT structure, the blade element momentum theory is adopted to model the loads. Assuming that the blade is composed of blade element, whose shape is obtained from a series of airfoil profiles and has its own aerodynamic characteristics, The loads on blade element are generated by air lift F L and drag F D and the mechanical analysis of it is shown in Fig. 2, where c a is the airfoil chord length; α is the angle of attack; ϕ is the inflow angle; V rel is the relative velocity of blade profile.
The lift and drag of blade element can be calculated as follows: where C L (α) and C D (α) are the air lift and drag coefficient, which are the nonlinear function of the attack angle; r b is the radial length of blade. According to Fig. 2, the resultant force of each blade element in the x direction is as follows: Hence, the total axial wind loads can be obtained by summing the loads of each blade element in x direction: where N b is the number of blades.

2) WAVE LOADS
For an offshore WT, wave forces are the dominant loads which is normally due to wind-generated random waves. The height of the sea at a certain point is a random time series which can be converted to a wave spectrum by Fourier transform. In this paper, Pierson-Moskowitz spectrum, which was mapped under long-term stable sea conditions in the Atlantic Ocean is adopted to describe the irregular waves. The energy spectrum density of Pierson-Moskowitz spectrum can be represented as follows: where H s is the significant wave height and T z is the average period of zero crossing. The spectrum can be transformed into the superposition of sine wave with different frequency and phase angle by IFFT. For each harmonic, the wave kinematics are modeled using Airy wave theory, which can solve the velocity and acceleration of water particle. The velocity of the water particle in the horizontal direction is given by [31].
where A m and f p are the amplitude and frequency of wave. z and h w are the depth and height of the sea. k w is the wavenumber.
The wave loads on a submerged vertical cylinder can be calculated adopting Morison's equation. The wave-induced horizontal force per unit length on the structure is the sum of a nonlinear drag component f d and a linear inertial component f i according to Morison's equation [32], which can be expressed as: where C d and C m are empirical drag and inertia coefficients; D c is the cylinder diameter and ρ w is the water density.
Hence, the wave loads on offshore WT can be obtained by integrating the infinitesimal element in (11) in the vertical direction.

3) ICE LOADS
When ice sheet hits the offshore WT, ice loads can be found in the form of creep failure, crushing failure and buckling failure. The influence of creep failure on WT is studied in this paper, because large creep deformations can develop over long periods of time due to the assumption that WT structure is stiff enough, which has low indentation velocities and low aspect ratios.
During the ice creep, the ice is in full contact with the WT and has a uniform pressure at the ice-structure interface. We consider ice as a viscous material flowing around WT structure with pressure p ice , which is shown in the front view and top view of Fig. 3.
In the creep model, the force of ice sheet on the WT structure has the Korzhavin's empirical expression [33]: where I is the indentation factor, which has the range of 1 to 3. k c is the contact coefficient, with the range of 0.3 (for non-simultaneous failure) to 1 (for simultaneous failure, such as creep) for small scale structures. Sanderson [34] considers that the contact coefficient k c must be very low (0.02-0.13) for full-scale structures. m s is the shape factor, and it is 0.9 for cylindrical structures and 1 for flat indenters. h i is the thickness of the ice sheet and σ is the uniaxial compressive strength of ice. The uniaxial compressive strength of ice depends on the strain rate. Michel and Toussaint [35] proposed the method, which adopted U/4D c as indentation strain rateε. For freshwater granular ice, we use the following equation: where R u = 8.314J·mol −1 K −1 is the universal gas constant; T is the temperature in kelvin; Q g is the activation energy; A g is a constant which depends only on crystal type.
Before the ice stress reaches the ''yield stress'', we assume ice is under elastic strain. The elastic strain when ice begins to ''yield'' can be calculated as: where E y is the Young's Modulus of ice. Then the ice loads increase gradually over time towards a peak value with the assumption of constant strain rate and VOLUME 7, 2019 then stays at a steady-state value. The time from the beginning to the stabilization of ice loads is as follows: Hence, the ice loads on offshore WT can be represented as follows: where t 0 is the start time of ice creep loads. Additionally, a module of FAST named IceDyn [36] has integrated several models of ice load, which can simulate proposed ice model combining offshore WT conveniently. The module of quasi-static ice loading on vertical structure whose sub-model is creep model is selected and the main parameters of the model are given in Table 2.

C. STRUCTURAL VIBRATION MODELING
The structure vibration of offshore WT usually reflects the magnitude of mechanical loads. However, the complexity of structure and the environmental uncertainty make it difficulty  building a general model. Hence, simplification of the WT structure should be the first task. Fig. 4 shows the idealization of vibration problem, which can also be found in [37]. The foundation can be modelled as four springs connected by a rigid base with a lumped mass m b because of the four-legged jacket on shallow foundation. Furthermore, the structure and WT tower can be modelled as a beam with a lumped mass m t at the tip. It is illustrated from Fig. 4 that degrees of freedom in the system are the movement of the rigid base, the bending of the tower and the rotation of the rigid base where x 1 , x 2 are the translations of the springs related to the small angle of rotation θ ; x 3 is the translation due to bending of the tower with the stiffness coefficient k t ; x g represents the translation of the centre of mass, which can be calculated as 0.5(x 1 + x 2 ).
For solving the dynamics problem of non-free particle system, the application of Lagrange's equation, which is represented as follows is more convenient than general kinetic equations.
where q k is the generalized degree of freedom; L is kinetic energy of the particle system; Q k is the generalized force corresponding to q k . The simplified offshore WT structure system is embodied into these parameters as follows: (20) where J g is the moment of inertia of the rigid base; U is potential energy of the system; x t is the displacement of tip of the tower.
After further simplification, (17) can be written in Matrix format, which is shown in (23) where mass and stiffness matrixes are shown in Appendix A. Hence, the vibration of tower top displacement is modeled as a second order forced vibration system, whose vibration frequency depends on the natural frequency and external forces. If external loads are not considered, natural frequencies of the system can be calculated easily. However, the time-varying and nonlinear environmental loads make the system more complex and aggravate the vibration of tower top displacement. Additionally, there is a theoretical basis for designing effective control strategies to alleviate the unfavorable ice loads on tower top displacement by reducing wind loads.

D. PERFORMANCE LOSSES OF WIND TURBINE DUE TO ICE LOADS
In this section, we carry out quantifiable analyses and present the performance losses of WT under different operating conditions due to ice loads. In order to illustrate the performance of generator power, we select Mean Absolute Error (MAE), which can accurately reflect the error and avoid error cancellation and Mean-Square Error (MSE), which can describe the degree of power variation as the quantitative indexes. Considering the significance of the vibration amplitude of TTD, mean and Standard Deviation (SD) of the displacement are selected to evaluate the performance of WT. Additionally, the Maximum Absolute Deviation (MAD) can intuitively reflect the boundary performance. Since the cut-in, cut-out and rated wind speed of NREL 5MW offshore WT are 3 m/s, 25 m/s and 11.4 m/s respectively, we selected the same intensity of turbulence IEC wind with speed of 8 m/s, 20 m/s as the representative of those operating conditions. However, wind around rated speed is not considered because switching of control strategy is not the point of this paper. Fig. 5 shows the wind speed curves of above two operating conditions, which were generated according to the standard of IEC and it can be implemented conveniently by a tool named Turbsim [38].
The sampling interval of the simulation are 0.01 seconds and the simulation time are 60 seconds where only the later 50 seconds are adopted to analyze because the former 10 seconds contain the start-up procedure. When wind speed is  below the rated speed, the control strategy we adopted in this section is optimal torque control [39], which is one of the Maximum Power Point Tracking (MPPT) algorithm to capture the maximum wind energy. Besides, PI controller was applied to adjust generator speed so that output power could track rated power if wind speed is above rated speed.

E. GENERATOR POWER
The generator power curves for the clean and iced case under three operating condition are shown in Fig. 6. The MAE, MSE and MAD of power are used to quantify the economy of WT, as shown in Table 3, where 1 and 0 of Flag ice represent the iced and clean case respectively.
It can be seen from Fig. 6(a) that there is a significant power loss and the MAE of power reduce 51 kW when the wind speed is below rated speed because of the ice loads. Since the pitch controller is running to ensure generator power tracking rated power when the wind speed is above rated speed, as shown in Fig. 6(b), the difference of generator power between clean and iced case is too small to take into account, where the difference of MAE is under 1kW.
Therefore, when the torque controller runs alone, the ice loads will directly reduce the generator power. However, when the pitch controller starts running, the performance losses of WT due to ice loads will not reflect in the generator power because the power losses caused by ice loads is equivalent to the effect of pitch controller. On the other hand, we can still design pitch controller to ensure the performance of generator power if the wind speed is above 11.4 m/s.

F. TOWER TOP DISPLACEMENT
The TTD are selected to represent the tower response to wind, wave and ice loads. The TTD include fore-aft and side-to-side displacement, but we just consider the fore-aft displacement since the side-to-side displacement is too small to take into account. The tower top fore-aft displacement (TTFD) curves for clean and iced cases under three operating conditions are shown in Fig. 7.
The mean, SD and MAD of TTFD are used to quantify the safety of WT, as shown in Table 4. Comparing to the clean WT, it is obvious that the ice loads results in greater TTFD and phase-shafting of time under the three operating conditions. Not only ice loads, but the variation of pitch angle can affect the TTFD. Therefore, the operation condition without pitch controller as seen in Fig. 7(a) can show the direct effect of ice loads on TTFD and there is a notable increase of the mean about 37.2% because of ice loads. It can be seen from Fig. 7(b) and Table 4 that although the mean of TTFD raises 41.2% with the wind speed of 20 m/s, pitch controller can degrade the influence of ice loads, where the SD under ice loads is not greater than the clean case, which improves the safety of WT. However, in order to design an excellent pitch controller, we must find a balance between the generator power and the displacement.

III. DESIGN OF OGS-PI PITCH ANGLE CONTROLLER
In this chapter, a switchable optimal gain-scheduled PI (OGS-PI) pitch angle controller is designed which can realize the automatic identification and the switching of the wind loads and wind-ice loads mode. Additionally, the generator-torque controller which is represented as follows is adopted to coordinate the pitch controller when wind speed  is above the rated speed but it does not the research emphasis in this paper.
The improvement of control strategy and these characteristics are shown in Fig. 8, which consists of control block diagram, optimization, strengths and weaknesses. The basic control strategy is a simple gain scheduling PI control strategy, which only adopts gain scheduling of pitch angle to reduce the power fluctuation. Although the basic strategy is easy implementation, the WT cannot obtain the best output because of the unoptimized parameters. Additionally, the mechanical loads caused by ice creep cannot be relieved applying this controller directly. Hence, PSO algorithm is applied for multi-objective optimization considering power tracking and mechanical loads, which can give the WT better performance and relieve the mechanical loads caused by ice creep. However, the improved controller only works well for the optimized scenario as the optimized parameters were generated in a certain scenario. In addition, it takes too much time for the optimization, which cannot guarantee real-time performance. In order to address these problems, the second improved controller OGS-PI is demonstrated in this section.
The block diagram of OGS-PI controller is shown in Fig. 9, where flag ice is the detection signal of ice loads which determines the selection of the gain functions. The PI functions can be described as: where δ ice is the selection function based on the signal of flag ice , which can be expressed as: The GK(β) is the gain scheduling 1, which overcomes the problem of a negative damping in the speed response because the generator torque drops with increasing speed error. The mathematical justification of the gain compensation can be found in [28]. The compensation function is written as: The output of OGS-PI controller can be expressed as: After the inverse Laplace transform and discretization: x where j represents the j-th interval of simulation and t j represents the corresponding real simulation time.
As shown in Fig. 9, the basis of the controller is the PI controller where Particle Swarm Optimization algorithm is applied to optimize the icing and ice-free scenes of PI parameters under different working conditions. The optimized parameters can meet the requirements of the objectives proposed in this paper under corresponding working condition. Then, Support Vector Regression (SVR) algorithm is adopted to non-linearly fit optimized wind-ice loads, which contribute to gain curves. After the offline optimization, the controller can adjust the PI parameters online according VOLUME 7, 2019 to the gain value under current operating condition mapped to the gain curve, which is recognized as gain scheduling 2.

A. OPTIMIZATION OF PI PARAMETER BASED ON PSO
Particle Swarm Optimization is a swarm intelligence optimization algorithm [40], which has been successfully applied in the fields of multivariate function optimization, neural network training and PID parameter tuning [41]. As the multi-objective optimization of the PI pitch angle controller is a PID parameter tuning problem, the algorithm is selected in this work.
In PSO, the solution of the problem corresponds to the position of a particle in the search space. Each particle has its own position, velocity and fitness value determined by the optimization objective function. In the search space of a D-dimensional target, the position of the i-th particle can be expressed as X i = (X i1 , X i2 , . . . , X iD ) T , and the flight speed can be expressed as V i = (V i1 , V i2 , . . . , V iD ) T . The position with the largest fitness value in the process of each particle optimization is called the individual optimal solution, denoted p i = (p i1 , p i2 , . . . , p iD ) T , and the position with the largest fitness value experienced by the whole particle population is called the global optimal solution, denoted p g = (p g1 , p g2 , . . . , p gD ) T . The initial population of particles is usually randomly generated within the allowed range, and then each particle evolves according to certain rules. The state update equation of evolution can be described as: where i represents the i-th particle; d represents the d-th dimension of the particle; k represents the k-th generation; c 1 and c 2 are acceleration factors, which are used to adjust the step size of particles flying towards their own and the global optimal solution direction, respectively. r 1 and r 2 are random numbers in [0, 1]; w is the inertia weight. The flow chart of the PSO employed in this work is presented in Fig. 10 and PSO is utilized to minimize output power fluctuations and reduce TTFD. The objective function and constraints for the PSO are defined as follows: where w 1 and w 2 are the weights of power deviation and the weights of TTFD deviation; P i and X i are deviation of the two parameters. P max and X max are the maximum operating values of the two parameter deviations. Equation (32) and (33) are implicit expressions of the output power  deviation and fore-aft displacement of tower top. Equation (34) and (35) describe the range and variety rate of pitch angle. Nine IEC turbulent wind with an average wind speed of 12m/s to 20m/s are adopted as operating conditions. Then, PSO algorithm is applied to optimize PI parameters in ice and ice-free scenes respectively. In order to illustrate the optimization effect, the optimization results under the ice scene with an average wind speed of 20m/s is selected to show. The basic parameters are set in Table 5. The selection of population size depends on the number of CPU cores where our servers enable parallel computing, which can run 24 simulations simultaneously. As numerical computation of fluid and structural dynamics equations takes much time, the parameters are optimized off-line, which costs about 7800 seconds per operating condition. The optimization process of PSO is shown in Fig. 10.
The black particles are initial population; green indicates the second generation; blue indicates the third generation; pink indicates the fourth generation and red indicates the fifth generation. The optimal value is the best of red population where P is 0.1335 and I is 0.008.
As indicated in Fig. 11, the particle swarm optimized the parameter step by step and finally found the best point, which demonstrated the swarm intelligence of PSO algorithm. In order to illustrate the impact of optimization on WT performance, Fig. 12 depicted the response of WT with PSO-PI controller compared with PI controller, which referred to the baseline-gain scheduled controller in [28]. The generator power of WT is exhibited by Fig. 12(a), which shows either the PSO-PI or PI controller can track the rated power with little deviations that can be ignored. Obviously, the figure in Fig. 12(b) shows that the TTFD of WT with PSO-PI controller has significantly fewer deviations compared with PI controller. The PI and PSO-PI parameters and the corresponding objective values are documented by Table 6, which shows the PSO optimal PI-based WT models have smaller TTFD deviation and fitness values than the PI-based systems. As the data shown, the performance of WT improves approximate 22.4% for using PSO. More specifically, although the MAE of power increases a little bit about 0.3514 kW, the SD of TTFD reduces by 0.078 m and the MAD of TTFD decreases by 0.2332 m.
Furthermore, dynamic response under random wind speed changing from 12m/s to 20m/s, including clean and ice scenarios are optimized by PSO algorithm, which costs about 140,400 seconds. The results are shown in Appendix B, which can be found that WT with ice loads has different  optimization results compared with clean scenario and the minimum improvement is more than 10%.

B. NONLINEAR REGRESSION OF PI PARAMETERS BASED ON SVR
In order to achieve the improvement of WT performance under variable scenarios, the nonlinear regression of optimized parameters based on SVR algorithm is proposed in this section.
As shown in Fig 9, the inputs of SVR are the ice signal and the nine data pair of wind speed, optimized parameters P and I. Additionally, analytic expression of the parameter P and I with respect to the wind speed and ice signal is the output of SVR. As the parameters under the clean and ice scenario have the same regression process, the regression of parameter P under certain scenario of ice is given as an example in this section. The flow chart of the adopted SVR algorithm is shown in Fig. 13.
The basic idea of support vector regression is to map the sample data to the high-dimensional feature space for linear fitting through a nonlinear mapping (v). Given a training sample: the fitting function is expressed as f p (v) = w· (v)+b, which minimizes the expected risk: where l(·) is the loss function. The nonlinear function fitting of support vector regression introduces an ε-insensitive loss function to modify distance, which means it can tolerate f P (v i ) and P i have a maximum deviation of ε. And ε is used to control the number of support vectors and increase the robustness of regression so that the solution of the nonlinear function is expressed as the following constrained optimization problem [42]: where C is the regularization constant to control the fitting accuracy. And the relaxation variable ξ is introduced in consideration of the fitting error exceeding the precision. Then the Lagrange multiplier method is used to solve this optimization problem, and the Lagrange multiplier is introduced as: to construct the following Lagrange function.
In order to obtain the minimum value of the above formula, we make the partial derivative of L to the independent variables zero respectively and obtain the solution. Then the original optimization problem can be transformed into the corresponding dual problem: Through substituting the solution into (38), we consider the feature mapping form (v) to get: Finally, we obtain the fitting function expression by assuming α is the optimal solution: is a kernel function that satisfies the Mercer condition, and the function can realize the nonlinearization of the algorithm without knowing the specific form of the nonlinear transformation.
The types of kernel functions are Linear Kernel Functions, Polynomial Kernel Functions, and Radial Basis Functions (RBF) [43]. The RBF used in this paper has the advantages of  simple form and good smoothness (the existence of any order derivative). It defines as: Considering that the selection of parameters C and ε will affect the fitting precision (the larger C is, the overfitting will occur; The smaller C is, the underfitting will occur), we adopt the grid search method in this paper to improve the precision of regression and find the optimal parameter value to obtain the best fitting effect. Having applying SVR four times under the different scenarios and parameters, the regression results of PI parameters are shown in Fig. 14. More importantly, the application of the regression results in OGS-PI controller is shown in Fig. 15. When the ice signal or wind speed changes, the OGS-PI controller can adjust the parameter real-time by calculate the regression function, which realizes performance compensation and improve the dynamic characteristics of the WT under ice creep scenario. Additionally, in order to verify the stability of the discrete gain-scheduling system, for all specific conditions, the optimal parameters of adjacent conditions are substituted into the operation for testing. All conditions have the simulation results with stability, which proves the robustness and stability of the system.

C. SIMULATION RESULTS
As variable pitch control is the main control method in this paper, the wind speed is selected in range from 11.4 m/s to 25 m/s. The simulation time is set to 70 seconds where the latter 60 seconds are adopted to analyze. In this simulation, PI, PSO-PI and OGS-PI are applied as the control strategy respectively and the results are compared together. The dynamic response under random wind speed with ice loads  are presented in Fig.16 and the detail data of the objective value are given in Table 8.
The comparison performs in the Table 8 shows that the OGS-PI has the cost function value 27.5% less than the PI and 6.6% less than the PSO-PI. And in the part of generator power, the OGS-PI has approximately zero deviation with rated power as shown in Fig. 16(b). The advantage of OGS-PI is particularly reflected between the 35 and 45 seconds.
During that time, Fig. 16(d) exhibits the output pitch angle of the three controllers where the pitch angle of OGS-PI is adjusted faster under the constraint of rate with less adjustment magnitude than the other two controllers as the wind speed changes rapidly. In other words, the OGS-PI can ensure the generator power exactly tracking rated power. More importantly, in the part of TTFD, the OGS-PI reduces SD of TTFD about 45.5% compared with PSO-PI and 77.1% less than the PI as shown in Fig. 16(c). Additionally, the real maximum TTFD deviation of OGS-PI decreases 42.4% compared with PI and 14.9% less than PSO-PI.
As a result, together with OGS-PI controller, it can weaken the performance degradation of WT due to ice loads. The performance improvements of power and TTFD are mainly attributed to the adaptability of OGS-PI controller to operating conditions and environmental scenarios while the PI and PSO controller can only work effectively for one specific condition.

IV. CONCLUSION AND FUTRUE WORK
In present study, based on the proposed WT model, environmental loads model and structural vibration model, the performance losses of WT due to creep ice loads are analyzed and demonstrated adequately under variable operating conditions. It should be noted that the vibration intensification of TTFD is the main influence of ice loads, which also affects the generator power due to the reduced windward area and the change of attack angle. More importantly, an Optimal Gain Scheduled PI (OGS-PI) pitch angle controller, including the application of PSO and SVR has been designed to regulate the WT to capture the rated wind power and reduce TTFD under ice loads when the wind speed exceeds the rated value. The PSO algorithm is applied repeatedly to optimize the parameter P and I of different operating conditions off-line, which can improve the performance of WT at least 10% under random wind with a certain mean speed. The SVR algorithm is used to fit optimized parameters so that the regression results can be adopted to schedule the gain of controller online. The simulation results from the high-fidelity simulator FAST prove that the OGS-PI based pitch angle controller performs better in constant generator power regulation and TTFD minimization, with stable actuator usage comparing with the PI and PSO-PI controller. However, it should be noted that wind in real life is more volatile and random than wind in simulation. Thus, this causes practical challenges for any controller that relies on wind speed measurement. First, the average wind speed is difficult to measure accurately where estimator or lidar often applied to improve the accuracy. Then, the robustness, optimality and stability cannot be guaranteed when applied in an actual scenario, which can be further improved through adopting the controller with predictive function such as the model predictive control (MPC). Additionally, ice loads have three types, including buckling failure, creep of which we focus on the type in this study, and crushing failure. Future work on rest type of the ice loads will be carried out with more accurate models.

APPENDIXES APPENDIX A MASS AND STIFFNESS MATRIXES
M w and K w as shown at the top of this page.

APPENDIX B
See Table 7.