MODELING, DESIGN AND CONTROL OF GLIDING ROBOTIC FISH
By
Feitian Zhang

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Electrical Engineering – Doctor of Philosophy
2014

ABSTRACT
MODELING, DESIGN AND CONTROL OF GLIDING ROBOTIC FISH
By
Feitian Zhang
Autonomous underwater robots have been studied by researchers for the past half century. In
particular, for the past two decades, due to the increasing demand for environmental sustainability,
significant attention has been paid to aquatic environmental monitoring using autonomous underwater robots. In this dissertation, a new type of underwater robots, gliding robotic fish, is proposed
for mobile sensing in versatile aquatic environments. Such a robot combines buoyancy-driven
gliding and fin-actuated swimming, inspired by underwater gliders and robotic fish, to realize both
energy-efficient locomotion and high maneuverability. Two prototypes, a preliminary miniature
underwater glider and a fully functioning gliding robotic fish, are presented. The actuation system
and the sensing system are introduced. Dynamic model of a gliding robotic fish is derived by integrating the dynamics of miniature underwater glider and the influence of an actively-controlled
tail. Hydrodynamic model is established where hydrodynamic forces and moments are dependent
on the angle of attack and the sideslip angle. Using the technique of computational fluid dynamics
(CFD) water-tunnel simulation is carried out for evaluating the hydrodynamic coefficients. Scaling
analysis is provided to shed light on the dimension design.
Two operational modes of gliding robotic fish, steady gliding in the sagittal plane and tailenabled spiraling in the three-dimensional space, are discussed. Steady-state equations for both
motions are derived and solved numerically. In particular, for spiral motion, recursive Newton’s
method is adopted and the region of convergence for this method is numerically examined. The
local asymptotic stability of the computed equilibria is established through checking the Jacobian

matrix, and the basins of attraction are further numerically explored. Simulation and experiments
are conducted to validate steady-state models and calculated equilibria for both motions.
Tail-enabled feedback control strategies are studied in both sagittal-plane glide stabilization
and three-dimensional heading maintenance. A passivity-based controller and a sliding mode controller are designed and tested in experiments for those two problems, respectively. In sagittalplane glide stabilization, a nonlinear observer is designed and implemented to estimate velocityrelated states. A three-dimensional curve tracking problem is also discussed and a two-degree-offreedom control scheme is proposed by integrating static inverse mapping and H∞ control technique. The differential geometric features, such as the torsion and curvature, are explored for
planning the trajectory.
Finally, the field tests with the lab-developed prototype of gliding robotic fish are conducted
in the Kalamazoo River, Michigan and the Wintergreen Lake, Michigan for detecting oil spill and
sampling harmful algal blooms, respectively. Both gliding and spiraling motions are tested in the
experiments as well as the fish-like swimming. The field test results are presented to show the
effectiveness of the designed robot in environmental monitoring tasks.

Dedicated to
my wife Mi Zhou
and our son
Andrew Anbang Zhang
with all my love.

iv

ACKNOWLEDGMENTS

First, I would like to express my warmest gratitude to my advisor, the director of Smart Microsystems Lab at Michigan State University, Prof. Xiaobo Tan. Five years ago, Prof. Tan provided
me a valuable opportunity of starting my Ph.D. journey here in the field of underwater robotics,
and now upon the completion of my Ph.D. training, with experience and knowledge gained during
the past five years, I am excitedly looking forward to my postdoctoral research life at University of
Maryland. I genuinely thank him for his advising in my research and mentoring in my career path.
I want to thank my academic committee members, Prof. Hassan Khalil, Prof. Guoming Zhu
and Prof. Mukherjee at Michigan State University, and our research collaborator, Prof. Fumin
Zhang at Georgia Institute of Technology, for their insightful comments in the area of control and
robotics. I also thank Prof. Elena Litchman at Michigan State University and Prof. Derek Paley at
University of Maryland for their generosity and kindness in providing experimental facilities and
testing fields for the project.
In Smart Microsystems Lab, I owe many people a debt of gratitude. The research in gliding
robotic fish needs diverse expertise and extensive collaboration. I want to thank John Thon and
Cody Thon for their continuous contribution in mechanical design and robot manufacturing, and
Jianxun Wang and Hong Lei for their help in circuit development, and Bin Tian, Osama En-Nasr,
Scott O’Connor, Suriya Madhan, Tingyuan Zhang and Yujie Hao for their help in conducting
experiments and collecting data.
I would also like to thank the other students in Smart Microsystems Lab and the faculty, students and post-docs in my department, whose comments and discussions have helped my research.
Thanks also to all of the administrators and staff members in the ECE department who have helped
me over the years.
v

I also want to acknowledge the support from the funding agencies NSF (IIS 0916720, ECCS
1050236, IIS 1319602, IIP 1343413, CCF 1331852) that makes this work possible.
Most of all, I want to thank my mother, Dongmei He, my father, Guangbo Zhang and my wife,
Mi Zhou for their unconditional love in both good times and bad times. I could not have finished
my Ph.D. program in the past five years without their understanding, support and encouragement.

vi

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

x

Chapter 1
Introduction . . . . . . . . . . . . . . . .
1.1 Technology in Aquatic Environmental Monitoring
1.2 Gliding Robotic Fish . . . . . . . . . . . . . . .
1.2.1 Design Concept . . . . . . . . . . . . . .
1.2.2 Motion and Control . . . . . . . . . . . .
1.3 Overview of Contributions . . . . . . . . . . . .

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

. 1
. 1
. 3
. 3
. 6
. 12

Chapter 2
Implementation of Gliding Robotic Fish
2.1 Actuation System . . . . . . . . . . . . . . . . .
2.2 Gliding Robotic Fish Components . . . . . . . .
2.3 Mechanical Design . . . . . . . . . . . . . . . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

14
14
16
19

Chapter 3
Dynamic Model of Gliding Robotic Fish . . . . . .
3.1 Full Dynamic Model . . . . . . . . . . . . . . . . . . . . .
3.2 Hydrodynamic Model . . . . . . . . . . . . . . . . . . . . .
3.2.1 Hydrodynamic Equations . . . . . . . . . . . . . . .
3.2.2 CFD-based Evaluation of Hydrodynamic Coefficients

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

22
22
26
26
27

Chapter 4
Steady Gliding in the Sagittal Plane . . . . . . .
4.1 Reduced Model in the Sagittal Plane . . . . . . . . . . .
4.2 Computation of Steady Gliding Path in the Sagittal Plane
4.3 Scaling Analysis . . . . . . . . . . . . . . . . . . . . . .
4.4 Experimental Results and Model Validation . . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

36
36
37
40
44

Chapter 5
Steady Spiral in Three-Dimensional Space . . . . . . . . . . . .
5.1 Steady-State Spiraling Equations . . . . . . . . . . . . . . . . . . . . . .
5.2 Computation of Spiral Path and Evaluation of Stability . . . . . . . . . .
5.2.1 Newton’s method for Solving the Steady-State Spiraling Equations
5.2.2 Region of Convergence for Newton’s Method . . . . . . . . . . .
5.2.3 Stability Analysis of Spiraling Motion . . . . . . . . . . . . . . .
5.2.4 Basins of Attraction for the Spiraling Dynamics . . . . . . . . . .
5.3 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

55
55
58
58
60
63
64
67

.
.
.
.
.

.
.
.
.
.

Chapter 6
Passivity-based Stabilization with a Whale-type Tail . . . . . . . . . . . . 71
6.1 Model of a Gliding Robotic Fish with a Whale-type Tail . . . . . . . . . . . . . . . 71

vii

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

72
75
78
78
81
85
89
97
97
98
101

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

115
115
116
121
129

Chapter 8
Three-Dimensional Curve Tracking . . . . . . . . . . . . . . . .
8.1 Three-Dimensional Steady Spiral and Its Differential Geometry Features .
8.2 Influence of Control Inputs on Spiral Trajectory . . . . . . . . . . . . . .
8.3 Two Degree-of-Freedom Control Design . . . . . . . . . . . . . . . . . .
8.3.1 Feedforward Control via Inverse Mapping of Steady Spiral Motion
8.3.2 Linearized Model Around the Steady Spiral Trajectory . . . . . .
8.3.3 H∞ Controller Design . . . . . . . . . . . . . . . . . . . . . . . .
8.4 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

133
133
136
143
145
147
148
151

6.2

6.3
6.4

6.1.1 Dynamic Model in the Sagittal Plane . . . . . . . . . . . . . . .
6.1.2 System Reduction via Singular Perturbation . . . . . . . . . . .
Passivity-based Controller Design . . . . . . . . . . . . . . . . . . . .
6.2.1 Passivity-based Controller for the Approximated Reduced Model
6.2.2 Stability Analysis for the Full Closed-loop System . . . . . . .
6.2.3 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . .
Observer Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.4.1 Filter Design . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.4.2 Open-loop Experiments . . . . . . . . . . . . . . . . . . . . . .
6.4.3 Closed-loop Experiments . . . . . . . . . . . . . . . . . . . . .

Chapter 7
Yaw Stabilization Using Sliding Mode Control
7.1 Problem Statement . . . . . . . . . . . . . . . . . . .
7.2 Sliding Mode Control for Yaw Stabilization . . . . . .
7.3 Simulation . . . . . . . . . . . . . . . . . . . . . . . .
7.4 Experimental Results . . . . . . . . . . . . . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

Chapter 9
Field Test Results for Environmental Monitoring . . . . . . . . . . . . . 159
9.1 Kalamazoo River Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159
9.2 Wintergreen Lake Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
Chapter 10 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . . . . 169
10.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169
10.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170
BIBLIOGRAPHY . . . . . . . . . . . . . .

viii

. . . . . . . . . . . . . . . . . . . 173

LIST OF TABLES

Table 2.1:

Selected components used in “Grace”. . . . . . . . . . . . . . . . . . . . 18

Table 4.1:

Computed steady gliding path under different values of the center of gravity zCG , the movable mass displacement r p1 , and the excess mass m0 , for
the two models shown in Fig. 3.8. . . . . . . . . . . . . . . . . . . . . . 38

Table 4.2:

Computed steady gliding path for the scaled models of the larger wing
prototype. In computation, r p = 5mm is used for the original scale model
(1:1) while the value is scaled linearly with dimension for other models. . 41

Table 4.3:

Comparison of our miniature glider with underwater gliders reported in
the literature. Velocity refers to the terminal velocity. The bottom row
shows the metrics of the glider reported in this dissertation (with the larger
wings). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

Table 5.1:

Parameters of the lab-developed underwater robot used in the steady-state
spiraling equations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

Table 5.2:

Computed spiraling steady states through Newton’s method. . . . . . . . 61

Table 7.1:

Parameters of gliding robotic fish “Grace”. . . . . . . . . . . . . . . . . . 121

Table 9.1:

Locations and sensor readings in Kalamazoo River test. . . . . . . . . . . 160

Table 9.2:

Locations and sensor readings in Wintergreen Lake test. . . . . . . . . . . 163

ix

LIST OF FIGURES

Figure 1.1:

Classic underwater gliders. (a) Seaglider [1]; (b) Spray [2]; (c) Slocum [3].

4

Figure 1.2:

Examples of robotic fish. (a) From Mechatronics Research laboratory at
Massachusetts Institute of Technology [4]; (b) From Smart Microsystems
Lab at Michigan State University. . . . . . . . . . . . . . . . . . . . . . .

5

Figure 1.3:

Schematics of working pattern “gliding in sagittal plane”. . . . . . . . . .

7

Figure 1.4:

Schematics of working pattern “spiraling in 3D space”. . . . . . . . . . .

8

Figure 2.1:

Protoype I, the miniature underwater glider operating in the swimming
pool. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

Figure 2.2:

Gliding robotic fish “Grace” gliding in the large indoor tank in the Smart
Microsystems Lab. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

Figure 2.3:

Gliding Rototic fish “Grace” swimming in Kalamazoo river. . . . . . . . 16

Figure 2.4:

The schematic of internal configuration for gliding robotic fish. . . . . . . 17

Figure 2.5:

Components on the outside of prototype “Grace”. . . . . . . . . . . . . . 18

Figure 2.6:

Components of the tail compartment for “Grace”. . . . . . . . . . . . . . 20

Figure 2.7:

Top cap interface of Grace to the outside when cap covered. . . . . . . . . 20

Figure 2.8:

Top cap interface of Grace to the outside when cap open. . . . . . . . . . 21

Figure 3.1:

The mass distribution of the gliding robotic fish (side view). . . . . . . . . 23

Figure 3.2:

Illustration of the reference frame and hydrodynamic forces. . . . . . . . 23

Figure 3.3:

The meshing used for the water tunnel simulation for miniature underwater glider prototype (generated with Gambit). . . . . . . . . . . . . . . . 28

Figure 3.4:

An example contour of the static pressure with the angle of attack and the
sideslip angle set as zero for miniature underwater glider prototype. . . . 29
x

Figure 3.5:

Data fitting for the drag force coefficient as a function of α . . . . . . . . . 31

Figure 3.6:

Data fitting for the lift force coefficient as a function of α . . . . . . . . . . 32

Figure 3.7:

Data fitting for the pitch moment coefficient as a function of α . . . . . . . 33

Figure 3.8:

Illustration of the glider body with two different wing designs. Area of
wings in (b) is half of that in (a). . . . . . . . . . . . . . . . . . . . . . . 33

Figure 3.9:

An example contour of the static pressure with the angle of attack and the
sideslip angle set as zero for gliding robotic fish prototype. . . . . . . . . 34

Figure 4.1:

The glider cost performance index with respect to model scales. . . . . . 42

Figure 4.2:

The horizontal velocity with respect to model scales. . . . . . . . . . . . 43

Figure 4.3:

Superimposed snapshot of the gliding experiments in an indoor tank. . . . 46

Figure 4.4:

Glide angle with respect to the movable mass displacement, with fixed
net buoyancy of -20 g, for the prototype with larger wings. . . . . . . . . 47

Figure 4.5:

Total glide speed with respect to the movable mass displacement, with
fixed net buoyancy of -20 g, for the prototype with larger wings. . . . . . 48

Figure 4.6:

Glide angle with respect to the net buoyancy, with fixed movable mass
displacement of 0.5 cm, for the prototype with larger wings. . . . . . . . 49

Figure 4.7:

Total glide speed with respect to the net buoyancy, with fixed movable
mass displacement of 0.5 cm, for the prototype with larger wings. . . . . 50

Figure 4.8:

Glide angle of the miniature glider with respect to the movable mass displacement, with fixed net buoyancy of -20 g, when using two different
sets of wings. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

Figure 4.9:

Total glide speed of the miniature glider with respect to the movable mass
displacement, with fixed net buoyancy of -20 g, when using two different
sets of wings. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

xi

Figure 5.1:

Convergence test results for Newton’s method with respect to initial conditions. Color yellow (light) means that convergence to the steady-state
spiraling equilibrium is achievable with the corresponding initial values;
color blue (dark) means that there is no convergent solution or the convergent solution is not at the steady-state spiraling equilibrium. In the
test, the used set of control inputs is r p1 = 5 mm, m0 = 30 g, δ =
45◦ ; and the corresponding equilibrium state values are θ = −42.1281◦,
φ = −34.2830◦, ω3i = 0.4229 rad/s, V = 0.2809 m/s, α = −0.9014◦ ,
β = 4.2414◦. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

Figure 5.2:

Convergence time in seconds for spiraling dynamics with respect to different initial values of states in roll angle φ , pitch angle θ , and translational velocity v1 along the Oxb direction, for the control inputs of
r p1 = 5 mm, m0 = 30 g and δ = 45◦ . The corresponding equilibrium
state values are θ = −42.1281◦ , φ = −34.2830◦ , and v1 = 0.2801 m/s. . 65

Figure 5.3:

Convergence time in seconds for spiraling dynamics with respect to different initial values of states in angular velocities, for the control inputs
of r p1 = 5 mm, m0 = 30 g and δ = 45◦ , displayed in orthogonal slice
planes. The corresponding equilibrium state values are ω1 = 0.2837
rad/s, ω2 = −0.1767 rad/s, and ω3 = 0.2592 rad/s. . . . . . . . . . . . . 66

Figure 5.4:

Snap shots of the robot spiraling in the experiment tank. . . . . . . . . . . 67

Figure 5.5:

Spiraling radius with respect to the tail angle, with fixed movable mass
displacement of 0.5 cm and fixed excess mass of 30 g. . . . . . . . . . . . 69

Figure 5.6:

Spiraling radius with respect to the excess mass , with fixed movable mass
displacement of 0.5 cm and fixed tail angle of 45◦ . . . . . . . . . . . . . . 70

Figure 6.1:

The schematic of a gliding robotic fish with forces and moments defined
in the corresponding coordinate frames (side view). . . . . . . . . . . . . 72

Figure 6.2:

Simulation results on the trajectories of the gliding angle θg for the openloop uδ = 0 and closed-loop (Kc = 2) cases, respectively. . . . . . . . . . 86

Figure 6.3:

The trajectory of control uδ for the closed-loop simulation with different
values for Kc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

Figure 6.4:

The trajectory of gliding angle θg for the closed-loop simulation with
different values for Kc . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

xii

Figure 6.5:

Simulated trajectory of the measured pitch angle which is corrupted with
noise. The measurement noise is modeled as Gaussian random process
with R = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

Figure 6.6:

Simulation results: the trajectories of the gliding angle θg of the real state
and nonlinear observer estimation with measurement noise. . . . . . . . . 96

Figure 6.7:

A snapshot of the open-loop experiment with fixed tail angle using gliding
robotic fish “Grace” in the lab tank. . . . . . . . . . . . . . . . . . . . . . 99

Figure 6.8:

The trajectories of pitch angle of the on-board sensor reading, observer
estimation and computer-based simulation result, in the open-loop experiments using gliding robotic fish “Grace”. (a) δ = 15◦ ; (b) δ = 30◦ ; (c)
δ = 45◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

Figure 6.9:

A snapshot of stabilization experiment using gliding robotic fish “Grace”
in NBRF, University of Maryland. . . . . . . . . . . . . . . . . . . . . . 102

Figure 6.10:

A snapshot of stabilization experiment in the view of the 12 underwater
cameras of the Qualysis motion capture system in NBRF, University of
Maryland. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

Figure 6.11:

The trajectories of pitch angle and gliding angle in the experiment of
gliding stabilization without feedback control. . . . . . . . . . . . . . . . 104

Figure 6.12:

The trajectory of predefined tail angle in the experiment of gliding stabilization without feedback control. . . . . . . . . . . . . . . . . . . . . . . 105

Figure 6.13:

The trajectories of pitch angle and gliding angle in the passivity-based
stabilization experiment with Kc = 3. . . . . . . . . . . . . . . . . . . . . 106

Figure 6.14:

The trajectory of tail angle in the passivity-based stabilization experiment
with Kc = 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

Figure 6.15:

The trajectories of pitch angle and gliding angle in the passivity-based
stabilization experiment with Kc = 1. . . . . . . . . . . . . . . . . . . . . 108

Figure 6.16:

The trajectory of tail angle in the passivity-based stabilization experiment
with Kc = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

xiii

Figure 6.17:

The trajectories in the comparative stabilization experiments with a proportional controller with KP = 2. (a) pitch angle and gliding angle; (b)
tail angle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

Figure 6.18:

The trajectories in the comparative stabilization experiments with a proportional controller with KP = 3. (a) pitch angle and gliding angle; (b)
tail angle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

Figure 6.19:

The trajectories in the comparative stabilization experiments with a PI
controller with KP = 2, KI = 1. (a) pitch angle and gliding angle; (b) tail
angle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

Figure 6.20:

The trajectories in the comparative stabilization experiments with a PI
controller with KP = 2, KI = 3. (a) pitch angle and gliding angle; (b) tail
angle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

Figure 7.1:

Yaw angle trajectory with respect to different controller parameters k1 . . . 122

Figure 7.2:

Tail angle trajectory with respect to different controller parameters k1 . . . 123

Figure 7.3:

Sideslip angle trajectory with respect to different controller parameters k1 . 124

Figure 7.4:

Yaw angle trajectory with respect to different controller parameters k2 . . . 125

Figure 7.5:

Tail angle trajectory with respect to different controller parameters k2 . . . 126

Figure 7.6:

Sideslip angle trajectory with respect to different controller parameters k2 . 127

Figure 7.7:

Trajectory of the gliding robotic fish for controller parameters k1 = 30,
k2 = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128

Figure 7.8:

Snapshots of controlled motion with yaw stabilization using sliding mode
controller, under parameters k1 = 50, k2 = 1. . . . . . . . . . . . . . . . . 130

Figure 7.9:

Yaw angle trajectory using “Grace” . . . . . . . . . . . . . . . . . . . . . 131

Figure 7.10:

Tail angle trajectory using “Grace” . . . . . . . . . . . . . . . . . . . . . 132

Figure 8.1:

Trajectory characteristic parameters at different tail angles with m0 = 30
g and r p = 0.5 cm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138

xiv

Figure 8.2:

Trajectory characteristic parameters at different net buoyancy with δ =
30◦ and r p = 0.5 cm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139

Figure 8.3:

Trajectory characteristic parameters at different displacements of movable mass with δ = 30◦ and m0 = 30 g. . . . . . . . . . . . . . . . . . . 140

Figure 8.4:

The gliding robotic fish “Grace” spiraling in Neutral Buoyancy Research
Facility, University of Maryland. . . . . . . . . . . . . . . . . . . . . . . 141

Figure 8.5:

Snapshots of spiral motion with 12 underwater cameras from different
angles of view using a Qualysis underwater motion capture system. . . . . 142

Figure 8.6:

Illustration of robot’s rigid body and coordinates in spiral motion using a
Qualysis underwater motion capture system. . . . . . . . . . . . . . . . . 142

Figure 8.7:

Model prediction and experimental results in spiral motion at different tail
angles. The other two control inputs are fixed m0 = 30 g, and r p = 0.5
cm. (a) curvature; (b) torsion. . . . . . . . . . . . . . . . . . . . . . . . . 144

Figure 8.8:

The control system diagram with combination of open-loop control and
closed-loop control. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149

Figure 8.9:

Transformed 2-DOF control configuration in H∞ control framework . . . . 149

Figure 8.10:

The simulation results of the geometric parameters when tracking a steady
spiral trajectory. (a) curvature; (b) torsion; (c) velocity. . . . . . . . . . . 154

Figure 8.11:

The simulation results of control inputs on when tracking a steady spiral trajectory. (a) displacement of movable mass; (b) tail angle; (c) net
buoyancy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155

Figure 8.12:

3D trajectory when tracking a steady spiral under the 2-DOF controller. . 156

Figure 8.13:

The simulation results of the geometric parameters when the reference
velocity changes as a sinusoid function with respect to time while curvature and torsion are kept constant. (a) curvature; (b) torsion; (c) velocity. . 157

Figure 8.14:

The simulation results of control inputs on when the reference velocity
changes as a sinusoid function with respect to time while curvature and
torsion are kept constant. (a) displacement of movable mass; (b) tail angle; (c) net buoyancy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158

xv

Figure 9.1:

Three selected sampling locations illustrated in Google Map. . . . . . . . 160

Figure 9.2:

Swimming trajectory of gliding robotic fish in an open water area at testing spot A. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161

Figure 9.3:

Swimming trajectory of gliding robotic fish in the woods at testing spot B. 161

Figure 9.4:

Swimming trajectory of gliding robotic fish under a bridge at testing spot
C. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162

Figure 9.5:

Selected sampling points in the Wintergreen Lake with surface swimming. 163

Figure 9.6:

Depth trajectory of steady glide when sampling water in Wintergreen
Lake, Michigan. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164

Figure 9.7:

Sensor readings of steady glide when sampling water in Wintergreen
Lake, Michigan. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165

Figure 9.8:

Yaw angle trajectory when sampling water in Wintergreen Lake, Michigan.165

Figure 9.9:

Pitch angle trajectory when sampling water in Wintergreen Lake, Michigan.166

Figure 9.10:

Depth trajectory of steady spiral when sampling water in Wintergreen
Lake, Michigan. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166

Figure 9.11:

Sensor readings of steady spiral when sampling water in Wintergreen
Lake, Michigan. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167

Figure 9.12:

Yaw angle trajectory of steady spiral when sampling water in Wintergreen
Lake, Michigan. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167

Figure 9.13:

Pitch angle trajectory of steady spiral when sampling water in Wintergreen Lake, Michigan. . . . . . . . . . . . . . . . . . . . . . . . . . . . 168

xvi

Chapter 1
Introduction
Gliding robotic fish, a new type of underwater robots, is proposed for monitoring aquatic environment in this dissertation. Such a robot combines the advantages of both underwater gliders and
robotic fish, and features long operation duration and high maneuverability.
In this introduction, I will first discuss the existing methods in aquatic environmental monitoring, and conduct a brief literature review on two important underwater robots, the underwater
glider and robotic fish. Afterwards, gliding robotic fish is introduced as a new type of underwater
robots, inspired by the above two well-known classes of robots. Locomotion mechanism and advantages of the gliding robotic fish are briefly discussed. The structure of this dissertation is then
clarified. At last, an overview of the contributions is presented.

1.1 Technology in Aquatic Environmental Monitoring
There is a growing interest in monitoring aquatic environments, due to the emerging problems
of environmental pollution and expanding demand for sustainable development. Such pollutions
involve various types contaminants, including industrial waste, chemicals, and bacteria, etc. The
pollution could happen in different aquatic environments, such as ponds, rivers, lakes, and even
the ocean. Due to ongoing industrialization and expanding exploration of aquatic resources, water
pollution problems have become increasingly frequent and severe, which has drawn a global attention [5–7]. In particular, the massive 2010 oil spill in the Gulf of Mexico has brought world-wide
1

