Tang Wei,Wu Jian,Shi Zhongke

School of Automation,Northwestern Polytechnical University,Xi’an 710072,China

Identification of reduced-order model for an aeroelastic system from flutter test data

Tang Wei*,Wu Jian,Shi Zhongke

School of Automation,Northwestern Polytechnical University,Xi’an 710072,China

Aeroelastic system;Flutter test;Maximum likelihood;Reduced-order model;Subspace identification

Recently,flutter active control using linear parameter varying(LPV)framework has attracted a lot of attention.LPV control synthesis usually generates controllers that are at least of the same order as the aeroelastic models.Therefore,the reduced-order model is required by synthesis for avoidance of large computation cost and high-order controller.This paper proposes a new procedure for generation of accurate reduced-order linear time-invariant(LTI)models by using system identification from flutter testing data.The proposed approach is in two steps.The well-known poly-reference least squares complex frequency(p-LSCF)algorithm is firstly employed for modal parameter identification from frequency response measurement.After parameter identification,the dominant physical modes are determined by clear stabilization diagrams and clustering technique.In the second step,with prior knowledge of physical poles,the improved frequencydomain maximum likelihood(ML)estimator is presented for building accurate reduced-order model.Before ML estimation,an improved subspace identification considering the poles constraint is also proposed for initializing the iterative procedure.Finally,the performance of the proposed procedure is validated by real flight flutter test data.

1.Introduction

Among the various phenomena that can affect the flight safety of an aircraft,flutter is one of the most feared since this dynamic instability caused by the interaction of structural,inertial and aerodynamic loads can lead to a catastrophic structural failure within a sufficiently short stretch of time.1,2It is therefore of the utmost importance to prevent flutter from happening and actively suppress it when possible.In the past decade,flutter active control using linear parameter varying(LPV)framework has received a lot of attention.3–7The class of LPV models has been demonstrated to be naturally suitable for describing the operating parameter dependence of the dynamic flutter models across the flight envelope.

Unfortunately,the high-fidelity LPV model derived in industrial context is based on structural finite elements and lifting-surface theory.The resulting LPV model has about hundreds of states and no specific structure.The high-ordermodels prohibit their immediate use in conjunction with modern control theory for controller synthesis.Hence,the reduced-ordermodels(ROM)ofaeroelasticsystem are needed.A variety of model order reduction techniques have been developed to reduce high-order aeroelastic model into a reduced-order state-space model while dominant behavior of the system is retained.8The effort in this field includes balanced realization and truncation9,Krylov-based projection10and hybrid SVD-Krylov approach11,and proper orthogonal decomposition(POD)based fluid ROM.12,13

However,the aforementioned analytical models cannot ensure high accuracy when modeling true aircraft.It is essential to incorporate flight data into the model development process because this data provides the only true indication of the actual aircraft dynamics.14The system identification technique provides another solution for LPV modeling of aeroelastic system.

In the literature,typically two different approaches can be distinguished in LPV identification:global and local methods.15,16The global method enables us to obtain an LPV model in a single experiment,but it requires that the schedule sequences(flight conditions)keep varying during the experiment,which is not realistic for flight flutter test because of high risk.On the contrary,the local method consists in performing several experiments on different stabilized test points,during which the scheduling parameters remain constant.Then the resulting set of LTI models are combined into an LPV model in next interpolation step.For local method,a straightforward approach of generating LTI model from the flight data is to identify a system model using mature system identification algorithms.17–19Direct application of these methods to aeroelastic systems rarely produces an accurate reduced-order model.One of the difficulties is that the flight data are typically of poor quality because of low signal to noise ratio(SNR).Moreover,these algorithms also face another challenge brought by reduced-order model identification,which requires the simplest mathematical representation that captures the dominant behavior of system.But it is not easy to make a suitable order selection for aeroelastic system,which is in essence an infinite dimensional system.Consequently,overfitting usually leads to high-order identified models containing not only the dominant but also minor dynamics of system.On the other hand,even with right reduced-order,most identification algorithms using curve fitting cannot guarantee an accurate dominant mode estimation with such a low order.

