# Robust Switching Control and Subspace Identification for Flutter of Flexible Wing

Authors: Yizhe Wang; Zhiwei Xu; Wei Li
Journal: Shock and Vibration
 (2018)

---

## Abstract

Active flutter suppression and subspace identification for a flexible wing model using micro fiber composite actuator were experimentally studied in a low speed wind tunnel. NACA0006 thin airfoil model was used for the experimental object to verify the performance of identification algorithm and designed controller. The equation of the fluid, vibration, and piezoelectric coupled motion was theoretically analyzed and experimentally identified under the open-loop and closed-loop condition by subspace method for controller design. A robust pole placement algorithm in terms of linear matrix inequality that accommodates the model uncertainty caused by identification deviation and flow speed variation was utilized to stabilize the divergent aeroelastic system. For further enlarging the flutter envelope, additional controllers were designed subject to the models beyond the flutter speed. Wind speed was measured online as the decision parameter of switching between the controllers. To ensure the stability of arbitrary switching, Common Lyapunov function method was applied to design the robust pole placement controllers for different models to ensure that the closed-loop system shared a common Lyapunov function. Wind tunnel result showed that the designed controllers could stabilize the time varying aeroelastic system over a wide range under arbitrary switching.


---

## Body

## 1. Introduction

Aeroelastic problem has been always attracting extensive attentions in past few decades. Along with the rapid development of aeronautical science, the requirements of endurance and maneuverability are constantly growing, leading to lightweight and flexible characteristics so that the aeroelastic problem is more noticeable. Flutter is a self-excited phenomenon that occurs when fluid couples energy into structure motion causing instability and excessive vibration. So flutter suppression is one of the major objective problems in aeroelastic control.Over the past decades, many research results have been obtained by using different actuator and algorithm to suppress flutter. In early 1980s, NASA Landley center launched active flexible wing plan (AFW) and Benchmark technology projects to suppress the wing model flutter under subsonic and transonic case by using flap [1, 2]. In recent years, a series of studies are conducted to show the efficiency of flap to suppress the bending torsion flutter and limit cycle oscillation [3–5].Piezoelectric material is also a major topic of discussion due to its advance of broadband frequency range and ease to be embedded. Guan et al. [6] conduct a flutter suppression experiment to the lifting surface by using distributed piezoelectric actuator. Since the driving force of piezoelectric ceramic is relatively small, advanced actuators derived from piezoelectric material are applied to the aeroelastic control. Ardelean et al. [7] designed a V-stack piezoelectric actuator to suppress the subsonic flutter of wing section which is capable of producing large force and stroke to drive the flap over a wide frequency range. Burnham and Pitt [8] conduct the research on buffeting alleviation for the vertical tail by using the macro fiber composite actuator, which is proved to be sufficiently forceful and less influential to the structure.In addition to the experimental research, control algorithm is also playing an important role in active flutter suppression. Due to the variation of flying condition and model error, uncertainty of system model is inevitable. In order to cope with the impact of uncertain, robust control algorithms are applied to ensure the stability and performance of the closed-loop system. Han et al. thought over the parameter uncertainty and structure uncertainty of aeroelastic system and designed the LQG and controller synthesis, respectively, by comparing two schemes implemented on a vertical thin plate, exhibiting a better performance and scope of action to the flutter control which improve the closed-loop critical speed [9]. Zeng et al. designed a controller for a semifuselage wing subjecting to the model uncertain and gust distribution based on ARMA model and the dynamic pressure could be improved in the presence of uncertainty [10]. Blue and Balas construct the LPV model by taking the dynamic pressure and Mach number as parameter perturbation, by separating the nominal model and uncertain model using linear fractional transform; the LPV controller was designed to suppress flutter of time varying system [11]. A unique controller could only be effective in a limit range. To enlarge the flutter envelope, Huang et al. employed the closed-loop identification to obtain the dynamic models beyond the critical speed and robust controllers are synthesized based on the corresponding model to suppress the wing model within the range between two identified models [12]. McEver et al. conducted Q parameterization method which make the frequency range selectable to identify the closed-loop model and suppressed the flutter by PPF controller [13].Although the corresponding controllers designed for different models are mentioned above, the switching process between each controller has not been fully considered. Even if every controller ensures the performance under the static condition, the stability and performance under arbitrary switching could not been guaranteed [14]. The system stability is ensured by satisfying the inequality AX+XAT<0 according to the theory of Lyapunov stability of linear system, that educes a Lyapunov function xTXx to be negative definite. However, for the switching system, the dynamic behavior is substantially nonlinear, which cannot be analyzed by the primitive theory of Lyapunov stability. The Lyapunov function of the switching system has a risk to increase at the time of switching. Therefore, it is essential to study the active flutter suppression under switching process.The stability of arbitrary switching control is so far attributed to the common Lyapunov function method [15, 16], multiple Lyapunov function method [17, 18], and dwell time method [19, 20]. The common Lyapunov method ensures that if each subsystem shares a common Lyapunov function, then the switching system could maintain the stability under arbitrary switching. Switching synthesis considers the design of controllers as a unified process, each controller is designed relatively to ensure the performance under random switching. The common Lyapunov function theory has been used for analyzing the individual linear model and partial linear model such as fuzzy T-S system [21]. Since the common Lyapunov function based controller design is essentially conservative [15], multiple Lyapunov function method is proposed. The switching stability can be ensured even if each subsystem has an individual Lyapunov function under certain conditions [22]. Dwell time method is emphasis on the critical switching time for system stability. It can be precisely computed based on exact dynamic model and the conservative of controller design is significantly decreased while the switching speed is less than dwell time. For slow switching system, several efficient methods are proposed based on dwell time theory [22–24].In this study, system analysis and controller design are all based on the common Lyapunov theory since the dynamic model of the aeroelastic system is time varying and uncertain so that the certain conditions based on exact models are not available. The contents of this paper are as follows: dynamic model of flutter is obtained by system identification, and then multicontrollers are designed by robust pole placement. Dynamic model beyond the critical wind speed could be obtained by using frequency domain subspace method under closed-loop conditions. The common Lyapunov function and liner matrix inequality are used to correlate with each controller to ensure the arbitrary switching stability. Wind tunnel test is implemented to indicate that the flutter of the wing under a specific wind speed can be suppressed by each robust controller and stabilization of arbitrary switching among controllers can be ensured perfectly. The main contributions of this research are that the switching stability is completely insured subjecting to the wind speed variation and the robust switching control method is effective significantly in suppressing the wing flutter.


## 2. Wind Tunnel Model and Aeroelastic Equation

### 2.1. Model Description

A wing model is designed to verify the control algorithm in wind tunnel test. The model is made of elastoplastic with elastic modulus 1.1 Gpa, density 1430kg/m3, and Poisson ratio 0.39 and the cross section of the wing model is NACA0006. The max thickness is 6% of the chord length. A lumped mass of 36 grams is mounted on the tip of the wing in order to reduce the flutter speed in wind tunnel test. The wing root is vertically fixed on the base of the closed wind tunnel. Real model and low speed wind tunnel are shown in Figure 1. Due to the limitation of maximum wind speed, wind tunnel test of flutter suppression is carried out within the flow speed range below 35 m/s.Figure 1
Wind tunnel mode.The configurations of actuator and sensor are shown in Figure2. MFC is adhered to the root of the wing model as the control actuator and accelerometer are mounted on the tip of the wing as the monitoring and feedback signal.Figure 2
Configuration of the wing model for wind tunnel test.Modal parameter of the wing can be obtained via FEM analysis or system identification test. The modal analysis is implemented in Nastran FEM software subjected to the first four modes. In system identification test, a low-pass white noise is applied to the MFC as input and the accelerometer signal is measured as output. The natural frequency of the wing model can be identified from the input and output signals via subspace method that will be mentioned later. Comparison between the numerical and experimental results is shown in Table1.Table 1
Comparison of numerical and experimental natural frequency.
Mode Numerical Experimental Mode shape 1 1.99 Hz 2.42 Hz First bending 2 8.60 Hz 11.46 Hz Second bending 3 23.47 Hz 25.93 Hz Third bending 4 59.24 Hz 67.11 Hz First torsionThe difference between the numerical and experimental results is contributed to the deviation of the material parameters and the neglection of MFC actuation. Analysis and test results later show that the flutter occurs in the third mode at wind speed 28 m/s.