attention and attracted intensive research into this urgent issue. The cleanup and recovery will take
many years and millions of dollars. Throughout this process, monitoring the water quality and detecting the remaining contaminants in the water are an essential task [8, 9]. Similar contamination
incidents include the 2010 oil spill in Kalamazoo River, Michigan, [10], 2010 Xingang Port oil
spill in China Yellow Sea [11], 2011 Nigeria oil spill [12], and 2012 Arthur Kill storage tank spill
in New Jersey [13].
Throughout the development in technology for aquatic environmental monitoring, a variety of
different water sampling methods have been studied and employed. Manual sampling, via boat/ship or with handheld devices, is still a common practice in aquatic environmental monitoring.
This approach is labor-intensive and has difficulty capturing dynamic and spatially distributed
phenomena of interest. An alternative is in-situ sensing with fixed or buoyed/moored sensors [14],
including vertical profilers that can move up and down along a water column [15–19]. However,
since buoyed sensors have little or no freedom to move around, it could take a prohibitive number
of them to capture spatially inhomogeneous information. The past decade or so has seen great
progress in the use of robotic technology in aquatic environmental sensing [20–36]. Predominant
examples of these technologies include remotely operated vehicles (ROVs) [20, 24, 27, 33], autonomous surface vehicles (ASVs) [22, 23, 25, 31, 35], propeller-powered autonomous underwater
vehicles (AUVs) [24, 30, 34, 36–39], and underwater gliders [1–3]. ROVs typically have limited
spatial access and autonomy due to their tethered nature, while the sampling space of ASVs is
limited to the two-dimensional (2D) water surface. AUVs, on the other hand, can operate freely
and autonomously in the 3D water body, but their high price tags (upward of $150K per vehicle)
presents a huge barrier to their deployment in large numbers for high-resolution spatiotemporal
coverage. Besides, there has been a growing interest in bio-inspired underwater robots [40–46],
such as robotic fish, which holds great potential in wide application in water monitoring. The
2

shortcoming of this kind of robots is the limited operational duration and thus area coverage, due
to the energy constraints and the requirement of constant actuation for propulsion..

1.2 Gliding Robotic Fish
1.2.1 Design Concept
The design concept of a gliding robotic fish comes from energy-efficient underwater gliders and
highly-maneuverable robotic fish.
Underwater seagliders are known for their great energy-efficiency and long-duration operation
in oceanographic applications [47]. An underwater glider utilizes its buoyancy and gravity to
enable motion without any additional propulsion, and adjusts its center of gravity to achieve certain
attitude, which results in glide and thus horizontal travel. Since energy is needed only for buoyancy
and center-of-gravity adjustment when switching the glide profile, underwater gliders are very
energy-efficient, as proven by the great success of the Seaglider [1], Spray [2] and Slocum [3]
(Fig. 1.1). The maneuverability of underwater gliders, however, is quite limited. The large size
(1– 2 m long), heavy weight (50 kg and above), and high cost [47] of these vehicles also impede
their adoption in the application of networked sensing and in versatile environments like ponds and
inland lakes.
On the other hand, over the past two decades, there has been significant interest in developing
robots that propel and maneuver themselves like real fish do. Often called robotic fish (Fig. 1.2),
they accomplish swimming by deforming the body and fin-like appendages [48–68]. In many
designs, a fish-like flapping tail is used to provide propulsion force and a biased tail angle is applied
to realize turning. Due to the similarity to the real fish and the fact that the rotation of the tail fin
is usually enabled by a motor that is easy to control and fast in response, robotic fish typically
3

(a)

(b)

(c)
Figure 1.1: Classic underwater gliders. (a) Seaglider [1]; (b) Spray [2]; (c) Slocum [3].

4

(a)

(b)
Figure 1.2: Examples of robotic fish. (a) From Mechatronics Research laboratory at Massachusetts
Institute of Technology [4]; (b) From Smart Microsystems Lab at Michigan State University.

5

have high maneuverability (e.g., small turning radius); However, as the forward propelling force is
generated from the flapping motion of the fins, such robots require constant actuation for swimming
and cannot work for extended periods of time without battery recharge.
In this dissertation, inspired by the design and merits of both underwater gliders and robotic
fish, a new type of underwater robots, gliding robotic fish, is proposed. Gliding robotic fish combines both mechanisms of gliding and swimming and features energy efficiency and high maneuverability at the same time. Such a robot would realize most of its locomotion through gliding like
underwater gliders, by utilizing its buoyancy and gravity to enable motion without any additional
propulsion, and adjusting its center of gravity to achieve a certain attitude. It would use actively
controlled fins to achieve high maneuverability, during turning and orientation maintenance. Of
course, fins can also provide additional propulsive power during locomotion, if needed.
The dimension of the gliding robotic fish is supposed to be smaller than traditional underwater
gliders, and one use of such a robot is to provide a mobile sensing platform in relatively shallow
waters, such as lakes, rivers, and even ponds, where larger gliders are not quite suitable due to
their large size and high cost. The small size and low cost of gliding robotic fish also facilitate the
research of networked sensing and operative control.

1.2.2 Motion and Control
For gliding robotic fish, there could be various interesting motions generated by integrating the
gliding and swimming mechanisms. In this dissertation, two main steady motion profiles are discussed as the regular working patterns for sampling water. One is steady-state gliding in the sagittal
plane, and the other is steady-state spiraling in the three-dimensional (3D) space.
The steady-state gliding is also the common operating mode for traditional underwater gliders, used for sampling the water field 1.3. In the zig-zag trajectory, the gliding robotic fish only
6

Figure 1.3: Schematics of working pattern “gliding in sagittal plane”.

7

Figure 1.4: Schematics of working pattern “spiraling in 3D space”.

8

consumes energy during the transition between descending and ascending. In the phase of steady
glide, balances of forces and moments are achieved. If we interpret the motion intuitively, the
net buoyancy of the robot enables a propelling force to counteract the hydrodynamic resistance.
Zero energy consumption in the steady gliding period makes the gliding robotic fish highly energyefficient.
The other motion, steady spiral in the three-dimentional space, is proposed as a novel method
of sampling a water column for gliding robotic fish 1.4. A water column is a conceptual narrow
volume (like a narrow cylinder) of water stretching vertically from the surface to the bottom. Water
column sampling is a routine surveying method in environmental studies to evaluate the stratification or mixing of water layers [69]. This motion is achieved by incorporating the steady gliding
with a deflected fish-like tail. The non-zero angled tail fin will introduce a turning moment to the
steady gliding robot and lead to a 3D spiraling motion. With the actively-controlled tail, the gliding
robotic fish is capable of spiral motions with tight turning radius. If needed, a gradually changing
tail angle will form spiral-in or spiral-out 3D trajectories.
For most of its operation time, a gliding robotic fish holds a preset steady gliding path without
any energy consumption. The gliding angle is calculated before deployment so that the robot will
follow a designed trajectory, which either minimizes the energy cost, or maximizes the field mapping capability [70]. However, a gliding robot fish is subject to many non-negligible uncertainties
from the aquatic environment (e.g. current disturbance), which results in additional energy cost because counteracting the deviation from the preset course requires re-calibration and control effort
to keep the robot on or bring it back to the designed path [71] [29] [72]. The stability property of
the steady gliding path and fast convergence to the path, are very important for the gliding robotic
fish to reduce the energy expenditure on path correction. The stabilization involves both sagittal
and lateral dynamics, corresponding to sagittal and lateral disturbances.
9

For the gliding angle stabilization in the sagittal plane, there are several related results in literature for underwater gliders. Most of them use the net buoyancy and internal mass displacement
as control variables. For instance, in [73] and [74] the authors presented an LQR (linear quadratic
regulator) controller and a PID controller, respectively, but both based on a linearized model. Although in the LQR method energy is used as a cost function, the approach does not consider the
additional cost for the course correction and guidance. Both LQR and PID controller can not
fully address the nonlinear gliding characteristics. Nonlinear controllers involving torque control
and buoyancy control are proposed by Bhatta [75], but those methods ignore the dynamics of the
control surface and also require full state feedback for the controller implementation, which will
increase the complexity of the software and hardware. A Lyapunov-based control design is reported in [76]; however, for the elevator control, it only deals with a fixed-value control input to
achieve a certain equilibrium gliding path. In this dissertation, a novel, passivity-based nonlinear
controller is proposed for the sagittal-plane stabilization problem using only a whale-like tail fin of
the gliding robotic fish. The singular perturbation results from [76] are utilized to reduce the fullorder system to a slow-mode second-order system. A passivity-based controller is designed based
on the approximated reduced system, and the controller is applied back into the original full-order
system. Through checking the Jacobian matrix, the local stability of the full-order closed-loop system is established given sufficient time-scale separation. Furthermore, a nonlinear observer is also
proposed to estimate velocity-related system states, which are used in controller implementation.
Both open-loop and closed-loop experiments are carried out to illustrate the effectiveness of the
designed controller and observer.
For the lateral motion, attitude stabilization is of great importance to underwater robots, which
are subjected to various environmental disturbances. There has been extensive work on pitch
motion stabilization in the longitudinal plane for underwater gliders as discussed in the previous
10

paragraph, where typically the absence of lateral motion is assumed [73, 76–78]. There has also
been some limited research into three-dimensional gliding involving lateral motion, most of which
focuses on the steady-state turning or spiraling [79–81]. However, little work has been reported on
the yaw angle stabilization where ambient flow disturbances could easily push the robot away from
its desired heading orientation. In this dissertation, the relative degree of the system dynamics is
first identified and error dynamics is then derived based on the system normal form. A sliding mode
controller is designed for the error dynamics to achieve attitude stabilization. Both simulation and
experiments are conducted to evaluate the effectiveness of the proposed control scheme.
Three-dimensional curve tracking is also crucial for underwater robots, including in particular, gliding robotic fish. For example, it is critical in sampling water columns, seeking pollutant
sources, and mapping the whole aquatic environment. Precise tracking control is very challenging
due to dynamic nonlinearity and strong coupling among multiple control inputs. As stated previously, there has been some limited literature that covers three-dimensional gliding involving lateral
motion, most of which focuses on the steady-state turning or spiraling [79–81]. However, it is still
an unexplored area of three-dimensional curve tracking control for buoyancy-driven underwater
robots. In this dissertation, I propose a novel two degree-of-freedom (DOF) control strategy for
gliding robotic fish to track three-dimensional curves based on the differential geometric features
of steady spirals. The control strategy includes a feedforward controller that is designed through
inverse mapping of steady spiral motion, and a feedback controller designed using the robust H∞
framework based on local, linearized dynamics. The effectiveness of the proposed 2-DOF control
scheme is demonstrated in simulation, where PI control and open-loop inverse mapping control
are also conducted for the purpose of comparison.
There are a number of challenges in aquatic environmental monitoring, especially in field experiments using autonomous underwater robots. Gliding robotic fish is proposed to carry out the
11

water sampling work with great energy efficiency, high maneuverability and adaptability to versatile environments. As the final test of the design concept and robot development, the gliding robotic
fish prototype is taken to Kalamazoo River, Michigan, to detect oil spill, and Wintergreen Lake,
Michigan, to sample harmful algal blooms. The external swappable environmental sensor reads
the field data of interest and the on-board micro-controller sends back the information wirelessly
to the base station. The GPS-based positioning system facilitates the water sampling task with
precise location information, where the motions of swimming, gliding and spiraling are tested and
integrated. The field test results are presented to demonstrate the functionality and usefulness of
the proposed gliding robotic fish as a novel platform for solving the real-world problem of aquatic
environmental monitoring.

1.3 Overview of Contributions
The contributions of this research reside mainly on dynamics analysis and tail-enabled control of
a new type of underwater robot, the gliding robotic fish. The details are as follows.
First, the design idea of combining buoyancy-driven propulsion and tail-enabled maneuverability is novel. The advantages of long operation duration and high maneuverability from underwater
gliders and robotic fish, make the gliding robotic fish a suitable underwater platform in shallow
water environmental sensing.
Second, the discussion about the two water-sampling working patterns, the sagittal-plane glide
and three-dimensional spiral, provides insights of the energy-efficient feature of gliding robotic
fish. Particularly, the tail-enabled steady spiral, is a novel motion proposed for sampling water
columns, and potential path planning tasks. The turning radius could be as tight as 0.5 m, providing
much higher sampling resolution compared to the 30-40 m radius from traditional underwater

12

gliders.
Third, the steady glide stabilization is usually realized through buoyancy control or massdistribution control in literature. In this dissertation an actively-control tail fin is adopted as the
control input to stabilize the steady gliding motion. For the first time in literature of the sagittalplane motion of buoyancy-driven underwater robots, a passivity-based feedback controller is designed to obtain a fast convergence speed with partial state feedback, which is convenient for
implementation. For the heading maintenance in lateral dynamics, a sliding mode controller is
proposed with only a requirement on the yaw angle feedback information. Moreover, both simulation and experimental results are presented to provide insight into the stabilization problem for
buoyancy-driven underwater robots.
Fourth, three-dimensional curve tracking for buoyancy-driven underwater robots is very challenging and little relevant work could be found. In this dissertation, A novel two-degree-offreedom control strategy, which consists of a feedforward controller designed through inverse mapping of steady spiral motion, and a feedback controller designed using the robust H∞ framework
based on local, linearized dynamics, is proposed for the curve tracking problem using the differential geometric features of steady spirals. The study of the geometric features of the spiral motion,
such as torsion and curvature, also shows a promising motion planning and navigation method.
Finally, the successful development of gliding robotic fish prototype and the field tests in the
Kalamazoo River and the Wintergreen Lake, proves the concept of the design, and demonstrates
the functions of the robot as a platform for aquatic environmental monitoring. This leads to more
research directions, such as adaptive sampling, networked control and cooperative control with
multiple agents, and it opens the door to further collaboration with biologists and sociologists to
solve bigger real-world water-pollution problems.

13

Chapter 2
Implementation of Gliding Robotic Fish
During the research of the gliding robotic fish, two prototypes have been developed. The first
one is a miniature underwater glider, reported in [82, 83], built to study the gliding component of
gliding robotic fish and help to evaluate the steady glide model in the sagittal plane (Fig. 2.1). The
second one is a fully functioning gliding robotic fish prototype, named “Grace” (Fig. 2.2) [84]. In
this section, I will focus on the implementation of the gliding robotic fish “Grace”, which naturally
explains the case for the miniature underwater glider.

2.1 Actuation System
Integrating an actively-controlled tail into the miniature underwater glider, the research team has
developed a fully functioning gliding robotic fish, named “Grace”. The robot has three actuation
systems for locomotion, including the buoyancy system, the mass distribution system, and the
actively-controlled tail fin system.
In the buoyancy system, water is pumped in and out of the robot’s body to change the net
buoyancy. When the robot is heavier than the water it displaces (negatively buoyant), the robot
will descend (Fig. 2.2); and when it is lighter than the water it displaces (positively buoyant), the
robot will ascend. The pumping system of “Grace” is enabled by a linear actuator with integrated
feedback, which allows the precise control of water volume despite the pressure differences at
different depths, while a DC pump is used for the miniature underwater glider prototype.
14

Figure 2.1: Protoype I, the miniature underwater glider operating in the swimming pool.

Figure 2.2: Gliding robotic fish “Grace” gliding in the large indoor tank in the Smart Microsystems
Lab.

15

PDF created with pdfFactory Pro trial version www.pdffactory.com

Figure 2.3: Gliding Rototic fish “Grace” swimming in Kalamazoo river.
For the mass distribution system, a linear actuator is used to push a mass (battery pack) back
and forth along a guiding rail to change the center of the mass, for the purpose of manipulating the
pitch angle, in both miniature underwater glider and “Grace”.
The fish-like tail fin system in “Grace” is driven by a servo motor through a chain transmission. In three-dimensional gliding, a deflected tail can be used to control the turning motion and
heading orientation. Like a real fish, the robot can also flap the tail to realize the swimming motion
(Fig. 2.3).

2.2 Gliding Robotic Fish Components
A schematic of selected components for the gliding robotic fish is shown in Fig. 2.4. In the figure,
the components of the three actuation systems can be identified. Besides, there are two physically
separated printed circuit boards (PCBs). One is the ctrl PCB containing the micro-controller and
navigational sensors, such as gyroscopes, accelerometers and a digital compass; the other is the
16

Figure 2.4: The schematic of internal configuration for gliding robotic fish.
driver PCB containing regulators and driver components for the actuators, including the linear
actuators and the servo motor. A pressure sensor is used to measure the depth of the current
location of the robot, with one port connecting to the ambient water.
There are some other components equipped on the shell of the gliding robotic fish, including
the wireless communication antenna, environmental sensors and the GPS receiver, as shown in
Fig. 2.5. “Grace” is equipped with a crude oil sensor and a temperature sensor, and has been
tested in the Kalamazoo River, Michigan, to sample the oil concentration near the site of a 2010 oil
spill (Fig. 2.3). The sensor can be easily swapped to measure other environmental processes, such
as chlorophyll, harmful algae, turbidity, rhodamine. The GPS unit is used to measure the global
position and provide the universal time, when the gliding robotic fish surfaces. Table 2.1 lists the
details of the used components mentioned above.

17

Figure 2.5: Components on the outside of prototype “Grace”.

Table 2.1: Selected components used in “Grace”.
Component name
1
Micro-controller
2
Battery
3
Linear actuator
4 Pump 1 (miniature underwater glider)
5
Pump 2 (“Grace”)
6
Servo motor
7
Pressure sensor
8
GPS
9
Gyro
10
Accelerometer+Compass
11
Wireless module
12
Wireless antenna
13
Crude oil sensor
14
Temperature sensor

Component model
Microchip dsPIC6014A
Batteryspace 18.5V Polymer-Li-Ion battery pack
Firgelli L16-140-63-12-P
Flight Works Model 300C
Servocity 180 lbs thrust linear actuator
Hitec Servo HS-7980TH
Honeywell 40PC100G2A
Garmin GPS 18x LVC
ST LPY503AL
ST LSM303DLH
XBee Pro 900 XSC RPSMA
900MHz Duck Antenna RP-SMA
Turner Designs Cyclops-7 Crude Oil Sensor
TMP36

18

2.3 Mechanical Design
The mechanical structure of the robot is designed considering the requirements of compactness,
low cost, and energy efficiency, and the limitation of having to house all necessary electrical components. The outer shell is designed to have a fish-like shape with a streamline profile for the
purpose of energy efficiency. The pair of wings is designed to provide enough hydrodynamic
forces and moments for a satisfactory glide. A guide rail is designed to help push the battery pack
using the linear actuator to change the mass distribution.
SolidWorks is used to draw schematics of most of the mechanical components. A computer
numerical control (CNC) machine is used to manufacture parts of the system, such as the mold of
the shell. 3D printer is used to print other parts of the system, such as the tail compartment.
The robot shell is made of carbon-fiber, which is strong enough to sustain the underwater
pressure. The tail compartment is made of composite polymer via 3D printing, which enables
rapid prototyping. The guide rail is made of stainless steel, which is rigid enough to provide a
straight traveling path for the battery pack.
Some physical parameters of “Grace” are as the follows. Length is 65 cm (body) / 90 cm
(total); width is 15 cm (body) / 75 cm (with wings); and height is 18 cm (body) / 34 cm (including
tail). Here, the term “total” includes the body, the tail, the antenna and the environmental sensor.
The wings are in a trapezoidal shape with a wingspan of 30 cm (one side) and an aspect ratio of
1.45. The weight is 9 kg in total. The tail with the servo compartment weighs 0.8 kg itself.
The tail is special in the design in that it is the only moving part seen from the outside, so it
must be designed properly to be waterproof, robust, and reliable with respect to different flapping
amplitude and frequency. The tail flapping motion is transformed from the rotation of a servo by a
chain system, in which two identical gears are used as shown in Fig. 2.6.

19

Figure 2.6: Components of the tail compartment for “Grace”.

Figure 2.7: Top cap interface of Grace to the outside when cap covered.

20

Figure 2.8: Top cap interface of Grace to the outside when cap open.
The robot must be waterproofed for all electrical/electronic components. In “Grace”, all electrical/electronic components are placed in the main body. In particular, the system switch, charging
port, and programming port are located under the top cap. A plunger-like cap covered by rubber is
used to seal the top interface, with two screw bolts to ensure the sealing, as shown in Figs. 2.7 and
2.8.
In addition, the tail is detachable from the main body, enabling us to replace the tail if something goes wrong. The detachable tail enforces the reliability of the sealing of the main body with
double O-rings. Furthermore, the tail can be configured to flap sideways (like a shark) or up and
down (like a whale).

21

Chapter 3
Dynamic Model of Gliding Robotic Fish

3.1 Full Dynamic Model
A gliding robotic fish is a hybrid of a miniature underwater glider and a robotic fish, and its modeling will need to incorporate the effects of both. The tail of gliding robotic fish is used to control
the lateral motion for underwater gliding, and used to propel the robot for swimming. We treat
it as a control surface and a source for external forces and moments. The robot is modeled as a
rigid-body system, including an internal movable mass for center-of-mass control and a water tank
for buoyancy adjustment [73, 82, 83]. On the other hand, the deflected tail provides external thrust
force and side force as well as the yaw moment.
Fig. 3.1 shows the mass distribution within the robot. The stationary body mass ms (excluding
the movable mass) has three components: hull mass mh (assumed to be uniformly distributed),
point mass mw accounting for nonuniform hull mass distribution with displacement r w with respect to the geometry center (GC), and ballast mass mb (water in the tank) at the GC, which is a
reasonable simplification since the effect on the center of gravity caused by the water in the tank
is negligible compared with the effect from the movable mass. The movable mass m,
¯ which is
located at r p with respect to the GC, provides a moment to the robot. The motion of the movable
mass is restricted to the longitudinal axis. The robot hull displaces a volume of fluid of mass m .
Let m0 = ms + m¯ − m represents the excess mass (negative net buoyancy). The robot will sink if
m0 > 0 and ascend if m0 < 0.
22

yb

O

xb

rp
zb

mb
mh

m
O

mw

rw

Figure 3.1: The mass distribution of the gliding robotic fish (side view).

A
z

x
y
D

-FS
O

L

xb
xv

yv

zv

į
zb

//zb

yb

Figure 3.2: Illustration of the reference frame and hydrodynamic forces.

23

The relevant coordinate reference frames are defined following the standard convention. The
body-fixed reference frame, denoted as Oxb yb zb and shown in Fig. 3.2, has its origin O at the
geometry center, so the origin will be the point of application for the buoyancy force. The Oxb
axis is along the body’s longitudinal axis pointing to the head; the Ozb axis is perpendicular to Oxb
axis in the sagittal plane of the robot pointing downward, and Oyb axis is automatically formed by
the right-hand orthonormal principle. In the inertial frame Axyz, Az axis is along gravity direction,
and Ax/Ay are defined in the horizontal plane, while the origin A is a fixed point in space.
As commonly used in the literature, R represents the rotation matrix from the body-fixed reference frame to the inertial frame. R is parameterized by three Euler angles: the roll angle φ , the
pitch angle θ , and the yaw angle ψ . Here


 cθ cφ sφ sθ cψ − cφ sψ cφ sθ cψ + sφ sψ


R =  cθ sψ cφ cψ + sφ sθ sψ −sφ cψ + cφ sθ sψ


−sθ
sφ cθ
cφ cθ
T

where s(·) is short for sin(·) and c for cos(·). Let v b =

v1 v2 v3









and ω b =

(3.1)

T

ω1 ω2 ω3

represent the translational velocity and angular velocity, respectively, expressed in the body-fixed
frame. The subscript b indicates that the vector is expressed in the body-fixed frame, and this
notation is applied throughout this dissertation.
We assume that the tail fin is rigid and pivots at the junction between the body and the tail about
the Ozb axis. The tail induces an external thrust force F t on the robot when it flaps. There are also
other hydrodynamic forces and moments generated because of the relative movement between the
tail and the surrounding water, like the side force and the yaw moment.
By extending the previous modeling work for underwater gliders [81, 82], we obtain the dy-

24

namic model for a gliding robotic fish with an actively-controlled tail fin as the following,

b˙ i = R v b

(3.2)

R˙ = R ωˆ b

(3.3)

RT k + F ext
v˙ b = M −1 M v b × ω b + m0 gR

(3.4)

ω˙ b = J −1 −J˙ ω b + J ω b × ω b + M v b × v b + T ext
+mw grr w × R T k + mgr
¯ r p × RT k

(3.5)

Here M = (ms + m)I
¯ I + M f = diag{m1 , m2 , m3 }, where I is the 3 × 3 identity matrix, and M f
is the added-mass matrix, which can be calculated via strip theory [85]. J = diag{J1, J2 , J3 } is
the sum of the inertia matrix due to the stationary mass distribution and the added inertia matrix
in water. In addition, k is the unit vector along the Az direction in the inertial frame, r w is a
constant vector, and r p is the controllable movable mass position vector, which has one degree of
T

freedom in the Oxb direction, r p =

r p1 0 0

T

. bi =

x y z

is the position vector of

the robot in the inertial reference frame, ωˆ b is the skew-symmetric matrix corresponding to ω b .
F ext stands for all external forces: the external thrust force F t induced by tail flapping, and the
external hydrodynamic forces (lift force, drag force and side force) acting on the gliding robotic
fish body, expressed in the body-fixed frame. Finally, T ext is the total hydrodynamic moment
caused by F ext .

25

3.2 Hydrodynamic Model
3.2.1 Hydrodynamic Equations
In order to model the hydrodynamics, we first introduce the velocity reference frame Oxv yv zv .
Oxv axis is along the direction of the velocity, and Ozv lies in the sagittal plane perpendicular to
Oxv . Rotation matrix R bv represents the rotation operation from the velocity reference frame to the
body-fixed frame:


 cα cβ −cα sβ −sα


Rbv =  sβ
cβ
0


sα cβ −sα sβ cα






,



(3.6)

where the angle of attack α = arctan (v3 /v1 ) and the sideslip angle β = arcsin(v2 / v b )
The hydrodynamic forces include the lift force L, the drag force D, and the side force FS ; the
hydrodynamic moments include the roll moment M1 , the pitch moment M2 , and the yaw moment
M3 . All of those forces and moments are defined in the velocity frame [86]. And if we further
assume that the tail is flapping relatively slowly and smoothly, usually true for the yaw control
motion during steady glide, the propelling force from the tail will be negligible compared to the
buoyancy-induced propelling force, which means F t = 0 . Then we will have the following relationship:
T