In this context,the present work focuses on the generation of accurate reduced-order LTI models based on state space form from flight testing data.In the proposed procedure,the poly-reference least squares complex frequency-domain(p-LSCF)method20,21is utilized firstly to provide extremely clear stabilization diagrams even with high levels of noise.The physical dominant modes can be easily determined and extracted by interpreting the clear stabilized diagrams with fuzzy clustering technique.Secondly,based on the preidentified poles,an accurate reduced-order model is established by the enhanced stochastic ML estimator in frequency domain.Meanwhile,the subspace identification algorithm with preknown poles is also proposed to give an initial estimation.The new procedure enables us to derive the reduced-order model straightforwardly from flight data without model order reduction processing.It guarantees the model simplicity and fidelity simultaneously.

The layout of the paper is as follows.The problem formulation is presented in Section 2.Section 3 presents a detailed introduction of the proposed identification procedure.The continuous-time domain model is obtained in Section 4 from the discrete-time model.Section 5 illustrates the efficiency of new algorithm by using real flight flutter testing data.Finally,the conclusion of the study is presented in Section 6 by summarizing the main results.

2.Problem statement

2.1.Flight flutter tests

Current flight flutter tests typically consist of flight under different conditions of airspeed and altitude while applying some form of excitation to the structure.One of the excitations widely used in industrial aeronautics is control surface oscillating.The sine sweep signal is artificially generated and fed into the flight control system to actuate the control surface.The produced oscillatory forces will excite the structure.Such exciters allow the flight-test engineer to explore the dynamics within the interesting frequency range.The sweep signalsrepresentative forthe force inputare typically recorded and dynamic responses are usually measured by accelerometers.

The flight testing data are typically characterized by high noise level due to atmospheric turbulence.The immeasurable turbulence can be viewed as the process noise(Fig.1).U0and Y0are true sweep and its accelerometer response free noise respectively.The measured output Ymis disturbed with noise.Nmis the measurement noise at output.Ngdenotes atmospheric turbulence.The Fourier spectrum of measured output Ymcan be written as

With Eq.(1),the cross power spectrum density(PSD)of Ymand U0is represented as

where SU0(k)is the auto PSD of U0.

2.2.Model description

2.2.1.Right matrix fraction description(MFD)

where N(zk)∈ Cno×ni,D(zk) ∈ Cni×niare the numerator matrix polynomial and denominator matrix polynomial

wherezk=ejωk,nis the polynomial order;the estimated polynomial coefficients Biand Aican be assembled in the following matrices:

2.2.2.State-space model

The estimated FRF data can also be modeled by means of state-space model in frequency domain:

where A ∈ Cn×n,B ∈ Cn×ni,C ∈ Cno×n,D ∈ Cno×niare the system matrices to be estimated,X(k)∈ Cn×niis the state,and W(k)∈ Cn×niand V(k)∈ Cno×niare the system and sensor noise matrices that are circular complex normally distributed according to

where vec(·)is the vector obtained by stacking all columns of matrix on top of one another,Inianni×niidentity matrix,and?the Kronecker operator.

Compared with MFD,the reduced model in state space form is more suitable for controller synthesis.Moreover,the parameterization of noise variance will help us find more accurate estimation in stochastic framework.

3.Proposed procedure

The aeroelastic system is typical distributed parameter system(DPS).The work herein is to find a lumped model with reduced-orderto approximatetheDPS havinginfinitedimensional time/space nature.On the other hand,the corrupted testing data often introduce some unstable spurious poles into models.The physical nature of aeroelastic system decides that it should be stable when the flutter testing is implemented on subcritical points.Therefore the key challenge behind this work is to determine and select dominant modes.To tackle such challenge,a new procedure is developed for aeroelastic model identification.The main idea is to perform the identification process in two steps.In the first step,as shown in Fig.2,the well-known p-LSCF algorithm is employed for modal parameter identification and model order determination.After parameter identi fication,the clear stabilization diagrams and clustering technique enable us to extract the physical modes automatically.In the second step,a modified ML estimation is presented for highly accurate model identi fication with pre-identi fied physical poles.A variant of classical subspace identification is also proposed for providing initial estimates.

3.1.Stable physical poles extraction from noisy flight data

3.1.1.Modal parameter identification by p-LSCF

