April 2006 

NASA/TM-2006-214296 

 

Exploratory Studies in Generalized Predictive 

Control for Active Gust Load Alleviation 

Raymond G. Kvaternik, Kenneth W. Eure, and Jer-Nan Juang 

Langley Research Center, Hampton, Virginia 
The NASA STI Program Office ... in Profile 

Since its founding, NASA has been dedicated 

to the advancement of aeronautics and space 

science. The NASA Scientific and Technical 

Information (STI) Program Office plays a key 

part in helping NASA maintain this important 

role. 

The NASA STI Program Office is operated by 

Langley Research Center, the lead center for 

NASA’s scientific and technical information. 

The NASA STI Program Office provides 

access to the NASA STI Database, the largest 

collection of aeronautical and space science 

STI in the world. The Program Office is also 

NASA’s institutional mechanism for 

disseminating the results of its research and 

development activities. These results are 

published by NASA in the NASA STI Report 

Series, which includes the following report 

types: 

• TECHNICAL PUBLICATION. Reports 

of completed research or a major 

significant phase of research that 

present the results of NASA programs 

and include extensive data or theoretical 

analysis. Includes compilations of 

significant scientific and technical data 

and information deemed to be of 

continuing reference value. NASA 

counterpart of peer-reviewed formal 

professional papers, but having less 

stringent limitations on manuscript 

length and extent of graphic 

presentations. 

• TECHNICAL MEMORANDUM. 

Scientific and technical findings that are 

preliminary or of specialized interest, 

e.g., quick release reports, working 

papers, and bibliographies that contain 

minimal annotation. Does not contain 

extensive analysis. 

• CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees. 

• CONFERENCE PUBLICATION. 

Collected papers from scientific and 

technical conferences, symposia, 

seminars, or other meetings sponsored 

or co-sponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and missions, 

often concerned with subjects having 

substantial public interest. 

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign 

scientific and technical material 

pertinent to NASA’s mission. 

Specialized services that complement the 

STI Program Office’s diverse offerings 

include creating custom thesauri, building 

customized databases, organizing and 

publishing research results ... even 

providing videos. 

For more information about the NASA STI 

Program Office, see the following: 

• Access the NASA STI Program Home 

Page at http://www.sti.nasa.gov

• E-mail your question via the Internet to 

help@sti.nasa.gov 

• Fax your question to the NASA STI 

Help Desk at (301) 621-0134 

• Phone the NASA STI Help Desk at (301) 

621-0390 

• Write to: 

 NASA STI Help Desk 

 NASA Center for AeroSpace Information 

 7121 Standard Drive 

 Hanover, MD 21076-1320
National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23681-2199

 

April 2006 

NASA/TM-2006-214296 

 

Exploratory Studies in Generalized Predictive 

Control for Active Gust Load Alleviation 

Raymond G. Kvaternik, Kenneth W. Eure, and Jer-Nan Juang 

Langley Research Center, Hampton, Virginia 
Available from: 

NASA Center for AeroSpace Information (CASI) National Technical Information Service (NTIS) 

7121 Standard Drive 5285 Port Royal Road 

Hanover, MD 21076-1320 Springfield, VA 22161-2171 

(301) 621-0390 (703) 605-6000 

The use of trademarks or names of manufacturers in this report is for accurate reporting and does not constitute an 

official endorsement, either expressed or implied, of such products or manufacturers by the National Aeronautics 

and Space Administration. 
iii 

Table of Contents 

Abstract............................................................................................................................................1 

Nomenclature...................................................................................................................................1 

Introduction......................................................................................................................................2 

Basic Features of GPC.....................................................................................................................3 

Pertinent Equations Underlying the GPC Method...........................................................................4 

 Form of System and Control Law Equations..........................................................................5 

 System Identification ..............................................................................................................5 

 Multi-Step Output Prediction Equation ..................................................................................7 

 Derivation of Control Law......................................................................................................9 

Selection of Control Model Parameters.........................................................................................10 

Practical Issues...............................................................................................................................10 

 Treatment of External Disturbances .....................................................................................11 

 Basic Signal Processing Guidelines......................................................................................11 

Results and Discussion ..................................................................................................................12 

 Comment on Codes Used for Simulation .............................................................................12 

 Numerical Simulations Using 3-DOF Math Model..............................................................12 

 Representative Results.................................................................................................14 

 Steady-State Disturbances ...................................................................................14 

 Non-Steady-State Disturbances ...........................................................................15 

 Summary of Observations....................................................................................15 

 Numerical Simulations Using Data from Wind-Tunnel Test of HiLDA Wing Model ........16 

 Test Conditions ............................................................................................................16 

 Test Procedures............................................................................................................17 

 Simulation Conditions and Parameter Settings............................................................17 

 Representative Results.................................................................................................18 

 Feedback with Embedded Feedforward...............................................................19 

 Feedback with Explicit Feedforward...................................................................20 

 Summary of Observations....................................................................................20 

Concluding Remarks......................................................................................................................21 

Appendix A ⎯ Treatment of Measurable Coherent Disturbances................................................22 

 Batch Identification of Disturbance Models........................................................................23 

 Adaptive Identification of Disturbance Models...................................................................23 

Appendix B ⎯ State-Space Equations for 3-DOF Mass-Spring-Dashpot System .......................34 

References......................................................................................................................................37 

Tables.............................................................................................................................................39 

Figures............................................................................................................................................40

1 

Abstract 

The results of numerical simulations aimed at assessing the efficacy of Generalized 

Predictive Control (GPC) for active gust load alleviation using trailing- and leadingedge control surfaces are presented. GPC is a linear, time-invariant, multi-input/multioutput predictive control method that uses an AutoRegressive with eXogenous input 

(ARX) model to characterize the system and to design the controller. The equations 

underlying the method are presented and discussed, including system identification, 

calculation of control law matrices, and calculation of the commands applied to the 

control effectors. Both embedded and explicit feedforward paths for incorporation of 

disturbance effects are addressed. Representative results from two types of simulations 

are presented. The first used a three-degree-of-freedom mathematical model of a massspring-dashpot system subject to a variety of user-defined external disturbances. The 

second used open-loop data obtained during a wind-tunnel test in which a wing model 

was subjected to sinusoidal vertical gust excitation; closed-loop behavior was simulated 

in post-test calculations. Results obtained from these simulations have been decidedly 

positive. In particular, the closed-loop simulations for the wing model showed reductions 

in wing root moments by factors as high as 1000, depending on whether the excitation is 

from a constant- or swept-frequency gust and on the direction of the response. 

Nomenclature 

AR AutoRegressive 

ARX AutoRegressive with eXogenous input 

d disturbance acting on plant 

dk, du known and unknown portions of disturbance 

GPC Generalized Predictive Control 

hp, hc prediction and control horizons 

k integer denoting time index 

A number of data points used for system identification 

m number of outputs 

nd order of AR disturbance model 

OMP Observer Markov Parameters 

p order of ARX model used for system identification 
2 

r total number of inputs (r = rc + rd) 

rc number of control inputs

rd number of disturbance inputs 

SID System IDentification 

T, A, B, Dp, Df coefficient matrices in multi-step output prediction equation 

u input 

uc control input 

uc(k) control input at time step k

u(k), y(k), d(k) vectors of inputs, outputs, and disturbances at k

th time step 

u(k-j), y(k-j), d(k-j) vectors of past inputs, outputs, and disturbances at (k-j)

th time step; j < k

u(k+j), y(k+j), d(k+j) vectors of future inputs, outputs, and disturbances at (k+j)

th time step; j ≥ 0 

y output 

Y matrix of observer Markov parameters 

χ

c

, αc, βc, δc

, γ

c

 control law gain matrices 

αi, βi and γi observer Markov parameters 

Introduction 

The Aeroelasticity Branch and the Army Research Laboratory’s Vehicle Technology 

Directorate at NASA Langley Research Center, in collaboration with Bell Helicopter Textron 

Inc, have recently completed a series of experimental evaluations of a predictive control 

technique known as Generalized Predictive Control (GPC) (ref. 1) to assess its potential for 

actively controlling the swashplate of tiltrotor aircraft to enhance aeroelastic stability in both 

helicopter and airplane modes of flight. The studies employed a 1/5-scale, semispan aeroelastic 

model of the V-22 that had been refurbished to serve as a tiltrotor research testbed for use in the 

Langley Transonic Dynamics Tunnel (TDT). The results of the exploratory experimental 

investigations conducted with this model using a multi-input/multi-output GPC-based method for 

active control are reported in references 2-5. The GPC algorithm employed was highly effective 

in increasing the stability of the critical (least stable) mode for all of the configurations and 

conditions tested. GPC has apparently not yet been employed to evaluate its potential for active 
3 

gust load alleviation (GLA) of fixed-wing aircraft using wing control surfaces. To this end, the 

intention was to upgrade and extend the GPC-based software system that was developed for the 

tiltrotor model and to conduct an investigation of the GPC method for active gust load alleviation 

as an adjunct to a planned test in the TDT of a large semi-span wing model of a SensorCraft 

concept being studied by Northrop Grumman as part of a High Lift over Drag Active (HiLDA) 

Wing Program (ref. 6). Although it was not possible to make all the necessary changes to the 

software in time for the subject test, appropriate data were taken during the test with the 

objective of conducting a post-test assessment and evaluation of the method. 

 

The purpose of this report is to present results of numerical simulations obtained from 

recently-developed MATLAB (ref. 7) programs intended to serve as the basis for the planned 

extensions to the tiltrotor GPC software system. First, a qualitative discussion of the essential 

features of GPC is given. Then, the equations underlying the method are presented and 

discussed, including system identification, calculation of control law matrices, and calculation of 

the commands applied to the control effectors. Procedures for incorporation of both known and 

unknown disturbance effects are addressed. Some general but relevant comments on the 

treatment of external disturbances, signal processing, and the codes developed and used for the 

simulations are then presented. Representative results from two types of simulations of openand closed-loop behavior are presented next. The first set of simulations was made using a 

three-degree-of-freedom mathematical model of a mass-spring-dashpot system subject to a 

variety of user-defined external disturbances. The second set of simulations was made using 

data obtained during the HiLDA wing test conducted in the TDT in September 2004 (ref. 6) in 

which the model was subjected to sinusoidal vertical gust excitation at constant and continuously 

varying frequencies. The particular results presented from each type of simulation were chosen 

to illustrate the potential gust-load-alleviation capability that can be obtained with a GPC-based 

active control system. 

Basic Features of GPC 

Predictive control refers to a strategy wherein the decision for the current control action is 

based on minimization of a quadratic objective function that involves a prediction of the system 

response at some number of time steps into the future. A variety of predictive controllers have 

been proposed (ref. 8). Among these, Generalized Predictive Control, which was introduced in 

1987 (ref. 9), has received notable attention by researchers. GPC is a linear, time-invariant, 

multi-input/multi-output predictive control method that uses an AutoRegressive with eXogenous 

input (ARX) model (i.e., a linear difference equation) to describe the input-output relationship of 

the system. The coefficient matrices of the ARX equation are determined by system 

identification. Such a model gives the current outputs as a linear combination of past inputs and 

outputs. The input-output equation is used to form a multi-step output prediction equation over a 

finite prediction horizon. An expression for the control to be imposed at the next time step is 

determined by minimizing the deviation of the predicted controlled plant outputs from the 

desired (or target) outputs, subject to a penalty on control effort. A version of the GPC 

procedure was developed at NASA Langley Research Center in 1997 for efficient computation 

and unknown disturbance rejection by Dr. Jer-Nan Juang and his associates. Their work resulted 

in a suite of MATLAB m-files that were collected into a predictive control toolbox that can be 
4 

used by researchers for GPC studies. A summary of the theory underlying their development is 

found in references 1 and 10-20, among others. 

The essential features of a GPC-based control system are depicted in figure 1. The system 

has rc control inputs u, m measured outputs y, and is subject to external disturbances d that may 

be regarded as having known (dk) and unknown (du) portions. Measurement noise (not shown) is 

also present. There are two fundamental steps involved in GPC implementation: (1) 

identification of the system, and (2) use of the identified model to design a controller. Any 

system identification technique that returns an ARX model for the system may be used. System 

identification is done on-line, but not necessarily at every time point, in the presence of the 

disturbances acting on the system. If all the external disturbances are unknown (or treated as 

unknown), the effects of the disturbance-induced output are automatically embedded in the 

identified model and there is no need for explicit inclusion of disturbance terms in the ARX 

model (ref. 1). This approach leads to a feedback controller with embedded (or implicit) 

feedforward. However, if any of the disturbances acting on the system are either known or 

measurable, that data can be used in a feedforward path simultaneously with the feedback data to 

enhance closed-loop performance, as suggested in figure 1. In this situation, the identified ARX 

model includes the effects of the disturbance-induced output associated with the unknown 

disturbances while an AutoRegressive (AR) model identified on-line can be used to include the 

effects of the measurable disturbances. This approach leads to a feedback controller with both 

embedded and explicit feedforward. In either case, the parameters of the identified model are 

used to compute the predictive control law. A random excitation uid is applied initially with uc

equal to zero to identify the open-loop system. Dither is added to the closed-loop control input

uc if it is necessary to re-identify the system while operating in the closed-loop mode. It should 

be pointed out that this method lends itself to implementation strategies that allow all 

computations to be done on-line and in real time, as was demonstrated in experimental 

investigations with the tiltrotor model. 

The GPC software system developed for the tiltrotor model identifies the plant and 

computes the control law gain matrices assuming that the aerodynamic forces (the disturbances) 

acting on the model are unknown. This is the case of feedback with embedded feedforward 