### 2.2. Dynamic Equation of Aeroelastic System

MFC actuator shown in Figure3 is internally arranged by a combination of piezoelectric fibers, which can be applied to the piezoelectric material by means of a cross electrode to produce a large driving force. The constitutive equations of piezoelectric materials can be expressed as follows:(1)ε=sEσ+dtE,D=dσ+εTE,where ε, σ, D, E are the strain, stress, electric displacement, and electric field, respectively, and sE, d, εT denote elastic compliance and piezoelectric constant.Figure 3
Macro fibers composite.The equivalent stressσs of the driven force as a result of applied voltage on the cross electrode can be expressed as follows:(2)σs=cεE=cdtE.And the moment applied by the MFC to the wing model can be expressed as follows:(3)Mt=∭z-σsdx-dy-dz-,where z denote the distance between stress point and middle layer of the airfoil. The integral domain includes all the areas covered by the piezoelectric material.To determine the aeroelastic equation, Hamilton principle is used derived from the virtual work of potential energy, kinetic energy, and generalized external force. The total potential energyU and kinetic energy T of the aeroelastic model can be given by(4)T=12∭μx-,y-,z-w2˙dx-dy-dz-,U=12∭σεdx-dy-dz-,where μ(x-,y-,z-) is density variable and is displacement variable. The generalized external force of the system is composed of MFC actuator force and aerodynamic force, which can be expressed as(5)W=∭cdtTE+Qxxkr,M∞wdx-dy-dz-,where Qxxkr,M∞ is unsteady aerodynamic matrix which can be computed via doublet lattice method [25]. kr,M∞ are reducing frequency and Mach number, respectively. By applying Hamilton principle, the variation formulation can be established by(6)∫δT-U+Wdt=0.In order to obtain the equations of motion in modal space, coordinate transformation of the displacement variable is carried out. Displacement of the wing model is expressed in terms of generalized coordinates:(7)w=∑i=1nφix-,y-,z-xi.Substituting (7) into (6) and performing the variation operation in terms of w, expression of the aeroelastic equation can be obtained as(8)Mx′¨+Cx′˙+Kx′=Fu+Qξξkr,M∞x′,where M, C, K are modal mass, modal damping, and modal stiffness matrix, namely,(9)M=mξ1mξ2⋱mξ1,C=cξ1cξ2⋱cξ1,K=kξ1kξ2⋱kξ1,mξi=∭μx-,y-,z-φξi2x-,y-,z-dx-dy-dz-,kki=∭σεφξi2x-,y-,z-dx-dy-dz-,u=∭φix-,y-,z-cdTEdx-dy-dz-,Qξξkr,M∞=∭Qxxkr,M∞φξi2x-,y-,z-dx-dy-dz-,and Qξξkr,M∞ is the unsteady aerodynamic matrix under modal coordinate and F(u) is the dynamic system from the control output to the generalized force applied to modal variable.


### 2.3. Identification of Aeroelastic System

In this study, frequency domain subspace identification method is used to obtain the open-loop and closed-loop aeroelastic system. The identification algorithm is implemented based on the frequency response data to estimate system matrix of state space model and it exhibits a reliable performance for open-loop system and closed-loop system identification [26].Computation simplicity is the major advantage over other method. Frequency domain estimation ensures that the identified data be selectable near the flutter frequency that provides an accurate information. Sufficient accuracy can be obtained if the outside disturbance is small enough. The state space description of multi-input and multi-output system in frequency domain can be expressed as:(10)ξkXk=AXk+BUk,Yk=CXk+DUk,where X(k)∈Cm×1, U(k)∈CNi×1, Y(k)∈CN0×1 are the discrete Fourier transformations of x, y, u. In this study, the aim of system identification is to obtain the system matrices A, B, C, D. y is the response of signal by the acceleration transducer adhered on the wing model. u is the output signal of the controller. ξk is the generalized operator. For the kth frequency point in frequency domain, ξk=j2πfsk/N denotes the corresponding value of circular frequency. fs is sampling frequency. By iteration of (10), the high order equation can be obtained by(11)ξkpYk=ξkp-1CξkXk+ξkDUk=ξkp-1CAXk+CBUk+ξkDUk=⋮CApCXk+CAp-1BUk+CAp-2BUkξk+⋯+CBUkξkp-1+DUkξkpUk.Rewrite p=1,2,…,r-1 into matrix form:(12)YkξkYk⋮ξkr-1kYk=OrXk+GUkξkUk⋮ξkr-1Uk,where Or is the generalized observable matrix, which can be expressed as(13)Or=CCA⋮CAr-1.G is the lower triangle Toeplitz matrix:(14)G=D0⋯0CBD⋯0⋮⋮⋮⋮CAr-2BCAr-3B⋯D.Several values in the frequency domain are selected and all satisfied (12), rewriting into matrix form:(15)Yv=OrX+GUv,where X=[X(1)X(2)⋯X(N)], N is the number of discrete points in frequency domain, and Yv, Uv are the Van Dermond block matrix expressed as:(16)Yv=Y1Y2⋯YNξkY1ξkY2⋯ξkYN⋮⋮⋮⋮ξkr-1ξkr-1Y2⋯ξkr-1YN,Uv=U1U2⋯UNξkU1ξkU1⋯ξkUN⋮⋮⋮⋮ξkr-1U1ξkr-1U1⋯ξkr-1UN.Separating (15) into real part and imaginary part, it can be expressed as(17)Yre=OrXre+GUre,where Yre=[Re⁡YvIm⁡Yv],  Ure=Re⁡UvIm⁡Uv.To estimateOr, QR decomposition is implemented to the matrix of input and output matrix in UreYreT frequency domain:(18)UreYre=R11T0R12TR22TQ1TQ2T.The singular value decomposition is implemented to matrix R22T:(19)R22T=UsU0Ss00S0VsTV0T.By applying  Yre orthogonal projection to the row space of  Ure, (20)YreUre,⊥=OrXreUre,⊥=R22TQ2T.Consequently we can get the estimation values of  Or as(21)Or^=Us.According to (21),(22)Or^1:Nor-1,  :A^=Or^No+1:Nor,  :.The least square method was utilized to estimate A and C; thus, namely, (23)A^=Or^1:Nor-1,  :+Or^No+1:Nor,  :,C^=O^1:No,  :.The transfer function of the control system from U(k) to Y(k) is(24)Yk=ICξI^-A-1D^B^Uk,and it can be written in the form of Kronecker product(25)Yk=Uk⊗IC^ξI-A^-1vecD^B^=ψkθ,where vec(∙) represents a column vector of the matrix which is arranged in order and then B^ and D^ can be obtained by using least square method.


## 2.1. Model Description