Determination of physical modes is a key step in the process of reduced-order modelidentification.Since the identified aeroelastic system is typical distributed parameter system with infinite dimensions in theory,classical model order selection criterion(e.g.AIC)is not applicable to such system.Moreover,heavy noise of flight data impairs the identification results.The mathematical poles resulting from noise usually lead to unstable identified model even at the subcritical points.For this reason,the stabilization diagram is commonly used for order selection and physical modes determination in practice.The well-known p-LSCF method is utilized to provide very clear stabilization diagram,which is helpful to make a selection of the dominant physical modes.

The linear least-squares(LS)equation error is obtained by right multiplication of Eq.(4)by D(zk)and taking the difference between the left and right part,resulting in

where W(zk)is weighted function.

The LS solution of matrix coefficients can be obtained by solving

where ?(·)is real part of a complex matrix and(·)Hconjugate transpose of matrix.The interesting problem is denominator matrix coefficients α from which poles and modal parameters can be computed.Consequently,the β coefficients are eliminated by

The so-called reduced normal equation can be obtained as follows:

This equation can be solved for the denominator coefficients matrix α in a LS sense.To obtain a clear stabilization diagram,A0=Iniare chosen.This selection principle23can ensure that the spurious modes are unstable and the physical modes are stable.For the constraint,the LS solution of α is obtained by

3.1.2.Establishment of stabilization diagrams

Establishment of a stabilization diagram needs the information of system poles from the denominator coefficients matrix α with an increasing model orders.The basic idea is that several runs of the complete pole identification process are made by using models of increasing order.Experience shows that in such analysis,the pole values of physical modes always appear at a nearly identical frequency,while spurious poles tend to scatter around the frequency range.The pole values from all these analyses at different orders can be combined in one single diagram with the natural frequency as horizontal axis and the solution order as vertical axis.In these diagrams,modes that appear in most of these models with consistent frequency,mode shape and damping are classified as stable and are likely to be physical.Modes that only appear in some models are considered spurious.

To obtain a clear stabilization diagram,most of spurious modes can be removed by using the so-called stabilization criteria.These criteria include two types:

Preliminary criteria.The frequenciesfiand damping ratios ξiare expected to be obtained in certain rangesfmin≤fi≤fmaxand ξmin≤ ξi≤ ξmax.Poles having frequenciesfiand damping ratios ξiout of these ranges will be removed.Only modes verifying these criteria are plotted in the stabilization diagram.The frequency range can be defined by pre-knowledge of interesting flutter modes,and the damping ratio ξiis allowed to vary from 0 to 0.2.

Stabilization criteria.Modes,for which the differences in modal parameters between two consecutive model orders are higher than certain threshold values,are not shown in the diagram.In the classical implementation of the stabilization diagram,the typical stability criteria values are as follows:εf=1%for frequency and εξ=5%for damping.

In fact,despite frequency and damping criteria,the modal assurance criteria(MAC)are also commonly used in stabilization diagram establishment.The experience shows that the frequency and damping stability criteria are enough to setup a clear plot from flutter testing data.Hence only the first two criteria are employed in our work.

3.1.3.Clustering of stabilized diagram

To make the order determination easier,another way for viewing stabilized diagram is used here,in which the identified stable poles are drawn in the damping ratio–frequency plane.Looking briefly at the graph shown in Fig.3 already reveals that the physical poles correspond to dense areas of data points.Inspired by this finding,fuzzy clustering has been introduced to interpret the stabilized diagram.The goal of cluster is to group the estimated modes that represent the same physical mode.

A variety of fuzzy clustering techniques have been introduced into this field,e.g.fuzzy C-means(FCM)clustering algorithm.24The clusters related to physical modes in Fig.3 show ‘vertical columned’features rather than spherical form.The research in Ref.25has demonstrated that the extended version of FCM,namely Gustafson–Kessel(GK)algorithm,has great advantages in such hyper ellipsoidal clustering.This algorithm is employed herein,and the clustering procedure works for grouping physical modes as follows:

(1)Normalization of the data sets.The magnitude of damping ratios is much less than frequencies.A direct application of geometric measures(distances)will implicitly assign bigger contributions to the frequencies than the damping.Therefore,the attributes should be dimensionless to reduce the effects of units of measurements on results of clustering.The values used herein are normalized between[0,1].

(2)Clustering of data set using GK algorithm.Before the