described above and represents the capability that was incorporated in the m-files of the 

aforementioned predictive control toolbox. A Langley team has recently identified, implemented 

(as MATLAB codes), and verified the changes needed to include a feedforward path in the GPCbased software system developed for the tiltrotor model. This extension was motivated by the 

fact that the data needed for a feedforward path would be available from the HiLDA wing test in 

the form of gust-induced flow-angle measurements made upstream of the model. It is not 

difficult to envision other active control situations in which one or more of the disturbances 

acting on a system are measurable (or known) and which could take advantage of feedforward 

control. 

Pertinent Equations Underlying the GPC Method 

 The key ingredients of the GPC method are system identification and output prediction 

based on the identified model; the control law is then derived by minimizing the error between 
5 

the predicted output and the desired output. The section begins by showing the general form of 

the system and control law equations. 

Form of System and Control Law Equations 

The relationship between the input and output discrete time histories of a linear, time 

invariant, multi-input/multi-output (MIMO) system is described by an ARX model that has the 

form 

 1 2 01

0 1

( ) ( -1) + ( -2) + + ( - ) + ( ) + ( -1) + + ( - )

 + ( ) + ( -1) + + ( - )

p p

p

yk yk yk ykp uk uk ukp

dk dk dkp

α α α ββ β

γγ γ

= " "

" (1)

This equation states that the current output y(k) at time step k may be estimated by using p sets of 

the previous output, input, and disturbance measurements, y(k-1),..., y(k-p), u(k-1),...,u(k-p), and 

d(k-1),...d(k-p); and the current input and disturbance measurements u(k) and d(k). The integer p

is called the order of the ARX model. The coefficient matrices αi, βi and γi are referred to as 

observer Markov parameters (OMP) or ARX parameters and are the quantities to be determined 

by the identification algorithm. Closed-loop performance is enhanced by performing the system 

identification in the presence of the external disturbances acting on the system, thereby ensuring 

that disturbance information will be incorporated in the system model. The goal of system 

identification is to determine the OMP based on input and output data. The OMP may be 

determined by any identification technique that returns an ARX model of the system. 

The ARX model is used to design the controller that, in the case of a regulator problem, 

has the general form given by 

 

c1 2 1 2

1 2 12

( ) ( -1) + ( -2) + + ( - ) + ( -1) + ( -2) + + ( - )

 + d( -1) + d( -2) + + d( - ) + d( ) + d( 1) + + d( 1)

c c c cc c

p p

c c c cc c

p p

u k yk yk ykp uk uk ukp

k k kp k k k p

α α α ββ β

δ δ δ γγ γ

=

+ +−

" "

" " (2)

Equation 2 indicates that the current control input uc(k) may be computed using p sets of the 

previous input and output measurements, u(k-1),..., u(k-p), d(k-1),..., d(k-p), and y(k-1),..., y(kp); and the future disturbances d(k), d(k+1), ..., d(k+p-1). The coefficient matrices αc

, βc, δc

, and 

γ

c

 are the control law gain matrices. If the p past disturbances d(k-1), d(k-1),..., d(k-p) are 

measurable (known) they can be inserted directly into equation 2. The future disturbances, d(k), 

d(k+1), ..., d(k+p-1), are typically unknown. However, if there is correlation in the disturbance 

signal (such as in the sinusoidal gust excitation of the HiLDA wing model) an AutoRegressive 

(AR) model of the disturbance can be determined and used to predict the future disturbances 

needed in equation 2. 

System Identification 

System identification (SID) in the presence of the operational disturbances is the first of 

the two major computational steps of the GPC method. The number of control inputs is rc, the 

number of disturbance inputs is rd, the number of measured outputs is m, and r = rc + rd. To 

initiate the identification process, the system is excited with band-limited white noise. These 

independent random excitations are applied to the rc inputs simultaneously and the m responses 
6 

are measured. The resulting input and output time histories (u and y) at A time points are then 

used to form the data matrices y and V in the equation 

 y = Y V (3)

where 

 (0) (1) (2) ( ) ( 1)

m

y y y y yp y ⎡ ⎤ ⎢ ⎥

⎣ ⎦

×

= −

A

" "A (4)

and 

[ ( )]

(0) (1) (2) ( ) ( 1)

(0) (1) (2) ( ) ( 1)

(0) (1) ( 1) ( 2) V=

(0) ( 2) ( 3)

(0) ( 1)

r r mp

d d d dp d

u u u up u

v v vp v

v vp v

v vp

⎡ ⎤ ⎢ ⎥

⎣ ⎦

++ ×

−

−

− −

− −

− −

A

" " A

A

" "A

" "A

%#" #

" A

 (5)

Equations 4 and 5 follow from writing the equations for a discrete-time state-space observer 

model of a linear time-invariant system at a sequence of time steps k = 0, 1, 2, ..., ( A -1) and 

grouping (assembling) them into a matrix form (ref. 11). Equations 4 and 5 are valid under the 

assumption that the system is initially at rest. For nonzero unknown initial conditions, the first p

columns of y and V must be deleted to eliminate the effect of the initial conditions (ref. 11). The 

vector ν(k) appearing in the data matrix V is formed from the vectors d(k), u(k) and y(k) 

according to 

 

( )1

( )

( ) ( )

( )

r m

d k

v k u k

y k

⎧ ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎭

+ ×

= (6) 

The order of the ARX model, p, and the number of time steps, A , are specified by the user. 

Some guidelines for their selection are given later. The sizes of key vectors and matrices are 

noted here and throughout the report as appropriate to facilitate computational considerations 

related to computer implementation. 

In forming the matrices given in equations 4 and 5, it has been assumed that a deadbeat 

observer has been added to the system (ref. 18). It is through this assumption that the matrix V

can be reduced to a size amenable for practical numerical computation of its pseudo-inverse. 

The system identification process yields observer Markov parameters rather than system Markov 

parameters because of the inclusion of the observer. A complete discussion of these aspects of 

the development may be found in reference 11. 

 In equation 3, Y is the matrix of observer Markov parameters that is to be identified and 

has the form 
7 

 

[ ( )]

0 01 1 1 2 2 23 3 3 p

m r r mp

Y p p γ β γ βαγ βαγ βα γ βα ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦

×++

= " (7)

The solution for Y is obtained by solving equation 3 for Y according to 

 -1 † y y = T T Y V V VV = ⎡ ⎤ ⎣ ⎦ (8) 

where † denotes the pseudo-inverse. If the product V VT is a well-conditioned matrix of 

reasonable size, the ordinary inverse can be taken as shown. Otherwise, a pseudo-inverse must 

be used. 

Multi-Step Output Prediction Equation 

 

 The one-step-ahead output prediction equation given in equation 1 is the starting point for 

deriving the multi-step output prediction equation that is needed for designing a GPC controller. 

Using equation 1, the output at time step k+j may be written in the form 

 

() () ()

1 2

(1) ( ) ( ) (j) ( )

0 12 o

(1) ( ) ( ) (j) ( )

0 12 o

( ) ( -1) + ( -2) + + ( - ) 

 + ( ) + ( 1) ( ) + ( -1) + ( 2) + ( - )

+ ( ) + ( 1) ( ) + ( -1) + ( 2) + (

jj j

p

jj j

o p

jj j

o p

yk j yk yk ykp

uk j uk j uk uk uk ukp

dk j dk j dk dk dk d

αα α

β β ββ β β

γ γ γγ γ γ

+ =

+ +−+ + − +

+ +− + + − +

"

" "

" " k p- )

 (9)

where the coefficient matrices are given by recursive expressions involving the matrices αi, βi, 

and γi appearing in the ARX equation (ref. 10). The system identification process described 

earlier determines the matrices αi, βi, and γi. Equation 9 shows that the output y(k+j) at time step 

k+j may be estimated by using p sets of the previous output and input measurements, y(k-1),..., 

y(k-p) and u(k-1),...,u(k-p), and the (unknown) current and future inputs u(k), u(k+1),...,u(k+j) 

and disturbances d(k), d(k+1),...,d(k+j). The GPC algorithm is based on system output 

predictions over a finite horizon hp known as the prediction horizon. To predict future plant 

outputs, some assumption needs to be made about future control inputs. In determining the 

future control inputs for GPC, it is assumed that control is applied over a finite horizon hc known 

as the control horizon that is equal to or less than the prediction horizon. Beyond the control 

horizon the control input is assumed to be zero. Letting j in equation 9 range over the set of 

values j = 1, 2, ..., hp-1, the resulting equations can be assembled into a multi-step output 

prediction equation having the form 

 

1 1 1 1

1 1

() () ( ) ( )

( ) ()

p cc p c c p

p c c

p d d p pp d d

p c h m hr h m pr pr pm h m pm h m h r

p f f h m pr pr hm hr hr

p p h h

p

y k u k ukp ykp

dkp dk

× × × × ×

× ×

× × × ×

= + −+ −

+ −+

TB A

D D (10)

The coefficient matrices T, A, B, Dp and Df are formed from combinations of the observer 

Markov parameters αi, βi, and γi. The quantity yhp(k) is the vector containing the future outputs, 

whereas uhc(k) is the vector containing the future control inputs yet to be determined. The 
8 

quantities up(k-p), yp(k-p), and dp(k-p) are vectors containing the previous p sets of control inputs, 

outputs, and disturbance inputs, respectively. The vector df (k) contains the future disturbance 

inputs, which are usually unknown. However, if the disturbance is a coherent signal, future 

values of the disturbance can be estimated using an AR model of the disturbance obtained in a 

separate identification computation using a portion of the measured disturbance time history. 

The expanded form of this multi-step output prediction equation for hc + hp is shown in equation 

11. 

( ) ( 1) ( 1)

0 1 10

(1)

( 1) ( 2)

( ) ( 1) (1)

( 1) ( 2) ( )

( )

( 1) ( )

( 1) ( 1)

( ) ( 1)

( 1) p p pc

qq q

o o

q q

oo o

q q

oo o

c

h h hh

p oo o

y k o

y k u k

u k yk q

yk q uk h

yk h

β β αβ

β β

ββ β

ββ β

ββ β

β

− −

− −

−

−− −

= +

⎧ ⎫ ⎡ ⎤

⎪ ⎪ ⎢ ⎥ + ⎪ ⎪

⎪ ⎪

⎪ ⎪ + ⎨ ⎬ + − =

⎪ ⎪ + ⎪ ⎪ + − ⎪ ⎪

⎪ ⎪ + − ⎩ ⎭ ⎣ ⎦

# # #%

" # "

# # #%#

"

( ) ( 1) ( 1) ( ) ( 1) ( 1)

1 11 1 11

(1) (1) (1) (1)

12 1

( 1) ( 1) ( 1) ( 1)

12 1

() () () ()

12 1

( 1) ( 1) ( 1)

12 1

12 1

pp p

qq q qq q

pp p pp p

p p

qq qq

p p

qq qq

p p

h

p

p

h h

p

h h c p

α α αα β β αβ

αα αα

α α

αα α

α α

αα αα

α αα

α

α

−− −−

−− −−

−

−

−

− −−

−

−

−− −

−

⎧ ⎫

⎪ ⎪

⎨ ⎬

⎪ ⎪

⎪ ⎪ ⎩ ⎭

=+ =+

+

≤

"

"

# #%# #

"

"

# #%# #

"

(1) (1) (1) (1)

12 1

( 1) ( 1) ( 1) ( 1)

12 1

() () () ()

12 1

( 1) ( 1) ( 1) ( 1)

1 2

1

1

2 1

( 1)

( 2)

( 1)

( )

p pp p

p p

qq qq

p p

qq qq

p p

h

p

p p

hh h

p

y k

y k

yk p

yk p

β

ββ ββ

ββ ββ

ββ ββ

β β

β β

β

β

−

−− −−

−

−

− −− −

−

⎡ ⎤ −

⎢ ⎥ ⎧ ⎫ −

⎪ ⎪ − ⎪ ⎪

⎨ ⎬ +

⎪ ⎪ − + ⎪ ⎪

⎪ ⎪ ⎩ ⎭ −

⎣ ⎦

"

"

# #%# #

# "

"

# #%# #

" ( 1)

( ) ( 1) ( 1) ( ) ( 1) ( 1)

1 1 1 0 1 10

(1) (1) (1) (1)

12 1

( 1) ( 1) ( 1) ( 1)

1 2

1 2

1

(

1

1

( 1)

( 2)

( 1)

( )

p h

p

q q q qq q

pp p

p p

qq qq

p

q

p p

p

u k

u k

uk p

uk p

β

γ γ αγ γ γ αγ

γγ γ γ

γ

γγ γ

γ γ

γ

γ

γ

−

− − − −

− −

−

−− −

−

−

−

⎡ ⎤

⎢ ⎥ ⎧ ⎫ −

⎪ ⎪ − ⎪ ⎪

⎨ ⎬

⎪ ⎪ − + ⎪ ⎪

⎪ ⎪ ⎩ ⎭ −

⎣ ⎦

=+ =+

+

#

"

"

# #%# #

"

(1)

( 1) ( 2)

) () () () ( ) ( 1) (1)

2 1

(2) (1)

( 1) ( 1) ( 1) ( 1) ( 1)

12 1

( 1)

( 2)

( 1)

( )

pp pp p

o o

q q

oo o

q qq q q

p p o o oo

oo o

hh hh h

o

o

p p

d k

d k

dk p

dk p

γ γ

γγ γ

γ γγ γ γ γγ

γγ γ

γγ γγ γ γ

γ

− −

−

−

−− −− −

−

⎡ ⎤

⎢ ⎥ ⎧ ⎫ −

⎪ ⎪ − ⎪ ⎪

⎨ ⎬ +