A wing model is designed to verify the control algorithm in wind tunnel test. The model is made of elastoplastic with elastic modulus 1.1 Gpa, density 1430kg/m3, and Poisson ratio 0.39 and the cross section of the wing model is NACA0006. The max thickness is 6% of the chord length. A lumped mass of 36 grams is mounted on the tip of the wing in order to reduce the flutter speed in wind tunnel test. The wing root is vertically fixed on the base of the closed wind tunnel. Real model and low speed wind tunnel are shown in Figure 1. Due to the limitation of maximum wind speed, wind tunnel test of flutter suppression is carried out within the flow speed range below 35 m/s.Figure 1
Wind tunnel mode.The configurations of actuator and sensor are shown in Figure2. MFC is adhered to the root of the wing model as the control actuator and accelerometer are mounted on the tip of the wing as the monitoring and feedback signal.Figure 2
Configuration of the wing model for wind tunnel test.Modal parameter of the wing can be obtained via FEM analysis or system identification test. The modal analysis is implemented in Nastran FEM software subjected to the first four modes. In system identification test, a low-pass white noise is applied to the MFC as input and the accelerometer signal is measured as output. The natural frequency of the wing model can be identified from the input and output signals via subspace method that will be mentioned later. Comparison between the numerical and experimental results is shown in Table1.Table 1
Comparison of numerical and experimental natural frequency.
Mode Numerical Experimental Mode shape 1 1.99 Hz 2.42 Hz First bending 2 8.60 Hz 11.46 Hz Second bending 3 23.47 Hz 25.93 Hz Third bending 4 59.24 Hz 67.11 Hz First torsionThe difference between the numerical and experimental results is contributed to the deviation of the material parameters and the neglection of MFC actuation. Analysis and test results later show that the flutter occurs in the third mode at wind speed 28 m/s.


## 2.2. Dynamic Equation of Aeroelastic System

MFC actuator shown in Figure3 is internally arranged by a combination of piezoelectric fibers, which can be applied to the piezoelectric material by means of a cross electrode to produce a large driving force. The constitutive equations of piezoelectric materials can be expressed as follows:(1)ε=sEσ+dtE,D=dσ+εTE,where ε, σ, D, E are the strain, stress, electric displacement, and electric field, respectively, and sE, d, εT denote elastic compliance and piezoelectric constant.Figure 3
Macro fibers composite.The equivalent stressσs of the driven force as a result of applied voltage on the cross electrode can be expressed as follows:(2)σs=cεE=cdtE.And the moment applied by the MFC to the wing model can be expressed as follows:(3)Mt=∭z-σsdx-dy-dz-,where z denote the distance between stress point and middle layer of the airfoil. The integral domain includes all the areas covered by the piezoelectric material.To determine the aeroelastic equation, Hamilton principle is used derived from the virtual work of potential energy, kinetic energy, and generalized external force. The total potential energyU and kinetic energy T of the aeroelastic model can be given by(4)T=12∭μx-,y-,z-w2˙dx-dy-dz-,U=12∭σεdx-dy-dz-,where μ(x-,y-,z-) is density variable and is displacement variable. The generalized external force of the system is composed of MFC actuator force and aerodynamic force, which can be expressed as(5)W=∭cdtTE+Qxxkr,M∞wdx-dy-dz-,where Qxxkr,M∞ is unsteady aerodynamic matrix which can be computed via doublet lattice method [25]. kr,M∞ are reducing frequency and Mach number, respectively. By applying Hamilton principle, the variation formulation can be established by(6)∫δT-U+Wdt=0.In order to obtain the equations of motion in modal space, coordinate transformation of the displacement variable is carried out. Displacement of the wing model is expressed in terms of generalized coordinates:(7)w=∑i=1nφix-,y-,z-xi.Substituting (7) into (6) and performing the variation operation in terms of w, expression of the aeroelastic equation can be obtained as(8)Mx′¨+Cx′˙+Kx′=Fu+Qξξkr,M∞x′,where M, C, K are modal mass, modal damping, and modal stiffness matrix, namely,(9)M=mξ1mξ2⋱mξ1,C=cξ1cξ2⋱cξ1,K=kξ1kξ2⋱kξ1,mξi=∭μx-,y-,z-φξi2x-,y-,z-dx-dy-dz-,kki=∭σεφξi2x-,y-,z-dx-dy-dz-,u=∭φix-,y-,z-cdTEdx-dy-dz-,Qξξkr,M∞=∭Qxxkr,M∞φξi2x-,y-,z-dx-dy-dz-,and Qξξkr,M∞ is the unsteady aerodynamic matrix under modal coordinate and F(u) is the dynamic system from the control output to the generalized force applied to modal variable.


## 2.3. Identification of Aeroelastic System

In this study, frequency domain subspace identification method is used to obtain the open-loop and closed-loop aeroelastic system. The identification algorithm is implemented based on the frequency response data to estimate system matrix of state space model and it exhibits a reliable performance for open-loop system and closed-loop system identification [26].Computation simplicity is the major advantage over other method. Frequency domain estimation ensures that the identified data be selectable near the flutter frequency that provides an accurate information. Sufficient accuracy can be obtained if the outside disturbance is small enough. The state space description of multi-input and multi-output system in frequency domain can be expressed as:(10)ξkXk=AXk+BUk,Yk=CXk+DUk,where X(k)∈Cm×1, U(k)∈CNi×1, Y(k)∈CN0×1 are the discrete Fourier transformations of x, y, u. In this study, the aim of system identification is to obtain the system matrices A, B, C, D. y is the response of signal by the acceleration transducer adhered on the wing model. u is the output signal of the controller. ξk is the generalized operator. For the kth frequency point in frequency domain, ξk=j2πfsk/N denotes the corresponding value of circular frequency. fs is sampling frequency. By iteration of (10), the high order equation can be obtained by(11)ξkpYk=ξkp-1CξkXk+ξkDUk=ξkp-1CAXk+CBUk+ξkDUk=⋮CApCXk+CAp-1BUk+CAp-2BUkξk+⋯+CBUkξkp-1+DUkξkpUk.Rewrite p=1,2,…,r-1 into matrix form:(12)YkξkYk⋮ξkr-1kYk=OrXk+GUkξkUk⋮ξkr-1Uk,where Or is the generalized observable matrix, which can be expressed as(13)Or=CCA⋮CAr-1.G is the lower triangle Toeplitz matrix:(14)G=D0⋯0CBD⋯0⋮⋮⋮⋮CAr-2BCAr-3B⋯D.Several values in the frequency domain are selected and all satisfied (12), rewriting into matrix form:(15)Yv=OrX+GUv,where X=[X(1)X(2)⋯X(N)], N is the number of discrete points in frequency domain, and Yv, Uv are the Van Dermond block matrix expressed as:(16)Yv=Y1Y2⋯YNξkY1ξkY2⋯ξkYN⋮⋮⋮⋮ξkr-1ξkr-1Y2⋯ξkr-1YN,Uv=U1U2⋯UNξkU1ξkU1⋯ξkUN⋮⋮⋮⋮ξkr-1U1ξkr-1U1⋯ξkr-1UN.Separating (15) into real part and imaginary part, it can be expressed as(17)Yre=OrXre+GUre,where Yre=[Re⁡YvIm⁡Yv],  Ure=Re⁡UvIm⁡Uv.To estimateOr, QR decomposition is implemented to the matrix of input and output matrix in UreYreT frequency domain:(18)UreYre=R11T0R12TR22TQ1TQ2T.The singular value decomposition is implemented to matrix R22T:(19)R22T=UsU0Ss00S0VsTV0T.By applying  Yre orthogonal projection to the row space of  Ure, (20)YreUre,⊥=OrXreUre,⊥=R22TQ2T.Consequently we can get the estimation values of  Or as(21)Or^=Us.According to (21),(22)Or^1:Nor-1,  :A^=Or^No+1:Nor,  :.The least square method was utilized to estimate A and C; thus, namely, (23)A^=Or^1:Nor-1,  :+Or^No+1:Nor,  :,C^=O^1:No,  :.The transfer function of the control system from U(k) to Y(k) is(24)Yk=ICξI^-A-1D^B^Uk,and it can be written in the form of Kronecker product(25)Yk=Uk⊗IC^ξI-A^-1vecD^B^=ψkθ,where vec(∙) represents a column vector of the matrix which is arranged in order and then B^ and D^ can be obtained by using least square method.


## 3. Controller Design

Flutter is caused by the instability of the aeroelastic system in general and the poles of the system should be located in the left half of complex plane for the linear stability system. Now, pole placement controllers are designed to make the closed-loop system poles in the left half of the complex plane beyond the flutter speed. Since aeroelastic system is time varying according to the change of the Mach number, attack angle, dynamic pressure, and other factors, single controllers cannot match well with the model and would be invalid if the factors vary in large scope. So single controller is limited in suppressing the wing flutter. In this research, several robust pole placement controllers are designed based on the aeroelastic model under several specific wind speed range; therefore each controller can theoretically place the closed-loop system poles within the designated area in complex plane.After the state space model of aeroelastic system is established, linear fractional model of uncertain system shown in Figure4 is taken into consideration.Figure 4
Linear fractional model of uncertain system.


### 3.1. System Modeling

The linear fractional model in form of additive uncertainty is consider as a feed forward dynamic uncertainty which is added to the model that represents the dynamic deviation caused by the model identification and wind speed variation. The state space with additive uncertainty can be expressed as(26)x˙=Ax+B1w+B2u,z=C1x+D11w+D12u,y=C2x+D21w+D22u,w=Δz,where Δ∈Cm×m represents an uncertain complex gain matrix satisfying σmax(Δ)≤γ-1 and γ is a prescribed constant which can be explained as the peak error of the model in frequency domain. A, B2, C2, D22 are the state space realization of P(s), respectively, shown in Figure 4. As the uncertainty Δ(s) is represented by the form of additive, it can be inferred that B1=0,C1=0,D11=0,D12=1,D21=1.Consider a dynamic feedback controller which can be expressed as(27)x˙k=Akxk+Bky,u=Ckxk+Dky.Substituting (27) into (26), the closed-loop system can be expressed as(28)x˙=AclΔx,where the closed-loop system of uncertain model can be expressed as(29)AclΔ=Acl+BclI-ΔDcl-1ΔCcl,Acl=A+B2DKC2B2CKBKC2AK,Bcl=B1+B2DKD21BKD21,Ccl=C1+D12DKC2D12CK,Dcl=D11+D12DKD21.The purpose of the controller design is to place the poles of the uncertain closed-loop system within a specific region.


### 3.2. Robust Pole Placement Controller Design

LMI Region. A LMI region is any subset D of the complex plane with characteristic function:(30)fD=L+zM+z-MT,where L, M are the Hermitian matrices. For the region of left half plane and Rez<α, the characteristic function can be expressed as(31)fD=z+z-+2α<0.The symbol “<0” denotes that fD is negative definite. Following content introduces the detail of realizing (31) for the closed-loop system.Definition 1.
The uncertain system (27) is robustly D-stable if the eigenvalues of A(Δ) lie in D for all admissible uncertainties Δ.Lemma 2 ([27], robust D-stability).
The sufficient and necessary conditions for the robustD-stability of the closed-loop system are that a positive definite symmetric matrix which satisfies with (32) is existence:(32)MDAc1,X∗∗M1⊗Bc1TX-γI∗M2⊗Cc1I⊗Dc1-γIMDA,X=L⊗X+M⊗XA+MT⊗XAT.Equation (32) is a nonlinear matrix inequality for X,Ak,Bk,Ck,Dk; by using variable substitution, it can be converted to the linear form.Theorem 3.
A output feedback controllerK(s) and a symmetric matrix X>0 are available such that (32) holds if and only if two symmetric matrices R,S and state space matrices Ak,Bk,Ck,Dk are existence, namely,(33)ΛR,S=RIIS>0,(34)L⊗ΛR,S+M⊗ΦAMT⊗ΦAT∗∗M1⊗ΦBT-γI∗M2⊗ΦCI⊗ΦD-γI,(35)ΦA=AR+B2CK¯A+B2DKC2AK¯SA+BK¯C2,ΦB=B1+B2DKD21SB1+BK¯D21,ΦC=C1R+D12CK¯C1+D12DKC2,ΦD=D11+D12DKD21.Proof.
The symmetric positive definite matrixX that satisfied (32) can be expressed as the following form:(36)X=SNNTW,X-1=RMMTZ.Two matrices Π1 and Π2 are defined by(37)Π1=RIMT0,Π2=IS0NT.And it can be verified that(38)XΠ2=Π1,I-RS=MTN.It has been proven [28] that the full rank matrices M, N are necessarily existence if X satisfied (32), then it can be inferred that Π1 and Π2 are invertible. Pre- and postmultiplied X by Π2T and Π2, this yields Π2TΠ1=R00S>0. Define the variable substitution equation:(39)DK=D-K,BK=N-1B-K-N-1SB2DK,CK=C-KMT-1-DKC2RMT-1,AK=N-1A-KMT-1-BKC2RMT-1-N-1SB2CK-N-1SA+B2DKC2RMT-1.
Pre-and postmultiply (32) by diag⁡(Π2T,I,I) and diag⁡(Π2,I,I) and pre-and postmultiply X by Π2T and Π2, substituting (39) and I-RS=MTN into the result, then (33) and (34) are satisfied:(40)Π2TIIL⊗X+M⊗XA+MT⊗XAT∗∗M1⊗Bc1TX-γI∗M2⊗Cc1I⊗Dc1-γIΠ2II=L⊗ΛR,S+M⊗ΦA+MT⊗ΦAT∗∗M1⊗ΦBT-γI∗M2⊗ΦCI⊗ΦD-γI<0.
SinceX>0, Π2TXΠ2=Λ(R,S)>0, and the necessity is proved.Similarly, pre- and postmultiply (33) by (Π2T)-1 and (Π2)-1 and pre- and postmultiply (34) by diag⁡((Π2T)-1,I,I) and diag⁡((Π2)-1,I,I): (41)Π2T-1IIL⊗ΛR,S+M⊗ΦA+MT⊗ΦAT∗∗M1⊗ΦBT-γI∗M2⊗ΦCI⊗ΦD-γIΠ2-1II=L⊗X+M⊗XA+MT⊗XAT∗∗M1⊗Bc1TX-γI∗M2⊗Cc1I⊗Dc1-γI<0,Π2T-1ΛR,SΠ2-1=X>0.If the above LMIs are feasible, the matrices Ak, Bk, Ck,Dk are derived as follows:1.
Calculate the matricesM, N based on MNT=I-RS.2.
Solve (39) for Ak, Bk, Ck, Dk.


### 3.3. Switching Controller Design

Equations (33), (34), and (39) give the synthesis process of robust controller. However, the robust controller will lose efficiency if the system uncertainty is larger than γ. A group of switching stabilized controllers K1s,K2s,…,Kns are synthesized according to the time varying aeroelastic model Gi(s) and its additive uncertainty Δi(s) under different wind speed, shown in Figure 5.Figure 5
Switching controller model of uncertain system.The uncertain switched system of aeroelastic model in state space form can be expressed as follows:(42)x˙=Aσx+Bσ1w+Bσ2u,z=Cσ1x+Dσ11w+Dσ12u,y=Cσ2x+Dσ21w+Dσ22u,w=Δσz,σ=1,2,…,n,where σ is the time varying switching signal in set 1,2,…,n that represents different model under the time variant condition. In this study, wind speed is considered as the characteristic condition which arbitrarily varied with time. The switching signal is determined by decision variable ψ which reflect the varying of wind speed. Each subsystem of (43) represents the state space model of different wind speed. A set of controllers are designed by the following theory to ensure the stability during arbitrary process.Theorem 4.
System equations (42) have a common Lyapunov function if there exists R,S,A¯σK,B¯σK,C¯σK,DσK.
That satisfied(43)ΛR,S=RIIS>0,(44)Lσ⊗ΛR,S+M⊗ΦAσ+MT⊗ΦAσT∗∗M1⊗ΦBσT-γσI∗M2⊗ΦCσI⊗ΦDσ-γσI,where(45)ΦAσ=AσR+Bσ2CσK¯Aσ+Bσ2DσKCσ2AσK¯SAσ+BσK¯Cσ2,ΦBσ=Bσ1+Bσ2DσKDσ21SBσ1+BσK¯Dσ21,ΦCσ=Cσ1R+Dσ12CσK¯Cσ1+Dσ12DσKCσ2,ΦDσ=Dσ11+Dσ12DσKDσ21.Proof.
It can be seen that if the closed-loop system (42) is satisfied with (43) and (44), there exists a common R,S and M,N which satisfy with (43). For each subsystem, pre- and postmultiply Λ(R,S) by (Π2T)-1 and Π2T, namely,(46)Π2T-1ΛR,SΠ2-1=X>0,Π2TIILσ⊗ΛR,S+M⊗ΦAσ+MT⊗ΦAσT∗∗M1⊗ΦBσT-γσI∗M2⊗ΦCσI⊗ΦDσ-γσIΠ2II=Lσ⊗X+M⊗XAσ+MT⊗XAσT∗∗M1⊗Bc1σTX-γσI∗M2⊗Cc1σI⊗Dc1σ-γσI<0.From Lemma 2, (47) is equivalent to the inequality AσX+(XAσ)T<-Lσ⊗X<0 which indicated that each subsystem σ shares a common Lyapunov function.


## 3.1. System Modeling

The linear fractional model in form of additive uncertainty is consider as a feed forward dynamic uncertainty which is added to the model that represents the dynamic deviation caused by the model identification and wind speed variation. The state space with additive uncertainty can be expressed as(26)x˙=Ax+B1w+B2u,z=C1x+D11w+D12u,y=C2x+D21w+D22u,w=Δz,where Δ∈Cm×m represents an uncertain complex gain matrix satisfying σmax(Δ)≤γ-1 and γ is a prescribed constant which can be explained as the peak error of the model in frequency domain. A, B2, C2, D22 are the state space realization of P(s), respectively, shown in Figure 4. As the uncertainty Δ(s) is represented by the form of additive, it can be inferred that B1=0,C1=0,D11=0,D12=1,D21=1.Consider a dynamic feedback controller which can be expressed as(27)x˙k=Akxk+Bky,u=Ckxk+Dky.Substituting (27) into (26), the closed-loop system can be expressed as(28)x˙=AclΔx,where the closed-loop system of uncertain model can be expressed as(29)AclΔ=Acl+BclI-ΔDcl-1ΔCcl,Acl=A+B2DKC2B2CKBKC2AK,Bcl=B1+B2DKD21BKD21,Ccl=C1+D12DKC2D12CK,Dcl=D11+D12DKD21.The purpose of the controller design is to place the poles of the uncertain closed-loop system within a specific region.


## 3.2. Robust Pole Placement Controller Design

LMI Region. A LMI region is any subset D of the complex plane with characteristic function:(30)fD=L+zM+z-MT,where L, M are the Hermitian matrices. For the region of left half plane and Rez<α, the characteristic function can be expressed as(31)fD=z+z-+2α<0.The symbol “<0” denotes that fD is negative definite. Following content introduces the detail of realizing (31) for the closed-loop system.Definition 1.
The uncertain system (27) is robustly D-stable if the eigenvalues of A(Δ) lie in D for all admissible uncertainties Δ.Lemma 2 ([27], robust D-stability).
The sufficient and necessary conditions for the robustD-stability of the closed-loop system are that a positive definite symmetric matrix which satisfies with (32) is existence:(32)MDAc1,X∗∗M1⊗Bc1TX-γI∗M2⊗Cc1I⊗Dc1-γIMDA,X=L⊗X+M⊗XA+MT⊗XAT.Equation (32) is a nonlinear matrix inequality for X,Ak,Bk,Ck,Dk; by using variable substitution, it can be converted to the linear form.Theorem 3.
A output feedback controllerK(s) and a symmetric matrix X>0 are available such that (32) holds if and only if two symmetric matrices R,S and state space matrices Ak,Bk,Ck,Dk are existence, namely,(33)ΛR,S=RIIS>0,(34)L⊗ΛR,S+M⊗ΦAMT⊗ΦAT∗∗M1⊗ΦBT-γI∗M2⊗ΦCI⊗ΦD-γI,(35)ΦA=AR+B2CK¯A+B2DKC2AK¯SA+BK¯C2,ΦB=B1+B2DKD21SB1+BK¯D21,ΦC=C1R+D12CK¯C1+D12DKC2,ΦD=D11+D12DKD21.Proof.
The symmetric positive definite matrixX that satisfied (32) can be expressed as the following form:(36)X=SNNTW,X-1=RMMTZ.Two matrices Π1 and Π2 are defined by(37)Π1=RIMT0,Π2=IS0NT.And it can be verified that(38)XΠ2=Π1,I-RS=MTN.It has been proven [28] that the full rank matrices M, N are necessarily existence if X satisfied (32), then it can be inferred that Π1 and Π2 are invertible. Pre- and postmultiplied X by Π2T and Π2, this yields Π2TΠ1=R00S>0. Define the variable substitution equation:(39)DK=D-K,BK=N-1B-K-N-1SB2DK,CK=C-KMT-1-DKC2RMT-1,AK=N-1A-KMT-1-BKC2RMT-1-N-1SB2CK-N-1SA+B2DKC2RMT-1.
Pre-and postmultiply (32) by diag⁡(Π2T,I,I) and diag⁡(Π2,I,I) and pre-and postmultiply X by Π2T and Π2, substituting (39) and I-RS=MTN into the result, then (33) and (34) are satisfied:(40)Π2TIIL⊗X+M⊗XA+MT⊗XAT∗∗M1⊗Bc1TX-γI∗M2⊗Cc1I⊗Dc1-γIΠ2II=L⊗ΛR,S+M⊗ΦA+MT⊗ΦAT∗∗M1⊗ΦBT-γI∗M2⊗ΦCI⊗ΦD-γI<0.
SinceX>0, Π2TXΠ2=Λ(R,S)>0, and the necessity is proved.Similarly, pre- and postmultiply (33) by (Π2T)-1 and (Π2)-1 and pre- and postmultiply (34) by diag⁡((Π2T)-1,I,I) and diag⁡((Π2)-1,I,I): (41)Π2T-1IIL⊗ΛR,S+M⊗ΦA+MT⊗ΦAT∗∗M1⊗ΦBT-γI∗M2⊗ΦCI⊗ΦD-γIΠ2-1II=L⊗X+M⊗XA+MT⊗XAT∗∗M1⊗Bc1TX-γI∗M2⊗Cc1I⊗Dc1-γI<0,Π2T-1ΛR,SΠ2-1=X>0.If the above LMIs are feasible, the matrices Ak, Bk, Ck,Dk are derived as follows:1.
Calculate the matricesM, N based on MNT=I-RS.2.
Solve (39) for Ak, Bk, Ck, Dk.


## 3.3. Switching Controller Design