clustering,the number of the clusters must be given by the user,but it is rarely known a priori,in this case it must be searched by using validity measures.Several scalar validity measures including Xie and Beni’s Index(XB),26Dunn’s Index(DI)27and Alternative Dunn’s Index(ADI)index are used to help us find the optimal number of clusters.Additionally,a modified Gustafson-Kessel algorithm28is utilized to group the data set.

(3)Cluster selection.After clustering,a decision should be made whether this cluster is selected to represent a physical pole based on the number and scattering of members of clusters.The contribution ratio is a measure of the number,and it is defined based on the number of members of each cluster as a percentage of the total number of points in the data set.The compactness is an index of the scatter of the members of the cluster,which is derived for each cluster by calculating the largest Euclidian distance of all the members to the cluster centers.

(4)Physical poles selection.By using above clustering analysis method,the cluster centers suggested by the fuzzy GK algorithms are selected as data points corresponding to physical poles.A pair of complex conjugate poles are given by

where ωi=2πfi,fiand ξiare the corresponding frequency and damping ratio of theith cluster center,and j denotes the imaginary number with j2=-1.

3.2.Frequency-domain subspace identi fication with known poles(SSKP)

With the selected physical poles,the next work is to identify the reduced-order state-space model(zeros identification).Subspace identification is a traditional method for statespace model identification. Classical frquency-domain subspace identification from FRF data just uses singular value decomposition to determine model order and provide an estimate of(A,B,C,D)without consideration of the poles constraint.Additionally,the estimated model derived by classical algorithm may contain unstable poles because of noisy data.This section presents a frequency-domain subspace algorithm with known poles(SSKP)to improve original algorithm.

3.2.1.Brief review of classical subspace algorithm29

Recursive use of Eq.(9)gives

Writing Eq.(22)forp=1，2，···，r-1,and stacking them on top of each other give

Collecting Eq.(23)fork=1，2，···，Mwith the numberMof discrete frequency lines in the interesting frequency band gives

Eq.(25)is converted into a real set of equations

where(·)re=[?(·)，?(·)]. ?(·)and ?(·)represent real and imaginary parts of a complex matrix.

Next,Π⊥denotes a matrix which projects onto the null space of UreM.Since UreMΠ⊥=0,GreMΠ⊥=OrXreMΠ⊥can be directly obtained.In practice,this step is performed by QR-factorization,followed by singular value decomposition.

Frequency-domain subspace algorithm uses singular value decomposition to judge model’s order.Assume that Σscontainsnlargest singular values,and then model order isn.

The extended observability matrix Oris estimated as follows:

The estimates of A and C from Orin a least-squares sense are given by

3.2.2.Proposed subspace identification

It is reasonable to assume that a real physical system has distinct poles.Hence,it is assumed that coefficient matrix A of state-space model has distinct eigenvalues λi(i=1，2，···，n)with corresponding linear independent eigenvector pi.

For real eigenvalue λi,from Section 3.1,there is eigenvalue equation

For aeroelastic system,complex conjugate eigenvalues pair λi,λHare commonly used.In this case,eigenvalue equation is

iwritten as

With Eqs.(36)and (37),matrix A hasa similarity transformation

whereΛ0=diag(λ1，λ2，···，λn)andP=[p1，p2，···，pn].

To estimate matrix A with fixed eigenvalues,it is clearly necessary to derive the similarity transformation matrix P.This may be simply achieved via the results of the following Lemma.

Lemma 1.For a state-space model(A,B,C,D)with fixed poles,ifλiis one of eigenvalues ofAandpiis the eigenvector corresponding toλi,then the eigenvector^piof matrixT-1ATcorresponding toλican be found by solving the following equation:

whereTis non-singular matrix,andT=XreMΠ⊥Q2VsΣ-1s.

Proof.As mentioned in Ref.29,the estimated^Orin Eq.(32)only represents the range space of Or. If we set T=XreMΠ⊥Q2VsΣ-1s,there is^Or=OrT.

Thus,combining Eqs.(24)and(32)yields

3.3.EM algorithm based ML estimation with known poles(MLKP)

Maximum likelihood(ML)estimation enjoys reputation in the field of system identi fication for its high accuracy.The main dif ficulty in identifying a system in state-space form is that the likelihood function is non-convex.The ML estimation based on EM algorithm hasshown greatadvantage over Gauss-Newton type search when one tackles such problem.30,31

