INTEGRATED MODELING AND CONTROL OF FLEXIBLE AIRCRAFT WINGS
By
Dagmara Anna Wehr

A THESIS
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Electrical Engineering – Master of Science
2014

ABSTRACT
INTEGRATED MODELING AND CONTROL OF FLEXIBLE AIRCRAFT WINGS
By
Dagmara Anna Wehr
Structural control for vibration reduction has important applications in many research
areas, including the effect of earthquakes on buildings and aerodynamic forces on aircraft
stability and performance. Both passive and active control techniques have been implemented,
with the best solution usually involving a passive approach followed by an active one. This
thesis presents an integrated modeling and controller design approach. Modal Cost Analysis
(MCA) and Output Covariance Constraint (OCC) control are used to reduce a high-order
aeroelastic wing model to establish the best controller for the reduced-order model, with a
constraint on the covariance of the vibration outputs. MCA seeks to keep the modes that have
the highest contribution to a given cost function. Using iterations on the two processes will
allow a lower-order controller to be designed and result in the same performance.
The OCC and MCA methods and their respective algorithms are presented, and an
approach to integrate the two procedures is given. NASA’s model used in this thesis is applied
to the MCA and OCC algorithms using MATLAB. A 40th-order wing model is derived. The
model reduction technique initially reduces the system to a 12th order one. A simulation of the
OCC algorithm is performed on the reduced-order model and applied to the full-order model.
The controller resulting in the best closed-loop performance is shown to significantly reduce the
vibrations due to wind. A corresponding weighting matrix used in OCC is then used for a
second round of MCA to further reduce the model to an 8th order model. A lower-order
controller designed for this second model is shown to similarly reduce the output vibrations.

ACKNOWLEDGEMENTS

Over the past two years, I have received support and encouragement from a number of
individuals. First, I would like to thank my advisor, Dr. Guoming (George) Zhu for his guidance
and support in my research and thesis. His guidance made this a thoughtful, manageable effort.
I would also like to thank Dr. Xiaobo Tan and Dr. Hassan Khalil for serving on my guidance
committee.
I would also like to thank Dr. Percy Pierre and Dr. Xiaobo Tan for their support in
helping me attend and be successful in my graduate studies.
Finally, I would like to thank my family for their support. They believed in me and made
me believe in myself during rough periods.

iii

TABLE OF CONTENTS

LIST OF TABLES ...........................................................................................................................v
LIST OF FIGURES ....................................................................................................................... vi
CHAPTER 1: INTRODUCTION ....................................................................................................1
1.1 Structural Control for Vibration Reduction .................................................................1
1.1.1 Passive Control .............................................................................................1
1.1.2 Active Control ...............................................................................................2
1.1.3 Combining Passive and Active Control ........................................................3
1.2 Related Research ..........................................................................................................5
1.2.1 Related Applications .....................................................................................5
1.2.2 Model Reduction Methods ............................................................................6
1.3 Objective and Problem Statement ................................................................................7
1.4 Thesis Outline ..............................................................................................................9
CHAPTER 2: REVIEW OF METHODS ......................................................................................11
2.1 Output Covariance Constraint (OCC) Control ..........................................................11
2.1.1 Output Covariance Constraint Control Problem .........................................11
2.1.2 OCC Algorithm ...........................................................................................13
2.2 Model Reduction Using Modal Cost Analysis ..........................................................14
2.3 Integration of Controller Design and Model Reduction ............................................15
CHAPTER 3: AIRCRAFT WING MODEL AND SIMULATION ..............................................17
3.1 Modeling Description ................................................................................................17
3.2 Model Reduction ........................................................................................................20
3.3 Simulation Study........................................................................................................22
3.3.1 Simulation Set-Up.......................................................................................22
3.3.2 Simulation Results ......................................................................................24
3.3.3 Iteration on OCC and MCA ........................................................................28
CHAPTER 4: CONCLUSION ......................................................................................................35
APPENDIX ....................................................................................................................................36
BIBLIOGRAPHY ..........................................................................................................................40

iv

LIST OF TABLES

Table 3.1. Modal Cost Analysis ....................................................................................................21
Table 3.2. Modal Cost Analysis for Q Loop 1 ..............................................................................30

v

LIST OF FIGURES

Figure 2.1.

Loop Diagram ......................................................................................................16

Figure 3.1.

Locations of Measured and Control Outputs ..........................................................20

Figure 3.2a. Impulse Response of Evaluation and Design Models: 6th Input to 4th Output ........22
Figure 3.2b. Impulse Response of Evaluation and Design Models: 6th Input to 5th Output ........22
Figure 3.3.

Output Variance vs. Iteration Number ....................................................................25

Figure 3.4.

Output Variance vs. Input Covariance ....................................................................26

Figure 3.5a. Impulse Response of the Open and Closed-Loop Systems from Wind Input to
5th Output .................................................................................................................26
Figure 3.5b. Impulse Response of the Open and Closed-Loop Systems from Wind Input to
7th Output .................................................................................................................27
Figure 3.6a. Open-Loop Wing Deformation with Time ..............................................................28
Figure 3.6b. Closed-Loop Wing Deformation with Time ...........................................................28
Figure 3.7a. Impulse Response of Evaluation and 8th Order Design Models: 6th Input to 4th
Output ......................................................................................................................30
Figure 3.7b. Impulse Response of Evaluation and 8th Order Design Models: 4th Input to 5th
Output ......................................................................................................................31
Figure 3.8a. Output Variance vs. Iteration Number for Q Loop 1 ..............................................32
Figure 3.8b. Output Variance vs. Input Covariance for Q Loop 1 ..............................................32
Figure 3.9a. Open vs. Closed-Loop Impulse Response of Wind to 5th Output for Q Loop 1 .....33
Figure 3.9b. Open vs. Closed-Loop Impulse Response of Wind to 7th Output for Q Loop 1 .....33
Figure 3.10. Closed-Loop Wing Deformation with Time for Q Loop 1 .......................................34

vi

CHAPTER 1: INTRODUCTION

1.1

Structural Control for Vibration Reduction
Structural control for vibration reduction has been employed in various fields including

civil and aerospace structures to protect against the possibly devastating effects of vibrations
caused by natural forces such as wind and earthquakes. With the invention of structures with
newer designs, using modern materials and shapes, comes the increased need for methods that
reduce vibrations caused by environmental or other forces. To achieve this, the control design
should consider both passive and active control. The following sections discuss different types
of control, as well as their traditional and recent applications.
1.1.1

Passive Control

Passive control can refer to two kinds of vibration suppression. One describes passive
actuators, traditionally used in building designs, and the other to “structural redesign”, in which
system parameters, such as those corresponding to damping, are modified to improve
performance [1], [2].
In the first case, passive control describes the use of devices that do not use external
energy sources to reduce vibrations of structures. As opposed to active controllers, they reduce
vibrations without any sensors or feedback. Passive controllers, such as tuned liquid dampers
and tuned mass dampers (TMDs), have been widely used for civil engineering applications,
especially as a means to counteract the effects of earthquakes. These devices have the ability to
control high-amplitude forces reliably and cost-effectively [1]. Purely passive TMDs, also
referred to as tuned vibration absorbers (TVAs), have been applied to tall buildings, bridges,
offshore platforms, and pipelines, each of which experience one of more forces from winds,

1

earthquakes, or sea waves. Aerospace applications include helicopter rotor hubs, and noise
radiation from structural resonance vibrations such as cabin noise [3]. [4] also discusses passive
control methods for various sources of aircraft vibration suppression. The author presents a
method consisting of “fuselage double-partition constructions with acoustic insulation coating”
for noise vibration reduction. He also describes actuators mounted on the fuselage skin or frame
to minimize vibrations [4].
In the second definition of passive control, system redesign is used in structures as a
means to minimize the active control needed [2]. Therefore, it involves system design rather
than direct control techniques.

Often, changes to the plant can improve the system’s

performance beyond what can be achieved by feedback control, and these changes can be
implemented more easily than designing a new controller. Plant and controller design are a
related process, because controllers are designed to modify the dynamics of a system which
could otherwise be modified by a change in plant parameters [2]. For this reason, [2] proposes
an iterative technique to simultaneously design the plant and controller. When the plant and
controller are redesigned in subsequent iterations, the goal is to maintain the desired performance
while minimizing the control energy [2]. This structural redesign approach is also applied in [5]
to a cantilever beam as a means to control vibrations of structures that resemble cantilevers.
1.1.2

Active Control

Active control is often necessary and allows going beyond passive control by efficiently
applying an appropriate control force at each moment to counteract structural vibrations, using
sensors and actuators. This type of control allows moderating multiple vibration modes at once
with one device. So, it is an effective method for controlling structures with several vibration
modes that contribute to unwanted structural dynamics. Furthermore, active controllers are

2

versatile and can be completely stable if collocated control is used. However, active control
techniques are expensive, as a result of a higher level of technology and required maintenance.
Additionally, often large control devices and power supply systems are needed [1]. Reliability
problems arise in civil structures if the power to a building is cut, which is possible during an
earthquake [6].
The use of active mass dampers for buildings is mentioned in [7]. [1] also applies active
controllers in the form of active mass dampers to a 3 story building model and a real footbridge,
as well as an active bracing system to the building model. Applications of active TMDs can
include classical control, fuzzy logic, and neural networks [3]. Active control in the form of
active velocity feedback is demonstrated for a cantilever beam in [8] to reduce its first 3
vibration modes. [9] describes the use of active control to reduce structural vibrations in a
blended wing-body type aircraft for ride comfort and improved handling qualities.
1.1.3

Combining Passive and Active Control

Due to the advantages and disadvantages of both passive and active controllers, the two
techniques are often combined in the form of structural redesign (Section 1.1.1), semi-active
control, or hybrid control. Use of passive techniques is frequently used to supplement the
necessity of active controllers. Therefore, when considering control techniques, passive control
should be considered first, then improved upon using active control approaches.
In an aircraft, aeroelasticity is determined by the aircraft’s structural, damping, and mass
characteristics, not only by the external forces. Therefore, considering passive control in the
sense of model redesign is extremely beneficial. This way, an efficient and simpler active
controller can be designed to an already properly-damped system. Since plant and controller
design problems are not independent, basic changes in the plant can also often result in great

3

improvements in performance and robustness over just using feedback control [2]. In [2], an
iterative technique for plant and controller redesign is proposed to make control of the system
easier.
Semi-active control refers to using a passive device that can be controlled actively [6].
Semi-active, sometimes called adaptive, control devices have advantages over fully active
systems. Passive controllers have their advantages, but can lose their effectiveness over time. It
is therefore useful to be able to modify, or adapt, these devices actively to changing operating
and environmental conditions. The disadvantage of semi-active systems lies in their many
components which may need continual maintenance [3]. The use of semi-active TMDs is
mentioned for several applications in [3]. Semi-active control is also studied in the context of
building vibration control as a solution that offers higher energy-efficiency than active devices
and is more effective at reducing seismic vibrations than passive devices [10]. In [10], a
controller is developed to adjust the damping of a type of semi-active control called a
magnetorheological damper, which changes its damping based on changing the viscosity of a
fluid using a magnetic field. [6] also mentions the use of semi-active control in the form of
hydraulic dampers, fluid viscous dampers, variable-orifice dampers, variable-friction dampers,
electrorheological, and magetorheological devices.

The main advantage of semi-active

controllers is that they run on battery power and are not destabilizing [6]. [7] mentions the use of
a semi-active hydraulic damper on a bridge.
Both [6] and [7] discuss hybrid control, such as a hybrid mass damper in [8], which uses
a combination of passive actuators and active control devices. Since a passive device is used,
less control energy is needed for the active component. Along with semi-active control, hybrid
control techniques are more reliable against power failures, as they can support some level of

4

damage prevention over purely active control systems in this case [7]. Many active TMDs are
actually hybrid, consisting of passive devices as well. In the context of a building, the active
device would be used for small vibrations due to wind, and passive for stronger reactions caused
by earthquakes [3].
1.2

Related Research
1.2.1

Related Applications

This research is closely based on the work in [11]. In [11], the authors simultaneously
consider modeling and controller design of flexible aircraft wings. Modal Cost Analysis (MCA)
is used to reduce the high-order finite-dimension state-space model for a desired output
covariance, with the output consisting of displacement and rates of points on the wing. A
reduced-order model is used to design a dynamic output feedback controller using Output
Covariance Constraint (OCC) control. The controller is then evaluated with the full-order model
to confirm handling performance requirements have been met, and the wing fluttering motion
has been almost completely eliminated [11].
Active flutter suppression of the NASA Benchmark Active Controls Technology (BACT)
wing using a gain-scheduled controller utilizing a linear fractional transformation is studied in
[12]. The goal of [12] is to increase disturbance rejection and the stable operating range. Since
the wing section varies with Mach and dynamic pressure changes, the controller is made as a
function of the two changing dynamics. As a comparison, a linear controller is also designed
using D-K iteration. The performance is analyzed with maximum singular value plots and time
simulations. The gain-scheduled controller is found to be stable over the specified operating
range, and have significant performance improvement [12].
Generalized Predictive Control (GPC) for a wind-tunnel wing model for flutter

5

suppression is presented in [13]. The control technique minimizes a cost function consisting of
the sum of the mean square output and the square rate of change of the control input. The
algorithm uses the model output to predict an arbitrary input’s effect on plant dynamics. A
BACT plant is reduced using an auto-regressive moving-average model. The cost function is
modified to account for the output being invariant to low-frequency inputs caused by a drift in
the control input. Simulations including the wind-tunnel test confirm flutter suppression robust
to a wide range of conditions and modeling errors [13].
1.2.2

Model Reduction Methods

There are several model reduction techniques that have been developed for largedimensioned systems. In [14], the author mentions several such methods appropriate to his
application of large space structures including the Internal Balancing Method, Component Cost
Analysis (CCA), Modal Cost Analysis (MCA), and the Canonical Correlation Analysis Method.
The Internal Balancing Method is based on measures of controllability and observability,
removing the least controllable and least observable parts to form a reduced-order model. CCA
and MCA (which is a special case of CCA), and their variations, Weighted Component Cost
Analysis (WCCA) and Weighted Modal Cost Analysis (WMCA), take into account the system’s
different components to reduce them separately. CCA considers any system as made up of
components, and reduces it based on the different components’ contributions to some quadratic
cost functional. MCA is a special case of CCA, which uses the system described by its modes.
The weighted versions of these two methods are used when input dynamics are considered,
instead of just white noise input. The Canonical Correlation Analysis Method also uses a
system’s components, as well as their interactions, to rank the components by their canonical
correlation coefficients and delete the smallest ones to form the reduced-order model [14]. The

6

methods applied in [14] are WMCA and a combination of CCA with Canonical Correlation
Analysis. MCA is also used in [15] and [16] for aerospace applications.
Several other model reduction methods are considered in [6] and [17], which are divided
into three types: nodal truncation, singular value decomposition (SVD)-based methods, and
Krylov-based methods.

Nodal truncation algorithms, including Guyan reduction, truncates

nodes (from finite element methods) based on their influence on the structure’s mass and
stiffness matrices. The standard approach used in industry for civil structure applications is
modal truncation, since higher modes usually have less influence on the system response than
lower modes. However, this method has high computational costs. The SVD-based methods,
including balanced model reduction, rely on the system’s Hankel singular values and balancing
of the system to determine the contribution of each state.
computationally demanding for large-scale systems.

However, this method is

Finally, the Krylov methods rely on

moment matching to iteratively match moments of original and lower-order models [6]. [6]
presents these methods in the context of civil structure modeling and control, and [17] mentions
an application to the International Space Station.
1.3

Objective and Problem Statement
With the development of modern materials, their use in new airplane wing designs has

been studied in recent years [11], [18]. Using materials that are more lightweight and flexible
leads to more energy-efficient aerodynamic designs, improving cruise efficiency and lowering
drag [11]–[13], [18]. Materials such as composites are less rigid than in previous-generation
designs, while preserving the necessary load-carrying capacity.

The Boeing 787 Dreamliner

aircraft is an example of the highly-flexible wing structure. The concept is that the flexible wing
surfaces are shaped during flight to change the angle of attack in a way that drag can be reduced,

7

resulting in lower fuel burn during cruise [18]. With increased structural flexibility, however,
come increased effects of aerodynamic forces on the aircraft’s stability and performance [11].
The objective of this work was to combine model reduction and controller design to
actively suppress the wing vibrations of a NASA aircraft model. Model reduction and controller
design was done together and iteratively to produce the best controller for the reduced-order
model, and again, the best reduced model for the chosen controller.

Output Covariance

Constraint (OCC) control was used for controller design, and Modal Cost Analysis (MCA) was
used for model reduction. It was expected that the controller designed using OCC for the design
model could significantly reduce the wing vibrations due to a random wind input. In addition,
through an iteration of model reduction and controller design, the closed-loop performance could
be matched with a lower-order controller than the initial reduced-order model.
MCA and OCC was an appropriate choice of model reduction and controller design
methods due to their interaction, to achieve the goal of better performance by coupling the two
design steps. MCA is effective for reducing large models that can be expressed by a sum of their
vibration modes. In MCA, a cost is found for each mode to an output covariance cost function,
and the lowest-cost modes are truncated. The performance requirements for aeroelastic aircraft
wings vibrating due to wind can also be expressed in terms of output covariance [11]. A
physical interpretation of OCC can be given by a bound on the peak magnitude of the time
response [15]. If the output covariance constraint is expressed by , and if there is a constraint
,
then the peak magnitude of the time response of the output
energy disturbances,

and , is bounded such that [15]:
.

8

in the presence of bounded

The OCC problem seeks to find a proper weighting matrix that leads to a controller which
minimizes the control effort while meeting constraints on output covariance [11].

This

weighting matrix is in turn used in MCA.
The work of this thesis is based on [11], which uses MCA and OCC to control vibrations
of a flexible wing structure. However, it goes a step further by iterating on the model reduction
and controller design steps for a second time, as mentioned previously.

By choosing an

appropriate weighting matrix found in the first iteration of the OCC control design, this
weighting is used to better determine the significant modes that contribute to the output. This
allows an even lower-order controller to be designed with very similar performance.
1.4

Thesis Outline
The rest of this thesis consists of three chapters. In Chapter 2, the methods used in the

research are presented. First, the definition and algorithm of the OCC problem is presented.
Next, MCA is discussed as a method for model reduction. Finally, the importance of coupling
model reduction and controller design is discussed.
Chapter 3 presents the application, mentioned above, of the methods given in Chapter 2.
The model of the aircraft is derived and the simulation parameters are given. Next, details of
how the reduced model was found using MCA is presented. A comparison of the full-order and
reduced-order model is shown with their impulse responses. Section 3.3 summarizes the steps
taken to carry out the OCC control design algorithm with the reduced-order model. The output
variances are plotted against the trace of the input (control) covariance, and a time response
comparing open and closed-loop systems is presented. Afterward, the iteration on MCA and
OCC is given, and the model is taken through a second round of model reduction and controller
design. Similar plots as in the first iteration are shown.