⎪ ⎪ − + ⎪ ⎪

⎪ ⎪ ⎩ ⎭ −

⎣ ⎦

# #%

# "

" "

# #%# # # #"

" ( 2) (3) (2) (1)

( )

( 1)

( 1)

p

p

h

o o o oo

d k

d k

dk h

γγ γγ −

⎡ ⎤

⎢ ⎥

⎧ ⎫

⎪ ⎪ +

⎨ ⎬

⎪ ⎪

⎪ ⎪ + − ⎩ ⎭

⎣ ⎦

#

"

 (11) 
9 

The observer Markov parameters αi , βi, and γi determined in system identification are used 

to form the first block row in each of the coefficient matrices T, A, B, Dp and Df in equation 10. 

These terms are colored red in equation 11. The terms in the remaining block rows are computed 

using the recursive relations indicated in the boxes above the coefficient matrices in equation 11 

(ref. 10). All terms in equation 11 are known, except for the hc sets of future commands, the hp

sets of predicted responses, and the hp sets of future disturbance inputs df (k). The goal of the 

GPC control algorithm is to determine the set of future commands u(k), u(k+1), ..., u(k+hc-1) 

that are required to achieve a desired predicted response y(k), y(k+1), ..., y(k+hp-1). It should be 

noted that the system Markov parameters, which are commonly used as the basis for identifying 

discrete-time models for linear dynamical systems, form the first block column in the matrix T; 

the remaining block columns are formed from subsets of these parameters. The system Markov 

parameters are the pulse response of a system and are unique for a given linear system. The 

discrete-time state-space matrices A, B, C, and D are embedded in the system Markov 

parameters. 

Derivation of Control Law 

 

 The predictive control law is obtained by minimizing the deviation of the predicted 

controlled response (as computed from the multi-step output prediction equation) from a desired 

response over a prediction horizon hp. To this end, one first defines an error function that is the 

difference between the desired (target) response yT (k) and the predicted response yhp (k): 

 

() ()

() () ( ) ( ) ( ) ()

T

T p ff p p

p

c

h

h p

k k

k k k p k p dk p dk

y y

yu u y

ε = −

= − − −− −− −− TB A D D (12)

An objective function J quadratic in the error and the unknown future controls is then formed: 

 c c

T T

h h J Q uR = + ε ε u (13)

Two weighting matrices are included in the objective function: R (symmetric and positive 

definite) is used to weight the control effort and stabilize the closed-loop system; Q (symmetric 

and positive semi-definite) is used to weight the relative importance of the differences between 

the target and predicted responses. Typically, R and Q are assumed to be diagonal with R having 

the same value wc along its diagonal and Q having the same value wr along its diagonal. 

Minimizing J with respect to uhc(k) and then solving for uhc(k) gives 

 ( ) †

1

() ( ) () ( ) ( ) ( ) ()

C c

T T

T p pp f f c p

h r

hu k R Q R yk uk p yk p dk p dk

×

=− + − + − + − + − + TT T B A D D

 (14) 

as the control sequence to be applied to the system over the next hc time steps. However, only 

the first rc values (corresponding to the first future time step) 
10 

 1

() () ( ) ( ) ( ) () c c

Tp p p f

c

cc c

r

c u k yk uk p yk p dk p dk χβ α δ γ

×

=− + − + − + − + (15) 

are applied to the rc control inputs, the remainder is discarded, and a new control sequence is 

calculated at the next time step. Equation 15 is the desired control law equation. The target 

response yT (k) is zero for a regulator problem and non-zero for a tracking problem. The 

coefficient matrices χ

c

, αc, βc, δc

, and γ

c

 are the control law gain matrices and have sizes rc × mhp, 

rc × pm, rc × prc, rc × prd, and rc × rdhp, respectively. 

Selection of Control Model Parameters 

An expression for estimating the order of the ARX model that is to be used for SID is 

given by 

 number of system states number of disturbance states p ceil

m

+

≥ ⎛ ⎞ ⎜ ⎟ ⎝ ⎠ (16) 

where ceil denotes rounding up the value of the quantity in parentheses to the next higher 

integer. The number of system states is typically chosen to be twice the number of significant 

structural modes; the number of disturbance states is set to twice the (estimated) number of 

frequencies in the disturbance; m is the number of output measurements. If measurement noise 

is of concern, the value of p so computed should be increased to allow for computational poles 

and zeros to improve system identification in the presence of noise. In practice, simply choosing 

p to be about five times the number of significant modes in the system is often adequate. The 

prediction and control horizons are set according to the relations 

 h p hh p cp ≥ ≤ (17)

By extending the prediction and control horizons to very large values, the GPC solution 

approaches that of the linear quadratic regulator (LQR). Thus, GPC approximates an optimal 

controller for large hp (ref. 9). 

Weighting matrices R and Q are used to weight the control effort and to weight the relative 

importance of the differences between the target and predicted responses, respectively. In 

practice, the relative values of the weighting matrices play the major role in the optimization 

process. Therefore, the weighting matrix Q is often set to an identity matrix. The control weight 

matrix R then becomes the only one needed to be tuned to produce an acceptable solution 

without going unstable. Reducing wc increases controller performance but may saturate the 

actuators or drive the control system unstable. 

Practical Issues 

 Two issues need to be addressed before applying the GPC method to a real problem: the 

treatment of external disturbances and basic signal processing guidelines. 
11 

Treatment of External Disturbances 

If the external disturbances acting on a system are unknown, the best that can be done is to 

identify the system in the presence of these disturbances so as to include their effects in the 

identified OMP. The control law obtained using the OMP with the embedded disturbance 

information is then used to compute the control inputs. This is the procedure that was employed 

successfully in the GPC-based active controls testing of the tiltrotor model. 

If any of the external disturbances acting on the system are known or can be measured a 

feedforward loop can be employed in combination with the feedback loop. There are two 

feedforward terms (e.g., see equation 10): a term involving past values dp and a term involving 

future values df. If only past values are available one simply discards the other term. If future 

values are also available from an AR prediction model then both terms are retained. If there is 

correlation in the disturbance signals and the character of those disturbances is not changing 

significantly with time an AR model can be determined using an initial portion of the disturbance 

data. This model can then be used to predict the future values needed in equation 15. If the 

disturbances are changing significantly with time (e.g., due to frequency or amplitude variations, 

changes in the pattern or type of disturbance, etc) but there is coherence in the signals, an AR 

model can still be determined. However, for the AR model to adapt to the changing conditions it 

would have to be updated (ideally at each time step) using some type of recursive least squares 

algorithm. If there is no coherence in the disturbance signals no predictions of future 

disturbances are possible and only the term involving past values of the disturbance can be 

retained. These aspects are addressed in Appendix A. 

Basic Signal Processing Guidelines 

 Several signal-processing considerations must be taken into account when using either 

measured data or simulated data in GPC-based active controls work. The measured response 

time histories must be passed through a low-pass filter with a cut-off frequency fc chosen so that 

the maximum frequency of interest is about 75% of fc. Assuming that fc is equal to the Nyquist 

frequency fN, the required sampling frequency fS should be at least twice this value of fN to 

prevent aliasing. However, if fS is made too large the low-frequency modes will be poorly 

identified due to a loss of frequency resolution. A sampling rate between 2 to 3 times fN is 

generally sufficient. Once the sampling frequency has been selected, the minimum number of 

data points required for system identification is at least 5-10 cycles of the lowest frequency mode 

in the measured response time histories. Normalization of the input and output data that is used 

for identification on the maximum actual or expected values of the data is often helpful 

numerically. 

 It should be noted that active controllers are usually designed using a reduced-order 

mathematical model that includes only the most important subset of modes of a system. The 

frequency band spanned by the retained system modes sets the effective bandwidth of the 

resulting reduced-order controller. Such controllers may destabilize higher-order modes of the 

system due to truncation and spill over effects. The modes having frequencies in the transition 

region between retained and truncated modes are particularly susceptible. While these aspects of 

controller design are dealt with only casually below (e.g., through the choice of p, sampling 
12 

frequency, decimation factor, frequency range of interest, filter cutoff frequency, and number of 

data points used for SID), it is recognized that such considerations often require substantial 

attention by controls engineers. 

Results and Discussion 

 Results from two types of closed-loop simulations are presented. One uses a 

mathematical model of a three-degree-of-freedom mass-spring-dashpot system subject to userspecified disturbance force time histories. The other uses gust response data obtained from the 

wind-tunnel test of the HiLDA wing model. The intent of these simulations is to demonstrate the 

efficacy of a GPC-based active control system for gust load alleviation and the improvement in 

performance associated with including a feedforward path if any of the external disturbances are 

known, measurable, or predictable. The results presented below have been selected primarily 

with these objectives in mind. 

Comment on Codes Used for Simulations 

Two MATLAB programs were developed and used for the simulation studies that 

provided the numerical results for this report. One is an extension of the three-degree-offreedom (3-DOF) code that was used extensively to develop and evaluate the computational 

procedures that underlie the GPC software system that was used for the active controls testing of 

the tiltrotor model. This program treats a simple system with known dynamic properties and 

user-defined external disturbances. It also lends itself to easy modification should the need arise 

to address either a what-if question or a computational issue. The other program is new and was 

written specifically to access and process the data sets that were obtained during the wind-tunnel 

test of the HiLDA wing model (ref. 6). For that reason, it is tailored to treat that system. 

Although the programs are considerably different if viewed from a coding perspective, they are 

based on the same theoretical and computational principles and are therefore equivalent. 

Numerical Simulations Using 3-DOF Math Model 

 The 3-DOF mass-spring-dashpot system used in the 3-DOF simulations is depicted in 

figure 2 and the corresponding equations of motion are summarized in Appendix B. The model 

has three masses, three springs, and three dashpots. Disturbance forces ud1, ud2, and ud3 and 

control forces uc1, uc2, and uc3 act on masses m1, m2, and m3, respectively. The degrees of 

freedom x1, x2, and x3 are the absolute displacements of the corresponding masses. For the 

simulations, the masses were set equal to 1.0 lb-sec2

/in, the spring rates to 1000 lb/in, and the 

damping coefficients to 0.01 lb/in/sec. This yields natural frequencies of 2.240, 6.276, and 9.069 

Hz and corresponding damping ratios of 7.0367e-005, 1.9716e-004, and 2.8491e-004 for the 

three natural modes of the system. Disturbances, control forces, and responses (displacements, 

velocities, or accelerations) can be specified and selected in any combination (see Appendix B). 

However, for the results to be shown here, controls were imposed at masses 1 and/or 3 (except 

for one case of three control inputs), disturbances were imposed at masses 1 and/or 2, and 

accelerations were measured at masses 1 and/or 2. Both steady-state (constant amplitude and 

frequency) and non-steady-state (varying amplitude and frequency) disturbances were imposed 

on the system. 
13 

 The steady-state disturbances (figure 3) were chosen to be sinusoidal with a frequency of 

6.276 Hz, which coincides with the second natural frequency of the system. The expressions 

used to represent these disturbances in the simulations were ud1 = 2*cos(2*pi*f*t) and ud2 = 

sin(2*pi*f*t). The non-steady-state disturbances had variable and/or discontinuous amplitude 

and frequency time histories. The first of these was composed of four variable-amplitude, 

constant-frequency segments (fig. 4a); the other represented a disturbance with increasing 

frequency and decreasing amplitude (fig. 4b). 

 The MATLAB function c2d was used to convert the continuous-time state-space models 

of the 3-DOF system (Appendix B) to discrete-time using t =0.05 as set below. The resulting 

discrete-time state-space models, along with their corresponding discretized disturbances, were 

then used in the MATLAB function dlsim to generate the open-loop time histories needed for 

SID. 

 The maximum frequency of interest for the 3-DOF system is just over 9 Hz. Based on 

the signal processing guidelines mentioned earlier, the sampling frequency fs was set to 20 Hz 

(t=0.05 sec) and the number of time steps A of input/output data was set to 300. The control 

horizon hc was set equal to the prediction horizon hp and hp was set equal to the order of the 

ARX model p. The value of p was set to 6, 8, or 10 depending on the control case being 

considered and the number of disturbances acting on the system. For the cases of feedback-only 

and feedback with feedforward p was set to 6. The use of feedback with embedded feedforward 

requires p to be increased by (at least) two for each of the dominant frequencies present in the 

external disturbances acting on the system. This expedient provides the additional degrees of 

freedom needed to ensure that the effects of the unknown periodic disturbances are embedded in 

the observer Markov parameters that are computed during system identification. Thus, p was set 

to 8 when using the single disturbances in figure 4 and to 10 when using the double disturbances 

in figure 3. Based on previous experience with the 3-DOF model, the response weights wr were 

set to 1.0 and the control weights wc were taken to be equal and set to 0.3. In practice, the 

response weights wr are usually set to unity and wc is varied to tune the controller to the nuances 

of the subject system. Tuning proceeds by initially setting the value of wc to a “conservative” 

value (e.g., 1, as was initially done for the 3-DOF system) and then tracking the value of the 

objective function J given in equation 13 as wc is incrementally reduced until either a minimum 

value of J is identified or the control system actuators become saturated. This process for 

selecting wc is analogous to the familiar procedure for selecting the controller gain in a singleinput/single-output system. Previous experience with the 3-DOF model indicates that the 

selected values are all reasonable but not optimum. No attempt was made to “optimize” (i.e., 

tune) the parameters p, hc, hp, wr, and wc to the nuances of the system under the different 

disturbances to elicit the best (lowest) possible closed-loop responses. In the simulations, the 

loop was closed after 100 time steps (5 sec) when using steady-state disturbances and 20 time 