3.3.1.Background of EM algorithm

In the ML framework,the following log-likelihood function is maximized:

The basic idea of the EM algorithm is to decompose the log-likelihood function as

By using the Jensen inequality,it is possible to proveH(θ，θi) ≤H(θi，θi).32Therefore,any value of θ for whichQ(θ，θi)≥Q(θi，θi)impliesl(θ) ≥l(θi).This suggests a strategy of maximizingQ(θ，θi),which increasesl(θ),and then setting θi+1equal to the maximum and repeating the process.That is well known as EM algorithm with following two steps:

(1)Estep

(2)Mstep

3.3.2.E step

To implement the EM algorithm,the first step is to derive a method for calculation ofQ(θ，^θi).This may be achieved via the results of the following Lemma.

3.3.3.M step

4.Continuous-time model identification

Although the discrete-time model leads to better condition than continuous-time model since powers of ejωk=ej2πk/Mform a natural orthogonal basis,the direct discrete-time model identification in narrow frequency range still suffers from illconditioned problem.For this reason,an alternative approach used in current work is to represent the aeroelastic system based on continuous-time model firstly,and then transform the continuous-time problem into discrete-time domain by using bilinear transformation,which is defined as

This nonlinear mapping enables us to perform wellconditioned system identification based on discrete-time model.Once the discrete-time model(A，B，C，D)is identified at the scaled frequency point,the continuous-time model can be derived later by use of the inverse mapping29:

5.Application to real flight flutter tests

This section intends to use real flight test to verify the effectiveness of the proposed approach.A large aircraft was excited by the outer aileron control signals with a sine sweep signal added in.The sine sweep generated signals are chosen as input,while the measured acceleration responses at various locations of the airplane:fuselage,engines and wings are chosen as outputs.Data from 4 accelerometers are analyzed.The whole test system can be regarded as a single input multiple output system(1 input and 4 outputs).The sample frequency for data acquisition is 256 Hz,and the data length is 4608.It needs to be noted that the interested modal parameters are distributed in a narrow frequency band(for confidentiality reasons,the specific frequency coordinates are not given in the figures of this section except Fig.4 with normalized frequency).

Table 1 Numerical value of validty measures.

Firstly,the flight testing data are handled with a normalized processing,and the FRF data are estimated by local polynomial method(LPM).33Secondly,the singular value decomposition(SVD)of subspace identification algorithm is used to determine model order.As shown in Fig.5,gap between the singular values indicates that the order of this system isn=4.

Meanwhile,a stabilized diagram is also drawn in Fig.3 to determine the correct model order.The p-LSCF algorithm is employed for model orders decreasing fromn=35 ton=10.From the stabilized diagram,we can see at least three vertical alignments of stable poles,which represent the first 3 modes of the aircraft in range of the narrow frequency band(the fourth mode in higher frequency is not clear).This implies that the system order is equal to or larger than 6,which is incompatible with the SVD result in Fig.5.Obviously,if we use singular values to determine the order,the model order will be set to a lower value,and at least one of these modes is not identified.

However,it needs to be noted that the fourth mode in higherfrequency isobscure,and thereisno obvious corresponding peak of FRF curve.To determine whether it represents a physical mode or not,all the stable damping poles are plotted in damping ratio-frequency plane(Fig.4).

Before cluster analysis,the optimal cluster number should be decided.Validity measures including XB,DI and ADI index are used with Gustafson-Kessel algorithm to help us find the optimal number of clusters.The values of the validity measures depending from the number of cluster are embraced in Table 1.The XB index reaches the local minimum at number=6.However,the DI and ADI indices are maximized at number=4.To determine the optimal number clearly,the clustering results for number=4 and 6 are shown in normalized frequency-damping ratio plane(Fig.4).It can be seen that some members originally involved in clusters C1 and C2 are partitioned into new clusters C5 and C6 when the number increases from 4 to 6.The other clusters like C3 and C4 are unchanged.

