Topics in the Structural Dynamics and Control of an Articulated Micro Air Vehicle

Material Information

Topics in the Structural Dynamics and Control of an Articulated Micro Air Vehicle
Regisford, Sean James Sidney, 1977-
Place of Publication:
[Gainesville, Fla.]
University of Florida
Publication Date:
Physical Description:
1 online resource (133 p.)

Thesis/Dissertation Information

Doctorate ( Ph.D.)
Degree Grantor:
University of Florida
Degree Disciplines:
Aerospace Engineering
Mechanical and Aerospace Engineering
Committee Chair:
Lind, Richard C.
Committee Co-Chair:
Kurdila, Andrew
Committee Members:
Ifju, Peter
Bloomquist, David G.
Graduation Date:


Subjects / Keywords:
Aircraft ( jstor )
Aircraft wings ( jstor )
Coordinate systems ( jstor )
Equations of motion ( jstor )
Fuselages ( jstor )
Kinetic energy ( jstor )
Matrices ( jstor )
Modeling ( jstor )
Optimal control ( jstor )
Simulations ( jstor )
Mechanical and Aerospace Engineering -- Dissertations, Academic -- UF
Electronic Thesis or Dissertation
bibliography ( marcgt )
theses ( marcgt )
Aerospace Engineering thesis, Ph.D.


This work presents a methodology to determine of accurate computer models of articulated micro air vehicles. Prior to this work, much of the innovation in design of a micro air vehicle was based on empirical data. This emphasis on empirical data necessitates a time intensive iterative process where a design is built and tested, then modifications are made and the process starts anew. This work looks to alleviate some of the load of this iteration by producing accurate models that can be inserted into the process at greatly reduced cost in terms of time. The modeling is accomplished by the application of Lagrange's equations to create a set of equations that incorporate both the rigid and flexible degrees of freedom associated with the materials that Micro Air Vehicles (MAVs) employ. These equations of motion also include the generalized forces that result from the body forces due to gravity and external forces acting as a result of the interaction between the structure and the air through which it passes. These equations are formulated in a manner that accomplishes the above while allowing for the reconfiguration of the system. The equations of motion are then transformed into a form amenable for use in the solution of an optimal control problem. ( en )
General Note:
In the series University of Florida Digital Collections.
General Note:
Includes vita.
Includes bibliographical references.
Source of Description:
Description based on online resource; title from PDF title page.
Source of Description:
This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Thesis (Ph.D.)--University of Florida, 2008.
Adviser: Lind, Richard C.
Co-adviser: Kurdila, Andrew.
Statement of Responsibility:
by Sean James Sidney Regisford

Record Information

Source Institution:
Rights Management:
Copyright by Sean James Sidney Regisford. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Embargo Date:
LD1780 2008 ( lcc )


This item has the following downloads:

Full Text

