Abstract
 Title of dissertation: MODELING AND SYSTEM
 IDENTIFICATION OF AN ORNITHOPTER
 FLIGHT DYNAMICS MODEL
 Jared Grauer
 Doctor of Philosophy, 2012
 Dissertation directed by: Langley Distinguished Professor
 James E. Hubbard, Jr.
 Department of Aerospace Engineering
 Ornithopters are robotic  ight vehicles that employ  apping wings to generate
 lift and thrust forces in a manner that mimics avian  yers. At the small scales and
 Reynolds numbers currently under investigation for miniature aircraft where viscous
 e ects deteriorate the performance of conventional aircraft, ornithopters achieve e -
 cient  ight by exploiting unsteady aerodynamic  ow  elds, making them well-suited
 for a variety of unmanned vehicle applications. Parsimonous dynamic models of these
 systems are requisite to augment stability and design autopilots for autonomous op-
 eration; however,  apping  ight is fundamentally di erent than other means of engi-
 neered  ight and requires a new standard model for describing the  ight dynamics.
 This dissertation presents an investigation into the  ight mechanics of an ornithopter
 and develops a dynamical model suitable for autopilot design for this class of system.
 A 1.22 m wing span ornithopter test vehicle was used to experimentally investigate
  apping wing  ight. Flight data, recorded in trimmed straight and level mean  ight
 using a custom avionics package, reported pitch rates and heave accelerations up to
 5.62 rad/s and 46.1 m/s2 in amplitude. Computer modeling of the vehicle geometry
 revealed a 0.03 m shift in the center of mass, up to a 53.6% change in the moments
 of inertia, and the generation of signi cant inertial forces. These  ndings justi ed
 a nonlinear multibody model of the vehicle dynamics, which was derived using the
 Boltzmann-Hamel equations. Models for the actuator dynamics, tail aerodynamics,
 and wing aerodynamics, di cult to obtain from  rst principles, were determined
 using system identi cation techniques with experimental data. A full nonlinear  ight
 dynamics model was developed and coded in both Matlab and Fortran programming
 languages.
 An optimization technique is introduced to  nd trim solutions, which are de ned
 as limit cycle oscillations in the state space. Numerical linearization about straight
 and level mean  ight resulted in both a canonical time-invariant model and a time-
 periodic model. The time-invariant model exhibited an unstable spiral mode, stable
 roll mode, stable dutch roll mode, a stable short period mode, and an unstable short
period mode. Floquet analysis on the identi ed time-periodic model resulted in an
 equivalent time-invariant model having an unstable second order and two stable  rst
 order modes, in both the longitudinal and lateral dynamics.
MODELING AND SYSTEM IDENTIFICATION OF AN
 ORNITHOPTER FLIGHT DYNAMICS MODEL
 by
 Jared A. Grauer
 Dissertation submitted to the Faculty of the Graduate School of the
 University of Maryland, College Park in partial ful llment
 of the requirements for the degree of
 Doctor of Philosophy
 2012
 Advisory Committee:
 Professor James E. Hubbard, Jr., Chair/Advisor
 Professor Amr Baz
 Professor J. Sean Humbert
 Professor Darryll Pines
 Professor Robert Sanner
c Copyright by
 Jared A. Grauer
 2012
Preface
 This research was funded in part by the National Aeronautics and Space Admin-
 istration (Nasa) Langley Distinguished Professor grant and by the Army Research
 Laboratory (ARL) Micro Autonomous Systems and Technology (Mast) Collabora-
 tive Technology Alliance (CTA).
 ii
Dedication
 For my family
 iii
Acknowledgments
 Working towards a graduate degree over a signi cant duration of time while dis-
 tributed among a consortium of institutions, one naturally accrues an unwieldy list
 of people to whom they are thankfully indebted. In this section I humbly attempt to
 express my gratitude to those who have most helped me in this journey.
 Firstly, I am extremely fortunate to have had Dr. James E. Hubbard Jr. as a
 mentor, advisor, teacher, and friend. His patience, enthusiasm, and outlook on life
 remain a constant inspiration to me, and his wisdom has strengthened me during
 the times when it was needed most. Furthermore, he has provided a great deal of
 academic encouragement and freedom in my work, as well as a nurturing and team-
 oriented environment in which to pursue it. My best wishes for him and his family.
 I have had a generous committee and a helpful faculty supporting me in this
 research. Dr. Amr Baz and Dr. Alison Flatau, in addition to their years of advising,
 have graciously provided to me coveted real estate in their o ces and laboratories.
 Dr. Darryll Pines has always met me with an enthusiastic, encouraging interest in
 my work. Dr. Eugene Morelli has spent numerous hours coaching me on system
 identi cation methods and discussing experimental results. Dr. Sean Humbert has
 treated me as a member of his own research group, supplying equipment vital to this
 work and  nding always the time to pour over experimental results and speculate
 on their meanings. Dr. Robert Sanner is responsible for awakening my interest in
 dynamics and control, and for helping me strive to do my best work. I very much
 appreciate the extraordinary amount of e ort he has put into my coursework and our
 research discussions. Dr. Ella Atkins has given me my most practical and frequently
 used engineering skills, and warmly remains the audience member at conferences who
 asks the hardest questions.
 Equally as helpful has been the support of lab mates and friends in the Morpheus
 Lab, the Autonomous Vehicle Lab, and the Rotorcraft Center. Morpheus members
 of \Gen 1" (Nelson Guerreiro, Robyn Harmon, Benjamin Nickless, Geo Slipher, and
 Sandra Ugrina) and \Gen 2" (Cornelia Altenbuchner, Alex Brown, and Aimy Wissa)
 have comprised what I can truly only describe as a wolf pack of a research group.
 The Autonomous Vehicle Lab has served as a second home for me, and I am most
 thankful to members Joseph Conroy, Imraan Faruque, and Evan Ulrich for their years
 of friendship and comradery. Among all the many Rotorcraft Center members who
 have enriched my graduate experience, I am most indebted to Brandon Bush, Chen
 Friedman, and Kan Yang for their inspiration across many dimensions, and to J urgen
 Rauleder and Zohreh Ghorbani for my most cherished of friendships.
 Providing perhaps a hidden layer of support has been the unconditional endorse-
 ment of my family. My father, Lawrence Grauer, has instilled in me the patience
 and work ethic needed to complete this degree, while my mother, Barbara Grauer,
 has fed the soul needed to endure this degree. A vital part of this journey has been
 the support of my sisters, Kathryn and Jennifer Grauer, who have provided caring
 ears and a welcomed comic relief. Grandparents, aunts, uncles, and cousins have also
 cheered me along the way, and I feel especially lucky to have such a large and strong
 family behind me.
 iv
Lastly, this section would belie its title without avowing the help of Shelby Dennis
 and her family over the past years. Shelby breathed a purpose and motivation into
 this work that alone I could not respire. Furthermore, the Tasmanian devil that often
 resulted from the occupational hazards of graduate school was gracefully met with a
 contagious warmth that quickly melted troubles into a calmed happiness.
 v
Contents
 List of Tables viii
 List of Figures ix
 Nomenclature x
 1 Introduction 1
 1.1 Background and Motivation . . . . . . . . . . . . . . . . . . . . . . . 1
 1.2 Bio-Inspired Flapping Wing Aircraft . . . . . . . . . . . . . . . . . . 2
 1.3 Flapping Wing Literature Review . . . . . . . . . . . . . . . . . . . . 3
 1.3.1 Vehicle Design . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
 1.3.2 Aerodynamics Modeling . . . . . . . . . . . . . . . . . . . . . 5
 1.3.3 Vehicle Dynamics Modeling . . . . . . . . . . . . . . . . . . . 6
 1.3.4 Feedback Control . . . . . . . . . . . . . . . . . . . . . . . . . 7
 1.4 Scope and Contributions of Current Research . . . . . . . . . . . . . 7
 2 Ornithopter Test Platform Characterizations 10
 2.1 Mathematical Representation of an Aircraft . . . . . . . . . . . . . . 10
 2.2 Ornithopter Aircraft Description . . . . . . . . . . . . . . . . . . . . . 12
 2.3 Measurements from Flight Data . . . . . . . . . . . . . . . . . . . . . 15
 2.3.1 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . . . 15
 2.3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
 2.4 Con guration-Dependent Mass Distribution . . . . . . . . . . . . . . 17
 2.5 Quasi-Hover Aerodynamics . . . . . . . . . . . . . . . . . . . . . . . . 22
 2.6 Implications for Flight Dynamics Modeling . . . . . . . . . . . . . . . 24
 2.7 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
 3 Rigid Multibody Vehicle Dynamics 27
 3.1 Model Con guration . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
 3.2 Kinematic Equations of Motion . . . . . . . . . . . . . . . . . . . . . 30
 3.3 Dynamic Equations of Motion . . . . . . . . . . . . . . . . . . . . . . 31
 3.3.1 Inertial E ects . . . . . . . . . . . . . . . . . . . . . . . . . . 32
 3.3.2 Nonlinear Coupling E ects . . . . . . . . . . . . . . . . . . . . 35
 3.3.3 Gravitational E ects . . . . . . . . . . . . . . . . . . . . . . . 35
 3.4 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
 vi
4 System Identi cation of Aerodynamic Models 37
 4.1 System Identi cation Method . . . . . . . . . . . . . . . . . . . . . . 37
 4.2 Tail Aerodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
 4.2.1 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . . . 38
 4.2.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
 4.3 Wing Aerodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
 4.3.1 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . . . 43
 4.3.2 Measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
 4.3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
 4.4 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
 5 Simulation Results 54
 5.1 Software Simulation Architecture . . . . . . . . . . . . . . . . . . . . 54
 5.2 Determining Trim Solutions . . . . . . . . . . . . . . . . . . . . . . . 56
 5.3 Numerical Linearization about Straight and Level Mean Flight . . . . 59
 5.3.1 Linear Time-Invariant Model . . . . . . . . . . . . . . . . . . 62
 5.3.2 Linear Time-Periodic Model . . . . . . . . . . . . . . . . . . . 64
 5.4 Modeling Implications for Control . . . . . . . . . . . . . . . . . . . . 70
 5.5 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
 6 Concluding Remarks 75
 6.1 Summary of Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
 6.2 Summary of Modeling Assumptions . . . . . . . . . . . . . . . . . . . 76
 6.3 Summary of Original Contributions . . . . . . . . . . . . . . . . . . . 76
 6.4 Recommendations for Future Research . . . . . . . . . . . . . . . . . 77
 A Field Calibration of Inertial Measurement Units 80
 A.1 Theory and Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
 A.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
 B Actuator Dynamics System Identi cation 85
 B.1 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
 B.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
 B.3 Coupling to Vehicle Dynamics . . . . . . . . . . . . . . . . . . . . . . 89
 C Equations of Motion for Single-Body Flight Vehicles 90
 D Linearization of a Conventional Aircraft Model 93
 vii
List of Tables
 2.1 Aircraft parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
 2.2 Avionics measurement speci cations . . . . . . . . . . . . . . . . . . 16
 3.1 Multibody model mass properties . . . . . . . . . . . . . . . . . . . . 28
 4.1 Tail aerodynamic parameters and standard errors . . . . . . . . . . . 42
 4.2 Marker position, rigid body position, and rigid body orientation stan-
 dard errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
 4.3 Wing aerodynamic parameters and standard errors . . . . . . . . . . 52
 5.1 Modal parameters of the decoupled linear time-invariant model . . . . 66
 5.2 Estimated parameters and standard errors for the decoupled linear
 time-periodic model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
 5.3 Modal parameters of the decoupled linear time-periodic model . . . . 71
 A.1 Measurement speci cations for avionics and visual positioning system 81
 A.2 IMU calibration estimates and standard errors . . . . . . . . . . . . . 84
 B.1 Actuator system identi cation measurement speci cations . . . . . . 85
 B.2 DC motor parameter estimates and standard errors . . . . . . . . . . 86
 B.3 Servo motor parameter estimates and standard errors . . . . . . . . . 89
 D.1 Modal parameters of the decoupled F-16 linear model . . . . . . . . . 95
 viii
List of Figures
 1.1 Flight control architecture block diagram . . . . . . . . . . . . . . . . 2
 1.2 Forces acting on an avian wing section, as viewed from a body- xed
 reference frame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
 2.1 An annotated generic aircraft . . . . . . . . . . . . . . . . . . . . . . 11
 2.2 Ornithopter  ight test platform . . . . . . . . . . . . . . . . . . . . . 13
 2.3 Ornithopter planform geometries . . . . . . . . . . . . . . . . . . . . 14
 2.4 Custom avionics board used to record  ight data . . . . . . . . . . . 15
 2.5 Ensemble averages of  ight data with two standard deviation error
 bounds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
 2.6 Power spectrum of  ight data . . . . . . . . . . . . . . . . . . . . . . 19
 2.7 Ornithopter mass distribution contributions . . . . . . . . . . . . . . 20
 2.8 Mass distribution due to wing position variation . . . . . . . . . . . . 21
 2.9 Moment of inertia rates due to  apping between 0 Hz and 10 Hz . . . 22
 2.10 Quasi-hover test experimental setup [17] . . . . . . . . . . . . . . . . 23
 2.11 Representative thrust and lift measurements while  apping at 4.7 Hz 24
 2.12 Ornithopter open-loop  ight dynamics model block diagram . . . . . 25
 3.1 Generic rigid body linkage i on chain j . . . . . . . . . . . . . . . . . 28
 3.2 Multibody ornithopter model schematic . . . . . . . . . . . . . . . . . 29
 4.1 Wind tunnel test experimental setup for tail aerodynamics modeling . 39
 4.2 Flow states from four concatenated wind tunnel tests used for modeling
 the tail aerodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
 4.3 Tail aerodynamic model  ts . . . . . . . . . . . . . . . . . . . . . . . 42
 4.4 Tail aerodynamics model prediction case . . . . . . . . . . . . . . . . 43
 4.5 Flight testing experimental setup for wing aerodynamics modeling . . 44
 4.6 Retro-re ective marker locations on the ornithopter . . . . . . . . . . 45
 4.7 Rigid body ornithopter geometry  t to retro-re ective marker data . 46
 4.8 Longitudinal state variable measurements from  ight data . . . . . . 47
 4.9 Force and moment magnitudes over a wing stroke . . . . . . . . . . . 48
 4.10 Scalar measurement distributions over a wing stroke . . . . . . . . . . 50
 4.11 Parameter estimation of wing aerodynamic models using equation-
 error in the time and frequency domains . . . . . . . . . . . . . . . . 51
 4.12 Wing aerodynamic model prediction case . . . . . . . . . . . . . . . . 52
 ix
5.1 Software simulation block diagram . . . . . . . . . . . . . . . . . . . 55
 5.2 Optimization results for  nding limit cycle trim trajectories . . . . . . 58
 5.3 Variations in trim characteristics with  apping frequency . . . . . . . 58
 5.4 Limit cycle oscillations for  apping frequencies between 4 Hz and 10 Hz 60
 5.5 Trimmed  ight trajectory solution . . . . . . . . . . . . . . . . . . . . 61
 5.6 Linear time-invariant model pole locations . . . . . . . . . . . . . . . 64
 5.7 Linear time-invariant model Eigenvector polar plots . . . . . . . . . . 65
 5.8 Numerically linearized (solid) and curve  tted (dashed) decoupled sys-
 tem matrix element time histories over one wing stroke . . . . . . . . 68
 5.9 Linear time-periodic model pole locations . . . . . . . . . . . . . . . . 70
 5.10 Periodic control gains designed for the LTP system using LQR . . . . 72
 A.1 IMU calibration method block diagram . . . . . . . . . . . . . . . . . 81
 A.2 Uncalibrated 16 bit IMU measurements from a roll-pitch-yaw maneuver 82
 A.3 Calibrated magnetometer and accelerometer signals with centered spheres 83
 A.4 Model  t to the attitude kinematic equation using equation-error . . 84
 B.1 Actuators and instrumented test stand . . . . . . . . . . . . . . . . . 86
 B.2 Measurements used for actuator system identi cation . . . . . . . . . 87
 B.3 Model  ts for actuator dynamics using equation-error and output-error
 in the time and frequency domains . . . . . . . . . . . . . . . . . . . 88
 D.1 Linearized F-16 model pole locations . . . . . . . . . . . . . . . . . . 94
 x
Nomenclature
 Throughout this document, scalar mathematical symbols are represented by lower
 case characters, vectors by lower case bold characters, and matrices by upper case
 bold characters. Di erentiation of a scalar quantity by a column vector results in a
 row vector.
 Roman Symbols
 A, B linear system matrices
 AR aspect ratio
 a(p;v) generalized aerodynamic forces
 a acceleration vector
 b wing span
 C(p;v) dynamic coupling matrix
 Cij, Kij center of mass location and reference frame
 Cx, Cy, Cz force coe cient components
 Cl, Cm, Cn moment coe cient components
 C controllability matrix
  c mean aerodynamic chord
 ex, ey, ez elementary unit vectors
 Ff:g Fourier transform
 f frequency
 ff , Tf  apping frequency and period
 g(p) generalized gravitational forces
 g acceleration vector due to gravity
 h magnetic  eld vector
 I inertia tensor
 I identity matrix
 =f:g, <f:g imaginary and real parts
 Jf force Jacobian matrix
 J attitude Jacobian matrix
 J( ) cost function
 j imaginary number
 L, M , N body- xed moment components
 l, r in-board and out-board position vectors
 M(p) generalized mass matrix
 m scalar mass
 np, nv number of position and velocity states
 xi
p generalized position vector
 p, q, r rotational body- xed velocity components
 Q dynamic pressure
 R rotation matrix
 R2 coe cient of determination
 Re Reynolds number
 r position vector
 S reference area
 S(:) skew operator
 T (p;v) kinetic energy
 t time
 U(p) potential energy
 u perturbation control input vector
 u, v, w translational body- xed velocity components
 V airspeed
 v generalized velocity vector
 W weighting matrix
 X, Y , Z body- xed force components
 x perturbation state vector
 x, y, z Cartesian position components
 y model output vector
 Z, z rotation axis matrix and vector
 z physical state vector
 Greek Symbols
  angle of attack
  sideslip angle
  Hamel coe cient matrix
  increment, perturbation
  w,  lon,  lat wing and tail control inputs
  ,  quaternion scalar and vector parts
  orientation vector
  vector of articulated joint angles
  , v Eigenvalue and Eigenvector
  physical control input vector
  translational velocity vector
  density
  ,  2 covariance matrix and scalar
  generalized force
  (t; t0) state transition matrix
  ,  kinematic Jacobian matrices
  parameter vector
  ,  ,  roll, pitch, and yaw Euler angles
 ! angular velocity vector
 ! radian frequency
 xii
Superscripts
  1 matrix inverse
 T matrix transpose
 _ time derivative
  ensemble average
  trim value
 ^ estimated value
 Subscripts
 ij linkage i on kinematic chain j
 m measured value
 Other
 k:kp p-norm
  Kronecker product
 Acronyms
 CAD computer aided design
 DC direct current
 EE equation-error
 IMU inertial measurement unit
 LTI linear time-invariant
 LTP linear time-periodic
 Mems micro electro-mechanical system
 ODE ordinary di erential equation
 OE output-error
 PSE predicted square error
 Quest quaternion estimator
 Sidpac system identi cation programs for aircraft
 Triad tri-axial attitude determination system
 UAV unmanned air vehicle
 xiii
Chapter 1
 Introduction
 1.1 Background and Motivation
 Unmanned air vehicles (UAVs), or aircraft without on-board pilots, are prolifer-
 ating into many sectors of society and have the potential to signi cantly impact daily
 life. Enabled by advances in electronics miniaturization and composite material man-
 ufacturing, hobbyists, commercial companies, defense contractors, universities, and
 government research laboratories are currently developing small unmanned platforms.
 UAVs will play an important role in the future because in addition to the relatively
 low cost, they can un-man aircraft to perform \dull, dirty, and dangerous" work.
 In the civilian sector, UAVs are currently used as toys, hobby aircraft, and research
  ight platforms. They have been employed in such tasks as airport wildlife control and
 population monitoring, and are expected to contribute to autonomous crop surveying,
 atmospheric weather monitoring, and search and rescue missions. In the military
 sector, these vehicles are employed in intelligence, surveillance, and reconnaissance
 missions. It is envisioned that these vehicles will perform missions such as chemical
 substance detection and autonomous perimeter surveillance, as well as achieve multi-
 mission capabilities including perch and stare, robust long-duration outdoor  ight,
 and agile  ight through cluttered indoor environments.
 To accomplish such mission pro les autonomously, a  ight control system, such
 as that shown Figure 1.1, must be implemented in order to stabilize the vehicle and
 achieve desired performance goals while attenuating the e ects of noise, disturbances,
 and model uncertainty. The aircraft plant model contains a set of equations that
 evolve the vehicle position and velocity state variables as a function of the control
 inputs and disturbances. Sensors on board the aircraft provide measurements, which
 are used in an estimator to predict the aircraft state. The tracking error is supplied
 to the controller, which is designed to shape the inputs to achieve  ight performance
 metrics.
 Fixed-wing aircraft remain the most prevalent type of UAV, as they are in general
 relatively simple aircraft and the dynamic models describing their motions are well
 known and understood. However, as the vehicle size decreases and viscous e ects
 become more pronounced,  xed-wing aircraft su er from decreased lift to drag ratios
 1
Control
 Laws
 Plant
 Dynamics
 Sensor
 Dynamics
 State
 Estimator
  
+
  
reference error input state
 outputstate
 disturbance
 noise
 Figure 1.1: Flight control architecture block diagram
 that degrade  ight performance [1]. Additionally, most  xed-wing aircraft require
 forward speed to generate aerodynamic forces, making it di cult for them to hover
 and  y slowly, which are abilities needed for indoor  ight.
 Miniature rotary-wing vehicles present an alternative to  xed-wing designs. These
 vehicles are highly maneuverable and have the ability to hover, making them suit-
 able for indoor reconnaissance and surveillance missions; however, rotorcraft have a
 number of drawbacks. As with  xed-wing aircraft, viscous e ects reduce the aerody-
 namic e ciency of the vehicle [2]. The added complexity of the rotor and swash plate
 systems increases the vehicle cost and maintenance. Furthermore, these aircraft are
 noisy, open-loop unstable, and potentially hazardous to bystanders and operators,
 making them di cult for human operators to  y without training and closed-loop
 control.
 Although both conventional  xed-wing and rotary-wing vehicles  ll niches in the
 design space of UAVs, a new vehicle design is needed to both  y in a robust manner
 outdoors and in an agile manner in cluttered environments.
 1.2 Bio-Inspired Flapping Wing Aircraft
 Miniature air vehicles that use  apping wings to generate aerodynamic forces and
 moments have advantages over both  xed- and rotary-wing vehicles. For vehicle sizes
 on the order of 0.1 kg and below,  apping wings are required to gain maneuverability
 and aerodynamic e ciency [3, 4, 5]. These vehicles are also safe for humans to operate,
 as they consist only of light weight material and do not have spinning rotors or fuel
 tanks, and they don a level of contextual camou age due to the bird-like appearance.
 It is envisioned that by mimicking the agility, maneuverability, and robustness of
 natural  yers, such as insects, bats, and birds, these vehicles will  ll the demand for
 multi-mission capable UAVs.
 2