To facilitate optimal number decision and physical poles extraction,the contribution ratio and compactness of each cluster are given in Table 2.The compactness is also graphically visualized in Fig.4 by a radius of a circle around each cluster.The smaller circle implies the better compactness of the cluster.As illustrated in Table 2,the new emerging clusters C5 and C6 have low contribution ratios for a few members included.Compared with the clusters C5 and C6,the clusters C1 and C2 remain the dominant number of cluster member.They have great contribution ratios so that the corresponding modes are distinguished as physical modes.Additionally,if carefully comparing the centers of clusters C1 and C2 in Fig.4(a)and(b),we can find that the centers are almost in the same positions,and the differences are minor.This indicates that the partition with less clusters has little influence on physical poles estimation.Therefore,in this example,the optimal number of cluster for FCM GK algorithm is 4.The cluster C3 is selected as representing a physical modes for large contribution ratio and good compactness,while the cluster C4 represents the spurious poles for little contribution ratio and bad compactness.

From the above analysis,three clusters including C1,C2 and C3 are selected as representing the physical modes so that the model order is determined to be 6.The poles of the system can be calculated by using Eq.(21)firstly,and then the proposed subspace algorithm(SSKP)and ML algorithm(MLKP)with fixed poles are used to identify the state-space model of aeroelastic system.The Figs.6–9 show the estimated FRF using two algorithms comparing with non-parametric FRF estimation by LPM from flight testing data.

Table 2 Numerical value of validity measures.

Comparing the results illustrated in Figs.6–9,we can find that both two algorithms provide good estimate.The results show that the reduced-order model can fully represent the three dominant modes.Meanwhile,the second algorithm based on ML yields better results because the model error is smaller.This indicates that the proposed procedure can ensure the model’s simplicity and fidelity simultaneously.

6.Conclusions

This paper presents a new reduced-order modeling approach for aeroelastic system by identifying it from flight test data.The p-LSCF algorithm and clustering algorithm are employed to extract the physical dominant modes firstly.Then two enhanced frequency algorithms are developed to ensure accurate reduced-order model identification based on the pre-known poles.The real application results illustrate its eff iciency.Future work will focus on the establishment of LPV model of the aircraft aeroelastic system with a set of reduced-order LTI models derived by the proposed approach.

Acknowledgements

This study was co-supported by the National Natural Science Foundation of China(Nos.61134004 and 61573289),Aeronautical Science Foundation of China(No.20140753010)and the FundamentalResearch Funds forthe Central Universities(No.3102015BJ004).

1.Mukhopadhyay V.Historical perspective on analysis and control of aeroelastic responses.J Guid Control Dynam2003;26(5):673–84.

2.Kekoe MW.A historical preview offlight fluttertesting.Washington D.C.:Dryden Flight Research Center;1995.Report No:4702.

3.Barker JM,Balas GJ.Comparing linear parameter-varying gainscheduled control techniques for active flutter suppression.J Guid Control Dynam2000;23(5):948–55.

4.Prime Z,Cazzolato B,Doolan C,Strganac T.Linear-parametervarying control of an improved three-degree-of-freedom aeroelastic model.J Guid Control Dynam2010;33(2):615–9.

5.Seiler P,Balas GJ.Packard A Linear parameter-varying control for the X-53 active aeroelastic wing.In:Mohammadpour J,Scherer WC,editors.Control of linear parameter varying systems with applications.New York:Springer Science+Business Media;2012.p.483–512.

6.Poussot-Vassal C,Roos C.Generation of a reduced-order LPV/LFT model from a set of large-scale MIMO LTI flexible aircraft models.Control Eng Pract2012;20(9):919–30.

7.Visser M,Navalkar ST,Wingerden VJW.LPV model identification for flutter prediction:A comparison of methods.IFAC PapersOnline2015;48(26):121–6.

8.Moreno CP,Seiler P,Balas G.Model reduction for aeroservoelastic systems.J Aircraft2014;51(1):280–90.

9.Lall S,Beck C.Error-bounds for balanced model-reduction of linear time-varying systems.IEEE Trans Autom Control2003;48(6):946–56.

10.Gugercin S,Antoulas AC,Beattie C.H2model reduction for large-scale linear dynamical systems.SIAM J Matrix Anal Appl2008;30:609–38.

11.Gugercin S.An iterative SVD-Krylov based method for model reduction of large-sale dynamical systems.Linear Algebra Appl2008;428(8–9):1964–86.

12.Thomas JP,Dowell EH,Hall KC.Three-dimensional transonic aeroelasticity using proper orthogonaldecomposition-based reduced-order models.J Aircraft2003;40(3):544–51.

