**1. Introduction**

Upsizing capacity of wind turbines to megawatt-class can increase wind energy capture and has grea<sup>t</sup> 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.



Data-driven modeling of input–output characteristics is a common way, including machinelearning [7–10], standard-model-set approximation [11–13] and subspace identification [14–16]. In [7–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–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 ge<sup>t</sup> 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–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 grea<sup>t</sup> 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–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–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*2/*H*∞ control for maximum power tracking and fatigue load alleviation of the wind turbine [2–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 ge<sup>t</sup> 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:


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.

 **Figure 1.** Modern variable-speed variable-pitch (VSVP) wind turbine 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 \mathcal{C}\_P V^3 \tag{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:

$$\mathbb{C}\_{\mathcal{P}} = f(V, \omega\_{\mathfrak{r}\_{\prime}} \beta) = f(\lambda, \beta) \tag{2}$$

where ωr is turbine rotor speed; λ = *R*<sup>ω</sup>r/*V* is tip-speed-ratio. Usually, *C*P can be used as fitted threedimensional surface or look-up table. Torque of turbine rotor is defined by:

$$T\_{\rm r} = \frac{1}{2} \rho \pi R^3 \mathcal{C}\_{\rm T}(\lambda, \beta) V^2 = \frac{1}{2} \rho \pi R^3 \frac{\mathcal{C}\_{\rm P}(\lambda, \beta)}{\lambda} V^2 = \frac{1}{2} \rho \pi R^5 \frac{\mathcal{C}\_{\rm P}(\lambda, \beta)}{\lambda^3} \omega\_{\rm r}^2 = K\_{\rm r} \omega\_{\rm r}^2 \tag{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 \mathbb{C}\_{\mathcal{F}}(\lambda\_\prime \beta) V^2 \tag{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{cases} \begin{aligned} \mathbf{\_{r}}\dot{\boldsymbol{\omega}}\_{\mathbf{r}} &= \mathbf{T\_{r}} - T\_{\mathbf{ls}} \\ \mathbf{\_{ls}} &= A\_{\text{stiff}}^{\text{ls}} (\boldsymbol{\delta\_{\mathbf{r}}} - \boldsymbol{\delta\_{\mathbf{ls}}}) + B\_{\text{damp}}^{\text{ls}} (\boldsymbol{\omega}\_{\mathbf{r}} - \boldsymbol{\omega}\_{\mathbf{ls}}) \\ &\boldsymbol{J\_{\mathbf{g}}}\dot{\boldsymbol{\omega}}\_{\mathbf{g}} = T\_{\mathbf{hs}} - T\_{\mathbf{g}} \end{aligned} \end{cases} \tag{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; *<sup>A</sup>*stifls and *<sup>B</sup>*dampls are equivalent stiffness and damping coefficients of the low-speed shaft; ωls and <sup>ω</sup>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)/d*t* = ωr, d(δg)/d*t* = <sup>ω</sup>g and δhs = <sup>δ</sup>g, ωhs = <sup>ω</sup>g, *<sup>N</sup>*gear = *T*ls/*T*hs = δhs/δls = <sup>ω</sup>hs/<sup>ω</sup>ls. *<sup>N</sup>*gear is the gearbox ratio. Then, the equivalent of Equation (5) is:

$$\begin{cases} \begin{aligned} T\_{\rm r} &= f\_{\rm r} \dot{\omega}\_{\rm r} + T\_{\rm shaf} \\ T\_{\rm shaf} &= A\_{\rm str} \Big( \delta\_{\rm r} - \frac{\delta\_{\rm g}}{N\_{\rm gur}} \Big) + B\_{\rm damp} \Big( \omega\_{\rm r} - \frac{\omega\_{\rm g}}{N\_{\rm gur}} \Big) \\ &- T\_{\rm g} = J\_{\rm g} \dot{\omega}\_{\rm g} - \frac{T\_{\rm shaf}}{N\_{\rm gur}} \end{aligned} \tag{6} \end{cases} \tag{6}$$

where *T*shaf = *T*ls, *A*stif = *<sup>A</sup>*stifls and *<sup>B</sup>*damp = *<sup>B</sup>*dampls; *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 *<sup>B</sup>*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.

**Figure 2.** Two-mass model.

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\_t \ddot{d} + D\_t \dot{d} + K\_t d = F\_t \tag{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.

**Figure 3.** VSVP control loops of modern wind turbine.

#### *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. *<sup>P</sup>*(*z*<sup>−</sup>1) is transfer function in forward channel; *<sup>C</sup>*(*z*<sup>−</sup>1) 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]:


4. If the controller switches among several regulation laws, closed-loop is structurally identifiable. For a multi-variable control system, it requires *l* ≥ 1 + *<sup>r</sup>*/*<sup>m</sup>*, where *l* is number of feedback controllers; *r* and *m* are input and output dimensions of closed-loop.

**Figure 4.** Typical closed-loop control system.

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.

**Figure 5.** Closed control loops below and above rated wind speed.

#### *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 coe fficient, linear regression and path analysis, are unsuitable while MI (mutual information) [29] becomes an e fficient 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(\mathbf{X}) = -\int p(\mathbf{x}) \log p(\mathbf{x}) d\mathbf{x} \tag{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(\mathbf{X}, \mathbf{Y}) = \iint p\_{\mathbf{X}\mathbf{Y}}(\mathbf{x}, \mathbf{y}) \log \frac{p\_{\mathbf{X}\mathbf{Y}}(\mathbf{x}, \mathbf{y})}{p\_{\mathbf{X}}(\mathbf{x}) p\_{\mathbf{Y}}(\mathbf{y})} \mathrm{d}\mathbf{x} d\mathbf{y} \tag{9}$$

where *pXY*(*<sup>x</sup>*, *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, *pXY*(*<sup>x</sup>*, *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\_{\rm min} = \sum\_{j=1}^{N\_{\rm Out}} \sum\_{i=1}^{N} \alpha\_j [y\_j(i) - \overline{y}\_j(i)]^2 \tag{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; %*y* is estimated output. To evaluate the fitting degree between estimated output and measured output, the fitting percent is used

$$\mathbf{r}\_{\text{FitPercent}} = \frac{\|\overline{y} - y\|}{\|y - \overline{y}\|} \times 100\% \tag{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 grea<sup>t</sup> 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 . ωr ≈ (<sup>ω</sup>r(*k* + 1) − <sup>ω</sup>r(*k*))/*T*, . <sup>ω</sup>g ≈ <sup>ω</sup>g(*<sup>k</sup>* + 1) − <sup>ω</sup>g(*k*) /*T* (*T* is sampling period). Taking *J*r, *J*g, *A*stif and *<sup>B</sup>*damp as identified parameters, the discrete model is:

$$
\begin{bmatrix} T\_{\mathbf{f}}(k) \\ T\_{\mathbf{f}}(k) \end{bmatrix} = \begin{bmatrix} \frac{\omega\_{\mathbf{f}}(k+1) - \omega\_{\mathbf{f}}(k)}{T} & 0 & \delta\_{\mathbf{f}}(k) - \frac{\delta\_{\mathbf{f}}(k)}{N\_{\text{gour}}} & \omega\_{\mathbf{f}}(k) - \frac{\omega\_{\mathbf{f}}(k)}{N\_{\text{gour}}} \\\ 0 & \frac{\omega\_{\mathbf{f}}(k+1) - \omega\_{\mathbf{f}}(k)}{T} & \frac{1}{N\_{\text{gour}}} \left(\delta\_{\mathbf{f}}(k) - \frac{\delta\_{\mathbf{f}}(k)}{N\_{\text{gour}}}\right) & \frac{1}{N\_{\text{gour}}} \left(\omega\_{\mathbf{f}}(k) - \frac{\omega\_{\mathbf{f}}(k)}{N\_{\text{gour}}}\right) \end{bmatrix} \begin{bmatrix} I\_{\mathbf{f}} \\ I\_{\mathbf{f}} \\ A\_{\text{\mathbf{diff}}} \\ B\_{\text{\mathbf{d}}\text{amp}} \end{bmatrix} \tag{12}
$$

Taking *M*t, *D*t and *K*t as identified parameters, discrete model of the tower model (Equation (7)) is

$$F\_{\mathbf{t}}(k) = \begin{bmatrix} \frac{d(k+1) - 2d(k) + d(k-1)}{T^2} \\ \frac{d(k+1) - d(k)}{T} \\ d(k) \end{bmatrix}^T \begin{bmatrix} M\_{\mathbf{t}} \\ D\_{\mathbf{t}} \\ K\_{\mathbf{t}} \end{bmatrix} \tag{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:

