Robust Switching Control and Subspace Identification for Flutter of Flexible

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.


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][4][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 closedloop 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  parameterization method which make the frequency range selectable to identify the closedloop 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  +   < 0 according to the theory of Lyapunov stability of linear system, that educes a Lyapunov function    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][23][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.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 Table 1.

Wind Tunnel Model and Aeroelastic Equation
The 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. 3 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: where , , ,  are the strain, stress, electric displacement, and electric field, respectively, and   , ,   denote elastic compliance and piezoelectric constant.

Dynamic Equation of Aeroelastic System. MFC actuator shown in Figure
The equivalent stress   of the driven force as a result of applied voltage on the cross electrode can be expressed as follows: And the moment applied by the MFC to the wing model can be expressed as follows: where  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 energy  and kinetic energy  of the aeroelastic model can be given by where (, , ) 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 where   (  ,  ∞ ) is unsteady aerodynamic matrix which can be computed via doublet lattice method [25].  ,  ∞ are reducing frequency and Mach number, respectively.By applying Hamilton principle, the variation formulation can be established by 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: Substituting ( 7) into ( 6) and performing the variation operation in terms of , expression of the aeroelastic equation can be obtained as where , ,  are modal mass, modal damping, and modal stiffness matrix, namely, = ∭ (, , )  2   (, , )   , and   (  ,  ∞ ) is the unsteady aerodynamic matrix under modal coordinate and () is the dynamic system from the control output to the generalized force applied to modal variable.

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 openloop 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: where () ∈  ×1 , () ∈    ×1 , () ∈   0 ×1 are the discrete Fourier transformations of , , .In this study, the aim of system identification is to obtain the system matrices , , , . is the response of signal by the acceleration transducer adhered on the wing model. is the output signal of the controller.  is the generalized operator.For the  th frequency point in frequency domain,   = 2  / denotes the corresponding value of circular frequency.  is sampling frequency.By iteration of ( 10), the high order equation can be obtained by Rewrite  = 1, 2, . . .,  − 1 into matrix form: where   is the generalized observable matrix, which can be expressed as is the lower triangle Toeplitz matrix: Several values in the frequency domain are selected and all satisfied (12), rewriting into matrix form: where  = [(1) (2) ⋅ ⋅ ⋅ ()],  is the number of discrete points in frequency domain, and  V ,  V are the Van Dermond block matrix expressed as: Separating ( 15) into real part and imaginary part, it can be expressed as where To estimate   ,  decomposition is implemented to the matrix of input and output matrix in [ re  re ]  frequency domain: The singular value decomposition is implemented to matrix   22 : By applying  re orthogonal projection to the row space of  re , Consequently we can get the estimation values of   as According to (21), The least square method was utilized to estimate  and ; thus, namely, The transfer function of the control system from () to () is and it can be written in the form of Kronecker product 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.

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 closedloop 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 Figure 4 is taken into consideration.
where Δ ∈  × 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.,  2 ,  2 ,  22 are the state space realization of (), respectively, shown in Figure 4.As the uncertainty Δ() is represented by the form of additive, it can be inferred that  1 = 0,  1 = 0,  11 = 0,  12 = 1,  21 = 1.
Consider a dynamic feedback controller which can be expressed as Substituting ( 27) into ( 26), the closed-loop system can be expressed as where the closed-loop system of uncertain model can be expressed as The purpose of the controller design is to place the poles of the uncertain closed-loop system within a specific region.

Robust Pole Placement Controller Design
LMI Region.A LMI region is any subset  of the complex plane with characteristic function: where ,  are the Hermitian matrices.For the region of left half plane and Re() < , the characteristic function can be expressed as The symbol "<0" denotes that   is negative definite.Following content introduces the detail of realizing (31) for the closed-loop system.
2. Solve (39) for   ,   ,   ,   .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  1 (),  2 (), . . .,   () are synthesized according to the time varying aeroelastic model   () and its additive uncertainty Δ  () under different wind speed, shown in Figure 5.The uncertain switched system of aeroelastic model in state space form can be expressed as follows:  where  is the time varying switching signal in set {1, 2, . . ., } 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.42) have a common Lyapunov function if there exists , ,   ,   ,   ,   .

Description of Experimental Setup of Wind Tunnel Test.
The aeroelastic control system, shown in Figure 6, 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.
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 Figure 7.

Model Identification of Aeroelastic System in Wind Tunnel.
From Section 2.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   () and self-power spectrum   () as follows: 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 Figure 8.
From Figure 8, 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: 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.

Robust Pole Placement Control Test in Wind Tunnel.
Vibration divergence of the wing model occurs when wind speed exceeds the critical value shown in Figure 9; however, it would not increase infinitely.Due to hardening characteristics of the structure, approximate stable vibration appears after a specific magnitude is attained.
The robust pole placement controllers are designed offline based on above three identified models, and each controller is solved by LMI.(49) Three identified state space models with different wind speed are as follows:

Shock and Vibration
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 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 Figure 10.
The closed-loop frequency response functions are shown in Figure 11.
The 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 Figure 12. Figure 13 shows the root locus of open-loop and closedloop 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.(52) 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: 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 Figure 15.
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 Figure 16 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%. 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 as  2 + 2   +  2  ; 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.
The uncertainty of the system is caused by two aspects: change of wind speed and the error of identification.The  ∞ 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  ∞ norm can be considered as less than this value while wind speed varies from 27 to 31 m/s.Similarly,  ∞ 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.

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

𝑥:
coordinate of wing model : coordinate of wing model : coordinate of wing model : S t i ff n e s sc o n s t a n t : Displacement of the wing model : S y m b o lo fv a r i a t i o n (, , ): Modal eigenvector (): Outputvectorofthestatespacemodel * : Abbreviation of symmetric matrix  max : Maxsingularvalueofsystemmatrix.

Figure 4 :
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 ẋ =  +  1  +  2 ,  =  1  +  11  +  12 ,  =  2  +  21  +  22 ,  = Δ,

Figure 5 :
Figure 5: Switching controller model of uncertain system.

Figure 6 :
Figure 6: Block diagram of flutter control system of wing model for wind tunnel tests.
FRF under the condition with 27 m/s wind speed FRF under the condition with 31 m/s wind speed Measured magnitude-frequency Identified magnitudeFRF under the condition with 33 m/s wind speed

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: speed < 28 m/s controller 1, 28 < speed < 32 m/s controller 2, 32 < speed < 35 m/s controller 3.

Figure 10 :
Figure 10: Time domain response of closed-loop system.

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,
loop Second mode open loop Third mode open loop First mode closed loop Second mode closed loop Third mode closed loop

Figure 13 :Figure 14 :Figure 15 :Figure 16 :
Figure 13: Comparison between the root locus of open-loop and closed-loop systems.

Table 1 :
Comparison of numerical and experimental natural frequency.
Theorem 3. A output feedback controller () and a symmetric matrix  > 0 are available such that (32) holds if and only if two symmetric matrices ,  and state space matrices   ,   ,   ,   are existence, namely,