A MIMO periodic ARX identification algorithm for the Floquet stability analysis of wind turbines

The paper presents a new stability analysis approach applicable to wind turbines. At first, a reduced order periodic model is identified from response time histories, and then stability is assessed using Floquet theory. The innovation of the proposed approach is in the ability of the algorithm to simultaneously consider multiple response time histories, for example in the form of measurements recorded both on the rotor and in the stand still system. As each different measurement carries a different informational content on the system, the simultaneous use of all available signals improves the quality and robustness of the analysis.


Introduction
The linear time periodic approximation for the stability analysis of wind turbines is now amply understood and adopted. Numerical methods as the implicit Floquet theory are able to lower its computational cost, but they require a linearization of the system. For this reason some authors preferred to resort to a system identification approach, which has also the advantage of being applicable to field data. In Ref. [1], the authors conducted the stability analysis of a wind turbine by using a SISO (single input, single output) Periodic ARX (autoregressive exogenous). This identification method successfully accomplished the task. However, it was also observed that the use of multiple inputs and outputs at once might improve the quality of the results, and at the same time simplify the procedure.
In this work an extension of the SISO-PARX to the MIMO (multiple input, multiple output) case is presented. The resulting identification method has guaranteed stability properties, and allows for the mitigation of mode contamination problems. The accuracy of the method has been tested against a simplified wind turbine model.
The task of performing the stability analysis of periodic systems by means of MIMO identifications has been previously tackled in various papers. In Ref. [2] the authors presented a subspace identification method and used it to identify the periodic modes of a simple helicopter model. Wang [3] applied the Partial Floquet analysis to the multibody model of a helicopter. Reference [4] developed a MIMO adaptive identification algorithm for general time-varying systems. The algorithm presented in the present paper is similar to the one discussed in Ref. [5], but differs in some aspects. First, the initial guess for the minimization problem is provided by the Equation-Error method, while in Ref. [5] a two stage least-squares method is used. Secondly, the convergence of the algorithm is enhanced by two constraints derived from Floquet theory. The paper starts by briefly reviewing Floquet theory in Sect. 2. The identification algorithm is illustrated in Sect. 3, by separating the discussion between the Equation-and the Output-Error. Section 4 illustrates the main characteristics of the proposed approach with reference to an analytical wind turbine model. Finally the conclusions are drawn in section 5.

Review of Floquet theory for discrete time systems
In this section we briefly review Floquet theory, while for an in depth explanation the reader is referred to Refs. [6,7]. A discrete time, strictly proper, Linear Time Periodic (LTP) system is written as Here x(k) ∈ R Ns is the state at time index k, u(k) ∈ R Nu is the input and y(k) ∈ R Ny the output. Matrices A(k), B(k), C(k) are periodic, with period K. The State Transition Matrix (STM) associated to this system obeys the following recursion and allows one to compute the state as x(k) = Φ(k, κ)x(κ). The monodromy matrix is defined as the STM after one period, i.e. Ψ(κ) = Φ(κ+K, κ). The eigenvalues θ j of the monodromy matrix are named characteristic multipliers, and do not depend on κ. The system is asymptotically stable if and only if the characteristic multipliers lie within the unit circle, i.e.
It can be proved that the observed STM can be decomposed over the modes j and harmonics n as where and ı = √ −1. The scalars η jn are called characteristic exponents, and play the same role of the eigenvalues of an LTI system. The complex vectors c jn , with N y elements, are named observed periodic mode shapes. Each mode j ∈ [1, N s ] of an LTP system is defined by K characteristic exponents and their associated periodic mode shapes. To measure the relative importance of each characteristic exponent η jn within its mode j, we can introduce the scalars which generalize the output-specific participation factors first introduced in Ref. [1] to the multiple outputs case. Within each mode, the characteristic exponent with the maximum participation is called principal harmonic, while all others are termed super-harmonics. It is important to recognize that the index n of the principal harmonic depends on C(k). The characteristic exponents are converted in continuous time [1] as being ∆t the time step. The natural frequencies ω jn and damping factors ξ jn are obtained from the continuous time characteristic exponents

The identification algorithm
The identification algorithm presented in this paper belongs to the class of the Prediction Error Methods (PEM), and it represents an extension to the MIMO case of the one presented in Ref. [1]. As its SISO version, it is composed of two parts: the Equation-Error (EE) and the Output-Error (OE). The model identification step and the subsequent Floquet analysis is performed by the OE, because its predictor is closer to the true system and it has a better treatment of measurement noise. However, the OE involves a nonlinear optimization, and we use the solution provided by the EE for its initial guess.
To ease the identification process, measurements were first resampled, in order to have an integer number of samples in each period. Next, they were also scaled with respect to their maximum absolute values, in order to improve the problem conditioning.