Avian  apping  ight is characterized by the up stroke and the down stroke of the
 wings. A cross-section of a seagull wing in forward  ight is illustrated in Figure 1.2,
 as viewed from a body- xed reference frame [6]. Lift and drag forces are generated
 over the wing, perpendicular and parallel to the local  ow direction, respectively, due
 to the forward motion of the vehicle and the plunging motion of the wing. During
 the down stroke, the resultant force is pointed such that a propulsive thrust force
 and a lift force are transmitted to the body, whereas the resultant force is oriented
 such that a lift force and a drag force are transmitted to the body during the up
 stroke. Typically the wing generates most of the lift force near the root of the wing,
 whereas the majority of the thrust is created on the outboard sections of the wing.
 Asymmetries between the up stroke and down stroke are required to produce a net
 lift and thrust force over the  apping cycle. Natural  iers achieve this e ect either
 by  ying with a positive mean angle of attack, modulating the up/down stroke ratio,
 or by changing the wing kinematics through morphing [7, 8].
 1.3 Flapping Wing Literature Review
 1.3.1 Vehicle Design
 The  rst recorded ornithopter vehicle design was DaVinci?s 1490 concept for a
 human-powered ornithopter, having membrane wings and a mechanically e cient
 transmission. In Germany during 1894, Lilienthal created the \Kleiner Schlagf ugelapparat,"
 a human powered  apping device inspired by his previous glider work and experi-
 ments with birds. During the late 19th century designers began experimenting with
 new forms of power generation, including Trouv e?s 1870  ight powered by gunpow-
 der charges and Hargrave?s steam-powered ornithopters in 1890. Beginning in 1942
 with Schmid?s design, ornithopters were powered using internal combustion engines.
 Man-powered ornithopter  ights extending the vehicle endurance time were achieved
 by Lippisch, Maule, Hartmann, Toporov, and Rousseau between 1929 and 1993. The
 most successful manned ornithopter has been that of the University of Toronto, which
 is able to take o under its own power [9].
 Unmanned ornithopter designs began with the rubber band powered models de-
 veloped by Jobert, P enaud, and Villeneuve in the 1870?s. Begining with Spence?s
 1961 ornithopter, designers incorporated smaller power sources and radio electronics
 into the vehicles. Similar ornithopters include the quarter scale model by DeLaurier
 in 1991 and the \Elektro Vogel" series of ornithopter started by R abinger in 1975 [10].
 Ornithopters for hobbyists have also developed, including the Kinkade and Cybird
 ornithopters, in the late 1990?s. Universities have designed ornithopters, including the
 University of Arizona ornithopter [11], \Microbat" at the California Institute of Tech-
 nology [12], the \Phoenix" ornithopter at the Massachusetts Institute of Technology
 [13], and the Morpheus ornithopter at the University of Maryland [14].
 3
plunge
 velocity
 forward velocity
 lift
 drag
 (a) down stroke
 plunge
 velocity
 forward velocity
 lift
 drag
 (b) up stroke
 Figure 1.2: Forces acting on an avian wing section, as viewed from a body- xed
 reference frame
 4
1.3.2 Aerodynamics Modeling
 It is di cult to model the aerodynamics of a  apping wing. One problem is that
 the Reynolds number for these vehicles is typically under 105, which encompasses
 a regime where aerodynamic models are currently innaccurate. Additionally, the
 thin membrane airfoils and light weight structures add complex aero-elastic couplings
 which require additional states or  nite-element models to capture. Furthermore, the
  apping motion of the wings generates complex, unsteady, three-dimensional  ow
  elds that are di cult to model.
 It is generally accepted in the literature that the aerodynamic lift force of  apping
 wings can be crudely, but analytically, represented using the quasi-steady model [15]
 CL = CL  + CL _ _ + CL _h
 _h (1.1)
 where the lift coe cient is a function of the instantaneous angle of attack  , its time
 derivative _ , and the local plunge velocity _h. The drag coe cient is often modeled
 using the drag polar
 CD = CD0 + CD  + CD 2 
2 (1.2)
 in order to capture the parasite drag and drag parabola. Often times these aero-
 dynamic models are posed using two dimensional airfoil sections in a blade element
 model [15, 16, 17, 18, 19, 20], where the local  ow velocities may be more accurately
 represented. Coe cients in these models are derived from  rst principles or measured
 from experiments.
 The most popular analytical model for avian  yers is the model developed by
 DeLaurier, who used a quasi-steady blade element model to capture wake/vortex
 interactions, post-stall phenomenon, and partial leading edge suction. Variants of
 this model are used widely in the literature [21, 22, 23, 16, 17, 24, 25].
 Several attempts to empirically measure aerodynamic models from data have also
 been conducted, mostly as wind tunnel tests where the ornithopter is mounted to a
 load cell measuring lift and drag forces, resulting in a look-up table of aerodynamic
 forces and moments [11, 26, 24]. This method constrains the ornithopter in transla-
 tion and rotation, thus altering the natural response of the vehicle and aerodynamic
  ow environment seen by the wings, and is therefore not representative of free  ight
 aerodynamics.
 Numerical tools have also been developed to examine  apping  ight. While these
 methods require more computation time and can not be used directly for feedback
 control, they provide insight into the fundamental  ight physics, and can be reduced
 to obtain low-order models. Panel methods solving three dimensional, unsteady, in-
 compressible  ows around avian-based  apping models in wind tunnels have shown
 good results [27]. Reynolds-averaged computational  uid dynamics solvers have also
 been used to illuminate the aerodynamic  ow  eld around  apping wings and con-
 tribute the the fundamental understanding of  apping  ight [28, 29].
 The  eld currently lacks an accurate and low-order model for predicting aero-
 dynamic  ows over wings, useful in a control-theoretic framework. This is in part
 due to the variety of ornithopter designs, having membrane wings, thin airfoils, and
 5
structural deformations. Another reason is the lack of experimental data available
 for  apping-wing vehicles in free  ight, due in part to the di culty in obtaining mea-
 surements. Recently developed conformal sensors show potential in measuring lift
 and drag forces in real-time and produce direct measurements of the aerodynamic
  ow  elds [30]. Additionally, non-invasive, vision-based measurements are becom-
 ing popular, which when combined with system identi cation methods can produce
 aerodynamic measurements.
 1.3.3 Vehicle Dynamics Modeling
 There are several types of models that describe the  ight dynamics of  apping-
 wing aircraft. These models are traditionally derived  rst as nonlinear models that
 evolve position and velocity states. For larger aircraft that have signi cant mass
 and inertia in the wings, and/or where the  apping frequencies are near the rigid
 body and control frequencies, multibody models are required. These models are
 more complicated and contain more states, but can capture the inertial e ects of the
  apping wings. Orlowski concluded, based upon a multibody model and empirical
 scaling laws observed in nature, that linear momentum e ects from  apping wings
 are always signi cant, but angular momentum e ects begin to attenuate for  apping
 frequencies above 40 Hz [20]. Rigid body dynamics for multibody ornithopters have
 been derived using a variety of methods [18, 31, 32, 33], and have also been included
 with  nite-element models to capture  uid/structure interactions [25].
 If the inertial e ects of wings are negligible and the wings  ap signi cantly faster
 than the dynamics of interest, then the conventional single-body equations of motion
 for aircraft [34, 35, 36] can be used to model the vehicle dynamics. These models
 are simpler, and the wing motion is either a kinematic state or a control input. This
 approach is common in the literature [37, 38, 39], even for larger vehicles where its
 validity is questionable [22, 11, 40].
 Linear perturbation models are often desired for stability analysis and controls
 synthesis. Linear time-invariant (LTI) models result in equilibrium points with Eigen-
 values and Eigenvectors, for which a large set of classical and modern control tools
 can be applied. Several authors have employed averaging theory to simplify the pe-
 riodic forcing of the wings [39, 41, 20]. Taylor et al. computed the state matrix using
 empirical  nite di erences on insects in a wind tunnel [42, 43]. Faruque and Hum-
 bert used system identi cation to obtain a low-order equivalent system for an insect
 [44, 45]. Kim et al. used a linear discrete-time model of hyper-states to map inputs
 to outputs [46]. Analyses have resulted in a wide variety of modal structures, most of
 which have unstable motions. Several authors have questioned the validity of these
 LTI models [42, 43, 46, 47] and have consequently developed more complex models
 [48].
 The higher  delity models result in trim conditions characterized as limit cycle
 oscillations in the state space, rather than  xed points. Stability has been inferred
 from simulation runs [48], and also proven analytically. For instance, Dietl and Garcia
 applied Floquet analysis using Poincare maps to create an analogous discrete-time
 system. Both their  ndings and a similar study by Bolender showed unstable trim
 6
solutions [22, 49].
 1.3.4 Feedback Control
 One of the goals of modeling the  apping wing vehicle is to be able to control the
 vehicle and complete mission scenarios autonomously. These goals range from robust
 way-point navigation in outdoor gusty environments, to agile  ight through cluttered
 indoor environments, to outdoor precision perching maneuvers.
 Smaller  apping wing vehicles  ap at frequencies too high to perform continuous
 control; instead, a control signal is computed at the beginning of each  ap cycle,
 in a sample and hold fashion. Several authors used variables parameterizing the
 stroke kinematics of each wing, as well as a mechanical bob weight, to achieve a fully
 actuated system capable of position tracking in simulation [19, 23, 37]. Deng et al.
 used averaging theory to obtain an analogous discrete-time, under-actuated system,
 for which stabilized hover was demonstrated in simulation [41].
 Control methods for larger vehicles resemble those of standard aircraft. Krashan-
 itsa et al. integrated a commercial  xed-wing autopilot into an ornithopter and was
 able to demonstrate navigation between two way point coordinates using tuned pro-
 portional control loops [11]. The work demonstrated successful closed-loop control,
 but also the limits of using conventional avionics and control without a model for
  apping-wing aircraft. Recently, Dietl and Garcia explored short comings of using
 classical control and suggested a discrete-time linear quadratic regulator control based
 on a linear time-periodic model [50]. Lee et al. has also recently explored using a bio-
 inspired nonlinear control for the tail to attenuate pitch oscillations and stabilize the
 vehicle [25]. Tedrake has suggested coupling a simple  ight dynamics model with an
 advanced machine learning technique to achieve agile  ight [51].
 1.4 Scope and Contributions of Current Research
 There is currently no  ight dynamics model of ornithopter  ight that is both
 accurate and amenable for control purposes. The goal of this work is to create such
 a model, with insights into design and control, so that the community may progress
 towards realizing an autonomous ornithopter  ight platform.
 In Chapter 2, standard nomenclature is discussed and the experimental  ight
 platform used in this study is introduced. As models and  ight data for this type
 of vehicle are not available in the literature, a series of experimental investigations
 into the characteristics of  ight are presented. The results suggest a nonlinear multi-
 body model is needed to capture the dynamics. Those nonlinear multibody vehicle
 dynamics are derived in Chapter 3. Energy methods are used determine equations
 of motion, which were then cast into a canonical form that is useful for simulation,
 system identi cation, and nonlinear feedback control design.
 Chapter 4 presents the system identi cation work performed to identify parts of
 the ornithopter  ight dynamics model that are not known a priori and are di cult
 to model analytically. A wind tunnel test is discussed where lift and drag models for
 7
the ornithopter tail were identi ed. Afterwards, an experiment is presented using a
 visual tracking system to obtain  ight test data, from which an aerodynamic model
 for the wings was identi ed.
 Chapter 5 presents  ight simulation results obtained using the  ight dynamics
 model developed in Chapters 3 and 4. A method for determining trim solutions is
 presented and is employed to trim the ornithopter model for straight and level mean
  ight. The model is then linearized about the trim trajectory, resulting in a canonical
 time-invariant model as well as a time-periodic model. Feedback control laws are then
 designed for the ornithopter longitudinal dynamics that could be implemented in an
 autopilot.
 The original contributions of this research to the state of the art are the following:
  Flight data suitable for modeling
 This work presents  ight data of an ornithopter, sampled at a bandwidth that
 illuminated the e ects of the  apping wings. Trimmed  ight data showed pitch
 oscillations up to 4.97 rad/s and heave accelerations up to 41.7 m/s2, which re-
 quire nonlinear multibody models to capture. Additionally, harmonic responses
 were seen in the data, indicating either structural modes or nonlinearities. High
 accelerations preclude the use of traditional attitude determination methods and
 require model-based techniques for state estimation.
  Modeling of mass distribution variation due to  apping wings
 The e ects of the  apping motion on the mass distribution of the ornithopter
 were investigated using a CAD model. The center of mass travels 10.3 cm in
 the vertical direction, the moments of inertia can vary up to 53.6%, and the
 inertia rates are signi cant. Multibody models are required to capture these
 variations.
  Derivation of rigid multibody vehicle dynamics
 It was determined that a three-body model was su cient to model the or-
 nithopter. Vehicle dynamics were derived using the Boltzmann-Hamel equa-
 tions, and were cast into a canonical form used for the nonlinear control of
 Euler-Lagrange systems.
  Tail aerodynamics system identi cation
 Wind tunnel tests were conducted on the ornithopter tail, having free stream
 velocities, angles of attack, and angles of sideslip ranging between 4.50 to 6.33
 m/s,  0:88 to  0:35 rad, and  0:57 rad, respectfully. System identi cation
 methods were applied to determine models of the aerodynamic coe cients.
  Flight testing and wing aerodynamics system identi cation
 Flight data, obtained using a visual tracking system, was presented. This data
 and the vehicle dynamics model were used to explore variations in angle of
 attack, Reynolds number, reduced frequency, velocity, and structural deforma-
 tions over the wing and throughout the wing stroke cycle. Additionally, force
 8
and moment sources were identi ed throughout the wing stroke. This data
 was also used to determine aerodynamic models for the wings using system
 identi cation techniques.
  Models for stability and control
 In addition to the full nonlinear model, a linear time-invariant model and a
 linear time-periodic model are presented. The time-invariant model is shown to
 not contain enough accuracy for control design; rather, the time-periodic model
 is needed.
 9
Chapter 2
 Ornithopter Test Platform
 Characterizations
 The form and structure of mathematical models describing conventional aircraft
  ight dynamics are well known [34, 52, 36]; however, it is not obvious whether these
 same models can be used to describe the  ight dynamics of ornithopters. This chapter
 presents experimental investigations into  apping wing  ight using an ornithopter test
 platform. Standard nomenclature for describing aircraft  ight dynamics is introduced,
 followed by a description of the ornithopter test platform. Afterwards,  ight data is
 presented which exhibits nonlinear phenomenon including fast and large amplitude
 motions, as well as limit cycle oscillations. The mass distribution of the ornithopter
 is then explored, and it is shown that the  apping wings induce signi cant changes
 in the center of mass location, inertia tensor, and inertia rates. An experiment is
 summarized where measured lift and thrust forces were modeled using quasi-steady
 aerodynamics. Finally, implications for  ight dynamics modeling are discussed and a
 new model for  apping-wing aircraft is outlined.
 2.1 Mathematical Representation of an Aircraft
 Consider the generic aircraft shown in Figure 2.1. An inertial reference frame
 KI = fexI ; eyI ; ezIg is  xed on the surface of the Earth at CI with unit vectors
 pointing North, East, and down, respectively. A body frame K0 = fex0; ey0; ez0g is
  xed to the aircraft center of mass C0, located at
 rI0;I =
 2
 4
 x
 y
 z
 3
 5 ; (2.1)
 with axes pointing out the nose, right wing, and underside of the aircraft. A wind
 frame KW is collocated with the body frame and is oriented into the oncoming wind.
 The attitude, or orientation, of the vehicle is speci ed by  I0;I . An Euler angle
 10
ex0ey0
 ez0
  0;I
 !0;I
 exIeyI
 ezI
 r0;I
 CI
 C0
 Figure 2.1: An annotated generic aircraft
 parametrization admits the three-element attitude vector
  I0;I =
 2
 4
  
 
 
3
 5 (2.2)
 with corresponding rotation matrix
 R0;I =
 2
 4
 1 0 0
 0 cos sin 
0  sin cos 
3
 5
 2
 4
 cos  0  sin  
0 1 0
 sin  0 cos  
3
 5
 2
 4
 cos sin 0
  sin cos 0
 0 0 1
 3
 5 (2.3)
 where  ,  , and  are the roll, pitch, and yaw angles. Alternatively, a quaternion
 parametrization admits the four-element attitude vector
  I0;I =
  
 
 
 
(2.4)
 and corresponding rotation matrix
 R0;I = ( 2   T )I + 2  T  2 S( ) (2.5)
 where  is the vector part and  is the scalar part of the quaternion, and where the
 skew operator is de ned
 S(x) =
 2
 4
 0  x3 x2
 x3 0  x1
  x2 x1 0
 3
 5 (2.6)
 11
for any three-element vector.
 The aircraft has body- xed translational velocity
  00;I =
 2
 4
 u
 v
 w
 3
 5 (2.7)
 and body- xed rotational velocity
 !00;I =
 2
 4
 p
 q
 r
 3
 5 : (2.8)
 Finally, the local airspeed, angle of attack, and sideslip angle
 V =
 p
 u2 + v2 + w2
  = arctan(w=u)
  = arcsin(v=V ) (2.9)
 are used to construct the rotation matrix
 RW;0 =
 2
 4
 cos cos  sin  sin cos  
 cos sin  cos   sin sin  
 sin 0 cos 
3
 5 (2.10)
 to rotate from the body frame to the wind frame.
 2.2 Ornithopter Aircraft Description
 The commercially available \Slow Hawk" ornithopter, by Sean Kinkade [53], was
 selected for this study due to its payload capacity, durability, and handling character-
 istics. A photograph of this ornithopter is shown in Figure 2.2 and relevant aircraft
 parameters are provided in Table 2.1.
 The fuselage is a 3 mm ply of carbon  ber and has a xed the electronic compo-
 nents. The wings are a rip-stop polyester membrane having 5 mm and 3 mm diameter
 carbon  ber spars along the leading edge and near the trailing edge, respectively. The
 tail is also a membrane, with small carbon  ber spars running from forward to aft.
 Planform geometries are illustrated in Figure 2.3.
 To  y the ornithopter, the aircraft is launched into the air and then piloted using
 standard 72 MHz hobby radio electronics. Power is delivered from a two-cell, 7:4 V,
 1320 mAh lithium-polymer battery to an electronic speed controller, which drives a
 brushless DC motor. This motor is attached to a gear box, four-bar mechanism, and
 leading edge wing spar, and drives the  apping motion of the wing. As the throttle
 is increased, the wings  ap faster, thereby increasing the lift and thrust. Power is
 also delivered to two servo motors which pitch and roll the tail in a serial fashion to
 create pitching and yawing moments on the aircraft, respectively. Flight duration for
 12
Figure 2.2: Ornithopter  ight test platform
 Table 2.1: Aircraft parameters
 Parameter Symbol Value Unit
 total mass m 0.45 kg
 wing span bw 1.22 m
 wing mean aerodynamic chord  cw 0.29 m
 wing area Sw 0.30 m2
 wing aspect ratio ARw 4.40 -
 tail span bt 0.20 m
 tail mean aerodynamic chord  ct 0.20 m
 tail area St 0.04 m2
 tail aspect ratio ARt 1.50 -
 13
-0.1
 +0.0
 +0.1
 -0.6 -0.4 -0.2 +0.0
 z
 [m
 ]
 x [m]
 (a) fuselage
 -0.4
 -0.2
 +0.0
 -0.6 -0.4 -0.2 +0.0 +0.2 +0.4 +0.6
 x
 [m
 ]
 y [m]
 (b) wings
 -0.3
 -0.2
 -0.1
 +0.0
 -0.2 -0.1 +0.0 +0.1 +0.2
 x
 [m
 ]
 y [m]
 (c) tail
 Figure 2.3: Ornithopter planform geometries
 14
(a) obverse (b) reverse
 Figure 2.4: Custom avionics board used to record  ight data
 this aircraft is approximately 15 minutes.
 2.3 Measurements from Flight Data
 Flight data can help guide the  ight dynamics modeling process; however, there
 was no adequate  ight data for ornithopter  ight available in the literature. As such,
 the ornithopter was out tted with an avionics unit to record vehicle pose and inertial
 sensor data during  ight. This section presents the experimental setup,  ight test,
 and data analysis. See reference [54] for a complete discussion.
 2.3.1 Experimental Setup
 Autopilot and stability augmentation systems commonly used on aircraft employ
 sensors including magnetometers, accelerometers, and gyroscopes to provide infor-
 mation on the vehicle motion. Commercially available packages for unmanned air-
 craft typically include sensors that have range and bandwidth too low to investigate
 ornithopter  ight. A custom avionics package was therefore developed to provide
 measurements in  ight from inertial and pose sensors.
 The avionics board, pictured in Figure 2.4 has a 40 cm by 70 cm footprint and
 a 27.9 g mass. The board was designed using CadSoft Eagle, manufactured by an
 external company, and populated using re- ow solder equipment. Power is supplied
 by a dedicated 2-cell, 7:4 V, 200 mAh lithium-polymer battery. Four voltage regu-
 lators deliver power at 3:3 V, 5 V, and 7 V for the on-board processing and sensors.
 Additional connections are included to reprogram the microprocessor and to connect
 to the RS-232 interface
 The MemSense Mag3d inertial measurement unit (IMU) was chosen for its small
 1.78 cm by 1.78 cm footprint and 10 g mass. This IMU provides analog voltage
 outputs,  ltered to 50 Hz, from orthogonal triads of magnetometers, gyroscopes,
 15
Table 2.2: Avionics measurement speci cations
 Measurement Range Resolution Error Unit
 time - 0:03 10 3 - s
 throttle input 0 to 1 1:53 10 3 1:42 10 3 -
 longitudinal input  0:87 to  0:31 8:52 10 6 6:85 10 3 rad
 lateral input  0:70 0:02 10 3 7:44 10 3 rad
 wing angle  3:14 0:10 10 3 20:8 10 3 rad
 longitudinal angle  0:87 to  0:31 9:37 10 6 25:1 10 3 rad
 lateral angle  0:70 0:02 10 3 55:1 10 3 rad
 magnetic  eld  3:17 0:05 10 3 5:83 10 3 -
 rotational velocity  5:24 0:08 10 3 0:14 10 0 rad/s
 linear acceleration  98:1 1:50 10 3 49:6 10 3 m/s2
 and accelerometers, which measure local magnetic  elds h00;I , body- xed rotational
 velocities !00;I , and local accelerations a
 0
 0;I , respectively. Additional sensors were
 installed on the ornithopter to measure temperature, power consumption, wing angle,
 pilot stick commands, and tail surface de ections.
 The board is controlled with a reprogrammable 40 MHz Microchip Pic18f8722
 microprocessor, programmed in the C computer language. During data collection,
 analog sensor outputs are digitized using 16 bit analog to digital converters and are
 recorded on a removable 1 GB memory card. Aggregate delays from data acquisition,
 computation, and memory write times resulted in a maximum sampling rate of 146
 Hz. Measurement speci cations are summarized in Table 2.2.
 The IMU was installed on the ornithopter at the wings level center of mass loca-
 tion, which mitigated DC motor magnetic  eld e ects and mechanical vibrations, as
 well as minimized accelerometer corrections due to rotation [55]. Sensors that mea-
 sured power consumption, wing angle, pilot stick inputs, and tail surface de ections
 were calibrated in the laboratory. Inertial sensors were calibrated at the time and
 place of  ight testing using the method presented in Appendix A.
 2.3.2 Results
 A  ight test was conducted where the ornithopter was manually piloted and
 trimmed for straight and level mean  ight. A 6.4 s segment of data containing 28 com-
 plete wing strokes was chosen for analysis; a con ned  ight volume and the presence
 of wind gusts prohibited longer segments of usable  ight data.
 The  ight data is periodic with the  apping period Tf . As such, ensemble averag-
 ing techniques were applied to compute mean wave forms and variances. For a noisy
 measurement of a Tf -periodic signal ym, the ensemble average and variance are
  ym(t) =
 1
 n
 n 1X
 k=0
 ym(t+ kTf ) (2.11)
 16
 2y(t) =
 1
 n
 n 1X
 k=0
 [ym(t+ kTf )  ym(t)]
 2 (2.12)
 where here n is the number of complete periods contained within the data.
 Ensemble averages are shown in Figure 2.5 with 2 error bounds. The wing
 angle  w, de ned as the roll angle of the right wing relative to the fuselage, oscillates
 with a 0.44 rad amplitude at 4.69 Hz. Magnetometer measurements, normalized
 to the strength of the magnetic  eld of the Earth, exhibit oscillations in the pitch
 axis. The gyroscopes show that while there is some roll and yaw oscillation, the
 dominant motion is in the pitch rate, which obtains a maximum value of +4:97 rad/s
 at 0:40Tf and a minimum value of  4:36 rad/s at 0:90Tf . Periodic forcing of the
 wings and gravitational forces create large heave motions, which are evident by the
 +41:7 m/s2 acceleration at 0:25Tf , and the  21:8 m/s2 acceleration at 0:68Tf . These
 fast and large amplitude motions observed in trimmed  ight require nonlinear models
 to capture. Moreover, the fact that these signals are periodic indicates the presence of
 limit cycle oscillations, which are a nonlinear phenomenon. Increases in error bounds
 at 0:75Tf in the rotational velocity and linear acceleration measurements may indicate
 the presence of unsteady aerodynamic forces in the transition from down stroke to
 up stroke.
 The frequency content of the  ight data was also examined using a high accuracy
 implementation of the chirp-z algorithm [56, 35] to compute the Fourier transform
 Ffy(t)g =
 Z t
 0
 y( )e j! d (2.13)
 using an arbitrarily speci ed set of frequencies !. Frequency transforms were com-
 puted using a 0.10 Hz spacing up to 50 Hz, and are shown in Figure 2.6. The sensors
 show power at harmonics of the  apping frequency, which indicate structural modes
 in the airframe and/or nonlinear e ects in the  ight dynamics.
 2.4 Con guration-Dependent Mass Distribution
 Flight dynamics models typically approximate aircraft as a single rigid body and
 neglect changes in mass distribution due to rotating machinery, control surface de ec-
 tions, and fuel burn. The ornithopter wings however have appreciable mass compared
 with the rest of the vehicle and  ap at frequencies near the expected rigid body modes
 of the system. In this section it is shown that the  apping wings create signi cant
 changes in the center of mass location, inertia tensor, and inertia rate, which cannot
 be captured by the standard aircraft  ight dynamics model.
 The ornithopter consists of three sections of parts: the fuselage, wings, and tail.
 For this ornithopter, the moving parts constitute a signi cant portion of the vehi-
 cle mass, as illustrated in Figure 2.7. Both a component build up method and a
 computer-aided design (CAD) program were employed to model the vehicle mass dis-
 tribution. The  rst method consisted of measuring part masses and geometries from
 which inertia tensors were then computed and transformed to the vehicle center of
 17