9

The third and final chapter presents the conclusions.

10

CHAPTER 2: REVIEW OF METHODS

2.1

Output Covariance Constraint (OCC) Control
2.1.1

Output Covariance Constraint Control Problem

Consider a continuous linear time-invariant system,
=

(2.1a)
(2.1b)
(2.1c)

,
with white noise signals

and

of intensities

and , respectively, and a strictly proper full-

.

(2.2a)
(2.2b)

order dynamic controller,

The closed-loop system that results from applying the controller (2.2) to the plant (2.1) is
(2.3a)
(2.3b)

,
where

(2.4a)
(2.4b)
(2.4c)

.
Substituting Equations (2.1) and (2.2) into (2.3), Equation (2.3) can be expressed as:

(2.5a)
.

(2.5b)

Then, let
.

11

(2.6)

If

represents the closed-loop controllability Gramian from the disturbance input

, then

satisfies
.

(2.7)

The following is the definition for the Output Covariance Constraint (OCC) problem: Find a
controller to minimize the cost function,
,
where

(2.8)

is a diagonal weighting matrix, subject to the system Equations (2.1) to (2.5) above,

and
,
where

(2.9)

is the dimension of the output y,
.

weighted

control

energy

),

, and

This means that the OCC problem is to minimize the

with

certain

constraints

on

the

block

diagonal

matrix

consisting of the output variances [19], [20].
Suppose that

is an optimal solution to the OCC problem. Then there exists a

semi-definite matrix,
,

(2.10)

such that
(2.11a)
(2.11b)
(2.11c)
(2.11d)
(2.11e)
(2.11f)
(2.11g)

,
where

is the closed-loop covariance of

12

and

is the covariance of

output

[19], [20].
2.1.2

OCC Algorithm

An algorithm for finding the optimal solution of the OCC problem discussed above is
described as follows [19], [20]:
Given

initial

, and tuned parameters

and

let

and proceed to 1).
1) Compute

and

2) Compute

by solving (2.11c) and (2.11e).

and

by solving
(2.12)
(2.13)

.
3) Compute

by solving
.

4) Compute

,

(2.14)

.

5) Let

and
.

If

, where

step 2 with

is defined below.

is some specified tolerance, stop. Otherwise, go back to

.

When the algorithm above ends, the dynamic controller
(2.15a)
(2.15b)

,
is an optimal solution to the OCC problem [19], [20].
The operator
[15]: Let

used in step 5 of the algorithm above is defined in the following way

represent a real, symmetric matrix that can be written in the Schur decomposition

form as follows:
13

,
where

and

respectively, and

(2.16)

are diagonal and consist of the positive and non-positive eigenvalues of

,

is a unitary orthogonal matrix. Then,

.

(2.17)

A MATLAB program used in Section 3.3 for the simulation implementing the OCC
algorithm is given in the Appendix.
2.2

Model Reduction Using Modal Cost Analysis
There are several ways to design a reduced-order controller for a high-order system, for

which a full-order controller has implementation limitations due to a computational limitation of
the real-time microcontroller. One method is to directly design a fixed-order controller. Another
approach is to design a full-order controller and then conduct a controller order reduction to get
the reduced-order controller. Lastly, the method used in this thesis is to first reduce the system
model, and then design a full-order controller for the reduced model [16]. One advantage of
using this approach is that the designed reduced-order control at least stabilizes the reducedorder model. However, the quality of the reduced-order controller is heavily dependent on the
accuracy of the model reduction. One such model-reduction approach presented here, Modal
Cost Analysis (MCA), can be used to reduce a high-order system model, especially for a
structure system model, when it is expressed in its modal coordinates. MCA takes into account
the modes of the system, finding and truncating those modes that have the lowest contribution to
the output. This method calculates each modal contribution,

, to the weighted cost function
,

where

is the expectation operator,

is the output vector,
14

(2.18)
is the number of modes in the

system, and

is the weighting matrix for modal cost analysis which could be from Equation

(2.10) [15]. The function corresponds to the weighted output variance with respect to the white
noise input [16].
Recall the system Equations (2.1). If the state-space parameters were transformed into
their modal coordinates,

, then the Lyapunov equation,
,

gives the solution , where

(2.19)

is as in Equation (2.6). Then, the modal costs can be calculated

by
,
where

(2.20)

is the number of modes. The resulting modes can be ranked as
,

where

now corresponds to the highest modal cost, not necessarily to the first mode. The

modes corresponding to the lowest modal costs

are then truncated to give a lower-order model

[11].
2.3

Integration of Controller Design and Model Reduction
When designing a controller for a reduced-order model, model reduction and controller

design is an integrated process. Using this method of designing a controller for a reduced-order
model, the neglected dynamics of the system will have an effect on how well the controller
works with the actual, full-order model [16]. Since the controller designed for the reduced
model may not be the best for the full-order model that it has to control, an iterative approach is
used in [15], [16] with an application to space structures. During this process, the best controller
is designed for a reduced-order model, after which the model reduction process is repeated to
produce the best model for the current controller [16].

15

The weighting matrix

introduced in Equation (2.10) of section 2.1.1 is the same one

used in both MCA and the beginning of OCC.

is a diagonal matrix, with relative higher-

valued elements emphasizing certain system outputs. In this way, the choice of
modal costs, and therefore the resulting model reduction.
calculations.

affects the

also affects the output variance

Combining the model reduction and controller design processes allows

to

represent the most significant outputs and retain the most significant modes [11].
A diagram of what will be referred to as the
2.1. As the figure shows, the chosen

Loop in Chapter 3 is presented in Figure

from the OCC controller design is used in another

iteration of model reduction using MCA.
A software package in MATLAB including the algorithms for the integration of model
reduction and controller design are used in [15], [16] and is described in [21].

Figure 2.1. Q Loop Diagram

16

CHAPTER 3: AIRCRAFT WING MODEL AND SIMULATION

3.1

Modeling Description
The aircraft wing model is described in [11] by combining the aeroelastic and rigid body

