Comparative Analysis of Identification Methods for Mechanical Dynamics of Large-Scale Wind Turbine

Abstract

With increasing size and flexibility of modern grid-connected wind turbines, advanced control algorithms are urgently needed, especially for multi-degree-of-freedom control of blade pitches and sizable rotor. However, complex dynamics of wind turbines are difficult to be modeled in a simplified state-space form for advanced control design considering stability. In this paper, grey-box parameter identification of critical mechanical models is systematically studied without excitation experiment, and applicabilities of different methods are compared from views of control design. Firstly, through mechanism analysis, the Hammerstein structure is adopted for mechanical-side modeling of wind turbines. Under closed-loop control across the whole wind speed range, structural identifiability of the drive-train model is analyzed in qualitation. Then, mutual information calculation among identified variables is used to quantitatively reveal the relationship between identification accuracy and variables’ relevance. Then, the methods such as subspace identification, recursive least square identification and optimal identification are compared for a two-mass model and tower model. At last, through the high-fidelity simulation demo of a 2 MW wind turbine in the GH Bladed software, multivariable datasets are produced for studying. The results show that the Hammerstein structure is effective for simplify the modeling process where closed-loop identification of a two-mass model without excitation experiment is feasible. Meanwhile, it is found that variables’ relevance has obvious influence on identification accuracy where mutual information is a good indicator. Higher mutual information often yields better accuracy. Additionally, three identification methods have diverse performance levels, showing their application potentials for different control design algorithms. In contrast, grey-box optimal parameter identification is the most promising for advanced control design considering stability, although its simplified representation of complex mechanical dynamics needs additional dynamic compensation which will be studied in future.

1. Introduction

Upsizing capacity of wind turbines to megawatt-class can increase wind energy capture and has great potential to reduce the LCOE (levelized cost of energy) per kilowatt hour for grid-connected wind power [1]. Yet, larger wind turbine causes higher fatigue load, greater rotational inertia and more challenging control difficulty [2,3]. Then, advanced control of the modern VSVP (variable-speed variable-pitch) wind turbine becomes very important, to fully utilize the multiple-degree-of-freedom control potentiality of blade pitches or sizable rotor while considering the multi-objectives such as maximizing conversion efficiency and alleviating fatigue load [4]. However, advanced control design usually depends heavily on a state-space model of the physical system [5,6]. Nowadays, it is still widely studied by industry and academia.

For wind turbine modeling, there are mainly three routes in the literature, shown in

Table 1. Summary of different identification methods.

Methods	Model Forms	Applications	References
Data-driven input-output modeling (black-box)	Machine-learning-based modeling, such as neural network and deep learning neural network.	Dynamic modeling, anomaly identification	[7,8,9,10]
	Standard-model-set-based modeling, such as ARX, ARMAX, BJ and OE with LS or PEM criterion	Dynamic modeling, control design	[11,12,13]
	Subspace identification	Dynamic modeling	[14,15,16]
Mechanism-oriented Modeling (white-box)	Complex mechanism model	High-fidelity simulation	[17,18,19]
	Simplified mechanism model	Control verification of theoretic algorithms	[20,21,22]
Combination-based modeling (grey-box)	RLS parameter identification	Dynamic modeling	[23]
	Optimization-based parameter identification	Dynamic modeling	[24,25]

.

Data-driven modeling of input–output characteristics is a common way, including machine- learning [7,8,9,10], standard-model-set approximation [11,12,13] and subspace identification [14,15,16]. In [7,8,9,10], machine learning algorithms, such as artificial neural network (ANN), support vector machine (SVM), random forest and deep neural network were useful for black-box modeling. Yet, trial-and-error often existed to select neural network structure and to tune parameters. If optimization is used, computation burden is non-negligible. In [11,12,13], model structures such as ARX (auto-regressive), ARMAX (auto-regressive moving average), BJ (Box-Jenkins) and OE (Output-Error) were adopted for identification with LS (least square) or PEM (prediction-error method) criterion, where PRBS (pseudo-random binary excitation signal) was often used as input excitation. To get numerical solution, it usually has requirements about forms and amplitudes of excitation signals, open-loop or closed-loop structure and sampling period, etc. In [14,15,16], subspace identification via MOESP (multivariable output error state space) and PBSIDopt (prediction-based subspace identification) were studied for wind turbine modeling. A state-space model could be obtained via subspace identification, but the reconstructed states did not have physical meanings. Besides, this method suffered great influence from excitation signals, control system structure and sampling period, etc. In general, the above black-box identification methods only focused on the external input–output characteristics of system and bypass the internal operation mechanism. Thus, these black-box models had poor interpretability, unsuitable for control design with stability.

From the first-principle, high-fidelity models of wind turbines were obtained undoubtedly for digital design and simulation [17,18,19]. Yet, the model structures were too complex to be used for control design. Relatively, control design models were often simplified in continuous or discrete state-space form, just reflecting leading dynamics of the system [20,21,22]. For an actual wind turbine with high-order dynamics, only leading dynamics represented by simplified models are concerned for control design while how to accurately identify their parameters is a still key problem.

Utilizing a simplified mechanism model of a wind turbine, only the unknown parameters need to be identified where grey-box parameter identification is useful to combine mechanism-oriented modeling and data-driven modeling. In [23], a normal five-order model of DFIG (double-fed induction generator) in d–q axis was adopted and an RLS (recursive least square) algorithm was used to identify model parameters. Through transforming DFIG model into ARX structure, a LS identification problem was built. Then, excitation signals, sampling period and control structure were also concerned. In contrast, optimal parameter identification provided a different way, where an optimization problem is formed to search parameters with optimal objective. In [24], Cava et.al. brought in an evolutionary multi-objective optimization problem to identify parameters of symbolic model under preselected nonlinear structure, where only tower and rotor speed dynamics versus wind speed, pitch angle and aerodynamic torque were studied. In [25], grey-box parameter identification was applied on five-order DFIG model and one-mass drive-train model via particle swarm optimization (PSO). Due to the model structure preselected, only optimal parameters estimation was needed. It was more insensitive to sampling period and more accommodative to closed-loop and process noises. Identification of complexity and difficulty were also greatly reduced. Additionally, value ranges of parameters based on their physical meaning could be set as constraints for optimal search and then to guarantee reasonableness of identified models. It was very important for advanced control design to be based on state-space models considering multi-objectives such as robust H_∞ control or mixed H₂/H_∞ control for maximum power tracking and fatigue load alleviation of the wind turbine [2,3,4]. Furthermore, intuitiveness of the control design model was helpful not only for steady-state performance but also for transient performance.

In summary, grey-box parameter identification is more advantageous in the three routes to get state-space models, so it is adopted in this paper. Different types of methods will be compared to show their diverse performances and features for advanced control design.

In spite of what control algorithms are used, critical equipment such as aerodynamic system, tower system and drive-train system play the dominant roles. With consideration of multi- objectives for control design such as dispatched-power-point-tracking and fatigue load alleviation, the tower model along fore–aft direction and the two-mass model of drive-train are generally used due to their appropriate representations to mechanical dynamics. However, identification research of them has not been conducted, especially for parameter identification of the two-mass model under closed-loop structure without excitation experiment. In this paper, it will be carefully discussed. As a result, the main contributions of the paper are as follows:

• Structural identifiability analysis of the two-mass model under general closed-loop control conditions across the whole wind speed range is presented and is useful to judge feasibility of parameter identification in the closed-loop.
• Influence of identified data to identification performance is discussed based on mutual information analysis of identified variables and is helpful to select great identified datasets.
• Grey-box identifications for mechanical dynamics of a wind turbine, including the two-mass model of drive-train and the two-order damping model of tower-top, are studied where subspace identification, RLS and optimal identification are compared under wind scenarios with different turbulence intensities.
• Identification performances and features of different methods are analyzed from views of control design, providing guiding opinions for further improvement.

The rest of this paper is organized as follows. Section 2 introduces basic knowledge of the VSVP wind turbine including rationality analysis of the Hammerstein structure and selected models. Section 3 analyzes structural identifiability of models in closed-loop and nonlinear correlations among identified variables. The executable identification procedure is proposed in Section 4. Simulation and comparative analysis are shown in Section 5. Section 6 concludes the paper.

2. Basic Knowledge of The Modern VSVP Wind Turbine

In this section, rationality of the Hammerstein structure is analyzed to simplify the identification process. Meanwhile, the selected mechanism models are introduced for parameter identification.

2.1. Rationality of Hammerstein Structure

In [14], wind turbine dynamics were approximated by the connection of a static nonlinear aerodynamic mapping and a linear time invariant mechanical subsystem—the so-called Hammerstein structure. In this paper, this structure is also adopted, shown in Figure 1. From Figure 1, the modern VSVP wind turbine usually has mechanical-side, including an aerodynamic system, drive-train system and tower system, and electrical-side, including an electrical system. The electrical-side has faster response than the mechanical-side. In this paper, only the mechanical-side identification is studied under the Hammerstein structure.

During wind energy capture, the aerodynamic system contributes the main nonlinearity. Response to varying wind speeds, the three flexible blades and their vibration modes yield both low and high frequency dynamics. However, due to the low-pass filtering effects of the wind turbine rotor, the variables such as rotor thrust, torque and speed often manifest themselves. Their spectrum characteristics mainly concentrate in the low-frequency part. Then, the aerodynamic system can be represented by its static nonlinearity where aerodynamic thrust and torque are used for modeling.

For drive-train with gearbox, the two-mass model is capable enough to represent the required leading dynamics for control design. For the tower system, the motion in side–side direction, caused by drive-train and generator dynamics, is often neglected. Only the motion in fore–aft direction is concerned and can be modeled by a mass-spring-damper. The two-mass model and the mass-spring-damper model are tightly coupled with aerodynamic torque and thrust, respectively. Under the Hammerstein structure, aerodynamic torque and thrust become known nonlinear inputs and then drive-train and tower systems are simplified into linear dynamic models.

The electrical subsystem mainly consists of an asynchronous generator and PWM (pulse width modulation) inverter where the generator produces the main nonlinearity. Taking DFIG for example, although its nonlinearity in the five-order model is greatly weakened under d–q rotating-coordinate system, structural nonlinearity between generator rotor speed and current still exists. However, in a sampling period of generator rotor speed, it can be seen to be fixed, relative to the fast change of the current. Then, the five-order model becomes a linear time invariant one. Considering the rapid response ability of the electrical system, it is often simplified as a first-order inertial process.

In summary, the Hammerstein structure is reasonable for wind turbine modeling and helpful to reduce modeling complexity. Considering the different time-scales, only parameter identification of mechanical-side including drive-train and tower systems will be studied in this paper.

2.2. Selected Subsystem Models

Under the Hammerstein structure, models of the mechanical-side with required leading dynamics for advanced control design are selected and introduced.

2.2.1. Aerodynamic Subsystem

In theory, captured wind power of a wind turbine rotor is defined by:

$$P=\frac{1}{2}\rho \pi R^2{C}_P{V}^3$$

(1)

where ρ is air density; R is rotor radius; V is upstream wind speed; C_P is the dimensionless power coefficient, adjusted by rotor speed ω_r and pitch angle β. The nonlinear relationship among them is:

$$C_{\mathrm{P}}=f\left(V,{\omega }_{\mathrm{r}},\beta \right)=f\left(\lambda ,\beta \right)$$

(2)

where ω_r is turbine rotor speed; λ = Rω_r/V is tip-speed-ratio. Usually, C_P can be used as fitted three- dimensional surface or look-up table. Torque of turbine rotor is defined by:

$$T_{\mathrm{r}}=\frac{1}{2}\rho \pi R^3{C}_{\mathrm{T}}\left(\lambda ,\beta \right)V^2=\frac{1}{2}\rho \pi R^3\frac{C_P\left(\lambda ,\beta \right)}{\lambda }{V}^2=\frac{1}{2}\rho \pi R^5\frac{C_P\left(\lambda ,\beta \right)}{{\lambda }^3}{\omega }_{\mathrm{r}}^2=K_{\mathrm{r}}{\omega }_{\mathrm{r}}^2$$

(3)

where C_T (λ,β) is torque coefficient; K_r is torque gain. When optimal λ and C_P are adopted, the optimal K_r can be generated. Aerodynamic thrust to tower is defined by:

$$F_t=\frac{1}{2}\rho \pi R^2{C}_{\mathrm{F}}\left(\lambda ,\beta \right)V^2$$

(4)

where C_F is thrust coefficient.

Equations (1)–(4) represent static characteristics, different from that based on BEM (blade element momentum) theory and AEC (aero-elastic code) [26]. Front inflows are often low-pass filtered even facing rapidly changing inflow angles and wake effects. For a larger size megawatt wind turbine, low-pass filtering effect is more significant in yielding relatively steady characteristic.

2.2.2. Drive-Train Subsystem

Drive-train mainly consists of a low-speed shaft, gearbox and high-speed shaft. The torsional stiffness of them can be considered optionally. Total torsional flexibility in the low-speed side is caused by turbine rotor, tower top, rotor hub, low-speed stage of gearbox, yaw bearing roll and low-speed shaft, etc. In the high-speed side, total torsional flexibility involves high-speed stage of gearbox, high-speed shaft and generator rotor, etc. Equivalent inertias of both sides depend on the drive-train structure, which can be measured via the inertial measurement unit. Usually, equivalent inertia of the low-speed side is ten times more than that of the high-speed side while gearbox inertia is less than that of the high-speed side. Thus, gearbox inertia is often merged into the high-speed side. Then, the main inertia sources are simplified into two parts, yielding the two-mass model. Considering torsional flexibility of the low-speed shaft while taking the high-speed shaft to be rigid, the two-mass model is shown in Figure 2. Torque of the turbine rotor is an input from the aerodynamic subsystem. It drives the turbine rotor to rotate. Due to the existence of torsional flexibility of the low-speed shaft, the transmitted torque of the low-speed shaft is different from torque of the turbine rotor, yielding a torsional angle by different rotation speeds and angle displacements of two-ends. Torque of the generator rotor is a reaction input to balance torque of the turbine rotor. Because the high-speed shaft is taken to be rigid, two-ends of the high-speed shaft have the same torque, rotation speed and angle displacement. The mathematical description of the two-mass model is represented by:

$$\{\begin{array}{c}{J}_{\mathrm{r}}{\dot{\omega }}_{\mathrm{r}}=T_{\mathrm{r}}-T_{\mathrm{ls}}\\ T_{\mathrm{ls}}=A_{\mathrm{stif}}^{\mathrm{ls}}\left({\delta }_{\mathrm{r}}-{\delta }_{\mathrm{ls}}\right)+B_{\mathrm{damp}}^{\mathrm{ls}}\left({\omega }_{\mathrm{r}}-{\omega }_{\mathrm{ls}}\right)\\ J_{\mathrm{g}}{\dot{\omega }}_{\mathrm{g}}=T_{\mathrm{hs}}-T_{\mathrm{g}}\end{array}$$

(5)

where J_r and J_g are equivalent inertias of the low-speed and high-speed sides; T_ls and T_hs are mechanical torques of the low-speed and high-speed shafts in the gear-box; T_g is generator reaction torque; A_stif^ls and B_damp^ls are equivalent stiffness and damping coefficients of the low-speed shaft; ω_ls and ω_g are speed of the low-speed shaft and generator rotor; δ_r, δ_ls, δ_g and δ_hs are angle displacements of the turbine rotor, low-speed shaft, generator rotor and high-speed shaft, respectively. Note that d(δ_r)/dt = ω_r, d(δ_g)/dt = ω_g and δ_hs = δ_g, ω_hs = ω_g, N_gear = T_ls/T_hs = δ_hs/δ_ls = ω_hs/ω_ls. N_gear is the gearbox ratio. Then, the equivalent of Equation (5) is:

$$\{\begin{array}{c}{T}_{\mathrm{r}}=J_{\mathrm{r}}{\dot{\omega }}_{\mathrm{r}}+T_{\mathrm{shaf}}\\ T_{\mathrm{shaf}}=A_{\mathrm{stif}}\left({\delta }_{\mathrm{r}}-\frac{{\delta }_{\mathrm{g}}}{N_{\mathrm{gear}}}\right)+B_{\mathrm{damp}}\left({\omega }_{\mathrm{r}}-\frac{{\omega }_{\mathrm{g}}}{N_{\mathrm{gear}}}\right)\\ -T_{\mathrm{g}}=J_{\mathrm{g}}{\dot{\omega }}_{\mathrm{g}}-\frac{T_{\mathrm{shaf}}}{N_{\mathrm{gear}}}\end{array}$$

(6)

where T_shaf = T_ls, A_stif = A_stif^ls and B_damp = B_damp^ls; T_shaf is internal shaft torque; J_r and J_g can be measured by inertial measurement unit. Torsional flexibility, represented by A_stif and B_damp, needs to be identified. Effective wind speed can be usually estimated via a LIDAR (Light Detection and Range) system in nacelle or soft sensing methods. Giyanani et al. [27] studied the estimation of effective wind speed using LIDAR data. Assume that C_P, C_T and C_F are known from design parameters. Then, aerodynamic torque T_r can be calculated.

In addition to the two-mass model, drive-train can also be modeled as one-mass model or three-mass model. If low-speed and high-speed shafts are both deemed to be rigid, the one-mass model can be obtained. If torsional flexibility is considered, the three-mass model can be obtained. Because the mechanical load is paid more and more attention, the one-mass model is too simplified. The two-mass model is sufficient to represent drive-train dynamic while the three-mass model appears to be complex and unnecessary [26].

2.2.3. Tower Subsystem

Tower dynamics mainly include two aspects, bending motions of fore–aft direction and side– side direction. The former is often caused by wind thrust on turbine rotor. Because the turbine rotor, linked to the low-speed shaft of drive-train, is mounted in the nacelle on a flexible tower top, effects of wind on the rotor can be transmitted to the tower top. Additionally, the latter is caused by coupling effects of drive-train and generator dynamics, etc.

Serving for control design, dominant tower dynamics are considered. In this case, only the tower’s first bending mode of fore–aft direction is modeled by an equivalent mass-spring-damper system, where tower torsion deformation, yawing effects and higher bending modes are neglected. As a result, a sufficiently simplified tower model can be obtained as follows:

$$M_{\mathrm{t}}\ddot{d}+D_{\mathrm{t}}\dot{d}+K_{\mathrm{t}}d=F_{\mathrm{t}}$$

(7)

where M_t, D_t and K_t are mass, damping and spring coefficients; d is fore–aft motion displacement.

3. Identifiability Analysis under Closed-Loop Condition

3.1. Control Strategy of Modern VSVP Wind Turbine

Below rated wind speed, OTC (optimal torque control) strategy is adopted by many industrial wind turbines and seen as a variable-speed controller. Generator torque demand is given by Equation (3) using optimal K_r. Measuring rotor speed, generator torque demand can be derived to control wind turbine operating around optimal operation points.

Above rated wind speed, generator rotor speed is limited below the maximum or rated value, which can be taken as set point. To track it, a variable-pitch controller is activated to regulate pitch angle. Due to the nonlinear dynamics of a wind turbine, in practice, a gain scheduling PI (proportional-integral) controller is generally used, changing with quasi-steady wind speeds.

Both the two control strategies and their switching mechanism are shown in Figure 3.

3.2. Structural Identifiability Analysis of Closed-Loop Control System

The modern wind turbine works in closed-loop which cannot be cut off during operation. As a result, closed-loop identification is necessary. There are mainly two ways for this—direct and indirect methods. Direct identification uses input–output data of controlled object even under closed-loop condition. The indirect method identifies the augmented closed-loop system and deduces the model of the controlled object based on the known controller structure and parameters. When the input–output data of the forward channel is measurable and disturbance signal exists in the feedback channel, the direct method should be preferentially considered as it is more convenient than indirect method.

A typical closed loop is shown in Figure 4. P(z⁻¹) is transfer function in forward channel; C(z⁻¹) is that in feedback channel; R(k) is set point signal; E(k) is deviation signal between set point and output; U(k) is controller output with noise; Z(k) is process output with noise; v(k) and w(k) are noise signals. A closed-loop system is identifiable if any one of the following conditions is satisfied [28]:

• If feedback channel is linear and invariant while no disturbance signal exists and set-point is constant, the identifiable condition is that cancellation between zeros and poles of closed-loop transfer function does not happen, caused by model structure of the feedback channel. Meanwhile, n_p ≥ n_b, n_q ≥ n_a−d_t where n_p and n_q are denominator and numerator order of C(z⁻¹); n_a and n_b are those of P(z⁻¹); d_t is time delay between output and input of the forward channel.
• If continuous excitation signal with enough order exists on the feedback channel and is irrelevant with the noise on the forward channel, a closed-loop system is structurally identifiable.
• If controller is time-varying or nonlinear, a closed-loop is structurally identifiable.
• If the controller switches among several regulation laws, closed-loop is structurally identifiable. For a multi-variable control system, it requires l ≥ 1 + r/m, where l is number of feedback controllers; r and m are input and output dimensions of closed-loop.

All the above conditions provide a decision support to determine whether the system can be successfully identified from a view of structure. This is a fundamental step before identification.

For a VSVP wind turbine, closed-loops below and above rated wind speed are shown in Figure 5. In Figure 5, the variable-speed controller is nonlinear and variant with rotor speed. The variable-pitch controller is usually a gain scheduling proportional-integral controller which is linear and variant. No continuous excitation signals exist on feedback channels. Both of them are single-variable control systems. When wind speed varies around the rated value, the two control strategies with their regulation laws also switch. Thus, condition 3 can be fulfilled and condition 4 can be partially fulfilled. As a result, the closed-loop of a modern VSVP wind turbine is structurally identifiable. This suggests that system or parameter identification can be executed for the controlled objects on forward channel under a closed-loop condition using direct or indirect methods.

3.3. Correlation Analysis of Identified Data

Structural identifiability analysis is qualitative, judging whether the closed-loop identification can succeed. This suggests that enough dynamic information may be contained in the data samples of system variables. In execution, appropriate data samples need to be selected according to their correlations from a relative point. Usually, correlations among system variables are nonlinear. Linear correlation analysis methods such as Person coefficient, linear regression and path analysis, are unsuitable while MI (mutual information) [29] becomes an efficient tool. It is defined by probability density of data without requirements to data distribution and both linear and nonlinear correlation can be analyzed. Firstly, define information entropy H(X) of a time series X as:

$$H\left(X\right)=-{\displaystyle \int p\left(x\right)\log p\left(x\right)}\mathrm{d}x$$

(8)

where H(X) is also called Shannon entropy [28]; p(X) is probability density function of X. Then, the numerical form of uncertainty degree of X can be used to represent its information content. Greater entropy value means stronger information uncertainty. On this basis, MI between X and Y is:

$$I\left(X,Y\right)={\displaystyle \iint p_{XY}\left(x,y\right)\log \frac{p_{XY}\left(x,y\right)}{p_X\left(x\right)p_Y\left(y\right)}}\mathrm{d}x\mathrm{d}y$$

(9)

where p_XY(x, y) is joint probability density function. If X and Y are independent, I(X, Y) = 0; if X and Y are highly dependent, I(X, Y) becomes greater.

In order to keep proper correlations among variables, MI values can be calculated for quantitative analysis. In this paper, p_XY(x, y) are calculated by kernel density estimation [30]. As a result, appropriate data can be determined for identification.

4. Execution of Identification

4.1. Data Acquisition and Preprocessing

For parameter identification, the sampling period is mainly concerned during data acquisition. It depends on maximum or cut-off frequency of the identified object. However, before identification, the frequencies are difficult to be determined, which can be estimated according to the power spectral densities of signals. Of course, trial and error can be used. To avoid missing important information, a smaller sampling period should be preferred to try.

For the data acquired from field, inappropriate low or high frequency components may exist, which affects the parameter identification effect. Low-frequency component mainly refers to slowly changing drift or trend characteristics. High-frequency component refers to interference noise. In this case, band-pass filter design is a feasible solution.

Additionally, if identification results are sensitive to the starting point of time, zero initialization of sampled data may be necessary but optional. Subtracting the mean value of several initial points by the acquired signal can realize zero initialization.

4.2. Optimization Criterion

The two-mass model of drive-train can be seen as a two-input two-output system. Inputs are aerodynamic torque and generator torque. Outputs are generator rotor speed and internal shaft speed. The tower model in fore–aft direction can be seen as a one-input one-output system. The input is aerodynamic thrust. The output is displacement of fore–aft direction. Then, parameter identification of multiple-input multiple-output system appears. Adopting the weighted loss function, optimization criterion in the LS form is defined by:

$$O_{\min }={\displaystyle \sum _{j=1}^{N_{\mathrm{Out}}}{\displaystyle \sum _{i=1}^N{\alpha }_j{\left[y_j\left(i\right)-{\tilde{y}}_j\left(i\right)\right]}^2}}$$

(10)

where N_Out is number of outputs; N is number of sampled data points; T_s is sampling period; α_j is weighting coefficient of each output; y is measured output; $\tilde{y}$ is estimated output. To evaluate the fitting degree between estimated output and measured output, the fitting percent is used

$${\mathrm{r}}_{\mathrm{FitPercent}}=\frac{\Vert \tilde{y}-y\Vert }{\Vert y-\overline{y}\Vert }\times 100\%$$

(11)

where ӯ represents average value. It is actually the normalized root mean squared error.

4.3. Brief Introduction of Identification Algorithms

Facing the advanced control design of the wind turbine in future, only the identification methods applicable for the state-space model are adopted. Among the black-box identification methods, subspace identification is very representative with great performance and convenient executability under state-space structure. Then, for grey-box identification methods using selected mechanism state-space models, there are mainly two categories: RLS and optimal parameter identification. As a result, the three methods are adopted for comparison in this paper.

4.3.1. Subspace Identification

From a statistical perspective, CVA (canonical variate analysis) method was early proposed by Larimore et al. [31]. Based on a geometric concept, Verhaegen et al. [32] proposed MOESP method. Using ARX estimation, SSARX (space state autoregressive exogenous) was given by Jansson et al. [33]. In the simulation, all three methods can be used as candidate methods. For parameter identification of the two-mass model or tower model, the subspace identification method with less MSE (mean squared error) and better fitting percent in Equation (11) can be selected to compare with the other classes of methods. They are applicable for a multi-input multi-output system using Equation (10) as objective.

4.3.2. Grey-Box Parameter Identification via Recursive Least Squares Algorithm

RLS algorithm is used for online parameter estimation of single-input single-output or multi-input single-output system. For the two-mass model (Equation (6)), approximate discretization is used where ${\dot{\omega }}_{\mathrm{r}}\approx \left({\omega }_{\mathrm{r}}\left(k+1\right)-{\omega }_{\mathrm{r}}\left(k\right)\right)/T$, ${\dot{\omega }}_{\mathrm{g}}\approx \left({\omega }_{\mathrm{g}}\left(k+1\right)-{\omega }_{\mathrm{g}}\left(k\right)\right)/T$ (T is sampling period). Taking J_r, J_g, A_stif and B_damp as identified parameters, the discrete model is:

$$\left[\begin{array}{c}{T}_{\mathrm{r}}\left(k\right)\\ T_{\mathrm{g}}\left(k\right)\end{array}\right]=\left[\begin{array}{cccc}\frac{{\omega }_{\mathrm{r}}\left(k+1\right)-{\omega }_{\mathrm{r}}\left(k\right)}{T}& 0& {\delta }_{\mathrm{r}}\left(k\right)-\frac{{\delta }_{\mathrm{g}}\left(k\right)}{N_{\mathrm{gear}}}& {\omega }_{\mathrm{r}}\left(k\right)-\frac{{\omega }_{\mathrm{g}}\left(k\right)}{N_{\mathrm{gear}}}\\ 0& \frac{{\omega }_{\mathrm{g}}\left(k+1\right)-{\omega }_{\mathrm{g}}\left(k\right)}{T}& \frac{1}{N_{\mathrm{gear}}}\left({\delta }_{\mathrm{r}}\left(k\right)-\frac{{\delta }_{\mathrm{g}}\left(k\right)}{N_{\mathrm{gear}}}\right)& \frac{1}{N_{\mathrm{gear}}}\left({\omega }_{\mathrm{r}}\left(k\right)-\frac{{\omega }_{\mathrm{g}}\left(k\right)}{N_{\mathrm{gear}}}\right)\end{array}\right]\left[\begin{array}{c}{J}_{\mathrm{r}}\\ J_{\mathrm{g}}\\ A_{\mathrm{stif}}\\ B_{\mathrm{damp}}\end{array}\right]$$

(12)

Taking M_t, D_t and K_t as identified parameters, discrete model of the tower model (Equation (7)) is

$$F_{\mathrm{t}}\left(k\right)={\left[\begin{array}{c}\frac{d\left(k+1\right)-2d\left(k\right)+d\left(k-1\right)}{T^2}\\ \frac{d\left(k+1\right)-d\left(k\right)}{T}\\ d\left(k\right)\end{array}\right]}^T\left[\begin{array}{c}{M}_{\mathrm{t}}\\ D_{\mathrm{t}}\\ K_{\mathrm{t}}\end{array}\right]$$

(13)

Based on the discrete model, the infinite-history recursive estimation algorithms [34] via forgetting factor are used for parameter identification.

4.3.3. Grey-Box Parameter Identification via Optimization Algorithm

Using Equations (6) and (7), discrete grey-box models can be established via zero-order holder. Then, optimization algorithms are adopted to optimize model parameters, using Equation (10) as loss function. The methods such as Gauss-Newton, Levenberg-Marquardt and trust-region- reflective are candidate optimization algorithms. Then, dominant dynamics of drive-train and tower systems can be approximated. Optimal parameter identification procedure mainly includes:

Step 1:

Set sampling period, acquire data samples and execute data preprocessing.

Step 2:

Set weighting coefficients of optimization criterion.

Step 3:

Estimate and set initial domains of identified parameters.

Step 4:

Identify unknown parameters using optimization algorithms.

Step 5:

Test identified dynamics of each input–output channel. If huge deviation happens, adjust certain steps and repeat Step 4. If required performance is fulfilled, the procedure is over.

5. Simulation

Based on the benchmark simulation demo for a 2 MW wind turbine in GH Bladed software, no excitation signals are applied on the controller or the output. Meanwhile, wind speeds with different turbulence intensities are produced to stimulate the benchmark demo where operation data under different wind scenarios can be acquired. The defaulted controllers in the benchmark demo are used as shown in Figure 3 and Figure 5. The controllers can arbitrarily switch with the varying wind speed.

5.1. Parameters Setting

GH Bladed is a high-fidelity software for wind turbine simulation. In this software, a benchmark simulation demo of a 2 MW wind turbine with gear-box and DFIG is adopted. It uses the controllers and control strategies shown in Figure 3 and Figure 5. The main parameters are listed as follows: rated capacity (2 MW), radius (40 m), hub height (61.5 m), gear-box ratio (83.33), cut-in wind speed (4 m/s), rated wind speed (10m/s), cut-out wind speed (25 m/s), generator inertia (60 kg·m²), stiffness coefficient (1.6 × 10⁸ Nm/rad), damping coefficient (2.5 × 10⁵ Nm·s/rad), variable speed controller (OTC), variable pitch controller (gain scheduling PI controller). Concretely, for this benchmark 2 MW wind turbine model, the horizontal-axis three blades are controlled integratively. For each blade, the blade length is 38.75 m where thickness to chord ratio, Reynolds number, pitching moment center and deployment angle are set as 21%, 2 × 10⁶, 25% and 0°. More information of the turbine blades such as blade structure and aerofoil parameters are shown in Figure A1,

Table A1. Design parameters of blades.

Distance from root (m)	0	3.44444	5.74074	9.18519	16.0741	26.4074	35.5926	38.75
Centre of mass (%)	50	38	29	29	29	29	29	29
Mass per unit length (kg/m)	1084.77	277.356	234.212	209.558	172.577	103.546	55.4713	24.6539
Flapwise stiffness (Nm/rad)	7.472 × 109	1.408 × 109	8.341 × 108	5.561 × 108	2.058 × 108	2.954 × 107	2.259 × 106	3127.98
Edgewise stiffness (Nm/rad)	7.472 × 109	2.085 × 109	1.425 × 109	1.286 × 109	5.648 × 108	1.216 × 108	2.433 × 107	8167.51

and

Table A2. Aerofoil parameters of blades.

Angle of incidence (deg)	−180	−141.27	−100.65	−19	0	21.61	115.85	150.54	180
Lift coefficient	−0.088	0.716	0.094	−0.354	0.449	0.806	−0.539	−0.674	−0.088
Drag coefficient	0.036	0.772	1.167	0.191	0.007	0.288	1.094	0.466	0.036
Pitch coefficient	−0.041	0.362	0.313	0.042	−0.079	−0.097	−0.363	−0.301	−0.041

in the ‘Appendix A’ part.

5.2. Scenarios Setting and Simulation

Under the closed-loop condition, four types of wind scenarios with different turbulence intensities are produced to stimulate the benchmark demo of a 2 MW wind turbine along the whole wind speed range. For each wind scenario, operation data are acquired for identification.

Wind speeds excite wind turbine dynamics. To evaluate volatility of wind speed, turbulence intensity is used, calculated as follows:

$$I_{\mathrm{tur}}=\frac{{\sigma }_{\mathrm{std}}}{V_{\mathrm{mean}}}$$

(14)

$${\sigma }_{\mathrm{std}}=\sqrt{\frac{1}{N_{\mathrm{WS}}-1}{\displaystyle \sum _{i=1}^{N_{\mathrm{WS}}}{\left(V_i-V_{\mathrm{mean}}\right)}^2}}$$

(15)

where V_mean is mean wind speed in time window T_WS with sampling period T_sp; N_WS = T_WS/T_SP is number of sampling points; σ_std is standard deviation; I_tur is turbulence intensity; V_i is sampled values. According to IEC 61400-1 [35], three grades of turbulence intensity such as A(0.16), B(0.14) and C(0.12) are given. In this section, using wind-generation function of GH Bladed [19], wind scenarios with different mean values and turbulence intensities are generated, shown in Figure 6. Equations (14) and (15) were used to calculate mean values and turbulence intensities with different time windows, shown in

Table 2. Characteristic parameters of wind scenarios.

Scenario Types	Mean Wind Speed (m/s)/Turbulence Intensity
	0–10 min	10–20 min	0–20 min
Type 1	7.02/0.15	6.30/0.13	6.66/0.15
Type 2	8.62/0.13	7.72/0.14	8.17/0.15
Type 3	12.31/0.13	11.68/0.18	12.00/0.16
Type 4	15.23/0.15	14.76/0.17	15.00/0.16

. Different wind scenarios’ settings have some differences but the whole wind speed range and the whole intensity range of turbulence can both be tested. Finally, the identified datasets with different operation information caused by different wind scenarios and switched controllers can be obtained for further grey-box parameter identification and validation.

Using wind speeds in Figure 6 as inputs of a 2 MW wind turbine model in GH Bladed, operation data can be acquired. Using the benchmark demo of a 2 MW wind turbine in the GH Bladed software, the virtual sensors and their measuring points are shown in Figure 7. For the identification of the two-mass model, the measuring points include ‘Nominal wind speed at hub position’, ‘Rotor speed’, ‘Rotor azimuth angle’, ‘Generator speed’, ‘Generator azimuthal position’, ‘Generator torque’, ‘Aerodynamic torque’ and ‘Low speed shaft torque’. For the identification of the tower model, the measuring points include ‘Nominal wind speed at hub position’, ‘Tower F_x–Tower station height = 60 m’, ‘Nacelle x-deflection’, ‘Nacelle x-velocity’ and ‘Nacelle x-acceleration’.

Then, MI values of two-mass model and tower model were calculated under different wind scenarios, shown in

Table 3. Mutual information (MI) values of two-mass model under different wind scenarios.

Scenarios	MI Values of Two-Mass Model
	Tr–Tg	Tr–ωr	Tr–ωg	Tr–Tshaf	Tg–ωr	Tg–ωg	Tg–Tshaf	ωr–ωg	ωr–Tshaf	ωg–Tshaf
Type 1	0.4856	0.4555	0.4552	0.5442	3.0705	3.1005	2.5336	5.5664	2.2144	2.2122
Type2	0.6511	0.5651	0.5643	0.7068	2.2504	2.2601	2.5607	5.2716	1.9310	1.9292
Type 3	1.3548	0.2447	0.2417	1.4188	0.2759	0.2751	3.2777	3.2968	0.2760	0.2746
Type 4	0.1381	0.0479	0.0467	0.2241	0.0590	0.0601	0.6571	2.2888	0.0656	0.0658
Scenarios	Normalized MI Values of Two-Mass Model
	Tr–Tg	Tr–ωr	Tr–ωg	Tr–Tshaf	Tg–ωr	Tg–ωg	Tg–Tshaf	ωr–ωg	ωr–Tshaf	ωg–Tshaf
Type 1	0.0872	0.0818	0.0818	0.0978	0.5516	0.5570	0.4552	1	0.3978	0.3974
Type2	0.1235	0.1072	0.1070	0.1341	0.4269	0.4287	0.4858	1	0.3663	0.3660
Type 3	0.4109	0.0742	0.0733	0.4304	0.0837	0.0834	0.9942	1	0.0837	0.0833
Type 4	0.0603	0.0209	0.0204	0.0979	0.0258	0.0263	0.2871	1	0.0287	0.0287

and

Table 4. MI values of tower model under different wind scenarios.

Scenarios		Type 1	Type 2	Type 3	Type 4
MI values	Ft–z	1.3200	1.3080	1.1340	1.0938

. In Table 3, divided by the maximum value of each row, normalized MI values of each row can be obtained.

According to MI values in Table 3, ω_r–ω_g has the highest correlation. Correlation of T_g–T_shaf is relatively high, reflecting the similar dynamics of T_g and T_shaf. Correlation of T_g–ω_r is similar to T_g–ω_g. Correlations of ω_r–T_shaf and ω_g–T_shaf are similar. Correlation of T_r–ω_r is similar to that of T_r–ω_g. Obviously, due to different control strategies and controllers being used below or above rated wind speed, MI values are different under different wind scenarios. From Figure 5, the two-mass model locates at the forward channel of a closed loop. Scenarios Type 1 and Type 2 work below rated wind speed where OTC controller acts. Scenarios Type 3 and Type 4 work above rated wind speed where variable-pitch controller acts. Depending on the acquired data, different data sets yield different MI values. Then, through analysis, connection of control loops and information contained in the acquired data can be preliminarily understood according to control strategies in Figure 5.

The tower model (Equation (7)) is an open-loop system. As a result, in Table 4, correlations between F_t and z are very similar under different wind scenarios.

From the view of basic principle, subspace identification is different from grey-box identification. For subspace identification, estimating the Kalman state vector from the input–output data is the first step and using the Kalman state vector to reconstruct the state space model is the second step. Among them, how to estimate the Kalman state vector from the input–output data is a critical step for subspace identification. For grey-box identification including RLS identification and optimal identification, estimating the model parameters from input–output data based on the selected state-space model structure is the first step and using the obtained state-space model to estimate the Kalman state vector is the second step. The above three identification methods have different basic principles.

Firstly, due to the identified two-mass model including no-self-balancing channel, which is unstable, exploration of different identification methods for direct identification of the two-mass model without excitation experiment under closed-loop condition is mainly studied as follows.

For two-mass model (Equation (6)), define x = [ω_r, ω_g, Δδ]^T, u = [T_r, T_g]^T, y = [ω_g, T_shaf]^T. Then, the state-space equation and transfer function matrix can be obtained, shown as following

$$M_{\text{Two-mass}}=\left[\begin{array}{cc}\frac{b_1{s}^2+b_{2}s+b_3}{s\left(s^2+a_{1}s+a_2\right)}& \frac{-\left(b_{4}s+b_5\right)}{s\left(s^2+a_{1}s+a_2\right)}\\ \frac{b_{6}s+b_7}{s^2+a_{1}s+a_2}& \frac{b_{8}s+b_9}{s^2+a_{1}s+a_2}\end{array}\right]$$

(16)

where a₁ = B_damp/J_r + B_damp/[J_g (N_gear)²], a₂ = A_stif/J_r + A_stif/[J_g (N_gear)²]; b₁ = 1/J_r, b₂ = B_damp/[J_rJ_g (N_gear)²], b₃ = A_stif/[J_rJ_g (N_gear)²]; b₄ = N_gear b₂, b₅ = N_gear b₃; b₆ = B_damp b₁, b₇ = A_stif b₁; b₈ = J_r b₄, b₉ = J_r b₅. Obviously, for the input–output channels from T_r to ω_g and T_g to ω_g, transfer functions with no-self-balancing ability are obtained. For the input–output channels from T_r to T_shaf and T_g to T_shaf, transfer functions with self-balancing ability are obtained. In execution, no-self-balancing object is difficult to be identified.

For the two-mass model, the measured data are obtained from operation of the two megawatts wind turbine in GH Bladed using wind speed inputs in Figure 6. Under the Hammerstein structure, inputs are turbine rotor torque, T_r, and generator reaction torque, T_g, for two-mass model. Outputs are generator rotor speed, ω_g, and internal shaft torque, T_shaf. Then, measured data can be used for identification. Utilizing the identified state-space model with T_r and T_g as inputs, the estimated generator rotor speed and internal shaft torque can be acquired. Herein, all the measured data for identification are produced by the wind turbine model in GH Bladed with high-order nonlinear characteristics. As a result, comparisons of the three identification methods under different wind scenarios are shown in

Table 5. Comparison of identification methods for two-mass model.

Scenarios	Methods	nx, Jr, Jg, Astif, Bdamp	MSE	Fit-Percent	Stability
Type 1	Subspace (MOESP)	nx = 10	7.964 × 107	−175.6%; 96.3%	−0.090 ± 0.951i; −0.284 ± 0.934i; −0.010 ± 0.627i; −0.866 ± 0.365i; 0.995; −0.437; Instability.
	RLS	3.544 × 106; 1; 2.388 × 105; 1.010 × 107	4.214 × 107 3.073 × 1011	−2.785 × 105% −130.064%	−1.458 × 103; −2.363 × 10−2; 0; Instability.
	Optimization	1.381 × 105; 31.8; 2.101 × 108; 1.655 × 105	8.932 × 108	−1.37 × 105%; 87.67%	−0.974±49.719i; −0. Stability.
Type 2	Subspace (MOESP)	nx = 10	8.883 × 107	53.76%; 89.39%	0.784 ± 0.521i; 0.259 ± 0.841i; −0.286 ± 0.824i; −0.712 ± 0.592i; −0.570; 0.994; Instability.
	RLS	5.776 × 106; 1; 1.028 × 105; −2.706 × 106	7.466 × 106; 5.547 × 1010	−206.69%; −165.14%	390.057; 0; 0.038; Instability.
	Optimization	2.063 × 105; 19.43; 2.062 × 108; 1.702 × 105	7.774 × 108	−1.69 × 104%; 72.67%	−1.043 ± 50.266i; −0; Stability.
Type 3	Subspace (MOESP)	nx = 10	3.502 × 107	90.11%; 89.53%	−0.878 ± 0.384i; −0.244 ± 0.939i; 0.435 ± 0.796i; 0.112 ± 0.974i; 0.997; 0.680; Instability.
	RLS	6.366 × 106; 1; 6.920 × 104; −4.531 × 106	3.187 × 106; 2.443 × 1010	−1.527 × 104% −176.42%	653.17; 0; 0.0152; Instability.
	Optimization	2.197 × 105; 20.06; 2.153 × 108; 1.901 × 105	4.559 × 108	−1.768 × 104%; 63.28%	−1.115 ± 50.241i; −0; Stability.
Type 4	Subspace (MOESP)	nx = 10	3.733 × 108	−449%; 70.33%	−0.789 ± 0.405i; −0.400 ± 0.791i; 0.084 ± 0.722i; 0.922 ± 0.362i; −0.845; 0.979; Instability.
	RLS	4.686 × 107; 1; 2.767 × 105; 1.620 × 107	4.650 × 107; 3.395 × 1011	−4.955 × 105%; −794.72%	−2.336 × 103; 0; −0.017; Instability.
	Optimization	2.578 × 105; 15.89; 1.954 × 108; 2.087 × 105	2.742 × 109	−2.734 × 105%; 33.22%	−1.351 ± 50.274i; −0; Stability.

and Figure 8, Figure 9, Figure 10 and Figure 11.

Subspace identification—a data-driven black-box identification method—is mainly realized via robust numerical solution with QR decomposition or singular value decomposition. All the subspace identification methods assume that input signal and process noise signal are independent with each other. It suggests that the system does not have a feedback loop and all poles strictly fall in the left-half-plane, and thus subspace identification can have unbiased estimation for open-loop system. Under closed-loop condition, input signal becomes dependent with process noise, so the closed-loop system identification via subspace identification is difficult. If directly using subspace identification for the closed-loop condition, a biased estimate or even error estimate may happen. Then, eliminating the dependence between input signal and process noise becomes a basic step when using subspace identification for a closed-loop system. In this case, if sufficient excitation can be designed and injected into the closed-loop system and higher signal-to-noise ratio can be yielded, better results can be obtained by subspace identification for the closed-loop system. Its identification accuracy highly depends on the design of the test signal. However, in this paper, system identification under closed-loop condition without excitation signal is mainly studied for the wind turbine model in GH Bladed with high-order nonlinear characteristics. In theory, the identified system will remain biased or with error if no-excitation signal is used under closed-loop condition. From Table 5 and Figure 8, Figure 9, Figure 10 and Figure 11, these methods all have limited identification performance. Through comparison, subspace identification via MOESP has the best performance, showing its good adaptability to closed-loops. Even for the mixed input–output channels, where T_g–ω_g and T_r–ω_g represent no-self-balancing channels and T_g–T_shaf and T_r–T_shaf represent self-balancing channels, subspace identification has good fit-percent using weighted loss function. However, reconstructed states of the discrete high-order system have no-physical-meanings and are disadvantage for control design to improve dynamic performance of the system. Besides, the identified system cannot guarantee its autonomous stability. As a result, subspace identification only yields a good black-box result from the view of fit-percent based on its high-order discrete model structure, whereas the identified system is still unbiased both in theory and in practice.

RLS identification—a black-box or grey-box identification method—is a generalization to LS identification and is mainly solved by numerical solution algorithm. Herein, the drive-train model structure has been selected, so a grey-box problem is formed. RLS identification can have unbiased estimation for an open-loop system without feedback loop where the process noise signal and input signal are independent. Yet, the drive-train model is a no-self-balancing and instable model where closed-loop identification is necessary. Under closed-loop condition, the system needs to be fully excited where a suitable excitation signal should be designed and injected into the system to reduce or even eliminate the dependence between process noise signal and input signal. In this case, the RLS algorithm may have a better identification result. Herein, the simplest identification condition is provided for RLS identification where no-excitation-signal is injected into the closed-loop of drive-train and the data yielded by the closed-loop system are directly acquired for identification. In theory, only unbiased or error estimate may be obtained. RLS identification is just applicable for single output system. Based on Equation (12), T_r is used to estimate parameters. Fit-percent and estimation error of the identification results for the two-mass model performs the worst, as shown in Table 5 and Figure 8, Figure 9, Figure 10 and Figure 11. Besides, for the identified parameters, negative values often appear. These parameters does not always keep consistent with their physical meanings and cannot guarantee autonomous stability of system. Using the estimated parameters, the simulated outputs deviate a lot to the measured outputs and error estimate happens in practice. As a result, for the direct identification of the two-mass model using operation data under the closed-loop condition, RLS identification is difficult to have proper results. It reflects that RLS method hugely relies on the acquired data and excitation experiment.

Grey-box identification utilizes the physical structure and input–output data to establish a system model where a priori physical knowledge and model uncertainty caused by process noise are both considered. It has obvious advantages for capturing the natural behavior of the actual system. Usually, the selected model structure needs include the noise terms to describe model uncertainty. Yet, many unknown noise signals are difficult to be estimated and improper estimation may yield bad results. Herein, the two-mass model is selected as the identified model structure whereas its noise terms are not added, so the approaching ability of the identified two-mass model may limited to the high-order nonlinearity. Then, based on the two-mass model structure, optimization search algorithm is used to find the parameters with the best estimation performance. Shown in Table 5 and Figure 8, Figure 9, Figure 10 and Figure 11, grey-box optimization identification has moderate performance. Although it is applicable for multi-output system, it performs poorly for the mixed channels. Because estimation bias always exists and even diverges for the no-self balancing channels, the fit-percent and estimation error have poor performance. To display the performance, single-step estimates of grey-box optimization identification are also shown in Figure 8a, Figure 9a, Figure 10a, and Figure 11a. It can be found that the single-step estimates are convergent and show their trend tracking abilities to the measured data. For the self-balancing channels, only limited estimation performance is obtained due to the simplified model structure without noise terms and weighted identification with the no-self-balancing channels. Even so, grey-box identification has obvious advantages such as better adaptability to identified data for direct identification under closed-loop condition and is a more convenient identification procedure. Besides, it has the most important advantage that the estimated parameters have physical meaning and they can guarantee autonomous stability of the identified system. As a result, even if the grey-box identification may have limited accuracy, it is still very useful for control design. Especially, through choosing better identified data, a better result of grey-box identification can be yielded. For wind scenarios with higher MI values, better identification performance can be obtained. In Type 1, Type 2 and Type 3, T_r–ω_g, T_g–ω_g, T_r–T_shaf and T_g–T_shaf have higher MI values than that of Type 4, so the identification results of Type 1, Type 2 and Type 3 have less MSE and better fit-percent.

In summary, for parameter identification of the two-mass model, subspace identification gets the best fit-percent. It is suitable for the predictive control design while unsuitable for many advanced optimal control design methods, such as H_∞ control and linear quadratic regulator control, etc. RLS identification is proved to be invalid for the two-mass model under the simple identification condition in this paper. Optimization identification has moderate performance. It is suitable for all kinds of control design algorithms. However, dynamics based on the simplified model structure need to be compensated. Besides, a novel method is needed to realize accurate identification for the no-self-balancing channel. Comparison analyses of above identification results for the two-mass model are briefly summarized and shown in

Table 6. Comparison analyses of identification results for two-mass model.

Identification Methods	Identification Condition Herein	Identification Results	Reason Analysis to Identification Results
Subspace (MOESP)	Direct identification under closed-loop condition without excitation signal	Valid: best fit-percent; instability	Black-box high-order model; closed-loop condition without excitation; no-self-balancing channel
RLS		Invalid: worst fit-percent; instability	Grey-box low-order model without process noise term; closed-loop condition without excitation; no-self-balancing channel
Optimization		Valid: moderate fit-percent; stability	Grey-box low-order model without process noise term; no-self-balancing channel

.

Then, open-loop identification of the tower model is studied as follows. It can directly show the identification performance of different methods.

For tower model (Equation (7)), define x = [d, $\dot{d}$]^T, u = F_t, y = d. Then, state-space equation and transfer function of this two-order spring-damping model can be obtained, shown as follows:

$$M_{\mathrm{Tower}}=\frac{1}{s\left(M_{\mathrm{t}}s+D_{\mathrm{t}}\right)+K_{\mathrm{t}}}$$

(17)

For the input–output channel from F_t to d, it has self-balancing ability.

Under open-loop condition, performance of the three identification methods are strongly dependent on that whether the identified data contain fully excited dynamic information. Herein, only a very simple identification condition is set without excitation signal, so the identification result only depends on the acquired operation data.

For tower model, comparisons of three identification methods are shown in

Table 7. Comparison of identification methods for tower model.

Scenarios	Methods	nx, Mt, Dt, Kt	MSE	Fit-Percent	Stability
Type 1	Subspace (CVA)	nx = 7	2.64 × 104	55.53%	0.698 ± 0.633i; −0.573 ± 0.620i; −0.071 ± 0.596i; 0.191; Instability.
	RLS	4.54 × 104; −1.68 × 104; 1.09 × 106.	8.70 × 10−5	74.55%	0.185 ± 4.896i; Instability.
	Optimization	1.16 × 104; 1 × 103; 1.099 × 106.	1.38 × 10−4	67.93%	−0.129 ± 9.728i; Stability.
Type 2	Subspace (CVA)	nx = 3	2.72 × 10−4	55.66%	−0.318 ± 0.589i; 0.535; Instability.
	RLS	2.88 × 104; −8.93 × 103; 1.13 × 106.	1.04 × 10−4	72.61%	0.155 ± 6.25i; Instability.
	Optimization	1.155 × 104; 1 × 103; 1.151 × 106	1.49 × 10−4	67.33%	−0.130 ± 9.98i; Stability.
Type 3	Subspace (CVA)	nx = 7	2.20 × 10−4	48.09%	−0.549 ± 0.575i; −0.907 ± 0.203i; 0.707 ± 0.482i; 0.921; Instability.
	RLS	4.14 × 104; −7.06 × 103; 1.18 × 106	9.24 × 10−5	66.52%	0.0852 ± 5.343i; Instability.
	Optimization	1.239 × 104; 1 × 103; 1.257 × 106	1.63 × 10−4	55.53%	−0.121 ± 10.072i; Stability.
Type 4	Subspace (CVA)	nx = 3	6.72 × 10−4	38.35%	0.623; −0.150 ± 0.347i; Instability.
	RLS	1.00 × 104; 2.12 × 103; 1.07 × 106	1.79 × 10−4	68.35%	−0.106 ± 10.358i; Stability.
	Optimization	1.116 × 104; 1 × 103; 1.127 × 106	3.07 × 10−4	58.48%	−0.134 ± 10.049i; Stability.

and Figure 12, Figure 13, Figure 14 and Figure 15. RLS identification and grey-box optimization identification both performs better than the subspace identification. It suggests the great advantage of priori model structure information to accurately approach the high-order nonlinear dynamics of tower system. Meanwhile, the arbitrarily acquired data may contain insufficient excited dynamic information which affect the estimation of Kalman state vector form the input–output data, so subspace identification performs the worst. It also suggests that subspace identification is a data-driven black-box identification method and its performance strongly depends on the acquired data and excitation signal.

Especially, RLS achieves the best fit-percent, while parameters incongruent with its physical meaning may be identified and they cannot guarantee autonomous stability of the system. Subspace identification via CVA has the worst fit-percent, which just uses input–output data and reconstructs the states with no-physical meanings. It also cannot guarantee autonomous stability of the identified system. In contrast, optimization identification has moderate fit-percent, but it is sensitive to the set of initial domains of parameters. Its advantage is that the identified parameters have clear physical meanings and they can guarantee autonomous stability of the identified system. However, the simplified model structure limits its representational ability to the nonlinear tower dynamics. For the subsequent applications, subspace identification is suitable to test identifying feasibility of acquired data in the early stage. RLS is suitable to determine the rough range of identified parameters in the middle stage. Then, optimization identification is suitable for the final stage to obtain parameters with physical meaning while guaranteeing autonomous stability of the open-loop tower system. Of course, if simplified mechanism model is adopted and the fit-percent does not perform well, reasonable dynamic compensation can be added. Besides, because the two-order spring-damping model of the tower system lies in an open-loop, similar MI values are obtained under different wind scenarios. As a result, identification performance under different wind scenarios are very similar, too. Comparison analyses of the above identification results for the tower model are briefly summarized and shown in

Table 8. Comparison analyses of identification results for tower model.

Identification Methods	Identification Condition Herein	Identification Results	Reason Analysis to Identification Results
Subspace (CVA)	Direct identification under open-loop condition without excitation signal	Valid: worst fit-percent; instability	Black-box high-order model; insufficient excitation of operation data
RLS		Valid: best fit-percent; instability	Grey-box low-order model without process noise term
Optimization		Valid: moderate fit-percent; stability	Grey-box low-order model without process noise term

.

The main purpose of this paper is to explore the identification performance of three representative methods under a simple condition with minimal complexity. At last, application potential and features of these identification methods are summarized and shown in Figure 16. Overall, using the minimal complexity for identification, grey-box optimization identification has the best adaptability to obtain a valid state-space model with physical meaning and guaranteed stability. It is very helpful for subsequent control design to improve both the stable and dynamic performance of the system. The work herein will be part of a hybrid modeling method in future where the biased estimation of the identified state-space model will be compensated by a non-parametric machine learning algorithm. As a result, the hybrid modeling method can efficiently balance modeling complexity and accuracy where the state-space model with high-precision approaching ability to high-order nonlinearity can be obtained.

6. Conclusions

In this paper, identification of wind turbine mechanical dynamics is studied under non-excitation condition. Identification performance under different wind scenarios is tested using three types of methods. For VSVP wind turbine, the drive-train subsystem is structurally identifiable under closed-loop condition and direct identification is feasible for the two-mass model. Through MI calculation, nonlinear correlations among identified variables can not only validate whether the linkage of these variables are consistent with the control loop but also reveal the relationship between identified data and identification performance. It can be found that higher correlations among identified variables can yield better identification performance. In contrast, state-space model from optimal identification can reflect the physical meaning of parameters and natural stability of the identified system which is important for advanced control algorithms. In summary, grey-box optimal identification shows its feasibility to identify complex wind turbine dynamics and its great potential in advanced control design. Additionally, the limitation of the simplified mechanism model to represent complex and practical dynamics should be paid attention. In future, dynamic compensation to the identified simple mechanism model based on machine-learning will be studied. It may balance modeling complexity and difficulty and would be attractive to the application of digital-twin modeling of wind turbines.