ROBOTIC FISH: DEVELOPMENT, MODELING, AND APPLICATION TO MOBILE
SENSING
By
Jianxun Wang

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

ABSTRACT
ROBOTIC FISH: DEVELOPMENT, MODELING, AND APPLICATION TO MOBILE
SENSING
By
Jianxun Wang
Robotic fish are underwater robots that emulate locomotion of live fish through actuated fin
and/or body movements. They are of increasing interest due to their potential applications such as
aquatic environmental monitoring and robot-animal interactions.
In this work, several bio-inspired robotic fish prototypes have been developed that make use
of periodic tail motions. A dynamic model for a tail-actuated robotic fish is presented by merging
rigid-body dynamics with Lighthill’s large-amplitude elongated-body theory. The model is validated with extensive experiments conducted on a robotic fish prototype. The role of incorporating
the body motion in evaluating the tail-generated hydrodynamic forces is assessed, which shows
that ignoring the body motion (as often done in the literature) results in significant overestimate of
the thrust force and robot speed. By exploiting the strong correlation between the angle of attack
and the tail-beat bias, a computationally efficient approach is further proposed to adapt the drag
coefficients of the robotic fish.
It has been recognized that the flexibility of the body and fin structures has a pronounced
impact on the swimming performance of biological and robotic fish. To analyze and utilize this
trait, a novel dynamic model is developed for a robotic fish propelled by a flexible tail actuated at
the base. The tail is modeled with multiple rigid segments connected in series through rotational
springs and dampers. For comparison, a model using linear beam theory is created to capture the
beam dynamics. Experimental result show that the two models have almost identical predictions
when the tail undergoes small deformation, but only the proposed multi-segment model matches

the experimental measurement closely for all tail motions.
Motivated by the need for system analysis and efficient control of robotic fish, averaging of
robots’ dynamics is of interest. For dynamic models of robotic fish, however, classical or geometric
averaging typically cannot produce an average model that is accurate and the in the meantime
amenable to analysis or control design. In this work, a novel averaging approach for tail-actuated
robotic fish dynamics is proposed. The approach consists of scaling the force and moment terms
and then conducting classical averaging. Numerical investigation reveals that the scaling function
for the force terms is a constant independent of tail-beat patterns, while the scaling function for the
moment term depends linearly on the tail-beat bias. Existence and local stability of the equilibria
for the average model are further analyzed. Finally, as an illustration of the utility of the average
model, a semi-analytical framework is presented for obtaining steady turning parameters.
Sampling and reconstruction of a physical field using mobile sensor networks have recently
received significant interest. In this work, an adaptive sampling framework is proposed to reconstruct aquatic environmental fields (e.g., temperature, or biomass of harmful algal blooms) using
schools of robotic sensor platforms. In particular, it is assumed that the field of interest can be
approximated by a low rank matrix, which is exploited for successive expansion of sampling area
and analytical reconstruction of the field. For comparison, an Augmented Lagrange Multiplier
optimization approach is also taken to complete the matrix reconstruction using a limited number
of samples. Simulation results show that the proposed approach is more computationally efficient
and requires shorter travel distances for the robots.

To Mom and Dad
for their love and encouragement, and for building the foundation that supported me throughout
this journey.

iv

ACKNOWLEDGMENTS

I am most grateful to my advisor, Prof. Xiaobo Tan, for his enthusiastic encouragement and
technical guidance during my doctoral study, and for pushing me to excel as an independent researcher. I will be forever grateful for this doctoral research assistant opportunity he has provided
me.
I would like to thank Prof. Hassan Khalil, Prof. Guoming Zhu, and Prof. Guoliang Xing, for
kindly joining my advisory committee and providing many insightful suggestions. In particular,
I thank Prof. Khalil and Prof. Zhu for leading me into the area of advanced controls. I am also
grateful to Prof. Xing for broadening my knowledge, and for aiding me in my research.
I would like to thank Prof. Philips McKinley, for taking the time to edit my work and providing
useful suggestions. I am also grateful for his advice and support on my job applications.
I would like to thank Prof. Hayder Radha for taking the time to review my research and offering
me comments. I also benefited a lot from his stochastic systems, signal compression classes, which
helped me accomplish the work in the last chapter.
I am grateful to my colleagues and friends at Michigan State University who have offered
help in various ways: Prof. Songlin Chen, Yu Wang, Dr. Rui Tan, Ruogu Zhou, Anthony Clark,
Mohammed Al-Qizwini, Dr. Nima Moghadam, Jianguo Zhao, Prof. Bingtuan Gao, Prof. Jing Xu,
Hong Lei, Feitian Zhang, Jun Zhang, Sanaz Behbahani, Prof. Yan Shen, Dr. Alex Esbrook, Dr.
Ahmad T. Abdulsadda, Cody Thon, Freddie Alequin-Ramos, Osama En-Nasr, Felix Adisaputra,
and many others. Special thanks are due to John Thon, for his help building the mechanical parts
of robotic fish platforms, which enables me to focus on other important issues in robotic fish
development.
I am also grateful to the technical and administrative staff of Electrical and Computer Engineerv

ing for their assistance during my study at Michigan State University, in particular, Brian Wright,
Gregg Mulder, Roxanne Peacock, Meagan Kroll, and Laurie Rashid.
I would like to acknowledge the financial support of my research by National Science Foundation (CNS-0751155, CCF-0820220, DBI-0939454, IIS-0916720, CNS-0915855, ECCS 1050236,
ECCS 1029683, CNS-1059373, CCF-1331852, IIS-1319602) and Office of Naval Research (Grant
N-000140810640) for their support of my research.
Last, but not least, I am foremost thankful for my family. I would like to thank my mother
Ruiru Du and my father Yongcheng Wang, for being there for me, for their everlasting support in
my endeavors. I would never have been able to accomplish what I have without them. I am deeply
grateful to my fiancee, Liyang Mao, for her love, support and encouragement.

vi

TABLE OF CONTENTS

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

x

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

xi

Chapter 1
Introduction . . . . . . . . . . . .
1.1 Robotic Fish as Mobile Sensing Platforms
1.2 Organization . . . . . . . . . . . . . . . .
1.3 Overview of Contributions . . . . . . . .
1.4 Publications . . . . . . . . . . . . . . . .
1.4.1 Journal Articles . . . . . . . . . .
1.4.2 Conference Proceedings . . . . .

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

1
1
4
5
6
6
7

Chapter 2
Development of Robotic Fish and Their Applications
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Network-Robotic-Fish . . . . . . . . . . . . . . . . . . . . .
2.2.1 Hardware structure . . . . . . . . . . . . . . . . . . .
2.2.2 Software architecture . . . . . . . . . . . . . . . . . .
2.2.3 System assembly . . . . . . . . . . . . . . . . . . . .
2.3 Gliding Robotic Fish . . . . . . . . . . . . . . . . . . . . . .
2.4 Predator Robotic Fish . . . . . . . . . . . . . . . . . . . . . .
2.5 Eduction Activities and Outreach Using Robotic Fish . . . . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

9
9
9
10
10
12
14
17
18

A Dynamic Model for Tail-actuated Robotic Fish with Drag Coeff cient
Adaptation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The Dynamic Model for Robotic Fish . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Rigid-body dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Evaluation of the hydrodynamic forces . . . . . . . . . . . . . . . . . . .
3.2.2.1 Lighthill’s large-amplitude elongated-body theory . . . . . . .
3.2.2.2 Hydrodynamic force due to an oscillating tail . . . . . . . . . .
3.2.2.3 Drag and lift on the body . . . . . . . . . . . . . . . . . . . . .
Experimental Model Validation . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Parameter identification . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.3 Results on model validation . . . . . . . . . . . . . . . . . . . . . . . .
3.3.4 Hydrodynamic forces with/without incorporating robotic fish body motion
Efficient Adaptation of Drag Coefficients . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.

20
20
22
23
25
25
27
29
30
30
32
34
36
39

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

Chapter 3
3.1
3.2

3.3

3.4

vii

3.5

3.4.1 Bias-based adaptation of drag coefficients . . . . . . . . . . . . . . . . . . 39
3.4.2 Elimination of the angle of attack from the lift force term . . . . . . . . . . 46
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

Chapter 4
Dynamic Modeling of Robotic Fish with a Base-actuated Flexible Tail
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 The Model Framework for Tail-actuated Robotic Fish . . . . . . . . . . . . . .
4.3 Dynamic Modeling of a Tail Using Euler-Bernoulli Beam Theory . . . . . . . .
4.4 Dynamic Modeling of a Tail Using Multi-segment Approximation . . . . . . .
4.5 Experimental Model Validation . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5.1 Parameter identification . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5.2 Model verification . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

52
52
55
56
62
65
66
70
74

Chapter 5
Control-oriented Averaging of Tail-actuated Robotic Fish Dynamics
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Dynamic Model of Robotic Fish . . . . . . . . . . . . . . . . . . . . . . . .
5.3 Averaging with Existing Methods . . . . . . . . . . . . . . . . . . . . . . . .
5.3.1 First-order Averaging for Systems with Oscillatory Inputs . . . . . . .
5.3.2 First-Order Averaging for Robotic Fish Model . . . . . . . . . . . . .
5.3.3 Classical Averaging . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4 Computationally Efficient Averaging with Scaled Forcing . . . . . . . . . .
5.4.1 Proposed Averaging Approach . . . . . . . . . . . . . . . . . . . . .
5.4.2 Numerical Evaluation of the Average Model . . . . . . . . . . . . . .
5.4.3 Experimental Evaluation of the Average Model . . . . . . . . . . . .
5.5 Analysis of the Average Model . . . . . . . . . . . . . . . . . . . . . . . . .
5.6 Application to Semi-analytical Modeling of Steady Turning . . . . . . . . . .
5.6.1 Semi-analytical Modeling Framework . . . . . . . . . . . . . . . . .
5.6.2 Experimental Evaluation . . . . . . . . . . . . . . . . . . . . . . . .
5.7 Nonlinear Controller Design . . . . . . . . . . . . . . . . . . . . . . . . . .
5.7.1 System Description . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.7.2 Hybrid Controller Design . . . . . . . . . . . . . . . . . . . . . . . .
5.7.3 Simulation and Experimental Results . . . . . . . . . . . . . . . . .
5.8 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

80
80
83
84
85
86
91
92
92
100
103
106
110
110
113
114
115
117
119
120

Sampling
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .

.
.
.
.
.
.
.
.
.

123
123
126
128
131
132
133
135
137

Chapter 6
6.1
6.2
6.3
6.4
6.5

Environmental Monitoring with Robotic Sensor Network:
and Reconstruction of Low-rank Matrix Fields . . . . . . . .
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . .
Matrix Completion Using the Augmented Lagrange Multiplier Method
Analytical Reconstruction via Successive Expansion of Sampling Area
Trace-Driven Simulations . . . . . . . . . . . . . . . . . . . . . . . .
6.5.1 Measurement of temperature field in a water tank . . . . . . .
6.5.2 Simulation results . . . . . . . . . . . . . . . . . . . . . . . .
6.5.3 Comparisons of energy cost for data sampling . . . . . . . . .
viii

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

6.6

6.7

Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.6.1 Impact of more accurate energy model . . . . . . . . . . . . . . . . . . .
6.6.2 Limitation of the proposed matrix completion algorithm via successive reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. 139
. 139
. 141
. 144

Chapter 7
Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . . . . 147
7.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
7.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
BIBLIOGRAPHY . . . . . . . . . . . . . . .

ix

. . . . . . . . . . . . . . . . . . . . . 151

LIST OF TABLES

Table 2.1

Selected components used in network robotic fish. . . . . . . . . . . . . . 14

Table 3.1

Parameter values for simulations. . . . . . . . . . . . . . . . . . . . . . . 33

Table 3.2

Identified drag force and moment coefficients for different tail beating
biases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

Table 5.1

Average of moment coefficients for different bias angles. . . . . . . . . . 100

Table 5.2

Surge velocities predicted by the original and average models, and the
relative error by the average model. . . . . . . . . . . . . . . . . . . . . . 101

Table 5.3

Best scaling coefficients for different bias angle. (mb = 0.622 kg in the
simulation) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

x

LIST OF FIGURES

Figure 1.1

Prototypes of robotic fish developed by Smart Microsystems Lab at Michigan State University [1]: (a) testing in a inland lake; (b) prototype for
dynamically monitoring the dissolved oxygen level in aquafarms . . . . .

3

Figure 2.1

Navigation board . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

Figure 2.2

Software architecture of network-robotic-fish. . . . . . . . . . . . . . . . . . 12

Figure 2.3

Three generations of robotic fish: (a) First generation of robotic fish; (b)
DO fish (second generation); (c) network-robotic-fish (third generation) . 13

Figure 2.4

Illustration of gliding robotic fish motion profiles [2]. (a) Rectilinear gliding motion; (b) 3D spiraling. . . . . . . . . . . . . . . . . . . . . . . . . 15

Figure 2.5

Gliding robotic fish prototype developed at Smart Microsystems Lab. [2] . . . . 16

Figure 2.6

Swimming trajectory of GRACE along the Kalamazoo River. [2] . . . . . . . . 16

Figure 2.7

Robot-fish interactions. (a) Predator robotic fish prototype; (b) live stickleback fish approaching the robot. . . . . . . . . . . . . . . . . . . . . . 17

Figure 2.8

Smart Microsystems Lab’s robotic fish exhibit at the 2010 USA Science
and Engineering Festival. . . . . . . . . . . . . . . . . . . . . . . . . . . 19

Figure 3.1

Top view of the tail-actuated robotic fish undergoing planar motion.

Figure 3.2

Illustration of the coordinate system for the spinal column of the elongated body
(top view). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

Figure 3.3

A free-swimming robotic fish prototype used for model validation. . . . . . . . 32

Figure 3.4

Comparison between model prediction and experimental measurement of (a)
turning radius and (b) turning period as a function of the tail beat bias α0 . The
amplitude and frequency of the tail beat are fixed at 15◦ and 0.9 Hz, respectively.

xi

. . . . . . 24

35

Figure 3.5

Comparison between model prediction and experimental measurement of (a)
turning radius and (b) turning period as a function of the tail beat amplitude.
The bias and frequency of the tail beat are fixed at 20◦ and 0.9 Hz, respectively. . 37

Figure 3.6

Comparison between model prediction and experimental measurement of (a)
turning radius and (b) turning period as a function of the tail beat frequency.
The bias and amplitude of the tail beat are fixed at 20◦ and 15◦ , respectively. . . 38

Figure 3.7