steps (1 sec) when using non-steady-state disturbances. 

 Results are presented below for three types of control cases: (1) feedback only, in which 

no account is taken of any external disturbances acting on the system during SID and controller 

design; (2) feedback with embedded feedforward, in which all the disturbances acting on the 
14 

system are assumed to be unknown while performing the SID and controller design; and (3) 

feedback plus feedforward, in which all of the external disturbances acting on the system are 

assumed to be measurable and are used during the SID and design of the feedforward controller. 

A second-order AR disturbance model (nd = 2) (see Appendix A) was used for predicting the 

future disturbances needed by the feedforward controller. This portion of the feedforward 

controller was adaptive and used RLS with a forgetting factor (ref. 21) of 0.9. Results are 

presented for both steady-state (constant amplitude and frequency) and non-steady-state (varying 

amplitude and frequency) disturbances. For ease of discussion, the second control case will be 

referred to simply as ‘embedded feedforward’ and the third control case will be referred to as 

‘explicit feedforward’. 

 While some of the results shown below indicate the effect on performance of increasing 

the number of control effectors, there is no discussion of trade-offs between the number of 

control inputs used and the costs associated with implementation or power requirements. These 

would be important considerations in any studies aimed at full-scale implementation but are 

beyond the scope of this report. 

Representative Results 

Steady-State Disturbances: Representative results obtained from simulations made using the 

steady-state disturbances shown in figure 3 acting on masses 1 and 2 (ud1 and ud2 in fig. 2) are 

given in figures 5-9. Loop closure (active control on) was taken to occur after five seconds. 

Figures 5-7 show results for feedback only, embedded feedforward, and explicit feedforward, 

respectively, for the case in which control is applied at mass 3 (uc3) and feedback is the 

acceleration of mass 1 ( 1 x ). In all three cases, the response was quickly reduced to very small 

values when the loop was closed and then maintained at those low values thereafter. The closedloop response levels for embedded feedforward (fig. 6) are essentially the same as those for 

explicit feedforward (fig. 7), but both are smaller than those for the feedback-only case of figure 

5. However, the feedback-only case also produces a considerably reduced response level. This 

result indicates that ignoring the disturbances during SID may produce an acceptable ARX 

model for designing a GPC controller for a system. The control inputs for all three cases are 

nearly the same. 

 The effect of applying control forces at two masses (uc1 and uc3) rather than just at one 

mass (uc3) as in figure 7 is illustrated in figure 8 for explicit feedforward. Comparing these 

results to those of figure 7 indicates that the use of an additional control force reduced both the 

closed-loop response and the magnitudes of the required control forces below the values 

associated with a single control, as might be expected. 

 Figure 9 shows results for explicit feedforward for the case of control inputs at masses 1 

and 3 (uc1 and uc3) as in figure 8 but with acceleration outputs at both masses 1 and 2 ( 1 x and 2 x ). 

Both responses are very quickly quashed after the loop is closed. There is a slight increase in the 

magnitudes of the required control forces because two responses are now being worked by the 

controller rather than just one as in figure 8.
15 

Non-Steady-State Disturbances: Typical results obtained from simulations made using the 

piecewise varying frequency and amplitude disturbance shown in figure 4a acting on mass 2 

(ud2) and taking the feedback to be the acceleration of mass 1 ( 1 x ) are given in figures 10-14 for 

cases involving one to three control inputs. Loop closure (active control on) was imposed after 

one second. All responses were reduced to varying extents depending on the number of control 

inputs. Figures 10 and 11 show results using embedded feedforward and explicit feedforward, 

respectively, for the case in which control is applied at mass 3 (uc3). The closed-loop response 

using explicit feedforward is slightly better than for embedded feedforward while the magnitudes 

of the control inputs are either slightly larger, smaller or about the same depending on the 

segment of the disturbance that is active. 

 Figures 12 and 13 show results for embedded feedforward and explicit feedforward, 

respectively, for the case of control inputs at masses 1 and 3 (uc1 and uc3). The closed-loop 

response is only slightly better using explicit feedforward (fig. 13) than embedded feedforward 

(fig. 12) while the required control inputs are about the same. Increasing the number of control 

inputs from one to two produces lower responses (compare figure 12 to figure 10 and figure 13 

to figure 11). 

 

 Figure 14 shows results using explicit feedforward for the case in which control inputs 

are at masses 1, 2 and 3 (uc1, uc2 and uc3). The closed-loop responses and controls in this case 

are improved relative to the two-control-input results shown in figure 13, as should be expected.

 Results from simulations made using the disturbance shown in figure 4b acting at mass 2 

(ud2) and taking the feedback to be the acceleration of mass 1 ( 1 x ) are given in figures 15-17. 

Loop closure was taken to occur after one second. Figures 15 and 16 show results for embedded 

feedforward and explicit feedforward, respectively, for the case of a control input at mass 3 (uc3). 

The response is reduced in both cases but is decidedly better using explicit feedforward with 

adaptive identification of the disturbance. Control inputs for the latter are slightly larger during 

the early (lower frequency) portion of the disturbance but lower during the latter (higher 

frequency) portion of the disturbance. Figure 17 shows results using explicit feedforward for the 

case of control inputs at masses 1 and 3 (uc1 and uc3). The closed-loop response in this case is 

lower than that for the case of one control in figure 16, and control inputs are slightly smaller 

during the early (lower frequency) portion of the disturbance but generally the same or larger 

during the latter (higher frequency) portion. 

Summary of Observations: Taken as a whole, the 3-DOF simulations indicate that GPC-based 

controllers are able to significantly reduce acceleration responses under a variety of imposed 

disturbance conditions and to do so very quickly after closing the loop. For the cases involving 

steady-state disturbances (figs. 5-9), the lightly-damped system was being excited at one of its 

natural frequencies and the open-loop responses are all diverging with time. Thus, the 

percentage reduction in response after closing the loop in these cases depends on when the loop 

is closed as the magnitude of the steady-state closed-loop response is essentially independent of 

loop closure time. For figures 5-8 (2 disturbances, 1 output, 1 control) the amplitude of the 

open-loop response at time of closure (5 sec) is about 70 in/sec2

. With respect to this value, the 

reduction in acceleration responses ranged from about 90% to essentially 100%. For the case of 

figure 9 (2 disturbances, 2 outputs, 2 controls) both 1x and 2 x were effectively reduced to zero. 
16 

The results in figures 10-17 are for disturbances with varying frequency and amplitude and can’t 

be quantified with a single number because the percentage reduction in response varies with the 

frequency and amplitude of the disturbance. However, depending on the position in the time 

history, the reduction varies from about 50% to essentially 100%. Given that the maximum 

values of the disturbance forces vary from one to two pounds in magnitude the required steadystate control forces in figures 5-17, which range from about 0.2 lbs to 2 lbs, are considered 

reasonable. 

 With respect to the effects on performance associated with varying the number of control 

inputs and response outputs, the 3-DOF simulation results indicate that increasing the number of 

controls rc for a fixed number of outputs m improves performance (both responses and control 

inputs are lower). Conversely, increasing the number of outputs m for fixed rc increases the 

required control inputs slightly but still quells all the responses. If both m and rc are increased 

simultaneously, all the responses are reduced with only a slight increase in the required control 

inputs. As mentioned earlier, no attempt was made to tune any of the GPC parameters to the 

nuances of the system under the different disturbances to obtain the best (lowest) possible 

closed-loop responses. However, the results (not shown) of some limited parametric variations 

indicate that increasing the order p of the ARX model and/or reducing the magnitude of the 

control weight wc can improve the performance of the controller in all of the 3-DOF simulations 

presented above. 

Numerical Simulations Using Data from Wind-Tunnel Test of HiLDA Wing Model 

 The HiLDA wing model, as mounted in the Langley Transonic Dynamics Tunnel (TDT), 

is shown in figure 18. The model has four trailing-edge control surfaces that span the lifting 

portion of the wing and one leading-edge control surface located near the wing tip. A simulated 

vertical gust field was applied to the test section flow in the form of a sinusoidal oscillation of 

the flow direction. This oscillating flow was generated by an airstream oscillator system 

consisting of a biplane arrangement of vanes located on the side walls of the tunnel about 70 feet 

upstream from the model (fig. 19). The gust-induced flow angle was measured by a vane-type 

flow direction transmitter that was mounted on the top of a seven-foot high post attached to the 

tunnel floor six feet ahead of the model (fig. 18). A schematic of the wing planform is depicted 

in figure 20, which shows the arrangement of the five control surfaces that were used as control 

effectors and the locations of the three strain gage sensors that provided the feedback response 

time histories for the simulations. The reader is referred to reference 6 for a description of the 

model. 

Test Conditions 

 

 The data used in the simulations presented below were obtained during a run in the TDT 

heavy gas (R134a) test medium at a dynamic pressure of 51 lb/ft2

. The corresponding Mach 

number and velocity were 0.25 and 138 ft/sec, respectively. Sinusoidal vertical gust excitations 

were induced in the test section by the vertical components of the tip vortices shed from the 

oscillating biplane vanes. The amplitude of oscillation of the airstream oscillator was set to ±6 

degrees. This amplitude was sufficient to produce gust angles as large as one degree at the 

location of the model at the lower airspeed/frequency combinations. The airstream oscillator 
17 

excitation signatures consisted of constant frequency dwells at the natural frequencies of the first 

and second vertical bending modes (f1v ≈ 2 Hz, f2v ≈ 8 Hz), and frequency sweeps from 0.5 Hz 

to 10 Hz in 30 seconds. For SID purposes random excitations of the five control surfaces were 

imposed simultaneously and consisted of uncorrelated, normally-distributed random numbers 

with mean zero and a variance chosen to limit the oscillation angles of the control surfaces to 

about ± 2 degrees. 

Test Procedures 

 The SID data for the case in which the gust excitation frequency is held constant were 

obtained as follows: (1) Set tunnel dynamic pressure to the desired value; (2) Set the frequency 

of the airstream oscillator to its desired value; (3) Impose a random excitation on each of the five 

control surfaces; and (4) Record the resulting input and output time histories for 35 seconds. For 

the case in which the gust frequency is swept from a low value to a high value the procedure 

was: (1) Set tunnel dynamic pressure to the desired value; (2) Set the frequency of the airstream 

oscillator to its starting value of 0.5 Hz; (3) Impose a random excitation on each of the five 

control surfaces; and (4) Sweep the gust frequency continuously from 0.5 Hz to 10 Hz in 30 

seconds while recording the resulting input and output time histories. 

Simulation Conditions and Parameter Settings 

 The original data were taken using a sampling rate of 500 Hz and recorded for 35 

seconds, yielding data records consisting of 17500 time samples. To exclude transients, the data 

outside time samples 1000 and 15000 were discarded. Since the wing modes to be controlled 

were less than 15 Hz, the truncated data were then decimated by a factor of 10 using the 

MATLAB command decimate, resulting in data representing a sampling rate of 50 Hz. Finally, 

the mean value of each measurement was subtracted from its time history. 

 The conditioned input-output data were used to identify the plant and to develop the 

controller. SID was accomplished using the command arx_ud from the predictive control 

toolbox. The outputs from strain-gage sensors SBI2, SBY, and STI2 shown in figure 20 were 

selected as the quantities to be minimized collectively (m = 3). All five control inputs, CTL1, 

CTL2, CTL3, CTL4, and CTL5, were utilized simultaneously as control effectors (rc = 5). The 

feedforward signal was taken from the flow direction transmitter. 

 System identification returned the ARX parameters of equation 1. The control law gain 

matrices in equation 15 were then obtained by processing the ARX parameters according to the 

computational methods described earlier in this report. The ARX parameters were used for the 

time history simulation of both the open- and closed-loop plant, as well as for the design of both 

the embedded feedforward and explicit feedforward controllers. Controller design was 

accomplished using a modified version of the m-file msop from the predictive control toolbox. 

The same controller was used for embedded feedforward and explicit feedforward simulations, 

the difference being that in the former case the feedforward controller path was not utilized. For 

closed-loop simulations employing embedded feedforward, only the terms involving the past 

control inputs up and feedback outputs yp in equation 15 were used to compute the control inputs 

uc. For cases employing explicit feedforward, all four terms were used. In order to use the df
18 

term of equation 15, some estimate must be made of future disturbance inputs given the history 

of past disturbance measurements. As discussed earlier, depending on the character of the 

disturbances acting on the system, either an adaptive or a non-adaptive AR model can be used to 

predict future disturbances. Here, an adaptive 6th order AR model was chosen to predict the 

future disturbances. A recursive least squares procedure with a forgetting factor of 0.97 was 

used to update the AR model parameters every time step, thus allowing for enhanced tracking of 

the swept sine gust. 

 A comment is in order here regarding some additional data processing that was necessary 

to allow a proper comparison of measured open-loop time histories with simulated closed-loop 

time histories. The desired open-loop time histories (quiescent controls during each of the three 

gust-excitation conditions) were not among the data sets collected during the test. The closest set 

of data was that taken for SID wherein the five control surfaces were randomly excited during 

gust excitation. Fortunately, it is possible to compute the desired open-loop responses using that 

data set. Essentially, one first performs an SID using as input the random excitations of the 

control surfaces and the gust vane measurements. This yields the OMP αi, βi and γi of equation 

1. One then performs a time-history simulation for y(k) using equation 1 but retaining only the 

αi and γi terms on the right-hand-side and using for d the same set of gust vane measurements 

used for the SID. The resulting time histories are the desired open-loop time histories. To 

maintain consistency, the same set of gust vane measurements is added to the plant response y at 

each time step during the closed-loop simulations. 

 The system model identified as described above was used in both the open-loop and 

