**Adaptive and Fault Tolerant Flight Control**

92 Automatic Flight Control Systems – Latest Developments

Russo, G. (2009). DTFT-1: Analysis of the first USV flight test, *Acta Astronautica*, Volume 65,

Ryali, V., Moudgalya, K.M. (2005). Practical stability analysis of uncertain nonlinear

Tancredi, U., Grassi, M., Corraro, F., Filippone, E., and Russo, M. (2008). A Novel Approach

Tancredi, U., Grassi, M., Corraro, F., and Filippone E. (2009). Robustness Analysis for

Tancredi, U., Grassi, M., Corraro, F., Vitale, A., & Filippone, E. (2011). A linear time varying

Van der Merwe, R., De Freitas, N., Doucet, A., and Wan, E. A. (2000). The Unscented Particle Filter, *Advances in Neural Information Processing Systems*, pp.584–590. Wang, Q., and Stengel, R. F. (2002). Robust control of nonlinear systems with parametric

*Control in Aerospace*, Vol. 1, No. 1, Paper 2, Retrieved from:

uncertainty, *Automatica*, Vol. 38, No. 9, pp. 1591-1599.

systems, *Proceedings of the National Conference on Control and Dynamic Systems*, I.I.T.

to Clearance of Flight Control Laws over Time Varying Trajectories", *Automatic*

Terminal Phases of Reentry Flight, *AIAA Journal of Guidance, Control, and Dynamics*,

approach for robustness analyses of a re-entry flight technology demonstrator, *Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace*

Issues 9-10, Pages 1196-1207, ISSN 0094-5765.

<http://www.aerospace.unibo.it/index.php?e=5>

Vol. 32, No. 5, pp. 1679 - 1683.

*Engineering* , in press.

Bombay.

**1. Introduction**

Active control techniques for the gust loads alleviation/flutter suppression have been investigated extensively in the last decades to control the aeroelastic response, and improve the handling qualities of the aircraft. Nonadaptive feedback control algorithms such as classical single input single output techniques (Schmidt & Chen, 1986), linear quadratic regulator (LQR) theory (Mahesh et al., 1981; Newsom, 1979), eigenspace techniques (Garrard & Liebst, 1985; Leibst et al., 1988), optimal control algorithm (Woods-Vedeler et al., 1995), *H*∞ robust control synthesis technique (Barker et al., 1999) are efficient methods for the gust loads alleviation/flutter suppression. However, because of the time varying characteristics of the aircraft dynamics due to the varying configurations and operational parameters, such as fuel consumption, air density, velocity, air turbulence, it is difficult to synthesize a unique control law to work effectively throughout the whole flight envelope. Therefore, a gain scheduling

**Adaptive Feedforward Control for**

Jie Zeng1, Raymond De Callafon2 and Martin J. Brenner3

**Gust Loads Alleviation**

1*ZONA Technology, Inc.* 2*University of California, San Diego*

3*NASA Dryden Flight Research Center*

**0**

**4**

*USA*

An alternative methodology is the feedforward and/or feedback adaptive control algorithms by which the control law can be updated at every time step (Andrighettoni & Mantegazza, 1998; Eversman & Roy, 1996; Wildschek et al., 2006). With the novel development of the airborne LIght Detection and Ranging (LIDAR) turbulence sensor available for an accurate vertical gust velocity measurement at a considerable distance ahead of the aircraft (Schmitt, Pistner, Zeller, Diehl & Navé, 2007), it becomes feasible to design an adaptive feedforward control to alleviate the structural loads induced by any turbulence and extend the life of the structure. The adaptive feedforward control algorithm developed in (Wildschek et al., 2006) showed promising results for vibration suppression of the first wing bending mode. However, an unavoidable constraint for the application of this methodology is the usage of a high order Finite Impulse Response (FIR) filter. As a result, an overwhelming computation effort was

In this chapter, an adaptive feedforward control algorithm where the feedforward filter is parameterized using orthonormal basis expansions along with a recursive least square algorithm with a variable forgetting factor is proposed for the feedforward compensation of gust loads. With the use of the orthonormal basis expansion, the prior flexible modes information of the aircraft dynamics can be incorporated to build the structure of the feedforward controller. With this strategy, the order of the feedforward filter to be estimated

technique is necessary to account for the time varying aircraft dynamics.

needed to suppress the structural vibration of the aircraft.

## **Adaptive Feedforward Control for Gust Loads Alleviation**

Jie Zeng1, Raymond De Callafon2 and Martin J. Brenner3 1*ZONA Technology, Inc.* 2*University of California, San Diego* 3*NASA Dryden Flight Research Center USA*

## **1. Introduction**

Active control techniques for the gust loads alleviation/flutter suppression have been investigated extensively in the last decades to control the aeroelastic response, and improve the handling qualities of the aircraft. Nonadaptive feedback control algorithms such as classical single input single output techniques (Schmidt & Chen, 1986), linear quadratic regulator (LQR) theory (Mahesh et al., 1981; Newsom, 1979), eigenspace techniques (Garrard & Liebst, 1985; Leibst et al., 1988), optimal control algorithm (Woods-Vedeler et al., 1995), *H*∞ robust control synthesis technique (Barker et al., 1999) are efficient methods for the gust loads alleviation/flutter suppression. However, because of the time varying characteristics of the aircraft dynamics due to the varying configurations and operational parameters, such as fuel consumption, air density, velocity, air turbulence, it is difficult to synthesize a unique control law to work effectively throughout the whole flight envelope. Therefore, a gain scheduling technique is necessary to account for the time varying aircraft dynamics.

An alternative methodology is the feedforward and/or feedback adaptive control algorithms by which the control law can be updated at every time step (Andrighettoni & Mantegazza, 1998; Eversman & Roy, 1996; Wildschek et al., 2006). With the novel development of the airborne LIght Detection and Ranging (LIDAR) turbulence sensor available for an accurate vertical gust velocity measurement at a considerable distance ahead of the aircraft (Schmitt, Pistner, Zeller, Diehl & Navé, 2007), it becomes feasible to design an adaptive feedforward control to alleviate the structural loads induced by any turbulence and extend the life of the structure. The adaptive feedforward control algorithm developed in (Wildschek et al., 2006) showed promising results for vibration suppression of the first wing bending mode. However, an unavoidable constraint for the application of this methodology is the usage of a high order Finite Impulse Response (FIR) filter. As a result, an overwhelming computation effort was needed to suppress the structural vibration of the aircraft.

In this chapter, an adaptive feedforward control algorithm where the feedforward filter is parameterized using orthonormal basis expansions along with a recursive least square algorithm with a variable forgetting factor is proposed for the feedforward compensation of gust loads. With the use of the orthonormal basis expansion, the prior flexible modes information of the aircraft dynamics can be incorporated to build the structure of the feedforward controller. With this strategy, the order of the feedforward filter to be estimated

Gust Loads Alleviation 3

Adaptive Feedforward Control for Gust Loads Alleviation 97

*w*(*t*) *d*(*t*) *+ e*(*t*)

*F G*

Fig. 1. Block Diagram of the Structural Vibration Control with Feedforward Compensation.

• A feedforward-looking measuring of 50 to 150 m to ensure that the measured air flow is

• The standard deviation of the wind speed measurement should be small, at least in the

A sensor system that meets these requirements is a so-called short pulse UV Doppler LIDAR, and was developed in (Schmitt, Pistner, Zeller, Diehl & Navé, 2007). This short pulse UV Doppler LIDAR was successfully applied to an Airbus 340 to measure the vertical gust speed (Schmitt, Rehm, Pistner, Diehl, Jenaro-Rabadan, Mirand & Reymond, 2007). The authors in (Schmitt, Rehm, Pistner, Diehl, Jenaro-Rabadan, Mirand & Reymond, 2007) claimed that the

Assuming a perfect gust perturbation signal can be measured via the LIDAR beam sensor,

In case the transfer functions in Eq. (1) are known, an ideal feedforward controller, *F*(*q*) =

*Fi*(*q*) = − *H*(*q*)

in case, *Fi*(*q*), is a stable and causal transfer function. The solution of *Fi*(*q*) in Eq. (2) assumes full knowledge of *G*(*q*) and *H*(*q*). Moreover, the filter, *Fi*(*q*), may not be a causal or stable filter due to the dynamics of *G*(*q*) and *H*(*q*) that dictate the solution of the feedforward controller, *Fi*(*q*). An approximation of the feedforward filter, *Fi*(*q*), can be made by an output-error based optimization that aims at finding the best causal and stable approximation, *F*(*q*), of the ideal

A direct adaptation of the feedforward controller *F*(*q*, *θ*) can be performed by considering the

• The sensor must be able to produce reliable signals in the absence of aerosols. • A good longitudinal resolution (the thickness of the air slice measured ahead).

system is ready to be used to design feedforward control for gust loads alleviation.

*u*(*t*) *y*(*t*)

*+*

*e*(*t*) = [*H*(*q*) + *G*(*q*)*F*(*q*)] *w*(*t*) (1)

*e*(*t*, *θ*) = *H*(*q*)*w*(*t*) + *F*(*q*, *θ*)*G*(*q*)*w*(*t*). (3)

*d*(*t*) := *H*(*q*)*w*(*t*), *u f*(*t*) := −*G*(*q*)*w*(*t*) (4)

*G*(*q*) (2)

*H*

*Gust perturbation.*

*Lidar Beam*

the one actually affecting the aerodynamics around the aircraft.

• The sensor must be able to measure the vertical wind speed.

that means, *n*(*t*) = *w*(*t*), the error signal, *e*(*t*), can be described by

*n*(*t*)

range of [2-4] m/s.

*Fi*(*q*), can be obtained by

feedforward controller in *Fi*(*q*) in Eq. (2).

parameterized error signal, *e*(*t*, *θ*),

Defining the signals

• A good temporal resolution.

can be largely reduced. As a result, the computation effort is greatly decreased, and the performance of the feedforward controller for gust loads alleviation will be enhanced. Furthermore, an FFT based PolyMAX identification method and the stabilization diagram program (Baldelli et al., 2009) are proposed to estimate the flexible modes of the aircraft dynamics.

The need for an integrated model of flight dynamics and aeroelasticity is brought about by the emerging design requirements for slender, more flexible and/or sizable aircraft such as the Oblique Flying Wing (OFW), HALE, Sensorcraft and morphing vehicles, etc. Furthermore, a desirable unified nonlinear simulator should be formulated in principle by using commonly agreeable terms from both the flight dynamics and aeroelasticity fields in a consistent manner.

A unified integration framework that blends flight dynamics and aeroelastic modeling approaches with wind-tunnel or flight-test data derived aerodynamic models has been developed in (Baldelli & Zeng, 2007). This framework considers innovative model updating techniques to upgrade the aerodynamic model with data coming from CFD/wind-tunnel tests for a rigid configuration or data estimated from actual flight tests when flexible configurations are considered.

Closely following the unified integration framework developed in (Baldelli & Zeng, 2007), an F/A-18 Active Aeroelastic Wing (AAW) aeroelastic model with gust perturbation is developed in this chapter, and this F/A-18 AAW aeroelastic model can be implemented as a test-bed for flight control system evaluation and/or feedback/feedforward controller design for gust loads alleviation/flutter suppression of the flexible aircraft.

The outline of the chapter is as follows. In Section II, a feedforward compensation framework is introduced. Section III presents the formulation of the orthonormal finite impulse filter structure. A brief description of a frequency domain PolyMAX identification method is presented in Section IV. In Section V, a recursive least square estimation method with variable forgetting factor is discussed. Section VI includes the development of a linear F/A-18AAW aeroelastic model and the application of the adaptive feedforward controller to F/A-18 AAW aeroelastic model.

## **2. Basic idea of the feedforward controller**

In order to analyze the design of the feedforward controller, *F*, consider the simplified block diagram of structural vibration control of a single input signal output (SISO) dynamic system depicted in Fig. 1. The gust perturbation, *w*(*t*), passes through the primary path, *H*, the body of the aircraft, to cause the structural vibrations. Mathematically, *H* can be characterized as the model/transfer function from the gust perturbation to the accelerometer sensor position. The gust perturbation, *w*(*t*), can be measured by the coherent LIDAR beam airborne wind sensor. The measured signal, *n*(*t*), is fed into the feedforward controller, *F*, to calculate the control surface demand, *u*(*t*), for vibration compensations. The structural vibrations are measured by the accelerometers providing the error signal, *e*(*t*). *G* is the model/transfer function from the corresponding control surface to the accelerometer sensor position, and which is so called the secondary path.

In order to apply the feedforward control algorithm for gust loads alleviation, developing a proper sensor to accurately measure the gust perturbation is crucial for the success of the feedforward control application. As mentioned in (Schmitt, Pistner, Zeller, Diehl & Navé, 2007), such a sensor should meet the following criteria:

2 Will-be-set-by-IN-TECH

can be largely reduced. As a result, the computation effort is greatly decreased, and the performance of the feedforward controller for gust loads alleviation will be enhanced. Furthermore, an FFT based PolyMAX identification method and the stabilization diagram program (Baldelli et al., 2009) are proposed to estimate the flexible modes of the aircraft

The need for an integrated model of flight dynamics and aeroelasticity is brought about by the emerging design requirements for slender, more flexible and/or sizable aircraft such as the Oblique Flying Wing (OFW), HALE, Sensorcraft and morphing vehicles, etc. Furthermore, a desirable unified nonlinear simulator should be formulated in principle by using commonly agreeable terms from both the flight dynamics and aeroelasticity fields in a consistent manner. A unified integration framework that blends flight dynamics and aeroelastic modeling approaches with wind-tunnel or flight-test data derived aerodynamic models has been developed in (Baldelli & Zeng, 2007). This framework considers innovative model updating techniques to upgrade the aerodynamic model with data coming from CFD/wind-tunnel tests for a rigid configuration or data estimated from actual flight tests when flexible configurations

Closely following the unified integration framework developed in (Baldelli & Zeng, 2007), an F/A-18 Active Aeroelastic Wing (AAW) aeroelastic model with gust perturbation is developed in this chapter, and this F/A-18 AAW aeroelastic model can be implemented as a test-bed for flight control system evaluation and/or feedback/feedforward controller design for gust