-1.0
 +0.0
 +1.0
 0.00 0.25 0.50 0.75 1.00
 ? w
 [ra
 d]
 t/Tf
 (a) wing angle
 -1.0
 +0.0
 +1.0
 h x
 [-]
 -1.0
 +0.0
 +1.0
 h y
 [-]
 -1.0
 +0.0
 +1.0
 0.00 0.25 0.50 0.75 1.00
 h z
 [-]
 t/Tf
 (b) magnetometer
 -9.0
 +0.0
 +9.0
 p
 [ra
 d/
 s]
 -9.0
 +0.0
 +9.0
 q
 [ra
 d/
 s]
 -9.0
 +0.0
 +9.0
 0.00 0.25 0.50 0.75 1.00
 r
 [ra
 d/
 s]
 t/Tf
 (c) gyroscope
 -50
 +0
 +50
 a
 x
 [m
 /s2
 ]
 -50
 +0
 +50
 a
 y
 [m
 /s2
 ]
 -50
 +0
 +50
 0.00 0.25 0.50 0.75 1.00
 a
 z
 [m
 /s2
 ]
 t/Tf
 (d) accelerometer
 Figure 2.5: Ensemble averages of  ight data with two standard deviation error bounds
 18
0.0
 0.5
 1.0
 0 10 20 30 40 50
 ? w
 [ra
 d]
 f [Hz]
 (a) wing angle
 +0.0
 +0.3
 +0.6
 h x
 [-]
 +0.0
 +0.3
 +0.6
 h y
 [-]
 +0.0
 +0.3
 +0.6
 0 10 20 30 40 50
 h z
 [-]
 f [Hz]
 (b) magnetometer
 +0.0
 +4.5
 +9.0
 p
 [ra
 d/
 s]
 +0.0
 +4.5
 +9.0
 q
 [ra
 d/
 s]
 +0.0
 +4.5
 +9.0
 0 10 20 30 40 50
 r
 [ra
 d/
 s]
 f [Hz]
 (c) gyroscope
 +0
 +25
 +50
 a
 x
 [m
 /s2
 ]
 +0
 +25
 +50
 a
 y
 [m
 /s2
 ]
 +0
 +25
 +50
 0 10 20 30 40 50
 a
 z
 [m
 /s2
 ]
 f [Hz]
 (d) accelerometer
 Figure 2.6: Power spectrum of  ight data
 19
fuselage
 75:6%
 wings
 19:5%
 tail 4.90%
 Figure 2.7: Ornithopter mass distribution contributions
 mass. In the second method, software was used to compute the inertia tensor from
 mass measurements and three-view photographs. Both methods produced consistent
 results, however the CAD results were used because they are generally considered to
 be more accurate.
 The e ect of wing position on the gross vehicle center of mass location, referenced
 to the fuselage center of mass, is shown in Figure 2.8(a). The symmetric motion of
 the wings creates a 10.3 cm (0.36  cw) variation in the vertical center of mass position.
 The vertical bias is due to the wings-level position being located above the fuselage
 center of mass. This change in center of mass location introduces complexity into the
 conventional aircraft model, as the body and stability axes are located at the center
 of mass.
 The e ect of wing position on the gross vehicle inertia tensor, referenced to the
 gross vehicle center of mass, is shown in Figure 2.8(b). The products of inertia
 are small and do not exhibit large changes. The roll inertia shows 32.2% variations
 due to the wing hinge location above the gross center of mass location. The pitch
 inertia is biased due the hinge o set and the moving center of mass, and experiences
 a 53.6% increase in value as the wings rotate. The yaw inertia is centered about
 the wings-level position and experiences a symmetric 52.6% change in inertia as the
 wings move. These variations introduce a position-dependent inertia matrix into the
 standard aircraft model.
 The time rate of change of the inertia tensor also contributes inertial forces to the
 equations of motion. Inertia rates resulting from a sinusoidal approximation to the
 wing trajectory
  w = 0:56 sin(2 ff t) 0:13 (2.14)
 are shown in Figure 2.9 as a function of the wing angle, where the  apping frequency
 20
-80.0
 -40.0
 +0.0
 +40.0
 +80.0
 -2 -1 +0 +1 +2
 ce
 n
 te
 r
 o
 fm
 as
 s
 po
 sit
 io
 n
 [m
 m
 ]
 ?w [rad]
 x
 y
 z
 (a) center of mass
 +0.0
 +6.0
 +12.0
 +18.0
 -2 -1 +0 +1 +2
 in
 er
 tia
 [g
 ?
 m
 2
 ]
 ?w [rad]
 Ixx Iyy Izz Ixy Ixz Iyz
 (b) inertia tensor
 Figure 2.8: Mass distribution due to wing position variation
 21
-50
 -25
 +0
 +25
 +50
 ? I x
 x
 [g
 ?
 m
 2
 /s]
 -200
 -100
 +0
 +100
 +200
 ? I y
 y
 [g
 ?
 m
 2
 /s]
 -200
 -100
 +0
 +100
 +200
 -0.8 -0.6 -0.4 -0.2 +0.0 +0.2 +0.4 +0.6 +0.8
 ? I z
 z
 [g
 ?
 m
 2
 /s]
 ?w [rad]
 0
 2
 4
 6
 8
 10
 Figure 2.9: Moment of inertia rates due to  apping between 0 Hz and 10 Hz
 ff = 1=Tf is varied between 0 and 10 Hz. The pitch inertia rate is smaller in amplitude
 because the pitch inertia experiences smaller changes, and is much larger above  0:5
 rad due to a greater slope. The pitch rate inertia and yaw rate inertia are null at
 approximately 0.1 rad and 0 rad, where the inertias reach a minimum and a maximum,
 respectively. The plots grow in magnitude for negative wing angles because the
 asymmetric nature of (2.14) puts more travel and faster velocity in that regime. The
 inertial rate e ects are on the same order as magnitude as the inertial e ects for the
 ornithopter and require multibody models to capture.
 2.5 Quasi-Hover Aerodynamics
 It is di cult to model aerodynamic forces and moments on an arbitrary  apping
 wing aircraft from  rst principles. To experimentally study the ornithopter aero-
 dynamics, force and wing shape measurements were collected while  apping in a
 quasi-hover state.
 The experimental setup is pictured in Figure 2.10 [17]. Lift and thrust measure-
 ments were recorded from a six-axis strain gauge transducer, and wing angle mea-
 surements form a potentiometer, at 1000 Hz using LabVIEW. The ornithopter wing
 was also  tted with one hundred 3 mm diameter retro-re ective spherical markers.
 A visual tracking system, consisting of six cameras and a software suite, was used to
 track the spatial position of each marker at a rate of 350 Hz, providing a discretized
 22
Figure 2.10: Quasi-hover test experimental setup [17]
 measurement of the wing shape. The ornithopter was  apped at a number of frequen-
 cies between 3 and 7 Hz while the aerodynamic forces, wing angle, and wing shape
 were recorded.
 A representative segment of lift and thrust data is shown in Figure 2.11, where
 the ornithopter was  apping at 4.7 Hz. Thrust forces range between 0.26 and 6.86 N,
 and oscillate at twice the  apping frequency, with the larger peak occurring at the
 transition between down stroke to up stroke and the smaller peak occurring at the
 transition between up stroke and down stroke. The lift force oscillates at the  apping
 frequency, and is zero mean with a 15.6 N amplitude.
 This experiment resulted in two models that used wing shape data as input to
 validate the aerodynamic forces. The  rst was a computational  uid dynamics model,
 which used a Reynolds-averaged, Navier-Stokes solver called Umturns [57] and a
 deforming grid algorithm to capture the wing deformations [28]. Good agreement
 was found at the low  apping frequencies; however at higher frequencies inaccuracies
 were found in the thrust model. The second analysis was a blade element model using
 aerodynamics from  rst principles [16, 17]. It was found that the forward section of
 the wing creates the lift while the aft section creates the thrust. A quasi-steady
 aerodynamics model could predict the lift force but not the thrust force.
 23
+0.0
 +2.0
 +4.0
 +6.0
 +8.0
 th
 ru
 st
 [N
 ]
 -16
 -8
 +0
 +8
 +16
 0.0 0.2 0.4 0.6 0.8 1.0
 lif
 t[
 N
 ]
 t [s]
 Figure 2.11: Representative thrust and lift measurements while  apping at 4.7 Hz
 2.6 Implications for Flight Dynamics Modeling
 The  ight data showed the presence of fast and large amplitude motions, as well as
 limit cycle oscillations. Additionally, the harmonic responses measured by the sensor
 are representative of nonlinear systems or structural vibrations. The large heave
 accelerations preclude the used of conventional attitude determination algorithms
 [58] and require model-based state estimators to estimate orientation.
 The  apping wings generate a moving center of mass, a variable inertia tensor,
 signi cant inertia rates, and a moving reference frame. These di culties can be
 modeled in a straight forward manner using multibody system techniques.
 The quasi-hover experiment provided a priori information on the thrust and lift
 waveforms that can be expected in free  ight. Modeling results showed that quasi-
 steady models may be su cient to model the lift force.
 A block diagram for the proposed ornithopter  ight dynamics model is shown in
 Figure 2.12. The pilot control inputs  are the  apping rate and the tail angles, which
 impact the aerodynamics. These forces and torques, along with gravitational e ects,
 force the rigid body dynamics, which evolve the generalized position p and velocity v
 states. The rigid body dynamics are nonlinear multibody equations of motion which
 will be developed in Chapter 3, along with the gravitational e ects. The aerodynamic
 e ects will be determined using system identi cation techniques on wind tunnel and
 free  ight data in Chapter 4. Additional actuator models for this ornithopter have
 been determined using system identi cation and are provided in Appendix B.
 24
Gr
 avitationa
 l
 E ect
 s
 Win
 g
 Aer
 o
 dynamic
 s
 T
 ai
 l
 Aer
 o
 d
 y
 namic
 s
 Nonlinea
 r
 Multi
 b
 o
 d
 y
 Kinematic
 s
 &
 Dynamic
 s
 a tai
 l(
 p
 ;v
 )
 a w
 in
 g
 (p
 ;v
 )
 g
 (p
 )
  
z
 =
  
p v
  
Figur
 e
 2.12
 :
 Ornithopte
 r
 op
 en-l
 oo
 p
  ig
 ht
 dynamic
 s
 m
 ode
 l
 bl
 oc
 k
 diagra
 m
 25
2.7 Chapter Summary
 This chapter presented an investigation into the  ight characteristics of an exper-
 imental  apping-wing ornithopter test platform. A  ight test was conducted where
 the ornithopter was  tted with a custom avionics package. Results showed that the
 ornithopter experiences harmonic responses, large amplitude pitch excursions, fast
 pitch rates, and large heave accelerations. Geometry modeling using a CAD software
 showed that as the ornithopter experiences large changes in the vertical position of
 the center of mass, inertia tensor, and inertia rates as the wings  ap. An experiment
 in a quasi-hover condition showed small deformations of the wing spars, complex
  ow structures, and the e ectiveness of a quasi-steady blade element model of the
 aerodynamics. Findings in this chapter suggest a nonlinear multibody model of the
 ornithopter vehicle dynamics and a quasi-steady model of the aerodynamics.
 26
Chapter 3
 Rigid Multibody Vehicle Dynamics
 Chapter 2 concluded with a recommendation to model ornithopter  ight dynam-
 ics using nonlinear multibody techniques. In this chapter, the nonlinear equations
 of motion are derived for a multibody ornithopter. The model con guration is  rst
 presented, where the ornithopter is simpli ed as a system of three rigid bodies. The
 kinematic equations are then derived for this system, followed by the dynamic equa-
 tions. Energy methods were employed, which resulted in a minimum state model,
 cast into a canonical form and having Lyapunov functions derived from scalar energy
 functions. At the conclusion of this chapter, only the aerodynamics model remains
 to be determined.
 3.1 Model Con guration
 The ornithopter is assumed to be a set of three rigid bodies, with one body
 allocated for the fuselage and one for each wing, all connected by revolute joints.
 Although susceptible to bending and torsion, the fuselage is assumed rigid for control
 modeling purposes. Similarly, the wings are  exible membranes and the aft section
 de ects to create thrust; however, the wings are assumed to be a rigid body for control
 modeling, and the aft section contains only 12% of the wing mass, which is 2.3% of
 the vehicle mass. The mass and inertia of the tail was absorbed into the fuselage, as
 the tail comprises only 4.9% of the vehicle mass, is actuated slowly (see Appendix B),
 and adds considerable complexity to the multibody model. Similarly, the four-bar
 linkage that  aps the wings was not modeled because although it does integrate the
  apping kinematics into the equations of motion, it does not contribute a signi cant
 source of momentum and adds considerable complexity. Momentum contributions
 from spinning gears were also neglected.
 A generic rigid body, linkage i along chain j, is shown in Figure 3.1. Each rigid
 body has at its center of mass Cij a mass mij and inertia tensor Iij. Position vectors
 lij and rij describe the locations of in-board and out-board joint locations, which
 rotate about vectors zij and z(i+1)j through angles  ij and  (i+1)j. Additionally, an
 arbitrary point P can be speci ed by the vector rP;ij.
 The ornithopter multibody con guration is shown in Figure 3.2 with mass and
 27
Cij
 zij z(i+1)j
  ij  (i+1)j
 lijrij
 rP;ij
 Figure 3.1: Generic rigid body linkage i on chain j
 Table 3.1: Multibody model mass properties
 Index Mass (kg) Inertia Tensor (kg m2)  10 6
 ij m Ixx Iyy Izz Ixy Ixz Iyz
 00 0.3412 +112:57 +3799:3 +3739:4  51:394 +5:3746  1:1234
 11 0.0414 +915:06 +422:48 +1337:5  166:07  0:0001  0:0000
 12 0.0414 +915:06 +422:48 +1337:5  166:07  0:0001  0:0000
 inertia properties listed in Table 3.1. The generalized position vector is
 p =
 2
 4
 rI0;I
  I0;I
  w
 3
 5 (3.1)
 where the  rst term is the inertial position of fuselage, the second term is a quaternion
 describing the orientation of the fuselage, and where the last term
  w =  11 =    12 (3.2)
 is the wing angle. The generalized velocity vector is
 v =
 2
 4
  00;I
 !00;I
 _ w
 3
 5 (3.3)
 where the  rst two terms are the body- xed translational and rotational velocities of
 the fuselage, and where the last term is the wing rate. This choice of state variables
 is the same as for conventional aircraft, except that an additional degree of freedom
 is provided for the wings.
 28
CI
 C0
 C11
 C12
 r0;I
 r01
 l11
 l12
 z11
  11
 z12
  12
 Figure 3.2: Multibody ornithopter model schematic
 29
3.2 Kinematic Equations of Motion
 The kinematic equations describe the evolution of the position state as a function
 of the velocity state. A generic point P on link i of chain j has the position vector
 rIP;I = r
 I
 0;I + R
 I;0r00j +
 i 1X
 k=1
 RI;kj(rkjkj  l
 kj
 kj) + R
 I;ij(rijP;ij  l
 ij
 ij) (3.4)
 expressed in the inertial reference frame. The rotational velocity of a linkage in the
 ornithopter system can be written in its body frame as
 !ijij;I = R
 ij;0!00;I +
 iX
 k=1
 zkjkj _ kj
 = Rij;0!00;I + Z
 ij
 ij
 _ (3.5)
 which forms a convenient way of writing the equations of motion. The translational
 velocity of a point P on linkage i of chain j is
  0P;I =  
0
 0;I +S(!
 0
 0;I)r
 0
 0j +
 i 1X
 k=1
 S(!0kj;I)R
 0;kj(rkjkj  l
 kj
 kj) +S(!
 0
 ij;I)R
 0;ij(rijp;ij l
 ij
 ij) (3.6)
 expressed in the fuselage frame of reference.
 The kinematic di erential equation relating the position and velocity states is
 _p =  v (3.7)
 where
  =
 2
 4
 RI;0 0 0
 0 JI;0 0
 0 0 I
 3
 5 (3.8)
 is a block-diagonal Jacobian matrix with an inverse relation  . The  rst entry is
 a rotation matrix, de ned using either (2.3) or (2.5), that rotates the translational
 velocities of the fuselage frame into the inertial coordinate frame. The second term
 is a rotational Jacobian matrix, de ned as either
 JI;0 =
 2
 4
 1 sin tan  cos tan  
0 cos  sin 
0 sin cos  cos cos  
3
 5 (3.9)
 or
 JI;0 =
 1
 2
  
 I + S( )
   T
  
(3.10)
 and relates body- xed velocities to attitude rates. The last entry is an identity matrix,
 re ecting the kinematically  at relationship of the joint angles and their rates.
 30
3.3 Dynamic Equations of Motion
 The dynamic equations describe the evolution of the velocity as a function of the
 states and controls. These equations are derived here using energy methods, which
 produce  rst-order state space equations with a minimal number of state variables.
 As aircraft dynamics are not typically derived using energy methods, Appendix C is
 included to provide detail and familiarity.
 The kinetic and potential energy, which may serve as future Lyapunov control
 functions [59], are
 T (p;v) =
 1
 2
 niX
 i=0
 njX
 j=0
 mij( 0ij;I)
 T ( 0ij;I) + (!
 ij
 ij;I)
 T Iijij(!
 ij
 ij;I) (3.11)
 U(p) =  
niX
 i=0
 njX
 j=0
 mij(rIij;I)
 TgI (3.12)
 where nj are the total number of kinematic chains, ni are the number of linkages on
 chain j, and g is the local acceleration due to gravity.
 The Boltzmann-Hamel equations [60, 61], a generalization of the Lagrange equa-
 tions which account for di erent reference frames,
 d
 dt
  
@T
 @v
  T
 +
  
nvX
 k=1
 @T
 @vk
  k
 !
 v   T
  
@T
 @p
  T
 =  (3.13)
 may then be applied, where potential energy sources are neglected and the Hamel
 coe cient matrices are de ned as
  k =  T k (3.14)
 with individual matrix elements
 f kgij =
 @ ki
 @pj
  
@ kj
 @pi
 : (3.15)
 For a given multibody system and choice of reference frames, the relationship between
 the velocities and quasi-velocities of the system is  xed, and so the Hamel coe cient
 matrices (3.14) are constant, skew symmetric matrices.
 The equations (3.13) can be manipulated into a canonical form used in the non-
 linear control of Euler-Lagrange systems [62, 63, 64, 59]. The kinetic energy of the
 system can be written
 T (p;v) =
 1
 2
 vTM(p)v (3.16)
 which is a quadratic form in the velocity. Inserting (3.16) into (3.13), di erentiating,
 31
and regrouping the terms, yields the reformulation
 M_v +
  
_M +
 nvX
 k=1
 @T
 @vk
  k  
1
 2
  T (I vT )
  
@M
 @p
  T
 !
 v =  (3.17)
 where  is a Kronecker product operator [65].
 Grouping terms in (3.17) and adding forces imparted on the vehicle by the external
 environment yields the canonical form
 M(p) _v + C(p;v)v + g(p) + a(p;v) =  (3.18)
 where M(p) is a generalized mass matrix, C(p;v) contains nonlinear coupling forces
 arising from centripetal and Coriolis accelerations, g(p) describes gravitational e ects,
 and a(p;v) describes aerodynamic e ects. With this choice of state variables, the
 forcing vector  corresponds to forces on the fuselage center of mass, torques on
 the fuselage center of mass, and torques on the wings. In the following sections these
 matrices and vectors are derived and discussed, with the exception of the aerodynamic
 terms, which are presented in Chapter 4. Properties of the canonical form are further
 discussed in Kelly et al. [63] and Lewis et al. [62].
 3.3.1 Inertial E ects
 The mass matrix M(p) describes the generalized mass and inertia of the system
 as a function of the vehicle pose. Contributions to this matrix can be determined for
 each rigid body by using algebraic manipulations to rearrange the kinetic co-energy
 (3.11) into the quadratic form (3.16), for which the mass matrix is then obvious.
 The mass matrix for the fuselage body is found by rearranging the energy and
 has the same form as that found in conventional aircraft equations of motion [34, 36].
 The mass matrix is
 M(p) =
 2
 6
 4
 m0I 0 0
 0 I00 0
 0 0 0
 3
 7
 5 (3.19)
 which is a diagonal matrix with the mass and inertia tensor along the block diagonal.
 The mass matrix for the  rst linkage (i = 1, e.g. the wings) is given in (3.20), and
 should a second linkage (i = 2, e.g. the four-bar mechanism, tail surface, etc.) need
 to be modeled, the mass matrix is given in (3.21). Each mass matrix has along its
 block diagonal the scalar mass, the inertia relative to the fuselage center of mass, and
 the inertia of the linkage. O -diagonal terms represent inertial couplings within the
 system. The total mass matrix of the system is a summation of the mass matrices of
 each rigid body.
 32
M
 (p
 )
 =
 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4
 m
 1j
 I
 m
 1j
 S
 (R
 0;
 1j
 l1
 j 1j
 )
  
m
 1j
 S
 (r
 0 0j
 )
 m
 1j
 S
 (R
 0;
 1j
 l1
 j 1j
 )Z
 0 1j
 :::
 m
 1j
  (
 r0 0
 j)
 T
 (r
 0 0j
 )
  
2(
 l1
 j 1j
 )T
 R
 1j
 ;0
 r0 0
 j
 +
 (l
 1j 1j
 )T
 (l
 1j 1j
 ) 
I
  
m
 1j
 (r
 0 0j
 )(
 r0 0
 j)
 T
 +
 m
 1j
 R
 0;
 1j
 l1
 j 1j
 (r
 0 0j
 )T
 +
 m
 1j
 r0 0
 j(
 l1
 j 1j
 )T
 R
 1j
 ;0
  
m
 1j
 R
 0;
 1j
 (l
 1j 1j
 )(
 l1
 j 1j
 )T
 R
 1j
 ;0
 +
 R
 0;
 1j
 I1
 j 1j
 R
 1j
 ;0
 m
 1j
  (
 l1
 j 1j
 )T
 (l
 1j 1j
 )
  
(l
 1j 1j
 )T
 R
 1j
 ;0
 r0 0
 j 
Z
 0 1j
 +
 m
 1j
 R
 0;
 1j
 l1
 j 1j
 (r
 0 0j
 )T
 Z
 0 1j
  
m
 1j
 R
 0;
 1j
 l1
 j 1j
 (l
 1j 1j
 )T
 R
 1j
 ;0
 Z
 0 1j
 +
 R
 0;
 1j
 I1
 j 1j
 Z
 0 1j
 :::
 :::
 m
 1j
 (Z
 0 1j
 )T
  (
 l1
 j 1j
 )T
 (l
 1j 1j
 ) 
Z
 0 1j
  
m
 1j
 (Z
 0 1j
 )T
 R
 0;
 1j
 l1
 j 1j
 (l
 1j 1j
 )T
 R
 1j
 ;0
 Z
 0 1j
 +
 (Z
 1j 1j
 )T
 I1
 j 1j
 Z
 1j 1j
 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5
 (3.20
 )
 33