13.Zhang W,Wang BB,Ye ZY,Quan JG.Efficient method for limit cycle flutter analysis based on nonlinear aerodynamic reducedorder models.AIAA J2012;50(5):1019–28.

14.Lind R,Brenner M.Incorporating flight data into a robust aeroelastic model.J Aircraft1998;35(3):470–7.

15.Lovera M,Bergamasco M,Casella F.LPV modeling and identification:an overview.In:Sename O,Ga′spa′r P,Bokor J,editors.Robustcontrolandlinearparametervarying approaches.Berlin:Springer;2013.p.3–24.

16.Wingerden JWV.Control of wind turbines with ’Smart’rotors:proof of concept&LPV subspace identification[dissertation].Delft:Delft University of Technology;2008.

17.Verboven P,Cauberghe B,Guillaume P,Vanlanduit S,Parloo E.Modal parameter estimation and monitoring for on-line flight flutter analysis.Mech Sys Sig Proc2004;18(3):587–610.

18.Miller DN,De Callafon RA,Brenner MJ.Covariance-based realization algorithm for the identification of aeroelastic dynamics.J Guid Control Dynam2012;35:1169–77.

19.Tang W,Shi ZK,Chen J.Aircraft flutter modal parameter identification using a numerically robust least-squares estimation frequency domain.Chin J Aeronaut2008;21(6):550–8.

20.Peeters B,Auweraer VD,Guillaume P,Schoukens J.The PolyMAX frequency-domain method:A new standard for modal parameter estimation.Shock Vib2004;11(3–4):395–409.

21.Cauberghe B.Applied frequency-domain system identification in the field ofexperimentaland operationalmodalanalysis[dissertation].Brussel:Vrije University;2004.

22.Pintelon R,Schoukens J.System identification:a frequency domain approach.1st ed.Piscataway:IEEE Press;2001.p.155.

23.Tang W,Qiao Q,Shi ZK.On the choice of parameter constraint for frequency–domain least squares identification.Acta Aeronaut Astronaut Sin2012;33(12):2253–60[Chinese].

24.Bezdek JC.Pattern recognition with fuzzy objective function.New York:Plenum Press;1981.

25.Scionti M,Lanslots JP.Stabilisation diagrams:pole identification using fuzzy clustering techniques.Adv Eng Softw2005;36:768–79.

26.Xie XL,Beni GA.Validity measure for fuzzy clustering.IEEE Trans Patt Analy Mach Intell1991;3(8):841–6.

27.Dunn JC.A fuzzy relative of the isodata process and its use in detecting compact well-separated clusters.J Cybernet1973;3(3):32–57.

28.Babuska R,VanderVeen VJ,Kaymak U.Improved co-variance estimation for Gustafson Kessel clustering.Proceedings of 2002 IEEE international conference on fuzzy systems;2002 May 12–17.Honolulu,HI.Piscataway(NJ):IEEE Press;2002.p.1081–5.

29.McKelvey T,Akc?ay H,Ljung L.Subspace-based multivariable system identification from frequency response data.IEEE Trans Autom Control1996;41(7):960–79.

30.Agu¨ero JC,Tang W,Yuz IJ,Delgado R,Goodwin GC.Dual time–frequency domain system identification.Automatica2012;48(12):3031–41.

31.Wills AG,Ninness B,Gibson SH.Maximum likelihood estimation of state space models from frequency domain data.IEEE Trans Autom Control2009;54(1):19–33.

32.Dempster AP,Laird NM,Rubin DB.Maximum likelihood from incomplete data via the EM algorithm.J R Stat Soc B(Meth)1977;39(1):1–38.

33.Ugryumova D,Pintelon R,Vandersteen G.Frequency response function estimation in the presence of missing output data.IEEE Trans Instrum Meas2015;64(2):541–53.

7 August 2016;revised 9 November 2016;accepted 14 November 2016

Available online 27 December 2016

?2016 Chinese Society of Aeronautics and Astronautics.Production and hosting by Elsevier Ltd.This is anopenaccessarticleundertheCCBY-NC-NDlicense(http://creativecommons.org/licenses/by-nc-nd/4.0/).

*Corresponding author.

E-mail address:tangwei@nwpu.edu.cn(W.Tang).

Peer review under responsibility of Editorial Committee of CJA.