The outline of the chapter is as follows. In Section II, a feedforward compensation framework is introduced. Section III presents the formulation of the orthonormal finite impulse filter structure. A brief description of a frequency domain PolyMAX identification method is presented in Section IV. In Section V, a recursive least square estimation method with variable forgetting factor is discussed. Section VI includes the development of a linear F/A-18AAW aeroelastic model and the application of the adaptive feedforward controller to F/A-18 AAW

In order to analyze the design of the feedforward controller, *F*, consider the simplified block diagram of structural vibration control of a single input signal output (SISO) dynamic system depicted in Fig. 1. The gust perturbation, *w*(*t*), passes through the primary path, *H*, the body of the aircraft, to cause the structural vibrations. Mathematically, *H* can be characterized as the model/transfer function from the gust perturbation to the accelerometer sensor position. The gust perturbation, *w*(*t*), can be measured by the coherent LIDAR beam airborne wind sensor. The measured signal, *n*(*t*), is fed into the feedforward controller, *F*, to calculate the control surface demand, *u*(*t*), for vibration compensations. The structural vibrations are measured by the accelerometers providing the error signal, *e*(*t*). *G* is the model/transfer function from the corresponding control surface to the accelerometer sensor position, and which is so called

In order to apply the feedforward control algorithm for gust loads alleviation, developing a proper sensor to accurately measure the gust perturbation is crucial for the success of the feedforward control application. As mentioned in (Schmitt, Pistner, Zeller, Diehl & Navé,

loads alleviation/flutter suppression of the flexible aircraft.

**2. Basic idea of the feedforward controller**

2007), such a sensor should meet the following criteria:

dynamics.

are considered.

aeroelastic model.

the secondary path.

Fig. 1. Block Diagram of the Structural Vibration Control with Feedforward Compensation.

A sensor system that meets these requirements is a so-called short pulse UV Doppler LIDAR, and was developed in (Schmitt, Pistner, Zeller, Diehl & Navé, 2007). This short pulse UV Doppler LIDAR was successfully applied to an Airbus 340 to measure the vertical gust speed (Schmitt, Rehm, Pistner, Diehl, Jenaro-Rabadan, Mirand & Reymond, 2007). The authors in (Schmitt, Rehm, Pistner, Diehl, Jenaro-Rabadan, Mirand & Reymond, 2007) claimed that the system is ready to be used to design feedforward control for gust loads alleviation.

Assuming a perfect gust perturbation signal can be measured via the LIDAR beam sensor, that means, *n*(*t*) = *w*(*t*), the error signal, *e*(*t*), can be described by

$$e(t) = \left[H(q) + G(q)F(q)\right]w(t)\tag{1}$$

In case the transfer functions in Eq. (1) are known, an ideal feedforward controller, *F*(*q*) = *Fi*(*q*), can be obtained by

$$F\_i(q) = -\frac{H(q)}{G(q)}\tag{2}$$

in case, *Fi*(*q*), is a stable and causal transfer function. The solution of *Fi*(*q*) in Eq. (2) assumes full knowledge of *G*(*q*) and *H*(*q*). Moreover, the filter, *Fi*(*q*), may not be a causal or stable filter due to the dynamics of *G*(*q*) and *H*(*q*) that dictate the solution of the feedforward controller, *Fi*(*q*). An approximation of the feedforward filter, *Fi*(*q*), can be made by an output-error based optimization that aims at finding the best causal and stable approximation, *F*(*q*), of the ideal feedforward controller in *Fi*(*q*) in Eq. (2).

A direct adaptation of the feedforward controller *F*(*q*, *θ*) can be performed by considering the parameterized error signal, *e*(*t*, *θ*),

$$e(t, \theta) = H(q)w(t) + F(q, \theta)G(q)w(t). \tag{3}$$

Defining the signals

$$d(t) := H(q)w(t),\ u\_f(t) := -G(q)w(t) \tag{4}$$

Gust Loads Alleviation 5

Adaptive Feedforward Control for Gust Loads Alleviation 99

*L*−1 ∑ *k*=0

where *q*−1 denotes the usual time shift operator, *q*−1*x*(*t*) = *x*(*t* − 1). Adaptive filter estimation using FIR filters converges to optimal and unbiased estimates irrespective of the coloring of the noise on the output data. However, a FIR filter is usually too simple to model complex system dynamics such as AE/ASE systems with many resonance modes being excited by atmospheric perturbations. As a result, many tapped delay coefficients of the FIR filter are required to approximate the optimal filter. Even though an IIR filter is appealing as an alternative, the inherent stability and bias estimation problems limit the use of an IIR filter

To improve the approximation properties of the adaptive filter, *F*, in Fig. 1, the linear combination of tapped delay functions, *q*−1, in the FIR filter of Eq. (8) can be generalized

> *L*−1 ∑ *k*=0

where *Bk*(*q*) are generalized (orthonormal) basis functions (Heuberger et al., 1995) that contain some *a-priori* knowledge on the desired filter dynamics. In other words, the orthonormal basis functions that are used in the parametrization of the ORTFIR filter will be tuned on the fly by

The application of orthonormal basis functions to parameterize and estimate dynamical systems has obtained extensive attention in recent years. Different constructions of the orthonormal basis structure has been reported in (Heuberger et al., 1995; Ninness & Gustafsson, 1997; Zeng & de Callafon, 2006). It is assumed that the pole locations are already known with the use of the standard open-loop prediction error system identification methods. Suppose the poles {*ξi*}*i*=1,2,···,*N* are known, an all pass function, *P*(*q*), can be created by these

> *N* ∏ *i*=1

Let (*A*, *B*, *C*, *D*) be a minimal balanced realization of an all pass function, *P*(*q*), define the input to state transfer function, *B*0(*q*)=(*q I* − *A*)−1*B*, then a set of functions, *Bi*(*q*), can be

*Bi*(*q*) = *B*0(*q*)*Pi*