M
 (p
 )
 =
 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4
 m
 2j
 I
  
m
 2j
 S
 (r
 0 0j
 )
  
m
 2j
 S
 (R
 0;
 1j
 (r
 1j 1j
  
l1
 j 1j
 ))
 +
 m
 2j
 S
 (R
 0;
 2j
 l2
 j 2j
 )
  
m
 2j
 S
 (R
 0;
 1j
 (r
 1j 1j
  
l1
 j 1j
 ))
 Z
 0 1j
 +
 m
 2j
 S
 (R
 0;
 2j
 l2
 j 2j
 )Z
 0 2j
 :::
 m
 2j
  (
 r0 0
 j)
 T
 (r
 0 0j
 ) 
I
  
m
 2j
 (r
 0 0j
 )(
 r0 0
 j)
 T
 +
 m
 2j
  (
 r0 0
 j)
 T
 R
 0;
 1j
 (r
 1j 1j
  
l1
 j 1j
 ) 
I
  
m
 2j
 r0 0
 j(
 r1
 j
 1j
  
l1
 j 1j
 )T
 R
 1j
 ;0
 +
 m
 2j
 r0 0
 j(
 l2
 j 2j
 )T
 R
 2j
 ;0
  
m
 2j
  (
 r0 0
 j)
 T
 R
 0;
 2j
 l2
 j 2j
  
I
 +
 m
 2j
  (
 r1
 j
 1j
  
l1
 j 1j
 )T
 R
 1j
 ;0
 r0 0
 j 
I
  
m
 2j
 R
 0;
 1j
 (r
 1j 1j
  
l1
 j 1j
 )(
 r0 0
 j)
 T
 +
 m
 2j
 (r
 1j 1j
  
l1
 j 1j
 )T
 (r
 1j 1j
  
l1
 j 1j
 )I
  
m
 2j
 R
 0;
 1j
 (r
 1j 1j
  
l1
 j 1j
 )(
 r1
 j
 1j
  
l1
 j 1j
 )T
 R
 1j
 ;0
 +
 m
 2j
 R
 0;
 1j
 (r
 1j 1j
  
l1
 j 1j
 )(
 l2
 j 2j
 )T
 R
 2j
 ;0
  
m
 2j
  (
 r1
 j
 1j
  
l1
 j 1j
 )T
 R
 1j
 ;2
 j l
 2j 2j
  
I
 +
 m
 2j
 R
 0;
 2j
 (l
 2j 2j
 )(
 r0 0
 j)
 T
  
m
 2j
  (
 l2
 j 2j
 )T
 R
 2j
 ;0
 r0 0
 j 
I
 +
 m
 2j
 R
 0;
 2j
 l2
 j 2j
 (r
 1j 1j
  
l1
 j 1j
 )T
 R
 1j
 ;0
  
m
 2j
  (
 l2
 j 2j
 )T
 R
 2j
 ;1
 j (
 r1
 j
 1j
  
l1
 j 1j
 ) 
I
 +
 m
 2j
  (
 l2
 j 2j
 )T
 l2
 j 2j
  
I
  
m
 2j
 R
 0;
 2j
 (l
 2j 2j
 )(
 l2
 j 2j
 )T
 R
 2j
 ;0
 +
 R
 0;
 2j
 I2
 j 2j
 R
 2j
 ;0
 m
 2j
  (
 r1
 j
 1j
  
l1
 j 1j
 )T
 R
 1j
 ;0
 r0 r
 0j
  
Z
 0 1j
  
m
 2j
 R
 0;
 1j
 (r
 1j 1j
  
l1
 j 1j
 )(
 r0 0
 j)
 T
 Z
 0 1j
 +
 m
 2j
  (
 r1
 j
 1j
  
l1
 j 1j
 )T
 (r
 1j 1j
  
l1
 j 1j
 ) 
Z
 0 1j
 +
 m
 2j
 R
 0;
 1j
 (r
 1j 1j
  
l1
 j 1j
 )(
 r1
 j
 1j
  
l1
 j 1j
 )T
 R
 1j
 ;0
 Z
 0 1j
 +
 m
 2j
 R
 0;
 1j
 (r
 1j 1j
  
l1
 j 1j
 )(
 l2
 j 2j
 )T
 R
 2j
 ;0
 Z
 0 1j
  
m
 2j
  (
 r1
 j
 1j
  
l1
 j 1j
 )T
 R
 1j
 ;2
 j l
 2j 2j
  
Z
 0 1j
 ;
 +
 m
 2j
 R
 0;
 2j
 (l
 2j 2j
 )(
 r0 0
 j)
 T
 Z
 0 2j
  
m
 2j
  (
 l2
 j 2j
 )T
 R
 2j
 ;0
 r0 0
 j 
Z
 0 2j
 +
 m
 2j
 R
 0;
 2j
 l2
 j 2j
 (r
 1j 1j
  
l1
 j 1j
 )T
 R
 1j
 ;0
 Z
 0 2j
  
m
 2j
  (
 l2
 j 2j
 )T
 R
 2j
 ;1
 j (
 r1
 j
 1j
  
l1
 j 1j
 ) 
Z
 0 2j
 +
 m
 2j
  (
 l2
 j 2j
 )T
 (l
 2j 2j
 ) 
Z
 0 2j
  
m
 2j
 R
 0;
 2j
 (l
 2j 2j
 )(
 l2
 j 2j
 )T
 R
 2j
 ;0
 Z
 0 2j
 R
 0;
 2j
 I2
 j 2j
 Z
 2j 2j
 :::
 :::
 m
 2j
 (Z
 0 1j
 )T
 (r
 1j 1j
  
l1
 j 1j
 )T
 (r
 1j 1j
  
l1
 j 1j
 )Z
 0 1j
  
m
 2j
 (Z
 0 1j
 )T
 R
 0;
 1j
 (r
 1j 1j
  
l1
 j 1j
 )(
 r1
 j
 1j
  
l1
 j 1j
 )T
 R
 1j
 ;0
 Z
 0 1j
 +
 m
 2j
 (Z
 0 2j
 )T
 R
 0;
 1j
 (r
 1j 1j
  
l1
 j 1j
 )(
 l2
 j 2j
 )T
 R
 2j
 ;0
 Z
 0 1j
  
m
 2j
 (Z
 0 2j
 )T
  (
 r1
 j
 1j
  
l1
 j 1j
 )T
 R
 1j
 ;2
 j l
 2j 2j
  
Z
 0 1j
 +
 m
 2j
 (Z
 0 1j
 )T
 R
 0;
 2j
 l2
 j 2j
 (r
 1j 1j
  
l1
 j 1j
 )T
 R
 1j
 ;0
 Z
 0 2j
  
m
 2j
 (Z
 0 1j
 )T
  (
 l2
 j 2j
 )T
 R
 2j
 ;1
 j (
 r1
 j
 1j
  
l1
 j 1j
 ) 
Z
 0 2j
 +
 m
 2j
 (Z
 0 2j
 )T
  (
 l2
 j 2j
 )T
 (l
 2j 2j
 ) 
Z
 0 2j
  
m
 2j
 (Z
 0 2j
 )T
 R
 0;
 2j
 (l
 2j 2j
 )(
 l2
 j 2j
 )T
 R
 2j
 ;0
 Z
 0 2j
 +
 (Z
 2j 2j
 )T
 I2
 j 2j
 Z
 2j 2j
 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5
 (3.21
 )
 34
3.3.2 Nonlinear Coupling E ects
 The nonlinear coupling matrix is a non-unique, square matrix that contains forces
 arising from Coriolis and centripetal accelerations, and is most easily derived from
 the mass matrix. It is apparent that one de nition of the coupling matrix is
 C(p;v) = _M +
 nvX
 k=1
 @T
 @vk
  k  
1
 2
  T (I vT )
  
@M
 @p
  T
 (3.22)
 from comparing (3.17) and (3.18). A property exploited in passivity-based control is
 that there exists a representation of the coupling matrix that satis es
 _M(p;v) 2C(p;v) (3.23)
 is a skew-symmetric matrix [64, 66]. The representation (3.22) does not satisfy this
 relation; rather, by extension of the discussion by Lewis et al. [62], the representation
 is
 C(p;v) =
 1
 2
 _M +
  
nvX
 k=1
 @T
 @vk
  k
 !
 +
 1
 2
  
@M
 @p
  
(v I)  
1
 2
  T (I vT )
  
@M
 @p
  T
 (3.24)
 so that
 _M 2C =  2
  
nvX
 k=1
 @T
 @vk
  k
 !
  
 
@M
 @p
  
(v  I) +  T (I vT )
  
@M
 @p
  T
 (3.25)
 which satis es the property (3.23) because the  rst term is a summation of skew
 Hamel coe cient matrices and the last two terms represent the skew decomposition
 of a matrix.
 Evaluating the coupling matrix (3.24) for a given ornithopter state is not a trivial
 task. The mass matrix derivatives, both respect to time and position, were derived.
 Some simpli cations of the equations arise, as the Kronecker products have many null
 entries and the mass matrix is only a function of the joint position states.
 3.3.3 Gravitational E ects
 Gravitation e ects depend on the orientation of the ornithopter, and are derived
 from the potential energy as
 g(p) =  T
  
@U
 @p
  T
 (3.26)
 however, this formulation leads to complicated expressions involving tensors. It is
 simpler in this case to include gravitational e ects as external forces applied at the
 35
centers of mass
 g(p) =
 niX
 i=0
 njX
 j=0
 Jf;ij
 2
 4
 0
 0
 mijgI
 3
 5 (3.27)
 where the force Jacobian matrices are
 Jf;0 =
 2
 4
 R0;I
 0
 0
 3
 5 (3.28)
 for the fuselage i = 0,
 Jf;1j =
 2
 4
 R0;I
 S(r00j  R
 0;1jl1j1j)
 (Z1j1j)
 TS( l1j1j)R
 1j;I
 3
 5 (3.29)
 for the  rst linkage i = 1 on a chain j, and
 Jf;2j =
 2
 4
 R0;I
 S(r00j + R
 0;1j(r1j1j  l
 1j
 1j) R
 0;2jl2j2j)R
 0;I
 (Z1j)TS((r
 1j
 1j  l
 1j
 1j) R
 1j;2jl2j2j)R
 1j;I + (Z2j2j)
 TS( l2j2j)R
 2j;I
 3
 5 (3.30)
 for the second linkage i = 2 on a chain j.
 3.4 Chapter Summary
 In this chapter, the nonlinear multibody equations of motion for an ornithopter
 were derived. A rigid body is allocated for the fuselage and each wing. Kinematic and
 dynamic equations of motion are presented. The Boltzmann-Hamel equations were
 used, which result in Lyapunov control functions and a minimal state space model
 that was cast into a canonical form. Gravitational e ects were also derived. At this
 point, the only unmodeled forces on the ornithopter are due to aerodynamics, which
 will be presented in Chapter 4.
 36
Chapter 4
 System Identi cation of
 Aerodynamic Models
 In Chapter 2 it was determined that a nonlinear multibody model was needed to
 describe the  ight dynamics of the ornithopter. In Chapter 3 the rigid body equa-
 tions of motion, as well as the gravitational e ects, were derived from  rst principles.
 The aerodynamics for an arbitrary  apping wing vehicle are di cult to model ana-
 lytically, and are determined in this work using system identi cation methods with
 experimental data. This chapter begins with a presentation of the system identi -
 cation methods employed herein. Afterwards, tail aerodynamic models are identi ed
 from wind tunnel data, and wing aerodynamic models are identi ed from  ight test
 data. This chapter concludes the nonlinear multibody model of the ornithopter  ight
 dynamics.
 4.1 System Identi cation Method
 System identi cation is the process of determining, using observations of inputs
 and outputs, a mathematical model that behaves similarly to the physical system
 under speci ed conditions. Applied to aircraft, system identi cation can be used to
 create  ight simulators, validate numerical codes, design control laws, assess handling
 qualities, perform model reduction, and gain physical insight into the system. In this
 work, system identi cation was employed to develop aerodynamic models for the
 ornithopter that could not be determined from  rst principles modeling. Several
 texts provide a thorough description of system identi cation techniques [67, 35, 68].
 The system identi cation process begins with an experiment where inputs and
 outputs are recorded while the system is su ciently excited. After any necessary
 calibration, signal processing, data compatibility and data collinearity analyses, the
 next step is model structure determination, which is the process of  nding an adequate
 functional representation relating inputs to outputs. While there are many techniques
 for this step, step-wise regression is most heavily relied on in this work [69, 35]. In this
 method, a pool of candidate regressors are selected and are sequentially included in or
 excluded from the model per engineering intuition and several statistical metrics. The
 37
process is iterative and interactive, and can be applied to linear or nonlinear problems,
 in the time or frequency domains, using any method of parameter estimation.
 Once the form of the model is known, the next task is to estimate the unknown
 parameters in the model. For the system
 _z = f(t; z; ) (4.1)
 with states z, inputs  , and measurements ym, the maximum likelihood cost function
 J( ) =
 1
 2
 nX
 i=1
 eTW 1e (4.2)
 where  is the vector of unknown parameters, e is the residual vector between
 the measurements and the outputs, and W is the inverse noise covariance matrix.
 Two simpli cations of the maximum likelihood estimator are the equation-error and
 output-error methods [70, 35]. In equation-error, process noise is assumed and the
 measurements are the state derivatives. This results in a linear estimation problem
 with an analytical solution. In output-error, measurement noise is assumed and the
 measurements are the states. This results in a nonlinear, iterative solution that re-
 quires initial guesses, but is considered more accurate than equation-error. Once
 parameter estimates are obtained, a new data set is used to test the predictive capa-
 bility of the identi ed model.
 All signal processing and system identi cation algorithms are implemented in this
 work using a Matlab toolbox called System IDenti cation Programs for AirCraft
 (Sidpac), developed at Nasa Langley Research Center and documented in Reference
 [35].
 4.2 Tail Aerodynamics
 This section describes the system identi cation of the ornithopter tail aerodynamic
 model using wind tunnel data. The tail aerodynamics are assumed to be steady due to
 limited pilot bandwidth [71] and actuator bandwidths (provided in Appendix B), and
 independent of the wing aerodynamics. These simplifying assumptions are common
 in the literature due to the complex nature of the interactions [25, 22, 18].
 4.2.1 Experimental Setup
 The ornithopter was tested in a free-jet wind tunnel, having a 0.56 m by 0.56 m
 test section. The wings were removed and the ornithopter was  tted on a custom
 test stand in the center of the test section, as shown in Figure 4.1. Free stream
 velocity measurements were obtained using a U-tube manometer and a barometer.
 The test stand included a load cell that produced six-axis force and torque data. Tail
 positions were measured using the internal potentiometers in the servo motors. The
 tail orientation was controlled using potentiometers and Pic18F452 microprocessors,
 as well as a signal generator.
 38
Figure 4.1: Wind tunnel test experimental setup for tail aerodynamics modeling
 39
+4.0
 +5.0
 +6.0
 +7.0
 +8.0
 V
 [m
 /s]
 -1.0
 -0.8
 -0.6
 -0.4
 -0.2
 ?
 [ra
 d]
 -0.8
 -0.4
 +0.0
 +0.4
 +0.8
 0 100 200 300 400 500 600 700 800
 ?
 [ra
 d]
 t [s]
 Figure 4.2: Flow states from four concatenated wind tunnel tests used for modeling
 the tail aerodynamics
 Four tests were selected for system identi cation of aerodynamic force and moment
 coe cients. Each test occurred at a di erent free stream velocity, and a potentiome-
 ter was used to step the angle of attack on the wings while the sideslip angle was
 cycled at 0.06 Hz using a signal generator. Aerodynamic quantities for these tests are
 concatenated and shown in Figure 4.2. Flow velocities ranged between 4.50 m/s and
 6.33 m/s; the lower bound was the slowest velocity with an acceptable signal to noise
 ratio, and the upper bound began to exhibit structural deformations and increases
 in actuator current consumption. The angle of attack range, -0.88 to -0.35 rad, and
 angle of sideslip range,  0:57 rad, comprise the full travel of the tail orientation. A
 frequency domain analysis showed that no signi cant unsteady e ects were present
 in the data due to the cycling of the sideslip angle, and the associated reduced fre-
 quency ranges between 0.006 and 0.009, which is a regime where unsteady e ects are
 not signi cant [15].
 40
4.2.2 Results
 Data were calibrated,  ltered to 1 Hz using a low-pass smoother, and then con-
 catenated together. Aerodynamic force and moment coe cients were computed as
 Cx = X=QSt
 Cy = Y=QSt
 Cz = Z=QSt
 Cl = L=QStbt
 Cm = M=QStct
 Cn = N=QStbt (4.3)
 where
 Q =
 1
 2
  V 2 (4.4)
 is the dynamic pressure, X, Y , and Z are the body axis forces, L, M , and N are the
 body axis moments, and tail geometry parameters are listed in Table 2.1.
 Model structure determination was performed using step-wise regression in the
 time domain using equation-error. The model structure
 Cx = Cx0 + Cx 2 
2
 Cy = Cy    
Cz = Cz0 + Cz  + Cz  
Cl = Cl  + Cl    
Cm = Cm0 + Cm  + Cm  + Cm 2 
2 + Cm V  V
 Cn = Cn  (4.5)
 was found to  t the data well. The model is nonlinear and expresses the aerodynamic
 force and moment coe cients in terms of the measured aerodynamic  ow variables.
 The lateral models are zero mean while the longitudinal models contain biases to
 account for pro le drag, as well as lift and pitching moment at zero angle of attack.
 The model for drag is the traditional drag polar. The yawing moment is dependent
 only on the sideslip angle. The remaining models contain a coupling between the angle
 of attack and sideslip angle, due to the serial con guration of the tail surface. The
 measured yawing moment is much larger than the rolling moment, which indicates
 the lateral control is primarily in yaw.
 Time domain equation-error was used to estimate model parameters, as the model
 equations precluded an output-error analysis, and testing at steady free stream ve-
 locities prohibited an accurate frequency domain analysis. Measured and modeled
 coe cients are shown in Figure 4.3, and estimated parameters with standard errors
 are given in Table 4.1. The Cz model had a lower  t value due to the residual seen in
 the third and fourth data sets, and the Cz and Cm 2 coe cients had slightly larger
 errors, but otherwise the models produced high  t values and parameter errors below
 10%, indicating good model  ts.
 41
-0.60
 -0.40
 -0.20
 C
 x
 -0.10
 +0.00
 +0.10
 C
 y
 +0.00
 +0.50
 +1.00
 C
 z
 -0.01
 +0.00
 +0.01
 C
 l
 +0.00
 +1.00
 +2.00
 C
 m
 -0.40
 +0.00
 +0.40
 0 100 200 300 400 500 600 700 800
 C
 n
 t [s]
 data model
 Figure 4.3: Tail aerodynamic model  ts
 Table 4.1: Tail aerodynamic parameters and standard errors
 Parameter Estimate Fit
   ^  ( ^) R2
 Cx0  0:3181 0:0058 0.88Cx 2  0:2310 0:0123
 Cy  +0:1153 0:0035 0.88
 Cz0 +0:3346 0:0243
 0.59Cz  0:2729 0:0364
 Cz +0:0884 0:0066
 Cl  0:0054 0:0005 0.93Cl   0:0161 0:0007
 Cm0  0:3486 0:0272
 0.96
 Cm  3:3182 0:1131
 Cm +0:0975 0:0067
 Cm 2  0:4184 0:0847
 Cm V +0:3053 0:0195
 Cn   0:5094 0:0141 0.97
 42
-0.60
 -0.40
 -0.20
 C
 x
 -0.10
 +0.00
 +0.10
 C
 y
 +0.00
 +0.50
 +1.00
 C
 z
 -0.01
 +0.00
 +0.01
 C
 l
 +0.00
 +1.00
 +2.00
 C
 m
 -0.40
 +0.00
 +0.40
 0 25 50 75 100 125 150 175 200
 C
 n
 t [s]
 data model
 Figure 4.4: Tail aerodynamics model prediction case
 A similar test, not used for identi cation and having a free stream velocity of
 4.86 m/s, was used to validate the identi ed model. Measured and predicted force
 and moment coe cients are shown in Figure 4.4. The rolling and pitching moment
 coe cients have small residuals at low angles of attack, but overall the models predict
 well and all have coe cients of determination above 0.81, indicating a good model.
 4.3 Wing Aerodynamics
 This section describes the system identi cation of the ornithopter wing aerody-
 namics model from  ight data. The experimental setup is  rst presented, followed by
 a discussion of the data reduction and kinematic  ight results, and lastly a presenta-
 tion of the system identi cation results.
 4.3.1 Experimental Setup
 It is di cult to obtain  ight data for system identi cation of an ornithopter.
 Avionics add weight and inertia to the ornithopter, and must sample faster than
 with conventional aircraft to measure the e ects of  apping. Additionally, the large
 heave accelerations prohibit an attitude estimation using conventional methods with
 magnetometers and accelerometers, and air data probes incur inaccuracies due to
 the fast pitching and heaving motions. Flexibility in the aircraft also adds errors to
 the responses at harmonics of the  apping frequency. To mitigate these problems, a
 visual tracking system was used to measure the spatial location of markers placed on
 43
Figure 4.5: Flight testing experimental setup for wing aerodynamics modeling
 the ornithopter. At the cost of restricting the capture volume for data collection, this
 is a non-invasive and highly accurate method that additionally averages out e ects
 due to  exibility.
 Flight tests consisted of  ying the ornithopter down a corridor, as shown in Fig-
 ure 4.5, while collecting data. This method was su cient for modeling straight and
 level  ight while minimizing atmospheric disturbances and facilitating the use of the
 visual tracking system. Eight cameras were placed at the corners of a capture volume
 10 m long by 6 m wide and 5 m tall. These cameras recorded at 500 Hz the spatial
 position of retro-re ective markers placed on the ornithopter fuselage, wings, and tail,
 as shown in Figure 4.6. Marker position measurements were  ltered to 250 Hz using
 a low-pass smoother.
 To compute the state vector, rigid bodies were  t to the marker data at each data
 frame. A simplex search algorithm was used  nd the center of mass position and
 44
(a) fuselage
 (b) wings (c) tail
 Figure 4.6: Retro-re ective marker locations on the ornithopter
 Table 4.2: Marker position, rigid body position, and rigid body orientation standard
 errors
 Rigid Body Marker Position Orientation
 Cij  (rIk;I) [m]  (r^
 I
 ij;I) [m]  ( ^
 I
 ij;I) [rad]
 fuselage 0.0009 0.0009 0.0056
 right wing 0.0017 0.0011 0.0041
 left wing 0.0018 0.0013 0.0037
 tail 0.0008 0.0006 0.0036
 orientation that minimized the cost function
 J( ) =
 1
 2
 eTWe (4.6)
 where  is a vector containing the rigid body position and orientation quaternion,
 W is a unity weighting matrix, and the error
 e = rIk;I  (r
 I
 ij;I + R
 I;ijrijk;ij) (4.7)
 is the di erence between the measured marker locations and those computed from
 the rigid body center of mass position, a rotation matrix derived from its orientation,
 and known positions of the markers relative to the center of mass. An example of the
  t is shown in Figure 4.7. Errors for the marker measurements, as well as the rigid
 body positions and orientations, are given in Table 4.2. From these estimates, the
 position states are assembled and subsequently smoothly di erentiated [35] to yield
 the velocity states. Previous studies have shown these measurements contain lower
 noise levels than measurements obtained using Mems gyroscopes and accelerometers
 [72].
 45