Simulated average thrust force with/without incorporating the body motion: (a)
Thrust force versus tail-beat amplitude (frequency fixed at 0.9 Hz); (b) thrust
force versus tail-beat frequency (amplitude fixed at 15◦ . . . . . . . . . . . . . 40

Figure 3.8

Simulated average velocity with/without incorporating the body motion: (a) Velocity versus tail-beat amplitude (frequency fixed at 0.9 Hz); (b) velocity versus
tail-beat frequency (amplitude fixed at 15◦ . . . . . . . . . . . . . . . . . . . 41

Figure 3.9

Simulation results on the dependence of the angle of attack on tail beat parameters, where β0 is the cycle-averaged absolute value of the angle attack: (a) β0
versus the tail beat bias α0 ; (b) β0 versus the tail beat amplitude αA ; (c) β0 versus
the tail beat frequency ωα . In (a), the amplitude and frequency are fixed at 15◦
and 0.9 Hz; in (b), the bias and the frequency are fixed at 20◦ and 0.9 Hz; in (c),
the bias and the amplitude are fixed at 20◦ and 15◦ . . . . . . . . . . . . . . . 43

Figure 3.10

Comparison between experimental measurement of (a) turning radius and (b)
turning period with those predicted by the model, where the drag coefficients are
adapted and the lift force is evaluated using (3.20). The amplitude and frequency
of the tail beat are fixed at 15◦ and 0.9 Hz, respectively. . . . . . . . . . . . . 45

Figure 3.11

Comparison between model predictions and experimental measurement for forward swimming (including transients): (a) time trajectory of X -coordinate of the
robot; (b) time trajectory of the Y -coordinate of the robot; (c) path of the robot
in the XY -plane. For both the experiment and simulations, the amplitude and
frequency of the tail beat are fixed at 15◦ and 0.9 Hz, respectively. . . . . . . . 47

Figure 3.12

Comparison between experimental measurement of (a) turning radius and (b)
turning period with those predicted by the model, where the drag coefficients
are not adapted and the lift force is evaluated using (3.34). The amplitude and
frequency of the tail beat are fixed at 15◦ and 0.9 Hz, respectively. . . . . . . . 49

Figure 3.13

Comparison between experimental measurement of (a) turning radius and (b)
turning period with those predicted by the model, where the drag coefficients are
adapted and the lift force is evaluated using (3.34). The amplitude and frequency
of the tail beat are fixed at 15◦ and 0.9 Hz, respectively. . . . . . . . . . . . . 50

xii

Figure 4.1

Schematic representation of the robotic fish in planar motion.

Figure 4.2

Passive flexible tail modeled by Euler-Bernoulli linear beam approach.

Figure 4.3

Passive flexible tail modeled by multiple rigid segments. . . . . . . . . . 63

Figure 4.4

A free-swimming robotic fish prototype used for model validation. . . . . . . . 66

Figure 4.5

Experimental setup for measuring the Young’s modulus of the flexible tail. . . . 68

Figure 4.6

Experimental results for measuring the Young’s modulus of the flexible tail. . . . 69

Figure 4.7

Displacement of the tail tip generated by the model using different number of
rigid segments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

Figure 4.8

Computation time needed to simulate the model using different numbers of rigid
segments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

Figure 4.9

Comparison between model predictions and experimental measurement of the
speed versus tail-beat frequency. The amplitude is fixed at 13.6◦ . . . . . . . . . 72

Figure 4.10

Comparison between model predictions and experimental measurement of the
speed versus tail-beat frequency, with the new tail. The amplitude is fixed at 13.6◦ . 73

Figure 4.11

Comparison between multi-segment model predictions and experimental measurement for forward swimming (including transients): (a) time trajectory of X coordinate of the robot; (b) time trajectory of the Y -coordinate of the robot; (c)
path of the robot in the XY -plane. For both the experiment and simulations, the
amplitude and frequency of the tail beat are fixed at 13.6◦ and 0.9 Hz, respectively. 76

Figure 4.12

Experimental setup to capture the dynamic shape of a flexible tail actuated at the
base. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

Figure 4.13

Comparison between experimental measurement of the time-dependent tail shape
with model predictions. The tail beats at 0.4 Hz with 0◦ bias and 14◦ amplitude.
The black solid line, blue dashed line with circles and the red dash-dot line imply the experimental measurement, predictions from multi-segment model and
Euler-Bernoulli beam model, respectively. . . . . . . . . . . . . . . . . . . . 78

xiii

. . . . . . 56
. 57

Figure 4.14

Comparison between experimental measurement of the time-dependent tail shape
with model predictions. The tail beats at 0.9 Hz with 0◦ bias and 13.6◦ amplitude.
The black solid line, blue dashed line with circles and the red dash-dot line imply the experimental measurement, predictions from multi-segment model and
Euler-Bernoulli beam model, respectively. . . . . . . . . . . . . . . . . . . . 79

Figure 5.1

The trajectories of robotic fish predicted by the original model and by the classical average model, respectively. . . . . . . . . . . . . . . . . . . . . . . . 92

Figure 5.2

Simulated trajectory of the center of the robotic fish when the tail beats at
1 Hz with bias 20◦ and amplitude 15◦ : (a) The global path; (b) close-up
view of the segment within the red rectangular region in (a). . . . . . . . . 95

Figure 5.3

Simulated heading angle of the robotic fish when the tail beats at 1 Hz
with bias 20◦ and amplitude 15◦ : (a) The heading angle trajectory; (b)
close-up view of the segment within the red rectangular region in (a). . . . 96

Figure 5.4

K f orce and Kmoment versus the tail-beat bias and amplitude. In all cases,
the tail beats at 1 Hz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

Figure 5.5

K f orce and Kmoment versus the tail-beat bias and frequency. In all cases,
the tail beats with amplitude 15◦ . . . . . . . . . . . . . . . . . . . . . . . 98

Figure 5.6

K f orce and Kmoment versus the tail-beat amplitude and frequency. In all
cases, the tail beats with bias 20◦ . . . . . . . . . . . . . . . . . . . . . . . 99

Figure 5.7

Comparison of steady-state turning radii and periods predicted by the original
and average models. In all cases, the tail-beat frequency is 1 Hz and amplitude
is 20◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

Figure 5.8

Simulation results from the original and average models (versus different tailbeat bias), when mb = 0.622 kg. In all cases, the tail-beat frequency is 1Hz and
amplitude is 20◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

Figure 5.9

Simulation results from the original and average models (versus different tailbeat amplitude), when mb = 0.622 kg. In all cases, the tail-beat frequency is
1Hz and bias is 20◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

Figure 5.10

Trajectory of the center of the robot when the tail-beat bias is updated by
a proportional feedback controller law (a) The simulated trajectory using
the original model; (b) The simulated trajectory using the average model;
(c) The control input to both the original and average models. . . . . . . . 105

xiv

Figure 5.11

Comparison between prediction from original dynamic model, average model
and experimental measurement for forward swimming. . . . . . . . . . . . . . 106

Figure 5.12

Comparison between predictions from the average model and experimental measurement for forward swimming (including transients): (a) time trajectory of
X -coordinate of the robot; (b) time trajectory of the Y -coordinate of the robot;
(c) path of the robot in the XY -plane. For both the experiments and the simulation, the amplitude and frequency of the tail beat are fixed at 15◦ and 0.9 Hz,
respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

Figure 5.13

Top view of the tail-actuated robotic fish undergoing “steady” turning motion.

Figure 5.14

Comparison between model predictions and experimental measurement of turning behavior versus tail-beat bias. . . . . . . . . . . . . . . . . . . . . . . . 114

Figure 5.15

Comparison between model predictions and experimental measurement of turning behavior versus tail-beat frequency. . . . . . . . . . . . . . . . . . . . . 114

Figure 5.16

Comparison between model predictions and experimental measurement of turning behavior versus tail-beat amplitude. . . . . . . . . . . . . . . . . . . . . 115

Figure 5.17

Target tracking control schematic diagram. . . . . . . . . . . . . . . . . . . . 116

Figure 5.18

Switching strategy between two controllers. . . . . . . . . . . . . . . . . . . 119

Figure 5.19

Simulation results with the target located at the origin and the robotic fish starting
from different initial conditions (denoted by solid blue squares). The locations
of controller switching are denotes by green diamonds. . . . . . . . . . . . . . 119

Figure 5.20

Experimental setup. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

Figure 5.21

Experimental results on tracking multiple target points continuously. . . . . . . 122

Figure 6.1

Oil spill in the Kalamazoo River. 2010 (Photo Credit: National Wildlife
Federation) (a) Oil map on the surface of Kalamazoo River; (b) discrete
uniform sampling; (c) random samplings . . . . . . . . . . . . . . . . . . 127

Figure 6.2

Sequential reconstruction of field of interest. (a) Successive expansion
of sampling area, where the new samples in each iteration are only taken
from the outer row and column; (b) Illustration of the sampling strategy
with the oil spill example. . . . . . . . . . . . . . . . . . . . . . . . . . . 133

xv

. 112

Figure 6.3

Sketch of the experimental setup for measuring the water-temperature distribution near the surface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134

Figure 6.4

Scalar field of water-temperature. . . . . . . . . . . . . . . . . . . . . . . . 135

Figure 6.5

Reconstruction of field of interest using MC-ALM method. (a) sampled
data; (b) reconstructed field. . . . . . . . . . . . . . . . . . . . . . . . . . 136

Figure 6.6

Reconstruction of field of interest using the proposed method.

Figure 6.7

Robot paths to cover all the sampling points using the multiple traveling salesman method. The red line and the blue denotes the paths of the two robots. . . . 139

Figure 6.8

Robot paths to cover all the sampling points using our proposed adaptive sampling algorithm. The red line and the blue denotes the paths of two robots. . . . 140

Figure 6.9

Scalar field described by a 20 × 20 matrix. . . . . . . . . . . . . . . . . . . . 143

Figure 6.10

Reconstruction of field of interest using Algorithm 2. . . . . . . . . . . . . . . 144

Figure 6.11

Reconstruction of field of interest using Algorithm 3. . . . . . . . . . . . . . . 145

Figure 6.12

Locations of samplings using Algorithm 3. . . . . . . . . . . . . . . . . . . . 146

xvi

. . . . . . . . . 137

Chapter 1
Introduction
Biomimetic systems have been receiving increasing attention from the robotics community, since
natural organisms can provide important insights into the theory and design of engineer systems.
For example, in the area of aquatic robots, the maneuverability and efficiency of live fish [3, 4]
has motivated significant scientific interest over the past two decades in developing, modeling and
controlling robotic fish [5–19]. In addition to providing platforms for underwater applications such
as environmental monitoring [1, 20–23], these robots offer a means to study the behavior of live
fish [24].
In this chapter, a review of robotic fish related state of arts is given in Section 1.1. The general
organization of this report and major contributions of this work are clarified in Section 1.2 and 1.3,
respectively. The author’s publications during his Ph.D research are listed at Section 1.4

1.1 Robotic Fish as Mobile Sensing Platforms
With five hundred million years of evolution, fish and other aquatic animals are endowed with
a variety of morphological and structural features that enable them to move through water with
speed, efficiency, and agility [3, 4]. The remarkable feats in biological swimming have stimulated extensive theoretical [18, 25, 26], experimental [27–31], and computational [32, 33] research
by biologists, mathematicians, and engineers, in an effort to understand and mimic locomotion,
maneuvering, and sensing mechanisms adopted by aquatic animals.
1

Over the past two decades, there has also been significant interest in developing underwater
robots that propel and maneuver themselves like real fish do [5, 7, 11, 12, 19, 24, 34–49]. Often termed robotic fish, these robots provide an experimental platform for studying fish swimming, and hold strong promise for a number of underwater applications such as environmental
monitoring [1]. Instead of using propellers, robotic fish accomplish swimming by deforming the
body and/or fin-like appendages, mostly functioning as caudal fins and sometimes as pectoral fins.
Body deformation and fin movements are typically achieved with motors. On the other hand, advances in smart materials have been explored to realize noiseless and compact actuation of robotic
fish [6, 9, 12, 45, 50–55]. Fig. 1.1a shows a robotic fish swimming in an inland lake. Fig. 1.1b
shows the close-up of another robot prototype, equipped with a dissolved oxygen (DO) sensor,
global position system (GPS), internal measurement unit (IMU), and other electronic components,
which is designed to monitor the DO level in aquafarms.
Aquatic environments have been facing various threats, and monitoring them and reconstructing aquatic process are of great importance to public health and biological balance. Manual sampling is still a common practice in environment monitoring, which is labor-intensive with difficulty
in capturing dynamic phenomena of interest. Deployment of in-situ sensing with fixed or buoyed
sensors or vertical profilers has been proposed [56–61]. However, since buoyed sensors have little freedom to move around, it cannot provide high spatial resolution to effectively adapt to the
destruction of parts of the network, and will require a large number of units for capturing the distributed information of dynamic phenomena. The past decade has also seen great progress in the
use of robotic technology, like autonomous underwater vehicles (AUVs) [62–68], and underwater
gliders [69–71] for aquatic environmental sensing [63, 66–68, 72–84]. However, it is difficult to
deploy many of them due to the high manufacturing and operational costs.
On the other hand, schools of autonomous robotic fish sensing platforms provide a competitive
2

(a)

(b)
Figure 1.1: Prototypes of robotic fish developed by Smart Microsystems Lab at Michigan State
University [1]: (a) testing in a inland lake; (b) prototype for dynamically monitoring the dissolved
oxygen level in aquafarms
alternative for the task of environment monitoring, with their low-cost, high-volume traits. However, the minimal human intervention and limited onboard resources prohibit the robots’ ability to
operate robustly in the unfriendly and unpredictable environment. This poses challenges across a
wide spectrum, ranging from locomotion and maneuvering mechanisms, energy-efficient deigns,
3

communication schemes, control and coordination, to sampling and reconstruction strategies.
This dissertation spans robotic fish development, dynamic modeling, averaging, nonlinear controller design, and an adaptive sampling framework to guide a group of robotic fish.

1.2 Organization
This dissertation is organized as follows. In Chapter 2, guidelines for designing robotic fish are
provided and a detailed description of one robotic fish system is explained as a representative
example. Then some other prototypes that the author made significant contributions to, e.g., a
gliding robotic fish and a robotic predator are briefly discussed. The applications of these robots
are also covered in this chapter. In Chapter 3, a complete dynamic model is presented for robotic
fish actuated by a rigid caudal fin, where Lighthill’s elongated-body theory is used to evaluate
the tail-generated hydrodynamic forces. It shows that, the drag force and moment coefficients are
strongly dependent on the tail-beat bias, and subsequently a novel scheme is proposed to adapt
these coefficients and lift force based on the tail movement pattern. In Chapter 4, a model for
robotic fish with a base-actuated flexible tail is developed. The tail is modeled as multiple rigid
segments connected by springs and dampers, which is compared with the model based on the
linear beam theory. In Chapter 5, a control-oriented averaging approach for tail-actuated robotic
fish dynamics is developed, which is applied to a hybrid controller for target control. In Chapter 6,
an adaptive sampling framework is developed for a group of robotic fish, where it is assumed that
the field can be represented by a low-rank matrix. Motivated by motion constraints of the robots,
the presented approach features successive expansion of the sampling area, where increasingly
larger sub-matrices of the original field are reconstructed analytically.

4

1.3 Overview of Contributions
The contributions of this dissertation are primarily on the complex dynamic modeling of tail-driven
robotic fish, its related averaging and control theory, and an adaptive sampling framework using a
group of robots. These contributions are further elaborated below.
First, the guidelines of robotic fish system design are provided and several types of robotic fish
prototypes are developed for aquatic environmental monitoring, robot-animal interaction study,
and education and outreach.
Second, a dynamic model for a tail-actuated robotic fish has been developed. In addition, a
computationally efficient approach is proposed and validated for adapting the drag coefficient of
the robotic fish based only on the tail-beat bias.
Third, to design and control robotic fish actuated with a flexible tail, it is essential to have
a faithful and efficient dynamic model. Existing studies typically adopt a linear Euler-Bernoulli
beam model is adopted to describe the beam dynamics; however, these models cannot accurately
capture the beam dynamics when the beam undergoes large deformations. A model consisting of
multiple, serially connected rigid segments is developed to describe the flexible beam, the effectiveness of which is validated with experimental results.
Fourth, a novel control-oriented data-driven averaging approach is proposed to approximate
the original dynamics of robotic fish. The merit of the proposed approach, is demonstrated in the
design of a hybrid controller for targeting tracking and in the application to study of steady turning.
Last, a school of robotic fish is proposed for adaptive sampling and reconstruction of aquatic
fields. Distinguished from most existing work, we propose to expand the sampling area successively and reconstruct the field iteratively via analytical solutions. The effectiveness of the approach is validated via extensive simulations based on real data traces of the temperature field in

5

a large water tank. It shows that the proposed adaptive sampling solution is more computationally
efficient and requires shorter travel distances for the robots than a competing matrix completionbased method.

1.4 Publications
1.4.1 Journal Articles
1. J. Wang, X. Tan, “Averaging of Tail-Actuated Robotic Fish Dynamics through Force and
Moment Scaling”, IEEE Transactions on Robotics, under review, 2014
2. J. Wang, P. K. McKinley, X. Tan, “Dynamic modeling of robotic fish with a base-actuated
flexible tail”, Journal of Dynamic Systems, Measurement, and Control, to appear, 2013
3. Y. Wang, R. Tan, G. Xing, X. Tan, J. Wang, R. Zhou, “Spatiotemporal aquatic field reconstruction using cyber-physical robotic sensor systems”, ACM Transaction on Sensor Networks, vol. 10, no. 57, 2014
4. Y. Wang, R. Tan, G. Xing, J. Wang, X. Tan, “Profiling aquatic diffusion process profiling
using robotic sensor networks”, IEEE Transactions on Mobile Computing, vol. 13, no. 4,
pp. 880-893, 2014
5. J. Wang, X. Tan, “A dynamic model for tail-actuated robotic fish with drag coefficient adaptation”, Mechatronics, vol. 23, no. 6, pp. 659-668, 2013

6

1.4.2 Conference Proceedings
1. C. Liu, S. Chen, J. Wang, X. Tan, “Study on shortest path planning for robotic fish actuated
by a caudal fin” [in Chinese], Proceedings of the 33rd Chinese Control Conference, Nanjing,
China, to appear, 2014
2. Y. Wang, R. Tan, G. Xing, J. Wang, X. Tan, X. Liu, “Aquatic debris monitoring using
smartphone-based robotic sensors”, Proceedings of the 13th ACM/IEEE International Conference on Information Processing in Sensor Networks, Berlin, Germany, pp. 13-24, 2014
[Acceptance rate: 21%; Runner-up for Best Paper Award]
3. F. Zhang, J. Wang, J. Thon, C. Thon, E. Litchman, X. Tan, “Gliding robotic fish for mobile sampling of aquatic environments” (Invited), Proceedings of 11th IEEE International
Conference on Networking, Sensing and Control, Miami, FL, pp. 167-172, 2014
4. S. Chen, J. Wang, X. Tan, “Target-tracking control design for a robotic fish with caudal fin,”
Proceedings of 32nd Chinese Control Conference, Xi’an, China, pp. 844-849, 2013
5. S. Chen, J. Wang, X. Tan, “Backstepping-based hybrid target tracking control for a carangiform robotic fish,” Proceedings of 2013 ASME Dynamic Systems and Control Conference,
Palo Alto, CA, Paper DSCC2013-3963, 2013
6. S. B. Behbahani, J. Wang, X. Tan, “A dynamic model for robotic fish with flexible pectoral
fins,” Proceedings of 2013 IEEE/ASME International Conference on Advanced Intelligent
Mechatronics, Wollongong, Australia, pp. 1552-1557, 2013
7. J. Wang, S. Chen, X. Tan, “Control-oriented averaging of tail-actuated robotic fish dynamics” (Invited), Proceedings of the 2013 American Control Conference, Washington, DC, pp.
591-596, 2013
7

8. Y. Wang, R. Tan, G. Xing, X. Tan, J. Wang, R. Zhou,“Spatiotemporal aquatic field reconstruction using robotic sensor swarm,” Proceedings of the 33rd IEEE Real-Time Systems
Symposium, San Juan, Puerto Rico, pp. 205-214, 2012 [Acceptance rate: 22%]
9. J. Wang, P. K. McKinley, X. Tan, “Dynamic modeling of robotic fish with a flexible caudal
fin” (Invited), Proceedings of the 2012 ASME Dynamic Systems and Control Conference,
Fort Lauderdale, FL, Paper DSCC2012-MOVIC2012-8695, 2012
10. A. J. Clark, J. M. Moore, J. Wang, X. Tan, P. K. McKinley, “Evolutionary design and experimental validation of a flexible caudal fin for robotic fish,” Proceedings of the 13th International Conference on the Simulation and Synthesis of Living Systems (Artificial Life 13),
East Lansing, MI, pp. 325-332, 2012 [Best Paper Award in Behavior and Intelligence]
11. Y. Wang, R. Tan, G. Xing, J. Wang, X. Tan, “Accuracy-aware aquatic diffusion process
profiling using robotic sensor networks”, Proceedings of the 11th ACM/IEEE Conference on
Information Processing in Sensor Networks, Beijing, China, pp. 281-292, 2012 [Acceptance
rate: 15%]
12. J. Wang, F. Alequin-Ramos, X. Tan, “Dynamic modeling of robotic fish and its experimental Validation” (Invited), Proceedings of the 2011 IEEE/RSJ International Conference on
Intelligent Robots and Systems, San Francisco, CA, pp. 588-594, 2011 [Acceptance rate:
32%]

8

Chapter 2
Development of Robotic Fish and Their
Applications
2.1 Introduction
Over the past two decades, there has been significant interest in the development of robotic fish.
These low-cost, small form-factor platforms present a variety of constraints for design. The main
design criteria in this work are to develop a robust system that is low-cost and easy to be manufactured, extended and serviced. The chapter is organized as follows. In Section 2.2, a detailed
description of the network-robotic-fish is presented. In Section 2.3, a gliding robotic fish prototype
and its application to environmental monitoring are described. In Section 2.4, a predator robotic
fish and its application to robot-fish interaction are discussed. The use of robotic fish in education
and outreach activities are briefly summarized in Section 2.5. These prototypes were developed
jointly with several other members of the Smart Microsystems Lab.

2.2 Network-Robotic-Fish
A robotic fish prototype, called network-robotic-fish in this dissertation, is first described. The
name comes from the intended purpose of this class of robotic fish for cooperative sensing as a
mobile network. The design of such robots is faced with several challenges. Simplicity and ease of
9

fabrication are critical to keeping the cost of the robot node low, however, system robustness and
modularity are also important considerations.

2.2.1 Hardware structure
To address the challenges discussed above, we have developed a two-layer hardware structure
for the robotic node with two micro-controllers. The navigation board is in charge of all the
navigation-related functions, including drivers of sensors and actuator, and their coordinations. On
the other hand, the processing board to the board takes care of higher-level information processing
and decision-making tasks. The processing board also deals with communication and coordination
among individual nodes. Such a design can significantly reduce the chances of conflict between
time-critical components, such as the wireless radio and the navigation controller. More importantly, the hardware results in good extensibility and robustness properties for the design. The
navigation board is relatively fixed in terms of low-level control of the robot and plenty of hardware interfaces are left for additional sensors. On the other hand, the processing board is designed
to accommodate different applications.
The processing board is a duplicate of a standard Telosb mote, which features a MSP430 microcontroller and CC2420 radio stack. The layout, as shown in Fig. 2.1, of the the navigation board
for network-robotic-fish accounts for interference isolation and minimization of total space.

2.2.2 Software architecture
The robotic fish software architecture, shown in Fig. 2.2, consists of the firmware layer, the nodelevel module controllers layer, and customizable head-level computation and user-level application
layer.

10

Gyro
y

Power

Compass
p

GPS

Microcontroller
SPI

Current

Temperature

HABs/Crude Oil

(a)
Switching Regulator

Linear Regulators
Servo connector
Inductor

(b)
Figure 2.1: Navigation board
The novelty of this design is that we implement the software using two self-governed systems
with an SPI interface between them. On the navigation board, a motion controller and sensor
controller are developed upon the drivers for the hardware. On the processing board, a storage
controller and network controller are implemented to control the data flow in the node and between
the nodes, respectively. Furthermore, a user-level applications layer is developed to adapt different
objectives.

11

Processing Board

Sensor Controller
(
(Coordination
d
&
Filtering)

Sto
orage
Con
ntroller

Motion
Controller

SPI Co
ontroller

Applications

SPI Co
ontroller

Contro
oller

App
plication

Navigation Board

802.15.4
Radio
Stack
Network
Layer

Hardw
ware

MAC
Layer

Servo
motor

Sensors

SPI

SPI

Flash
Storage

Radio

Figure 2.2: Software architecture of network-robotic-fish.

2.2.3 System assembly
One of the most time-consuming tasks in robotic fish development is the fabrication of shell with
a streamlined shape with good hydrodynamic performance. In the first two generations of robotic
fish, shown in Figs. 2.3(a) and 2.3(b) we have manually manufactured the shells. The mold has
to be redesigned every time when evolving the robot’s body. Another disadvantage of such an
approach is that the body is not accessible1 , which creates significant problems when it comes to
diagnosis and service of the robot.
To accelerate the assembling process, an Objet Connex350 Multi-Material 3D printer is utilized
to design and develop the robot shell. An important feature of this 3D printed platform is the selfcontainment of the system structure: the waterproof servo can be directly mounted to the bracket
at the end portion; the navigation board is placed on the stand in the shell; extra capsules are left
1 The only way to access the robot’s body is to cut off the shell and repair it afterwards, which

is time-consuming and reduces the reliability.

12

(a)

Zigbee antenna
GPS
Charger &
programming port
Servomotor
(in silicone)

DO sensor

Heat sink

(b)
GPS holder

Zigbee antenna

Waterproof
servomotor

Pectoral fin
(for balance)

Temperature
sensor

(c)
Figure 2.3: Three generations of robotic fish: (a) First generation of robotic fish; (b) DO fish
(second generation); (c) network-robotic-fish (third generation)
13

Table 2.1: Selected components used in network robotic fish.

1
2
3
4
5
6
7
8
9

Component name
Micro-controller ()
Micro-controller (on TelosB)
Wireless module (on TelosB)
Battery
Servo motor
GPS
Gyro
Accelerometer+Compass
Temperature sensor

Component model
Microchip dsPIC6014A
TI MSP430
TI CC2420
Batteryspace 7.4V (13.32Wh) Li-ion Polymer Battery Pack
Hitec Servo HS-5646WP (waterproof)
Garmin GPS 18x LVC
ST LPY503AL
ST LPY503AL
TMP36

for the balance purposes. Any of these components can be easily detached and replaced within a
hour or so, which highly improves the flexibility and extensibility of the system design, and also
cuts the cost.
Table 2.1 lists the details of components on network robotic fish mentioned before in this
section.

2.3 Gliding Robotic Fish
Gliding robotic fish is a new type of underwater vehicles, which is a combination of underwater
glider and robotic fish. Other than the swimming motion by flapping the tail as robotic fish do, the
robot has a mechanism of gliding by utilizing the difference between its gravity and buoyancy.
Fig. 2.4 illustrates the two gliding-based locomotion modes for the gliding robotic fish. The
robot achieves a zig-zag motion by periodically varying its weight and the center of mass, as shown
in Fig. 2.4(a). It can also realize a novel 3D spiral motion, shown in Fig. 2.4(b) by deflecting the
tail during gliding.
We have successfully built a gliding robotic fish prototype, called GRACE (Gliding Robot

14

(a)

(b)
Figure 2.4: Illustration of gliding robotic fish motion profiles [2]. (a) Rectilinear gliding motion;
(b) 3D spiraling.
ACE) [2, 85], as shown in Fig. 2.5.
In November, 2012, we took GRACE, which was equipped with a crude oil sensor , to Kalamazoo River2 , MI, to detect the crude oil content. In this test, GRACE swam along the river by
flapping its tail as robotic fish and collected the data. The data obtained from the sensor were
2 Kalamazoo River was polluted by an oil spill in July 2010, when a pipeline burst and oil

flowed into Talmadge Creek, a tributary of the Kalamazoo River.
15

Figure 2.5: Gliding robotic fish prototype developed at Smart Microsystems Lab. [2]
transmitted wirelessly to the base station with the location information from the GPS unit. Fig. 2.6
shows one sampling trajectory of GRACE.

Figure 2.6: Swimming trajectory of GRACE along the Kalamazoo River. [2]

Results from this section have appeared in the publication [2]. I was mainly responsible for
system integration of GRACE and troubleshooting during the tests.
16

2.4 Predator Robotic Fish
We are also interested in using our robotic fish to study behaviors of live fish. Collaborating with
Dr. Boughman’s group at Michigan State University, we have developed a predator robotic fish to
study the response behavior of stickleback fish as shown in Fig. 2.7.

(a)

(b)
Figure 2.7: Robot-fish interactions. (a) Predator robotic fish prototype; (b) live stickleback fish
approaching the robot.

A C# code is implemented to capture the signal from a joystick and transfer it to the Xbee

17

module3 , as shown in Fig.2.7(a). Users can control the robot’s tail motion (i.e., beat-frequency
and beat-bias) by operating the joystick. The robot is shaped and painted to mimic a trout, which
serves as a predator to stickleback fish. Fig. 2.7(b) shows one frame of the video taken when a
stickleback swam up tp the robot.

2.5 Eduction Activities and Outreach Using Robotic Fish
Aside from their vast potentials for research, robotic fish offer rich opportunities for reaching out to
students and the general public. Bio-inspired robots are valuable and effective means for attracting
students to the areas of science, technology, engineering and math (STEM) [86–88].
Robotic fish have been used as educational tools all over the world. For instance, the carp-like
robots in [8, 89] have attracted thousands of youth during their exhibits at London Aquarium [89].
We have presented our robotic fish at various outreach programs, including the first and second
USA Science and Engineering Festivals, in 2010 and 2012, respectively. At these events, the
visitors got a chance to not only see the robotic fish in action, but also interacted with our lab
members to understand the technology and its applications, as shown in Fig. 2.8.

3 http://www.digi.com/xbee/.

18

(a)

(b)
Figure 2.8: Smart Microsystems Lab’s robotic fish exhibit at the 2010 USA Science and Engineering Festival.

19

Chapter 3
A Dynamic Model for Tail-actuated Robotic
Fish with Drag Coeff cient Adaptation
3.1 Introduction
Dynamic modeling [11, 12, 18, 19, 45, 47, 90–97], trajectory planning, and control [11, 38, 98–103]
of robotic fish have received extensive attention. A major challenge in the modeling of robotic fish
lies in properly capturing the fluid-structure interactions and the resulting force and moment on the
robot. While computational fluid dynamics (CFD) modeling [32, 33, 104, 105] can provide highfidelity representation of such interactions, it is not amenable to control design. A few alternatives
are available. The first is to apply quasi-steady lift and drag models from airfoil theory to the
body and fin surfaces of underwater robots [11, 90, 93, 94]. Another approach is to assume perfect
fluids (irrotational potential flow) and exploit the symmetry to obtain a finite-dimensional model
[95, 102], but some coefficients of such models may have to be obtained numerically [95]. Effects
of vorticity can be accommodated by assuming, for example, vortices periodically shed from the
tail fin [18, 92].
The third approach to the evaluation of hydrodynamic forces is to use Lighthill’s elongatedbody theory [25]. This theory captures the reactive force between body and fluids with an addedmass effect, and it approximates the effect of wake dynamics by considering the momentum balance in a half-space control volume, which contains the deforming body and is bounded by a
20

vertical plane attached to the trailing edge. Consequently, the resulting force and moment only
depend on the configuration of the (robotic) fish. Lighthill’s model strikes a sound balance between fidelity and simplicity, and its effectiveness in robotic fish modeling has been demonstrated
in [45, 106].
The contributions of the work reported in this chapter are two-fold. First, we present a dynamic
model for a tail-actuated robotic fish by combining rigid-body dynamics with Lighthill’s largeamplitude elongated-body theory, and validate the model with extensive experiments conducted
on a robotic fish prototype. Specifically, we consider rigid-body dynamics for the robot body
and adopt Lighthill’s theory to capture the tail-actuation effect. To our best knowledge, such an
approach has not been reported in the literature. Existing work [45, 106] involving Lighthill’s
theory was focused on steady-state motions, where the dynamics of the robot was not considered.
Chen et al. applied a similar approach to evaluate the thrust produced by undulatory body motion
of a robot tuna [107], but the analysis therein dealt with only point-mass dynamics for the robot,
where forward-swimming was studied via simulation. Note that rigid body dynamics allows one
to capture both body translation and rotation, and is significantly more complex and relevant than
point-mass dynamics in terms of capturing the true behavior of the robot. Using the proposed
model, we investigate the role of incorporating the body motion in evaluating the tail-produced
hydrodynamic force, and show that ignoring the body motion (as often done in the literature [12,
14, 45]) results in significant overestimate of the thrust force and robot speed.
The second contribution is a computationally efficient approach to the adaptation of drag coefficients for robotic fish. These coefficients play a significant role in the modeling of underwater
vehicles and are often functions of the angle of attack [108, 109]. For a carangiform robotic fish,
however, its body and/or fins constantly undergo large deformation or rigid-body motion, and
consequently, the angle of attack oscillates due to the body recoil. For these reasons, modeling
21

variable drag coefficients is challenging, and existing work on robotic fish modeling has predominantly adopted constant drag coefficients [12, 93, 110]. In this work, we propose a novel approach
to adapting the drag coefficients. For ease of presentation, we focus on a robotic fish actuated
by an oscillating tail, where the tail beat bias, frequency, and amplitude are the three kinematic
parameters one can change. We first study the effect of each parameter on the average of the angle
of attack, and find that the bias plays a much more significant role than the other two. We then
propose determining the drag coefficients as simple polynomial functions of the tail beat bias. The
effectiveness of the adaptation scheme is demonstrated through the comparison of simulation and
experimental results, for both steady turning and transient dynamics in forward swimming.
The remainder of the chapter is organized as follows. First, the dynamic model is described
in Section 3.2. Experimental validation of the model is presented in Section 3.3. The drag coefficient adaptation scheme is discussed in Section 3.4. Finally, concluding remarks are provided in
Section 3.5.

3.2 The Dynamic Model for Robotic Fish
The robot is considered to consist of two parts, a body that does not deform and a moving or
deforming tail. The motion of the robot body is governed by rigid-body dynamics with the addedmass effect incorporated. Lighthill’s large-amplitude elongated-body theory is utilized to evaluate
the hydrodynamic forces generated by the actuating tail (as opposed to the whole robot). We
assume that the depth of the tail does not change abruptly along the length direction, thus meeting
the “elongated body” requirement [25]. Although the proposed modeling approach is applicable
for a general actuating tail, for ease of presentation, we focus on the case where the tail is rigid
and actuated at its base. As mentioned in Section 3.1, one contribution of this work lies in the

22

combination of the two elements, rigid-body dynamics and Lighthill’s theory, for the dynamic
modeling of robotic fish. While both elements are well known, for completeness of the chapter,
brief review of each is provided in Sections 3.2.1 and 3.2.2.1, respectively.

3.2.1 Rigid-body dynamics
Fig. 3.1 shows the top view of the robotic fish, consisting of a rigid body and a tail. We follow [12]
closely in describing the dynamics of the rigid body. [X ,Y, Z] denote the inertial coordinates while
[x, y, z] denote the body-fixed coordinates with unit vectors [x,
ˆ y,
ˆ zˆ]. We denote by mˆ and nˆ the
unit vectors parallel and perpendicular to the tail, respectively. We assume that both the body
and the tail are neutrally buoyant, and that, for the body, the center of gravity coincides with
the center of geometry at point C. The velocity at C expressed in the body-fixed coordinates
VC = [VCx ,VCy ,VCz ]T comprises surge (VCx ), sway (VCy ) and heave (VCz ) components. In addition,
the angular velocity ω = [ωx , ωy , ωz ]T comprises roll (ωx ), pitch (ωy ) and yaw (ωz ) expressed
in the body-fixed coordinates. We use α to denote the tail deflection angle with respect to the
negative x−axis, and β to denote the angle of attack, formed by the direction of VC with respect to
the x−axis. Finally, ψ denotes the heading angle, formed by the x−axis relative to the X −axis.
The linear momentum P and angular momentum H of the body in the body-fixed coordinates
are expressed as

P =M · VC + DT · ω ,

(3.1)

H =D · VC + J · ω ,

(3.2)

where M and J are the mass and inertia matrices, respectively, and D is the Coriolis and centripetal
matrix. For a rigid body in an inviscid fluid, Kirchhoff’s equations of motion in the body-fixed
23

Y
yˆ

MD

y

Zz
FL

FD

&
VC

x
B

E
\

Fish body
z

xˆ

C

A

D

Z

D

mˆ

Tail
nˆ

X

x

Figure 3.1: Top view of the tail-actuated robotic fish undergoing planar motion.
frame are [111, 112]

P˙ =P × ω + F,

(3.3)

H˙ =H × ω + P × VC + M,

(3.4)

where F = [Fx , Fy , Fz ]T , M = [Mx , My , Mz ]T denote the external forces and moments about the
body center C, respectively, and “×” denotes the vector product. The assumption of neural buoyancy implies that these forces and moments will only come from the hydrodynamic interactions
between the robot and the fluid. We further assume that the body is symmetric with respect to
the xz−plane and the tail moves in the xy−plane. Consequently, the heave velocity VCz , the roll
rate ωz , and the pitch rate ωy are all equal to zero, in which case the system has three degrees of
freedom, namely, surge (VCx ), sway (VCy ) and yaw (ωz ). We further assume that the inertial cou-

24

pling between the surge, sway and yaw motions is negligible [12], implying D = 0. Under these
assumptions, following [111], (3.3) and (3.4) can be simplified as

(mb − XV˙ )V˙Cx =(mb −YV˙ )VCy ωz + Fx ,
Cx
Cy

(3.5)

(mb −YV˙ )V˙Cy = − (mb − XV˙ )VCx ωz + Fy ,
Cy
Cx

(3.6)

(Jbz − Nω˙ z )ω˙ z =(YV˙

Cy

− XV˙ )VCxVCy + Mz ,
Cx

(3.7)

where mb is the mass of the body, and Jbz is the inertia of the body about the z−axis. XV˙ , YV˙
Cx
Cy
and Nω˙ z are the hydrodynamic derivatives that represent the effect of added mass/inertia on the
body.

3.2.2 Evaluation of the hydrodynamic forces
To complete the model (3.5) – (3.7), we need to evaluate the hydrodynamic force (Fx , Fy ) and
moment Mz , which are due to both the tail actuation and the interaction of the body itself with the
fluid. We first review Lighthill’s large-amplitude elongated-body theory, and use it to evaluate the
forces exerted on the tail (and transmitted to the rigid body) due to tail-fluid interactions. The drag
and lift forces on the robotic fish body will then be described.

3.2.2.1 Lighthill’s large-amplitude elongated-body theory
Consider an elongated “body” moving in the XY -plane in water; here (and throughout Section 3.2.2.1)
the term “body” could mean a live fish, a robotic fish, or in the context of our work, a flapping tail.
As illustrated in Fig. 3.2, a frame of reference is chosen such that the water far from the body
is at rest. The spinal column (analogy from a live fish) of the elongated body is assumed to be
inextensible and is parameterized by τ , with τ = 0 denoting the anterior end of the body and τ = L
25

denoting the posterior end. The coordinates (X (τ ,t),Y(τ ,t)), 0 ≤ τ ≤ L, denote the trajectory of
each point τ on the spinal column at time t, which could be due to the body undulation/oscillation
or the resulting translational/rotational motion of the body.

nˆ

Y

mˆ
vA
W

Z

v||

0
Spinal
column

x

W

L

X

Figure 3.2: Illustration of the coordinate system for the spinal column of the elongated body (top view).
Following [25], given (X (τ ,t),Y(τ ,t)), the hydrodynamic reactive force density at each point

τ < L, due to the added-mass effect, is




d
 f X (τ ) 
ˆ
f (τ ) = 
 = −m (v⊥ n),
dt
fY (τ )

(3.8)

and at τ = L, there is a concentrated force




1 2
 FLX 
FL = 
.
 = − mv⊥ mˆ + mv⊥ v nˆ
2
τ
=L
FLY

(3.9)

The two terms in (3.9) account for the pressure force acting on Π and the force due to convection of
momentum out of V across Π, respectively, where Π is the plane at τ = L that is perpendicular to
m,
ˆ and V is the half-space bounded by Π that includes the fish but excludes the wake [25]. In (3.8)
26

and (3.9), m denotes the virtual mass per unit length and can be approximated by 14 πρ d 2 , where

ρ is the density of water and d is the depth of fish cross-section (in −Z direction) at τ (and thus
can vary with τ ). As illustrated in Fig. 3.2, mˆ = (∂ X /∂ τ , ∂ Y /∂ τ )T and nˆ = (−∂ Y /∂ τ , ∂ X /∂ τ )T
represent the unit vectors tangential and perpendicular to the spinal column, respectively, and v
and v⊥ represent the components of the velocity v = (∂ X /∂ t, ∂ Y /∂ t)T at τ in mˆ and nˆ directions,
respectively:

∂ X ∂ X ∂Y ∂Y
+
,
∂t ∂τ
∂t ∂τ
∂ X ∂Y ∂Y ∂ X
v⊥ = < v, nˆ >= −
+
,
∂t ∂τ
∂t ∂τ
v = < v, mˆ >=

(3.10)
(3.11)

where < ·, · > denotes the inner product of vectors. Note that the actual mass of the considered
elongated-body is negligible (relative to the virtual mass) when the width of the body cross section
is much thinner than the depth [25]. The latter typically holds true for the tail of a robotic fish.

3.2.2.2 Hydrodynamic force due to an oscillating tail
Following Lighthill’s large-amplitude elongated-body theory, we need to know the motion of every
point along the tail over time, to evaluated the tail actuation-induced hydrodynamic force.
Refer to Fig. 3.1. The velocity VA at the body-tail joint A can be expressed as

VA = VC − ωz cy,
ˆ

(3.12)

where c is the distance from the body center to the joint. The velocity at each point τ along the tail
is
ˆ
Vτ = VC − ωz cyˆ + (α˙ + ωz )τ n,

27

(3.13)

where α˙ denotes the angular velocity of the tail with respect to the body. Note here that τ = 0
would represent the point A. In (3.13), Vτ is calculated incorporating the robot body motion (VC
and ωz ), which is often ignored in the literature [12, 14, 45]. In Section 3.3.4, we will discuss the
role of incorporating the body motion in the appropriate evaluation of hydrodynamic forces.
ˆ Using (3.8), we
Let Vτ⊥ be the projection of Vτ along the nˆ direction, i.e., Vτ⊥ =< Vτ , nˆ > n.
evaluate the hydrodynamic force density at τ as

f (τ ) = −m

dVτ⊥
dt

.

(3.14)

Then the hydrodynamic force on the tail due to the added-mass effect is

Fa =

L
0

ˆ
f (τ )dτ = Fam mˆ + Fan n,

(3.15)

where L is the length of the tail, and Fam and Fan are the components of Fa along mˆ and n,
ˆ respectively.
From (3.9), the concentrated force acting at the tail tip is evaluated as Fc = Fc1 mˆ + Fc2 n,
ˆ where
1
Fc1 = − mVL2 ,
⊥
2
Fc2 =mVL⊥ VL .

The total hydrodynamic force acting on the tail is then

Fh = Fhx · xˆ + Fhy · y,
ˆ

28

where

Fhx = − (Fam + Fc1 ) cos α + (Fan + Fc2 ) sin α ,

(3.16)

Fhy = − (Fam + Fc1 ) sin α − (Fan + Fc2 ) cos α .

(3.17)

The hydrodynamic force-induced moment relative to the center C of the robot body is

Mh =

L
0

rCτ × f (τ ) d τ + rCL × Fc ,

(3.18)

where rCτ denotes the vector from the body center C to the point τ on the tail, i.e., rCτ = −(c +

τ cos α ) xˆ − (τ sin α ) y.
ˆ It is clear that all vectors in (3.18) lie in the xy-plane, and thus the only
non-zero component of Mh is along the z−direction, which we denote as Mhz.
3.2.2.3 Drag and lift on the body
As shown in Fig. 3.1, besides the hydrodynamic force and moment transmitted from the tail, the
robot body also experiences drag force FD , lift force FL , and drag moment MD [11, 12, 93]:
1
FD = ρ |VC |2 SCD ,
2

(3.19)

1
FL = ρ |VC |2 SCLβ β ,
2

(3.20)

MD = −KD ωz2 sgn(ωz ).

(3.21)

where S is a suitably defined reference surface area for the robot body. CD and KD are the drag
force coefficient and drag moment coefficient, respectively. CLβ is the lift coefficient per angle of
attack β . The drag coefficients, as we will discuss later, are generally dependent on the angle of
29

attack β of the robot body.
Finally, by adding the hydrodynamic forces and moments from the tail and directly on the body,
we obtain Fx , Fy , and Mz in (3.5) – (3.7) as

Fx =Fhx − FD cos β + FL sin β ,

(3.22)

Fy =Fhy − FD sin β − FL cos β ,

(3.23)

Mz =Mhz + MD .

(3.24)

Note that, Fhx , Fhy and Mhz contain the derivative of VCx , VCy , and ωz since the robotic fish body’s
motion is incorporated in the evaluation of hydrodynamic forces at the tail section. With fundamental transformations, Eq. (5.1) could be re-written with VCx , VCy and ωz on the right side of the
equations and their time derivative on the left side to avoid the numerical loop in the simulation.

3.3 Experimental Model Validation
3.3.1 Experimental setup
To validate the proposed dynamic model in Section 3.2, experiments have been conducted with
a free-swimming robotic fish prototype as shown in Fig. 3.3. The robot has a simple mechanism
for actuation, and it satisfies most of the assumptions used in the modeling work, which facilitates
model validation. Through a chain transmission mechanism, a servomotor (HS-5085MG from
Hitec) was used to control the angular position of the tail shaft and thus the tail deflection angle

α . The tail was a rectangular carbon-fiber foil, which was 8 cm long, 2.5 cm wide (deep), and
1.1 mm thick. With an onboard microcontroller, the servomotor was programmed to rotate the tail

30

according to

α (t) = α0 + αA sin(ωα t + φα )

(3.25)

with α0 , αA and ωα denoting the bias, amplitude, and frequency of the tail beat, respectively; and

φα denoting the phase that indicates the initial condition. When the bias α0 = 0, the robot performs
forward swimming. When α0 = 0, the trajectory of the robot converges to a circular orbit. For
model validation, we measured the turning radius (radius of the circular orbit) and turning period
(time for completing a turn) during steady turning, under different combinations of α0 , αA , and

ωα . The experiments were conducted in an indoor kiddie pool (as shown in Fig. 7 of [106]), and
the turning radius and period were measured with a timer and yardstick, respectively. Specifically,
for each experiment, we would first wait until the robot entered steady turning, and then place
a yardstick (suspended slightly higher than the top of the robot) along the line connecting the
visually determined leftmost and rightmost points of the trajectory. Afterwards (while the robot
continued turning) we would mark down the points A and B where the robot entered and left,
respectively, the half-space defined by the yardstick. These two points were typically close to
the visually determined left/rightmost points mentioned earlier. The distance between A and B
and the difference between the times when the robot hit these two points (captured by a timer)
were used to determine the turning radius and period, respectively. With any prescribed open-loop
control pattern, the observed circular orbit shifted slightly over time, which could be due to the
flow disturbance created by the swimming robot in the relatively small pool (diameter 1.75 m).
To mitigate the impact of such shift on the measurement error, for each pattern, experiments were
repeated for 6 times (restarting the robot from rest every time), from which the mean and standard
deviation of the turning radius/period were evaluated. We note that, while the turning experiments
conducted in this work did offer valuable insight on how tail-beat parameters could impact turning

31

behavior, their primary use is in providing data for model validation under a wide range of openloop inputs. This is different from [113], where the authors focused on turning control strategies
for multi-link robotic fish.

Servo motor

Chain transmission

Tail
Tail shaft

Figure 3.3: A free-swimming robotic fish prototype used for model validation.

3.3.2 Parameter identif cation
The parameters in simulation are measured directly or calculated based on measurement as shown
1 m (2c)2 . The added masses
in Table 3.1. The body inertia about the z−axis is evaluated as Jbz = 12
b

and inertia are calculated by approximating the robot body as a prolate spheroid [12, 112].
The drag and lift coefficients CD , CLβ , and KD can be potentially obtained by performing tests
in a water tunnel or from the CFD simulation, while in this work, we choose a more efficient way to
empirically identify these coefficients. In particular, we have tuned these parameters to match the
turning radius, turning period, and the average of |β | obtained in simulation with the experimental
measurement under a particular tail beat pattern, with bias α0 = 20◦ , amplitude αA = 15◦ , and
frequency ωA = 1.8π rad/s (0.9 Hz). The resulting coefficients are CD = 0.2231, CLβ = 4.50, and
32

Table 3.1: Parameter values for simulations.
Component
Body

Parameter
mass (mb )
distance from body center to the shaft (c)
nominal body surface area (S)
inertia (Jbz )
−XV˙
Cx
−YV˙

Value
0.311 kg
0.07 m
0.0108 m2
5.0797 × 10−4 kg·m2
0.0621 kg
0.2299 kg

−Nω˙ z
tail length (L)
tail depth (d)
water density (ρ )

1.0413 × 10−4 kg·m2
0.08 m
d = 0.025m
1 × 103 kg/m3

Cy

Tail

KD = 11.29 × 10−4 kg·m2 . These parameters are then used in independent model validation for
all other tail beat patterns.
For a fully submerged body taking a shape of a prolate spheroid, which we will assume in this
work, the added-mass terms can be evaluated analytically. Let the profile of the body be described
by
x2 y2 z2
+ +
= 1,
a2 b2 b2

(3.26)

where a and b are the semi-axis length. In particular, considering the robotic fish body volume and
geometric distribution, the values of a and b are set to 7 and 3.3 cm, respectively. Then the added
mass terms are given by

XV˙

= −k1 m f

YV˙

= −k2 m f ,

Cx

Cy

(3.27)

Nω˙ z = −k3 J f z

where positive constant k1 , k2 and k3 are Lamb’s k-factors dependent on the geometry of the body,

33

defined using two constant α0 and β0 :

α0
2(1 − e2 ) 1 1 + e
( ln
− e),
, α0 =
2 − α0
2 1−e
e3
β0
1 1 − e2 1 + e
k2 =
,β = −
ln
,
2 − β0 0 e2
1−e
2e3
e4 (β0 − α0 )
k3 =
,
(2 − e2 )[2e2 − (2 − e2 )(β0 − α0 )]
k1 =

where e =

(3.28)

1 − ( ba )2 , designates the eccentricity of the ellipsoid. Furthermore, in Eq. (3.27), the

mass of the displaced fluid is evaluated using m f = 4/3ρπ ab2, which is 0.3193 kg and the moment
of inertia of the fluid’s spheroid about the z−axis is J f z = 1/5m f (a2 + b2 ).

3.3.3 Results on model validation
Figs. 3.4 – 3.6 show the comparison between the model prediction and experimental measurement
of the turning radius and turning period, when we vary the tail beat bias α0 , amplitude αA , and
frequency ωα , respectively, while holding the other two parameters constant.
Experimental results in Fig. 3.4 show that both the turning radius and the turning period drop
monotonically with an increasing tail beat bias; in other words, the robot makes tighter and quicker
turns as the bias is increased. For these experiments, the tail beat amplitude and frequency are held
fixed at 15◦ and 0.9 Hz, respectively. From Fig. 3.4, these trends are well captured by the model
prediction, which closely matches the experimental measurement.
Fig. 3.5 shows how the turning radius and period depend on the tail beat amplitude, where the
bias and frequency are set to be 20◦ and 0.9 Hz, respectively. The model predicts that the turning
radius is nearly constant (about 27.5 cm) over the tested range for the tail beat amplitude. The
mean of the measured radii for each given amplitude stays within the band [26, 30.5] cm, which is

34

Turning radius (cm)

30

Experiment
Model

25
20
15
10
5
0

20

25

30

35

40

Bias (°)

Turning period (s)

(a)

Experiment
Model

60

40

20

0

20

25

30

35

40

Bias (°)
(b)
Figure 3.4: Comparison between model prediction and experimental measurement of (a) turning radius

and (b) turning period as a function of the tail beat bias α0 . The amplitude and frequency of the tail beat are
fixed at 15◦ and 0.9 Hz, respectively.

35

close to the predicted value. The turning period is predicted by the model to decrease (i.e., faster
turning) as the tail beat amplitude increases, which is supported by the experimental data both
qualitatively and quantitatively.
The impact of the tail beat frequency on the turning behavior is shown in Fig. 3.6, where the
chosen tail beat bias and amplitude are 20◦ and 15◦ , respectively. The model predicts that the
turning radius has negligible dependence on the tail frequency, while the turning period drops with
the frequency. As seen in Fig. 3.6, both trends are well supported by the experimental data.

3.3.4 Hydrodynamic forces with/without incorporating robotic f sh body motion
It is often tempting to consider the body as anchored when evaluating fin-produced hydrodynamic
forces [12, 14, 45], since this would greatly simplify the model. It is of interest to assess the validity of such simplifications for a free-swimming robot where the robot body undergoes movement.
We will use the proposed model in Section 3.2, which has proven effective as shown in the previous subsection, to evaluate the necessity of incorporating the body motion, by comparing the
simulation results when such motion is and is not considered.
For this study, we take the tail beat bias α0 to be 0, and thus the robot will perform forward
swimming at the steady state. We will evaluate the impact of incorporating (or neglecting) body
motion on the computed average thrust (over a tail-beat cycle) and the robot swimming speed. The
thrust force is evaluated via

Fthrust = Fhx cos(ψ − ψ0 ) − Fhy sin(ψ − ψ0 ),

(3.29)

where we recall from Fig. 3.1 that ψ denotes the heading of the robotic fish with respect to the
36

Turning radius (cm)

Experiment
Model

40
30
20
10
0

15

20

25

30

Amplitude (°)

Turning period (s)

(a)

Experiment
Model

60

40

20

0

15

20

25

30

Amplitude (°)
(b)
Figure 3.5: Comparison between model prediction and experimental measurement of (a) turning radius
and (b) turning period as a function of the tail beat amplitude. The bias and frequency of the tail beat are
fixed at 20◦ and 0.9 Hz, respectively.

37

Turning radius (cm)

40

Experiment
Model

30
20
10
0

1

1.5

2

2.5

3

Frequency (Hz)

Turning period (s)

(a)

Experiment
Model

60

40

20

0

1

1.5

2

2.5

3

Frequency (Hz)
(b)
Figure 3.6: Comparison between model prediction and experimental measurement of (a) turning radius
and (b) turning period as a function of the tail beat frequency. The bias and amplitude of the tail beat are
fixed at 20◦ and 15◦ , respectively.

38

X −axis in the inertial frame. In (3.29), ψ0 denotes the average of ψ over a tail-beat cycle at the
steady state, while Fhx and Fhy are evaluated using (3.16) and (3.17). Ignoring the body motion in
hydrodynamic force evaluation implies that the velocity Vτ at point τ on the tail is changed from
(3.13) to
˙ n.
ˆ
Vτ = ατ

(3.30)

Fig. 3.7 compares the calculated average thrusts under different tail-beat amplitudes and frequencies, when the body motion is and is not incorporated in evaluating Vτ , respectively. It shows
clearly that ignoring the body motion would result in significant over-estimate of the thrust. For
example, when the tail beats at 0.9 Hz with an amplitude of 15◦ , the average thrust is 1.188 mN
when the body motion is ignored, which is 57% larger than the value of 0.755 mN when the body
motion is considered. Fig. 3.8 further compares the computed average swimming speeds under
the two modeling schemes, from which we can see that ignoring the body motion results in an
over-estimate of steady-state speed by about 25%. These results illustrate the profound impact of
the body motion on the tail-generated hydrodynamic force.

3.4 Eff cient Adaptation of Drag Coeff cients
3.4.1 Bias-based adaptation of drag coeff cients
The drag force and moment coefficients CD and KD in (3.19) and (3.21) have been assumed to be
constants. Most (if not all) modeling work on robotic fish has assumed constant drag coefficients
[12, 93, 94, 110]. However, it is known in fluid mechanics that the angle of attack influences the
drag on a body that is moving through a fluid [11, 94, 108, 109]. For robotic fish, one challenge is
that the angle of attack oscillates constantly because of the recoil motion of the body in response to

39

Average of thrust force (N)

−3

6

x 10

Considering body motion
Ignoring body motion

5
4
3
2
1
0

15

20

25

30

Amplitude (°)

Average of thrust force (N)

(a)

0.015

Considering body motion
Ignoring body motion

0.01

0.005

0

1

1.5

2

2.5

3

Frequency (Hz)
(b)
Figure 3.7: Simulated average thrust force with/without incorporating the body motion: (a) Thrust force
versus tail-beat amplitude (frequency fixed at 0.9 Hz); (b) thrust force versus tail-beat frequency (amplitude
fixed at 15◦ .

40

Average of velocity (m/s)

Considering body motion
Ignoring body motion

0.08
0.06
0.04
0.02
0

15

20

25

30

Amplitude (°)

Average of velocity (m/s)

(a)

Considering body motion
Ignoring body motion
0.1

0.05

0

1

1.5

2

2.5

3

Frequency (Hz)
(b)
Figure 3.8: Simulated average velocity with/without incorporating the body motion: (a) Velocity versus
tail-beat amplitude (frequency fixed at 0.9 Hz); (b) velocity versus tail-beat frequency (amplitude fixed at
15◦ .

41

tail movements. Furthermore, evaluation of the angle of attack involves a highly nonlinear function
of the velocity components, β = arctan(VCy /VCx ). In this section, we present an approach to the
efficient adaptation of drag coefficients for robotic fish, to maintain or even improve model fidelity
without the need to explicitly evaluate the angle of attack. The word “efficient” here refers to
computational efficiency (or low modeling complexity), as we will show in the later part of this
section that such a varying drag effect can be well captured with a simple function of the tail-beat
bias.
At the core of our approach is the exploitation of the close correlation between the angle of
attack and the tail beat bias. Recall the form of the tail beat angle α in (3.25). While the bias α0 ,
amplitude αA , and frequency ωα all could be varied over time (for example, when one controls
the robot to follow some arbitrary path), such variations typically take place at a much slower time
scale than the tail beat motion – a tail beat pattern will be repeated for at least several cycles before
it is updated. Figs. 3.4–3.6 suggest that the turning radius strongly depends on the bias α0 , but
is nearly independent of the amplitude αA and frequency ωα . The turning radius is intimately
related to the angle of attack. To explicitly reveal the relationships between the angle of attack and
the tail-beat pattern parameters, we have used the model identified in Section 4.5 to numerically
evaluate the average β0 of the absolute value |β | of the angle of attack over each beat cycle, i.e.,

β0 =

ωα t0 +2π /ωα
|β (t)|dt.
2π t0

(3.31)

From Fig. 3.9, β0 shows nearly linear dependence on the tail-beat bias angle, with little dependence
on the amplitude or the frequency.
Based on the observations from Fig. 3.9 and on the relationships between the drag coefficients
and the angle of attack [108], we propose adapting these coefficients using the tail-beat bias α0

42

0.2

β0

0.15
0.1
0.05
0

20

25

30

35

40

Bias (°)

(a)
0.2

β0

0.15
0.1
0.05
0

15

20

25

30

Amplitude (°)

(b)
0.2

β0

0.15
0.1
0.05
0

1

1.5

2

2.5

3

Frequency (Hz)

(c)
Figure 3.9: Simulation results on the dependence of the angle of attack on tail beat parameters, where β0

is the cycle-averaged absolute value of the angle attack: (a) β0 versus the tail beat bias α0 ; (b) β0 versus
the tail beat amplitude αA ; (c) β0 versus the tail beat frequency ωα . In (a), the amplitude and frequency are
fixed at 15◦ and 0.9 Hz; in (b), the bias and the frequency are fixed at 20◦ and 0.9 Hz; in (c), the bias and the
amplitude are fixed at 20◦ and 15◦ .

43

(note that the underlying assumption here is that, the impact of the angle of attack, β , on drag can
be well captured through the cycle-averaged absolute value of β , β0 ):

′ +C ′ α 2 ,
CD =CD0
D1 0

(3.32)

′ + K′ α ,
KD =KD0
D1 0

(3.33)

′ , C ′ , K ′ and K ′ are constants that are empirically identified as follows. First, as
where CD0
D1
D0
D1

described in Section 3.3.2, we seek values of CD , KD , and CLβ such that the simulated turning
radius and period best match the experimental results when the tail beats at 0.9 Hz with amplitude
of 15◦ and bias of 20◦ . We then maintain the value of CLβ , and find the best parameter estimates
for CD and KD when the bias is 30◦ and 40◦ , respectively. We have deliberately left out the cases
of 25◦ and 35◦ biases at this stage, so that those data can be used in independent model validation.
Table 3.2 lists the identified CD and KD identified for the three biases. We then use the values
′ , C ′ , K ′ and K ′ in (3.32) and (3.33) through the least-squarein Table 3.2 to identify CD0
D1
D0
D0
′ = 0.1936, C ′ = 0.1412/rad2 , K ′ = 9.655 ×
error fitting. The resulting parameters are CD0
D1
D0
′ = 4.383 × 10−4 kg·m2/rad, respectively.
10−4 kg·m2 and KD1

Fig. 3.10 shows the comparison between the measured turning radii and periods and the values
predicted by the model with drag coefficient adaptation. We can see that, with the bias-based
adaptation of drag coefficients, the match between the model prediction and the experimental data
is further improved over that in Fig. 3.4.
Table 3.2: Identified drag force and moment coefficients for different tail beating biases.
Bias α0 (◦ )
CD
KD (kg·m2 )
20
0.2231 1.129 × 10−3
30
0.2112 1.174 × 10−3
40
0.2712 1.282 × 10−3

44

Turning radius (cm)

30

Experiment
Model

25
20
15
10
5
0

20

25

30

35

40

Bias (°)

Turning period (s)

(a)

Experiment
Model

60

40

20

0

20

25

30

35

40

Bias (°)
(b)
Figure 3.10: Comparison between experimental measurement of (a) turning radius and (b) turning period
with those predicted by the model, where the drag coefficients are adapted and the lift force is evaluated
using (3.20). The amplitude and frequency of the tail beat are fixed at 15◦ and 0.9 Hz, respectively.

45

To further validate the approach for drag coefficient adaptation, we have conducted additional
experiments involving transients in forward swimming. A yellow LED powered by a button cell
was taped to the top of the robot as a marker. The weight of the LED/cell add-on was 3 g, so its
impact on the total mass (over 300 g) and inertia of the robot was negligible. Starting from at rest,
the robot was programmed to swim forward (with the tail beat bias α0 = 0◦ ) with the amplitude
and frequency of the tail-beat equal to 15◦ and 0.9 Hz, respectively. We video-recorded the motion
of the LED (and thus of the robot) with the ambient light turned off, and extracted the robot’s
time-trajectory by applying background subtraction to the images. The dynamic model has been
simulated with the adapted and non-adapted drag coefficients, respectively. For the model with
non-adapted drag coefficients, CD = 0.2231, KD = 11.29 × 10−4 kg·m2 , which are identical to the
coefficients obtained from the turning experiments as presented in Section 3.3.2. For the model
with adapted drag coefficients, CD = 0.1936, KD = 9.655 × 10−4 kg·m2 , which are evaluated
based on (3.32) and (3.33), with α0 = 0. Fig. 3.11(a) and (b) compare the model predictions of the
X /Y -coordinate time-trajectories of the robot with the experimental measurement, and Fig. 3.11(c)
provides the comparison on the robot path in the XY -plane. It can be seen that the model with
adapted drag coefficients outperforms the one with non-adapted coefficients, and is able to capture
well both the transient and steady-state behaviors of the robot.

3.4.2 Elimination of the angle of attack from the lift force term
Note that, while the drag coefficients in (3.32) and (3.33) do not involve the angle of attack β ,
the evaluation of the lift force for the robot body in (3.20) still requires knowing β . In order to
eliminate β from the overall model and thus reduce the model complexity, in this subsection we

46

Experiment
Adapted model
Non−adapted model

1

X (m)

0.8
0.6
0.4
0.2
0
0

10

20

30

Time (s)

(a)

Experiment
Adapted model
Non−adapted model

0.25

Y (m)

0.2
0.15
0.1
0.05
0
0

10

20

30

Time (s)

(b)
0.4

Experiment
Adapted model
Non−adapted model

Y (m)

0.3
0.2
0.1
0
0

0.2

0.4

0.6

0.8

X (m)

(c)
Figure 3.11: Comparison between model predictions and experimental measurement for forward swimming (including transients): (a) time trajectory of X -coordinate of the robot; (b) time trajectory of the
Y -coordinate of the robot; (c) path of the robot in the XY -plane. For both the experiment and simulations,
the amplitude and frequency of the tail beat are fixed at 15◦ and 0.9 Hz, respectively.

47

explore replacing β in (3.20) with the tail-beat bias α0 :
1
FL = ρ |VC |2 SCL′ β α0 ,
2

(3.34)

where, CL′ β is the new lift coefficient associated with the tail bias. To identify CL′ β , we tune CL′ β ,
CD , and KD so that the simulated turning radius, turning period, and average of |β | match those
obtained in experiments when the tail beats with bias α0 = 20◦ , amplitude αA = 15◦ , and frequency

ωA = 1.8π rad/s. The estimates are CL′ β = 0.935, CD = 0.2131, and KD = 1.144 × 10−3 kg·m2.
With CL′ β fixed, we can then follow the similar steps as in Section 3.4.1 to determine the drag
coefficient adaptation laws (3.32) and (3.33). Note that the four parameters in (3.32) and (3.33)
are identified twice since we want to compare the two cases when the evaluation of the lift force
involves the angle of attack β or its surrogate α0 . In practice, one only needs to identify these
parameters once given the choice of (3.20) or (3.34) for the lift force evaluation.
Fig. 3.12 and Fig. 3.13 show the model prediction performance when drag coefficients are fixed
and adapted, respectively, where the lift force is evaluated using (3.34) instead of using (3.20) as
for Figs. 3.4-3.6 and 3.10. From Fig. 3.12, the model with the modified lift term (3.34) is still able
to capture the trends of how the turning radius and period vary with the tail bias, although the match
has degraded from that seen in Fig. 3.4. From Fig. 3.13, one can see that, with the drag adaption,
significantly better match between the model prediction and experimental data is achieved.

3.5 Conclusion
In this work, we have presented a complete dynamic model for robotic fish actuated by a caudal
fin, where Lighthill’s elongated-body theory is used to evaluate the tail-generated hydrodynamic

48

Turning radius (cm)

30

Experiment
Model

25
20
15
10
5
0

20

25

30

35

40

Bias (°)

Turning period (s)

(a)

Experiment
Model

60

40

20

0

20

25

30

35

40

Bias (°)
(b)
Figure 3.12: Comparison between experimental measurement of (a) turning radius and (b) turning period
with those predicted by the model, where the drag coefficients are not adapted and the lift force is evaluated
using (3.34). The amplitude and frequency of the tail beat are fixed at 15◦ and 0.9 Hz, respectively.

49

Turning radius (cm)

30

Experiment
Model

25
20
15
10
5
0

20

25

30

35

40

Bias (°)

Turning period (s)

(a)

Experiment
Model

60

40

20

0

20

25

30

35

40

Bias (°)
(b)
Figure 3.13: Comparison between experimental measurement of (a) turning radius and (b) turning period
with those predicted by the model, where the drag coefficients are adapted and the lift force is evaluated
using (3.34). The amplitude and frequency of the tail beat are fixed at 15◦ and 0.9 Hz, respectively.

50

force. For experimental validation purpose, we have considered a rigid tail with biased sinusoid
oscillation and treated the turning radius and turning period as outputs in the comparison between
model prediction and experimental results. Additional model validation has also been conducted
with experimental data of forward swimming including the transients. We have found that, incorporating the body motion in the evaluation of hydrodynamic forces plays a significant role in
capturing the dynamics. We have also shown that, the drag force and moment coefficients are
strongly dependent on the tail-beat bias, and subsequently we have proposed a novel scheme to
adapt these coefficients and lift force based on the tail movement pattern. The effectiveness of this
approach has been further demonstrated in experiments.

51

Chapter 4
Dynamic Modeling of Robotic Fish with a
Base-actuated Flexible Tail
4.1 Introduction
Numerous designs of actuation mechanisms have been proposed for robotic fish [5, 11, 12, 19,
24, 35–39, 41, 43–46, 49, 114]. A typical approach is to use multiple links actuated separately or
jointly to deform the body itself [5, 8, 16, 19, 115, 116], which requires multiple actuators and/or
complex transmission mechanisms. Another class of designs involves an oscillating caudal fin
(e.g., [6,11,12,44,45] ), sometimes in conjunction with pectoral fins [11]. Among various designs,
tail actuation is especially attractive since it is easy to realize, enables both forward swimming
and turning maneuvers, and leaves the majority of the body free of moving parts. The latter is
important when the robot is used in applications such as mobile sensing, where the body space can
be maximally used to house sensors and electronics [1]. Many tail-actuation designs involve rigid,
oscillating plates [11,14,15,93]. This design lends itself to simple construction and tractable analysis. However, it has been recognized that the flexibility of body and fin structures has a pronounced
impact on the swimming performance of biological and robotic fish [18, 117, 118]. Flexible caudal
fins can be realized by motor-driven compliant beams or plates [19, 44], or directly through soft
actuation materials such as ionic polymer-metal composites (IPMCs) [6, 9, 12, 45]. While these
active materials possess intriguing properties, the thrust they can produce is still relatively weak
52

and the long-term repeatability of their actuation behavior is yet to be established. Therefore,
base-actuated soft, passive structures remain a competitive option for flexible tails.
To design and control robotic fish actuated with a flexible tail, it is essential to have a faithful
and efficient dynamic model. Modeling of a flexible beam attached to a moving base in the air has
been studied by a number of researchers [119–121]. However, in this work, a major challenge in the
modeling lies in properly capturing the fluid-structure interactions and the resulting force and moment on the robot. Most existing work on modeling has dealt with rigid fins [11, 14, 15, 90, 92, 93],
although modeling of flexible fins has been conducted recently by several groups [12, 19, 45, 122].
In [19], Alvarado et al. focused on designing compliant bodies to achieve biomimetic locomotion
efficiency and maneuverability. To demonstrate the results, they also compared the experimental
results with the linear beam model. Researchers have also developed models that capture the actuation physics of an anchored IPMC beam and the complex hydrodynamic interactions between
IPMC and fluid [12, 45]. In [122], Kopman et al. described a modular biomimetic robotic fish
developed for educational activities, along with a modeling framework for predicting the robot’s
static thrust production. In all these studies, a linear Euler-Bernoulli beam model was adopted
to describe the beam dynamics. A disadvantage of this approach is that, under large oscillations,
these models do not accurately capture the beam dynamics [123–125] and consequently the hydrodynamics.
In this work, we take a significant departure from the beam model approach by approximating a flexible tail with multiple rigid segments connected in series through torsional springs and
dampers. A similar multi-segment approach has been used to model a beam under large deformation in air [125,126] and model a flexible beam in flow sensing [127]. However, the multi-segment
beam model has not been explored in the modeling of a base-actuated flexible beam subject to
complex hydrodynamic interactions. For comparison, we also present a model based on linear
53

beam theory, as widely adopted in other studies [12, 19, 45, 122]. Note that, for ease of presentation, a tail with uniform height is considered in this work. Despite its simple appearance, this
case captures the key essential challenges (nonlinear beam dynamics under large deformation,
coupled with hydrodynamic interactions) that would be present for cases with more general tail
shapes. We evaluate the hydrodynamic force on the actuating tail using Lighthill’s large-amplitude
elongated-body theory [25], because it has a sound balance between fidelity and simplicity and its
effectiveness in robotic fish modeling has been demonstrated in our prior work [15, 45]. A weakness of Lighthill’s theory is that it neglects the effect of vortex wake on the pressure distribution
on the body and thus neglects the effect of vortex shedding on the thrust production [128]. Incorporating such vortex-shedding effects (for example, using the vortex ring panel method [129]),
however, typically requires CFD simulation and is not amenable to the development of a model
that is suitable for robot optimization and control, which is the goal of this work.
To compare and validate the models, we conduct extensive experiments with a robotic fish
prototype. We find that both models have similar predictions that are close to the experimental
measurements when the flexible beam is under small deformations. However, the discrepancy
between the predictions becomes greater as the beam experiences larger deformations. In particular, we show that the multi-segment model is able to predict precisely the transient trajectory and
steady-state speed of the free-swimming robot, as well as the dynamic shape of the tail, under a
wide range of actuation inputs. In comparison, the linear beam model-based approach is only able
to capture the robot speed and tail shape when the tail undergoes relatively small deformations.
The remainder of the chapter is organized as follows. The overall dynamic model for a tailactuated robotic fish is first reviewed in Section 4.2, where the hydrodynamic forces generated
by the tail are evaluated by Lighthill’s theory. The use of Lighthill’s theory requires knowing
the tail shape as a function of time, which is the main focus of this chapter. The Euler-Bernoulli
54

beam-based model and the multi-segment model for the based-actuated flexible tail are developed
in Sections 4.3 and 4.4, respectively. Experimental validation and comparison of the models are
presented in Section 4.5. Finally, concluding remarks are provided in Section 4.6.

4.2 The Model Framework for Tail-actuated Robotic Fish
The robot is assumed to comprise two parts, a rigid body and a flexible tail. The motion of the robot
body is governed by rigid-body dynamics with the added-mass effect incorporated. Lighthill’s
large-amplitude elongated-body theory is adopted to evaluate the hydrodynamic forces generated
by the flexible tail’s motion. We assume that the height of the tail does not vary abruptly along its
length, thus meeting the “elongated body” requirement [25].
Fig. 4.1 shows a schemata of the robotic fish, with [X ,Y, Z] denoting the inertial coordinates,
and [x, y, z] denoting the body-fixed coordinates with unit vectors [x,
ˆ y,
ˆ zˆ]. We denote by mˆ and
nˆ the unit vectors parallel and perpendicular to the tail, respectively. We assume that both the
body and the tail are neutrally buoyant, and that the center of gravity of the body coincides with
the center of geometry at point C. The velocity at C expressed in the body-fixed coordinates
VC = [VCx ,VCy ,VCz ]T comprises surge (VCx ), sway (VCy ) and heave (VCz ) components. In addition,
the angular velocity ω = [ωx , ωy , ωz ]T comprises roll (ωx ), pitch (ωy ) and yaw (ωz ) expressed in
the body-fixed coordinates. We use α to denote the tail deflection angle (the tangential direction
of the flexible tail at the base) with respect to the negative x−axis, and β to denote the angle of
attack, formed by the direction of VC with respect to the x−axis. Finally, ψ denotes the heading
angle, formed by the x−axis relative to the X −axis.
We note that the actual mass of the tail is negligible (relative to the virtual mass) when the
width of the tail is much thinner than the height [25], which is typically true for the tail of a robotic

55

Y

yˆ

MD

y

Zz
FL

FD

xˆ

&
VC

x

B

E
\

Fish body
z

C
A

Tail

D

nˆ

Z

D

mˆ

X

x

Figure 4.1: Schematic representation of the robotic fish in planar motion.
fish. For example, in the model validation part of this work, we use a flexible tail with a height of
2.5 cm and a thickness of 0.3 mm, which indicates that mw is 490.9 g per meter. In contrast, the
physical tail, which has a length of 8 cm, weighs 0.54 g, resulting in a mass of 6.75 g per meter,
which is 1.38% of the virtual mass.

4.3 Dynamic Modeling of a Tail Using Euler-Bernoulli Beam
Theory
The application of Lighthill’s theory requires knowing the shape trajectory of the base-actuated
flexible tail. We propose a multi-segment approach to the modeling of the tail that undergoes
large deformations in Section 4.4. But in this section we first develop a comparative approach
based on linear beam theory, which has been widely used in the relevant literature. We focus on
56

beam modeling with the robot anchored, as typically adopted in the literature [19, 122]. Fig. 4.2
illustrates this approach; the dashed line represents a rotating frame, which coincides with the line
tangential to the tail at the base, and the blue solid curve represents the shape of the flexible caudal
fin. The transverse displacement of points on the flexible tail relative to the rotating base, due
to the beam’s vibration, is given by w(τ ,t). We note that the tail is assumed to be under small
deformation in accordance with the Euler-Bernoulli theory. In other words, the motion direction
of each point along the flexible tail is perpendicular to the dashed line.

y

Y

x
B

Fish body
z

x

C
A
D (t )

Caudal Fin
w(W , t )

D

Z

X

x

Figure 4.2: Passive flexible tail modeled by Euler-Bernoulli linear beam approach.

The tail is excited by the oscillating support at the point A, with:

∂ w1 (τ ,t)
∂ w2 (τ ,t)
=
= α (t).
∂τ
∂τ
τ =0
τ =0
57

(4.1)

Here w1 (τ ,t) denotes the beam displacement in the inertial frame and w2 (τ ,t) accounts for the
rigid-body (rotating base) motion:

w1 (τ ,t) = w2 (τ ,t) + w(τ ,t),
w2 (τ ,t) = τ · α (t).

(4.2)

The forced underwater vibration of the flexible beam is described by the following equation
[122, 130]

∂2
∂ 2 w1 (τ ,t)
∂ 2 w1 (τ ,t)
EI(
τ
)
+
m(
τ
)
= p(τ ,t),
∂ τ2
∂ τ2
∂ t2

(4.3)

where E denotes the Young’s modulus of the flexible tail, I denotes the area moment of inertia,
m(τ ) denotes the mass of the tail per unit length and p(τ ,t) denotes the transverse loading on the
tail which is caused by the interactions between the tail and the surrounding aquatic environment
that can be evaluated by Lighthill’s theory. Unlike [19, 122], in this work we also consider distributed viscous damping introduced by the internal resistance opposing the strain velocity [130],
which leads to

2
3
∂ 2 EI(τ ) ∂ w1 (τ ,t) + κ ∂ w1 (τ ,t)
∂ τ2
∂ τ2
∂ τ 2∂ t

+ m(τ )

∂ 2 w1 (τ ,t)
∂ t2

= p(τ ,t),

(4.4)

where κ is the stiffness proportionality factor for Rayleigh damping.
Similar to [14, 15], we assume that the tail itself has negligible mass compared to the added
mass effects. Then following Lighthill’s theory, we obtain:

EI

∂ 5 w(τ ,t)
∂ 4 w(τ ,t)
∂ 2 w(τ ,t)
+
κ
+
m
= Peff (τ ,t)
w
∂ τ4
∂ τ 4∂ t
∂ t2
58

(4.5)

in which

Peff (τ ,t) =

∂ 4 w2 (τ ,t)
∂ 5 w2 (τ ,t)
+
κ
∂ τ4
∂ τ 4∂ t
∂ 2 w2 (τ ,t)
−mw
∂ t2
−EI

represents the effective distributed dynamic loading caused by the prescribed support excitations,
and mw is virtual mass per unit length defined previously, accounting for the added mass effect.
Since the first two terms of Peff (τ ,t) are equal to zero as w2 (τ ,t) = τ · α (t), the tail beam can be
modeled as
EI

∂ 4 w(τ ,t)
∂ 2 w(τ ,t)
∂ 5 w(τ ,t)
+
κ
+
m
= −mw τ · α¨ (t).
w
∂ τ4
∂ τ 4∂ t
∂ t2

(4.6)

Using the modal analysis method, we can express the solution of (4.6) as the sum of an infinite
number of modes:
w(τ ,t) =

∞

∑ φi(τ )ηi(t),

(4.7)

i=1

where φi (τ ) is the beam shape for the ith mode and ηi (t) is the corresponding generalized coordinate. There is no concentrated force in the nˆ direction with the anchored body assumption, which
implies that the boundary conditions for w(τ ,t) is the same as those for a cantilever beam:

w(0,t) = 0,
EI

∂ 2w
(L,t) = 0,
∂ τ2

∂w
(0,t) = 0,
∂τ
∂ 3w
EI 3 (L,t) = 0.
∂τ
(4.8)

59

The mode shape φi (τ ) takes the form

φi (τ ) = (cos βi τ − cosh βi τ ) − δi (sin βiτ − sinh βi τ ),

(4.9)

where βi can be obtained by solving

1 + cos βi L cosh βi L = 0,

and

δi =

cos βi L + cosh βi L
.
sin βi L + sinh βi L

With the damping ratio for the ith mode

ξi =

κωi
,
2

where ωi is the natural frequency for the ith mode,

ωi = βi2

EI
,
mw

(4.10)

ηi (t) can be solved from
d2 ηi (t)
dη (t)
Q (t)
+ 2ξi ωi i + ωi2 ηi (t) = i ,
2
dt
Mi
dt

i = 1, 2 . . .

(4.11)

where
Qi (t) =

L
0

φi (τ )Peff(τ ,t)dτ ,

60

(4.12)

and
Mi =

L
0

φ i ( τ ) 2 mw d τ .

(4.13)

When the tail shaft oscillates sinusoidally, we can obtain a closed-form solution for ηi (t). Consider
in particular the following form for the base angle

α (t) = α0 + αA sin(ωα t + φα ),

(4.14)

where α0 , αA , ωα , and φα denote the bias, amplitude, frequency, and initial phase of the tail beat,
respectively. Then we have
mw αA ωα2 0L τφi (τ )dτ
ηi (t) =
sin(ωα t + φα + φi ),
MiZi ωα

(4.15)

where
1
(2ωi ξi )2 + 2 (ωα2 − ωi2 )2
ωα

Zi =

is the magnitude of the impedance function and

φi = arctan

2ωα ωi ξi
ωα2 − ωi2

is the phase lag of the oscillation relative to α (t). The total hydrodynamic force acting on the tail
is then
L

F=

0

f (τ )d τ −

1
mw v2 mˆ
,
⊥
2
τ =L

(4.16)

2
where f (τ ) = −mw (τ α¨ + ∂ 2w )nˆ − α˙ (τ α˙ + ∂∂wt )mˆ . The corresponding moment relative to the

∂t

61

center C of the robot body is
L

M=

0

rCτ × f (τ )d τ + rCL ×

−

1
mw v2 mˆ
,
⊥
2
τ =L

(4.17)

where rCτ and rCL denote the vectors from the body center C to the point τ and L on the tail,
respectively.

4.4 Dynamic Modeling of a Tail Using Multi-segment Approximation
A fundamental underlying assumption in establishing the model described in the previous section
is that the flexible tail is under small deformation. However, based upon observation of the dynamic
shape of the tail, this assumption does not always hold, especially when the tail undergoes large
angular displacement and/or high-frequency oscillations. In this section, we propose a novel model
by representing the tail as multiple rigid elements connected in series through torsional springs and
dampers, to capture the large deformation of the beam. N rigid segments, with equal length of l,
are used to represent the tail, as illustrated in Fig. 4.3. Each segment is linked with its neighboring
segments through joints modeled by a torsional spring KS and a viscous damper KD .
Following [125], we can evaluate the stiffness of each torsional spring as

KS =

Edh3
,
12l

(4.18)

with h denoting the thickness of tail. KD can be evaluated as KD = κ KS , where κ is the proportional
constant as defined in the previous section. Following Lighthill’s large-amplitude elongated-body

62

y

Zz

Y

&
VC

x

E

Fish body

\

z

B

C

D1
Di

A0
A1

Ai 1

nˆi

Ai

mˆ i

AN 1

Caudal Fin

Z

x

AN D

X

Figure 4.3: Passive flexible tail modeled by multiple rigid segments.
theory, we need to compute the motion of every point along the tail over time, in order to evaluate
the tail actuation-induced hydrodynamic force. Consequently, we need to know the joint angles

αi , made by the ith link with respect to the negative x-axis (as illustrated in Fig. 4.3), so that (3.8)
and (3.9) can be applied to evaluate the hydrodynamic forces on the tail section. The displacement
rτi of every point τi on the ith segment in the inertial frame can be described as:
i−1

rτi = l ·

∑ mˆ k

+ τi mˆ i ,

(4.19)

k=1

then the velocity perpendicular to the ith segment is
i−1

vτi ⊥ = l ·

∑ α˙ k cos(αk − αi)

k=1

63

+ τi α˙ i .

(4.20)

Therefore, the force density acting on the ith segment can be evaluated as

fτi = −mw

d
(v nˆ )
dt τi ⊥ i
i−1

= −mw l

∑

α¨ k cos(αi − αk )

k=1

−α˙ k sin(αi − αk )(α˙ i − α˙ k ) + τi α¨ i nˆ i
+mw α˙ i l

i−1

∑

α˙ k cos(αi − αk ) + α˙ i τi mˆ i .

(4.21)

k=1

The total reactive force Fi on the ith segment, and the moment Mi relative to point Ai−1 can be
evaluated using
l

Fi =

fτi d τ ,

l

Mi =

τi mˆ i × fτi d τ .

(4.22)

Defining FAi and MAi to be the force and the moment exerted by the (i + 1)th segment on the ith
segment, respectively, we can express the interactions between adjacent segments as,

FAi = Fi+1 + FAi+1 ,

(4.23)

MAi = Mi+1 + MAi+1 + l mˆ i × FAi+1 .

(4.24)

For the last segment of the tail, as illustrated in Fig. 4.3, the reactive force and moment are
N
1
2
FAN = − mw l ∑ α˙ k cos(αN − αk ) mˆ N
2
k=1
N

+mw l

∑ α˙ k cos(αN − αk )

·

k=1
N

l

∑ α˙ k sin(αN − αk ) nˆN

(4.25)

k=1

MAN = 0.

(4.26)

64

Defining M(S+D) the moment produced by the spring and damper at joint Ai , then the moment
i

balance equation implies, for i = 1, 2, · · · , N − 1,

MAi = M(S+D) ,

(4.27)

M(S+D) = KS (αi+1 − αi ) + KD (α˙ i+1 − α˙ i ) zˆ.

(4.28)

i

where
i

Eq. (4.27) has (N − 1) scalar equations involving (N − 1) unknown variables α2 , · · · , αN , which is
solvable.
FA0 is the force that the flexible tail exerts on the robotic fish body, which can be written as
ˆ where Fhx and Fhy are the components of FA along xˆ and y,
ˆ respectively. The
FA0 = Fhx · xˆ + Fhy · y,
moment to the center of the body caused by the oscillation of the tail can be evaluated as

MC = MA0 − cxˆ × FA0 .

(4.29)

It is clear that MC is along the z direction, which we denote as Mhz .

4.5 Experimental Model Validation
To evaluate the dynamic models described in Sections 4.3 and 4.4, we have conducted experiments
with the free-swimming robotic fish prototype shown in Fig. 4.4. The robot has a simple mechanism for actuation, and it satisfies most of the assumptions used in the modeling work, which
facilitates model validation. In particular, the height of the tail does not change abruptly along
the length direction, thus meeting the “elongated body” requirement. The thickness of the tail is

65

much less than its height, so that the actual mass of the tail is much less than the virtual mass as
assumed. Through a chain transmission mechanism, a servomotor (HS-5085MG from Hitec) is
able to control the angular position of the tail shaft and thus the tail deflection angle α , i.e. sinusoidal motion, accurately. On the other hand, in the modeling work, we assume that the robot is
anchored, as typically adopted in the literature. This assumption, of course, does not hold fully during free-swimming experiments, which might explain the slight discrepancies between the model
predictions and experimental measurement in those experiments. The tail was a rectangular plastic
slice, which was 8 cm long, 2.5 cm wide (deep), and 0.3 mm thick.

Programming port

Servo motor

Chain transmission

Tail shaft

Flexible Tail

Figure 4.4: A free-swimming robotic fish prototype used for model validation.

4.5.1 Parameter identif cation
The same experimental prototype (with a different tail) was used in [14] to validate a dynamic
model for a robotic fish with a rigid tail. The following parameters for the robot were identified
in [14]: c = 0.07 m, mb = 0.311 kg, ρ = 1000 kg/m3, S = 0.0108 m2 , Jbz = mb (2c)2 /12. The added
masses and inertias are calculated by approximating the robot body as a prolate spheroid [12,112]:

66

−XV˙

Cx

= 0.0621 kg, −YV˙

Cy

= 0.2299 kg, and −Nω˙ z = 1.0413 × 10−4 kg·m2 .

The drag and lift coefficients, CD , CL , and CM are also identified empirically for the robotic
fish with a rigid tail and then used in validating the flexible tail models. In particular, we have
tuned these parameters to match the turning radius, turning period, and the average of |β | obtained
in simulation, with the experimental measurements recorded for a particular tail beat pattern. We
set amplitude αA = 13.6◦, frequency ωα = 1.8π rad/s (0.9 Hz) and evaluate the suitable drag and
lift coefficients for the tail-beat bias α0 equal to 20◦ , 30◦ and 40◦ as [15]. Using least-square-error
fitting, we obtain CD = 0.276, CL = 4.5, and CM = 7.4 × 10−4 kg·m2, when α0 = 0.
For measuring the Young’s modulus, we set up the experiment shown in Fig. 4.5, and evaluate
E using:
E=

L3b PL
,
3Ib wL

(4.30)

where Lb and Ib are the length and the area of moment inertia of the testing beam, respectively. PL
is the load at the tip end, and wL is the end’s displacement. A plastic beam (the same material as
the flexible tail) is clamped at the upper surface of a rectangular block, which is fixed on a piston
that can move up and down along the stand. A custom LabVIEW (2011 SP1) virtual instrument
graphical user interface (GUI) is developed to perform the data acquisition through a dSPACE
system (RTI 1104, dSPACE). The force exerted on the beam is captured by a load cell (GS-10,
Transducer Techniques) with a custom-made amplifier circuit. The displacement of the beam tip
is measured with a laser sensor (OADM 20I6441/S14F, Baumer Electric). Prior to experiment, the
system is calibrated with a weight applied to the load cell. Three sets of data are collected to calculate the Young’s modulus, and in each set, a least square error method is adopted to approximate
the slope between the force and displacement. Fig. 4.6 shows results for three tests, indicating that
E ∈ [1.41, 1.51] GPa. In this work, we takes the average E = 1.48 GPa.

67

Laser displacement sensor

Beam
Load cell
Stand
Clamp
Figure 4.5: Experimental setup for measuring the Young’s modulus of the flexible tail.
When modeling a tail using the multi-segment approximation, the number of rigid segments
affects both modeling accuracy and computational complexity [127]. Specifically, a higher number
of elements results in a more accurate model, but is also more computation-intensive. Fig. 4.7
shows the simulated responses of the beam tip displacement (relative to the x axis along the body)
to a sinusoidal base excitation, when different numbers of segments are used to model the flexible
beam. The properties of the beam used in the simulation are same as those of the beam identified

68

−3

x 10

First Test
Second Test
Third Test

Force (N)

9.5

9

8.5

8

6.2

6.4

6.6

6.8

Tip Displacement (m)

7
−3

x 10

Figure 4.6: Experimental results for measuring the Young’s modulus of the flexible tail.
from the experiments. It can be observed that, the beam tip displacements gradually converge to
each other when the number of segments increases. Fig. 4.8 shows the running time needed for
the simulation illustrated in Fig. 4.7, which was conducted with Matlab/Simulink on a desktop
PC (Dell Vostro 460 with 3.1 GHz Intel i5-2400 CPU and 4 GB memory). In particular, for all
cases, the fundamental sampling time in Simulink is fixed at 0.000167 s, the simulation time is
set to be 10 s in total, and the elapsed CPU time is obtained using the Matlab macro cputime.
We can see that the computation time increases rapidly with the number of rigid segments used
in the simulation. A five-segment approximation achieves a sound tradeoff between the modeling
accuracy and computational efficiency, and therefore is adopted in this study.
According to Eq. (4.18), we have KS = 5.2 × 10−3 N·m that is used in the model using multisegment approximation. The values of κ are identified empirically by fitting the data of forward

69

Five Segments
Four Segments
Three Segments

Tip Displacement (m)

0.05
0.04
0.03
0.02
0.01
0
−0.01
−0.02
−0.03
0

0.5

1

1.5

2

2.5

3

Time (s)
Figure 4.7: Displacement of the tail tip generated by the model using different number of rigid segments.
speed versus beating frequency for both models. In particular, we choose the value of κ such
that the simulated steady-state speed of the robot, when α0 = 0◦ , matches the experimental results
when the tail beats at 0.45 Hz, 0.9 Hz and 1.35 Hz, respectively (for which cases the flexible tail
undergoes relatively small deformations). This results in κ = 0.0846 s for the Euler-Bernoulli
beam model and κ = 0.0731 s for the multi-segment approximation model. We take the average

κ = 0.079 s, which leads to KD = 4.1 × 10−3 N·m·s, using KD = κ KS . These parameters are then
applied in the simulation of all other cases.

4.5.2 Model verif cation
Fig. 4.9 compares the two models and experimental measurements in terms of the steady-state
speed of the robotic fish with respect to the actuation frequency. The trends are similar for the two

70

12

CPU time (s)

10
8
6
4
2
0
2

3

4

5

6

Number of segments)
Figure 4.8: Computation time needed to simulate the model using different numbers of rigid segments.
model predictions. In particular, the robot’s speed increases with the frequency up to a threshold
value (1.35 Hz in the simulated case) and then starts to drop. While the predictions from both
models match the experiments relatively well (due to the tuning process for parameter κ as described in the previous subsection), the discrepancy between the predictions becomes larger as the
actuation frequency increases. The flexible tail is subjected to the added mass effect as explained
in Eq. (3.8). Under a high frequency excitation, the lateral hydrodynamic loading on the tail increases and causes larger deformation of the tail. By using linear Euler-Bernoulli beam theory, we
assume that the tail is under small deformation, and thus the motion direction of each point along
the flexible tail is perpendicular to the dashed line in Fig. 4.2. The latter assumption no longer
holds when the excitation frequency gets high due to the large deformation, which explains the
poor prediction performance of the linear beam model at relatively high frequencies. On the other

71

hand, the model using the multi-segment approximation matches the experimental measurement
closely throughout the actuation frequency range used in the experiments.

12

Experiment
Multi−segment Model
Beam Model

Speed (cm/s)

10
8
6
4
2
0
0

0.5

1

1.5

2

2.5

3

Tail−beat frequency (Hz)
Figure 4.9: Comparison between model predictions and experimental measurement of the speed versus
tail-beat frequency. The amplitude is fixed at 13.6◦ .

To further compare the two models under different conditions, a second experiment has been
conducted with another tail of different dimensions, which is 2.34 cm high and 9.8 cm long. Following Eq. (4.18) and the linear relationship KD = κ KS , the values of KS and KD are updated
from the original KS = 5.2 × 10−3 N·m and KD = 4.1 × 10−3 N·m·s to KS = 4.0 × 10−3 N·m and
KD = 3.13 × 10−3 N·m·s, respectively. Fig. 4.10 shows the comparison between the two models
and the experimental measurements for the robot with the new tail. Consistent with the results in
Fig. 4.9, one can see that the model using the multi-segment approximation matches the experimental measurement closely for all actuation frequencies, while the linear beam model performs
well only for low frequency actuation.
72

12

Experiment
Multi−segment Model
Beam Model

Speed (cm/s)

10
8
6
4
2
0
0

0.5

1

1.5

2

2.5

3

Tail−beat frequency (Hz)
Figure 4.10: Comparison between model predictions and experimental measurement of the speed versus
tail-beat frequency, with the new tail. The amplitude is fixed at 13.6◦ .

To further assess the proposed multi-segment model, we conduct an experiment involving transients in forward motion, where the robot started from at rest. Fig. 4.11(a) and (b) compare the
model predictions of the X /Y -coordinate time-trajectories of the robot with the experimental measurements, and Fig. 4.11(c) compares the predicted and actual robot paths in the XY -plane. It can
be seen that the proposed model is able to capture well both the transient and steady-state behaviors
of the robot.
We have conducted additional experiments to compare the time-dependent shape of the flexible
tail with those predicted by the models. As can be seen in Fig. 4.12, the free swimming robot shown
in Fig. 4.4 is fixed by a bracket and set to oscillate the tail. A Casio Exilim EX-FH25 high-speed
camera with a frame rate of 120 frames/s is used to record the tail’s motion from above.
Figs. 4.13 and 4.14 compare the measured time-dependent tail shape and those predicted by
73

the two models for 0.4 Hz and 0.9 Hz, respectively. To save space, we show every ninth frame
for half a period of the tail oscillation. It can be seen that when the tail beats at the relative low
frequency, 0.4 Hz, both models produce very close approximations of the tail shape. However,
at 0.9 Hz, when the beam deformation is bigger, the multi-segment modeling approach produces
much more precise predictions of the beam shape than the linear beam theory-based approach.

4.6 Conclusions
In this chapter, we have developed a model for robotic fish with a base-actuated flexible tail.
The tail is modeled as multiple rigid segments connected by springs and dampers, and Lighthill’s
elongated-body theory is used to evaluate the tail-generated hydrodynamic forces. For comparison,
we have also constructed a model using linear beam theory. We compare predictions of both models to experimental results with a robotic fish, in terms of steady-state cruising speeds and dynamic
tail shapes under different actuation frequencies. From these results, we conclude that when the
tail is excited under a relatively low frequency, and consequently experiences small deformation,
both models produce similar predictions that are close to experimental measurements. However,
when the actuation frequency increases, the two models differ, and the model using multi-segment
approximation is able to predict much better the robot speed and the tail shape. Additional experimental results also indicate that the proposed model is capable of capturing the transient dynamics
of the robot.
In summary, the work presented in this chapter provides a computationally efficient and accurate model for capturing large tail deformation and the resulting hydrodynamic force for a tailactuated robotic fish. The model will facilitate effective tail design optimization and controller
development for such robots. While one could also use nonlinear beam models and computational

74

fluid dynamics (CFD) to achieve faithful modeling of the tail, the latter approach would be much
more computationally expensive, and would be difficult to integrate with the robot dynamics for
controller design.
This work can be extended in several directions. First, although this chapter has focused on
the case of a rectangular tail for ease of presentation, the mathematical derivation itself is general
and the approach can be extended to tails of other shapes (e.g., trapezoidal) following a similar
treatment, where the properties of each segment (mass, inertia) and each joint (spring and damper
constants) will depend on the local shape. Second, we are particularly interested in using the proposed model to understand the effect of tail shape and stiffness properties on the robot’s locomotion
performance, and exploit such understanding for design optimization. Finally, we will utilize the
proposed model to design controllers for flexible tail-actuated robotic fish.

75

Experiment
Model

X (m)

1
0.8
0.6
0.4
0.2
0
0

5

10

15

20

Time (s)
(a)

Experiment
Model
Y (m)

0.2

0.1

0
0

5

10

15

20

Time (s)
(b)
0.2

Experiment
Model

Y (m)

0.15
0.1
0.05
0
0

0.2

0.4

0.6

X (m)
(c)
Figure 4.11: Comparison between multi-segment model predictions and experimental measurement for
forward swimming (including transients): (a) time trajectory of X -coordinate of the robot; (b) time trajectory of the Y -coordinate of the robot; (c) path of the robot in the XY -plane. For both the experiment and
simulations, the amplitude and frequency of the tail beat are fixed at 13.6◦ and 0.9 Hz, respectively.

76

Casio HighͲspeed Camera

Flexible Tail

Bracket
Robotic Fish

Figure 4.12: Experimental setup to capture the dynamic shape of a flexible tail actuated at the base.

77

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4.13: Comparison between experimental measurement of the time-dependent tail shape with model
predictions. The tail beats at 0.4 Hz with 0◦ bias and 14◦ amplitude. The black solid line, blue dashed line
with circles and the red dash-dot line imply the experimental measurement, predictions from multi-segment
model and Euler-Bernoulli beam model, respectively.

78

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4.14: Comparison between experimental measurement of the time-dependent tail shape with model
predictions. The tail beats at 0.9 Hz with 0◦ bias and 13.6◦ amplitude. The black solid line, blue dashed line
with circles and the red dash-dot line imply the experimental measurement, predictions from multi-segment
model and Euler-Bernoulli beam model, respectively.

79

Chapter 5
Control-oriented Averaging of Tail-actuated
Robotic Fish Dynamics
5.1 Introduction
Much of the work in the area of robotic fish has focused on the design ( [5, 6, 8, 10, 131]) and
dynamic modeling [11, 14, 15, 45, 92, 95, 96, 132]. Due to the rhythmic nature of body and fin
movements of robotic fish, averaging represents a promising approach in obtaining models that
simplify motion analysis and control design [133]. While classical averaging [134], which simply
integrates the original vector field over one period, is easy to implement, it cannot be applied
directly to robotic fish dynamics, because the latter typically involves large-amplitude oscillatory
inputs and fails to satisfy the slow-dynamics assumption required by classical averaging theory.
Geometric averaging [135, 136], on the other hand, can handle systems with highly oscillatory
inputs, but as we show in this chapter, for an original dynamic model that is reasonably accurate
(and complex), the resulting average model is overly complicated and has limited use for system
analysis and control design.
In this chapter, a novel averaging approach is proposed for robotic fish dynamics, with a focus
on planar motions of a rigid tail-actuated robotic fish. Note that among various actuation mechanisms, actuation with an oscillating caudal fin is particularly attractive (e.g., [6,11,12,44,45]). This
is because the latter is easy to realize, enables both forward swimming and turning maneuvers, and
80

leaves the majority of the robot body free of moving parts. We adopt a model where Lighthill’
large-amplitude elongated-body theory [25] is used to evaluate the tail-generated hydrodynamic
force, since such a model strikes a sound balance between fidelity and simplicity [15].
The proposed averaging approach consists of properly scaling the force and moment terms in
the dynamic model and then conducting classical averaging. To identify the scaling functions, a
large number of tail-beat patterns are used in simulation of the original dynamics; for each tail-beat
pattern (bias, amplitude, and frequency), the scaling coefficients for the force and moment terms
are determined such that the corresponding average model will result in turning radius and turning
period at the steady state that best match those of the original dynamics. It is found that, for a
given set of parameters for the robotic fish model, its force scaling function is a constant (thus
independent of the tail-beat pattern) while its moment scaling function is linear in the tail-beat bias
but independent of tail-beat frequency and amplitude.
The proposed averaging model is evaluated extensively with simulation and experimental results, including the comparison of predictions of turning period/radius and cruising speed at the
steady state, by the original model and the average model, respectively. Also included is the validation of the average models ability to predict transient behavior of the robot. Most of the simulation
and experimental scenarios used in the evaluation have not been used in identifying the scaling
coefficients, which provides strong support for the effectiveness of the average model.
To gain analytical insight into the obtained average model, we analyze the existence and stability of its equilibrium points. In particular, with several reasonable assumptions, we establish the
existence and uniqueness of equilibrium points for the (approximate) average model under a given
tail-beat pattern. Furthermore, through linearization analysis, we show that any equilibrium of the
average model is locally asymptotically stable.
Finally, to illustrate the utility of the proposed average model in motion analysis, we present
81

a semi-analytical framework for modeling the steady turning of robotic fish. Based upon the proposed average model, the force and moment balance equations during steady turning are set up
explicitly in terms of three unknowns, the turning radius, angular speed (equivalent to turning period information) and the angle of attack. These algebraic equations can then be solved to get the
turning-related parameters for any given tail-beat pattern. Such a framework is different from existing work (e.g. [14,15,137,138]), where dynamical equations are used for simulation and turning
parameters can only be found by observing the simulated trajectories. We note that the average
model also holds potential for controller design; for example, Chen et al. has used a preliminary
form of this model in backstepping-based control design for a robotic fish [139].
The remainder of the chapter is organized as follows. In Section 5.2, we summarize the dynamic model for a tail-actuated carangiform robotic fish. In Section 5.3 we review the first-order
geometric averaging method and the classical averaging method, and discuss why they are inadequate for averaging the robotic fish dynamics. The proposed averaging method is developed and
evaluated in Section 5.4. Analysis of the resulting average model is presented in Section 5.5. The
semi-analytical framework for turning is described in Section 5.6. A nonlinear controller design for
the robotic fish to achieve trajectory tracking using the average model is presented in Section 5.7.
Finally, concluding remarks are provided in Section 5.8.

82

5.2 Dynamic Model of Robotic Fish
As described in Section 3.2 of Chapter 3, the equation of planar motion for a rigid body in an
inviscid fluid is governed by the Kirchhoff’s equations

(mb − XV˙ )V˙Cx = (mb −YV˙ )VCy ωz + Fx ,
Cx
Cy

(5.1)

(mb −YV˙ )V˙Cy = −(mb − XV˙ )VCx ωz + Fy ,
Cy
Cx

(5.2)

(Jbz − Nω˙ z )ω˙ z = (YV˙ − XV˙ )VCxVCy + Mz ,
Cy
Cx

(5.3)

where Fx , Fy , Mz denote the external forces and moments about the body center, which are caused
by the tail motion and the interaction between the body itself and the surrounding fluid, and can be
evaluated by














Fx =< ( 0L f (ζ )dζ + FL ), xˆ > −FD cos β + FL sin β ,
Fy =< ( 0L f (ζ )dζ + FL ), yˆ > −FD sin β − FL cos β ,

(5.4)

Mz = 0L rCζ × f (ζ )dζ + rCL × FL + MD ,

where rCζ denotes the vector from the body center to the point ζ on the tail, i.e., rCζ = −(c +

ζ cos α )xˆ − (ζ sin α )y,
ˆ and c is the distance from the body center to the joint. “×” represents the
cross product of vectors.
Finally, the dynamic model for the robotic fish can be summarized as



V˙Cx = f1 (VCx ,VCy , ωz ) + f4 (t),




V˙Cy = f2 (VCx ,VCy , ωz ) + f5 (t),





 ω˙ = f (ω ) + f (t),
z
3 z
6

83

(5.5)

where



f1 (VCx ,VCy , ωz ) =
















f2 (VCx ,VCy , ωz ) =








m2
c1
m1 VCy ωz − m1 VCx

2 +V 2
VCx
Cy
V

c

2 +V 2 arctan( Cy ),
+ m2 VCy VCx
V
Cy
1

Cx

m

c

2 +V 2
− m1 VCx ωz − m1 VCy VCx
Cy
2

c
− m2 VCx
2

2

V

2 +V 2 arctan( Cy ),
VCx
V
Cy
Cx

f3 (ωz ) = (m1 − m2 )VCxVCy − c4 ωz2 sgn(ωz ),



























c

f4 (t) = − m3 α¨ sin α ,
1
c

f5 (t) = + m3 α¨ cos α ,
2
f6 (t) = −c5 α¨ cos α − c6 α¨ ,

with c1 = 1/2ρ SCD , c2 = 1/2ρ SCL, c3 = 1/2mL2 , c4 = 1/(J3 )KD , c5 = 1/(2J3)L2 mc and c6 =
1/(3J3)L3 m, m1 = mb − XV˙ , m2 = mb −YV˙ , J3 = Jbz − Nω˙ z , where L is the tail length and α¨
Cx

Cy

denotes the second-derivative of the tail angle α with respect to time.

5.3 Averaging with Existing Methods
Since the tail-actuation for robotic fish is typically periodic, averaging is a useful tool for gaining
insight into the effect of the input parameters (such as tail beat bias, period, and frequency) on
the dynamics [11, 133], and for designing the controller [11, 139]. In this section we review some
fundamental results on averaging theory [134, 136], apply them to the robotic fish model, and
discuss their limitations.

84

5.3.1 First-order Averaging for Systems with Oscillatory Inputs
To address the averaging of the dynamics for a robotic fish, one can consider first-order averaging
for systems with oscillatory control input [136]. In particular, consider the following dynamics

x˙ = f (x) + (1/ε )g(x,t/ε ),

(5.6)

where ε > 0 is a small number, and (1/ε )g(x,t/ε ) is a T −periodic, zero mean function in t/ε to
indicate a high-amplitude and large-frequency oscillatory input. With the transformation of time,
t/ε → τ , one can obtain
dx
= ε f (x) + g(x, τ ).
dτ

(5.7)

The term ε f (x) can be treated as a perturbation to the primary vector field g(x, τ ). Define the
following

g

F(y, τ ) =

(Φ0,τ )∗ f (y),

(5.8)

Fav (y) =

1 T
F(y, τ )d τ ,
T 0

(5.9)

g

g

g

where Φ0,τ is the flow of the vector field g, and (Φ0,τ )∗ f is the pull-back of f by Φ0,τ . According
to the Variation of Constants formula, the solution x(τ ) is obtained by

g

(5.10)

dy/d τ = ε F(y, τ ),

(5.11)

x(τ ) = Φ0,τ (y(τ )),

where y(τ ) is the solution to the system

85

with y(0) = x0 . An approximation to the solution y can be obtained from the average equation

dz/d τ = ε Fav (z).

(5.12)

It can be shown that [136] z(τ ) − y(τ ) = O(ε ) on the time scale 1. Moreover, if the origin is
an asymptotically stable equilibrium point for the average dynamics, then z(τ ) − y(τ ) = O(ε ) as

ε → 0 for all τ ∈ R+ and the differential equation (5.11) possesses a unique periodic orbit that is
locally asymptotically stable and belongs to an O(ε )-neighborhood of the origin.
The pull-back used in Eq. (5.8) can be calculated as

g

(Φ0,t )∗ f
∞

= f+

∑

(5.13)
t

k=1 0

···

sk−1
0

(adg(s ) · · · adg(s ) f )dsk · · · ds1 ,
k

1

where adX Y = [X ,Y ] is the Lie bracket of X and Y , or the Lie derivative of Y with respect to X .
The convergence of the infinite sum in (5.13) is generally not guaranteed unless the series contains
only a finite number of non-vanishing terms.

5.3.2 First-Order Averaging for Robotic Fish Model
In [133], Morgansen et al. applied the first and higher-order averaging method to a simplified
carangiform robotic fish model and demonstrated its effectiveness with trajectory stabilization controller design. We now apply the first-order averaging analysis to the robotic fish model presented
in Section 5.2.

86

We consider the following periodic pattern for the tail deflection angle:

α (t) = α0 + αA sin(ωα t),

(5.14)

where α0 , αA and ωα denote the bias, amplitude, and frequency of the tail beat, respectively.
Define ε = 1/ωα and τ = t/ε , which allows us to express

α (t) = α0 + αA sin(τ ),

(5.15)

1
α¨ (t) = − ωα αA sin(τ ).
ε

(5.16)

Define the states x1 = u, x2 = v, x3 = ω . The system (5.5) can be rewritten in the form of (5.6)
1 t
x(t)
˙ = f (x) + g( ),
ε ε

(5.17)

where g(·) is zero mean over a period of T . To facilitate the computation of the average, we
first use the second-order Taylor series expansion to approximate sin α ≈ α − 61 α 3 and cos α ≈
1 − 21 α 2 . The terms like arctan(VCy /VCx ) in Eq. (5.5) impede the averaging process due to its
highly nonlinear nature. To address this challenge, we further simplify (5.5) with the assumptions
that VCx ≫ VCy and VCx > 0. In other words, the forward velocity is positive, and much larger than
the sideway velocity. These assumptions are reasonable, and they imply that the terms

2 +V 2
VCx
Cy

and arctan(VCy /VCx ) can be approximated by VCx and VCy /VCx , respectively. We also neglect the
coupling term (m1 − m2 )VCxVCy , when evaluating the total moment in Eq. (10), based upon the
observation that it is usually less than 2% of the moment Mz that is caused by external forces. With

87

the aforementioned simplifications, the system (5.5) can be re-written as:
 

 

x˙1 
 f1 (x) g1 (t)
 

 

1 t
 

 

x˙2  = f (x) + g( ) =  f2 (x) + g2 (t)
 

 

ε ε
 

 

x˙3
f3 (x)
g3 (t)

(5.18)

where

k
k
m2
x2 x3 − 1 x21 + 2 x22 + k4 ,
m1
m1
m1
(k + k2 )
m
f2 (x) = − 1 x1 x3 − 1
x1 x2 + k5 ,
m2
m2
f1 (x) =

f3 (x) = −k3 x23 sgn(x3 ) − k6 ,
g1 (t) =

k7 (k8 sin ωα t − k9 cos 2ωα t − k10 sin3 ωα t
+k11 cos 2ωα t − k12 cos 4ωα t) ,

g2 (t) = −k13 (k14 sin ωα t + k15 cos 2ωα t − k16 sin3 ωα t),
g3 (t) =

k17 (k18 sin ωα t + k19 cos 2ωα t − k20 sin3 ωα t)
+k21 sin ωα t .

K
m L2 α 2 ω 2 (3α − 3 α 2 α − 3 α 3 ), k = m L2 α α 2 ω 2 ,
with k1 = 21 ρ SCD , k2 = 12 ρ SCL , k3 = JD , k4 = 12m
A 2 0 A 8 A
0 A α
5 4m2
A α
3
1
2

L mcα α 2 ω 2 , k = m L2 α ω 2 , k = 6α − α 3 , k = 3α − 3 α 2 α , k = 3α α 2 , k =
k6 = 4J
0 A α 7
0
0 A 11
A α 8
A 2 0 A 10
12m1
0 9
3
1 3
1 3
m 2
L2
2
2
2
2
2 αA , k12 = 8 αA , k13 = 4m2 L αA ωα , k14 = 2 − α0 , k15 = α0 αA , k16 = αA , k17 = 4J3 mcαA ωα ,
3
k18 = 2 − α02 , k19 = α0 αA , k20 = αA2 , and k21 = L3Jm αA ωα2 ; note that g1 (t), g2 (t) and g3 (t) are
3

T −periodic, zero mean functions.
For the approximate model (5.18), the pull-back operation can be obtained directly using its
original definition rather than through (5.13), since the g field does not involve the state vari-

88

g

g

ables. It can be easily shown that (Φ0,τ )∗ f = f (Φ0,τ ). With the definition of pull-back function
g

(Φ0,τ )∗ f , we get
k
m2
(y2 + h2 )(y3 + h3 ) − 1 (y1 + h1 )2
m1
m1
k2
+ (y2 + h2 )2 + k4 ,
m1
(k + k2 )
y˙2 = − 1
(y1 + h1 )(y2 + h2 )
m2
m
+k5 − 1 (y1 + h1 )(y3 + h3 ),
m2

y˙1 =

y˙3 = −k3 (y3 + h3 )2 sgn(y3 + h3 ) − k6 ,

(5.19)

(5.20)
(5.21)

where h1 (t) = 0t g1 (τ )d τ , h2 (t) = 0t g2 (τ )d τ , h3 (t) = 0t g3 (τ )d τ , respectively. The term (y3 +
h3 )2 sgn(y3 + h3 ) in (5.21) will make the averaging overly complicated.
1 k k , and k = k k −
Define k22 = k17 k18 − 34 k17 k20 + k21 , k23 = 12 k17 k19 , k24 = 12
17 20
17 18
25
3 k k + k + 1 k k . We have
21 12 17 20
4 17 20

− k22 cos ωα t + k23 sin 2ωα t − k24 cos 3ωα t + k25
ωα
1
≈
− k22 cos ωα t + k22 ,
ωα

h3 (t) =

(5.22)

since the size of k22 is typically two orders of magnitude bigger than those of k23 , k24 , |k22 − k25 |.

89

With such an approximation, the average dynamics (5.12) is given below,

z˙1 = z2 z3 + h¯ 2 z3 + h¯ 3 z2 + h¯ 23 − k1 z21
−2k1 h¯ 1 z1 − k1 h¯ 11 + k2 z22 + 2k2 h¯ 2 z2 ,

(5.23)

z˙2 = −z1 z3 − h¯ 1 z3 + h¯ 3 z1 + h¯ 13 − (k1 + k2 )z1 z2 + k5
−(k1 + k2 )h¯ 1 z2 − (k1 + k2 )h¯ 2 z1 − (k1 + k2 )h¯ 12 ,
z˙3 = −k3 {π − 2 arccos(1 +

(5.24)

ωα
)/π }
k22 z3

k22 2k22 1
ωα 2
+
1 − (1 +
z )
ωα
π ωα
k22 3
ωα
arccos(1 +
z )
k22 3
k2
1
ωα
− 2 22
z ))
sin(2 arccos(1 +
ωα π 4ωα
k22 3
ωα 2
ωα
3
1 − (1 +
arccos(1 +
z3 ) +
z )
k22
2ωα
k22 3

−2k3 z3
−

1
ωα

−k3
−
+

2
ωα

2
3k22

(5.25)

2ωα2

with h¯ 1 = T1 0T h1 dt, h¯ 2 = T1 0T h2 dt, h¯ 3 = T1 0T h3 dt, h¯ 11 = T1 0T h21 dt, h¯ 12 = T1 0T h1 h2 dt,
g
h¯ 13 = T1 0T h1 h3 dt, and h¯ 23 = T1 0T h2 h3 dt. Then the first-order average model is x(t) = Φt z(t) ,

according to Variation of constants formula [136].
With the tail-beat bias, frequency, and amplitude set to be 20◦ , 1 Hz, and 15◦ , respectively,
we compare the turning radii and periods predicted by the model (5.18) and the corresponding
average model. The parameters used in the simulation are based on those identified for the robotic
fish as discussed later in the next section. The results are (0.26776 m, 59.7734 s) and (0.270276 m,
60.0873 s), respectively, which indicates that the first-order averaging is able to capture the original
dynamics reasonably well. However, the resulting average model (5.23)-(5.25) is overly complex

90

to be useful for analysis or controller design, even though it is derived based on simplifications of
the original dynamics.

5.3.3 Classical Averaging
Considering its simplicity, we have also examined the possible use of the classical averaging
method. The classical averaging can be applied to a slow system of the form

x˙ = ε f (t, x),

x(0) = x0 ,

(5.26)

where ε is a small positive number and f (t, x) is T −periodic in t. The associated averaging system
is
y˙ = ε fav (y),

y(0) = x0 ,

(5.27)

1 T
f (τ , y)d τ .
T 0

(5.28)

where
fav (y) =

Theorem 1 [134] There exists a positive ε0 , such that for all 0 < ε ≤ ε0 ,
1)

x(t) − y(t) = O(ε ) as ε → 0 on the time scale 1/ε and

2)

if the origin is asymptotically stable for fav , then x(t) − y(t) = O(ε ) as ε → 0 for all t ∈ R+

and the system (5.26) possesses a unique periodic orbit which is locally asymptotically stable and
belongs to an O(ε ) neighborhood of the origin.
This averaging method is easy to implement; however, it cannot be directly applied to the
robotic fish dynamics as the latter is not a slow system that satisfies Eq. (5.26). To illustrate this
point, we have conducted simulation comparison of the original dynamic model and the averaged
model with classical averaging, for the tail-beat patterns of α0 = 20◦ , ωα = 2π rad/s, and αA = 15◦ .
91

As shown in Fig. 5.1, the averaged model predicts a motion orbit that is significantly different from
that predicted by the original dynamics.

0

Original Model
Averaged Model

Y (m)

−0.1
−0.2
−0.3
−0.4
−0.5
−0.2

−0.1

0

0.1

0.2

0.3

0.4

X (m)
Figure 5.1: The trajectories of robotic fish predicted by the original model and by the classical average
model, respectively.

5.4 Computationally Eff cient Averaging with Scaled Forcing
5.4.1 Proposed Averaging Approach
As demonstrated in the previous section, although the first-order averaging method for systems
with oscillatory inputs can approximate well the original dynamics for robotic fish, the resulting
model is overly complex for controller design and other applications. On the other hand, classical
averaging is easy to implement with simple expressions, but it cannot be applied to the robotic fish
model directly since the latter does not satisfy the slow dynamics requirement. Motivated by these
92

observations, we propose scaling the original forcing terms by some functions that are potentially
dependent on the tail-beat parameters, α0 , αA and ωα , and then conduct the classical averaging
over the modified dynamics. Specifically, we modify the original system (5.5) to



x˙1 = f1 (x) + K f orce · f4 (t),




x˙2 = f2 (x) + K f orce · f5 (t),





 x˙ = f (x) + K
moment · f6 (t),
3
3

(5.29)

where K f orce , Kmoment are scaling functions dependent on tail-beat parameters. Note that the physical meanings of f4 (t) and f5 (t) are the projections of instantaneous hydrodynamic forces along
the longitudinal and transverse directions in the body-fixed coordinates, so the scaling functions
for them are the same. Using the second-order Taylor series expansion to approximate the sin α
and cos α terms in (5.29), we obtain the averaged system as


m L2 ω 2 α 2 (3 − 3 α 2 − 3 α 2 ) ,

x˙1 = f1 (x) + K f orce 12m

α A

2 0
8 A

1

m L2 ω 2 α 2 α ,
x˙2 = f2 (x) + K f orce 4m
α A 0

2




m 2 2 2
 x˙ = f (x) + K
moment (− 4J L cωα αA α0 ) .
3
3

(5.30)

3

We have conducted extensive simulations with different sets of α0 , αA , and ωα , to identify the
corresponding values of K f orce and Kmoment that would result in good match between the average
model (5.30) and the original dynamics (5.5).
In the simulation, we use the parameters identified for a free-swimming robotic fish prototype
as shown in Fig. 3.3, which has been reported in Chapter 3. The tail deflection angle α is controlled
by a servomotor through the tail shaft. The hydrodynamic coefficients for the body are CD = 0.386,
CL = 4.50, and KD = 7.82 × 10−4 kg·m2 .

93

For a tail-beat pattern (5.14) with α0 = 0, the trajectory of the robotic fish will converge to
a circular orbit. For example, Figs. 5.2 and 5.3 show the simulated trajectory and heading of the
robotic fish, respectively, based on the original dynamics, when the tail beats at ωA = 2π rad/s with
bias α0 = 20◦ and amplitude αA = 15◦ . The trajectory is close to a circle with oscillatory ripples,
as seen in the close-up in Fig. 5.2b. The average heading changes continuously with similar smallamplitude oscillations as seen in Fig. 5.3b. From the simulation, we can extract the turning radius
and turning period (time taken to complete one turn) of the robot at the steady state.
With the average model (5.30), the trajectory of the robot converges to a circle free of ripples.
Analysis of the stability properties for the average model will be presented in Section 5.5. For
each set of (α0 , αA , ωα ), we have conducted blanket-search of parameters K f orce and Kmoment ,
such that the resulting average model would produce the best match in turning radius and period
with the original dynamics.
In the simulation, we have used five tail-beat frequencies: 1, 2, 3, 4, 5 Hz; five bias angles:
10, 20, 30, 40, 50 degrees; and five amplitudes: 10, 15, 20, 25, 30 degrees, and find the optimal
coefficients for all 5 × 5 × 5 = 125 combinations of these parameters. Fig. 5.4 shows the optimal
coefficients for different combinations of tail-beat bias and amplitude when the tail beats at 1 Hz.
The results for other frequencies are nearly identical to those in Fig. 5.4, and are thus not presented
in the interest of brevity. From Fig. 5.4(a), we can see that the optimal force coefficient is nearly
constant for different biases and amplitudes. On the other hand, as shown in Fig. 5.4(b), the
optimal moment coefficient shows linear dependence on the bias angle and barely changes with
the amplitude.
Fig. 5.5 shows how the scaling coefficients vary with the bias and the frequency for a given
amplitude (15◦ ); nearly identical results hold for the other amplitudes simulated. Again, K f orce
is nearly constant for different combinations of tail-beat bias and frequency, and Kmoment shows
94

0

Y (m)

−0.1
−0.2
−0.3
−0.4
−0.5
−0.2

0

0.2

0.4

X (m)
(a)

Y (m)

−0.21

−0.22

−0.23

−0.24
0.35

0.36

0.37

0.38

X (m)
(b)
Figure 5.2: Simulated trajectory of the center of the robotic fish when the tail beats at 1 Hz with
bias 20◦ and amplitude 15◦ : (a) The global path; (b) close-up view of the segment within the red
rectangular region in (a).
linear dependence on the bias and independence of the frequency. Fig. 5.6 shows the obtained
coefficients versus the frequency and the amplitude, for a fixed bias of 20◦ . In consistency with
Figs. 5.4 and 5.5, we see that both coefficients are independent of the frequency and the amplitude.
In addition, in all cases, the values of K f orce are very close, as are the slopes of Kmoment with
95

Heading angle (rad)

2
0
−2
−4
−6
−8
0

20

40

60

80

Time (s)
(a)

Heading angle (rad)

−3
−3.2
−3.4
−3.6
−3.8
−4
35

40

45

Time (s)
(b)
Figure 5.3: Simulated heading angle of the robotic fish when the tail beats at 1 Hz with bias 20◦
and amplitude 15◦ : (a) The heading angle trajectory; (b) close-up view of the segment within the
red rectangular region in (a).
respect to the bias.
Based upon above observations, we get the force scaling coefficient by taking the average of

96

Force coefficient

0.88

1

0.86
0.5
0.84
0
50

0.82
30
20

Bias (°)

0

10

0.8

Amplitude (°)

Moment coefficient

(a)

0.4
0.35
0.4
0.3
0.25

0.2

0.2
0
50

0.15
30
20

Bias (°)

0

10

0.1

Amplitude (°)

(b)
Figure 5.4: K f orce and Kmoment versus the tail-beat bias and amplitude. In all cases, the tail beats
at 1 Hz.
all corresponding coefficients obtained in 125 individual simulations, which gives

K f orce = 0.8485.

For each tail beating bias, we take the average of the optimal Kmoment for 25 different combinations
97

Force coefficient

0.84

1

0.83
0.5
0.82
0
50

0.81
5
0.8

Bias (°)

0 0

Frequency (Hz)
(a)

Moment coefficient

0.4
0.35
0.4

0.3
0.25

0.2

0.2
0
50

0.15
5

Bias (°)

0 0

0.1

Frequency (Hz)
(b)

Figure 5.5: K f orce and Kmoment versus the tail-beat bias and frequency. In all cases, the tail beats
with amplitude 15◦ .
of amplitude and frequency, as summarized in Table 5.1. Then the least-square-error fitting is used
to identify the relationship between Kmoment and the bias α0 , resulting in

Kmoment = −0.0074 + 0.4831α0,

98

Force coefficient

0.87

1

0.86
0.5

0.85
0.84

0

0.83

5

30
20

Freq (Hz)

0

10

0.82

Amplitude (°)

(a)

Moment coefficient

0.17
0.4
0.165
0.2
0.16
0
5

30
20

Freq (Hz)

0

10

0.155

Amplitude (°)

(b)
Figure 5.6: K f orce and Kmoment versus the tail-beat amplitude and frequency. In all cases, the tail
beats with bias 20◦ .
where α0 is expressed in radians.

99

Table 5.1: Average of moment coefficients for different bias angles.
Bias (◦ )
10
20
30
40
50

Average of moment coefficient
0.0796
0.1596
0.2424
0.33
0.416

5.4.2 Numerical Evaluation of the Average Model
In order to evaluate the proposed averaging approach, we simulate the average model (5.30) with
the obtained scaling functions. Fig. 5.7 shows the comparison of the turning radii and turning
periods predicated by the original model and the average model, respectively, under different tailbeat biases, where the tail-beat frequency and amplitude are fixed at 1 Hz and 20◦ , respectively.
We have deliberatively chosen cases of 25◦ and 35◦ biases at this stage (and the latter part in this
section), which are not in the aforementioned 125 motion combinations, to provide independent
validation for the proposed averaging approach. It can be seen that the average model is able to
capture the original dynamics in the steady turning motion with high fidelity.
We then validate the effectiveness of the proposed averaging approach by comparing the cruising speeds predicted by the original and average models, when the tail beat bias α0 = 0. Note that

α0 = 0 has not been used in obtaining the scaling coefficients K f orce and Kmoment , and thus this
new set of simulation provides a meaningful test for the averaging approach. For this simulation,
we randomly pick some tail-beat frequencies and amplitudes, and the simulated velocities using
the original and average models are shown in Table 5.2. It can be seen that the prediction error by
the average model is less than 2% for all cases.
To further evaluate the generality of the proposed averaging method, we compare the original and average models for new sets of parameters. In particular, we double the weight of the
robot (and the inertia of robot’s body, correspondingly) from mb = 0.311 kg to mb = 0.622 kg,
100

35

Dynamic model
Average model

70

25

Turning period (s)

Turning radius (cm)

30

Dynamic model
Average model

80

20

15

10

60
50
40
30
20

5

0

10

20

25

30

35

0

40

Bias (°)

20

25

30

35

40

Bias (°)

Figure 5.7: Comparison of steady-state turning radii and periods predicted by the original and average
models. In all cases, the tail-beat frequency is 1 Hz and amplitude is 20◦ .

Table 5.2: Surge velocities predicted by the original and average models, and the relative error by
the average model.
Frequency (Hz)
3.1181
1.0955
1.8308
1.1385
1.2914
3.4704
3.0845
1.9513
3.8507
1.1033

Amplitude (◦ )
18.7749
17.6312
25.3103
25.9040
13.7375
19.7953
18.9117
22.9263
24.1873
25.0937

Original Averaged Evel (%)
0.112114 0.114017
1.70
0.036991 0.037648
1.78
0.088677 0.089746
1.21
0.056431 0.057084
1.16
0.033976 0.034661
2.0
0.131556 0.133695
1.62
0.111713 0.113599
1.69
0.085645 0.086837
1.39
0.178275 0.180581
1.29
0.052984 0.053632
1.22

while keeping the other parameters unchanged (note that the added mass/inertia are determined by
robot’s dimensions and independent of robot’s mass/inertia). Then following the same procedure
as described in the early part of this section, one could identify the best scaling functions for the

101

system, by matching the turning parameters (radius and period) from the original dynamic model
and from the average model. We fix the tail-beat amplitude and frequency to be 20◦ and 1 Hz,
respectively, and vary the tail-beat bias. We note that, instead of blanket-searching of parameters
by running extensive simulations, the observations of behaviors for K f orce and Kmoment can be exploited to find these parameters in a significantly more efficient manner. In particular, we find the
best scaling coefficients for the system when the bias is 10◦ , 20◦ and 30◦ , respectively, as shown
in Table 5.3.
Table 5.3: Best scaling coefficients for different bias angle. (mb = 0.622 kg in the simulation)
Bias (◦ )
10
20
30

Force coefficients Moment coefficient
0.8780
0.2500
0.8980
0.4650
0.8980
0.610

The force coefficient is then chosen to be K ′f orce = 0.891 by averaging the coefficients. The
′
moment coefficients can be obtained by a linear interpolation Kmoment
= 0.0817 + 1.0313α0. Then

the simulations are conducted with the updated force and moment coefficients and the results are
shown in Figures 5.8 and 5.9. We can see that the predictions by average model match well those
predicted by the original dynamic model.
To investigate the performance of the average model when the tail-beat parameters vary over
time, in particular, under a feedback control context, we design a simple proportional controller for
the original system. In the simulations, the proportional controller depends on the present error,
which is defined as the difference between the orientation of the target relative to the robot center
and the heading of the robot. The tail is programmed to beat at 1 Hz, with an amplitude of 15◦ ,
while the bias is updated by the proportional controller every five seconds. Fig. 5.10a shows the
simulation trajectory of the robot, where it starts at the origin from rest and moves to the target
location (2, 2). The corresponding bias obtained from the feedback is shown in Fig. 5.10c. The
102

40

60

Dynamic model
Average model

35

Dynamic model
Average model

Turning period (s)

Turning radius (cm)

50
30
25
20
15

40

30

20

10
10
5
0

10

15

20

25

0

30

Bias (°)

10

15

20

25

30

Bias (°)

Figure 5.8: Simulation results from the original and average models (versus different tail-beat bias), when
mb = 0.622 kg. In all cases, the tail-beat frequency is 1Hz and amplitude is 20◦ .

time trajectory of this bias is then applied to the averaged system, which predicts the robot trajectory as shown in Fig. 5.10b. With the same bias input, the terminal point with the original model
is (1.9176, 2.2430) while the robot trajectory with the average model stops at (2.0849, 2.2803)
after 200 s of simulation time. By comparing Fig. 5.10a and 5.10b, we can see that there is some
modest error, which can be explained by that the average model does not account for the transients
when the tail-beat parameters change. However, overall the average model is able to predict the
trajectory of the robot.

5.4.3 Experimental Evaluation of the Average Model
For model validation purpose, we have conducted experiments with the robot in forward swimming. We note that this forward swimming case (α0 = 0) was not used in obtaining the scaling
103

25

Dynamic model
Average model
20

50

Turning period (s)

Turning radius (cm)

Dynamic model
Average model

60

15

10

40

30

20

5
10

0

10

20

30

40

0

50

Amplitude (°)

10

20

30

40

50

Amplitude (°)

Figure 5.9: Simulation results from the original and average models (versus different tail-beat amplitude),
when mb = 0.622 kg. In all cases, the tail-beat frequency is 1Hz and bias is 20◦ .

functions, and thus this experimental result provides a significant validation for the averaging approach. Fig. 5.11 shows the experimental and simulation results by comparing the time that the
robotic fish travels 50 cm long after it reaches steady speed, and we can see that the average model
is able to not only capture the original dynamics but also predict the robot’s motion with high
fidelity.
Previous discussions have exclusively focused on the steady motion behavior of a robotic fish.
It is of interest to investigate the performance of the average model in predicting the transient
dynamics. For this purpose, we have conducted experiments to record the transient trajectories
of the robot when it is started from at rest. A yellow LED is taped to the top of the robot as
a marker and its motion has been video-recorded with the ambient light turned off. The robot’s
time trajectory, as shown in Fig. 5.12, is then extracted by applying background subtraction to the
104

2

1.5

1.5

Y (m)

Y (m)

2

1
0.5
0

1
0.5

0

0.5

1

1.5

2

0

2.5

0

0.5

1

1.5

X (m)

X (m)

(a)

(b)

2

2.5

Control input (rad)

1
0.5
0
−0.5
−1
0

50

100

150

200

Time (s)

(c)
Figure 5.10: Trajectory of the center of the robot when the tail-beat bias is updated by a proportional feedback controller law (a) The simulated trajectory using the original model; (b) The
simulated trajectory using the average model; (c) The control input to both the original and average
models.

images. Fig. 5.12(a) and (b) compare the model predictions of time-trajectories of the robot with
the experimental measurements in X and Y coordinates, respectively. Fig. 5.12(c) provides the
comparison of the trajectories in the XY -plane. It can be seen that the developed average model is
able to capture not only the steady motion but also the transients.

105

25

Experiment
Original model
Average model

Time (s)

20

15

10

5

0

1

1.5

2

2.5

3

Frequency (Hz)

Figure 5.11: Comparison between prediction from original dynamic model, average model and experimental measurement for forward swimming.

5.5 Analysis of the Average Model
It is of interest to analyze the average model (5.30) in terms of its stability properties. In particular,
for a given tail-beat pattern, it is important to know whether the model admits an equilibrium, and
if so, whether the equilibrium is unique and stable.
The equilibrium points can be obtained by setting x˙ = 0. The highly nonlinear terms like
arctan(x2 /x1 ) make it impossible to get an analytical solution. However, if we make the same
assumptions as in Section 5.3.2, x1 ≫ x2 and x1 > 0, then we can obtain unique, analytical solutions. In particular, with these assumptions, we can approximate arctan(x2 /x1 ) ≈ (x2 /x1 ) and

106

0.2

0.8

Y (m)

X (m)

Experiment
Average Model

Experiment
Average Model

1

0.6
0.4

0.1

0.2
0
0

10

20

0
0

30

10

Time (s)
(a)

30

Experiment
Average Model

0.2

Y (m)

20

Time (s)
(b)

0.15
0.1
0.05
0
0

0.2

0.4

0.6

X (m)
(c)
Figure 5.12: Comparison between predictions from the average model and experimental measurement for
forward swimming (including transients): (a) time trajectory of X -coordinate of the robot; (b) time trajectory
of the Y -coordinate of the robot; (c) path of the robot in the XY -plane. For both the experiments and the
simulation, the amplitude and frequency of the tail beat are fixed at 15◦ and 0.9 Hz, respectively.

107

(x21 + x22 ) ≈ x21 . The equilibrium point is then solved as
K f orce mL2 ωα2 αA2 (3 − 32 α02 − 83 αA2 )
6ρ SCD

x1 =
x2 =

6mρ SCD L2 ωα2 αA2 α02

K f orce
2ρ S(CD +CL )

K f orce (3 − 23 α02 − 38 αA2 )

2m1
ρ S(CD +CL )

Kmoment mL2 cωα2 αA2 α0
4KD

+
x3 = −

Kmoment mL2 cωα2 αA2 α0
4KD

Next we analyze the stability property of the equilibrium for the model (5.30). Due to the
highly nonlinear nature of the model, analysis of global asymptotic stability is difficult. Instead,
we focus on the local stability in the neighborhood of the equilibrium. We rewrite Eq. (5.30) in a
compact form
x˙ = f(x, u),
where x = [x1 , x2 , x3 ]T and u = [α0 , αA , ωα ]T . We then linearize the system




 a11 a12 a13 

∂f 


A=
=  a21 a22 a23 

∂x 


a31 a32 a33

108

(5.31)

to get the Jacobian matrix A, the elements of which are

c
a11 = − 1
m1

x21

c
+ 2
x2 + x2 m1
1

x1 x2

x
arctan( 2 )
x1
x2 + x2

2

1

2

x2
c
c
− 1 x21 + x22 − 2 x21 + x22 2 2 2 ,
m1
m1
x1 + x2
x
x
c
c
m2
x
1 2
a12 =
x3 − 1
+ 2 x21 + x22 arctan( 2 )
m1
m1 x2 + x2 m1
x1
c
+ 2
m1

1
2
2
x2
x
c
arctan( 2 ) + 2
x1
m1
x21 + x22

m2
x ,
m1 2
c
m
a21 = − 1 x3 − 1
m2
m2

x x
x21 + x22 2 1 2 2 ,
x1 + x2

a13 =

c
− 2
m2
c
a22 = − 1
m2

x
x21 + x22 arctan( 2 )
x1

x1 x2

c
− 2
x2 + x2 m2
1

2

x21

x
c
arctan( 2 ) + 2
x1
m2
x2 + x2
1

2

x22

c
− 2
x2 + x2 m2
1

2

c
c
− 1 x21 + x22 − 2
m2
m2
m
a23 = − 1 x1 ,
m2

x x
x21 + x22 2 1 2 2 ,
x1 + x2

x1 x2

x
arctan( 2 )
x1
x2 + x2
1

2

x2
x21 + x22 2 1 2 ,
x1 + x2

a31 = 0,
a32 = 0,
a33 = −2c4 x3 sgn(x3 ).

By examining the Hurwitz property of A, we can know whether the averaged system is asymptotically stable at its equilibrium points. Note that a33 is negative except when x3 = 0. Since both a31

109

and a32 are 0, we only need to check the Hurwitz property of a two-two matrix




 a11 a12 
A′ = 
,
a21 a22
the characteristic equation of which is

λ 2 − (a11 + a22 )λ + (a11 a22 − a12 a21 ) = 0.

In order to show that A′ is a Hurwitz matrix, we need to have a11 + a22 < 0 and a11 a22 − a12 a21 >
0. It is easy to see that a11 + a22 = −(3c1 + c2 ) x21 + x22 , and thus is less than 0 except when
m
m
x1 = x2 = 0. By canceling the terms m2 x3 and − m1 x3 in view of the negligible coupling terms in
1

2

x

(8)-(9), we have a11 a22 −a12 a21 = 1/(m1 m2 ) (2c21 +2c1 c2 )(x21 +x22 ) +2c22 (x21 +x22 ) arctan2 ( x2 ) ,
1

which is always positive.
From the above analysis, we show that the linearization matrix A is Hurwitz at an equilibrium
point independent of its specific values, which proves the local asymptotic stability of the average
model given a particular set of control input (α0 , αA , ωα ).

5.6 Application to Semi-analytical Modeling of Steady Turning
5.6.1 Semi-analytical Modeling Framework
The proposed averaging approach is instrumental in several ways. For example, it facilitates
model-based controller design, as was demonstrated in [139], where a backstepping approach is
developed for tracking control of robotic fish. In this section, we present another application of the
average model, where it is used to gain analytical insight into steady turning of robotic fish.
110

Section 5.2 provides a complete dynamic model for a tail-actuated robotic fish, and one can
extract turning information by observation from simulated trajectories as in [14]. However, these
simulations are time-consuming and cannot provide analytical insights into how the turning parameters are related to system parameters or tail-beat parameters. In our group’s prior work [106],
a modeling framework was proposed for computation of steady-state turning motion given asymmetric, periodic body/tail deformation of a robotic fish. In that approach, it is postulated that the
two key parameters of turning motion, the radius and the period, can be obtained by solving the
implicit force and moment balance equations for the averaged, steady-state motion. However, the
modeling framework there adopts two major assumptions: 1) the body direction is aligned with
the motion direction, i.e. the angle of attack remains to be zero during the robot’s motion; 2) the
classical averaging method can be directly applied to robotic fish dynamics (which, as we show in
Section III-C, is generally not true). In this section, we apply the average model obtained with the
proposed approach in this chapter, to develop an semi-analytical modeling framework understanding the turning behavior of robotic fish.
Fig. 5.13 illustrates a robotic fish undergoing a steady turning motion. [X, Y, Z] and [x, y, z]
denote the aforementioned inertial coordinates and body-fixed coordinates, respectively. [x′ , y′ , z′ ]
denote the moving coordinates attached to the fish body, where x′ points to the tangential direction
of the circular path. α is the aforementioned oscillating input (tail angle).
With the averaging method developed in Section IV, the time-explicit oscillating terms are
replaced with constant inputs, particularly, the tail-beat parameters, α0 , αA and ωα in (5.30). Then
the robotic fish can be treated as a rigid body under force and moment balances. The averaged