*k* (1/*q*)

The construction of the orthonormal basis function is illustrated in Fig. 3. It should be noted that if *B*0(*q*), includes all information of a dynamical system, then only one parameter, *β*0, needs to be estimated to approximate this dynamic system. It means that the parameters

d*q q* =

 1 − *ξ*∗ *i q q* − *ξi*

*I i* = *k*

*F*(*q*, *θ*) =

taking full advantage of the modal information embedded in the flight data.

*P*(*q*) =

*Bi*(*q*)*BT*

1 2*πj*  *βkq*−*k* (8)

*βkBk*(*q*) (9)

(*q*) (11)

0 *i* �= *k* (12)

(10)

Generally, the discrete time linear time invariant (LTI) FIR filter, *F*(*q*), can be presented as:

*F*(*q*) =

for adaptive filtering in aeroservoelastic systems.

Construction of the Orthonormal Basis Sets

to the following form:

poles, and is given as

and *Bi*(*q*) has orthogonal property

obtained via

where *d*(*t*) can actually be measured, and *u f*(*t*) is called filtered input signal, Eq. (3) is reduced to

$$e(t, \theta) = d(t) - F(q\_\prime \theta) \mu\_f(t) \tag{5}$$

for which the minimization

min *θ* 1 *N N* ∑ *t*=1 *e* 2(*t*, *θ*) (6)

to compute the optimal feedforward filter, *F*(*q*, *θ*), is a standard output-error (OE) minimization problem in a prediction error framework (Ljung, 1999).

The minimization of Eq. (6) for lim*N*→∞ can be rewritten into the frequency domain expression

$$\min\_{\theta} \int\_{\pi}^{-\pi} \left| H(e^{jw}) + G(e^{jw}) F(e^{jw}, \theta) \right|^2 \tag{7}$$

using Parceval's theorem (Ljung, 1999). It can be observed that the standard output-error (OE) minimization problem in Eq. (6) can be used to compute the optimal feedforward filter *F*(*q*, *θ*), provided *d*(*t*) and *u f*(*t*) in Eq. (4) are available.

Fig. 2. Block Diagram of the Structural Vibration Control with Adaptive Feedforward Compensation.

For a proper derivation of the adaptation of the feedforward filter, *F*, an approximate model, *G*ˆ, of the control path, *G*, is required to create the filtered signal, *u f*(*t*), for adaptive filtering purpose. The adaptation of the feedforward filter is illustrated in Fig. 2. The filtered signal, *u f*(*t*), and the error signal, *e*(*t*), are used for the computation of the coefficients of the feedforward filter by the adaptive filtering. Thus, the coefficients of the feedforward filter, *F*, can be updated at each time constant for structural vibration reduction.

#### **3. ORThonormal Finite Impulse Response (ORTFIR) filter structure**

In general, the feedforward filter, *F*, in Fig. 1 can be realized by adopting both the finite impulse response (FIR) structure as well as the infinite impulse response (IIR) structure. Because the FIR filter incorporates only zeros, it is always stable and it will provide a linear phase response. It is the most popular adaptive filter widely used in adaptive filtering. 4 Will-be-set-by-IN-TECH

where *d*(*t*) can actually be measured, and *u f*(*t*) is called filtered input signal, Eq. (3) is reduced

*N* ∑ *t*=1 *e*

to compute the optimal feedforward filter, *F*(*q*, *θ*), is a standard output-error (OE)

The minimization of Eq. (6) for lim*N*→∞ can be rewritten into the frequency domain expression

using Parceval's theorem (Ljung, 1999). It can be observed that the standard output-error (OE) minimization problem in Eq. (6) can be used to compute the optimal feedforward filter *F*(*q*, *θ*),

*w*(*t*) *d*(*t*) *+ e*(*t*)

*F G*

For a proper derivation of the adaptation of the feedforward filter, *F*, an approximate model, *G*ˆ, of the control path, *G*, is required to create the filtered signal, *u f*(*t*), for adaptive filtering purpose. The adaptation of the feedforward filter is illustrated in Fig. 2. The filtered signal, *u f*(*t*), and the error signal, *e*(*t*), are used for the computation of the coefficients of the feedforward filter by the adaptive filtering. Thus, the coefficients of the feedforward filter,

In general, the feedforward filter, *F*, in Fig. 1 can be realized by adopting both the finite impulse response (FIR) structure as well as the infinite impulse response (IIR) structure. Because the FIR filter incorporates only zeros, it is always stable and it will provide a linear phase response. It is the most popular adaptive filter widely used in adaptive filtering.

*u*(*t*) *y*(*t*)

*Adaptive Filtering* *+*

*H*

Fig. 2. Block Diagram of the Structural Vibration Control with Adaptive Feedforward

*F*, can be updated at each time constant for structural vibration reduction.

**3. ORThonormal Finite Impulse Response (ORTFIR) filter structure**

|*H*(*ejw*) + *G*(*ejw*)*F*(*ejw*, *θ*)|

min *θ*

minimization problem in a prediction error framework (Ljung, 1999).

 −*π π*

min *θ*

provided *d*(*t*) and *u f*(*t*) in Eq. (4) are available.

*Gust perturbation.*

*Lidar Beam*

*n*(*t*)

1 *N*

*e*(*t*, *θ*) = *d*(*t*) − *F*(*q*, *θ*)*u f*(*t*) (5)

2(*t*, *θ*) (6)

2 (7)

to

for which the minimization

Compensation.

Generally, the discrete time linear time invariant (LTI) FIR filter, *F*(*q*), can be presented as:

$$F(q) = \sum\_{k=0}^{L-1} \beta\_k q^{-k} \tag{8}$$

where *q*−1 denotes the usual time shift operator, *q*−1*x*(*t*) = *x*(*t* − 1). Adaptive filter estimation using FIR filters converges to optimal and unbiased estimates irrespective of the coloring of the noise on the output data. However, a FIR filter is usually too simple to model complex system dynamics such as AE/ASE systems with many resonance modes being excited by atmospheric perturbations. As a result, many tapped delay coefficients of the FIR filter are required to approximate the optimal filter. Even though an IIR filter is appealing as an alternative, the inherent stability and bias estimation problems limit the use of an IIR filter for adaptive filtering in aeroservoelastic systems.

To improve the approximation properties of the adaptive filter, *F*, in Fig. 1, the linear combination of tapped delay functions, *q*−1, in the FIR filter of Eq. (8) can be generalized to the following form:

$$F(q, \theta) = \sum\_{k=0}^{L-1} \beta\_k B\_k(q) \tag{9}$$

where *Bk*(*q*) are generalized (orthonormal) basis functions (Heuberger et al., 1995) that contain some *a-priori* knowledge on the desired filter dynamics. In other words, the orthonormal basis functions that are used in the parametrization of the ORTFIR filter will be tuned on the fly by taking full advantage of the modal information embedded in the flight data.

#### Construction of the Orthonormal Basis Sets

The application of orthonormal basis functions to parameterize and estimate dynamical systems has obtained extensive attention in recent years. Different constructions of the orthonormal basis structure has been reported in (Heuberger et al., 1995; Ninness & Gustafsson, 1997; Zeng & de Callafon, 2006). It is assumed that the pole locations are already known with the use of the standard open-loop prediction error system identification methods. Suppose the poles {*ξi*}*i*=1,2,···,*N* are known, an all pass function, *P*(*q*), can be created by these poles, and is given as

$$P(q) = \prod\_{i=1}^{N} \left[ \frac{1 - \mathfrak{F}\_i^\* q}{q - \mathfrak{F}\_i} \right] \tag{10}$$

Let (*A*, *B*, *C*, *D*) be a minimal balanced realization of an all pass function, *P*(*q*), define the input to state transfer function, *B*0(*q*)=(*q I* − *A*)−1*B*, then a set of functions, *Bi*(*q*), can be obtained via

$$B\_{\dot{1}}(q) = B\_0(q)P^{\dot{1}}(q) \tag{11}$$

and *Bi*(*q*) has orthogonal property

$$\frac{1}{2\pi \text{j}} \oint B\_i(q) B\_k^T(1/q) \frac{\mathbf{d}q}{q} = \begin{cases} \text{I i } i = k\\ \text{0 i } i \neq k \end{cases} \tag{12}$$

The construction of the orthonormal basis function is illustrated in Fig. 3. It should be noted that if *B*0(*q*), includes all information of a dynamical system, then only one parameter, *β*0, needs to be estimated to approximate this dynamic system. It means that the parameters

Gust Loads Alleviation 7

Adaptive Feedforward Control for Gust Loads Alleviation 101

*m*1 = *m*2 = *m*3 = *m*4 = 1 *k*1 = *k*3 = *k*5 = 1750 *k*2 = *k*4 = 2000 *c*1 = *c*3 = *c*5 = 0.7 *c*2 = *c*4 = 0.8

A mathematical model of this lumped system can be easily derived with the use of Newton's second law. The natural frequencies and damping ratios of this lumped system are also obtained. For simplicity purposes, all the units of this 4-DOF lumped system are omitted. This mathematical model is applied in this case example as the real model, and an FIR model and an ORTFIR model will be implemented to approximate this real model, respectively. To facilitate the model estimation using input and output data of the 4-DOF lumped parameter system, a band limited white noise (zero mean) is injected to the 4-DOF lumped parameter system, and an additional band limited white noise (zero mean) is added to the output response to simulate the measurement noise. With the collected input/output data, an FIR filter with varying order is applied to fit the real model, the variance of the simulation error (the difference of the measured and the simulated output) is used to indicate the performance of the FIR filter. Furthermore, the PolyMAX identification method described in Section 4 is applied in this case example to estimate the four physical modes. These estimated physical modes (shown in Table 2, Section 4) are used for the basis function generation of the ORTFIR filter. Finally, an eight order ORTFIR filter is applied to approximate the physical system. The estimation results are shown in Table 1. From Table 1, it is clearly seen that with FIR filter, the optimal FIR filter will be 400th order, with the smallest simulation error 36.18. However, with the simplest eight order ORTFIR filter, the variance of the simulation error is only 13.18, which is almost three times smaller than that of 400th order FIR filter. Fig. 5 compared the model estimation results using 50th/400th order FIR filters and 8th order ORTFIR filter. From Fig. 5, it is observed that with 50th order FIR filter, the essential dynamics of the physical system can hardly be catched. With 400th order FIR filter, even though the physical system can be correctly approximated, the estimated model has evident variation, especially in the high frequency range. On the other hand, with 8th order ORTFIR filter, the physical system can be perfectly approximated, no visible variation of the estimated model was found in a

> Order 1000 500 400 300 200 100 50 8 Filter Type FIR FIR FIR FIR FIR FIR FIR ORTFIR

Variance of Simulation Error 42.53 36.46 36.18 36.70 40.28 52.91 70.10 13.18

A rather general frequency-domain identification method using the standard least squares estimator algorithm is introduced and applied to extract the modal characteristics of a dynamic system from a set of measured data. Consider a set of noisy complex Frequency Response Functions (FRF) measurement data, *G*(*ωj*), (*j* = 1, ··· , *N*). The approximation of

*E*(*ωj*) = *G*(*ωj*) − *P*(*ωj*) *j* = 1, ··· , *N* (14)

the data by a model *P*(*ω*) is addressed by considering the following additive error,

Table 1. Model Estimation Results Using FIR Fiter and ORTFIR Filter.

**4. Modal parameters estimation-frequency PolyMAX identification**

wider frequency range.

(13)

estimated will directly depend on the *a-priori* system information injected into the basis functions, *Bi*(*q*). An important property and advantage of the ORTFIR filter is the knowledge

Fig. 3. ORTFIR Filter Topology.

of the (desired) dynamical behavior can be incorporated throughout the basis functions, *Bi*(*q*). As a result, an accurate description of the filter to be estimated can be achieved by a relatively small number of coefficients.

Case Example: Illustration of the Advantage of Using ORTFIR Filter Over FIR Filter

A 4 degrees-of-freedom (DOF) lumped parameter system is considered to demonstrate the advantage of using ORTFIR filter over FIR filter. An illustration of this 4-DOF lumped parameter system is shown in Fig. 4, where *ki* and *ci*(*i* = 1, ··· , 5) indicate the system stiffness and damping, respectively and *mi*(*i* = 1, ··· , 4) are the masses. The nominal values of these parameters are given as,

Fig. 4. Lumped Parameter System.

6 Will-be-set-by-IN-TECH

estimated will directly depend on the *a-priori* system information injected into the basis functions, *Bi*(*q*). An important property and advantage of the ORTFIR filter is the knowledge

P(q) ... P(q)

of the (desired) dynamical behavior can be incorporated throughout the basis functions, *Bi*(*q*). As a result, an accurate description of the filter to be estimated can be achieved by a relatively

A 4 degrees-of-freedom (DOF) lumped parameter system is considered to demonstrate the advantage of using ORTFIR filter over FIR filter. An illustration of this 4-DOF lumped parameter system is shown in Fig. 4, where *ki* and *ci*(*i* = 1, ··· , 5) indicate the system stiffness and damping, respectively and *mi*(*i* = 1, ··· , 4) are the masses. The nominal values of these

✲

✲

+

β*N*−1 ✻

✻

✲ +

y(t)

✻

B0

✲

β1

✻

+

✲

+ ✻

B0

✻

Case Example: Illustration of the Advantage of Using ORTFIR Filter Over FIR Filter

✻

✲ v(t)

B0

Fig. 3. ORTFIR Filter Topology.

small number of coefficients.

parameters are given as,

Fig. 4. Lumped Parameter System.

β0

✻

✲

...

$$\begin{array}{l}m\_1 = m\_2 = m\_3 = m\_4 = 1\\k\_1 = k\_3 = k\_5 = 1750\\k\_2 = k\_4 = 2000\\c\_1 = c\_3 = c\_5 = 0.7\\c\_2 = c\_4 = 0.8\end{array} \tag{13}$$

A mathematical model of this lumped system can be easily derived with the use of Newton's second law. The natural frequencies and damping ratios of this lumped system are also obtained. For simplicity purposes, all the units of this 4-DOF lumped system are omitted. This mathematical model is applied in this case example as the real model, and an FIR model and an ORTFIR model will be implemented to approximate this real model, respectively. To facilitate the model estimation using input and output data of the 4-DOF lumped parameter system, a band limited white noise (zero mean) is injected to the 4-DOF lumped parameter system, and an additional band limited white noise (zero mean) is added to the output response to simulate the measurement noise. With the collected input/output data, an FIR filter with varying order is applied to fit the real model, the variance of the simulation error (the difference of the measured and the simulated output) is used to indicate the performance of the FIR filter. Furthermore, the PolyMAX identification method described in Section 4 is applied in this case example to estimate the four physical modes. These estimated physical modes (shown in Table 2, Section 4) are used for the basis function generation of the ORTFIR filter. Finally, an eight order ORTFIR filter is applied to approximate the physical system. The estimation results are shown in Table 1. From Table 1, it is clearly seen that with FIR filter, the optimal FIR filter will be 400th order, with the smallest simulation error 36.18. However, with the simplest eight order ORTFIR filter, the variance of the simulation error is only 13.18, which is almost three times smaller than that of 400th order FIR filter. Fig. 5 compared the model estimation results using 50th/400th order FIR filters and 8th order ORTFIR filter. From Fig. 5, it is observed that with 50th order FIR filter, the essential dynamics of the physical system can hardly be catched. With 400th order FIR filter, even though the physical system can be correctly approximated, the estimated model has evident variation, especially in the high frequency range. On the other hand, with 8th order ORTFIR filter, the physical system can be perfectly approximated, no visible variation of the estimated model was found in a wider frequency range.

Table 1. Model Estimation Results Using FIR Fiter and ORTFIR Filter.

#### **4. Modal parameters estimation-frequency PolyMAX identification**

A rather general frequency-domain identification method using the standard least squares estimator algorithm is introduced and applied to extract the modal characteristics of a dynamic system from a set of measured data. Consider a set of noisy complex Frequency Response Functions (FRF) measurement data, *G*(*ωj*), (*j* = 1, ··· , *N*). The approximation of the data by a model *P*(*ω*) is addressed by considering the following additive error,

$$E(\omega\_j) = G(\omega\_j) - P(\omega\_j) \quad j = 1, \dots, N \tag{14}$$

Gust Loads Alleviation 9

Adaptive Feedforward Control for Gust Loads Alleviation 103

rule, unstable poles are not considered in the plot. Physical poles will appear as *stable* poles, independent of the number of the assumed model order. On the other hand, *mathematical* poles that intent to model the noise embedded in the data, will change with the assumed

As an example, the 4-DOF lumped system used in the Section 3 is applied in this section to demonstrate the PolyMAX identification method. Fig. 6 (a) depicts the stabilization diagram for the 4-DOF lumped system where four physical modes can be easily appreciated with the parameter constraint of *Ana* = *I*1. The estimation results can be easily extracted through the access of the stabilization diagram, and they are shown in Table 2. However, with in *An*0 = *I*1, the physical modes are difficult to extract from the stabilization diagram, because with this parameter constraint, all the mathematical poles are also estimated as the stable poles. This phenomenon is illustrated in Fig. 6 (b). In Fig. 6 (a) and (b), the solid curve indicates the frequency response function (FRF) estimate from the input and output time domain data, the dotted curve indicates the estimated model with the highest order 50. The right ordinate axis is the magnitude in dB, and which is used to present in the frequency domain the magnitude of the FRF and the estimated 50th order model. The markers indicate different damping values of each estimated stable poles displayed in the stabilization diagram. The detailed meanings

In Table 2, the second column indicates the frequency and damping of the true modes. The third column presents the estimated frequency and damping of the true modes using the proposed PolyMAX identification method and stabilization diagram. Comparing the estimated modes and real modes (calculated from the mathematical equation of motion of 4-DOF lumped system) in Table 2, it is obvious to see that the frequency, *fi*, and damping,

> *f*1(Hz) 4.1866 4.1868 *f*2(Hz) 7.8648 7.8630 *f*3(Hz) 11.3191 11.3303 *f*4(Hz) 13.1320 13.1320 *ζ*1(%) 0.5261 0.5622 *ζ*2(%) 0.9883 1.1035 *ζ*3(%) 1.4224 1.4836 *ζ*4(%) 1.6502 1.6333

Range of Damping Ratio Marker Sign 0 < *ζ* < 0.1% + 0.1% < *ζ* < 1% × 1% < *ζ* < 2% ∗ 2% < *ζ* < 4% ♦ 4% < *ζ* < 6% ∇ 6% < *ζ* �

Real Modes Identified Modes

*ζi*(*i* = 1, 2, 3, 4), of these four physical modes are estimated consistently.

Table 2. Comparison Between the Identified Modes and Real Modes.

Table 3. Damping Markers in the Stabilization Diagram

model order.

of these damping markers are presented in Table 3.

Fig. 5. Comparison of the Model Estimation Using the FIR and ORTFIR Filters.

Then, it is assumed that the model, *P*(*ω*), can be represented by a right polynomial fraction matrix given by,

$$P(\omega) = [B(\omega)][A(\omega)]^{-1} \tag{15}$$

where *P*(*ω*) ∈ C*p*×*m* is the FRF matrix with *p* outputs and *m* inputs, *B*(*ω*) ∈ C*p*×*m* is the numerator matrix polynomial, and *A*(*ω*) ∈ C*m*×*m* is the denominator matrix polynomial.

The matrix polynomial *B*(*ω*) is parameterized by

$$B(\omega) = \sum\_{k=0}^{n\_b} B\_k \mathfrak{f}\_k(\omega) \tag{16}$$

where *Bk* ∈ R*p*×*m*, and *nb* is the number of non-zero matrix coefficients in *B*(*ω*), or the order of *B*(*ω*). *ξk*(*ω*) are the polynomial basis functions. For continuous time model, *ξk*(*ω*) = −*iωk*. For discrete time model, *ξk*(*ω*) = *e*−*iωkT* (*T* is the sampling time).

The matrix polynomial *A*(*ω*) is parameterized by

$$A(\omega) = \sum\_{k=0}^{n\_d} A\_k \xi\_k(\omega) \tag{17}$$

where *Ak* ∈ R*m*×*m*, and *na* is the number of non-zero matrix coefficients in *A*(*ω*).

Assuming that the coefficients of the denominator, *A*(*ω*), are [*A*0, *A*1, ··· , *Ana* ], then, a constraint, *A*0 = *Im*, is set to obtain a stable model to fit the measured frequency-domain data. Here, a constraint, *Ana* = *Im*, is adopted to extract physical modes from the measured frequency-domain data (Cauberghe et al., 2004).

With, *Ana* = *Im*, the poles of the estimated model are separated into *stable physical* poles and *unstable mathematical* poles, from which a very clean stabilization diagram can be obtained, and the *physical* modal parameters of the real system can be estimated from a quick evaluation of the generated stabilization diagram (Cauberghe et al., 2005).

The stabilization diagram assumes an increasing model order (number of poles noted in the left ordinate axis), and it indicates where on the frequency axis the poles are located. As a 8 Will-be-set-by-IN-TECH

True 8th Order Model 50th Order FIR 400th Order FIR 8th Order ORTFIR

*P*(*ω*)=[*B*(*ω*)][*A*(*ω*)]−1 (15)

*Bk ξk* (*ω*) (16)

*Ak ξk* (*ω*) (17)

100 101 −50

Fig. 5. Comparison of the Model Estimation Using the FIR and ORTFIR Filters.

*B*(*ω*) =

*A*(*ω*) =

where *Ak* ∈ R*m*×*m*, and *na* is the number of non-zero matrix coefficients in *A*(*ω*).

For discrete time model, *ξk*(*ω*) = *e*−*iωkT* (*T* is the sampling time).

of the generated stabilization diagram (Cauberghe et al., 2005).

Then, it is assumed that the model, *P*(*ω*), can be represented by a right polynomial fraction

where *P*(*ω*) ∈ C*p*×*m* is the FRF matrix with *p* outputs and *m* inputs, *B*(*ω*) ∈ C*p*×*m* is the numerator matrix polynomial, and *A*(*ω*) ∈ C*m*×*m* is the denominator matrix polynomial.

> *nb* ∑ *k*=0

where *Bk* ∈ R*p*×*m*, and *nb* is the number of non-zero matrix coefficients in *B*(*ω*), or the order of *B*(*ω*). *ξk*(*ω*) are the polynomial basis functions. For continuous time model, *ξk*(*ω*) = −*iωk*.

> *na* ∑ *k*=0

Assuming that the coefficients of the denominator, *A*(*ω*), are [*A*0, *A*1, ··· , *Ana* ], then, a constraint, *A*0 = *Im*, is set to obtain a stable model to fit the measured frequency-domain data. Here, a constraint, *Ana* = *Im*, is adopted to extract physical modes from the measured

With, *Ana* = *Im*, the poles of the estimated model are separated into *stable physical* poles and *unstable mathematical* poles, from which a very clean stabilization diagram can be obtained, and the *physical* modal parameters of the real system can be estimated from a quick evaluation

The stabilization diagram assumes an increasing model order (number of poles noted in the left ordinate axis), and it indicates where on the frequency axis the poles are located. As a

−40 −30 −20 −10 0 10 20 30 40 50

The matrix polynomial *B*(*ω*) is parameterized by

The matrix polynomial *A*(*ω*) is parameterized by

frequency-domain data (Cauberghe et al., 2004).

matrix given by,

Magnitude [dB]

rule, unstable poles are not considered in the plot. Physical poles will appear as *stable* poles, independent of the number of the assumed model order. On the other hand, *mathematical* poles that intent to model the noise embedded in the data, will change with the assumed model order.

As an example, the 4-DOF lumped system used in the Section 3 is applied in this section to demonstrate the PolyMAX identification method. Fig. 6 (a) depicts the stabilization diagram for the 4-DOF lumped system where four physical modes can be easily appreciated with the parameter constraint of *Ana* = *I*1. The estimation results can be easily extracted through the access of the stabilization diagram, and they are shown in Table 2. However, with in *An*0 = *I*1, the physical modes are difficult to extract from the stabilization diagram, because with this parameter constraint, all the mathematical poles are also estimated as the stable poles. This phenomenon is illustrated in Fig. 6 (b). In Fig. 6 (a) and (b), the solid curve indicates the frequency response function (FRF) estimate from the input and output time domain data, the dotted curve indicates the estimated model with the highest order 50. The right ordinate axis is the magnitude in dB, and which is used to present in the frequency domain the magnitude of the FRF and the estimated 50th order model. The markers indicate different damping values of each estimated stable poles displayed in the stabilization diagram. The detailed meanings of these damping markers are presented in Table 3.

In Table 2, the second column indicates the frequency and damping of the true modes. The third column presents the estimated frequency and damping of the true modes using the proposed PolyMAX identification method and stabilization diagram. Comparing the estimated modes and real modes (calculated from the mathematical equation of motion of 4-DOF lumped system) in Table 2, it is obvious to see that the frequency, *fi*, and damping, *ζi*(*i* = 1, 2, 3, 4), of these four physical modes are estimated consistently.

Gust Loads Alleviation 11

Adaptive Feedforward Control for Gust Loads Alleviation 105

*θT*(*t* − 1)*φ*(*t*) (20)

*P*(*t* − 1) (22)

*θ*(*t*)*Tφ*(*i*)]2 (24)

*λ*1(*t*) + *φT*(*t*)*P*(*t* − 1)*φ*(*t*) (21)

*ξ*(*t*) = *y*(*t*) − ˆ

2. Update the inverse correlation matrix, *P*(*t*), and the forgetting factor, *λ*(*t*):

*P*(*t*) = *λ*(*t*)−1

infinity, *λ* reaches its minimum value. The RLS minimization is posted as:

*t* ∑ *i*=1

estimation when the parameters, *θ*, are still far from the optimal value.

*b* + Ω*b* × *Vb* − *Rbg*(*E*) [0, 0, *g*]

*λ*(*i*)*t*−*i*

*J*(*t*) =

**6. Application to closed loop F/A-18 AAW linear model**

**6.1 Linear aeroelastic solver formulation approach**

*m V*˙

from inertial to body-axes, (*E* = [ *φ*, *θ*, *ψ* ]).

*k*(*t*) = *P*(*t* − 1)*φ*(*t*)

1 − *k*(*t*)*φT*(*t*)

*λ*(*t*) = *λmin* + (1 − *λmin*) · 2−*L*(*t*)

where *ρ* is a design parameter which controls the change rate and the width of a unity zone, *ξ*(*t*) is the estimation error which is calculated via Eq. (20). *L*(*t*) is defines as the nearest integer of *ρ* · *ξ*(*t*)2 at each instant time step. *λmin* defines the lower bound of the *λ*.

In Eq. (23), it is shown that when the estimation error, *ξ*(*t*) and *L*(*t*) is small, 2−*L*(*t*) −→ 1, and *λ*(*t*) −→ 1 at an exponential rate, and this rate is controller by *ρ*. When *ξ*(*t*) increases to

By choosing the variable forgetting factor indicated in Eq. (23), the fast decrease of the inverse correlation matrix, *P*(*t*), can be avoided at the beginning of the estimation. In general this will result in an accelerated convergence by maintaining a high adaptation at the beginning of the

A unified aeroelastic formulation to take into account the influence of aeroelastic effects on the flight dynamic behavior of the whole aircraft has been developed in (Baldelli et al., 2006). A general formulation of a flexible aircraft with respect to a body-fixed reference system driven by aerodynamic, thrust, and gravity (*g*) forces and moments can be defined as:

where, *m* and *J* are the air vehicle mass and inertia tensor, and *Rbg*(*E*) is the rotation mapping

Eq. (25) is driven by the forces and moments on its right hand side, where *FA* and *MA* are the external aerodynamic forces and moments on the air vehicle. *FA* and *MA* are a function of the aerodynamic flight states (*V α*, *β*, *E*,..., *etc*.), Mach number, body angular rates (Ω*b*), and control surface deflections and are usually obtained by wind-tunnel or flight tests. In either case, the quasi-steady influence of the deformed air-vehicle is included by considering flexible-to-rigid ratios or Parameter Identification (PID) techniques, (Morelli, 1995; Morelli & Klein, 1997). *Fδ* and *Mδ* are the aerodynamic forces and moments from the control surfaces

*T*

= *FA* + *Fδ* + *FT* + Δ*F*

*J*Ω˙ *b* + Ω*b* × *J*Ω*b* = *MA* + *Mδ* + *MT* + Δ*M* (25)

[*y*(*i*) − ˆ

*L*(*t*) = *round*(*ρ* · *ξ*(*t*)2) (23)

(b) Stabilization Diagram with *An*0 = *I*1

Fig. 6. Illustration of Stabilization Diagram Using PolyMAX Identification Method (for the Meaning of the Markers, See Table 3).

#### **5. Recursive least square adaptive algorithm**

The adaptive algorithm to be implemented is the Recursive Least Square (RLS) algorithm (Haykin, 2002). Given the input and output data can be written in regressor form:

$$y(t) = \boldsymbol{\phi}^T(t)\boldsymbol{\theta} + \boldsymbol{e}(t), \; \boldsymbol{\theta} = \begin{bmatrix} \boldsymbol{\beta}\_{0\prime}\boldsymbol{\beta}\_{1\prime}\ldots\boldsymbol{\beta}\_{L-1} \end{bmatrix}^T \tag{18}$$

where *φT*(*t*)=[*uT* 0 (*t*), ..., *uT L*−1(*t*)] is the available input data vector, *θ* is the parameter vector to be estimated of the ORTFIR feedforward controller, and *e*(*t*) is the residue error. The parameters, *θ*, can be identified with the available input-output data up to time *t* by a standard RLS algorithm. It is well known that the RLS algorithm at the steady-state operation exhibits a windup problem if the forgetting factor remains constant, which will deteriorate the estimation results. As a result, a variable forgetting factor (Park, 1991), is sought to prevent this problem from occurring. The parameters, *θ*, can be estimated by the RLS algorithm using a variable forgetting factor through a two steps approach at each sample time:

1. Compute the gain vector, *k*(*t*), and the parameters, ˆ *θ*(*t*), at the current sample time as:

$$
\hat{\theta}(t) = \hat{\theta}(t-1) + k(t)\hat{\xi}^T(t) \tag{19}
$$

10 Will-be-set-by-IN-TECH

0 10 20 30 40 50 30

(a) Stabilization Diagram with *Ana* = *I*1

0 10 20 30 40 50 30

(b) Stabilization Diagram with *An*0 = *I*1

The adaptive algorithm to be implemented is the Recursive Least Square (RLS) algorithm

*y*(*t*) = *φT*(*t*)*θ* + *e*(*t*), *θ* = [*β*0, *β*1, ..., *βL*−1]

vector to be estimated of the ORTFIR feedforward controller, and *e*(*t*) is the residue error. The parameters, *θ*, can be identified with the available input-output data up to time *t* by a standard RLS algorithm. It is well known that the RLS algorithm at the steady-state operation exhibits a windup problem if the forgetting factor remains constant, which will deteriorate the estimation results. As a result, a variable forgetting factor (Park, 1991), is sought to prevent this problem from occurring. The parameters, *θ*, can be estimated by the RLS algorithm using

(Haykin, 2002). Given the input and output data can be written in regressor form:

a variable forgetting factor through a two steps approach at each sample time:

ˆ *θ*(*t*) = ˆ

Fig. 6. Illustration of Stabilization Diagram Using PolyMAX Identification Method (for the

Frequency [Hz]

Frequency [Hz]

−44

−57

*L*−1(*t*)] is the available input data vector, *θ* is the parameter

−39.4

−21.8

Magnitude [dB]

*T* (18)

*θ*(*t*), at the current sample time as:

*θ*(*t* − 1) + *k*(*t*)*ξT*(*t*) (19)

−4.2

13.4

31

−29.2

−14.4

Magnitude [dB]

0.4

15.2

30

40

Model Order

40

Model Order

**5. Recursive least square adaptive algorithm**

0 (*t*), ..., *uT*

1. Compute the gain vector, *k*(*t*), and the parameters, ˆ

Meaning of the Markers, See Table 3).

where *φT*(*t*)=[*uT*

50

50

$$\dot{\xi}(t) = y(t) - \dot{\theta}^T(t-1)\phi(t) \tag{20}$$

$$k(t) = \frac{P(t-1)\phi(t)}{\lambda\_1(t) + \phi^T(t)P(t-1)\phi(t)}\tag{21}$$

2. Update the inverse correlation matrix, *P*(*t*), and the forgetting factor, *λ*(*t*):

$$P(t) = \lambda(t)^{-1} \left[ 1 - k(t) \phi^T(t) \right] P(t-1) \tag{22}$$

$$\begin{array}{l} \lambda(t) = \lambda\_{\text{min}} + (1 - \lambda\_{\text{min}}) \cdot 2^{-L(t)}\\ L(t) = \quad \text{round}(\rho \cdot \xi(t)^2) \end{array} \tag{23}$$

where *ρ* is a design parameter which controls the change rate and the width of a unity zone, *ξ*(*t*) is the estimation error which is calculated via Eq. (20). *L*(*t*) is defines as the nearest integer of *ρ* · *ξ*(*t*)2 at each instant time step. *λmin* defines the lower bound of the *λ*.

In Eq. (23), it is shown that when the estimation error, *ξ*(*t*) and *L*(*t*) is small, 2−*L*(*t*) −→ 1, and *λ*(*t*) −→ 1 at an exponential rate, and this rate is controller by *ρ*. When *ξ*(*t*) increases to infinity, *λ* reaches its minimum value. The RLS minimization is posted as:

$$f(t) = \sum\_{i=1}^{t} \lambda(i)^{t-i} [y(i) - \hat{\theta}(t)^T \phi(i)]^2 \tag{24}$$

By choosing the variable forgetting factor indicated in Eq. (23), the fast decrease of the inverse correlation matrix, *P*(*t*), can be avoided at the beginning of the estimation. In general this will result in an accelerated convergence by maintaining a high adaptation at the beginning of the estimation when the parameters, *θ*, are still far from the optimal value.

#### **6. Application to closed loop F/A-18 AAW linear model**

#### **6.1 Linear aeroelastic solver formulation approach**

A unified aeroelastic formulation to take into account the influence of aeroelastic effects on the flight dynamic behavior of the whole aircraft has been developed in (Baldelli et al., 2006).

A general formulation of a flexible aircraft with respect to a body-fixed reference system driven by aerodynamic, thrust, and gravity (*g*) forces and moments can be defined as:

$$m\left[\dot{V}\_b + \Omega\_b \times V\_b - \mathcal{R}\_b \underline{\boldsymbol{\beta}}(E) \begin{bmatrix} \boldsymbol{0}, \boldsymbol{0}, \boldsymbol{\varrho} \end{bmatrix}^T\right] = F\_A + F\_\delta + F\_T + \Delta\_F$$

$$J\dot{\Omega}\_b + \Omega\_b \times f\Omega\_b = M\_A + M\_\delta + M\_T + \Delta\_M \tag{25}$$

where, *m* and *J* are the air vehicle mass and inertia tensor, and *Rbg*(*E*) is the rotation mapping from inertial to body-axes, (*E* = [ *φ*, *θ*, *ψ* ]).

Eq. (25) is driven by the forces and moments on its right hand side, where *FA* and *MA* are the external aerodynamic forces and moments on the air vehicle. *FA* and *MA* are a function of the aerodynamic flight states (*V α*, *β*, *E*,..., *etc*.), Mach number, body angular rates (Ω*b*), and control surface deflections and are usually obtained by wind-tunnel or flight tests. In either case, the quasi-steady influence of the deformed air-vehicle is included by considering flexible-to-rigid ratios or Parameter Identification (PID) techniques, (Morelli, 1995; Morelli & Klein, 1997). *Fδ* and *Mδ* are the aerodynamic forces and moments from the control surfaces

Gust Loads Alleviation 13

Adaptive Feedforward Control for Gust Loads Alleviation 107

Usually, the aerodynamic force coefficient matrix, *Q*(*s*), is approximated using the Rational

*L*2 *V*2 *s*

In this formulation, the [*Ai*] coefficient matrices represent the quasi-steady aerodynamic forces, and the remanent terms are used to model the flow unsteadiness by Padé

Using the Minimum-State approach during the RFA implemented in the ZAERO/ASE module (Karpel, 1992), and due to the performed similarity transformation the aero-lag terms are

By including Eqs. (30) and (32) into Eq. (28), the aeroelastic EoM can now be easily partitioned in accordance with the airframe degrees of freedom, elastic dynamics, aerodynamic lag terms,

> *ξas* 000 *Arr*0 *Arr*1 *Are*0 *Are*1 *ArL* 0 0 0 *I* 0 *Aer*0 *Aer*1 *Aee*0 *Aee*1 *AeL EL*∗ *ELr* 0 *ELe*

⎤ ⎥ ⎥ ⎥ ⎥ ⎦

⎧ ⎨ ⎩

Now, the aeroelastic incremental loads, Δ*F* and Δ*M*, should be implemented in a way to allow a seamless integration between the nonlinear flight dynamics and the linear aeroelastic EoMs. In fact, this can be easily achieved in accordance with the partitions showed in Eq. (34) between rigid, elastic and aerodynamic lag dynamics. Hence, the aeroelastic incremental

*δu δu*˙ *δu*¨

*V L R*

> ⎫ ⎬ ⎭ +

*ξas* and *EL*∗ are coupling matrices due to the similarity transformation executed.

*EL*∗ *ELr ELe ELδ*

⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪⎪⎪⎪⎨

*ξas* ˙ *ξas ηe η*˙*e xL*

⎫ ⎪⎪⎪⎪⎬

⎪⎪⎪⎪⎭

⎤ ⎥ ⎥ ⎥ ⎥ ⎦

� *wG w*˙ *G* �

0 0 *Brw*1 *Brw*2 0 0 *Bew*1 *Bew*2 0 *ErG*

⎪⎪⎪⎪⎩

⎡ ⎢ ⎢ ⎢ ⎢ ⎣ � ⎧ ⎪⎪⎨

⎪⎪⎩

*ξas* ˙ *ξas η*˙*e δu*˙

� (33)

⎫ ⎪⎪⎬

⎪⎪⎭

2 + [*D*]

*Ar Ae Aδ*

� *sI* − *V L* [*R*]

*Mrr* = *diag*(*mI*3, *J*) (29)

�−1

� (31)

[*E*]*s* (30)

(32)

(34)

off-diagonal terms contain the products of inertia,

Function Approximation (RFA) approach as

approximation.

computed as,

where *Aξas*, *A* ˙

*Q*(*s*)=[*A*0]+[*A*1]

{*x*˙*L*} = *V*

*zL* = �

a set of control inputs and gust perturbation as:

⎫ ⎪⎪⎪⎪⎬

⎪⎪⎪⎪⎭

=

⎡ ⎢ ⎢ ⎢ ⎢ ⎣

+

⎡ ⎢ ⎢ ⎢ ⎢ ⎣

⎧ ⎪⎪⎪⎪⎨

˙ *ξas* ¨ *ξas η*˙*e η*¨*e x*˙ *L*

⎪⎪⎪⎪⎩

*L V s* + [*A*2]

where *i* = *r*, *e*, *δ* are the airframe, elastic and control related states.

where the [*Ai*], *i* = 0, 1, 2, [*D*] and [*E*] matrices are column partitioned as,

*L* [*R*] {*xL*} + �

*D* � � *xL*

*Aξas A* ˙

000 *Br*0 *Br*1 *Br*2 000 *Be*0 *Be*1 *Be*2 0 *ELδ* 0

[*Ai*] = �

commanded by the flight control system and pilot inputs while *FT* and *MT* includes the thrust loads.

In addition, Δ*F* and Δ*M* are the *aeroelastic incremental loads* due to the structural deformation. Usually, these loads are assumed to be quasi-statics and can be computed by a static aeroelastic analysis. However, this quasi-static assumption may not be sufficient for a highly reconfigurable and flexible aircraft like the new generation of Morphing UAV, HALEs, etc., where the interaction between the dynamic structural deformation due to unsteady flow and rigid body motion can play an important role.

During the integration process, the aeroelastic equations of motion underwent two similarity transformation steps, so the generalized coordinates related with the six rigid-body modes originally defined in principle axes are mapped into the airframe states (stability-axes definition). Specifically, for symmetric maneuvers the transformation matrix, [*TA*]*long* (Baldelli et al., 2006), reads as:

$$\begin{Bmatrix} T\_x \\ T\_z \\ R\_y \\ T\_x \\ \dot{T}\_x \\ \dot{R}\_y \end{Bmatrix} = \underbrace{\begin{bmatrix} -1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -V & V & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix}}\_{\text{\(}T\_A\)\_{\text{long}}} \begin{Bmatrix} \chi \\ \mu \\ h \\ \alpha \\ \theta \\ q \end{Bmatrix} \tag{26}$$

For an anti-symmetric maneuver, [*TA*]*lat* is,

$$\begin{Bmatrix} T\_y \\ R\_x \\ R\_z \\ \dot{T}\_y \\ \dot{R}\_x \\ \dot{R}\_z \end{Bmatrix} = \underbrace{\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 & 0 & -1 \\ 0 & V & 0 & 0 & 0 & V \\ 0 & 0 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & -1 & 0 & 0 \end{bmatrix}}\_{[T\_A]\_{lat}} \begin{Bmatrix} \Psi \\ \beta \\ p \\ r \\ \Phi \\ \Psi \end{Bmatrix} \tag{27}$$

For an asymmetric maneuver, the matrix [*TA*] ∈ **R**12×12 will be composed by the proper allocation of the elements that form the rows and columns of the [*TA*]*long* and [*TA*]*lat* matrices.

In this new coordinate system, the linear aeroelastic Equations of Motion (EoM) are:

$$\begin{aligned} \left\{ \left[ \begin{matrix} M\_{rr} & 0\\ 0 & M\_{\text{er}} \end{matrix} \right] s^2 + \begin{bmatrix} \mathbf{C}\_{rr} & 0\\ 0 & \mathbf{C}\_{\text{er}} \end{bmatrix} s + \begin{bmatrix} K\_{rr} & 0\\ 0 & K\_{\text{er}} \end{bmatrix} \right\} \begin{Bmatrix} \xi\_{\text{as}}\\ \eta\_{\text{er}} \end{Bmatrix} = \\\ q\_{\text{os}} \left\{ \left[ \begin{matrix} Q\_{rr}(s) & Q\_{rr}(s) \\ Q\_{\text{er}}(s) & Q\_{\text{er}}(s) \end{matrix} \right] \begin{Bmatrix} \xi\_{\text{as}}\\ \eta\_{\text{er}} \end{Bmatrix} \right\} + \begin{bmatrix} Q\_{r\delta}(s) \\ Q\_{\text{er}}(s) \end{bmatrix} \delta\_{u} + \begin{bmatrix} 1 & Q\_{r\mathcal{G}}(s) \\ V & Q\_{\text{er}}(s) \end{bmatrix} w\_{\text{G}} \right\} \end{aligned} \tag{28}$$

where *wG* is the gust input; the elastic generalized coordinates, *ηe*, input, and *δu*, vectors are,

$$\begin{aligned} \eta\_{\varepsilon}^{T} &= \begin{bmatrix} \eta\_{\varepsilon\_{1'}} \dots \eta\_{\varepsilon\_{Ne}} \end{bmatrix}^{T} \\ \delta\_{\mu}^{T} &= \begin{bmatrix} \delta\_{elec\nu} \ \delta\_{ail\nu} \ \delta\_{rud\nu} \dots \end{bmatrix}^{T} \end{aligned}$$

It should be noted that the equations are only coupled via external forces and moments. In addition, after the transformation is applied the generalized mass matrix of the finite element model, it is no longer necessarily diagonal. In fact, the sub-matrix, *Mrr*, associated with the rigid body modes is identical to the mass matrix in the flight dynamics equation; i.e. the

12 Will-be-set-by-IN-TECH

commanded by the flight control system and pilot inputs while *FT* and *MT* includes the thrust

In addition, Δ*F* and Δ*M* are the *aeroelastic incremental loads* due to the structural deformation. Usually, these loads are assumed to be quasi-statics and can be computed by a static aeroelastic analysis. However, this quasi-static assumption may not be sufficient for a highly reconfigurable and flexible aircraft like the new generation of Morphing UAV, HALEs, etc., where the interaction between the dynamic structural deformation due to unsteady flow and

During the integration process, the aeroelastic equations of motion underwent two similarity transformation steps, so the generalized coordinates related with the six rigid-body modes originally defined in principle axes are mapped into the airframe states (stability-axes definition). Specifically, for symmetric maneuvers the transformation matrix, [*TA*]*long* (Baldelli

> −1 0 0 0 00 0 0 1 0 00 0 0 0 0 10 0 −10 0 0 0 0 00 −*V V* 0 0 0 0 0 01

� �� � [*TA*]*long*

10 0 0 0 0 00 0 0 −1 0 00 0 0 0 −1 0 *V* 000 *V* 0 0 −10 0 0 00 0 −10 0

� �� � [*TA*]*lat*

For an asymmetric maneuver, the matrix [*TA*] ∈ **R**12×12 will be composed by the proper allocation of the elements that form the rows and columns of the [*TA*]*long* and [*TA*]*lat* matrices.

> � � *ξas ηe* � + � *Qrδ*(*s*) *Qeδ*(*s*)

where *wG* is the gust input; the elastic generalized coordinates, *ηe*, input, and *δu*, vectors are,

*ηe*1 , ..., *ηeNe*

It should be noted that the equations are only coupled via external forces and moments. In addition, after the transformation is applied the generalized mass matrix of the finite element model, it is no longer necessarily diagonal. In fact, the sub-matrix, *Mrr*, associated with the rigid body modes is identical to the mass matrix in the flight dynamics equation; i.e. the

In this new coordinate system, the linear aeroelastic Equations of Motion (EoM) are:

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪⎪⎪⎪⎪⎨ *x u h α θ q* ⎫ ⎪⎪⎪⎪⎪⎬

⎪⎪⎪⎪⎪⎭

(26)

(27)

(28)

⎪⎪⎪⎪⎪⎩

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪⎪⎪⎪⎪⎨ *y β p r φ ψ* ⎫ ⎪⎪⎪⎪⎪⎬

⎪⎪⎪⎪⎪⎭

⎪⎪⎪⎪⎪⎩

�� � *ξas ηe* � =

�*T*

*δelev*, *δail*, *δrud*, ... �*T*

� *δu* + � 1 *V*

*QrG*(*s*) *QeG*(*s*) � *wG* �

loads.

et al., 2006), reads as:

rigid body motion can play an important role.

For an anti-symmetric maneuver, [*TA*]*lat* is,

�� *Mrr* 0 0 *Mee* � *s* 2 + � *Crr* 0 0 *Cee* � *s* + � *Krr* 0 0 *Kee*

*q*∞

⎧ ⎪⎪⎪⎪⎪⎨

*Tx Tz Ry T*˙ *x T*˙ *z R*˙ *y* ⎫ ⎪⎪⎪⎪⎪⎬ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

⎪⎪⎪⎪⎪⎭ =

⎪⎪⎪⎪⎪⎩

⎧ ⎪⎪⎪⎪⎪⎨

*Ty Rx Rz T*˙ *y R*˙ *x R*˙ *z* ⎫ ⎪⎪⎪⎪⎪⎬ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

⎪⎪⎪⎪⎪⎭ =

�� *Qrr*(*s*) *Qre*(*s*) *Qer*(*s*) *Qee*(*s*)

> *ηT e* = �

*δT u* = �

⎪⎪⎪⎪⎪⎩

off-diagonal terms contain the products of inertia,

$$M\_{\mathcal{II}} = \text{diag}(mI\_{\mathfrak{J}'}I) \tag{29}$$

Usually, the aerodynamic force coefficient matrix, *Q*(*s*), is approximated using the Rational Function Approximation (RFA) approach as

$$Q(s) = [A\_0] + [A\_1]\frac{L}{V}s + [A\_2]\frac{L^2}{V^2}s^2 + [D]\left(sI - \frac{V}{L}[R]\right)^{-1}[E]s \tag{30}$$

where the [*Ai*], *i* = 0, 1, 2, [*D*] and [*E*] matrices are column partitioned as,

$$\begin{bmatrix} A\_{l} \end{bmatrix} = \begin{bmatrix} A\_{r} \ A\_{\ell} \ A\_{\delta} \end{bmatrix} \tag{31}$$

where *i* = *r*, *e*, *δ* are the airframe, elastic and control related states.

In this formulation, the [*Ai*] coefficient matrices represent the quasi-steady aerodynamic forces, and the remanent terms are used to model the flow unsteadiness by Padé approximation.

Using the Minimum-State approach during the RFA implemented in the ZAERO/ASE module (Karpel, 1992), and due to the performed similarity transformation the aero-lag terms are computed as,

$$\begin{Bmatrix} \{\dot{\mathbf{x}}\_{L}\} = \frac{V}{L} \begin{bmatrix} \mathbf{R} \end{bmatrix} \begin{Bmatrix} \mathbf{x}\_{L} \end{Bmatrix} + \begin{bmatrix} \mathbf{E}\_{L\*} \ \mathbf{E}\_{L\*} \ \mathbf{E}\_{L\varepsilon} \ \mathbf{E}\_{L\delta} \end{Bmatrix} \begin{Bmatrix} \begin{Bmatrix} \mathbf{\tilde{\xi}}\_{as} \\ \mathbf{\tilde{\xi}}\_{as} \\ \dot{\mathbf{\tilde{\eta}}\_{\varepsilon}} \\ \boldsymbol{\delta}\_{\dot{u}} \end{Bmatrix} \end{Bmatrix} \tag{32}$$

$$\mathbf{z}\_{\rm L} = \begin{bmatrix} D \end{bmatrix} \begin{Bmatrix} \mathbf{x}\_{\rm L} \end{Bmatrix} \tag{33}$$

By including Eqs. (30) and (32) into Eq. (28), the aeroelastic EoM can now be easily partitioned in accordance with the airframe degrees of freedom, elastic dynamics, aerodynamic lag terms, a set of control inputs and gust perturbation as:

⎧ ⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎩ ˙ *ξas* ¨ *ξas η*˙*e η*¨*e x*˙ *L* ⎫ ⎪⎪⎪⎪⎬ ⎪⎪⎪⎪⎭ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ *Aξas A* ˙ *ξas* 000 *Arr*0 *Arr*1 *Are*0 *Are*1 *ArL* 0 0 0 *I* 0 *Aer*0 *Aer*1 *Aee*0 *Aee*1 *AeL EL*∗ *ELr* 0 *ELe V L R* ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎩ *ξas* ˙ *ξas ηe η*˙*e xL* ⎫ ⎪⎪⎪⎪⎬ ⎪⎪⎪⎪⎭ + ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 000 *Br*0 *Br*1 *Br*2 000 *Be*0 *Be*1 *Be*2 0 *ELδ* 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎨ ⎩ *δu δu*˙ *δu*¨ ⎫ ⎬ ⎭ + ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 *Brw*1 *Brw*2 0 0 *Bew*1 *Bew*2 0 *ErG* ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ � *wG w*˙ *G* � (34)

where *Aξas*, *A* ˙ *ξas* and *EL*∗ are coupling matrices due to the similarity transformation executed. Now, the aeroelastic incremental loads, Δ*F* and Δ*M*, should be implemented in a way to allow a seamless integration between the nonlinear flight dynamics and the linear aeroelastic EoMs. In fact, this can be easily achieved in accordance with the partitions showed in Eq. (34) between rigid, elastic and aerodynamic lag dynamics. Hence, the aeroelastic incremental

Gust Loads Alleviation 15

Adaptive Feedforward Control for Gust Loads Alleviation 109

The previous equation is used to estimate the elastic and lag states as a function of the

The quasi-static elastic deformation , *ηe*0, is computed by static residualization of the elastic modes, that is the *η*˙*e* = *η*¨*e* = *xLe* = 0 condition needs to be fulfilled. Therefore, the quasi-static

2. First-order differential Eq. (40) to compute the generalized coordinates related vectors *η*¯*e*, *η*˙*e*, and *η*¨*e*, as well as the aerodynamic lag terms related with the elastic modes, *xLe* . 3. Algebraic Eq. (41) to estimate the quasi-static deformation vector *ηe*0 at that specific flight

Fig. 7 illustrates the interconnection of the F/A-18 AAW 6-DOF Dynamics subsystem and the Incremental Aeroelastic Solver, Control Surface Mixer and Control

> \*HQHUDO B&RRUGLQDWH HWD BHHH

Fig. 7. Addition of the Incremental Aeroelastic Loads Solver to the Nonlinear Rigid-Body

In order to demonstrate the proposed feedforward filter design algorithm, a simplified closed loop F/A-18 AAW linear simulink model with gust excitation is developed/implemented for the evaluation purposes. This high-fidelity aeroelastic model was developed using the

HWD BH

,QFUHPHQWDO)RUFHV 0RPHQWV6ROYHU

)\$ \$\$:5LJLG '2)

) \$ \$\$:2XWSXWV

ERROHDQ HQDEOH[

&RQWURO(IIHFWRUV

'HOWD)

'HOWD0

*x*˙*e* = 0 =⇒ *xe* = −*A*−1(*B*1 *δU* + *B*2 *δxi* + *B*3 *w*) (41)

*T* ), and incremental airframe states, (*δT*

*ξ* =

)\$ \$\$:2XWSXWV

*u* , *δT u*˙ , *δT u*¨ ]

*U* = [ *δT*

and from *xe* the quasi-static elastic influence vector, *ηe*0, can be recovered.

1. Algebraic Eq. (35) to compute the incremental aeroelastic loads, Δ*F* and Δ*M*.

 (40) as:

[ *δT ξas*, *δT* ˙ *ξas* ]

condition.

&RQWURO6XUIDFH

loads are computed similarly to the approximation given by Eq. (30),

$$\begin{Bmatrix} \Delta\_F\\ \Delta\_M \end{Bmatrix} = q\_{\infty} \left\{ A\_{0\epsilon} \vec{\eta}\_{\ell} + A\_{1\epsilon\epsilon} \frac{L}{V} \vec{\eta}\_{\ell} + A\_{2\epsilon\epsilon} \frac{L^2}{V^2} \vec{\eta}\_{\ell} + D\_{\ell\epsilon} \mathbf{x}\_{L\epsilon} \right\} \tag{35}$$

Clearly to implement this algebraic equation, the generalized coordinate, *η*¯*e* = *ηe* − *ηe*0 , its rate, *η*˙*e*, and acceleration, *η*¨*e*, vectors as well as the aerodynamic lag terms related with the elastic modes, *xLe* , must be estimated at each time iteration.

Note that the Minimum-State method is formulated to only use a single set of lag states, *xL*, in Eq. (32). Therefore, the following augmented equation is devised to decouple the generalized coordinates' aero lag terms from the airframe states, *ξas* and ˙ *ξas*, related ones,

$$
\begin{Bmatrix}
\dot{\mathbf{x}}\_{L\_{\text{os}}} \\
\dot{\mathbf{x}}\_{L\_{\text{e}}}
\end{Bmatrix} = \begin{bmatrix}
\frac{V}{L}R & 0 \\
0 & \frac{V}{L}R
\end{bmatrix} \begin{Bmatrix}
\mathbf{x}\_{L\_{\text{os}}} \\
\mathbf{x}\_{L\_{\text{e}}}
\end{Bmatrix} + \begin{bmatrix}
E\_{L\star} \ E\_{L\prime} \ 0 \ E\_{L\prime} \\
0 & 0 \ 0 \ E\_{L\prime}
\end{bmatrix} \begin{Bmatrix}
\dot{\mathbf{x}}\_{\delta\dot{s}} \\
\dot{\mathbf{y}}\_{\delta\dot{s}} \\
\dot{\mathbf{y}}\_{\delta\varepsilon}
\end{Bmatrix}
$$

$$
+ \begin{bmatrix}
0 \ 0.5 E\_{L\delta} \ 0 \\
0 \ 0.5 E\_{L\delta} \ 0
\end{bmatrix} \begin{Bmatrix}
\delta\_{\dot{u}} \\
\delta\_{\dot{u}}
\end{Bmatrix} + \begin{bmatrix}
0 \ 0.5 E\_{rG} \\
0 \ 0.5 E\_{rG}
\end{bmatrix} \begin{Bmatrix}
w\_{G} \\
\dot{w}\_{G}
\end{Bmatrix}
\tag{36}
$$

$$
\begin{Bmatrix}
\mathbf{z}\_{L\_{\text{as}}} \\
z\_{L\_{\text{e}}}
\end{Bmatrix} = \begin{bmatrix}
D \ D
\end{bmatrix} \begin{Bmatrix}
\mathbf{x}\_{L\_{\text{as}}} \\
\mathbf{x}\_{L\_{\text{e}}}
\end{Bmatrix}
$$

In this way, only elastic lag terms are considered to avoid any possible coupling with the rigid-body airframe related states (i.e. *ξas* and *ξas*˙ ). Now, the following differential equation is obtained combining the lower partition of Eq. (34) with the new devised Eq. (36).

$$
\underbrace{\begin{aligned}
\begin{Bmatrix}\dot{\eta}\_{\ell}\\ \dot{\eta}\_{\ell}\\ \dot{x}\_{L\_{\text{os}}}\\ \dot{x}\_{L\_{\text{os}}}\end{Bmatrix}}\_{\dot{\mathbf{x}}\_{\mathcal{E}}} &= \underbrace{\begin{bmatrix}0&I&0&0\\ \frac{A\_{\text{er}\_{0}}}{0}\underbrace{A\_{\text{er}\_{1}}}\_{L\_{\text{er}}}\underbrace{A\_{\text{el}}}\_{L\_{\text{el}}}\underbrace{A\_{\text{el}}}\_{L\_{\text{el}}}\\ 0&\underbrace{E\_{L\varepsilon}\left[\frac{\mathbf{V}}{\mathcal{E}}\,\mathbf{0}\right]}\_{\mathcal{E}}\underbrace{\begin{bmatrix}\eta\_{\varepsilon}\\ \dot{\eta}\_{\varepsilon}\\ \mathbf{x}\_{L\_{\text{e}}}\end{bmatrix}}\_{\mathcal{X}\_{\mathcal{E}}} &\underbrace{\begin{bmatrix}\eta\_{\varepsilon}\\ \dot{\eta}\_{\varepsilon}\\ \mathbf{x}\_{L\_{\text{e}}}\end{bmatrix}}\_{\mathcal{X}\_{\mathcal{E}}} + \underbrace{\begin{bmatrix}0&0&0\\ \frac{B\_{\text{e}0}}{0}\underbrace{B\_{\text{el}}}\_{L\_{\text{el}}}&\mathbf{0}\end{bmatrix}}\_{\mathcal{B}\_{1}}\underbrace{\begin{Bmatrix}\delta\_{\text{u}}\\ \delta\_{\text{i}}\\ \delta\_{\text{u}}\end{Bmatrix}}\_{\delta\_{\text{II}}}\\ + \underbrace{\begin{bmatrix}0&0\\ \frac{A\_{\text{e}0}}{0}\underbrace{A\_{\text{e}0}}\_{L\_{\text{e}}}\end{bmatrix}}\_{\mathcal{S}\_{\text{Z}}} + \underbrace{\begin{bmatrix}0&0\\ \frac{B\_{\text{e}0}}{0}\left(0.5\mathbf{E}\_{\text{e}}\mathbf{u}\right)}\_{\mathcal{E}}\end{bmatrix}}\_{\mathcal{B}\_{3}}\end{aligned}} \tag{37}$$

where *δξas* and *δu* are defined as the incremental airframe states and inputs (perturbation from trim values),

$$
\delta\_{\tilde{\xi}\_{as}} = \tilde{\xi}\_{as} - \tilde{\xi}\_{as\vert\_0} \tag{38}
$$

$$
\delta\_{\boldsymbol{\mu}} = \boldsymbol{\mu} - \boldsymbol{\mu}\_{|\boldsymbol{0}} \tag{39}
$$

*ξas*| 0 and *u*|0 being the airframe state and input vectors computed at some specific trim condition *ab-initio* of the simulation run. Using a short notation form, Eq. (37) can be expressed as,

$$\dot{\mathbf{x}}\_{\ell} = A \,\mathbf{x}\_{\ell} + B\_1 \,\delta\_{\ell I} + B\_2 \,\delta\_{\tilde{\xi}} + B\_3 \,\mathbf{w} \tag{40}$$

14 Will-be-set-by-IN-TECH

*L*

Clearly to implement this algebraic equation, the generalized coordinate, *η*¯*e* = *ηe* − *ηe*0 , its rate, *η*˙*e*, and acceleration, *η*¨*e*, vectors as well as the aerodynamic lag terms related with the

Note that the Minimum-State method is formulated to only use a single set of lag states, *xL*, in Eq. (32). Therefore, the following augmented equation is devised to decouple the generalized

> � + �

> > *δu δu*˙ *δu*¨

In this way, only elastic lag terms are considered to avoid any possible coupling with the rigid-body airframe related states (i.e. *ξas* and *ξas*˙ ). Now, the following differential equation

> ⎧ ⎪⎪⎨

*ηe η*˙*e xLas xLe*

� �� � *xe*

⎡ ⎢ ⎢ ⎣ ⎫ ⎪⎪⎬

⎪⎪⎭

+

0 0 *Bew*1 *Bew*2 0 0.5*ErG* 0 0.5*ErG*

� �� � *B*3

⎡ ⎢ ⎢ ⎣

⎪⎪⎩

+

where *δξas* and *δu* are defined as the incremental airframe states and inputs (perturbation from

*δξas* = *ξas* − *ξas*|

0 and *u*|0 being the airframe state and input vectors computed at some specific trim condition *ab-initio* of the simulation run. Using a short notation form, Eq. (37) can be expressed

⎫ ⎬ ⎭ + �

� ⎧ ⎨ ⎩

�

*V η*˙*e* + *A*2*re*

*L*2

*EL*∗ *ELr* 0 *ELe* 0 00 *ELe*

> 0 0.5*ErG* 0 0.5*ErG*

*V*2 *η*¨*e* + *DrexLe*

*ξas*, related ones,

⎪⎪⎩

� � *wG w*˙ *G* �

000 *Be*0 *Be*1 *Be*2 0 0.5*ELδ* 0 0 0.5*ELδ* 0

� �� � *B*1

> � *wG w*˙ *G* �

� �� � *w*

*δu* = *u* − *u*|0 (39)

*x*˙*e* = *A xe* + *B*1 *δU* + *B*2 *δξ* + *B*3 *w* (40)

0 (38)

⎤ ⎥ ⎥ ⎦ *ξas ξas*˙ *ηe η*˙*e*

⎫ ⎪⎪⎬

⎪⎪⎭

⎤ ⎥ ⎥ ⎦

⎧ ⎨ ⎩

*δu δu*˙ *δu*¨

⎫ ⎬ ⎭ � �� � *δU*

� ⎧ ⎪⎪⎨ �

(35)

(36)

(37)

loads are computed similarly to the approximation given by Eq. (30),

elastic modes, *xLe* , must be estimated at each time iteration.

coordinates' aero lag terms from the airframe states, *ξas* and ˙

+ � *A*0*reη*¯*e* + *A*1*re*

� � *xLas xLe*

0 0.5*ELδ* 0 0 0.5*ELδ* 0

*D D* � � *xLas*

0 *I* 0 0 *Aee*0 *Aee*1 *AeL AeL*

0 *ELe* 0 *V*

� �� � *A*

> 0 0 *Aer*0 *Aer*1 *EL*∗ *ELr* 0 0

� �� � *B*2

*V c*¯ *R* 0

> ⎤ ⎥ ⎥ ⎦

0 *ELe*

*xLe*

is obtained combining the lower partition of Eq. (34) with the new devised Eq. (36).

⎤ ⎥ ⎥ ⎦

*c*¯ *R*

� *δξas δ* ˙*ξas* �

� �� � *δξ*

� Δ*F* Δ*M* � = *q*∞ �

� *x*˙*Las x*˙*Le*

� *zLas zLe* � = �

⎧ ⎪⎪⎨

*η*˙*e η*¨*e x*˙*Las x*˙*Le*

� �� � *x*˙*e*

⎫ ⎪⎪⎬

⎪⎪⎭

=

⎡ ⎢ ⎢ ⎣

+

⎡ ⎢ ⎢ ⎣

⎪⎪⎩

trim values),

*ξas*|

as,

� = � *V L R* 0 0 *V L R* The previous equation is used to estimate the elastic and lag states as a function of the incremental control input, (*δT U* = [ *δT u* , *δT u*˙ , *δT u*¨ ] *T* ), and incremental airframe states, (*δT ξ* = [ *δT ξas*, *δT* ˙ *ξas* ] *T* ) at each time iteration.

The quasi-static elastic deformation , *ηe*0, is computed by static residualization of the elastic modes, that is the *η*˙*e* = *η*¨*e* = *xLe* = 0 condition needs to be fulfilled. Therefore, the quasi-static elastic influence is estimated from Eq. (40) as:

$$\mathfrak{x}\_{\varepsilon} = 0 \Longrightarrow \mathfrak{x}\_{\varepsilon} = -A^{-1}(B\_1 \delta\_{\mathrm{II}} + B\_2 \delta\_{\mathrm{x}} i + B\_3 \,\mathrm{w}) \tag{41}$$

and from *xe* the quasi-static elastic influence vector, *ηe*0, can be recovered.

In summary, the linear aeroelastic solver will be built based on:

Fig. 7 illustrates the interconnection of the F/A-18 AAW 6-DOF Dynamics subsystem and the Incremental Aeroelastic Solver, Control Surface Mixer and Control Command Transform blocks.

Fig. 7. Addition of the Incremental Aeroelastic Loads Solver to the Nonlinear Rigid-Body 6-DOF Subsystem.

#### **6.2 Closed loop F/A-18 AAW linear model with gust excitation**

In order to demonstrate the proposed feedforward filter design algorithm, a simplified closed loop F/A-18 AAW linear simulink model with gust excitation is developed/implemented for the evaluation purposes. This high-fidelity aeroelastic model was developed using the following elements:

Gust Loads Alleviation 17

Adaptive Feedforward Control for Gust Loads Alleviation 111

u(t)

Basis function

In Out RLS Taps1

1

1

Gain -K-

Fig. 8.

*G*ˆ(*q*), is shown in Fig. 9.

Table 4.

[FilterF ]

T

e(t)

[error]

2 − *Nzcg*

PolyMAX Identification

6-DOF Linear FA 18 AAW Dynamics

error\_Nz

[FilterF ]

The error signal, *e*(*t*), can be selected as the vertical accelerometer reading at the aircraft left/right wing folder positions, i.e., *Nzkm*023*R* or *Nzkm*023*L*. An alternative choice could be

*Nzkm*023*R* + *Nzkm*023*L*

In this paper, Eq. (43) is served as a feedback signal for feedforward filter design purpose. The advantage of choosing Eq. (43) is that the rigid body dynamics can be partly removed, and

Upon initialization of the feedforward controller, a 20th order ORTFIR model, *G*ˆ(*q*), was estimated in order to create the filtered signal, *u*ˆ *f*(*t*). The amplitude bode plot of the estimated,

The modes used to build the orthonormal basis, *Bi*(*q*) are extracted from the stabilization diagram in Fig. 10. From Fig. 10, five elastic modes can be extracted, and they are shown in

> Mode Number Frequency [Hz] Damping *ζ* [%] 5.9246 4.5311 10.083 4.182 13.602 10.072 18.377 2.7409 21.569 2.5183

Table 4. Estimated Modes of Feedforward Filter Using FFT-PolyMAX Method.

 PROPRIETARY**

> lat\_stk\_cmmd Lat\_dyn Lon\_dyn Surf Deflection lon\_stk\_cmmd

Lateral Output

Lonitudinal Output

Flight variables Visualization

Generic Longitudinal Output 2

(43)

Generic Lateral Output 1

For continuous vertical gust perturbation, a low pass filter followed by a Dryden vertical velocity shaping filter is used to shape the power of the gust perturbation. The low pass filter is used to obtain the derivative of the gust perturbation. The low pass filter is given as *TLPF*(*s*) = *a s*+*a* where *a* = 200*π rad*/*s* is chosen in the remainder of the section.

The Dryden vertical velocity shaping filter is given as

$$T\_{\mathcal{S}}(\mathbf{s}) = \sigma\_{\mathcal{w}\_G} \frac{\sqrt{3}\tau\_{\mathcal{S}}^{-1/2}\mathbf{s} + \tau\_{\mathcal{S}}^{-3/2}}{(\mathbf{s} + 1/\tau\_{\mathcal{S}})^2}$$

where *τg* = *Lg*/*V*, and *Lg* = 1750 *f t*, *V* is the aircraft body axis velocity; *σwG* = 100 *f t*/*s*.

For a more detailed development of the F/A-18 AAW simulink model with gust excitation, please refer to the NASA SBIR Phase I final report (Zeng & de Callafon, 2008). The implementation of the adaptive feedforward control algorithm to the linearized F/A-18 AAW simulink model is illustrated in Fig. 8. It should be noted that during the simulation study considered in this paper, the dynamics of airborne LIDAR turbulence sensor has not been considered. We assumed that a perfect gust perturbation can be measured by the airborne LIDAR turbulence sensor, i.e., the sensor dynamics has an ideal constant dynamics of 1. However, the practical effects of the airborne LIDAR turbulence sensor on the performance of the feedforward controller has to be addressed in the future study.

#### **6.3 Implementation of the adaptive feedforward control**

The construction of the feedforward controller can be separated into two steps, initialization and the recursive estimation of the filter. In the initialization step, a secondary path transfer function, *G*(*q*), is estimated, which is done by performing an experiment using an external signal injected into the left and right trailing edge flaps as the excitation signal, and the error signal, *e*(*t*), as the output signal. Since, *G*ˆ(*q*), is only used for filtering purposes, a high order model can be estimated to provide an accurate reconstruction of the filtered input, *u*ˆ *f*(*t*), via

$$\text{if } \mathfrak{a}\_f(t) = \hat{G}(q)w(t) \tag{42}$$

as described in Eq. (4).

To facilitate the use of the ORTFIR filter, a set of modal parameters need to be extracted to build the ORTFIR filter, using the frequency domain PolyMAX identification methodology in Section 4. With *u*ˆ *f*(*t*) given in Eq. (42) and *d*(*t*) = *H*(*q*)*w*(*t*) in place, the modal parameters can be easily estimated using the PolyMAX method. With the signal, *d*(*t*), *u*ˆ *f*(*t*) and the basis function, *Bi*(*q*), a recursive minimization of the feedforward filter is done via the recursive least squares minimization technique described in Section 5.

16 Will-be-set-by-IN-TECH

• Aerodynamic Forces and Moments subsystem using the set of non-dimensional stability and control derivatives obtained through a set of AAW parameter identification flight tests. • An incremental aeroelastic load solver including gust excitation generated by the ZAERO

For continuous vertical gust perturbation, a low pass filter followed by a Dryden vertical velocity shaping filter is used to shape the power of the gust perturbation. The low pass filter is used to obtain the derivative of the gust perturbation. The low pass filter is given as

*s*+*a* where *a* = 200*π rad*/*s* is chosen in the remainder of the section.

where *τg* = *Lg*/*V*, and *Lg* = 1750 *f t*, *V* is the aircraft body axis velocity; *σwG* = 100 *f t*/*s*.

For a more detailed development of the F/A-18 AAW simulink model with gust excitation, please refer to the NASA SBIR Phase I final report (Zeng & de Callafon, 2008). The implementation of the adaptive feedforward control algorithm to the linearized F/A-18 AAW simulink model is illustrated in Fig. 8. It should be noted that during the simulation study considered in this paper, the dynamics of airborne LIDAR turbulence sensor has not been considered. We assumed that a perfect gust perturbation can be measured by the airborne LIDAR turbulence sensor, i.e., the sensor dynamics has an ideal constant dynamics of 1. However, the practical effects of the airborne LIDAR turbulence sensor on the performance of

The construction of the feedforward controller can be separated into two steps, initialization and the recursive estimation of the filter. In the initialization step, a secondary path transfer function, *G*(*q*), is estimated, which is done by performing an experiment using an external signal injected into the left and right trailing edge flaps as the excitation signal, and the error signal, *e*(*t*), as the output signal. Since, *G*ˆ(*q*), is only used for filtering purposes, a high order model can be estimated to provide an accurate reconstruction of the filtered input, *u*ˆ *f*(*t*), via

To facilitate the use of the ORTFIR filter, a set of modal parameters need to be extracted to build the ORTFIR filter, using the frequency domain PolyMAX identification methodology in Section 4. With *u*ˆ *f*(*t*) given in Eq. (42) and *d*(*t*) = *H*(*q*)*w*(*t*) in place, the modal parameters can be easily estimated using the PolyMAX method. With the signal, *d*(*t*), *u*ˆ *f*(*t*) and the basis function, *Bi*(*q*), a recursive minimization of the feedforward filter is done via the recursive

*u*ˆ *f*(*t*) = *G*ˆ(*q*)*w*(*t*) (42)

√3*τ*−1/2

*g s* + *τ*−3/2 *g*

(*s* + 1/*τg*)2

• Six-degree-of-freedom solver using Euler angles subsystem.

The Dryden vertical velocity shaping filter is given as

software system using rational function approximation techniques.

*Tg*(*s*) = *σwG*

the feedforward controller has to be addressed in the future study.

**6.3 Implementation of the adaptive feedforward control**

least squares minimization technique described in Section 5.

• The AAW flight control system.

• Actuators and sensors.

*TLPF*(*s*) = *a*

as described in Eq. (4).

Fig. 8. Closed Loop F/A-18 AAW Linear Simulink Model.

The error signal, *e*(*t*), can be selected as the vertical accelerometer reading at the aircraft left/right wing folder positions, i.e., *Nzkm*023*R* or *Nzkm*023*L*. An alternative choice could be

$$e(t) = \left[\frac{Nz\_{km023R} + Nz\_{km023L}}{2} - Nz\_{\mathcal{G}}\right] \tag{43}$$

In this paper, Eq. (43) is served as a feedback signal for feedforward filter design purpose. The advantage of choosing Eq. (43) is that the rigid body dynamics can be partly removed, and the vertical wing bending is still observed.

Upon initialization of the feedforward controller, a 20th order ORTFIR model, *G*ˆ(*q*), was estimated in order to create the filtered signal, *u*ˆ *f*(*t*). The amplitude bode plot of the estimated, *G*ˆ(*q*), is shown in Fig. 9.

The modes used to build the orthonormal basis, *Bi*(*q*) are extracted from the stabilization diagram in Fig. 10. From Fig. 10, five elastic modes can be extracted, and they are shown in Table 4.

Table 4. Estimated Modes of Feedforward Filter Using FFT-PolyMAX Method.

Gust Loads Alleviation 19

Adaptive Feedforward Control for Gust Loads Alleviation 113

Without Feedforward Control With Feedforward FIR With Feedforward ORTFIR

10−1 100 101 102 10−6

(a) Frequency Spectrum Plot of the *Nzkm*023*R* .

101 10−6

observed in regards to *Nzkm*023*L*, and which is shown in Fig. 12 (a) and (b).

Frequency (Hz)

Without Feedforward Control With Feedforward FIR With Feedforward ORTFIR

Frequency (Hz)

(b) Zoomed Frequency Spectrum Plot of the *Nzkm*023*R*.

zoomed-in plot of Fig. 11 (a) in the frequency range of [4 30] Hz. It is clearly seen that with the ORTFIR filter, a better magnitude reduction of auto spectrum of *Nzkm*023*R* can be obtained in most of the frequency range compared to the FIR filter. Similar performance could also be

The corresponding time responses are illustrated in Fig. 13 and Fig. 14. Fig. 13 (b) and Fig. 14 (b) are the zoomed-in plots of Fig. 13 (a) and Fig. 14 (a), respectively. These time responses clearly demonstrate that with the adaptive feedforward controller using ORTFIR filter, a better structural vibration reduction can be obtained. From these figures, it is clearly demonstrated

Fig. 11. Spectral Content Estimates of the *Nzkm*023*R* Without Control (Solid), With Control Using 20th Order FIR Filter (Dashed), and Using 20th Order ORTFIR Filter (Dotted).

10−4

10−4

10−2

Magnitude [Db]

100

10−2

Magnitude [Db]

100

Fig. 9. Bode Plot of the Estimated 20th Order Secondary Path Model, *G*ˆ(*q*).

Fig. 10. Stabilization Diagram.

For implementation purposes, only *L* = 2 parameters in the ORTFIR filter are estimated. With a 10th order basis, *Bi*(*q*), this amounts to 20th order ORTFIR filter. To evaluate the performance of the proposed ORTFIR filter for feedforward compensation, a 20th order Finite Impulse Response Filter is also designed to reduce the vertical wing vibration.

For a clear performance comparison between FIR filter and ORTFIR filter, the frequency response of the *Nzkm*023*R*, and *Nzkm*023*L* using FIR filter and ORTFIR filter are plotted in Fig. 11 and Fig. 12, respectively. The solid line in Fig. 11 (a) is the auto spectrum of the accelerometer measurement *Nzkm*023*R* without feedforward controller integrated in the system; the dashed line in Fig. 11(a) indicates the auto spectrum of the accelerometer measurement *Nzkm*023*R* with the adaptive feedforward controller using FIR filter added in the system; the dotted line in Fig. 11 shows the auto spectrum of the accelerometer measurement *Nzkm*023*R* with the adaptive feedforward controller using ORTFIR filter added in the system. Fig. 11 (b) is the

18 Will-be-set-by-IN-TECH

Frequency Transfer Function Estimated Second Path Model

0 10 20 30 40 50 60 70 80 90 100 −40

0 10 20 30 40 50 60 70 80 90 100 −200

Fig. 9. Bode Plot of the Estimated 20th Order Secondary Path Model, *G*ˆ(*q*).

Frequency [Hz]

0 10 20 30 40 50 40

For implementation purposes, only *L* = 2 parameters in the ORTFIR filter are estimated. With a 10th order basis, *Bi*(*q*), this amounts to 20th order ORTFIR filter. To evaluate the performance of the proposed ORTFIR filter for feedforward compensation, a 20th order Finite

For a clear performance comparison between FIR filter and ORTFIR filter, the frequency response of the *Nzkm*023*R*, and *Nzkm*023*L* using FIR filter and ORTFIR filter are plotted in Fig. 11 and Fig. 12, respectively. The solid line in Fig. 11 (a) is the auto spectrum of the accelerometer measurement *Nzkm*023*R* without feedforward controller integrated in the system; the dashed line in Fig. 11(a) indicates the auto spectrum of the accelerometer measurement *Nzkm*023*R* with the adaptive feedforward controller using FIR filter added in the system; the dotted line in Fig. 11 shows the auto spectrum of the accelerometer measurement *Nzkm*023*R* with the adaptive feedforward controller using ORTFIR filter added in the system. Fig. 11 (b) is the

Impulse Response Filter is also designed to reduce the vertical wing vibration.

Frequency [Hz]

−69

−52.4

−35.8

−19.2

Magnitude [dB]

−2.6

14

Frequency [Hz]

50

Model Order

Fig. 10. Stabilization Diagram.

60

Phase [deg]

Magnitude [Db]

(a) Frequency Spectrum Plot of the *Nzkm*023*R* .

Fig. 11. Spectral Content Estimates of the *Nzkm*023*R* Without Control (Solid), With Control Using 20th Order FIR Filter (Dashed), and Using 20th Order ORTFIR Filter (Dotted).

zoomed-in plot of Fig. 11 (a) in the frequency range of [4 30] Hz. It is clearly seen that with the ORTFIR filter, a better magnitude reduction of auto spectrum of *Nzkm*023*R* can be obtained in most of the frequency range compared to the FIR filter. Similar performance could also be observed in regards to *Nzkm*023*L*, and which is shown in Fig. 12 (a) and (b).

The corresponding time responses are illustrated in Fig. 13 and Fig. 14. Fig. 13 (b) and Fig. 14 (b) are the zoomed-in plots of Fig. 13 (a) and Fig. 14 (a), respectively. These time responses clearly demonstrate that with the adaptive feedforward controller using ORTFIR filter, a better structural vibration reduction can be obtained. From these figures, it is clearly demonstrated

Gust Loads Alleviation 21

Adaptive Feedforward Control for Gust Loads Alleviation 115

Without Feedforward Control With Feedforward FIR With Feedforward ORTFIR

0 10 20 30

Time (s)

Without Feedforward Control With Feedforward FIR With Feedforward ORTFIR

9 9.2 9.4 9.6 9.8 10

Time (s)

(b) Zoomed Time Domain Response Plot of the *Nzkm*023*L*.

Fig. 13. Time domain Response of the *Nzkm*023*R* Without Control (Solid), With Control Using

20th Order FIR Filter (Dashed), and Using 20th Order ORTFIR Filter (Dotted).

(a) Time Domain Response Plot of the *Nzkm*023*L* .

−2

5

−5

0

N

(g)

Z

KM023R

−1

0

N

(g)

Z

KM023R

1

2

3

(b) Zoomed Frequency Spectrum Plot of the *Nzkm*023*L*.

Fig. 12. Spectral Content Estimates of the *Nzkm*023*L* Without Control (Solid), With Control Using 20th Order FIR Filter (Dashed), and Using 20th Order ORTFIR Filter (Dotted).

that both FIR filter and ORTFIR filter are efficient to reduce the normal acceleration at the left wing folder position and right wing folder position. With the use of the both the ORTFIR filter and FIR filter, the spectral content of the *Nzkm*023*R* and *Nzkm*023*L* have been reduced significantly in the frequency range from 2Hz to 20Hz. However, with the use of the ORTFIR filter, more efficient vibration reduction performances are expected compared to the FIR filter. 20 Will-be-set-by-IN-TECH

Without Feedforward Control

With Feedforward FIR With Feedforward ORTFIR

10−1 100 101 102 10−6

(a) Frequency Spectrum Plot of the *Nzkm*023*L* .

101 10−6

Frequency (Hz)

Without Feedforward Control

With Feedforward FIR With Feedforward ORTFIR

Frequency (Hz)

(b) Zoomed Frequency Spectrum Plot of the *Nzkm*023*L*.

that both FIR filter and ORTFIR filter are efficient to reduce the normal acceleration at the left wing folder position and right wing folder position. With the use of the both the ORTFIR filter and FIR filter, the spectral content of the *Nzkm*023*R* and *Nzkm*023*L* have been reduced significantly in the frequency range from 2Hz to 20Hz. However, with the use of the ORTFIR filter, more efficient vibration reduction performances are expected compared to the FIR filter.

Fig. 12. Spectral Content Estimates of the *Nzkm*023*L* Without Control (Solid), With Control Using 20th Order FIR Filter (Dashed), and Using 20th Order ORTFIR Filter (Dotted).

10−4

10−4

10−2

Magnitude [Db]

100

10−2

Magnitude [Db]

100

(b) Zoomed Time Domain Response Plot of the *Nzkm*023*L*.

Fig. 13. Time domain Response of the *Nzkm*023*R* Without Control (Solid), With Control Using 20th Order FIR Filter (Dashed), and Using 20th Order ORTFIR Filter (Dotted).

Gust Loads Alleviation 23

Adaptive Feedforward Control for Gust Loads Alleviation 117

In this chapter, an adaptive feedforward control methodology has been proposed for the active control of gust loads alleviation using an ORTFIR filter. The ORTFIR filter has the same linear parameter structure as a taped delay FIR filter that is favorable for (recursive) estimation. The advantage of using the ORTFIR filter is that it allows the inclusion of prior knowledge of the flexible mode information of the aircraft dynamics in the parametrization of the filter for better

In addition, by combining the flight dynamics model for rigid body dynamics and an aeroelastic solver for aeroelastic incremental loads to accurately mimic in-flight recorded dynamic behavior of the air vehicle, a unified integration framework that blends flight

The proposed methodology in this chapter is implemented on an F/A-18 AAW aeroelastic model developed with the unified integration framework. The feedforward filter is updated via the recursive least square technique with the variable forgetting factor at each time step. Compared with a traditional FIR filter and evaluated on the basis of the simulation data from the F/A-18 AAW aeroelastic model, it demonstrates that applying the adaptive feedforward controller using the ORTFIR filter yields a better performance of the gust loads alleviation of

Andrighettoni, M. & Mantegazza, P. (1998). Multi-Input/Multi-output adaptive active flutter

Baldelli, D. H., Chen, P. C. & Panza, J. (2006). Unified aeroelastic and flight dynamics formulation via rational function approximations, *Journal of Aircraft* 43(3): 763–772. Baldelli, D. H. & Zeng, J. (2007). Unified nonlinear flight dynamics aeroelastic solver tool,

Baldelli, D. H., Zeng, J., Lind, R. & Harris, C. (2009). Flutter-prediction tool for flight-test-based

Barker, J. M., Balas, G. J. & Blue, P. A. (1999). Active flutter suppression via gain-scheduled

Cauberghe, B., Guillaume, P., Verboven, P., Parloo, E. & Vanlanduit, S. (2004). A poly-reference

Cauberghe, B., Guillaume, P., Verboven, P., Vanlanduit, S. & Parloo, E. (2005). On the influence

Eversman, W. & Roy, I. D. (1996). Active flutter suppresion using MIMO adaptive LMS

Garrard, W. L. & Liebst, B. S. (1985). Active flutter suppression using eigenspace and linear

*Material Conference and Exhibit*, Salt Lake, UT. AIAA-1996-1345.

Haykin, S. (2002). *Adaptive Filter Theory*, Prentice Hall, Englewood Cliffs, NJ.

aeroelastic parameter-varying models, *Journal of Guidance, Control, and Dynamics*

linear fractional control, *Proceedings of the American Control Conference*, San Diego,

implementation of the maximum likelihood complex frequency-domain estimator and some industrial applications, *Proceedings of the 22nd International Modal Analysis*

of the parameter constraint on the stability of the poles and the discrimination capabilities of the stabilization diagrams, *Mechanical Systems and Signal Processing*

control, *37th AIAA / ASME / ASCE / AHS / ASC Structures, Structural Dynamics, and*

quadratic design techniques, *Journal of Guidance, Control, and Dynamics* 8(3): 304–311.

suppression for a wing model, *Journal of Aircraft* 35(3): 462–469.

*Technical Report SBIR Phase I NNL07AA85P*, ZONA Technology Inc.

dynamics and aeroelastic model is developed to facilitate the pre-flight simulation.

**7. Conclusions**

the aircraft.

**8. References**

32(1): 158–171.

19(5): 989–1014.

California, pp. 4014–4018.

*Conference*, Dearborn, US.

accuracy of the feedforward filter.

(b) Zoomed Time Domain Response Plot of the *Nzkm*023*L*.

Fig. 14. Time Domain Response of the *Nzkm*023*L* Without Control (Solid) and Using 20th Order FIR Filter (Dashed), and With Control Using 20th Order ORTFIR Filter (Dotted).

## **7. Conclusions**

22 Will-be-set-by-IN-TECH

Without Feedforward Control

With Feedforward FIR With Feedforward ORTFIR

0 10 20 30

Time (s)

Without Feedforward Control

With Feedforward FIR With Feedforward ORTFIR

9 9.2 9.4 9.6 9.8 10

Time (s)

(b) Zoomed Time Domain Response Plot of the *Nzkm*023*L*.

Fig. 14. Time Domain Response of the *Nzkm*023*L* Without Control (Solid) and Using 20th Order FIR Filter (Dashed), and With Control Using 20th Order ORTFIR Filter (Dotted).

(a) Time Domain Response Plot of the *Nzkm*023*L* .

−2

5

−5

0

N

(g)

Z

KM023L

−1

0

N

(g)

Z

KM023L

1

2

3

In this chapter, an adaptive feedforward control methodology has been proposed for the active control of gust loads alleviation using an ORTFIR filter. The ORTFIR filter has the same linear parameter structure as a taped delay FIR filter that is favorable for (recursive) estimation. The advantage of using the ORTFIR filter is that it allows the inclusion of prior knowledge of the flexible mode information of the aircraft dynamics in the parametrization of the filter for better accuracy of the feedforward filter.

In addition, by combining the flight dynamics model for rigid body dynamics and an aeroelastic solver for aeroelastic incremental loads to accurately mimic in-flight recorded dynamic behavior of the air vehicle, a unified integration framework that blends flight dynamics and aeroelastic model is developed to facilitate the pre-flight simulation.

The proposed methodology in this chapter is implemented on an F/A-18 AAW aeroelastic model developed with the unified integration framework. The feedforward filter is updated via the recursive least square technique with the variable forgetting factor at each time step. Compared with a traditional FIR filter and evaluated on the basis of the simulation data from the F/A-18 AAW aeroelastic model, it demonstrates that applying the adaptive feedforward controller using the ORTFIR filter yields a better performance of the gust loads alleviation of the aircraft.

## **8. References**

**0**

**5**

*Canada*

**Fault Tolerant Flight Control Techniques with**

**Application to a Quadrotor UAV Testbed**

*Department of Mechanical and Industrial Engineering, Concordia University*

Unmanned Aerial Vehicles (UAVs) are gaining more and more attention during the last few years due to their important contributions and cost-effective applications in several tasks such as surveillance, search and rescue missions, geographic studies, as well as various military and security applications. Due to the requirements of autonomous flight under different flight conditions without a pilot onboard, control of UAV flight is much more challenging compared with manned aerial vehicles since all operations have to be carried out by the automated flight control, navigation and guidance algorithms embedded on the onboard flight microcomputer/microcontroller or with limited interference by a ground pilot

As an example of UAV systems, the quadrotor helicopter is relatively a simple, affordable and easy to fly system and thus it has been widely used to develop, implement and test-fly methods in control, fault diagnosis, fault tolerant control as well as multi-agent based technologies in formation flight, cooperative control, distributed control, mobile wireless networks and communications. Some theoretical works consider the problems of control (Dierks & Jagannathan, 2008), formation flight (Dierks & Jagannathan, 2009) and fault diagnosis (Nguyen et al., 2009; Rafaralahy et al., 2008) of the quadrotor UAV. However, few research laboratories are carrying out advanced theoretical and experimental works on the system. Among others, one may cite for example, the UAV SWARM health management project of the Aerospace Controls Laboratory at MIT (SWARM, 2011), the Stanford Testbed of Autonomous Rotorcraft for Multi-Agent Control (STARMAC) project (STARMAC, 2011) and

A team of researchers at the Department of Mechanical and Industrial Engineering of Concordia University, with the financial support from NSERC (Natural Sciences and Engineering Research Council of Canada) through a Strategic Project Grant and a Discovery Project Grant and three Canadian-based industrial partners (Quanser Inc., Opal-RT Technologies Inc., and Numerica Technologies Inc.) as well as Defence Research and Development Canada (DRDC) and Laval University, have been working on a research and development project on fault-tolerant and cooperative control of multiple UAVs since 2007 (see the Networked Autonomous Vehicles laboratory (NAV, 2011)). In addition to the work that has been carried out for the multi-vehicles case, many fault-tolerant control (FTC) strategies have been developed and applied to the single vehicle quadrotor helicopter UAV system. The objective is to consider actuator faults and to propose FTC methods to

the Micro Autonomous Systems Technologies (MAST) project (MAST, 2011).

**1. Introduction**

if needed.

Youmin Zhang and Abbas Chamseddine