closed-loop simulations for the cases of both embedded feedforward and explicit feedforward 

controllers. The control inputs uc(k) were set to zero for the open-loop simulation. For the 

closed-loop simulation the control signal computed using equation 15 was used. The control 

input for embedded feedforward used only the past outputs yp and past control inputs up; no use 

was made of the past and future feedforward terms (dp and df). For the explicit feedforward 

controller the computed control signal uc(k) was obtained using all four terms on the right-handside of equation 15: the past outputs, past control inputs, past measured feedforward values, and 

future predicted feedforward values. The same weights that were used for the 3-DOF 

simulations (wr= 1 and wc= 0.3) were tried for the HiLDA wing data and found to be acceptable. 

The value of p was estimated using equation 16 and set to 10. Its value was not varied with 

controller type as for the 3-DOF model. The controller was turned on five seconds into the 

simulations. 

Representative Results 

Figures 21 through 62 show simulated open- and closed-loop responses of wing root 

vertical bending (SBI2), fore-and-aft bending (SBY), and torsion (STI2) to sinusoidal vertical 

gust excitation held constant at the natural frequencies of the first and second vertical bending 

modes and while sweeping the frequency of the gust from 0.5 Hz to 10 Hz. Bending and torsion 

moment response time histories, power spectral densities of bending and torsion moment time 

histories, and control angle input time histories are shown for each of the three gust excitation 

conditions for two control cases: (1) Feedback with embedded feedforward; and (2) Feedback 

with explicit feedforward and using an adaptive AR model for predicting the future disturbances. 
19 

Feedback with Embedded Feedforward: Figures 21-27 indicate controller performance in 

reducing wing root moments due to gust excitation at the natural frequency of the first vertical 

bending mode (2 Hz). The SBI2 closed-loop time history (fig. 21) shows a dramatic reduction of 

that modal response almost immediately after the controller is turned on. The corresponding 

SBY and STI2 closed-loop time histories (figs. 22-23) are also significantly, and promptly, 

reduced. The PSD plots (figs. 24-26) indicate that the primary effect has been to reduce the 

response at the frequency of excitation (the first mode), as would be expected. However, the 

closed-loop spectra are not flat indicating that there is residual correlation in the sensor signals. 

It may not be possible to remove all correlation in regulated systems due to stability issues, the 

presence of unobservable and/or uncontrollable modes, or limits on actuator authority. Based on 

experience with GPC applied to other systems, it is known that tuning of the GPC parameters 

(horizons, control penalty, order) oftentimes leads to enhanced performance. However, no 

adjustments were made in these parameters to optimize performance in the results shown here or 

below. Figure 27 shows the commanded control deflection angles to be small once transients 

decay. In a physical control system, the overshoot occurring at loop closure would be 

suppressed by low-pass filtering the control signal and/or placing upper and lower bounds on the 

allowable control amplitudes in both the hardware and software. The roll-off occurring in the 

magnitudes of the PSD curves at frequencies above about 20 Hz is due to the filtering action of 

the function decimate that was used to decimate the original wind-tunnel data. 

 Figures 28-34 summarize controller performance due to gust excitation of the wing at its 

second vertical bending natural frequency (8 Hz). The response time histories (figs. 28-30) show 

a dramatic and essentially immediate drop in amplitude once the control is applied. The PSD 

plots (figs. 31-33) show a relatively flat controlled spectrum. Again, the primary effect has been 

to reduce the response at the frequency of excitation (the second mode), as would be expected. 

In this case the response in that mode is virtually eliminated. Additionally, the SBY responses at 

5 Hz (the first fore-and-aft bending mode) and at 16 Hz are greatly reduced. Figure 34 shows the 

control effort to be small once transients decay. In practice, the overshoot occurring at 5 seconds 

would be suppressed by low pass filtering the control signal and/or placing upper and lower 

bounds on the allowable control amplitudes. 

 Figures 35-41 show controller performance for the case in which the gust frequency is 

swept from 0.5 Hz to 10 Hz. In comparison to the constant-frequency gust excitation cases of 

figures 21-34, controller performance is only fair. This is due to the fact that SID in the presence 

of a varying-frequency gust leads to an internal (i.e., embedded) disturbance model that is 

effectively an average model over the bandwidth of the gust. This leads to a controller that can 

not effectively reduce the responses at all the frequencies contained in the gust as in the constantfrequency gust cases. One possible solution to this problem is to use a GPC controller in which 

all the coefficient matrices χ

c

, αc, βc, δc

, and γ

c

 in equation 15 periodically adapt (see, for 

example, ref. 15). Figure 41 shows that the control deflection angles in this case are large 

compared to those required for a constant-frequency gust excitation. It is thought that this is 

because the responses to the disturbance input (gust vane measurement added to y) were larger 

than the responses to the random disturbances input by the control surfaces during SID. There 

was thus a greater correlation in the gust vane disturbance input-to-output path than in the 

control input-to-output path. A more accurate mapping from control input u to system output y
20 

would produce a better model of the plant. This could be achieved by either increasing the 

amplitude of the random noise applied to the control inputs for SID or reducing the level of the 

disturbance (gust vane measurement added to y) applied to the system during SID. 

Feedback with Explicit Feedforward: Figures 42-48 show the performance of the explicit 

feedforward controller for the case in which the gust frequency is set to that of the first vertical 

bending mode of the model. Recursive least squares is employed to update the coefficients in 

the AR model used to estimate the future values of the gust disturbance df given the past 

disturbance values dp, as described in Appendix A. In comparison to the embedded feedforward 

controller (figs. 21-27), the explicit feedforward controller exhibits better performance as can be 

seen by comparing time histories, PSDs, and required control angles. The steady-state control 

angles are all small after the transients have decayed. The improved performance is due to the 

incorporation of the measured past disturbance values dp into the controller as well as the ability 

of the future component of the feedforward path to tune itself to the changing conditions at each 

time step, thereby providing improved estimates of future disturbance values df. 

 Figures 49-55 are for the case of explicit feedforward with a gust frequency set to that of 

the second vertical bending mode of the model. The adaptive part of the controller is the matrix 

Dfp in equation A4 that is used to predict the future disturbance values df (k) needed in the last 

term in equation 15. Enhanced performance over the embedded feedforward case is observed as 

it was for the first vertical bending mode. The control angles are again all small once the 

transients decay. It is of interest to note that the closed-loop spectra (figs. 53-55) are effectively 

flat, indicating that the controller is approximately optimal. 

 Figures 56-62 show performance of the explicit feedforward controller for the case in 

which the gust frequency is swept from 0.5 Hz to 10 Hz. Improvement over the performance of 

the embedded feedforward controller (figs. 35-41) is due to the incorporation of a feedforward 

path (the last two terms in eq. 15) and the ability of the explicit feedforward component (the last 

term in eq. 15) to track changes in the gust frequency. The closed-loop spectra (figs. 60-62) are 

reasonably flat, indicating that the controller is near optimal. The required control angles are 

quite large compared to that required for quelling the responses for constant-frequency gust 

excitations. As discussed earlier in the case of the embedded feedforward controller, this may be 

due to the system identification being compromised by the magnitude of the disturbance signal 

relative to the noise input on the control channels during identification. Again, it is suspected 

that if the controller were fully adaptive (αc

, βc, δc

, and γ

c

 in equation 15 all adapt at each time 

step) performance would be much better and require much smaller control angles. 

Summary of Observations: Taken as a whole, the simulated closed-loop behavior of the HiLDA 

wing model during constant- and swept-frequency vertical gust excitation indicates that GPCbased controllers are able to significantly reduce wing bending and torsion moments and to do so 

rapidly after closing the loop. Tables 1 and 2 summarize the overall factors by which the SBI2, 

SBY, and STI2 wing root moments were reduced for each of the three gusts for the cases of 

embedded feedforward and explicit feedforward, respectively. For the case of embedded 

feedforward (table 1), the reductions range from a factor of four to as high as 667, depending on 

whether the excitation is from a constant- or swept-frequency gust and on the direction of the 

response. For the case of explicit feedforward (table 2), the reductions range from a factor of 
21 

four to as high as 1000, again depending on whether the excitation is from a constant- or sweptfrequency gust and on the direction of the response. The control angles required to effect these 

reductions ranged from one to three degrees for fixed-frequency gust excitations and from 10 to 

20 degrees for swept-frequency gusts. Viewed as a whole, these results are decidedly positive in 

regard to establishing the credibility of GPC as a serious candidate approach for active GLA. 

Additional wind-tunnel testing of the HiLDA wing model is planned to provide additional 

validation of the method. 

Concluding Remarks 

The results of numerical simulations aimed at assessing and evaluating the efficacy of a 

predictive control method known as Generalized Predictive Control (GPC) for active gust load 

alleviation (GLA) using trailing- and leading-edge control surfaces have been summarized. GPC 

is a linear, time-invariant, multi-input/multi-output predictive control method that uses an 

AutoRegressive with Exogenous input (ARX) model to describe the input-output relationship of 

the system. The coefficient matrices of the ARX equation are determined using system 

identification techniques. The input-output equation is used to form a multi-step output 

prediction equation over a finite prediction horizon. The control to be imposed at the next time 

step is determined by minimizing the deviation of the predicted controlled plant outputs from the 

desired (or target) outputs, subject to a penalty on control effort. The method lends itself to an 

implementation that allows all computations to be done on-line and in real time. Following a 

qualitative discussion of the essential features of GPC, the equations underlying the method were 

presented and discussed, including system identification, calculation of control law matrices, and 

calculation of the commands applied to the control effectors. Both embedded and explicit 

feedforward paths for incorporation of disturbance effects were addressed. 

Representative results from two types of simulations were presented. The first used a 

three-degree-of-freedom mathematical model of a mass-spring-dashpot system subject to a 

variety of user-defined external disturbances. The second used open-loop data obtained in a 

wind-tunnel while subjecting a wing model to constant-frequency and swept-frequency 

sinusoidal vertical gust excitations. Closed-loop behavior of the model was simulated in posttest calculations. Results obtained from these simulations have been decidedly positive. In 

particular, the results of all the closed-loop simulations conducted for the wing model showed 

reductions in wing root moments ranging from factors of four to as much as 1000, depending on 

whether the excitation is from a constant- or swept-frequency gust and on the direction of the 

response. The control angles required to achieve these reductions ranged from one to three 

degrees for fixed-frequency gust excitations and from 10 to 20 degrees for swept-frequency 

gusts. The GPC-based controller was evaluated for both embedded feedforward and explicit 

feedforward modes of operation. The effectiveness of the explicit feedforward controller using 

an adaptive AutoRegressive (AR) model for predicting future disturbances was particularly 

notable, both in improving performance and in tracking the variable disturbance frequency. 

Based on the results of these studies, it appears that GPC is a strong candidate approach for 

active GLA and warrants further investigation in the subject area. Additional wind-tunnel testing 

of the HiLDA wing model is planned to provide additional validation of the method. 
22 

Appendix A

Treatment of Measurable Coherent Disturbances 

 If the disturbance is measurable and coherent, an AutoRegressive (AR) model of the 

disturbance giving the disturbance at a future time step k can be written in the form 

 1 2 ( ) ( 1) ( 2) ( ) dn d dk dk dk dk n = −+ − + + − η η η " (A1) 

where the coefficient matrices ηi are rd × rd , the disturbance vectors d are rd × 1, rd is the number 

of disturbances, and nd is the order of the AR model. The coefficient matrices ηi in the model 

can be determined using established system identification techniques. The disturbance d(k+j) at 

an arbitrary future time step k+j can be written as 

 1 2 ( ) ( 1) ( 2) ( ) dn d dk j dk j dk j dk j n + = +−+ +− + + +− η η η " (A2) 

Letting j range over the set of values 0 to hp-1 one obtains the set of equations 

 

1 2

1 2

1 2

1 23

1 2

( 1

( 1) ( 2) ( )

( 1) ( 1)

( 2)

( 3) ( 3)

( 1) ( 2)

( 2)

( 2

( )

( )

( )

(

)

( 1

( 3

(

(

)

)

) 1) )

d

d

d

d

d

n d

n d

n d

n d

p p pn

dk dk dk n

dk dk n

dk n

dk dk n

dk h dk h dk

d k

d k d k

d k

d k

d k

d k

h dk

d k d k

h

η η η

ηη η

ηη η

ηηη η

ηη η

= −+ − + + −

= + −+ + − +

= + ++ − +

+= + + ++ −

+

+ +

+ −= + − + + −+

+

+

+ +

+

"

"

"

"

"

#

1) p d − − n

 (A3) 

The first equation in this set is the AR equation given in (A1) and succeeding equations are given 

using previously computed future disturbances (indicated in color) and measured past 

disturbances according to the pattern shown. The shifting pattern of the computations indicated 

in these equations is easily implemented in MATLAB. The results of these shifting calculations 

can be assembled into a matrix Dfp relating future disturbances to past disturbances according to 

the equation 

 

{ } [ ]{ }

1 1

() ( )

pd pd d d

f fp p

h r h r pr pr

dk D dk p

×× ×

= − (A4) 

and the term Df df (k) in equation 10 can be replaced by Df Dfp dp(k-p). 
23 

Batch Identification of Disturbance Models 

 Comparisons of several imposed/measured disturbance time histories with those 

predicted using a non-adaptive AR model are summarized in figures A1-A5. The AR model was 

identified using a small block of “measured” disturbance data. These results are intended to 

provide an indication of the ability of AR models to adequately represent correlated disturbances 

for GPC-based active control systems. 

 The imposed disturbance time history used for figures A1 and A2 is composed of four 

equal-length time segments, each containing 100 data points, for a total time period of 20 