Equations (33), (34), and (39) give the synthesis process of robust controller. However, the robust controller will lose efficiency if the system uncertainty is larger than γ. A group of switching stabilized controllers K1s,K2s,…,Kns are synthesized according to the time varying aeroelastic model Gi(s) and its additive uncertainty Δi(s) under different wind speed, shown in Figure 5.Figure 5
Switching controller model of uncertain system.The uncertain switched system of aeroelastic model in state space form can be expressed as follows:(42)x˙=Aσx+Bσ1w+Bσ2u,z=Cσ1x+Dσ11w+Dσ12u,y=Cσ2x+Dσ21w+Dσ22u,w=Δσz,σ=1,2,…,n,where σ is the time varying switching signal in set 1,2,…,n that represents different model under the time variant condition. In this study, wind speed is considered as the characteristic condition which arbitrarily varied with time. The switching signal is determined by decision variable ψ which reflect the varying of wind speed. Each subsystem of (43) represents the state space model of different wind speed. A set of controllers are designed by the following theory to ensure the stability during arbitrary process.Theorem 4.
System equations (42) have a common Lyapunov function if there exists R,S,A¯σK,B¯σK,C¯σK,DσK.
That satisfied(43)ΛR,S=RIIS>0,(44)Lσ⊗ΛR,S+M⊗ΦAσ+MT⊗ΦAσT∗∗M1⊗ΦBσT-γσI∗M2⊗ΦCσI⊗ΦDσ-γσI,where(45)ΦAσ=AσR+Bσ2CσK¯Aσ+Bσ2DσKCσ2AσK¯SAσ+BσK¯Cσ2,ΦBσ=Bσ1+Bσ2DσKDσ21SBσ1+BσK¯Dσ21,ΦCσ=Cσ1R+Dσ12CσK¯Cσ1+Dσ12DσKCσ2,ΦDσ=Dσ11+Dσ12DσKDσ21.Proof.
It can be seen that if the closed-loop system (42) is satisfied with (43) and (44), there exists a common R,S and M,N which satisfy with (43). For each subsystem, pre- and postmultiply Λ(R,S) by (Π2T)-1 and Π2T, namely,(46)Π2T-1ΛR,SΠ2-1=X>0,Π2TIILσ⊗ΛR,S+M⊗ΦAσ+MT⊗ΦAσT∗∗M1⊗ΦBσT-γσI∗M2⊗ΦCσI⊗ΦDσ-γσIΠ2II=Lσ⊗X+M⊗XAσ+MT⊗XAσT∗∗M1⊗Bc1σTX-γσI∗M2⊗Cc1σI⊗Dc1σ-γσI<0.From Lemma 2, (47) is equivalent to the inequality AσX+(XAσ)T<-Lσ⊗X<0 which indicated that each subsystem σ shares a common Lyapunov function.


## 4. Wind Tunnel Test of the Wing Flutter Suppression

### 4.1. Description of Experimental Setup of Wind Tunnel Test

The aeroelastic control system, shown in Figure6, consists of embedded controller, accelerometer, wind speed transducer, signal conditioner, power amplifier, and MFC actuator. The acceleration response signal is acquired by the sensor and fed back to the controller. NI PXIe-1071 controller is used as the measured signal processing and control algorithm implementation. PXIe-6356 16 bit input and output module which is integrated with the controller performs a high speed and precision measurement that ensure supporting the control performance. Real-time performance of the control system can be obtained by the NI LabVIEW software and high calculation of the controller. The control and signal processing algorithm are programmed on LabVIEW language. Control signal derived from the control algorithm operation is applied to MFC actuator via the power amplifier. The controller is designed at wind speeds 27, 31, and 33 m/s, respectively, and discretized at the sampling time 0.002 s. Experimental data are recorded and displayed on the platform to verify the control performance.Figure 6
Block diagram of flutter control system of wing model for wind tunnel tests.In order to verify the dynamic characteristic and stability under arbitrary switching of the closed-loop system, wind speed is manually regulated by the knob on the control cabinet of wind tunnel, shown in Figure7.Figure 7 (a)
Experimental apparatus of wind tunnel test (b)
Control cabinet of wind tunnel speed


### 4.2. Model Identification of Aeroelastic System in Wind Tunnel

From Section2.3 we know that the state space model of aeroelastic system can be obtained by frequency domain subspace method. A 25 Hz to 45 Hz band pass signal was used as the input signal which is applied to the MFC actuator via power amplifier and the accelerometer response signal is used as the output signal of the model. In this study, frequency response function of the input and output signals is estimated as the identified data by the ratio of cross power spectrum pyu(ω) and self-power spectrum puu(ω) as follows:(47)Gjω=pyuωpuuω.Power spectrum can be estimated based on the time domain data via Welch method. Three frequency response models are estimated under the condition with wind speeds 27, 31, and 33 m/s, respectively. The identification results of three frequency response function of input and output signals are shown in Figure8.Figure 8
FRF of identified aeroelastic model.
(a)
FRF under the condition with 27 m/s wind speed (b)
FRF under the condition with 31 m/s wind speed (c)
FRF under the condition with 33 m/s wind speedFrom Figure8, we find that the major energy of flutter is concentrated near the range of flutter frequency. To simplify the process of the controller model, all modes except for the flutter mode will be neglected in identification process. A second-order state space model is used to approximate the aeroelastic system by mean of optimal fitting in frequency domain. The identified models in transfer function forms of three different wind speeds by subspace method are as follows:(48)G1s=1.867s+925.2s2+4.621s+38999v=27m/s,G2s=4.555s+981.9s2-3.696s+36293v=31m/s,G3s=-0.315s2+5.158s-11220s2-6.373s+31685v=33m/s.From the comparison of frequency response function between identified and test models, we find that the modulus error and phase error caused by identification deviation and variation of flow speed are existence. It will be compensated in the process of robust controller design.


### 4.3. Robust Pole Placement Control Test in Wind Tunnel

Vibration divergence of the wing model occurs when wind speed exceeds the critical value shown in Figure9; however, it would not increase infinitely. Due to hardening characteristics of the structure, approximate stable vibration appears after a specific magnitude is attained.Figure 9
Divergence signal.The robust pole placement controllers are designed offline based on above three identified models, and each controller is solved by LMI.By assigning the pole position to a specific region and estimating the uncertainty according to the frequency response, the robust pole placement controller under the condition witch 27 m/s, 31 m/s and 33 m/s is designed by setting(49)λ1=0.6,α1=1.5,λ2=0.8,α2=1.5,λ3=0.8,α3=1.5.Three identified state space models with different wind speed are as follows:(50)A1=33.83167.4-240.9-39.44,B1=0.12960.1119,C1=0.06356-0.1518,v<27m/s,A2=6.298366.5-486.5-166.3,B2=2.17434.84,C2=10.78-89.63,v<31m/s,A3=20.91153.5-250-76.18,B3=2.5987.322,C3=3.271-9.249,v<33m/s.Before flutter suppression experiment, the designed controller should be discretized. Zero-order hold transformation is used to convert the analog controller to the form of difference equation. The discrete transfer functions of designed controller in three wind speed are(51)K1z=-1.831×10-6+1.308×10-5z-1+1.491×10-5z-21-1.841z-1+0.989z-2v<27m/s,K2z=-0.994+0.206z-1+1.211z-21-0.589z-1-0.173z-2v<31m/s,K3z=-0.043+0.022z-1+0.065z-21-1.764z-1+0.899z-2v<33m/s.To demonstrate the effectiveness of designed controller, three tests are implemented with 28 m/s, 31 m/s, and 33 m/s single wind speed. During the test, controllers are switched off at first and then vibration is converged quickly after the controllers were switched on, shown in Figure10.Figure 10
Time domain response of closed-loop system.
(a)
Time domain response with 28 m/s wind speed (b)
Time domain response with 31 m/s wind speed (c)
Time domain response with 33 m/s wind speedThe closed-loop frequency response functions are shown in Figure11.Figure 11
FRF of closed-loop system.
(a)
FRF under the condition with 28 m/s wind speed (b)
FRF under the condition with 31 m/s wind speed (c)
FRF under the condition with 33 m/s wind speedThe comparisons between the power spectrum of the controlled and uncontrolled response signal illustrate that the major vibration energy of flutter and turbulence is efficiently alleviated by the robust controllers shown in Figure12.Figure 12
Comparison of power spectrum with different speed.
(a)
Power spectrums with 28 m/s wind speed (b)
Power spectrums with 28 m/s wind speed (c)
Power spectrums with 33 m/s wind speedFigure13 shows the root locus of open-loop and closed-loop aeroelastic system with the variation of wind speed identified by subspace method. We find that the aeroelastic system is stabilized by the controller over the flutter speed and all poles of closed-loop system located in stabilization area.Figure 13
Comparison between the root locus of open-loop and closed-loop systems.