fuselage
 right wing
 left wing
 tail
 +3.6
 +3.8
 +4.0
 +4.2
 +4.4
 exI [m]
 -0.6
 -0.3
 +0.0
 +0.3
 +0.6
 eyI [m]
 -1.1
 -1.0
 -0.9
 -0.8
 ezI [m]
 Figure 4.7: Rigid body ornithopter geometry  t to retro-re ective marker data
 4.3.2 Measurements
 Longitudinal state variables for a representative segment of  ight data, where the
 ornithopter was  apping at 5.91 Hz, is shown in Figure 4.8. The periodic wing forcing
 induces a heave velocity between  1:1 m/s and +0:90 m/s, which results in a 0.08
 m oscillation about the nominal altitude. Smoothly di erentiating the heave velocity
 shows that at the fuselage center of mass, the ornithopter experiences translational
 accelerations between  4:0 and +4:9 times the that of gravity. The wing beats also
 induce a pitch rate between  1:4 and +1:6 rad/s, resulting in a 0.02 rad oscillation
 about the mean pitch angle. These observations are consistent with the inertial
 measurements presented in Chapter 2. The longitudinal tail de ection, although
 commanded to a trim condition, oscillates due to the vehicle and actuator dynamics.
 The down stroke phase of the  apping cycle has a range of  0:67 to +0:39 rad, which
 occurs between 0:25Tf and 0:78Tf , yielding a 1.18 down stroke to up stroke ratio.
 Using the identi ed model structure (3.18), the  ight data can be partitioned into
 contributions arising from inertial, dynamic coupling, gravitational, and aerodynamic
 sources. Magnitudes of force and moment contributions are shown in this manner in
 Figure 4.9 as polar plots, where the angle represents the location within the  ap-
 ping cycle and the radius denotes the magnitude. Figure 4.9(a) shows that forces
 arising from inertial and wing aerodynamic contributions dominate just before and
 after transitions between up and down stroke. Figure 4.9(b) shows that the nonlin-
 ear coupling and wing aerodynamic terms contribute most heavily to the moments
 throughout the wing stroke, but diminish at the transitions between up and down
 stroke, which is when the wing changes direction.
 In a similar fashion, kinematic quantities can be displayed through the wing stroke,
 as in Figure 4.10, where the angular coordinate denotes the position in the  apping
 cycle, and where the radial component denotes the position along the wing span.
 Figure 4.10(a) shows estimates of the main wing spar bending deformation, found by
 interpolating three marker locations along the wing span. Maximum tip de ections
 of  0:02 m, which are relatively small, are found 0:13Tf after the transitions between
 up and down strokes. Figure 4.10(b) shows the velocity magnitude, which has a
 46
x [m
 ]
 u [m
 /s]
 z [m
 ]
 w [m
 /s]
 ? [ra
 d] q
 [ra
 d/
 s]
 ? w [ra
 d]
 ? ? w [ra
 d/
 s]
 ? l
 o
 n
 [ra
 d]
 t [s]
 ? ? l
 o
 n
 [ra
 d/
 s]
 t [s]
 +0.00
 +3.00
 +6.00
 +8.0
 +9.0
 +10.0
 -0.95
 -0.90
 -0.85
 -1.50
 +0.00
 +1.50
 -0.05
 +0.00
 +0.05
 -2.00
 +0.00
 +2.00
 -0.80
 +0.00
 +0.80
 -20.0
 +0.0
 +20.0
 -0.20
 -0.10
 +0.00
 0.00 0.10 0.20 0.30 0.40 0.50
 -10.0
 +0.0
 +10.0
 0.00 0.10 0.20 0.30 0.40 0.50
 (a) time domain
 x [m
 ]
 u [m
 /s]
 z [m
 ]
 w [m
 /s]
 ? [ra
 d] q
 [ra
 d/
 s]
 ? w [ra
 d]
 ? ? w [ra
 d/
 s]
 ? l
 o
 n
 [ra
 d]
 f [Hz]
 ? ? l
 o
 n
 [ra
 d/
 s]
 f [Hz]
 +0.00
 +0.02
 +0.04
 +0.00
 +0.02
 +0.04
 +0.00
 +0.01
 +0.02
 +0.00
 +0.20
 +0.40
 +0.00
 +0.01
 +0.02
 +0.00
 +0.20
 +0.40
 +0.00
 +0.10
 +0.20
 +0.00
 +3.00
 +6.00
 +0.00
 +0.01
 +0.02
 0 10 20 30 40 50
 +0.00
 +1.00
 +2.00
 0 10 20 30 40 50
 (b) frequency domain
 Figure 4.8: Longitudinal state variable measurements from  ight data
 47
30
 20
 10
 0
 10
 20
 30
 20 10 0 10 20 30 40
 inertia
 coupling
 gravity
 tail aero
 wing aero
 (a) body force magnitude (N)
 4
 2
 0
 2
 4
 4 2 0 2 4
 inertia
 coupling
 gravity
 tail aero
 wing aero
 (b) body moment magnitude (Nm)
 Figure 4.9: Force and moment magnitudes over a wing stroke
 48
minimum value of 8.9 m/s due to forward velocity, and a maximum value at the
 wing tips of 11.0 m/s due to the  apping motion. The angle of attack distribution
 is shown in Figure 4.10(c), which varies between  0:94 rad and +0:86 rad. The
 Reynolds number is shown in Figure 4.10(d) and varies between 19,000 and 232,000
 for this  ight, which encompasses the transition region where steady lift to drag ratios
 decrease [1]. Reduced frequency measurements in Figure 4.10(e) show that this  ow
 regime may be classi ed as unsteady [15].
 4.3.3 Results
 System identi cation was performed to determine a model for the longitudinal
 force, heave force, and pitching moment generated by the wings. Rearranging (3.18),
 the generalized aerodynamic force due to each wing is
 awing(p;v) =   [M(p) _v + C(p;v)v + g(p) + atail(p;v)] (4.8)
 where atail(p;v) is the aerodynamic contributions of the tail identi ed in the previous
 section. The  rst, third, and  fth elements of this vector represent the longitudinal
 force, heave force, and pitching moment, respectively, and are normalized in a similar
 fashion as (4.3). Other forces and moments were not identi ed because the excitation
 during  ight testing was only in the longitudinal axis.
 Model structure determination was performed in the time domain with equation-
 error and step-wise regression. The resulting model structure was
 Cx = Cx0 + Cx w  w + Cxqq + Cx _ w
 _ w + Cx 2w  
2
 w
 Cz = Cz0 + Czqq + Cz _ w
 _ w
 Cm = Cm0 + Cm _ w
 _ w (4.9)
 which is a nonlinear expansion in terms of the state variables. Model structures were
 also found using orthogonal regressors [73] and a traditional quasi-steady model [33],
 but yielded less accurate results. The thrust force includes a polar to capture the
 e ects of drag, and also depends on the pitch rate and wing motion. The lift is a
 function of the pitch rate and wing motion, similar to a quasi-steady model. The
 pitching moment variation is due only to the wing motion.
 As many state variables are periodic with harmonics of the  apping frequency, a
 collinearity analysis was performed on the model regressors. One method examined
 the correlation matrix between the regressors, which had maximum value of 0.57;
 because this value is less than 0.9, guidelines indicate that collinearity is not an
 issue [35, 74]. Additionally, an Eigensystem analysis was performed, which resulted
 in condition numbers under 3.70; because these numbers are much less than 100,
 guidelines also indicate here that collinearity is not an issue [35, 75].
 Equation-error was used to  nd the parameters in (4.9) that  t the aerodynamic
 coe cients found using (4.8). Time domain and frequency domain  ts are shown in
 Figure 4.11 and parameter estimates are provided in Table 4.3. The time domain
 49
(0/4)Tf
 (1/4)Tf
 (2/4)Tf
 (3/4)Tf
 -0.02
 -0.01
 +0.00
 +0.01
 (a) wing deformations (m)
 (0/4)Tf
 (1/4)Tf
 (2/4)Tf
 (3/4)Tf
 10
 11
 12
 13
 14
 (b) velocity magnitude (m/s)
 (0/4)Tf
 (1/4)Tf
 (2/4)Tf
 (3/4)Tf
 -0.8
 -0.4
 +0.0
 +0.4
 +0.8
 (c) angle of attack (rad)
 (0/4)Tf
 (1/4)Tf
 (2/4)Tf
 (3/4)Tf
 4e+04
 8e+04
 1e+05
 2e+05
 2e+05
 (d) Reynolds number
 (0/4)Tf
 (1/4)Tf
 (2/4)Tf
 (3/4)Tf
 +0.1
 +0.2
 +0.3
 +0.4
 +0.5
 +0.6
 +0.7
 (e) reduced frequency
 Figure 4.10: Scalar measurement distributions over a wing stroke
 50
-0.5
 +0.0
 +0.5
 +1.0
 C
 x
 0.00
 0.02
 0.04
 0.06
 -4.0
 -2.0
 +0.0
 +2.0
 +4.0
 C
 z
 0.00
 0.20
 0.40
 0.60
 -4.0
 -2.0
 +0.0
 +2.0
 +4.0
 0.0 0.1 0.2 0.3 0.4 0.5
 C
 m
 t [s]
 0.00
 0.10
 0.20
 0.30
 0.40
 0.0 4.0 8.0 12.0 16.0
 f [Hz]
 data
 model
 Figure 4.11: Parameter estimation of wing aerodynamic models using equation-error
 in the time and frequency domains
 results for the longitudinal force coe cient had larger error bounds and a lower  t
 value than the frequency domain results. While frequency domain modeling is often
 superior to time domain, this may also be due to low excitation, a short data record,
 or a large number of regressors. Heave force results were similar to the longitudinal
 force results. Moment coe cient results were very similar and accurate.
 The identi ed models were used to predict the aerodynamic forces and moments
 in a second set of  ight data. These results are shown in Figure 4.12. The heave
 force and the pitching moment were well predicted by both models. The longitudinal
 force model, however, appears to be over-parametrized; longer data records with more
 airspeed excitation, e.g. from  ap variation or during push-over/pull-up maneuvers,
 are required to improve this model. The frequency domain model had higher  t values
 and was selected as the wing aerodynamic model.
 4.4 Chapter Summary
 Previously in Chapter 3, the rigid body dynamics of the multibody ornithopter sys-
 tem were derived analytically. This chapter applied system identi cation techniques
 to extract models for the aerodynamic forces on the tail and wings from experimental
 data.
 A wind tunnel test was conducted to determine the aerodynamic forces resulting
 from de ections of the tail  ap. Simple models based on angle of attack and side slip
 51
Table 4.3: Wing aerodynamic parameters and standard errors
 Parameter Time Estimate Frequency Estimate Time Fit Frequency Fit
   ^  ( ^)  ^  ( ^) R2 R2
 Cx0 +0:2262 0:1068 +0:0634 0:0145
 Cx w  0:1010 0:0896  0:0568 0:0071
 Cxq  0:1913 0:1168  0:3647 0:0068 0.66 0.80
 Cx _ w +0:0127 0:0030 +0:0205 0:0002
 Cx 2w  1:4455 0:2236  1:1843 0:0165
 Cz0  0:9222 0:2052  0:9655 0:0628
 Czq  1:8822 0:4566  2:4701 0:0594 0.73 0.86
 Cz _ w  0:0627 0:0197  0:0619 0:0018
 Cm0 +0:0955 0:0172 +0:0942 0:0154 0.96 0.96Cm _ w +0:0862 0:0038 +0:0836 0:0005
 -0.5
 +0.0
 +0.5
 +1.0
 C
 x
 -6.0
 -3.0
 +0.0
 +3.0
 +6.0
 C
 z
 -2.0
 +0.0
 +2.0
 +4.0
 0.00 0.10 0.20 0.30 0.40 0.50 0.60
 C
 m
 t [s]
 data
 time domain
 frequency domain
 Figure 4.12: Wing aerodynamic model prediction case
 52
angle were able to describe the aerodynamic forces generated.
 A  ight test was conducted to observe the ornithopter under free  ight condi-
 tions. A visual tracking system was employed to measure the spatial position of
 retro-re ective markers placed on the ornithopter, from which the state vector was
 computed. Several kinematic variables were presented for the ornithopter over the
 course of a wing stroke from  ight data, providing insight into the  ight physics.
 System identi cation yielded a low-order model for the aerodynamic forces generated
 by the  apping wings.
 This chapter and the previous chapter fully describe the nonlinear, multibody
 model of the ornithopter  ight dynamics, determined using a combination of analyt-
 ical and experimental methods.
 53
Chapter 5
 Simulation Results
 In the previous chapters, the equations of motion were determined for the or-
 nithopter. In this chapter, those equations are programmed into a nonlinear  ight
 simulator for the purposes of simulation and model simpli cation. An overview of
 the simulation environment is  rst presented. Afterwards, a method is discussed for
 trimming the model for straight and level mean  ight. Numerical methods are then
 employed to numerically linearize the model about straight and level mean  ight, re-
 sulting in both a conventional time-invariant model as well as a time-periodic model.
 5.1 Software Simulation Architecture
 The ornithopter  ight dynamics (3.7) and (3.18) were written as the  rst-order
 state space model
 _z = f(z; ) (5.1)
 and were coded in both Fortran 95 and Matlab, according to the block diagram
 shown in Figure 5.1. The Fortran code ran very quickly; however, at the expense of
 computation time, the Matlab environment provided a wealth of internal functions
 and tools such as plotting and optimization routines. Simulation of a single wing  ap
 using Matlab used approximately 0.94 seconds of computation time on a dual-core
 laptop, which was su cient for this work.
 A main script is  rst run, which initializes an ornithopter con guration  le, sim-
 ulation parameters, control laws, and an initial state vector. This information is
 given to an ordinary di erential equation (ODE) solver, and returns the state vector
 time history. Data is then plotted, animated, analyzed, and saved. The ODE solver
 numerically integrates the equations of motion from an initial state and returns the
 state vector time history. In the Fortran implementation of the code, an explicit,
 fourth-order, Runge-Kutta method was employed, where inversion of the mass matrix
 was accomplished using a Cholesky decomposition and back-substitution. Solution
 time histories converged using a 10 4 s  xed time step interval. The Matlab imple-
 mentation of the code utilized the internal solver ode15s, which was chosen because
 it is an adaptive time step method for accurate solutions, implicit solver for saving
 computation time, and valid for sti ODE problems.
 54
Main
 Script
 ODE
 Solver
 Equations of Motion
 _z = f(z; )
 con guratio
 n
 kinematic
 s
 mas
 s
 matri
 x
 M
 (p
 )
 couplin
 g
 matri
 x
 C
 (p
 ;v
 )
 gr
 avi
 ty
 g
 (p
 )
 aer
 odynamic
 s
 a(
 p
 ;v
 )
 co
 ntro
 l
  
initial conditions
 z(t0)
 state time history
 z
 state vector
 z(ti)
 state derivative
 _z(ti)
 Figure 5.1: Software simulation block diagram
 55
The ODE solver calls the equations of motion  le and computes the current state
 derivative. Within the routine, several subroutines are called. The  rst is a con gu-
 ration  le, speci c to each ornithopter, that loads mass and inertia properties, linkage
 lengths, and actuator constants into the workspace. Next is a kinematics subroutine,
 which computes the necessary rotation matrices and kinematic variables. Lastly, sub-
 routines are called that assemble the mass matrix, coupling matrix, gravitational and
 aerodynamic forces, and control inputs.
 Several checks were performed to validate the  ight simulation. First, the hand-
 derived equations of motion were checked against output of the software algebra
 package Mathematica. Second, simulations from random initial conditions without
 external forces resulted in state trajectories having constant kinetic energy. Third,
 simulations from random initial conditions with only conservative forces resulted in
 state trajectories having constant total energy, again within numerical accuracy.
 5.2 Determining Trim Solutions
 Conventional aircraft have trim solutions characterized by  xed points in the state
 space that correspond to a constant state z and control   . Nonlinear, iterative
 solvers are traditionally used to numerically solve
 _z = f(z ;  ) = 0 (5.2)
 for the trim states and controls. Flapping-wing vehicles, on the other hand, have
 Tf -periodic forcings which manifest as limit cycles, or trim conditions that are closed
 trajectories in the state space. Stable limit cycles are relatively easy to  nd because
 state trajectories beginning within the basin of attraction converge on the limit cycle
 as time progresses. While Lee et al. [23] developed a simulation that exhibits a stable
 limit cycle behavior, the majority of work in the literature consists of models that
 yield instabilities in forward  ight [22, 49, 42, 76, 20].
 Simulations of the the ornithopter dynamics model, during straight and level mean
  ight, exhibited divergent  ight trajectories, suggesting an unstable limit cycle. As
 limit cycles are periodic orbits de ned as
 z(t) = z(t+ Tf ) (5.3)
 the cost function
 J( ) =
 1
 2
 eTWe (5.4)
 can be minimized to  nd the the initial state and control to start a trajectory on this
 56
limit cycle. In this cost function, the error vector
 e =
 2
 6
 6
 6
 6
 4
 z(Tf ) z(0)
  (Tf )  (0)
 u(Tf ) u(0)
 w(Tf ) w(0)
 q(Tf ) q(0)
 3
 7
 7
 7
 7
 5
 (5.5)
 is the di erence between the longitudinal states after one  ap cycle, and the unknown
 parameter vector
  =
 2
 6
 6
 6
 6
 4
  (0)
 u(0)
 w(0)
 q(0)
   lon
 3
 7
 7
 7
 7
 5
 (5.6)
 consists of the initial longitudinal states and longitudinal trim control setting. The
 simplex search optimization method was employed to minimize the cost function (5.4),
 where the weighting matrix
 W = diagf1; 10000; 10; 100; 100g (5.7)
 accounted for the relative magnitudes in the error vector and produced adequate re-
 sults. Figure 5.2 shows an example of the errors and the total cost for a solution,
 which required approximately 200 iterations before convergence, although solutions
 were close after about 50 iterations. This is a simple method to trim to the or-
 nithopter; other more complicated methods in the literature include multiple shoot-
 ing algorithms [22] and monitoring mean values of  ight variables after an initial
 launching force vector is supplied [25].
 The ornithopter was trimmed for a range of  apping frequencies between 4 Hz and
 10 Hz. Below this range the ornithopter could not be trimmed, and above this range
 the dynamics model loses validity. Flight variables are shown as a function of  apping
 frequency in Figure 5.3. The longitudinal tail de ection becomes more negative as
 the wings  ap faster in order to counter an increasing pitching moment. Although
 this de ection adds more drag, even more thrust is produced as the  apping rate
 increases, which in turn increases the mean forward speed. The faster wing movements
 also generate a lower mean pitch angle and oscillation, as the inertial e ects cause a
 more rapid pitching motion that has less time to travel. Interestingly, the pitch rate
 amplitude remains approximately constant with  apping frequency because the tail is
 adjusted to compensate for the mean change in the pitching moment, which has also
 been reported by Lee et al. [25]. In addition, the altitude oscillation decreases and
 the heave velocity oscillation increases with increasing  apping frequency. Similar to
 the pitching dynamics, this is due to a faster excitation, behind which the response
 lags.
 The main features of the trim condition consist of the heave dynamics, involving
 the altitude and heave velocity, and the pitch dynamics, involving the pitch angle and
 57
e
 1 [m
 ]
 e
 2 [ra
 d]
 e
 3
 [m
 /s]
 e
 4
 [m
 /s]
 e
 5
 [ra
 d/
 s]
 J
 ?
 10
 3
 iteration number
 0.0
 0.4
 0.8
 0.0
 0.3
 0.6
 0.0
 1.5
 3.0
 0.0
 1.5
 3.0
 0.0
 5.0
 10.0
 0.0
 4.0
 8.0
 0 25 50 75 100 125 150 175 200
 Figure 5.2: Optimization results for  nding limit cycle trim trajectories
 ? l
 o
 n
 [ra
 d]
 u? [m
 /s]
 ?z [m
 ]
 ?w [m
 /s]
 ? ? [ra
 d]
 ?? [ra
 d]
 ?q [ra
 d/
 s]
 ff [Hz]
 -1.5
 -1.0
 -0.5
 +7.5
 +8.0
 +8.5
 +0.0
 +0.2
 +0.4
 +0.0
 +2.0
 +4.0
 +0.0
 +0.1
 +0.2
 +0.0
 +0.3
 +0.6
 +12
 +14
 +16
 3 4 5 6 7 8 9 10 11
 Figure 5.3: Variations in trim characteristics with  apping frequency
 58
pitch rate, which are shown in Figure 5.4. The characteristics shown in Figure 5.3
 are evident in these plots, as well as a change in shape, or bifurcation, of the limit
 cycles around 7 Hz. A similar phenomenon has been reported by Lee et al. [25], who
 claimed that simpli ed models may be possible for these higher  apping frequencies.
 The remainder of this chapter will deal with a single  apping frequency of 5.91
 Hz, to correspond with the  ight test used for system identi cation of the aerody-
 namics model. The trimmed  ight trajectory for this  apping frequency is shown in
 Figure 5.5, which has similar trajectories as those reported in the literature [50, 25].
 Due to the instability of the trim solution, integration over more than approximately
  ve wing strokes will tend to diverge from the limit cycle due to the accumulation
 of numerical errors. The trimmed model resulted in lower longitudinal velocity mean
 and oscillations than measured in the visual tracking  ight test, which is due to the
 di culty in modeling the thrust and drag of the wings. The altitude and heave oscil-
 lations, as well as the pitch oscillations and rates were consistent with those measured
 in  ight tests. The simulation data also has much lower frequency content than the
  ight data, due to the unmodeled  uid/structure interactions, actuator dynamics,
 and simpli ed aerodynamics.
 5.3 Numerical Linearization about Straight and
 Level Mean Flight
 Within the range of its validity, a linearized model of the  ight dynamics facilities
 the use of modern linear analysis tools, which can provide physical insight into the
 system and be used to design control laws to guarantee performance metrics. Ad-
 ditionally, linear models can be incorporated with standard commercial autopilots.
 Towards these goals, a linear model of the ornithopter  ight dynamics
 _x = Ax + Bu (5.8)
 is desired, where
 x = z z 
u =     (5.9)
 59
w
 [m
 /s]
 z [m]
 -2.0
 -1.0
 +0.0
 +1.0
 +2.0
 -0.2 -0.1 +0.0 +0.1 +0.2
 4
 5
 6
 7 8 9 10
 (a) heave dynamics
 q
 [ra
 d/
 s]
 ? [rad]
 -10.0
 -5.0
 +0.0
 +5.0
 +10.0
 -0.3 +0.0 +0.3 +0.6
 4 7 8 9 10
 (b) pitch dynamics
 Figure 5.4: Limit cycle oscillations for  apping frequencies between 4 Hz and 10 Hz
 60
x
 [m
 ]
 u
 [m
 /s]
 z
 [m
 ]
 w
 [m
 /s]
 ?
 [ra
 d]
 q
 [ra
 d/
 s]
 ? w
 [ra
 d]
 t [s]
 ? ? w
 [ra
 d/
 s]
 t [s]
 +0.00
 +0.70
 +1.40
 +7.00
 +8.00
 +9.00
 -0.10
 +0.00
 +0.10
 -1.00
 +0.00
 +1.00
 -0.25
 +0.00
 +0.25
 -10.0
 +0.0
 +10.0
 -1.00
 +0.00
 +1.00
 0.00 0.05 0.10 0.15 0.20
 -30.0
 +0.0
 +30.0
 0.00 0.05 0.10 0.15 0.20
 Figure 5.5: Trimmed  ight trajectory solution
 are state and control perturbations away from the trim values. To parallel the models
 currently employed in autopilots, the reduced state and control vectors are
 z =
 2
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 4
  
u
 w
 q
  
v
 p
 r
 3
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 5
  =
  
 lon
  lat
  
: (5.10)
 Linear models may be determined using a variety of methods. The nonlinear model
 may be linearized analytically using Jacobian matrices, but this method is prohibitive
 for complex and large-order systems. Alternatively, system identi cation experiments
 can be performed on the nonlinear simulation, but this method may require several
 iterations or long records of data to obtain adequate information content. In this work,
 the nonlinear model is linearized numerically, using  nite di erences, which requires
 little analysis and can be automated. A nonlinear system (5.1) can be numerically
 linearized about some point or trajectory in the state space z and some control   .
 61
Let
  1 = f(z +  iei;  )
  2 = f(z   iei;  )
  3 = f(z + 2 iei;  )
  4 = f(z  2 iei;  ) (5.11)
 where  i is a small perturbation and ei is an elementary vector. Then the ith column
 of the system state matrix A is computed
 8( 1   2) ( 3   4)
 12 i
 (5.12)
 and is accurate to O( 3i ). A similar process, where instead each control variable
 is perturbed, is performed to  nd columns of the system control matrix B. The
 perturbations  i are typically decreased iteratively until the system matrices converge
 [36].
 5.3.1 Linear Time-Invariant Model
 Conventional aircraft  ight dynamics models used for control design are typically
 linear time-invariant (LTI) models, where the system matrices are constant in time.
 The ornithopter  ight dynamics can be cast into an LTI model by using averaging
 theory [77, 39, 37, 20]. The cycle-averaged state derivatives (5.11) are then
  1 = (1=Tf )
 Z Tf
 0
 f(z +  iei;  )dt
  2 = (1=Tf )
 Z Tf
 0
 f(z   iei;  )dt
  3 = (1=Tf )
 Z Tf
 0
 f(z + 2 iei;  )dt
  4 = (1=Tf )
 Z Tf
 0
 f(z  2 iei;  )dt (5.13)
 and the remainder of the process is the same. Typically when this technique is applied,
 a time scale separation exists in the dynamics. Since the  apping frequency is close
 to the expected rigid body modes of the ornithopter, this method is questionable,
 as is the practice of using commercial autopilots, based on LTI models, with an
 62