seconds (4*100*0.05 = 20 sec, where t = 0.05 sec). Each segment has a different amplitude 

and frequency. The predictions in figure A1 were based on a model that was identified using 

only the first 50 data points (first 50*0.05 = 2.5 sec of data) while the predictions in figure A2 

were based on a model that used all 400 data points. A 6th order (six-term) AR model (nd = 6) 

was assumed for both. These results suggest that a disturbance that has considerable variation in 

its amplitude and/or frequency can be represented rather well by a low-order model based on a 

single identification using a modest-length block of data if there is correlation in the disturbance. 

This appears to be so even if there are abrupt changes in the character of the disturbance (and 

hence in the coefficient matrices ηi in equation A1). 

 Figures A3-A5 show comparisons of gust angles measured during the HiLDA test at a 

dynamic pressure of 51 lb/ft2

 with predictions obtained using a six-term AR model of the 

disturbance. The data were obtained at a sampling rate of 500 Hz but were decimated by a factor 

of ten to reduce the effective sampling rate to 50 Hz. Hence, the effective sampling time is t = 

0.02 sec. The decimated time history contained 1300 points. The first 500 data points (first 

500*0.02 = 10 sec) of this decimated time history were used to identify the coefficient matrices 

ηi of the AR disturbance model. Figure A3 shows the measured versus predicted time histories 

for the case in which the gust frequency was held constant near the first vertical bending mode 

(about 2 Hz). Similar results for a constant-frequency excitation near the second vertical 

bending mode (about 8 Hz) are given in figure A4. Figure A5 is for the case in which the gust 

frequency was swept from 0 to 10 Hz in 30 seconds. In spite of the fact that the gusts are not 

pure sinusoids and are contaminated by higher harmonics and noise, there appears to be 

sufficient coherence in the signals to allow use of non-adaptive methods to identify disturbance 

models that produce reasonable results, particularly at the lower frequencies. 

Adaptive Identification of Disturbance Models 

 Sometimes it is necessary to update the model of a system or a disturbance on-line and in 

real time while a system is in operation. Methods using models that are adjusted on-line using 

measured input/output data up to the current time are called adaptive. To be adaptive in real time 

requires that the measurements obtained at any time point can be processed during one sampling 

period. In such instances, recourse must be made to recursive identification techniques, which 

process data sequentially as they become available. A variety of recursive estimation methods 

are described in the literature (see, for example, refs. 11 and 21). Among these, the classical 
24 

recursive least squares (RLS) method is probably the most popular and simplest. For models 

whose assumed order is relatively small (such as the disturbance models of interest here), it is the 

method of choice. However, there are several fast versions of the method available that would 

be better suited to large models. 

 

 Adaptive identification of a disturbance model would update the rd × rd coefficient 

matrices ηi in equation A1 at each time step. This updated model would then be used at each 

time step to compute the vector df (k) in equation A4 and then df (k) used in equations 10, 11, 12, 

14, and 15 as appropriate. 

 

 Although the results of figures A1-A5 indicate that a non-adaptive identification scheme 

may be adequate for a wide range of disturbance time histories, there could be times when an 

adaptive scheme may prove advantageous or necessary. Figures A6-A9 show comparisons of 

several imposed and measured disturbance time histories with those predicted using an AR 

model of the disturbance that is updated at each time step using the classical recursive least 

squares technique. Figures A6 and A7 show results using 2-term and 6-term AR models, 

respectively. The imposed disturbance time history used is the same one that was employed in 

figures A1 and A2. Increasing the order of the AR model gives a predicted disturbance that is in 

better agreement with the imposed disturbance, except when transitioning from one type of 

disturbance to another where the differences increase away from transition points. 

 Figures A8 and A9 show results for the case of random disturbances using 6-term and 

16-term AR models, respectively. No agreement is evident (or should be expected) because 

there is no correlation in random signals. 
25 

0 2 4 6 8 10 12 14 16 18 20 −1.5

−1

−0.5

0

0.5

1

1.5

Time, sec

Force, lb

Imposed

Predicted

Figure A1.- Imposed versus predicted disturbance time histories using a 6-term AR model 

identified using the first 50 data points (first 2.5 sec of data). 
26 

0 2 4 6 8 10 12 14 16 18 20 −1.5

−1

−0.5

0

0.5

1

1.5

Time, sec

Force, lb

Imposed

Predicted

Figure A2.- Imposed versus predicted disturbance time histories using a 6-term AR model 

identified using all 400 data points. 
27 

0 5 10 15 20 25 30 −2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Time, sec

Gust angle, deg

Measured

Predicted

20 20.5 21 21.5 22 22.5 23 23.5 24 24.5 25 −2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Time, sec

Gust angle, deg

Measured

Predicted

Figure A3.- Measured versus predicted gust time histories using a 6-term AR model identified 

using the first 500 data points (out of 1300) from HiLDA wing test. Airstream was excited 

at 2 Hz. 

 Expanded view of curves 

 between 20 and 25 sec 
28 

0 5 10 15 20 25 30 −1.5

−1

−0.5

0

0.5

1

Time, sec

Gust angle, deg

Measured

Predicted

20 20.5 21 21.5 22 22.5 23 23.5 24 24.5 25 −1.5

−1

−0.5

0

0.5

1

Time, sec

Gust angle, deg

Measured

Predicted

Figure A4.- Measured versus predicted gust time histories using a 6-term AR model identified 

using the first 500 data points (out of 1300) from HiLDA wing test. Airstream was excited 

at 8 Hz. 

 Expanded view of curves 

 between 20 and 25 sec 
29 

0 5 10 15 20 25 30 −2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Time, sec

Gust angle, deg

Measured

Predicted

20 20.5 21 21.5 22 22.5 23 23.5 24 24.5 25 −1.5

−1

−0.5

0

0.5

1

1.5

Time, sec

Gust angle, deg

Measured

Predicted

Figure A5.- Measured versus predicted gust time histories using a 6-term AR model obtained 

using the first 500 data points (out of 1300) from HiLDA wing test. Airstream was excited by 

sweeping gust frequency from 0.5 to 10 Hz. 

 

 Expanded view of curves 

 between 20 and 25 sec 
30 

0 2 4 6 8 10 12 14 16 18 20 −1.5

−1

−0.5

0

0.5

1

1.5

Time, sec

Force, lb

Imposed

Predicted

Figure A6.- Imposed versus predicted disturbance time histories using a 2-term AR model 

updated at each time step using RLS. 
31 

0 2 4 6 8 10 12 14 16 18 20 −1.5

−1

−0.5

0

0.5

1

1.5

Time, sec

Force, lb

Imposed

Predicted

Figure A7.- Imposed versus predicted disturbance time histories using a 6-term AR model 

updated at each time step using RLS. 
32 

0 2 4 6 8 10 12 14 16 18 20 −6

−4

−2

0

2

4

6

8

Time, sec

Force, lb

Imposed

Predicted

Figure A8.- Imposed versus predicted disturbance time histories using a 6-term AR model 

updated at each time step using RLS. 
33 

0 2 4 6 8 10 12 14 16 18 20 −10

−5

0

5

10

15

20

Time, sec

Force, lb

Imposed

Predicted

Figure A9.- Imposed versus predicted disturbance time histories using a 16-term AR model 

updated at each time step using RLS. 
34 

Appendix B 

State-Space Equations for 3-DOF Mass-Spring-Dashpot System 

 

 The second-order matrix equations of motion for the 3-degree-of-freedom mass-springdashpot system shown in figure 2 are given by 

 

1 1 12 2 1 12 2 1 1

2 2 2 23 3 2 2 23 3 2 2

33 3 3 3 3 3 3 3

00 0 0

0 0

00 0 0

100100

010010

001001

m x cc c x kk k x f

m x c cc c x k kk k x f

mx c c x k k x f

u

⎡ ⎤⎧ ⎫ ⎡ ⎤⎧ ⎫ ⎡ ⎤⎧ ⎫ ⎧ ⎫ +− +−

⎢ ⎥⎢ ⎥⎢ ⎥ ⎪⎪ ⎪⎪ ⎪⎪⎪⎪ ⎨⎬ ⎨⎬ ⎨⎬⎨⎬ +− + − +− + − =

⎪⎪ ⎪⎪ ⎪⎪⎪⎪ ⎣ ⎦⎩ ⎭ ⎣ ⎦⎩ ⎭ ⎣ ⎦⎩ ⎭ ⎩ ⎭ − −

⎡ ⎤

⎢ ⎥ =

⎣ ⎦

 

 

 

[ ][ ] { }

{ } [ ]{ }

1

2

3

2

1

2

3

ud uc

d

ud

ud ud

I I Bu

uc uc

uc

uc

⎧ ⎫

⎪ ⎪

⎪ ⎪

⎪ ⎪ ⎧ ⎫ = = ⎡ ⎤ ⎨⎬ ⎨ ⎬ ⎣ ⎦ ⎪ ⎪ ⎩ ⎭

⎪ ⎪

⎪ ⎪ ⎪ ⎪ ⎩ ⎭

(B1) 

where the external forces acting on the system are divided into disturbance forces udi and control 

forces uci for convenience in controls simulation work. The matrix [B2] is an input influence 

matrix indicating the locations and types of force inputs. Equation B1 can be written compactly 

as 

 [M ]{  x Cx Kx f B u } + + == [ ]{ } [ ]{ } { } [ 2 ]{ } (B2) 

 The corresponding first-order equations describing the continuous-time state-space model 

for the system are given by the state equation 

 {X} = + [AX Bu c c ]{ } [ ]{ }  (B3) 

and the output equation 

 

 {y C X Du } = + [ c c ]{ } [ ]{ } (B4) 

where the state matrix [Ac], state matrix [Bc], and state vector {X} are given by 

 [ ] [ ] [ ] [ ] [ ] { } { }

{ } 11 1

2

0 0

; ; c c

I x

A BX

MK MC MB x −− −

⎡ ⎤ ⎡⎤ ⎧ ⎫ = == ⎢ ⎥ ⎢⎥ ⎨ ⎬ − − ⎡ ⎤⎡ ⎤ ⎡ ⎤ ⎣ ⎦ ⎣⎦ ⎣ ⎦⎣ ⎦ ⎣ ⎦ ⎩ ⎭  (B5) 
35 

The outputs (measured physical quantities) are related to the state variables and the inputs by the 

output equation which can be written as 

 { } [ ] { }

{ } dv a [ ]{ } x

y CC Cx

x

⎧ ⎫

= + ⎨ ⎬ ⎩ ⎭

  (B6) 

where Cd, Cv, and Ca are 3×3 output influence matrices for displacement, velocity, and 

acceleration, respectively, and have the form 

 [ ] [] [] 1 11

22 2

33 3

; ;

d va

d d vv aa

dv a

c cc

C c Cc Cc

ccc

⎡ ⎤⎡ ⎤⎡ ⎤

⎢ ⎥⎢ ⎥⎢ ⎥ = ==

⎣ ⎦⎣ ⎦⎣ ⎦

 (B7) 

These matrices describe the relationship between the vectors {x}, , {x x   } { } and the measurement 

vector {y}. If the measurement unit is a physical unit (e.g., displacement, velocity, or 

acceleration) the corresponding term is unity; otherwise there is a conversion factor from 

physical units to the measurement units. 

Solving equation B2 for {x} and substituting into equation B6 gives 

 [ ] [ ] { }

{ } [ ] {} {} {} 11 1

dv a 2

x

y C C C MC x MK x MB u

x

−− − ⎧ ⎫

= +− − + ⎨ ⎬ ⎡ ⎤ ⎡⎤⎡ ⎤⎡ ⎤ ⎣ ⎦ ⎣⎦⎣ ⎦⎣ ⎦ ⎩ ⎭

  (B8a) 

or 

 { } { }

{ } { } 11 1

da va a 2

x

y C CM K C CM C CM B u

x

−− − ⎧ ⎫

=− − + ⎡ ⎤⎡ ⎤ ⎣ ⎦⎣ ⎦ ⎨ ⎬ ⎩ ⎭  (B8b) 

or 

 { } { } { } 1 11

da va a 2 y C CM K C CM C X CM B u − −− =− − + ⎡ ⎤⎡ ⎤ ⎣ ⎦⎣ ⎦ (B8c) 

which can be written compactly as 

 {y C X Du } = + [ c c ]{ } [ ]{ } (B9) 

where [Cc ] and [Dc ] are the output and direct transmission matrices and are given by 
36 

 

[ ]

[ ]

1 1

1

2

c da va

c a

C C CM K C CM C

D CM B

− −

−

=− − ⎡ ⎤ ⎣ ⎦

= ⎡ ⎤ ⎣ ⎦

 (B10) 

The direct transmission matrix[Dc ]is zero if accelerations are not included in the output 

measurements {y} (in which case [Ca ] = [0] ), or if accelerations are measured but there are no 

external forces (disturbance and control) acting at those locations (in which case some of the 

diagonal elements of Iud and Iuc will be zero). 

 Special cases of the general equations given above are easily obtained. For example, the 

number, type, and combination of input and output quantities can be varied by deleting 

appropriate rows of the output equation and/or zeroing appropriate diagonal elements in Iud, Iuc, 

Cd, Cv, and Ca. 
37 

References 

1. Juang, J.-N., and Eure, K. W.: Predictive Feedback and Feedforward Control for Systems 

with Unknown Disturbances, NASA/TM-1998-208744, December 1998. 

2. Kvaternik, R. G.; Juang, J-N.; and Bennett, R. L.: Exploratory Studies in Generalized 

Predictive Control for Active Aeroelastic Control of Tiltrotor Aircraft. NASA TM-2000- 

