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,zyH2CB3, ,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(B12)
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]
(B13)
Examining a representative term:
/ BBEf [a12] [T (a)] [a1] [a3] [T (3)] []32T [N3] dV
rB1 Hl,r, Czsra12,stT (o),tu 21,vut23,vwT (/),wx N3,yz dV +
TlCB2,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 +
TH2CI 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
Jv
rH1CB2,vC21,vw x2 21 12 21 r CB2,z dV
12 L21 12:21
rH1CB2,sa21,st a,ut 2U a,uv ;'21,wva23,w8 X (),xy (' TH2CB3,z dV +
Jv
12 21 12 21 r CB3 dV +
THlCB2,s21,st 2 autl 221 1 ., .23,wxT ( y 2 f dop,z dV +
JW
12 21 12 21 3 dV
rHlCB2,rt21,rs ats atuO21,vua23,vwT (/) ,w (' qf, N3 dV +
12 21 1221 f
JV
12:21 12:21 f
rB2H2,va21,vw axw a; 1, .,_ (YsrM _H2,z dV +
rB2H2,a21,st 21 12 21 T f, dV +
23v
B2H2,v 21,w 2 1 12 21 yB2H2,z
FB2H2,sCt21,st l2 ut 1a i lur 21wQS,wvT2 (33) 0' IrH2CB3,z dV +
rB2H2,sa21,st2 a,ut ( 21 .i ., 23,w (),xy op,z dV +
JV
12:21 WH12:2,v dV3
rB2H2,ra21,rs al2 ,tsl2 atuaO21,vua23,vw (/),wx ,z N3,y dV +
JV
1221 12:21
rH2CB3,sa32,stT ((j))ut a 21,vw 1, 1 :yrHlCB2,j z \ dV +
Jv
H2CB3,sT32,st ( 21,w l2 21 2 21, :yB2H2,z dV +
SrT2C ,, ,T(23 ,r)g 1s23,srQ 21,st2 nut uau~n (21,wt23,wa ,y H2CB3,z dV
TH2C ) a 21,st lt12 1 12 21 T. .. 23,ws ( op dV +
J7 V
rH2C, 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 :yTHCB2,z uop,s dV +
v/
a32,stT (3)t ,. '21i,vw2 2 1~~ l iyrB2H2,z op,s dV +
Jv
1221 12 21 T io
' ,rg 23,sr21,stl a,utl La,uv( 21,wvt ;23,wx y itH2CB3,z Hopp dV +
S, (),r a23,sra21,stl aut auva21,wva23,w X (3),xy Uop,pUop,z dV +
Jv
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, BH2,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
JJV
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 +
V+
V+
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 finaltime problems, fixed finalstate problems and
minimumtime 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.
CHAPTER 1
INTRODUCTION
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,
inthefield 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
where
/fvHTB,dV f T B BdV
[M] d p [B] [B dV p v 2 (7112)
jv fv HTBdV fvB T 2dV
Looking at a representative term of Equation 7112 and using the indicial notation,
V B BI dV
0 1 01 f
TBIH1,x 1,yx 1l,yzrBlHl,z dV +
rB1H1,w 0 ,xw 1 T (2) ,yz rH1CB2,z dV +
TB1Hw ol ,Xw
rB1 HI,w (dl, W I U (
B1Hl,w ,XW 1,1y ( 0,y, Uz dV +
rH1CB2,wT (2),w 01, 0 1,yz B1H,z dV +
JV
rH1CB2,vT (2), 0W w (1 2, 0) rHiCB2, dV +
JV
rHICB2,vT (2), 0W w 1,yT (2, yz) Uop u dV +
2 W, O ,y 12 2 9 H1B1Hlz oopw dV +
T (0 2) 0;J 1 T (02, yz) rHICB2,z j u^ dV +
T ( 02),W U 1 0, ,
Jv
0T 1 1 yT (12, yZ) U opUoz dV (7113)
Jv
It is important to note that the integrals left in Equation 7113 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 + B1H1 + HICB2 + B2H2 + H2CB3 + op
(7114)
where rB2H2 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 bodyfixed coordinate system
associated with body 2, rH2CB3 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 bodyfixed 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 743 produces:
[0] q + [t]t 0 (744)
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 744
reduces to:
[ q],q 0 (745)
Taking the second total derivative of Equation 743 or the first total derivative of
Equation 744 produces,
[1q = [] 2 [qt ([q q (746)
which for autonomous systems reduces to
[q _ ([:I^_q ") ]q (7 47)
Any system that simultaneously satisfies Equation 744 and Equation 746 also
satisfies Equation 743. Thus, using the accelerations of the generalized coordinates,
Equation 746 can be augmented to the system equations of motion to form the matrix
equation
[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 righthand side of Equation 746.
The constraint matrix for the threebar 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
lowpressure 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
CHAPTER 7
NUMERICAL EXAMPLETHREE BAR MANIPULATOR
The equations of motion that were generated in C'! ipter 3 are also applicable to other
multibody systems. One of the most common examples of a multibody 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 73, 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 73. At the
initial time all of the bodyfixed 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
system.
Figure 11: 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 (312)
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+ (313)
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. 313 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 (315)
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)
Thus,
S B RCB_ + [ t] o+ [0w] + 0[w ] (317)
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:
4
min J J E fdt
i=1
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 514and Equation 515.
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
LegendreGauss 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 11. As opposed to more traditional aircraft, the MAV presented in
figure 11 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[25], 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 =
rH1CB2,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)
Finally,
S[N2 [N2] dV N2,nmN2p dV (A14)
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 bodyfixed 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 (793)
where RCBI is written in terms of the inertial coordinate system and u is written
in terms of the bodyfixed coordinate frame associated with body 1. Taking the first
derivative of Equation 793, with respect to the inertial coordinate frame, yields:
Ed E (LRC)) +B E B' x (794)
but the revolute joint between the ground and body 1 prevents the rigidbody translation
of body 1, thus RCBI = 0. Therefore,
v [o01] tu (795)
where [0'1] is the skew symmetric matrix formed from EUB1. However,
[OL1] 0, [ (796)
Thus,
v = 01 [0)] Ut (797)
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 Gullwing MAV illustrated in
Fig. 11. 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 (7131)
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) (7132)
where
X = qs (7133)
F (X, u) }I (7134)
[ Tv + .U J
CHAPTER 2
LITERATURE REVIEW
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 lifttodrag 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 1ir 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 reattachment 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,000100,000.[12'13] At Reynolds'
numbers less than 200,000 the lifttodrag 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 7132, 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
solution.
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 + (712)
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 712 with respect to the inertial frame
produces:
Ed Ed E dU E Bd
t (_p) Ed (RCB1) l (OP) + x op
= RCBI + [ ] IP (713)
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)
Thus,
p RCBi + 1i [o01] rp + 01 [o0L)w] ]U, + i [0 ] (715)
Rewriting Equation 715 so all components are written in terms of the inertial frame
produces:
p, = RCB + 1 [Ci] [,1 u] P 01 [CI] [00 )L] U + [CI] [0
where b2 is a unit basis vector of the bodyfixed coordinate frame of body 2 that is
perpendicular to the bodyfixed 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 764 in terms of the inertial coordinate system produces:
6 b IT[Cl]T[2]b b 0 (765)
Taking the partial derivative of Equation 765 with respect to the generalized coordinates
yields:
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
(766)
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 (767)
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 bodyfixed 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 bodyfixed
coordinate system of body 3, and uH2 is the vector from the origin of the bodyfixed
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 bodyfixed coordinate system of body 2.
Writing all components of Equation 767 in terms of the inertial coordinate system yields:
7 = RCB3 + [C3] 32 CB2 [C2 22 = 0 (768)
Rewrite Eq. 317 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 ]
Then
VP= [Ci]T [B] R a
where
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]
2
(325)
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
(318)
(319)
(320)
(321)
T P{CB1
f p [C,] [B] dV
f p [B]T [B] dV
(322)
(323)
(3 24)
RCB1
o_
q 
APPENDIX A
MASS MATRIX FOR INBOARD WING SEGMENT
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
(A1)
Examining the integrals term by term and using the indicial notation:
S/lp[I]dV= 2pI,mdV (A2)
[CI] [B] [Ci] B2 B3 (A3)
Thus examining a representative term,
j [CI] B dV
C1,mn ,nprB1H1,p / dV +
Jv
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 (A5)
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)
SBBIdV JBB,2dV f BBdV _
BIOGRAPHICAL SKETCH
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. 5559 is required. The difficulty in
solving those equations lies in getting a closedform 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 infinitedimensional optimal control problem into a
finitedimensional 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 finitedimensional 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'70731 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 nonlinear programming problem, direct methods for the
solution of the optimal control problem have wellknown 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.[7678]
Taking the partial derivative of Equation 753 with respect to the generalized coordinates
produces:
0q= O eT [Ci] [0] b 4 [CI] ] [ ]bl 4 i, o [ ,] bl O 0 0 0 O 0 0 0
q (7 54)
(754)
Similarly,
S3 = e b' = 0
(755)
Or presenting all quantities in terms of the inertial frame:
3 = [CI] b'= 0
(756)
Taking the partial derivative of Equation 756 with respect to the generalized coordinates
produces:
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
(757)
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
(758)
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 bodyfixed 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 bodyfixed
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
(759)
For the remaining terms replace o0D, with O^1, and 0:D, where appropriate.
S[BB] DDidV j BB DD dV
< V v 2DD
BBTDD1
Examining a representative term:
SBBTDD, dV
rB1H1,ul ,vual2,vT (a),w 12;?1 a21,zyrH1CB2, V
rB1 H1,uOcv1,vua12,vt (a),., 12 21
BlHIu0 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+
rB1H1,r Wo,srZa12,stT (,12a) tu WLu 21,uwv23,tuT (,3)xy a32,zy uop,z dV +
rBlHl,q rqa12,rT (a),st 121tua21 ,vua23,wT (3) ,a z f N3,yz dV +
rH1CB2,ra21,rsT (a),ts l2,uto,0w ,o l2,uT (a),wx 12 ?21 ,
0 a21,zyrH1CB2,z dV d
rH1CB2,ra21,rsT (a),t ail2,ut O l,vual2,vTuT (a),wu 1 ,21,zyB2H2,z
0 a21,zyrB23H2,z dV af
rH1CB2,na21,npT (a)qp al2,rq0 ,sr a12,stT (a) ,tu 121,az21,v23,tu T (/3) ,y a32,zyrH2C3,z dV +
rH1CB2,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
B2H2,rO21,rT (a) ,ts al2,ut o ,vua12,vT (a),wxz 1,zyB22,z d
0 a 21,zyrB2H2,z dV d
rB2H2,na21,npT (a),qp a 2,rqO ,,s 12,stT (),tu 12w1 21,u23,xT (3), 32,zH2CB3,z dV +
rB2H2,n21,npT (a),qp a2,rq w,sral2,stT(a),tu 122,va21,wa23,wxT (3),xy 32,zy op,z dV +
rB2H2,ma21,mnT (a)pn al2,qp ,rqal2,rsT (a),st 1221 ta21,vua23,vwT (/), N3, d
rH2CB3,na32,npT (3),qp a23,rqa21,rsT (a),t, a12,ut ,vua12,vwT (a),wx 12 ,21, yfH1CB2,z d +
rH2CB3,na32,npT (),qp a23,rqa21,rsT (a),t al2,ut0 o,val2,wT (a),wx 12 21,zyB2H2,z V +
rH2CB3,ja32,jkT (),mk 23,nm 21,npT (a),qp 2,r ,sl2,stT (a),t 12w ,u21,w 23,aT (3)y 32,z
rH2CB3,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
rH2CB3,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,zyrH1CB2,z op,n +
a32,npT (3)qp a23,rqa21,rsT (a),ts l12,ut 0 cvu12,uT (a) ,x 12 ?1 ,t21,zyrB2H2,z op,n dV
v
(B9)
yrH2CB3,z ; dV
y op,z dV +
,,zfN
APPENDIX B
MASS MATRIX FOR OUTBOARD WING SEGMENT
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 (B1)
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 (B2)
where [I] is the identity matrix.
/ [Ci] [BB] dV = [Ci] BB, BB2 B 3 dV (B3)
Examining a representative term,
J [C11] BB dV =
Jv
C1,xy ryzBlH,z dV +
Cl,uv Lo l ... .T(a),_, ,irHICB2,z dV +
JV
C1,uv1aO ,vwO2,wxT (C),xy (,_ ,rB2H2,z dV +
1,rs o,sta12,tuT (a), 21,wva23,wxT (3)y ,rH2CB3,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 (B4)
Jv
For the remaining terms replace 011 with 0o1i and o1;, where appropriate.
J [C'] DD, dV
Cl,uval2,vwT ({a),w 12Z1, 2 '., :yrHICB2,z dV +
Jv
40
20
time
Figure 78: Time history of the angular rates as determined by integration.
20
40
0 0.2 0.4 0.6 0.8 1
time
Figure 78: 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
XCBI
YCBI
ZCBI
0
q
acl
01
a1
a2
02
For the purposes of this numerical example, the invariant integrals were carried out using
fourpoint quadrature along with an isoparametric element used for each body. Thus,
4
/ pdV = 1 det J p t (61)
i=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 (62)
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 + (1i) 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
lifttodrag 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
overarcing 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 Fluidstructure interaction is not only limited to aerospace interests. In
work by Kral and Kreuzer, the modeling of fluidstructure interaction was addressed as it
impacts offshore platforms.[35]
In order to capture this interaction, the choice of modeling becomes important. There
are several vi to describe any given system. Standard rigidbody 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 (377)
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 (378)
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 + (379)
but
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 (381)
Rewriting Eq. 381 in terms of the B1 frame yields:
6r [C[i] [B] { e C (382)
S[J 60
Cl,uvla2,vwT (c),wx 1 21, .1 :yrB2H2,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
Jv
BB BB BBT Bi
IBB BB, BB Bi
BBTBB2 BB Bi
,rH2CB3,z dV +
J Up,z dV +
Jv
q j N3,yz dV
Jv
TrH2CB3,z dV +
J V
op,z dV
fz N3,yz dV
BB BB3
1 3
BB BB3 dV
BBTBB3
Examining a representative term:
J BTBB, dV
JV1 1
rBIH1,x021 01 BlHl,z dV
rBlH1,uO ,vu0l ',val2,wxT(a),y a21,zyrH1CB2,z dV +
01 01 dV
rBlH 1,u0 ,^vu0 v"vwal2,wxa T(),xy a21,zyrB2H2,z dV +
rBlHl,r ,sro W1,stal2,uT (a),uv 21,wv23,wxT ( y3),xy 32,zy Uop,z dV +
0 1 01 dV
rBlHl,qOW,rqo wC,rsal2,stT(a) ,tu21,vua23,vwT (3),a N3,yzdV+
rH1CB2,ua21,uvT (a),w 2,xw 01 01 rBlHl,z dV+
rH1CB2,2,r2,rsT (a),ts al2,utO0 w,vuwOvwal2,w (a),xy a21,zyrH1CB2,z dV
rHlCB2,ra21,rsT (a),ts al2,ut ,uOvwal2,wx (a),y 21,zyrB2H2,z f dV +
rH1CB2,n~21,npT (),qp al2,rq0 s ,stO 2,tuT (a),uv 2,w v23,wx )T )xy a32,zyrH2CB3,z dV
rH1CB2,2,na2npT(a)qp al2,rq o ,sr0o ,st2,tuT (a),v 21,wva23,wxT (3),xy 32,zy op,z dV +
(B5)
(B6)
(B7)
B2
B2
V2
CHAPTER 4
DETERMINATION OF AEROELASTIC FORCES
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 3D
Navierstokes 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 fluidstructure 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 opensource
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 singlei 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 illbehaved.
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
equations.
19.7 19.9
Figure 61: Geometric model of gullwing 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 multibody formulation of the equations
of motion. The multibody 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 multibody
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 timedependent. 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
maneuver.
indicial notation convention is used. Thus,
f p [I] dV plm dV (326)
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
(327)
Looking at the first term of Eq. 327
Sp [CI] [0)] uOdV C Ci,imo, f dV (328)
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 (329)
_, 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 (330)
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 + LB1HI+ 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 (3t11)
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
SDDDD2 dV =
rHlCB2,sa21,tl2 ,ut W '1 a23,vwT ( ,w 2332 I rH2CB3, dV +
Hl12CB2,21 21 t T tvu3,vw 23 32 / opz dV +
B2CB2,sa2l,st ,l2 ut21 ,v ,vw ,w2 2 ,rI H CB3, dV +
Jv
TB2H2,sa21,st l2 ,uta21,v 23,v (/3),w 23 32 dV +
Jv
,,2: 3 j N3,, dV
TB2H2,r'21,rs1 221 23,u ( 2332 vT ,(,)z N dV+
Jv
rH2CL ,rq 3) a23,sr 21,stl221 ,r ) 23 32, ,H2CB3,. dV +
12:21 T23
rH2CL ,q a23,sr21,st5 ut221 u23,vwT ( ,w 23 32 UopI, dV +
J V
rH2CL 1 T(3) 232, '22, 1, 0 .,' 23 132 , 3 N3,, dV +
1221 T(0w 23 32 /pd
J V
,T q () 23,sra 21,st2 a ,ut a21,vu T23,vw ( 2332 I rH2CB3,z op dV
Jv
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
rH1CB2,sa21,t ...T (/3),w 3,, dV
JV
Jv
rB2H2,C 21,t w L,ut a2,1,23,vw (21 ),wx tl2 ( N3,v 3, dV +
H2CT (q
1221 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 (B17)
Jv
/ DD DD2 dV
Jv
23r32 23 32 B, d
rH2CB3,v<132,vw W3 w rH2CB3,z dV +
23L32 23 32 3fN, dV
TH2CB3,v1'32,vw 3, w L .q pz +
2332 2332 /i
TH2CB3,uQa32,uv23 3223 32 N3,uy dV +
'32,vw 3 3,xw r IH2CB3,z Uop,v dV +
Jv
3 23 32 23 32 dV +
a32,vw W3 w J
Q32,uv23,wv 2 32, , ,fz Uop,uN3,yz dV +
Jv
f ..... /32, 3 TH2CB3,z 3,vu dV +
3 2332 2332 [
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
2332 f
Jv
0H2CB3,vw 32,vw op 3yz dV +
JV
qf w23u32 2,vw L,, N3,vuN3,yz dV (B19)
Jv
Finally,
S[N3T [N3] dV N3,yN3,yz dV (B20)
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 KarushKuhnTucker 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 62 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 63 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 64 and Figure 65. 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 nonzero 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:
CB1
vPp [C]T [B] D, [a12] [T (a)] [221]T [N2] (3 48)
f2
Now the kinetic energy of body 2 can be written.
T j= tp(v v) dV = pvvPdV (349)
Substituting Eq. 348 into Eq. 349 yields
1
T = 4 [M]q (350)
2
where
RCB1
S< > (351)
.2
[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(352)
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 bodyfixed coordinate frame associated
with body 3. Taking the first time derivative of Equation 7115 with respect to the
inertial coordinate system yields:
Sd B + d E X BH HCB2 E B HCB2
dt dt C d\B1t (rBlHI1 H1+ dt CB2I + Wl
2 d B2 B3 d E EB3
S (rB2H2) X LB2H2 + (rH2CB3) X H2CB3
d dtC
B3 d (3 x (7115)
However,
E B2 E B1 +B B2 (7 116)
EB2 E + B1 B2 + B2B3 (7117)
Inserting Equation 7116 and Equation 7117 into Equation 7115 and factoring the skew
symmetric matrices produces:
vp 1 ([ 1 rBH, + [01 ] [T (2)] rH1CB2 + [0)j1] [T (02)] rB2H2
[o0 ] [T (02) [T (3)] rH2CB3 + [O:1] [T (02)] [T (03)] p) +
2 ([T (2)] [12] H1rCB2 + [T (02)1 [1] rB2H2+
[T (02)] [1W2] [T (03)] rH2CB3 + [T (02)] [1)2 [T (03)] p)
<3 ([T (02)] [T (03)] [11 TrH2CB3 [T (2)] [T (03)] [1] o) (7 118)
Define
BB1 = [ ro H)i +] BH1 [T (2) H1CB2 + [0)Wj1 [T (2)] r2H2 +
[O] [T (02)] [T (03)] rH2CB3 + [0J w [T (02)] [T (03)] p (7119)
BR2 [T (02)] [1] rHICB2 + [T (2)] [1] W r B2H2 +
[T (2)1 [1w] [T (3)] rH2CB3 + [T (2)1 [1W] [T (03)] (7 120)
[T(02) 621 21 _O
stated as:
"B2 [T()r3 (88)
where
1 0 0
[T(3)] = 0 cos 3 sin 3 (789)
0 sin 3 cos 03
and rB3 is a vector written in terms of the bodyfixed coordinate system associated with
body 3 and rB2 is the same vector written in terms of the bodyfixed 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 bodyfixed coordinate frame associated
with the first body and the inertial frame is written as:
E 1 13 (790)
where b is the unit vector of the bodyfixed 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 bodyfixed coordinate frame associated with the second body and the
bodyfixed coordinate frame associated with the first body is written as:
B1 B2 2b (791)
where b is the unit vector of the bodyfixed 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 bodyfixed coordinate frame associated with the third body and
the bodyfixed coordinate frame associated with the second body is written as:
B2 B3 3_ (792)
w1
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0 0.2 0.4 0.6 0.8
time (sec)
Figure 71: Norm of the constraint vector
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0.4 0.6
time (sec)
Figure 72: Norm of the constraint vector
Taking the partial derivative of Equation 769 with respect to the generalized coordinates
produces:
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
(769)
Looking at the constraints on the relative rotational motion between body 2 and body 3
imposed by the revolute joint:
s = b b 0 (770)
where b is the unit basis vector of the bodyfixed 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 bodyfixed coordinate frame in body 3 that is perpendicular to
the unit basis vector of the bodyfixed 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 (771)
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 772)
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 (773)
where b3 is a unit basis vector of the bodyfixed coordinate frame in body 3 that is
perpendicular to the unit basis vector of the bodyfixed 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 773 in terms of the inertial system yields:
9 = b2T [2 [C3] b 1 0 (774)
Taking the partial derivative of Equation 774 with respect to the generalized coordinates
produces:
6i = 0 00
0 0U i [2[ ( 2\:( 3 3
o or ]bi L1]T [C2] T [C3]_b7
(775)
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 748 are ready for integration. There are
several methods to integrate these equations. The simplest involves direct integration.
Define
[M] [D]T
[MA] (776)
[][ q [0]
q x (777)
A
F } (7 78)
QC
Thus, Equation 748 is written as:
[Mx] qA F (779)
4 [Ci] [Oi[0]
C0 [i [0 (386)
TC1], [ [l [0:]
[T / =[T (Q)] [i2w:]
Thus
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 ( HCB2+
[Ci] [OLO] [ [a [a2i]T U2 + [Ci] [01;] [ai2 [T (a)] [0a2lT [N2] q) 2 +
OcI] [0L ]]Tb [CIBH1 1 [a121 [T (a)] [a2 H1CB2
[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 (387)
Rewriting Eq. 387 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:
6RCB1
6a
6q 2
ff
Performing the dot product of the forces acting on body 2 with the corresponding virtual
displacement vector of the form of Eq. 388 and comparing the result with Eq. 377
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
TOPICS IN THE STRUCTURAL DYNAMICS AND CONTROL OF AN
ARTICULATED MICRO AIR VEHICLE
By
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.1745.
[92] Benson, D., "A Gauss Pseudospectral Transcription for Optimal Control", Ph.D.
Thesis, Massachustts Institute of Technology 2004.
REFERENCES
[1] Shabana, A.A., D;in,.i of Multibody S';,. n New York NY: John Wiley & Sons,
1989.
[2] Meirovitch, L. and Stemple, T., "Hybrid State Equations of Motion for Flexible
Bodies in Terms of QuasiCoordi ii. ,Journal of Guidance, Control, and D;.i,,.i
ics, Vol. 18, NO. 4, 1995, pp.678688.
[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/CR2003211748, 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.714727.
[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,
pp.814828.
[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.13421353.
[7] GadelHak, M., i\! roAirVehicles: Can They be Controlled Better?", Journal of
Aircraft, Vol.38, No.3, 2001, pp.419429.
[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.425465.
[9] Mueller, T.J. and DeLaurier, J.D., "Aerodynamics of Small \. !hi !. Annual
Review of Fluid Mechanics, Vol.35, 2003, pp.89111.
[10] Muti Lin, J.C. and Pauley, L.L., "LowReynoldsNumber Separation on an
Airfoil",AIAA Journal, Vol.34, No.8, 1996, pp.15701577.
[11] Lutz, T., Wurz, W., and Wagner, S., "Numerical Optimization and WindTunnel
Testing of Low ReynoldsNumber 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.,
"FlexibleWingBased Micro Air \. !i 1!, American Institute of Aeronautics and
Astronautics Paper No. AIAA20020705, 2002.
[ ] [12] [T (a)] [a ] rB2H2 +
[o0w] [012] [T (a)] [a2] [a23] [T (03)] [T H2CB3 +
[Owl] [ 12] [T (a)] [a2l]T [a23] [T (03)] [a32]T +
[O] la2 [T [(a)] [2T rB2H2 +
[0O] [012] [T (a)] [a21l] K ] [T (8)] r rH2CB3 +
[ 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 rB2H2+
[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) (363)
Define
BB = [o] rBlHi + [w] [a12] [T (a)] [21] TarHCB2 +
[o1] [al2] [T (a)] [a21]T rB2H2 +
[0 [a 2] [T] (a)] [a21]T [a23] [T (3)] [a3]T 2CB3 +
[oL] [2 [T0 (a)] [a21T [a23] [T (/3)] [a3a]
[O [C 12] [T2] (a)] [a21]T [23] [T (3)] [a32]T [N3] q 3 (364)
BB2 [0 B1] _rB1Hi + [01w] [a2] [T (a)] [a21]T H1CB2 +
Since, [4], has full row rank, [MA] is invertible and thus
q = [A]1F (780)
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 ( (781)
Thus,
X = (782)
S[]1 Q T ( U (7 82)
where u represents a control input to the system. Figure 71 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 779. 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. 83851 This objective is reached
by making the constraint equations analogous to the standard secondorder ordinary
differential equation of control theory. This analog is represented by:
K + 2aK + 32 0 (783)
and 0 D, where appropriate.
Sv BDidV
S[B]T DIdV = f B DidV (A8)
fS B DidV
Examining a representative term,
SBRTD dV =
v
01 12(21
rBiH1,i L ,o, T ,km W ,, TH1CB2,p dV +
rB1Hl,i 3, ji(' m a12 mn, (I / opp dV +
JV
01 1221 2dV
TBlHs L jT (a), a,krn21,nmrfp N2 dV +
TH1CB2,wa21,wrT (a,, a12,is T (a) ,km 12nn( rHlCB2,p jdV +
JV
TH1CB2,wa21,wrT (a),sr 12,iso i T (a) ,k n 12 n I dV +
SJvann
Jv
0 12 1 u dV +2
T21,wrT (a),,, a12,is () ,km mn(' rH1CB2,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 +
Jv
0( 1 )121 /
q2 .1 .,rT (a),sr o 12,is 0 ji T (a) ,kn n ( I f,nrH1CB2 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 (A9)
JV
CHAPTER 6
NUMERICAL EXAMPLEARTICULATED MICRO AIR VEHICLE
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. 61 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 (799)
2
where q = i and using the indicial notation,
M 0oijy 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 + LB1HI H1CB2 + p (7101)
,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 bodyfixed coordinate frame
associated with body 1, rHICB2 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 bodyfixed 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 bodyfixed coordinate frame
associated with body 2. Taking the first time derivative of Equation 7101 with respect to
the inertial coordinate frame yields:
E d d Bid + EwLB1 Hir
dt (E ) E (CB1) dt B1H) + BIH
B2d (rH12) + E xB2 + B2 d (B) + Ep7102)
but
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, rB1H1 is in terms of the B1 basis, rHICB2 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 BH1 + 2 H1CB) +
(L X _) H +B2 + x u^ + B d[ ([N,] q +
EwB2 H _ B2 d ( B+ B + 22)
E B2 x [N2 (332)
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 B1H1 E X H1CB2 d \ J /
(333)
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 rewriting
all cross products in terms of the same basis yields:
EB1
VP RCBI + V1 x B1 Hi + EB x [a121 [T (a)1 [ITi LHICB2 +
2 >21x + EHi CB +EWBI x [a12] [T (a)I [acl+ u +
J1221 [a21iT [N2] (335)
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] BIHi + [LO1] [Cai2] [T (a)] [ailT HICB2
rH1CB2,uT (02),vu 0w 2 ,y T y H2CB3,z dV +
rHCB2, 2) 0w1 0 1 2 w dV +
JV
rB2H2,uT 2 ( w 01,y0 2,yzHjB1Hl,z dV+
Jv
rB2H2,vT (02), w, 0 (/2),yz (tjH1CB2,z dV +
6 6 1,yzrB H
rB2H2,vT ( T2),W O1,W0 w T (yT(02),yr rHCCB2,z dV +
,x r HfJV
rH2H2,T 02),Wv T ( L 0 1 0w1 T +2V
S 01 0 dV +
6JJV
rB2H2,uT (02),v 01, 0 1,wT 0y T (3 ),y rH2CB3,z dV dV +
6 6 1 ,JJ v
TH2B3), 3~ ,u 30w 1 0 1( T (2),yz THCB32,z dV +
TH22CB3,T (3),t (3) 01 0 1 T3,yz H2,CB3,z d dV
Jv
rH2CB3,vT (33) ,w T (23), o X ,yz BlHlz op,vz dV +
JVv
rH2CB3,uT (3),v T 101 T ( 2 ,yz H CB2,z rHop, dV +
T T1Jv
rH2CB3,T ( vvu T v O 1 2 T2,yz B2H2,z op,u dV +
rH2 3,T (3 3),,.t 2 ( v), O)1,W 2) T (y3),yz TH2C3,z op,t dV +
; x 01 01 T U V+
T ( 3),t T 2) 0 1,w I, ,y T ( ,yz op,topz dV (7126)
Jv
Here again, the remaining integrals in Equation 7126 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 73: Three bar manipulator.
[14] Hall, J., Mohseni, K., Lawrence, D., and Geuzaine, P., l is i:, i of Variable
WingSweep for Applications in Micro Air V\ !i 1!, American Institute of
Aeronautics and Astronautics Paper No.AIAA20057171, 2005.
[15] Ramamurti, R. and Sandberg, W., "Simulation of the Dynamics of MicroAir
V\ !i, !I 38th Aerospace Sciences Meeting 6 Exhibit Reno, Nevada, Paper
No.AIAA20000896, 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,
pp.12571265.
[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.AIAA20041674,
2004.
[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. AIAA20024875, 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.AIAA20014005, 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.424426.
[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.24922494.
[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,
pp.865873.
[24] R in, :lr'vani, D.S., Lesieutre, G.A., Frecker, M., and Bharti, S., "Aircraft Structural
Morphing Using TendonActuated Compliant Cellular Ti, . Journal of Aircraft,
Vol.42, No.6, 2005, pp.16151621.
[25] Bae, J., Seigler, T.M., and Inman, D.J., "Aerodynamic and Static Aeroelastic
('C!i lteristics of VariableSpan Morphini Journal of Aircraft, Vol. 42, No.2, 2005,
pp.528534.
cb2CO2 sQ2C02 + c2AsO02s 2 sQ2Qs2 + cP2s020c 2
[C21 _C2 _( C_ + ',_ ( ', + .', (75)
s02 c02Se 2 c02CQ2
( 3 3 1 (' '103'3 ,' 3S 3 3 + (' '03C 3
[C3] 3 ( ,' .)3 + ,, + ,' )3 (76)
s03 c03s 3 c03sc3
where cQi is cos Qi and sQ is sin Qi. As all three bodies have similar relations between
their respective bodyfixed frame and the inertial frame, they have similar relations for
the angular velocity of those bodyfixed frames with respect to the inertial frame. These
relations are:
E LB (4j +1O1) bt + (Oico + iCOIS1) bl + (Oisoi + QicHicOi) b (77)
E1B2 (= 2 .) + ( i022 + 't0'C0252) b2 + (02SQ2 + '_C 02C62) b6 (78)
E B3 (43 03) b3 + (33CQ3 + 0* '3 b) 3 + (O3S3 + 0c) 3 (79)
where the angular velocities are written in terms of the basis vectors of the bodyfixed
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 (711)
Jv 2
of the force application on body 3 is written as:
S CBI + BHI + rH1CB2 + B2H2 + rH2CB3 + op + (389)
Rewriting Eq. 389 so that all components are in terms of the inertial frame yields
rp = RCB1 + [Ci] rB1H1 + [Ci] [a2] [T (a)] [i] rHICB2 +
[Oi] [l12] [T (0)] [a21]T B2H2 +
[C] [la121 [T (a)] f ]T [, ] [T (/3) [2 rH2CB3 +
[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]lH1)
I6 + +
a ([C1] [a12] [T (a)] [a2]T H1 B2 ) a ([Cl] [1121 [T (a)] [a21]T H10B2)
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+ OS +
a ([c1] [a12] [T (a)] [021]T [a23 Er (8)] [a32] TH2CB3)
a ([C lj [a12] [T (a)] [a2l]T [a231 [T (/3)1 [0 32]TTr /2C 3)
a0
9 ([C1] [012] [T (a)] [a21]T [023] [T (3)] [a32 ]T H2CB3)
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 +
a8
a ([C1] [a121 [T (a)] [2j]T [a231 T ()] [a32]T_]a2 0S)
ap
a ((c)] ]a12] 1] (a)] [a2ljT []a13 ]T (8)] []32]T_ !u )
a60 +
o~ ~ ~ ~ ~ ~~~~~~~S ([j[11I 1[lT[1I /) T. 6
de+
generalized coordinate. The preceding information is sufficient to define all of the
quantities on the left hand side of the Equation 510. Next the appropriate quantities
on the right hand side of the equations of motion must be defined. The gravitational
forces are applied to each multibody 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
S
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,,
as:
u Q] T Q = U (7127)
where
Ci
U = C2 (7128)
d (T[[
d ([ 1)D + t a (T (7 129)
 dt 2 O"qS sys 8
ss 2 (7130)
03
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
CHAPTER 5
OPTIMAL CONTROL STRATEGY
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.[5867,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 (51)
to
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 Ndimensional vector of states and u (t) is an Mdimensional vector
of control inputs. Equation 52 shows that the states x is a function of the control u.
Similarly,
E B3 E _B1 J12 21 J23 32 (358)
Therefore,
Vp C BI + EB X rB1H1 + (EI + J2 1) X rHCB2 +
E B J2 )2 X rB2 H2 E J21 J 23 + 32) x H2CB3 +
E B J2 J 2 J23 j32 3
(EBL1 + L V3+ )x P + [N3] _j +
(E B + ) J 21 + J23a2) x [N3] (359)
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 H1CB2 +
J121 X [a2 21] rH1CB2 + E 1 x [O12] [T (a)] [2i]T B2H2
P2 21 x [a221]T B2H2 + EB X [Ca12 [ (C [1T [o [T (/3) [] rH2CB3a
2 2 x [a21]T [(o ]i [T (/3)] [] H2CB3 + 33 X [( T rH2CB3 +
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 (360)
Introducing skew symmetric matrices to simplify the cross product expressions and
rewriting 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 rHCB2 [0: 1] [a12] [T (a)] [a2TrB2H2 +
[a 12 [T ( 12)21] [21 ]TB2H2 +
[Owj1] [I12 [T (I)] [a21]T [T (/3)] IrH2CB3 +
[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. A96347120912, July 1996, pp.4714793.
[40] de Jalon, J.G., "TwentyFive Years of Natural Coordi ii., ,Multibody System
D i,.i, Vol. 18, 2007, pp.1533.
[41] Leger, M. and McPhee, J., "Selection of Modeling Coordinates for Forward Dynamic
Multibody Simulations", Multibody System Dn,,ni,. Vol. 18, 2007, pp.277297.
[42] Anderson, K.S., and Duan, S., "Highly Parallelizable Loworder Dynamics
Simulation Algorithm for MultiRigidBody Systems", Journal of Guidance, Control,
and D.,u,,. Vol.23, No.2, 2000, pp.355364.
[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.791822.
[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. AIAA984927, 1998.
[46] Schaefer, E., Loesch, M., Rulka, W., "RealTiem Simulation of Flexible Space
Manipulator Dyir 1i.1 American Institute of Aeronautics and Astronautics, Paper
No. AIAA9737062, 1997.
[47] Kielau, G. and Maiber, P., Nuinholonomic Multibody Dyir 1,i11 Multibody System
D i, ,. . Vol. 9, 2003, pp.213236.
[48] Aoustin, Y. and Formalsky, A., "On the Feedforward Torques and Reference
Trajectory for Flexible TwoLink Arm", Multibody System Di',,,ii. Vol. 3, 1999,
pp.241265.
[49] Junkins, J.L., "Adventures on the Interface of Dynamics and Control", AIAA, 35th
Aerospace Sciences Meeting 6 Exhibit, Reno, NV, Jan.69, 1997.
[50] Rives, J., Portigliotti, S., and Leveugle, T., "Atmospheric ReEntry Demonstrator
Descent and Recovery SubSystem Qualification Test Flight Dynamics and
Trajectory Modeling during D. .i. i American Institue of Aeronautics and As
tronautics, Paper No. AIAA971441, 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.13921397.
[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 rHCB2 + [C [[0 '][a [o] T (a)] [ fT rB2H2 +
[Ci] [o)] [012] [T (a)] [21T [ [T (3)] [(I 2CB3 +
[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 rB2H2+
[Ci] [)] [121 [T (a)] 21] [a21] ] [T ()],, ]T r 2CB3 +
[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; (393)
Using the vectors and matrices defined in the derivation of the contribution to the kinetic
energy from body 3, and rewriting the virtual displacement of a point on body 3 in terms
of the B1 basis yields:
6RCBl
60
,p [C 1] [BB]DDDD2[a2 [T (a)] [21]T [o ][ (3)] [32]T [N3] 6a
6q3
(f 9
(394)
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 threedimensional 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 bodyfixed coordinate system and these systems are all related to a globalinertial
coordinate system that is used to describe the motion of the system as a whole. These
coordinate systems are illustrated in 74. Mathematically, the bodyfixed 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' (71)
e2 [C2] b2 (72)
e3 [C3] b3 (73)
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 bodyfixed reference frame of bodies one, two, and three
respectively. Also, the transformation matrices are constructed using a 123 Euler angle.
Thus,
cylcOl syicl + ceplsOlspl sQlsol + cep1sO1cp1
[Ci] = syic1i c ico1 + ssO1isk c isoi + s isO1ici (74)
so1 cOisQi c1icQi
[52] Tadikonda, S.S.K., "Articulated, Flexible Multibody Dynamics Modeling:
Geostationary Operational Environmental Satellite8 Case Study", Journal of
Guidance, Control, and D;.,,.i,, .:. Vol.20, No.2, 1997, pp.276283.
[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.830837.
[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.665672.
[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.3744.
[56] Dias, J.M.P. and Pereira, M.S., "Sensitivity Analysis of RigidFlexible Multibody
Systems", Multibody System D m;,,i. Vol. 1, 1977, pp.303322.
[57] Ortiz, J.L. and Barhorst, A.A., "Modeling FluidStructure Interaction", Journal of
Guidance, Control, and D ;,,in. Vol.20, No.6, 1997, pp.12211228.
[58] Paiewonsky, B., "Optimal Control: A Review of Theory and Practice", AIAA
Journal, Vol.3, No.ll, 1965, pp.19852006.
[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.148156.
[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.428434.
[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,
pp.231239.
[62] Seywald, H. and Kumar, R.R., "A Method for Automatic Costate Calculation",
AIAA, Guidance, Navigation, and Control Conference, San Diego, CA, Paper No.
AIAA963699, 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.12521261.
[64] Venkataraman, P., "Dynamics and Optimal Control", AIAA Guidance, Navigation,
and Control Conference and Exhibit, Dever, CO, Paper No. AIAA20004560, 2000.
quadrature used here,
1 1
1 1
I = v< 3 > v=< > (63)
1 1
V3 V3
1 1
Similarly,
4
/ V = p opp det J t (64)
where op,i is the position of the quadrature point. Lastly,
4
/ /'I" Up,pdV p op, Up det J t (65)
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
100
40
20
20
20
40L I
0 0.2 0.4 0.6 0.8 1
time
Figure 74: Time history of the angles as determined by GPOCS.
150
100
50
50
0
0.2 0.4 0.6 0.8 1
time
Figure 75: Time history of the angular rates as determined by GPOCS.
a ([C1] [E12] [T (a)] [ 23] [] (3)] [( 328a]T ,)
ab
a ([C]] [a12] [T (a) 21]T [] [a ] (8)] [ ]T o') +
a ([c1] [al2] [T (a)] [al]T [023] [T (3)] [ ]T32] ) +
Sa1
a0
a ([c1] [a12] [T (a)] [a1]T [23] (3)] [ [3]T [N ] qf)
ao
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)
aa
a ([C1] [a12] [T (a)] [a2]T [oa23] [T (3)] [a32]T [N3] q3)
= [c ] [oo ]
a1
[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 (391)
f_
but
a(ciE [Cl] [0:l]
ac ]o [OCi [0:] (392)
a o(a [T (a)] [i2Bl3]
aor8)] [ [T\ () [23:32]
Thus,
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 rrH2CB3 +
[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 rH2CB3 +
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 (513)
Thus the equations of motion can be written as:
S=f (x, ) (514)
where
f (515)
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 (516)
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 516. 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 4spacecraft 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 secondorder equations of motion into a firstorder form that is
standard for optimal control problems. The secondorder equations of motion are written
as:
[M] Q +[K] q = Q +Q +Q +u (510)
T grav aero prop
Where
RCB1
al
02
q= < > (511)
02
/32
'f
q 4
[Q=[ }q+ ': (512)
.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 spaceplatformbased
mobile manipulator system.[541 Attia uses a multibody formulation to model a pendulum
with a rotating end.[551
The benefit of using a multibody dynamics formulation is that the individual parts
can be modeled in completely different viv. 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 flexiblerigid 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
ConfigurationsPast and Future", Journal of Aircraft, Vol. 40, No.6, 2003,
pp.10471065.
[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. AIAA981927, 1998.
[28] Reich, G.W., Raveh, D.E., and Zink, P.S., "Application of ActiveAeroelasticWing
Technology to a JoinedWing S. i .~ I i Journal of Aircraft, Vol 41, No.3, 2004,
pp.594602.
[29] Zheng, L. and Ramaprian, B.R., "A PiezoElectrically Actuated Wing for a Micro
Air \. !i. !, American Institute of Aeronautics and Astronautics Paper No.
AIAA20022302, 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.9499.
[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,
pp.131137.
[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.105121.
[35] Kral, R. and Kreuzer, E., \!liil lhody Systems and FluidStructure Interactions
with Application to Floating Structures", Multibody System D;.inii. Vol. 8, 1999,
pp.6583.
[36] Shabana, A.A., "Flexible Multibody Dynamics: Review of Past and Recent
D. !,pin, iii Multibody System D;.,iii. Vol. 1, 1997, pp.189222.
[37] Schiehlen, W., "Research Trends in Multibody System Dyri 1i1i1, Multibody System
D;.ir,. Vol. 18, 2007, pp.313.
[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.3558.
bodies that comprise the system. Thus the system generalized coordinates can be written
as:
RCB1
01
01
RCB2
q (738)
sys
02
62
RCB3
03
03
Thus the system kinetic energy takes on the form of
[M1]
1.T
Ty f :, _V] qsys (739)
Substituting Equation 739 into Equation 710 results in:
[M] q = : + Q (740)
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.
rH1CB2,ma21,mnT (a)p a12,qp0 crq0 W,rsa2,sT (a), 21,vuta23,vwT (3) 32,yqf, N3,yz dV +
rB2H2,va21,uvT(a) z a)12,xw I I rH dV +
rB2H22,ra rsT (a)t a12,utO 1 0 w 2, (, 21,zyHC2,z dV +
rB2H2,ra21,rsT (a) ,t al2,ut Dvu ,vwa2,wT (a)y 21,zyB2H2,z dV +
1 0 rf
rB2H2n2,n, npT (a),qp a12,rq0 esrO wstal2,tuT (a), a21,wva23,wxT (3), a32,zyrB2CB3, dV +
rB2H2,na21,npT(a) qp a12,rq W0 srO, wstal2,tuT (a), a21,wva23,wxT (3)xy a32,zy p,z dV +
rB2H2,ma21,mnT (a), pn a2,qp wrq O0w ,rsal2,stT (a) t a21,vvu23,vwT ( ) ,w a f N3,yz dV +
rH2CB3,ra32,rsT (3),ts 23,uta21,uvT (a),w 12,xw 0 rB1Bll HI dV +
rH2CB3,nOa32,npT (3),qp a23,rqa21,rsT (a),t a12,ut0o 1,vuO 'w 2,wxT (a),xy Z21,zyrHCB2,z dV +
rH2CB3,na32,np qp 3rqa21,rsT (a),ts a12,ut0, ,vwl2,wxT (a),y 1,zyrB2H2,z dV +
rH2CB3,ja32,jkT (),mk a23,nma21,npT (a)q 12,rq ,sr 0 W,st 12,tuT (a),u a21,wva23,wxT (3),xy Ca32,zyrH2C B3,z dV
rH2CB3,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 +
rH2CB3,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 01 rB Hl,z v op,r d +
a32,npT (),qp a23,rqa21,rsT (a),ts al2,utO,vuO1 ,,vwal2,wT (a),y 21,yrHCB2,z op,n dV +
a32,npT (3),qp 23,rqa2l,rsT (a),ts a12,ut vu0 vw 2 T (a)y a21,zyrB2H2,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,zyrB2CB3,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 rBHlHI,z N3,rq dV +
q9f,mO32,np (Q),qp a23,rqa21,rsT (a),ts al2,rqOw9vuo vwal12,wxT (a),y a21,zyrHlCB2, N3,m dV +
q3 002 N
,ma32,npT (),qp a23,rqa21,rsT (a)ts a2,rq1 01,vu vwl2,wxT(a),xy 21,zyrB2H2,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,zyrB2C'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)
CHAPTER 8
CONCLUSIONS
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 multibody 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 multibody equation of motion formulation on a
model that simulates the configuration changes possible on MAVs. The threebar 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 bodyfixed 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 (784)
where e, is a vector written in terms of the inertial coordinate frame and b{ is a vector
written in terms of the bodyfixed coordinate frame associated withe body 1. Here,
cos 1 sin Q1 0
[C] = sin 1 cos 1 0 (785)
0 0 1
In the alternative formulation illustrated in C'! plter 3 showed joint frames that were used
as the basis for relating the bodyfixed coordinate system of one body to the bodyfixed
coordinate system of an adjoining body. For this example these joint frames are not
necessary as the bodyfixed 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 (786)
where
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 bodyfixed coordinate system associated with
body 2 and rB1 is the same vector written in terms of the bodyfixed coordinate system
associated with body 1. Similarly, the relationship between body 2 and body 3 can be
LIST OF FIGURES
page
Figu
11
31
32
64 Position of MAV from integrated optimal control. ............... 72
65 Morphing angles from integrated optimal control. ............. .. .. 72
71 Norm of the constraint vector ............... ........ .. 103
72 Norm of the constraint vector ............... ........ .. 103
73 Three bar manipulator. ............... ........... 104
74 Time history of the angles as determined by GPOCS. . . 105
75 Time history of the angular rates as determined by GPOCS. . .... 105
76 Time history of the control torques. ............... ... .. 106
77 Time history of the angles as determined by integration. . . 106
78 Time history of the angular rates as determined by integration. ........ 107
re
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 gullwing 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 (741)
Ta[M]
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. 73; one between the ground and body 1 that allows
rotation about the bi axis of the bodyfixed coordinate system associated with body 1,
one between body 1 and body 2 that allows rotation about the b2 axis of the bodyfixed
coordinate system associated with body 2, and one between body 2 and body 3 that
allows rotation about the b3 axis of the bodyfixed coordinate system associated with body
3.
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 (743)
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 undercambered 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 : tan1 (U)
angle of sideslip : sin1 (v)
where
u XCBI
v = [C]T YCB1
w ZCBI
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 Zaxes 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 multibody 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 vv a multibody system can be analyzed. The two most
common techniques are a NewtonEuler formulation and a Lagrangian formation. Work by
Shabana, Schiehlen, and Yoo demonstrate some of the other popular choices for composing
the equations of motion.[3638] The NewtonEuler 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 multibody 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.[4246] 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, twolink 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 reentry 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
time
Figure 64: Position of MAV from integrated optimal control.
10
0
0.5 1 1.5 2 2.5
time
Figure 65: 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,zyrH2CB3,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,zyf1CB2,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,zyfH2CB3,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 (B10)
For the remaining terms, replace o01, with o0), and a1I where appropriate.
BB DD2
S[BBT DD2dV = BBTDD2 dV (B11)
BBT DD
Examining a representative term:
/BB DD, dV =
rBlHl,ro0's,ral2,stT (a),tua21,vuat23,vwT(3) ? a32,zyfH2CB3,z dV +
rBlHI,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+
rH1CB2,na21,npT(a)qp al2,rqO wsral2,stT(a)tu a21,vua23,vwT (3),wx 23 ,32,zyrH2CB3, dV
J V
rH1CB2,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 +
rB2H2,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
J 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 +
01 V
a32,jkT (P)k a23,nma21,npT (a) p a12,rq osral2,stT (a),t 21,vu~a23,vwT (3),ax 23 a32, yrH2CB3,z uop,j dV +
,q3 ,z H CB f oV d
Figure 31: Articulated MAV.
Figure 32: Coordinate frames of the fuselage, inboard wing segment, and outboard wing
segment.
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 MooneyRivlin
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 lowdrag and highlift 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 33: Joint coordinate frames between fuselage and inboard wing segment.
Figure 34: 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 shortcomings 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
reconfigurable MAV. Instead, more robust measures such as the use of a fully 3D
NavierStokes CFD software would be more appropriate. However, the use of such
a software on a reconfigurable MAV is a research topic in itself. Also, the current
stateoftheart 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
lifttodrag ratio of the proposed second generation vehicle for the University of Colorado's
MAV.[14]
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 bodyfixed coordinate frame associated with body 1:
[01] BH + [1 1] [T ( 2)] rHCB2 + [T (2)1 [2 rHCB2
+ [O' ] [T (42)] U + [T (02)] [12] u (7 104)
The skew symmetric matrices can be factored to yield:
v = di ([",] rB1H1 + [0 ] [T ()1 rHICB2 + [w1] [T (2)] Uop)
'2 ([T (2)] [1 2] rH1CB2 + [T (02)] [1 _2] P) (7 105)
Let
B L = '[0] rBiHi + [BH1 [T (0H2)1 rH1CB2 + [01,] [T (02)] Uo (7 106)
B2 [T (2)] [w122J rHICB2 + [T (_2)] [1W2] U (7 107)
[B] [ B B2 ] (7108)
Thus the velocity of a random point on body 2 is written as:
v, = [B] (7109)
2
Inserting Equation 7109 into the contribution to the system's kinetic energy from body 2
yields:
T = 1 '2 } 2Pp [B] [B] dV (7110)
Thus, the contribution to the kinetic energy from body 2 assumes the familiar quadratic
form:
T =4[M] (7111)
2
If the roots of Equation 783 lie in the lefthalf plane, the solution is guaranteed to
converge to zero. Thus, with the proper choice of a and 3 in Equation 783 the constraint
vector and the first and second time derivatives of the constraint vector can be made to
converge to zero.[8688] Standard application of the Baumgarte stabilization method
assumes 3 = a, with 15 < a < 50.[89] Figure 72 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 FiniteTime Feedback Controllers for
Nonlinear Systems with Terminal Constraints", Journal of Guidance, Control, and
D;i,.i, Vol.29, No.4, 2006, pp.921928.
[66] Betts, J.T., "Survey of Numerical Methods for Trajectory Optimization", Journal of
Guidance, Control, and D.,i,.. Vol.21, No.2, 1998, pp.193207.
[67] Hull, D.G., "Conversion of Optimal Control Problems into Parameter Optimization
PiN1.!. in Journal of Guidance, Control, and D.in,,. Vol.20, No.l, 1997,
pp.5760.
[68] Thompson, R.C., Junkins, J.L., and Vadali, S.R., \N. iMinimum Time OpenLoop
Slewing of Flexible \; I.! I. Journal of Guidance, Control, and D;;,i.. Vol.12,
No.l, 1989, pp.8288.
[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.358360.
[70] Bass, M., von Stryk, O., Bulirsch, R., and Schmidt, G., "Towards Hybrid Optimal
Control", atAutomatisierungstechnik, Vol.48, No.9, pp.448459, 2000.
[71] von Stryk, O., and Bulirsch, R., "Direct and Indirect Methods for Trajectory
Optimization", Annals of Operations Research, Vol.37, pp.357373, 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 CoastIsland 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.817823.
[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 Paper20076405.
[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 GaussLobatto Quadrature Method for Solving Optimal Control
PiI. in Australian and New Zealand Industrial and Applied Mathematics
Journal, Vol.47, 2006, pp.C101C115.
[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,
pp.270277.
150
100
50 L
1 AC C1
100 C
C3
150  A
0 0.2 0.4 0.6 0.8 1
time
Figure 76: Time history of the control torques.
60
40
20
20
40
0
0
0.2 0.4 0.6 0.8 1
time
Figure 77: Time history of the angles as determined by integration.
Taking the partial derivative of Equation 759 with respect to the generalized coordinates
yields:
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
(760)
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 (761)
which can be written in terms of the inertial coordinate system as:
)5 = h [C]T [C2] b2 0 (762)
where b1 is the unit basis vector of the bodyfixed 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 bodyfixed coordinate frame of body 2 that is perpendicular to
the bodyfixed 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 762 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 (763)
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 (764)

PAGE 1
1
PAGE 2
2
PAGE 3
3
PAGE 4
Iwouldliketothankallofthememberofmycommitteeforthesupport,direction,andconstructivecriticismtheyprovidedinthecompletionofthisresearch.Also,IthankmycolleaguesintheFlightControlLabfortheirsacricesmadebysettingasidetheirownresearchtoaidmeinmywork. 4
PAGE 5
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.2ShortComings ................................. 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
PAGE 6
........................... 100 8CONCLUSIONS ................................... 108 APPENDIX AMASSMATRIXFORINBOARDWINGSEGMENT ............... 111 BMASSMATRIXFOROUTBOARDWINGSEGMENT ............. 116 REFERENCES ....................................... 125 BIOGRAPHICALSKETCH ................................ 133 6
PAGE 7
Figure page 11MorphingMAV. ................................... 16 31ArticulatedMAV. ................................... 51 32Coordinateframesofthefuselage,inboardwingsegment,andoutboardwingsegment. ........................................ 51 33Jointcoordinateframesbetweenfuselageandinboardwingsegment. ...... 52 34Jointcoordinateframesbetweeninboardandoutboardwingsegments. ..... 52 61GeometricmodelofgullwingMAV. ........................ 70 62PositionofMAVwhenusingtheoptimalcontrolsolutionasdeterminedbyGPOCS. 71 63MorphingangleswhenusingtheoptimalcontrolsolutionasdeterminedbyGPOCS. 71 64PositionofMAVfromintegratedoptimalcontrol. ................. 72 65Morphinganglesfromintegratedoptimalcontrol. ................. 72 71Normoftheconstraintvector ............................ 103 72Normoftheconstraintvector ............................ 103 73Threebarmanipulator. ................................ 104 74TimehistoryoftheanglesasdeterminedbyGPOCS. ............... 105 75TimehistoryoftheangularratesasdeterminedbyGPOCS. ........... 105 76Timehistoryofthecontroltorques. ......................... 106 77Timehistoryoftheanglesasdeterminedbyintegration. ............. 106 78Timehistoryoftheangularratesasdeterminedbyintegration. ......... 107 7
PAGE 8
Thisworkpresentsamethodologytodetermineofaccuratecomputermodelsofarticulatedmicroairvehicles.Priortothiswork,muchoftheinnovationindesignofamicroairvehiclewasbasedonempiricaldata.Thisemphasisonempiricaldatanecessitatesatimeintensiveiterativeprocesswhereadesignisbuiltandtested,thenmodicationsaremadeandtheprocessstartsanew.Thisworklookstoalleviatesomeoftheloadofthisiterationbyproducingaccuratemodelsthatcanbeinsertedintotheprocessatgreatlyreducedcostintermsoftime.ThemodelingisaccomplishedbytheapplicationofLagrange'sequationstocreateasetofequationsthatincorporateboththerigidandexibledegreesoffreedomassociatedwiththematerialsthatMicroAirVehicles(MAVs)employ.Theseequationsofmotionalsoincludethegeneralizedforcesthatresultfromthebodyforcesduetogravityandexternalforcesactingasaresultoftheinteractionbetweenthestructureandtheairthroughwhichitpasses.Theseequationsareformulatedinamannerthataccomplishestheabovewhileallowingfortherecongurationofthesystem.Theequationsofmotionarethentransformedintoaformamenableforuseinthesolutionofanoptimalcontrolproblem. 8
PAGE 9
OneofthechallengesinutilizingtheMAVtechnologyisthepropermodelingofvehicles.ThechallengeinmodelingofMAVsisexempliedbythesmallsizesofthevehiclesandtheightregimeinwhichtheyoperate.Whereconventionalaircrafthavebeenstudiedingreatdetailandawealthofknowledgeexistsontheirightregime,thereismuchlessinformationforMAVs.TheinteractionoftheMAVcontrolsurfacesandthelowReynolds'owregimeinwhichtheyoperateisofgreatinterest.ConstructionmaterialsofconventionalaircraftarenotviableforMAVs.Thesmallsizedictatestheuseoflightermaterialbuttheextremeeectthattheightregimehasonthevehiclecallsforbuildingmaterialtoalsobemorefunctionalthanthematerialsinconventionalaircraft.Asaresult,compositestructuresarecommonplaceinMAVconstruction.ThisuseofcompositematerialfurtherstressestheneedfornewmodelingtechniqueswheninvestigatingMAVdynamics. However,anexcitingnewopportunityisprovidedbytheMAVtechnologies.ThesensitivityoftheMAVstotheirenvironmentsallowsforagreaterbenetfromdierentcontrolstrategiesthanconventionalaircraftwouldallow.Forexample,whilemorphingexistsonconventionalaircraft,theextentoftheeectofmorphingoftenpalesincomparisontothecomplexityofthemorphingsystem.OnMAVs,however,morphing 9
PAGE 10
10
PAGE 11
11
PAGE 12
11 .Asopposedtomoretraditionalaircraft,theMAVpresentedingure 11 haswingsthatarearticulatedtoallowfordierentanglesbetweentheinboardwingandthefuselageandalsobetweentheoutboardwingandtheinboardwing.Ratherthanattempttomodelthecompletevehicleinasingleformulation,amultibodyapproach,asdemonstratedbyShabana[ Asidefromthemultibodyapproachtomodeling,thereisanothermatterthataddstothecomplexityoftheformulationoftheequationsofmotion.Thevehiclewillbemodeledwithbothrigidbodyandexibledegreesoffreedom.ForMAVs,thefuselagecanbemodeledasarigidbodyasthereisnotmuchdeformationofthefuselage.Thewingstructure,howeverishighlyexible.Theaeroelasticloadsonthewingscoupledwiththeirlightweightandlowstiness,incomparisontothefuselage,canleadtolargedeformations.Thus,asrstproposedbyMeirovitch[ 12
PAGE 13
Thelevelofdelityofthemodelisalsoaconcern.Intheworkpresentedhere,idealizationswillbemadetosimplifythederivationoftheequationsofmotion.Asanexample,themassesofhingesandactuatorsthatallowforarticulationwillbeignored.Theequationsofmotionthemselveswillalsobeformulatedinaminimalbasisinordertoeliminateconstraintreactionsbetweenbodies.Thesesimplicationsgreatlyreducethecomplexityofthenalsetofequations,withoutneedlesslyoversimplifyinganysimulationsofthevehicleresultingfromtheequationsofmotionderived. Thenalstepinthederivationoftheequationsofmotionistheformulationofthegeneralizedforces.Aerodynamiceectsaccountformostoftheforcesthattheexibleportionsofthemorphingmicroairvehicleexperiences.Theexactnatureoftheseforceshavetobeidealizedforthescopeofthiswork.Theactualforceswhichthevehicleissubjectedisacontinuumthatexistsovertheentiresystem.However,theseforcesandthemomentstheyproducecannotbeaccuratelymeasuredovertheentiresystem.Atthispointintheresearch,itisassumedthataerodynamicforceswillberepresentedbystaticpressuresataparticulartimeatdiscretepointsalongthevehicle.Includedinthisforcemodelingistheassumptionofhowtoconvertpressuredataintoaerodynamicforces.Ultimately,thepressuremeasurementswillbeenoughtoapproximatelyobtaintheactualaerodynamicforcesandmomentsthatcharacterizeofthemorphingMicroAirVehicleatanypointintimeandbasedonthecongurationofthecraft. Thesecondstepinthesolutionoftheproblemstatementiscastingtheequationsofmotioninaformthatisconvenientforcontrol.Fromtheassumptionofaniteelement 13
PAGE 14
Inthefollowingchapters,thecourseofthisresearchwillbefurtherexplored.Chapter 2 presentsabrieflookatthemajortopicsinvolvedinthisresearchandsomeoftheworkthathasbeendoneinthoseareas.Chapter 3 presentsadetailedderivationoftheequationsofmotion.Chapter 4 explainstherepresentationoftheforcesusedintheequationsofmotion.Chapter 5 thencaststheproblemasacontrolsproblemalongwithadetaileddiscussionoftheassumptions,goals,andlimitationsassociatedwithdoingso.Chapter 6 thenappliesthederivationoftheequationsofmotionandthecastingofsaidequationsasacontrolproblemtoanexamplebasedonGullwingMAVillustratedinFig. 11 .Chapter 7 presentsanotherexamplebasedonathreebarmanipulator.Finally,Chapter 8 exploressomeconclusionsofthework. 14
PAGE 15
Thesecondarycontributionofthisworkisinthecontrolstrategiesusedonmicroairvehicles.Theequationsofmotionasderivedinthisresearcharehighlynonlinearandtimedependent.Standardcontrolstrategieswouldchooseanoperatingpointtolinearizeabout,thencreateacontrolstrategybasedonthatpoint,andnallytunetheperformanceforconditionsthatdon'tmatchthelinearizationpoint.Thistypeofcontrolstrategy,evenwhenusingmultiplelinearizationpointsandswitchingcontrollers,doesnotachievethefullperformancepotentialforamorphingMAV.Morphingallowsforaninnitenumberofcongurationswithassociatedperformance.Thus,acontrolstrategythattakesallofthecongurationsintoconsiderationishighlydesirable.Theequationsofmotionasderivedinthisworkcanthenbecastintoanoptimalcontrolproblem.Thesolutionoftheoptimalcontrolproblemnotonlyndsacontrolstrategyforthevehicle,butdependingonthecostfunctionalchosencanpresentaformatthatallowsthemorphingofthevehicletobeusedastheonlymeansbywhichtoobtainsomedesiredmaneuver. 15
PAGE 16
MorphingMAV. 16
PAGE 17
BelowachordReynolds'numberof50,000,thefreeshearlayerafterseparationdoesnottransitiontoturbulentowsoonenoughtoreattachtotheairfoilsurface.[ 17
PAGE 18
TheaerodynamicsofMAVsinlowReynolds'numberregimediersgreatlyfromtheaerodynamicsofminivehicles.[ ThereareseveraldierentclassesofMAVs.InapaperbyMuellerandDelaurier,severaldierentMAVcongurationsarepresented.[ 18
PAGE 19
TheUniversityofFloridahasdevelopedseveralsuccessfulMAVplatforms.TherearesmallMAVsthathaveamaximaldimensionof6in.andtherearelargerMAVsthathaveupto24in.wingspan.BoththelargerandsmallerMAVsdevelopedbytheUniversityofFloridarelyonexiblemembranewingsasopposedtorigidwings.ThedevelopmentofexiblewingsattheUniversityofFloridawasbasedonbiologicalinspiration.BirdsandbatsthatalsoyinthesamelowReynolds'numberregimeasMAVshavethinandundercamberedwings.[
PAGE 20
20
PAGE 21
Severalinitiativeshavebeenstartedbyorganizationstoinvestigatemorphing.AccordingtotheprogramatDARPA,amorphingvehiclewillchangeitsshapesubstantiallytoadapttothemissionenvironment,providesuperiorsystemcapabilitynotpossiblewithoutreconguration,anduseintegrateddesignofmaterials,distributedactuators,eectorsandmechanismstorecongureinight.[ 21
PAGE 22
InworkbySaundersetal.conductedexperimentsona16%scalemodelofanF18wingthatusedshapememoryalloytoeectasmoothdeformationofthetraillingedgecontrolsurfaces.[ InpapersbyAbdulrahimetal.,studieswereperformedtoevaluatetheuseoftorquerods,toeectwingtwist,onthecontrolauthorityofaMAV.[ MAVstypicallyhavewingswithlowaspectratios(AR).[ 22
PAGE 23
23
PAGE 24
24
PAGE 26
31 representsgraphicalcomputermodelofthearticulateMAV.Thisvehicleiscomprisedofacorebody,thefuselage,andfourwingsegments.Thefuselageandeachofthefourwingsegmentshavetheirowncoordinatessystem.Theequationsofmotionforthissystemaredesired.TheseequationsneedtobeabletodescribethecompletelocationandorientationoftheMAVinsomeinertialframe.Sincethegoalistobeabletowritetheequationsofmotionofthisvehiclewithrespecttoaninertialframe,thetransformationmatricesmustbedenedforeachbodythatrelatetheinertialframetothebodyxedframe. Forthepurposesofthederivationoftheequationsofmotion,thevechicleisidealizedasisillustratedinFigure 32 .Thevehicleissymmetricaboutthe^b11^b13plane,thusonlyasinglesideisconsideredintheidealization.Inthefullsimulationofthesystem,asymmetricmodesdooccurandthecontributionfromtheotherhalfoftheMAVcanbeignoredforthepurposesofthederivationofthequationsofmotionwithoutlossofgenerality.Thefollowingderivationswillonlycontainthefuselage(body1),therightinboardwingsegment(body2),andtherightoutboardwingsegment(body3)andtheassociatedcoordinateframes,B1;B2;andB3respectively.Eachbodyxedcoordinateframe,B1;B2;andB3,originateatthecenterofmassofthierrespectivebody.Therelationshipsbetweenthebodyxedcoordinatesystemsofthesebodiesandtheinertialcoordinatesystem,E,areasfollows: ^e ^e ^e 26
PAGE 27
[C1]=266664c1c1s1c1+c1s1s1s1s1+c1s1c1s1c1c1c1+s1s1s1c1s1+s1s1c1s1c1s1c1c1377775(3{4) wherec1iscos1ands1issin1. Inadditiontotheinertialcoordinatesystemandthebodyxedcoordinatesystemsassociatedwitheachbody,thereareintermediatecooridinateframesthatarealignedwiththejointsbetweenbodies.ThesejointcooridnatesystemsareillustratedinFigure 33 andFigure 34 .Therelationshipbetweenthebodiesandthejointcoordinateframesare: 27
PAGE 28
^j ^j 28
PAGE 29
Fortheminimalsetformulationallbodieswillbereferencedtoasinglebody.Inthiscase,allbodieswillbereferencedtotheB1frameassociatedwiththefuselage.TheequationsofmotionwillbederivedusingLagrange'sEquations.Lagrange'sequationsare: dt@T whereTisthekineticenergy,Visthepotentialenergy,RistheRayleighdissipationfunction,q 29
PAGE 30
2v 2v wheretheintegralistakenoverthevolumeofthebodyandv whereR 3{13 yields dtr dt(R dtu dtu SincethefuselageisrigidthepointuopdoesnotchangeintheB1basis.Thus, where[0~!1]isaskewsymmetricmatrixcomposedofthecomponentsofE!B1.Theskewsymmetricmatrixcanbedecomposedas: Thus, 30
PAGE 31
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
PAGE 32
whichisaconstantmatrix. LookingatthersttermofEq. 3{27 whereRVuop;ndVisconstant.Forthesecondandthirdtermssubstitutein0~!1and0~!1whereappropriate. Lookingatthetermintherstrowandrstcolumnoftheresultingmatrixyields whereRVuop;iuop;ndVisconstant.Fortheremainingtermssubstitutein0~!1and0~!1whereappropriate. Nextcreatethekineticenergyoftheinboardwing.Thepositionvectortoanarbitrarypointontheinboardwingisasfollows: wherer 32
PAGE 33
dtr dtr dtr dtu dt[N2]q Byvirtueofsomeofthevectorsbeingconstantintheframethatthederivativeistakenin dt[N2]q Throughtheuseofsimpleangularrotationsbetweenframes, However,B1! .Thus,separatingtheangularvelocitytermsandrewritingallcrossproductsintermsofthesamebasisyields: Let[0~!1]beaskewsymmetricmatrixcomposedofthecomponentsofE!B1and[12~!21]beaskewsymmetricmatrixcomposedofthecomponentsofJ12!J21.Thus, 33
PAGE 34
rewritingallEq. 3{36 intermsoftheB1coordinatesystemyields However,theskewsymmetricmatricescanbefactoredasfollows: Thus, 34
PAGE 35
UsingEquations 3{41 through 3{44 ,Equation 3{40 canbewrittenas: Dene [B]=B 35
PAGE 36
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
PAGE 37
wherer dtr dtr dtr dtr dtr dtu dt[N3]q AswiththeangularvelocityoftheB1framewithrespecttotheinertialframe,theangularvelocityoftheB2framewithrespecttotheinertialframecanbewrittenasthesumofsimpleangularrotations.Thus, However, (3{56) Thus, 37
PAGE 38
Therefore, Expandingthecrossproductsandexpressingthecomponentsofthecrossproductsinthesamebasesyields IntroducingskewsymmetricmatricestosimplifythecrossproductexpressionsandrewritingthewholeequationtobeintermsoftheB1basisyields: 38
PAGE 39
where[0~!1]istheskewsymmetricmatrixcomposedfromthecomponentsofE! Thus, 39
PAGE 40
Dene 40
PAGE 41
[BB]=BB 41
PAGE 42
InsertingEq. 3{70 intotheexpressionforthekineticenergyyields: 2_q where _q and Nowthatthekineticenergyforeachbodyhasbeenderived.Theycanbecombinedintoanexpressionforthekineticenergyfortheentiresystem.Theglobalkineticenergyis 42
PAGE 43
2_q where _q ThenextpieceofinformationneededforLagrange'sequationisthepotentialenergyofthesystem.Aswiththekineticenergyofthesystem,thepotentialenergyofthesystemcanbedecomposedintothepotentialenergyoftheindividualbodiesthatmakeupthesystem.Forthepurposesoftheworkpresentedinthisdocument,thepotentialenergyduetochangesinelevationwillbeignored.Thus,thepotentialenergyofthesystemwillbearesultoftheexibledegreesoffreedom.Asbody1isconsideredrigid,ithasnocontributiontothepotentialenergyofthesystem.Thecontributionstothepotentialenergyofthesystemfrombodies2and3havethesameform.Thearetheresultoftheassumptionoftheniteelementsusedtorepresentthebodies.Thesymbolicformofthepotentialenergycontributionsare: 2f" where"isavectorofthestrainsand[E 43
PAGE 44
AssumethereissomeforceF whereu but Thusforapointonbody1: RewritingEq. 3{81 intermsoftheB1frameyields: 44
PAGE 45
3{82 andthencomparingtheexpressiontoEq. 3{77 ,thecontributionstothegeneralizedforcesfromtheforcesactingonbody1canbederived. Thederivationofthecontributiontothegeneralizedforcesfromtheforcesactingonbody2proceedswithndinganequationforthevirtualdisplacementofapointonbody2.Thepositionvectorofapointonbody2is: WritingEq. 3{83 intermsoftheinertialframeproduces: ThevirtualdisplacementofEq. 3{84 is 45
PAGE 46
Thus RewritingEq. 3{87 intermsoftheB1basisandusingthematricesintroducedinthederivationofthecontributiontothesystemkineticenergyfrombody2: q Performingthedotproductoftheforcesactingonbody2withthecorrespondingvirtualdisplacementvectoroftheformofEq. 3{88 andcomparingtheresultwithEq. 3{77 producesthecontributionstothegeneralizedforcesfromtheforcesactingonbody2. Thederivationofthecontributiontothegeneralizedforcesfromtheforcesactingonbody3proceedsfromtheprincipleofvirtualwork.First,thepositionvectortothepoint 46
PAGE 47
RewritingEq. 3{89 sothatallcomponentsareintermsoftheinertialframeyields Takingthevirtualdisplacementproceedsas:
PAGE 48
Thus, 48
PAGE 49
Usingthevectorsandmatricesdenedinthederivationofthecontributiontothekineticenergyfrombody3,andrewritingthevirtualdisplacementofapointonbody3intermsoftheB1basisyields: q 49
PAGE 50
50
PAGE 51
ArticulatedMAV. Figure32: Coordinateframesofthefuselage,inboardwingsegment,andoutboardwingsegment. 51
PAGE 52
Jointcoordinateframesbetweenfuselageandinboardwingsegment. Figure34: Jointcoordinateframesbetweeninboardandoutboardwingsegments. 52
PAGE 53
Sincethewingstructuresaremodeledasexiblebodies,thegenerationofaeroelasticforcesbecomesimportant.Theaeroelasticforcesusedinthemodelwillinuencetheshapeofthewings,duetotheirexibility,andtheshapeofthewingswillinuencetheaerodynamicforces.Therearemanydierentschemesavailableforthegenerationofaeroelasticforces.Theseschemesrangefrompotentialowtheorytofull3DNavierstokescomputationaluiddynamicssimulations.Thegenerationofaeroelasticforcesisasucientresearchtopicinitsownright.WorkbyOrtizandBarhorstdemonstratesanumberofdierenttechniquesinthegenerationofaerodynamicforcesforuseinaeroelastictystudies.[ 3 53
PAGE 54
AVL'sowanalysisallowsforthecalculationofseveraldierentparameters,suchasstabilityderivativesandhingemoments,baseduponseveralinputsthatdenetheightcondition.Theseinputsincludepitchangle,sideslipangle,pitchrate,rollrate,yawrate,controlsurfacedeectionangle,andwindvelocity.Foragivensetofinputs,AVLdisplaystheforcesandmomentsactinguponabodycenteredcoordinatesystem.Thiscoordinatesystemhasthexaxisrunningalongthefuselageofthemodelandtheyaxisalongtherightwingofthemodel.Theindividualelementforcesand/orstripforcescanbewrittentolesasrequestedbytheuser.AVLalsocanperformseveraldierentoperationsontheinputedmodelsuchasgeometryplotting,eigenmodecalculation,andTretzPlaneplotting. Thus,theaerodynamicforcesandmomentsproducedbyagivenAVLrunassumesthattheightconditionsasdenedbythegeometryandtheinputparametersdonotchange 54
PAGE 55
Inaddition,AVLisonlyapplicableforlowMachnumberow.AnAVLanalysiscanbecompletedwithouttheuseofafuselageinthegeometrydenition.Whenafuselageisincluded,itistreatedasablubodyandasaresulttheforcescalculatedthatactonthefuselagearenotquiteasaccurateasamorerobustanalysistoolwouldgenerate. TheseshortcomingsindicatethattheuseofAVLisnotoptimalforthegenerationofaerodynamicforcesinthesimulationofanunsteadyprocessofaeroelasticityofarecongurableMAV.Instead,morerobustmeasuressuchastheuseofafully3DNavierStokesCFDsoftwarewouldbemoreappropriate.However,theuseofsuchasoftwareonarecongurableMAVisaresearchtopicinitself.Also,thecurrentstateoftheartincomputationaluiddynamicsisverycomputationallyexpensivetogenerateasolutionfortheforcesactingovergivengeometry.Firstanappropriatediscretizationofthegeometryandtheuidsurroundingthegeometrymustbegenerated.Oncethediscretizationhastakenplace,thesolutionoftheaerodynamicforcesisaniterativeprocessandmaytakealargeamountoftimetoconvergetoasinglesolution.Thisproblemiscompoundedwhenusedinanaeroelasticsimulationasthecongurationofthevehiclechangesateverytimestepandthecomputationalcostneededtogeneratetheaerodynamicforcesmustalsobepaidateverytimestep.Incontrast,AVLtakesamatterofsecondstoconvergetoasolutionforagivengeometryinput,duetothesimplicityofthelatticevortexmethod.Thusasacomputationallyinexpensivesourceofaerodynamicforces,AVLisusedeventhoughproperapplicationoftheprocesswouldchooseamorerobustmethod. 55
PAGE 56
Inthepreviouschapter,thecontributionstothekineticenergyandthepotentialenergyfromthevariousbodiesthatcomposethesystemwerederived.Thegeneralizedformofthecontributionstothegeneralizedforcesfromthevariousforcesactingonthesebodieswasalsoderived.Thegoalofthissectionistocombinethesecomponentsintoasystemofequationsthatgovernthemotionoftheentirevehicleandusethisformulationinthederivationofacontrolstrategythatwillbeusedtodrivethevehicletoadesiredstate.Anoptimalcontrolstrategycanbeusedtoaddressanumberofcontrolobjectives,fromtrackingtoregulation. Thereareseveraldierentwaysinwhichtheoptimalproblemcanbesolved.However,thesemethodsgenerallyfallwithintwoclassications.Theseclassicationsaredirectmethodsandindirectmethods.Muchworkhasbeendoneintheuseofindirectmethodsinthesolutionofoptimalcontrolproblems.[ ;u Where(x _x wherex 5{2 showsthatthestatesx 56
PAGE 57
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.Theseequationsformatwopointboundaryvalueproblem,thesolutionofwhichdeterminestheoptimalcontrolu 57
PAGE 58
5{5 5{9 isrequired.Thedicultyinsolvingthoseequationsliesingettingaclosedformsolutionofthecostateequations.Theadvantagesofindirectmethodsisthatthesolutionofthecostateequationsservesasagoodindicatoroftheoptimalityofthesolution.Numerically,indirectmethodsaretypicallysolvedusingshootingormultipleshootingmethods.Theseshootingmethodshavemuchincommonwithstandardrootndingtechniques. Directmethods,alsoknownasdirecttranscription,solvetheoptimalcontrolproblemsbyconvertingtheoptimalcontrolproblemintoanonlinearprogrammingproblem.Thetranscriptionprocessconvertstheinnitedimensionaloptimalcontrolproblemintoanitedimensionaloptimizationproblem.Thedynamicsoftheproblemarediscretizedatanitesetofpointsandthecostfunctionalisreplacedwithacostfunctionthatcanbeevaluatedatthosediscretepoints.ThisnitedimensionalproblemisknownasaNonLinearProgram(NLP).Directmethodshavetheadvantageofnotrequiringthecostateequationstobeexplicitlyderived.Inaddition,thedirectmethodallowsforfasterconvergenceofthesolution.Asaresult,directmethodsareusedquiteofteninpractice,asillustratedinworkbyBass,vonStryk,Hu,Frazzoli,andHull.[
PAGE 59
3 beginswiththeconversionofthesecondorderequationsofmotionintoarstorderformthatisstandardforoptimalcontrolproblems.Thesecondorderequationsofmotionarewrittenas: [M]q Where 1122q 28>>>><>>>>:_q 59
PAGE 60
Thustheequationsofmotioncanbewrittenas: _x ;u where Nowthattheequationsofmotionarewritteninrstorderform,itisnecessarytoexplicitlydenethecontrolproblem.Anexampleofagoalofthecontrolproblemistodrivethevehicletoadesiredpositionusingthecontroltorqueslocatedathingeaxestochangethecongurationofthevehicle.Thischangeofthecongurationwilltheninuencetheaerodynamicsofthevehicle.Writingthisinmathematicaltermsproducesthefollowingobjectivefunction: ThiscostfunctionseekstominimizethedistancebetweenthenalglobalpositionoftheMAVandsomedesirednalglobalposition.TheoptimalcontrolisthenthecontrolthatminimizesEquation 5{16 .Theobjectivefunctionissubjecttophysicalconstraintsonthesystemandmustalsosatisfytheequationsofmotion.Anexampleofphysicalconstraintsistolimittherangeofmotionontheinboardandoutboardsegmentsto90degwithrespecttotheinboardbody(fortheinboardsegmenttheinboardbodyisthefuselageand 60
PAGE 61
minJ=jXCB1XDj2+jYCB1YDj2+jZCB1ZDj2subjectto_x ;u 6 presentsanumericalexampleofanarticulatedMAV.Includedinthisexampleistheformationofanoptimalcontrolproblem. 61
PAGE 62
ThederivationsoftheequationsofmotionthatwerederivedinChapter 3 wereappliedtoanumericalexampletoshowthattheequationscanthenbecastintoaformamenabletoanoptimalcontrolproblem.ThisexampleisimplementedinMATLAB.Therststepinformingthisexampleisgeneratingageometricmodel.Thegeometricmodelisconstructedtoexhibittheprominentfeaturesoftheactualmodel,suchasthearticulationofthewing,butwithoutthecomplexitiesassociatewiththefullmodel.Issuessuchasjointfrictionandcompositematerialpropertiesareignoredwithoutlossofgenerality. 61 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
PAGE 63
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
PAGE 64
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
PAGE 65
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
PAGE 66
5{10 .Nexttheappropriatequantitiesontherighthandsideoftheequationsofmotionmustbedened.Thegravitationalforcesareappliedtoeachmultibodysegmentusingthemassesofeachsegmentandlocalaccelerationofgravity.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
PAGE 67
67
PAGE 68
;u 5{14 andEquation 5{15 Inthedeterminationofthesolutionoftheoptimalcontrolproblem,therearemanydierentcongurationsthatwouldbeselected.Suchalargedesignspaceincreasesthecomputationaleortneededtosolvetheproblem.Toalleviatesomeofthiscomputationalburden,thisexampleisfurtherconstrainedtofollowsymmetricdeectionsonly.Thisconstraintremovestwodegreesoffreedomandtwocontrolinputs,associatedwiththedeectionoftheoppositewing. Adirectmethodwasusedinthesolutionofthisoptimalcontrolproblem.ThismethodusedaGaussPseudospectralmethodinthesolution.ThismethodusestwentyLegendreGausspointsinthediscretizationrepresentthesolutionoftheoptimalcontrolproblem.Thenaltimeforthesystemwasfreetobedeterminedbytheoptimizationprocedurethatminimizedcontroleort.ThisoptimalcontrolproblemwasimplementedintheMATLABprogramusedinthesolutionofstandardoptimalcontrolproblems.[ 68
PAGE 69
Figure 62 showstheglobalpositionofthevehiclegeneratedbytheoptimizationprogramafteritwasunabletoimprovethecostfunctionoverthepreviousiteration'scalculatedcontrolvector.Similarly,Figure 63 showsthemorphanglesgeneratedbytheoptimizationprogram.Theplotsproducedfromtheoptimizationsoftwareseemstoindicatethedesireddecreaseinaltitudeof5mwiththeconstraintofthemorphanglesreturningtozerodeectionatthetimethevehiclecompletesthedecreaseinaltitude.However,theseresultsdonotreectthesuccessfulsolutionoftheoptimalcontrolproblem.ThisfactisillustratedbyFigure 64 andFigure 65 .Theseguresshowthepositionofthevehicleandthemorphanglesafterintegratingtheequationsofmotionusingtheoptimalcontrolasdeterminedbytheoptimizationsoftware.Asillustratedinthegures,notonlydoesthevehiclenotcompletethe5mdecreaseinaltitude,butthemorphanglesarealsononzeroattheterminationofwhattheoptimizationprogramsstatesasthemaneuvertime. 69
PAGE 70
Figure61: GeometricmodelofgullwingMAV. 70
PAGE 71
PositionofMAVwhenusingtheoptimalcontrolsolutionasdeterminedbyGPOCS. Figure63: MorphingangleswhenusingtheoptimalcontrolsolutionasdeterminedbyGPOCS. 71
PAGE 72
PositionofMAVfromintegratedoptimalcontrol. Figure65: Morphinganglesfromintegratedoptimalcontrol. 72
PAGE 73
TheequationsofmotionthatweregeneratedinChapter 3 arealsoapplicabletoothermultibodysystems.Oneofthemostcommonexamplesofamultibodysystemisthethreebarmanipulator.Thissystemconsistsofthreebarsthatareconnectedviasimple,singledegreeoffreedomlinkages.References[ 90 ]and[ 91 ]havesectionsinwhichthiscommonexampleisstudied. 73 ,therstbarisconnectedtoaninertialframeviaarevolutejointthatallowsthebartorotateabouttheaxisalignedwiththebar.Thesecondbarisconnectedtotherstbarbyarevolutejointthatallowsthebartopivotaboutanaxisthatisperpendiculartotheaxisalignedwiththebarandperpendiculartotheplaneillustratedinthegure.Likewise,thethirdbarisconnectedtothesecondbarbyarevolutejointthatallowsthebartopivotaboutanaxisthatisperpendiculartotheaxisalignedwiththebarandperpendiculartotheplaneillustratedinFigure 73 .Attheinitialtimeallofthebodyxedcoordinateframesarealigned. Inthisexamplesomeassumptionsaremadetofacilitatetheexecutions.Primaryamongtheseassumptionsisthemodelingeachbarasarigidbody.ThisassumptioneliminatestheneedforaniteelementrepresentationofthebarsandthuseliminatestheglobalsystemstinessmatrixandthetermsintheapplicationofLagrange'sequationsinthederivationoftheequationsofmotionthatresultfromtheexibledegreesoffreedom.Therevolutejointsconnectingthebarsareassumedtobeideal.Thus,theyhavenomassandonlyallowforasingledegreeoffreedombetweenthebodiestheyconnect.Alsoincludedintheassumptionofidealjointsisthelackoffrictioninthejoints.Itisfurtherassumedthatthebarsareallofonemeterinlengthandarealsoidealizedtohaveaunitarydensityandazeroradius.Finally,itisassumedthatnoexternalforcesactonthesystem. 73
PAGE 74
3 .Inusingthesetwoformulationsthestrengthsandweaknessescanbecontrasted. 74 .Mathematically,thebodyxedcoordinatesystemsarerelatedtotheinertialcoordinatesystembytheuseofatransformationmatrix.Eachtransformationmatrix,isitselfafunctionofthreeangles.Thus, wheree [C1]=266664c1c1s1c1+c1s1s1s1s1+c1s1c1s1c1c1c1+s1s1s1c1s1+s1s1c1s1c1s1c1c1377775(7{4) 74
PAGE 75
[C3]=266664c3c3s3c3+c3s3s3s3s3+c3s3c3s3c3c3c3+s3s3s3c3s3+s3s3c3s3c3s3c3c3377775(7{6) wherec1iscos1ands1issin1.Asallthreebodieshavesimilarrelationsbetweentheirrespectivebodyxedframeandtheinertialframe,theyhavesimilarrelationsfortheangularvelocityofthosebodyxedframeswithrespecttotheinertialframe.Theserelationsare: wheretheangularvelocitiesarewrittenintermsofthebasisvectorsofthebodyxedcoordinatesystem. ThederivationoftheequationsofmotionstartswithapplicationofLagrange'sequations.Forarigidbody,Lagrange'sequationstakestheformof: dt@T WhereTisthekineticenergyofthesystem,q 2v 75
PAGE 76
whereR 7{12 withrespecttotheinertialframeproduces: dtr dt(R dtu where[0~!1]isaskewsymmetricmatrixformedfromthecomponentsoftheangularvelocityofbody1withrespecttotheinertialframe.However, Thus, RewritingEquation 7{15 soallcomponentsarewrittenintermsoftheinertialframeproduces: 76
PAGE 77
Let Therefore SubstitutingEquation 7{21 intoEquation 7{11 produces: 2RTCB1_1_1_1266666664RV[I]dVRVB 7{22 areexpandedas: 77
PAGE 78
ItisimportanttonotethattheintegralsontherighthandsideofEquation 7{23 throughEquation 7{35 areallinvariantintegrals.Forbody1,thegeneralizedcoordinatescanbedenedas: ThusEquation 7{22 takesonthefamiliarquadraticformof: 2_q Thederivationstothekineticenergyofthesystemfrombodies2and3proceedinthesamemannerasthederivationofthecontributionfrombody1illustratedintheabovewith0~!1x1beingreplacedwith1~!2x2forbody2and2~!3x3forbody3.Thesystem'sgeneralizedcoordinatesarecomposedofthegeneralizedcoordinatesoftheindividual 78
PAGE 79
Thusthesystemkineticenergytakesontheformof 2_q SubstitutingEquation 7{39 intoEquation 7{10 resultsin: [M]q whereQ 79
PAGE 80
Thethreebodiesarenotcompletelyfreetomoveintheinertialframe.Throughtheuseoftherevolutejoints,thebodiesareconstrainedtomoverelativetoeachother.Sinceinthestandardderivationoftheequationsofmotionthebodiesarenotdenedrelativetoeachother,anothertermappearsinLagrange'sequationstoaccountfortheconstraints.Lagrange'sequationsbecome: dt@T where[Cq]isthematrixthatresultsfromtakingthederivativeoftheconstraintmatrixwithrespecttothegeneralizedcoordinatesalsoknownasthesystemJacobianmatrix.Arevolutejointonlyallowsforasingledegreeoffreedombetweentwobodies.Thus,eachrevolutejointhasveconstraintequationsassociatedwithit.Thissystemhasthreerevolutejoints,asillustratedinFig. 73 ;onebetweenthegroundandbody1thatallowsrotationabouttheb13axisofthebodyxedcoordinatesystemassociatedwithbody1,onebetweenbody1andbody2thatallowsrotationabouttheb21axisofthebodyxedcoordinatesystemassociatedwithbody2,andonebetweenbody2andbody3thatallowsrotationabouttheb31axisofthebodyxedcoordinatesystemassociatedwithbody3. Theconstraintmatrixisdevelopedbysettingtheconstraintequationsthatgovernthemotionofthebodiesequaltozero.Thus,forthesystem: []=0 (7{43) 80
PAGE 81
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 Theconstraintmatrixforthethreebarmanipulatorisderivedinthefollowing.Therevolutejointbetweenthegroundandbody1restrictsthepositionoftheendofbody1tobecoincidentwiththepointonthegroundthatcorrespondstotherevolutejointlocation. 81
PAGE 82
(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
PAGE 83
7{53 withrespecttothegeneralizedcoordinatesproduces: Similarly, 3=e Orpresentingallquantitiesintermsoftheinertialframe: 3=e TakingthepartialderivativeofEquation 7{56 withrespecttothegeneralizedcoordinatesproduces: Fortherevolutejointbetweenbody1andbody2theconstraintontherelativetranslationalmotionbetweenthebodiesisexpressedas: (7{58) whereR (7{59) 83
PAGE 84
7{59 withrespecttothegeneralizedcoordinatesyields: 0 0 Nowlookingattheconstraintsthesecondrevolutejointimposesontheangulardegreesoffreedombetweenbody1andbody2 5=b whichcanbewrittenintermsoftheinertialcoordinatesystemas: 5=b whereb 7{62 is: Similarly,thesecondconstraintonrotationaldegreesoffreedomimposedbytherevolutejointbetweenbody1andbody2isillustratedin: 6=b 84
PAGE 85
7{64 intermsoftheinertialcoordinatesystemproduces: 6=b TakingthepartialderivativeofEquation 7{65 withrespecttothegeneralizedcoordinatesyields: Theconstraintsontherelativetranslationalmotionbetweenbody2andbody3causedbythethirdrevolutejoint,areillustratedin (7{67) whereR 7{67 intermsoftheinertialcoordinatesystemyields: (7{68) 85
PAGE 86
7{69 withrespecttothegeneralizedcoordinatesproduces: 0 0 Lookingattheconstraintsontherelativerotationalmotionbetweenbody2andbody3imposedbytherevolutejoint: 8=b whereb 8=b ThecontributiontotheJacobianmatrixfromthisconstraintis: Similarly,thesecondrotationalconstraintimposedbytherevolutejointbetweenbody2andbody3isexempliedin: 9=b 86
PAGE 87
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
PAGE 88
whichgivesthesimultaneoussolutionforthesecondtimederivativeofthegeneralizedcoordinatesandtheLagrangemultipliers.Thesolutionofthesecondtimederivativeofthegeneralizedcoordinatescanthenbeusedinarstordersystemtosolveforthegeneralizedcoordinatesatthenexttimestep.Towritetheequationsofmotionasarstordersystem,rstlet: Thus, _X whereu 71 demonstratesthesolutionoftheequationsofmotiontosimulatethesystemafteronesecondinwhichacontrolisappliedto1.Thereisaknownproblemwiththedirectintegrationoftheequationsofmotion.ThisproblemisthatthesolutionoftheLagrangemultipliersisnotexact.Asthesimulationtimegrows,theseerrorsalsogrow.Asseeninthegure,thenormofthevectorthatrepresentstheconstraintvectorforthesystemgrowsasthesimulationprogresses.ThereisastabilizationmethodthatreducestheerrorsassociatedwiththesolutionoftheaugmentedsystemillustratedinEquation 7{79 .ThisprocesswasrstpublishedintheworkbyBaumgarte[ +2_ +2 =0 (7{83) 88
PAGE 89
7{83 lieinthelefthalfplane,thesolutionisguaranteedtoconvergetozero.Thus,withtheproperchoiceofandinEquation 7{83 theconstraintvectorandtherstandsecondtimederivativesoftheconstraintvectorcanbemadetoconvergetozero.[ 72 showsthecomparisonbetweenthenormoftheconstraintvectorwithoutstabilizationandwithstabilization.Asillustratedinthegure,withoutconstraintstabilizationtheerrorintheconstraintvectorgrowsexponentiallywithtime.Stabilizationhelpstoreducetheerrorsaccumulatedandtherateatwhichtheyaccumulate.However,asillustratedthesimultaneoussolutionoftheLagrangemultipliersandtheaccelerationsthatisnecessarywhenusingthestandardformulationoftheequationsofmotionstillleadstoviolationsoftheconstraintsonthesystem.Forthisreason,formingtheequationsofmotionwithouttheneedtocoupleconstraintsequationsisbenecialintheoverallsimulationofthesystem. 3 allowsfortheeliminationofsomeofthesediculties.Thenumberofgeneralizedcoordinatesismuchlesswiththisalternativeformulation.Sinceallbodiesinthesystemarereferencedtoa\central"body,onlythefullcoordinatesassociatedwiththisbodyarenecessary.Thosecoordinatesarethenaugmentedwithcoordinatesthatrelatethedierencesinorientationbetweensubsequentbodiesandthe\central"body. Thederivationoftheequationsofmotion,usingthealternativeformulation,forthethreebarsystembeginswiththeestablishmentofcoordinateframes.Asinthestandard 89
PAGE 90
wheree [C]=266664cos1sin10sin1cos10001377775(7{85) InthealternativeformulationillustratedinChapter 3 showedjointframesthatwereusedasthebasisforrelatingthebodyxedcoordinatesystemofonebodytothebodyxedcoordinatesystemofanadjoiningbody.Forthisexamplethesejointframesarenotnecessaryasthebodyxedcoordinateframesarealignedwiththejointaxisoftherevolutejointsforallbodies.Thustherelationshipbetweentherstbodyandthesecondbodycanbestatedas: where [T(2)]=2666641000cos2sin20sin2cos2377775(7{87) andr 90
PAGE 91
where [T(3)]=2666641000cos3sin30sin3cos3377775(7{89) andr whereb whereb 91
PAGE 92
whereRCB1iswrittenintermsoftheinertialcoordinatesystemandu 7{93 ,withrespecttotheinertialcoordinateframe,yields: dtr dt(R dtu buttherevolutejointbetweenthegroundandbody1preventstherigidbodytranslationofbody1,thusR where[0~!1]istheskewsymmetricmatrixformedfromE! Thus, Thusthecontributiontothekineticenergyfrombody1canbewrittenas: 2ZVv 2ZVu 92
PAGE 93
2_qM_q(7{99) where_q=_1andusingtheindicialnotation, Thederivationofthecontributiontothekineticenergyofthesystemfrombody2beginswiththedenitionofthevectorfromtheoriginoftheinertialsystemtoarandompointonbody2.Thisvectoriswrittenas: ,wherer 7{101 withrespecttotheinertialcoordinateframeyields: dtr dt(R dtr dtr dtu but 93
PAGE 94
Theskewsymmetricmatricescanbefactoredtoyield: Let [B]=B Thusthevelocityofarandompointonbody2iswrittenas: InsertingEquation 7{109 intothecontributiontothesystem'skineticenergyfrombody2yields: 2ZV[B]T[B]dV8><>:_1_29>=>;(7{110) Thus,thecontributiontothekineticenergyfrombody2assumesthefamiliarquadraticform: 2_q 94
PAGE 95
[M]=ZV[B]T[B]dV=264RVB LookingatarepresentativetermofEquation 7{112 andusingtheindicialnotation, ItisimportanttonotethattheintegralsleftinEquation 7{113 areconstantsintegrals. Thederivationofthecontributiontothesystem'skineticenergyfrombody3beginswiththeformationofapositionvectorfromtheoriginoftheinertialsystemtoarandompointonbody3.Thisvectorisstatedas: wherer 95
PAGE 96
7{115 withrespecttotheinertialcoordinatesystemyields: dtr dt(R dtr dtr dtr dtr dtu However, InsertingEquation 7{116 andEquation 7{117 intoEquation 7{115 andfactoringtheskewsymmetricmatricesproduces: Dene 96
PAGE 97
Therefore, InsertingEquation 7{122 intothecontributiontothesystem'skineticenergyfrombody3yieldsthefamiliarquadraticform: 2_q where _q=8>>>><>>>>:_1_2_39>>>>=>>>>;(7{124) [M]=266664RVBB ExaminingarepresentativetermofEquation 7{125 andusingtheindicialnotationyields: 97
PAGE 98
Hereagain,theremainingintegralsinEquation 7{126 areinvariantoverthebody. Combiningthekineticenergycontributionsfrombody1,body2,andbody3yieldsthekineticenergyoftheentiresystem.Then,insertingthiskineticenergytermintotheLagrange'sequationsproducestheequationsofmotionofthesystemwhichcanbestated 98
PAGE 99
[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.TheoptimalcontrolproblemdenedintheabovewastranscribedintoaNonLinearProgram(NLP)usingaGaussPseudospectralmethod.Figures 74 and 75 showthetimehistoryoftheanglesandtheangularratesundertheinuenceoftheoptimalcontrolillustratedingure 76 .Thesegureswereproducedbytheoptimizationsoftware,GPOCS,aftertheoptimizationwasterminatedbecauseofalackofimprovementintheminimizationofthecostfunctionoverpreviousiterations.Althoughtheequationsofmotionwereformulatedinawaythatwouldallowforthesolutionoftheoptimalcontrolproblem,theresultsdonotdemonstrateasuccessfulsolution.AsillustratedinFigure 77 andFigure 78 ,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 Figure72: Normoftheconstraintvector 103
PAGE 104
Threebarmanipulator. 104
PAGE 105
TimehistoryoftheanglesasdeterminedbyGPOCS. Figure75: TimehistoryoftheangularratesasdeterminedbyGPOCS. 105
PAGE 106
Timehistoryofthecontroltorques. Figure77: Timehistoryoftheanglesasdeterminedbyintegration. 106
PAGE 107
Timehistoryoftheangularratesasdeterminedbyintegration. 107
PAGE 108
ThisresearchwasconcernedwiththeemergenceofMicroAirVehiclesandthevariousnewopportunitiesthattheyprovide.OneofthemostcompellingistheabilitytousemorphingonMAVs.ThesmallstructureandlightweightmaterialsallowsformorphingstrategiesthatcancompletelyalterthecongurationofMAVs.Thesecongurationchangescanthenbeusedtoexploitightregimesandmaneuversthatarealientoconventionalaircraft.However,withthesenewopportunitiescomenewobstaclesthatmustbeovercome.ThechiefamongtheseobstaclesisthelackofawaytocharacterizetheperformanceofMAVsbasedonconventionalaircraftmodeling.Newmodelingtechniquesarerequiredthatallowforthe,insomecases,completerecongurationofthevehiclebecauseintheightregimethatMAVstypicallyoperateinthechangeshavegreatimpactontheperformance.Inanattempttocreatenewmodelingtechniques,thisstudyoutlinedthederivationofequationsofmotionbasedonmultibodydynamicsthatallowsforthelargecongurationchangespossiblewithmorphing.Theseequationsofmotionalsoincorporatethecouplingofrigidbodyandexiblebodydynamicsintoasinglesetofequations.Alsopresentedwastheabilitytotaketheequationsofmotionasderivedandcastthemintoaformamenableforanoptimalcontrolproblem.ThenumericalexampledemonstratedtheimplementationofthemultibodyequationofmotionformulationonamodelthatsimulatesthecongurationchangespossibleonMAVs.Thethreebarexampleshowsthatthemethodsderivedareapplicabletotraditionaldynamicsystemsaswellandthroughtheoptimalcontrolproblem,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,forwhichitiswellsuited.Perhapstherichestresearchtopicthatextendsfortheworkpresentedhereliesintheareaofoptimalcontrol.Apossibletopicforresearchliesintheuseofsomesymbolicmanipulationtocastthehighlystatedependentequationsofmotionintoaformthatallowsfortheanalyticalsolutionofthecostateequationsandthusanapplicationofindirectmethods.AsdemonstratedinChapter 6 andChapter 7 directmethodsdonotguaranteethesuccessfulsolutionoftheoptimalcontrolproblem.Moreresearchcanbedonetoeithertransformtheequationsofmotion,asderivedinthisresearch,intoaformthattheoptimizationsoftwareGPOCScansuccessfullyndanoptimumforortondadirectmethodthatismoreappropriate. TherearealsotopicsforfurtherresearchintheapplicationofthisworktoMAVsinspecic.Determiningtheinteractionbetweenthestructuresandtheuidthroughwhich 109
PAGE 110
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
Itisimportanttonotethatalloftheintegralsontherighthandsideoftheaboveequationareinvariant.FortheremainingtermsofRV[B]T[B]dVreplace0~!1;with0~!1;
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
Examiningarepresentativeterm:
PAGE 119
Examiningarepresentativeterm:
PAGE 120
Examiningarepresentativeterm:
PAGE 121
121
PAGE 122
122
PAGE 123
123
PAGE 124
Finally, Notealloftheintegralsleftintheindicialnotationareconstant. 124
PAGE 125
[1] Shabana,A.A.,DynamicsofMultibodySystemsNewYorkNY:JohnWiley&Sons,1989. [2] Meirovitch,L.andStemple,T.,\HybridStateEquationsofMotionforFlexibleBodiesinTermsofQuasiCoordinates",JournalofGuidance,Control,andDynamics,Vol.18,NO.4,1995,pp.678688. [3] Meirovitch,L.andTuzcu,I.,\IntegratedApproachtotheDynamicsandControlofManeuveringFlexibleAircraft",NASATechnicalPaper,PaperNo.NASA/CR2003211748,2003. [4] Meirovitch,L.andTuzcu,I.,\UniedTheoryfortheDynamicsandControlofManeuveringFlexibleAircraft",AIAAJournal,Vol.42,No.4,April2004,pp.714727. [5] Meirovitch,L.andTuzcu,I.,\TimeSimulationsoftheResponseofManeuveringFlexibleAircraft",JournalofGuidance,Control,andDynamics,Vol.27,No.5,2004,pp.814828. [6] Kwak,M.K.andMeirovitch,L.,\NewApproachtotheManeuveringandControlofFlexibleMultibodySystems",JournalofGuidance,Control,andDynamics,Vol.15,No.6,1992,pp.13421353. [7] GadelHak,M.,\MicroAirVehicles:CanTheybeControlledBetter?",JournalofAircraft,Vol.38,No.3,2001,pp.419429. [8] Lian,Y.,ShyyW.,Viieru,D.,andZhang,B.,\MembraneWingAerodynamicsforMicroAirVehicles",ProgressinAerospaceSciences,Vol.39,2003,pp.425465. [9] Mueller,T.J.andDeLaurier,J.D.,\AerodynamicsofSmallVehicles",AnnualReviewofFluidMechanics,Vol.35,2003,pp.89111. [10] MutiLin,J.C.andPauley,L.L.,\LowReynoldsNumberSeparationonanAirfoil",AIAAJournal,Vol.34,No.8,1996,pp.15701577. [11] Lutz,T.,Wurz,W.,andWagner,S.,\NumericalOptimizationandWindTunnelTestingofLowReynoldsNumberAirfoils",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.,\FlexibleWingBasedMicroAirVehicles"AmericanInstituteofAeronauticsandAstronauticsPaperNo.AIAA20020705,2002. 125
PAGE 126
Hall,J.,Mohseni,K.,Lawrence,D.,andGeuzaine,P.,\InvestigationofVariableWingSweepforApplicationsinMicroAirVehicles",AmericanInstituteofAeronauticsandAstronauticsPaperNo.AIAA20057171,2005. [15] Ramamurti,R.andSandberg,W.,\SimulationoftheDynamicsofMicroAirVehicles",38thAerospaceSciencesMeeting&ExhibitReno,Nevada,PaperNo.AIAA20000896,2000. [16] Raney,D.L.andSlominski,E.C.,\MechanizationandControlConceptsforBiologicallyInspiredMicroAirVehicles",JournalofAircraft,Vol.41,No.6,2004,pp.12571265. [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.AIAA20041674,2004. [19] Waszak,M.R.,Davidson,J.B.,andIfju,P.G.,\SimulationandFlightControlofanAeroelasticFixedWingMicroAirVehicle",AIAAAtmosphericFlightMechanicsConference,Monterey,CA,PaperNo.AIAA20024875,2002. [20] Waszak,M.R.,Jenkins,L.N.,andIfju,P.,\StabilityandControlPropertiesofanAeroelasticFixedWingMicroAirVehicle",AIAAAtmosphericFlightMechanicsConference,Montreal,Canada,PaperNo.AIAA20014005,2005. [21] Lian,Y.,Shyy,W.,andHaftka,R.T.,\ShapeOptimizationofaMembraneWingforMicroAirVehicles",AIAAJournal,Vol.42,No.2,2003,pp.424426. [22] Lian,Y.,Shyy,W.,Ifju,P.G.,andVerron,E.,\MembraneWingModelofMicroAirVehicles",AIAAJournal,Vol.41,No.12,2003,pp.24922494. [23] Lian,Y.andShyy,W.,\NumericalSimulationsofMembraneWingAerodynamicsforMicroAirVehicleApplications",JournalofAircraft,Vol.42,No.4,2005,pp.865873. [24] Ramrakhyani,D.S.,Lesieutre,G.A.,Frecker,M.,andBharti,S.,\AircraftStructuralMorphingUsingTendonActuatedCompliantCellularTrusses",JournalofAircraft,Vol.42,No.6,2005,pp.16151621. [25] Bae,J.,Seigler,T.M.,andInman,D.J.,\AerodynamicandStaticAeroelasticCharacteristicsofVariableSpanMorphing",JournalofAircraft,Vol.42,No.2,2005,pp.528534. 126
PAGE 127
Livne,E.andWeisshaar,T.A.,\AeroelasticityofNonConventionalAirplaneCongurations{PastandFuture",JournalofAircraft,Vol.40,No.6,2003,pp.10471065. [27] Wlezien,R.W.,Horner,G.C.,McGowan,A.R.,Padula,S.L.,Scott,M.A.,Silcox,R.J.,andSimpson,J.O.,\TheAircraftMorphingProgram",AmericanInstituteofAeronauticsandAstronauticsPaperNo.AIAA981927,1998. [28] Reich,G.W.,Raveh,D.E.,andZink,P.S.,\ApplicationofActiveAeroelasticWingTechnologytoaJoinedWingSensorcraft",JournalofAircraft,Vol41,No.3,2004,pp.594602. [29] Zheng,L.andRamaprian,B.R.,\APiezoElectricallyActuatedWingforaMicroAirVehicle",AmericanInstituteofAeronauticsandAstronauticsPaperNo.AIAA20022302,2002. [30] Saunders,B.,Eastep,F.E.,andForster,E.,\AerodynamicandAeroelasticCharacteristicsofWingswithConformalControlSurfaces",JournalofAircraft,Vol.40,No.1,2003,pp.9499. [31] Abdulrahim,M.,Garcia,H.,andLind,R.,\FlightCharacteristicsofShapingtheMembraneWingofaMicroAirVehicle",JournalofAircraft,Vol.42,No.1,2005,pp.131137. [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.105121. [35] Kral,R.andKreuzer,E.,\MultibodySystemsandFluidStructureInteractionswithApplicationtoFloatingStructures",MultibodySystemDynamics,Vol.8,1999,pp.6583. [36] Shabana,A.A.,\FlexibleMultibodyDynamics:ReviewofPastandRecentDevelopment",MultibodySystemDynamics,Vol.1,1997,pp.189222. [37] Schiehlen,W.,\ResearchTrendsinMultibodySystemDynamics",MultibodySystemDynamics,Vol.18,2007,pp.313. [38] Yoo,W.,Kim,K.,Kim,H.,andSohn,J.,\DevelopmentsofMultibodySystemDynamics:ComputerSimulationsandExperiments",MultibodySystemDynamics,Vol.18,2007,pp.3558. 127
PAGE 128
Boutin,B.A.andMisra,A.K.,\DynamicsofManipulatingTrussStructures",CollectionofTechnicalPapers,AIAA/AASAstrodynamicsConference,SanDiego,CA,PaperNo.A96347120912,July1996,pp.4714793. [40] deJalon,J.G.,\TwentyFiveYearsofNaturalCoordinates",MultibodySystemDynamics,Vol.18,2007,pp.1533. [41] Leger,M.andMcPhee,J.,\SelectionofModelingCoordinatesforForwardDynamicMultibodySimulations",MultibodySystemDynamics,Vol.18,2007,pp.277297. [42] Anderson,K.S.,andDuan,S.,\HighlyParallelizableLoworderDynamicsSimulationAlgorithmforMultiRigidBodySystems",JournalofGuidance,Control,andDynamics,Vol.23,No.2,2000,pp.355364. [43] Banerjee,A.K.,\ContributionsofMultibodyDynamicstoSpaceFlight:ABriefReview",JournalofGuidance,Control,andDynamics,Vol.26,No.3,2003,pp.385394. [44] ChenS.,HansenJ.M.,Tortorelli,D.A.,\UnconditionallyEnergyStableImplicitTimeIntegration:ApplicationtoMultibodySystemAnalysisandDesign",InternationalJournalforNumericalMethodsinEngineering,Vol.48,2000,pp.791822. [45] Chen,S.,Tortorelli,D.A.,\OptimizationofMultibodySystemsViaanEnergyConservingMinimalCoordinateFormulation",AmericanInstituteofAeronauticsandAstronautics,PaperNo.AIAA984927,1998. [46] Schaefer,E.,Loesch,M.,Rulka,W.,\RealTiemSimulationofFlexibleSpaceManipulatorDynamics",AmericanInstituteofAeronauticsandAstronautics,PaperNo.AIAA9737062,1997. [47] Kielau,G.andMaiber,P.,\NonholonomicMultibodyDynamics",MultibodySystemDynamics,Vol.9,2003,pp.213236. [48] Aoustin,Y.andFormalsky,A.,\OntheFeedforwardTorquesandReferenceTrajectoryforFlexibleTwoLinkArm",MultibodySystemDynamics,Vol.3,1999,pp.241265. [49] Junkins,J.L.,\AdventuresontheInterfaceofDynamicsandControl",AIAA,35thAerospaceSciencesMeeting&Exhibit,Reno,NV,Jan.69,1997. [50] Rives,J.,Portigliotti,S.,andLeveugle,T.,\AtmosphericReEntryDemonstratorDescentandRecoverySubSystemQualicationTestFlightDynamicsandTrajectoryModelingduringDescent",AmericanInstitueofAeronauticsandAstronautics,PaperNo.AIAA971441,1997. [51] Tadikonda,S.S.K.,\ModelingofTranslationalMotionBetweenTwoFlexibleBodiesConnectedviaThreePoints",JournalofGuidance,Control,andDynamics,Vol.18,No.6,1995,pp.13921397. 128
PAGE 129
Tadikonda,S.S.K.,\Articulated,FlexibleMultibodyDynamicsModeling:GeostationaryOperationalEnvironmentalSatellite8CaseStudy",JournalofGuidance,Control,andDynamcis,Vol.20,No.2,1997,pp.276283. [53] Jain,A.,andRodriguez,G.,\RecursiveDynamicsAlgorithmforMultibodySystemswithPrescribedMotion",JournalofGuidance,Control,andDynamics,Vol.16,No.5,1993,pp.830837. [54] Pradhan,S.,Modi,V.J.,andMisra,A.K.,\OrderNFormulationforFlexibleMultibodySystemsinTreeTopology:LagrangianApproach",JournalofGuidance,Control,andDynamics,Vol.20,No.4,1997,pp.665672. [55] Attia,H.,\DynamicSimulationofConstrainedMechanicalSystemsUsingRecursiveProjectionAlgorithm",JournalofBrazilianSocietyofMechanicalSciencesandEngineering,Vol.28,2006,pp.3744. [56] Dias,J.M.P.andPereira,M.S.,\SensitivityAnalysisofRigidFlexibleMultibodySystems",MultibodySystemDynamics,Vol.1,1977,pp.303322. [57] Ortiz,J.L.andBarhorst,A.A.,\ModelingFluidStructureInteraction",JournalofGuidance,Control,andDynamics,Vol.20,No.6,1997,pp.12211228. [58] Paiewonsky,B.,\OptimalControl:AReviewofTheoryandPractice",AIAAJournal,Vol.3,No.11,1965,pp.19852006. [59] Hodges,D.H.andBless,R.R.,\WeakHamiltonianFiniteElementMethodforOptimalControlProblems",JournalofGuidance,Control,andDynamics,Vol.14,No.1,1991,pp.148156. [60] Schaub,H.andJunkins,J.L.,\NewPenaltyFunctionsandOptimalControlFormulationforSpacecraftAttitudeControlProblems",JournalofGuidance,Control,andDynamics,Vol.20,No.3,1997,pp.428434. [61] Seywald,H.andKumar,R.R.,\FiniteDierenceSchemeforAutomaticCostateCalculation",JournalofGuidance,Control,andDynamics,Vol.19,No.1,1996,pp.231239. [62] Seywald,H.andKumar,R.R.,\AMethodforAutomaticCostateCalculation",AIAA,Guidance,Navigation,andControlConference,SanDiego,CA,PaperNo.AIAA963699,1996. [63] Seywald,H.andKumar,R.R.,\MethodforAutomaticCostateCalculation",JournalofGuidance,Control,andDynamics,Vol.19,No.6,1996,pp.12521261. [64] Venkataraman,P.,\DynamicsandOptimalControl",AIAAGuidance,Navigation,andControlConferenceandExhibit,Dever,CO,PaperNo.AIAA20004560,2000. 129
PAGE 130
Vadali,S.R.andSharma,R.,\OptimalFiniteTimeFeedbackControllersforNonlinearSystemswithTerminalConstraints",JournalofGuidance,Control,andDynamics,Vol.29,No.4,2006,pp.921928. [66] Betts,J.T.,\SurveyofNumericalMethodsforTrajectoryOptimization",JournalofGuidance,Control,andDynamics,Vol.21,No.2,1998,pp.193207. [67] Hull,D.G.,\ConversionofOptimalControlProblemsintoParameterOptimizationProblems",JournalofGuidance,Control,andDynamics,Vol.20,No.1,1997,pp.5760. [68] Thompson,R.C.,Junkins,J.L.,andVadali,S.R.,\NearMinimumTimeOpenLoopSlewingofFlexibleVehicles",JournalofGuidance,Control,andDynamics,Vol.12,No.1,1989,pp.8288. [69] Warner,M.S.,andHodges,D.H.,\TreatmentofControlConstraintsinFiniteElementSolutionofOptimalControlProblems",JournalofGuidance,Control,andDynamics,Vol.22,No.2,1999,pp.358360. [70] Bass,M.,vonStryk,O.,Bulirsch,R.,andSchmidt,G.,\TowardsHybridOptimalControl",atAutomatisierungstechnik,Vol.48,No.9,pp.448459,2000. [71] vonStryk,O.,andBulirsch,R.,\DirectandIndirectMethodsforTrajectoryOptimization",AnnalsofOperationsResearch,Vol.37,pp.357373,1992. [72] Hu,G.S.,Ong,C.J.,andTeo,C.L.,\DirectCollocationandNonlinearProgrammingforOptimalControlProblemUsinganEnhancedTranscribingScheme",Proceedingofthe1999IEEEInternationalSymposiumonComputerAidedControlSystemDesign,KohalaCoastIslandofHawaii,Hawaii,USA,1999. [73] Frazzoli,E.andBulle,F.,\OnQuantizationandOptimalControlofDynamcalSystemswithSymmetries",Proceedingofthe41stIEEEConferenceonDecisionandControl,LasVegas,NV,2002,pp.817823. [74] Huntington,G.,Benson,D.,andRao,A.,\AComparisonofAccuracyandComputationalEciencyofThreePseudospectralMethods"2007AIAANavigationandControlConverenceHiltonHead,SC,2007,AIAAPaper20076405. [75] Huntington,G.andRao,A.,\AComparisonofLocalandGlobalOrthogonalCollocationMethodsforSolvingOptimalControlProblems",2007AmericanControlConference,NewYork,NY,2007. [76] Williams,P.,\AGaussLobattoQuadratureMethodforSolvingOptimalControlProblems",AustralianandNewZealandIndustrialandAppliedMathematicsJournal,Vol.47,2006,pp.C101C115. [77] Fahroo,F.andRoss,I.M.,\CostateEstimationbyaLegendrePseudospectralMethod",JournalofGuidance,Control,andDynamics,Vol.24,No.2,2001,pp.270277. 130
PAGE 131
Herman,A.L.,andConway,B.A.,\DirectOptimizationUsingCollocationBasedonHighOrderGaussLobattoQuadratureRules",JournalofGuidance,Control,andDynamics,Vol.19,No.3,1996,pp.592599. [79] Lo,J.,Huang,G.,andMetaxas,D.,\HumanMotionPlanningBasedonRecursiveDynamicsandOptimalControlTheory",MultibodySystemDynamics,Vol.1,1997,pp.189222. [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.116. [83] Weijia,Z.,Zhenkuan,P.,andYibing,W.,\AnAutomaticConstrainViolationStabilizationMethodforDierential/AlgebraicEquationsofMotioninMultibodySystemDynamics",AppliedMathematicsandMechanics,Vol.21,2000,pp.103108. [84] Ascher,U.,Chin,H.,Petzold,L.,andReich,S.,\StabilizationofConstrainedMechanicalSystemswihtDAEsandInvariantManifolds",JournalofMechanicsofStructuresandMachines,Vol.23,1995,pp.135158. [85] Lin,S.andHong,M.,\StabilizationMethodfortheNumericalIntegrationofControlledMultibodyMechanicalSystem:AHybridIntegrationApproach",TheJapanSocietyofMechanicalEngineersInternationalJournal,Vol.44,2001,pp.7988. [86] Thomas,S.,\DynamicsofSpacecraftandManipulators",SIMULATION,Vol.57,1991,pp.5672. [87] Lin,S.andHuang,J.,\StabilizationofBaumgarte'sMethodUsingtheRungeKuttaApproach",JournalofMechanicalDesign,Vol.124,2002,pp.633641. [88] Anderson,K.andCritchley,J.,\Improved'OrderN'PerformanceAlgorithmfortheSimulationofConstrainedMultiRigidBodyDynamicSystems",MultibodySystemsDynamics,Vol.9,2003,pp.185212. [89] Simeon,B.\OnLagrangeMultipliersinFlexibleMultibodyDynamics",ComputerMethodsAppliedMechanicalEngineering,Vol.195,2006,pp.69937005. [90] Bottasso,C.L.andCroce,A.,\OntheSolutionofInverseDynamicsandTrajectoryOptimizationProblemsforMultibodySystems",MultibodySystemDynamics,Vol.11,2004,pp.122. 131
PAGE 132
Bottasso,C.L.andCroce,A.,\OptimalControlofMultibodySystemsUsinganEnergyPreservingDirectTranscriptionMethod",MultibodySystemDynamics,Vol.12,2004,pp.1745. [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