ornithopter. The LTI model resulting from numerical linearization with averaging is
 2
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 4
  _ 
 _u
  _w
  _q
  _ 
 _v
  _p
  _r
 3
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 5
 =
 2
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 4
  0:69 +0:01  0:83 +0:84  0:04  0:01 +0:06  0:03
  7:38  0:51 +2:42  2:59 +0:06  0:01  0:21 +0:09
 +0:33  0:72  0:59 +0:41 +0:04 +0:03  0:01  0:03
  2:64  2:03  12:4 +0:74 +0:26 +0:12 +0:66  0:62
  0:00  0:00  0:00  0:00  0:50  0:12 +0:87 +0:42
  0:00 +0:00  0:00  0:01 +8:65  0:60 +1:65  6:27
 +0:00  0:01  0:08 +0:00 +4:79 +1:87  0:14  3:27
  0:01 +0:01  0:00  0:00  0:13 +0:06 +1:67 +0:49
 3
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 5
 2
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 4
   
 u
  w
  q
   
 v
  p
  r
 3
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 5
 +
 2
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 4
  3:10 +0:27
 +10:0  0:82
  2:09 +0:32
  44:5 +3:37
  0:02 +0:05
  0:03  3:20
  0:30  2:91
  0:02 +4:84
 3
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 5
 "
   lon
   lat
 #
 : (5.14)
 where states represent perturbations from an averaged trim value. While linearized
 lateral dynamics are presented in this chapter for completeness, their validity is ques-
 tionable, as lateral motions were not excited during the  ight test used to model the
 wing aerodynamics in Chapter 4.
 There are several points of interest in these matrices as compared to those of a
  xed-wing aircraft, for example the F-16 model in Appendix D. This ordering of the
 state vector illustrates that the state matrix decouples in much the same manner as
 conventional aircraft in straight and level  ight. The  rst and second 4  4 block-
 diagonal matrices represents the longitudinal and lateral dynamics, respectively. Here
 however, the pitch rate dynamics however still have signi cant coupling with the
 lateral dynamics, which is not common to  xed-wing aircraft and is due to the wings
  apping. The block-diagonal matrices are also more heavily populated than models
 for  xed-wing aircraft: typically the equations for the pitch angle and bank angle
 derivatives involve only the pitch rate and roll rate; however, for the ornithopter,
 they involve many states. The equations for the forward and side velocity typically
 involve the acceleration of gravity; here the signs are correct but the value is slightly
 lower than expected.
 The control matrix is also signi cantly populated and decouples in a similar man-
 ner as the state matrix. A positive longitudinal tail angle perturbation creates less
 drag and negative lift on the tail, resulting in a positive forward acceleration, a nega-
 tive heave acceleration, and a negative pitching moment. A positive lateral tail angle
 perturbation creates a negative side force, a negative roll rate, and a positive yaw
 rate. The controllability matrix for this system has full rank and is fully controllable.
 The modes of this system, described by Eigenvalues and Eigenvectors of the state
 matrix, are shown in Figures 5.6 and 5.7 and tabulated in Table 5.1. The longitudinal
 63
=
 (?
 )
 <(?)
 stable short period
 unstable short period
 spiral
 dutch roll
 roll
 -6
 -4
 -2
 0
 2
 4
 6
 -4 -2 0 2 4
 Figure 5.6: Linear time-invariant model pole locations
 dynamics include two oscillatory modes, one unstable and one stable, occurring at
 approximately the same frequency. The unstable mode is under-damped, while the
 stable mode is highly damped. The Eigenvectors of these modes indicate that they
 both are short period type motions, involving a fast and strong interplay of the pitch
 rate and forward velocity. Unlike conventional aircraft, a phugoid mode does not
 appear. The classical spiral, dutch roll, and roll modes are present in the decoupled
 lateral dynamics. The spiral mode is  rst order and unstable, involving mostly the
 sidewards velocity. The dutch roll mode is unstable and is strongly dependent on the
 sideways velocity. The roll mode is a  rst order mode, strongly dependent on the roll
 rate.
 5.3.2 Linear Time-Periodic Model
 There is not a time scale separation between the periodicity of the wings  ap-
 ping and the rigid body modes found in the LTI model, which calls into question
 the validity of using LTI models. A linear time-periodic (LTP) model results from
 computing (5.11) at each time step over the wing stroke, so that the system matrices
 vary with time. Elements of the system matrices are shown for the longitudinal and
 lateral decoupled systems in Figures 5.8(a) and 5.8(b), respectively.
 From observations of the frequency content, a functional representation of the
 matrix elements was obtained using the model
 y =  0 +  1 sin(2 ff t) +  2 sin(4 ff t) +  3 cos(2 ff t) (5.15)
 64
(0/2)pi
 (1/2)pi
 (2/2)pi
 (3/2)pi
 ??
 ?u
 ?w
 ?q
 (a) stable short period
 (0/2)pi
 (1/2)pi
 (2/2)pi
 (3/2)pi
 ??
 ?u
 ?w
 ?q
 (b) unstable short period
 (0/2)pi
 (1/2)pi
 (2/2)pi
 (3/2)pi
 ??
 ?v
 ?p
 ?r
 (c) spiral
 (0/2)pi
 (1/2)pi
 (2/2)pi
 (3/2)pi
 ???v
 ?p ?r
 (d) dutch roll
 (0/2)pi
 (1/2)pi
 (2/2)pi
 (3/2)pi
 ?? ?v?p
 ?r
 (e) roll
 Figure 5.7: Linear time-invariant model Eigenvector polar plots
 65
Table 5.1: Modal parameters of the decoupled linear time-invariant model
 Mode Eigenvalue Eigenvector Damping Frequency
 Ratio [rad/s]
  i vi  i !i
 stable
 short period
  2:07 1:57j
 2
 6
 6
 4
 +0:17 0:19j
  0:56 0:32j
  0:11 0:00j
  0:72
 3
 7
 7
 5 +0:80 2:60
 unstable
 short period
 +1:55 2:75j
 2
 6
 6
 4
  0:14 0:06j
 +0:74
  0:18 0:13j
  0:36 0:50j
 3
 7
 7
 5  0:49 3:16
 spiral  0:50
 2
 6
 6
 4
 +0:58
  0:18
  0:36
 +0:71
 3
 7
 7
 5 +1:00 0:50
 dutch roll +0:98 1:68j
 2
 6
 6
 4
 +0:07 0:25j
 +0:49 0:23j
 +0:58
 +0:18 0:53j
 3
 7
 7
 5  0:51 1:94
 roll  2:23
 2
 6
 6
 4
 +0:26
 +0:65
  0:61
 +0:37
 3
 7
 7
 5 +1:00 2:23
 66
and applying an equation-error analysis, where f gi are unknown, constant coe -
 cients. The functions consists of a bias and harmonic functions that are Tf -periodic
 and capture the major features of the responses. Simulation data and curve  ts are
 shown in Figure 5.8, whereas parameter estimates, standard error bounds, and coe -
 cients of determination for each matrix element are provided in Table 5.2. Generally
 the function (5.15)  ts all the elements well with high coe cients of determination
 and relatively low standard errors. Whereas the curve  ts for the lateral dynamics
 retain most of the frequency content, some elements in the longitudinal model, for
 instance a41 and a43, contain frequency content up to 10 times the  apping frequency
 and is not captured by the curve  t. This higher frequency content in the longitudinal
 dynamics is due to the longitudinal  apping of the wings.
 For a linear time varying system with x 2 <n 1, controllability at time t requires
 that the matrix
 C(t) =
  
M0(t) M1(t) ::: Mn 1(t)
  
(5.16)
 have rank n, where in general
 Mk+1(t) =  A(t)Mk(t) + _Mk(t) (5.17)
 M0(t) = B(t) (5.18)
 for k 2 [0; n 2]. Employing (5.15) and the identi ed parameters, computation of the
 controllability matrix for a wing stroke shows that both the longitudinal and lateral
 decoupled dynamics are fully controllable for all time.
 Given the functional representation of the LTP system (5.15), the state transition
 matrix can be computed by assembling the responses to initial conditions xi(0) = ei
 as a fundamental matrix
  (Tf ; 0) =
  
x1(Tf ) x2(Tf ) x3(Tf ) x4(Tf )
  
(5.19)
 stemming from a Floquet decomposition [78]. The decoupled longitudinal and lateral
 systems result in the analogous LTI system matrices
 Alon =
 2
 6
 6
 4
 +0:75 +0:06  0:01 +0:07
  0:96 +0:73  0:21  0:15
 +0:02 +0:07 +1:05 +0:03
  2:16 +0:16  1:04 +0:20
 3
 7
 7
 5 (5.20)
 Alat =
 2
 6
 6
 4
 +0:53  0:06 +0:05 +0:40
 +1:05 +1:00 +0:21  0:84
  1:95 +0:69 +0:25 +1:48
  1:20  0:38  0:10 +2:08
 3
 7
 7
 5 (5.21)
 from which the stability can be ascertained using a discrete Eigenvalue analysis.
 Poles are are tabulated in Table 5.3 and shown in Figure 5.9. The longitudinal
 dynamics exhibit a stable oscillatory mode, a stable subsidence mode, and an unstable
 divergence mode. The oscillatory mode looks most like a short period mode. The
 67
a
 1
 1
 a
 1
 2
 a
 1
 3
 a
 1
 4
 b 1
 a
 2
 1
 a
 2
 2
 a
 2
 3
 a
 2
 4
 b 2
 a
 31
 a
 32
 a
 33
 a
 34 b 3
 a
 41
 t [s]
 a
 42
 t [s]
 a
 43
 t [s]
 a
 44
 t [s]
 b 4
 t [s]
 -4.00
 -2.00
 +0.00
 -1.50
 +0.00
 +1.50
 -1.00
 +0.00
 +1.00
 +0.60
 +0.90
 +1.20
 -1.00
 +0.00
 +1.00
 +2.00
 -12.0
 -8.0
 -4.0
 -12.0
 +0.0
 +12.0
 -15.0
 +0.0
 +15.0
 -3.00
 +0.00
 +3.00
 -12.0
 -6.0
 +0.0
 -15.0
 +0.0
 +15.0
 -30.0
 -15.0
 +0.0
 -25.0
 +0.0
 +25.0
 -2.00
 +0.00
 +2.00
 -4.00
 +0.00
 +4.00
 -100
 -50
 +0
 0.0 0.1 0.2
 -80.0
 +0.0
 +80.0
 0.0 0.1 0.2
 -60.0
 +0.0
 +60.0
 0.0 0.1 0.2
 -25.0
 +0.0
 +25.0
 0.0 0.1 0.2
 -60.0
 +0.0
 +60.0
 0.0 0.1 0.2
 (a) decoupled longitudinal dynamics
 a
 1
 1
 a
 1
 2
 a
 1
 3
 a
 1
 4
 b 1
 a
 2
 1
 a
 2
 2
 a
 2
 3
 a
 2
 4
 b 2
 a
 31
 a
 32
 a
 33
 a
 34 b 3
 a
 41
 t (s)
 a
 42
 t (s)
 a
 43
 t (s)
 a
 44
 t (s)
 b 4
 t (s)
 -3.00
 +0.00
 +3.00
 -2.00
 +0.00
 +2.00
 +0.00
 +1.00
 +2.00
 +0.00
 +1.00
 +2.00
 +0.00
 +0.30
 +0.60
 +0.0
 +6.0
 +12.0
 -5.00
 +0.00
 +5.00
 +0.00
 +1.00
 +2.00
 -8.00
 -4.00
 +0.00
 -1.00
 +0.00
 +1.00
 -80.0
 -40.0
 +0.0
 -60.0
 +0.0
 +60.0
 -20.0
 -10.0
 +0.0
 +0.0
 +25.0
 +50.0
 +0.0
 +15.0
 +30.0
 -6.00
 -3.00
 +0.00
 0.0 0.1 0.2
 -4.00
 +0.00
 +4.00
 0.0 0.1 0.2
 -6.00
 +0.00
 +6.00
 0.0 0.1 0.2
 +0.0
 +7.0
 +14.0
 0.0 0.1 0.2
 +0.0
 +7.0
 +14.0
 0.0 0.1 0.2
 (b) decoupled lateral dynamics
 Figure 5.8: Numerically linearized (solid) and curve  tted (dashed) decoupled system
 matrix element time histories over one wing stroke
 68
Table 5.2: Estimated parameters and standard errors for the decoupled linear time-
 periodic model
 element  ^0  s( ^0)  ^1 + s( ^1)  ^2  s( ^2)  ^3  s( ^3) R2
 a11  0:44 0:05 +0:19 0:07 +0:46 0:07 +0:31 0:07 0.28
 a21  9:03 0:10  2:31 0:15  0:52 0:15  0:25 0:14 0.62
 a31  0:18 0:30  0:70 0:43  1:23 0:43  0:67 0:43 0.07
 a41  18:7 1:71 +4:23 2:43 +12:8 2:43 +14:0 2:42 0.28
 a12 +0:13 0:02  0:62 0:03 +0:05 0:03  0:44 0:02 0.85
 a22  1:24 0:07 +4:91 0:10  0:83 0:10 +1:75 0:10 0.95
 a32  0:39 0:24 +7:63 0:34  2:23 0:34 +0:86 0:34 0.77
 a42 +1:89 1:62 +14:3 2:29  17:1 2:29  27:5 2:28 0.59
 a13 +0:21 0:02  0:40 0:03 +0:34 0:03  0:07 0:03 0.66
 a23  1:53 0:09 +8:02 0:12  1:83 0:12 +3:03 0:12 0.97
 a33  0:12 0:27 +3:92 0:38  1:83 0:38 +0:03 0:38 0.43
 a43  11:3 1:30 +6:82 1:85  8:52 1:85  14:0 1:84 0.36
 a14 +0:79 0:00 +0:19 0:01  0:04 0:01 +0:12 0:01 0.87
 a24  1:14 0:02  0:60 0:03 +0:04 0:03  0:81 0:03 0.87
 a34 +0:35 0:02  0:90 0:03 +0:07 0:03 +0:51 0:03 0.90
 a44  5:50 0:30  0:21 0:43  2:10 0:43 +7:76 0:42 0.68
 b1 +0:55 0:02  0:90 0:03 +0:19 0:03  0:17 0:03 0.88
 b2  5:55 0:09 +3:71 0:13 +0:11 0:13 +3:69 0:12 0.91
 b3 +0:78 0:03 +0:77 0:04 +0:06 0:04  2:05 0:04 0.94
 b4  15:0 0:86  12:6 1:22 +15:2 1:22  36:5 1:21 0.88
 a11  0:90 0:01 +1:33 0:01  0:55 0:01 +0:66 0:01 0.99
 a21 +8:21 0:07 +1:99 0:09 +0:35 0:09 +0:17 0:09 0.74
 a31  28:0 0:66  0:61 0:94 +1:13 0:94 +33:5 0:93 0.89
 a41  1:43 0:03 +1:83 0:05  0:78 0:05 +1:45 0:05 0.94
 a12  0:22 0:01 +1:23 0:01  0:07 0:01 +0:49 0:01 0.98
 a22  0:76 0:06 +2:06 0:09 +0:49 0:09 +0:07 0:09 0.78
 a32 +4:41 0:71  19:6 1:00 +8:24 1:00 +46:6 1:00 0.94
 a42 +0:04 0:07 +1:30 0:10 +0:20 0:10 +0:96 0:10 0.61
 a13 +0:89 0:00 +0:27 0:00  0:06 0:00 +0:09 0:00 1.00
 a23 +1:24 0:02 +0:20 0:02 +0:15 0:02 +0:00 0:02 0.45
 a33  5:94 0:12 +1:00 0:18 +0:82 0:18 +8:62 0:18 0.94
 a43 +1:47 0:04  3:38 0:06 +0:59 0:06  0:61 0:06 0.96
 a14 +0:67 0:01  0:61 0:01 +0:06 0:01  0:36 0:01 0.97
 a24  5:99 0:04  1:16 0:05  0:58 0:05 +0:36 0:05 0.80
 a34 +20:9 0:43 +4:60 0:61 +0:02 0:61  22:6 0:61 0.90
 a44 +1:79 0:26  0:76 0:36 +4:68 0:36  3:30 0:36 0.60
 b1 +0:25 0:00  0:20 0:01 +0:01 0:01  0:12 0:01 0.91
 b2  0:24 0:01  0:26 0:02  0:26 0:02 +0:28 0:02 0.85
 b3 +9:27 0:13  0:36 0:19 +0:73 0:19  8:58 0:19 0.93
 b4 +6:64 0:12 +0:17 0:17 +2:46 0:17  1:95 0:16 0.69
 69
=
 (?
 )
 <(?)
 ?1,2
 ?3
 ?4
 ?5
 ?6
 ?7
 ?8
 -1
 -0.5
 0
 0.5
 1
 0 0.5 1 1.5
 Figure 5.9: Linear time-periodic model pole locations
  rst order modes involve most of the state variables and in that regard are consistent
 with previous  ndings [42, 22, 50], although these modes are much slower. The
 lateral dynamics contain two subsidence modes and two divergence modes. One of
 the subsidence appears to be a roll mode, and one of the divergence appears as a
 spiral mode.
 5.4 Modeling Implications for Control
 In this chapter, three models have been presented: a nonlinear multibody model,
 a LTP model for straight and level mean  ight, and a LTI model for straight and
 level mean  ight.
 The nonlinear model contains the full set of modeled dynamics, including the mass
 variation, inertial couplings, and nonlinear aerodynamics, from which trim solutions
 and simpli ed models can be determined. This model is of the highest  delity, and
 could be used to explore perching maneuvers and agile  ight. Although it is di cult
 to design controllers for nonlinear systems, this model has been cast into a canonical
 form for Euler-Lagrange systems so that techniques such as passivity control can
 be applied in a straight forward manner [64, 62, 63, 66]. Additionally, techniques
 such as sliding mode or nonlinear damping could be used to attenuate the pitch
 oscillation, so that video images are more stable and electronic components experience
 less vibration; similar control strategies have shown promise [24, 20, 33]. The heave
 motion, however, cannot be attenuated without morphing the wing structure to more
 70
Table 5.3: Modal parameters of the decoupled linear time-periodic model
 Eigenvalue Magnitude Eigenvector Damping Frequency
 Ratio [rad/s]
  i k ik2 vi  i !i
 +1:04 1.04
 2
 6
 6
 4
 +0:13
  0:17
  0:78
 +0:59
 3
 7
 7
 5  1:00 0.21
 +0:60 0.60
 2
 6
 6
 4
  0:20
  0:47
 +0:03
 +0:86
 3
 7
 7
 5 +1:00 3.01
 +0:54 0:41j 0.68
 2
 6
 6
 4
  0:09 0:16j
 +0:36 0:07j
  0:07 0:04j
 +0:91
 3
 7
 7
 5 +0:52 4.45
 +1:04 1.04
 2
 6
 6
 4
 +0:59
 +0:19
 +0:14
 +0:77
 3
 7
 7
 5  1:00 0.23
 +0:70 0.70
 2
 6
 6
 4
 +0:41
  0:03
  0:86
 +0:29
 3
 7
 7
 5 +1:00 2.11
 +1:84 1.84
 2
 6
 6
 4
 +0:28
  0:42
 +0:24
 +0:83
 3
 7
 7
 5  1:00 3.60
 +0:28 0.28
 2
 6
 6
 4
  0:13
  0:18
 +0:97
  0:07
 3
 7
 7
 5 +1:00 7.52
 71
k 1
 1
 k 1
 2
 k 1
 3
 k 1
 4
 t [s]
 -10
 -5
 +0
 +5
 +10
 -12
 -6
 +0
 +6
 +12
 -2.0
 -1.0
 +0.0
 +1.0
 +2.0
 -4.0
 -2.0
 +0.0
 +2.0
 +4.0
 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00
 Figure 5.10: Periodic control gains designed for the LTP system using LQR
 gracefully generate the lift and thrust forces.
 The LTP model describes the dynamics of perturbations about a trimmed straight
 and level mean  ight trajectory. This model is still complex and must be determined
 and scheduled for each  ight condition. However, controllers are simple enough for
 implementation on microprocessors, the dynamics still re ect the true motions, and
 analytical tools can be applied, e.g. the controllability result. This model could be
 used for an accurate model for navigation, or for inner-loop stabilization. Straight
 forward synthesis tools, such as a periodic Linear Quadratic Regulator (LQR) result,
 also exist [79]. For instance, solving LQR problem, using unity weighting matrices,
 results in the periodic controller gains shown in Figure 5.10. For an impulse applied
 at any point within the  apping cycle, the state responses settle within 0.6 seconds.
 Dietl and Garcia [50] have also shown that a discrete-time, in nite-horizon version
 of this control works well in simulation and suggests it is practical for  ight control,
 despite control computations, controller gain memory requirements, and knowledge
 of the wing angle.
 The simplest model presented is the LTI model, which describes the stroke-
 averaged dynamics of perturbations away from trimmed  ight trajectories. This
 model enables a large set of analysis and synthesis tools. However, models must still
 be determined and scheduled for each  ight condition, but more importantly these
 models average the signi cant e ects of the wings  apping. For that reason, it should
 not be trusted for fast or accurate tasks, and should be relegated to slow, forgiving
 tasks such as navigation. Krashanitsa et al. were able to demonstrate navigation us-
 ing a heuristically-tuned proportional/derivative control law [11]; however, the  ight
 72
trajectories exhibited a large amount of variation between laps. Dietl and Garcia
 noted that although classical controllers are easy to implement, the system dynam-
 ics are periodically changing, so that these controllers are sometimes improving the
  ight performance, and sometimes degrading the  ight performance. To control the
 pitch angle based on the identi ed LTI model, the proportional/integral/derivative
 controller
  lon(t) =  6:67  (t) + 0:27
 Z t
 0
   ( )d  2:03 _ (t) (5.22)
 can be used, but only predicts a settling time of 32.7 seconds. Applying the LQR
 controller, again solved with unity weighting matrices,
  lon(t) =  0:23  (t) + 0:14 u(t) 0:02 w(t) 0:15 q(t) (5.23)
 can more e ectively capture cross couplings and predicts settling times under three
 seconds. Another possibility is to use this simple model in conjunction with a robust
 control scheme, such as H1 or  -synthesis, or rather a robust adaptive control law
 such as L1, to learn how to compensate for the changing system dynamics [51].
 5.5 Chapter Summary
 This chapter presented an investigation into the dynamics of the identi ed or-
 nithopter model using simulation tools. The model was programmed into both For-
 tran and Matlab simulation environments, which o ered trade o s between com-
 putational speed and the availability of analysis tools. A time step of 10 4 s was
 determined necessary for convergence with a fourth order Runge-Kutta routine.
 The periodic forcing of the wings admits trimmed  ight trajectories that are limit
 cycle oscillations in the state space. Simulations from initial conditions suggest that
 these limit cycles are unstable, and a simple method was introduced for locating initial
 conditions and trimmed control settings that would put the ornithopter state on the
 limit cycle. This method was used to trim the ornithopter for straight and level mean
  ight for a variety of  apping frequencies. As the  apping frequency increased, the
 tail angle became more negative, the mean forward speed increased, the mean pitch
 angle and pitch oscillations decreased, the altitude variation decreased and the heave
 velocity oscillations increased. Additionally a bifurcation is evident around 7 Hz, due
 to a shape change in the pitch dynamics and heave dynamics limit cycles.
 Linear time-periodic and linear time-invariant models were numerically deter-
 mined from the nonlinear model. The LTI model exhibited a stable and unstable
 short period mode, while the LTP model showed a stable short period mode, a stable
 subsidence mode, and an unstable divergence mode.
 Finally, implications for control strategies were discussed. For high precision or
 rapid maneuvering control, the nonlinear model should be used. Although not a
 straight forward process, the model has been cast into a canonical form for nonlin-
 ear control. The LTP model is suitable for accurate regulation of trimmed  ight
 conditions, and the LTI model is only suitable for simple tasks such as heading or
 73
