# Implicit Floquet analysis of wind turbines using tangent matrices of a non-linear aeroelastic code P. F. Skjoldan1 and M. H. Hansen2 1 Loads, Aerodynamic and Control, Siemens Wind Power A/S, DK-2630 Taastrup, Denmark 2 Wind Energy Division, National Laboratory for Sustainable Energy, Risø DTU, DK-4000 Roskilde, Denmark # ABSTRACT The aeroelastic code BHawC for calculation of the dynamic response of a wind turbine uses a non-linear finite element formulation. Most wind turbine stability tools for calculation of the aeroelastic modes are, however, based on separate linearized models. This paper presents an approach to modal analysis where the linear structural model is extracted directly from BHawC using the tangent system matrices when the turbine is in a steady state. A purely structural modal analysis of the periodic system for an isotropic rotor operating at a stationary steady state was performed by eigenvalue analysis after describing the rotor degrees of freedom in the inertial frame with the Coleman transformation. For general anisotropic systems, implicit Floquet analysis, which is less computationally intensive than classical Floquet analysis, was used to extract the least damped modes. Both methods were applied to a model of a three-bladed $2.3\;\mathrm{MW}$ Siemens wind turbine model. Frequencies matched individually and with a modal identification on time simulations with the non-linear model. The implicit Floquet analysis performed for an anisotropic system in a periodic steady state showed that the response of a single mode contains multiple harmonic components differing in frequency by the rotor speed. Copyright $\copyright$ 2011 John Wiley & Sons, Ltd. # KEYWORDS modal analysis; Floquet analysis; rotor dynamics # Correspondence P. F. Skjoldan, Loads, Aerodynamic and Control, Siemens Wind Power A/S, Dybendalsvænget 3, DK-2630 Taastrup, Denmark. E-mail: peter.skjoldan@siemens.com Received 26 June 2010; Revised 7 October 2010; Accepted 12 February 2011 # 1. INTRODUCTION Today, advanced non-linear finite element codes1–3 are routinely used for load calculations on wind turbines. Most wind turbine stability tools for calculation of the aeroelastic modes are, however, based on separate linearized models. Stability analysis can be divided into three steps: first, a calculation of the steady state; then, a linearization of the equations of motion about the steady state and last, a modal analysis to extract modal frequencies, damping and mode shapes. This paper presents an approach to structural modal analysis applicable to any periodic steady state where the linearization is obtained directly from the non-linear wind turbine aeroelastic code BHawC.3 The equations of motion for a wind turbine operating at a constant mean rotor speed contain periodic coefficients, preventing direct eigenvalue analysis of the system. Most recent wind turbine stability tools $4\mathrm{-}7$ incorporate the Coleman transformation, also known as the multiblade coordinate transformation, which describes the rotor degrees of freedom in the inertial frame. This transformation eliminates the periodic coefficients if the system is isotropic, i.e. the rotor consists of identical symmetrically mounted blades, and the environment conditions are symmetric. Floquet analysis is, however, applicable to anisotropic systems and any periodic steady state. It requires integration of the equations of motion over a period of rotor rotation, as many times as there are state variables in the system. Because of the computational burden of this approach, it has only been applied to reduce or simplify wind turbine models with a limited number of degrees of freedom.8–10 One way to reduce the computation time is to use the Fast Floquet Theory11 where only one third of the integrations are necessary for a three-bladed isotropic rotor. Another way is to use implicit Floquet analysis12 where the least damped modes can be extracted after a limited number of integrations. Stol et al.13 compare the Floquet analysis with the Coleman transformation approach applied to a periodic steady state, where the remaining periodic coefficients are eliminated by averaging and find small differences in modal frequencies and damping, concluding that it is not necessary to use Floquet analysis. Another approach to modal analysis is system identification,14–16 which operates on the response from numerical simulations or measurements, and no knowledge of the system equations is needed to extract the modal properties. The accuracy of the methods is, however, limited and depends on the chosen excitation. In this paper, tangent matrices for mass, damping and stiffness are extracted from the aeroelastic code BHawC. If the system is isotropic and the steady state is stationary, the Coleman transformation is applied before extracting the modal parameters by eigenvalue analysis. For an anisotropic system, implicit Floquet analysis is used for the modal analysis. When the system is isotropic, the response of a single mode contains a single harmonic component for tower degrees of freedom and up to three components for the blades. The response of a single mode in the anisotropic system on both blades and tower contains multiple harmonic components differing in frequency by the rotor speed. Section 2 of this paper describes the BHawC model, and Section 3 explains the methods for modal analysis, the Coleman transformation approach, the implicit Floquet analysis and also the partial Floquet analysis, a system identification technique. In Section 4, the methods are applied to a model of a wind turbine. Section 5 discusses the approaches, and Section 6 concludes the paper. # 2. STRUCTURAL MODEL The BHawC wind turbine aeroelastic code3 is based on a structural finite element model sketched in Figure 1, where the main structural parts, tower, nacelle, shaft, hub and blades, are modelled as two-node 12-degrees of freedom Timoshenko beam elements. The code uses a corotational formulation, where each element has its own coordinate system that rotates with the element. The elastic deformation is described in the element frame, whereas the movement of the element coordinate system accounts for rigid body motion. In this way, a geometrically non-linear model is obtained using linear finite elements. The configuration of the system, defined by nodal positions $\pmb{p}$ and orientations $\pmb q$ , nodal velocities $\dot{\pmb u}$ (of both positions and orientations) and nodal accelerations $\ddot{u}$ , must satisfy the equilibrium equation given in global coordinates as $$ f_{\mathrm{iner}}(\boldsymbol{p},\boldsymbol{q},\dot{\boldsymbol{u}},\ddot{\boldsymbol{u}})+f_{\mathrm{damp}}(\boldsymbol{q},\dot{\boldsymbol{u}})+f_{\mathrm{int}}(\boldsymbol{p},\boldsymbol{q})=f_{\mathrm{ext}} $$ where $f_{\mathrm{iner}},f_{\mathrm{damp}},f_{\mathrm{int}}$ and $\pmb{f}_{\mathrm{ext}}$ are the inertial, damping, internal and external force vectors, respectively, and $\dot{\mathrm{O}}=\mathrm{d}/\mathrm{d}t$ denotes a time derivative. The inertial forces depend on the acceleration of the masses, the damping forces are given by viscous damping, the internal forces are due to elastic forces and the external forces contain the aerodynamic forces.17 To find ![](images/8c6afcabcd9b46b013d0b9eaf1928f2d43001587902377b212d0d820be50f7d8.jpg) Figure 1. Sketch of the BHawC model substructures. this equilibrium configuration, increments of the positions and the orientations $\delta\pmb{u}$ , the velocities $\delta\dot{\pmb{u}}$ and the accelerations $\delta\ddot{\pmb{u}}$ are obtained using Newton–Raphson iteration with the tangent relation obtained from the variation of Equation (1) as $$ \mathbf{M}(q)\delta{\ddot{u}}+\mathbf{C}(q,{\dot{u}})\delta{\dot{u}}+\mathbf{K}(p,q,{\dot{u}},{\ddot{u}})\delta u=r $$ where M, C and $\mathbf{K}$ are the tangent mass, damping/gyroscopic and stiffness matrices, respectively, and $r=f_{\mathrm{ext}}+f_{\mathrm{iner}}+$ $f_{\mathrm{damp}}\!-\!f_{\mathrm{int}}$ is the residual. The stiffness matrix is composed of constitutive, geometric and inertial stiffness. The orientation $\pmb q$ of the nodes is described by quaternions, also known as the Euler parameters,18 a general four-parameter representation equivalent to a triad, which for node number $i$ is updated as $$ \pmb q_{i}:=q u a t(\delta\pmb u_{i,\mathrm{rot}})*\pmb q_{i} $$ where $\delta\mathbf{\boldsymbol{u}}_{i,\mathrm{rot}}$ contains three rotations that are assumed infinitesimal and thus commute and where this rotation pseudovector is transformed by the function termed quat into a quaternion, which is used to update the nodal quaternion $\pmb q_{i}$ employing the special quaternion product denoted by $^*$ , which maintains the unity of the quaternion. The nodal positions $\pmb{p}$ , the nodal velocities $\dot{\pmb u}$ and the accelerations $\ddot{u}$ are updated by regular addition of the positional part of $\delta\pmb{u},\,\delta\dot{\pmb{u}}$ and $\delta\ddot{\pmb{u}}$ , respectively. All components in $\pmb{p}$ , $\pmb q$ and $\delta\pmb{u}$ are absolute and described in a global frame. The present work considers small perturbations in position and orientation ${\bf\delta y}$ , velocity $\dot{\mathbf{y}}$ and acceleration $\ddot{\mathbf{y}}$ to a steady state with constant mean rotor speed $\varOmega$ defined by $(p_{\mathrm{ss}},q_{\mathrm{ss}},\dot{\pmb u}_{\mathrm{ss}},\ddot{\pmb u}_{\mathrm{ss}})$ , the steady state positions, orientations, velocities and accelerations, respectively, all periodic with the rotor period $T=2\pi/\varOmega$ . The linearized equations of motion are obtained from equation (2) at $r\approx\theta$ as $$ {\bf M}({q}_{\mathrm{ss}})\ddot{\boldsymbol{y}}+{\bf C}({q}_{\mathrm{ss}},\dot{\boldsymbol{u}}_{\mathrm{ss}})\dot{\boldsymbol{y}}+{\bf K}({p}_{\mathrm{ss}},{q}_{\mathrm{ss}},\dot{\boldsymbol{u}}_{\mathrm{ss}},\ddot{\boldsymbol{u}}_{\mathrm{ss}}){\boldsymbol{y}}=\boldsymbol{\theta} $$ where the matrices $\mathbf{M}$ , $\mathbf{C}$ and $\mathbf{K}$ are the $T$ -periodic tangent system matrices that are employed in the modal analysis described in the next section. # 3. METHODS In this section, the four methods for modal analysis of structures with rotors are presented. # 3.1. Coleman approach The Coleman transformation requires identical degrees of freedom on each blade, and therefore, the equations of motion (equation (4)) in global coordinates were first transformed into substructure coordinates $y_{\mathrm{T}}$ . The transformation is $$ \begin{array}{r l}&{\boldsymbol{y}=\mathrm{\mathbf{T}}\boldsymbol{y}_{\mathrm{T}}}\\ &{\mathbf{T}=\mathbf{diag}(\mathbf{I}_{N_{s}},\mathbf{T}_{\mathrm{r}},\mathbf{T}_{\mathrm{b1}},\mathbf{T}_{\mathrm{b2}},\mathbf{T}_{\mathrm{b3}})}\end{array} $$ where $\mathbf{T}$ is a block diagonal time-variant matrix composed of the identity matrix $\mathbf{I}_{N_{\mathrm{s}}}$ sized by the number of degrees of freedom of the tower, the nacelle and the drivetrain, $\mathbf{T_{r}}$ transforms the degrees of freedom on the shaft and the hub into a hub centre frame and $\mathrm{T}_{\mathfrak{b}j}$ transforms the degrees of freedom on blade number $j=1,2,3$ into a local frame for blade $j$ . The triads were obtained in the periodic steady state, and thus, $\mathbf{T}$ is $T$ -periodic. The time-variant transformation into inertial frame coordinates $z$ is $$ \begin{array}{r l}&{{\mathbf y}_{\mathrm{T}}={\mathbf B}\,z}\\ &{{\mathbf B}={\textbf d i a g}({\mathbf I}_{N_{\mathrm{s}}},{\mathbf B}_{\mathrm{r}},{\mathbf B}_{\mathrm{b}})}\end{array} $$ where $\mathbf{B}_{\mathrm{r}}$ is a simple rotational transformation of the shaft and the hub and $\mathbf{B}_{\mathrm{b}}$ is the Coleman transformation introducing multiblade coordinates for a three-bladed rotor11,19 as $$ \mathbf{B}_{\mathrm{b}}=\left[\begin{array}{l l l}{\mathbf{I}_{N_{\mathrm{b}}}}&{\mathbf{I}_{N_{\mathrm{b}}}\cos\psi_{1}}&{\mathbf{I}_{N_{\mathrm{b}}}\sin\psi_{1}}\\ {\mathbf{I}_{N_{\mathrm{b}}}}&{\mathbf{I}_{N_{\mathrm{b}}}\cos\psi_{2}}&{\mathbf{I}_{N_{\mathrm{b}}}\sin\psi_{2}}\\ {\mathbf{I}_{N_{\mathrm{b}}}}&{\mathbf{I}_{N_{\mathrm{b}}}\cos\psi_{3}}&{\mathbf{I}_{N_{\mathrm{b}}}\sin\psi_{3}}\end{array}\right] $$ where $\psi_{j}=\varOmega t+2\pi(j-1)/3$ is the mean azimuth angle to blade number $j$ and $N_{\mathrm{b}}$ is the number of degrees of freedom on each blade. The inertial frame coordinate vector $$ \boldsymbol{z}=\{y_{\mathrm{s}}^{\mathrm{T}}\,z_{\mathrm{r}}^{\mathrm{T}}\,a_{0}^{\mathrm{T}}\,a_{1}^{\mathrm{T}}\,b_{1}^{\mathrm{T}}\}^{\mathrm{T}} $$ contains the untransformed coordinates for tower, nacelle and drivetrain ${\mathfrak{y}}_{\mathrm{s}}$ , the coordinates for shaft and hub $z_{\mathrm{r}}$ measured in a non-rotating frame aligned with the hub and the multiblade symmetric coordinates $\pmb{a}_{0}$ , cosine coordinates $\pmb{a}_{1}$ and sine coordinates $\pmb{b}_{1}$ . The details on how multiblade coordinates describe the motion of a wind turbine rotor in the inertial frame are discussed by Hansen.20,21 The Coleman transformed equations were obtained by first inserting equation (5) into equation (4), then converting it to first order form and lastly introducing the inertial frame transformation in equation (6) as ${\displaystyle y_{\mathrm{T}2}=\mathrm{diag}({\bf B},{\bf B})z_{2}}$ where ${y}_{\mathrm{T}2}=\{{y}^{\mathrm{T}}\,{\dot{{y}}}^{\mathrm{T}}\}^{\mathrm{T}}$ and $z_{2}=\dot{\{z^{\mathrm{T}}\,\tilde{z}^{\mathrm{T}}\}}^{\mathrm{T}}$ are the state v ectors in substructure and ine rtial frames, respectively, with $\tilde{z}=\dot{z}+\bar{\omega}z$ and the constant matrix $\bar{\boldsymbol{\omega}}=\mathbf{B}^{-1}\mathbf{B}$ . The result is $$ \begin{array}{r l}&{\dot{z}_{2}=\mathbf{A}_{\mathrm{B}}z_{2}}\\ &{\mathbf{A}_{\mathrm{B}}=\left[\mathbf{-}\mathbf{\bar{\omega}}-\mathbf{\bar{\omega}}\mathbf{\bar{\omega}}_{\mathrm{KB}}\quad\mathbf{-M}_{\mathrm{B}}^{-1}\mathbf{C}_{\mathrm{B}}-\mathbf{\bar{\omega}}\right]}\end{array} $$ where $\mathbf{A}_{\mathrm{B}}$ is the Coleman transformed system matrix and $$ \begin{array}{r l}&{\mathbf{M}_{\mathrm{B}}=\mathbf{B}^{-1}\mathbf{T}^{\mathrm{T}}\mathbf{M}\mathbf{T}\,\mathbf{B}}\\ &{\mathbf{C}_{\mathrm{B}}=\mathbf{B}^{-1}\mathbf{T}^{\mathrm{T}}(\mathbf{C}\,\mathbf{T}+2\,\mathbf{M}\,\dot{\mathbf{T}})\mathbf{B}}\\ &{\mathbf{K}_{\mathrm{B}}=\mathbf{B}^{-1}\mathbf{T}^{\mathrm{T}}(\mathbf{K}\,\mathbf{T}+\mathbf{C}\,\dot{\mathbf{T}}+\mathbf{M}\,\ddot{\mathbf{T}})\mathbf{B}}\end{array} $$ are the Coleman transformed mass, damping/gyroscopic and stiffness matrices, respectively. If the system is isotropic, then $\mathbf{A}_{\mathrm{B}}$ is time-invariant, and a transient solution of equation (9) is $$ z_{2}=\mathrm{e}^{\mathbf{A}_{\mathrm{B}}t}z_{2}(0)=\mathbf{V}\mathrm{e}^{\mathbf{A}t}q(0) $$ where $\mathbf{A}$ is a diagonal matrix containing the eigenvalues of $\mathbf{A}_{\mathrm{B}}$ , $\mathbf{V}$ contains the corresponding eigenvectors as columns and $\pmb q(0)=\mathbf V^{-1}z_{2}(0)$ are the initial conditions in modal c oordinates. It is assumed th at all eigenvect ors are linearly independen t . The blade motion given in the inertial frame in equation (11) can be transformed back into the rotating frame using equation (6) as21 $$ \Gamma,i k=\mathrm{e}^{\sigma_{k}t}\left(A_{0,i k}\cos(\omega_{k}t+\varphi_{0,i k})+A_{\mathrm{BW},i k}\cos\left((\omega_{k}+\Omega)t+\varphi_{j}+\varphi_{\mathrm{BW},i k}\right)+A_{\mathrm{FW},i k}\cos\left((\omega_{k}-\Omega)t-\varphi_{j}+\varphi_{\mathrm{GW},i k}\right)\right), $$ where $\varphi_{j}=2\pi(j-1)/3$ and $\sigma_{k}$ and $\omega_{k}$ pare the modal damping and frequency of mode number $k$ , respectively, given by the eigenvalue $\lambda_{k}=\sigma_{k}+\mathrm{i}\omega_{k}$ with $\mathrm{i}=\sqrt{-1}$ . The amplitudes for degree of freedom number $i$ were determined from the components of the eigenvector $\nu_{k}$ gi ven in multiblade coordinates of equation (8) as $A_{0,i k}=|a_{0,i k}|$ and $$ \begin{array}{r l}&{A_{\mathrm{BW},i k}=\frac{1}{2}\big((\mathrm{Re}\,(a_{1,i k})+\mathrm{Im}\,(b_{1,i k}))^{2}+(\mathrm{Re}\,(b_{1,i k})-\mathrm{Im}\,(a_{1,i k}))^{2}\big)^{1/2}}\\ &{A_{\mathrm{FW},i k}=\frac{1}{2}\big((\mathrm{Re}\,(a_{1,i k})-\mathrm{Im}\,(b_{1,i k}))^{2}+(\mathrm{Re}\,(b_{1,i k})+\mathrm{Im}\,(a_{1,i k}))^{2}\big)^{1/2}}\end{array} $$ where the subscripts 0, BW and FW denote symmetric, backward whirling and forward whirling motion, respectively. # 3.2. Classical Floquet analysis Floquet analysis enables the solution of the periodic equations of motion directly without an explicit transformation. Equation (4) is written in first order form $$ \begin{array}{r l}&{\dot{\boldsymbol{y}}_{2}=\mathbf{A}\boldsymbol{y}_{2}}\\ &{\mathbf{A}=\left[\mathbf{-M}^{-1}\mathbf{K}\quad\mathbf{-M}^{-1}\mathbf{C}\right]}\end{array} $$ where ${\mathbf y}_{2}=\{{\mathbf y}^{{\mathrm T}}\,\dot{{\mathbf y}}^{{\mathrm T}}\}^{{\mathrm T}}$ is the state vector and $\mathbf{A}$ is the $T$ -periodic system matrix. Floquet theory22 states that the solution to equation (15) is of the form $$ \mathbf{\boldsymbol{y}}_{2}=\mathbf{\boldsymbol{U}}\mathbf{\boldsymbol{e}}^{\mathbf{\boldsymbol{\Lambda}}t}\mathbf{\boldsymbol{U}}^{-1}(0)\mathbf{\boldsymbol{y}}_{2}(0) $$ where $\mathbf{U}$ is a $T$ -periodic matrix and $\mathbf{A}$ is a diagonal matrix. One way to construct this solution is to form a fundamental solution to equation (15) as $$ \displaystyle\varphi=\bigl[\varphi_{1}\quad\varphi_{2}\quad.\ .\quad\varphi_{N}\bigr] $$ over one period, $t\ \in\ [0;T]$ , where $N$ is the number of state variables, such that $\dot{\varphi}\;=\;{\bf A}\varphi$ . The monodromy matrix defined as $$ \mathbf{C}=\boldsymbol{\varphi}^{-1}(0)\boldsymbol{\varphi}(T) $$ contains all modal properties, which can be extracted from the eigenvalue decomposition $$ \mathbf{C}=\mathbf{V}\mathbf{J}\mathbf{V}^{-1} $$ where $\mathbf{V}$ contains the column eigenvectors $\nu_{k}$ of $\mathbf{C}$ , which are all assumed to be linearly independent and $\mathbf{J}$ is a diagonal matrix containing the eigenvalues $\rho_{k}$ of $\mathbf{C}$ , called the characteristic multipliers. The characteristic exponents $\lambda_{k}=\sigma_{k}+\mathrm{i}\omega_{k}$ contain the frequency $\omega_{k}$ and damping $\sigma_{k}$ and are related to the characteristic multipliers as $\rho_{k}=\exp(\lambda_{k}T)$ . Because the complex logarithm is not unique, the frequency is not determined uniquely, and the principal frequency $\omega_{\mathrm{p},k}$ and the damping $\sigma_{k}$ are defined from the characteristic multipliers as $$ \begin{array}{c}{\displaystyle\sigma_{k}=\frac{1}{T}\ln(\vert\rho_{k}\vert)}\\ {\displaystyle\omega_{\mathrm{p},k}=\frac{1}{T}\arg(\rho_{k})}\end{array} $$ where $\arg(\rho_{k})\in]-\pi;\pi]$ is implied, resulting in $\omega_{\mathrm{p},k}\in\big[-\Omega/2;\Omega/2\big]$ . Any integer multiple of the rotor speed can be added to the principal frequency to obtain a more physically meaningful frequency23,24 $$ \omega_{k}=\omega_{\mathrm{p},k}+j_{k}\Omega $$ a choice that also affects the periodic modal matrix $\mathbf{U}$ in equation (16). This matrix $\mathbf{U}$ contains the periodic mode shapes uk and is given as24 $$ \pmb{u}_{k}=\varphi\nu_{k}\mathrm{e}^{-\lambda_{k}t} $$ where the real part of $\lambda_{k}$ is given by equation (20) and the imaginary part of $\lambda_{k}$ is defined by equation (21) by selecting $j_{k}$ such that $\pmb{u}_{k}$ is as constant as possible for degrees of freedXom measured in the inertial frame. Introducing the Fourier transform of the periodic mode shape $$ {\pmb u}_{k}=\sum_{j=-\infty}^{\infty}u_{j k}\mathrm{e}^{\mathrm{i}j\Omega t} $$ the transient solution in equation (16) can be written as a sum of harmonic components $$ y_{2}=\sum_{k=1}^{N}\sum_{j=-\infty}^{\infty}\mathcal{U}_{j k}\mathrm{e}^{(\sigma_{k}+\mathrm{i}(\omega_{k}+j\Omega))t}q_{k}(0) $$ where $\begin{array}{r}{\pmb q(0)=\mathbf{U}^{-1}(0)\mathbf{y}_{2}(0).}\end{array}$ . Note that equation (12) is a special case of this expression for $j=-1,0,1$ # 3.3. Implicit Floquet analysis The implicit Floquet method is here described based on the detailed description in Bauchau and Nikishkov,12 which focuses on computation of the characteristic multipliers from the state transition matrix $\Phi(T,0)$ . It can be defined in classical Floquet theory as $$ \boldsymbol{\varphi}(T)=\boldsymbol{\Phi}(T,0)\,\boldsymbol{\varphi}(0) $$ Using equation (18), the relationship between the state transition and monodromy matrices is derived as $$ \Phi(T,0)=\varphi(0){\bf C}\,\varphi^{-1}(0) $$ showing that $\Phi(T,0)$ and C have identical eigenvalues (characteristic multipliers), and their eigenvectors are related as $\pmb{\nu}_{k}=\pmb{\varphi}^{-1}(0)\pmb{w}_{k}$ , where $w_{k}$ represents the eigenvectors of $\Phi(T,0)$ . The key feature of the state transition matrix is that it defines the solution $y_{2}(T)=\Phi(T,0)y_{2}(0)$ for a time integration of the system equations (equation (15)) over one period $T$ with initial conditions ${\mathfrak{y}}_{2}(0)$ . Hence, without knowing the state transition matrix, it is possible to obtain the product of it with an arbitrary vector (the initial state vector) by integration of equation (15) over one period. The Arnoldi algorithm25 is a method to approximate the eigenvalues and the eigenvectors of a matrix, say $\Phi(T,0)$ , using only the matrix multiplication with $\Phi(T,0)$ to construct an $m$ -sized subspace $$ \mathbf{P}=\left[p_{1}\quad p_{2}\quad\ldots\quad p_{m}\right] $$ that satisfies the orthonormality condition $$ \mathbf{P}^{\mathrm{T}}\mathbf{P}=\mathbf{I} $$ and where the eigenvalues $\tilde{\rho}_{k}$ of the subspace projected state transition matrix $$ \mathbf{H}=\mathbf{P}^{\mathrm{T}}\Phi(T,0)\mathbf{P} $$ converge towards the eigenvalues $\rho_{k}$ of $\Phi(T,0)$ with the largest modulus as the size $m$ of the subspace increases. The subspace eigenvectors $\tilde{w}_{k}$ of $\mathbf{H}$ projected back to the full state space converge towards the eigenvectors $w_{k}$ of $\Phi(T,0)$ , i.e. $w_{k}\approx\mathbf{P}\tilde{w}_{k}$ . The Arnoldi algorithm proceeds as follows: Choose an arbitrary vector $\pmb{p}_{1}$ with $|p_{1}|=1$ for $n=1,2,\ldots,m$ $\pmb{a}:=\Phi(T,0)p_{n}$ (integration of equation (15) over $t\in[0;T])$ $\begin{array}{l}{b:=a}\\ {\mathrm{for~}j=1,2,\dotsc,n}\\ {\quad h_{j,n}:=p_{j}^{\operatorname{T}}a}\\ {\quad b:=b-h_{j,n}p_{j}}\end{array}$ end if $n