aircraft dynamics with wind turbulence and actuator dynamics.

The aircraft dynamics are

described by the following equations:
(3.1a)
,
where

(3.1b)

are the bending and torsional displacements of the 20

aeroelastic modes,

are the rigid body states (angle of attack and pitch rate),

elevator deflection angle,

is the flap deflection, and

is the

is the slat deflection. Equation (3.1)

can be written in state-space form as follows:
(3.2)
where

,

,

(3.3a)

.

(3.3b)

A turbulence model is then added to the above to include the effect of random wind gust
dynamics. It is assumed to be in the following form:

,
with states

, random wind gust

modeled by zero-mean white noise with intensity
17

(3.4a)
(3.4b)
, and

the total random wind input to the wing model. Combining the turbulence model with the
wing dynamics model gives the following representation:
.
The input

that comes from the actuators follows a system:
(3.6a)
(3.6b)

,
with states

(3.5)

and control command .

To obtain the complete open-loop state-space system, Equations (3.4-3.6) are combined
to give
,
where

(3.7)

,

(3.8a)

(3.8b)

.

Then, the performance outputs

and measurement outputs

,
where

(3.8c)

are given by
(3.9a)
(3.9b)

is the zero-mean white noise of the sensors, with intensity . Combining the previous

two equations with Equation (3.7) gives the open-loop state-space representation of the whole,
full-order system, which was provided by NASA. This system will be referred to as the
evaluation model:

18

(3.10a)
(3.10b)
(3.10c)
where

and
. The first two states of

represent angle of

attack and pitch rate, the next 20 states the 10 bending and 10 torsional states, the following 20
states the rates of the 10 bending and 10 torsional states, and lastly the 3 turbulence states.
However, to simplify the simulation, the rigid body mode and turbulence states were removed
after transforming the model to its modal coordinates (discussed in Section 3.2), leaving only the
40 aeroelastic states. This simplified model assumes all forces are applied to the wings.
After simplification, the control outputs are

where

bending displacement measurements (in feet) at 5 locations and
displacement measurements (in radians) at those locations.

are

are the torsional

The measurement outputs

are the bending and torsional deflections along with their
corresponding rates at 10 locations. Figure 3.1 below, provided by NASA, shows the locations
of the measurements and control outputs – 10 equally spaced points along the wing with the last
one being on the wingtip for the measurements, and the last 5 of these coinciding with the
control outputs [11]. The controls include an elevator, 6 flaps, and 6 slats. The random wind
gust
5% of

is

ft/s, which gives

. Therefore,

.

is the random measurement noise, chosen as

for each measurement.

19

Figure 3.1. Locations of Measured and Control Outputs
3.2

Model Reduction
To use MCA (described in Chapter 2), the model in Equation (3.10) is transformed into

modal coordinates using a transformation matrix . The system matrices are then described by
(3.11a)
(3.11b)
(3.11c)
(3.11d)
(3.11e)
with

block diagonal – the first 3x3 block corresponding to turbulence, the next 2x2 block

corresponding to the rigid body mode, and lastly the 20 aeroelastic modes. So, for the simplified
model, the first 5x5 block was removed.
The algorithm described in Chapter 2 was used to perform MCA on the full-order model
in Equation (3.10) in modal coordinates. The weighting matrix

was chosen to be equal to the

identity matrix. 6 modes with the highest modal costs were kept, producing a 12th order reducedorder model, which will be referred to as the design model:

.
Table 3.1 below gives the modal costs of the evaluation model [21].

20

(3.12a)
(3.12b)
(3.12c)

Table 3.1. Modal Cost Analysis

Figures 3.2a and 3.2b show a comparison of the impulse responses from the sixth input to
the fourth output and from the fourth input to the fifth output, respectively. The plots show that
the responses after model reduction vary only slightly in the first half second, and are nearly
identical afterward. Therefore, the 12th order design model is a good lower-order representation
of the full-order evaluation model.

21

100
evaluation model
design model

80
60

Amplitude

40
20
0
-20
-40
-60
-80
0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Time (sec)

Figure 3.2a. Impulse Response of Evaluation and Design Models: 6th Input to 4th
Output
300
evaluation model
design model

200

Amplitude

100
0
-100
-200
-300
0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

Time (sec)

Figure 3.2b. Impulse Response of Evaluation and Design Models: 4th Input to 5th
Output
3.3

Simulation Study
3.3.1

Simulation Set-up

The evaluation model described in Section 3.1 and the design model described in 3.2 was
input into a MATLAB script that performed the OCC algorithm (given in the Appendix). First,

22

the variables that do not depend on the weighting matrix

,

and

from Equations (2.11c and

2.11e), were calculated. The program includes two loops: the outer loop iterates on the output
covariance constraint

, which is decreased, and an inner loop that corresponds to the OCC

algorithm presented in Section 2.1.2.
At the beginning of the outer loop,

is calculated based on the following equation:
,

where

,

of the outer loop, and

is initially 0.95 and is decreased by 0.05 each iteration

is solved from the following equation:
,

where

(3.13)