altitude-hold navigation.
 74
Chapter 6
 Concluding Remarks
 6.1 Summary of Work
 Flight dynamics models are needed for  apping-wing aircraft, so that the dy-
 namics can be better understood and  ight controllers can be designed to perform
 mission scenarios ranging from robust outdoor navigation to indoor agile  ight and
 perching. This work contributes a nonlinear model of the  ight dynamics, as well as
 an investigation into linearizing the model and the rami cations on control design.
 Chapter 2 presented characterizations of an ornithopter  ight test platform. The
 aircraft design and operation is  rst discussed. Then an experiment is presented
 where the ornithopter was  tted with a custom avionics package to measure state
 variables, in real time, during trimmed straight and level mean  ight. Afterwards, the
 variation in the mass distribution is examined using computer aided design software.
 An experiment is then summarized where aerodynamic and con guration data were
 measured in a constrained quasi-hover condition. The implications for  ight dynamics
 modeling are then discussed.
 The results of Chapter 2 suggest a nonlinear multibody model of the ornithopter
 vehicle dynamics. Chapter 3 begins with a review of the available tools for obtaining
 multibody models. Based on this review, a Lagrangian type model was selected and
 the derivation is presented using the Boltzmann-Hamel equations. The equations of
 motion are then cast into a canonical form for nonlinear control and the model is
 discussed.
 Chapter 4 presents the system identi cation work completed to identify parts
 of the ornithopter  ight dynamics model not known a priori. A review of system
 identi cation techniques and methods is presented. Afterwards, a wind tunnel test is
 discussed where lift and drag values were calculated for the ornithopter tail. Finally
 an experiment is presented using a visual tracking system to obtain  ight test data
 with which an aerodynamic model for the wings was identi ed.
 Chapter 5 presents  ight simulation results obtained using the  ight dynamics
 model developed in Chapters 3 and 4. A method for  nding trim solutions is presented
 and was employed to trim the ornithopter model for straight and level mean  ight.
 The model was then linearized about the trim trajectory, resulting in a canonical
 75
time-invariant model as well as a time-periodic model.
 6.2 Summary of Modeling Assumptions
  Vehicle dynamics are approximated by three interconnected bodies
 The ornithopter multibody model consisted of a fuselage and two wings, which
 comprise 95.1% of the mass and captures the inertial e ects of the  apping
 wings. Extra bodies for the tail surface were not modeled, as the tail has only
 4.9% of the mass, moves slowly relative to the  apping wings, and does not
 create a signi cant change in the mass distribution. Instead, the tail mass and
 inertia was absorbed into the fuselage. Similarly, no additional bene t was
 observed from modeling the four-bar mechanism used to  ap the wings.
  Structural dynamics in the vehicle dynamics are ignored
 Structural deformations have been observed in the fuselage, wing spars, and
 wing sail membrane; however, these bodies are assumed to be rigid in order to
 develop models for control. Torsion and lateral bending modes in the fuselage
 are at higher frequencies than used in control. Additionally, tip displacements of
 the wing spars are relatively small, less than 8% of the wingspan [17]. Finally,
 although the aft section of the wings de ect to produce thrust, this portion
 contains less than 2.3% of the vehicle mass and does not signi cantly vary the
 mass distribution.
  Actuator dynamics are negligible
 It is assumed that the pilot or autopilot directly controls the  apping frequency
 and orientation of the tail. This bypasses the actuator dynamics, presented in
 Appendix B, and thus any transmission delays, time lags, backlash, mechanical
 slop, and current/torque saturation limits. These models can be easily incorpo-
 rated to increase model  delity. Momentum sources from the spinning motors
 and gears are assumed negligible.
  Tail aerodynamics are steady and do not interact with the wings
 The downwash from the wings onto the tail is assumed negligible due to the
 placement of the tail between and above the wings. Additionally, limited actu-
 ator and pilot bandwidths led to steady models of the tail aerodynamics. The
 full e ect of downwash from the wings, unsteady aerodynamic e ects from the
 body pitching motion, and  exibility in mechanisms is not currently known.
 6.3 Summary of Original Contributions
  Flight data suitable for modeling
 This work presented  ight data of an ornithopter, sampled at a bandwidth
 that illuminated the e ects of the  apping wings. Trimmed  ight data showed
 76
pitch oscillations up to 4.97 rad/s and heave accelerations up to 41.7 m/s2,
 which require nonlinear multibody models to capture. Additionally, harmonic
 responses were seen in the data, indicating either structural modes or nonlinear-
 ities. High accelerations preclude the use of traditional attitude determination
 methods and require model-based techniques.
  Modeling of mass distribution variation due to wings  apping
 The e ects of the  apping motion on the mass distribution of the ornithopter
 was investigated using a CAD model. The center of mass travels 10.3 cm in
 the vertical direction, the moments of inertia can vary up to 53.6%, and the
 inertia rates are signi cant. Multibody models are required to capture these
 variations.
  Derivation of rigid multibody vehicle dynamics
 It was determined that a three-body model was su cient to model the or-
 nithopter. Vehicle dynamics were derived using the Boltzmann-Hamel equa-
 tions, and were cast into a canonical form used for the nonlinear control of
 Euler-Lagrange systems.
  Tail aerodynamics system identi cation
 Wind tunnel tests were conducted on the ornithopter tail, having free stream
 velocities, angles of attack, and angles of sideslip ranging between 4.50 to 6.33
 m/s,  0:88 to  0:35 rad, and  0:57 rad, respectfully. System identi cation
 methods were applied to determine a models for the aerodynamic coe cients.
  Flight testing and wing aerodynamics system identi cation
 The ornithopter  ight data, obtained using a visual tracking system, was pre-
 sented. This data and the vehicle dynamics model were used to explore vari-
 ations in angle of attack, Reynolds number, reduced frequency, velocity, and
 structural deformations over the wing and throughout the wing stroke cycle.
 Additionally, force and moment sources were identi ed throughout the wing
 stroke. Finally, this data was used to determine a wing aerodynamics model
 obtained from  ight data using system identi cation methods.
  Models for stability and control
 In addition to the full nonlinear model, a linear time-invariant model and a
 linear time-periodic model was presented, with modeling implications and ram-
 i cations for control.
 6.4 Recommendations for Future Research
  Aerodynamic model expansion
 The aerodynamic model presented in Chapter 4 is valid for a relatively small
 region of the  ight envelope. The tail was tested in steady conditions for  ow
 77
velocities between 4.50 and 6.33 m/s, angles of attack between  0:88 and  0:35
 rad, and angles of sideslip between  0:57 rad. The model could be improved
 by increasing these ranges and be looking at the dynamic e ects of actuating
 the tail.
 The  ight test was conducted for a  apping frequency of 5.91 Hz forward speeds
 between 9 and 10 m/s. For these tests the tail was  xed, and so the only
 dynamics excited were a result of the  apping wings. More testing should
 be performed where the lateral dynamics are excited, the throttle changed,
 and the tail excited, in order to expand the aerodynamic model of the wings.
 Additionally, longer sets of data are needed to excite the forward speed dynamics
 in order to model the thrust force more accurately. In order to complete these
 tasks, a larger  ight space is needed. This could be accomplished using a large
 building instrumented with numerous tracking cameras, or perhaps by  ying
 outdoors with more instrumentation; however, problems still arise when using
 air data booms and trying to solve the attitude problem due to the pitching
 and heaving motion.
  Wing/tail aerodynamic interactions
 In this work, it was assumed that the wing and tail aerodynamics do not in-
 teract. Although the tail is centered between and located above the wings, it
 is currently unknown how much downwash from the wing impinges on the tail.
 Flow visualization or particle image velocimetry experiments would provide this
 information, which would then be used to update the aerodynamic model struc-
 ture of the tail. Additionally, the  ight test mentioned above where the wings
 and tail were both excited could also provide this insight.
  Aeroelastic e ects
 It is known that aeroelastic e ects play a role in the  ight dynamics of or-
 nithopters; however, this work focus on models for  ight control design and
 neglected the aeroelastic e ects. The nonlinear model  delity could be im-
 proved by incorporating  nite element models to capture the  uid/structure
 interactions, as in Lee et al. [25], so that the e ect of aeroelasticity could be
 assessed.
  Control law design
 Control laws can now be designed from the models presented in this work. It
 remains then to implement these in  ight controllers and evaluate their perfor-
 mance. This capability would allow direct comparison of controllers developed
 from the LTI and LTP models, and would facilitate the development of perching
 and aggressive maneuvering experiments.
  Integrated platform for morphing
 Provided a  ight dynamics model and a stabilized  ight platform, work on
 morphing wing structures could progress. Morphing would allow more degrees
 78
of freedom, with which the ornithopter could perform more e cient and agile
  ight.
 79
Appendix A
 Field Calibration of Inertial
 Measurement Units
 Inertial measurement units (IMUs) consist of orthogonal triads of magnetometers,
 gyroscopes, and accelerometers, and are used in autopilots to provide information on
 the vehicle pose and motion. The IMUs used on unmanned aircraft su er from
 the low signal-to-noise ratios and strong temperature dependencies typical of Micro-
 electromechcanical systems (Mems) devices, requiring tedious laboratory testing or
 on-board state estimators for calibration. In this appendix, a method is presented for
 estimating IMU calibration constants, based on three successive least-squares compu-
 tations. Memory and computation requirements are low enough for implementation
 on 8 bit microprocessors, and experimental data shows adequate calibration results.
 A.1 Theory and Method
 The IMU provides measurements of the magnetic  eld h0d, linear acceleration a
 0
 d,
 and rotational velocity !0d, corrected to the aircraft body frame [55, 35]. The subscript
 d is used to denote measurements in raw digital format, e.g. an integer between 0 and
 65535. Assuming a constant temperature, as unmanned aircraft have limited  ight
 durations and altitude ceilings, the calibration model is
 y0 =  y0d + b + n (A.1)
 where y is the measurement,  = diagf x;  y;  zg is a matrix of scale factors, b is a
 vector of biases, and n is a vector of additive measurement noise.
 A block diagram of the calibration process is shown in Figure A.1. Data is collected
 from a maneuver where quick angular deviations are made to excite each sensor
 channel without introducing large translational accelerations. Existing least-square
 calibration algorithms can be used to calibrate the magnetometer and accelerometer
 [80, 81], which can then estimate the vehicle attitude [58, 82].
 Substituting the sensor model (A.1), the rotational kinematic equation of motion
 is
 _ I =  
 
 !!0d + b! + n!
  
(A.2)
 80
Inertial
 Measurement
 Unit
 Magnetometer
 Calibration
 Accelerometer
 Calibration
 Attitude
 Estimation
 Gyroscope
 Calibration
 h0d
 a0d
 !0d
 h0
 a0
  I
 !0
 Figure A.1: IMU calibration method block diagram
 Table A.1: Measurement speci cations for avionics and visual positioning system
 Sensor Symbol Range Resolution Error Units
 time t - 25:6 10 6 - s
 magnetometer h0  3:17 48:3 10 6 5:83 10 3 -
 gyroscope !0  5:24 0:16 10 3 0:14 10+0 rad/s
 accelerometer a0  98:1 3:88 10 3 0:05 10+0 m/s2
 attitude  I - - 0:03 10 3 rad
 which is a least-squares problem that solves for the gyroscope calibration parameters
 using the attitude estimates in a certainty-equivalence fashion. Using a quaternion
 to represent attitude, (A.2) can be solved quickly and accurately by using only the
 scalar portion of the quaternion.
 A.2 Results
 The avionics package shown in Figure 2.4, was used to record measurements from
 an IMU, through a 16 bit analog to digital converter, at 185 Hz. A visual tracking
 system was employed to obtain high accuracy estimates of the avionics orientation, de-
 rived from spatial measurements of ten retro-re ective markers placed on the avionics
 board. Hardware speci cations for the IMU and the visual tracking system are given
 in Table A.1. The maneuver chosen for analysis was a sequence of rotational doublets
 along each sensing axes, as shown by the digitized measurements in Figure A.2.
 The least-squares routine developed by Gebre-Egziabher [80] was used to cali-
 brate the magnetometers and accelerometers. Field domain calibrations are shown
 in Figure A.3. Vehicle attitude was then computed using the Quest algorithm, and
 equation-error was used to match the quaternion, as shown in Figure A.4. Estimated
 parameters and standard errors are given in Table A.2. Magnetometer scale factors
 81
+20
 +30
 +40
 +50
 h
 0d
 x
 [co
 u
 nts]
 ?
 10
 3
 +20
 +30
 +40
 +50
 h
 0d
 y
 [co
 u
 nts]
 ?
 10
 3
 +20
 +30
 +40
 +50
 h
 0d
 z
 [co
 u
 nts]
 ?
 10
 3
 t[s]
 (a
 )
 magnetomete
 r
 +20
 +30
 +40
 +50
 a
 0d
 x
 [co
 u
 nts]
 ?
 10
 3
 +20
 +30
 +40
 +50
 a
 0d
 y
 [co
 u
 nts]
 ?
 10
 3
 +20
 +30
 +40
 +50
 a
 0d
 z
 [co
 u
 nts]
 ?
 10
 3
 t[s]
 (b
 )
 acceleromete
 r
 +20
 +30
 +40
 +50
 ?
 0d
 x
 [co
 u
 nts]
 ?
 10
 3
 +20
 +30
 +40
 +50
 ?
 0d
 y
 [co
 u
 nts]
 ?
 10
 3
 +20
 +30
 +40
 +50
 ?
 0d
 z
 [co
 u
 nts]
 ?
 10
 3
 t[s]
 (c
 )
 gyrosco
 p
 e
 Figur
 e
 A.2
 :
 Uncalibrate
 d
 16
 bi
 t
 IM
 U
 measureme
 nt
 s
 fro
 m
 a
 roll-pit
 ch-
 ya
 w
 maneu
 ve
 r
 82
-1.0
 -0.5
 +0.0
 +0.5
 +1.0
 h0x
 [?]
 -1.0
 -0.5
 +0.0
 +0.5
 +1.0
 h0y
 [?]
 -1.0
 -0.5
 +0.0
 +0.5
 +1.0
 h0z
 [?]
 (a) magnetometer
 -10
 -5
 +0
 +5
 +10
 a0x
 [m/s2]
 -10
 -5
 +0
 +5
 +10
 a0y
 [m/s2]
 -10
 -5
 +0
 +5
 +10
 a0z
 [m/s2]
 (b) accelerometer
 Figure A.3: Calibrated magnetometer and accelerometer signals with centered spheres
 had a high error bound, perhaps due to the small amount of data used, but otherwise
 all estimated parameters had less than a 10% predicted error, indicating an accurate
 set of measurements.
 The method presented in this section is a batch method, executed after the data
 has been collected. Alternatively, the least-squares routines could be implemented in
 a recursive fashion, as in [35, 83], so that the calibration constants and their associated
 error bounds evolve in real-time as more data becomes available. Estimating parame-
 ters in real-time would also allow for variations in the calibration due to temperature
 or  ight condition changes over long periods of time.
 83
-2
 -1
 +0
 +1
 +2
 0 2 4 6 8 10
 ? ?
 [ra
 d/
 s]
 t [s]
 data
 model
 Figure A.4: Model  t to the attitude kinematic equation using equation-error
 Table A.2: IMU calibration estimates and standard errors
 Parameter Equation-Error Multiplier
   ^  ( ^)
  hx +0:16 2:05 10 3
  hy +0:14 1:89 10 3
  hz +0:16 2:04 10 3
 bhx  6:36 0:40 10+0
 bhy  4:98 0:36 10+0
 bhz +5:47 0:37 10+0
   x +0:50 0:02 10+0
   y +0:34 0:01 10+0
   z +0:37 0:02 10+0
 b x  15:6 3:17 10+0
 b y  11:7 2:74 10+0
 b z  12:3 2:81 10+0
  !x +0:27 0:01 10 3
  !y +0:37 0:04 10 3
  !z +0:33 0:02 10 3
 b!x +8:58 0:46 10+0
 b!y  12:5 1:20 10+0
 b!z  10:2 0:76 10+0
 84
Appendix B
 Actuator Dynamics System
 Identi cation
 The ornithopter has on board a Feigao GF2030 brushless motor to drive the
  apping motion, and two Hitec HS-56 servo motors to orient the tail. Due to time
 lags observed in  ight, models for the actuator dynamics were identi ed to determine
 actuator bandwidths and to enable higher  delity models of the  ight dynamics. This
 section presents the system identi cation experiments and modeling results.
 B.1 Experimental Setup
 The actuators were mounted to a reaction torque sensor for measuring the gen-
 erated torque. The DC motor was  tted with a shaft encoder for a rotary speed
 measurement, while the servo motor internal potentiometer was measured for a ro-
 tary position measurement. A Pic18F452 microprocessor was used to generate input
 waveforms and to measure data at 500 Hz using 16 bit analog to digital converters.
 The experimental setup is shown in Figure B.1 and measurement speci cations are
 provided in Table B.1.
 The data used for system identi cation are shown in Figure B.2. Both actuators
 were excited with a digital waveform that had an increasing and decreasing frequency.
 Table B.1: Actuator system identi cation measurement speci cations
 Measurement Range Resolution Error Unit
 time - 0:05 10 3 - s
 throttle input 0.00 to 1.00 0:25 10 3 0:07 10 3 -
 angle input  0:70 10+0 0:35 10 3 0:16 10 3 rad
 velocity  3:77 10+3 8:47 10 3 0:03 10+0 rad/s
 angle  0:70 10+0 0:16 10 3 0:13 10+0 rad
 torque  0:35 10+0 0:77 10 6 0:90 10 3 N m
 85
(a) DC motor (b) servo motor
 Figure B.1: Actuators and instrumented test stand
 Table B.2: DC motor parameter estimates and standard errors
 Method Mass Damping Sti ness Fit
 m 10 3 c 10 3 k  10+0 R2
 time equation-error +0:28 0:01 +0:92 0:06  64:1 3:08 0.69
 frequency equation-error +0:23 0:01 +0:89 0:02  58:4 1:30 0.78
 time output-error +0:22 0:01 +1:01 0:03  64:6 2:07 0.98
 frequency output-error +0:20 0:01 +0:92 0:04  58:0 2:98 0.99
 B.2 Results
 DC and servo motors are typically linear transducers. Step-wise regression in the
 time domain using equation-error showed that a second order system
 m q + c _q + kq +  = b (B.1)
 was su cient to model the systems. The mass encompasses the shaft inertia; damping
 describes e ects of friction, back electro-mechanical force, current/torque coupling,
 and electrical resistance; sti ness models the position feedback used in the servo
 motor;  is the shaft torque; and  is the pilot stick input. Higher  delity models
 such as Coulomb friction and third-order dynamics [84, 85, 62] made no signi cant
 increase in model  delity.
 Equation-error and output-error parameter identi cation methods were applied in
 both the time and frequency domains, as shown in Figure B.3. Estimated parameters
 are provided in Table B.2 and Table B.3. All models had coe cients of determina-
 tion of 0.70 or above, and output-error parameter estimates were within statistical
 agreement. The DC motor has a bandwidth of 0.48 Hz and the servo motor 3.20 Hz.
 86
+0.12
 +0.13
 +0.14
 +0.15
 +0.16
 ?
 [-]
 -4.00
 -3.50
 -3.00
 -2.50
 -2.00
 q?
 [ra
 d/
 s]
 ?
 10
 3
 -2.00
 -1.00
 +0.00
 +1.00
 +2.00
 0 5 10 15 20
 ?
 [N
 ?
 m
 ]?
 10
 ?
 3
 t [s]
 (a) DC motor
 -1.0
 -0.5
 +0.0
 +0.5
 +1.0
 ?
 [ra
 d]
 -1.0
 -0.5
 +0.0
 +0.5
 +1.0
 q
 [ra
 d]
 -20
 -10
 +0
 +10
 +20
 0 2 4 6 8 10 12
 ?
 [N
 ?
 m
 ]?
 10
 ?
 3
 t [s]
 (b) servo motor
 Figure B.2: Measurements used for actuator system identi cation
 87
-10
 -5
 +0
 +5
 +10
 q?
 [ra
 d/
 s2
 ]?
 10
 3
 0.0
 2.5
 5.0
 7.5
 10.0
 q?
 [ra
 d/
 s2
 ]?
 10
 3
 -4.0
 -3.5
 -3.0
 -2.5
 -2.0
 0 5 10 15 20
 q?
 [ra
 d/
 s]
 ?
 10
 3
 t [s]
 0.0
 2.0
 4.0
 6.0
 0 2 4 6
 q?
 [ra
 d/
 s]
 ?
 10
 3
 f [Hz]
 data
 model
 (a) DC motor
 -0.6
 -0.3
 +0.0
 +0.3
 +0.6
 q?
 [ra
 d/
 s2
 ]?
 10
 3
 0.0
 0.1
 0.2
 0.3
 0.4
 q?
 [ra
 d/
 s2
 ]?
 10
 3
 -1.0
 -0.5
 +0.0
 +0.5
 +1.0
 0 2 4 6 8 10 12
 q
 [ra
 d]
 t [s]
 0.0
 0.5
 1.0
 1.5
 2.0
 0 2 4 6 8 10 12
 q
 [ra
 d/
 s]
 f [Hz]
 data
 model
 (b) servo motor
 Figure B.3: Model  ts for actuator dynamics using equation-error and output-error
 in the time and frequency domains
 88
Table B.3: Servo motor parameter estimates and standard errors
 Method Mass Damping Sti ness Fit
 m 10 3 c 10 3 k  10+0 R2
 time equation-error +0:01 0:00 +0:15 0:01 +2:22 0:27 0.81
 frequency equation-error +0:03 0:00 +0:15 0:01 +2:15 0:03 0.96
 time output-error +0:04 0:00 +0:53 0:04 +7:12 0:82 0.98
 frequency output-error +0:04 0:00 +0:57 0:01 +5:58 0:14 0.97
 B.3 Coupling to Vehicle Dynamics
 The rigid body dynamic equations (3.18) can be written as
 Mb _v + Cbv + g + a =  : (B.2)
 The actuator dynamics (B.1) can be written in matrix form, in terms of the or-
 nithopter state variables as
 Ma _v + Cav + Ga + Kap = Ba (B.3)
 where Ga is a matrix of gear ratios, relating the actuator outputs to the ornithopter
 joint angles. Substituting (B.2) into (B.3) and rearranging yields
 (Ma + GaMb) _v + (Ca + GaCb)v + (Gag + Gaa + Kap) = Ba (B.4)
 which is the same canonical form as the equations of motion described in Chapter 3.
 89
Appendix C
 Equations of Motion for
 Single-Body Flight Vehicles
 Many  ight vehicles, including conventional  xed-wing aircraft, rotary-wing air-
 craft, and spacecraft, can often be modeled as a single body. Derivation of the
 equations of motion is most commonly performed using the Newton-Euler approach
 [34, 35]; however, to provide familiarity with the energy approach presented in Chap-
 ter 3, the equations of motion for a single-body vehicle are derived here using the
 Boltzmann-Hamel equations.
 Consider the single-body vehicle in Figure 2.1, having position and velocity state
 variables
 p =
  
rI0;I
  I0;I
  
v =
  
 00;I
 !00;I
  
: (C.1)
 Written in terms of these state variables, the kinetic energy
 T (p;v) =
 1
 2
 ( 00;I)
 T (m0I)( 00;I) +
 1
 2
 (!00;I)
 T I00(!
 0
 0;I) (C.2)
 is the summation of the system energy. Using algebraic manipulations, (C.2) can be
 rearranged into the quadratic form
 T (p;v) =
 1
 2
 vTMv (C.3)
 to admit the generalized mass matrix
 M =
  
m0I 0
 0 I00
  
: (C.4)
 As the mass matrix does not depend on the position states of the aircraft, the
 coupling matrix given by (3.22) reduces to
 C(p;v) = _M +
 nvX
 k=1
 @T
 @vk
  k (C.5)
 90