111

y

yc
Fish body

Tail

E

xc
x

FD
FL

Y

D

z

Z

x

X

Figure 5.13: Top view of the tail-actuated robotic fish undergoing “steady” turning motion.
forces along and perpendicular to the body direction are


 Fx = K f orce m L2 ωα2 α 2 (3 − 3 α 2 − 3 α 2 ) ,
A
12m1
2 0
8 A

 F =K
m 2 2 2
y
f orce 4m L ωα αA α0 .

(5.32)

2

Then the forces along the radial and tangential directions (with respect to the steady circular
path) are


 F ′ = Fx cos β + Fy sin β ,
x

 F = −F sin β + F cos β .
x
y
y′

(5.33)

As the robot is under steady turning, we have the following force and moment balance equa-

112

tions:









Fx′ − FD = 0,

Fy′ − FL + mb ω 2 R = 0,






Mz − MD = 0,

(5.34)

where Mz is the moment relative to the center of mass caused by the tail motion, which is equal to
Kmoment (− 4Jm L2 cωα2 αA2 α0 ) . The unknowns in (5.34) are angular speed ω (which is related to
bz
turning period Tp by Tp = 2ωπ ), turning radius R, and the angle of attack β . Therefore, these turning
parameters can be directly obtained by solving the algebraic equation (5.34), instead of simulating
the differential equation (5.30) and observing its steady-state orbit, as done in Section IV-B and
IV-C. Note that Eq. (5.34) can only be numerically solved (hence explaining our term “semianalytical”), due to its transcendental nature. In particular, we limit the values of β ∈ [0◦ , 20◦] and
use the fsolve command in Matlab to get the solution.

5.6.2 Experimental Evaluation
To validate the steady turning framework, we compare the model predictions and experimental
results. The experiments were conducted in our previous work [15]. As shown in Figures 5.14 5.16, the steady turning framework can capture the original dynamic model very well and predict
the robot’s turning behavior in all cases. Here the turning radius and turning period under “the
average model” are obtained by solving the algebraic equations discussed in Section VI-A, instead
of through integration of the differential equations as in Section IV. We note that the tail-beat
motion here is numerically different (for example, tail-beat frequency is 0.9 Hz) from those in
obtaining the scaling function in Section 5.4, and thus these results provide strong support for the
proposed approach.

113

Turning radius (cm)

30

20

10

0

20

25

30

35

100

Experiment
Dynamic model
Average model

80

Turning period (s)

Experiment
Dynamic model
Average model

40

60
40
20
0

40

20

25

Bias (°)
(a)

30

35

40

Bias (°)
(b)

Figure 5.14: Comparison between model predictions and experimental measurement of turning behavior

Experiment
Dynamic model
Average model

40
30
20
10
0

1

2

100

Turning period (s)

Turning radius (cm)

versus tail-beat bias.

80
60
40
20
0

3

Frequency (Hz)
(a)

Experiment
Dynamic model
Average model

1

2

3