Equation-Error
In the EE framework, measures are assumed to be generated by the following PARX process (see [1,8]) being y(k) the vector of measures at time index k, and u(k) the vector of inputs. A i (k) is the matrix of Auto-Regressive coefficients at delay i. Similarly, B j (k) is the matrix of eXogenous coefficients. Noise is modeled by the sequence e(k), which is assumed to be white, Gaussian and with periodic variance. By realizing Eq. (9) in state-space form, and proceeding as explained in Ref. [9], it can be proved that the optimal one step ahead predictor associated to Eq. (9) iŝ beingŷ(k) the predicted output. In order to reduce the number of unknowns, we resort to a parsimonious representation of the coefficient matrices. A i (k) and B j (k) are approximated with truncated Fourier series where ψ(k) is the rotor azimuth angle. The AR parameters are collected in the following block row vector Similarly, the X parameters are collected inB, with elements from . All the parameters are then collected in one block row vector Two column vectors for the AR and X basis are formed as ϕ a (k) = 1 cos(ψ(k)) sin(ψ(k)) · · · cos(N Fa ψ(k)) sin(N Fa ψ(k)) T , (14a) By means of the Kronecker product, denoted with ⊗, one can define the regressor vectors which in turn are used to express prediction (10) as a matrix-vector product By collecting the measures from q = max(N a , N b )+1 to N in matrixȳ and forming the regressor matrix Υȳ all predicted outputs can be computed. The parameters are then obtained by solvinḡ in a least-squares sense.

Output-Error
In the OE framework, measures are assumed to be generated by the following PARX process [1,8] With derivations similar to the ones for the EE case, it can be proved that the optimal one step ahead predictor associated to Eq. (19) iŝ By defining the regressor vectorsα the prediction at the current time index can be efficiently computed aŝ and its covariance matrix is given by The PEM aims at estimating parameters Θ by minimizing a scalar function of R N defined as Parameters Θ are computed by solving a nonlinear minimization problem, where the initial guess Θ 0 is provided by the EE. To improve the conditioning of the problem, the parameters are scaled with respect to Γ i,j = 10 log 10 |Θ 0 i,j | .
It is well known that cost functions based on prediction errors may be affected by local minima. To improve convergence towards the global minimum, while at the same time augmenting the flexibility of the method, a two-pronged approach has been adopted. Firstly, two constraints are added to the cost function, and, secondly, multiple local optimizations are performed by random perturbations of the initial guess.

Stability constraint
The minimization of V N (Θ) does not guarantee that the identified system is asymptotically stable, even if measures having a decaying trend are used. This feature can be enforced by realizing the AR sequence in state-space form, and applying Floquet theory.
By following a procedure similar to the one exposed in Ref. [1], one can show that the OE process (19) is equivalent to the following system where The monodromy matrix associated with A(k) can be formally computed as The constraint simply requires that the characteristic multipliers lie within the unit circle, and it is written as The idea of adding a stability constraint has been already proposed in Ref. [10], although for the LTI case and only constraining the largest eigenvalue. However, since the max function is discontinuous, the approach is not suited for a gradient-based minimization. The present expression of the constraint is instead inspired by the one used in Ref. [11] for the Periodic SISO-ARMAX predictor.

Expected participation constraint
Systems having whirling modes, as wind turbines, are particularly prone to mode contamination problems caused by nearly overlapping superharmonics [12]. Allen et al. determined the contamination by looking at the observed periodic mode shapes. However this problem can also be examined by looking at the much simpler outputspecific participation factors. Simplified wind turbine models, as well as previous experience on high-fidelity models and field data, may give a qualitative idea on the behavior of the outputspecific participation factors. To guide the optimization routine away from un-physical solutions, one may impose lower and upper bounds to selected output-specific participation factors. By denoting with j the mode, n the super-harmonic and y the output channel, the constraint is written as φ y jn min ≤ φ y jn ≤ φ y jn max .