where the number of velocity states nv is six. Partial di erentiation of (C.3) with
 respect to the velocity states reveals that the derivative term in (C.5) refers to kth
 entry of the (nv  1) row vector vTM. The second term in (C.5) is the set of Hamel
 coe cient matrices  k, which can be computed using (3.8), (3.14), and (3.15) as
  1 =
 2
 6
 6
 6
 6
 6
 6
 4
 0 0 0 0 0 0
 0 0 0 0 0 1
 0 0 0 0  1 0
 0 0 0 0 0 0
 0 0 1 0 0 0
 0  1 0 0 0 0
 3
 7
 7
 7
 7
 7
 7
 5
  2 =
 2
 6
 6
 6
 6
 6
 6
 4
 0 0 0 0 0  1
 0 0 0 0 0 0
 0 0 0 1 0 0
 0 0  1 0 0 0
 0 0 0 0 0 0
 1 0 0 0 0 0
 3
 7
 7
 7
 7
 7
 7
 5
  3 =
 2
 6
 6
 6
 6
 6
 6
 4
 0 0 0 0 0 0
 0 0 0 0 0 1
 0 0 0 0  1 0
 0 0 0 0 0 0
 0 0 1 0 0 0
 0  1 0 0 0 0
 3
 7
 7
 7
 7
 7
 7
 5
  4 =
 2
 6
 6
 6
 6
 6
 6
 4
 0 0 0 0 0  1
 0 0 0 0 0 0
 0 0 0 1 0 0
 0 0  1 0 0 0
 0 0 0 0 0 0
 1 0 0 0 0 0
 3
 7
 7
 7
 7
 7
 7
 5
  5 =
 2
 6
 6
 6
 6
 6
 6
 4
 0 0 0 0 0 0
 0 0 0 0 0 1
 0 0 0 0  1 0
 0 0 0 0 0 0
 0 0 1 0 0 0
 0  1 0 0 0 0
 3
 7
 7
 7
 7
 7
 7
 5
  6 =
 2
 6
 6
 6
 6
 6
 6
 4
 0 0 0 0 0  1
 0 0 0 0 0 0
 0 0 0 1 0 0
 0 0  1 0 0 0
 0 0 0 0 0 0
 1 0 0 0 0 0
 3
 7
 7
 7
 7
 7
 7
 5
 (C.6)
 regardless of the attitude parametrization. Applying these terms yields the dynamic
 coupling matrix
 C(p;v) =
  
_m0I 0
 0 _I00
  
+
  
0  m0S( 00;I)
  m0S( 00;I)  S(I
 0
 0!
 0
 0;I)
  
(C.7)
 for a single body vehicle. Furthermore, if the vehicle can be assumed rigid so that it
 neither changes mass nor inertia, the mass rate matrix in (C.5) and (C.7) is null.
 The remaining forces on the vehicle are those imparted by the environment through
 which it passes. The e ects of gravity can be modeled by  rst writing the potential
 energy of the aircraft
 U(p) =  m0(rI0;I)
 TgI (C.8)
 and then using (3.26) to form the generalized gravitational forces vector
 g(p) =
  
m0R0;IgI
 0
  
: (C.9)
 Combining terms and adding aerodynamic contributions a(p;v), the full nonlinear
 91
dynamics of the rigid body are
 M(p) _v + C(p;v)v + g(p) + a(p;v) =  (C.10)
 where  is a vector of exogenous forces and torques. This parametrization of the
 dynamics satis es the property used in passivity control design that _M  2C is a
 skew-symmetric matrix [66, 59, 64].
 92
Appendix D
 Linearization of a Conventional
 Aircraft Model
 This section presents the linearization of a conventional aircraft model to provide
 a baseline for comparison with the results presented in Chapter 5. The aircraft chosen
 was the single-seat  ghter F-16 Fighting Falcon, which has a nominal 9279 kg mass
 and 9.96 m wing span. A nonlinear dynamics model is contained within Sidpac, and
 is documented in references [35, 36, 86].
 The aircraft was set at a 3048 m altitude and a 0.09 rad pitch angle. A modi ed
 Newton-Raphson algorithm was then used to solve (5.2) for the remaining trim states,
 which were a 137 m/s airspeed, a 19.5% throttle setting, and a  0:08 rad elevator
 de ection. The center of mass is at the mean aerodynamic quarter-chord location to
 synthetically increase open loop stability.
 Numerical linearization using  nite di erences resulted in the time-invariant model
 2
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 4
  _ 
 _u
  _w
  _q
  _ 
 _v
  _p
  _r
 3
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 5
 =
 2
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 4
 +0:00 +0:00 +0:00 +1:00 +0:00 +0:00 +0:00 +0:00
  9:77  0:01 +0:06  11:5 +0:00 +0:00 +0:00 +0:00
 +0:00  0:08  0:68 +127: +0:00  0:00 +0:00 +0:00
 +0:00 +0:00  0:03  1:15 +0:00 +0:00 +0:00 +0:00
 +0:00 +0:00 +0:00  0:00 +0:00 +0:00 +1:00 +0:00
 +0:00  0:00  0:00 +0:00 +9:77  0:21 +12:1  136:
 +0:00  0:00 +0:00  0:00 +0:00  0:16  2:37 +0:60
 +0:00  0:00 +0:00  0:00 +0:00 +0:05  0:05  0:35
 3
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 5
 2
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 4
   
 u
  w
  q
   
 v
  p
  r
 3
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 5
 +
 2
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 4
 +0:00 +0:00 +0:00
 +2:08 +0:00 +0:00
  11:2 +0:00 +0:00
  6:43 +0:00 +0:00
 +0:00 +0:00 +0:00
 +0:00 +1:54 +4:20
 +0:00  25:3 +4:28
 +0:00  1:11  2:28
 3
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 5
 2
 6
 4
   e
   a
   r
 3
 7
 5 (D.1)
 93
=
 (?
 )
 <(?)
 phugoid
 short period
 spiral
 dutch roll
 roll
 -3
 -2
 -1
 0
 1
 2
 3
 -3 -2 -1 0
 Figure D.1: Linearized F-16 model pole locations
 where  e is the elevator de ection,  a is the aileron de ection, and  r is the rudder
 de ection.
 In this  ight condition, the block-diagonal structure of the system matrices results
 in a decoupling of the longitudinal and lateral dynamics. Kinematic integrations are
 seen in the equations for the pitch angle and bank angle. The local acceleration due
 to gravity, 9.77 m/s2, appears in the equations for longitudinal and lateral velocity.
 The elevator primarily a ects the pitch rate, but also the longitudinal and heave
 velocities. The aileron and rudder primarily a ect the roll and yaw rate, respectively,
 but also couple into the other lateral states.
 Poles are shown in Figure D.1, and modal parameters of the decoupled subsystems
 are provided in Table D.1. The longitudinal dynamics have two oscillatory, under-
 damped modes. The phugoid mode is slower and involves primarily the forward
 velocity, whereas the short period mode is faster and involves primarily the heave
 velocity. The lateral dynamics comprise three modes. The spiral mode is slow,  rst-
 order mode that involves primarily the bank angle and side velocity. The roll mode
 is fast,  rst-order mode that involves mostly the roll rate. The dutch roll mode is a
 fast, second-order mode using mostly the side velocity.
 94
Table D.1: Modal parameters of the decoupled F-16 linear model
 Mode Eigenvalue Eigenvector Damping Frequency
 Ratio [rad/s]
  i vi  i !i
 phugoid  0:01 0:09j
 2
 6
 6
 4
  0:00 0:01j
 +0:99
 +0:05 0:00j
 +0:00 0:00j
 3
 7
 7
 5 +0:11 0:09
 short period  0:91 1:84j
 2
 6
 6
 4
 +0:01 0:00j
  0:07 0:03j
 +0:99
  0:00 0:01j
 3
 7
 7
 5 +0:44 2:06
 spiral  0:01
 2
 6
 6
 4
 +0:89
 +0:44
  0:02
 +0:06
 3
 7
 7
 5 +1:00 0:01
 dutch roll  0:35 2:82
 2
 6
 6
 4
 +0:01 0:01j
 +0:99
  0:03 0:04j
  0:00 0:02j
 3
 7
 7
 5 +0:12 2:84
 roll  2:22
 2
 6
 6
 4
  0:34
  0:56
 +0:75
 +0:03
 3
 7
 7
 5 +1:00 2:22
 95
Bibliography
 [1] T. Mueller and J. DeLaurier, An Overview of Micro Air Vehicle Aerodynamics,
 vol. 195 of Progress in Aeronautics and Astronautics, chapter 1, pp. 1{10, Aiaa,
 2001.
 [2] F. Bohorquez, D. Pines, and P. Samuel, \Small rotor design optimization using
 blade element momentum theory and hover tests", Journal of Aircraft, vol. 47,
 no. 1, pp. 268{283, January-February 2010.
 [3] U. Pesavento and Z. Wang, \Flapping wing  ight can save aerodynamic power
 compared to steady  ight", Physical Review Letters, vol. 103, no. 11, pp. 1{4,
 September 2009.
 [4] R. Malhan, M. Benedict, and I. Chopra, \Experimental investigation of a  ap-
 ping wing concept in hoverand forward  ight for micro air vehicle applications",
 Phoenix, AZ, May 2010, American Helicopter Society, 66th Annual Forum.
 [5] P. Seshadri, M. Benedict, and I. Chopra, \Experimental investigation of an
 insect-based  apping wing hovering micro air vehicle", San Francisco, CA, Jan-
 uary 2010, American Helicopter Society, Aeromechanics Specialists? Conference.
 [6] T. Liu, K. Kuykendoll, R. Rhew, and S. Jones, \Avian wings", Portland, OR,
 June 2004, Aiaa, number 2004-2186 in Aerodynamic Measurement Technology
 and Ground Testing Conference.
 [7] J. Videler, Avian Flight, Oxford Ornithology Series. Oxford University Press,
 2005.
 [8] D. Alexander, Nature?s Flyers: Birds, Insects, and the Biomechanics of Flight,
 The Johns Hopkins University Press, 2002.
 [9] J. Delaurier, \An ornithopter wing design", Canadian Aeronautics and Space
 Journal, vol. 40, no. 1, pp. 10{18, March 1994.
 [10] H. R abinger, \Die entwicklung der ev-schlag  ugelmodelle", Tech. Rep.,
 N urnberg, Germany, 2002.
 [11] R. Krashanitsa, D. Silin, S. Shkarayev, and G. Abate, \Flight dynamics of a
  apping-wing air vehicle", International Journal of Micro Air Vehicles, vol. 1,
 no. 1, pp. 35{49, 2009.
 96
[12] M. Keennon and J. Grasmeyer, \Development of the black widow and microbat
 mavs and a vision for the future of mav design", Dayton, OH, July 2003, Aiaa,
 number 2003-2901 in International Air and Space Symposium and Exposition
 conference.
 [13] Z. Jackowski, \Design and construction of an autonomous ornithopter", Master?s
 thesis, Massachusetts Institute of Technology, 2009.
 [14] E. Avadikian and J. Hubbard, \Ornithopter development and optimization
 guide", Tech. Rep., University of Maryland, 2010.
 [15] J. Leishman, Principles of Helicopter Aerodynamics, Cambridge University
 Press, 2nd edition, 2006.
 [16] R. Harmon, J. Grauer, J. Hubbard, J. Conroy, S. Humbert, J. Sitaraman, and
 B. Roget, \Experimental determination of ornithopter membrane wing shapes
 used for simple aerodynamic modeling", Aiaa, 2008, number 2008-6237 in Ap-
 plied Aerodynamics Conference.
 [17] R. Harmon, \Aerodynamic modeling of a  apping membrane wing using motion
 tracking experiments", Master?s thesis, University of Maryland, 2008.
 [18] T. Rashid, \The  ight dynamics of a full-scale ornithopter", Master?s thesis,
 University of Toronto, 1995.
 [19] C. Orlowski, A. Girard, and W. Shyy, \Open loop pitch control of a  apping
 wing micro-air vehicle using a tail and control mass", Baltimore, MD, June 2010,
 Ieee, American Controls Conference.
 [20] C. Orlowski, Flapping Wing Micro Air Vehicles: An Analysis of the Importance
 of the Mass of the Wings to Flight Dynamics, Stability, and Control, PhD thesis,
 University of Michigan, 2011.
 [21] J. Delaurier, \An aerodynamic model for  apping-wing  ight", Royal Aeronau-
 tical Society, pp. 125{130, 1993.
 [22] J. Dietl and E. Garcia, \Stability in ornithopter longitudinal  ight dynamics",
 Journal of Guidance, Control, and Dynamics, vol. 31, no. 4, pp. 1157{1162,
 July-August 2008.
 [23] J. Lee, J. Kim, D. Kim, and J. Han, \Longitudinal  ight dynamics of bio-inspired
 ornithopter considering  uid-structure interaction", Toronto, Ontario Canada,
 August 2010, Aiaa, Atmospheric Flight Mechanics Conference.
 [24] D. Kim, S. Lee, and J. Han, \Improved aerodynamic model for e cient analysis
 of  apping-wing  ight", Aiaa Journal, vol. 49, no. 4, pp. 868{872, April 2011.
 [25] J. Lee, J. Kim, D. Kim, and J. Han, \Longitudinal  ight dynamics of bioin-
 spired ornithopter considering  uid-structure interaction", Journal of Guidance,
 Control, and Dynamics, vol. 34, no. 3, pp. 667{677, May-June 2011.
 97
[26] D. Kim, H. Kim, J. Han, and K. Kwon, \Experimental investigation on the aero-
 dynamic characteristics of a bio-mimetic  apping wing with macro- ber compos-
 ites", Journal of Intelligent Material Systems and Structures, vol. 19, no. 3, pp.
 423{431, 2008.
 [27] M. Vest and J. Katz, \Unsteady aerodynamic model of  apping wings", Aiaa
 Journal, vol. 34, no. 7, pp. 1435{1440, July 1996.
 [28] B. Roget, J. Sitaraman, R. Harmon, J. Grauer, J. Hubbard, and S. Humbert,
 \Computational study of  exible wing ornithopter  ight", Journal of Aircraft,
 vol. 46, no. 6, pp. 2016{2031, November-December 2009.
 [29] W. Shyy, Y. Lian, J. Tang, D. Viieru, and H. Liu, Aerodynamics of Low Reynolds
 Number Flyers, Cambridge Aerospace Series. Cambridge University Press, 2007.
 [30] D. Keats-Pullen, J. Hubbard, and N. Guerreiro, \Spatial weighting of smart ma-
 terials for real-time measurement of aerodynamic forces", Journal of Intelligent
 Material Systems and Structures, vol. 19, no. 7, pp. 837{844, July 2007.
 [31] M. Bolender, \Rigid multi-body equations-of-motion for  apping wing mavs
 using kane?s equations", Chicago, IL, August 2009, Aiaa, Guidance, Navigation,
 and Control Conference.
 [32] K. Sibilski, \Dynamics of micro-air-vehicle with  apping wings", Acta Polytech-
 nia, vol. 44, no. 2, pp. 15{21, 2004.
 [33] J. Grauer and J. Hubbard, \Modeling of ornithopter  ight dynamics for state
 estimation and control", Baltimore, MD, June 2010, Ieee, American Controls
 Conference.
 [34] D. McRuer, I. Ashkenas, and D. Graham, Aircraft Dynamics and Automatic
 Control, Princeton University Press, 1973.
 [35] V. Klein and E. Morelli, Aircraft System Identi cation: Theory and Practice,
 Aiaa, Aiaa Education Series, 2006.
 [36] B. Stevens and F. Lewis, Aircraft Control and Simulation, Wiley, 2003.
 [37] D. Doman, M. Oppenheimer, and D. Sigthorsson, \Wingbeat shape modulation
 for  apping-wing micro-air-vehicle control during hover", Journal of Guidance,
 Control, and Dynamics, vol. 33, no. 3, pp. 724{739, May-June 2010.
 [38] W. Dickinson, A. Straw, and M. Dickinson, \Integrative model of drosophila
  ight", Aiaa Journal, vol. 46, no. 9, pp. 2150{2164, September 2008.
 [39] X. Deng, L. Schenato, W. Wu, and S. Sastry, \Flapping- ight for biomimetic
 robotic insects: Part i - system modeling", Ieee Transactions on Robotics, vol.
 22, no. 4, pp. 776{788, August 2006.
 98
[40] S. Chung and M. Dorothy, \Neurobiologically inspired control of engineered
  apping  ight", Journal of Guidance, Control, and Dynamics, vol. 33, no. 2, pp.
 440{453, March-April 2010.
 [41] X. Deng, L. Schenato, and S. Sastry, \Flapping- ight for biomimetic robotic
 insects: Part ii -  ight control design", Ieee Transactions on Robotics, vol. 22,
 no. 4, pp. 789{803, August 2006.
 [42] G. Taylor and A. Thomas, \Dynamic  ight stability in the desert locust schis-
 tocerca gregaria", Journal of Experimental Biology, pp. 2803{2829, 2003.
 [43] G. Taylor, R. Bomphrey, and J. Hoen, \Insect  ight dynamics and control",
 Reno, NV, Jan. 2006, Aiaa, number 2006-32 in Aerospace Sciences Meeting and
 Exhibit.
 [44] I. Faruque and J. Humbert, \Dipteran insect  ight dynamics: Part i - longitu-
 dinal motion about hover", Journal of Theoretical Biology, pp. 1{15, 2010.
 [45] I. Faruque, Control-Oriented Reduced Order Modeling Dipteran Flapping Flight,
 PhD thesis, University of Maryland, 2011.
 [46] J. Kim, J. Lee, and J. Han, \Limit-cycle oscillation suppression of ornithopter
 longitudinal  ight dynamics", Portland, OR, Aug. 2011, Aiaa, number 2011-
 6404 in Guidance, Navigation, and Control Conference.
 [47] J. Lee, J. Kim, and J. Han, \Cycle-averaged mean and zero-mean time-varying
  ight dynamics of the ornithopters at trim", Smart Systems and Structures Lab,
 KAIST, 2010.
 [48] G. Taylor and R. Zbikowski, \Nonlinear time-periodic models of the longitudinal
  ight dynamics of desert locusts shistocerca gregaria", Journal of the Royal
 Society Interface, vol. 5, pp. 197{221, May 2005.
 [49] M. Bolender, \Open-loop stability of  apping  ight in hover", Toronto, Ontario
 Canada, August 2010, Aiaa, Guidance, Navigation, and Control Conference.
 [50] J. Dietl and E. Garcia, \Ornithopter control with periodic in nite horizon con-
 trollers", Journal of Guidance, Control, and Dynamics, vol. 34, no. 5, pp. 1412{
 1422, Sep.-Oct. 2011.
 [51] R. Tedrake, Z. Jackowski, R. Cory, J. Roberts, and W. Hoburg, \Learning to  y
 like a bird", under review, 2010.
 [52] B. Etkin and L. Reid, Dynamics of Flight: Stability and Control, Wiley, 1996.
 [53] S. Kinkade, \Hobby technik", www.hobbytechnik.com, 2010.
 [54] J. Grauer and J. Hubbard, \Inertial measurements from  ight data of a  apping-
 wing ornithopter", Journal of Guidance, Control, and Dynamics, vol. 32, no. 1,
 pp. 326{331, Jan.-Feb. 2009.
 99
[55] T. Gainer and S. Ho man, \Summary of transformation equations and equations
 of motion used in free- ight and wind tunnel data reduction and analysis", Tech.
 Rep. SP-3070, Nasa, 1972.
 [56] E. Morelli, \High accuracy evaluation of the  nite fourier transform using sam-
 pled data", Tech. Rep. 110340, Nasa, June 1997.
 [57] J. Sitaraman and J. Baeder, \On the  eld velocity approach and geometric
 conservation law for unsteady  ow simulations", 2003, number 2003-3835 in
 Applied Aerodynamics Conference.
 [58] M. Shuster and S. Oh, \Three-axis attitude determination from vector observa-
 tions", Journal of Guidance and Control, vol. 4, no. 1, pp. 70{77, 1981.
 [59] J. Slotine and W. Li, Applied Nonlinear Control, Prentice Hall, 1991.
 [60] E. Whittaker, A Treatise on the Analytical Dynamics of Particles and Rigid
 Bodies, Cambridge University Press, 1917.
 [61] D. Greenwood, Advanced Dynamics, Cambridge University Press, 2003.
 [62] F. Lewis, D. Dawson, and C. Abdallah, Robot Manipulator Control: Theory and
 Practice, Marcel Dekker, 2004.
 [63] R. Kelly, V. Santib a nez, and A. Lor  a, Control of Robot Manipulators in Joint
 Space, Springer, 2005.
 [64] R. Ortega and A. Loria, Passivity-Based Control of Euler-Lagrange Systems,
 Springer, 1998.
 [65] J. Brewer, \Kronecker products and matrix calculus in system theory", Ieee
 Transactions on Circuits and Systems, vol. 25, no. 9, pp. 772{781, September
 1978.
 [66] I. Fantoni and R. Lozano, Nonlinear Control for Underactuated Mechanical
 Systems, Springer, 2002.
 [67] L. Ljung, System Identi cation: Theory for the User, Information and System
 Sciences Series. Prentice Hall, 1999.
 [68] M. Tischler and R. Remple, Aircraft and Rotorcraft System Identi cation: Engi-
 neering Methods with Flight Test Examples, Aiaa Eduation Series. Aiaa, 2006.
 [69] V. Klein, J. Batterson, and P. Murphy, \Determination of airplane model struc-
 ture from  ight data using modi ed stepwise regression", Tech. Rep. TP-1916,
 Nasa, October 1981.
 [70] R. Maine and K. Ili , \Application of parameter estimation to aircraft stability
 and control: The output-error approach", Tech. Rep. TM 1168, Nasa, 1986.
 100
[71] Anonymous, \Flying qualities of piloted aircraft", Tech. Rep. MIL-STD-1797,
 US Department of Defense, 1995.
 [72] J. Grauer, J. Conroy, J. Hubbard, J. Humbert, and D. Pines, \System identi ca-
 tion of a miniature helicopter", Journal of Aircraft, vol. 46, no. 4, pp. 1260{1269,
 July-Aug. 2009.
 [73] J. Grauer, E. Ulrich, J. Hubbard, D. Pines, and J. Humbert, \Model structure
 determination of an ornithopter aerodynamics model from  ight data", Orlando,
 FL, January 2010, Aiaa, Aerospace Sciences Meeting.
 [74] D. Montgomery, E. Peck, and G. Vining, Introduction to Linear Regression
 Analysis, Wiley, 2001.
 [75] D. Belsley, E. Kuh, and R. Welsh, Regression Diagnostics: Identifying In uential
 Data and Sources of Collinearity, Wiley, 1980.
 [76] Y. Xiong and M. Sun, \Stabilization control of a bumblebee in hovering and
 forward  ight", Acta Mech Sin, vol. 25, pp. 13{21, 2008.
 [77] H. Khalil, Nonlinear Systems, Prentice Hall, 2011.
 [78] W. Rugh, Linear System Theory, Information and system sciences. Prentice
 Hall, second edition, 1996.
 [79] P. Dorato, C. Abdallah, and V. Cerone, Linear Quadratic Control: An Intro-
 duction, Krieger, 2000.
 [80] D. Gebre-Egziabher, G. Elkaim, D. Powell, and B. Parkinson, \Calibration
 of strapdown magnetometers in magnetic  eld domain", Journal of Aerospace
 Engineering, vol. 19, no. 2, pp. 87{102, April 2006.
 [81] Q. Li, \Least squares ellipsoid speci c  tting", 2004, Ieee paper No. 0-7695-
 2078-2.
 [82] J. Crassidis, F. Markley, and Y. Cheng, \Survey of nonlinear attitude estimation
 methods", Journal of Guidance, Control, and Dynamics, vol. 30, no. 1, pp. 12{
 28, January-February 2007.
 [83] E. Morelli, \Real-time dynamic modeling: Data information requirements and
  ight test results", Journal of Aircraft, vol. 46, no. 6, pp. 1894{1905, Nov.-Dec.
 2009.
 [84] L. Meirovitch, Fundamentals of Vibrations, McGraw-Hill, 2001.
 [85] T. Tarn, A. Bejczy, X. Yun, and Z. Li, \E ect of motor dynamics on nonlinear
 feedback robot arm control", Ieee Transactions on Robotics and Automation,
 vol. 7, no. 1, pp. 114{122, February 1991.
 101
[86] L. Nguyen, M. Ogburn, W. Gilbert, K. Kibler, P. Brown, and P. Deal, \Sim-
 ulator study of stall/post-stall characteristics of a  ghter airplane with relaxed
 longitudinal static stability", Tech. Rep. TP 1538, Nasa, 1979.
 102