Frequency (Hz)
(b)

Figure 5.15: Comparison between model predictions and experimental measurement of turning behavior
versus tail-beat frequency.

5.7 Nonlinear Controller Design
One promising application of robotic fish is their use as mobile sensing platforms in aquatic environment to monitor water quality and track oil spills. To execute such tasks, robotic fish will need
114

100

Experiment
Dynamic model
Average model

40

Turning period (s)

Turning radius (cm)

50

30
20
10
0

Experiment
Dynamic model
Average model

80
60
40
20

15

20

25

0

30

Amplitude (°)
(a)

15

20

25

30

Amplitude (°)
(b)

Figure 5.16: Comparison between model predictions and experimental measurement of turning behavior
versus tail-beat amplitude.

be able to track a given trajectory or a set of target points. A fish-like swimming robot is a system
that is difficult to control due to its nonlinear, nonholonomic nature.
Most prior work on the control of carangiform robotic fish proposed to use model-free controllers, such as fuzzy logic control [140, 141]. These approaches are useful for engineering application, but they are not amenable to stability analysis, which is fundamental to understanding the
closed-loop system’s behavior. In this work, we investigate target tracking control problem for a
tail-actuated robotic fish, based on the averaged robot dynamic model, as reported in Section 5.4 in
this Chapter. To simplify the problem and emphasize the controller design, the tail-beat frequency
and amplitude are fixed and the only control input is the tail-beat bias α0 .