F ext = R bv

(3.7)

−D FS −L
T

T ext = R bv

M1 M2 M3

(3.8)

The hydrodynamic forces and moments are dependent on the angle of attack α , the sideslip

26

angle β , the velocity magnitude V [87–90], and the tail angle δ :

α α 2 +C δ δ 2 )
D = 1/2ρ V 2 S(CD0 +CD
D
β
S

(3.9)

FS = 1/2ρ V 2 S(CF β +CFδ δ )

(3.10)

L = 1/2ρ V 2 S(CL0 +CLα α )

(3.11)

S

β

M1 = 1/2ρ V 2 S(CM β + Kq1 ω1 )

(3.12)

α α +K ω )
M2 = 1/2ρ V 2 S(CM0 +CM
q2 2

(3.13)

R

P

β

δ δ)
M3 = 1/2ρ V 2 S(CM β + Kq3 ω3 +CM
Y

Y

(3.14)

where ρ is the density of water and S is the characteristic area of the gliding robotic fish. The tail
angle δ is defined as the angle between the longitudinal axis Oxb and the center line of the tail
projected into the Oxb yb plane, with Ozb axis as the positive direction. Kq1 , Kq2 , Kq3 are rotation
damping coefficients. All other constants with ‘C’ in their notations are hydrodynamic coefficients,
whose values can be evaluated through CFD-based water tunnel simulation [82].

3.2.2 CFD-based Evaluation of Hydrodynamic Coefficients
In this dissertation we use CFD simulation to obtain the hydrodynamic coefficients for any given
gliding robotic fish body geometry. Experimental methods like towing experiments can be further
used to verify the CFD results [91], [92].
We first look into the influence of the angle of attack on the hydrodynamic forces and moments.
We simulate the steady gliding motion when the sideslip angle and the tail angle (no activelycontrolled tail for miniature underwater) are both zero in order to eliminate their influence. In the
simulation the gliding robotic fish model is created in SolidWorks 2009. We use Gambit 2.3.16 to

27

Figure 3.3: The meshing used for the water tunnel simulation for miniature underwater glider
prototype (generated with Gambit).

28

Figure 3.4: An example contour of the static pressure with the angle of attack and the sideslip
angle set as zero for miniature underwater glider prototype.

29

create a mesh file (Fig. 3.3), which contains the shape and size information of the gliding robotic
fish, and then use Fluent 6.2.16 to simulate the flow and pressure distribution around the robot
body (Fig. 3.4), which is placed in a water channel. With a virtual water tunnel as the simulation
workspace, the boundary conditions are set to different inlet velocities, and under such different
boundary conditions, CFD simulations are run with different angles of attack. With given values
for the characteristic area and length, we form a table for the convergent coefficients obtained from
CFD simulation for lift force, drag force and pitch moment. From the convergent results, lift, drag
and pitch moment functions are approximated with polynomials in α by curve fitting. When tail
angle δ and sideslip angle β are fixed at zero, drag coefficient CD , lift coefficient CL , and pitch
moment coefficient CM2 can be expressed as
α α 2,
CD = CD0 +CD

CL = CL0 +CLα α ,
α α,
CM2 = CM0 +CM
P

Data fitting is conducted to compute the drag, lift, and pitch moment coefficients, as a function
of the angle of attack α . For example, Fig. 3.5, Fig. 3.6 and Fig. 3.7 show the fitted function for the
drag, lift and pitch moment coefficients, respectively, for the miniature underwater glider model.
The constants in the lift, drag, and pitch moment functions are estimated to be:
α = 17.5948,
CD0 = 0.45275, CD

CL0 = 0.074606, CLα = 19.5777,
α = 0.5665.
CM0 = 0.0075719, CM
P

As a comparative trial, another set of wings is used with the same wingspan but doubled aspect
30

1.2

Drag coefficient C

D

1

0.8

0.6

0.4

0.2

0
−10

−5

0

Angle of attack α(°)

5

Figure 3.5: Data fitting for the drag force coefficient as a function of α .

31

10

4
3

Lift coefficient C

L

2
1
0
−1
−2
−3
−4
−10

−5

0

Angle of attack α(°)

5

Figure 3.6: Data fitting for the lift force coefficient as a function of α .

32

10

Pitch moment coefficient C

M

0.15

0.1

0.05

0

−0.05

−0.1
−10

−5

0

Angle of attack α(°)

5

10

Figure 3.7: Data fitting for the pitch moment coefficient as a function of α .

(a) larger wings

(b) smaller wings

Figure 3.8: Illustration of the glider body with two different wing designs. Area of wings in (b) is
half of that in (a).
33

Figure 3.9: An example contour of the static pressure with the angle of attack and the sideslip
angle set as zero for gliding robotic fish prototype.

ratio (i.e., chord length is half), while the body is left unchanged. Fig. 3.8 illustrates the two models
created in SolidWorks that have different sets of wings. The hydrodynamic coefficients relevant to
the angle of attack for the one with smaller wings are:
α = 10.298,
CD0 = 0.44724, CD

CL0 = 0.054273, CLα = 11.5545,
α = 0.2903.
CM0 = 0.0062683, CM
P

With similar CFD water-tunnel simulation, we can obtain the hydrodynamic coefficients re-

34

garding the sideslip angle when we set the angle of attack and the tail angle to zero. And the
hydrodynamic coefficients involving the tail angle can be obtained similarly by setting the other
two hydrodynamic angles to zero (3.9). Here we want to point out that in this dissertation we
ignore the coupling effect of those three angles (the sideslip angle, the angle of attack and the tail
angle) on the dynamics of the glider.

35

Chapter 4
Steady Gliding in the Sagittal Plane

4.1 Reduced Model in the Sagittal Plane
Steady-state gliding in the sagittal plane is one of most important working patterns for the gliding
robotic fish. When the robot motion is restricted to the sagittal plane, we have


 cos θ 0


R=
0
1


− sin θ 0


 0 




ω b =  ω2  ,




0







sin θ 
 x 
 

 

, b i =  0 ,
0 
 

 

cos θ
z




 r p1 






r p =  0 , r w = 






0





 v1 




v b =  0 ,




v3

0 


, δ = 0.
0 


rw3

Here we assume the point mass mw is just below the center of geometry O by rw3 as such bottomheavy design is desirable for stability concern and also achievable with manufacture. Plugging the

36

hydrodynamic forces and moment into the robot dynamics equations, we get the following model:

v˙1 = (m1 + m)
¯ −1 − (m3 + m)
¯ v3 ω2 − m0 g sin θ + L sin α − D cos α

(4.1)

v˙3 = (m3 + m)
¯ −1 (m1 + m)
¯ v1 ω2 + m0 g cos θ − L cos α − sin α

(4.2)

ω˙ 2 = J2−1 M2 + (m3 − m1 )v1 v3 − mw grw3 sin θ − mgr
¯ p1 cos θ

(4.3)

x˙ = v1 cos θ + v3 sin θ

(4.4)

z˙ = −v1 sin θ + v3 cos θ

(4.5)

θ˙ = ω2

(4.6)

During steady glide, the angular velocity is zero, while the velocity stays unchanged. The
control r p1 and m0 are constant, which means that the position of the movable mass is fixed with
respect to the origin O and the pumping rate is zero. So the steady motion can be described by:

0 = −m0 g sin θ + L sin α − D cos α

(4.7)

0 = m0 g cos θ − L cos α − D sin α

(4.8)

0 = M2 + (m3 − m1 )v1 v3 − mw grw3 sin θ − mgr
¯ p1 cos θ

(4.9)

The solution to the above equation gives us the steady gliding path.

4.2 Computation of Steady Gliding Path in the Sagittal Plane
With the hydrodynamic parameters obtained from CFD simulation, let us take a look at the solution
of the steady gliding equations (4.7)–(4.9). These equations are highly nonlinear due to the terms
involving trigonometric functions and inverse trigonometric functions in the state. When the angle

37

Table 4.1: Computed steady gliding path under different values of the center of gravity zCG , the
movable mass displacement r p1 , and the excess mass m0 , for the two models shown in Fig. 3.8.
(a) larger wings
zCG (cm)
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.2
0.2
0.2

r p1 (cm)
0.3
0.5
0.7
0.3
0.5
0.7
0.3
0.5
0.7
0.3
0.5
0.7

m0 (g)
10
10
10
30
30
30
50
50
50
10
10
10

V, α , θg (m/s,◦,◦ )
(0.1129 3.0470 -29.5522)
(0.1366 1.6543 -43.0404)
(0.1485 1.0936 -52.7389)
(0.1766 3.9483 -25.0827)
(0.2245 2.0106 -38.4395)
(0.2495 1.3300 -48.7594)
(0.2245 4.0967 -24.5276)
(0.2846 2.1371 -37.0375)
(0.3174 1.3980 -47.0516)
(0.0856 5.8988 -20.2069)
(0.1084 3.3827 -27.6331)
(0.1240 2.3211 -35.1835)

(b) smaller wings
zCG (cm)
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.2
0.2
0.2

r p1 (cm)
0.3
0.5
0.7
0.3
0.5
0.7
0.3
0.5
0.7
0.3
0.5
0.7

m0 (g)
10
10
10
30
30
30
50
50
50
10
10
10

38

V, α , θg (m/s,◦,◦ )
(0.1221 4.0658 -37.0187)
(0.1396 2.4940 -48.3662)
(0.1486 1.7575 -56.4820)
(0.2260 3.2732 -41.9002)
(0.2477 2.2108 -51.2303)
(0.2598 1.6385 -58.0061)
(0.3105 2.5525 -47.8110)
(0.3290 1.8747 -55.0414)
(0.3401 1.4598 -60.4123)
(0.0949 7.7979 -26.1314)
(0.1136 5.0110 -32.7957)
(0.1265 3.6361 -39.4833)

of attack is small enough, we can use the approximation sin α ≈ α and cos α ≈ 1, and derive an
approximate analytical solution for the desired control r p and m0 in order to achieve some given
steady states [74]. However, here we are interested in the calculation of the steady gliding states
themselves under a fixed control. Unfortunately, there are no feasible analytical solutions for this
problem. So with Matlab command solve(), we numerically solve equations (4.7)–(4.9) to get the
velocity v , pitch angle θ , and glide angle θg for a given movable mass displacement r p and net
buoyancy m0 , under different conditions for r w , the location of nonuniform stationary mass. There
is only one feasible solution for each pair of (rr p , m0 ). Other solutions are rejected based on their
physical interpretations.
Table 4.1 shows scan results where the steady gliding path is presented with different sets of
center of mass distribution, location of movable mass, and net buoyancy, for both wing designs
illustrated in Section 3.2.2. The gliding angle θg = θ − α is the angle between Oxv and Ax with
gliding up as positive. zCG stands for the center of gravity expressed in the z-axis coordinate of the
body-fixed frame and there is a bijective function from rw3 to zCG :

zCG =

mw
r .
m w3

(4.10)

Here we ignore the influence of the excess mass m0 on zCG , which is really small compared to
that of mw . For example, m0 is generally around 30 g while mw is up to several kilograms. From
the results, we can see that different wing designs lead to different static gliding profiles. For
example, the larger wings result in shallower gliding paths (longer horizontal travel) but slower
total speed compared to the smaller wings, given the same set of control inputs. Since the wings
can be easily replaced in our design, we can potentially tailor the wing designs, while leaving the
glider body and its inside intact, to accommodate the requirements of different applications. On

39

the other hand, the results in the table show that, for a fixed wing design, the speed is influenced by
both the excess mass m0 and the pair (rr p , zCG ) while the pitch angle is affected mainly by the pair
(rr p , zCG ). Therefore, the center of mass plays an important role in determining the steady gliding
attitude. In particular, if we compare the cases where the values of zCG are different but the other
parameters are the same, we find that smaller zCG results in higher speed and larger glide angle.
This observation has been used in our design – by making zCG small, we can achieve desired glide
angle with very small displacement of the movable mass.

4.3 Scaling Analysis
We study the larger-wing glider model (Fig. 3.8a) at different scales and introduce a new cost
performance index, which reflects the horizontal travel distance per unit energy consumption. For
one dive (descent and ascent), the horizontal travel distance Dd is approximated as
V h
Dd = Vhtd = 2 h ,
Vv

(4.11)

where Vh and Vv are the steady-state horizontal speed and vertical speed, respectively, td is the
travel time for one dive, and h is the vertical travel depth. The energy consumption in one dive
comes from two sources, the pump actuation and the movable mass actuation, while the energy
consumed for the latter is negligible compared to that for pumping since the pump needs to overcome large pressure when the glider switches to ascent from descent. So the energy consumption
per dive Ed can be approximated as

Ed = ρ gh0 S p l p + ρ g(h0 + h)S p l p .

40

(4.12)

Table 4.2: Computed steady gliding path for the scaled models of the larger wing prototype. In
computation, r p = 5mm is used for the original scale model (1:1) while the value is scaled linearly
with dimension for other models.
scale

mass (kg) m0 (kg) Vh (m/s) Vv (m/s)
0.25:1
1
0.0075
0.063
0.018
0.5:1
2
0.015
0.11
0.039
1:1
4
0.03
0.19
0.094
2:1
8
0.06
0.28
0.207
4:1
16
0.12
0.39
0.377
8:1
32
0.24
0.54
0.574

Vh
-1
Vv m0 (kg )

glide ratio

488.35
203.25
74.55
28.30
12.01
5.72

3.5
2.82
2.02
1.35
1.03
0.94

Here, ρ is the water density, h0 is the equivalent water depth of the atmosphere pressure, S p is
the cross-section area of the pump tank inlet (and outlet) and l p represents the length of the water
column if the water pumped in each cycle is placed in a cylindrical container with cross-section
area S p . Noting the net buoyancy m0 = 12 ρ S p l p , we further simplify the energy consumption
per dive to Ed = 2m0 g(2h0 + h). Then we have the horizontal travel distance per unit energy
consumption

V
1
Dd
.
= h
Ed
Vv m0 1 + 2 h0

(4.13)

h

For a specific task, the depth is fixed and we have

Dd
V
∝ h ,
Ed Vv m0

(4.14)

which we define as the cost performance index τ .
We now conduct scaling analysis to examine how the cost performance metric evolves with
the dimensions (and consequently the weight) of the glider. CFD simulation shows that the drag
coefficient CD and lift coefficient CL stay almost the same at different scales we considered (from
0.25:1 to 8:1), while the pitch moment coefficient CM scales linearly with the characteristic dimen41

500
r = 4mm
p

450

r = 5mm
p

r = 6mm
p

350
300

0

−1

V / V / m (kg )

400

v

250

h

200
150
100
50
0
0.125

0.25

0.5

1

2

4

8

16

model scale in weight
Figure 4.1: The glider cost performance index with respect to model scales.

42

0.7
r = 4mm
p

r = 5mm

0.6

p

r = 6mm
p

0.4

h

V (m/s)

0.5

0.3
0.2
0.1
0
0.125

0.25

0.5

1

2

4

8

model scale in weight
Figure 4.2: The horizontal velocity with respect to model scales.

43

16

sion l of the glider. All related masses of the glider will scale as l 3 , including the movable mass
m¯ and the negative net buoyancy m0 . Taking the total length of the glider as l, the scale 1:1 would
imply l = 50 cm. By plugging those new parameter values into equations (4.7)–(4.9), we can solve
the glide path for the scaled model.
Table 4.2 shows the glide paths for glider models at different scales. Fig. 4.1 shows the relaV

h and the scale. The results show us that with a
tionship between the cost performance index Vv m
0

larger body the glider has a smaller glide ratio and a smaller cost performance index value, thus
consuming more energy for a given horizontal travel distance. This is consistent with the fact that
a larger glider needs to pump more water for a proper net buoyancy to provide the propelling force,
which is also the main energy consumption source. However, a larger-scale glider is able to achieve
faster horizontal speed as shown in Fig. 4.2. There is a trade-off between the achieved horizontal
speed and the horizontal distance coverage per unit energy cost, when selecting the optimal scale
for the glider. Other factors, like the dimension and the mass of the sensors and actuation devices,
should be also taken into account in the design process.

4.4 Experimental Results and Model Validation
With the developed miniature glider prototype, underwater gliding experiments are conducted in
order to validate the model. Most experiments are conducted in a large water tank that measures
15-foot long, 10-foot wide, and 4-foot deep.
For the steady gliding in the sagittal plane, we use the zero-angle tail. First, we set the initial
net buoyancy (negatively buoyant) and the linear actuator position to desired values. Then the
glider is released on the water surface with a stopwatch started, and the glider dives down until
it reaches the programmed depth indicated by the pressure sensor. Then it pumps water out and

44

resets its attitude to glide up. Here we focus on the gliding-down section and record the whole
gliding process with an underwater video camera fixed inside the water tank, pointing directly to
the glide motion plane (Fig. 4.3). Then we conduct video post-processing to extract the steady
gliding data, including the operating depth, horizontal travel distance and the time spent. We first
carry out one series of experiments with net buoyancy fixed at −20 g, and vary the linear actuator
positions. Then in another series of experiments, we vary the net buoyancy while holding linear
actuator position at 0.5 cm forward. For each setting of net buoyancy and linear actuator position,
we repeat the experiments 10 times and evaluate both the means and standard deviations for the
measured variables. Most experiments are conducted with the larger set of wings, but we have also
experimented with the smaller set of wings to further validate the model.
There are several things that must be given careful treatment. First, although the pressure
sensor can provide accurate depth measurement when the glider is at rest, the measurement is
subject to larger error when the glider is moving through the water. Therefore, we have chosen
to measure the actual depth each time. Second, we need the value of the center of gravity of the
glider zCG to obtain model-based predictions. We calculate it by hanging the glider up at different
points on the glider with strings and then taking the intersection of different hanging strings. This
value is further fine-tuned with collected tank test data. Third, to deal with the image distortion,
we have made a grid board for calibration. The board measures 2.5 m by 1.5 m, with grid cell size
of 10 cm by 10 cm. We fix the grid board in the glide motion plane and take its images. Then with
post-processing techniques, those grid images are incorporated into the gliding videos, to facilitate
the extraction of glide paths. Fourth, there is a period of transients before the glider reaches steady
gliding. To minimize the error introduced by the transients, we have chosen to use only the data
from last three seconds of descent in each run to compute the steady glide parameters.
Fig. 4.4 and Fig. 4.5 show the comparison between model predictions and experiment results
45

Figure 4.3: Superimposed snapshot of the gliding experiments in an indoor tank.

46

0
Experimental measurement
Model prediction

−5
−10
−15

θg (°)

−20
−25
−30
−35
−40
−45
−50
2

3

4

5

6

7

8

rp (mm)
Figure 4.4: Glide angle with respect to the movable mass displacement, with fixed net buoyancy
of -20 g, for the prototype with larger wings.

47

20
18
16

V (cm/s)

14
12
10
8
6
4
Experimental measurement
Model prediction

2
0
2

3

4

5

6

7

8

rp (mm)
Figure 4.5: Total glide speed with respect to the movable mass displacement, with fixed net buoyancy of -20 g, for the prototype with larger wings.

48

0
Experimental measurement
Model prediction

−5
−10

θg (°)

−15
−20
−25
−30
−35
−40
−26

−24

−22

−20

−18

−16

−14

m0 (g)
Figure 4.6: Glide angle with respect to the net buoyancy, with fixed movable mass displacement
of 0.5 cm, for the prototype with larger wings.

49

20

V (cm/s)

15

10

5
Experimental measurement
Model prediction
0
−26

−24

−22

−20

−18

−16

−14

m0 (g)
Figure 4.7: Total glide speed with respect to the net buoyancy, with fixed movable mass displacement of 0.5 cm, for the prototype with larger wings.

50

when we vary the movable mass position while holding the net buoyancy fixed, and Fig. 4.6 and
Fig. 4.7 show the results when the net buoyancy is changing while the movable mass location is
fixed. From the comparison results, we can see that the velocity and glide angle calculated from
the model match the experimental data reasonably well. In particular, the model has predicted well
the trends of how glide speed and angle vary with the center of the gravity and the net buoyancy.
We note that there are some non-negligible factors contributing to errors in the measurement.
When the glider starts gliding from rest, it is accelerating rather than steadily gliding for a few
seconds. We have already tried to remove the accelerating section from the data in the video processing process; however, it is difficult to completely eliminate that effect, especially considering
the relatively shallowness of the test tank. This effect would be reduced with deeper gliding; however, conducting precise measurement in a deep water body itself presents challenges. In addition,
flow disturbances in the tank influence the experiment results as well. So with these uncertainties,
we consider the match between our experimental results and the model predictions in Figs. 4.4 –
4.7 satisfactory.
We have further compared the glide performance when the glider is equipped with the larger
and smaller wings, respectively. Fig. 4.8 and Fig. 4.9 show the glide angle and the glide speed,
respectively, as the movable mass displacement is varied while the net buoyancy is held fixed at
−20 g. Both model predictions and experimental measurement are shown in the figures, and they
match well for both sets of wings. From the results, we can see that with smaller wings, the
glider tends to have a deeper gliding profile but higher speed, which matches our model predictions. These results further validate our derived model, and prove the effectiveness of our design
method with CFD-based evaluation of hydrodynamic coefficients in Section 3.2.2. Meanwhile, the
results also indicate that we can realize various gliding performance and meet different mission
requirements by replacing the wings, which are designed to be easy to change.
51

0
Experimental measurement (smaller wings)
Model prediction (smaller wings)
Experimental measurement (larger wings)
Model prediction (larger wings)

−5
−10
−15

θg (°)

−20
−25
−30
−35
−40
−45
−50
2.5

3

3.5

4

4.5

5

5.5

rp (mm)
Figure 4.8: Glide angle of the miniature glider with respect to the movable mass displacement,
with fixed net buoyancy of -20 g, when using two different sets of wings.

52

20
18
16

V (cm/s)

14
12
10
8
6
Experimental measurement (smaller wings)
Model prediction (smaller wings)
Experimental measurement (larger wings)
Model prediction (larger wings)

4
2
0
2.5

3

3.5

4

4.5

5

5.5

rp (mm)
Figure 4.9: Total glide speed of the miniature glider with respect to the movable mass displacement, with fixed net buoyancy of -20 g, when using two different sets of wings.

53

Table 4.3: Comparison of our miniature glider with underwater gliders reported in the literature.
Velocity refers to the terminal velocity. The bottom row shows the metrics of the glider reported
in this dissertation (with the larger wings).
Name
Length Weight
Slocum Electric Glider [47] 1.8 m
52 kg
Spray Glider [47]
2m
51 kg
Seaglider [47]
1.8 m
52 kg
ROGUE Lab Glider [73]
–
11.2 kg
ALBAC Glider [93]
1.4 m
45 kg
USM Glider [93]
1.3 m
–
ANT Littoral Glider [93]
2m
120 kg
Miniature Glider (this work) 0.5 m
4.2 kg

NetBuoyancy
520 g
900 g
840 g
360 g
–
–
–
20 g

Velocity
0.4m/s
0.45 m/s
0.45 m/s
–
0.51 m/s
–
1.03 m/s
0.17 m/s

The size of our miniature underwater glider is much smaller compared to traditional underwater
gliders, since we expect the future gliding robotic fish to be lightweight, easy to carry, and operate
in relatively shallow waters such as inland lakes and ponds. Table 4.3 shows some metrics of
existing gliders and our miniature glider prototype. From Table 4.3, we can see that the ratio
between net buoyancy and total weight of our miniature glider is similar to those of traditional
underwater gliders, at the order of 1 %, while the velocity magnitude is less than 1/3 of those typical
gliders. However, the cost performance index, the achieved horizontal travel distance per unit
energy consumption, by our glider is over 18 times bigger than those of reported gliders, assuming
they have similar gliding angles. All the above unique features of our miniature underwater glider
are consistent with the scaling analysis in Section 4.3, which could be used to design the glider
based on the specific applications.

54

Chapter 5
Steady Spiral in Three-Dimensional Space

5.1 Steady-State Spiraling Equations
If control inputs are fixed with nonzero tail angle, we can treat the influence of the tail on the
hydrodynamic forces and moments as the effects of increased hydrodynamic angles (α , β ), and
we know that the gliding robotic fish will perform three-dimensional steady spiraling motion (
[94], [81]), where the yaw angle ψ changes at constant rate while the roll angle φ and pitch angle

θ are constants. Then R T k is constant since






 0   − sin θ
  
  
T
T
R k = R  0  =  sin φ cos θ
  
  
1
cos φ cos θ









(5.1)

Taking time derivative of R T k , we have

RT k ) = 0,
ω b × (R

(5.2)

so the angular velocity has only one degree of freedom with ω3i in Oz axis in the inertial frame.
Then
RT k )
ω b = ω3i (R

55

(5.3)

The translational velocity in the body-fixed frame
T

v b = R bv

(5.4)

V 0 0

There are two important parameters in the spiraling motion: the turning radius R and the vertical speed Vvertical. By projecting the total velocity into the horizontal plane and vertical direction,
we have
T

Vvertical = R bv
R =

V 0 0

RT k

(5.5)

V 2 −Vv2 /ω3i

(5.6)

The steady-state spiraling equations are obtained by setting time derivatives to zero for the
robot’s dynamics:

R T k + F ext
0 = M v b × ω b + m0 gR

(5.7)

¯ r p × RT k
0 = J ω b × ω b + M v b × v b + T ext + mw grr w × R T k + mgr

(5.8)

From equations (3.1), (3.6), (5.3), (5.4) and the above steady-state spiraling equations, we know
there are six independent states for describing the steady spiral motion:
with

m0 r p1 δ

θ φ ω3i V α β

as the three control inputs. Expanding equations (5.7) and (5.8), then trans-

forming the original states to the above six independent states, we can obtain the nonlinear steady-

56

state spiraling equations as in (5.9) - (5.14).

α α 2 +C δ δ 2 )cα cβ
0 = m2 sβ V cφ cθ ω3i − m3 sα cβ V sφ cθ ω3i − m0 gsθ − 1/2ρ V 2 S(CD0 +CD
D
β

δ δ )cα sβ + 1/2ρ V 2 S(C +C α α )sα
−1/2ρ V 2 S(CSF β +CSF
L0
L

(5.9)

α α 2 +C δ δ 2 )sβ
0 = −m3 sα cβ V sθ ω3i − m1 cα cβ V cφ cθ ω3i − 1/2ρ V 2 S(CD0 +CD
D
β

δ δ )cβ + m gsφ cθ
+1/2ρ V 2 S(CSF β +CSF
0

(5.10)

α α 2 +C δ δ 2 )sα cβ
0 = m1 cα cβ V sφ cθ ω3i + m2 sβ V sθ ω3i + m0 gcφ cθ − 1/2ρ V 2 S(CD0 +CD
D
β

δ δ )sα sβ − 1/2ρ V 2 S(C +C α α )cα
−1/2ρ V 2 S(CSF β +CSF
L0
L

(5.11)
β

2 + (m − m )sβ sα cβ V 2 + 1/2ρ V 2 S(C β − K sθ ω )cα cβ
0 = (J2 − J3 )sφ cθ cφ cθ ω3i
2
3
q1
3i
M
R

α α + K sφ cθ ω )cα sβ − m gr sφ cθ
−1/2ρ V 2 S(CM0 +CM
w w
q2
3i
P

β

δ δ )sα
−1/2ρ V 2 S(CM β + Kq3 cφ cθ ω3i +CM

(5.12)

Y

Y

2 + (m − m )cα cβ sα cβ V 2 − m gr sθ − mgr
¯ p1 cφ cθ
0 = (J1 − J3 )sθ cφ cθ ω3i
w w
3
1

β

α α + K sφ cθ ω )cβ
+1/2ρ V 2 S(CM β − Kq1 sθ ω3i )sβ + 1/2ρ V 2 S(CM0 +CM
q2
3i
P

R

(5.13)

β

2
0 = (J2 − J1 )sθ sφ cθ ω3i
+ (m1 − m2 )cα cβ sβ V 2 + 1/2ρ V 2 S(CM β − Kq1 sθ ω3i )sα cβ
R

α α + K sφ cθ ω )sα sβ + mgr
−1/2ρ V 2 S(CM0 +CM
¯ p1 sφ cθ
q2
3i
P

β

δ δ )cα
+1/2ρ V 2 S(CM β + Kq3 cφ cθ ω3i +CM

(5.14)

Y

Y

Here, we assume the mass matrix and inertia matrix have the following form:


0
 m1 0


M =  0 m2 0


0
0 m3





 J1 0 0





 J =  0 J2 0




0 0 J3

57









5.2 Computation of Spiral Path and Evaluation of Stability
5.2.1 Newton’s method for Solving the Steady-State Spiraling Equations
The steady-state spiraling equations are highly nonlinear due to the terms involving trigonometric
functions and inverse trigonometric functions. Given the angle of attack α , the sideslip angle β ,
and the velocity magnitude V , a recursive algorithm based on fixed-point iteration could potentially
be applied to solve the equations for the other system states and control inputs [80]. However, we
are more interested in the converse problem of how to calculate steady-state solutions given fixed
control inputs, which are more useful for path planning and control purposes. Unfortunately, this
problem does not admit analytical solutions and the convergence condition for the corresponding
fixed-point problem is not satisfied. In the following we apply Newton’s method to solve the
problem.
Let x = [ θ φ ω3i V α β ]T be the six states that we want to solve for steady-state spiral
gliding equations. And let u = [ m0 r p1 δ ]T be the three control inputs. For convenience of
presentation, we write the governing equations in a compact form

0 = f (xx, u ) = [ fi (xx, u )]i=1,···,6

(5.15)

For example, f1 is the right hand side of equation (5.9).
The iterative algorithm for Newton’s method reads [95]

xˆ i+1 = xˆ i − J −1 (ˆx i , u ) f (ˆx i , u )

58

(5.16)

Here xˆ i is the ith-step iteration for the steady states, and J(xx, u) is the Jacobian matrix of f (xx, u)

J(xx, u ) =

∂f
∂ fi
=
∂x
∂ x j 6×6

(5.17)

The first row elements of the Jacobian matrix are given in equations (5.18) - (5.23) while the
others are omitted for succinct presentation, which can be calculated similarly.

∂ f1 /∂ x1 = −m2 sβ V cφ sθ ω3i + m3 sα cβ V sφ sθ ω3i − m0 gcθ

(5.18)

∂ f1 /∂ x2 = −m2 sβ V sφ cθ ω3i − m3 sα cβ V cφ cθ ω3i

(5.19)

∂ f1 /∂ x3 = m2 sβ V cφ cθ − m3 sα cβ V sφ cθ

(5.20)

α α 2 +C δ δ 2 )cα cβ
∂ f1 /∂ x4 = m2 sβ cφ cθ ω3i − m3 sα cβ sφ cθ ω3i − m0 gsθ − ρ V S(CD0 +CD
D
β

δ δ )cα sβ + ρ V S(C +C α α )sα
−ρ V S(CSF β +CSF
L0
L

(5.21)

α α cα cβ + 1/2ρ V 2 S(C +C α α 2 +C δ δ 2 )sα cβ
∂ f1 /∂ x5 = −m3 cα cβ V sφ cθ ω3i − ρ V 2 SCD
D0
D
D
β

δ δ )sα sβ + 1/2ρ V 2 SC α sα
+1/2ρ V 2 S(CSF β +CSF
L

+1/2ρ V 2 S(CL0 +CLα α )cα

(5.22)

α α 2 +C δ δ 2 )cα sβ
∂ f1 /∂ x6 = m2 cβ V cφ cθ ω3i + m3 sα sβ V sφ cθ ω3i + 1/2ρ V 2 S(CD0 +CD
D
β

β

δ δ )cα cβ
−1/2ρ V 2 SCSF cα sβ − 1/2ρ V 2 S(CSF β +CSF

(5.23)

Based on the parameters of the miniature underwater glider prototype as listed in Table 5.1,
Newton’s iterative formula is used to solve the steady-state spiraling equations. Characteristic
parameters for steady spiraling motion, including the turning radius and ascending/descending
speed, are obtained with different inputs as shown in Table 5.2. To apply Newton’s method, the

59

Table 5.1: Parameters of the lab-developed underwater robot used in the steady-state spiraling
equations.
Parameter
m1
m3
CD0

Value
3.88 kg
5.32 kg
0.45

CF

-2 rad−1

J1
J3

0.075
0.8 kg·m2
0.08 kg·m2

β
S
CL0

β
CM
R
β
CM
Y
Kq1

Kq3

Parameter
m2
m¯
α
CD

Value
9.9 kg
0.8 kg
17.59 rad−2

CFδ

1.5 rad−1

CM0

19.58 rad−1
0.05 kg·m2
0.0076 m

S
CLα
J2

-0.3 m/rad

α
CM

0.57 m/rad

5 m/rad
-0.1 m·s/rad
-0.1 m·s/rad

δ
CM
Y
Kq2
S

-0.2 m/rad
-0.5 m·s/rad
0.012 m2

P

initial values of states for iteration are chosen to be θ = −10◦ , φ = −10◦ , ω3i = 0.1 rad/s, V = 0.3
m/s, α = 0◦ , β = 0◦ . From the calculated results, we can see that a small turning radius requires
a large tail angle, a large displacement of movable mass, and a small net buoyancy, while a low
descending or ascending speed demands a small tail angle, a small displacement of movable mass,
and a medium net buoyancy.

5.2.2 Region of Convergence for Newton’s Method
For Newton’s method, the choice of the initial condition is important to the convergence of the
algorithm. Here, we numerically explore the region of convergence. For a fixed set of control
inputs, e.g., r p1 = 5 mm, m0 = 30 g and δ = 30◦ , we carry out the convergence test by running
the algorithm starting from different initial values of the states, and record whether a given initial
condition leads to convergence. Fig. 5.1a shows the convergence test results for different initial
conditions of the roll angle φ , pitch angle θ and spiraling speed V while the initial values of the

60

Table 5.2: Computed spiraling steady states through Newton’s method.
m0 (g)
25
25
25
25
25
10
15
20
30
35
40
25
25
25
25
25
25

r p1 (cm)
0.3
0.4
0.5
0.6
0.7
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5

δ (◦ )
(θ , φ , ω3i ,V, α , β ) (◦ ,◦ ,rad/s,m/s,◦,◦ )
(Vvertical, R) (m/s,m)
45
(-44.5, -31.0, 0.425, 0.264, -0.914, 4.10)
(0.182, 0.450)
45
(-46.8, -36.6, 0.448, 0.267, -1.32, 4.52)
(0.190, 0.417)
45
(-48.3, -40.6, 0.464, 0.268, -1.61, 4.87)
(0.195, 0.396)
45
(-49.3, -43.8, 0.476, 0.267, -1.84, 5.18)
(0.197, 0.380)
45
(-50.2, -46.5, 0.486, 0.267, -2.04, 5.48)
(0.211, 0.338)
45
(-70.8, -49.3, 0.589, 0.184, -3.64, 7.36)
(0.169, 0.121)
45
(-63.5, -52.7, 0.571, 0.218, -3.30, 6.98)
(0.189, 0.190)
45
(-55.5, -47.8, 0.517, 0.247, -2.46, 5.85)
(0.197, 0.287)
45
(-42.1, -34.3, 0.423, 0.281, -0.901, 4.24)
(0.185, 0.500)
45
(-36.9, -29.3, 0.392, 0.289, -0.306, 3.85)
(0.172, 0.591)
45
(-32.3, -25.3, 0.368, 0.293, 0.224, 3.60)
(0.157, 0.670)
30
(-37.6, -11.9, 0.235, 0.242, 0.854, 2.19)
(0.151, 0.806)
35 ( -43.4, -20.7, 0.311, 0.258, 0.0698, 2.87)
(0.178, 0.602)
40
( -46.8, -31.2, 0.389, 0.266, -0.761, 3.77)
(0.192, 0.474)
50
( -49.2, -48.8, 0.537, 0.264, -2.54, 6.19)
(0.192, 0.337)
55
(-51.1, -56.4, 0.615, 0.257, -3.62, 7.86)
(0.190, 0.283)
60
( -55.0, -63.8, 0.705, 0.247, -4.95, 10.0)
(0.189, 0.225)

other three states are set as α = 0◦ , β = 0◦ , ω3i = 0.1 rad/s; Fig. 5.1b shows the results for different
initial conditons of the angle of attack α , sideslip angle β and the angular speed ω3i with the initial
values of the other three states set as φ = −10◦ , θ = −10◦ ,V = 0.3 m/s. From the results, we
see that a small roll angle, a small pitch angle and a large velocity in the reasonable range will
lead to convergence; and the signs of the angle of attack and sideslip angle are very important
to the convergence of the solution. These observations offer insight into how to properly choose
the initial conditions when running the Newton’s method to obtain the steady spiraling path; for
example, one may want to select zero degree for the initial values of the angle of attack and sideslip
angle when having no idea about the signs of those two variables.

61

1
0.8

V

0.6
0.4
0.2
0
0
−20

0
−20

−40

−40

−60

θ (°)

−60
−80

−80

φ (°)

(a) Convergence with respect to the roll angle, pitch angle and
spiraling speed

ω3i (rad/s)

0.8
0.6
0.4
0.2
0
10
10
5

0

0
−5

β (°)

−10

−10

α (°)

(b) convergence with respect to the angle of attack, sideslip angle
and angular speed
Figure 5.1: Convergence test results for Newton’s method with respect to initial conditions. Color
yellow (light) means that convergence to the steady-state spiraling equilibrium is achievable with
the corresponding initial values; color blue (dark) means that there is no convergent solution or
the convergent solution is not at the steady-state spiraling equilibrium. In the test, the used set
of control inputs is r p1 = 5 mm, m0 = 30 g, δ = 45◦ ; and the corresponding equilibrium state
values are θ = −42.1281◦, φ = −34.2830◦ , ω3i = 0.4229 rad/s, V = 0.2809 m/s, α = −0.9014◦ ,
β = 4.2414◦.
62

5.2.3 Stability Analysis of Spiraling Motion
It is of interest to understand the stability of the spiral motion under a given set of control inputs.
Global stability analysis, however, is very difficult if not impossible due to the highly nonlinear
dynamics of the system. In this subsection we perform local stability analysis for the steady spiraling motion obtained from (5.9)–(5.14). A solution to these equations can be considered a relative
equilibrium of the system since it is independent of the coordinates of the robot in the inertial
frame. We denote with Jd (xxd , u ) the Jacobian matrix for the dynamics (3.4) and (3.5), where in a
compact form the system state vector is represented as xd = [ v b ω b ]T and the system dynamics
as x˙ d = f d (xxd , u ). As the (relative) equilibrium point is computed using a different set of system
states x = [ θ φ ω3i V α β ]T , and the Jacobian matrix J(xx, u ) for the steady-state equations has been evaluated with those states (Section 5.2.1), we can evaluate Jd through the chain
rule



Jd (xxd , u) = 

M −1
0

J



0 
dxxd −1
x
u
 J(x , )
dxx
−1



(5.24)

0
0
0
cα cβ −V sα cβ −V cα sβ




0
0
0
sβ
0
V cβ




0
0
0
sα cβ V cα cβ −V sα sβ
dxxd

where
= hspiral(xx)6×6 = 

dxx
 −cθ ω3i
0
0
0
0
−sθ



 −sφ sθ ω3i cφ cθ ω3i sφ cθ
0
0
0


0
0
0
−cφ sθ ω3i −sφ cθ ω3i cφ cθ





















So the value of linearization matrix Jd at the equilibrium point can be obtained by just plugging
the steady state values xe computed in Section 5.2.1 into the above equation, and by checking the
63

Hurwitz property of the linearization matrix Jd , we can understand the local stability property of
the steady-state spiraling motion.
We test the listed steady-state spiral motions in Table 5.2 for local stability. For example, for
the steady spiral corresponding to the control input set m0 = 25 g, r p1 = 0.3 cm, and δ = 45◦ , the
eigenvalues of Jacobian matrix Jd are −0.91 ± 5.02i, −6.69, −2.09, −0.42, −0.090, which shows
exponential stability. We find that all spirals in Table 5.2 have a Hurwitz Jacobian matrix and thus
the equilibrium of each spiral is locally asymptotically stable.

5.2.4 Basins of Attraction for the Spiraling Dynamics
The analysis in the previous subsection suggests that the relative equilibria associated with steady
spiraling are locally stable. It is of interest to gain some insight into the sizes of basins of attraction
for those equilibria. In this subsection, we run the dynamics simulation starting from different
initial conditions for a given fixed control input, and then record whether each initial state configuration will lead to convergence to steady spiraling, and if yes, what is the approximate time
it takes to converge. Since one cannot visualize a state space of more than three dimensions, we
have chosen to visualize the basin of attraction in three-dimensional subspaces of the original state
space.
Fig. 5.2 shows the simulation results of convergence time based on the parameters of our experimental prototype with respect to three states φ , θ , v1 . To obtain the results shown in this figure,
the control inputs are given as r p1 = 5 mm, m0 = 30 g and δ = 45◦ . Following this simulation
method, we can get the basins of attraction with convergence time for any other set of control
inputs. Similarly, we can obtain the convergence test results, shown in Fig. 5.3, when we vary
the initial conditions for the angular velocities in the body-fixed frame. From the results, it seems
that the basin of attractions for the spiraling dynamics is very large, which means that, starting
64

70

v1

1
0.8

60

0.6

50

0.4

40

0.2
30
0
0

20
−20
−40
−60

θ (°)

−80

−80

−60

−40

−20

0

10
0

φ (°)

(a) Display in orthogonal slice planes

300
1

250

0.8
200

v1

0.6
0.4

150

0.2
100

0
100

−100
50

50
0

0
−50

φ (°)

−100

100

θ (°)

0

(b) Display in a half sphere surface
Figure 5.2: Convergence time in seconds for spiraling dynamics with respect to different initial
values of states in roll angle φ , pitch angle θ , and translational velocity v1 along the Oxb direction,
for the control inputs of r p1 = 5 mm, m0 = 30 g and δ = 45◦ . The corresponding equilibrium state
values are θ = −42.1281◦, φ = −34.2830◦, and v1 = 0.2801 m/s.
65

65
60

0.5

ω3 (rad/s)

55
50
0

45
40
35

−0.5
0.5

30
0.5
0

ω2 (rad/s)

0
−0.5

−0.5

25
20

ω1 (rad/s)

Figure 5.3: Convergence time in seconds for spiraling dynamics with respect to different initial
values of states in angular velocities, for the control inputs of r p1 = 5 mm, m0 = 30 g and δ = 45◦ ,
displayed in orthogonal slice planes. The corresponding equilibrium state values are ω1 = 0.2837
rad/s, ω2 = −0.1767 rad/s, and ω3 = 0.2592 rad/s.

66













Figure 5.4: Snap shots of the robot spiraling in the experiment tank.
from almost every state configuration in a reasonable state value range, the glider is able to achieve
the steady spiraling motion eventually. However, we also notice that the convergence time varies
significantly with the starting condition. When the pitch angle and roll angle are negative, and
the speed is neither too high nor too low, the convergence time is relatively short. This provides
us with some ideas about when to switch to a desired gliding profile and how long we expect for
the transient period. We also notice that among all three angular velocity states, only the initial
condition of ω1 takes a noticeable influence on the convergence time of the glider dynamics. This
is consistent with the slow dynamics of the rotation motion in Oxb direction due to the enhanced
inertia from the large wings.

5.3 Experimental Results
With the miniature underwater glider prototype featuring a swappable tail fin, experiments are
conducted in order to confirm the spiraling motion and validate the derived mathematical model.

67

The experiments are performed in a large water tank that measures 15-foot long, 10-foot wide,
and 4-foot deep, as shown in Fig. 5.4. We set the net buoyancy (negatively buoyant), the linear
actuator position and the tail angle to fixed values. Then the glider is released on the water surface
and enters into the spiraling mode. Cameras are set to record the videos in both top view and
side view. The turning radii of the spirals are extracted after video processing. The comparison
between model predictions and experiment results on turning radius for different tail angles and
different excess masses are shown in Fig. 5.5 and Fig. 5.6, respectively. From the results, we can
see that the turning radius of the spiral is smaller with a larger deflected tail angle and a smaller
net buoyancy. The error bar in the figures shows the average value and standard deviation of the
turning radius out of ten repeated experiments. The model prediction shares the same trend with
the experimental results, and generally speaking, the match between those two are good.
There are some factors contributing to the measurement errors. First, the scales of camera
images are different at different distances. Here, an average scale is used in the information acquisition during video processing. A grid board is used for calibration, captured with the camera
at the average distance. Second, there are some initial transient processes, which is difficult to
be completely separate from the steady spiraling period. Experimental environment with deeper
water will effectively reduce the influence of initial transient; however, the complexity of experiments setup will be increased as a result. The environmental disturbances such as currents will
also affect the experimental results. So with these uncertainties, we consider the match between
our experimental results and the model predictions satisfactory.

68

1
Experimental Data
Model prediction

0.9
0.8

R (m)

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
25

30

35

40

45

50

55

60

65

δ (°)
Figure 5.5: Spiraling radius with respect to the tail angle, with fixed movable mass displacement
of 0.5 cm and fixed excess mass of 30 g.

69

0.6

R (m)

0.5
0.4
0.3
0.2
Experimental Data
Model prediction

0.1
0
24

26

28

30

32

34

36

m0 (g)
Figure 5.6: Spiraling radius with respect to the excess mass , with fixed movable mass displacement
of 0.5 cm and fixed tail angle of 45◦ .

70

Chapter 6
Passivity-based Stabilization with a
Whale-type Tail
In this chapter, the sagittal-plane stabilization problem of a gliding robotic fish with a whale-type
tail is discussed. The dynamics of an underwater glider in the sagittal plane is first reviewed and
separated into the slow dynamics and fast dynamics based on singular perturbation analysis. In
Section 6.2, a passivity-based nonlinear controller for the approximated reduced model is proposed, and stability analysis for the full system is conducted via linearization. Simulation results
are presented to show the effectiveness of the designed controller. A nonlinear observer is then
proposed in Section 6.3 for the implementation of the control strategy. In Section 6.4, both openloop and closed-loop experimental results are presented using a gliding robotic fish to illustrate the
effectiveness of both the controller and the observer.

6.1 Model of a Gliding Robotic Fish with a Whale-type Tail
In this section, the effect of a whale-type tail of a gliding robotic fish is considered as a control
surface. A dynamic model for underwater gliders in the sagittal plane is first reviewed [76, 78],
which is an invariant plane for such robots. Then with singular perturbation analysis, the model is
separated into two subsystems, a fast-mode system and a slow-mode system, for controller design
in Section 6.2.
71

xb

L

D
į

M2
O

Mį

Į

ș

xv
șg

Fį

m0g

zv zb

Figure 6.1: The schematic of a gliding robotic fish with forces and moments defined in the corresponding coordinate frames (side view).

6.1.1 Dynamic Model in the Sagittal Plane
This dissertation is focused on the motion in the robot’s invariant plane. As shown in Fig. 6.1, a
gliding robotic fish is modeled as a rigid-body system, with an ellipsoidal shape and fixed wings
as typically reported in the underwater glider literature. The relevant coordinate reference frames
are defined as follows: the body-fixed reference frame, denoted as Oxb yb zb , has its origin O at
the geometric center. The Oxb axis is along the body’s longitudinal axis pointing to the head;
the Ozb axis is perpendicular to Oxb axis in the sagittal plane of the robot pointing downwards,
and Oyb axis will be automatically formed by the right-hand orthonormal principle. In the velocity
reference frame Oxv yv zv , Oxv axis is along the direction of velocity whose magnitude is V , and Ozv
lies in the sagittal plane perpendicular to Oxv . In the inertial frame Axyz (not shown in Fig. 6.1), Az
axis points in the gravitational acceleration direction, and Ax is defined in the same direction as the
intersection line of the horizontal plane and the sagittal plane, while the origin A is a fixed point
72

in space. There are three angles defined following the standard convention in marine applications,
based on the aforementioned reference frames. The pitch angle θ is the angle between Oxb and Ax
with nose up as positive; the gliding angle θg is the angle between Oxv and Ax with gliding up as
positive; the angle of attack α is the angle from Oxb axis to Oxv axis with rotation axis Oyb .
We define the sum of the mass of the gliding robotic fish, mg , and the added mass in Oxb
direction as m1 , and similarly, the sum of mg and the added mass in Ozb direction as m3 . The robot
displaces a volume of water of mass mw . Let m0 = mg − mw represent the excess mass (negative
buoyancy).
The forces acting on the robot body include gravitational force, buoyancy force, hydrodynamic
forces (lift and drag) and control force. Due to the symmetric shape of the robot, the center of
buoyancy will be through the origin O. The assumptions in [76] are taken that the movable mass
is fixed at the origin O (during steady gliding), with the stationary mass distributed uniformly, and
the added masses are equally valued (m1 = m3 = m). Then the center of gravity will coincide with
the center of buoyancy at the origin. The force pair, gravitational force and buoyancy force, act like
one force of excess mass m0 at the origin O in Az direction. The hydrodynamic lift force L is along
negative Ozv axis, while the drag force D is along negative Oxv direction. The control force Fδ
is in Ozb direction, exerted by the control surface (e.g., a whale fluke-type tail) traveling through
the fluid medium, which is essentially another hydrodynamic force. The control surface angle δ
is defined as the angle between the control surface plane and the Oxb yb plane. The hydrodynamic

73

forces are dependent on the angle of attack and the velocity as follows:

L = (KL0 + KL α )V 2

(6.1)

D = KD0 + KD α 2 V 2

(6.2)

Fδ = KF V 2 uδ
δ

(6.3)

where KL0 , KL are lift coefficients, and KD0 , KD are drag coefficients. uδ is the effective angle of
attack that the control surface contributes to the gliding robotic fish. There is a linear relationship
between uδ and the control surface angle δ , uδ = Kuδ δ , where Kuδ is a scale constant. KF is the
δ

coupling factor that describes the additional force that the control surface induces.
There are two moments about the Oyb axis, which rotate the robot to a specific attitude. One is
the hydrodynamic pitch moment M2 , and the other is the control moment Mδ . They are modeled
as

M2 = KM0 + KM α + Kq2 ω2 V 2

(6.4)

Mδ = −KM uδ V 2

(6.5)

where KM0 and KM are pitch moment coefficients, Kq is the pitching damping coefficient and ω2
is the angular velocity for the pitch.
By applying Newton’s second law and the moment of momentum equation, the gliding robotic

74

fish dynamics are obtained as

1
V˙ = −
m g sin θg + D − Fδ sin α
m1 0
1
−m0 g cos θg + L + Fδ cos α
θ˙g =
m1V
1
α˙ = ω2 −
−m0 g cos θg + L + Fδ cos α
m1V
1
KM0 + KM α + Kq2 ω2 − KM uδ V 2
ω˙ 2 =
J2

(6.6)
(6.7)
(6.8)
(6.9)

where J2 is the total inertia about Oyb axis, consisting of stationary mass inertia and added inertia
in water, and g represents the gravitational acceleration.
For the open-loop system (i.e., uδ = 0), the steady gliding profile can be obtained from (6.6)(6.9). The state variables at the equilibrium have the following relationships

−KDe
,
θge = arctan
KLe

K
αe = − M0 ,
KM

ω2e = 0,



Ve = 

|m0 | g
2 + K2
KD
Le
e

1

2



where KDe = KD0 + KD αe2 , KLe = KL0 + KL αe .

6.1.2 System Reduction via Singular Perturbation
Bhatta and Leonard [76] have shown with singular perturbation analysis that for the open-loop
system, the dynamic model can be reduced to a second-order system with good approximation,

75

and the corresponding non-dimensional full state model is:
dV¯
1
m0 g sin θ¯g + θge + D − Fδ sin (α¯ + αe )
=−
dtn
KDe Ve2
d θ¯g
1
=
−m0 g cos θ¯g + θge + L + Fδ cos (α¯ + αe )
2
dtn
KDe Ve (1 + V¯ )

ε1
ε2

(6.10)
(6.11)

d α¯
1
−m0 g cos θ¯g + θge + L + Fδ cos (α¯ + αe )
= ω¯ 2 − ε1
2
dtn
KDe Ve (1 + V¯ )

d ω¯ 2
= − (α¯ + ω¯ 2 − uδ ) (1 + V¯ )2
dtn

(6.12)
(6.13)

where the new state variables are defined as

V −Ve
,
V¯ =
Ve

θ¯g = θg − θge ,

α¯ = α − αe ,

ω¯ 2 =

Kq
ω
KM 2

ε1 =

Kq 1
KM τs

the non-dimensional time tn and some related constants are defined as

τs =

m3
,
KDeVe

ε2 = −

J2 1
,
KqVe2 τs

tn = t/τs,

For the new state model, the hydrodynamic forces and moment can be described as

D = KD0 + KD (α¯ + αe )2 Ve2 (1 + V¯ )2

(6.14)

L = (KL0 + KL (α¯ + αe ))Ve2 (1 + V¯ )2

(6.15)

M2 = KM0 + KM α + Kq ω¯ 2 Ve2 (1 + V¯ )2

(6.16)

Fδ = KF uδ Ve2 (1 + V¯ )2
δ

(6.17)

76

The system can be further written in a compact form
dξ
= f (ξ , η , uδ )
dtn
dη
= Ag (ξ , η , ε , uδ )
µ
dtn

(6.18)
(6.19)

where,




¯
 V 

ξ =
, η = 
θ¯g

µ
 ε1 0
A=
0 εµ
2



α¯ 
,
ω¯ 2


,









 f1 
 g1 
f =
, g = 
,
f2
g2


 ε1 
ε =
 , µ = max(ε1 , ε2 ),
ε2

and f and g are defined accordingly based on (6.10) - (6.13).
From singular perturbation analysis, by setting µ = 0, we arrive at ω¯ 2 = 0 and α¯ = uδ . Plugging
those two fast-mode states into the other two state equations, the reduced model for the full system
is obtained. Now we further set α¯ = 0 in the reduced model for design convenience, since α¯ is
relatively small in value. Then the approximation of the reduced model can be expressed as
dξ
= f (ξ , 0, uδ )
dtn

(6.20)

and this second-order system will be used in the controller design. The effectiveness of designing the controller based on the approximated reduced system for the original full system, will be
demonstrated in next section.

77

6.2 Passivity-based Controller Design
6.2.1 Passivity-based Controller for the Approximated Reduced Model
The open-loop reduced model (6.20) with uδ = 0 has an exponentially stable equilibrium point
at the origin, which can be proven by Lyapunov analysis with the following positive definite Lyapunov function [76]

Φ=

1
2
− (1 + V¯ ) cos θ¯g + (1 + V¯ )3
3
3

(6.21)

f (ξ , 0, 0) ≤ −b1 ξ with b1 > 0.
and ∂∂ Φ
ξ
Now the objective is to design a feedback controller to stabilize the origin of the approximated
reduced model, which also provides a faster convergence speed. The approximated reduced system
is linear in control
dξ
= f (ξ , 0, 0) + gr (ξ ) uδ
dtn

(6.22)

where



2
¯
 KF (1 + V ) sin αe /KDe 
gr (ξ ) =  δ

¯
KF (1 + V ) cos αe /KDe

(6.23)

δ

For passivity-based controller design, an output yr needs to be defined for the approximated
reduced system, to make the system passive [96]. The output is chosen as

yr =

∂Φ
gr ( ξ )
∂ξ

78

(6.24)

where

∂Φ
=
∂ξ
=

∂Φ
∂ V¯

∂Φ
∂ θ¯g

− cos θ¯g + (1 + V¯ )2 (1 + V¯ ) sin θ¯g

(6.25)

We check the following expression for the approximated reduced model,
dΦ ∂ Φ
=
( f (ξ , 0, 0) + gr (ξ ) uδ )
dtn
∂ξ
It is known that

∂Φ
f (ξ , 0, 0) ≤ 0
∂ξ
So
dΦ
≤ uδ yr
dtn
Then by the definition of a passive system, the following system

dξ



= f (ξ , 0, uδ )
dtn
∂Φ


 yr =
gr ( ξ )
∂ξ

(6.26)

is passive. Let control uδ for system (6.22) be

uδ = −φ (yr )

for some function φ , where yr uδ = −yr φ (yr ) ≤ 0.

79

(6.27)

Now we take Φ in (6.21) as the Lyapunov function for the closed-loop system (6.26). Then
dΦ ∂ Φ
( f (ξ , 0, 0) + gr (ξ ) uδ )
=
dtn
∂ξ
∂Φ
∂Φ
=
f (ξ , 0, 0) +
gr ( ξ ) uδ
∂ξ
∂ξ
≤ −b1 ξ + yr uδ

For the robot controller design, there is limitation on the magnitude of the control variable uδ ,
so in this dissertation, we take

φ (yr ) =

1
arctan(yr )
Kc

(6.28)

where Kc is the control parameter that is used to limit the control output magnitude. We then have
dΦ
1
≤ −b1 ξ − yr arctan(yr )
dtn
Kc

(6.29)

which proves the asymptotic stability of the origin. Furthermore, the additional negative term
−

1
yr arctan(yr ) in the derivative of Lyapunov function provides an extra stabilization advantage.
Kc

With that term, the Lyapunov function will converge to zero more quickly, which results in a faster
convergence speed. That would help the robot to return to its steady gliding path with less time.

80

6.2.2 Stability Analysis for the Full Closed-loop System
From Eqs. (6.23) - (6.25), the designed controller uδ in (6.27) is dependent on the reduced model
states as

uδ =

1
1
(K cos αe sin θ¯g (1 + V¯ )2 + KF sin αe (1 + V¯ )2 ((1 + V¯ )2 − cos θ¯g )))
arctan(
δ
Kc
KDe Fδ
(6.30)

If applying the above controller to the full-order system (6.10) - (6.13), we will have the closedloop system, expressed in the compact form as

dz
= h (z, uδ (ξ ))
dtn

(6.31)

T

where the full system state vector z =

ξ η

.

It is challenging to establish the global stability of the origin in (6.31), so it is focused on the
local stability in this dissertation. The linearization matrix is defined as

A=

∂h
= ai j 4×4
∂z 0

81

(6.32)

The calculated elements in the Jacobian matrix A are

a11 = −2 −

2KF 2 sin2 αe
δ

Kc KDe 2

,

2
m g cos θge KFδ sin αe cos αe
a12 = − 0
,
−
KDeVe2
Kc KDe 2
2K αe
a13 = − D ,
KDe
2
m g cos θge KL αe 2KFδ sin αe cos αe
a21 = 0
,
+
−
KDe
KDeVe2
Kc KDe 2
2
2
m0 g sin θge KFδ cos αe
−
a22 =
,
KDeVe2
Kc KDe 2

a23 = KL /KDe ,
a14 = 0,

a24 = 0,

a31 = a21 ,

a32 = a22 ,

a33 = a23 ,

1 2KFδ sin αe
,
ε2 Kc KDe
1 KFδ cos αe
,
a42 = −
ε2 Kc KDe
a41 = −

a34 = 1/ε1 ,

a43 = −1/ε2 ,

a44 = −1/ε2 .

By examining the A matrix for the Hurwitz property, we would know whether the closed-loop
system with the passivity-based controller is asymptotically stable at its equilibrium. However, it
is technically difficult to check this 4-by-4 matrix directly unless we use a numerical approach.
But due to the time-scale separation property of system (6.31), there is an alternative way to check
the stability, by checking two 2-by-2 matrices on their Hurwitz property [97].
First, we break the matrix A into four blocks using four 2-by-2 matrices A11 , A12 , A21, A22 .

82

Here,




A12 
 A11
A=

1A
1A
µ 21 µ 22
Following [97], instead of checking matrix A, we only need to check 2-by-2 matrices A22 and
A11 − A12 A−1
22 A21 on whether they are Hurwitz, for proving the stability of the full-order system
when µ is sufficiently small. Plugging the calculated elements of linearization matrix A into those
four matrices, we have
KL
µ
 µK
ε1

De
A22 = 

µ
µ
−
−
ε2
ε2







,






 b1 b2 
A11 − A12 A−1
A
=


21
22
b3 b4

where

b1 = −2 −

2KF 2 sin2 αe
δ

2
Kc KDe

+

4µ 2 KF αe sin αe KD
δ

2
ε1 ε2 Kc KDe

2
2
m0 g cos θge KFδ sin αe cos αe 2µ KFδ KD αe cos αe
−
+
b2 = −
2
2
KDeVe2
Kc KDe
ε1 ε2 Kc KDe

b3 =

2
2
m0 g cos θge KL αe 2KFδ sin αe cos αe 2µ KFδ KL sin αe
+
−
−
2
2
KDe
KDeVe2
Kc KDe
ε1 ε2 Kc KDe

2
2
m0 g sin θge KFδ cos αe µ KFδ KL cos αe
b4 =
−
−
2
2
KDeVe2
Kc KDe
ε1 Kc KDe

For matrix A22 , the characteristic equation is

λ 2 + k1 λ + k0 = 0

83

where

k1 =

µ
K
−µ L ,
ε2
KDe

k0 =

µ2
K µ
−µ L
ε1 ε2
KDe ε2

Because µ ∼ O(εi ), i.e., µ is in the same order of the time scale εi , and they are both sufficiently
small, it can be shown easily that the coefficients k1 and k0 are both positive, which means A22 is
Hurwitz.
Now for matrix A11 − A12 A−1
22 A21 , let the coefficients of the characteristic equation be l1 and
l0 , defined similar to k1 and k0 . It is easy to get

l1 = −(b1 + b4 ),

l 0 = b1 b4 − b2 b3

(6.33)

Here we exploit information about the gliding robotic fish system and set a range for those
parameters. From the gliding robotic fish application perspective, we have

sin αe ≈ αe ,
KLe ∼ O(10),

cos(αe ) ≈ 1,
KL ∼ O(100),

0 ≤ KF ∼ O(10),
δ

KD ∼ O(100),

Ve ≤ 1,

KDe ∼ O(10),

Kuδ ∼ O(1),

0.25 ≤ |θge | ≤ 0.75

where metric units are applied to all above parameter values. Besides, we bound the control parameter Kc in an open set (1, 10) to restrain the magnitude of the control uδ from 9◦ to 90◦ , which is
consistent with the application constraints. Then we find that for each element in the 2-by-2 matrix
A11 − A12 A−1
22 A21 , there is a dominant term, which determines the sign of that element, with all
other minor terms only influencing the value. We express the approximation using the dominant

84

terms for the matrix elements

b1 ≈ −2,

b2 ≈ −

m0 g cos θge
,
KDeVe2

b3 ≈

m0 g cos θge
,
KDeVe2

b4 ≈

m0 g sin θge
KDeVe2

Since m0 and θge are always opposite in sign, b4 < 0. We also notice that b2 and b3 are opposite
in sign. So from equation (6.33), we have l1 > 0 and l0 > 0. With all characteristic equation
coefficients positive, the matrix A11 − A12 A−1
22 A21 is Hurwitz.
From the above analysis, we have shown that the original 4-by-4 linearization matrix A is
Hurwitz for sufficiently small µ , which proves that the passivity-based controller derived from the
approximated reduced model also stabilizes the full-order system. Furthermore, by the fact that
the controller is designed based on the approximated reduce system through passivity analysis,
we conjecture that this controller will be similarly effective for the full system as it does for the
approximated reduced system. In particular, we anticipate that the controller will provide a faster
convergence speed than the open-loop controller uδ = 0, due to the additional negative term it
introduced into (6.29).

6.2.3 Simulation Results
We apply the passivity-based controller to the full dynamics model, and use Matlab Simulink
to simulate the controller performance. The robot parameters we used are: m = 10 kg, J2 =
0.08 kgm2 , KL0 = 0 kg/m, KL = 303.6 kg/m, KD0 = 3.15 kg/m, KD = 282.8 kg/m, Kq = −0.8 kgs,
KM0 = 0.39 kg, KM = −14.7 kg, δ = 29.5, m0 = 0.05 kg. The equilibrium point is Ve = 0.24 m/s,

θge = −22.5◦ , αe = 1.52◦ and ω2e = 0 rad/s. Suppose that we have a current disturbance that
makes the robot deviate off its steady gliding path. From that point we want the robot to return to its
equilibrium gliding profile. The initial states are given as V0 = 0.2 m/s, θg0 = −35◦ , α0 = 1◦ and
85

−18
−20
−22
−24
Closed loop
Open loop

g

θ (°)

−26
−28
−30
−32
−34
−36
−38
0

10

20

30

40

50

60

t (s)
Figure 6.2: Simulation results on the trajectories of the gliding angle θg for the open-loop uδ = 0
and closed-loop (Kc = 2) cases, respectively.

ω20 = 0 rad/s. In simulation, we also consider the dynamics of the actuator for moving the control
surface, approximated by a first-order system with a time constant of 10 ms. The simulation time
is 60 seconds.
Fig. 6.2 shows that the passivity-based controller designed for the reduced model works for
the original full-order system, not only stabilizing the steady gliding equilibrium but also speeding
up the convergence process as we expected from the analysis. Figs. 6.3-6.4 show the influences
of the control parameter Kc on the control output and the glide angle transients. It can be seen
that with a smaller Kc , the system converges faster but requires larger initial control effort. With
the arctangent function in (6.28), the tunable parameter Kc makes it easy to balance between the

86

25

20

Kc = 2

δ

°

u ()

15
Kc = 5

10

Kc = 8
5

0

−5
0

5

10

15

20

25

30

t (s)
Figure 6.3: The trajectory of control uδ for the closed-loop simulation with different values for Kc .

87

−20
−22
−24

Kc = 2

−28

Kc = 5

g

θ (°)

−26

−30

Kc = 8

−32
−34
−36
−38
0

5

10

15

20

25

30

t (s)
Figure 6.4: The trajectory of gliding angle θg for the closed-loop simulation with different values
for Kc .

88

control effort and the convergence speed.

6.3 Observer Design
In this section, we propose a nonlinear model-based observer to estimate the velocity V and the
gliding angle θg (see (6.30)), using only the pitch angle which can be measured from onboard
sensors. The local stability of the observer is analyzed with a constructed Lyapunov function [98].
A nonlinear observer is proposed in order to estimate system states with a relatively large convergence region. The observer gain structure is selected to be linear, to enable efficient computation
and onboard implementation in experiments. The gain is obtained by solving the Riccati equation
as in the Extended Kalman Filter, in order to take the robustness to the measurement noise into the
design consideration. The nonlinear observer can be expressed as

x˙ˆ = f (ˆx , uδ ) + Ho (θ − θˆ )

(6.34)

θˆ = h(ˆx )

(6.35)

Here, xˆ = [ Vˆ θˆg αˆ ωˆ 2 ]T is the estimated system state, f (ˆx , uδ ) is the observer dynamics, as
described in Eqs. (6.6)–(6.9), h(ˆx ) = θˆg + αˆ is the output function, and Ho is the observer gain.
Let Q be a 4-by-4 symmetric positive definite matrix, which denotes the process noise covariance for the state dynamics. Let R be a constant, representing the measurement noise variance. Let
P be the solution of the Riccati equation

Ao P + PATo + Q − PCoT R−1Co P = 0

89

(6.36)

where Ao and Co are the linearization matrix of the system,


1 ∂D
−
m1 ∂ V





1 ∂L


m1V ∂ V

Ao = 

1 ∂L

−

m1V ∂ V


 2V (KM0 + KM α )
J2
Co =

m0 g cos θg
−
m1

1 ∂D
−
m1 ∂ α

m0 g sin θg
m1V

1 ∂L
m1V ∂ α

−

m0 g sin θg
m1V
0

−

1 ∂L
m1V ∂ α
KMV 2
J2

0








0





1



2
Kq2V 
J2
e

0 1 1 0

Here, [·]e means that matrix elements are evaluated at the equilibrium point.

∂L
= KLV 2 ,
∂α

∂L
= 2V (KL0 + KL α ),
∂V

∂D
= 2KDV 2 α ,
∂α

∂D
= 2V (KD0 + KD α 2 ).
∂V

The existence of a positive definite matrix P is guaranteed by the observability of (Ao , Co ). If
P exists, the observer gain Ho can be designed as

Ho = PCoT R−1

(6.37)

By tuning the value of R, we can adjust the robustness of the observer to different levels of measurement noise.
The stability of the designed nonlinear observer is analyzed below. We first define the dynamics

90

system states as x = [ V θg α ω2 ]T . Then the system dynamics can be expressed as

x˙ = f (xx, uδ )

(6.38)

θ = h(xx)

(6.39)

Define the state estimation error as x˜ = x − xˆ , and then from Eqs. (6.34) (6.35) (6.38) (6.39) we
have the estimation error dynamics

x˙˜ = f (xx, uδ ) − f (ˆx , uδ ) − Ho (h(xx) − h(ˆx ))

(6.40)

Taylor series of functions f and h are taken at the equilibrium point, and we have

x˙˜ = (Ao − HoCo )˜x + η (˜x ,t)

(6.41)

Here η (˜x ,t) represents the sum of Taylor series terms that contain second-order and higher-order
x˜ . Furthermore, η (0,t) = 0. There exist positive constants c0 and k1 such that when estimation
error |˜x | < c0 and x is bounded, we have

η (˜x,t) ≤ k1 x˜ 2

(6.42)

We define a Lyapunov function V (˜x ) for the system (6.41)

V (˜x ) = x˜ T P−1 x˜

91

(6.43)

We take the time derivative of V (˜x)

T
V˙ (˜x ) = x˜ T P−1 x˙˜ + x˙˜ P−1 x˜

(6.44)

From Eqs. (6.36), (6.37), (6.41), (6.44), we have

V˙ (˜x ) = −˜x T (C0T R−1C0 + P−1 QP−1 )˜x + 2˜x T P−1 η (˜x ,t)

(6.45)

Here (C0T R−1C0 + P−1 QP−1 ) ≥ σmin 1 4×4 is positive definite, and σmin is a finite scaler. From
Eq. (6.42), the time derivative of Lyapunov function is bounded

V˙ (˜x ) ≤ −c1 x˜ 2 + c2 x˜ 3

(6.46)

Here c1 and c2 are positive scalar coefficients that depend on the noise signals and the system
matrices. The coefficients can be selected as

c1 = σmin ,

c2 = 2k1 P−1

(6.47)

It can be easily shown from Eq. (6.46) that the state estimation error x˜ dynamics is locally
exponentially stable. Furthermore, the local exponential stability of the full closed-loop system
can be verified by the separation principle [96] with the proven local exponential stability property
of the passivity-based controller and the nonlinear observer.
We also consider the influence of the system noise and measurement noise on the observer

92

stability. The noise-corrupted system dynamics can be expressed as

x˙ = f (xx, uδ ) + ν1

(6.48)

θ = h(xx) + ν2

(6.49)

where ν1 and ν2 represent the process noise and measurement noise, respectively. The state estimation error dynamics for such a system is

x˙˜ = f (xx, uδ ) − f (ˆx, uδ ) − Ho (h(xx) − h(ˆx)) + ν1 − Ho ν2

(6.50)

Taylor series of functions f and h are taken at the equilibrium point, and we have

x˙˜ = (Ao − HoCo )˜x + η (˜x ,t) + ξ (t)

(6.51)

Here ξ (t) = ν1 − Ho ν2 . Assuming that x is bounded and the noise signals ν1 and ν2 are also
bounded, there exist constant k1 and k2 such that

η (˜x ,t) ≤ k1 x˜ 2 ,

ξ (t) ≤ k2

(6.52)

We take the same Lyapunov function V (˜x ) for the system (6.51)

V (˜x ) = x˜ T P−1 x˜

93

(6.53)

From Eqs. (6.36), (6.37), (6.44), (6.51), the time derivative of V (˜x) is

V˙ (˜x ) = −˜x T (C0T R−1C0 + P−1 QP−1 )˜x + 2˜x T P−1 (η (˜x ,t) + ξ (t))

(6.54)

Because of the boundedness properties of the noise signal, the time derivative of Lyapunov
function is bounded

V˙ (˜x ) ≤ −c1 x˜ 2 + c2 x˜ 3 + c3 x˜

(6.55)

Here c1 , c2 and c3 are positive scalar coefficients that depend on the noise signals and the system
matrices. The coefficients can be selected as

c1 = σmin ,

c2 = 2k1 P−1 ,

c3 = 2k2 P−1

(6.56)

It can be shown from Eq. (6.55) that the state estimation error will be bounded

|˜x | ≤

c1
2c2

c
when 4c2 c3 − c21 ≤ 0 and 1 < c0
2c2

(6.57)

The condition inequalities in Eq. (6.57) can be checked for a given system with known system
matrices and the noise signal characteristics. We speculate that the stabilization output of the full
closed-loop system will also be bounded around the nominal value, based on the boundedness of
the state estimation error and the exponential stability of the designed passivity-based controller.
Simulation is carried out to examine the performance of the designed observer. The closed-loop
system with the passivity-based controller for stabilization is simulated with full state feedback.
The observer runs in parallel with the closed-loop dynamics with the measured pitch angle as the

94

−10

−15

°

θ+υ ( )

−20

−25

−30

−35

−40
0

5

10

15

20

25

30

t (s)
Figure 6.5: Simulated trajectory of the measured pitch angle which is corrupted with noise. The
measurement noise is modeled as Gaussian random process with R = 0.1.
observer input. In simulation, measurement noise υ (t), added to the system output θ (t), is modeled
as a Gaussian random process with zero-mean, E{υ } = 0, and variance R, E{υ (t)υ (τ )} = Rδ (t −

τ ). Fig. 6.5 shows the noise-corrupted system output, the measured pitch angle with the variance of
measurement noise R = 0.1. Fig. 6.6 shows both the real and estimated gliding angle trajectories.
Gliding angle θg is used here to evaluate the effectiveness of state estimation because it is the
signature variable in the sagittal-plane glide motion. Control parameter Kc = 2 is used. The
simulation results show that the proposed observer is able to estimate the system states with good
robustness to the measurement noise.

95

−22
−24
Estimated State
Real State

−26

g

θ (°)

−28
−30
−32
−34
−36
−38
0

5

10

15

20

25

30

t (s)
Figure 6.6: Simulation results: the trajectories of the gliding angle θg of the real state and nonlinear
observer estimation with measurement noise.

96

6.4 Experimental Results
In order to test the effectiveness of the proposed passivity-based controller (Section 6.2.1) and the
nonlinear observer (Section 6.3), we conduct both open-loop and closed-loop experiments using a
gliding robotic fish prototype “Grace”. The tail fin system in “Grace” is driven by a servo motor
(Hitec Servo HS-7980TH) through a chain transmission. In experiments, the tail is adjusted up or
down like a whale fluke to modulate the glide motion. A microcontroller (dsPIC 30F6014A) in
the robot runs the glide control program and provides storage memory for the measurement data.
“Grace” is also equipped with the inertial measurement units including gyros (ST LPY503AL),
accelerometers and a digital compass (ST LSM303DLMTR), which are used to measure the robot’s
attitude.

6.4.1 Filter Design
The IMU sensor output is corrupted with a high-frequency noise. The system output, namely,
the pitch angle, is computed from the accelerometer sensor output. In this dissertation, in order
to smooth the measured pitch angle, a second-order Butterworth digital filter is adopted. The
discrete-time transfer function of the filter can be expressed as
a + a1 z−1 + a2 z−2
H(z) = 0
1 + b1 z−1 + b2 z−2

(6.58)

where a0 , a1 , a2 , b1 , b2 are filter coefficients.
Let fr denote the ratio between the sampling frequency fs and the cutoff frequency fc

fr =

97

fs
fc

(6.59)

For the Butterworth filter, the above coefficients can be determined as the follows

ω ′2
a0 = a2 = c
c
a1 = 2a0
2(ωc′2 − 1)
b1 =
c
1 − 2 cos(π /4)ωc′ + ωc′2
b2 =
c
where ωc′ = tan(π / fr ), and c = 1 + 2 cos(π /4)ωc′ + ωc′2 .
In this dissertation, the sampling frequency is 50 Hz, and the cutoff frequency is designed to
be 10 Hz. The low-pass Butterworth filter for the pitch angle can be described using the following
recursive difference equation,

θ (n) =a0 θs (n) + a1θs (n − 1) + a2θs (n − 2) − b1 θ (n − 1) − b2 θ (n − 2)

(6.60)

where θ (n) and θs (n) are the Butterworth filter’s output and input (pitch angle) at the nth step,
respectively. In this dissertation, a0 = a2 = 0.20644, a1 = 0.41289, b1 = −0.3702, b2 = 0.19597.

6.4.2 Open-loop Experiments
Open-loop experiments are first conducted in a large indoor tank, which measures 5 m long, 3.3 m
wide, and 1.3 m deep, using the gliding robotic fish prototype named “Grace” (Fig. 6.7). The
robot is first released from the water surface with a fixed tail angle δ , and then water is pumped
into the robot’s body until the net buoyancy m0 reaches 50 g. Two seconds after the gliding
down motion is initiated, the observer is initialized in the microcontroller. During the following
period, the robot records the pitch angle readings from both onboard sensors and the observer. The
98

Tail

Wing

Figure 6.7: A snapshot of the open-loop experiment with fixed tail angle using gliding robotic fish
“Grace” in the lab tank.
readings are further compared together with the pure computer-based Matlab simulation result,
which is obtained by running the same observer dynamics continuously in Matlab based on the
sensor output history recorded onboard. Fig. 6.8 shows the experimental results on the pitch angle

θ for different fixed tail angles δ = 15◦ , 30◦ , 45◦ , where the measured pitch angle is compared
with the values derived from the state estimate with Eq. (6.35). First, we can observe that the
computer-based Matlab simulation of the nonlinear observer produces estimate trajectories almost
identical to those from the onboard observer, confirming that the microcontroller is able to execute
the observer with little loss in accuracy. Second, the match between the observer estimation and
sensor output further validates the design of the nonlinear observer.

99

Sensor feedback
Simulation result
Onboard observer

°

θ()

−15
−20
−25
−30
0

1

2

3

4

3

4

3

4

t (s)
(a)

Sensor feedback
Simulation result
Onboard observer

°

θ()

−10
−20
−30
0

1

2

t (s)
(b)

Sensor feedback
Simulation result
Onboard observer

°

θ()

0
−10
−20
−30
−40
0

1

2

t (s)
(c)
Figure 6.8: The trajectories of pitch angle of the on-board sensor reading, observer estimation and
computer-based simulation result, in the open-loop experiments using gliding robotic fish “Grace”.
(a) δ = 15◦ ; (b) δ = 30◦ ; (c) δ = 45◦ .
100

6.4.3 Closed-loop Experiments
In this subsection, experimental results on sagittal-plane stabilization using the passivity-based
controller are presented to verify the effectiveness of both the proposed nonlinear controller and
the nonlinear observer. Experiments are conducted in the Neutral Buoyancy Research Facility
(NBRF) at the University of Maryland, College Park, which is a water tank measuring 15 m across
and 7.5 m deep.
In the experiments, the robot is released from the water surface with a net buoyancy m0 = 50 g.
The tail angle δ is initially set to 60◦ . Then the robot enters gliding down motion. When the robot
reaches a preset depth of 1 m, the tail flaps to δ = 0◦ to provide an initial perturbation for the
stabilization process. Then the designed passivity-based controller (6.30) kicks in to stabilize the
system, with the state estimation using the nonlinear observer (Fig. 6.9). A Qualysis underwater
motion capture system, which consists of 12 underwater cameras and motion tracking software, is
used to record the whole stabilization process and analyze the motion afterwards (Fig. 6.10).
Fig. 6.11 shows the experimental results on three types of pitch angle trajectories including
onboard sensor measurements, onboard observer estimation, and motion capture system output,
together with the onboard gliding angle estimation, when no feedback control exists, with tail angle trajectory designed as in Fig. 6.12. Figs. 6.13–6.16 show the experimental results for passivitybased stabilization for two different values of the controller gain. From the experimental results,
we observe a good match among the pitch angle estimation, the on-board sensor reading, and the
motion capture system output. The results further verify the proposed nonlinear observer design.
Besides, the gliding angle converges to the equilibrium point around −23.5◦ . The results confirm that the passivity-based controller effectively stabilizes the sagittal-plane glide motion and
speeds up the convergence process, comparing with the glide angle trajectory in the open-loop

101

Figure 6.9: A snapshot of stabilization experiment using gliding robotic fish “Grace” in NBRF,
University of Maryland.

102

Figure 6.10: A snapshot of stabilization experiment in the view of the 12 underwater cameras of
the Qualysis motion capture system in NBRF, University of Maryland.

103

−16
On−board sensor θ
Qualysis system θ
On−board observer θ
On−board observer θ

−18
−20

g

°

angle ( )

−22
−24
−26
−28
−30
−32
−34
−36
0

2

4

6

8

10

12

t (s)
Figure 6.11: The trajectories of pitch angle and gliding angle in the experiment of gliding stabilization without feedback control.
case (Fig. 6.11). It also shows that with a smaller controller gain Kc , the arising time of the pitch
angle is shorter (1 second in Fig. 6.15 vs 2 seconds in Fig. 6.13) and the control output magnitude
is larger (25◦ in Fig. 6.16 vs 15◦ in Fig. 6.14), which is consistent with and complementary to
the simulation findings in Section 6.2.3. This provides insight into the control parameter design in
order to balance between convergence time and control efforts.
In experiments, we also implemented proportional controller and PI controller for glide stabilization for the purpose of comparison. Experimental results using a proportional controller and a
PI controller are shown in Figs. 6.17–6.18 and Figs. 6.19–6.20, respectively. We varied the controller parameters, the proportional gain KP and the integral gain KI in order to obtain a better

104

0
−10

°

δ()

−20
−30
−40
−50
−60
0

2

4

6

8

10

12

t (s)
Figure 6.12: The trajectory of predefined tail angle in the experiment of gliding stabilization without feedback control.

105

−14
On−board sensor θ
Qualysis system θ
On−board observer θ
On−board observer θ

−16
−18

g

°

angle ( )

−20
−22
−24
−26
−28
−30
−32
0

2

4

6

8

10

12

t (s)
Figure 6.13: The trajectories of pitch angle and gliding angle in the passivity-based stabilization
experiment with Kc = 3.

106

20
10
0

°

δ()

−10
−20
−30
−40
−50
−60
0

2

4

6

8

10

12

t (s)
Figure 6.14: The trajectory of tail angle in the passivity-based stabilization experiment with Kc = 3.

107

−16
−18
−20

°

angle ( )

−22
−24
−26
−28
On−board sensor θ
Qualysis system θ
On−board observer θ
On−board observer θg

−30
−32
−34
0

2

4

6

8

10

12

t (s)
Figure 6.15: The trajectories of pitch angle and gliding angle in the passivity-based stabilization
experiment with Kc = 1.

108

20
10
0

°

δ()

−10
−20
−30
−40
−50
−60
0

2

4

6

8

10

12

t (s)
Figure 6.16: The trajectory of tail angle in the passivity-based stabilization experiment with Kc = 1.

109

performance. We find that with a larger KP and KI , the arising time of the gliding angle trajectory is reduced, however, with the cost of larger overshoot and oscillation amplitude, which shows
similar tradeoff in tuning the passivity-based controller gain. Furthermore, in the comparison of
the experimental results between using P/PI control and passivity-based control, the differences
in arising time, percent overshoot, and oscillation time show that passivity-based controller has
an overall better performance than the P/PI controller, especially in terms of the balance between
convergence speed and control effort.

110

−14
−16
−18

°

angle ( )

−20
−22
−24
−26

On−board sensor θ
Qualysis system θ
On−board observer θ
On−board observer θg

−28
−30
−32
0

2

4

6

8

10

12

8

10

12

t (s)

(a)
20
10
0

°

δ()

−10
−20
−30
−40
−50
−60
0

2

4

6

t (s)

(b)
Figure 6.17: The trajectories in the comparative stabilization experiments with a proportional controller with KP = 2. (a) pitch angle and gliding angle; (b) tail angle.

111

−16
−18

angle (°)

−20
−22
−24
−26
On−board sensor θ
Qualysis system θ
On−board observer θ
On−board observer θ

−28
−30

g

−32
0

2

4

6

8

10

12

8

10

12

t (s)

(a)
20
10
0

°

δ()

−10
−20
−30
−40
−50
−60
0

2

4

6

t (s)

(b)
Figure 6.18: The trajectories in the comparative stabilization experiments with a proportional controller with KP = 3. (a) pitch angle and gliding angle; (b) tail angle.

112

−16
−18
−20

°

angle ( )

−22
−24
−26
−28
On−board sensor θ
Qualysis system θ
On−board observer θ
On−board observer θg

−30
−32
−34
0

2

4

6

8

10

12

8

10

12

t (s)

(a)
10
0

−20

°

δ()

−10

−30
−40
−50
−60
0

2

4

6

t (s)

(b)
Figure 6.19: The trajectories in the comparative stabilization experiments with a PI controller with
KP = 2, KI = 1. (a) pitch angle and gliding angle; (b) tail angle.

113

−10
On−board sensor θ
Qualysis system θ
On−board observer θ
On−board observer θg

°

angle ( )

−15

−20

−25

−30

−35
0

2

4

6

8

10

12

8

10

12

t (s)

(a)
30
20
10

°

δ()

0
−10
−20
−30
−40
−50
−60
0

2

4

6

t (s)

(b)
Figure 6.20: The trajectories in the comparative stabilization experiments with a PI controller with
KP = 2, KI = 3. (a) pitch angle and gliding angle; (b) tail angle.

114

Chapter 7
Yaw Stabilization Using Sliding Mode
Control

7.1 Problem Statement
Steady gliding motion is the most commonly used profile for underwater gliders, providing the
capability of sampling water in the field while saving energy at the same time. Setting the right
hand side of Eqs. (3.2)–(3.5) to zero, we can solve those equations for the steady glide path given
a fixed movable mass displacement r p and excess mass m0 , with zero tail angle [80]. Due to the
existence of ambient currents or disturbances, the robot is susceptible to yaw deviation from its
desired direction, beside the sagittal-plane perturbation discussed in Chapter 6, which makes yaw
angle stabilization very important.
For succinctness, we first rewrite the system dynamics Eqs. (3.2)–(3.5) in a compact form

x˙ = f (xx) + ∆1 (xx) + g(xx)(u + ∆2 (t, x, u))

(7.1)

y = h(xx) = ψ

(7.2)

where x is the system state, x = [ φ θ ψ v1 v2 v3 ω1 ω2 ω3 ]T , u = δ is the tail angle,
which is the control input in the current setting, and ∆1 (xx) and ∆2 (t, x, u) represent uncertainties.

115

The system output is chosen to be the yaw angle ψ . The function g(xx) is dependent on the state


0 3×1



1
2 δ
 −
 2m1 ρ V SCSF cos α sin β


1
δ cos β

ρ V 2 SCSF

2m
2


1
2 δ
g(xx) = 
 − 2m ρ V SCSF sin α sin β

3

1

δ sin α
−
ρ V 2 SCM

Y
2J

1


0



1
δ cos α
ρ V 2 SCM
Y
2J3


























(7.3)

The yaw motion stabilization problem is how to design the state-feedback controller for the tail
angle δ , to stabilize the yaw angle ψ to a desired value r in the presence of disturbances.

7.2 Sliding Mode Control for Yaw Stabilization
Sliding model control is a practical nonlinear control method, especially for robust stabilization
[96]. In this dissertation, we adopt sliding mode control for tail-enabled yaw motion stabilization
of gliding robotic fish. Basic control design procedure follows the approach described in [96]. We
further construct a simplified controller requiring only partial state feedback information based on
the derived controller [99].
In order to obtain the relative degree of the system, we take the time derivatives of h(xx)

˙ x) = ψ˙ = sin φ sec θ ω2 + cos φ sec θ ω3 = L f h(xx)
h(x

(7.4)

¨ x) = ψ¨ = L2 h(xx) + L∆ L f h(xx) + Lg L f h(xx)(u + ∆2 (t, x, u))
h(x
f
1

(7.5)

116

where L f h(·) represents the Lie derivative of function h(·) with respect to the vector field f (·) [96],
and L2f h(xx) is equal to L f L f h(xx).
˙ x) does not depend on control input u and h(x
¨ x) does, imply that
The results that h(x

Lg h(xx) = 0

(7.6)

Lg L f h(xx) = 0

(7.7)

so the relative degree of the system ρsys = 2.
From Frobenius Theorem [100], there exists a transform function T (xx), which converts the
original system to the normal form with system states [η ξ ]T .







 p1 (xx)

..
 
 

.

η   Φ(xx)  
=
=
 
 
 p7 (xx)

ξ
Ψ(xx)

 h(xx)


L f h(xx)










 = T (xx)







(7.8)

where

∂ pi
g(xx) = 0, for i = 1, 2, . . ., 7
∂x
ξ˙1 = ξ2

(7.9)
(7.10)

Let r denote the reference trajectory for the yaw angle, which would be a constant number in

117

the stabilization problem. Take




 r 
R= 
r˙

(7.11)

The yaw error vector e is




 ξ1 − r 
e = ξ −R = 

ξ2 − r˙

(7.12)

Then the error dynamics is expressed as

η˙ = f0 (η , ξ )

(7.13)

e˙1 = e2

(7.14)

e˙2 = L2f h(xx) + L∆ L f h(xx)
1

+ Lg L f h(xx)(u + ∆2 (t, x, u)) − r¨

(7.15)

Assume that η˙ = f0 (η , ξ ) is bounded-input-bounded-state stable with ξ as the input. Then we
design a sliding manifold

s = e2 + k0 e1

(7.16)

where k0 is a positive constant.
Sliding mode control can be taken as the following to cancel the known terms as in feedback

118

linearization,

u=−

1
k e + L2f h(xx) − r¨ + ν
Lg L f h(xx) 0 2

(7.17)

where ν is the switching component, and

∂ ψ˙
f
∂x
∂ ψ˙
g
Lg L f h(xx) =
∂x

L2f h(xx) =

(7.18)
(7.19)

cos φ sec θ ω2



 sec θ tan θ (cos φ ω2 + cos φ ω3 )

∂ ψ˙ 

=
05 × 1
∂x 



sin φ sec θ


cos φ sec θ

T















(7.20)

Or we can take the controller as the pure switching component,

u=ν

(7.21)

Then in either case, the s-equation
˙
can be written as

s˙ = Lg L f h(xx)ν + ∆(t, x , ν )

(7.22)

Suppose that the uncertainty satisfies the following inequality,
∆(t, x , ν )
≤ ρ (xx) + κ0 |ν |, 0 ≤ κ0 ≤ 1
Lg L f h(xx)
119

(7.23)

where ρ (xx) represents the upper bound of the uncertainty related to the system states.
Design the switching component

ν = −γ (xx )sat(s/ε )

(7.24)

where sat(s/ε ) is a high-slope saturation function with a small constant ε , used to reduce chattering, and γ (xx) ≥ ρ (xx)/(1 − κ0) + γ0 with constant γ0 to deal with the non-vanishing disturbance
∆(t, x, ν ) if that is the case.
In this dissertation, we choose

k
γ (xx) = k1 x − xe 22 + k3

(7.25)

where x e is the system equilibrium point, which can be calculated given a steady gliding profile
[82], and k1 , k2 , k3 are tunable controller parameters, determined by the uncertainty type and also
capable of adjusting closed-loop dynamics performance.
Based on the fact that the yaw angle ψ is the state we really care about, we further simplify the
sliding mode controller (Eq. (7.21),(7.24),(7.25)) to

k
u = −(k1 ψ − ψe 22 + k3 )sat(s/ε )

(7.26)

where the sliding mode controller only requires the information of yaw angle, making it easy
to implement. The effect of this simplification on stability can be compensated by increasing
controller parameters k1 , k2 and especially the non-zero constant k3 . Although this will increase
the tracking error in general, the controller implementation becomes much more simpler. The
effectiveness of the designed controller will be evaluated in simulation and experiments in the
120

Table 7.1: Parameters of gliding robotic fish “Grace”.
Parameter
m1
m3
CD0

Value
8.0 kg
10.8 kg
0.45

CF

-2 rad−1

J1
J3

0.075
1.27 kg·m2
0.13 kg·m2

CM

-0.3 m/rad

α
CM

0.71 m/rad

CM
Y
Kq1
Kq3

5 m/rad
-0.16 m·s/rad
-0.16 m·s/rad

δ
CM
Y
Kq2
S

-0.2 m/rad
-0.80 m·s/rad
0.019 m2

β
S
CL0

β

β

R

Parameter
m2
m¯
α
CD

Value
19.8 kg
1.6 kg
17.59 rad−2

CFδ

1.5 rad−1

CM0

19.58 rad−1
0.08 kg·m2
0.0076 m

S
CLα
J2

P

following sections.

7.3 Simulation
To evaluate the designed sliding mode controller, simulation is carried out in Matlab. The parameters used in the simulation is based on the gliding robotic fish prototype “Grace”, obtained via
scaling analysis on our previously developed prototype in [82] (Tab. 7.1).
The initial state values for the simulation are

φ = 0,
v3 = 0,

θ = −30◦ ,
ω1 = 0,

ψ = 30◦ ,

v1 = 0.27 m/s,

ω2 = 0,

v2 = 0,

ω3 = 0.

The controller parameters used in the simulation are

k0 = 10,

k1 = 10/30/50,

k2 = 0.8/1/1.2,

121

k3 = 0.01.

Yaw angle ψ (°)

35
30

k = 10

25

k = 30

1
1

k = 50

20

1

15
10
5
0
−5
0

20

40

60

80

100

Time (s)
Figure 7.1: Yaw angle trajectory with respect to different controller parameters k1 .

122

5

Tail angle δ (°)

0
−5
−10
−15
k = 10
1

−20

k = 30
1

k = 50

−25
−30
0

1

20

40

60

80

100

Time (s)
Figure 7.2: Tail angle trajectory with respect to different controller parameters k1 .

123

Sideslip angle β (°)

1

0

−1

−2

k = 10
1

k = 30
1

−3

k = 50
1

−4
0

20

40

60

80

100

Time (s)
Figure 7.3: Sideslip angle trajectory with respect to different controller parameters k1 .

124

Yaw angle ψ (°)

35
30

k = 0.80

25

k = 1.0

2
2

k = 1.20

20

2

15
10
5
0
−5
0

20

40

60

80

100

Time (s)
Figure 7.4: Yaw angle trajectory with respect to different controller parameters k2 .
In Figs. 7.1–7.3, the trajectories of yaw angle ψ , controller command δ , and sideslip angle

β are shown for varying controller parameter k1 while Figs. 7.4–7.6 shows the simulation results
when varying controller parameter k2 . From the results, we can see that under proper controller
parameter setting, the sliding mode controller is able to regulate the yaw angle, which is deviated
from the desired orientation, back to the original, zero angle within a relative short time. Consequently, the trajectory of the glider is adjusted to the desired path, with the heading orientation
being zero degree, as shown in Fig. 7.7. From the comparison under different controller parameters, we find that k1 and k2 control the balance between response speed and controller effort. With
larger k1 and smaller k2 , the system responses faster and requires bigger control output amplitude.

125

5

Tail angle δ (°)

0

−5

−10

k = 0.80
2

k = 1.0
2

−15

k = 1.20
2

−20
0

20

40

60

80

100

Time (s)
Figure 7.5: Tail angle trajectory with respect to different controller parameters k2 .

126

0.5

Sideslip angle β (°)

0
−0.5
−1
k = 0.80

−1.5

2

k = 1.0
2

−2
−2.5
0

k = 1.20
2

20

40

60

80

100

Time (s)
Figure 7.6: Sideslip angle trajectory with respect to different controller parameters k2 .

127

0
−2
end point
start point

−4
Z (m)

−6
−8
−10
−12
−14
1
0.5
0

Y (m)

0

20

10
X (m)

(a) Three-dimensional view

1

Y (m)

0.8
0.6
0.4
end point
start point

0.2
0
0

5

10
15
X (m)

20

25

(b) Top view for X-Y plane
Figure 7.7: Trajectory of the gliding robotic fish for controller parameters k1 = 30, k2 = 1.
128

Regarding parameter k3 , as in the sliding mode design principle, it should balance the steady-state
error and uncertainty tolerance capability. With a larger k3 , the controller is able to work under
larger uncertainty while leading to bigger steady-state error.

7.4 Experimental Results
We implement the sliding mode controller on our lab-developed gliding robotic fish “Grace”, and
conduct yaw angle stabilization experiments in a large indoor water tank that measures 15-foot
long, 10-foot wide, and 4-foot deep. Yaw angle is calculated in the on-board micro-controller
using real-time sensor feedback from compass, accelerometer, and gyro. The sampling time for
sensor reading is 100 ms. A discrete-time low pass filter is applied directly after the raw data to
smooth out high frequency noise, with the difference equation

o(n) = a × i(n) + (1 − a) × o(n − 1)

(7.27)

where i(n) and o(n) are the low-pass filter’s input and output at time n, respectively. Filter parameter a ∈ (0, 1) is selected to be 0.5 in the experiment.
In the experiments, we release the gliding robotic fish on the water surface with a deviated
yaw angle. Then the robot starts to pump water in and translate the movable mass forward to
the set value, m0 = 40 g, and r p = 5 mm. After the robot is fully submerged and enters gliding
down motion, the sliding mode controller is applied to regulate the yaw angle to the desired value,
defined as zero degree. The robot motion in the whole process is recorded with a top-view camera
hung on a guiding rail and a side-view underwater camera placed in the tank (Fig. 7.8).
The sensor reading is stored in the on-board memory and sent back to the laptop through

129

T=4s

T=6s

T=10s

T=14s

(a) Top view

T=4s

T=8s

T=12s

T=14s

(b) Side view
Figure 7.8: Snapshots of controlled motion with yaw stabilization using sliding mode controller,
under parameters k1 = 50, k2 = 1.

130

0

Yaw angle ψ (°)

−10
−20
−30
−40
−50
−60
0

5

10

15

20

Time (s)
Figure 7.9: Yaw angle trajectory using “Grace”
wireless communication when the robot surfaces. Due to the range limitation of the tank, we
conduct the experiment for about 15 seconds. That corresponds to the most of transient process,
leaving the period that is close to the steady-state motion unrecorded.
In Figs. 7.9 and 7.10, we plot experimental results for both yaw angle ψ and tail angle δ . From
the experiment, we can see that the proposed sliding mode controller is able to regulate the deviated
yaw angle to the desired set value, with reasonable tail angle amplitude. There exist sensor noise
and tail rotation dynamics, which contributes to the unsmooth yaw angle readings. More advanced
sensor and filtering method could be used to solve the problem.

131

60

Tail angle δ (°)

50
40
30
20
10
0
0

5

10

15

Time (s)
Figure 7.10: Tail angle trajectory using “Grace”

132

20

Chapter 8
Three-Dimensional Curve Tracking

8.1 Three-Dimensional Steady Spiral and Its Differential Geometry Features
The three-dimensional motion control for gliding robotic fish, in terms of curve tracking, is very
challenging because the influences of the control inputs on the robot’s locomotion are strongly
nonlinear and coupled. It is more convenient to look into the influence of control inputs on the
robot’s differential geometry features, such as curvature and torsion, because we can examine the
relationship between those geometric characteristic parameters and the control inputs by studying
the steady-state spiral motions of the robot.
We decompose an arbitrary three-dimensional curve into a set of continuously evolving spirals.
In this way, at any point of the space curve, there is an imaginary matching spiral curve with the
same curvature and torsion. With this interpretation, instead of using the Euclidean positions, we
will explore the task of three-dimensional curve tracking via designing and following continuously
evolving spirals from the point of view of differential geometry [101].
First, Let us review the results of the steady spiral motion discussed in Chapter 4.
There are three control variables available to manipulate the robot’s motion profile: the excess
mass m0 , the position of the movable mass r p , and the tail angle δ . From [80] and [81], we know
that when all three controls are fixed at non-zero values, the gliding robotic fish will perform three-

133

dimensional spiraling motion and finally enter a steady spiral, where the yaw angle ψ changes at a
constant rate while the roll angle φ and pitch angle θ remain constant.
The dynamics of the spiral motion, derived from (3.2)-(3.5), can be presented in a compact
form as
x˙ s = f (xxs , u ) = [ fi (xxs , u )]8×1

(8.1)

T

where x s =

φ θ v1 v2 v3 ω1 ω2 ω3

T

, and control inputs u =

r p m0 δ

.

The steady-state spiraling equations can be obtained by setting time derivatives to zero in (8.1)

0 = f (xxs , u ) = [ fi (xxs , u )]8×1

(8.2)

In a steady spiral, R T k is constant since






 0   − sin θ
  
  
R T k = R T  0  =  sin φ cos θ
  
  
cos φ cos θ
1









(8.3)

The angular velocity has only one degree of freedom with ω3i in Ox axis in the inertial frame

RT k )
ω b = ω3i (R

(8.4)

Therefore, in the system of algebraic equations (8.2), there are nine independent variables (including control inputs) for describing the steady spiral motion: ( φ θ ω3i r p m0 δ V α β )T .
Hereafter we will use a state transformation on linear velocity variables for the sake of calculation

134

convenience





vb = 











v1 
 V   V cos α cos β


 


 
R
= bv  0  = 
v2 
V sin β


 


 
v3
0
V sin α cos β









(8.5)

In the elementary differential geometry, a three-dimensional curve is captured by its curvature
and torsion. The curvature κ is the amount by which a geometric object deviates from being flat,
or the degree by which a geometric object bends, while torsion τ measures the departure of a curve
from a plane, or how sharply a curve twists. Any time-trajectory of a smooth space curve can
be completely described mathematically using curvature, torsion and velocity with Frenet-Serret
formulas [102].
The geometric parameters (curvature, torsion, velocity) of a steady spiral can be expressed as

r

κ =

r2 + c2
c
r2 + c2

τ =

(8.6)
(8.7)

where r is the steady spiral radius, and 2π c is the steady spiral pitch, or the vertical separation
between two steady spirals. Furthermore,

Vh
ω3i
Vv
c =
ω3i
r =

(8.8)
(8.9)

where Vh and Vv are the horizontal velocity and vertical velocity, respectively, of the steady spiral
motion.

135

We also have

Vh2 +Vv2 = V 2

(8.10)

From (8.6)-(8.10), the angular velocity ω3i and the vertical velocity Vv can be described by the
three geometric parameters, κ , τ , and V

ω3i = V

κ2 + τ2

(8.11)

Vv = V

τ2
τ2 + κ2

(8.12)

8.2 Influence of Control Inputs on Spiral Trajectory
As for the spiral motion, the three control inputs, including the movable mass displacement r p , the
net buoyancy m0 , and the tail angle δ , have different influences on both the steady-state motion
profile and the transient dynamics. In this dissertation, we focus on the influences of control inputs
on the steady-state spiral trajectory characteristics, which provide useful insight for path planning
in three-dimensional curve tracking.
We study the relationship between three system control inputs and three trajectory characteristic parameters, including curvature κ , torsion τ and total speed V , which are used to completely
describe any three-dimensional trajectory. With the system dynamic model in Chapter 3 and system parameters as in “Grace”, we conduct simulation with different sets of values of system control
inputs, and then record the corresponding steady-state spiral paths.
Fig. 8.1 shows the relationship between tail angle δ and the three trajectory characteristic
parameters, while the net buoyancy m0 and the displacement of movable mass r p are fixed at

136

30 g and 0.5 cm, respectively. Simulation results of varying m0 and r p are shown in Figs. 8.2–8.3.
From those figures, we see that all control inputs have significant influence on the motion profile,
although the degree of influence varies. For example, δ and r p have greater influence on κ and

τ than m0 . Most relationships show monotonic trends, while non-monotonic relationships appear
between κ /τ and m0 (Fig.8.2a). The simulation results of the influences of control inputs on robot’s
spiral trajectory, provide insight into the capability of three-dimensional maneuvering as well as
the feedback controller design for three-dimensional curve tracking.
In order to verify the relationship between the control inputs and the trajectory characteristic
parameters of steady-state spirals, experiments are conducted using the prototype “Grace”. The
experiments are carried out in the Neutral Buoyancy Research Facility (NBRF), University of
Maryland. The water tank measures 50 feet across and 25 feet deep.
In experiments, “Grace” are remotely controlled via XbeePro communication to perform spiral
motions with different control input values (Fig. 8.4). The whole spiral process is recorded using
a Qualysis underwater motion capture system. The motion capture system features 12 underwater
cameras around the water tank, with 8 at a shallower depth and 4 at a deeper depth. Each camera captures the spiral motion from a different angle of view (Fig. 8.5). The robot is equipped
with five markers, which the motion capture system uses to identify the rigid body. Some of the
robot’s states, such as linear and angular positions can be measured and outputted using the system
(Fig. 8.6), and other system states, including translational and angular velocities can be estimated
from those measurements.
Fig. 8.7a and Fig. 8.7b show the comparison results between the model prediction and the
experimental results on the spiral curvature and torsion when varying δ from 20◦ to 50◦ with
the m0 and r p fixed at 30 g and 0.5 cm, respectively. The results on the spiral total speed is not
presented as the influence of the tail angle on that variable is not very obvious as shown in Fig. 8.1b.
137

2

2

1

0
10

1

20

30

°

40

50

0
60

40

50

60

τ

κ

curvature
torsion

δ( )
(a)
0.072
0.071
0.07

V(m/s)

0.069
0.068
0.067
0.066
0.065
0.064
0.063
0.062
10

20

30

°

δ( )

(b)
Figure 8.1: Trajectory characteristic parameters at different tail angles with m0 = 30 g and r p = 0.5
cm.

138

1

1.4
curvature
torsion

0.8

1

0.6

τ

κ

1.2

0.8
10

20

30

40

50

60

70

80

0.4
90

m (g)
0

(a)
0.08

0.075

V(m/s)

0.07

0.065

0.06

0.055

0.05
0

20

40

60

80

100

m0(g)

(b)
Figure 8.2: Trajectory characteristic parameters at different net buoyancy with δ = 30◦ and r p =
0.5 cm.

139

2

1

1

0
2

0.5

4

6

8

10

12

14

0
16

12

14

16

τ

κ

curvature
torsion

r (mm)
p

(a)
0.14
0.13
0.12

V(m/s)

0.11
0.1
0.09
0.08
0.07
0.06
0.05
0.04
2

4

6

8

10

rp(mm)

(b)
Figure 8.3: Trajectory characteristic parameters at different displacements of movable mass with
δ = 30◦ and m0 = 30 g.

140

Figure 8.4: The gliding robotic fish “Grace” spiraling in Neutral Buoyancy Research Facility,
University of Maryland.

141

Figure 8.5: Snapshots of spiral motion with 12 underwater cameras from different angles of view
using a Qualysis underwater motion capture system.

Figure 8.6: Illustration of robot’s rigid body and coordinates in spiral motion using a Qualysis
underwater motion capture system.

142

The spiral experiments at each set of control input values are conducted five times. The mean and
standard deviation of κ and τ are provided with the error bar in the figures. Considering the current
disturbance in the water tank due to the boundary effects and constantly active filtering system, the
match between the model prediction and the experimental results is reasonable well, which further
validates the derived system model.
The experiments for varying m0 and r p are not carried out and presented here because the speed
and the gliding angle could be much increased so that the robot will bump into the metal frame
located in the center of the water tank (setup for other experiments and not removable currently).

8.3 Two Degree-of-Freedom Control Design
In this section, we propose a 2-DOF control strategy for the curve tracking problem where an
inverse mapping is used as a feedforward controller while a robust H∞ controller is used as the
feedback controller, which is designed based on the linearized model around the steady spiral
trajectory. The open-loop feedforward controller obtained from inverse mapping of steady spiral
motion serves as a driving force pushing the robot towards the desired steady spiral. However, the
convergence time is relatively long, which will be demonstrated in simulation later. The feedback
H∞ controller aims to speed up the convergence and enhance the performance robustness. The idea
of the 2-DOF controller is that with the feedforward inverse mapping, the dynamic nonlinearity is
reduced so that a feedback H∞ controller can be designed based on the linearized model to achieve
improved transient performance.

143

2

κ

1.5

1

0.5
10

20

30

40

°

50

60

70

50

60

70

δ()

(a)
1.6
1.4
1.2

τ

1
0.8
0.6
0.4
0.2
0
10

20

30

40

°

δ()

(b)
Figure 8.7: Model prediction and experimental results in spiral motion at different tail angles. The
other two control inputs are fixed m0 = 30 g, and r p = 0.5 cm. (a) curvature; (b) torsion.

144

8.3.1 Feedforward Control via Inverse Mapping of Steady Spiral Motion
In this subsection, we study the feedforward control for 3D curve tracking for gliding robotic fish.
Based on the fact that a three-dimensional curve can be decomposed into a set of continuously
evolving spirals, we propose a 3D curve tracking method by tracking geometric characteristics
of these spirals instead of following Euclidean positions. An open-loop feedforward controller
is designed using inverse mapping of the steady spiral motion. This inverse controller will also
become part of the proposed 2-DOF controller.
We will calculate the desired control inputs for a given steady spiral profile, which is parameterized by curvature, torsion and velocity. This inverse mapping solution can be adopted as an
open-loop feedforward controller for the 3D curve tracking problem.
With (8.4), it can be shown that the first two equations in (8.2) always hold, thus redundant.
With the value of Vv known from (8.12), we have one more constraint equation




V 



 T
R k)
Vv = R bv  0  (R




0

(8.13)

Given κ , τ and V , we can calculate the value of the angular velocity ω3i from (8.11). Knowing
the values of V and ω3i , there are seven unknown variables left out of nine independent states
for the steady spiral motion ( φ θ ω3i r p m0 δ V α β )T . Correspondingly, there are
seven independent algebraic equations from (8.2) and (8.12). The inverse mapping problem is then
formulated as
0 = g (xx) = [gi(xx)]7×1

(8.14)

where x = ( φ θ α β r p m0 δ )T . The expansion of gi (xx) is shown in (8.15)-(8.21).
145

α α 2 +C δ δ 2 )cα cβ
0 = m2 sβ V cφ cθ ω3i − m3 sα cβ V sφ cθ ω3i − m0 gsθ − 1/2ρ V 2 S(CD0 +CD
D
β

δ δ )cα sβ + 1/2ρ V 2 S(C +C α α )sα
− 1/2ρ V 2 S(CSF β +CSF
L0
L

(8.15)

α α 2 +C δ δ 2 )sβ
0 = −m3 sα cβ V sθ ω3i − m1 cα cβ V cφ cθ ω3i + m0 gsφ cθ − 1/2ρ V 2 S(CD0 +CD
D
β

δ δ )cβ
+ 1/2ρ V 2 S(CSF β +CSF

(8.16)

α α 2 +C δ δ 2 )sα cβ
0 = m1 cα cβ V sφ cθ ω3i + m2 sβ V sθ ω3i + m0 gcφ cθ − 1/2ρ V 2 S(CD0 +CD
D
β

δ δ )sα sβ − 1/2ρ V 2 S(C +C α α )cα
− 1/2ρ V 2 S(CSF β +CSF
L0
L

(8.17)
β

2 + (m − m )sβ sα cβ V 2 + 1/2ρ V 2 S(C β − K sθ ω )cα cβ
0 = (J2 − J3 )sφ cθ cφ cθ ω3i
2
3
q1
3i
M
R

α α + K sφ cθ ω )cα sβ − m gr sφ cθ
− 1/2ρ V 2 S(CM0 +CM
w w
q2
3i
P

β

δ δ )sα
− 1/2ρ V 2 S(CM β + Kq3 cφ cθ ω3i +CM

(8.18)

Y

Y

2 + (m − m )cα cβ sα cβ V 2 − m gr sθ − mgr
¯ p cφ cθ
0 = (J1 − J3 )sθ cφ cθ ω3i
w w
3
1

β

α α + K sφ cθ ω )cβ
+ 1/2ρ V 2 S(CM β − Kq1 sθ ω3i )sβ + 1/2ρ V 2 S(CM0 +CM
q2
3i
P

R

(8.19)

β

2
0 = (J2 − J1 )sθ sφ cθ ω3i
+ (m1 − m2 )cα cβ sβ V 2 + 1/2ρ V 2 S(CM β − Kq1 sθ ω3i )sα cβ
R

α α + K sφ cθ ω )sα sβ + mgr
− 1/2ρ V 2 S(CM0 +CM
¯ p sφ cθ
q2
3i
P

β

δ δ )cα
+ 1/2ρ V 2 S(CM β + Kq3 cφ cθ ω3i +CM
Y

Y

0 = Vv /V + cα cβ sθ − sβ cθ sφ − sα cβ cθ cφ

(8.20)
(8.21)

Unfortunately, there is no closed-form solution to this system of equations. In this dissertation,
we use a Newton’s method to find solutions recursively, which provide the desired open-loop

146

control inputs. The iterative algorithm for Newton’s method reads [95]

xˆ i+1 = xˆ i − J −1 (ˆx i )gg(ˆx i )

(8.37)

Here xˆ i is the ith-step iteration for x, and J(xx) is the Jacobian matrix of g(xx)

J(xx, u ) =

∂g
=
∂x

∂ gi
∂ x j 7×7

(8.38)

8.3.2 Linearized Model Around the Steady Spiral Trajectory
The linearized model can be obtained by linearization of the spiral dynamics around the equilibrium spiral trajectory. Recall the spiral dynamics (8.1)

x˙s = f (xxs , u ) = [ fi (xxs , u )]8×1

where x s = ( φ θ v1 v2 v3 ω1 ω2 ω3 )T , and u = ( r p m0 δ )T . We define transformed system states z = ( φ θ V α β ω1 ω2 ω3 )T for the convenience of computation
of the Jacobian matrix J(xxs , u ). The linearized system matrices are

A = [J(xxs , u )]e =
B=

∂f
∂f
=
∂x e
∂z

∂ x −1
∂z

(8.39)
e

∂f
∂u e

(8.40)

Here [·]e means the matrix elements are evaluated at the equilibrium point.
We define the linearized system output as y = ( φ θ v1 )T , linear in system states. The

147

linearized system output matrices are

C = [11 3×3 0 3×5]

(8.41)

D = 0 3×3

(8.42)

The objective of the controller design for the linearized model is to drive the three selected system
outputs to the desired values that are computed via inverse mapping (Section 8.3.1) at an improved
speed.

8.3.3 H∞ Controller Design
A 2-DOF control is adopted for the 3D curve tracking problem to increase the system bandwidth.
The 2-DOF control system configuration is shown in Fig. 8.8. The transfer function G(s) represents the spiral dynamics system. K f b (s) is the feedback controller. r = ( φr θr v1r )T is the
perturbation of the reference signal from the nominal values. u = ( r p m0 δ )T is the control
input to the plant. e = r − y represents the tracking error. We (s) is the user-defined weighting
function to impose the requirements for the tracking bandwidth and tracking error amplitude. The
state-space realization of We (s) is as follows,

x˙ w = Aw x w + Bw u w

(8.43)

z w = Cw x w + Dw u w

(8.44)

Wu = diag(wu1 , wu2 , wu3 ) is the weighting function to help control the magnitude of system control
inputs, and z u = Wu u .
The tracking performance can be characterized by the tracking error e . Meanwhile, the control

148

Figure 8.8: The control system diagram with combination of open-loop control and closed-loop
control.

Figure 8.9: Transformed 2-DOF control configuration in H∞ control framework .
effort can be characterized by the control input u which is desired to be small for the consideration
of energy consumption. The objective of the feedback control design is to minimize those two
signals e and u.
This optimization problem with the feedback control for the linearized model can be transformed into an H∞ robust control framework as shown in Fig. 8.9. In this H∞ control system

149

configuration,

K(s) = K f b (s)

(8.45)

z = ( z w z u )T

(8.46)

n =r−y

(8.47)

The interconnected system




 Ap Bp 
P=

Cp D p

(8.48)

where



Ap = 






Cp = 







B 
 0
Bp = 

Bw −Bw D


Wu 
 0




D p =  Dw −Dw D 




1
−D

0 

−BwC Aw

0
0 


−DwC Cw 


−C
0
A



The design objective is then to minimizing the H∞ norm of the transfer function from r to z

min Trz (s) ∞
K

(8.49)

To solve the above H∞ optimization problem, we adopt the command hinfmix(·) in Matlab LMI
toolbox, the output of which provides the feedback controller K(s).

150

8.4 Simulation Results
Given desired trajectories of curvature, torsion and velocity, simulation is conducted to test the
effectiveness of the proposed 2-DOF control algorithm. For the purpose of comparison, we have
also conducted simulation with pure inverse mapping control (see Section 8.3.1) and ProportionalIntegral (PI) control. The PI controller is designed as

δ = KPδ ∆κ + KIδ
rp

rp

r p = KP ∆τ + KI

m
m
m0 = KP 0 ∆V + KI 0

∆κ
∆τ
∆V

(8.50)
(8.51)
(8.52)

where ∆ stands for the difference between the desired value and actual value of the variable that
follows. The particular form of the PI controller, where one control input is only dependent on
the error feedback from one geometric parameter, is adopted for design convenience and based on
the observed influences (in simulation) of the control inputs on the geometric parameters, where it
appears that each control input has more pronounced impact on one of the geometric parameters
than other two inputs. The PI controller coefficients are designed as KPδ = 0.1, KIδ = 0.01, KP =
rp

rp
m
m
0.05, KI = 0.005, KP 0 = 0.1, KI 0 = 0.01. There are three control inputs and three geometric

variables to track, so the strong coupling between the control inputs makes the parameter tuning
quite challenging. The PI control parameters are tuned in simulation in order to obtain the best
tracking performance.
The model parameters used in simulation are based on the lab-developed gliding robotic fish

151

prototype “Grace” as in Table 7.1. The initial values of system states used in simulation are

θ = −7.2◦

φ =0

v3 = 0.04◦

ω1 = 0

v1 = 0.1 m/s

v2 = 0

ω2 = 0

ω3 = 0

This represents a sagittal-plane glide motion.
The weighting functions for the feedback H∞ control design are chosen as

Aw = diag(−500, −100, −200)

Bw = diag(9, 9, 9)

Cw = diag(−1, −1, −1)

Dw = diag(0.1, 0.1, 0.1)

Wu = diag(0.2, 0.05, 0.2)

The solution to the H∞ optimization problem is a eleventh-order linear system for the controller K(s). Through model reduction techniques by investigating the dominant singular values,
a seventh-order controller is used in simulation. Saturation is also imposed to restrict the control
inputs to the feasible range of actuators.
Figures 8.10 and 8.11 show the simulation results of the tracking performance of three geometric parameters and the control efforts of three control inputs, respectively, for tracking a steady
spiral trajectory with constant geometric parameters and Fig. 8.12 shows the tracking trajectory
in the 3D view under the proposed 2-DOF controller. From the simulation results, we see that
both PI controller and feedforward inverse mapping controller are able to stabilize the system to
the desired steady-spiral trajectory. However, both have convergence times between 30 and 40
seconds. There is also noticeable steady-state tracking error with the PI controller. With the proposed 2-DOF controller, the system tracking performance is improved significantly. Convergence

152

time decreases to less than 10 seconds with smaller steady-state errors. Meanwhile, the 2DOF
controller uses more control effort than the open-loop feedforward strategy. Figs. 8.13 and 8.14
show the simulation results of the tracking performance of three geometric parameters and the control efforts of three control inputs, respectively, when the reference velocity changes as a sinusoid
function with respect to time while the curvature and torsion are kept constant. There is a large
time/phase delay in the tracking of velocity for the feedforward control, which is expected due to
the observed slow convergence speed. Besides the time delay, there is also significant tracking error with the PI control. But with the 2-DOF controller, the tracking performance shows significant
improvement in terms of the time/phase delay and the tracking error. Besides, in Figs. 8.13a and
8.13b, the variables fluctuate even though the reference is a constant. This shows coupling among
control inputs on the spiral geometric features.

153

0.4
Reference to Track
PI Control

0.3

Inversion control

κ

2−DOF Control
0.2

0.1

0

0

5

10

15

20

25

30

35

40

45

50

t(s)

(a)
0.4
Reference to Track
PI Control

0.3

Inversion control

τ

2−DOF Control
0.2

0.1

0

0

5

10

15

20

25

30

35

40

45

50

t(s)

(b)
0.2

V (m/s)

0.18
0.16
0.14

Reference to Track
PI Control

0.12

Inversion control
2−DOF Control

0.1

0

5

10

15

20

25

30

35

40

45

50

t(s)

(c)
Figure 8.10: The simulation results of the geometric parameters when tracking a steady spiral
trajectory. (a) curvature; (b) torsion; (c) velocity.

154

50
PI Control

rp (mm)

40

Inversion control
2−DOF Control

30
20
10
0

0

5

10

15

t(s)

(a)
50
PI Control

δ (°)

40

Inversion control
2−DOF Control

30
20
10
0

0

5

10

15

t(s)

(b)
90

m0 (g)

80
70
60

PI Control
Inversion control

50
40

2−DOF Control
0

5

10

15

t(s)

(c)
Figure 8.11: The simulation results of control inputs on when tracking a steady spiral trajectory.
(a) displacement of movable mass; (b) tail angle; (c) net buoyancy.

155

end point
start point
0

Z (m)

−10
−20
−30
−40
0

−2

−4

−6

Y (m)

−4

−8

−2

0

2

X (m)

Figure 8.12: 3D trajectory when tracking a steady spiral under the 2-DOF controller.

156

0.17
0.16

κ

0.15
0.14
Reference to Track
0.13

PI Control
Inversion control

0.12

2−DOF Control
0.11
50

60

70

80

90

100

110

120

130

140

150

110

120

130

140

150

t(s)

(a)
0.115
Reference to Track
0.11

PI Control
Inversion control

0.105

τ

2−DOF Control
0.1
0.095
0.09
0.085
50

60

70

80

90

100

t(s)

(b)
0.22
Reference to Track

V (m/s)

0.21

PI Control
Inversion control

0.2

2−DOF Control
0.19
0.18
0.17
0.16
50

60

70

80

90

100

110

120

130

140

150

t(s)

(c)
Figure 8.13: The simulation results of the geometric parameters when the reference velocity
changes as a sinusoid function with respect to time while curvature and torsion are kept constant.
(a) curvature; (b) torsion; (c) velocity.

157

19
PI Control

rp (mm)

18

Inversion control
2−DOF Control

17
16
15
14
50

60

70

80

90

100

110

120

130

140

150

110

120

130

140

150

110

120

130

140

150

t(s)

(a)
25
PI Control

24

Inversion control

δ (°)

23

2−DOF Control

22
21
20
19
50

60

70

80

90

100

t(s)

(b)
90
85

m0 (g)

80
75
PI Control

70

Inversion control
65
60
50

2−DOF Control
60

70

80

90

100

t(s)

(c)
Figure 8.14: The simulation results of control inputs on when the reference velocity changes as a
sinusoid function with respect to time while curvature and torsion are kept constant. (a) displacement of movable mass; (b) tail angle; (c) net buoyancy.

158

Chapter 9
Field Test Results for Environmental
Monitoring
In this chapter, we present some preliminary field testing results for oil spill detection in the Kalamazoo River, Michigan, and algae bloom monitoring in the Wintergreen Lake, Michigan, where
basic functions of gliding robotic fish, the swimming, the gliding, and the spiraling are tested in
the field. It is found that the crude oil sensor readings were different at the three selected locations
along Kalamazoo River near the spill site, with slightly higher values at the downstream sites.
For the algae concentration, the sampling results were consistent with readings from a traditional
hand-held device [84].

9.1 Kalamazoo River Test
“Grace” was deployed in Kalamazoo River in November, 2012, to detect the crude oil content
near the 2010 Enbridge oil spill site near Marshall, Michigan. A Turner Designs Cyclops-7 crude
oil sensor was used to sample the water. We tested three locations along the Kalamazoo River
(Fig. 9.1). The spot A was downstream of the oil spill spot which was an open river area (Fig. 9.2);
the spot B was upstream which was in the woods (Figs. 9.3); and the spot C was also upstream
which was under a bridge (Fig. 9.4). “Grace” swam in the Kalamazoo river by flapping its tail
fin, collecting data in the selected locations. The readings from the sensor are recorded onboard
159

flow direction

oil spill location

B

A

C

Kalamazoo River

Figure 9.1: Three selected sampling locations illustrated in Google Map.
Table 9.1: Locations and sensor readings in Kalamazoo River test.
Spot
GPS Coordinate
Sensor Reading
A
(42.258655, -84.99888333)
3.33 V
B
(42.26156167, -84.953565)
3.04 V
C
(42.2616, -84.85534667)
2.87 V
and transmitted wirelessly to the base station onshore. The GPS position and the average sensor
reading at each spot are shown in Table 9.1.
The results show that the sensor readings were slightly higher at the downstream sites of the
spill spot than the upstream. There is a relationship between the sensor’s output (in voltage) and
the real concentration of crude oil. In this dissertation, we provide only the raw readings from the
crude oil sensor. The detection spectrum of the sensor is relatively wide so that it may not precisely
reflect the actual crude oil concentration.

160

Figure 9.2: Swimming trajectory of gliding robotic fish in an open water area at testing spot A.

Figure 9.3: Swimming trajectory of gliding robotic fish in the woods at testing spot B.

161

Figure 9.4: Swimming trajectory of gliding robotic fish under a bridge at testing spot C.

9.2 Wintergreen Lake Test
Gliding robotic fish “Grace” was also used to sample the Wintergreen Lake, Michigan for the bluegreen algae (cyanobacteria) concentration in July, 2013. The sensor we used was Turner Designs
Cyclops-7 Freshwater Blue-Green Algae Sensor. The sensor readings in this field test are compared
with those from a manually deployed profiler with the blue-green algae sensor (Hydrolab) that has
been traditionally used in the sampling of harmful algae.
Surface swimming was tested first at Wintergreen Lake using “Grace” (Fig. 9.5), with sensor
readings recorded for selected locations. Table 9.2 shows the GPS position and the sensor outputs
of algae concentration at those points.
Since algae concentration varies at different depths, we ran a number of steady glide and steady
spiral motions using the robot to detect the algae bloom in the three-dimensional water space.
Figs. 9.6–9.9 show the results of steady glide tests. The start (submerging) and the end (surfacing)

162

E

F

D

C

A
B

Figure 9.5: Selected sampling points in the Wintergreen Lake with surface swimming.

Table 9.2: Locations and sensor readings in Wintergreen Lake test.
Spot Sensor Reading
A
0.252979 V
B
0.336768 V
C
0.289233 V
D
0.263452 V
E
0.269897 V
F
0.422168 V

Hydrolab Chl
2.19
2.45
2.33
3.06
2.60
3.89

163

GPS Coordinate
(42.3987, -85.38334)
(42.39851, -85.38359)
(42.39874, -85.38367)
(42.39904, -85.38356)
(42.39919, -85.3831)
(42.39938, -85.38271)

3.5
3

depth (m)

2.5
2
1.5
1
0.5
0
0

50

100

150

t(s)
Figure 9.6: Depth trajectory of steady glide when sampling water in Wintergreen Lake, Michigan.
points are (42.398228, -85.384644), and (42.398361, -85.384636), respectively, with a travel distance of 14 meters. The readings of Chlorophyll from the Hydrolab device at this location are 3.11
at 1.26 m, 5.69 at 2.45 m, and 8.73 at 3.55 m. Figs 9.10–9.13 show the results for the steady spiral
test with the spiral center at (42.398338, -85.384727). The readings of Chlorophyll from the Hydrolab device at this location are 2.89 at 0.97 m, 6.59 at 3.48 m, and 10.26 at 4.05 m. The readings
from the robot were found to be consistent in general with those from the Hydrolab device.

164

1

sensor reading (V)

0.9
0.8
0.7
0.6
0.5
0.4
0

50

100

150

t(s)
Figure 9.7: Sensor readings of steady glide when sampling water in Wintergreen Lake, Michigan.

180
160
140

°

yaw ( )

120
100
80
60
40
20
0
0

50

100

150

t(s)
Figure 9.8: Yaw angle trajectory when sampling water in Wintergreen Lake, Michigan.

165

30

°

pitch ( )

20

10

0

−10

−20
0

50

100

150

t(s)
Figure 9.9: Pitch angle trajectory when sampling water in Wintergreen Lake, Michigan.

4
3.5

depth (m)

3
2.5
2
1.5
1
0.5
0
0

50

100

150

200

t(s)
Figure 9.10: Depth trajectory of steady spiral when sampling water in Wintergreen Lake, Michigan.

166

1.6

sensor reading (V)

1.4
1.2
1
0.8
0.6
0.4
0.2
0

50

100

150

200

t(s)
Figure 9.11: Sensor readings of steady spiral when sampling water in Wintergreen Lake, Michigan.

200
150

50

°

yaw ( )

100

0
−50
−100
−150
−200
0

50

100

150

200

t(s)
Figure 9.12: Yaw angle trajectory of steady spiral when sampling water in Wintergreen Lake,
Michigan.

167

30
20

°

pitch ( )

10
0
−10
−20
−30
0

50

100

150

200

t(s)
Figure 9.13: Pitch angle trajectory of steady spiral when sampling water in Wintergreen Lake,
Michigan.

168

Chapter 10
Conclusions and Future Work

10.1 Conclusions
We reported a new type of underwater robots, gliding robotic fish, designed for aquatic environmental monitoring. Combing design features of buoyancy-driven propulsion of underwater gliders
and tail-fin actuation of robotic fish, the gliding robotic fish shows a great potential in shallow
water sampling with energy efficiency and high maneuverability.
We first introduced two lab-developed prototypes, a miniature underwater glider and a gliding
robotic fish “Grace”. The actuation systems were discussed including the buoyancy system, mass
distribution system and tail-fin system, to explain the locomotion mechanism of the robot. The
sensor system was introduced, including the inertial measurement units and environmental sensing
units. The mechanical design of the robot was also discussed.
Dynamic model of gliding robotic fish was derived and further reduced according to two special steady-state motions, the sagittal-plane glide and 3D steady sprial. Solutions of the robot’s
gliding path were provided for both motions. Particularly, for steady spiral, we adopted Newton’s
method and numerically explored the basins of attraction for this recursive algorithm. Experiments
were conducted on the miniature underwater glider prototype. Comparison between model prediction and experimental results of the glide path was carried out to validate the derived model and
proposed recursive solving method.
Stabilization problems were studied, for both sagittal-plane and lateral motions. For the sagittal169

plane gliding path stabilization, a passivity-based controller using a whale-like tail fin was proposed with partial state feedback. A nonlinear observer was designed to estimate the velocityrelated system states in controller implementation. Both simulation and experiments were conducted to evaluate the effectiveness of the designed controller and observer. For lateral motion
stabilization, a sliding mode control scheme was adopted with a fish-like tail as the system control
input, which depended only on the pitch angle information. Both simulation and experimental
results were presented.
A 2-DOF control strategy was proposed for three-dimensional curve tracking problem for
gliding robotic fish, by following the trajectory characteristic parameters (curvature, torsion and
speed). We investigated the steady spiral motion and its geometric characteristics. A feedforward
controller was designed first via inverse mapping of steady spiral motion. A 2-DOF control design
was then proposed, which included an inverse mapping feedforward controller and a robust H∞
feedback controller, designed based on the linearized model. Simulation was conducted with comparison to PI control and pure feedforward inverse mapping control. The simulation results were
presented to verify the effectiveness of the proposed 2-DOF control design.
At last, full functions of the developed gliding robotic fish prototype were tested in field. The
results of field tests, including crude oil detection in Kalamazoo River and algae sampling in Wintergreen Lake, were presented.

10.2 Future Work
First, in both the lab/pool experiments and the field tests, significant drifting of the robot was
observed that was caused by the currents and waves of its aquatic environment, e.g., the water jets
shooting from the side wall of the swimming pool. Thus an environmental current flow estimation

170

or prediction is desirable in disturbance rejection and precise locomotion. The information of the
surrounding environments will benefit the controller design and the path following.
Second, a comprehensive path planning and path following control strategy is needed in water
sampling tasks. Take the algae sampling in Wintergreen Lake as an example. It is of interest to
develop adaptive sampling strategies and explore the combination of elementary energy-efficient
motions (steady gliding and steady spiraling), to cover the whole lake with high sampling resolution and and low energy consumption.

171

BIBLIOGRAPHY

172

BIBLIOGRAPHY
[1] C. C. Ericksen, T. J. Osse, R. D. Light, T. Wen, T. W. Lehman, P. L. Sabin, J. W. Ballard, and
A. M. Chiodi, “Seaglider: A long-range autonomous underwater vehicle for oceanographic
research,” IEEE Journal of Oceanic Engineering, vol. 26, pp. 424–436, 2001.
[2] J. Sherman, R. E. Davis, W. B. Owens, and J. Valdes, “The autonomous underwater glider
“Spray”,” IEEE Journal of Oceanic Engineering, vol. 26, pp. 437–446, 2001.
[3] D. C. Webb, P. J. Simonetti, and C. P. Jones, “SLOCUM, an underwater glider propelled by
environmental energy,” IEEE Journal of Oceanic Engineering, vol. 26, no. 4, pp. 447–452,
2001.
[4] S. American, “Researchers spawn a new breed of robotic fish.” [Online]. Available:
http://www.scientificamerican.com/article/robotic-fish/
[5] R. P. Schwarzenbach, T. Egli, T. B. Hofstetter, U. Von Gunten, and B. Wehrli, “Global
water pollution and human health,” Annual Review of Environment and Resources, vol. 35,
pp. 109–136, 2010.
[6] N. F. Gray et al., Water technology: an introduction for environmental scientists and engineers. IWA Publishing, 2010, no. Ed. 3.
[7] R. Rajagopalan, Environmental studies: from crisis to cure. Oxford University Press, 2011.
[8] United States National Oceanographic and Atmospheric Administration, “Erma gulf
response.” [Online]. Available: http://response.restoration.noaa.gov/maps-and-spatial-data/
environmental-response-management-application-erma/erma-gulf-response.html
[9] Office of the Maritime Administrator, “Deepwater horizon marine casualty investigation
report,” August 2011, Retrieved on February 25, 2013.
[10] USEPA, “EPA Response to Enbridge Spill in Michigan.” [Online]. Available: http:
//www.epa.gov/enbridgespill/documents.html
[11] Wikipedia, “Xingang port oil spill.” [Online]. Available: http://en.wikipedia.org/wiki/
Xingang Port oil spill

173

[12] United Nations Environment Programme, “United Nations report on Nigeria oil
spills.” [Online]. Available: http://www.shell.com.ng/environment-society/our-response/
unep-response.html
[13] Wikipedia, “Oil Spill.” [Online]. Available: http://en.wikipedia.org/wiki/Oil spill
[14] S. A. Ruberg, R. W. Muzzi, S. B. Brandt, J. C. Lane, T. C. Miller, J. J. Gray, S. A. Constant,
and E. J. Downing, “A wireless internet-based observatory: The Real-time Coastal Observation Network (ReCON),” in Proceedings of Marine Technology Society/IEEE Oceans 2007
Conference, 2007, 6 pp.
[15] H. Hongji, A. Kaneko, and K. Kawatate, “Self-governing profiling system,” Cont. Shelf
Res., vol. 7, pp. 1257–1265, 1987.
[16] R. E. Davis, D. C. Webb, L. A. Regier, and J. Dufour, “The autonomous Lagrangian circulation explorer (ALACE),” Journal of Atmospheric and Oceanic Technology, vol. 9, pp.
264–285, 1992.
[17] R. A. Luettich, “PSWIMS, a profiling instrument system for remote physical and chemical
measurements in shallow-water,” Estuaries, vol. 16, pp. 190–197, 1993.
[18] K. W. Doherty, D. E. Frye, S. P. Liberatore, and J. M. Toole, “A moored profiling instrument,” J Atmos Ocean Technol., vol. 16, pp. 1816–1829, 1999.
[19] J. V. Reynolds-Fleming, J. G. Fleming, and R. A. Luettich, “Portable autonomous vertical
profiler for estuarine applications,” Estuaries, vol. 25, pp. 142–147, 2002.
[20] J. B. Newman and B. H. Robison, “Development of a dedicated ROV for ocean science,”
Marine Technology Society Journal, vol. 26, pp. 46–53, 1994.
[21] J. Lee, M. Roh, K. Kim, and D. Lee, “Design of autonomous underwater vehicles for cage
aquafarms,” in Proceedings of the 2007 IEEE Intelligent Vehicles Symposium, Istanbul,
Turkey, 2007, pp. 938–943.
[22] ScienceDaily, “Rensselaer researchers experiment with solar underwater robots,” December
2004. [Online]. Available: http://www.sciencedaily.com/releases/2004/12/041212081548.
htm
[23] G. S. Sukhatme, A. Dhariwal, B. Zhang, C. Oberg, B. Stauffer, and D. A. Caron, “Design and development of a wireless robotic networked aquatic microbial observing system,”
Environmental Engineering Science, vol. 24, no. 2, pp. 205–215, 2007.
174

[24] D. R. Yoerger, A. M. Bradley, M. Jakuba, C. R. German, T. Shank, and M. Tivey, “Autonomous and remotely operated vehicle technology for hydrothermal vent discovery, exploration, and sampling,” Oceanography, vol. 20, no. 1, pp. 152–161, 2007.
[25] J. R. Higinbotham, J. R. Moisan, M. Schirtzinger, C.and Linkswiler, and P. Yungel, J.and Orton, “Update on the development and testing of a new long duration solar powered autonomous surface vehicle,” in Proceedings of the 2008 OCEANS Conference, Quebec City,
QC, 2008, 10 pp.
[26] M. Stealey, A. Singh, M. Batalin, B. Jordan, and W. Kaiser, “NIMS-AQ: A novel system
for autonomous sensing of aquatic environments,” in Proceedings of the 2008 IEEE International Conference on Robotics and Automation, Pasadena, CA, 2008, 621-628.
[27] A. Goldstein and S. Bentley, “Use of highly portable micro-sized remotely operated vehicles
for environmental monitoring and mapping,” in Proceedings of OCEANS 2010 Conference,
Seattle, WA, 2010, 6 pp.
[28] J. Das, H. Heidarsson, A. Pereira, F. Arrichiello, I. Cetnic, L. Darjany, M.-E. Garneau,
M. Howard, C. Oberg, M. Ragan, E. Seubert, E. Smith, B. Stauffer, A. Schnetzer, G. ToroFarmer, D. Caron, B. Jones, and G. Sukhatme, “USC CINAPS builds bridges,” IEEE Robot.
Automat. Mag., vol. 17, no. 1, pp. 20–30, 2010.
[29] R. N. Smith, Y. Chao, P. P. Li, D. A. Caron, B. H. Jones, and G. S. Sukhatme, “Planning and
implementing trajectories for autonomous underwater vehicles to track evolving ocean processes based on predictions from a regional ocean model,” International Journal of Robotics
Research, vol. 29, no. 12, pp. 1475–1497, 2010.
[30] F. S. Hover, R. M. Eustice, A. Kim, B. Englot, H. Johannsson, M. Kaess, and J. J. Leonard,
“Advanced perception, navigation and planning for autonomous in-water ship hull inspection,” The International Journal of Robotics Research, vol. 31, no. 12, pp. 1445–1464, 2012.
[31] M. Dunbabin and A. Grinham, “Experimental evaluation of an Autonomous Surface Vehicle
for water quality and greenhouse gas emission monitoring,” in Proceedings of the 2010
IEEE International Conference on Robotics and Automation, Anchorage, AK, 2010, pp.
5268–5274.
[32] L. Elkins, D. Sellers, and W. R. Monach, “The Autonomous Maritime Navigation (AMN)
project: Field tests, autonomous and cooperative behaviors, data fusion, sensors, and vehicles,” J. Field Robot., vol. 27, no. 6, pp. 790–818, 2010.
[33] S. J. M. Hughes, D. O. B. Jones, C.Hauton, A. R. Gates, and L. E.Hawkins, “An assessment
of drilling disturbance on echinus acutus var. norvegicus based on in-situ observations and
175

experiments using a remotely operated vehicle (ROV),” J. Exp. Marine Biol. Ecol., vol. 395,
no. 1-2, pp. 37–47, 2010.
[34] S. B. Williams, O. Pizarro, J. Webster, R. Beaman, I. Mahon, M. Johnson-Roberson, and
T. Bridge, “AUV-assisted surveying of drowned reefs on the shelf edge of the Great Barrier
Reef, Australia,” J. Field Robot., vol. 27, no. 5, pp. 675–697, 2010.
[35] G. Hitz, F. Polerleau, M. Gameau, C. Pradalier, T. Posch, J. Pernthaler, and R. Y. Siegwart,
“Autonomous inland water monitoring: Design and application of a surface vehicle,” IEEE
Robot. Automat. Mag, vol. 19, no. 1, pp. 62–72, 2012.
[36] S. B. Williams, O. R. Pizarro, M. V. Jakuba, C. R. Johnson, N. S. Barrett, R. C. Babcock,
G. A. Kendrick, P. D. Steinberg, A. J. Heyward, P. J. Doherty, I. Mahon, M. JohnsonRoberson, D. Steinberg, and A. Friedman, “Monitoring of benthic reference sites: Using an
autonomous underwater vehicle,” IEEE Robot. Automat. Mag., vol. 19, no. 1, pp. 73–84,
2012.
[37] P. R. Bandyopadhyay, “Trends in biorobotic autonomous undersea vehicles,” IEEE Journal
of Oceanic Engineering, vol. 30, no. 1, pp. 109–139, 2005.
[38] B. Bingham, B. Foley, H. Singh, R. Camilli, K. Delaporta, R. Eustice, A. Mallios, D. Mindell, C. Roman, and D. Sakellariou, “Robotic tools for deep water archaeology: Surveying
an ancient shipwreck with an autonomous underwater vehicle,” J. Field Robot., vol. 27,
no. 6, pp. 702–717, 2010.
[39] L. L. Whitcomb, M. V. Jakuba, J. C. Kinsey, S. C. Martin, S. E. Webster, J. C. Howland,
C. L. Taylor, D. Gomez-Ibanez, and D. Yoerger, “Navigation and control of the Nereus
hybrid underwater vehicle for global ocean science to 10,903 m depth: Preliminary results,”
in Proceedings of the 2010 IEEE International Conference on Robotics and Automation,
Anchorage, AK, 2010, 594-600.
[40] G. V. Lauder and E. G. Drucker, “Morphology and experimental hydrodynamics of fish fin
control surfaces,” IEEE Journal of Oceanic Engineering, vol. 29, no. 3, pp. 556–571, 2004.
[41] F. E. Fish and G. V. Lauder, “Passive and active flow control by swimming fishes and mammals,” Annual Review of Fluid Mechanics, vol. 38, pp. 193–224, 2006.
[42] A. Punning, M. Anton, M. Kruusmaa, and A. Aabloo, “A biologically inspired ray-like underwater robot with electroactive polymer pectoral fins,” in International IEEE Conference
on Mechatronics and Robotics, vol. 2004, 2004, pp. 241–245.

176

[43] R. Pfeifer, M. Lungarella, and F. Iida, “Self-organization, embodiment, and biologically
inspired robotics,” science, vol. 318, no. 5853, pp. 1088–1093, 2007.
[44] R. Altendorfer, N. Moore, H. Komsuoglu, M. Buehler, H. Brown Jr, D. McMordie,
U. Saranli, R. Full, and D. E. Koditschek, “Rhex: A biologically inspired hexapod runner,”
Autonomous Robots, vol. 11, no. 3, pp. 207–213, 2001.
[45] K. Melsaac and J. P. Ostrowski, “A geometric approach to anguilliform locomotion: modelling of an underwater eel robot,” in Robotics and Automation, 1999. Proceedings. 1999
IEEE International Conference on, vol. 4. IEEE, 1999, pp. 2843–2848.
[46] E. Garcia, M. A. Jimenez, P. G. De Santos, and M. Armada, “The evolution of robotics
research,” Robotics & Automation Magazine, IEEE, vol. 14, no. 1, pp. 90–103, 2007.
[47] D. L. Rudnick, R. E. Davis, C. C. Eriksen, D. M. Fratantoni, and M. J. Perry, “Underwater
gliders for ocean research,” Marine Technology Society Jouranl, vol. 38, no. 2, pp. 73–84,
2004.
[48] M. S. Triantafyllou and G. S. Triantafyllou, “An efficient swimming machine,” Scientific
America, vol. 273, no. 3, pp. 64–70, 1995.
[49] N. Kato, “Control performance in the horizontal plane of a fish robot with mechanical pectoral fins,” IEEE Journal of Oceanic Engineering, vol. 25, no. 1, pp. 121–129, 2000.
[50] J. Anderson and N. Chhabra, “Maneuvering and stability performance of a robotic tuna,”
Integrative and Comparative Biology, vol. 42, no. 1, pp. 118–126, 2002.
[51] P. R. Bandyopadhyay, “Maneuvering hydrodynamics of fish and small underwater vehicles,”
Integrative and Comparative Biology, vol. 42, pp. 102–117, 2002.
[52] J. Yu, M. Tan, S. Wang, and E. Chen, “Development of a biomimetic robotic fish and its control algorithm,” IEEE Transactions on Systems, Man, and Cybernetics, Part B: Cybernetics,
vol. 34, no. 4, pp. 1798–1810, 2004.
[53] J. H. Long, A. C. Lammert, C. A. Pell, M. Kemp, J. A. Strother, H. C. Crenshaw, and
M. J. McHenry, “A navigational primitive: Biorobotic implementation of cycloptic helical
klinotaxis in planar motion,” IEEE Journal of Oceanic Engineering, vol. 29, no. 3, pp. 795–
806, 2004.
[54] P. V. Alvarado and K. Youcef-Toumi, “Design of machines with compliant bodies for
biomimetic locomotion in liquid environments,” Journal of Dynamic Systems, Measurement, and Control, vol. 128, pp. 3–13, 2006.
177

[55] J. Liu and H. Hu, “Biologically inspired behavior design for autonomous robotic fish,” International Journal of Automation and Computing, vol. 4, pp. 336–347, 2006.
[56] M. Epstein, J. E. Colgate, and M. A. MacIver, “Generating thrust with a biologicallyinspired robotic ribbon fin,” in Proceedings of the 2006 IEEE/RSJ International Conference
on Intelligent Robots and Systems, Beijing, China, 2006, pp. 2412–2417.
[57] M. Krieg and K. Mohseni, “Dynamic modeling and control of biologically inspired vortex
ring thrusters for underwater robot locomotion,” IEEE Transactions on Robotics, vol. 26,
no. 3, pp. 542–554, 2010.
[58] K. A. Morgansen, B. I. Triplett, and D. J. Klein, “Geometric methods for modeling and control of free-swimming fin-actuated underwater vehicles,” IEEE Transactions on Robotics,
vol. 23, no. 6, pp. 1184–1199, 2007.
[59] G. V. Lauder, E. J. Anderson, J. Tangorra, and P. G. A. Madden, “Fish biorobotics: Kinematics and hydrodynamics of self-propulsion,” Journal of Experimental Biology, vol. 210,
pp. 2767–2780, 2007.
[60] P. Kodati, J. Hinkle, A. Winn, and X. Deng, “Microautonomous robotic ostraciiform
(MARCO): Hydrodynamics, design, and fabrication,” IEEE Transactions on Robotics,
vol. 24, no. 1, pp. 105–117, 2008.
[61] K. H. Low and C. W. Chong, “Parametric study of the swimming performance of a fish
robot propelled by a flexible caudal fin,” Bioinsp. Biomim., vol. 5, p. 046002, 2010.
[62] Z. C. S. Shatara and X. Tan, “Modeling of biomimetic robotic fish propelled by an ionic
polymer-metal composite caudal fin,” IEEE/ASME Transactions on Mechatronics, vol. 15,
no. 3, pp. 448–459, 2010.
[63] M. Aureli, V. Kopman, and M. Porfiri, “Free-locomotion of underwater vehicles actuated
by ionic polymer metal composites,” IEEE/ASME Transactions on Mechatronics, vol. 15,
no. 4, pp. 603–614, 2010.
[64] J. Yu, R. Ding, Q. Yang, M. Tan, W. Wang, and J. Zhang, “On a bio-inspired amphibious
robot capable of multimodal motion,” IEEE/ASME Transactions on Mechatronics, vol. PP,
no. 99, pp. 1–10, 2011.
[65] P. C. Strefling, A. M. Hellum, and R. Mukherjee, “Modeling, simulation, and performance of a synergistically propelled ichthyoid,” IEEE/ASME Transactions on Mechatronics,
vol. 17, no. 1, pp. 36–45, 2012.
178

[66] F. Liu, K.-M. Lee, and C.-J. Yang, “Hydrodynamics of an undulating fin for a wave-like
locomotion system design,” Mechatronics, IEEE/ASME Transactions on, vol. 17, no. 3, pp.
554–562, 2012.
[67] C. Zhou and K. H. Low, “Design and locomotion control of a biomimetic underwater vehicle
with fin propulsion,” IEEE/ASME Transactions on Mechatronics, vol. 17, no. 1, pp. 25–35,
2012.
[68] S. Marras and M. Porfiri, “Fish and robots swimming together: attraction towards the robot
demands biomimetic locomotion,” Journal of The Royal Society Interface, vol. 9, no. 73,
pp. 1856–1868, 2012.
[69] M. Schwartz, Encyclopedia of coastal science.

Kluwer Academic Pub, 2005, vol. 24.

[70] R. E. Davis, N. E. Leonard, and D. M. Fratantoni, “Routing strategies for underwater gliders,” Deep Sea Research Part II: Topical Studies in Oceanography, vol. 56, no. 3-5, pp. 173
– 187, 2009.
[71] D. Blidberg, “The development of autonomous underwater vehicles (AUV); a brief summary,” in IEEE ICRA, Seoul, Korea, 2001.
[72] N. Mahmoudian, “Efficient motion planning and control for underwater gliders,” Ph.D. dissertation, Virginia Polytechnic Institute and State University, 2009.
[73] N. E. Leonard and J. G. Graver, “Model-based feedback control of autonomous underwater
gliders,” IEEE Journal of Oceanic Engineering, vol. 26, no. 4, pp. 633–645, 2001.
[74] N. Mahmoudian and C. A. Woolsey, “Underwater glider motion control,” in Proceedings of
IEEE Conference on Decision and Control, Cancun, Mexico, 2008, pp. 552–557.
[75] P. Bhatta, “Nonlinear stability and control of gliding vehicles,” Ph.D. dissertation, Princeton
University, 2006.
[76] P. Bhatta and N. Leonard, “Nonlinear gliding stability and control for vehicles with hydrodynamic forcing,” Automatica, vol. 44, no. 5, pp. 1240–1250, 2008.
[77] J. G. Graver, “Underwater gliders: Dynamics, control, and design,” Ph.D. dissertation,
Princeton University, 2005.

179

[78] F. Zhang, X. Tan, and H. K.Khalil, “Passivity-based controller design for stabilization of
underwater gliders,” in Proceedings of the 2012 American Control Conference, Montreal,
Canada, 2012, pp. 5408–5413.
[79] N. Mahmoudian, J. Geisbert, and C. Woolsey, “Approximate analytical turning conditions
for underwater gliders: Implications for motion control and path planning,” IEEE Journal
of Oceanic Engineering, vol. 35, no. 1, pp. 131–143, 2010.
[80] S. Zhang, J. Yu, A. Zhang, and F. Zhang, “Spiraling motion of underwater gliders: Modeling, analysis, and experimental results,” Ocean Engineering, vol. 60, pp. 1–13, 2013.
[81] F. Zhang, F. Zhang, and X. Tan, “Steady spiraling motion of gliding robotic fish,” in 2012
IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS). IEEE, 2012,
pp. 1754–1759.
[82] F. Zhang, J. Thon, C. Thon, and X. Tan, “Miniature underwater glider: Design, modeling,
and experimental results,” in Proceedings of the 2012 IEEE International Conference on
Robotics and Automation, St. Paul, MN, 2012, pp. 4904–4910.
[83] ——, “Miniature underwater glider: Design and experimental results,” IEEE/ASME Transactions on Mechatronics, 2014, vol. 19, no. 1, pp. 394-399.
[84] F. Zhang, J. Wang, J. Thon, C. Thon, E. Litchman, and X. Tan, “Gliding robotic fish for
mobile sampling of aquatic environments,” in Proceedings of 11th IEEE International Conference on Networking, Sensing and Control, Miami, FL, 2014, accepted.
[85] J. Milgram, “Strip theory for underwater vehicles in water of finite depth,” Journal of Engineering Mathematics, vol. 58, no. 1, pp. 31–50, 2007.
[86] R. L. Panton, Incompressible Flow. New York: Wiley, 2005.
[87] D. T. Greenwood, Principles of Dynamics. Englewood Cliffs, NJ.: Prentice Hall, 1988.
[88] J. D. Anderson, Aircraft Performance and Design. New York: McGraw-Hill, 1998.
[89] F. O. Smetana, Introductory aerodynamics and hydrodynamics of wings and bodies: a
software-based approach. AIAA, 1997, vol. 1.
[90] R. Von Mises, Theory of flight. Dover Pubns, 1959.
[91] I. G. Currie, Fundamental Mechanics of Fluids. New York: McGraw-Hill, 1974.
180