210552, October 2000. 

3. Kvaternik, R. G.; Piatak, D. J.; Nixon, M. W.; Langston, C. W.; Singleton, J. D.; Bennett, R. 

L.; and Brown, R. K.: An Experimental Evaluation of Generalized Predictive Control for 

Tiltrotor Aeroelastic Stability Augmentation in Airplane Mode of Flight. J. Amer. Hel. Soc., 

Vol. 27, No. 3, July 2002, pp. 198-208. 

4. Nixon, M. W.; Langston, C. W.; Singleton, J. D.; Piatak, D. J.; Kvaternik, R. G.; Corso, L. 

M.; and Brown, R. K.: Hover Test of a Soft-Inplane Gimballed Tiltrotor Model, J. Amer. 

Helicopter Soc., Vol. 48, Jan 2003, pp. 63-66. 

5. Nixon, M. W.; Langston, C. W.; Singleton, J. D.; Piatak, D. J.; Kvaternik, R. G.; Corso, L. 

M.; and Brown, R.: Aeroelastic Stability of a Four-Bladed Semi-Articulated Soft-Inplane 

Tiltrotor Model. Presented at the AHS 59th Annual Forum and Technology Display, May 6-8, 

2003, Phoenix, Arizona. In American Helicopter Society 59th Annual Forum Proceedings. 

6. Vartio, E.; Shimko, A.; Tilmann, C. P.; and Frick, P. M.: Structural Model Control and Gust 

Load Alleviation for a SensorCraft Concept. Presented at 46th AIAA/ASME/ASCE/AHS/ASC 

Structures, Structural Dynamics, and Materials Conference, Austin, TX, April 2005. Paper No. 

AIAA-2005-1946. 

7. MATLAB Reference Guide. The MathWorks, Inc., Natick, MA, August 1992. 

8. Clark, D. W. (Ed): Advances in Model-Based Predictive Control. Oxford University Press, 

New York, 1994. 

9. Clarke, D. W.; Mohtadi, C.; and Tuffs, P. S.: Generalized Predictive Control – Parts I and II. 

Automatica, Vol. 23, No. 2, 1987, pp. 137-160. 

10. Juang, J-N., and Phan, M. Q.: Identification and Control of Mechanical Systems. Cambridge 

University Press, New York, NY, 2001, ISBN 0-521-78355-0. 

11. Juang, J-N.: Applied System Identification. Prentice-Hall, Englewood Cliffs, NJ, 1994. 

12. Eure, K. W.; and Juang, J.-N.: Broadband Noise Control Using Predictive Techniques. 

NASA TM 110320, January 1997. 
38 

13. Juang, J.-N.: State-Space System Realization with Input- and Output-Data Correlation. 

NASA TP 3622, April 1997. 

14. Juang, J.-N.; and Phan, M. Q.: Deadbeat Predictive Controllers. NASA TM 112862, 

 May 1997. 

15. Eure, K. W.: Adaptive Predictive Feedback Techniques for Vibration Control. Ph.D. 

Dissertation, Virginia Polytechnic Institute and State University, May 1998. 

16. Phan, M. Q.; and Juang, J.-N.: Predictive Controllers for Feedback Stabilization. 

Journal of Guidance, Control, and Dynamics, Vol. 21, No. 5, Sept.-Oct. 1998, pp. 747-753. 

17. Juang, J.-N.; and Eure, K. W.: Predictive Feedback and Feedforward Control for Systems 

with Unknown Disturbances. NASA/TM-1998-208744, December 1998. 

18. Phan, M. Q.; Horta, L. G.; Juang, J.-N.; and Longman, R. W.: Identification of Linear 

Systems by an Asymptotically Stable Observer. NASA TP 3164, June 1992. 

19. Juang, J.-N.; and Phan, M. Q.: Recursive Deadbeat Controller Design. NASA TM 112863, 

May 1997 (also, Journal of Guidance, Control, and Dynamics, Vol. 21, No. 5, Sept.-Oct. 1998, 

pp. 747-753). 

20. Phan, M. Q.; and Juang, J-N: Predictive Feedback Controllers for Stabilization of Linear 

Multivariable Systems. AIAA Guidance, Navigation and Control Conference, San Diego, CA, 

July 29-31, 1996. 

21. Åström, K. J.; and Wittenmark, B.: Adaptive Control. Addison-Wesley Publishing 

Company, New York, 1995. 
39 

Table 1.- Reduction Factors for HiLDA Wing Model Moment Responses 

Using Embedded Feedforward

Gust Frequency SBI2 SBY STI2 

2 Hz 667 10 34 

8 Hz 75 33 30 

Swept from 0.5 – 10 Hz 31 7 4 

Table 2.- Reduction Factors for HiLDA Wing Model Moment Responses 

Using Explicit Feedforward 

Gust Frequency SBI2 SBY STI2 

2 Hz 1000 30 170 

8 Hz 500 100 60 

Swept from 0.5 – 10 Hz 625 40 4 
40 

Figure 1.- Block diagram of GPC identification and control procedure. 

Figure 2.- Three-degree-of-freedom mass-spring-dashpot system. 
41 

0 0.5 1 1.5 2 2.5 3 −2

−1

0

1

2

Time, sec

Force, lb

Disturbance at mass 1

Disturbance at mass 2

Figure 3.- Steady-state sinusoidal disturbances imposed on 3-DOF math model 

at masses 1 and 2 (ud1 and ud2). 

0 5 10 15 −2

−1

0

1

2

Time, sec

Force, lb

(a) Piecewise varying frequency and amplitude 

0 5 10 15 −1

−0.5

0

0.5

1

Time, sec

Force, lb

(b) Continuously varying frequency and amplitude 

Figure 4.- Non-steady-state disturbances imposed on 3-DOF math model at mass 2 (ud2). 
42 

0 5 10 15 −40

−20

0

20

40

60

80

Time, sec

Force, lb

0 5 10 15 −300

−200

−100

0

100

200

300

Time, sec

Acceleration, in/sec2

Open Loop

Closed Loop

Control Input 

Responses 

Figure 5.- Feedback only for 3 DOF system with steady-state disturbances (Control: uc3; 

Output: 1 x ; Disturbance: ud1 and ud2). 
43 

0 5 10 15 −40

−20

0

20

40

60

80

Time, sec

Force, lb

0 5 10 15 −300

−200

−100

0

100

200

300

Time, sec

Acceleration, in/sec2

Open Loop

Closed Loop

Control Input 

Responses 

Figure 6.- Feedback with embedded feedforward for 3 DOF system with steady-state 

disturbances (Control: uc3; Output: 1 x ; Disturbance: ud1 and ud2). 
44 

0 5 10 15 −40

−20

0

20

40

60

80

Time, sec

Force, lb

0 5 10 15 −300

−200

−100

0

100

200

300

Time, sec

Acceleration, in/sec2

Open Loop

Closed Loop

Control Input 

Responses 

Figure 7.- Feedback with explicit feedforward of both steady-state disturbances for 3 DOF 

system (Control: uc3; Output: 1 x ; Disturbance: ud1 and ud2). 
45 

0 5 10 15 −50

−40

−30

−20

−10

0

10

20

Time, sec

Force, lb

0 5 10 15 −300

−200

−100

0

100

200

300

Time, sec

Acceleration, in/sec2

uc1

uc3

Open Loop

Closed Loop

Responses 

Control Inputs 

Figure 8.- Feedback with explicit feedforward of both steady-state disturbances for 3 DOF 

system (Control: uc1 and uc3; Output: 1 x ; Disturbance: ud1 and ud2). 
46 

0 5 10 15 −50

−25

0

25

50

Time, sec

Force, lb

0 5 10 15 −300

−150

0

150

300

Time, sec

Acceleration, in/sec2

0 5 10 15 −300

−150

0

150

300

Time, sec

Acceleration, in/sec2

Open Loop

Closed Loop

Open Loop

Closed Loop

uc1

uc3

Control Inputs 

x2 acceleration response 

x1 acceleration response 

Figure 9.- Feedback with explicit feedforward of both steady-state disturbances for 3 DOF 

system (Control: uc1 and uc3; Output: 1 x and 2 x ; Disturbance: ud1 and ud2) 
47 

0 5 10 15 −2

−1

0

1

2

Time, sec

Force, lb

0 5 10 15 −3

−2

−1

0

1

2

3

Time, sec

Acceleration, in/sec2

Open Loop

Closed Loop

Responses 

Control Input 

Figure 10.- Feedback with embedded feedforward of disturbance with piecewise varying 

amplitude and frequency (Control: uc3; Output: 1 x ; Disturbance: ud2). 

 
48 

0 5 10 15 −2

−1

0

1

2

Time, sec

Force, lb

0 5 10 15 −3

−2

−1

0

1

2

3

Time, sec

Acceleration, in/sec2

Open Loop

Closed Loop

Responses 

Control Input 

Figure 11.- Feedback with explicit feedforward of disturbance with piecewise varying amplitude 

and frequency (Control: uc3; Output: 1 x ; Disturbance: ud2). 
49 

0 5 10 15 −1

−0.5

0

0.5

1

Time, sec

Force, lb

uc1

uc3

0 5 10 15 −3

−2

−1

0

1

2

3

Time, sec

Acceleration, in/sec2

Open Loop

Closed Loop Responses 

Control Inputs 

Figure 12.- Feedback with embedded feedforward of disturbance with piecewise varying 

amplitude and frequency (Control: uc1 and uc3; Output: 1 x ; Disturbance: ud2). 
50 

0 5 10 15 −1

−0.5

0

0.5

1

Time, sec

Force, lb

uc1

uc3

0 5 10 15 −3

−2

−1

0

1

2

3

Time, sec

Acceleration, in/sec2

Open Loop

Closed Loop Responses 

Control Inputs 

Figure 13.- Feedback with explicit feedforward of disturbance with piecewise varying amplitude 

and frequency (Control: uc1 and uc3; Output: 1 x ; Disturbance: ud2). 
51 

0 5 10 15 −1

−0.5

0

0.5

1

Time, sec

Force, lb

0 5 10 15 −3

−2

−1

0

1

2

3

Time, sec

Acceleration, in/sec2

Open Loop

Closed Loop

uc1

uc3

uc3

Control Inputs 

Responses 

Figure 14.- Feedback with explicit feedforward of disturbance with piecewise varying amplitude 

and frequency (Control: uc1, uc2, and uc3; Output: 1 x ; Disturbance: ud2). 
52 

0 5 10 15 −1

−0.5

0

0.5

1

Time, sec

Force, lb

0 5 10 15 −3

−2

−1

0

1

2

3

Time, sec

Acceleration, in/sec2

Open Loop

Closed Loop

Control Input 

Responses 

Figure 15.- Feedback with embedded feedforward of disturbance with continuously varying 

frequency and amplitude (Control: uc3; Output: 1 x ; Disturbance: ud2). 
53 

0 5 10 15 −1

−0.5

0

0.5

1

Time, sec

Force, lb

0 5 10 15 −3

−2

−1

0

1

2

3

Time, sec

Acceleration, in/sec2

Open Loop

Closed Loop

Control Input 

Responses 

Figure 16.- Feedback with explicit feedforward of disturbance with continuously varying 

frequency and amplitude acting at mass 2 (Control: uc3; Output: 1 x ; Disturbance ud2). 
54 

0 5 10 15 −1

−0.5

0

0.5

1

Time, sec

Force, lb

uc1

uc3

0 5 10 15 −3

−2

−1

0

1

2

3

Time, sec

Acceleration, in/sec2

Open Loop

Closed Loop

Control Inputs 

Responses 

Figure 17.- Feedback with explicit feedforward of disturbance with continuously varying 

frequency and amplitude (Control: uc1 and uc3; Output: 1 x ; Disturbance: ud2). 
55 

Figure 18.- HiLDA wing model in Langley Transonic Dynamics Tunnel. View looking

upstream.

Flow angle 

sensor 
56 

Figure 19.- Biplane vanes of TDT airstream oscillator system. View looking downstream toward

model in test section.

Biplane vanes of airstream oscillator system 
57 

Figure 20.- Schematic of HiLDA model wing planform showing arrangement of control surfaces 

and locations of strain-gage response sensors used for simulations. 
58 

0 5 10 15 20 25 30

-2500

-2000

-1500

-1000

-500

0

500

1000

1500

2000

2500

Time, sec

SBI2, in-lb

Open Loop

Closed Loop

Figure 21.- Time histories of open- and closed-loop vertical bending moment (SBI2) during 

 gust excitation at 2 Hz. Embedded feedforward. 
59 

0 5 10 15 20 25 30

-100

-50

0

50

100

150

200

250

Time, sec

SBY, in-lb

Open Loop

Closed Loop

Figure 22.- Time histories of open- and closed-loop fore-and-aft bending moment (SBY) during 

gust excitation at 2 Hz. Embedded feedforward. 
60 

0 5 10 15 20 25 30

-200

-150

-100

-50

0

50

100

150

200

Time, sec

STI2, in-lb

Open Loop

Closed Loop

Figure 23.- Time histories of open- and closed-loop torsional moment (STI2) during 

gust excitation at 2 Hz. Embedded feedforward. 
61 

0 5 10 15 20 25

10-4

10-2

100

102

104

106

108

1010

Frequency, Hz

SBI2, (in-lb)2/Hz

Open Loop

Closed Loop

Figure 24.- Open- and closed-loop power spectral densities for vertical bending moment (SBI2) 

during gust excitation at 2 Hz. Embedded feedforward. 
62 

0 5 10 15 20 25

10-4

10-2

100

102

104

106

Frequency, Hz

SBY, (in-lb)2/Hz

Open Loop

Closed Loop