5.7.1 System Description
Fig. 5.17 shows the target tracking control schematic diagram. Here, (xs , ys ) and (x, y) denote the
locations of the target and the robot’s center of mass, respectively, relative to the inertial coordinate
115

frame.

( xs , ys )

Y

U
T

V

\

r
I

( x, y )

D
X

O



Figure 5.17: Target tracking control schematic diagram.

Define r as the distance between (xs , ys ) and (x, y), and θ as the orientation error, specified by
the angle between the robotic fish heading direction and the line connecting the body center to the
target. The pair (r, θ ) can be considered as the target location in the body-fixed polar coordinates.
We take θ and r as auxiliary state variables, which are expressed as


 r=

x2e + y2e ,

(5.35)


 θ = ψ − ϕ,

where ye = ys − y, xe = xs − x, ϕ = arctan(ye /xe ). The target tracking is achieved if r converges
to zero, which is the goal of the controller’s design in this chapter. The dynamic equations for the

116

two new variables θ and r can be written as


 r˙ = −VCx cos θ +VCy sin θ ,


 θ˙ = VCx sin θ + VCy cos θ + ω .
z
r
r

(5.36)

5.7.2 Hybrid Controller Design

In order to drive r to zero (or the neighborhood of zero), consider a Lyapunov function

1
V1 = r2 .
2

(5.37)

However, we cannot show V˙1 negative no matter what kind of control input is provided. Define a
second Lyapunov equation as
1
V2 = V1 + z21 ,
2

(5.38)

where z1 = kr r − VCx cos θ + VCy sin θ , and kr is a positive coefficient. We check the value of V˙2
and it turns out V˙2 would be negative if there exists an α0 satisfying f (α0 ) = 0, where

f (α0 ) = Aα02 + Bα0 +C = 0

117

(5.39)

and,
1
A = c¯1 k f ωα2 αA2 cos θ ,
2
B = c¯1 k f ωα2 αA2 sin θ ,
1
C = c¯1 k f ωα2 αA2 cos θ (1 − αA2 )
8
+ (k1 + kr )(−VCx cos θ +VCy sin θ ) + (k1 kr )r
− cos θ [VCy ω + fu (VCx ,VCy )] + sin θ [−VCx ω + fv (VCx ,VCy )]
VCy
V
+ (VCx sin θ +VCy cos θ )( Cx sin θ +
+ ω ).
r
r
We can show that there exists a solution of Eq. (5.39) when orientation error |θ | < 45◦ , which
brings in the second controller design.
As shown in [15], the robotic fish converges to a circle when the bias is nonzero, and the larger
the bias, the smaller the radius of the circle will be. If the orientation angle error is not in the
interval where the approaching controller works, we set the bias to the allowable maximum, which
indicates the most effective turning. The sign of the bias will be determined by the sign of θ . The
turning controller can be described as

α0 = α0max sgn(θ ),

(5.40)

where α0max is a positive constant which denotes the allowable maximal tail-beat bias.
To avoid, possible frequent switching between two controllers, a hysteretic switching strategy
is adopted by introducing two thresholds θ1 and θ2 , as shown in Fig. 5.18.

118

Turning
Controller

Approaching
Controller

T1

20$

T2

40$

Figure 5.18: Switching strategy between two controllers.

5.7.3 Simulation and Experimental Results
Extensive simulations have been conducted to verify the effectiveness of the proposed method.
Fig. 5.19 shows the results of target tracking from different initial locations but the same initial
orientations.

Figure 5.19: Simulation results with the target located at the origin and the robotic fish starting from
different initial conditions (denoted by solid blue squares). The locations of controller switching are denotes
by green diamonds.

119