a32,jkT (P)k a23,nma21,npT ()qp al2,rq 0 or12,tT (a), a21,va23,vwT (3, 23 ? ia32,zy Uopjop, dV +
a I (ai 23 & /
a32,ijT (/)k a23,mka21,mn. T (a)jp p 12,qp q 1al2,rs~T (a),, 21,ut~ 23,uvT (/) 23 ,z UoN dV +

q ,iJ 32,jkT (3)nk 23,nm. .21,npT ( p 012,rq ,sr12,stT (a),tu 21,vua 23,vwT (3),w 23 -? 32,zyH2-CB3, ,i dV
f ,3232,jyrT 0)B3,qP a32,iq d
q9 ,i2,jfT (32),mk 23,nm. .21,npT(oq) pa 1l2,rq0 C ,srQ12,stT (a),tu 21,vua23,vwT (),w 23r -2 a 32,zy N ,jiUop,z dV +
q ,i"32,ijT ()kj a23,mka21,mnT (a),pn 12,qp0 0 ,ql2,rsT (a),st a21,ut 23,uT (/3), 23w 2 wz N3,i ,yz dV(B-12)

For remaining terms replace o0D, with OD1, and O" where appropriate.

BB_ [a12] [T (a)] [a2]T [a23] [T (3)] [32 ]T [N3]
S[BB] [a12] [T (a)] [a2] [023] [T (3)] [a32] [N3] dV BB2 [012] [T (a)] [a21T [23] [T (P)] [32]T [N3] dV
BB3 [a12] [T (a)] [.1]T [.23 [T (3)] [.32] [N3]
Examining a representative term:

/ BBEf [a12] [T (a)] [a1] [a3] [T (3)] []32T [N3] dV

rB-1 Hl,r, Czsra12,stT (o),tu 21,vut23,vwT (/),wx N3,yz dV +
Tl-CB2,n21QnpT (a) O2,rq ,sr 12,stT () ,t 21,vua23,vwT ( N3,yx d +
rB2 H2O21p (O),20 Jq N3,0 dV1
.1 V
'B2H2l 2I,2npT (),p a12,rq OC,sra12,stT (o),tu 21,vu023,vwT (3),wx N3,yx dV +

TH2-CI T( k 3,nmL i21,npT (a o l2,rq 0 ,sr1,sstT (),tu 21,vu~23,vwT (),wx N3,yz dV +

T(/),k Oa23,nma21,npT (o)qp a 12,rqO,sra12,stT (a),tu 21,va23vwT(/),wx, V iN3,yz dV

qf,i T (/),mk 023,nm 21,npT (o),qp 2,rqO~,srL12,stT (),tu a21,vua23,vwT (), N3,jN3,yz (L14)
.7 V)

For the remaining terms, replace o0D, with 0o, and 0:D where appropriate.

/DD DD, dV
rH1-CB2,vC21,vw x2 21 12 21 r -CB2,z dV

12 L21 12:21

rH1-CB2,sa21,st a,ut 2U a,uv ;'21,wva23,w8 X (),xy (' TH2-CB3,z dV +
12 21 12 21 r CB3 dV +
TH-lCB2,s21,st 2 autl 221 1 ., .23,wxT ( y 2- f dop,z dV +

12 21 12 -21 3 dV
rHl-CB2,rt21,rs ats atuO21,vua23,vwT (/) ,w (' qf, N3 dV +

12 -21 12-21 f

12:21 12:21 f
rB2-H2,va21,vw axw a; 1, .-,_ (YsrM _H2,z dV +

rB2-H2,a21,st 21 12 21 T f, dV +
B2-H2,v 21,w 2 1 12 21 yB2-H2,z
FB2-H2,sCt21,st l2 ut 1a i lur 21wQS,wvT2 (33) 0' IrH2-CB3,z dV +

rB2-H2,sa21,st2 a,ut ( 21 .i ., 23,w (),xy op,z dV +
12-:21 WH12:2,v dV3
rB2-H2,ra21,rs al2 ,tsl2 atuaO21,vua23,vw (/),wx ,z N3,y dV +

12-21 12:21
rH2-CB3,sa32,stT ((j))ut a 21,vw 1, 1 :yrHl-CB2,j z \ dV +
H2-CB3,sT32,st ( 21,w l2 21 2 21, :yB2-H2,z dV +

SrT2-C ,, ,T(23 ,r)g 1s23,srQ 21,st2 nut uau~n (21,wt23,wa ,y H2-CB3,z dV

TH2-C ) a 21,st lt12 1 12 21 T. .. 23,ws ( op dV +
J7 V

rH2-C, T (0),, a,2rqa3,r21 ts2 ltu 21,uavT23,vw 2,z 1 2 J N3,yz d

Oq32,stT f t 23,vut 21,vw aw "12 21 r :yTH-CB2,z uop,s dV +

a32,stT (3)t ,. '21i,vw2 2 1~~ l iyr-B2-H2,z op,s dV +
12-21 12 21 T io
' ,rg 23,sr21,stl a,utl La,uv( 21,wvt ;23,wx y itH2-CB3,z Hopp dV +

S-, (),r a23,sra21,stl aut auva21,wva23,w X (3),xy Uop,pUop,z dV +

a3,p 12:21 12 .;1 .r (I3) o N 3,, r dV +

qf,. T (),ut a23,vua21,vwl 21 12 21 H

q T ,ut o (1212 2 21, B-H2,z N3,sr dV +
q' o T (/3 ),a23,ra21,tvw 1axw wv23,w ,y o N3,py dH2,z

qg ,0 ),Tr(3) 23,sr a2,stl2 2 21 .,wva23,wx 0,y o N3,pnop,z dV +

3 1221 1221 q 3 N
q (,32 p T qp ( 23) ,rp '1, s 2 I, rS a ts 2 21 ..IT ( ),w z i3,z N 3,y d
Lo' ~I r

dV +



W 15)

for the outboard segment the inboard body is the inboard segment). Thus written in the

standard optimal control formulation, the problem becomes

min J = IXcB XD2 + YCB1 YD2 + ZCB1 ZD 2

f(x, u)
ai1 < 900
subject to I/1 < 900

la2 < 900

12a < 900

The above optimal control problem is to be performed in some finite, but unspecified,

time. There are many different problems that fit into the category of an optimal control

problem. Examples include fixed final-time problems, fixed final-state problems and

minimum-time problems.[68] The objective function determines the goal of the control and

many different objective functions can be used for the same system. Chapter 6 presents a

numerical example of an articulated MAV. Included in this example is the formation of an

optimal control problem.


1.1 Motivation

In recent years there has been a surge of interest in the development of Micro Air

Vehicles(\ 1AVs). These MAVs are typically classified as flying vehicles that have a

maximal dimension that ranges from 6in to 24in. There are many different configurations

that qualify as MAVs. They vary based on propulsion type, wing configuration, and flight

regime. There are equally as many different mission profiles as there are configurations

of MAVs. Some of the proposed uses of MAVs are urban observation, chemical sensing,

in-the-field intelligence ,ill i._, :. and small p loadd delivery.

One of the challenges in utilizing the MAV technology is the proper modeling of

vehicles. The challenge in modeling of MAVs is exemplified by the small sizes of the

vehicles and the flight regime in which they operate. Where conventional aircraft have

been studied in great detail and a wealth of knowledge exists on their flight regime, there

is much less information for MAVs. The interaction of the MAV control surfaces and

the low Reynolds' flow regime in which they operate is of great interest. Construction

materials of conventional aircraft are not viable for MAVs. The small size dictates the

use of lighter material but the extreme effect that the flight regime has on the vehicle

calls for building material to also be more functional than the materials in conventional

aircraft. As a result, composite structures are common place in MAV construction. This

use of composite material further stresses the need for new modeling techniques when

investigating MAV dynamics.

However, an exciting new opportunity is provided by the MAV technologies. The

sensitivity of the MAVs to their environments allows for a greater benefit from different

control strategies than conventional aircraft would allow. For example, while morphing

exists on conventional aircraft, the extent of the effect of morphing often pales in

comparison to the complexity of the morphing system. On MAVs, however, morphing

/fvHTB,dV f T B BdV
[M] d p [B] [B dV p v --2 (7-112)
jv fv HTBdV fvB T 2dV

Looking at a representative term of Equation 7-112 and using the indicial notation,


0 1 01 f
TBI-H1,x 1,yx 1l,yzrBl-Hl,z dV +

rB1-H1,w 0 ,xw 1 T (2) ,yz rH1-CB2,z dV +
TB1-Hw ol ,Xw

rB1- HI,w (-dl, W I U (
B1-Hl,w ,XW 1,1y ( 0,y, Uz dV +

rH1-CB2,wT (2),w 01, 0 1,yz B1-H,z dV +

rH1-CB2,vT (2), 0W w (-1 2, 0-) rHi-CB2, dV +
rHI-CB2,vT (2), 0W w 1,yT (2, yz) Uop u dV +

2 W, O ,y 12 2 9 H1-B1-Hlz oopw dV +

T (0 2) 0;J 1 T (02, yz) rHI-CB2,z j u^ dV +
T ( 02),W U 1 0, ,
0T 1 1 yT (12, yZ) U opUoz dV (7-113)

It is important to note that the integrals left in Equation 7-113 are constants integrals.

The derivation of the contribution to the system's kinetic energy from body 3 begins

with the formation of a position vector from the origin of the inertial system to a random

point on body 3. This vector is stated as:

rp RCB1 + B1-H1 + HI-CB2 + B2-H2 + H2-CB3 + op


where rB2-H2 is the vector from the center of mass of the second body to the revolute

joint between body 2 and body 3, written in terms of the body-fixed coordinate system

associated with body 2, rH2-CB3 is the vector from the revolute joint between body 2 and

body 3 to the center of mass of body 3, written in terms of the body-fixed coordinate

frame associated with body 3, and up if the vector from the center of mass of body 3 to a

where [4] is the system constraint matrix and in general is a function of the generalized
coordinates and time. Taking the total derivative of Equation 7-43 produces:

[0] q + [t]t 0 (7-44)

where [+]i is the partial derivative of the constraint matrix with respect to the generalized
coordinates, also known as the system Jacobian matrix, and [l]t is the partial derivative
of the constraint matrix with respect to time. For autonomous systems, Equation 7-44
reduces to:

[ q],q 0 (7-45)

Taking the second total derivative of Equation 7-43 or the first total derivative of
Equation 7-44 produces,

[1q = [] 2 [qt ([q q (7-46)

which for autonomous systems reduces to

[q- _- ([:I^_q ") ]q (7 47)

Any system that simultaneously satisfies Equation 7-44 and Equation 7-46 also
satisfies Equation 7-43. Thus, using the accelerations of the generalized coordinates,
Equation 7-46 can be augmented to the system equations of motion to form the matrix
[M] [ ] J Q + QT [K] q(74)

[K]4 [0] q A Q,
where Q is used as the collection of the terms on the right-hand side of Equation 7-46.
The constraint matrix for the three-bar manipulator is derived in the following. The
revolute joint between the ground and body 1 restricts the position of the end of body 1 to
be coincident with the point on the ground that corresponds to the revolute joint location.

the velocity of the turbulence in the atmosphere is comparable to flight speeds of the

vehicles.[12] Thus there can be large variations in the airspeed over the wings. In fact,

wind speed can change the flight Reynold's number by more than S:l', [13]

MAV wings typically have a low aspect ratio. It has been found that in low aspect

ratio wings, tip vortices dominate the flow.[78]1 The tip vortices have two 1n! i i effects on

the lift characteristics of the wing. Tip vortices on low aspect ratio wings lower the lift

force by reducing the effective angle of attack and provide addition lift force by creating a

low-pressure zone on the upper surface of the wing.[8] Wings with aspect ratios less than 2

are difficult to stall, but are susceptible to roll instabilities. These flight characteristics are

typical of thin flexible wing MAVs.

The aerodynamics of MAVs in low Reynolds' number regime differs greatly from the

aerodynamics of mini vehicles.[15] Ramamurti et al. performed numerical simulations on

the N i.v Research Laboratory's Micro Tactical Expendable (\MITE) MAV.[151 The MITE

is a fixed wing MAV with a wingspan of 6in. and a span to wing chord ratio of 1.25. They

used a finite element based incompressible flow solver and used an Arbitrary Lagrangian

Eulerian (ALE) formulation of the governing equations of motion. Propeller effects

were also included through the use of an actuator disk model. However, this study only

investigated the flow around the MITE with no consideration to structural deformations or

aeroelastic effects.

There are several different classes of MAVs. In a paper by Mueller and Delaurier,

several different MAV configurations are presented.[91 It is evident from these MAVs that

there is a wide range of acceptable designs that result in acceptable mission performance.

These designs even include flapping flight ornithopters that try to take advantage of the

increased Reynolds' number without increased forward flight speed that is associated

with flapping flight in birds and insects. Flapping flight is being investigated as a

inspirational point for flapping wing MAVs. In a paper by Raney and Slominski the

flapping flight of hummingbirds is studied.[16] The focus of the study is the control over


The equations of motion that were generated in C'! ipter 3 are also applicable to other

multi-body systems. One of the most common examples of a multi-body system is the

three bar manipulator. This system consists of three bars that are connected via simple,

single degree of freedom linkages. References [90] and [91] have sections in which this

common example is studied.

7.1 Model Description

As illustrated in Figure 7-3, the first bar is connected to an inertial frame via a

revolute joint that allows the bar to rotate about the axis aligned with the bar. The

second bar is connected to the first bar by a revolute joint that allows the bar to pivot

about an axis that is perpendicular to the axis aligned with the bar and perpendicular to

the plane illustrated in the figure. Likewise, the third bar is connected to the second bar

by a revolute joint that allows the bar to pivot about an axis that is perpendicular to the

axis aligned with the bar and perpendicular to the plane illustrated in Figure 7-3. At the

initial time all of the body-fixed coordinate frames are aligned.

In this example some assumptions are made to facilitate the executions. Primary

among these assumptions is the modeling each bar as a rigid body. This assumption

eliminates the need for a finite element representation of the bars and thus eliminates the

global system stiffness matrix and the terms in the application of Lagrange's equations in

the derivation of the equations of motion that result from the flexible degrees of freedom.

The revolute joints connecting the bars are assumed to be ideal. Thus, they have no

mass and only allow for a single degree of freedom between the bodies they connect. Also

included in the assumption of ideal joints is the lack of friction in the joints. It is further

assumed that the bars are all of one meter in length and are also idealized to have a

unitary density and a zero radius. Finally, it is assumed that no external forces act on the


Figure 1-1: Morphing MAV.

derivation of the kinetic energy of the base body(the fuselage) in this minimal formulation.
The kinetic energy is given by:

T p(v, v,) dV = pvPTvdV (3-12)

where the integral is taken over the volume of the body and vp is the velocity of a random
point on the body. To find the velocity of a point on the body, the first time derivative,
taken with respect to the inertial frame, of the position vector of that point is calculated.
A random point on the fuselage is written as:

rp RCB l+ (3-13)

where RCBI is the position vector from the origin of the inertial frame, E frame, to the
origin of the B1 frame located at the center of mass of the fuselage. From this point on,

RCBI is written in the E frame and u, is written in terms of the B1 basis. Taking the
first time derivative of Eq. 3-13 yields

Ed Ed Ed p B d (U + EwBi x (314)
dt (P) dt (RCB) + dt(mr) RCBI + t-}

Since the fuselage is rigid the point uop does not change in the B1 basis. Thus,

P RCBI + [01] U (3-15)

where [0c1] is a skew symmetric matrix composed of the components of EwB The skew
symmetric matrix can be decomposed as:

[Ol] 4 [ol] + [o4] + 4 [O;] (3 16)


S B RCB_ + [ t] o+ [0w] + 0[w ] (3-17)
Ep C 10'6] UP+ 0')]--0 :)'l U0

segment and the fuselage and the angle between the outboard wing segment and the

inboard wing segment. Mathematically, the optimal control problem for this example is:
min J J E fdt

ail < 900

3i11 < 900
subject to ca2 < 900

I32 < 900

a[[(tf)l = (tf), = a (tf)2 (tf)2 = 00

ZCB,f ZCBI,o = -5 m

where Ci are the control torques applied at the joints between the inboard wing and the

fuselage and between the outboard wing and the inboard wing. Also, the system dynamics

are the ones illustrated in Equation 5-14and Equation 5-15.

In the determination of the solution of the optimal control problem, there are many

different configurations that would be selected. Such a large design space increases the

computational effort needed to solve the problem. To alleviate some of this computational

burden, this example is further constrained to follow symmetric deflections only. This

constraint removes two degrees of freedom and two control inputs, associated with the

deflection of the opposite wing.

A direct method was used in the solution of this optimal control problem. This

method used a Gauss Pseudospectral method in the solution. This method uses twenty

Legendre-Gauss points in the discretization represent the solution of the optimal control

problem. The final time for the system was free to be determined by the optimization

procedure that minimized control effort. This optimal control problem was implemented

in the MATLAB program used in the solution of standard optimal control problems.[92]

As with all direct methods, the cost functional is transcribed into cost function that

* Finally, once a control strategy is formulated the system is ready for simulation with
the implementation of the control. The main criterion of concern in this phase is the
performance of the control algorithm. The response of the system will be compared to
the expected response to determine the fitness of the controller.

1.4 Research Plan

The first step in the solution of the problem statement is the derivation of the

equations of motion for a morphing Micro Air Vehicle. To this extent, Lagrange's

equations will be used. There are several considerations that must be accounted for

first. The overall physical structure of the morphing micro air vehicle is complicated as

can be seen in figure 1-1. As opposed to more traditional aircraft, the MAV presented in

figure 1-1 has wings that are articulated to allow for different angles between the inboard

wing and the fuselage and also between the outboard wing and the inboard wing. Rather

than attempt to model the complete vehicle in a single formulation, a multibody approach,

as demonstrated by Shabana[11, will be used. This multibody approach allows for each

segment of the wing to be treated as a body. Thus, when modeling the system as a whole,

the individual bodies can be composed to represent the system. If the system is modeled

as a single body, each configuration would have to be modeled separately. In addition, the

simulation of the system would then need to connect the separate models based on the

current configuration in some cohesive manner, which is not intuitively a simple procedure.

Aside from the multibody approach to modeling, there is another matter that adds

to the complexity of the formulation of the equations of motion. The vehicle will be

modeled with both rigid body and flexible degrees of freedom. For MAVs, the fuselage can

be modeled as a rigid body as there is not much deformation of the fuselage. The wing

structure, however is highly flexible. The aeroelastic loads on the wings coupled with their

light weight and low stiffness, in comparison to the fuselage,can lead to large deformations.

Thus, as first proposed by Meirovitch[2--5], a hybrid formulation is necessary. The rigid

body motions of the vehicle will result in ordinary differential equations while the flexible

motions typically result in partial differential equations[6]. The ordinary differential

J DY [a12] [T (a)] [a2i]T [N2] dV
D [al 0 21 2 V dV =
rH1-CB2,ji(' 2I j mkn21m N2p dV +
12W21 j J dV +
('-1 -1 aTk (21,m Uop,jN2,np dV +

q .., 12 nck221ln N2,jiN2,p dV (A 13)

S[N2 [N2] dV N2,nmN2p dV (A-14)

It is important to note that all integrals left when written in the indicial notation are
constant integrals and only need be performed once.

where b3 is the unit vector of the body-fixed frame associated with body 3 aligned with
the joint axis of the revolute joint between body 1 and body 2. The derivation of the
equations of motion using the alternative formulation starts with the application of
Lagrange's equations. Looking at the kinetic energy contribution of body 1 to the system,
the position vector to a point on body 1 must be written as:

S= RCBI + op (7-93)

where RCBI is written in terms of the inertial coordinate system and u is written
in terms of the body-fixed coordinate frame associated with body 1. Taking the first
derivative of Equation 7-93, with respect to the inertial coordinate frame, yields:

Ed E (LRC)) +B E B' x (7-94)

but the revolute joint between the ground and body 1 prevents the rigid-body translation
of body 1, thus RCBI = 0. Therefore,

v [o01] tu (7-95)

where [0'1] is the skew symmetric matrix formed from EUB1. However,

[OL1] 0, [ (7-96)


v = 01 [0)] Ut (7-97)

Thus the contribution to the kinetic energy from body 1 can be written as:

= ~2P1 p V dV
t 9P T [ ]T' [04] u dV # (798)
2 Jv OP 6

description of the flexible degrees of freedom, the system equations of motion are second

order ordinary differential equations. Classic control would require these equations to

be cast as a system of first order ordinary differential equations. This system of first

order differential equations is a function of the states and the control. Here, the states

are defined as the generalized coordinates used to describe the motion of the vehicle and

the control are a series of actuation torques applied at the hinge joints of the vehicle that

control the configuration. Once the equations are cast in a form amenable to a control

problem, it is important to differentiate the form of the control problem to be solved.

This task is accomplished by the formulation of a cost functional which quantifies the

effectiveness of the control. In the work presented here, the purpose of the control is to

drive the vehicle to a final state in a finite amount of time. For the purposes of this study,

the cost functional must reflect the distance between the final position of the vehicle and

the desired final position. Here no distinction is made about the orientation of the MAV

when it obtains this final position. The final state is chosen based simply on the instant

the MAV reaches a desired altitude.

In the following chapters, the course of this research will be further explored. Chapter

2 presents a brief look at the major topics involved in this research and some of the

work that has been done in those areas. C'!i lpter 3 presents a detailed derivation of the

equations of motion. C'! Ilpter 4 explains the representation of the forces used in the

equations of motion. C'! Ilpter 5 then casts the problem as a controls problem along with

a detailed discussion of the assumptions, goals, and limitations associated with doing so.

C'! Ilpter 6 then applies the derivation of the equations of motion and the casting of said

equations as a control problem to an example based on Gull-wing MAV illustrated in

Fig. 1-1. C'! Ilpter 7 presents another example based on a three bar manipulator. Finally,

('!C lpter 8 explores some conclusions of the work.

addressed in the same manner as the preceding bodies. This fact allows for a simpler

application of the equations of motion. In the alternative formulation, each new body to

the system must be referenced to the "central" body and the contribution to the kinetic

energy of the system must be calculated anew. The various advantages associated with the

alternative formulation of the equations of motion make it a likely candidate for forming

an optimum control problem.

7.4 Optimum Control Problem

Work by Bottasso and Croce formulated an optimum control problem that sought

to minimize the control effort needed to effect the bars into a final state within one

second.[90] This work will also attempt to minimize the control effort needed to reach the

same final state in one second. To this end, a cost functional must be stated. In this work

the objective function is chosen as:

tf 3
IicCdt (7-131)
i 1

The equations of motion must now be put into a form amenable to the optimal control

problem. Thus the equations of motion become:

= F(X,u) (7-132)


X = qs (7-133)

F (X, u) }I (7-134)
[ Tv + .U J


2.1 Low Reynolds' Number Flow

Low Reynolds' number flow is typically regarded as flow that has a Reynolds' number

of 106 104. In low Reynolds' number flow, flow separation and low lift-to-drag ratios are

common phenomenon.[7',8 Flow separation, transition from laminar flow to turbulent flow,

and reattachment of the flow to the wing can all occur in short distance. The distance

that this occurs in had a large effect on the performance of the lifting surface.

Below a chord Reynolds' number of 50,000, the free shear 1i-r after separation does

not transition to turbulent flow soon enough to reattach to the airfoil surface.[8,91 Then,

when the separation point reaches the leading edge of the airfoil surface, the airfoil stalls.

At Reynolds' numbers greater than 50,000, the separated flow transitions to turbulent

flow and if the adverse pressure gradient is not too severe, the flow reattaches. The area

of the airfoil between the separation point and the reattachment point is typically referred

to as the separation bubble. The separation bubble can occupy between 15'. and 40'. of

the airfoil surface.[9,101 Work done by Lutz et. al has attempted to better qualify wing

performance in low Reynolds' number flow in the development of a predictive tool.[11]

After the flow separates from the airfoil surface, there is an increase in drag and a

decrease in lift resulting in decreased airfoil performance. In the separation bubble, there

is a region of uniform static pressure downstream of the separation point that ends with

an abrupt rise in pressure near the re-attachment point.[10]

2.2 Micro Air Vehicles

The Defense Advanced Research Project's Agency (DARPA) defines MAVs as

vehicles that operate in the Reynolds' number range of 10,000-100,000.[12'13] At Reynolds'

numbers less than 200,000 the lift-to-drag ratio for MAVs is significantly lower than

that of conventional aircraft.[141 Vehicles that operate in this range typically have small

mass moments and inertia that make stability and control difficult. Also, at this scale,

instabilities. There are many quantities in the equations of motion that are dependent

on the current state of the system. As a result the functional evaluation of the equations

when written in first order form, as illustrated by Equation 7-132, is computationally

intensive. The optimization software attempts to use a global polynomial to convert

the continuous equations of motion into a discrete set of equations. At each point in

the discretization, the states and the control are being solved. Thus, the number of

discretization points coupled with the nature of the evaluation of the equations of motion

numerical difficulties are possible. The optimization software used in this work was not

able to overcome these numerical instabilities and was unable to converge to the optimum


where vp is the velocity of a point on the body. This velocity is obtained by taking the
first time derivative, with respect to the inertial frame, of the position vector of a point on
the body. The position vector to a point on body 1 is written as:

Lp RcB + (7-12)

where RCBI is the position vector from the origin of the inertial frame to the center of
mass of body 1, written in terms of the inertial frame, and u, is the position vector from
the center of mass of body 1 to the point, written in terms of the basis vectors of body
1. Taking the first time derivative of Equation 7-12 with respect to the inertial frame

Ed Ed E dU E Bd
t (_p) Ed (RCB1) l (OP) + x op
= RCBI + [ ] IP (7-13)

where [0Z1] is a skew symmetric matrix formed from the components of the angular
velocity of body 1 with respect to the inertial frame. However,

[0w"] [01] uo + 01 [o0,J] + 41 [o~ ] (7 14)


p RCBi + 1i [o01] rp + 01 [o0L)w] ]U, + i [0 ] (7-15)

Re-writing Equation 7-15 so all components are written in terms of the inertial frame

p, = RCB + 1 [Ci] [,1 u] P 01 [CI] [00 )L] U + [CI] [0

where b2 is a unit basis vector of the body-fixed coordinate frame of body 2 that is

perpendicular to the body-fixed unit basis vector in body 2 that is aligned with the axis

of rotation of the revolute joint between body 1 and body 2. Writing all components of

Equation 7-64 in terms of the inertial coordinate system produces:

6 b IT[Cl]T[2]b b 0 (7-65)

Taking the partial derivative of Equation 7-65 with respect to the generalized coordinates

06 oT r c_ [o T U
9 1 [,]T [_lT [C2 b3 b1T ) ]T [Cu] ['02]Jb2
L O ] -1 03
bT [0o:, [CI ] [7 2 b 2 bT 1T [ \ 1 1 ] [C 2] [ 1, 2
6iT o ^ [Ci] [C2] b 0 T [Ci] [C2 1
T [C] [C2] T [: 10T
,b)T [271T [I2,1 L1w] b2 bIT [CI]T [671 [1w2 ] [ b2 t

0 0 0

The constraints on the relative translational motion between body 2 and body 3 caused by

the third revolute joint, are illustrated in

47 RCB3 + H2 -CB2 -H2 0 (7-67)

where RCB3 is the vector from the origin of the inertial system to the center of mass

of body 3 written in terms of the inertial system basis vectors, u 2 is the vector from

the origin of the body-fixed coordinate system in body 3 to the position of the revolute

joint between body 2 and body 3 written in terms of the basis vectors of the body-fixed

coordinate system of body 3, and uH2 is the vector from the origin of the body-fixed

coordinate system of body 2 to the position of the revolute joint between body 2 and body

3 written in terms of the basis vectors of the body-fixed coordinate system of body 2.

Writing all components of Equation 7-67 in terms of the inertial coordinate system yields:

7 = RCB3 + [C3] 32 CB2 [C2 22 = 0 (7-68)

Re-write Eq. 3-17 so that all terms are in the B1 basis:

Ev, = [C1]T RCBI + [0 ] uP + [;] + [o4] uOP

[B] [- [ ] p U

[OL:Dl] [0:l] o ]


VP= [Ci]T [B] R a


0 0

Thus the kinetic energy for the fuselage is

-T [ [Ci] [B] R_ }dV
S[B] L J eT \

L p [I] dV
f p [f]T [CI]T dV

where [I] is the identity matrix and let

Then the kinetic energy for the fuselage is

T = t ] [M]


where [M] is the mass matrix for the fuselage. Next examine the integrals to be performed
in the computation of the mass matrix for the fuselage. To facilitate this examination the






f p [C,] [B] dV

f p [B]T [B] dV



(3 24)


q -

This section describes the formation of the mass matrix for body 2. The use of
indicial notation will facilitate the derivation of terms to be integrated. The mass matrix
for body 2 is:
[I] [C ] [B] [Ci]D, [C1 [a C [T(a)] [a21T [N2
] = [B]T [1]T [B]T [B] [B]T D [B]T [a12] [T ()] [ ]T [N21 ]
2 DT [01]T DT[B] DTDI DT [-12] [T () [21T [N2
[N2]T [21 [T (o)]T [a12]T [C ]T [N ]T [a21 [T ()]T [l2 ]T [B] [N2]T [21 [T ()]T [al2]T [N2T [N2
Examining the integrals term by term and using the indicial notation:

S/lp[I]dV= 2pI,mdV (A2)

[CI] [B] [Ci] B2 B3 (A-3)

Thus examining a representative term,

j [CI] B dV

C1,mn ,nprB1-H1,p / dV +

0 T ( 1), a2,nmq p j 2,n dV (A4)

where fJ dV, J uop,p dV, and fV N2,Tp dV are invariants. For the terms fV [CI] B2 dV
and fv [Ci] B3 dV replace 01ow with O1, and o; respectively.

/j [C1[] 2] [T (a)] [21]T [N2] dV Cl, ._..T (a),km a21,m j N2,Tp dV (A-5)

fv B BjdV fv B B2dV fv B B3dV
[ I 2 V -3
J [B] [B] dV = f BTBjdV fv BB2dV f, BB3dV (A 6)


Sean James Sidney Regisford was born in St. Georges, Grenada, on July 25, 1977. In

1979 he and his family moved to Florida. He graduated from Miramar High School in 1995

and then attended the University of Florida, studying aerospace engineering. In 2000, he

graduated from the University of Florida with a bachelor's degree in aerospace engineering

and then was accepted to the graduate program at the University of Florida. In 2002,

he received his Master of Science degree in aerospace engineering and then continued his

studies in the Ph.D. program. Upon completion of the degree Doctor of Philosophy, he

hopes to obtain employment in industry

problem. The use of the Calculus of Variation in concert with Pontryagin's Minimum

Principle states that the optimal control is the control that minimizes the Hamiltonian.

In the indirect method, the solution of Equations. 5-5-5-9 is required. The difficulty in

solving those equations lies in getting a closed-form solution of the costate equations.

The advantages of indirect methods is that the solution of the costate equations serves

as a good indicator of the optimality of the solution. Numerically, indirect methods are

typically solved using shooting or multiple shooting methods. These shooting methods

have much in common with standard root finding techniques.

Direct methods, also known as direct transcription, solve the optimal control problems

by converting the optimal control problem into a nonlinear programming problem. The

transcription process converts the infinite-dimensional optimal control problem into a

finite-dimensional optimization problem. The dynamics of the problem are discretized

at a finite set of points and the cost functional is replaced with a cost function that

can be evaluated at those discrete points. This finite-dimensional problem is known as

a Non Linear Program(NLP). Direct methods have the advantage of not requiring the

costate equations to be explicitly derived. In addition, the direct method allows for faster

convergence of the solution. As a result, direct methods are used quite often in practice, as

illustrated in work by Bass, von Stryk, Hu, Frazzoli, and Hull.[67'70--731 A disadvantage

of using direct methods is that the optimality of the solution is not guaranteed as the

costate equations are not solved. Direct methods have a tendency to get stuck in local

minima. Once converted to a non-linear programming problem, direct methods for the

solution of the optimal control problem have well-known methods. The solution often

depends on the discretization chosen. Standard choices include the Lobatto methods that

include both endpoints in the discretization, Radau methods that include only a single

endpoint in the discretization, and Gauss methods that do not include either endpoints

in the discretization.[74'751 An example of the use of Lobatto discretization is provided in

studies by Williams, Fahroo, and Herman.[76--78]

Taking the partial derivative of Equation 7-53 with respect to the generalized coordinates


0q= O eT [Ci] [0] b 4 [CI] ] [ ]bl 4 i, o [ ,] bl O 0 0 0 O 0 0 0
q (7 54)


S3 = e b' = 0


Or presenting all quantities in terms of the inertial frame:

3 -= [CI] b'= 0


Taking the partial derivative of Equation 7-56 with respect to the generalized coordinates


aT= O e [Ci] [0)1] b T [C,] [0 )1] bl JT [Ci] [0"1] bi O 0 0 0 OT 0 0 0

For the revolute joint between body 1 and body 2 the constraint on the relative

translational motion between the bodies is expressed as:

14 -RCB2 + m1 CBI U 1


where RCB2 is the position vector from the origin of the inertial frame to the center of

mass of body 2 written in terms of the basis vectors of the inertial frame, uH is the

position vector from the center of mass of body 2 to the revolute joint between body 1 and

body 2 written in terms of the basis vectors of the body-fixed coordinate system attached

to body 2, and u~, is the position vector from the center of mass of body 1 to the revolute

joint between body 1 and body 2 written in terms of the basis vectors of the body-fixed

coordinate system attached to body 1. Expressing all components in terms of the inertial

frame yields:

(4 RCB2 + [C2]1 2I- RCB1 [C1] U1 = 0


For the remaining terms replace o0D, with O^1, and 0:D, where appropriate.

S[BB] DDidV j BB DD dV-
< V v 2DD

Examining a representative term:


rB1-H1,ul -,vual2,vT (a),w 12;?1 a21,zyrH1CB2, V
rB1 H1,uOc-v1,vua12,vt (a),., 12 -21
-Bl-HIu0 wal2, ,wa ,( l21,zyrB2H2,z dV +

Bl- H1,r c,sroal2,stT (o),tu 2 1,uva21,wuv23,wxaT (3),xy a32,zyrH2 CB3,z dV+

rB1-H1,r Wo,srZa12,stT (,12a) tu WLu 21,uwv23,tuT (,3)xy a32,zy uop,z dV +

rBl-Hl,q rqa12,rT (a),st 121tua21 ,vua23,wT (3) ,a z f N3,yz dV +
rH1-CB2,ra21,rsT (a),ts l2,uto,0w ,o l2,uT (a),wx 12 -?21 ,
0 a21,zyrH1-CB2,z dV d

rH1-CB2,ra21,rsT (a),t ail2,ut O l,vual2,vTuT (a),wu 1 ,21,zyB2-H2,z
0 a21,zyrB23H2,z dV af

rH1-CB2,na21,npT (a)qp al2,rq0 ,sr a12,stT (a) ,tu 121,az21,v23,tu T (/3) ,y a32,zyrH2-C3,z dV +
rH1-CB2,n21,npT(a) pal2,rqO ,sr212,stT (a),tu 12 %v 21,wt23,T (/3),y 32,zy op,z

rH1lCB2,ma21,nT (,pn a) 12,qpw O ,rqal2,rsT(a),st 12W1,tua21,vua23,vT (P),a V N3,yz dV +
rB2 H2,r,21,rsT (a),ts a2,utoW ,vual2,vwT (a),w 12 21,yH1CB2, V

B2-H2,rO21,rT (a) ,ts al2,ut o ,vua12,vT (a),wxz 1,zyB2-2,z d
0 a 21,zyrB2-H2,z dV d

rB2-H2,na21,npT (a),qp a 2,rqO ,,s 12,stT (),tu 12w1 21,u23,xT (3), 32,zH2-CB3,z dV +

rB2-H2,n21,npT (a),qp a2,rq w,sral2,stT(a),tu 122,va21,wa23,wxT (3),xy 32,zy op,z dV +

rB2-H2,ma21,mnT (a)pn al2,qp ,rqal2,rsT (a),st 1221 ta21,vua23,vwT (/), N3, d

rH2-CB3,na32,npT (3),qp a23,rqa21,rsT (a),t, a12,ut ,vua12,vwT (a),wx 12 ,21, yfH1-CB2,z d +

rH2-CB3,na32,npT (),qp a23,rqa21,rsT (a),t al2,ut0 o,val2,wT (a),wx 12 21,zyB2-H2,z V +

rH2-CB3,ja32,jkT (),mk 23,nm 21,npT (a),qp 2,r ,sl2,stT (a),t 12w ,u21,w 23,aT (3)y 32,z

rH2-CB3,ja32,jkT (3),mik 23,nm 21,npT (a) qp a2,rq c,srt12,stT (a),tu 121 ,uva21,wvOu23,wT (/3),xy a32,z

rH2-CB3,ia32,ijT (),kj a23,mka21,mnT (a) ,pn a12,qpO w~q,rq2,sT (a) ,t 12 lwtut21,vul 23,wT (3),x

o32,npT (3)p 23,rqO21,rsT (a),ts 12,u ,vu12,U T (a),wx 12 a21,zyrH1-CB2,z op,n +

a32,npT (3)qp a23,rqa21,rsT (a),ts l12,ut 0 cvu12,uT (a) ,x 12 ?1 ,t21,zyrB2-H2,z op,n dV


yrH2-CB3,z ; dV

y op,z dV +


This section details the formation of the mass matrix for body 3. The mass matrix for
body 3 is as follows:
[I] [CI] [BB] [CI] DD1 [C1 DD2 [C1 [12] [T ()] [a21]T [a23 [T (3)] [a32T [N3
[BB]T [BB] [BB]T DD, [BB]T DD2 [BB]T [012] [T (a)] [a21]T [-23] [T (3)] [.32T [N3
[M] = J P DD DD DDTDD2 DD [1a2] [T (a)] [a21] [23] [T (3)] [32 [3 dV (B-1)
sym DD2DD2 DDn [o12] [T (a)] [o21] [23] [T (3)] [32 [N3
N3]T [N3
Indicial notation will be used to examine the details of the integrals involved in the
preceding equation. Examining the first component:

j [I] dV = I,,,p dV (B-2)

where [I] is the identity matrix.

/ [Ci] [BB] dV -= [Ci] BB, BB2 B 3 dV (B-3)

Examining a representative term,

J [C11] BB dV =
C1,xy ryzBl-H,z dV +

Cl,uv Lo l ... .T(a),-_, ,irHI-CB2,z dV +
C1,uv1aO ,vwO2,wxT (C),xy (,_ ,rB2-H2,z dV +

1,rs o,sta12,tuT (a), 21,wva23,wxT (3)y ,rH2-CB3,z dV +

C1,rs Ostl12,tuT (a),uv o21,,wva23,wxT (3),y Uop,z dV +

Ci,pq1 ... rsT (a),st a21,uta23,T (/)w ,w ( ,z N3,y dV (B-4)

For the remaining terms replace 011 with 0o1i and o1;, where appropriate.

J [C'] DD, dV

Cl,uval2,vwT ({a),w 12Z1, 2 '., :yrHI-CB2,z dV +




Figure 7-8: Time history of the angular rates as determined by integration.


0 0.2 0.4 0.6 0.8 1

Figure 7-8: Time history of the angular rates as determined by integration.

example, the invariants involved in the generation of the mass matrix of the system, as

outlined in Appendix A, are reduced in number to only being fy dV, f Uopp dV, and

v uop,pUop, dV. Similarly, the generalized coordinates are reduced to









For the purposes of this numerical example, the invariant integrals were carried out using

four-point quadrature along with an isoparametric element used for each body. Thus,

/ pdV = 1 det J p t (6-1)

where t is the thickness of the element, and det J is the determinant of the Jacobian.
Using an isoparametric element, the Jacobian becomes

Ji (6-2)
S(- (1 y )) xi (1 )2 (1 iy)x 3 (1 t)x4 4 ( (1 u)) 1 (1 'u 2 (1 )3 (1 W) 4(62)
( (1 -l i)) (1+ i)2 + (1+i) 3 + (1-i) 4 1 (1 -(i)) 1 (1+i) 2 + (l+i) 3 + ( Y4i) 4

where xl and yl are the elements x and y coordinates of the first node of the element

and rli and ,i are the integration points in the transformed space. For the four point

attack. It is known that as the aspect ratio is increased, induced drag decreases.8 1 Thus,

in variable span morphing, the aspect ratio is increased in an attempt to obtain a better

lift-to-drag ratio.

2.4 Modeling

There are several concerns that must be addressed when modeling the behavior

of a Micro Air Vehicle in flight. As demonstrated in the preceding section, there is a

plethora of research into the individual components of the problem, such as research

in low Reynolds' number flow, research into various morphing strategies, and research

into the design and manufacturing of Micro Air Vehicles. Things become more involved

when the interaction of the various disciplines becomes the primary focus of the research.

The interaction of the vehicle and the fluid through which it interacts then becomes the

over-arcing mechanism that is needed to tie the various disciplines together. For a flexible

body, the forces produced by the medium through which the body moves influences

the deformation of the body which then in turn influences the forces generated by the

medium. This interaction between the body and the medium is know as aeroelasticity in

aircraft.[32'331 Aeroelasticity is also known to have an affect on the controls and thermal

effects of the craft, thus its study is of crucial importance in the proper modeling of an

air vehicle.[341 Fluid-structure interaction is not only limited to aerospace interests. In

work by Kral and Kreuzer, the modeling of fluid-structure interaction was addressed as it

impacts offshore platforms.[35]

In order to capture this interaction, the choice of modeling becomes important. There

are several v--i- to describe any given system. Standard rigid-body dynamics generally

treats the system as a single entity and is concerned about the overall motions of the

vehicle. This approach is clearly not acceptable to demonstrate aeroelasticity. On the

other end of the spectrum of modeling techniques are Finite Element techniques that

discretize the continuum of the vehicle in an attempt to capture the flexibility located in

true systems. These techniques, however, are usually more concerned with the deformation

body can be separated into components derived from the individual bodies that compose
the system. The principle of virtual work states

6W = Qi 6q = F, bri (3-77)
i i
Assume there is some force Fp applied at some point pr on body 1. The principle of
virtual work can be used to find the generalized force that results from the force Fp being
applied to the body. Write the position vector from the origin of the inertial frame to the
point pr on body 1 as:

rp r RCB1 + [Ci] ,b (3-78)

where ub, is the position vector from the origin of body 1 to the point pr. Taking the
virtual variation of the position vector yields:
0 ([C1] U,) 0 ([Ci] U,) 0 ([Ci] p) (
6pr, = 6CB + (3-79)

4, [CI] [0:i ]
9([C1]up) [Oil [0w 1] (3 -80)
-() 0:41ou
a ( C i lu o ') [ C_ I [ O i

Thus for a point on body 1:

6r, = [I] [C ] [B] 6 i (3-81)

Re-writing Eq. 3-81 in terms of the B1 frame yields:

6r [C[i] [B] { e C (3-82)
S[J 60

Cl,uvla2,vwT (c),wx 1 21, .1 :yrB2-H2,z dV +

Cl,rsO1'l2,stT ((i),tu 12 u a21, wva23,wxT ( y3),y -

Cl,rst'Il2,stT () ,tu a;Q ,uva21,wva23,wxs (/3),y ('

Cl,qrQa2,rsT (C),st 12atua21,vu23,vwT (3),w

S[Ci] DD2 dV

Cl,rscl2,stT (C),tu a21,vu'23,vwT (),w 23 32

Cl,rs'Q2,stT ('),tu o21, .. ),T L(23)32 ,

Cl,qr1l2,rsT (a'),st 21,uto23,uvT (3),V 23;32 -

S[BB] [BB] dV




,rH2-CB3,z dV +

J Up,z dV +
q j N3,yz dV

TrH2-CB3,z dV +

op,z dV

fz N3,yz dV

1 3


Examining a representative term:

JV1 1
rBI-H1,x021 0-1 Bl-Hl,z dV

rBl-H1,uO ,vu0l ',val2,wxT(a),y a21,zyrH1-CB2,z dV +

01 01 dV
rBl-H 1,u0 ,^vu0 v"vwal2,wxa T(),xy a21,zyrB2-H2,z dV +

rBl-Hl,r ,sro W1,stal2,uT (a),uv 21,wv23,wxT ( y3),xy 32,zy Uop,z dV +
0 1 01 dV
rBl-Hl,qOW,rqo wC,rsal2,stT(a) ,tu21,vua23,vwT (3),a N3,yzdV+

rH1-CB2,ua21,uvT (a),w 2,xw 01 0-1 rBl-Hl,z dV+

rH1-CB2,2,r2,rsT (a),ts al2,utO0 w,vuwOvwal2,w (a),xy a21,zyrH1-CB2,z dV

rHl-CB2,ra21,rsT (a),ts al2,ut ,uOvwal2,wx (a),y 21,zyrB2-H2,z f dV +
rH1-CB2,n~21,npT (),qp al2,rq0 s ,stO 2,tuT (a),uv 2,w v23,wx )T )xy a32,zyrH2-CB3,z dV

rH1-CB2,2,na2npT(a)qp al2,rq o ,sr0o ,st2,tuT (a),v 21,wva23,wxT (3),xy 32,zy op,z dV +








Since the wing structures are modeled as flexible bodies, the generation of aeroelastic

forces becomes important. The aeroelastic forces used in the model will influence the

shape of the wings, due to their flexibility, and the shape of the wings will influence

the aerodynamic forces. There are many different schemes available for the generation

of aeroelastic forces. These schemes range from potential flow theory to full 3-D

Navier-stokes computational fluid dynamics simulations. The generation of aeroelastic

forces is a sufficient research topic in its own right. Work by Ortiz and Barhorst

demonstrates a number of different techniques in the generation of aerodynamic forces

for use in aeroelasticty studies.[571 However, the generation of these forces is not the

focus of the work presented in this document. Thus, it is necessary to find a method to

generate the forces that is robust enough to influence the fluid-structure interaction but

simple enough to be integrated easily into the form of the equations of motion presented in

C'!h pter 3.

4.1 Athena Vortex Lattice

Athena Vortex Lattice (AVL) is software developed by Mark Drela of the Massachusetts

Institute of Technology and Harold Youngren of Aerocraft, Inc. It is an open-source

program that uses lattice vortex method to generate an aerodynamic simulation of a

given aircraft configuration. AVL is best suited for aerodynamic configurations which

consist mainly of thin lifting surfaces at small angles of attack and sideslip, thus making

it suitable for aerodynamic analysis for MAVs. These surfaces and their trailing wakes

are represented as single-i- r vortex sheets that are discretized into horseshoe vortex

filaments. AVL also has the capability to also model slender bodies such as fuselages and

nacelles via source+doublet filaments; however, force and moment predictions on these

surfaces should be used with caution.

The optimization program failing to converge to an optimum is a difficulty that direct

methods are known to exhibit. In particular, the adjustments to the equations of motion

that were necessary to enable execution in GPOCS also renders the problem ill-behaved.

Typically direct methods work best when the equations of motion can be discretized in

time to still form a smooth function. In the case of this example, however, this smoothness

is not guaranteed because of the nature of the generation of the aerodynamic forces.

Even if the equations still used a direct call of AVL in the generation of the forces,

configurations that the optimization uses as trial candidates would violate the validity

of those forces. In turn, these erroneous forces fracture the smoothness of the discretized


19.7 19.9

Figure 6-1: Geometric model of gull-wing MAV.

1.5 Contribution

The major contribution of this work to the pantheon of information on the dynamics

and control of Micro Air Vehicles is the use of the multi-body formulation of the equations

of motion. The multi-body formulation presented in this document allows for the

determination of the contribution of the individual airplane components. Standard

modeling choices are either concerned with the motion of the overall system, or with

the minutia of the flexible bodies without any consideration as to how small scale

deformations effect the total motion of the system as a whole. The use of multi-body

dynamics for the modeling of MAVs will bridge this existing gap. The ability to see

impact that typically small scale deformations have on the overall system meshes well with

the configuration changes that are possible through morphing on MAVs.

The secondary contribution of this work is in the control strategies used on micro

air vehicles. The equations of motion as derived in this research are highly nonlinear

and time-dependent. Standard control strategies would choose an operating point to

linearize about, then create a control strategy based on that point, and finally tune the

performance for conditions that don't match the linearization point. This type of control

strategy, even when using multiple linearization points and switching controllers, does

not achieve the full performance potential for a morphing MAV. Morphing allows for an

infinite number of configurations with associated performance. Thus, a control strategy

that takes all of the configurations into consideration is highly desirable. The equations

of motion as derived in this work can then be cast into an optimal control problem.

The solution of the optimal control problem not only finds a control strategy for the

vehicle, but depending on the cost functional chosen can present a format that allows the

morphing of the vehicle to be used as the only means by which to obtain some desired


indicial notation convention is used. Thus,

f p [I] dV plm dV (3-26)
which is a constant matrix.

Sp [C] [B] dV = fV p [Ci] [0o] udV f p [CI] [0:o] uodV fy p [CI] [o ] updV
Looking at the first term of Eq. 3-27

Sp [CI] [0)] uOdV C Ci,imo, f dV (3-28)

where fy, i dV is constant. For the second and third terms substitute in 0ow and 0:,
where appropriate.

Sp[B T[B] dV= p U IT[oI [[O [ ] [Ow [L] [ ] dV (3-29)
_, P [0 : 1

Looking at the term in the first row and first column of the resulting matrix yields

Sp [O]T LO ] dV = ,mionm jI pUopu dV (3-30)

where fS puop,iuop,TdV is constant. For the remaining terms substitute in 0ow and 01;
where appropriate.
Next create the kinetic energy of the inboard wing. The position vector to an
arbitrary point on the inboard wing is as follows:

p RCB1 + LB1-HI+ H1- CB2 + p+ [N2 q2 (3 31)

where rBIH1 is the vector from the center of mass of the fuselage to the hinge point
between the fuselage and the inboard wing, _HICB2 is the vector from the hinge point
between the fuselage and the inboard wing to the origin of the coordinate system fixed in

To my mother Maria, who has never allowed me to give up on my studies when my

motivation was lacking.

where [T (a)] is the transformation matrix between the joint coordinate system, aligned

with the joint between body 1 and body 2, fixed in body 1 and the joint coordinate

system, aligned with the joint between body 1 and body 2, fixed in body 2. This

transformation matrix is a function of a, which is the angle between the fuselage and

the inboard wing segment. Similarly, [T (/3)] is the transformation matrix between the

joint coordinate system, aligned with the joint between body 2 and body 3, fixed in body

2 and the joint coordinate system, aligned with the joint between body 2 and body 3,

fixed in body 3. This transformation matrix is a function of 3, which is the angle between

the inboard wing segment and the outboard wing segment. The fuselage is considered

rigid and the wing segments are considered flexible bodies. The inboard wing segments

are hinged to the fuselage. Likewise the outboard segments are hinged to the inboard

segments. In this work, these hinges are idealized as having no mass and being a single

point that is coincident in both bodies that the hinge connects. Also, the hinges or joints

only allow a single rotational degree of freedom between the bodies that are connected by

the hinge.

For the minimal set formulation all bodies will be referenced to a single body. In

this case, all bodies will be referenced to the B1 frame associated with the fuselage. The

equations of motion will be derived using Lagrange's Equations. Lagrange's equations are:

d 9T\ 9T 891V 98T
d(O T+ = 7+ Q (3-t11)
dt 6) 6q 6q 6q -

where T is the kinetic energy, V is the potential energy, R7 is the Rayleigh dissipation

function, q are the generalized coordinates, q are the first time derivative of the generalized

coordinates taken with respect to the inertial frame, and Q are the generalized forces

acting on the system. The derivation of the system equations of motion will start with

the derivation of the kinetic energy. The kinetic energy of the system can be written as

the sum of the kinetic energy of the bodies that compose the system, since the kinetic

energy is a scalar function. Thus, the derivation of the kinetic energy will start with the


rHl-CB2,sa21,tl2 ,ut W '1- a23,vwT ( ,w 2332 I rH2CB3, dV +
Hl12CB2,21 21 t T tvu3,vw 23 32 / opz dV +

B2-CB2,sa2l,st ,l2 ut21 ,v ,vw ,w2 2 ,rI H CB3, dV +

TB2-H2,sa21,st l2 ,uta21,v 23,v (/3),w 23 32 dV +
,,2: 3 j N3,, dV
TB2-H2,r'21,rs1 221 23,u ( 2332 vT ,(,)z N- dV+
rH2-CL ,rq 3) a23,sr 21,stl221 ,r ) 23 32, ,H2-CB3,. dV +
12:21 T23
rH2-CL ,q a23,sr21,st5 ut221 u23,vwT ( ,w 23 32 UopI, dV +
rH2-CL 1 T(3) 232, '22, 1, 0 .,'- 23 132 -, 3 N3,, dV +

1221 T(0w 23 32 /pd
,T q () 23,sra 21,st2 a ,ut a21,vu T23,vw ( 2332 I rH2-CB3,z op dV
S ,( 3),r 23,sra 21,st2 1uta 21,vua 23,vw T (jw 23; 32 op, o dV +

a32,npT ( qp a23,r I 12213 IT (3) 2332 ,z J opN3,yz dV +

qn ,T ( 3),rq 23,sra21,st ta2l. i- ..T (3),w 23032, N3,p opz dV +

qfma32npT )qp 23, l221 ., 3) 2332, ,q N3,nmN3,y i 16)

/ DDT [al2] [T (a)] [a21]T [a23] [T (/3)] [a32]T [N3] dV

rH1-CB2,sa21,t ...T (/3),w 3,, dV

rB2-H2,-C 21,t w L,ut a2,1,23,vw (21 ),wx tl2 ( N3,v 3, dV +
H2-CT (q

12-21 T '0),w,
S, ,T (3)q, a23,sa21,stl Lut21,vu23,vw (s3) UoppN3,yz dV +

q T (!),r 23,sr'21,stl2 -I N23, 3pN3,y dV (B-17)

/ DD DD2 dV
23r32 23 -32 B, d
rH2-CB3,v<132,vw W3 w rH2-CB3,z dV +

23L32 23 32 3fN, dV
TH2-CB3,v1'32,vw 3, w L .q pz +

2332 23-32 /i
TH2-CB3,uQa32,uv23 3223 32 N3,uy dV +

'32,vw 3 3,xw r IH2-CB3,z Uop,v dV +
3 23 32 23 32 dV +
a32,vw W3 w J

Q32,uv23,wv 2 32, ,- ,fz Uop,uN3,yz dV +

f ..... /32, 3 TH2-CB3,z 3,vu dV +

3 23-32 23-32 [-
qfua32,vw 3,xw LL) N3,vuUopz dV +

qf,t 32,uv /3,wv .." ,'f ,z 3,utN3,yz dV (B

/ DD [ 2] [T ()] [2 [23] [T ()] [32 ] 3 idV
23-32 f

0H2-CB3,vw 32,vw op 3yz dV +
qf w23u32 2,vw L,, N3,vuN3,yz dV (B-19)


S[N3T [N3] dV N3,yN3,yz dV (B-20)
Note all of the intervals left in the indicial notation are constant.
Note all of the integrals left in the indicial notation are constant.

can be evaluated based on the states and control inputs. This function, along with the

constraints and the equations of motion that must be satisfied, is discretized in time. At

each discretization point, the states and Karush-Kuhn-Tucker multipliers are solved. As

a result, the equations of motion must be simultaneously evaluated at each discretization

point. The equations of motion for the articulated MAV example illustrated in this

chapter, however must be integrated in time, as written. Moreover, the generation of the

generalized forces from the aerodynamic forces requires a execution of AVL based on the

current state of the vehicle. In the integration of the equations of motion, this process is

the largest computational expense. To alleviate some of this expense, the generation of

the aerodynamic forces for a range of values that the vehicle may occupy were determined

and the results tabulated. The data is then available for interpolation to determine the

aerodynamic forces without the use of AVL for each new state the vehicle attains during

the optimization process.

Figure 6-2 shows the global position of the vehicle generated by the optimization

program after it was unable to improve the cost function over the previous iteration's

calculated control vector. Similarly, Figure 6-3 shows the morph angles generated by

the optimization program. The plots produced from the optimization software seems

to indicate the desired decrease in altitude of 5m with the constraint of the morph

angles returning to zero deflection at the time the vehicle completes the decrease in

altitude. However, these results do not reflect the successful solution of the optimal control

problem. This fact is illustrated by Figure 6-4 and Figure 6-5. These figures show the

position of the vehicle and the morph angles after integrating the equations of motion

using the optimal control as determined by the optimization software. As illustrated in the

figures, not only does the vehicle not complete the 5m decrease in altitude, but the morph

angles are also non-zero at the termination of what the optimization programs states as

the maneuver time.

Therefore the velocity of a random point on body 1 is written as:


vPp [C]T [B] D, [a12] [T (a)] [221]T [N2] (3 48)


Now the kinetic energy of body 2 can be written.

T j= tp(v v) dV = pvvPdV (3-49)

Substituting Eq. 3-48 into Eq. 3-49 yields

T = 4 [M]q (3-50)


S< > (3-51)


[I] [CI] [B] [Ci] D [C1] [a12] [T (a)] [a21] [N2]

f 1 [B]T [B] [B]T D [B]T [a12 [T (a)] [a21] [N2(3-52)
J 2 DTD D [a12] [T (a)] [a21] [N2]

,/,/II [N2]T [N2]

Please refer to Appendix A to see the details of the integral involved in the formation of

the mass matrix for body 2.

As was done with body 1 and body 2, the contribution to the system equations of

motion from body 3 will now be derived. Thus the equations of motion for body 3 proceed

from application of Lagrange's equations. This process begins with the derivation of the

random point on body 3, written in terms of the body-fixed coordinate frame associated
with body 3. Taking the first time derivative of Equation 7-115 with respect to the
inertial coordinate system yields:

Sd B + d E X B-H H-CB2 E B H-CB2
dt dt -C d\B1t (rBlHI1 H1+ dt CB2I + Wl
2 d B2 B3 d E EB3
S (rB2-H2) X LB2-H2 + (rH2-CB3) X H2-CB3
d dtC
B3 d (3 x (7-115)

E B2 E B1 +B B2 (7 116)

EB2 E + B1 B2 + B2B3 (7117)

Inserting Equation 7-116 and Equation 7-117 into Equation 7-115 and factoring the skew
symmetric matrices produces:

vp 1 ([ 1 rB-H, + [01 ] [T (2)] rH1CB2 + [0)j1] [T (02)] rB2-H2

[o0 ] [T (02) [T (3)] rH2-CB3 + [O:1] [T (02)] [T (03)] p) +

2 ([T (2)] [12] H1r-CB2 + [T (02)1 [1] rB2-H2+

[T (02)] [1W2] [T (03)] -rH2CB3 + [T (02)] [1)2 [T (03)] p)

<3 ([T (02)] [T (03)] [11 TrH2-CB3 [T (2)] [T (03)] [1] o) (7 118)


BB1 = [ ro H)i +] BH1 [T (2) H1-CB2 + [0)Wj1 [T (2)] r2-H2 +

[O] [T (02)] [T (03)] rH2-CB3 + [0J w [T (02)] [T (03)] p (7-119)

BR2 -[T (02)] [1] rHICB2 + [T (2)] [1] W r B2-H2 +

[T (2)1 [1w] [T (3)] rH2-CB3 + [T (2)1 [1W] [T (03)] (7 -120)
[T(02) 621 21 _O

stated as:

"B2 [T()r3 (88)

1 0 0

[T(3)] = 0 cos 3 -sin 3 (7-89)

0 sin 3 cos 03

and rB3 is a vector written in terms of the body-fixed coordinate system associated with

body 3 and rB2 is the same vector written in terms of the body-fixed coordinate system

associated with body 2. The relationships between the bodies in the system are completed

with the definition of angular velocities between the coordinate systems that comprise the

complete system. The angular velocity between the body-fixed coordinate frame associated

with the first body and the inertial frame is written as:

E 1 13 (7-90)

where b is the unit vector of the body-fixed frame associated with body 1 aligned with

the joint axis of the revolute joint between the ground and body 1. The angular velocity

between the body-fixed coordinate frame associated with the second body and the

body-fixed coordinate frame associated with the first body is written as:

B1 B2 2b (7-91)

where b is the unit vector of the body-fixed frame associated with body 2 aligned

with the joint axis of the revolute joint between body 1 and body 2. The angular

velocity between the body-fixed coordinate frame associated with the third body and

the body-fixed coordinate frame associated with the second body is written as:

B2 B3 3_ (7-92)








0 0.2 0.4 0.6 0.8
time (sec)

Figure 7-1: Norm of the constraint vector








0.4 0.6
time (sec)

Figure 7-2: Norm of the constraint vector

Taking the partial derivative of Equation 7-69 with respect to the generalized coordinates


S[0] 0 0 0 [I] -[2 [1] 2] UH2

[C 2 [1W ] 1uH2 [C2] [1:] )2u [I] [3] [2:,3] u1 2 [03] [2 2 [C3] [2
Looking at the constraints on the relative rotational motion between body 2 and body 3

imposed by the revolute joint:

s = b b 0 (7-70)

where b is the unit basis vector of the body-fixed coordinate frame in body 2 that is

aligned with the axis of rotation of the revolute joint between body 2 and body 3 and b& is

a unit basis vector of the body-fixed coordinate frame in body 3 that is perpendicular to

the unit basis vector of the body-fixed coordinate frame in body 3 that is aligned with the

axis of rotation of the revolute joint between body 2 and body 3. Writing all components
in terms of the inertial coordinate system yields:

s =b2 [C2T [3] b 0 (7-71)

The contribution to the Jacobian matrix from this constraint is:

=[ o0 00

0 T 2T [ 3
o o0 b CT [l)2]T [2]T [3] b (7 72)
S2 3(2 7-72)
2T [1 3 2T 1 2
2T ]T [3]1 [2:C3] b~3 T ]T [3]1 [2:33 b 3 2T [0 2T [3]1 [2:,33] b
3 [22 3 32 2 3 2 32 3 2

Similarly, the second rotational constraint imposed by the revolute joint between body 2

and body 3 is exemplified in:
9 b2 b = 0 (7-73)

where b3 is a unit basis vector of the body-fixed coordinate frame in body 3 that is

perpendicular to the unit basis vector of the body-fixed coordinate frame in body 3 that is

aligned with the axis or rotation of the revolute joint between body 2 and body 3. Writing
all components of Equation 7-73 in terms of the inertial system yields:

9 = b2T [2 [C3] b 1 0 (7-74)

Taking the partial derivative of Equation 7-74 with respect to the generalized coordinates

6i = 0 00

0 0U i [2[ ( 2\:( 3 3
o or ]bi L1]T [C2] T [C3]_b7
2T [1:21T T 3 23 T
bTL ] [2]T [C3] b bi 2 [ [LC]T [ C]T []_b 0T
0 2 33 1 21 3] -3
2b [0C2T [03] [23] b3 b2T [021T [03] [2: ,3 b 2T b [C iT [03] [2:3] b ]
3 [C2 3 2 3 1 032 3 1 3 2 3

Using the system Jacobian matrix and the system constraint matrix, derived in the above,
the system of equations illustrated in Equation 7-48 are ready for integration. There are

several methods to integrate these equations. The simplest involves direct integration.
[M] [D]T
[MA] (7-76)
[][ q [0]

q x (777)

F } (7 78)
Thus, Equation 7-48 is written as:

[Mx] qA F (7-79)

4 [Ci] [Oi[0]

C-0 [i [0 (3-86)
TC1], [ [l [0:]
[T / =[T (Q)] [i2w:]


1rp _RcBi + (C] [oi ] LH, +[CI] [O 1] [12] [T (a)] [a2iT T H CB2
[Ci] [ i] l [1 T ()]l [oi21iTH [O [0w] [ai] [T (,)] [ai]T [N2]j ) 6p +
([ci] [o] 1iIH + [oI] 0 [012 [T (a)1 [[a i]T H C B2+
[Ci] [01;] [a12] [T (a)] [a2iT U + [Ci,] [01;'] [ai2 [T (a)] [a2iT [N2];) 2o +
([C] [oL []] [BHi + [T0() [[ i [T ( H-CB2+

[Ci] [OLO] [ [a [a2i]T U2 + [Ci] [01;] [ai2 [T (a)] [0a2lT [N2] q) 2 +
OcI] [0L ]]Tb [CIB-H1 1 [a121 [T (a)] [a2 -H1-CB2

[C] [ (] 1 12 [T ( [a21i] + [CBi ; 12] [T (a)] [a12 ] [2 21] o+
I 1 2 irll ~ u 12L ; 121 2 2. 1 T r2 12L ; 21
[Ci]i] ai2l [T (a)] [i2i] [a2i]T [[N2] ) i6 + [i] [T (a)] [a2iT [N2] 6q2 (3-87)

Re-writing Eq. 3-87 in terms of the B1 basis and using the matrices introduced in the
derivation of the contribution to the system kinetic energy from body 2:


6q 2
Performing the dot product of the forces acting on body 2 with the corresponding virtual
displacement vector of the form of Eq. 3-88 and comparing the result with Eq. 3-77
produces the contributions to the generalized forces from the forces acting on body 2.
The derivation of the contribution to the generalized forces from the forces acting on
body 3 proceeds from the principle of virtual work. First, the position vector to the point

Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy



Sean James Sidney Regisford

May 2008

('C! i: Richard Lind, Jr.
Cochair: Andrew Kurdila
Major: Aerospace Engineering

This work presents a methodology to determine of accurate computer models of

articulated micro air vehicles. Prior to this work, much of the innovation in design

of a micro air vehicle was based on empirical data. This emphasis on empirical data

necessitates a time intensive iterative process where a design is built and tested, then

modifications are made and the process starts anew. This work looks to alleviate some

of the load of this iteration by producing accurate models that can be inserted into the

process at greatly reduced cost in terms of time. The modeling is accomplished by the

application of Lagrange's equations to create a set of equations that incorporate both the

rigid and flexible degrees of freedom associated with the materials that Micro Air Vehicles

( I!AVs) employ. These equations of motion also include the generalized forces that result

from the body forces due to gravity and external forces acting as a result of the interaction

between the structure and the air through which it passes. These equations are formulated

in a manner that accomplishes the above while allowing for the reconfiguration of the

system. The equations of motion are then transformed into a form amenable for use in the

solution of an optimal control problem.

[91] Bottasso, C.L. and Croce, A., "Optimal Control of Multibody Systems Using an
Energy Preserving Direct Transcription Method", Multibody System D.;;,i,.-
Vol.12, 2004, pp.17-45.

[92] Benson, D., "A Gauss Pseudospectral Transcription for Optimal Control", Ph.D.
Thesis, Massachustts Institute of Technology 2004.


[1] Shabana, A.A., D;in,.i of Multibody S';,-. n- New York NY: John Wiley & Sons,

[2] Meirovitch, L. and Stemple, T., "Hybrid State Equations of Motion for Flexible
Bodies in Terms of Quasi-Coordi ii. ,Journal of Guidance, Control, and D;.i,,.i-
ics, Vol. 18, NO. 4, 1995, pp.678-688.

[3] Meirovitch, L. and Tuzcu, I., liI,. ,i I1,ed Approach to the Dynamics and
Control of Maneuvering Flexible Aircraft",NASA Technical Paper, Paper No.
NASA/CR-2003-211748, 2003.

[4] Meirovitch, L. and Tuzcu, I., "Unified Theory for the Dynamics and Control of
Maneuvering Flexible Aircraft", AIAA Journal, Vol.42, No.4, April 2004,pp.714-727.

[5] Meirovitch, L. and Tuzcu, I., "Time Simulations of the Response of Maneuvering
Flexible Aircraft", Journal of Guidance, Control, and Diimi. Vol.27, No.5, 2004,

[6] Kwak, M.K. and Meirovitch, L., \. v- Approach to the Maneuvering and Control of
Flexible Multibody Systems", Journal of Guidance, Control, and Diimi. Vol.15,
No.6, 1992, pp.1342-1353.

[7] Gad-el-Hak, M., i\! ro-Air-Vehicles: Can They be Controlled Better?", Journal of
Aircraft, Vol.38, No.3, 2001, pp.419-429.

[8] Lian, Y., Shyy W., Viieru, D., and Z!:i .1- B., i\. i:ioi ,i,. Wing Aerodynamics for
Micro Air \. !i. 1! Progress in Aerospace Sciences, Vol.39, 2003, pp.425-465.

[9] Mueller, T.J. and DeLaurier, J.D., "Aerodynamics of Small \. !hi !. Annual
Review of Fluid Mechanics, Vol.35, 2003, pp.89-111.

[10] Muti Lin, J.C. and Pauley, L.L., "Low-Reynolds-Number Separation on an
Airfoil",AIAA Journal, Vol.34, No.8, 1996, pp.1570-1577.

[11] Lutz, T., Wurz, W., and Wagner, S., "Numerical Optimization and Wind-Tunnel
Testing of Low Reynolds-Number Airfoils", Conference on Fixed, Fir,,.:,j and
Rotary Wing Vehicles at very Low R. ;,.,....1- Numbers Notre Dame, Indiana 2000.

[12] Ifju, P.G., Ettinger, S., Jenkins, D.A., and Martinez, L., "Composite Materials
for Micro Air V\ I, !. S .. .. I'; for the Advancement of Materials and Process
Engineering Annual Conference, May 2001.

[13] Ifju, P.G., Jenkins, D.A., Ettinger, S., Lian, Y., Shyy, W., and Waszak, M.R.,
"Flexible-Wing-Based Micro Air \. !i 1!, American Institute of Aeronautics and
Astronautics Paper No. AIAA2002-0705, 2002.

[ ] [12] [T (a)] [a ] rB2-H2 +
[o0w] [012] [T (a)] [a2] [a23] [T (03)] [T H2-CB3 +
[Owl] [ 12] [T (a)] [a2l]T [a23] [T (03)] [a32]T +

[O] la2 [T [(a)] [2T rB2-H2 +
[0O] [012] [T (a)] [a21l] K ] [T (8)] r rH2-CB3 +
[ l]a2 [T (a)] [a21]T ,.] [T )] T OP
[0o] [012] [T (a)] [a21]T,. ] [T (3)] ]T [N3] q3 +

S([a121 [T (a)] [ 12 21] [21] CB2 + [Ha12] [T(a)] [1221] [a2]T rB2-H2+
[a12] [T (a)] [12L21] [a21 [ [T (3)] [( ] rH2 -CB3 +
[a12 T )] [ 1221] [ 21 [a23] [T (3)] [a32] o +
[a12] [T (a)] [12 21] [a21T [a23] [ (03)] [a32] [N3] q)
) ([a12] [T (a)] [a21] [a23] [T (3)] [2332] [32]T rH2CB3+
[012] [T(a)] [I21] [a23 [T(03)] [2332] [a32 ]T
[a12] [T (a)] [a2] [a23] [T (3)] [ 323 ] a32 ]T [N3] qf) (3-63)


BB = [o] rBl-Hi + [w] [a12] [T (a)] [21] TarH-CB2 +
[o1] [al2] [T (a)] [a21]T rB2-H2 +
[0 [a 2] [T] (a)] [a21]T [a23] [T (3)] [a3]T 2-CB3 +
[oL] [2 [T0 (a)] [a21T [a23] [T (/3)] [a3a]
[O [C 12] [T2] (a)] [a21]T [23] [T (3)] [a32]T [N3] q 3 (3-64)

BB2 [0 B1] _rB1-Hi + [01w] [a2] [T (a)] [a21]T H1-CB2 +

Since, [4], has full row rank, [MA] is invertible and thus

q = [A]-1F (7-80)

which gives the simultaneous solution for the second time derivative of the generalized

coordinates and the Lagrange multipliers. The solution of the second time derivative

of the generalized coordinates can then be used in a first order system to solve for the

generalized coordinates at the next time step. To write the equations of motion as a first

order system, first let:

X= q ( (7-81)


X = (7-82)
S[]-1 Q T ( U (7 82)

where u represents a control input to the system. Figure 7-1 demonstrates the solution

of the equations of motion to simulate the system after one second in which a control is

applied to bi. There is a known problem with the direct integration of the equations of

motion. This problem is that the solution of the Lagrange multipliers is not exact. As the

simulation time grows, these errors also grow. As seen in the figure, the norm of the vector

that represents the constraint vector for the system grows as the simulation progresses.
There is a stabilization method that reduces the errors associated with the solution of

the augmented system illustrated in Equation 7-79. This process was first published in

the work by Baumgarte[82]. In the stabilization method, the objective is to force the

constraint vector, the first time derivative of the constraint vector and the second time

derivative of the constraint vector to zero simultaneously. 83--851 This objective is reached

by making the constraint equations analogous to the standard second-order ordinary

differential equation of control theory. This analog is represented by:

K + 2aK + 32 0 (7-83)

and 0 D, where appropriate.

Sv BDidV
S[B]T DIdV = f B DidV (A-8)
fS B DidV

Examining a representative term,

0-1 12(21
rBi-H1,i L ,o, T ,km W ,, TH1-CB2,p dV +

rB1-Hl,i 3, ji(' m a12 mn, (I / opp dV +
01 1221 2dV
TBl-Hs L jT (a), a,krn21,nmrfp N2 dV +

TH1-CB2,wa21,wrT (a,, a12,is T (a) ,km 12nn( rHl-CB2,p jdV +
TH1-CB2,wa21,wrT (a),sr 12,iso i T (a) ,k n 12 n I dV +

0 12 1 u dV +2
T21,wrT (a),,, a12,is () ,km mn(' rH1-CB2,p o

a21,Twr (a),, a12, is0 1 i1 T () -k 1 opp dV +
6 aT JV

T21,vwT, (a) 12,sr1 i2, ,T (a), j km 2nmqk p j opNp dV +
0( 1 )121 /
q2 .1 .,rT (a),sr o 12,is 0 ji T (a) ,kn n ( I f,nrH1-CB2 N V ,w dV +
22212 f1 2 dV)
q ... T (a), 2,is 0 1 .T ( a) 12 21 2,w opp dV +

q .wT (a) a12,sr 0 1 (aJk12 1 km lm N2,v 2,np dV (A-9)


The derivations of the equations of motion that were derived in ('i! iter 3 were

applied to a numerical example to show that the equations can then be cast into a form

amenable to an optimal control problem. This example is implemented in MATLAB. The

first step in forming this example is generating a geometric model. The geometric model is

constructed to exhibit the prominent features of the actual model, such as the articulation

of the wing,but without the complexities associate with the full model. Issues such as joint

friction and composite material properties are ignored without loss of generality.

6.1 Model Parameters

As illustrated in Fig. 6-1 the fuselage is modeled as a rectangular prism and the wing

segments as well as the horizontal and vertical tail are modeled as rectangular plates.

Although the idealized model is a simplification of the actual articulated MAV, the lengths

and masses of the idealization was chosen to reflect the actual characteristics of the full

vehicle. The model fuselage has a length of 0.3048 meters and a width and height of

0.0636 meters. The inboard wing segments have a span of 0.2032 meters and a chord

length of 0.0762 meters. The outboard wing segments have a span of 0.1016 meters and a

chord length of 0.0762 meters. The horizontal tail segments have a span of 0.127 meters

and a chord of 0.06858 meters. The vertical tail segment has a span of 0.1016 meters and

a chord of 0.06858 meters. All rectangular segments have a thickness of 0.012192 meters.

The density of all components is assumed to be 137.47k to give the total system a mass

of 0.29031 kg.

To further simplify the numerical example, flexibility effects are ignored. Eliminating

flexibility has the effect of reducing the generalized coordinates as well as simplifying

the generation of the matrices involved in the equations of motion. This simplification

does not invalidate the observations that result from the analysis of the idealized model.

Rather, the exclusion of flexibility effects greatly reduces computational effort. For

which can be written in the familiar quadratic form:

T M4 (7-99)

where q = i and using the indicial notation,

M 0oi-jy j, pp, u dV (7 100)

The derivation of the contribution to the kinetic energy of the system from body 2 begins

with the definition of the vector from the origin of the inertial system to a random point

on body 2. This vector is written as:

Lp = RCBl + LB1-HI -H1-CB2 + p (7-101)

,where TBI_-H is the vector from the center of mass of body 1 to the revolute joint

between body 1 and body 2, written in terms of the body-fixed coordinate frame

associated with body 1, rHI-CB2 is the vector from the revolute joint between body 1

and body 2 to the center of mass of body 2, written in terms of the body-fixed coordinate

frame associated with body 2, and u, is the vector from the center of mass of body 2

to the random point on body 2, written in terms of the body-fixed coordinate frame

associated with body 2. Taking the first time derivative of Equation 7-101 with respect to

the inertial coordinate frame yields:

E d d Bid + EwLB1 Hir
dt (E ) E (CB1) dt B1-H) + BI-H
B2d (rH1-2) + E xB2 + B2 d (B) + Ep7-102)

E B2 B1 + B1 B2 (7 103)

the inboard wing centered at the center of mass of the inboard wing, u is the vector from
the center of mass of the inboard wing to the undeformed point p, and [N2] q2 is the vector
from the undeformed point p to the deformed point p. For the following derivation, RCBI
is in terms of the inertial frame, rB1-H1 is in terms of the B1 basis, rHI-CB2 is in terms of
the B2 basis, u is in terms of the B2 basis, and [N2] 92 is in terms of the B2 basis. Thus

Ed () RCBI + B ( -Hi) + B X B-H1 + 2 H1-CB) +
(L X _) -H- +B2 + x u^ + B d[ ([N,] q +
EwB2 H _- B2 d ( B+ B + 22)

E B2 x [N2 (3-32)

By virtue of some of the vectors being constant in the frame that the derivative is taken in

E RCB +W1 EB2 E BCB2+ B2 d ([N12] q) + B2 X [12]
CB XB B1-H1 E X H1-CB2 d \ J /
Through the use of simple angular rotations between frames,

E ,B2 E LLB1 BI J1 J12 ,J21 + J21 L B2 (3 34)

However, B1 U2 j21 UB2 0. Thus, separating the angular velocity terms and re-writing
all cross products in terms of the same basis yields:

VP RCBI + V1 x B1 Hi + EB x [a121 [T (a)1 [ITi LHI-CB2 +
2 >21x + EHi CB +EWB-I x [a12] [T (a)I [acl+ u +

J1221 [a21iT [N2] (3-35)

Let [OZ1] be a skew symmetric matrix composed of the components of EWB1 and [12 21] be
a skew symmetric matrix composed of the components of J12w2 Thus,

R CBI + [Oi] BI-Hi + [LO1] [Cai2] [T (a)] [ailT HI-CB2

rH1-CB2,uT (02),vu 0w 2 ,y T y H2-CB3,z dV +

rH-CB2, 2) 0w-1 0- 1 2 w dV +
rB2-H2,uT 2 ( w 01,y0 -2,yzHjB1-Hl,z dV+
rB2-H2,vT (02), w, 0 (/2),yz (tjH1-CB2,z dV +
6 6 1,yzrB -H
rB2-H2,vT ( T2),W O1,W0 w T (yT(02),yr rHC-CB2,z dV +
,x r HfJV

rH2-H2,T 02),Wv T ( L 0 1 0w1 T +2V

S 0-1 0 dV +

rB2-H2,uT (02),v 01, 0 1,wT 0y T (3 ),y rH2-CB3,z dV dV +
6 6 1 ,JJ v

TH2-B3), 3~ ,u 30w 1 0- 1( T (2),yz TH-CB32,z dV +

TH2-2CB3,T (3),t (3) 0-1 0 -1 T3,yz H2,CB3,z d dV

rH2-CB3,vT (33) ,w T (23), o X ,yz BlHlz op,vz dV +

rH2-CB3,uT (3),v T 101 T ( 2 ,yz H -CB2,z rHop, dV +
T T1Jv

rH2-CB3,T ( vvu T v O 1 2 T2,yz B2-H2,z op,u dV +
rH2 3,T (3 3),,.t 2 ( v), O)1,W 2) T (y3),yz TH2-C3,z op,t dV +
; x 0-1 0-1 T U V+

T ( 3),t T 2) 0 1,w I, ,y T ( ,yz op,topz dV (7-126)

Here again, the remaining integrals in Equation 7-126 are invariant over the body.
Combining the kinetic energy contributions from body 1, body 2, and body 3 yields
the kinetic energy of the entire system. Then, inserting this kinetic energy term into the
Lagrange's equations produces the equations of motion of the system which can be stated

Figure 7-3: Three bar manipulator.

[14] Hall, J., Mohseni, K., Lawrence, D., and Geuzaine, P., l -is i:, i of Variable
Wing-Sweep for Applications in Micro Air V\ !i 1!, American Institute of
Aeronautics and Astronautics Paper No.AIAA2005-7171, 2005.

[15] Ramamurti, R. and Sandberg, W., "Simulation of the Dynamics of Micro-Air
V\ !i-, !I 38th Aerospace Sciences Meeting 6 Exhibit Reno, Nevada, Paper
No.AIAA2000-0896, 2000.

[16] Raney, D.L. and Slominski, E.C., \ I. Ii ,i~. ,ii ,i and Control Concepts for
Biologically Inspired Micro Air V\ 1!. Journal of Aircraft, Vol.41, No.6, 2004,

[17] Schenato, L., Campolo, D., and Sastry, S., "Controllability Issues in Flapping Flight
for Biomimetic Micro Air Vehicles(i\ AVs)" ,Proceedings of IEEE Conference on
Decision and Control, 2003.

[18] Abdulrahim, M., Garcia, H., Ivey, G.F., and Lind, R., "Flight Testing a Micro Air
Vehicle Using Morphing for Aeroservoelastic Control", AIAA Structures, Structural
DIm.II:. and Materials Conference, Palm Springs, CA, Paper No.AIAA2004-1674,

[19] Waszak, M.R., Davidson, J.B., and Ifju, P.G., "Simulation and Flight Control of an
Aeroelastic Fixed Wing Micro Air Vehicle", AIAA Atmospheric Flight Mechanics
Conference, Monterey, CA, Paper No. AIAA2002-4875, 2002.

[20] Waszak, M.R., Jenkins, L.N., and Ifju, P., "Stability and Control Properties of an
Aeroelastic Fixed Wing Micro Air Vehicle", AIAA Atmospheric Flight Mechanics
Conference, Montreal, Canada, Paper No.AIAA2001-4005, 2005.

[21] Lian, Y., Shyy, W., and Haftka, R.T., "Shape Optimization of a Membrane Wing for
Micro Air V\ I!.. AIAA Journal, Vol.42, No.2, 2003, pp.424-426.

[22] Lian, Y., Shyy, W., Ifju, P.G., and Verron, E., \. il.i ,i. Wing Model of Micro Air
V\ !. !, AIAA Journal, Vol.41, No.12, 2003, pp.2492-2494.

[23] Lian, Y. and Shyy, W., "Numerical Simulations of Membrane Wing Aerodynamics
for Micro Air Vehicle Applications", Journal of Aircraft, Vol.42, No.4, 2005,

[24] R in, :lr'vani, D.S., Lesieutre, G.A., Frecker, M., and Bharti, S., "Aircraft Structural
Morphing Using Tendon-Actuated Compliant Cellular Ti, --. Journal of Aircraft,
Vol.42, No.6, 2005, pp.1615-1621.

[25] Bae, J., Seigler, T.M., and Inman, D.J., "Aerodynamic and Static Aeroelastic
('C!i lteristics of Variable-Span Morphini Journal of Aircraft, Vol. 42, No.2, 2005,

cb2CO2 -sQ2C02 + c2AsO02s 2 sQ2Qs2 + cP2s020c 2
[C21 _-C2 _( C_ + ',_ ( ', + .-', (7-5)
-s02 c02Se 2 c02CQ2

( 3 3 -1 (' '103'3 ,' 3S 3 3 + (' '03C 3
[C3] 3 ( ,' .)3 + ,, + ,' )3 (7-6)
-s03 c03s 3 c03sc3
where cQi is cos Qi and sQ is sin Qi. As all three bodies have similar relations between
their respective body-fixed frame and the inertial frame, they have similar relations for
the angular velocity of those body-fixed frames with respect to the inertial frame. These
relations are:

E LB (4j +1O1) bt + (Oico + iCOIS1) bl + (-Oisoi + QicHicOi) b (7-7)

E1B2 (= 2 .) + ( i022 + 't0'C0252) b2 + (02SQ2 + '_C 02C62) b6 (7-8)

E B3 (43- 03) b3 + (33CQ3 + 0* '3 b) 3 + (-O3S3 + 0c) 3 (7-9)

where the angular velocities are written in terms of the basis vectors of the body-fixed
coordinate system.
The derivation of the equations of motion starts with application of Lagrange's
equations. For a rigid body, Lagrange's equations takes the form of:

dt ) Oq (7 10)

Where T is the kinetic energy of the system, q are the generalized coordinates, q are the
first time derivative of the generalized coordinates, and QE are the generalized forces due
to external forces exerted on the system. The kinetic energy of the system is the sum of
the kinetic energies of the parts that comprise the system. Looking at bar 1, then kinetic
energy is:
T = j p (vp. ) dV (7-11)
Jv 2

of the force application on body 3 is written as:

S CBI + B-HI + rH1-CB2 + -B2-H2 + rH2-CB3 + op + (389)

Re-writing Eq. 3-89 so that all components are in terms of the inertial frame yields

rp = RCB1 + [Ci] rB1-H1 + [Ci] [a2] [T (a)] [i] rHI-CB2 +

[Oi] [l12] [T (0)] [a21]T B2-H2 +

[C] [la121 [T (a)] f ]T [, ] [T (/3) [2 rH2-CB3 +

[CI] [a121 [T (a)] [c21T [o ] [ () [32 ]To +

[Ci] [a121 [T ( [21 [ [ ( [2 [ 3)] [a32]T [N3] q3f (3 90)

Taking the virtual displacement proceeds as:

6r 6RCI cB+ a([C1]rBI H1)
a([C1]LB1 -H1) + a ([CI]l-H1)
--I--6 + +
a ([C1] [a12] [T (a)] [a2]T H1- B2 ) a ([Cl] [1121 [T (a)] [a21]T H1-0B2)
p+ -60+

a ([C1] [012] [T (a)] [a21]T _B2H 2) a ([C1] [a12] [T (a)] [21]T -r B2 2 )
69----------' Sp+ 6-----------(a +

a ([C1] [012] [T (a)] 2 ([C]] [12] [[2 ()] [a [22]T] 2B2 H2)

S([Cl1 [12] [T (a)] [21] 2+ O-S +

a ([c1] [a12] [T (a)] [021]T [a23 Er (8)] [a32] TH2-CB3)
a ([C lj [a12] [T (a)] [a2l]T [a231 [T (/3)1 [0 32]TTr /2-C 3)

9 ([C1] [012] [T (a)] [a21]T [023] [T (3)] [a32 ]T H2-CB3)
a ([C] [012] [T (a)] [021]T [0231 ]T (/)] [a32]1 TH_2 op 3)

a ([C1] [012] [T (a)] [a2]T [023] [T ()] [032]T H2 -03)
68 +

a ([C1] [a121 [T (a)] [2j]T [a231 T ()] [a32]T_]a2 0S)
a ((c)] ]a12] 1] (a)] [a2ljT []a13 ]T (8)] []32]T_ !u )

a--60 +
o~ ~ ~ ~ ~ ~~~~~~~S ([j[11I 1[lT[1I /) T. 6

generalized coordinate. The preceding information is sufficient to define all of the

quantities on the left hand side of the Equation 5-10. Next the appropriate quantities

on the right hand side of the equations of motion must be defined. The gravitational

forces are applied to each multi-body segment using the masses of each segment and local

acceleration of gravity. These gravitational forces are applied to the center of gravity of

the respective segments. To define the propulsive forces, AVL is executed using a geometry

input file that corresponds to no deflection of both the inboard and outboard segments. In

addition, the velocity is defined in AVL to be 15I. This velocity is chosen to be associated
with the MAV flying in a trim condition in straight and level flight. Thus, the coefficients

generated by this configuration should produce a lift force that exactly balances the weight

force and a drag force that balances the thrust force. The magnitude of this drag force

is then used as the magnitude of the propulsive force. Following the derivation of the

generalized forces as outlined in ('!i lpter 3 the thrust force is associated with acting on

the fuselage. The direction of the force is assumed to act through the center of gravity

of the entire system. The center of gravity is found to l. i lie in the fuselage in all

configurations of the vehicle that are admissible as defined by the constraints on the angles

set forth in ('!i lpter 5. For open loop simulation, the angles of the inboard and outboard

wing segments are constrained to be zero.

Now all of the necessary components for an open loop simulation of the vehicle are

present. It was chosen to use a Euler integration scheme to integrate the equations of

motion. The Euler integration scheme has the advantage of only needing information from

the current time step in order to calculate the values at the next time step. However, care

must be taken in the determination of the time step to ensure convergence of the solution.

To ensure convergence, the time step must be as follows:

At < At,,


u Q] T Q = U (7-127)

U = C2 (7-128)

d (T[[
d ([ 1)D + t a (T (7 129)
-- dt 2 O"qS sys -8

ss 2 (7-130)

The alternative formulation of the equations of motion illustrated in the above results in

a system that has only three generalized coordinates. These coordinates are the angles

that represent the relative motion of the bodies in the system with respect to each other.

In the standard formulation, the number of generalized coordinates is greatly increased

due to the necessity to completely describe the motion of each body instead of the relative

motion of the bodies with respect to each other which is the result of the alternative

formulation. The reduced number of generalized coordinates has the primary advantage

of the reduction in computational cost when integrating the equations of motion. As

more bodies comprise the system, the reduction in computational cost by using the

alternative formulation of the equations of motion becomes more evident. Also included

in the computational savings the alternative formulation demonstrates is the ease of

integration of the equations of motion in comparison to the integration of the equations

of motion for the standard formulation. The equations of motion as derived using the

alternative formulation can be directly integrated without having to first simultaneously

solve for Lagrange constraints, as is necessary in the standard formulation of the equations

of motion. However, the standard formulation has the advantage of each new body is


In the previous chapter, the contributions to the kinetic energy and the potential

energy from the various bodies that compose the system were derived. The generalized

form of the contributions to the generalized forces from the various forces acting on these

bodies was also derived. The goal of this section is to combine these components into a

system of equations that govern the motion of the entire vehicle and use this formulation

in the derivation of a control strategy that will be used to drive the vehicle to a desired

state. An optimal control strategy can be used to address a number of control objectives,

from tracking to regulation.

There are several different v--,v in which the optimal problem can be solved.

However, these methods generally fall within two classifications. These classifications

are direct methods and indirect methods. Much work has been done in the use of indirect

methods in the solution of optimal control problems.[58--67,69] Indirect methods use the

calculus of variation in the derivation of the optimal control. First, the objective function

is written as:

J (x, U) (x (tf), tf)+/ (x (T) u (T) ,) d (5-1)
Where Q (x (tf) t ) is some metric of the states at the final time, and f, (x (r), u (r), r) dr

is the Lagrangian function that represents a metric of accumulated errors over the time

span. The solution of the optimal control problem is concerned with minimizing this cost

functional. The solution that minimizes this cost functional, however, must also satisfy

a set of equations of motion. For typical optimal control problems, these equations of

motion take the form of:
(t) f (x (t) ,u(t) ,t) (52)

x (0) = 0

where x (t) is an N-dimensional vector of states and u (t) is an M-dimensional vector

of control inputs. Equation 5-2 shows that the states x is a function of the control u.

E B3 E _B1 J12 21 J23 32 (3-58)


Vp C BI + EB X rB1-H1 + (EI + J2 1) X rH-CB2 +
E B J2 )2 X rB2- H2 E J21 J 23 + 32) x H2-CB3 +
E B J2 J 2 J23 j32 3
(EBL1 + L V3+ )x P + [N3] _j +
(E B + ) J 21 + J23a2) x [N3] (3-59)

Expanding the cross products and expressing the components of the cross products in the
same bases yields

Up RCB1 + X x B1H1i + x [a12] [T (a)] [21]T H1-CB2 +
J121 X [a2 21] rH1-CB2 + E 1 x [O12] [T (a)] [2i]T- B2-H2

P2 21 x [a221]T B2H2 + EB X [Ca12 [ (C [1T [o [T (/3) [] rH2-CB3a

2 2 x [a21]T [(o ]i [T (/3)] [] H2-CB3 + 33 X [( T rH2-CB3 +

EB x [Q12] [T (a)] [fT[ ] [TQ3)] [o .f 1
12 )21 x [ a21 ] [T [(/03) [32 ]T + J23W32 x [ _2] T +
[12 3.,B 2J21 -23 + 32 j2

12 .j21 x ]2 2 33 J23 .j 32 + 933 60)
[N3} f + EB1 x [ai12 [T (Ia)] [a2ll [a231 [T (3)] [321T [N3] qf

f f
J12 1 X [c21]T [ad23] [T (/3)] [Q32]T [N3] q + 23 3 x [c32]T [N3] q (3-60)

Introducing skew symmetric matrices to simplify the cross product expressions and
re-writing the whole equation to be in terms of the B1 basis yields:

Up [CR1T CB1 + [1)] rBH1 i 1w1] [12 [T (a)] [a21] rTHICB2 +

[012] [T (a)] [1221] [21a T rH-CB2 [0: 1] [a12] [T (a)] [a2T-rB2-H2 +

[a 12 [T ( 12)21] [21 ]TB2-H2 +

[Owj1] [I12 [T (I)] [a21]T -[T (/3)] IrH2-CB3 +

[39] Boutin, B.A. and Misra, A.K., "Dynamics of Manipulating Truss Structures",
Collection of Technical Papers, AIAA/AAS Astr '.;i,,:. Conference, San Diego,
CA, Paper No. A96-34712-09-12, July 1996, pp.471-4793.

[40] de Jalon, J.G., "Twenty-Five Years of Natural Coordi ii., ,Multibody System
D i,.i, Vol. 18, 2007, pp.15-33.

[41] Leger, M. and McPhee, J., "Selection of Modeling Coordinates for Forward Dynamic
Multibody Simulations", Multibody System Dn,,ni,.- Vol. 18, 2007, pp.277-297.

[42] Anderson, K.S., and Duan, S., "Highly Parallelizable Low-order Dynamics
Simulation Algorithm for Multi-Rigid-Body Systems", Journal of Guidance, Control,
and D.,u,,. Vol.23, No.2, 2000, pp.355-364.

[43] Banerjee, A.K., "Contributions of Multibody Dynamics to Space Flight: A Brief
Revi. Journal of Guidance, Control, and Diiii,. Vol.26, No.3, 2003,
p11' ;-'-394.

[44] C('! i, S., Hansen J.M., Tortorelli, D.A., "Unconditionally Energy Stable Implicit
Time Integration: Application to Multibody System Analysis and D. -1 i1 Interna-
tional Journal for Numerical Methods in Engineering, Vol.48, 2000, pp.791-822.

[45] C('1!, S., Tortorelli, D.A., "Optimization of Multibody Systems Via an Energy
Conserving Minimal Coordinate Formulation", American Institute of Aeronautics
and Astronautics, Paper No. AIAA-98-4927, 1998.

[46] Schaefer, E., Loesch, M., Rulka, W., "Real-Tiem Simulation of Flexible Space
Manipulator Dyir 1i.1 American Institute of Aeronautics and Astronautics, Paper
No. AIAA97-37062, 1997.

[47] Kielau, G. and Maiber, P., Nuinholonomic Multibody Dyir 1,i11 Multibody System
D i, ,. .- Vol. 9, 2003, pp.213-236.

[48] Aoustin, Y. and Formalsky, A., "On the Feedforward Torques and Reference
Trajectory for Flexible Two-Link Arm", Multibody System Di',,,ii. -Vol. 3, 1999,

[49] Junkins, J.L., "Adventures on the Interface of Dynamics and Control", AIAA, 35th
Aerospace Sciences Meeting 6 Exhibit, Reno, NV, Jan.6-9, 1997.

[50] Rives, J., Portigliotti, S., and Leveugle, T., "Atmospheric Re-Entry Demonstrator
Descent and Recovery Sub-System Qualification Test Flight Dynamics and
Trajectory Modeling during D. -.i. i American Institue of Aeronautics and As-
tronautics, Paper No. AIAA-97-1441, 1997.

[51] Tadikonda, S.S.K., "Modeling of Translational Motion Between Two Flexible Bodies
Connected via Three P.i l- Journal of Guidance, Control, and Di,,,ii. Vol.18,
No.6, 1995, pp.1392-1397.

[Ci] [ow] [a12] [T (a)] [a21]T [d23] [T (03)] [a32]T Uop+
[C ] [ 2] [a12T (a)] [I21]T [23] [T (3)] [a32]T [N3] q) 60 +
([C,] [0:J] rBIH1+
[Ci] [0w ] [a2] [0 T (a)] [ai]T rH-CB2 + [C [[0 '][a [o] T (a)] [ fT rB2-H2 +
[Ci] [o)] [012] [T (a)] [21T [ [T (3)] [(I 2-CB3 +
[Ci] [o ] [a12] [T (a)] [a21] [o I [T (3)] [ o(I ] +
[Ci] [o)] [a12] [T (a)] [a21T [ IT ()] [( J [N3] _) q '3 +
([c,] [,12] [T (a)] [12 1] [ 17,T .H -B2+ [C] [a12 [T(a)] [121] [ T rB2-H2+
[Ci] [)] [121 [T (a)] 21] [a21] ] [T ()],, ]T r 2-CB3 +
[Ci] [] [T0 (a)] [1221] [a21i]TK [T (j3)] [, ]jT +
[Ci] [a12] [T (a)] [12a1] [a2,]T K:. [T (03)] : jT [N3] q) 6a +
([Ci] [a12] [T (a)] [a21] [a23 [T ()] 2] [a32]T r2CB3+
[Ci] [a12] [T (a)] [a21] [a23] [T (3)] [23w2] [a32]TU +
[C] [a12 T ( [a [023] [T (03)] [232] [a32]T [N3] q3) 6, +
[Ci] [a12] [T (a)] [a21 [23] [T (3)] [23 2] [a32]T [N3] 6q; (3-93)

Using the vectors and matrices defined in the derivation of the contribution to the kinetic
energy from body 3, and re-writing the virtual displacement of a point on body 3 in terms
of the B1 basis yields:

,p [C 1] [BB]DDDD2[a2 [T (a)] [21]T [o ][ (3)] [32]T [N3] 6a

(f 9

Based on these assumptions, the derivation of a numerical example is divided into

two sections that follow. First the equations of motion are derived and a simulation

conducted based on the standard multibody dynamics procedure, where each rigid body

in three-dimensional space has six degrees of freedom. Then the equations of motion are

derived in the same manner as illustrated in C'! lpter 3. In using these two formulations

the strengths and weaknesses can be contrasted.

7.2 Standard Formulation

In the following, bars one, two, and three are each assumed to be rigid bodies and

as such have six degrees of freedom associated with their movement. Each body will

have a body-fixed coordinate system and these systems are all related to a global-inertial

coordinate system that is used to describe the motion of the system as a whole. These

coordinate systems are illustrated in 7-4. Mathematically, the body-fixed coordinate

systems are related to the inertial coordinate system by the use of a transformation

matrix. Each transformation matrix, is itself a function of three angles. Thus,

e, = [CI] b' (7-1)

e2 [C2] b2 (7-2)

e3 [C3] b3 (7-3)

where e,, e2, and e3 are the inertial frame version of the vectors b], b, and b3 which are

vectors written in terms of the body-fixed reference frame of bodies one, two, and three

respectively. Also, the transformation matrices are constructed using a 1-2-3 Euler angle.

cylcOl -syicl + ceplsOlspl sQlsol + cep1sO1cp1

[Ci] = syic1i c ico1 + ssO1isk -c isoi + s isO1ici (7-4)

-so1 cOisQi c1icQi

[52] Tadikonda, S.S.K., "Articulated, Flexible Multibody Dynamics Modeling:
Geostationary Operational Environmental Satellite-8 Case Study", Journal of
Guidance, Control, and D;.,,.i,, .:. Vol.20, No.2, 1997, pp.276-283.

[53] Jain, A., and Rodriguez, G., "Recursive Dynamics Algorithm for Multibody Systems
with Prescribed Motion", Journal of Guidance, Control, and D;i.,. Vol.16,
No.5, 1993, pp.830-837.

[54] Pradhan, S., Modi, V.J., and Misra, A.K., "Order N Formulation for Flexible
Multibody Systems in Tree Topology: Lagrangian Appr.o !I Journal of Guidance,
Control, and D; ,an.. Vol.20, No.4, 1997, pp.665-672.

[55] Attia, H., "Dynamic Simulation of Constrained Mechanical Systems Using Recursive
Projection Al, i..i I ii Journal of Brazilian S... :,. ;.i of Mechanical Sciences and
Engineering, Vol.28, 2006, pp.37-44.

[56] Dias, J.M.P. and Pereira, M.S., "Sensitivity Analysis of Rigid-Flexible Multibody
Systems", Multibody System D m;,,i. Vol. 1, 1977, pp.303-322.

[57] Ortiz, J.L. and Barhorst, A.A., "Modeling Fluid-Structure Interaction", Journal of
Guidance, Control, and D ;,,in. Vol.20, No.6, 1997, pp.1221-1228.

[58] Paiewonsky, B., "Optimal Control: A Review of Theory and Practice", AIAA
Journal, Vol.3, No.ll, 1965, pp.1985-2006.

[59] Hodges, D.H. and Bless, R.R., \V. I: Hamiltonian Finite Element Method for
Optimal Control Pil!. !ii- Journal of Guidance, Control, and D.;in..ii Vol.14,
No.l, 1991, pp.148-156.

[60] Schaub, H. and Junkins, J.L., X. .- Penalty Functions and Optimal Control
Formulation for Spacecraft Attitude Control P..I ii,- Journal of Guidance,
Control, and D; ,an.. Vol.20, No.3, 1997, pp.428-434.

[61] Seywald, H. and Kumar, R.R., "Finite Difference Scheme for Automatic Costate
Calculation", Journal of Guidance, Control, and D,,iiii.. Vol.19, No.l, 1996,

[62] Seywald, H. and Kumar, R.R., "A Method for Automatic Costate Calculation",
AIAA, Guidance, Navigation, and Control Conference, San Diego, CA, Paper No.
AIAA-96-3699, 1996.

[63] Seywald, H. and Kumar, R.R., \I. I !.)d for Automatic Costate Calculation",
Journal of Guidance, Control, and D;,i,.ii. Vol.19, No.6, 1996, pp.1252-1261.

[64] Venkataraman, P., "Dynamics and Optimal Control", AIAA Guidance, Navigation,
and Control Conference and Exhibit, Dever, CO, Paper No. AIAA-2000-4560, 2000.

quadrature used here,
1 1

1 1
I = v< 3 > v=< > (6-3)
1 1
V3 V3
1 1

/ V = p opp det J t (6-4)

where op,i is the position of the quadrature point. Lastly,


/ /'I" Up,pdV p op, Up det J t (6-5)
i= 1

Now that the invariants are calculated, the current configuration of the model is used

to determine the various quantities in the equations of motion. For the purposes of this

numerical example, the configuration dependence is simplified as much as possible. For

example, in the initial configuration the coordinate systems associated with each body

is aligned with the global coordinate system. Furthermore, the wing segment coordinate

frames are aligned with the respective joint coordinate frame located on that body. Thus,

the only change in coordinate frames are determined by the change in coordinate frames

between the fuselage and the global coordinate system, determined by the angles b, 0,

and 4, and the changes between the joint coordinate frame on one body and the joint

coordinate frame on the connected body, determined by the angles al, Pi, a2, and /32.

Once the appropriate transformation matrices are formed, this information along with

the relevant position vectors and the calculated invariants are used in formation of the

mass matrix of each individual segment. This process is expanded in detail in Appendices

A and B. These individual mass matrices are then combined to form the system mass

matrix. Similarly, the partial and time derivatives of the various matrices, position vectors,

and invariants are used to form various time and partial derivatives of the system mass





-40L I
0 0.2 0.4 0.6 0.8 1

Figure 7-4: Time history of the angles as determined by GPOCS.





0.2 0.4 0.6 0.8 1

Figure 7-5: Time history of the angular rates as determined by GPOCS.

a ([C1] [E12] [T (a)] [ 23] [] (3)] [( 328a]T ,)
a ([C]] [a12] [T (a) 21]T [] [a ] (8)] [- ]T o') +

a ([c1] [al2] [T (a)] [al]T [023] [T (3)] [ ]T32] ) +
a ([c1] [a12] [T (a)] [a1]T [23] (3)] [ [3]T [N ] qf)
a ([C1 []a12] [T (a)] [Ea2lT []23] [T (8)] [Ea32]T [N3] q )

a ([C1] [a12] [T (a)] [a2]T [023] [T (/)] [a32]T [N3] q3)

a ([Cl] [a12] [T (a)] [a2jT [a23] [T (/3)] [32]T [N3] q3)
a ([C1] [a12] [T (a)] [a2]T [oa23] [T (3)] [a32]T [N3] q3)

= [c ] [oo ]

[c9] [o (] [ a] T ( a)] [c21 [.._ lT (3)] [p. _] g -}
[c[] [o C] [1 T (a)I [a a231 T ()] [a32 0 T [N3] 6q (3-91)

a(ciE [Cl] [0:l]

ac ]o [OCi [0:] (3-92)
a o(a [T (a)] [i2Bl3]
aor8)] [ [T\ () [23:32]


r 6RCB + (CII [Ci[ rBi_-Hi+

[Ci [0w] [al2l [T (a)] [a2 r rHICB2 + il ['] [122l [T (a)] [a2i]-rB2H2 +

[CII [O;] [l2l [T (a)i [a2ilr I [T (3)l [o rrH2-CB3 +

[CII [O;] [12l [T(a)] [a2ilr I [T( _]3)l [o .^ +

[Cil [O1] [al2l [T (a)] [ai il] i [T (3)J [_I 7 [N3] q) )6+

([Ci] [0wL] rBiHi+

[Ci] [0w ] [al2l [T (a)]l [ ]Tr rHICB2 + [CI] [0;] [al2l [T (ca)l [i] -rB2H2 +

[Ci] [Ow;] [al2l [T (a)l [a2il] [a23] [T (3)] [a32 rH2-CB3 +

and [M] = M (q), u is a vector of control torques, and Q = Q (q). Where Q o

are aerodynamic forces, Q are gravitational forces, and Q are propulsive forces. As
-grav --prop
written, the equations of motion are a system of second order differential equations. For

control purposes it is preferable to have a system of first order equations. Define

x q (5-13)

Thus the equations of motion can be written as:

S=f (x, ) (5-14)


f (5-15)
S [M]- Qr + [M]-l r [M]- + [M]- [K] q +[M]-

Now that the equations of motion are written in first order form, it is necessary to

explicitly define the control problem. An example of a goal of the control problem is to

drive the vehicle to a desired position using the control torques located at hinge axes

to change the configuration of the vehicle. This change of the configuration will then

influence the aerodynamics of the vehicle. Writing this in mathematical terms produces

the following objective function:

J = IXcB XD + YCB1 YD2 + ZCB1 ZD 2 (5-16)

This cost function seeks to minimize the distance between the final global position of the

MAV and some desired final global position. The optimal control is then the control that

minimizes Equation 5-16. The objective function is subject to physical constraints on the

system and must also satisfy the equations of motion. An example of physical constraints

is to limit the range of motion on the inboard and outboard segments to 90 deg with

respect to the inboard body( for the inboard segment the inboard body is the fuselage and

The optimal control problem most often sees implementation in the solution of

inverse kinematic problems. In these problem, the forces and torques that are required

for the system to achieve some final condition from a known initial condition are sought

as the solution of the optimal control problem. Work by Lo et. al posed an optimal

control problem to find the torques necessary for a human animation to follow prescribed

motions[791. In works by Huntington, a Gauss Pseudospectral method is used to determine

the sequence of maneuvers for a 4-spacecraft formation the results in minimum fuel

consumption. [80,81]

The optimal control problem for the articulated MAV introduced in C'! lpter 3 begins

with the conversion of the second-order equations of motion into a first-order form that is

standard for optimal control problems. The second-order equations of motion are written


[M] Q +[K] q = Q +Q +Q +u (5-10)
T grav -aero -prop




q= < > (5-11)



q 4

[Q=-[ }q+ ': (5-12)
.Tr[M] J
-l aq.

satellite.[51,52] Jain and Rodriguez extended a recursive multibody dynamics approach

to a study of the cassini satellite and a dynamic forward simulator.[531 Pradhan et. al,

applied a multibody formulation to the equations of motion for a space-platform-based

mobile manipulator system.[541 Attia uses a multibody formulation to model a pendulum

with a rotating end.[551

The benefit of using a multi-body dynamics formulation is that the individual parts

can be modeled in completely different v--iv. However, when combined into the global

system the separate modeling concerns are coupled. If a section is modeled as rigid and

another as flexible then the global system matrix will have coupling between the rigid

and flexible components. The sensitivity of the flexible-rigid interaction can also become

important, as demonstrated in the study by Dias and Pereira.[56]

[26] Livne, E. and Weisshaar, T.A., "Aeroelasticity of Non Conventional Airplane
Configurations-Past and Future", Journal of Aircraft, Vol. 40, No.6, 2003,

[27] Wlezien, R.W., Horner, G.C., McGowan, A.R., Padula, S.L., Scott, M.A., Silcox,
R.J., and Simpson, J.O., "The Aircraft Morphing P1iJ nti ,American Institute of
Aeronautics and Astronautics Paper No. AIAA-98-1927, 1998.

[28] Reich, G.W., Raveh, D.E., and Zink, P.S., "Application of Active-Aeroelastic-Wing
Technology to a Joined-Wing S. i .~ I i Journal of Aircraft, Vol 41, No.3, 2004,

[29] Zheng, L. and Ramaprian, B.R., "A Piezo-Electrically Actuated Wing for a Micro
Air \. !i. !, American Institute of Aeronautics and Astronautics Paper No.
AIAA2002-2302, 2002.

[30] Saunders, B., Eastep, F.E., and Forster, E., "Aerodynamic and Aeroelastic
C'! i I :teristics of Wings with Conformal Control Su I ...- ,Journal of Aircraft,
Vol. 40, No.l, 2003, pp.94-99.

[31] Abdulrahim, M., Garcia, H., and Lind, R., "Flight C'!i :'teristics of Shaping the
Membrane Wing of a Micro Air \. !i. I. Journal of Aircraft, Vol.42, No.l, 2005,

[32] Bisplinghoff, R.L. and Ashley, H., Principles of Ae ,.. .r1.i: .:;' New York, NY: Dover
Publications, Inc. 1975.

[33] Bisplinghoff, R.L., Ashley, H., and Halfman, R.L., A( ,.., I.1.:..i: /; Mineola, NY: Dover
Publications, Inc. 1996.

[34] Friedmann, P.P., "Renaissance of Aeroelasticity and Its Future", Journal of Aircraft,
Vol.36, No.l, 1999, pp.105-121.

[35] Kral, R. and Kreuzer, E., \!liil lhody Systems and Fluid-Structure Interactions
with Application to Floating Structures", Multibody System D;.inii. Vol. 8, 1999,

[36] Shabana, A.A., "Flexible Multibody Dynamics: Review of Past and Recent
D. !,pin, iii Multibody System D;.,iii-. Vol. 1, 1997, pp.189-222.

[37] Schiehlen, W., "Research Trends in Multibody System Dyri 1i1i1, Multibody System
D;.ir,. Vol. 18, 2007, pp.3-13.

[38] Yoo, W., Kim, K., Kim, H., and Sohn, J., "Developments of Multibody System
Dynamics: Computer Simulations and Exp. iin,. iii- Multibody System D;,iii,.
Vol. 18, 2007, pp.35-58.

bodies that comprise the system. Thus the system generalized coordinates can be written





q (7-38)





Thus the system kinetic energy takes on the form of

Ty f :-, _V] qsys (7-39)

Substituting Equation 7-39 into Equation 7-10 results in:

[M] q = : + Q (7-40)

where QE is a vector of generalized forces due to external loads on the system and Q

is a vector of "forces" due to the partial derivative of the kinetic energy with respect to

the generalized coordinates and the first time derivative of the partial derivative of the

kinetic energy with respect to the first time derivative of the generalized coordinates.

rH1-CB2,ma21,mnT (a)p a12,qp0 crq0 W,rsa2,sT (a), 21,vuta23,vwT (3) 32,yqf, N3,yz dV +

rB2-H2,va21,uvT(a) z a)12,xw I I r-H dV +

rB2H22,ra rsT (a)t a12,utO 1 0 w 2, (, 21,zyH-C2,z dV +

rB2-H2,ra21,rsT (a) ,t al2,ut Dvu ,vwa2,wT (a)y 21,zyB2H2,z dV +
1 0 rf

rB2-H2n2,n, npT (a),qp a12,rq0 esrO wstal2,tuT (a), a21,wva23,wxT (3), a32,zyrB2-CB3, dV +

rB2-H2,na21,npT(a) qp a12,rq W0 srO, wstal2,tuT (a), a21,wva23,wxT (3)xy a32,zy p,z dV +

rB2-H2,ma21,mnT (a), pn a2,qp wrq O0w ,rsal2,stT (a) t a21,vvu23,vwT ( -) ,w a f N3,yz dV +

rH2-CB3,ra32,rsT (3),ts 23,uta21,uvT (a),w 12,xw 0 rB1Bl-l HI dV +

rH2-CB3,nOa32,npT (3),qp a23,rqa21,rsT (a),t a12,ut0o 1,vuO 'w 2,wxT (a),xy Z21,zyrH-CB2,z dV +

rH2-CB3,na32,np qp 3rqa21,rsT (a),ts a12,ut0, ,vwl2,wxT (a),y 1,zyrB2-H2,z dV +

rH2-CB3,ja32,jkT (),mk a23,nma21,npT (a)q 12,rq ,sr 0 W,st 12,tuT (a),u a21,wva23,wxT (3),xy Ca32,zyrH2-C B3,z dV

rH2-CB3,ja32,jkT (3),mk a23,nm 21,npT (a),qp a12,rq0 ,sr lst12,tuT (a), 21,w 23,wT (),y 3,z op,z dV +

rH2-CB3,ia32,ijT (),kj a23,mka21,mnT (a)pn a12,qp CWrq 0,rs a2,stT (a)tu a21,vuC223,vwT (3),x a f N3,y dV +

a32,npT (3),ts a23,uta21,uvT (a), a2,xw0 a 0-1 rB Hl,z v op,r d +

a32,npT (),qp a23,rqa21,rsT (a),ts al2,utO,vuO1 ,,vwal2,wT (a),y 21,yrH-CB2,z op,n dV +

a32,npT (3),qp 23,rqa2l,rsT (a),ts a12,ut vu0 vw 2 T (a)y a21,zyrB2-H2,z op,n dV +

a32,jT (),mk a23,nma21,npT (a),qp a12,rq0 0,sr, Wsta2,tT (a),v a21,wv(23,wxT (f3),xy 32,zyrB2-CB3,z v Uop,j dV

a32,jkT (),m a23,nma21,npT (a)qp al12,rq c ,sr ,ast2l2,tuT (a) 21, a wvea23,wxT (/3),xy 32,zy op,jUop,z dV +

a32,ijT (3),j a23,mka21,mnT (a)p a2,qpO erq0o,rs1a2,stT (a),ta21,vua23,vwT (3),wxa Uop,jN3,yz dV +

qqa32,rsT (3),ts a23,uta21,uvT (a),wv a12,w- rBHl-HI,z N3,rq dV +

q9f,mO32,np (Q),qp a23,rqa21,rsT (a),ts al2,rqOw9vuo vwal12,wxT (a),y a21,zyrHl-CB2, N3,m dV +
q3 002 N
,ma32,npT (),qp a23,rqa21,rsT (a)ts a2,rq1 01,vu vwl2,wxT(a),xy 21,zyrB2-H2,z N,m dV +

9,ia32,jkT (),mk a23,nma21,npT (a),p a,rq0 ,srO,sta12,tuT (a),uv21,wva23,wxT (3),xy 32,zyrB2-C'B3,z / N3, dV +

Vq,i"32,jkT (),mk a23,nm 21,npT ( 12,rqp ,srOe ,sta12,tuT (a),uv 21,wva23,wT (3) 32,zy N3,ijUop,z dV +
qf a 32,ijT ( ,j 23,mk 21,n ,qp O ,rqa ,s 2,st 21 u 23,v T ,v 3,ij 3,y
qJ3, h32,ijT (),kj d23,mkd21,mnT (a) p ai12,qp 011l" ~v",sal2,stT (a) alvua23,vuua N3,iN3,y JX (B 8)


This research was concerned with the emergence of Micro Air Vehicles and the various

new opportunities that they provide. One of the most compelling is the ability to use

morphing on MAVs. The small structure and light weight materials allows for morphing

strategies that can completely alter the configuration of MAVs. These configuration

changes can then be used to exploit flight regimes and maneuvers that are alien to

conventional aircraft. However, with these new opportunities come new obstacles that

must be overcome. The chief among these obstacles is the lack of a way to characterize the

performance of MAVs based on conventional aircraft modeling. New modeling techniques

are required that allow for the, in some cases, complete reconfiguration of the vehicle

because in the flight regime that MAVs typically operate in the changes have great impact

on the performance. In an attempt to create new modeling techniques, this study outlined

the derivation of equations of motion based on multi-body dynamics that allows for the

large configuration changes possible with morphing. These equations of motion also

incorporate the coupling of rigid body and flexible body dynamics into a single set of

equations. Also presented was the ability to take the equations of motion as derived and

cast them into a form amenable for an optimal control problem. The numerical example

demonstrated the implementation of the multi-body equation of motion formulation on a

model that simulates the configuration changes possible on MAVs. The three-bar example

shows that the methods derived are applicable to traditional dynamic systems as well and

through the optimal control problem, allows for new control strategies.

However, there are challenges introduce by this research that need to be overcome

as the techniques mature. Primarily among these are the numerical difficulties associated

with the solution of the Optimal Control Problem. With systems that incorporate rigid

and flexible modes, there are a large number of degrees of freedom. It is well known that

as the degrees of freedom increase, the solution time for direct methods also increases. As,

formulation, there is an inertial system and a body-fixed coordinate located at the center

of mass of each respective body in the system. The first body that is attached to the

ground is the "central" body in this system. This body is only allowed to rotate about the

axis of the revolute joint that connects this body to the ground. Thus, the relationship

between the first body and the inertial frame is represented by:

e = [C] b (7-84)

where e, is a vector written in terms of the inertial coordinate frame and b{ is a vector

written in terms of the body-fixed coordinate frame associated withe body 1. Here,

cos 1 -sin Q1 0

[C] = sin 1 cos 1 0 (7-85)

0 0 1

In the alternative formulation illustrated in C'! plter 3 showed joint frames that were used

as the basis for relating the body-fixed coordinate system of one body to the body-fixed

coordinate system of an adjoining body. For this example these joint frames are not

necessary as the body-fixed coordinate frames are aligned with the joint axis of the

revolute joints for all bodies. Thus the relationship between the first body and the second

body can be stated as:

r 1 [T(02)]rB (7-86)

1 0 0

[T(40)] 0 cos 2 -sin 2 (7 87)

0 sin 02 cos 2

and rB2 is a vector written in terms of the body-fixed coordinate system associated with

body 2 and rB1 is the same vector written in terms of the body-fixed coordinate system

associated with body 1. Similarly, the relationship between body 2 and body 3 can be







6-4 Position of MAV from integrated optimal control. ............... 72

6-5 Morphing angles from integrated optimal control. ............. .. .. 72

7-1 Norm of the constraint vector ............... ........ .. 103

7-2 Norm of the constraint vector ............... ........ .. 103

7-3 Three bar manipulator. ............... ........... 104

7-4 Time history of the angles as determined by GPOCS. . . 105

7-5 Time history of the angular rates as determined by GPOCS. . .... 105

7-6 Time history of the control torques. ............... ... .. 106

7-7 Time history of the angles as determined by integration. . . 106

7-8 Time history of the angular rates as determined by integration. ........ 107


M orphing M AV. .................. ... ........... 16

Articulated M AV. .................. .... .......... 51

Coordinate frames of the fuselage, inboard wing segment, and outboard wing
segm ent. . . . . . . ........... .. .. 51

Joint coordinate frames between fuselage and inboard wing segment. . 52

Joint coordinate frames between inboard and outboard wing segments. . 52

Geometric model of gull-wing MAV. .................. ..... 70

Position of MAV when using the optimal control solution as determined by GPOCS. 71

Morphing angles when using the optimal control solution as determined by GPOCS. 71

Mathematically, Q takes the form:

*T O[M]
*T 9[M]
Q LM < q q2q (7-41)


The three bodies are not completely free to move in the inertial frame. Through the use

of the revolute joints, the bodies are constrained to move relative to each other. Since in

the standard derivation of the equations of motion the bodies are not defined relative to

each other, another term appears in Lagrange's equations to account for the constraints.

Lagrange's equations become:

d (9T -9T )
dt ) +[% A (742)

where [Cq] is the matrix that results from taking the derivative of the constraint matrix

with respect to the generalized coordinates also known as the system Jacobian matrix.

A revolute joint only allows for a single degree of freedom between two bodies. Thus,

each revolute joint has five constraint equations associated with it. This system has three

revolute joints, as illustrated in Fig. 7-3; one between the ground and body 1 that allows

rotation about the bi axis of the body-fixed coordinate system associated with body 1,

one between body 1 and body 2 that allows rotation about the b2 axis of the body-fixed

coordinate system associated with body 2, and one between body 2 and body 3 that

allows rotation about the b3 axis of the body-fixed coordinate system associated with body


The constraint matrix is developed by setting the constraint equations that govern the

motion of the bodies equal to zero. Thus,for the system:

[=] 0 (7-43)

more greatly impacts the vehicle due to slight changes in the configuration of the MAV

having a large effect on its performance as a result of the high sensitivity that exists

between the MAV and the regime it operates in. Typically, motions induced by the flow

regime are on the order of the dimensions of the MAVs, themselves. Thus, changes in

the flow regime, an example of which is gust loading, can have radical effects on the

configuration of the MAV. In a synergistic manner, these configurations changes then

lead to differing responses of the overall vehicle in the same flight regime. Thus the

challenge to researchers in the study of MAVs is how to incorporate the flexibility of the

composite structures, the nature of the low Reynolds' number flow, and the interaction

of the vehicle and the regime into a framework that allows for successful modeling of the

MAV. Once a model exists that addresses the challenges of the flexibility of the structure

and the interaction of the structure and the medium through which it passes, it then

becomes possible to use this model as a predictive tool in the development of further

MAV technology. For example, a proper MAV model can be adjusted to reflect design

changes and then simulate the adjusted model to determine the performance as a result

of the design change or a control strategy can be applied to the model with the goal of

determining the effectiveness of the control strategy. Also control strategies can be used

on the model A proper investigation into the modeling of MAVs must first begin with an

understanding of the various fundamental topics that are concern MAVs.

1.2 Framework

The work presented in this work shows the evolution of modeling and control of a

morphing micro air vehicle. Micro air vehicles are a relatively new technology. There has

not been a comprehensive study on the modeling of a MAV that allows for a large degree

of articulation to allow for morphing. It is the goal of this research to bridge this gap

in the lexicon of vehicle dynamics to allow for a complete model of a highly articulated

MAV to enable future control studies of such vehicles. MAVs present a good platform

for a study of dynamics and control strategies due to the ri ir, interesting challenges

wingbeat kinematics. The hummingbird wing is modeled using simplified rigid structures.

This model is then implemented in SIMULINK in an attempt to reproduce wing tip

trajectories. A mechanical model is then constructed based on the simplified rigid model

that is used to verify the trajectories predicted by the simulation. This modeling attempt

neglects flexibility effects and does not demonstrate the mimicking of lift and thrust

characteristics of actual humming birds that would be necessary for a flapping wing

MAV. In a study by Schenato et al., the controllability of a flapping flight MAV was

investigated.[17 The study used simplified dynamics of an insect flight. Flexibility effects

were not taken into account.

The University of Florida has developed several successful MAV platforms. There are

small MAVs that have a maximal dimension of 6 in. and there are larger MAVs that have

up to 24 in. wingspan. Both the larger and smaller MAVs developed by the University

of Florida rely on flexible membrane wings as opposed to rigid wings. The development

of flexible wings at the University of Florida was based on biological inspiration. Birds

and bats that also fly in the same low Reynolds' number regime as MAVs have thin and

undercambered wings.12,13]

The larger UF MAVs generally consist of under-cambered wings composed of a

plastic membrane over composite battens.[18 As a result of this construction, ailerons are

not easily mounted on the wings. Thus control authority on these vehicles are provided

by rudder and elevator. The wings on these MAVs are highly flexible. This flexibility

results in a passive mechanism known as adaptive washout, which changes the angle of

attack along the wing and reduces sensitivity to external disturbances. [19 The adaptive

washout demonstrated by the flexible UF MAV wing is produced by the extension of the

membrane and twisting of the framework.[12] In nature, adaptive washout can be seen

when the wingtips of birds flair to accommodate sudden changes in airspeed. Experiments

performed on the 6 in. UF MAVs show that the primary mode of deformation occurs at

around 140Hz.[20]

matrix that are then utilized in the evaluation of the equations of motion as derived from

Lagrange's equations (see ('!Ci pter 3 for further detail).

For the purposes of calculating aerodynamic forces on this example, an AVL input

file is created that mimics the geometry of the example. However, in order for this

configuration to produce aerodynamic lift forces, the wing must have camber. To this

end, in the .avl input file, the wing segments are made to have a camber profile. AVL also

requires input parameters that define the runcase. These inputs are the velocity of the

vehicle, the angle of attack of the vehicle, and the angle of sideslip of the vehicle. These

quantities are derived from the states of the generalized coordinates as follows:

velocity : /U2 + v2 + w2

angle of attack : tan-1 (U)

angle of sideslip : sin1 (v)


v = [C]T YCB1


V = u2 + v2 + w2

The results of the AVL execution are force coefficients, based on the entire aircraft

configuration, along the axes of a coordinate system fixed in the fuselage. These

coefficients are then transformed into forces directed along the X, Y, and Z-axes of

the fuselage's coordinate system. AVL also produces moment coefficients but for the

purposes of this numerical example, the moments produced from the moment coefficients

are ignored.

At this point, the geometry of the example is fully defined. From this geometry,

the system matrices can be formed using the derivations in ('!i lpter 3, Appendix A,

and Appendix B. These system matrices are the mass matrix, the first time derivative

of the mass matrix, and the partial derivative of the mass matrix with respect to each

of a part rather than the overall motion of the system. Here the aeroelastic properties

can be investigated, but generally at the loss of the overall motions of the structure as

a whole. Using multi-body dynamics, the system as a whole is studied but with great

care as to the contribution from the individual parts that comprise the system. Thus,

the level of detail for each body can be chosen to determine the deformation of the

body while still connecting the body to the overall motions of the system as a whole.

However, there are a number of v--v a multi-body system can be analyzed. The two most

common techniques are a Newton-Euler formulation and a Lagrangian formation. Work by

Shabana, Schiehlen, and Yoo demonstrate some of the other popular choices for composing

the equations of motion.[36--38] The Newton-Euler formulation has the advantage of

producing equations of motion that are computationally efficient to solve. The Lagrangian

approach has the advantage of simpler derivation of the equations of motion. In work

done by Boutin and Misra, a low Earth orbit satellite was modeled as a collection of truss

structures using a multi-body dynamics approach.[39]

A paper by de Jalon presents a study of the different choices of natural coordinates

that have been made in the last 25 years of dynamic studies.[401 Similarly, in work by

Leger the choice of coordinate systems determines the speed of subsequent dynamic

simulations.[41] Included in the coordinate system effecting simulation time is the nature

of the time integration algorithms.[42--46] Kielau and Maiber presented a study in

which various coordinate systems and equations formulations are used on nonholonomic

systems.[47 Aoustin and Fromalsky presented work in determining feedforward simulation

of a flexible, two-link manipulator.[481

Multibody dynamics can be used to model a number of complex systems. Junkins

presented a study of orbital dynamics using a nonlinear multibody dynamics approach.[491

Rives et. al modeled the descent phase of re-entry vehicle.[50] Boutin and Misra modeled

a struss structure using multibody dynamics.[39] Tadikonda used a multibody approach

in the modeling of translation motion between two flexible bodies and a geostationary

-3 0 I
0 0.5 1 1.5 2 2.5 3

Figure 6-4: Position of MAV from integrated optimal control.


0.5 1 1.5 2 2.5

Figure 6-5: Morphing angles from integrated optimal control.

a32,jkT (3), a23,nm 21,npT (a)p al2,~rq 01 a2,tT (a)t 12 L 2v21 ,wva23,wxT (3),y a32,zyrH2-CB3,z f dV +

a32,jkT (P), a23,nma21,npT (a) a2,rq0 o ,- a2,stT(a),t 12 ,e21, va1 23,wvT (3),y 32,zy / op,jop,z dV +

a32,ijT (0),kj a23,mka21,mnT (a),pn al2,qp rqa12,rsT (a) ,s 12 1 tu~21,vu~a23,vu T (3) x z Uop,i3,yz dV +

q 3 0 12 21
9 ,m_32,npT (),qp a23,rqa21,rsT (a),ts a12,ut0 W ,vu~12,avwT (a),w 12 ,21,zyf1-CB2,z N3,nm +

9f,mJ 32,npT (/),p a23,rqa21,rsT (a),ts a12,utC ,~ a12,avT (a),wx 12 -2?1 21,zyB z m dV
q 32,jk ,mk ,23,nm21,npT (),p a12,rq D ,s.ro12,stT(a),tu 121uv21,wva23,wxT (),xy 32,zyfH2-CB3,z N3,ji d +
fua2laua23,axT(8),y a32,y0H2 qB3,

q ,iO32,jkT ( ,k 23, 21,npT (a) a,rq ~s 12stT (a),tu 12 Ce21 a 23,wxT () y 32,zy N3,jiop,z dV +

q3,ha32,ijT (3),kj 23,mk21I,mnT (a),pn a12,qp0 wrq12,rsT (a),st 12 1tu21, a23,wT 3) a N3,ihN I (B-10)

For the remaining terms, replace o01, with o0), and a1I where appropriate.

S[BBT DD2dV = BBTDD2 dV (B-11)


Examining a representative term:

/BB DD, dV =

rBl-Hl,ro0's,ral2,stT (a),tua21,vuat23,vwT(3) ? a32,zyfH2-CB3,z dV +

rBl-HI,r w,sral12,stT (a)t a21,vua23,vwT (S3), 23 32, / opz

rB1 -H,q ,rqa12,rsT (a),s a21,uta23,uvT (3), 23. a N3,yz dV+

rH1-CB2,na21,npT(a)qp al2,rqO wsral2,stT(a)tu a21,vua23,vwT (3),wx 23 ,32,zyrH2-CB3, dV
rH1-CB2,na21,npT(a),p al2,rq0 w ,sra12,sT (a),tu a21,vu~23,uvT (3) ,wx 23 zy op,z V +

rH1 CB2,ma21,mnT (a),pn a12,qp 0w ,rqa12,rsT (a),st a21,tu~23u (3) 23W2z 3 a z N3,y dV +

rB2-H2,na21,np Ta) qa12,rq0 wlsral2,stT(a),uta21,vua 23, 23 ? 3 32,zyrH2 CB3,z dV

B2 H2,na21,npT (a),qp a 2,rq0 sa12,stT (a),tu a21,vua23,vwT (),wx 23 32,zy Uop,z dV +

9 V
rB2H2-,mB21,i T 2 12,qp 2,qa12,rsT (a)'2 a t21 .a23,uvT () 2 a,~ z N3,yz dV +

rH2 CB3,ja32,jkT ()mk a23,nma21,npT (a),qp l12,rq W ,sr 12,stT (a),u a21 ,vu23,vwT (w3),wx 23, ,a32,zyH2 Bo3,z d V

rH2 CB3,ja32,jkT (63)mk a23 ,nma 21,npT (a) ,qp a 12,rq 0 ,sr~a12,stT (a),tu a21,vua23,vwu T () ,wx 23 a32,z ?op dV +
013 2332 y dV
rH2 CB3,ia32,ijT (3)kj a23,mka21,mn.T (a),pn a12,qp OW^ qal2,rsT (a),st 21,uta23,uvT (3), 23x32 ., / N3, dV +
0-1 V
a32,jkT (P)k a23,nma21,npT (a) p a12,rq osral2,stT (a),t 21,vu~a23,vwT (3),ax 23- a32, yrH2-CB3,z uop,j dV +
,q3 ,z H CB f oV d

Figure 3-1: Articulated MAV.

Figure 3-2: Coordinate frames of the fuselage, inboard wing segment, and outboard wing

In work performed by Lian et al. numerical simulations were performed on a 6in.

membrane wing model of the UF MAV configuration.[21,221 A full CFD/CSD interaction

simulation was run. The wing was modeled as a membrane using a Mooney-Rivlin

formulation. While the simulation did have idealized flexibility effects, there was not an

accounting of the overall motion of the wing. Only the deformations of the wing were used

in the computation of the effect of the wing shape on the aerodynamics for a particular

run. Numerically, the vibration of the wing was found to be on the order of 102Hz.[22,23]

The University of Colorado has developed a MAV known as C \! AV.141 The C \! AV,

much like the MAVs developed at the University of Florida has a low aspect ratio, flexible

membrane wing that is supported by carbon fiber ribs and battens.

2.3 Morphing Aircraft

Morphing is defined as a change in the shape or structure of an aircraft surface

in a continuous manner. Morphing has been used in aircraft since the powered flight's

inception. The Wright Flyer morphed the wings by pulling on cables in order to control

the flight.[24]

The goal of morphing a wing is to adapt to multiple flight regimes to obtain better

performance. The morphing of the wing can be accomplished by camber change, wing

twist, wing sweep change, and wing span change.[25] The goal of camber change is to

obtain the desired lift on a wing without the use of conventional control surfaces. Wing

twist is used in obtaining a low-drag and high-lift configurations. Wing sweep is used

to decrease drag as speed is increased. In wing span morphing, the wing aspect ratio

and wing area are altered to effect drag of the wing and consequently the range of the

vehicle. These are just simple examples of the type of performance changes that can be

accomplished by morphing. Camber change, wing twist, wing sweep change, and wing

span change can also be used to obtain other performance characteristics. An example of

this is the use of wing sweep in tailless configurations to increase stability.

Figure 3-3: Joint coordinate frames between fuselage and inboard wing segment.

Figure 3-4: Joint coordinate frames between inboard and outboard wing segments.

very much. While this assumption is valid for steady state applications, they are less than

optimal for any geometry undergoing transient behavior.

In addition, AVL is only applicable for low Mach number flow. An AVL analysis can

be completed without the use of a fuselage in the geometry definition. When a fuselage is

included, it is treated as a bluff body and as a result the forces calculated that act on the

fuselage are not quite as accurate as a more robust analysis tool would generate.

These short-comings indicate that the use of AVL is not optimal for the generation

of aerodynamic forces in the simulation of an unsteady process of aeroelasticity of a

re-configurable MAV. Instead, more robust measures such as the use of a fully 3D

Navier-Stokes CFD software would be more appropriate. However, the use of such

a software on a re-configurable MAV is a research topic in itself. Also, the current

state-of-the-art in computational fluid dynamics is very computationally expensive

to generate a solution for the forces acting over given geometry. First an appropriate

discretization of the geometry and the fluid surrounding the geometry must be generated.

Once the discretization has taken place, the solution of the aerodynamic forces is an

iterative process and may take a large amount of time to converge to a single solution.

This problem is compounded when used in an aeroelastic simulation as the configuration

of the vehicle changes at every time step and the computational cost needed to generate

the aerodynamic forces must also be paid at every time step. In contrast, AVL takes

a matter of seconds to converge to a solution for a given geometry input, due to the

simplicity of the lattice vortex method. Thus as a computationally inexpensive source

of aerodynamic forces, AVL is used even though proper application of the process would

choose a more robust method.

In a study by Bae et al. variable wing span morphing was tested on a model of a long

range cruise missile. Both the static aerodynamics and the aeroelastic characteristics of a

variable span wing are investigated. The aeroelastic analysis uses a panel code together

with a wingbox structural model. While this investigation does include the effects of

flexibility on the aerodynamics and vice versa, it does not dynamically simulate the

varying of the span in flight.

Several initiatives have been started by organizations to investigate morphing.

According to the program at DARPA, a morphing vehicle will change its shape substantially

to adapt to the mission environment, provide superior system capability not possible

without reconfiguration, and use integrated design of materials, distributed actuators,

effectors and mechanisms to reconfigure in flight.[26] The Aircraft Morphing Program

developed by the Airframe Systems Program Office of the NASA Office of Aeronautics

and Space Transportation Technology is designed to develop smart devices for airframe

applications using active component technologies.[27]

In research conducted by Reich, experiments were performed in which multiple

control surfaces were actuated in order to effect a change in the camber of the wing of a

sensor craft and thereby alter the orientation of the built in antennae.[281 Hall et al. have

conducted a numerical simulation to determine the effectiveness of variable sweep on the

lift-to-drag ratio of the proposed second generation vehicle for the University of Colorado's


In a study performed by Zheng and Ramaprian, a piezoelectric actuator was used

to change the lift coefficient on a MAV without changing the angle of attack of the

wing.[29] In the first part of the study, an interactive structural and aeroelastic analysis is

performed to determine the relationship between the DC voltage applied, the freestream

velocity and the angle of attack of the wing. To this end, the wing cross section was

modeled as a composite beam. A 2D airfoil viscous flow solver is used to solve for the flow

over the airfoil. The deflections from the piezoelectric actuator are then superimposed

Thus using the skew symmetric matrices to represent the vector cross product and writing
all components in terms of the body-fixed coordinate frame associated with body 1:

[01] -B-H + [1 1] [T ( 2)] rH-CB2 + [T (2)1 [2 rH-CB2

+ [O' ] [T (42)] U + [T (02)] [12] u (7 104)

The skew symmetric matrices can be factored to yield:

v = di ([",] rB1-H1 + [0 ] [T ()1 rHI-CB2 + [w1] [T (2)] Uop)

'2 ([T (2)] [1 2] rH1-CB2 + [T (02)] [1 _2] P) (7 105)

B L = '[0] rBiHi + [B-H1 [T (0H2)1 rH1-CB2 + [01,] [T (02)] Uo (7 106)

B2 [T (2)] [w122J rHI-CB2 + [T (_2)] [1W2] U (7 107)

[B] [ B B2 ] (7-108)

Thus the velocity of a random point on body 2 is written as:

v, = [B] (7-109)
Inserting Equation 7-109 into the contribution to the system's kinetic energy from body 2

T = 1 '2 } 2Pp [B] [B] dV (7-110)

Thus, the contribution to the kinetic energy from body 2 assumes the familiar quadratic
T =-4[M] (7-111)

If the roots of Equation 7-83 lie in the left-half plane, the solution is guaranteed to

converge to zero. Thus, with the proper choice of a and 3 in Equation 7-83 the constraint

vector and the first and second time derivatives of the constraint vector can be made to

converge to zero.[86--88] Standard application of the Baumgarte stabilization method

assumes 3 = a, with 15 < a < 50.[89] Figure 7-2 shows the comparison between the

norm of the constraint vector without stabilization and with stabilization. As illustrated

in the figure, without constraint stabilization the error in the constraint vector grows

exponentially with time. Stabilization helps to reduce the errors accumulated and the

rate at which they accumulate. However, as illustrated the simultaneous solution of the

Lagrange multipliers and the accelerations that is necessary when using the standard

formulation of the equations of motion still leads to violations of the constraints on the

system. For this reason, forming the equations of motion without the need to couple

constraints equations is beneficial in the overall simulation of the system.

7.3 Alternative Formulation

As illustrated in the above example, the standard formulation of the equations of

motion of a multibody system leads to a complex system. There are a large number

of generalized coordinates that increases with the number of bodies in the system.

Also, the correct integration of the equations of motion is difficult due to the difficulties

associated with the solution of the Lagrange multipliers. Even with constraint stabilization

techniques, the process is susceptible to error. Application of the formulation illustrated

in C'! lpter 3 allows for the elimination of some of these difficulties. The number of

generalized coordinates is much less with this alternative formulation. Since all bodies in

the system are referenced to a "central" body, only the full coordinates associated with

this body are necessary. Those coordinates are then augmented with coordinates that

relate the differences in orientation between subsequent bodies and the "central" body.

The derivation of the equations of motion, using the alternative formulation, for the

three bar system begins with the establishment of coordinate frames. As in the standard

[65] Vadali, S.R. and Sharma, R., "Optimal Finite-Time Feedback Controllers for
Nonlinear Systems with Terminal Constraints", Journal of Guidance, Control, and
D;i,.i, Vol.29, No.4, 2006, pp.921-928.

[66] Betts, J.T., "Survey of Numerical Methods for Trajectory Optimization", Journal of
Guidance, Control, and D.,i,.. Vol.21, No.2, 1998, pp.193-207.

[67] Hull, D.G., "Conversion of Optimal Control Problems into Parameter Optimization
PiN1.!. in- Journal of Guidance, Control, and,,. Vol.20, No.l, 1997,

[68] Thompson, R.C., Junkins, J.L., and Vadali, S.R., \N. i-Minimum Time Open-Loop
Slewing of Flexible \; I.! I. Journal of Guidance, Control, and D;;,i.. Vol.12,
No.l, 1989, pp.82-88.

[69] Warner, M.S., and Hodges, D.H., "Treatment of Control Constraints in Finite
Element Solution of Optimal Control Pi.I!. in- Journal of Guidance, Control, and
D;i,.. .- Vol.22, No.2, 1999, pp.358-360.

[70] Bass, M., von Stryk, O., Bulirsch, R., and Schmidt, G., "Towards Hybrid Optimal
Control", at-Automatisierungstechnik, Vol.48, No.9, pp.448-459, 2000.

[71] von Stryk, O., and Bulirsch, R., "Direct and Indirect Methods for Trajectory
Optimization", Annals of Operations Research, Vol.37, pp.357-373, 1992.

[72] Hu, G.S., Ong, C.J., and Teo, C.L., "Direct Collocation and Nonlinear Programming
for Optimal Control Problem Using an Enhanced Transcribing Scheme", Proceeding
of the 1999 IEEE International Symposium on Computer Aided Control System
Design, Kohala Coast-Island of Hawaii, Hawaii, USA, 1999.

[73] Frazzoli, E. and Bulle, F., "On Quantization and Optimal Control of Dynamcal
Systems with Symmetries", Proceeding of the l1st IEEE Conference on Decision and
Control, Las Vegas, NV, 2002, pp.817-823.

[74] Huntington, G., Benson, D., and Rao, A., "A Comparison of Accuracy and
Computational Efficiency of Three Pseudospectral Methods" 2007 AIAA Navi-
gation and Control Converence Hilton Head, SC, 2007, AIAA Paper2007-6405.

[75] Huntington, G. and Rao, A., "A Comparison of Local and Global Orthogonal
Collocation Methods for Solving Optimal Control Pi"!.. in 2007 American
Control Conference, New York, NY, 2007.

[76] Williams, P., "A Gauss-Lobatto Quadrature Method for Solving Optimal Control
PiI. in- Australian and New Zealand Industrial and Applied Mathematics
Journal, Vol.47, 2006, pp.C101-C115.

[77] Fahroo, F. and Ross, I.M., "Costate Estimation by a Legendre Pseudospectral
Method", Journal of Guidance, Control, and D;in, Vol.24, No.2, 2001,



-50 L

-1-- -A--C- C1
-100 C
-150 -- A
0 0.2 0.4 0.6 0.8 1

Figure 7-6: Time history of the control torques.






0.2 0.4 0.6 0.8 1

Figure 7-7: Time history of the angles as determined by integration.

Taking the partial derivative of Equation 7-59 with respect to the generalized coordinates

S -[I] [Ci] [o:)] U [Ci] [0:)1,1] U11 [CI [0:1)] U1 [I] [C02] [1w ] -

[C2] [1 2 u1 [021 [[w2] 2H1 [0] 0 0 0
Now looking at the constraints the second revolute joint imposes on the angular degrees of
freedom between body and body 2

5 = b- b- 2 0 (7-61)

which can be written in terms of the inertial coordinate system as:

)5 -= h [C]T [C2] b2 0 (7-62)

where b1 is the unit basis vector of the body-fixed coordinate frame of body 1 that is
aligned with the axis of rotation of the revolute joint between body 1 and body 2 and b2 is

a unit basis vector of the body-fixed coordinate frame of body 2 that is perpendicular to
the body-fixed unit basis vector in body 2 that is aligned with the axis of rotation of the
revolute joint between body 1 and body 2. The contribution to the Jacobian matrix from
Equation 7-62 is:

9 0 b1 [o C) [C]]T [[2C 62 1 0 iT 2
IJT [0:C)l [71]T [(172 0T hIT [r71]T [C 2] [1:2 h2
]T [o [C ] [C2] b 0 bT b C2
T 2 1 [ 12 (7-63)
1 02 a2 IT 1 2 0
b1T [l]T [2] L J 2 [] ~i]T [C21 [iw ^] h2 -T

0 0 0

Similarly, the second constraint on rotational degrees of freedom imposed by the revolute
joint between body 1 and body 2 is illustrated in:

)6 = b b2 0 (7-64)

Full Text








Iwouldliketothankallofthememberofmycommitteeforthesupport,direction,andconstructivecriticismtheyprovidedinthecompletionofthisresearch.Also,IthankmycolleaguesintheFlightControlLabfortheirsacricesmadebysettingasidetheirownresearchtoaidmeinmywork. 4


page ACKNOWLEDGMENTS ................................. 4 LISTOFFIGURES .................................... 7 ABSTRACT ........................................ 8 CHAPTER 1INTRODUCTION .................................. 9 1.1Motivation .................................... 9 1.2Framework .................................... 10 1.3ProblemStatement ............................... 11 1.4ResearchPlan .................................. 12 1.5Contribution ................................... 15 2LITERATUREREVIEW .............................. 17 2.1LowReynolds'NumberFlow .......................... 17 2.2MicroAirVehicles ............................... 17 2.3MorphingAircraft ............................... 20 2.4Modeling ..................................... 23 3DERIVATIONOFEQUATIONSOFMOTION .................. 26 3.1MinimalSetofGeneralizedCoordinates ................... 26 3.2GeneralizedForces ............................... 43 3.3Summary .................................... 50 4DETERMINATIONOFAEROELASTICFORCES ................ 53 4.1AthenaVortexLattice ............................. 53 4.2Short-Comings ................................. 54 5OPTIMALCONTROLSTRATEGY ........................ 56 6NUMERICALEXAMPLE{ARTICULATEDMICROAIRVEHICLE ...... 62 6.1ModelParameters ................................ 62 6.2Results ...................................... 67 7NUMERICALEXAMPLE{THREEBARMANIPULATOR ........... 73 7.1ModelDescription ............................... 73 7.2StandardFormulation ............................. 74 7.3AlternativeFormulation ............................ 89 5


........................... 100 8CONCLUSIONS ................................... 108 APPENDIX AMASSMATRIXFORINBOARDWINGSEGMENT ............... 111 BMASSMATRIXFOROUTBOARDWINGSEGMENT ............. 116 REFERENCES ....................................... 125 BIOGRAPHICALSKETCH ................................ 133 6


Figure page 1-1MorphingMAV. ................................... 16 3-1ArticulatedMAV. ................................... 51 3-2Coordinateframesofthefuselage,inboardwingsegment,andoutboardwingsegment. ........................................ 51 3-3Jointcoordinateframesbetweenfuselageandinboardwingsegment. ...... 52 3-4Jointcoordinateframesbetweeninboardandoutboardwingsegments. ..... 52 6-1Geometricmodelofgull-wingMAV. ........................ 70 6-2PositionofMAVwhenusingtheoptimalcontrolsolutionasdeterminedbyGPOCS. 71 6-3MorphingangleswhenusingtheoptimalcontrolsolutionasdeterminedbyGPOCS. 71 6-4PositionofMAVfromintegratedoptimalcontrol. ................. 72 6-5Morphinganglesfromintegratedoptimalcontrol. ................. 72 7-1Normoftheconstraintvector ............................ 103 7-2Normoftheconstraintvector ............................ 103 7-3Threebarmanipulator. ................................ 104 7-4TimehistoryoftheanglesasdeterminedbyGPOCS. ............... 105 7-5TimehistoryoftheangularratesasdeterminedbyGPOCS. ........... 105 7-6Timehistoryofthecontroltorques. ......................... 106 7-7Timehistoryoftheanglesasdeterminedbyintegration. ............. 106 7-8Timehistoryoftheangularratesasdeterminedbyintegration. ......... 107 7


Thisworkpresentsamethodologytodetermineofaccuratecomputermodelsofarticulatedmicroairvehicles.Priortothiswork,muchoftheinnovationindesignofamicroairvehiclewasbasedonempiricaldata.Thisemphasisonempiricaldatanecessitatesatimeintensiveiterativeprocesswhereadesignisbuiltandtested,thenmodicationsaremadeandtheprocessstartsanew.Thisworklookstoalleviatesomeoftheloadofthisiterationbyproducingaccuratemodelsthatcanbeinsertedintotheprocessatgreatlyreducedcostintermsoftime.ThemodelingisaccomplishedbytheapplicationofLagrange'sequationstocreateasetofequationsthatincorporateboththerigidandexibledegreesoffreedomassociatedwiththematerialsthatMicroAirVehicles(MAVs)employ.Theseequationsofmotionalsoincludethegeneralizedforcesthatresultfromthebodyforcesduetogravityandexternalforcesactingasaresultoftheinteractionbetweenthestructureandtheairthroughwhichitpasses.Theseequationsareformulatedinamannerthataccomplishestheabovewhileallowingfortherecongurationofthesystem.Theequationsofmotionarethentransformedintoaformamenableforuseinthesolutionofanoptimalcontrolproblem. 8


OneofthechallengesinutilizingtheMAVtechnologyisthepropermodelingofvehicles.ThechallengeinmodelingofMAVsisexempliedbythesmallsizesofthevehiclesandtheightregimeinwhichtheyoperate.Whereconventionalaircrafthavebeenstudiedingreatdetailandawealthofknowledgeexistsontheirightregime,thereismuchlessinformationforMAVs.TheinteractionoftheMAVcontrolsurfacesandthelowReynolds'owregimeinwhichtheyoperateisofgreatinterest.ConstructionmaterialsofconventionalaircraftarenotviableforMAVs.Thesmallsizedictatestheuseoflightermaterialbuttheextremeeectthattheightregimehasonthevehiclecallsforbuildingmaterialtoalsobemorefunctionalthanthematerialsinconventionalaircraft.Asaresult,compositestructuresarecommonplaceinMAVconstruction.ThisuseofcompositematerialfurtherstressestheneedfornewmodelingtechniqueswheninvestigatingMAVdynamics. However,anexcitingnewopportunityisprovidedbytheMAVtechnologies.ThesensitivityoftheMAVstotheirenvironmentsallowsforagreaterbenetfromdierentcontrolstrategiesthanconventionalaircraftwouldallow.Forexample,whilemorphingexistsonconventionalaircraft,theextentoftheeectofmorphingoftenpalesincomparisontothecomplexityofthemorphingsystem.OnMAVs,however,morphing 9






1-1 .Asopposedtomoretraditionalaircraft,theMAVpresentedingure 1-1 haswingsthatarearticulatedtoallowfordierentanglesbetweentheinboardwingandthefuselageandalsobetweentheoutboardwingandtheinboardwing.Ratherthanattempttomodelthecompletevehicleinasingleformulation,amultibodyapproach,asdemonstratedbyShabana[ Asidefromthemultibodyapproachtomodeling,thereisanothermatterthataddstothecomplexityoftheformulationoftheequationsofmotion.Thevehiclewillbemodeledwithbothrigidbodyandexibledegreesoffreedom.ForMAVs,thefuselagecanbemodeledasarigidbodyasthereisnotmuchdeformationofthefuselage.Thewingstructure,howeverishighlyexible.Theaeroelasticloadsonthewingscoupledwiththeirlightweightandlowstiness,incomparisontothefuselage,canleadtolargedeformations.Thus,asrstproposedbyMeirovitch[ 12


Thelevelofdelityofthemodelisalsoaconcern.Intheworkpresentedhere,idealizationswillbemadetosimplifythederivationoftheequationsofmotion.Asanexample,themassesofhingesandactuatorsthatallowforarticulationwillbeignored.Theequationsofmotionthemselveswillalsobeformulatedinaminimalbasisinordertoeliminateconstraintreactionsbetweenbodies.Thesesimplicationsgreatlyreducethecomplexityofthenalsetofequations,withoutneedlesslyoversimplifyinganysimulationsofthevehicleresultingfromtheequationsofmotionderived. Thenalstepinthederivationoftheequationsofmotionistheformulationofthegeneralizedforces.Aerodynamiceectsaccountformostoftheforcesthattheexibleportionsofthemorphingmicroairvehicleexperiences.Theexactnatureoftheseforceshavetobeidealizedforthescopeofthiswork.Theactualforceswhichthevehicleissubjectedisacontinuumthatexistsovertheentiresystem.However,theseforcesandthemomentstheyproducecannotbeaccuratelymeasuredovertheentiresystem.Atthispointintheresearch,itisassumedthataerodynamicforceswillberepresentedbystaticpressuresataparticulartimeatdiscretepointsalongthevehicle.Includedinthisforcemodelingistheassumptionofhowtoconvertpressuredataintoaerodynamicforces.Ultimately,thepressuremeasurementswillbeenoughtoapproximatelyobtaintheactualaerodynamicforcesandmomentsthatcharacterizeofthemorphingMicroAirVehicleatanypointintimeandbasedonthecongurationofthecraft. Thesecondstepinthesolutionoftheproblemstatementiscastingtheequationsofmotioninaformthatisconvenientforcontrol.Fromtheassumptionofaniteelement 13


Inthefollowingchapters,thecourseofthisresearchwillbefurtherexplored.Chapter 2 presentsabrieflookatthemajortopicsinvolvedinthisresearchandsomeoftheworkthathasbeendoneinthoseareas.Chapter 3 presentsadetailedderivationoftheequationsofmotion.Chapter 4 explainstherepresentationoftheforcesusedintheequationsofmotion.Chapter 5 thencaststheproblemasacontrolsproblemalongwithadetaileddiscussionoftheassumptions,goals,andlimitationsassociatedwithdoingso.Chapter 6 thenappliesthederivationoftheequationsofmotionandthecastingofsaidequationsasacontrolproblemtoanexamplebasedonGull-wingMAVillustratedinFig. 1-1 .Chapter 7 presentsanotherexamplebasedonathreebarmanipulator.Finally,Chapter 8 exploressomeconclusionsofthework. 14


Thesecondarycontributionofthisworkisinthecontrolstrategiesusedonmicroairvehicles.Theequationsofmotionasderivedinthisresearcharehighlynonlinearandtime-dependent.Standardcontrolstrategieswouldchooseanoperatingpointtolinearizeabout,thencreateacontrolstrategybasedonthatpoint,andnallytunetheperformanceforconditionsthatdon'tmatchthelinearizationpoint.Thistypeofcontrolstrategy,evenwhenusingmultiplelinearizationpointsandswitchingcontrollers,doesnotachievethefullperformancepotentialforamorphingMAV.Morphingallowsforaninnitenumberofcongurationswithassociatedperformance.Thus,acontrolstrategythattakesallofthecongurationsintoconsiderationishighlydesirable.Theequationsofmotionasderivedinthisworkcanthenbecastintoanoptimalcontrolproblem.Thesolutionoftheoptimalcontrolproblemnotonlyndsacontrolstrategyforthevehicle,butdependingonthecostfunctionalchosencanpresentaformatthatallowsthemorphingofthevehicletobeusedastheonlymeansbywhichtoobtainsomedesiredmaneuver. 15


MorphingMAV. 16


BelowachordReynolds'numberof50,000,thefreeshearlayerafterseparationdoesnottransitiontoturbulentowsoonenoughtoreattachtotheairfoilsurface.[ 17


TheaerodynamicsofMAVsinlowReynolds'numberregimediersgreatlyfromtheaerodynamicsofminivehicles.[ ThereareseveraldierentclassesofMAVs.InapaperbyMuellerandDelaurier,severaldierentMAVcongurationsarepresented.[ 18






Severalinitiativeshavebeenstartedbyorganizationstoinvestigatemorphing.AccordingtotheprogramatDARPA,amorphingvehiclewillchangeitsshapesubstantiallytoadapttothemissionenvironment,providesuperiorsystemcapabilitynotpossiblewithoutreconguration,anduseintegrateddesignofmaterials,distributedactuators,eectorsandmechanismstorecongureinight.[ 21


InworkbySaundersetal.conductedexperimentsona16%scalemodelofanF-18wingthatusedshapememoryalloytoeectasmoothdeformationofthetraillingedgecontrolsurfaces.[ InpapersbyAbdulrahimetal.,studieswereperformedtoevaluatetheuseoftorquerods,toeectwingtwist,onthecontrolauthorityofaMAV.[ MAVstypicallyhavewingswithlowaspectratios(AR).[ 22






3-1 representsgraphicalcomputermodelofthearticulateMAV.Thisvehicleiscomprisedofacorebody,thefuselage,andfourwingsegments.Thefuselageandeachofthefourwingsegmentshavetheirowncoordinatessystem.Theequationsofmotionforthissystemaredesired.TheseequationsneedtobeabletodescribethecompletelocationandorientationoftheMAVinsomeinertialframe.Sincethegoalistobeabletowritetheequationsofmotionofthisvehiclewithrespecttoaninertialframe,thetransformationmatricesmustbedenedforeachbodythatrelatetheinertialframetothebody-xedframe. Forthepurposesofthederivationoftheequationsofmotion,thevechicleisidealizedasisillustratedinFigure 3-2 .Thevehicleissymmetricaboutthe^b11^b13plane,thusonlyasinglesideisconsideredintheidealization.Inthefullsimulationofthesystem,asymmetricmodesdooccurandthecontributionfromtheotherhalfoftheMAVcanbeignoredforthepurposesofthederivationofthequationsofmotionwithoutlossofgenerality.Thefollowingderivationswillonlycontainthefuselage(body1),therightinboardwingsegment(body2),andtherightoutboardwingsegment(body3)andtheassociatedcoordinateframes,B1;B2;andB3respectively.Eachbody-xedcoordinateframe,B1;B2;andB3,originateatthecenterofmassofthierrespectivebody.Therelationshipsbetweenthebody-xedcoordinatesystemsofthesebodiesandtheinertialcoordinatesystem,E,areasfollows: ^e ^e ^e 26


[C1]=266664c1c1s1c1+c1s1s1s1s1+c1s1c1s1c1c1c1+s1s1s1c1s1+s1s1c1s1c1s1c1c1377775(3{4) wherec1iscos1ands1issin1. Inadditiontotheinertialcoordinatesystemandthebody-xedcoordinatesystemsassociatedwitheachbody,thereareintermediatecooridinateframesthatarealignedwiththejointsbetweenbodies.ThesejointcooridnatesystemsareillustratedinFigure 3-3 andFigure 3-4 .Therelationshipbetweenthebodiesandthejointcoordinateframesare: 27


^j ^j 28


Fortheminimalsetformulationallbodieswillbereferencedtoasinglebody.Inthiscase,allbodieswillbereferencedtotheB1frameassociatedwiththefuselage.TheequationsofmotionwillbederivedusingLagrange'sEquations.Lagrange'sequationsare: dt@T whereTisthekineticenergy,Visthepotentialenergy,RistheRayleighdissipationfunction,q 29


2v 2v wheretheintegralistakenoverthevolumeofthebodyandv whereR 3{13 yields dtr dt(R dtu dtu SincethefuselageisrigidthepointuopdoesnotchangeintheB1basis.Thus, where[0~!1]isaskewsymmetricmatrixcomposedofthecomponentsofE!B1.Theskewsymmetricmatrixcanbedecomposedas: Thus, 30


3{17 sothatalltermsareintheB1basis: Let [B]=0~!1u Then where Thusthekineticenergyforthefuselageis 2ZV_R Let [M]=264RV[I]dVRV[C1][B]dVRV[B]T[C1]TdVRV[B]T[B]dV375(3{23) where[I]istheidentitymatrixandlet _q Thenthekineticenergyforthefuselageis 2_q where[M]isthemassmatrixforthefuselage.Nextexaminetheintegralstobeperformedinthecomputationofthemassmatrixforthefuselage.Tofacilitatethisexaminationthe 31


whichisaconstantmatrix. LookingatthersttermofEq. 3{27 whereRVuop;ndVisconstant.Forthesecondandthirdtermssubstitutein0~!1and0~!1whereappropriate. Lookingatthetermintherstrowandrstcolumnoftheresultingmatrixyields whereRVuop;iuop;ndVisconstant.Fortheremainingtermssubstitutein0~!1and0~!1whereappropriate. Nextcreatethekineticenergyoftheinboardwing.Thepositionvectortoanarbitrarypointontheinboardwingisasfollows: wherer 32


dtr dtr dtr dtu dt[N2]q Byvirtueofsomeofthevectorsbeingconstantintheframethatthederivativeistakenin dt[N2]q Throughtheuseofsimpleangularrotationsbetweenframes, However,B1! .Thus,separatingtheangularvelocitytermsandre-writingallcrossproductsintermsofthesamebasisyields: Let[0~!1]beaskewsymmetricmatrixcomposedofthecomponentsofE!B1and[12~!21]beaskewsymmetricmatrixcomposedofthecomponentsofJ12!J21.Thus, 33


re-writingallEq. 3{36 intermsoftheB1coordinatesystemyields However,theskewsymmetricmatricescanbefactoredasfollows: Thus, 34


UsingEquations 3{41 through 3{44 ,Equation 3{40 canbewrittenas: Dene [B]=B 35


Nowthekineticenergyofbody2canbewritten. 2v 2v SubstitutingEq. 3{48 intoEq. 3{49 yields 2_q where _q 2266666664[I][C1][B][C1]D D PleaserefertoAppendixAtoseethedetailsoftheintegralinvolvedintheformationofthemassmatrixforbody2. Aswasdonewithbody1andbody2,thecontributiontothesystemequationsofmotionfrombody3willnowbederived.Thustheequationsofmotionforbody3proceedfromapplicationofLagrange'sequations.Thisprocessbeginswiththederivationofthe 36


wherer dtr dtr dtr dtr dtr dtu dt[N3]q AswiththeangularvelocityoftheB1framewithrespecttotheinertialframe,theangularvelocityoftheB2framewithrespecttotheinertialframecanbewrittenasthesumofsimpleangularrotations.Thus, However, (3{56) Thus, 37


Therefore, Expandingthecrossproductsandexpressingthecomponentsofthecrossproductsinthesamebasesyields Introducingskewsymmetricmatricestosimplifythecrossproductexpressionsandre-writingthewholeequationtobeintermsoftheB1basisyields: 38


where[0~!1]istheskewsymmetricmatrixcomposedfromthecomponentsofE! Thus, 39


Dene 40


[BB]=BB 41


InsertingEq. 3{70 intotheexpressionforthekineticenergyyields: 2_q where _q and Nowthatthekineticenergyforeachbodyhasbeenderived.Theycanbecombinedintoanexpressionforthekineticenergyfortheentiresystem.Theglobalkineticenergyis 42


2_q where _q ThenextpieceofinformationneededforLagrange'sequationisthepotentialenergyofthesystem.Aswiththekineticenergyofthesystem,thepotentialenergyofthesystemcanbedecomposedintothepotentialenergyoftheindividualbodiesthatmakeupthesystem.Forthepurposesoftheworkpresentedinthisdocument,thepotentialenergyduetochangesinelevationwillbeignored.Thus,thepotentialenergyofthesystemwillbearesultoftheexibledegreesoffreedom.Asbody1isconsideredrigid,ithasnocontributiontothepotentialenergyofthesystem.Thecontributionstothepotentialenergyofthesystemfrombodies2and3havethesameform.Thearetheresultoftheassumptionoftheniteelementsusedtorepresentthebodies.Thesymbolicformofthepotentialenergycontributionsare: 2f" where"isavectorofthestrainsand[E 43


AssumethereissomeforceF whereu but Thusforapointonbody1: Re-writingEq. 3{81 intermsoftheB1frameyields: 44


3{82 andthencomparingtheexpressiontoEq. 3{77 ,thecontributionstothegeneralizedforcesfromtheforcesactingonbody1canbederived. Thederivationofthecontributiontothegeneralizedforcesfromtheforcesactingonbody2proceedswithndinganequationforthevirtualdisplacementofapointonbody2.Thepositionvectorofapointonbody2is: WritingEq. 3{83 intermsoftheinertialframeproduces: ThevirtualdisplacementofEq. 3{84 is 45


Thus Re-writingEq. 3{87 intermsoftheB1basisandusingthematricesintroducedinthederivationofthecontributiontothesystemkineticenergyfrombody2: q Performingthedotproductoftheforcesactingonbody2withthecorrespondingvirtualdisplacementvectoroftheformofEq. 3{88 andcomparingtheresultwithEq. 3{77 producesthecontributionstothegeneralizedforcesfromtheforcesactingonbody2. Thederivationofthecontributiontothegeneralizedforcesfromtheforcesactingonbody3proceedsfromtheprincipleofvirtualwork.First,thepositionvectortothepoint 46


Re-writingEq. 3{89 sothatallcomponentsareintermsoftheinertialframeyields Takingthevirtualdisplacementproceedsas:


Thus, 48


Usingthevectorsandmatricesdenedinthederivationofthecontributiontothekineticenergyfrombody3,andre-writingthevirtualdisplacementofapointonbody3intermsoftheB1basisyields: q 49




ArticulatedMAV. Figure3-2: Coordinateframesofthefuselage,inboardwingsegment,andoutboardwingsegment. 51


Jointcoordinateframesbetweenfuselageandinboardwingsegment. Figure3-4: Jointcoordinateframesbetweeninboardandoutboardwingsegments. 52


Sincethewingstructuresaremodeledasexiblebodies,thegenerationofaeroelasticforcesbecomesimportant.Theaeroelasticforcesusedinthemodelwillinuencetheshapeofthewings,duetotheirexibility,andtheshapeofthewingswillinuencetheaerodynamicforces.Therearemanydierentschemesavailableforthegenerationofaeroelasticforces.Theseschemesrangefrompotentialowtheorytofull3-DNavier-stokescomputationaluiddynamicssimulations.Thegenerationofaeroelasticforcesisasucientresearchtopicinitsownright.WorkbyOrtizandBarhorstdemonstratesanumberofdierenttechniquesinthegenerationofaerodynamicforcesforuseinaeroelastictystudies.[ 3 53


AVL'sowanalysisallowsforthecalculationofseveraldierentparameters,suchasstabilityderivativesandhingemoments,baseduponseveralinputsthatdenetheightcondition.Theseinputsincludepitchangle-,sideslipangle-,pitchrate,rollrate,yawrate,controlsurfacedeectionangle,andwindvelocity.Foragivensetofinputs,AVLdisplaystheforcesandmomentsactinguponabodycenteredcoordinatesystem.Thiscoordinatesystemhasthex-axisrunningalongthefuselageofthemodelandthey-axisalongtherightwingofthemodel.Theindividualelementforcesand/orstripforcescanbewrittentolesasrequestedbytheuser.AVLalsocanperformseveraldierentoperationsontheinputedmodelsuchasgeometryplotting,eigenmodecalculation,andTretzPlaneplotting. Thus,theaerodynamicforcesandmomentsproducedbyagivenAVLrunassumesthattheightconditionsasdenedbythegeometryandtheinputparametersdonotchange 54


Inaddition,AVLisonlyapplicableforlowMachnumberow.AnAVLanalysiscanbecompletedwithouttheuseofafuselageinthegeometrydenition.Whenafuselageisincluded,itistreatedasablubodyandasaresulttheforcescalculatedthatactonthefuselagearenotquiteasaccurateasamorerobustanalysistoolwouldgenerate. Theseshort-comingsindicatethattheuseofAVLisnotoptimalforthegenerationofaerodynamicforcesinthesimulationofanunsteadyprocessofaeroelasticityofare-congurableMAV.Instead,morerobustmeasuressuchastheuseofafully3DNavier-StokesCFDsoftwarewouldbemoreappropriate.However,theuseofsuchasoftwareonare-congurableMAVisaresearchtopicinitself.Also,thecurrentstate-of-the-artincomputationaluiddynamicsisverycomputationallyexpensivetogenerateasolutionfortheforcesactingovergivengeometry.Firstanappropriatediscretizationofthegeometryandtheuidsurroundingthegeometrymustbegenerated.Oncethediscretizationhastakenplace,thesolutionoftheaerodynamicforcesisaniterativeprocessandmaytakealargeamountoftimetoconvergetoasinglesolution.Thisproblemiscompoundedwhenusedinanaeroelasticsimulationasthecongurationofthevehiclechangesateverytimestepandthecomputationalcostneededtogeneratetheaerodynamicforcesmustalsobepaidateverytimestep.Incontrast,AVLtakesamatterofsecondstoconvergetoasolutionforagivengeometryinput,duetothesimplicityofthelatticevortexmethod.Thusasacomputationallyinexpensivesourceofaerodynamicforces,AVLisusedeventhoughproperapplicationoftheprocesswouldchooseamorerobustmethod. 55


Inthepreviouschapter,thecontributionstothekineticenergyandthepotentialenergyfromthevariousbodiesthatcomposethesystemwerederived.Thegeneralizedformofthecontributionstothegeneralizedforcesfromthevariousforcesactingonthesebodieswasalsoderived.Thegoalofthissectionistocombinethesecomponentsintoasystemofequationsthatgovernthemotionoftheentirevehicleandusethisformulationinthederivationofacontrolstrategythatwillbeusedtodrivethevehicletoadesiredstate.Anoptimalcontrolstrategycanbeusedtoaddressanumberofcontrolobjectives,fromtrackingtoregulation. Thereareseveraldierentwaysinwhichtheoptimalproblemcanbesolved.However,thesemethodsgenerallyfallwithintwoclassications.Theseclassicationsaredirectmethodsandindirectmethods.Muchworkhasbeendoneintheuseofindirectmethodsinthesolutionofoptimalcontrolproblems.[ ;u Where(x _x wherex 5{2 showsthatthestatesx 56


5{1 isalsoafunctionofthecontrol.Optimalcontroltheorystatesthattheoptimalcontrol,u foreveryadmissiblecontrol,u ;u ; ;u where()areknownasthecostatesandaretheLagrangemultipliersoftheboundaryconditions.TheaugmentationofthecostfunctionalisstandardproceduretoconvertEquation 5{1 subjecttotheconstraintsofEquation 5{2 intoanunconstrainedoptimizationproblem.NowdenetheHamiltonianfunctionas: ;u ;)=L(x ;u ;u Thus,theaugmentedcostfunctionalbecomes: Afterapplicationofthecalculusofvariationtotheaugmentedcostfunctionalthefollowing: _x (5{9) WhereEquation 5{7 representsthestatedynamics,Equation 5{8 representsthecostatedynamics,andEquation 5{9 representstheconditionsforoptimalcontrol.Theseequationsformatwo-pointboundaryvalueproblem,thesolutionofwhichdeterminestheoptimalcontrolu 57


5{5 5{9 isrequired.Thedicultyinsolvingthoseequationsliesingettingaclosed-formsolutionofthecostateequations.Theadvantagesofindirectmethodsisthatthesolutionofthecostateequationsservesasagoodindicatoroftheoptimalityofthesolution.Numerically,indirectmethodsaretypicallysolvedusingshootingormultipleshootingmethods.Theseshootingmethodshavemuchincommonwithstandardrootndingtechniques. Directmethods,alsoknownasdirecttranscription,solvetheoptimalcontrolproblemsbyconvertingtheoptimalcontrolproblemintoanonlinearprogrammingproblem.Thetranscriptionprocessconvertstheinnite-dimensionaloptimalcontrolproblemintoanite-dimensionaloptimizationproblem.Thedynamicsoftheproblemarediscretizedatanitesetofpointsandthecostfunctionalisreplacedwithacostfunctionthatcanbeevaluatedatthosediscretepoints.Thisnite-dimensionalproblemisknownasaNonLinearProgram(NLP).Directmethodshavetheadvantageofnotrequiringthecostateequationstobeexplicitlyderived.Inaddition,thedirectmethodallowsforfasterconvergenceofthesolution.Asaresult,directmethodsareusedquiteofteninpractice,asillustratedinworkbyBass,vonStryk,Hu,Frazzoli,andHull.[


3 beginswiththeconversionofthesecond-orderequationsofmotionintoarst-orderformthatisstandardforoptimalcontrolproblems.Thesecond-orderequationsofmotionarewrittenas: [M]q Where 1122q 28>>>><>>>>:_q 59


Thustheequationsofmotioncanbewrittenas: _x ;u where Nowthattheequationsofmotionarewritteninrstorderform,itisnecessarytoexplicitlydenethecontrolproblem.Anexampleofagoalofthecontrolproblemistodrivethevehicletoadesiredpositionusingthecontroltorqueslocatedathingeaxestochangethecongurationofthevehicle.Thischangeofthecongurationwilltheninuencetheaerodynamicsofthevehicle.Writingthisinmathematicaltermsproducesthefollowingobjectivefunction: ThiscostfunctionseekstominimizethedistancebetweenthenalglobalpositionoftheMAVandsomedesirednalglobalposition.TheoptimalcontrolisthenthecontrolthatminimizesEquation 5{16 .Theobjectivefunctionissubjecttophysicalconstraintsonthesystemandmustalsosatisfytheequationsofmotion.Anexampleofphysicalconstraintsistolimittherangeofmotionontheinboardandoutboardsegmentsto90degwithrespecttotheinboardbody(fortheinboardsegmenttheinboardbodyisthefuselageand 60


minJ=jXCB1XDj2+jYCB1YDj2+jZCB1ZDj2subjectto_x ;u 6 presentsanumericalexampleofanarticulatedMAV.Includedinthisexampleistheformationofanoptimalcontrolproblem. 61


ThederivationsoftheequationsofmotionthatwerederivedinChapter 3 wereappliedtoanumericalexampletoshowthattheequationscanthenbecastintoaformamenabletoanoptimalcontrolproblem.ThisexampleisimplementedinMATLAB.Therststepinformingthisexampleisgeneratingageometricmodel.Thegeometricmodelisconstructedtoexhibittheprominentfeaturesoftheactualmodel,suchasthearticulationofthewing,butwithoutthecomplexitiesassociatewiththefullmodel.Issuessuchasjointfrictionandcompositematerialpropertiesareignoredwithoutlossofgenerality. 6-1 thefuselageismodeledasarectangularprismandthewingsegmentsaswellasthehorizontalandverticaltailaremodeledasrectangularplates.AlthoughtheidealizedmodelisasimplicationoftheactualarticulatedMAV,thelengthsandmassesoftheidealizationwaschosentoreecttheactualcharacteristicsofthefullvehicle.Themodelfuselagehasalengthof0.3048metersandawidthandheightof0.0636meters.Theinboardwingsegmentshaveaspanof0.2032metersandachordlengthof0.0762meters.Theoutboardwingsegmentshaveaspanof0.1016metersandachordlengthof0.0762meters.Thehorizontaltailsegmentshaveaspanof0.127metersandachordof0.06858meters.Theverticaltailsegmenthasaspanof0.1016metersandachordof0.06858meters.Allrectangularsegmentshaveathicknessof0.012192meters.Thedensityofallcomponentsisassumedtobe137:47kg m3togivethetotalsystemamassof0.29031kg. Tofurthersimplifythenumericalexample,exibilityeectsareignored.Eliminatingexibilityhastheeectofreducingthegeneralizedcoordinatesaswellassimplifyingthegenerationofthematricesinvolvedintheequationsofmotion.Thissimplicationdoesnotinvalidatetheobservationsthatresultfromtheanalysisoftheidealizedmodel.Rather,theexclusionofexibilityeectsgreatlyreducescomputationaleort.For 62


A ,arereducedinnumbertoonlybeingRVdV;RVuop;pdV;andRVuop;puop;zdV.Similarly,thegeneralizedcoordinatesarereducedtoq wheretisthethicknessoftheelement,anddetJisthedeterminantoftheJacobian.Usinganisoparametricelement,theJacobianbecomes 4((1i))x1+(1i)x2+(1+i)x3(1+ti)x41 4((1i))y1+(1i)y2+(1+i)y3(1+i)y41 4((1i))x1(1+i)x2+(1+i)x3+(1i)x41 4((1i))y1(1+i)y2+(1+i)y3+(1i)y435(6{2) 63


Similarly, whereuop;iisthepositionofthequadraturepoint.Lastly, Nowthattheinvariantsarecalculated,thecurrentcongurationofthemodelisusedtodeterminethevariousquantitiesintheequationsofmotion.Forthepurposesofthisnumericalexample,thecongurationdependenceissimpliedasmuchaspossible.Forexample,intheinitialcongurationthecoordinatesystemsassociatedwitheachbodyisalignedwiththeglobalcoordinatesystem.Furthermore,thewingsegmentcoordinateframesarealignedwiththerespectivejointcoordinateframelocatedonthatbody.Thus,theonlychangeincoordinateframesaredeterminedbythechangeincoordinateframesbetweenthefuselageandtheglobalcoordinatesystem,determinedbytheangles,,and,andthechangesbetweenthejointcoordinateframeononebodyandthejointcoordinateframeontheconnectedbody,determinedbytheangles1,1,2,and2.Oncetheappropriatetransformationmatricesareformed,thisinformationalongwiththerelevantpositionvectorsandthecalculatedinvariantsareusedinformationofthemassmatrixofeachindividualsegment.ThisprocessisexpandedindetailinAppendices A and B .Theseindividualmassmatricesarethencombinedtoformthesystemmassmatrix.Similarly,thepartialandtimederivativesofthevariousmatrices,positionvectors,andinvariantsareusedtoformvarioustimeandpartialderivativesofthesystemmass 64


3 forfurtherdetail). Forthepurposesofcalculatingaerodynamicforcesonthisexample,anAVLinputleiscreatedthatmimicsthegeometryoftheexample.However,inorderforthiscongurationtoproduceaerodynamicliftforces,thewingmusthavecamber.Tothisend,inthe.avlinputle,thewingsegmentsaremadetohaveacamberprole.AVLalsorequiresinputparametersthatdenetheruncase.Theseinputsarethevelocityofthevehicle,theangleofattackofthevehicle,andtheangleofsideslipofthevehicle.Thesequantitiesarederivedfromthestatesofthegeneralizedcoordinatesasfollows:velocity:p uangleofsideslip:sin1v V Atthispoint,thegeometryoftheexampleisfullydened.Fromthisgeometry,thesystemmatricescanbeformedusingthederivationsinChapter 3 ,Appendix A ,andAppendix B .Thesesystemmatricesarethemassmatrix,thersttimederivativeofthemassmatrix,andthepartialderivativeofthemassmatrixwithrespecttoeach 65


5{10 .Nexttheappropriatequantitiesontherighthandsideoftheequationsofmotionmustbedened.Thegravitationalforcesareappliedtoeachmulti-bodysegmentusingthemassesofeachsegmentandlocalaccelerationofgravity.Thesegravitationalforcesareappliedtothecenterofgravityoftherespectivesegments.Todenethepropulsiveforces,AVLisexecutedusingageometryinputlethatcorrespondstonodeectionofboththeinboardandoutboardsegments.Inaddition,thevelocityisdenedinAVLtobe15m s.ThisvelocityischosentobeassociatedwiththeMAVyinginatrimconditioninstraightandlevelight.Thus,thecoecientsgeneratedbythiscongurationshouldproducealiftforcethatexactlybalancestheweightforceandadragforcethatbalancesthethrustforce.Themagnitudeofthisdragforceisthenusedasthemagnitudeofthepropulsiveforce.FollowingthederivationofthegeneralizedforcesasoutlinedinChapter 3 thethrustforceisassociatedwithactingonthefuselage.Thedirectionoftheforceisassumedtoactthroughthecenterofgravityoftheentiresystem.ThecenterofgravityisfoundtoalwayslieinthefuselageinallcongurationsofthevehiclethatareadmissibleasdenedbytheconstraintsontheanglessetforthinChapter 5 .Foropenloopsimulation,theanglesoftheinboardandoutboardwingsegmentsareconstrainedtobezero. Nowallofthenecessarycomponentsforanopenloopsimulationofthevehiclearepresent.ItwaschosentouseaEulerintegrationschemetointegratetheequationsofmotion.TheEulerintegrationschemehastheadvantageofonlyneedinginformationfromthecurrenttimestepinordertocalculatethevaluesatthenexttimestep.However,caremustbetakeninthedeterminationofthetimesteptoensureconvergenceofthesolution.Toensureconvergence,thetimestepmustbeasfollows:ttcr




;u 5{14 andEquation 5{15 Inthedeterminationofthesolutionoftheoptimalcontrolproblem,therearemanydierentcongurationsthatwouldbeselected.Suchalargedesignspaceincreasesthecomputationaleortneededtosolvetheproblem.Toalleviatesomeofthiscomputationalburden,thisexampleisfurtherconstrainedtofollowsymmetricdeectionsonly.Thisconstraintremovestwodegreesoffreedomandtwocontrolinputs,associatedwiththedeectionoftheoppositewing. Adirectmethodwasusedinthesolutionofthisoptimalcontrolproblem.ThismethodusedaGaussPseudospectralmethodinthesolution.ThismethodusestwentyLegendre-Gausspointsinthediscretizationrepresentthesolutionoftheoptimalcontrolproblem.Thenaltimeforthesystemwasfreetobedeterminedbytheoptimizationprocedurethatminimizedcontroleort.ThisoptimalcontrolproblemwasimplementedintheMATLABprogramusedinthesolutionofstandardoptimalcontrolproblems.[ 68


Figure 6-2 showstheglobalpositionofthevehiclegeneratedbytheoptimizationprogramafteritwasunabletoimprovethecostfunctionoverthepreviousiteration'scalculatedcontrolvector.Similarly,Figure 6-3 showsthemorphanglesgeneratedbytheoptimizationprogram.Theplotsproducedfromtheoptimizationsoftwareseemstoindicatethedesireddecreaseinaltitudeof5mwiththeconstraintofthemorphanglesreturningtozerodeectionatthetimethevehiclecompletesthedecreaseinaltitude.However,theseresultsdonotreectthesuccessfulsolutionoftheoptimalcontrolproblem.ThisfactisillustratedbyFigure 6-4 andFigure 6-5 .Theseguresshowthepositionofthevehicleandthemorphanglesafterintegratingtheequationsofmotionusingtheoptimalcontrolasdeterminedbytheoptimizationsoftware.Asillustratedinthegures,notonlydoesthevehiclenotcompletethe5mdecreaseinaltitude,butthemorphanglesarealsonon-zeroattheterminationofwhattheoptimizationprogramsstatesasthemaneuvertime. 69


Figure6-1: Geometricmodelofgull-wingMAV. 70


PositionofMAVwhenusingtheoptimalcontrolsolutionasdeterminedbyGPOCS. Figure6-3: MorphingangleswhenusingtheoptimalcontrolsolutionasdeterminedbyGPOCS. 71


PositionofMAVfromintegratedoptimalcontrol. Figure6-5: Morphinganglesfromintegratedoptimalcontrol. 72


TheequationsofmotionthatweregeneratedinChapter 3 arealsoapplicabletoothermulti-bodysystems.Oneofthemostcommonexamplesofamulti-bodysystemisthethreebarmanipulator.Thissystemconsistsofthreebarsthatareconnectedviasimple,singledegreeoffreedomlinkages.References[ 90 ]and[ 91 ]havesectionsinwhichthiscommonexampleisstudied. 7-3 ,therstbarisconnectedtoaninertialframeviaarevolutejointthatallowsthebartorotateabouttheaxisalignedwiththebar.Thesecondbarisconnectedtotherstbarbyarevolutejointthatallowsthebartopivotaboutanaxisthatisperpendiculartotheaxisalignedwiththebarandperpendiculartotheplaneillustratedinthegure.Likewise,thethirdbarisconnectedtothesecondbarbyarevolutejointthatallowsthebartopivotaboutanaxisthatisperpendiculartotheaxisalignedwiththebarandperpendiculartotheplaneillustratedinFigure 7-3 .Attheinitialtimeallofthebody-xedcoordinateframesarealigned. Inthisexamplesomeassumptionsaremadetofacilitatetheexecutions.Primaryamongtheseassumptionsisthemodelingeachbarasarigidbody.ThisassumptioneliminatestheneedforaniteelementrepresentationofthebarsandthuseliminatestheglobalsystemstinessmatrixandthetermsintheapplicationofLagrange'sequationsinthederivationoftheequationsofmotionthatresultfromtheexibledegreesoffreedom.Therevolutejointsconnectingthebarsareassumedtobeideal.Thus,theyhavenomassandonlyallowforasingledegreeoffreedombetweenthebodiestheyconnect.Alsoincludedintheassumptionofidealjointsisthelackoffrictioninthejoints.Itisfurtherassumedthatthebarsareallofonemeterinlengthandarealsoidealizedtohaveaunitarydensityandazeroradius.Finally,itisassumedthatnoexternalforcesactonthesystem. 73


3 .Inusingthesetwoformulationsthestrengthsandweaknessescanbecontrasted. 7-4 .Mathematically,thebody-xedcoordinatesystemsarerelatedtotheinertialcoordinatesystembytheuseofatransformationmatrix.Eachtransformationmatrix,isitselfafunctionofthreeangles.Thus, wheree [C1]=266664c1c1s1c1+c1s1s1s1s1+c1s1c1s1c1c1c1+s1s1s1c1s1+s1s1c1s1c1s1c1c1377775(7{4) 74


[C3]=266664c3c3s3c3+c3s3s3s3s3+c3s3c3s3c3c3c3+s3s3s3c3s3+s3s3c3s3c3s3c3c3377775(7{6) wherec1iscos1ands1issin1.Asallthreebodieshavesimilarrelationsbetweentheirrespectivebody-xedframeandtheinertialframe,theyhavesimilarrelationsfortheangularvelocityofthosebody-xedframeswithrespecttotheinertialframe.Theserelationsare: wheretheangularvelocitiesarewrittenintermsofthebasisvectorsofthebody-xedcoordinatesystem. ThederivationoftheequationsofmotionstartswithapplicationofLagrange'sequations.Forarigidbody,Lagrange'sequationstakestheformof: dt@T WhereTisthekineticenergyofthesystem,q 2v 75


whereR 7{12 withrespecttotheinertialframeproduces: dtr dt(R dtu where[0~!1]isaskewsymmetricmatrixformedfromthecomponentsoftheangularvelocityofbody1withrespecttotheinertialframe.However, Thus, Re-writingEquation 7{15 soallcomponentsarewrittenintermsoftheinertialframeproduces: 76


Let Therefore SubstitutingEquation 7{21 intoEquation 7{11 produces: 2RTCB1_1_1_1266666664RV[I]dVRVB 7{22 areexpandedas: 77


ItisimportanttonotethattheintegralsontherighthandsideofEquation 7{23 throughEquation 7{35 areallinvariantintegrals.Forbody1,thegeneralizedcoordinatescanbedenedas: ThusEquation 7{22 takesonthefamiliarquadraticformof: 2_q Thederivationstothekineticenergyofthesystemfrombodies2and3proceedinthesamemannerasthederivationofthecontributionfrombody1illustratedintheabovewith0~!1x1beingreplacedwith1~!2x2forbody2and2~!3x3forbody3.Thesystem'sgeneralizedcoordinatesarecomposedofthegeneralizedcoordinatesoftheindividual 78


Thusthesystemkineticenergytakesontheformof 2_q SubstitutingEquation 7{39 intoEquation 7{10 resultsin: [M]q whereQ 79


Thethreebodiesarenotcompletelyfreetomoveintheinertialframe.Throughtheuseoftherevolutejoints,thebodiesareconstrainedtomoverelativetoeachother.Sinceinthestandardderivationoftheequationsofmotionthebodiesarenotdenedrelativetoeachother,anothertermappearsinLagrange'sequationstoaccountfortheconstraints.Lagrange'sequationsbecome: dt@T where[Cq]isthematrixthatresultsfromtakingthederivativeoftheconstraintmatrixwithrespecttothegeneralizedcoordinatesalsoknownasthesystemJacobianmatrix.Arevolutejointonlyallowsforasingledegreeoffreedombetweentwobodies.Thus,eachrevolutejointhasveconstraintequationsassociatedwithit.Thissystemhasthreerevolutejoints,asillustratedinFig. 7-3 ;onebetweenthegroundandbody1thatallowsrotationabouttheb13axisofthebody-xedcoordinatesystemassociatedwithbody1,onebetweenbody1andbody2thatallowsrotationabouttheb21axisofthebody-xedcoordinatesystemassociatedwithbody2,andonebetweenbody2andbody3thatallowsrotationabouttheb31axisofthebody-xedcoordinatesystemassociatedwithbody3. Theconstraintmatrixisdevelopedbysettingtheconstraintequationsthatgovernthemotionofthebodiesequaltozero.Thus,forthesystem: []=0 (7{43) 80


7{43 produces: []q_q (7{44) where[]qisthepartialderivativeoftheconstraintmatrixwithrespecttothegeneralizedcoordinates,alsoknownasthesystemJacobianmatrix,and[]tisthepartialderivativeoftheconstraintmatrixwithrespecttotime.Forautonomoussystems,Equation 7{44 reducesto: []q_q (7{45) TakingthesecondtotalderivativeofEquation 7{43 orthersttotalderivativeofEquation 7{44 produces, []q t[]q whichforautonomoussystemsreducesto []q AnysystemthatsimultaneouslysatisesEquation 7{44 andEquation 7{46 alsosatisesEquation 7{43 .Thus,usingtheaccelerationsofthegeneralizedcoordinates,Equation 7{46 canbeaugmentedtothesystemequationsofmotiontoformthematrixequation Q whereQ 7{46 Theconstraintmatrixforthethree-barmanipulatorisderivedinthefollowing.Therevolutejointbetweenthegroundandbody1restrictsthepositionoftheendofbody1tobecoincidentwiththepointonthegroundthatcorrespondstotherevolutejointlocation. 81


(7{49) whereu 7{49 intermsoftheinertialframeproduces: TakingthepartialderivativeofEquation 7{50 withrespecttothegeneralizedcoordinates,asstatedinEquation 7{38 ,yields: 0 0 [0]0 0 0 Therevolutejointalsointroducesconstraintsthateliminatetwooftherotationaldegreesoffreedom.Thustwovectorsalignedwiththeaxisofrotationbetweenthebodies,onethatisxedintherstbodyandthesecondbeingxedinthesecondbody,willalwaysbeparallel.Thisstatementisequivalenttosayingthatavectoralignedwiththeaxisofrotationandxedintherstbodywillalwaysbeperpendiculartoavectorthatisxedinthesecondbodythatisperpendiculartotheaxisofrotation.Mathematicallythesestatements,astheyrelatetotherevolutejointbetweenthegroundandbody1canbeexpressedas: 2=e Orpresentingallquantitiesintermsoftheinertialframe: 2=e 82


7{53 withrespecttothegeneralizedcoordinatesproduces: Similarly, 3=e Orpresentingallquantitiesintermsoftheinertialframe: 3=e TakingthepartialderivativeofEquation 7{56 withrespecttothegeneralizedcoordinatesproduces: Fortherevolutejointbetweenbody1andbody2theconstraintontherelativetranslationalmotionbetweenthebodiesisexpressedas: (7{58) whereR (7{59) 83


7{59 withrespecttothegeneralizedcoordinatesyields: 0 0 Nowlookingattheconstraintsthesecondrevolutejointimposesontheangulardegreesoffreedombetweenbody1andbody2 5=b whichcanbewrittenintermsoftheinertialcoordinatesystemas: 5=b whereb 7{62 is: Similarly,thesecondconstraintonrotationaldegreesoffreedomimposedbytherevolutejointbetweenbody1andbody2isillustratedin: 6=b 84


7{64 intermsoftheinertialcoordinatesystemproduces: 6=b TakingthepartialderivativeofEquation 7{65 withrespecttothegeneralizedcoordinatesyields: Theconstraintsontherelativetranslationalmotionbetweenbody2andbody3causedbythethirdrevolutejoint,areillustratedin (7{67) whereR 7{67 intermsoftheinertialcoordinatesystemyields: (7{68) 85


7{69 withrespecttothegeneralizedcoordinatesproduces: 0 0 Lookingattheconstraintsontherelativerotationalmotionbetweenbody2andbody3imposedbytherevolutejoint: 8=b whereb 8=b ThecontributiontotheJacobianmatrixfromthisconstraintis: Similarly,thesecondrotationalconstraintimposedbytherevolutejointbetweenbody2andbody3isexempliedin: 9=b 86


7{73 intermsoftheinertialsystemyields: 9=b TakingthepartialderivativeofEquation 7{74 withrespecttothegeneralizedcoordinatesproduces: UsingthesystemJacobianmatrixandthesystemconstraintmatrix,derivedintheabove,thesystemofequationsillustratedinEquation 7{48 arereadyforintegration.Thereareseveralmethodstointegratetheseequations.Thesimplestinvolvesdirectintegration.Dene [M]=264[M][]Tq Thus,Equation 7{48 iswrittenas: [M]q 87


whichgivesthesimultaneoussolutionforthesecondtimederivativeofthegeneralizedcoordinatesandtheLagrangemultipliers.Thesolutionofthesecondtimederivativeofthegeneralizedcoordinatescanthenbeusedinarstordersystemtosolveforthegeneralizedcoordinatesatthenexttimestep.Towritetheequationsofmotionasarstordersystem,rstlet: Thus, _X whereu 7-1 demonstratesthesolutionoftheequationsofmotiontosimulatethesystemafteronesecondinwhichacontrolisappliedto1.Thereisaknownproblemwiththedirectintegrationoftheequationsofmotion.ThisproblemisthatthesolutionoftheLagrangemultipliersisnotexact.Asthesimulationtimegrows,theseerrorsalsogrow.Asseeninthegure,thenormofthevectorthatrepresentstheconstraintvectorforthesystemgrowsasthesimulationprogresses.ThereisastabilizationmethodthatreducestheerrorsassociatedwiththesolutionoftheaugmentedsystemillustratedinEquation 7{79 .ThisprocesswasrstpublishedintheworkbyBaumgarte[ +2_ +2 =0 (7{83) 88


7{83 lieintheleft-halfplane,thesolutionisguaranteedtoconvergetozero.Thus,withtheproperchoiceofandinEquation 7{83 theconstraintvectorandtherstandsecondtimederivativesoftheconstraintvectorcanbemadetoconvergetozero.[ 7-2 showsthecomparisonbetweenthenormoftheconstraintvectorwithoutstabilizationandwithstabilization.Asillustratedinthegure,withoutconstraintstabilizationtheerrorintheconstraintvectorgrowsexponentiallywithtime.Stabilizationhelpstoreducetheerrorsaccumulatedandtherateatwhichtheyaccumulate.However,asillustratedthesimultaneoussolutionoftheLagrangemultipliersandtheaccelerationsthatisnecessarywhenusingthestandardformulationoftheequationsofmotionstillleadstoviolationsoftheconstraintsonthesystem.Forthisreason,formingtheequationsofmotionwithouttheneedtocoupleconstraintsequationsisbenecialintheoverallsimulationofthesystem. 3 allowsfortheeliminationofsomeofthesediculties.Thenumberofgeneralizedcoordinatesismuchlesswiththisalternativeformulation.Sinceallbodiesinthesystemarereferencedtoa\central"body,onlythefullcoordinatesassociatedwiththisbodyarenecessary.Thosecoordinatesarethenaugmentedwithcoordinatesthatrelatethedierencesinorientationbetweensubsequentbodiesandthe\central"body. Thederivationoftheequationsofmotion,usingthealternativeformulation,forthethreebarsystembeginswiththeestablishmentofcoordinateframes.Asinthestandard 89


wheree [C]=266664cos1sin10sin1cos10001377775(7{85) InthealternativeformulationillustratedinChapter 3 showedjointframesthatwereusedasthebasisforrelatingthebody-xedcoordinatesystemofonebodytothebody-xedcoordinatesystemofanadjoiningbody.Forthisexamplethesejointframesarenotnecessaryasthebody-xedcoordinateframesarealignedwiththejointaxisoftherevolutejointsforallbodies.Thustherelationshipbetweentherstbodyandthesecondbodycanbestatedas: where [T(2)]=2666641000cos2sin20sin2cos2377775(7{87) andr 90


where [T(3)]=2666641000cos3sin30sin3cos3377775(7{89) andr whereb whereb 91


whereRCB1iswrittenintermsoftheinertialcoordinatesystemandu 7{93 ,withrespecttotheinertialcoordinateframe,yields: dtr dt(R dtu buttherevolutejointbetweenthegroundandbody1preventstherigid-bodytranslationofbody1,thusR where[0~!1]istheskewsymmetricmatrixformedfromE! Thus, Thusthecontributiontothekineticenergyfrombody1canbewrittenas: 2ZVv 2ZVu 92


2_qM_q(7{99) where_q=_1andusingtheindicialnotation, Thederivationofthecontributiontothekineticenergyofthesystemfrombody2beginswiththedenitionofthevectorfromtheoriginoftheinertialsystemtoarandompointonbody2.Thisvectoriswrittenas: ,wherer 7{101 withrespecttotheinertialcoordinateframeyields: dtr dt(R dtr dtr dtu but 93


Theskewsymmetricmatricescanbefactoredtoyield: Let [B]=B Thusthevelocityofarandompointonbody2iswrittenas: InsertingEquation 7{109 intothecontributiontothesystem'skineticenergyfrombody2yields: 2ZV[B]T[B]dV8><>:_1_29>=>;(7{110) Thus,thecontributiontothekineticenergyfrombody2assumesthefamiliarquadraticform: 2_q 94


[M]=ZV[B]T[B]dV=264RVB LookingatarepresentativetermofEquation 7{112 andusingtheindicialnotation, ItisimportanttonotethattheintegralsleftinEquation 7{113 areconstantsintegrals. Thederivationofthecontributiontothesystem'skineticenergyfrombody3beginswiththeformationofapositionvectorfromtheoriginoftheinertialsystemtoarandompointonbody3.Thisvectorisstatedas: wherer 95


7{115 withrespecttotheinertialcoordinatesystemyields: dtr dt(R dtr dtr dtr dtr dtu However, InsertingEquation 7{116 andEquation 7{117 intoEquation 7{115 andfactoringtheskewsymmetricmatricesproduces: Dene 96


Therefore, InsertingEquation 7{122 intothecontributiontothesystem'skineticenergyfrombody3yieldsthefamiliarquadraticform: 2_q where _q=8>>>><>>>>:_1_2_39>>>>=>>>>;(7{124) [M]=266664RVBB ExaminingarepresentativetermofEquation 7{125 andusingtheindicialnotationyields: 97


Hereagain,theremainingintegralsinEquation 7{126 areinvariantoverthebody. Combiningthekineticenergycontributionsfrombody1,body2,andbody3yieldsthekineticenergyoftheentiresystem.Then,insertingthiskineticenergytermintotheLagrange'sequationsproducestheequationsofmotionofthesystemwhichcanbestated 98


[Msys]q where dt([Msys])+1 2@ @q Thealternativeformulationoftheequationsofmotionillustratedintheaboveresultsinasystemthathasonlythreegeneralizedcoordinates.Thesecoordinatesaretheanglesthatrepresenttherelativemotionofthebodiesinthesystemwithrespecttoeachother.Inthestandardformulation,thenumberofgeneralizedcoordinatesisgreatlyincreasedduetothenecessitytocompletelydescribethemotionofeachbodyinsteadoftherelativemotionofthebodieswithrespecttoeachotherwhichistheresultofthealternativeformulation.Thereducednumberofgeneralizedcoordinateshastheprimaryadvantageofthereductionincomputationalcostwhenintegratingtheequationsofmotion.Asmorebodiescomprisethesystem,thereductionincomputationalcostbyusingthealternativeformulationoftheequationsofmotionbecomesmoreevident.Alsoincludedinthecomputationalsavingsthealternativeformulationdemonstratesistheeaseofintegrationoftheequationsofmotionincomparisontotheintegrationoftheequationsofmotionforthestandardformulation.TheequationsofmotionasderivedusingthealternativeformulationcanbedirectlyintegratedwithouthavingtorstsimultaneouslysolveforLagrangeconstraints,asisnecessaryinthestandardformulationoftheequationsofmotion.However,thestandardformulationhastheadvantageofeachnewbodyis 99

PAGE 100

Theequationsofmotionmustnowbeputintoaformamenabletotheoptimalcontrolproblem.Thustheequationsofmotionbecome: _X ;u where ;u 100

PAGE 101

minJ=tfRt0u ;u Forthisexamplethenaltimeisonesecond.TheoptimalcontrolproblemdenedintheabovewastranscribedintoaNon-LinearProgram(NLP)usingaGaussPseudospectralmethod.Figures 7-4 and 7-5 showthetimehistoryoftheanglesandtheangularratesundertheinuenceoftheoptimalcontrolillustratedingure 7-6 .Thesegureswereproducedbytheoptimizationsoftware,GPOCS,aftertheoptimizationwasterminatedbecauseofalackofimprovementintheminimizationofthecostfunctionoverpreviousiterations.Althoughtheequationsofmotionwereformulatedinawaythatwouldallowforthesolutionoftheoptimalcontrolproblem,theresultsdonotdemonstrateasuccessfulsolution.AsillustratedinFigure 7-7 andFigure 7-8 ,usingthecontrolgeneratedbytheoptimizationsoftware,theintegratedequationsofmotiondonotreachthedesirednalangles.Moreover,theangularratesattheterminaltimearenonzero.SincetheGaussPseudospectralmethoddiscretizestheproblemintime,itisverydependentonthenumberofdiscretetimesatwhichthesolutionstrategyseeksanoptimum.Anincreaseinthenumberofdiscretizationpointswouldseemtoleadtoagreaterprobabilityofndinganoptimum.Theequationsofmotionastheyarewritten,however,requireevaluationateachdiscretizationpointusinginformationfromthepreviousdiscretizationpoint.Thusseekingasimultaneoussolutionatallthediscretizationpointsissusceptibletonumerical 101

PAGE 102

7{132 ,iscomputationallyintensive.Theoptimizationsoftwareattemptstouseaglobalpolynomialtoconvertthecontinuousequationsofmotionintoadiscretesetofequations.Ateachpointinthediscretization,thestatesandthecontrolarebeingsolved.Thus,thenumberofdiscretizationpointscoupledwiththenatureoftheevaluationoftheequationsofmotionnumericaldicultiesarepossible.Theoptimizationsoftwareusedinthisworkwasnotabletoovercomethesenumericalinstabilitiesandwasunabletoconvergetotheoptimumsolution. 102

PAGE 103

Normoftheconstraintvector Figure7-2: Normoftheconstraintvector 103

PAGE 104

Threebarmanipulator. 104

PAGE 105

TimehistoryoftheanglesasdeterminedbyGPOCS. Figure7-5: TimehistoryoftheangularratesasdeterminedbyGPOCS. 105

PAGE 106

Timehistoryofthecontroltorques. Figure7-7: Timehistoryoftheanglesasdeterminedbyintegration. 106

PAGE 107

Timehistoryoftheangularratesasdeterminedbyintegration. 107

PAGE 108

ThisresearchwasconcernedwiththeemergenceofMicroAirVehiclesandthevariousnewopportunitiesthattheyprovide.OneofthemostcompellingistheabilitytousemorphingonMAVs.ThesmallstructureandlightweightmaterialsallowsformorphingstrategiesthatcancompletelyalterthecongurationofMAVs.Thesecongurationchangescanthenbeusedtoexploitightregimesandmaneuversthatarealientoconventionalaircraft.However,withthesenewopportunitiescomenewobstaclesthatmustbeovercome.ThechiefamongtheseobstaclesisthelackofawaytocharacterizetheperformanceofMAVsbasedonconventionalaircraftmodeling.Newmodelingtechniquesarerequiredthatallowforthe,insomecases,completerecongurationofthevehiclebecauseintheightregimethatMAVstypicallyoperateinthechangeshavegreatimpactontheperformance.Inanattempttocreatenewmodelingtechniques,thisstudyoutlinedthederivationofequationsofmotionbasedonmulti-bodydynamicsthatallowsforthelargecongurationchangespossiblewithmorphing.Theseequationsofmotionalsoincorporatethecouplingofrigidbodyandexiblebodydynamicsintoasinglesetofequations.Alsopresentedwastheabilitytotaketheequationsofmotionasderivedandcastthemintoaformamenableforanoptimalcontrolproblem.Thenumericalexampledemonstratedtheimplementationofthemulti-bodyequationofmotionformulationonamodelthatsimulatesthecongurationchangespossibleonMAVs.Thethree-barexampleshowsthatthemethodsderivedareapplicabletotraditionaldynamicsystemsaswellandthroughtheoptimalcontrolproblem,allowsfornewcontrolstrategies. However,therearechallengesintroducebythisresearchthatneedtobeovercomeasthetechniquesmature.PrimarilyamongthesearethenumericaldicultiesassociatedwiththesolutionoftheOptimalControlProblem.Withsystemsthatincorporaterigidandexiblemodes,therearealargenumberofdegreesoffreedom.Itiswellknownthatasthedegreesoffreedomincrease,thesolutiontimefordirectmethodsalsoincreases.As, 108

PAGE 109

Thisresearchoersmanyavenuesforfurtherresearchinthemanyareasinvolvedinitsspan.Forexample,thecontinuedapplicationofamultibodyformulationprovidesmanyopportunitiesforfurtherresearch.Oneofthebenetsofthemultibodyformulationisthatitallowseachbodytobemodeledwithadierentlevelofdelity.Thisfactisespeciallyimportantwhenconsideringexiblestructures.TheFiniteElementMethodisthemostlogicalchoicefortherepresentationofexibledegreesoffreedom.Researchstillneedstobedoneintohowvariouslevelsofdelityinmodelingtheexibledegreesoffreedomaectstheoverallsimulationofthedynamicsofthesystem.Themultibodyformulationoftheequationsofmotion,asillustratedinthisresearch,leadstohighlynonlinearequationsthathaveastrongdependenceonthestateofthesystem.Theequationspresentedinthiswork,however,didnotincludethenonlineareectsthatexistinrealsystems.Anotherresearchopportunityexistsintheinvestigationofeectssuchasjointfrictionandhysteresisincontrolactuatorshasonthegoverningequations.Theseresearchtopicswouldbenecessaryforapplicationsofthisresearchintotheeldofrobotics,forwhichitiswell-suited.Perhapstherichestresearchtopicthatextendsfortheworkpresentedhereliesintheareaofoptimalcontrol.Apossibletopicforresearchliesintheuseofsomesymbolicmanipulationtocastthehighlystatedependentequationsofmotionintoaformthatallowsfortheanalyticalsolutionofthecostateequationsandthusanapplicationofindirectmethods.AsdemonstratedinChapter 6 andChapter 7 directmethodsdonotguaranteethesuccessfulsolutionoftheoptimalcontrolproblem.Moreresearchcanbedonetoeithertransformtheequationsofmotion,asderivedinthisresearch,intoaformthattheoptimizationsoftwareGPOCScansuccessfullyndanoptimumforortondadirectmethodthatismoreappropriate. TherearealsotopicsforfurtherresearchintheapplicationofthisworktoMAVsinspecic.Determiningtheinteractionbetweenthestructuresandtheuidthroughwhich 109

PAGE 110


PAGE 111

Thissectiondescribestheformationofthemassmatrixforbody2.Theuseofindicialnotationwillfacilitatethederivationoftermstobeintegrated.Themassmatrixforbody2is: 226666664[I][C1][B][C1]D 2[I]dV=ZV1 2I;mndV(A{2) [C1][B]=[C1]B Thusexaminingarepresentativeterm, whereRVdV;RVuop;pdV;andRVN2;npdVareinvariants.ForthetermsRV[C1]B 111

PAGE 112


PAGE 113

Examiningarepresentativeterm, 113

PAGE 114

Examiningarepresentativeterm, Fortheremainingterms,replace0~!1;with0~!1;and0~!1;whereappropriate. dV=rH1CB2;j21;jk12~!21;mk12~!21;mn21;pnrH1CB2;pZVdV+rH1CB2;j21;jk12~!21;mk12~!21;mn21;pnZVuop;pdV+rH1CB2;i21;ij12~!21;kj12~!21;km21;nmq2f;pZVN2;npdV+21;jk12~!21;mk12~!21;mn21;pnrH1CB2;pZVuop;jdV+21;jk12~!21;mk12~!21;mn21;pnZVuop;juop;pdV+21;ij12~!21;kj12~!21;km21;nmq2f;pZVuop;iN2;npdV+q2f;i21;jk12~!21;mk12~!21;mn21;pnrH1CB2;pZVN2;jidV+q2f;i21;jk12~!21;mk12~!21;mn21;pnZVN2;jiuop;pdV+q2f;s21;ij12~!21;kj12~!21;km21;nmq2f;pZVN2;isN2;npdV 114

PAGE 115

Finally, Itisimportanttonotethatallintegralsleftwhenwrittenintheindicialnotationareconstantintegralsandonlyneedbeperformedonce. 115

PAGE 116

Thissectiondetailstheformationofthemassmatrixforbody3.Themassmatrixforbody3isasfollows: 226666666664[I][C1][BB][C1]DD where[I]istheidentitymatrix. Examiningarepresentativeterm, Fortheremainingtermsreplace0~!1;with0~!1;and0~!1;whereappropriate. 116

PAGE 117


PAGE 119


PAGE 120


PAGE 121


PAGE 122


PAGE 123


PAGE 124

Finally, Notealloftheintegralsleftintheindicialnotationareconstant. 124

PAGE 125

[1] Shabana,A.A.,DynamicsofMultibodySystemsNewYorkNY:JohnWiley&Sons,1989. [2] Meirovitch,L.andStemple,T.,\HybridStateEquationsofMotionforFlexibleBodiesinTermsofQuasi-Coordinates",JournalofGuidance,Control,andDynam-ics,Vol.18,NO.4,1995,pp.678-688. [3] Meirovitch,L.andTuzcu,I.,\IntegratedApproachtotheDynamicsandControlofManeuveringFlexibleAircraft",NASATechnicalPaper,PaperNo.NASA/CR-2003-211748,2003. [4] Meirovitch,L.andTuzcu,I.,\UniedTheoryfortheDynamicsandControlofManeuveringFlexibleAircraft",AIAAJournal,Vol.42,No.4,April2004,pp.714-727. [5] Meirovitch,L.andTuzcu,I.,\TimeSimulationsoftheResponseofManeuveringFlexibleAircraft",JournalofGuidance,Control,andDynamics,Vol.27,No.5,2004,pp.814-828. [6] Kwak,M.K.andMeirovitch,L.,\NewApproachtotheManeuveringandControlofFlexibleMultibodySystems",JournalofGuidance,Control,andDynamics,Vol.15,No.6,1992,pp.1342-1353. [7] Gad-el-Hak,M.,\Micro-Air-Vehicles:CanTheybeControlledBetter?",JournalofAircraft,Vol.38,No.3,2001,pp.419-429. [8] Lian,Y.,ShyyW.,Viieru,D.,andZhang,B.,\MembraneWingAerodynamicsforMicroAirVehicles",ProgressinAerospaceSciences,Vol.39,2003,pp.425-465. [9] Mueller,T.J.andDeLaurier,J.D.,\AerodynamicsofSmallVehicles",AnnualReviewofFluidMechanics,Vol.35,2003,pp.89-111. [10] MutiLin,J.C.andPauley,L.L.,\Low-Reynolds-NumberSeparationonanAirfoil",AIAAJournal,Vol.34,No.8,1996,pp.1570-1577. [11] Lutz,T.,Wurz,W.,andWagner,S.,\NumericalOptimizationandWind-TunnelTestingofLowReynolds-NumberAirfoils",ConferenceonFixed,FlappingandRotaryWingVehiclesatveryLowReynoldsNumbersNotreDame,Indiana2000. [12] Ifju,P.G.,Ettinger,S.,Jenkins,D.A.,andMartinez,L.,\CompositeMaterialsforMicroAirVehicles",SocietyfortheAdvancementofMaterialsandProcessEngineeringAnnualConference,May2001. [13] Ifju,P.G.,Jenkins,D.A.,Ettinger,S.,Lian,Y.,Shyy,W.,andWaszak,M.R.,\Flexible-Wing-BasedMicroAirVehicles"AmericanInstituteofAeronauticsandAstronauticsPaperNo.AIAA2002-0705,2002. 125

PAGE 126

Hall,J.,Mohseni,K.,Lawrence,D.,andGeuzaine,P.,\InvestigationofVariableWing-SweepforApplicationsinMicroAirVehicles",AmericanInstituteofAeronauticsandAstronauticsPaperNo.AIAA2005-7171,2005. [15] Ramamurti,R.andSandberg,W.,\SimulationoftheDynamicsofMicro-AirVehicles",38thAerospaceSciencesMeeting&ExhibitReno,Nevada,PaperNo.AIAA2000-0896,2000. [16] Raney,D.L.andSlominski,E.C.,\MechanizationandControlConceptsforBiologicallyInspiredMicroAirVehicles",JournalofAircraft,Vol.41,No.6,2004,pp.1257-1265. [17] Schenato,L.,Campolo,D.,andSastry,S.,\ControllabilityIssuesinFlappingFlightforBiomimeticMicroAirVehicles(MAVs)",ProceedingsofIEEEConferenceonDecisionandControl,2003. [18] Abdulrahim,M.,Garcia,H.,Ivey,G.F.,andLind,R.,\FlightTestingaMicroAirVehicleUsingMorphingforAeroservoelasticControl",AIAAStructures,StructuralDynamics,andMaterialsConference,PalmSprings,CA,PaperNo.AIAA2004-1674,2004. [19] Waszak,M.R.,Davidson,J.B.,andIfju,P.G.,\SimulationandFlightControlofanAeroelasticFixedWingMicroAirVehicle",AIAAAtmosphericFlightMechanicsConference,Monterey,CA,PaperNo.AIAA2002-4875,2002. [20] Waszak,M.R.,Jenkins,L.N.,andIfju,P.,\StabilityandControlPropertiesofanAeroelasticFixedWingMicroAirVehicle",AIAAAtmosphericFlightMechanicsConference,Montreal,Canada,PaperNo.AIAA2001-4005,2005. [21] Lian,Y.,Shyy,W.,andHaftka,R.T.,\ShapeOptimizationofaMembraneWingforMicroAirVehicles",AIAAJournal,Vol.42,No.2,2003,pp.424-426. [22] Lian,Y.,Shyy,W.,Ifju,P.G.,andVerron,E.,\MembraneWingModelofMicroAirVehicles",AIAAJournal,Vol.41,No.12,2003,pp.2492-2494. [23] Lian,Y.andShyy,W.,\NumericalSimulationsofMembraneWingAerodynamicsforMicroAirVehicleApplications",JournalofAircraft,Vol.42,No.4,2005,pp.865-873. [24] Ramrakhyani,D.S.,Lesieutre,G.A.,Frecker,M.,andBharti,S.,\AircraftStructuralMorphingUsingTendon-ActuatedCompliantCellularTrusses",JournalofAircraft,Vol.42,No.6,2005,pp.1615-1621. [25] Bae,J.,Seigler,T.M.,andInman,D.J.,\AerodynamicandStaticAeroelasticCharacteristicsofVariable-SpanMorphing",JournalofAircraft,Vol.42,No.2,2005,pp.528-534. 126

PAGE 127

Livne,E.andWeisshaar,T.A.,\AeroelasticityofNonConventionalAirplaneCongurations{PastandFuture",JournalofAircraft,Vol.40,No.6,2003,pp.1047-1065. [27] Wlezien,R.W.,Horner,G.C.,McGowan,A.R.,Padula,S.L.,Scott,M.A.,Silcox,R.J.,andSimpson,J.O.,\TheAircraftMorphingProgram",AmericanInstituteofAeronauticsandAstronauticsPaperNo.AIAA-98-1927,1998. [28] Reich,G.W.,Raveh,D.E.,andZink,P.S.,\ApplicationofActive-Aeroelastic-WingTechnologytoaJoined-WingSensorcraft",JournalofAircraft,Vol41,No.3,2004,pp.594-602. [29] Zheng,L.andRamaprian,B.R.,\APiezo-ElectricallyActuatedWingforaMicroAirVehicle",AmericanInstituteofAeronauticsandAstronauticsPaperNo.AIAA2002-2302,2002. [30] Saunders,B.,Eastep,F.E.,andForster,E.,\AerodynamicandAeroelasticCharacteristicsofWingswithConformalControlSurfaces",JournalofAircraft,Vol.40,No.1,2003,pp.94-99. [31] Abdulrahim,M.,Garcia,H.,andLind,R.,\FlightCharacteristicsofShapingtheMembraneWingofaMicroAirVehicle",JournalofAircraft,Vol.42,No.1,2005,pp.131-137. [32] Bisplingho,R.L.andAshley,H.,PrinciplesofAeroelasticity,NewYork,NY:DoverPublications,Inc.1975. [33] Bisplingho,R.L.,Ashley,H.,andHalfman,R.L.,AeroelasticityMineola,NY:DoverPublications,Inc.1996. [34] Friedmann,P.P.,\RenaissanceofAeroelasticityandItsFuture",JournalofAircraft,Vol.36,No.1,1999,pp.105-121. [35] Kral,R.andKreuzer,E.,\MultibodySystemsandFluid-StructureInteractionswithApplicationtoFloatingStructures",MultibodySystemDynamics,Vol.8,1999,pp.65-83. [36] Shabana,A.A.,\FlexibleMultibodyDynamics:ReviewofPastandRecentDevelopment",MultibodySystemDynamics,Vol.1,1997,pp.189-222. [37] Schiehlen,W.,\ResearchTrendsinMultibodySystemDynamics",MultibodySystemDynamics,Vol.18,2007,pp.3-13. [38] Yoo,W.,Kim,K.,Kim,H.,andSohn,J.,\DevelopmentsofMultibodySystemDynamics:ComputerSimulationsandExperiments",MultibodySystemDynamics,Vol.18,2007,pp.35-58. 127

PAGE 128

Boutin,B.A.andMisra,A.K.,\DynamicsofManipulatingTrussStructures",CollectionofTechnicalPapers,AIAA/AASAstrodynamicsConference,SanDiego,CA,PaperNo.A96-34712-09-12,July1996,pp.471-4793. [40] deJalon,J.G.,\Twenty-FiveYearsofNaturalCoordinates",MultibodySystemDynamics,Vol.18,2007,pp.15-33. [41] Leger,M.andMcPhee,J.,\SelectionofModelingCoordinatesforForwardDynamicMultibodySimulations",MultibodySystemDynamics,Vol.18,2007,pp.277-297. [42] Anderson,K.S.,andDuan,S.,\HighlyParallelizableLow-orderDynamicsSimulationAlgorithmforMulti-Rigid-BodySystems",JournalofGuidance,Control,andDynamics,Vol.23,No.2,2000,pp.355-364. [43] Banerjee,A.K.,\ContributionsofMultibodyDynamicstoSpaceFlight:ABriefReview",JournalofGuidance,Control,andDynamics,Vol.26,No.3,2003,pp.385-394. [44] ChenS.,HansenJ.M.,Tortorelli,D.A.,\UnconditionallyEnergyStableImplicitTimeIntegration:ApplicationtoMultibodySystemAnalysisandDesign",Interna-tionalJournalforNumericalMethodsinEngineering,Vol.48,2000,pp.791-822. [45] Chen,S.,Tortorelli,D.A.,\OptimizationofMultibodySystemsViaanEnergyConservingMinimalCoordinateFormulation",AmericanInstituteofAeronauticsandAstronautics,PaperNo.AIAA-98-4927,1998. [46] Schaefer,E.,Loesch,M.,Rulka,W.,\Real-TiemSimulationofFlexibleSpaceManipulatorDynamics",AmericanInstituteofAeronauticsandAstronautics,PaperNo.AIAA97-37062,1997. [47] Kielau,G.andMaiber,P.,\NonholonomicMultibodyDynamics",MultibodySystemDynamics,Vol.9,2003,pp.213-236. [48] Aoustin,Y.andFormalsky,A.,\OntheFeedforwardTorquesandReferenceTrajectoryforFlexibleTwo-LinkArm",MultibodySystemDynamics,Vol.3,1999,pp.241-265. [49] Junkins,J.L.,\AdventuresontheInterfaceofDynamicsandControl",AIAA,35thAerospaceSciencesMeeting&Exhibit,Reno,NV,Jan.6-9,1997. [50] Rives,J.,Portigliotti,S.,andLeveugle,T.,\AtmosphericRe-EntryDemonstratorDescentandRecoverySub-SystemQualicationTestFlightDynamicsandTrajectoryModelingduringDescent",AmericanInstitueofAeronauticsandAs-tronautics,PaperNo.AIAA-97-1441,1997. [51] Tadikonda,S.S.K.,\ModelingofTranslationalMotionBetweenTwoFlexibleBodiesConnectedviaThreePoints",JournalofGuidance,Control,andDynamics,Vol.18,No.6,1995,pp.1392-1397. 128

PAGE 129

Tadikonda,S.S.K.,\Articulated,FlexibleMultibodyDynamicsModeling:GeostationaryOperationalEnvironmentalSatellite-8CaseStudy",JournalofGuidance,Control,andDynamcis,Vol.20,No.2,1997,pp.276-283. [53] Jain,A.,andRodriguez,G.,\RecursiveDynamicsAlgorithmforMultibodySystemswithPrescribedMotion",JournalofGuidance,Control,andDynamics,Vol.16,No.5,1993,pp.830-837. [54] Pradhan,S.,Modi,V.J.,andMisra,A.K.,\OrderNFormulationforFlexibleMultibodySystemsinTreeTopology:LagrangianApproach",JournalofGuidance,Control,andDynamics,Vol.20,No.4,1997,pp.665-672. [55] Attia,H.,\DynamicSimulationofConstrainedMechanicalSystemsUsingRecursiveProjectionAlgorithm",JournalofBrazilianSocietyofMechanicalSciencesandEngineering,Vol.28,2006,pp.37-44. [56] Dias,J.M.P.andPereira,M.S.,\SensitivityAnalysisofRigid-FlexibleMultibodySystems",MultibodySystemDynamics,Vol.1,1977,pp.303-322. [57] Ortiz,J.L.andBarhorst,A.A.,\ModelingFluid-StructureInteraction",JournalofGuidance,Control,andDynamics,Vol.20,No.6,1997,pp.1221-1228. [58] Paiewonsky,B.,\OptimalControl:AReviewofTheoryandPractice",AIAAJournal,Vol.3,No.11,1965,pp.1985-2006. [59] Hodges,D.H.andBless,R.R.,\WeakHamiltonianFiniteElementMethodforOptimalControlProblems",JournalofGuidance,Control,andDynamics,Vol.14,No.1,1991,pp.148-156. [60] Schaub,H.andJunkins,J.L.,\NewPenaltyFunctionsandOptimalControlFormulationforSpacecraftAttitudeControlProblems",JournalofGuidance,Control,andDynamics,Vol.20,No.3,1997,pp.428-434. [61] Seywald,H.andKumar,R.R.,\FiniteDierenceSchemeforAutomaticCostateCalculation",JournalofGuidance,Control,andDynamics,Vol.19,No.1,1996,pp.231-239. [62] Seywald,H.andKumar,R.R.,\AMethodforAutomaticCostateCalculation",AIAA,Guidance,Navigation,andControlConference,SanDiego,CA,PaperNo.AIAA-96-3699,1996. [63] Seywald,H.andKumar,R.R.,\MethodforAutomaticCostateCalculation",JournalofGuidance,Control,andDynamics,Vol.19,No.6,1996,pp.1252-1261. [64] Venkataraman,P.,\DynamicsandOptimalControl",AIAAGuidance,Navigation,andControlConferenceandExhibit,Dever,CO,PaperNo.AIAA-2000-4560,2000. 129

PAGE 130

Vadali,S.R.andSharma,R.,\OptimalFinite-TimeFeedbackControllersforNonlinearSystemswithTerminalConstraints",JournalofGuidance,Control,andDynamics,Vol.29,No.4,2006,pp.921-928. [66] Betts,J.T.,\SurveyofNumericalMethodsforTrajectoryOptimization",JournalofGuidance,Control,andDynamics,Vol.21,No.2,1998,pp.193-207. [67] Hull,D.G.,\ConversionofOptimalControlProblemsintoParameterOptimizationProblems",JournalofGuidance,Control,andDynamics,Vol.20,No.1,1997,pp.57-60. [68] Thompson,R.C.,Junkins,J.L.,andVadali,S.R.,\Near-MinimumTimeOpen-LoopSlewingofFlexibleVehicles",JournalofGuidance,Control,andDynamics,Vol.12,No.1,1989,pp.82-88. [69] Warner,M.S.,andHodges,D.H.,\TreatmentofControlConstraintsinFiniteElementSolutionofOptimalControlProblems",JournalofGuidance,Control,andDynamics,Vol.22,No.2,1999,pp.358-360. [70] Bass,M.,vonStryk,O.,Bulirsch,R.,andSchmidt,G.,\TowardsHybridOptimalControl",at-Automatisierungstechnik,Vol.48,No.9,pp.448-459,2000. [71] vonStryk,O.,andBulirsch,R.,\DirectandIndirectMethodsforTrajectoryOptimization",AnnalsofOperationsResearch,Vol.37,pp.357-373,1992. [72] Hu,G.S.,Ong,C.J.,andTeo,C.L.,\DirectCollocationandNonlinearProgrammingforOptimalControlProblemUsinganEnhancedTranscribingScheme",Proceedingofthe1999IEEEInternationalSymposiumonComputerAidedControlSystemDesign,KohalaCoast-IslandofHawaii,Hawaii,USA,1999. [73] Frazzoli,E.andBulle,F.,\OnQuantizationandOptimalControlofDynamcalSystemswithSymmetries",Proceedingofthe41stIEEEConferenceonDecisionandControl,LasVegas,NV,2002,pp.817-823. [74] Huntington,G.,Benson,D.,andRao,A.,\AComparisonofAccuracyandComputationalEciencyofThreePseudospectralMethods"2007AIAANavi-gationandControlConverenceHiltonHead,SC,2007,AIAAPaper2007-6405. [75] Huntington,G.andRao,A.,\AComparisonofLocalandGlobalOrthogonalCollocationMethodsforSolvingOptimalControlProblems",2007AmericanControlConference,NewYork,NY,2007. [76] Williams,P.,\AGauss-LobattoQuadratureMethodforSolvingOptimalControlProblems",AustralianandNewZealandIndustrialandAppliedMathematicsJournal,Vol.47,2006,pp.C101-C115. [77] Fahroo,F.andRoss,I.M.,\CostateEstimationbyaLegendrePseudospectralMethod",JournalofGuidance,Control,andDynamics,Vol.24,No.2,2001,pp.270-277. 130

PAGE 131

Herman,A.L.,andConway,B.A.,\DirectOptimizationUsingCollocationBasedonHigh-OrderGauss-LobattoQuadratureRules",JournalofGuidance,Control,andDynamics,Vol.19,No.3,1996,pp.592-599. [79] Lo,J.,Huang,G.,andMetaxas,D.,\HumanMotionPlanningBasedonRecursiveDynamicsandOptimalControlTheory",MultibodySystemDynamics,Vol.1,1997,pp.189-222. [80] Huntington,G.andRao,A.,\OptimalCongurationofSpacecraftFormationsViaaGaussPseudospectralMethod",15thSpaceFlightMechanicsMeeting,CopperMountain,CO,2005. [81] Huntington,G.andRao,A.,\OptimalRecongurationofaTetrahedralFormationViaaGaussPseudospectralMethod",2005AstrodynamicsSpecialistConference,LakeTahoe,CA,2005. [82] Baumgarte,J.,\StabilizationofConstraintsonIntegralsofMotion",ComputerMethodsinAppliedMechanicsandEngineering,Vol.1,1972,pp.1-16. [83] Weijia,Z.,Zhenkuan,P.,andYibing,W.,\AnAutomaticConstrainViolationStabilizationMethodforDierential/AlgebraicEquationsofMotioninMultibodySystemDynamics",AppliedMathematicsandMechanics,Vol.21,2000,pp.103-108. [84] Ascher,U.,Chin,H.,Petzold,L.,andReich,S.,\StabilizationofConstrainedMechanicalSystemswihtDAEsandInvariantManifolds",JournalofMechanicsofStructuresandMachines,Vol.23,1995,pp.135-158. [85] Lin,S.andHong,M.,\StabilizationMethodfortheNumericalIntegrationofControlledMultibodyMechanicalSystem:AHybridIntegrationApproach",TheJapanSocietyofMechanicalEngineersInternationalJournal,Vol.44,2001,pp.79-88. [86] Thomas,S.,\DynamicsofSpacecraftandManipulators",SIMULATION,Vol.57,1991,pp.56-72. [87] Lin,S.andHuang,J.,\StabilizationofBaumgarte'sMethodUsingtheRunge-KuttaApproach",JournalofMechanicalDesign,Vol.124,2002,pp.633-641. [88] Anderson,K.andCritchley,J.,\Improved'Order-N'PerformanceAlgorithmfortheSimulationofConstrainedMulti-Rigid-BodyDynamicSystems",MultibodySystemsDynamics,Vol.9,2003,pp.185-212. [89] Simeon,B.\OnLagrangeMultipliersinFlexibleMultibodyDynamics",ComputerMethodsAppliedMechanicalEngineering,Vol.195,2006,pp.6993-7005. [90] Bottasso,C.L.andCroce,A.,\OntheSolutionofInverseDynamicsandTrajectoryOptimizationProblemsforMultibodySystems",MultibodySystemDynamics,Vol.11,2004,pp.1-22. 131

PAGE 132

Bottasso,C.L.andCroce,A.,\OptimalControlofMultibodySystemsUsinganEnergyPreservingDirectTranscriptionMethod",MultibodySystemDynamics,Vol.12,2004,pp.17-45. [92] Benson,D.,\AGaussPseudospectralTranscriptionforOptimalControl",Ph.D.Thesis,MassachusttsInstituteofTechnology2004. 132

PAGE 133

SeanJamesSidneyRegisfordwasborninSt.Georges,Grenada,onJuly25,1977.In1979heandhisfamilymovedtoFlorida.HegraduatedfromMiramarHighSchoolin1995andthenattendedtheUniversityofFlorida,studyingaerospaceengineering.In2000,hegraduatedfromtheUniversityofFloridawithabachelor'sdegreeinaerospaceengineeringandthenwasacceptedtothegraduateprogramattheUniversityofFlorida.In2002,hereceivedhisMasterofSciencedegreeinaerospaceengineeringandthencontinuedhisstudiesinthePh.D.program.UponcompletionofthedegreeDoctorofPhilosophy,hehopestoobtainemploymentinindustry 133