[92] K. Fujii, “Progress and future prospects of CFD in aerospace wind tunnel and beyond,” Prog
Aerospace Sci, vol. 41, pp. 455–470, 2005.
[93] The Autonomous Underwater Vehicle Applications Center (AUVAC), “AUVs database in
AUVAC.” [Online]. Available: http://auvac.org/explore-database/simple-search
[94] S. Zhang, J. Yu, A. Zhang, and F. Zhang, “Steady three dimensional gliding motion of an
underwater glider,” in 2011 IEEE International Conference on Robotics and Automation
(ICRA). IEEE, 2011, pp. 2392–2397.
[95] C. Kelley, Solving Nonlinear Equations with Newton’s Method.
Mathematics, 2003.

Society for Industrial

[96] H. Khalil, Nonlinear systems. Upper Saddle River, N.J.: Prentice Hall, 2002.
[97] P. Kokotovi´c, H. Khalil, and J. O’Reilly, Singular perturbation methods in control: analysis
and design. Society for Industrial Mathematics, 1999.
[98] F. Zhang and X. Tan, “Nonlinear observer design for stabilization of gliding robotic fish,”
in Proceedings of the 2014 American Control Conference, Portland, OR, 2014, accepted.
[99] ——, “Gliding robotic fish and its tail-enabled yaw motion stabilization using sliding mode
control,” in Proceedings of the 2013 ASME Dynamic Systems and Control Conference, Palo
Alto, CA, 2013, paper DSCC2013-4015 (10 pp).
[100] F. W. Warner, Foundations of differentiable manifolds and Lie groups.
vol. 94.

Springer, 1983,

[101] F. Zhang and X. Tan, “Three-dimensional spiral tracking control for gliding robotic fish,”
in Proceedings of the 53rd IEEE Conference on Decision and Control, Los Angeles, CA,
2014, accepted.
[102] A. Pressley, Elementary differential geometry. Springer, 2010.

181