To further verify the proposed method, experiments have been conducted and the setup is
shown in Fig. 5.20. Two LED makers are fixed on the fish body and a camera is installed on
the ceiling above a large indoor water tank. A computer which is connected to the camera is used
to process the acquired images of the two LED markers, to obtain the state estimates for the robot.
The robotic fish can communicate with the computer via Zigbee. The proposed control scheme
is implemented with Visual C++ development tools. The control input is calculated based on the
state information obtained by camera every oscillating period and is then sent to robotic fish via the
Zigbee wireless module. We select several target points randomly in the field of view of camera
and the target point is updated when r reaches 0.02 m. Fig. 5.21 shows final experimental results
which validate the effectiveness of the proposed hybrid controller.
Results from this section have been published in Proceeds of 2013 ASME Dynamics Systems
and Control Conference [139]. Beyond my advisor, the other author of this work Songlin Chen
is an associated professor at Harbin Institute of Technology, Harbin, China, and visited our lab at
2012. Dr. Chen derived most of the back-stepping based hybrid controller and I was in charge of
experimental validation and part of the theoretical derivations.

5.8 Conclusions
In this chapter, a novel averaging method has been presented for dynamics of a carangiform robotic
fish model actuated through a tail. The merit of the approach lies in its simplicity, with the resulting average model amenable to analysis and control design. In particular, the proposed method
builds upon classical averaging, but with the forcing terms in the original dynamic model scaled
appropriately. It is found that the scaling functions take simple forms – constant for the force scaling while linear in bias for the moment scaling. The validity of the approach has been evaluated

120

Camera

Robot

Telosb
(wireless comm.)
Figure 5.20: Experimental setup.
with simulation and experimental results, where steady-state and transient behaviors predicted by
the average model and the original model, respectively, are compared under various scenarios not
used in obtaining the scaling functions. We have further conducted analysis of the average model
in terms of the existence and stability of its equilibria. As an application of the proposed average
model, we have demonstrated a new framework for studying steady turning behavior of a robotic
fish without having to simulate the differential equations for the original dynamic model or the
121

Trajectory
Target

1

Y(m)

0.8

1st Target

4th Target

2nd Target

0.6

0.4

Initail position

0.2
3rd Target

0
0

0.2

0.4

0.6

0.8
X(m)

1

1.2

1.4

Figure 5.21: Experimental results on tracking multiple target points continuously.
average dynamic model.
In the future work, we will explore analytical insight as to why the scaling functions take the
specific simple forms. We will also seek to extend the approach to the dynamic model of robotic
fish actuated with a flexible tail [132], as well as explore the validity of the approach for a wider
class of systems with oscillatory inputs.

122

Chapter 6
Environmental Monitoring with Robotic
Sensor Network: Sampling and
Reconstruction of Low-rank Matrix Fields
6.1 Introduction
Monitoring aquatic environments is of growing interest to public health, marine biology and
ecosystem sustainability. In recent years, aquatic environments have been faced with a variety
of threats from climate change and industrial pollutions [142–144]. For instance, there was an
oil spill in Kalamazoo River, Michigan in 2010, when a pipeline burst and oil flowed into the
river [145]. In another incident, more than 300,000 gallons of diesel fuel spilled in Arthur Kill,
New Jersey after a storage tank was ruptured [146]. During and after these accidents, monitoring
the water quality and detecting pollutants in the water is a vital task [147, 148].
Manual sampling using hand-held devices, is still the most common practice in monitoring
aquatic environment. This approach is labor-intensive and has difficulty in capturing large-scale
spatially distributed phenomena of interest. An alternative approach is in-situ sensing on fixed
buoys [56]. However, since buoyed sensors cannot move around, they have limited adaptability
in monitoring and forecasting dynamic aquatic processes. There is a growing interest in using

123

underwater sensor platforms like autonomous underwater vehicles (AUVs) [74, 149, 150] and sea
gliders [151] to monitor the environment. However, it is not applicable to deploy many of them due
to their high manufacturing and operational cost. In this chapter, we propose to use inexpensive,
highly adaptive and maneuverable robotic fish platforms to sample and reconstruct spatiotemporal
aquatic fields.
Sampling and reconstruction of physical fields using mobile sensor networks have received
significant interest. Early work has dealt with optimal placement of static sensors to reduce the
uncertainty [152], or guidance of one mobile sensor to reduce the field reconstruction error, given
the measurements of the a static sensor network [153]. Field reconstruction using mobile sensor
networks has been studied recently by several groups [154–156].
In this work, we investigate data-adaptive path planning schemes for wireless networks of
robotic sensor platforms. Traditional sampling theory deals with data collection processes that are
independent of the target map to be estimated, aside from some possible prior knowledge and/or
assumptions. On the other hand, sequential, adaptive data collection strategies [157, 158] have
gained increasingly interests due to the fact that the existing observations can be used to guide the
future actions to achieve more efficient solutions.
Field reconstruction using robotic sensor networks has been extensively studied, for fields modeled by Gaussian processes [154, 157–159]. In [157], a data-driven adaptive sampling algorithm is
presented. The next sample locations are selected such that the closest neighbors are able to communicate with each other while the distances between the robots are maximized at the same time.
It also has other rules such as taking samples in area where present sampling density is lower and
avoiding areas where measurements. In [158], the authors present an adaptive sampling algorithm
for multiple autonomous underwater vehicles by assuming that field variables vary linearly with
the location. A cost function is proposed to measure the dissimilarity between current sampling
124

and the next sampling, and the algorithm selects the new sampling locations to minimize the cost
function. In [21], researchers present an adaptive sampling framework to scan and reconstruct
a Gaussian temporal-spatial field using robotic fish sensor networks. In particular, these robots
are controlled to work collaboratively in the form of a swarm to sense the environment in a series of carefully chose rendezvous regions. While all these approaches have their own merits, in
real-world applications, the dynamics of the processes of interest could be very complicated and
sometimes kernel modeling-based techniques may not provide adequate performance.
On the other hand, Compressive Sensing (CS) [160, 161] techniques have been applied to field
monitoring, as they allow the reconstruction of a field with a subset of samples and hence save
time and energy. In [162], trajectory planning algorithms based on the standard CS paradigm are
presented for robots with a variety of mobility patterns, including Random Walk, Random Traveling Salesman, and Randomized Boustrophedon. The target is to minimize the cost function (i.e.,
traveling distance) by comparing these motion patterns and recovering the data using CS. Similar
work has been reported in [163, 164], where the authors assume that the data has a sparse structure
and come up with an optimal sampling solution by leveraging CS techniques. The major drawback
of these techniques is that they do not take into account the underlying data structure obtained
or estimated using real-time measurements, which makes them not amendable to implementing
adaptive sampling strategies.
This work presents a novel adaptive sampling framework to guide mobile robots to conduct
samplings more efficiently in an environmental field of interest, while maintaining the accuracy
of the recovered data. We consider a two-dimensional space (for example, the horizontal plane at
a certain depth), where the field can be represented by a low-rank matrix. The idea is to successively span the sampling area and reconstruct the field analytically. The methods are illustrated in
simulations based on the temperature data collected in a large indoor water tank. The results show
125

that the proposed adaptive sampling scheme is more computationally efficient and requires shorter
travel distances than an existing low-rank matrix completion scheme [165].
The remainder of the chapter is organized as follows. First, the problem formulation is described in Section 6.2. An matrix completion approach using Augmented Lagrange Method (MCALM) is reviewed in Section 6.3. Our proposed adaptive sampling and reconstruction framework
is discussed in Section 6.4. In Section 6.5, simulation results using real trace-driven data are presented to compare the performances of our proposed approach and the MC-ALM method. Discussions of weaknesses of the proposed method and potential solutions are presented in Section 6.6.
Finally, concluding remarks are provided in Section 6.7.

6.2 Problem Formulation
The goal of this work is to develop an active mobile sensing approaching using a group of mobile
robots (for example, robotic fish), to recover an unknown field with minimal energy consumption.
Fig. 6.1a shows the field of oil distribution on the the surface of Kalamazoo River following the
2010 spill. Traditional sampling methods suggest uniform sampling according to the Shannon
sampling theorem to recover the (discretized) field in Fig. 6.1b. While the latter approach is able
to provide the data field with high resolution, it could waste a significant amount of sensing resources, because many measurements collected might not provide any additional information. CS
for sparse matrices and Matrix Completion (MC) [166] for low-rank matrices are two approaches
for recovering the field of interest by collecting only a few random measurements, as illustrated
in Fig. 6.1c. MC is able to reconstruct the data with good fidelity, however, random samples do
not accommodate the mobility constraints of the robots and makes them unsuitable to implement
adaptive sampling strategies.

126

(a)

(b)

(c)
Figure 6.1: Oil spill in the Kalamazoo River. 2010 (Photo Credit: National Wildlife Federation) (a)
Oil map on the surface of Kalamazoo River; (b) discrete uniform sampling; (c) random samplings

127

In this work, we assume that the original field can be represented with a matrix as illustrated in
Fig. 6.1b, where the field values within each grid are identical. Matrix completion is a very useful
tool in many applications, e.g., machine learning [167], control [168], and computer vision [169].
It is often reasonable to assume that the matrix to recover is of low rank. The problem is formulated
as: given a number of mobile robotic sensing platforms (i.e. robotic fish [1,15]), how to drive those
robots to sample the low-rank matrix and recover the missing entries with high accuracy and low
energy consumptions.

6.3 Matrix Completion Using the Augmented Lagrange Multiplier Method
In this section, for the purpose of comparison, we first provide a summary of an existing algorithm that solves the matrix completion problem using Augmented Lagrange Multiplier (MCALM) [165].
The general matrix completion problem can be expressed as

min A ∗ ,
A

subject to Ai j = Di j ,

∀(i, j) ∈ Ω

(6.1)

where · ∗ denotes the nuclear norm of a matrix (sum of its singular values), Ω is the set of indices
of samples, and D denotes the original low-rank matrix. It has been shown that most matrices A
of low rank r can be perfectly recovered by solving the following optimization problem [166],
provided the number p of samples satisfies p ≥ Crn6/5 ln n for some constant C. This bound has
then been improved by the work of others [170, 171]. In [165], Lin et al. proposed to use MCALM algorithm to solve a robust principle component analysis problem and extend it to complete
128

low-rank matrices. The problem is re-formulated as

min A ∗ ,
A

subject to A + E = D,

πΩ (E) = 0,

(6.2)

where πΩ : Rm×n → Rm×n is a linear operator that keeps the entries in Ω unchanged and sets
those outside Ω zeros. Since the error matrix E will compensate for the unknown entries of D, the
unknown entries of D are simply set as zeros. Then the augmented Lagrangian function of (6.2) is
given by [165]:

L(A, E,Y, µ ) = A ∗ + < Y, D − A − E > +

where

µ
D − A − E 2F ,
2

(6.3)

· F is the Frobenius norm, µ is a positive scalar, and Y denotes the optimal Lagrange

multiplier. Following [165], the MC-ALM approach is described in Algorithm 1, where E is
updated with the constraint πΩ (E) = 0 enforced when minimizing L(A, E,Y, µ ).

129

Algorithm 1: Matrix Completion via ALM Method [165]
Inputs: Observation samples Di, j , (i, j) ∈ Ω, of matrix D ∈ Rm×n .
1

Y0 = 0; E0 = 0; µ0 > 0; ρ > 1; k = 0;

2

while not converged do

3

// Lines 4-5 solve Ak+1 = arg minA L(A, Ek ,Yk , µk ).

4

(U, S,V ) = svd(D − Ek + µk−1Yk );

5

Ak+1 = U T −1 [S]V T .
µ
k

6

// Line 7 solves Ek+1 = arg minπ (E)=0 L(Ak+1 , E,Yk , µk ).
Ω

7

Ek+1 = πΩ (D − Ak+1 + µk−1Yk ).

8

Yk+1 = Yk + µk (D − Ak+1 − Ek+1 ).

9

Update µk to µk+1 .

10

k ← k + 1.

11

end while

Output: (Ak , Ek ).
Tε [x] is the soft-thresholding operator:



x − ε,
if x > ε ,




Tε [x] =
x + ε,
if x < −ε ,





 0,
otherwise.

130

6.4 Analytical Reconstruction via Successive Expansion of Sampling Area
A fundamental requirement for the low-rank MC algorithm described in the previous section is
that the original matrix is sampled randomly. Such a requirement is difficult to satisfy for mobile
robots due to their mobility constraints. In this section, we propose a new strategy for sampling
and field-reconstruction by successively growing sub-matrices of the original matrix.
The algorithm is illustrated in Fig. 6.2. For each iterations, new samples (indicated by blue
crosses in Fig. 6.2a) are taken to recover the new row and new column. For instance, consider the
case where a matrix M ∈ R(k+1)×(k+1) needs to be recovered fully based on its known sub-matrix
M ′ ∈ Rk×k and some new samples in its outer row and column. For ease of explanation, we assume
that the rank of M is equal to the rank of M ′ , which is denoted by r, where r ≤ k. M ′ thus has r
linearly independent columns that can be used as a basis B for representing any other column v .
In particular, we can represent v through the set of coefficients with respect to the basis B:
 
c1 
 
 
c2 
 
RepB (v) =  
 .. 
.
 
 
cr
where B =< β1 , ..., βr > and v = c1 β1 + c2 β2 + · · · + cr βr . In M, those basis column vectors are
also independent from each other and the rank of M is r, thus the columns corresponding to B will
form a basis for M ∈ R(k+1)×(k+1) and can be used to recover other columns with RepB (v). The
same discussion applies to the rows. The sampling and reconstruction algorithm is summarized in

131

Algorithm 2.
Algorithm 2: Matrix Completion via Analytical Reconstruction
Inputs: Rank r of the original matrix r
1

Sample the first r × r matrix by initializing k = r

2

while the dimension of the recovered matrix is less than that of the original field, do

3

// Lines 4-12 reconstruct M

4

Find r linearly independent columns as the basis Bcol

5

Take new samples (on the outer row) in those r columns to form the new basis B′col

6

Express the other columns i as a linear combination of Bcol with (ci1 , . . . cir )T

7

Recover the outer row by multiplying coordinates with new basis B′col

8

Find r linearly independent rows as the basis Brow

9

Take new samples (on the outer column) in those r rows to form the new basis B′row

10

Express the other rows j (including the outer row) as a linear combination of
Brow with (c j1 , . . .c jr )T

11

Recover the outer column by multiplying coordinates with new basis B′row

12

Update k to (k + 1)

13

end while

Output: Recovered matrix.

6.5 Trace-Driven Simulations
We evaluate the performance of the MC-ALM method and our proposed algorithms through tracedriven simulations. 1
1 Trace driven simulation is getting the actual data from a real experiment and feeding it to the

simulation

132

y

y

t1

t2

t3

x

tk

ti

x
(a)

(b)
Figure 6.2: Sequential reconstruction of field of interest. (a) Successive expansion of sampling
area, where the new samples in each iteration are only taken from the outer row and column; (b)
Illustration of the sampling strategy with the oil spill example.

6.5.1 Measurement of temperature f eld in a water tank
Experiments are first conducted to measure the distribution of water surface temperature in a water
tank, as illustrated in Fig. 6.3. A 1000W water heater (the red dot point in Fig. 6.3) is fixed in the
tank. The heater is powered on five hours before the measurements and remains on through the
whole experiment, so that the distribution of temperature reaches its equilibrium and is considered
to be time-invariant during the sampling.

133

Tank area

23 inches

Area of Interest

23 inches

Heat source
Figure 6.3: Sketch of the experimental setup for measuring the water-temperature distribution near the
surface.

The 2-D scalar water-temperature field is then described by a matrix A ∈ R24×24 with singular
value decomposition:
A = U SV T ,

S = diag(σ1 , σ2 , . . . , σ24 ),

where the singular values are ordered in a decreasing order. To test the matrix completion algorithms, we approximate A with a low rank matrix A˜ by setting

˜ T,
A˜ = U SV

S˜ = diag(σ1 , σ2 , . . . , σ6 , 0, . . . , 0),

Fig. 6.4 shows the sensing field A˜ with pixels being scaled between 0 and 255.

134

Figure 6.4: Scalar field of water-temperature.

6.5.2 Simulation results
In this subsection we evaluate the performance of MC-ALM and our proposed algorithms using
the collected temperature data.
For the MC-ALM method, 50% measurements are randomly taken to recover the field. Fig. 6.5
shows the random sampling locations and the reconstructed field using those samples. It shows
that the MC-ALM method is able to reconstruct the field with good fidelity. The simulation is
conducted with Matlab on a desktop PC (Dell Vostro 460 with 3.1 GHz Intel i5-2400 CPU and 4
GB memory), and the elapsed CPU time is 0.4 s, obtained using the Matlab macro cputime.
For our proposed iterative adaptive sampling approach, 6 new samples are taken to recover
135

(a)

(b)
Figure 6.5: Reconstruction of field of interest using MC-ALM method. (a) sampled data; (b)
reconstructed field.
the every new row or new column in each iteration, and thus 43.75% measurements are taken
sequentially. The reconstructed data is identical to the original field as shown in Fig. 6.6 and the
136

CPU time for calculation is 0.04 s, which is 10% of that using the ALM-MC method.

Figure 6.6: Reconstruction of field of interest using the proposed method.

6.5.3 Comparisons of energy cost for data sampling
In order to evaluate the energy cost of each sampling strategy, consider a simplified energy model
for robotic fish:
n

E=

∑

j=1

mj

κ1 ∑ Li + κ2 m j

(6.4)

i=1

where E is the expended energy, κ1 and κ2 are coefficients that imply the linear motion cost
and measurement cost, respectively. Li is the moving distance of step i, and m j is the number

137

of measurements of the jth robotic moving platform. Here we neglect other costs, e.g., wireless
communication cost, cost to change the motion direction of the robot, etc., to simplify the problem.
By setting κ1 and κ2 with different values, we can simulate different scenarios. For instance, if

κ1 ≫ κ2 , it implies the movement cost is very expensive while the measurement cost is cheap and
negligible, and then the sampling strategy should drive the robots to have shorter travels. On the
other hand, if κ1 ≪ κ2 , the sampling strategy should minimize the number of samplings to reduce
the total cost. In this work, we primarily consider the first scenario, because that locomotion is the
biggest source of energy expenditure for autonomous robotic fish [1].
For the MC-ALM approach, the question is given a set of sampling locations, what is the
shortest possible route that visits every sampling point and returns to the initial point. In the
simulation, we set n = 2, to denote that the sampling algorithm is conducted with 2 robots. It is
a NP-complete multiple traveling salesman problem (MTSP) and its solution is shown in Fig. 6.4.
The total traveling distance of the two robots is 382 inches after about 50,000 iterations using a
genetic algorithm.
On the other hand, the movement planning of robots using our proposed adaptive sampling
framework is a natural result when determining the next sampling points at every time step, and
the results are shown in Fig. 6.8. The total traveling distance (including back to the initial locations
- a direct line back to the initial point) is 290 inches, which is 24% less than that using the MCALM approach.

138

Total Distance = 382.0066, Iteration = 48827
25

20

15

10

5

0

0

5

10

15

20

25

Figure 6.7: Robot paths to cover all the sampling points using the multiple traveling salesman method. The
red line and the blue denotes the paths of the two robots.

6.6 Discussions
6.6.1 Impact of more accurate energy model
The current work uses the distance to approximate the movement energy cost of the robotic fish
system. In reality, the cost is dependent on not only the distance but also the effort for changing
the direction. Compared with Eq. (6.5), a more accurate energy model is given below:
n

E=

∑

j=1

mj

mj

i=1

i=1

κ1 ∑ Li + κ2 m j + κ3 ∑ n(∆θi )2 ,

139

(6.5)

25

20

15

10

5

0

0

5

10

15

20

25

Figure 6.8: Robot paths to cover all the sampling points using our proposed adaptive sampling algorithm.
The red line and the blue denotes the paths of two robots.

where κ3 is the coefficient to denote the direction change cost and ∆θi is the magnitude of the
direction change in step i. Using this new model for the cost function, the path planning for the
robotic sensor networks using our proposed successive sampling method is no longer intuitive and
needed to be carefully addressed.

140

6.6.2 Limitation of the proposed matrix completion algorithm via successive
reconstruction
If the rank of the initial sampling matrix is less than that of the full field, then in each iteration, the
rank of the new sub-matrix might be different from that of its immediate predecessor. The proposed
algorithms in Section 6.4 will not work for this situation. In this part, we improve the proposed
method as explained in Algorithm 3, with the cost of taking more samplings. We note that the
extra samplings are picked randomly and there is no guarantee that a new linearly independent row
or column (if existing) would be caught using Algorithm 3.

141

Algorithm 3: Matrix Completion via Analytical Reconstruction
Inputs: Rank r of the original matrix r
1

Sample the first r × r matrix by initializing k = r

2

while the dimension of the recovered matrix is less than that of the original field do

3

// Lines 4-16 reconstruct M

4

Calculate the rank r′ of M ′

5

If r = r′ , go to Algorithm 2, else, continue to line 6

6

Find r′ linearly independent columns as the basis Bcol

7

Take samples (on the outer row) in r′ columns to form the new basis B′col and
another column v′col randomly

8

If v′col is independent of B′col , sample every new points, and go to line 11

9

Express other columns i as a linear combination of the basis vectors Bcol with (ci1 , . . .cir )T

10

Recover the outer row by multiplying coordinates with new basis B′col

11

Find r′ linearly independent rows as the basis Brow

12

Take samples (on the outer column) in r′ rows to form the new basis B′row and
another row v′row randomly

13

If v′row is independent of B′row , sample every new points, and go to step 16

14

Express other rows j (including the outer row) as a linear combination of the
basis vectors Brow with (c j1 , . . . c jr )T

15

Recover the outer column by multiplying coordinates with new basis B′row

16

Update k to (k + 1)

17

end while

Output: Recovered matrix.

142

To demonstrate the effectiveness of Algorithm 3, we generate a 20 × 20 matrix (shown in Fig. 6.9)
with rank 8, and its first 10 × 10 sub-matrix has rank 7 and the first 11 × 11 sub-matrix has rank 8.

Figure 6.9: Scalar field described by a 20 × 20 matrix.

Using Algorithm 2, the reconstructed matrix is shown in Fig. 6.10, which has a big error compared with the original scalar field.
Fig. 6.11 shows the reconstructed field using Algorithm 3, which is able to perfectly recover
the original data (matrix). Fig. 6.12 shows the sampling points when using Algorithm 3 to guide
the robotic networks with the red stars denoting the additional samplings that needed to be taken
compared with Algorithm 2.

143

Figure 6.10: Reconstruction of field of interest using Algorithm 2.

6.7 Conclusions
In this work, an adaptive sampling framework is proposed to reconstruct aquatic environmental
fields using a group of robotic sensor platforms. In particular, it is assumed that the field of interest
can be approximated by a low rank matrix, and then an iterative matrix completion approach via
analytical solutions is proposed and explored. For comparison, an MC approach using Augmented
Lagrange Multiplier optimization is also presented. Simulation results show that our proposed
approach is more computationally efficient and requires shorter travel distances for the robots.
The future work can be extended in two directions. First, it will be of interest to explore ignoring weakly linearly independent columns in the rank computation and analyze the propagation of
144

Figure 6.11: Reconstruction of field of interest using Algorithm 3.
reconstruction error caused by this approximation, and the performance of the proposed analytical
matrix reconstruction algorithm. Second, we will explore the validity of the approach using the
robotic fish prototypes developed at the Smart Microsystems Lab.

145

20
18
16
14

y

12
10
8
6
4
2
0

0

5

10

15

x
Figure 6.12: Locations of samplings using Algorithm 3.

146

20

Chapter 7
Conclusions and Future Work
7.1 Conclusions
In this work, the design, development, modeling, control of robotic fish and an adaptive sampling
framework have been studied.
First, the guidelines for robotic fish system design have been provided and demonstrated by
presenting several robot prototypes. Applications with these robots, including aquatic environmental monitoring, robot-animal interaction study, education and outreach, are also presented.
Second, a complete dynamic model for robotic fish actuated by a caudal fin has been developed,
where Lighthill’s elongated-body theory is used to evaluate the tail-generated hydrodynamic force.
It shows that incorporating the body motion in the evaluation of hydrodynamic forces, which is
often ignored in literatures, plays a significant role in capturing the dynamics. Besides, simulation
studies show that the drag force and moment coefficients are strongly dependent on the tail-beat
bias, and subsequently a novel scheme is developed to adapt these coefficients and lift force using
the parameters of tail movement.
Third, a model for robotic fish with a base-actuated flexible tail has been developed. The tail is
modeled as multiple rigid segments connected by springs and dampers, and Lighthill’s elongatedbody theory is used to evaluate the tail-generated hydrodynamic forces. Extensive experimental
results have been collected, and for comparison, a model using linear beam theory is also constructed. From these results, it can be concluded that when the tail is excited under a relatively
147

low frequency, and consequently experiences small deformation, both models produce similar predictions that are close to experimental measurements. However, when the actuation frequency
increases, the two models differ, and the model using multi-segment approximation is able to predict much better the robot speed and the tail shape.
Fourth, a novel control-oriented data-driven averaging approach has been studied that can approximate well the original dynamics to provide better insight for the controller design. To demonstrate the merit of the proposed averaging approach, a semi-analytical framework for studying the
robot’s turning motion and a hybrid controller to make the robot follow a set of waypoints have
been developed based upon the averaged model.
Last, an adaptive sampling framework for a group of robotic sensing platforms has been developed for the mapping of aquatic fields, where the original field is assumed to be represented by
a low-rank matrix. The proposed method successfully expands the sampling area and computes
the missing data points analytically. Real trace-driven simulation results show that the proposed
approach is more computationally efficient and requires shorter travel distances than an established
matrix completion scheme.

7.2 Future work
The future work can be extended in several directions. First, for the averaging of robotic fish dynamics, it is of interest to explore analytical insight as to why the scaling functions take the specific
simple forms. Second, We will also seek to extend the approach to the dynamic model of robotic
fish actuated with a flexible tail, as well as explore the validity of the approach for a wider class
of systems with oscillatory inputs. Third, for the environmental reconstruction using the proposed
successively spanning the sampling area strategy, it will be of interest to explore ignoring weakly

148

linearly independent columns in the rank computation and analyze the propagation of reconstruction error caused by this approximation, and the performance of the proposed analytical matrix
reconstruction algorithm. Forth, we will explore the validity of the approach using the robotic fish
prototypes developed at the Smart Microsystems Lab.

149

BIBLIOGRAPHY

150

BIBLIOGRAPHY
[1] X. Tan, “Autonomous robotic fish as mobile sensor platforms: Challenges and potential
solutions,” Marine Technology Society Journal, vol. 45, no. 4, pp. 31–40, 2011.
[2] F. Zhang, J. Wang, J. Thon, C. Thon, E. Litchman, and X. Tan, “Gliding robotic fish for
mobile sampling of aquatic environment,” in Proceedings of the 11th IEEE International
Conference on Networking, Sensing and Control, Miami, FL, 2014, pp. 167–172.
[3] 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.
[4] 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.
[5] M. S. Triantafyllou and G. S. Triantafyllou, “An efficient swimming machine,” Scientific
America, vol. 273, no. 3, pp. 64–70, 1995.
[6] S. Guo, T. Fukuda, and K. Asaka, “A new type of fish-like underwater microrobot,”
IEEE/ASME Transactions on Mechatronics, vol. 8, no. 1, pp. 136–141, 2003.
[7] 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.
[8] H. Hu, J. Liu, I. Dukes, and G. Francis, “Design of 3D swim patterns for autonomous robotic
fish,” in Proceedings of the 2006 IEEE/RSJ International Conference on Intelligent Robots
and Systems, Beijing, China, 2006, pp. 2406–2411.
[9] B. Kim, D. Kim, J. Jung, and J. Park, “A biomimetic undulatory tadpole robot using ionic
polymer-metal composite actuators,” Smart Materials and Structures, vol. 14, pp. 1579–
1585, 2005.
[10] X. Tan, D. Kim, N. Usher, D. Laboy, J. Jackson, A. Kapetanovic, J. Rapai, B. Sabadus, and
X. Zhou, “An autonomous robotic fish for mobile sensing,” in Proceedings of the IEEE/RSJ
International Conference on Intelligent Robots and Systems, Beijing, China, 2006, pp.
5424–5429.