(

and

(3.14)

is solved by
.

This loop continues as long as

(3.15)

is above 0. The purpose of this loop is to find a sequence of

OCC controllers with control effort from low to high.
The inner loop goes through the OCC algorithm, calculating

and

. For each

inner loop iteration, the convergence condition given in Step 3 of the OCC algorithm is tested
chosen as 10-6. If it fails, a new

with

Step 5. For this model, the values
all

is calculated according to the equation for
and

conditions besides the last one corresponding to

with this updated

in

were used, resulting in convergence for
. Then, the inner loop continues

value. If the condition in Step 3 passes (convergence), then the inner loop is

exited.
The block diagonal matrix consisting of the output variances for the evaluation model is
then computed:
,
23

(3.16)

where

,

, and the closed-loop state

is the

solution to
,

(3.17)

where
(3.18)
and

and

following from Equation (2.5) in Chapter 2 as:
,

(3.19)

where
(3.20a)
.
For plotting the results, the input covariance

(3.20b)

, which corresponds to the control energy, is

also calculated:
,
where

and

of the outer loop,
3.3.2

(3.21)

is the value from the last iteration of the inner loop. At the end

is decremented, and the loop iterates, with the last iteration at

.

Simulation Results

The following two plots were generated by the MATLAB program described in Section
3.3.1. Figure 3.3 shows the output variances with respect to the design model and the output
variances with respect to the evaluation model, against iteration of the outer loop. The 18
iterations that converged are included (excluding when

), and the 18 points each

correspond to a designed controller. The design model is plotted with a solid line and the
evaluation model in a dashed line. However, the overlap of the design and evaluation model

24

variances in most of the plot results in the appearance of only the design model being plotted.
The variances of each of the 10 outputs are plotted separately. Therefore, the plot includes 20
lines, 10 with respect to each of the models. Only 5 outputs can be easily seen however, because
the remaining ones are much closer to 0.
4

Output Variances

3.5
3
2.5
2
1.5
1
0.5
0
0

2

4

6

8

10

12

14

16

18

Iteration

Figure 3.3. Output Variances vs. Iteration Number
Figure 3.4 shows the same output variances, but plotted against the input covariance,
described in Equation (3.21). Again, each point corresponds to a controller designed for a given
output covariance constraint. The shape of the plot and the spacing of the points reveal the
tradeoff between output variance and control energy used to achieve that output variance. As
output variance of each output decreases, input covariance increases.

This means that as

performance increases, control energy increases, and at a faster rate. The ideal controller would
then be chosen by considering the desired performance and level of control, with the best
controller generally corresponding to one of the points in the lower left side of the plot (lower
output variance and input covariance).

25

4

Output Variances

3.5
3
2.5
2
1.5
1
0.5
0
0

0.5

1

1.5

2

2.5

3

3.5
-3
x 10

Input Covariance

Figure 3.4. Output Variances vs. Input Covariance
Figures 3.5a and 3.5b show a comparison of the open-loop and closed-loop impulse
response from the 16th iteration of the outer loop of the controller design. This controller was
chosen based on a reasonable tradeoff between performance and control energy. The response is
from the wind gust input to the fifth and seventh outputs, respectively. The plots verify that the
designed controller reduces the amplitude of the aeroelastic modes to nearly zero, such that the
amplitude cannot easily be seen on the plot, compared to the much larger open loop response.
10
closed loop
open loop

8

Amplitude

6
4
2
0
-2
-4
-6
0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Time (sec)

Figure 3.5a. Impulse Response of the Open and Closed-Loop Systems from Wind
Input to 5th Output
26

0.25
closed loop
open loop

0.2

Amplitude

0.15
0.1
0.05
0
-0.05
-0.1
0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Time (sec)

Figure 3.5b. Impulse Response of the Open and Closed-Loop Systems from
Wind Input to 7th Output
Figures 3.6a and 3.6b show the bending displacements of each measured point on the
wing for the open and closed-loop situations, respectively. The controller used in the closedloop case is again that from the 16th iteration of the outer loop of the controller design (also used
in Section 3.3.3 for the next round of MCA and OCC). If Figure 3.1 corresponds to the side
view of the wing with no forces applied, the following plots show the shape of the wing for 3
time steps to an impulse response on the wind input. As in Figure 3.1, the rightmost point
corresponds to that of the wing tip, and the 5 rightmost points are the control outputs. The lines
in Figure 3.6b are nearly horizontal on the same scale, showing that the bending vibrations have
been significantly reduced.

27

0.2

Displacement

0
t=1 s
t=3 s
t=5 s

-0.2
-0.4
-0.6
-0.8
-1
1

2

3

4
5
6
7
Measurement Point Along Wing

8

9

10

8

9

10

Figure 3.6a. Open-Loop Wing Deformation with Time
0.2

Displacement

0
-0.2
-0.4

t=1 s
t=3 s
t=5 s

-0.6
-0.8
-1
1

2

3

4
5
6
7
Measurement Point Along Wing

Figure 3.6b. Closed-Loop Wing Deformation with Time
3.3.3

Iteration on OCC and MCA

In the previous section, an optimal controller was found for the 12th order design model,
and was shown to significantly reduce the amplitude of the impulse response with the evaluation
model. The next step was to iterate once again on the model reduction algorithm and controller
design to achieve a similar performance output with a lower-order controller. The goal of this
section is to describe the iteration process and show the results of the second simulation.
28

18 weighting matrices,

, were found after completing the OCC algorithm described in

Section 3.3.1. As discussed in the last section, one of the latter
controller, was chosen. As mentioned in Section 2.2, this
of the system. Therefore, a modified version of this

values, that from the 16th

relates to the most significant modes

was used in the second round of model

reduction, which will be referred to as Q Loop 1 (Q Loop 0 is the initial model reduction and
controller design. See Figure 2.1). In order to capture the effect of the new

, it was multiplied

by a large scalar (specifically, 108), so that the effect of the 10 control outputs could be seen over
the weighting of the measured outputs, which remained unchanged in the second iteration of
MCA (Q Loop 1). Using this scaled

value, the same 6 modes had the highest modal costs.

However, the scaling shows that the relative difference between the 4 highest-cost modes and the
two immediately lower ones increased. This indicates that the new

value further allowed the

model to be reduced to its 4 most significant modes, resulting in an 8th order system. The results
of MCA are shown below in Table 3.2. The table shows the larger difference between modes 5
and 6 as compared to 2, 4, 3, and 1.

29

Table 3.2 Modal Cost Analysis for Q Loop 1

Again, two impulse responses were plotted, as in Figures 3.2a and 3.2b, this time to
compare the evaluation model to the 8th order reduced model. The result is plotted in Figures
3.7a and 3.7b, and shows the validity of this further reduced model.
100
evaluation model
design model

80
60
Amplitude

40
20
0
-20
-40
-60
-80
0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Time (sec)

Figure 3.7a. Impulse Response of Evaluation and 8th Order Design Models: 6th
Input to 4th Output

30

300
evaluation model
design model

200

Amplitude

100
0
-100
-200
-300
0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Time (sec)

Figure 3.7b. Impulse Response of Evaluation and 8th Order Design Models: 4th
Input to 5th Output
The OCC algorithm was then run using the original evaluation model and the new design
model. Instead of beginning the algorithm with

equal to identity, the

value that was found

before scaling for the model reduction was used to start. The same plots as in the first

loop

were produced, which are shown below in Figures 3.8a and 3.8b. The dashed lines correspond to
the evaluation model. The plots show that the output variance values decrease as before.
Additionally, for all but the topmost output, the variance values are below or very close to those
in the first OCC iteration, Q Loop 0.

31

4

Output Variances

3.5
3
2.5
2
1.5
1
0.5
0
0

2

4

6

8

10

12

14

16

18

Iteration

Figure 3.8a. Output Variances vs. Iteration Number for Q Loop 1
4

Output Variances

3.5
3
2.5
2
1.5
1
0.5
0
0

0.5

1

1.5

2

2.5

Input Covariance

3
-3

x 10

Figure 3.8b. Output Variances vs. Input Covariance for Q Loop 1
Lastly, Figures 3.9a and 3.9b show the open-loop vs. closed-loop responses made again
with the 16th controller of Q Loop 1. The plots show that the vibrations are once again reduced
significantly.

32

10
closed loop
open loop

8

Amplitude

6
4
2
0
-2
-4
-6
0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Time (sec)

Figure 3.9a. Open vs. Closed-Loop Impulse Response of Wind to 5th Output for Q
Loop 1
0.25
closed loop
open loop

0.2

Amplitude

0.15
0.1
0.05
0
-0.05
-0.1
0

1

2

3

4

5

6

7

8

Time (sec)

Figure 3.9b. Open vs. Closed-Loop Impulse Response of Wind to 7th Output for Q
Loop 1
Figure 3.10 corresponds to Figure 3.6b, and shows the bending displacement of the whole
wing at 3 time points using the 16th controller designed in Q Loop 1. By comparing to the openloop plot in Figure 3.6a, it can be seen that the bending displacements are reduced equally
effectively using the lower-order 8th order controller as with the 12th order controller from Q
Loop 0.

33

0.2

Displacement

0
t=1 s
t=3 s
t=5 s

-0.2
-0.4
-0.6
-0.8
-1
1

2

3

4
5
6
7
Measurement Point Along Wing

8

3.10. Closed-Loop Wing Deformation with Time for Q Loop 1

34

9

10

CHAPTER 4: CONCLUSION

In this work, control for vibration reduction is applied to modern, elastic aircraft wings.
With increased structural flexibility comes a greater interaction with aerodynamic forces, which
increases the need and challenge for an effective controller design. The aircraft physical system
is often very high-order and must be reduced. This work presented an integrated controller
design and model reduction method to realize the natural coupling of the two steps. MCA was
used to first reduce the order of the model to a certain level that would allow a practical
controller to be designed and implemented. After that, the OCC control design algorithm was
implemented to produce a controller for the reduced-order model, which was verified for the
full-order model. A second iteration on the model reduction and controller design process used
the weighting matrix found in the first iteration to show that an even lower-order model can be
used in the OCC algorithm, and a new controller was again verified for the full-order model.
Using plots of the variance values of the bending displacement and rates outputs, as well
as a comparison between open and closed-loop impulse responses to a random wind input, it was
shown that the controllers designed in both Q-loop iterations were successful in significantly
reducing wing vibrations, with the lower-order controller performing equally effectively.

35

APPENDIX

36

OCC PROGRAM CODE

%
load design model
['load cdmodel_' int2str(Q_loop) '.mat'];
eval(ans);
%
load evaluation model
load cemodel_0.mat
%
load some parameters
if Q_loop==0
load weight.mat
else
['load weight' int2str(Q_loop) '.mat'];
eval(ans);
end
%
for subsequent Q loop iterations
if Q_loop>0
Q=Qred; %
Q is that from the last Q loop iteration
end
%
define more parameters
epsi = 1*10^-6;
Beta = 0.8; %
0<Beta<1
alpha = .2; %
alpha>0
gamma = 0.95;
count = 1; %
gamma iteration
N = 500;
num_iter(1:19)=0;
Q1 = Q;
%
ititialize to 0
Yd_i = zeros(10,19);
Ye_i = zeros(10,19);
Yd_norm = zeros(1,18);
Ye_norm = zeros(1,18);
[sr,sc]=size(EMd0);
[X_tilde,L,G] =
care(EAd0',EMd0',EDd0*W*EDd0',V,zeros(sc,40),eye(length(EAd0)));
F = X_tilde*EMd0'/(V);
X_OL = lyap(EAd0,EDd0*W*EDd0');
Yb_OL=diag(diag(ECd0*X_OL*ECd0'));
while gamma >= 0

37

Yb_bar = eye(10)*norm(gamma*Yb_OL) + diag(diag(ECd0*X_tilde*ECd0'));
%
recalculate Q with new Yb_bar (after initial iteration when we start
with Q=I
if count ~= 1
M = Q + alpha*(Yb-Yb_bar);
%
Schur decomposition on M:
[U,T] = schur(M);
%
then make negative diag == 0
for k=1:10
if T(k,k)<0
T(k,k)=0;
end
end
P = U*T*U';
Q = Beta*Q+(1-Beta)*P;
eval(['Q' int2str(count) '=Q;']);
end
for j=1:N %
however long you want to allow to converge (get below
epsi)
num_iter(count) = num_iter(count)+1;
%
for each changing Q:
[sr,sc]=size(EBd0);
[K,blah,blahh] =
care(EAd0,EBd0,ECd0'*Q*ECd0,R,zeros(sr,13),eye(length(EAd0)));
G = -inv(R)*EBd0'*K;
X = lyap(EAd0+EBd0*G,F*V*F');

%
for design model
Yb = diag(diag(ECd0*(X_tilde+X)*ECd0'));
test = norm((Yb-Yb_bar)*Q);
if test < epsi
disp(['Converged after ', num2str(num_iter(count)), ' iterations
when gamma = ', num2str(gamma)]);
break % BREAK 1
elseif j==N
disp(['Did not converge for gamma = ', num2str(gamma)]);
break % BREAK 1
end
%
Else update Q
M = Q + alpha*(Yb-Yb_bar);
%
Schur decomposition on M:
[U,T] = schur(M);
%
then make negative diag == 0
for k=1:10
if T(k,k)<0
T(k,k)=0;
end
end
38

P = U*T*U';
Q = Beta*Q+(1-Beta)*P;
Qnorm(j,count) = norm(Q); % supposed to go up
eval(['Q' int2str(count) '=Q;']); % could be outside loop
end
%
BREAK 1 goes to here
%
For design model
Yd_i(:,count) = diag(ECd0*(X_tilde+X)*ECd0'); % gives column vector
Yd_norm(count) = norm(diag(diag(ECd0*(X_tilde+X)*ECd0')));
%
For evaluation model
%
Output Covariance Ye for evaluation model
Ac = EAd0+EBd0*G-F*EMd0;
A_CL = [EA0 EB0*G ; F*EM0 Ac];
[sr,sc]=size(ECd0);
C_CL = [EC0 zeros(10,sc)];
D_CL = [ED0 zeros(length(ED0),40) ; zeros(length(EAd0),1) F];
W_bar = [W zeros(1,40) ; zeros(40,1) V];
Xe_CL = lyap(A_CL,D_CL*W_bar*D_CL');
%
Save every closed-loop system
['save CLsys' int2str(Q_loop) '_' int2str(count) ' A_CL C_CL D_CL'];
eval(ans);
%
Input Covariance U_CL
Cu_CL = [zeros(13,length(EA0)) G];
U_CL = Cu_CL*Xe_CL*Cu_CL';
U_CL_plot(count) = trace(U_CL);
Ye_i(:,count) = diag(C_CL*(Xe_CL)*C_CL'); % gives column vector
Ye_norm(count) = norm(diag(diag(C_CL*(Xe_CL)*C_CL')));
gamma = gamma - 0.05;
count = count+1; % gamma changes
end
%
Plot results
figure
plot(Yd_i(:,1:(end-1))','.-')
hold on
plot(Ye_i(:,1:(end-1))','--.')
xlabel('Iteration')
ylabel('Output Variances')
figure
plot(U_CL_plot(1:(end-1)),Yd_i(:,1:(end-1))','.-')
hold on
plot(U_CL_plot(1:(end-1)),Ye_i(:,1:(end-1))','--.')
xlabel('Input Covariance')
ylabel('Output Variances')

39

BIBLIOGRAPHY

40

BIBLIOGRAPHY
[1] C. Moutinho, A. Cunha, and E. Caetano, “Implementation of active and passive control
systems in laboratorial and real structures,” S3T, 2010.
[2] K. M. Grigoriadis, G. Zhu, and R. E. Skelton, “Optimal Redesign of Linear Systems,” J.
Dyn. Syst. Meas. Control, vol. 118, no. 3, pp. 598–605, Sep. 1996.
[3] J. Q. Sun, M. R. Jolly, and M. A. Norris, “Passive, Adaptive and Active Tuned Vibration
Absorbers—A Survey,” J. Mech. Des., vol. 117, no. B, pp. 234–242, Jun. 1995.
[4] P. Gardonio, “Review of Active Techniques for Aerospace Vibro-Acoustic Control,” J.
Aircr., vol. 39, no. 2, pp. 206–214, Apr. 2002.
[5] K. Adachi, K. Sakamoto, and T. Iwatsubo, “Redesign of closed loop system for integrated
design of structure and its vibration control system,” in Proceedings of the 1999 IEEE
International Conference on Control Applications, 1999, 1999, vol. 1, pp. 80–85 vol. 1.
[6] E. Gildin, “Model and controller reduction of large-scale structures based on projection
methods,” The University of Texas at Austin, Austin, Texas, 2006.
[7] T. T. Soong and B. F. Spencer Jr, “ACTIVE, SEMI-ACTIVE AND HYBRID CONTROL
OF STRUCTURES,” in Proceedings of the Twelfth World Conference on Earthquake
Engineering, Auckland, New Zealand, 2000.
[8] B. J. Grobbelaar, G. D. Wood, and C. P. Constancon, “Application of Modal Control to
Flexible Structures,” NO J., pp. 18–24, Jul. 1992.
[9] “Active Control of Flexible 2020 Aircraft - Transport - Research & Innovation - European
Commission,” European Commission. [Online]. Available:
http://ec.europa.eu/research/transport/projects/items/acfa_2020_en.htm. [Accessed: 13-Jul2014].
[10] C. M. D. Wilson and M. M. Abdullah, “Structural Vibration Reduction Using Fuzzy
Control of Magnetorheological Dampers,” in Structures Congress 2005, American Society
of Civil Engineers, pp. 1–12.
[11] S. S.-M. Swei, G. G. Zhu, and N. T. Nguyen, “Integrated Model Reduction and Control of
Aircraft with Flexible Wings,” presented at the AIAA Guidance, Navigation, and Control
(GNC) Conference, Boston, MA.
[12] J. M. Barker, G. J. Balas, and P. . Blue, “Active flutter suppression via gain-scheduled
linear fractional control,” in American Control Conference, 1999. Proceedings of the 1999,
1999, vol. 6, pp. 4014–4018 vol.6.
[13] P. Haley and D. Soloway, “Generalized predictive control for active flutter suppression,”
IEEE Control Syst., vol. 17, no. 4, pp. 64–70, Aug. 1997.

41

[14] J. H. Kim, “Model reduction, sensor/actuator selection and control system redesign for
large flexible structures,” Ph.D., Purdue University, United States -- Indiana, 1991.
[15] G. Zhu, K. M. Grigoriadis, and R. E. Skelton, “Covariance Control Design for Hubble
Space Telescope,” J. Guid. Control Dyn., vol. 18, no. 2, pp. 230–236, Apr. 1995.
[16] G. Zhu and R. E. Skelton, “Integrated Modeling and Control for the Large Spacecraft
Control Laboratory Experiment Facility,” J. Guid. Control Dyn., vol. 17, no. 3, pp. 442–
450, Jun. 1994.
[17] A. C. Antoulasг and D. C. Sorensen, “Approximation of large-scale dynamical systems: An
overview,” 2001.
[18] N. Nguyen and J. Urnes, “Aeroelastic Modeling of Elastically Shaped Aircraft Concept via
Wing Shaping Control for Drag Reduction,” in AIAA Atmospheric Flight Mechanics
Conference, American Institute of Aeronautics and Astronautics.
[19] G. Zhu, M. A. Rotea, and R. Skelton, “A convergent algorithm for the output covariance
constraint control problem,” SIAM J. Control Optim., vol. 35, no. 1, pp. 341–361, 1997.
[20] G. Zhu, “L_2 AND L_INFINITY MULTIOBJECTIVE CONTROL FOR LINEAR
SYSTEMS,” Purdue University, West Lafayette, IN, 1992.
[21] G. Zhu and R. E. Skelton, “INTEGRATED MODELING AND CONTROL SOFTWARE,”
Dynamic Systems Research, INC, West Lafayette, IN, Sep. 1992.

42