### 4.4. Switching Control Test in Wind Tunnel

The switching control signal is obtained by the wind pressure transducer placed on the wall of closed tunnel. Switching rule is expressed by the following expression:(52)speed<28m/scontroller1,28<speed<32m/scontroller2,32<speed<35m/scontroller3.By setting the parametersλ1=0.6, α1=1.5, λ2=0.8, α2=2, λ3=0.8, α3=2 of the switching control and solving the LMI of each subsystem simultaneously, three controllers can share a common Lyapunov function. Then, discretized transfer functions of switching controllers are as follows:(53)Kσ1z=-7.345×10-7+4.672×10-6z-1+5.406×10-6z-21-1.841z-1+0.989z-2v<27m/s,Kσ2z=-0.039+0.019z-1+0.059z-21-1.776z-1+0.912z-2v<31m/s,Kσ3z=-0.071+0.032z-1+0.102z-21-1.739z-1+0.874z-2v<33m/s.By comparison of the damping ratios between open loop and closed loop, it shows that the damping ratio of aeroelastic system is improved by the switching controllers.The wind speed of wind tunnel is varied from the range of 28 m/s to 35 m/s which was controlled by the knob of the panel. The controller is switched dynamically over three candidate controllers and the time domain of the switching signal is shown in Figure15.Since the wind pressure signal was disturbed by high frequency noise generated by wind tunnel, there is a serious risk that the switching signal could be changed quickly. The existence of common Lyapunov function for each subsystem ensures the unconditional stability of arbitrary switching. From Figure16 we find that the switching signal keeps stability while the wind speed changes from 27 m/s to 35 m/s. The vibration amplitude is suppressed about 80%.


### 4.5. Discussion

Subspace identification is used to model the aeroelastic system in this research. Aeroelastic model under different wind speeds 27 m/s, 31 m/s, and 33 m/s is identified though experiments. The model under 27 m/s is below the critical speed and controller designed based on this model can efficiently enlarge the critical speed and increase the Signal-to-Noise ratio for the identification under 31 m/s, and it has the same effect under 31 m/s. When the flow speed is closed to the critical speed, flutter mode plays a dominant role of vibration, so a two-order model is used to approximate the unstable mode. The instability of vibration can be represented by “negative damping” in vibration. The characteristic equation of second-order transform function is expressed ass2+2ξωns+ωn2; as the damping coefficient ξ is negative, the system is unstable. The identification results in Figure 14 show that the root locus is always in the left half part of the complex plane below the critical speed, the damping coefficient is changed from 0.13 to 0.01 and will always be positive. While the flow speed is higher than the critical speed, the root locus locates on the right half part and the damping coefficient is negative. The identification result shows an accurate reflection to the aeroelastic instability.Figure 14
Damping ratio of open and closed-loop systems.Figure 15
Time domain signal of switching control.Figure 16
Switching control test results.The uncertainty of the system is caused by two aspects: change of wind speed and the error of identification. TheH∞ norm of the uncertain dynamic for controller design is based on the measured frequency response function. The modulus deviation of the FRF under 27 and 31 is less than 0.6, so the H∞ norm can be considered as less than this value while wind speed varies from 27 to 31 m/s. Similarly, H∞ norm of uncertainty under the wind speed variation from 31 to 33 m/s is set to be 0.8.In practice, if each controller is individually designed, a global Lyapunov function would not be shared by each closed-loop system so that the convergence can only be satisfied under no switching or slow switching condition. At switching time, total energy of the system is still possible to increase. If the flow speed variation and switching time satisfy a specific condition, vibration is still possible to divergence. The major advantage of common Lyapunov function method is absolutely alleviating the risk of switching instability and ensures the safety of flutter control.


## 4.1. Description of Experimental Setup of Wind Tunnel Test

The aeroelastic control system, shown in Figure6, consists of embedded controller, accelerometer, wind speed transducer, signal conditioner, power amplifier, and MFC actuator. The acceleration response signal is acquired by the sensor and fed back to the controller. NI PXIe-1071 controller is used as the measured signal processing and control algorithm implementation. PXIe-6356 16 bit input and output module which is integrated with the controller performs a high speed and precision measurement that ensure supporting the control performance. Real-time performance of the control system can be obtained by the NI LabVIEW software and high calculation of the controller. The control and signal processing algorithm are programmed on LabVIEW language. Control signal derived from the control algorithm operation is applied to MFC actuator via the power amplifier. The controller is designed at wind speeds 27, 31, and 33 m/s, respectively, and discretized at the sampling time 0.002 s. Experimental data are recorded and displayed on the platform to verify the control performance.Figure 6
Block diagram of flutter control system of wing model for wind tunnel tests.In order to verify the dynamic characteristic and stability under arbitrary switching of the closed-loop system, wind speed is manually regulated by the knob on the control cabinet of wind tunnel, shown in Figure7.Figure 7 (a)
Experimental apparatus of wind tunnel test (b)
Control cabinet of wind tunnel speed


## 4.2. Model Identification of Aeroelastic System in Wind Tunnel

From Section2.3 we know that the state space model of aeroelastic system can be obtained by frequency domain subspace method. A 25 Hz to 45 Hz band pass signal was used as the input signal which is applied to the MFC actuator via power amplifier and the accelerometer response signal is used as the output signal of the model. In this study, frequency response function of the input and output signals is estimated as the identified data by the ratio of cross power spectrum pyu(ω) and self-power spectrum puu(ω) as follows:(47)Gjω=pyuωpuuω.Power spectrum can be estimated based on the time domain data via Welch method. Three frequency response models are estimated under the condition with wind speeds 27, 31, and 33 m/s, respectively. The identification results of three frequency response function of input and output signals are shown in Figure8.Figure 8
FRF of identified aeroelastic model.
(a)
FRF under the condition with 27 m/s wind speed (b)
FRF under the condition with 31 m/s wind speed (c)
FRF under the condition with 33 m/s wind speedFrom Figure8, we find that the major energy of flutter is concentrated near the range of flutter frequency. To simplify the process of the controller model, all modes except for the flutter mode will be neglected in identification process. A second-order state space model is used to approximate the aeroelastic system by mean of optimal fitting in frequency domain. The identified models in transfer function forms of three different wind speeds by subspace method are as follows:(48)G1s=1.867s+925.2s2+4.621s+38999v=27m/s,G2s=4.555s+981.9s2-3.696s+36293v=31m/s,G3s=-0.315s2+5.158s-11220s2-6.373s+31685v=33m/s.From the comparison of frequency response function between identified and test models, we find that the modulus error and phase error caused by identification deviation and variation of flow speed are existence. It will be compensated in the process of robust controller design.


## 4.3. Robust Pole Placement Control Test in Wind Tunnel