Figure 25.- Open- and closed-loop power spectral densities for fore-and-aft bending moment 

(SBY) during gust excitation at 2 Hz. Embedded feedforward. 
63 

0 5 10 15 20 25

10-4

10-2

100

102

104

106

Frequency, Hz

STI2, (in-lb)2/Hz

Open Loop

Closed Loop

Figure 26.- Open- and closed-loop power spectral densities for torsional moment (STI2) 

during gust excitation at 2 Hz. Embedded feedforward. 
64 

Figure 27.- Time histories of closed-loop control inputs (CTL1-CTL5) during gust excitation 

 at 2 Hz. Embedded feedforward. 

0 5 10 15 20 25 30

-5

-4

-3

-2

-1

0

1

2

3

4

5

Time, sec

Control Surface Deflection Angle, deg

CTL1

CTL2

CTL3

CTL4

CTL5
65 

0 5 10 15 20 25 30

-400

-300

-200

-100

0

100

200

300

400

Time, sec

SBI2, in-lb

Open Loop

Closed Loop

Figure 28.- Time histories of open- and closed-loop vertical bending moment (SBI2) during 

 gust excitation at 8 Hz. Embedded feedforward. 
66 

0 5 10 15 20 25 30

-80

-60

-40

-20

0

20

40

60

80

Time, sec

SBY, in-lb

Open Loop

Closed Loop

Figure 29.- Time histories of open- and closed-loop fore-and-aft bending moment (SBY) 

 during gust excitation at 8 Hz. Embedded feedforward. 
67 

0 5 10 15 20 25 30

-100

-80

-60

-40

-20

0

20

40

60

80

100

Time, sec

STI2, in-lb

Open Loop

Closed Loop

Figure 30.- Time histories of open- and closed-loop torsional moment (STI2) during 

 gust excitation at 8 Hz. Embedded feedforward. 
68 

0 5 10 15 20 25

10-6

10-4

10-2

100

102

104

106

108

Frequency, Hz

SBI2, (in-lb)2/Hz

Open Loop

Closed Loop

Figure 31.- Open- and closed-loop power spectral densities for vertical bending moment (SBI2) 

during gust excitation at 8 Hz. Embedded feedforward. 
69 

0 5 10 15 20 25

10-6

10-4

10-2

100

102

104

106

Frequency, Hz

SBY, (in-lb)2/Hz

Open Loop

Closed Loop

Figure 32.- Open- and closed-loop power spectral densities for fore-and-aft bending moment 

(SBY) during gust excitation at 8 Hz. Embedded feedforward. 
70 

0 5 10 15 20 25

10-4

10-2

100

102

104

106

Frequency, Hz

STI2, (in-lb)2/Hz

Open Loop

Closed Loop

Figure 33.- Open- and closed-loop power spectral densities for torsional moment (STI2) 

 during gust excitation at 8 Hz. Embedded feedforward. 
71 

Figure 34.- Time histories of closed-loop control inputs (CTL1-CTL5) during gust 

excitation at 8 Hz. Embedded feedforward. 

0 5 10 15 20 25 30

-5

-4

-3

-2

-1

0

1

2

3

4

5

Time, sec

Control Surface Deflection Angle, deg

CTL1

CTL2

CTL3

CTL4

CTL5
72 

0 5 10 15 20 25 30

-2500

-2000

-1500

-1000

-500

0

500

1000

1500

2000

Time, sec

SBI2, in-lb

Open Loop

Closed Loop

Figure 35.- Time histories of open- and closed-loop vertical bending moment (SBI2) while 

sweeping gust frequency from 0.5 to 10 Hz. Embedded feedforward. 
73 

0 5 10 15 20 25 30

-150

-100

-50

0

50

100

150

Time, sec

SBY, in-lb

Open Loop

Closed Loop

Figure 36.- Time histories of open- and closed-loop fore-and-aft bending moment (SBY) while 

sweeping gust frequency from 0.5 to 10 Hz. Embedded feedforward. 
74 

0 5 10 15 20 25 30

-250

-200

-150

-100

-50

0

50

100

150

200

250

Time, sec

STI2, in-lb

Open Loop

Closed Loop

Figure 37.- Time histories of open- and closed-loop torsional moment (STI2) while sweeping 

gust frequency from 0.5 to 10 Hz. Embedded feedforward. 
75 

0 5 10 15 20 25

10-6

10-4

10-2

100

102

104

106

108

Frequency, Hz

SBI2, (in-lb)2/Hz

Open Loop

Closed Loop

Figure 38.- Open- and closed-loop power spectral densities for vertical bending moment (SBI2) 

while sweeping gust frequency from 0.5 to 10 Hz. Embedded feedforward. 
76 

0 5 10 15 20 25

10-4

10-2

100

102

104

106

Frequency, Hz

SBY, (in-lb)2/Hz

Open Loop

Closed Loop

Figure 39.- Open- and closed-loop power spectral densities for fore-and-aft bending moment 

(SBY) while sweeping gust frequency from 0.5 to 10 Hz. Embedded feedforward. 
77 

0 5 10 15 20 25

10-4

10-2

100

102

104

106

Frequency, Hz

STI2, (in-lb)2/Hz

Open Loop

Closed Loop

Figure 40.- Open- and closed-loop power spectral densities for torsional moment (STI2) while 

sweeping gust frequency from 0.5 to 10 Hz. Embedded feedforward. 
78 

0 5 10 15 20 25 30

-30

-20

-10

0

10

20

30

Time, sec

Control Surface Deflection Angle, deg

CTL1

CTL2

CTL3

CTL4

CTL5

Figure 41.- Time histories of closed-loop control inputs (CTL1-CTL5) while sweeping 

gust frequency from 0.5 to 10 Hz. Embedded feedforward. 
79 

0 5 10 15 20 25 30

-2500

-2000

-1500

-1000

-500

0

500

1000

1500

2000

2500

Time, sec

SBI2, in-lb

Open Loop

Closed Loop

Figure 42.- Time histories of open- and closed-loop vertical bending moment (SBI2) during 

gust excitation at 2 Hz. Explicit feedforward. 
80 

0 5 10 15 20 25 30

-100

-50

0

50

100

150

200

250

Time, sec

SBY, in-lb

Open Loop

Closed Loop

Figure 43.- Time histories of open- and closed-loop fore-and-aft bending moment (SBY) during 

gust excitation at 2 Hz. Explicit feedforward. 
81 

0 5 10 15 20 25 30

-200

-150

-100

-50

0

50

100

150

200

Time, sec

STI2, in-lb

Open Loop

Closed Loop

Figure 44.- Time histories of open- and closed-loop torsional moment (STI2) during 

gust excitation at 2 Hz. Explicit feedforward. 
82 

0 5 10 15 20 25

10-5

100

105

1010

Frequency, Hz

SBI2, (in-lb)2/Hz

Open Loop

Closed Loop

Figure 45.- Open- and closed-loop power spectral densities for vertical bending moment (SBI2) 

during gust excitation at 2 Hz. Explicit feedforward. 
83 

0 5 10 15 20 25

10-4

10-2

100

102

104

106

Frequency, Hz

SBY, (in-lb)2/Hz

Open Loop

Closed Loop

Figure 46.- Open- and closed-loop power spectral densities for fore-and-aft bending moment 

(SBY) during gust excitation at 2 Hz. Explicit feedforward. 
84 

0 5 10 15 20 25

10-4

10-2

100

102

104

106

Frequency, Hz

STI2, (in-lb)2/Hz

Open Loop

Closed Loop

Figure 47.- Open- and closed-loop power spectral densities for torsional moment (STI2) 

during gust excitation at 2 Hz. Explicit feedforward. 
85 

Figure 48.- Time histories of closed-loop control inputs (CTL1-CTL5) during gust 

excitation at 2 Hz. Explicit feedforward. 

0 5 10 15 20 25 30

-5

-4

-3

-2

-1

0

1

2

3

4

5

Time, sec

Control Surface Deflection Angle, deg

CTL1

CTL2

CTL3

CTL4

CTL5
86 

0 5 10 15 20 25 30

-400

-300

-200

-100

0

100

200

300

400

Time, sec

SBI2, in-lb

Open Loop

Closed Loop

Figure 49.- Time histories of open- and closed-loop vertical bending moment (SBI2) during 

 gust excitation at 8 Hz. Explicit feedforward. 
87 

0 5 10 15 20 25 30

-80

-60

-40

-20

0

20

40

60

80

Time, sec

SBY, in-lb

Open Loop

Closed Loop

Figure 50.- Time histories of open- and closed-loop fore-and-aft bending moment (SBY) during 

gust excitation at 8 Hz. Explicit feedforward. 
88 

0 5 10 15 20 25 30

-100

-80

-60

-40

-20

0

20

40

60

80

100

Time, sec

STI2, in-lb

Open Loop

Closed Loop

Figure 51.- Time histories of open- and closed-loop torsional moment (STI2) during 

gust excitation at 8 Hz. Explicit feedforward. 
89 

0 5 10 15 20 25

10-6

10-4

10-2

100

102

104

106

108

Frequency, Hz

SBI2, (in-lb)2/Hz

Open Loop

Closed Loop

Figure 52.- Open- and closed-loop power spectral densities for vertical bending moment (SBI2) 

during gust excitation at 8 Hz. Explicit feedforward. 
90 

0 5 10 15 20 25

10-6

10-4

10-2

100

102

104

106

Frequency, Hz

SBY, (in-lb)2/Hz

Open Loop

Closed Loop

Figure 53.- Open- and closed-loop power spectral densities for fore-and-aft bending moment 

(SBY) during gust excitation at 8 Hz. Explicit feedforward. 
91 

0 5 10 15 20 25

10-4

10-2

100

102

104

106

Frequency, Hz

STI2, (in-lb)2/Hz

Open Loop

Closed Loop

Figure 54.- Open- and closed-loop power spectral densities for torsional moment (STI2) 

during gust excitation at 8 Hz. Explicit feedforward. 
92 

Figure 55.- Time histories of closed-loop control inputs (CTL1-CTL5) during gust 

excitation at 8 Hz. Explicit feedforward. 

0 5 10 15 20 25 30

-5

-4

-3

-2

-1

0

1

2

3

4

5

Time, sec

Control Surface Deflection Angle, deg

CTL1

CTL2

CTL3

CTL4

CTL5
93 

0 5 10 15 20 25 30

-2500

-2000

-1500

-1000

-500

0

500

1000

1500

2000

Time, sec

SBI2, in-lb

Open Loop

Closed Loop

Figure 56.- Time histories of open- and closed-loop vertical bending moment (SBI2) while 

sweeping gust frequency from 0.5 to 10 Hz. Explicit feedforward. 
94 

0 5 10 15 20 25 30

-150

-100

-50

0

50

100

150

Time, sec

SBY, in-lb

Open Loop

Closed Loop

Figure 57.- Time histories of open- and closed-loop fore-and-aft bending moment (SBY) 

while sweeping gust frequency from 0.5 to 10 Hz. Explicit feedforward. 
95 

0 5 10 15 20 25 30

-250

-200

-150

-100

-50

0

50

100

150

200

250

Time, sec

STI2, in-lb

Open Loop

Closed Loop

Figure 58.- Time histories of open- and closed-loop torsional moment (STI2) while sweeping 

gust frequency from 0.5 to 10 Hz. Explicit feedforward. 
96 

0 5 10 15 20 25

10-6

10-4

10-2

100

102

104

106

108

Frequency, Hz

SBI2, (in-lb)2/Hz

Open Loop

Closed Loop

Figure 59.- Open- and closed-loop power spectral densities for vertical bending moment (SBI2) 

while sweeping gust frequency from 0.5 to 10 Hz. Explicit feedforward. 
97 

0 5 10 15 20 25

10-5

100

105

Frequency, Hz

SBY, (in-lb)2/Hz

Open Loop

Closed Loop

Figure 60.- Open- and closed-loop power spectral densities for fore-and-aft bending moment 

(SBY) while sweeping gust frequency from 0.5 to 10 Hz. Explicit feedforward. 
98 

0 5 10 15 20 25

10-4

10-2

100

102

104

106

Frequency, Hz

STI2, (in-lb)2/Hz

Open Loop

Closed Loop

Figure 61.- Open- and closed-loop power spectral densities for torsional moment (STI2) while 

sweeping gust frequency from 0.5 to 10 Hz. Explicit feedforward. 
99 

0 5 10 15 20 25 30

-30

-20

-10

0

10

20

30

Time, sec

Control Surface Deflection Angle, deg

CTL1

CTL2

CTL3

CTL4

CTL5

Figure 62.- Time histories of closed-loop control inputs (CTL1-CTL5) while sweeping 

gust frequency from 0.5 to 10 Hz. Explicit feedforward. 
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing

data sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or 

any other aspect of this collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate 

for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that 

notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently 

valid OMB control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY) 2. REPORT TYPE 3. DATES COVERED (From - To) 

4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

5c. PROGRAM ELEMENT NUMBER 

6. AUTHOR(S) 5d. PROJECT NUMBER 

5e. TASK NUMBER 

5f. WORK UNIT NUMBER 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION 

REPORT NUMBER 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSORING/MONITOR'S ACRONYM(S) 

11. SPONSORING/MONITORING

REPORT NUMBER 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

15. SUBJECT TERMS 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 

ABSTRACT 

18. NUMBER 

OF 

PAGES 

19a. NAME OF RESPONSIBLE PERSON 

a. REPORT b. ABSTRACT c. THIS PAGE 

19b. TELEPHONE NUMBER (Include area code) 

Standard Form 298 (Rev. 8-98) Prescribed by ANSI Std. Z39-18