151

[11] 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.
[12] 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.
[13] K. H. Low, “Locomotion and depth control of robotic fish with modular undulating fins,”
International Journal of Automation and Computing, vol. 4, pp. 348–357, 2006.
[14] J. Wang, F. Alequin-Ramos, and X. Tan, “Dynamic modeling of robotic fish and its experimental validation,” in Proceedings of the 2011 IEEE/RSJ International Conference on
Intelligent Robots and Systems, San Francisco, CA, 2011, pp. 588–594.
[15] J. Wang and X. Tan, “A dynamic model for tail-actuated robotic fish with drag coefficient
adaptation,” Mechatronics, vol. 23, no. 6, pp. 659–668, 2013.
[16] 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.
[17] Z. G. Zhang, N. Yamashita, M. Gondo, A. Yamamoto, and T. Higuchi, “Electrostatically
actuated robotic fish: Design and control for high-mobility open-loop swimming,” IEEE
Transaction on Robotics, vol. 24, no. 1, pp. 118–129, 2008.
[18] E. Kanso and P. K. Newton, “Passive locomotion via normal-mode coupling in a submerged
spring-mass system,” J. Fluid Mech, vol. 641, pp. 205–215, 2009.
[19] 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.
[20] Y. Wang, R. Tan, G. Xing, J. Wang, and X. Tan, “Profiling aquatic diffusion process profiling
using robotic sensor networks,” IEEE Transactions on Mobile Computing, vol. 13, no. 4,
2014.
[21] Y. Wang, R. Tan, G. Xing, X. Tan, J. Wang, and R. Zhou, “Spatiotemporal aquatic field
reconstruction using cyber-physical robotic sensor systems,” ACM Transactions on Sensor
Networks, vol. 10, no. 57, 2014.

152

[22] Y. Wang, R. Tan, G. Xing, J. Wang, and X. Tan, “Accuracy-aware aquatic diffusion process
profiling using robotic sensor networks,” in 11th ACM/IEEE Conference on Information
Processing in Sensor Networks, Beijing, China, 2012, pp. 281–292.
[23] Y. Wang, R. Tan, G. Xing, X. Tan, J. Wang, and R. Zhou, “Spatiotemporal aquatic field
reconstruction using robotic sensor swarm,” in Proceedings of the 33st IEEE Real-Time
Systems Symposium (RTSS), San Juan, Puerto Rico, 2012, pp. 205–214.
[24] 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.
[25] M. J. Lighthill, “Large-amplitude elongated-body theory of fish locomotion,” Proceedings
of the Royal Society of London B, vol. 179, pp. 125–138, 1971.
[26] T. Y. Wu, “Hydromechanics of swimming propulsion. Part 1. Swimming of a twodimensional flexible plate at variable forward speeds in an inviscid fluid,” Journal of Fluid
Mechanics, vol. 46, pp. 337–355, 1971.
[27] B. Ahlborn, S. Chapman, R. Stafford, R. W. Blake, and D. G. Harper, “Experimental simulation of the thrust phase of fast-start swimming of fish,” Journal of Experimental Biology,
vol. 200, pp. 2301–2312, 2007.
[28] I. L. Y. Spierts and J. L. Van Leeuwen, “Kinematics and muscle dynamics of C- and S-starts
of carp (cyprinus carpio l.),” Journal of Experimental Biology, vol. 202, pp. 393–406, 1999.
[29] U. K. Muller, E. J. Stamhuis, and J. J. Videler, “Riding the waves: The role of the body wave
in undulatory fish swimming,” Integrative and Comparative Biology, vol. 42, pp. 981–987,
2002.
[30] M. S. Triantafyllou, D. K. P. Yue, and G. S. Triantafyllou, “Hydrodynamics of fishlike
swimming,” Annu. Rev. Fluid Mech., vol. 32, pp. 33–53, 2000.
[31] M. S. Triantafyllou, A. H. Techet, and F. S. Hover, “Review of experimental work in
biomimetic foils,” IEEE Journal of Oceanic Engineering, vol. 29, no. 3, pp. 585–594, 2004.
[32] H. Liu, R. J. Wassersug, and K. Kawachi, “A computational fluid dynamics study of tadpole
swimming,” Journal of Experimental Biology, vol. 199, pp. 1245–1260, 1996.
[33] R. Mittal, “Computational modeling in biohydrodynamics: Trends, challenges, and recent
advances,” IEEE Journal of Oceanic Engineering, vol. 29, no. 3, pp. 595–604, 2004.
153

[34] D. Byun, J. Choi, K. Cha, J. Park, and S. Park, “Swimming microrobot actuated by two
pairs of helmholtz coils system,” Mechatronics, vol. 21, no. 1, pp. 357–364, 2011.
[35] 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.
[36] 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.
[37] P. R. Bandyopadhyay, “Maneuvering hydrodynamics of fish and small underwater vehicles,”
Integrative and Comparative Biology, vol. 42, pp. 102–117, 2002.
[38] 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.
[39] 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.
[40] 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.
[41] 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.
[42] 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.
[43] 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.
[44] 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.
[45] 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.
154

[46] 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.
[47] 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.
[48] 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.
[49] 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.
[50] J. W. Paquette and K. J. Kim, “Ionomeric electroactive polymer artificial muscle for naval
applications,” IEEE Journal of Oceanic Engineering, vol. 29, no. 3, pp. 729–738, 2004.
[51] T. Tao, Y.-C. Liang, and M. Taya, “Bio-inspired actuating system for swimming using shape
memory alloy composites,” International Journal of Automation and Computing, vol. 4, pp.
366–373, 2006.
[52] J. Tangorra, P. Anquetil, T. Fofonoff, A. Chen, M. Del Zio, and I. Hunter, “The application
of conducting polymers to a biorobotic fin propulsor,” Bioinspiration & Biomimetics, vol. 2,
pp. S6–S17, 2007.
[53] T. Wiguna, S. Heo, H. C. Park, and N. S. Goo, “Design and experimental parameteric study
of a fish robot actuated by piezoelectric actuators,” Journal of Intelligent Material Systems
and Structures, vol. 20, pp. 751–758, 2009.
[54] Y. Zhang, M. Cong, D. Guo, and D. Wang, “Design optimization of a bidirectional microswimming robot using giant magnetostrictive thin films,” IEEE/ASME Transactions on
Mechatronics, vol. 14, no. 4, pp. 493–503, 2009.
[55] P. Valdastri, E. Sinibaldi, S. Caccavaro, G. Tortora, A. Menciassi, and P. Dario, “A novel
magnetic actuation system for miniature swimming robots,” IEEE Transactions on Robotics,
vol. 27, no. 4, pp. 769–779, 2011.
[56] 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 Observa155

tion Network (ReCON),” in Proceedings of Marine Technology Society/IEEE Oceans 2007
Conference, 2007, 6 pp.
[57] 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.
[58] 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.
[59] H. Hongji, A. Kaneko, and K. Kawatate, “Self-governing profiling system,” Cont. Shelf
Res., vol. 7, pp. 1257–1265, 1987.
[60] 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.
[61] R. A. Luettich, “PSWIMS, a profiling instrument system for remote physical and chemical
measurements in shallow-water,” Estuaries, vol. 16, pp. 190–197, 1993.
[62] P. R. Bandyopadhyay, “Trends in biorobotic autonomous undersea vehicles,” IEEE Journal
of Oceanic Engineering, vol. 30, no. 1, pp. 109–139, 2005.
[63] 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.
[64] 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.
[65] 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.
[66] 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.

156

[67] 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.
[68] 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.
[69] 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.
[70] 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.
[71] 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.
[72] 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.
[73] 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.
[74] ScienceDaily, “Rensselaer researchers experiment with solar underwater robots,” December
2004. [Online]. Available: http://www.sciencedaily.com/releases/2004/12/041212081548.
htm
[75] 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.
[76] 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.

157

[77] 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.
[78] 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.
[79] 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.
[80] 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.
[81] 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.
[82] 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.
[83] 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
experiments using a remotely operated vehicle (ROV),” J. Exp. Marine Biol. Ecol., vol. 395,
no. 1-2, pp. 37–47, 2010.
[84] 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.
[85] F. Zhang and X. Tan, “Gliding robotic fish and its tail-enabled yaw motion stabilization using sliding mode control,” in Proceedings of the 2013 ASME Dynamic Systems and Control
Conferece, Palo Alto, CA, 2013, paper DSCC2013-4015 (10 pp).

158

[86] K. Cannon, M. Lapoint, N. Bird, K. Panciera, H. Veeraraghavan, N. Papanikolopoulos, and
M. Gini, “Using robots to raise interest in technology among underrepresented groups,”
Robotics Automation Magazine, IEEE, vol. 14, no. 2, pp. 73–81, 2007.
[87] N. Abaid, V. Kopman, and M. Porfiri, “An attraction toward engineering careers: The story
of a brooklyn outreach program for k-12 students,” Robotics Automation Magazine, IEEE,
vol. 20, no. 2, pp. 31–39, 2013.
[88] E. Matson, S. DeLoach, and R. Pauly, “Building interest in math and science for rural and
underserved elementary school children using robots,” Journal of STEM Education, vol. 5,
no. 3, pp. 35–46, 2004.
[89] H. Hu, “Biologically inspired design of autonomous robotic fish at essex,” in Proceedings
of the IEEE SMC UK-RI Chapter Conference, Sheffield, UK, 2006, pp. 1–6.
[90] K. A. Harper, M. D. Berkemeier, and S. Grace, “Modeling the dynamics of spring-driven
oscillating-foil propulsion,” IEEE Journal of Oceanic Engineering, vol. 23, no. 3, pp. 285–
296, 1998.
[91] S. D. Kelly, “The mechanics and control of robotic locomotion with applications to aquatic
vehicles,” Ph.D. dissertation, California Institute of Technology, 1998.
[92] S. D. Kelly and R. M. Murray, “Modelling efficient pisciform swimming for control,” International Journal of Robust and Nonlinear Control, vol. 10, pp. 217–241, 2000.
[93] R. Mason and J. W. Burdick, “Experiments in carangiform robotic fish locomotion,” in
Proceedings of the 2000 IEEE International Conference on Robotics & Automation, San
Francisco, CA, 2000, pp. 428–435.
[94] R. Mason, “Fluid locomotion and trajectory planning for shape-changing robots,” Ph.D.
dissertation, California Institute of Technology, 2003.
[95] E. Kanso, J. E. Marsden, C. W. Rowley, and J. B. Melli-Huber, “Locomotion of articulated
bodies in a perfect fluid,” Nonlinear Science, vol. 15, pp. 255–289, 2005.
[96] F. Boyer, M. Porez, A. Leroyer, and M. Visonneau, “Fast dynamics of an eel-like robot
c Comparisons with Navier-Stokers simulations,” IEEE Transactions on Robotics, vol. 24,
no. 6, pp. 1274–1288, 2008.
[97] Y. Or, S. Zhang, and R. M. Murray, “Dynamics and stability of low-reynolds-number swimming near a wall,” SIAM J. Applied Dynamical Systems, vol. 10, no. 3, pp. 1013–1041,
2011.
159

[98] K. A. McIsaac and J. P. Ostrowski, “Motion planning for anguilliform locomotion,” IEEE
Transactions on Robotics and Automation, vol. 19, no. 4, pp. 637–652, 2003.
[99] J. E. Colgate and K. M. Lynch, “Mechanics and control of swimming: A review,” IEEE
Journal of Oceanic Engineering, vol. 29, no. 3, pp. 660–673, 2004.
[100] J. Guo and Y.-J. Joeng, “Guidance and control of a biomimetic autonomous underwater
vehicle using body-fin propulsion,” Proceedings of the Institution of Mechanical Engineers,
Part M: Journal of Engineering for the Maritime Environment, vol. 218, pp. 93–111, 2004.
[101] S. Saimek and P. Y. Li, “Motion planning and control of a swimming machine,” International Journal of Robotics Research, vol. 23, no. 1, pp. 27–53, 2004.
[102] J. B. Melli, C. W. Rowley, and D. S. Rufat, “Motion planning for an articulated body in a
perfect planar fluid,” SIAM J. Applied Dynamical Systems, vol. 5, no. 4, pp. 650–669, 2006.
[103] J. Cochran, A. Siranosian, N. Ghods, and M. Krstic, “3-D source seeking for underactuated
vehicles without position measurement,” IEEE Transactions on Robotics, vol. 25, no. 1, pp.
117–129, 2009.
[104] Y.-H. Zhang, J.-H. He, J. Yang, and S.-W. Zhang, “A computational fluid dynamics (CFD)
analysis of an undulatory mechanical fin driven by shape memory alloy,” International Journal of Automation and Computing, vol. 4, pp. 374–381, 2006.
[105] M. Anton, Z. Chen, M. Kruusmaa, and X. Tan, “Analytical and computational modeling of
robotic fish propelled by soft actuation material-based active joints,” in Proceedings of the
2009 IEEE/RSJ International Conference on Intelligent Robots and Systems, St. Louis, MO,
2009, pp. 2126–2131.
[106] X. Tan, M. Carpenter, J. Thon, and F. Alequin-Ramos, “Analytical modeling and experimental studies of robotic fish turning,” in Proceedings of the IEEE International Conference on
Robotics and Automation, Anchorage, AK, 2010, pp. 102–108.
[107] H. Chen, C. Zhu, X. Yin, X. Xing, and G. Cheng, “Hydrodynamic analysis and simulation
of a swimming bionic robot tuna,” Journal of Hydrodynamics, Ser. B, vol. 19, no. 4, pp.
412–420, 2009.
[108] 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, Oct. 2001.

160

[109] P. Bhatta and N. E. Leonard, “A Lyapunov function for vehicles with lift and drag: stability
of gliding,” in Proceedings of the 43rd IEEE Conference on Decision and Control, Dec.
2004, pp. 4101–4106 Vol.4.
[110] K. A. Morgansen, V. Duindam, R. J. Mason, J. W. Burdick, and R. M. Murray, “Nonlinear
control methods for planar carangiform robot fish locomotion,” in Proceedings of the 2001
IEEE International Conference on Robotics & Automation, Seoul, Korea, 2001, pp. 427–
434.
[111] T. I. Fossen, Guidance and Control of Ocean Vehicles. New York: Wiley, 1994.
[112] H. N. Arafat, D. J. Stilwell, and W. L. Neu, “Development of a dynamic model of a small
high-speed autonomous underwater vehicle,” in Proceedings of Oceans, Boston, MA, 2006,
pp. 1–6.
[113] J. Yu, L. Liu, L. Wang, M. Tan, and D. Xu, “Turning control of a multilink biomimetic
robotic fish,” IEEE Transactions on Robotics, vol. 24, no. 1, pp. 201–206, 2008.
[114] F. Liu, K.-M. Lee, and C.-J. Yang, “Hydrodynamics of an undulating fin for a wave-like
locomotion system design,” IEEE/ASME Transactions on Mechatronics, vol. 17, no. 3, pp.
554–562, 2012.
[115] S. Lee, J. Park, and C. Han, “Optimal control of a mackerel-mimicking robot for energy
efficient trajectory tracking,” Journal of Bionic Engineering, vol. 4, no. 4, pp. 209–215,
2007.
[116] D. S. Barrett, M. S. Triantafyllou, D. K. P. Yue, M. A. Grosenbaugh, and M. J. Wolfgang,
“Drag-reduction in fish-like locomotion,” Journal of Fluid Mechanics, vol. 392, pp. 183–
212, 1999.
[117] J. H. Long, Jr., T. J. Koob, K. Irving, K. Combie, V. Engel, H. Livingston, A. Lammert, and
J. Schumacher, “Biomimetic evolutionary analysis: testing the adaptive value of vertebrate
tail stiffness in autonomous swimming robots,” Journal of Experimental Biology, vol. 209,
no. 23, pp. 4732–4746, 2006.
[118] J. L. Tangorra, G. V. Lauder, I. W. Hunter, R. Mittal, P. G. A. Madden, and M. Bozkurttas,
“The effect of fin ray flexural rigidity on the propulsive forces generated by a biorobotic fish
pectoral fin,” Journal of Experimental Biology, vol. 213, pp. 4043–4054, 2010.
[119] T. R. Kane, R. Ryan, and A. K. Banerjee, “Dynamics of a cantilever beam attached to a
moving base,” Journal of Guidance, Control, and Dynamics, vol. 10, no. 2, pp. 139–151,
1987.
161

[120] T. P. Mitchell and J. C. Bruch, “Free vibrations of a flexible arm attached to a compliant
finite hub,” Journal of Vibration, Acoustics Stress and Reliability in Design, vol. 1, no. 1,
pp. 118–120, 1988.
[121] K. H. Low, “Eigen-analysis of a tip-loaded beam attached to a rotating joint,” Journal of
Vibration and Acoustics, vol. 112, no. 4, pp. 497–500, 1990.
[122] V. Kopman and M. Porfiri, “Design, modeling, and characterization of a miniature robotic
fish for research and education in biomimetics and bioinspiration,” IEEE Transactions on
Mechatronics, vol. 18, no. 2, pp. 471–483, 2013.
[123] H. Yang, J. Hong, and Z. Yu, “Dynamics modelling of a flexible hub-beam system with a
tip mass,” Journal of Sound and Vibration, vol. 266, no. 4, pp. 759–774, 2003.
[124] P. Shi, J. McPhee, and G. R. Heppler, “A deformation field for Euler-Bernoulli beams with
applications to flexible multibody dynamics,” Multibody System Dynamics, vol. 5, no. 1, pp.
79–104, 2001.
[125] A. K. Banerjee and S. Nagarajan, “Efficient simulation of large overall motion of beams
undergoing large deflection,” Multibody System Dynamics, vol. 1, no. 1, pp. 113–126, 1997.
[126] N. Wereley, G. Wang, and A. Chaudhuri, “Demonstration of uniform cantileverd beam beding vibration using a pair of piezoelectric actuators,” J. Intell. Mater. Syst. Struct., vol. 22,
no. 4, pp. 307–316, 2011.
[127] X. Chen, G. Zhu, X. Yang, D. L. S. Hung, and X. Tan, “Model-based estimation of flow
characteristics using an ionic polymer-metal composite beam,” IEEE/ASME Transactions
on Mechatronics, vol. 18, no. 3, pp. 932–943, 2013.
[128] T. J. Pedley and S. J. Hill, “Kinematic condition for maximizing the thrust of a robotic fish
using a compliant caudal fin,” The Journal of Experimental Biology, vol. 202, pp. 3431–
3438, 1999.
[129] J. Cheng, L. Zhuang, and B. Tong, “Analysis of swimming three-dimensional waving
plates,” Journal of Fluid Mechanics, vol. 232, pp. 341–355, 1991.
[130] R. W. Clough and J. Penzien, Dynamics of Structures, 3rd ed. 1995 University Ave, Berkeley, CA: Computers & Structures Inc, 2003.
[131] V. Kopman and M. Porfiri, “A miniature and low-cost robotic fish for ethorobotics research
and engineering education: Iłbioinspired design,” in Proceedings of the ASME 2011 Dynamic Systems and Control Conference, Oct. 2011, pp. 209–216.
162

[132] J. Wang, P. K. McKinley, and X. Tan, “Dynamic modeling of robotic fish with a flexible
caudal fin,” in Proceedings of the 2012 ASME Dynamic Systems and Control Conference,
Fort Lauderdale, FL, 2012, dSCC2012-MOVIC2012-8695.
[133] K. A. Morgansen, P. A. Vela, and J. W. Burdick, “Trajectory stabilization for a planar
carangiform robot fish,” in Proceedings of the 2002 IEEE International Conference on
Robotics and Automation, Washington, DC, May 2002, pp. 756–762.
[134] J. A. Sanders, F. Verhulst, and J. Murdock, Averaging Methods in Nonlinear Dynamical
Systems. New York: Springer, 2007.
[135] P. A. Vela, K. A. Morgansen, and J. W. Burdick, “Second order averaging methods for
oscillatory control of underactuated mechanical systems,” in Proceedings of the American
Control Conference, Anchorage, AK, May 2002, pp. 576–581.
[136] F. Bullo and A. D. Lewis, Geometric Control of Mechanical Systems: Modeling, Analysis,
and Design for Simple Mechanical Control Systems. New York: Springer, 2004.
[137] J. Liu and H. Hu, “Mimicry of sharp turning behaviours in a robotic fish,” in Proceedings
of the IEEE International Conference on Robotics and Automation, Barcelona, Spain, 2005,
pp. 3318–3323.
[138] J. Yu, L. Liu, and L. Wang, “Dynamics and control of turning manerver for biomimetic
robotic fish,” in Proceedings of 2006 the IEEE/RSJ International Conference on Intelligent
Robots and Systems, Beijing, China, 2006, pp. 3318–3323.
[139] S. Chen, J. Wang, and X. Tan, “Backstepping-based hybrid target tracking control for a
carangiform robotic fish,” in Proceedings of the 2013 ASME Dynamic Systems and Control
Conferece, Palo Alto, CA, 2013, paper DSCC2013-3963 (9 pp).
[140] Y. Hu, W. Zhao, L. Wang, and Y. Jia, “Underwater target following with a vision-based
autonomous robotic fish,” in Proceedings of the 2009 American Control Conference, St.
Louis, MO, 2009, pp. 5265–5270.
[141] J. Yu, S. Wang, and M. Tan, “Basic motion control of a free-swimming biomimetic robot
fish,” in Proceedings of the 42nd IEEE International Conference on Decision and Control,
Maui, HI, 2003, pp. 1268–1273.
[142] 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.
163

[143] N. F. Gray et al., Water technology: an introduction for environmental scientists and engineers. IWA Publishing, 2010, no. Ed. 3.
[144] R. Rajagopalan, Environmental studies: from crisis to cure. Oxford University Press, 2011.
[145] USEPA, “EPA Response to Enbridge Spill in Michigan.” [Online]. Available: http:
//www.epa.gov/enbridgespill/documents.html
[146] Wikipedia, “Oil Spill.” [Online]. Available: http://en.wikipedia.org/wiki/Oil spill
[147] 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
[148] Office of the Maritime Administrator, “Deepwater horizon marine casualty investigation
report,” August 2011, Retrieved on February 25, 2013.
[149] B. Chen and D. Pompili, “Team formation and steeing algorithms for underwater gliaders
using acoustic communication,” Computer Communications, vol. 35, no. 9, pp. 1017–1028,
2012.
[150] B. Chen, P. Pandey, and D. Pompili, “A distributed adaptive sampling soluting using autonomous underwater vehicles,” in Proceedings of the Seventh ACM International Conference on Underwater Networks and Systems, New York, USA, 2012, article No. 29.
[151] 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.
[152] A. Krause, C. Guestrin, A. Gupta, and J. Kleinberg, “Near-optimal sensor placements: Maximizing information while minimizing communication cost,” in The 5th International Conference on Information Processing in Sensor Networks (IPSN). ACM/IEEE, 2006, pp.
2–10.
[153] B. Zhang and G. Sukhatme, “Adaptive sampling for estimating a scalar field using a robotic
boat and a sensor network,” in IEEE International Conference on Robotics and Automation
(ICRA). IEEE, 2007, pp. 3673–3680.
[154] Y. Xu, J. Choi, and S. Oh, “Mobile sensor network navigation using gaussian processes
with truncated observations,” IEEE Transactions on Robotics, vol. 27, no. 6, pp. 1118–1131,
2011.
164

[155] J. Cort´es, “Distributed kriged kalman filter for spatial estimation,” IEEE Transactions on
Automatic Control, vol. 54, no. 12, pp. 2816–2827, 2009.
[156] K. Low, J. Dolan, and P. Khosla, “Adaptive multi-robot wide-area exploration and mapping,”
in The 7th International Conference on Autonomous Agents and Multiagent Systems, 2008,
pp. 23–30.
[157] A. Munafo, E. Simetti, A. Turetta, A. Caiti, and G. Casalino, “Autonomous underwater
vehicle teams for adaptive ocean sampling: A data-driven approach,” Ocean Dynamics,
vol. 61, no. 11, pp. 1981–1994, 2011.
[158] D. Popa, A. Sanderson, R. Komerska, S. Mupparapu, D. Blidberg, and S. Chappel, “Adaptive sampling algorithms for multiple autonomous underwater vehicles,” in Proceedings of
IEEE INternational Conference on Autonomous Underwater Vehicles, Maine, USA, 2004,
pp. 108–118.
[159] Y. Xu and J. Choi, “Adaptive sampling for learning gaussian processes using mobile sensor
networks,” Sensors, vol. 11, no. 3, pp. 3051–3066, 2011.
[160] D. L. Donoho, “Compressed sensing,” IEEE Transactions on Information Theory, vol. 52,
no. 4, pp. 1289–1306, 2006.
[161] E. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Transactions on Information
Theory, vol. 52, no. 2, pp. 489–509, 2006.
[162] R. Hummel, S. Poduri, F. Hover, U. Mitra, and G. Sukahatme, “Mission design for compressive sensing with mobile robots,” in Proceedings of 2011 IEEE International Conference on
Robotics and Automation, Shanghai, China, 2011, pp. 2362–2367.
[163] S. Huang and J. Tan, “Adaptive sampling using mobile sensor networks,” in Proceedings of
2012 IEEE International Conference on Robotics and Automation, Minnesota, USA, 2012,
pp. 657–662.
[164] ——, “Compressive mobile sensing in robotic mapping,” in Proceedings of 2009 the
IEEE/RSJ International Conference on Intelligent Robots and Systems, Saint Louis, USA,
2009, pp. 3070–3075.
[165] Z. Lin, M. Chen, and Y. Ma, “The augmented lagrange multiplier method for exact recovery
of corrupted low-rank matrices,” 2009, uIUC Technical Report UILU-ENG-09-2215, 2009.

165

[166] E. Candes and B. Recht, “Exat matrix completion via convex optimization,” Foundations of
Computational Mathematics, vol. 9, pp. 717–772, 2009.
[167] A. Argyriou, T. Evgeniou, and M. Pontil, “Convex multi-task feature learning,” Machine
Learning, vol. 73, no. 3, pp. 243–272, 2008.
[168] M. Mesbahi and G. Papavassilopoulos, “On the rank minimization problem over a positive
semidefinite linear matrix inequality,” IEEE Transactions on Automatic Control, vol. 42,
no. 2, pp. 239–243, 1997.
[169] C. Tomasi and T. Kanade, “Shape and motion from image streams under orthography: a
factorization method,” International Journal of Computer Vision, vol. 9, no. 2, pp. 137–154,
1992.
[170] J. Cai, E. J. Candes, and Z. Shen, “A singular value thresholding algorithm for matrix completion,” SIAM Journal on Optimization, vol. 20, no. 4, pp. 1956–1982, 2010.
[171] K. Toh and S. Yun, “An accelerated proximal gradient algorithm for nuclear norm regularized least squares problems,” Pacific Journal of Optimization, vol. 6, pp. 615–640, 2010.

166