Vibration divergence of the wing model occurs when wind speed exceeds the critical value shown in Figure9; however, it would not increase infinitely. Due to hardening characteristics of the structure, approximate stable vibration appears after a specific magnitude is attained.Figure 9
Divergence signal.The robust pole placement controllers are designed offline based on above three identified models, and each controller is solved by LMI.By assigning the pole position to a specific region and estimating the uncertainty according to the frequency response, the robust pole placement controller under the condition witch 27 m/s, 31 m/s and 33 m/s is designed by setting(49)λ1=0.6,α1=1.5,λ2=0.8,α2=1.5,λ3=0.8,α3=1.5.Three identified state space models with different wind speed are as follows:(50)A1=33.83167.4-240.9-39.44,B1=0.12960.1119,C1=0.06356-0.1518,v<27m/s,A2=6.298366.5-486.5-166.3,B2=2.17434.84,C2=10.78-89.63,v<31m/s,A3=20.91153.5-250-76.18,B3=2.5987.322,C3=3.271-9.249,v<33m/s.Before flutter suppression experiment, the designed controller should be discretized. Zero-order hold transformation is used to convert the analog controller to the form of difference equation. The discrete transfer functions of designed controller in three wind speed are(51)K1z=-1.831×10-6+1.308×10-5z-1+1.491×10-5z-21-1.841z-1+0.989z-2v<27m/s,K2z=-0.994+0.206z-1+1.211z-21-0.589z-1-0.173z-2v<31m/s,K3z=-0.043+0.022z-1+0.065z-21-1.764z-1+0.899z-2v<33m/s.To demonstrate the effectiveness of designed controller, three tests are implemented with 28 m/s, 31 m/s, and 33 m/s single wind speed. During the test, controllers are switched off at first and then vibration is converged quickly after the controllers were switched on, shown in Figure10.Figure 10
Time domain response of closed-loop system.
(a)
Time domain response with 28 m/s wind speed (b)
Time domain response with 31 m/s wind speed (c)
Time domain response with 33 m/s wind speedThe closed-loop frequency response functions are shown in Figure11.Figure 11
FRF of closed-loop system.
(a)
FRF under the condition with 28 m/s wind speed (b)
FRF under the condition with 31 m/s wind speed (c)
FRF under the condition with 33 m/s wind speedThe comparisons between the power spectrum of the controlled and uncontrolled response signal illustrate that the major vibration energy of flutter and turbulence is efficiently alleviated by the robust controllers shown in Figure12.Figure 12
Comparison of power spectrum with different speed.
(a)
Power spectrums with 28 m/s wind speed (b)
Power spectrums with 28 m/s wind speed (c)
Power spectrums with 33 m/s wind speedFigure13 shows the root locus of open-loop and closed-loop aeroelastic system with the variation of wind speed identified by subspace method. We find that the aeroelastic system is stabilized by the controller over the flutter speed and all poles of closed-loop system located in stabilization area.Figure 13
Comparison between the root locus of open-loop and closed-loop systems.


## 4.4. Switching Control Test in Wind Tunnel

The switching control signal is obtained by the wind pressure transducer placed on the wall of closed tunnel. Switching rule is expressed by the following expression:(52)speed<28m/scontroller1,28<speed<32m/scontroller2,32<speed<35m/scontroller3.By setting the parametersλ1=0.6, α1=1.5, λ2=0.8, α2=2, λ3=0.8, α3=2 of the switching control and solving the LMI of each subsystem simultaneously, three controllers can share a common Lyapunov function. Then, discretized transfer functions of switching controllers are as follows:(53)Kσ1z=-7.345×10-7+4.672×10-6z-1+5.406×10-6z-21-1.841z-1+0.989z-2v<27m/s,Kσ2z=-0.039+0.019z-1+0.059z-21-1.776z-1+0.912z-2v<31m/s,Kσ3z=-0.071+0.032z-1+0.102z-21-1.739z-1+0.874z-2v<33m/s.By comparison of the damping ratios between open loop and closed loop, it shows that the damping ratio of aeroelastic system is improved by the switching controllers.The wind speed of wind tunnel is varied from the range of 28 m/s to 35 m/s which was controlled by the knob of the panel. The controller is switched dynamically over three candidate controllers and the time domain of the switching signal is shown in Figure15.Since the wind pressure signal was disturbed by high frequency noise generated by wind tunnel, there is a serious risk that the switching signal could be changed quickly. The existence of common Lyapunov function for each subsystem ensures the unconditional stability of arbitrary switching. From Figure16 we find that the switching signal keeps stability while the wind speed changes from 27 m/s to 35 m/s. The vibration amplitude is suppressed about 80%.


## 4.5. Discussion

Subspace identification is used to model the aeroelastic system in this research. Aeroelastic model under different wind speeds 27 m/s, 31 m/s, and 33 m/s is identified though experiments. The model under 27 m/s is below the critical speed and controller designed based on this model can efficiently enlarge the critical speed and increase the Signal-to-Noise ratio for the identification under 31 m/s, and it has the same effect under 31 m/s. When the flow speed is closed to the critical speed, flutter mode plays a dominant role of vibration, so a two-order model is used to approximate the unstable mode. The instability of vibration can be represented by “negative damping” in vibration. The characteristic equation of second-order transform function is expressed ass2+2ξωns+ωn2; as the damping coefficient ξ is negative, the system is unstable. The identification results in Figure 14 show that the root locus is always in the left half part of the complex plane below the critical speed, the damping coefficient is changed from 0.13 to 0.01 and will always be positive. While the flow speed is higher than the critical speed, the root locus locates on the right half part and the damping coefficient is negative. The identification result shows an accurate reflection to the aeroelastic instability.Figure 14
Damping ratio of open and closed-loop systems.Figure 15
Time domain signal of switching control.Figure 16
Switching control test results.The uncertainty of the system is caused by two aspects: change of wind speed and the error of identification. TheH∞ norm of the uncertain dynamic for controller design is based on the measured frequency response function. The modulus deviation of the FRF under 27 and 31 is less than 0.6, so the H∞ norm can be considered as less than this value while wind speed varies from 27 to 31 m/s. Similarly, H∞ norm of uncertainty under the wind speed variation from 31 to 33 m/s is set to be 0.8.In practice, if each controller is individually designed, a global Lyapunov function would not be shared by each closed-loop system so that the convergence can only be satisfied under no switching or slow switching condition. At switching time, total energy of the system is still possible to increase. If the flow speed variation and switching time satisfy a specific condition, vibration is still possible to divergence. The major advantage of common Lyapunov function method is absolutely alleviating the risk of switching instability and ensures the safety of flutter control.


## 5. Conclusion

In this paper, macro fiber composite is used to suppress the flutter of wing model under the variation of wind speed. The frequency subspace identification and robust pole placement switching controller design method are discussed and the efficiency of this method is validated though flutter suppression test. The main contributions of the paper are the following two aspects:(1)
The subspace identification method in frequency domain is utilized to obtain the dynamic model under open and closed-loop condition. Identified model effectively reflects the instability of aeroelastic system. The robust pole placement controllers are designed based on the identified model and solved by the LMIs. Theoretical analysis and wind tunnel test are implemented to prove that the time varying system within a range of uncertainty can be stabilized by the designed controllers.(2)
To enlarge the stability envelope, multiple controllers are applied to suppress the flutter under different wind speed range. Three switching controllers are designed by the common Lyapunov method which ensure that the closed-loop system is stable under arbitrary switching. Theoretical analysis is implemented to prove that the designed controller can stabilize the time varying system within a certain range of uncertainty. To verify the effective of control algorithm, the switching control test is implemented in wind tunnel. Wind speed is manually adjusted to simulate the random change of flow. Experimental result shows that the switching control algorithm can successfully suppress the flutter over a wide range of speed under the arbitrary switching.



---

*Source: 9734684-2018-02-26.xml*