Constrained minimization problem
The parameters are identified by solving the following constrained minimization problem The process performs a full (discrete time) Floquet analysis at each call of the nonlinear constraints function. However, due to the small number of states of the identified system, this does not result in excessive computational costs. The characteristic multipliers are the result of an eigenvalue computation, and consequently their order can not be preserved from one call of the routine to the next, thus challenging the convergence of the solver. To reorder the characteristic multipliers, the Modal Assurance Criterion (MAC) is used. The eigenvectors of the monodromy matrix are first computed for the initial guess, and then updated at the end of each iteration. These eigenvectors, denoted s old i , are considered as references for the calculation of the constraint gradients by forward differences. At each call of the nonlinear constraint routine a new set of eigenvectors s new j is computed, and then used to get the MAC matrix .
The characteristic multipliers are reordered by starting from the most consistent eigenvectors, and then moving towards the least consistent ones. A similar issue arises for the constraints on the participations. The user might want to specify bounds to the participations of only some of the modes, harmonics and output channels. To get the participations for the ith output channel it is sufficient to select the following output matrix where 1 i is a row of zeros with N y elements and the sole element i equal to 1. When dealing with wind turbines, the output channels might be in a rotating or in a nonrotating frame. This often implies that the principal harmonic of a mode is displaced by one position when the mode is observed in a rotating instead of a non-rotating frame. To overcome this difficulty, one may define an output matrix C MBC (k) that provides as principal harmonic the one that would have been obtained by the Multi Blade Coordinate (MBC) LTI approximation. This matrix should have the same pattern shown in Eq. (28), and hence it is defined by appropriately weighting the output channels being w i ∈ [0, 1] the weights. To ensure the order of the modes, it is sufficient to compute the participations also with respect to this accessory C MBC (k) output matrix, and then to sort the modes in ascending order of principal frequency. This way the modes will be sorted as in a standard Campbell diagram.

Application to a simplified wind turbine rotor
The performance of the proposed identification method was tested on the simplified three-bladed wind turbine model described in Ref. [1]. Each blade is modeled by two rigid bodies, joined by an equivalent hinge that allows for their edgewise motion. The hub is a point mass, constrained to move only in the side-side direction. Springs and dampers are included in each hinge, as well as at the hub. The rotor is forced to rotate at a constant angular speed, and it is subject only to gravitational loads. The model is tuned in order to represent the dynamics of a 6 MW wind turbine, and it represents the following modes: tower side-side, collective edgewise, backward and forward whirling edgewise.
To draw the Campbell diagrams in the partial load region II, multiple simulations were started from suitable initial conditions. After some test with different model order and output, it has been resolved to select as output signals the three lag angles ζ i , and the hub side-side displacement y H . The original system has eight states, and therefore N a = 2. A spectral analysis of the A(t) model matrix revealed that N Fa = 3 is an appropriate choice. The identification algorithm has been implemented in MATLAB ® . The execution time is of the order of the minutes on standard hardware, and the number of iterations increases with the system period.
The full (continuous time) Floquet analysis performed in Ref. [1] has shown the rich periodic content of the system, in turn indicating appropriate lower boundaries for the participation factors, shown in tables 1a and 1b at rated angular speed. Table 1: Expected lower boundaries of the output-specific participation factors, at rated angular speed (Coll.: collective; BW: backward; FW: forward). Super-harmonics are counted from the principal one, labeled n = 0. By a trial and error approach, weights for C MBC (k) were chosen as w ζ 1 = 0.3 and w y H = 1, while the others were set to zero. The motivation behind this choice is that the periodic mode shapes of the tower and whirling modes are best observed in the non-rotating frame, while the collective edgewise mode is best observed in the rotating frame. Experience as shown that the constraints to the participations should be added one at a time, and only if mode contamination problems arise. In particular, this constraint should not be used to address problems due to lack of identifiability. Figures 1a and 1b show an excerpt of the identification results at rated angular speed. The excellent superposition of the measured and predicted outputs denotes the quality of the results. Figures 2a to 2d show the Campbell diagrams of the identified system. An examination of the figures reveals that, as previously observed in Ref. [1], frequencies are well estimated, but accuracy degrades for the damping and the participation factors. The proximity of the collective edgewise and whirling modes at low angular speeds diminishes the identifiability of the system,

Conclusions and future outlooks
The proposed MIMO-PARX approach has been shown to correctly identify all the modes of the system with good accuracy, over the entire region II. By adding constraints to the participation factors, it was possible to eliminate small contamination problems among periodic modes. This algorithm significantly reduces user workload with respect to its SISO counterpart, since less identifications are needed to obtain the system Campbell diagram. It has however been noticed that the quality of the results strongly depends on the number and type of measures used, as well as on the selected model order. Future improvements include the addition of the Moving Average term to the model, to better account for atmospheric turbulence. It would be also advisable to implement a different method for choosing the initial guess, for example by using a two stage least-squares method, and compare it to the EE approach used here. Further studies will regard the application to multibody wind turbine systems, to understand the effect of unmodeled dynamics.