Three-Dimensional Modeling of a Robotic Fish Based on Real Carp Locomotion

Abstract

This work focuses on developing a complete non-linear dynamic model comprising entirely kinematic and hydrodynamic effects of Carangiform locomotion based on the Lagrange approach by adapting the parameters and behaviors of a real carp. In order to imitate biological features, swimming patterns of a real carp for forward, turning and up-down motions are analyzed by using the Kineova 8.20 software. The proportional optimum link lengths according to actual size, swimming speed, flapping frequency, proportional physical parameters and different swimming motions of the real carp are investigated with the designed robotic fish model. Three-dimensional (3D) locomotion is evaluated by tracking two trajectories in a MATLAB environment. A Reaching Law Control (RLC) approach for inner loop (Euler angles-speed control) and a guidance system for the outer loop (orientation control) are proposed to provide an effective closed-loop control performance. In order to illustrate the 3D performance of the proposed closed loop control system in a virtual reality platform, the designed robotic fish model is also implemented using the Virtual Reality Modeling Language (VRML). Simulation and experimental analysis show that the proposed model gives us significant key solutions to design a fish-like robotic prototype.

1. Introduction

In recent years, biologically inspired behavior-based systems have been more and more popular topic in underwater vehicles. The goal of the bio-inspired approach is to adapt the biological features of underwater creatures such as fish to Autonomous Underwater Vehicle (AUV) designs and also imitate the aquatic locomotion abilities of them. Biomimetic modeling also provides enhanced performance and increased efficiency, maneuverability and acceleration for novel AUV designs [1,2,3,4,5]. As it is known, many kinds of fish perform high-efficiently locomotion and maneuvering in the water. Propulsion efficiencies of the rotary propeller AUVs are limited below 70%, while the swimming mechanism of a real fish is 20% more efficient than rotary propellers [6]. In addition, the turning radius of propellers is big and speed performances are low. Therefore, this kind of propulsion is more noisy and ineffective than bio-inspired systems [6,7]. These advantages have great benefit for marine applications and fish-like robots have evolved to provide versatile solutions for a wide variety of underwater applications, such as undersea investigation, pollution detection, deep sea monitoring and mapping, military operations and protections, etc. [8,9,10,11].

An elaborate observation of the whole biological structure of a fish is necessary to identify the principles of a bio-inspired approach to the fish. In this way, significant limitations and constraints of any biological feature are determined before modeling the robotic fish. A bio-inspired robotic fish is defined as an aquatic vehicle, which is propelled by undulatory or oscillatory body motions [10,12]. In nature, most fish exhibit various swimming behaviors by bending their bodies and/or their caudal fins (BCF). Alternatively, some fish use their median and/or pectoral fins (MPF) [13]. Commonly, more than 85% of fish swim by BCF locomotion and on average 15% of fish swim in the world by MPF. Carangiform locomotion, the most common type of BCF, generates undulatory swimming patterns, high speed performance and low noise [11,13]. Thus, biomimetic Carangiform design is an appropriate solution for AUVs [14,15,16,17].

There have been some studies about robotic fish modeling in the literature. The first mathematical models including fish swimming were proposed by Wu and Lighthill in 1960s. Wu suggested simplified two-dimensional waving plate theory and Lighthill defined the body-traveling wave, which takes a sinusoidal form of undulatory motion [1,18]. Sfakiotakis explained the fish-swimming modes for oscillatory and undulatory motions of fins and bodies [14]. McIsaac and Ostrowski established a model for Anguilliform locomotion to give swimming gaits [19], while Morgansen et al. derived 3D models for pectoral fins and tails [20]. Also, Yu et al. applied the Central Pattern Generator (CPG)-based locomotor controller to their dynamic model. CPG is a neural circuit diagram inspired from biological musculoskeletal systems to perform gait generation of bio-inspired robots [21,22,23]. Liu and Hu built a simulator for the autonomous Carangiform motion of a four-joint robotic fish, but up-down motions were not modeled [24]. Yu et al. also modeled up-down motions [25,26]. Unlike these, Zhou et al. presented a dynamic backward analysis for a biomimetic Carangiform robotic fish [27]. Furthermore, Nakashima and Ohgishi presented a dynamic model of the double-jointed robotic fish to examine self-exciting conditions and propulsive characteristics [28]. Kim et al. presented the four-link motion mechanism of fish [29]. Masoomi and colleagues modeled a four degrees of freedom (4-DoF) mathematical model including hydrodynamic forces inspired by tuna fish [10]. A multi-DoF robotic fish model was suggested by Suebsaiprom and Lin to determine the position and orientation of the rigid body. They used a state feedback controller in the simulation [15,30]. According to the above literature survey, a 3D dynamic model including all kinematic and hydrodynamic effects by analyzing realistic parameters of a real fish would be a significant reference while developing a realistic prototype.

The contribution of this paper is to develop a complete dynamic model including kinematic and hydrodynamic effects of Carangiform locomotion based on the Lagrange approach by adapting the parameters and behaviors obtained from real carp. This model provides appropriate solutions for biomimetic design of fish-like motions with 3D carp gait patterns unlike the dynamic robotic fish models available in the literature with non-realistic parameters. In order to mimic real fish-swimming abilities, the designed robotic fish model has been focused on the proportional optimum link lengths according to actual size, swimming speed, flapping frequency, proportional physical parameters and different swimming motions of the real carp. 3D locomotion is evaluated by applying various trajectories in different MATLAB environments. The closed loop performance is provided by using the Reaching Law Control (RLC) approach in the inner loop and orientation control is also achieved with a guidance system in the outer loop. These simulations and experimental analyses show that the proposed model gives us significant key solutions to design a fish-like robotic prototype.

The rest of this paper is organized as follows: in the next section, Carangiform locomotion is described and a complete dynamic model, including kinematic and hydrodynamic effects, is given. Then the maneuverability of the robotic fish model is performed and simulation results are given. Finally, conclusions are presented.

2. Three Dimensional Dynamic Model and Motion Equations of the Robotic Fish

The fish model analyzed in this study is inspired by a Carangiform-mode fish based on the BCF locomotion. A Carangiform mode exhibits more significant swimming actions such as fast start, rapid turning, C-shape turning, cruise swimming and high accelerations [14]. A Carangiform fish swims on the lateral sinusoidal motion, which increases the amplitude from nose to caudal fin in a spine. This kinematic motion can be described using the body-traveling wave function, suggested by Lighthill [18,31]. However, this function only includes the simple kinematic equality for the fish tail and there is no whole dynamic effect to understand the fish locomotion. Based on the biological features of Carangiform locomotion, the non-linear dynamic model of the robotic fish should be derived, including 6-DoF motion. This study focuses on deriving a complete non-linear mathematical model of a Carangiform robotic fish by analyzing real carp locomotion. The dynamic model is performed by using realistic parameters from 3D swimming patterns of the carp. The proposed non-linear mathematical model of the robotic fish includes the body and tail kinematics, hydrodynamic effects acting on each link and rigid body motions. This model behaves like a real fish with complete non-linear derivations. In this way, biomimetic modeling of the robotic fish prototype can be easily implemented in future work.

2.1. Carp Locomotion

In this study, forward, turning and ascending-descending swimming patterns are investigated with top- and side-view cameras for 3D records in a 120 cm length, 80 cm width and 60 cm depth test tank to analyze the maneuverability of real carp. The resolution of the high-definition video camera is 1920 × 1080 p, the data rate is 19,949 kb/s, the total bit rate is 20,205 kb/s and frame rate is 29 fps, respectively. Stabilized swimming over a distance of 60 cm is selected for testing. Snapshots are also captured and examined by using the Kinovea 8.20 software. The sampling time is set to 0.05 s.

The main idea for kinematic modeling of the robotic fish is to assume the robotic fish to have a multi-link mechanism. The equivalent model of the robotic fish and description of the one link are given in Figure 1.

The designed robotic fish consists of five links. We will define $[l_1,\dots ,l_i]\to (i=1,2,\dots ,N)$ as one body link, three active tail links and one passive oscillating caudal fin link. $N=5$ is the number of links. ${\theta }_1$ is the angle of head during swimming, ${\theta }_k(k=2,3,4)$ is the each active link angle of the propulsion mechanism and ${\theta }_5$ is angle of the caudal fin. $m_1$ is the body mass and each link mass of the propulsion mechanism equals from $m_2$ to $m_5$. $(x_i^g,y_i^g)$ and $(x_i,y_i)$ are link coordinates of the center of mass and end points of the link coordinates, respectively. Also, ${\phi }_i={\theta }_i+{\theta }_{i-1}$ is the orientation angle of the joints in the horizontal plane. It can be seen from Figure 2b that each link is homogeneous and the middle points of the links are regarded as the center of mass. $a_i$ is the distance from joint to center of the link mass and ${\tau }_i$ is the turning moment of the each joint.

Motion of a real fish depends on the body-traveling wave. In order to imitate the body-traveling wave, there are two general methods: sine-based joint location or intersection methods [19,32], and CPGs [21,33]. In this paper, the improved intersection method with the Big Bang–Big Crunch optimization algorithm, which is also used in our previous work [32], is applied to the link points of the real carp. In this method, locations of links are determined within a scanning area to achieve the minimum error of the fitting body-travelling wave. Thus, optimization criteria is chosen as error of the envelop area with optimum link lengths [32,34]. In the analysis of the carp, five points are defined to verify optimized link lengths on forward and turning patterns obtained from the top view of the carp. The locations of the points and optimized link lengths on forward and the turning patterns are shown in Figure 2.

In the experiments, five points are marked on the real fish in order to define link locations according to optimized link lengths of the robotic fish. Locations of the points are tested by applying links between every two points to ensure a minimum area.

As seen from Figure 2, light blue line is the anterior body, orange, pink and deep blue lines are active links and the green line is the passive caudal fin, respectively.

2.2. Swimming Motion of the Robotic Fish

Motions of the carp obtained from experimental analysis according to all joints for forward and turning swimming modes are given in Figure 3.

The average speed of a carp ranges from 0 to 2.5 m/s in the ichthyology. In the experimental studies, speeds of forward swimming are recorded for some different ranges from 0.105 m/s, 0.125 m/s, 0.144 m/s, 0.233 m/s to 0.41 m/s and turning speed is obtained as 0.24 m/s, approximately. In order to determine flapping frequency of the carp, motion of the caudal fin is analyzed as shown in Figure 4 during different periods. The average recorded flapping frequencies of the carp are measured as 1.818 Hz, 2 Hz, 2.439 Hz and 2.857 Hz, approximately. In order to calculate one sample flapping frequency for 2.857 Hz, minimum and maximum angles of the caudal fin from the center of gravity are measured and time of the fin motion is determined. In this sample, motion time for one period is measured as 350 ms.

Angles of the active links are also shown in Figure 5 for forward and turning swimming modes during one period.

Motion of the robotic fish is realized by the powerful spine that generates the lateral sinusoidal undulations, while the anterior body is relatively rigid and determines the swimming direction. Active tail links are actuated by servomotors and the passive caudal fin is connected to the fourth link with a peduncle. It is noted that the thrust force is generated by the caudal fin.

In order to obtain real fish-like motions, the dynamic model of the robotic fish is derived by using following second order form:

$$M({\theta }_i)\ddot{\theta }+C({\theta }_i,{\dot{\theta }}_i){\dot{\theta }}_i+B({\dot{\theta }}_i)+K({\theta }_i)={\tau }_i$$

(1)

$M({\theta }_i)$ is the moment of inertia, $C({\theta }_i,{\dot{\theta }}_i)$ is the Coriolis matrix, $B({\dot{\theta }}_i)$ is the damping coefficient matrix, $K({\theta }_i)$ is the spring coefficient matrix and ${\tau }_i$ is the force vector that includes hydrodynamic forces.

2.3. Hydrodynamic Forces Acting on the Tail

As shown in Figure 6, hydrodynamic forces depend on the tail motion and there are five main hydrodynamic forces acting on the tail. Inertial fluid force is $F_V$, lift force is $F_J$, thrust force is $F_F$, lateral force is $F_C$ and drag force is $F_D$.

It is considered that tail fin of the robotic fish flaps in constant flow representing by $U_m$ so that the inertial and lift forces influence the caudal fin. Afterwards, thrust force and lateral force components are calculated by using inertial fluid and lift forces. The value of constant flow is determined as 0.08 m/s to reduce the turbulent effects [35]. According to Figure 6, $F_V$ is a proportional force to the acceleration and it is calculated by Equation (2) and the lift force is also vertical to the flow and it is defined as Equation (3):

$$F_V=\pi \rho L{C}^2\dot{U}\sin (\alpha )+\pi \rho L{C}^{2}U\dot{\alpha }\cos (\alpha )$$

(2)

$$F_J=2\pi \rho LC{U}^2\sin (\alpha )\cos (\alpha )$$

(3)

here, the chord length is 2C, the span of the caudal fin is L, density of the water is $\rho$, relative speed at the center of caudal fin is $U$ and the angle of attack is α. Relative speed and angle of attack variables are time dependent and the first derivative of these should be obtained, smoothly. These forces are separated into the thrust force on the x-axis and lateral force on the y-axis components:

$$F_F=(F_V-F_J)\sin (2\pi -{\displaystyle \sum _{i=1}^N{\theta }_i})$$

(4)

$$F_C=(F_V-F_J)\cos (2\pi -{\displaystyle \sum _{i=1}^N{\theta }_i})$$

(5)

If motion of the robotic fish is considered in the x-axis direction, the relative speed on the center of caudal fin in the y-axis direction is expressed by the following equation:

$$\begin{array}{cc}\hfill u_f=& l_1{\dot{\theta }}_1\cos ({\theta }_1)+l_2({\dot{\theta }}_1+{\dot{\theta }}_2)\cos ({\theta }_1+{\theta }_2)\hfill \\ & +l_3({\dot{\theta }}_1+{\dot{\theta }}_2+{\dot{\theta }}_3)\cos ({\theta }_1+{\theta }_2+{\theta }_3)\hfill \\ & +l_4({\dot{\theta }}_1+{\dot{\theta }}_2+{\dot{\theta }}_3+{\dot{\theta }}_4)\cos ({\theta }_1+{\theta }_2+{\theta }_3+{\theta }_4)\hfill \\ & +a_5({\dot{\theta }}_1+{\dot{\theta }}_2+{\dot{\theta }}_3+{\dot{\theta }}_4+{\dot{\theta }}_5)\cos ({\theta }_1+{\theta }_2+{\theta }_3+{\theta }_4+{\theta }_5)\hfill \end{array}$$

(6)

Relative speed on the center of caudal fin $U$ is determined by $\sqrt{{U_m}^2+{u_f}^2}$, because $U_m$ and $u_f$ are perpendicular to each other. The thrust force $F_F$ pushes the robotic fish in the forward direction. The drag force $F_D$ acts on the fish body in the flow direction by the friction. The definition of the drag force is presented in Equation (7):

$$F_D=\frac{1}{2}\rho C_D{V}_f^{2}S$$

(7)

where, $C_D$ is the drag coefficient, $V_f$ is the relative speed of the robotic fish and S is the area of the main body. The head of the robotic fish is assumed as conic and the value of drag coefficient is defined as 0.5. Both hydrodynamic and external forces act on the robotic fish tail. These forces are called with ${\tau }_i={[{\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4,{\tau }_5]}^T$ which affect to the ith link, and can be expressed by:

$${\mathrm{if}\{\begin{array}{c}active\_link\to {\tau }_i=F_C{l}_i\cos ({\displaystyle \sum _{i=1}^N{\theta }_i})-F_F{l}_i\sin ({\displaystyle \sum _{i=1}^N{\theta }_i})+T_{i-1}^u\\ passive\_link\to {\tau }_i=F_C{l}_i\cos ({\displaystyle \sum _{i=1}^N{\theta }_i})-F_F{l}_i\sin ({\displaystyle \sum _{i=1}^N{\theta }_i})\end{array}|}_{last\_link\to l_i=a_i}$$

(8)

The designed robotic fish is driven by three active joints that generate the thrust force. These joints are actuated by three DC servo motors with input torques $T_i^u=A_{\max }\sin (2\pi f_{T}t\pm {\Delta }_i)$. $A_{\max }$ is the amplitude of input torque, $f_T$ is the frequency of input torque and $\pm {\Delta }_i$ is the phase angle between the two input torques. By using Lagrange’s equation, the dynamic model of the robotic fish is described briefly in Equation (9):

$$\left[\begin{array}{ccccc}{M}_{11}& M_{12}& M_{13}& M_{14}& M_{15}\\ M_{21}& M_{22}& M_{23}& M_{24}& M_{25}\\ M_{31}& M_{32}& M_{33}& M_{34}& M_{35}\\ M_{41}& M_{42}& M_{43}& M_{44}& M_{45}\\ M_{51}& M_{52}& M_{53}& M_{54}& M_{55}\end{array}\right]\left[\begin{array}{c}{\ddot{\theta }}_1\\ {\ddot{\theta }}_2\\ {\ddot{\theta }}_3\\ {\ddot{\theta }}_4\\ {\ddot{\theta }}_5\end{array}\right]=\left[\begin{array}{c}{N}_1\\ N_2\\ N_3\\ N_4\\ N_5\end{array}\right]$$

(9)

$M_{ij}$ is the moment element of the inertia matrix and each value of ${\theta }_i$ can be obtained from Equation (9).

2.4. Modeling of the Up-Down Motion Mechanism

In nature, up-down motions such as ascending and descending are another important criterion to exhibit high swimming performance. In order to generate ascending and descending motions, a fish changes its center of mass (COM) position and uses the air bladder. Figure 7 gives the descending pattern of the carp during 1.05 s. In this analysis, pitch angle of the carp is measured as 35°. This performance is the high motion ability for a real fish and similar values of this angle are validated by using the derived model.

In Figure 8a,b, the model of the up-down motion mechanism and flow chart, which represents the working principle of the mechanism, are given, respectively. An up-down mechanism is designed and modeled by using the Lagrange approach to generate and control pitch and roll motions. This mechanism is placed inside the anterior robotic fish body. While the robotic fish is swimming under water, swimming speed is not zero and the depth position of the robotic fish in the vertical plane can be controlled. Control of the orientation angles are performed by changing COM position. There are two masses $(m_x,m_y)$ moving on the shaft. Shafts are connected to the pulley sets with DC motors. By moving the masses along the x and y-axis, the up-down mechanism changes the pitch and roll angles.

When mass positions change, pitch and roll torques are generated. These mentioned torques are considered in Equations (10) and (11), respectively.

$${\tau }_x=g((l_x\sin (\varphi )\sin (\theta )-h_x\sin (\varphi )\cos (\theta ))m_x+(l_y\cos (\varphi )-h_y\sin (\varphi )\cos (\theta ))m_y)$$

(10)

$${\tau }_y=g((l_x\cos (\theta )-h_x\sin (\theta ))m_x+(h_y\sin (\theta ))m_y)$$

(11)

Here, g is the acceleration of the gravity, $(l_x,l_y)$ are the distance from the origin and $(h_x,h_y)$ are the heights from the axis.

2.5. Three Dimensional Model Equations of Motion

Based on the ichthyology of carp, a multi-link and autonomous-swimming biomimetic robotic fish prototype is designed. The prototype consists of a rigid body, three active tail links and a passive caudal fin. In order to analyze the guidance states of the robotic fish, 6-DoF body motion equations are presented in this paper. These equalities define the kinematic and dynamic behaviors of the robotic fish in two coordinate frames and also include the hydrodynamic forces and moments acting on the robotic fish body. The moving coordinate frame is fixed to the body and it is expressed as a body-fixed reference frame. The origin of the body-fixed frame is pointed with the center of gravity. Also, the motion of the body-fixed frame depends on the inertial reference frame that is named the earth-fixed frame. The 6-DoF motion of the robotic fish prototype in the Earth-Fixed frame with translation and rotation relations is shown in Figure 9.

The generalized speed and positions of the robotic fish are completely given by the state variable vector $s={\left[\begin{array}{cccccccccccc}u& v& w& p& q& r& x& y& z& \varphi & \theta & \psi \end{array}\right]}^T$. The linear speeds $(u,v,w)$ are surge, sway and heave, respectively. The angular speeds $(p,q,r)$ are roll, pitch and yaw speeds, respectively. Surge, sway and heave positions $(x,y,z)$ are determined with the earth-fixed frame. Roll, pitch and yaw angles $(\varphi ,\theta ,\psi )$ are the rigid body orientations with respect to the earth-fixed frame.

The kinetic translation between the body-fixed and the earth-fixed frames for the robotic fish motion is expressed as:

$$\dot{\eta }=\left[\begin{array}{cc}{R}_{\Theta }(\eta )& 0_{3\times 3}\\ 0_{3\times 3}& T_{\Theta }(\eta )\end{array}\right]v$$

(12)

where, $\eta ={\left[\begin{array}{cccccc}x& y& z& \varphi & \theta & \psi \end{array}\right]}^T$ and $v={\left[\begin{array}{cccccc}u& v& w& p& q& r\end{array}\right]}^T$ are the generalized positions and speeds, respectively. $\eta$ is relative with the body-fixed frame and speed vector depends on the earth-fixed frame. $R_{\Theta }(\eta )$ is the rotation matrix and $T_{\Theta }(\eta )$ is the transformation matrix. In this way, the 3D equation of the rigid body motion is given by:

$$M\dot{v}+C(v)v+D(v)v+g(\eta )={\tau }_f$$

(13)

where, $M=M_{RB}+M_A$ is the mass matrix, $M_{RB}$ is the rigid body mass matrix and $M_A$ is the added mass matrix. $C(v)=C_{RB}(v)+C_A(v)$ is the Coriolis and centripetal matrix, $C_{RB}(v)$ is the rigid body Coriolis matrix and $C_A(v)$ is the hydrodynamic Coriolis matrix. $D(v)$ is the hydrodynamic damping matrix including linear and non-linear damping effects. ${\tau }_f$ is the vector that includes the hydrodynamic forces. Rigid body mass matrix and added mass matrix are defined as:

$$M_{RB}=\left[\begin{array}{cc}m{I}_{3\times 3}& -mS(r_G)\\ mS(r_G)& I_0\end{array}\right]$$

(14)

$$M_A=\left[\begin{array}{cccccc}{X}_{\dot{u}}& 0& 0& 0& 0& 0\\ 0& Y_{\dot{v}}& 0& 0& 0& Y_{\dot{r}}\\ 0& 0& Z_{\dot{w}}& 0& Z_{\dot{q}}& 0\\ 0& 0& 0& K_{\dot{p}}& 0& 0\\ 0& 0& M_{\dot{w}}& 0& M_{\dot{q}}& 0\\ 0& N_{\dot{v}}& 0& 0& 0& N_{\dot{r}}\end{array}\right]$$

(15)

In these equations, m is the total body mass. The position of the center of gravity is equal to $r_G={\left[\begin{array}{ccc}{x}_G& y_G& z_G\end{array}\right]}^T$. The moments of inertia are defined as $I_0$. Asymmetric matrix S is used for easier computation. $X_{\dot{u}}$ and so forth are the zero-frequency added mass coefficients. $C_{RB}(v)$ matrix indicates the Coriolis matrix ${\left[\begin{array}{ccc}p& q& r\end{array}\right]}^T\times {\left[\begin{array}{ccc}u& v& w\end{array}\right]}^T$ and it is expressed by:

$$C_{RB}(v)=\left[\begin{array}{cc}{0}_{3\times 3}& -S\left(v_{1}m{I}_{3\times 3}-v_{2}mS(r_G)\right)\\ -S\left(v_{1}m{I}_{3\times 3}+v_2-mS(r_G)\right)& -S\left(v_{1}mS(r_G)+v_2{I}_0\right)\end{array}\right]$$

(16)

$C_A(v)$ matrix contains Munk moments and it is given by:

$$C_A(v)=\left[\begin{array}{cccccc}0& 0& 0& 0& -Z_{\dot{w}}w& Y_{\dot{v}}v\\ 0& 0& 0& Z_{\dot{w}}w& 0& -X_{\dot{u}}u\\ 0& 0& 0& -Y_{\dot{v}}v& X_{\dot{u}}u& 0\\ 0& -Z_{\dot{w}}w& Y_{\dot{v}}v& 0& -N_{\dot{r}}r& M_{\dot{q}}q\\ Z_{\dot{w}}w& 0& -X_{\dot{u}}u& N_{\dot{r}}r& 0& -K_{\dot{p}}p\\ -Y_{\dot{v}}v& X_{\dot{u}}u& 0& -M_{\dot{q}}q& K_{\dot{p}}p& 0\end{array}\right]$$

(17)

Hydrodynamic damping effects of the robotic fish are separated into linear and non-linear terms. These matrixes can be defined as below equations:

$$D_L=\left[\begin{array}{cccccc}{X}_u& 0& 0& 0& 0& 0\\ 0& Y_v& 0& 0& 0& -Y_r\\ 0& 0& Z_w& 0& Z_q& 0\\ 0& 0& 0& K_p& 0& 0\\ 0& 0& M_w& 0& M_q& 0\\ 0& -N_v& 0& 0& 0& N_r\end{array}\right]$$

(18)

$$D_{NL}=\left[\begin{array}{cccccc}{X}_{|u|u}\left|u\right|& 0& 0& 0& 0& 0\\ 0& Y_{|v|v}\left|v\right|& 0& 0& 0& -Y_{|v|r}\left|v\right|\\ 0& 0& Z_{|w|w}\left|w\right|& 0& Y_{|w|q}\left|w\right|& 0\\ 0& 0& 0& K_{|p|p}\left|p\right|& 0& 0\\ 0& 0& M_{|q|w}\left|q\right|& 0& M_{|q|q}\left|q\right|& 0\\ 0& -N_{|r|v}\left|r\right|& 0& 0& 0& N_{|r|r}\left|r\right|\end{array}\right]$$

(19)

$D_L$ and $D_{NL}$ are the diagonal matrixes and these matrixes depend on the body shape and speed of the robotic fish, respectively. $X_u$, $X_{|u|u}$ and so forth terms are linear and non-linear positive damping coefficients. In marine applications, underwater vehicles move at low speed and adding mass terms can be neglected. However, these effects are used to obtain complete motion model.

$W=mg$ is the gravity force and $B=\rho g\Delta$ is the buoyancy force. Here, $\Delta$ is the submerged volume. Thus, hydrostatic effect $g(\eta )$ is calculated as:

$$g(\eta )=\left[\begin{array}{c}(W-B)s(\theta )\\ -(W-B)c(\theta )s(\varphi )\\ -(W-B)c(\theta )c(\varphi )\\ -(y_{g}W-y_{b}B)c(\theta )c(\varphi )+(z_{g}W-z_{b}B)c(\theta )s(\varphi )\\ -(z_{g}W-z_{b}B)s(\theta )+(x_{g}W-x_{b}B)c(\theta )c(\varphi )\\ -(x_{g}W-x_{b}B)c(\theta )s(\varphi )-(y_{g}W-y_{b}B)s(\theta )\end{array}\right]$$

(20)

Note that (s(), c(), t()) are equal to (sin(), cos(), tan()). $r_B={\left[\begin{array}{ccc}{x}_B& y_B& z_B\end{array}\right]}^T$ is the position vector from the body origin to center of buoyancy. The force vector ${\tau }_f$ includes the hydrodynamic forces and moments acting on the fish body and it is expressed as:

$${\tau }_f={\left[\begin{array}{ccc}\begin{array}{cccc}{F}_x& F_y& F_z& {\tau }_x\end{array}& {\tau }_y& {\tau }_z\end{array}\right]}^T$$

(21)

where, $(F_x,F_y,F_z)$ represent the surge, sway and heave forces. $({\tau }_x,{\tau }_y,{\tau }_z)$ are the yaw, pitch and roll moments, respectively. It is considered that $F_x=F_F\cos ({\theta }_2+{\theta }_3+{\theta }_4)$. $F_y$ and $F_z$ are chosen as zero, and ${\tau }_z=\left((l_1+l_2)\cos ({\theta }_2)+l_3\cos ({\theta }_2+{\theta }_3)+l_4\cos ({\theta }_2+{\theta }_3+{\theta }_4)\right)F_F\sin ({\theta }_2+{\theta }_3+\theta )$. Finally, 6-DoF motion block diagram is given in Figure 10 that includes the fish body and fish tail equations both body-fixed frame and earth-fixed frame.

By solving these equations, whole 3D non-linear dynamic model of the robotic fish is obtained.

3. Implementations of the Fish-Like Motion

In order to mimic real fish swimming abilities, the robotic fish model focuses on five special points according to the real carp:

• Proportional optimum link lengths according to actual size,
• Flapping frequency,
• Swimming speed performance,
• Proportional physical parameters according to the carp,
• Trajectory tracking.

As can be observed from the experiments, proportional link lengths generate sinusoidal angles while the robotic fish propels itself. This situation causes the sinusoidal undulation motions. Secondly, obtained flapping frequency from the observations is an appropriate value and it is applied to the robotic fish model. Forward and turning speeds of the carp are also examined in the simulations. Finally, guidance trajectory tracking is performed to validate the dynamic model.

3.1. Ability of the Fish-Like Motion

In this study, numerical analyses are performed in a MATLAB environment. The flapping frequency of the carp is observed from the experiments and it is determined as 2.857 Hz. This value is also applied to the robotic fish model to validate the ability of the fish-like motion. In Figure 11, free-swimming motion of the robotic fish is illustrated during an 80 s simulation time.

Figure 11 shows open loop free-swimming behavior on a two-dimensional (2D) reference frame. The start position of the robotic fish is X = 0 m, Y = 0 m. After 80 s simulation time, the robotic fish reaches the final position at X = 0.23 m, Y = 2.40 m. The speed of the robotic fish is measured as 0.184 m/s and 0.474 m/s for forward and turning swimming motions, respectively. Figure 12 shows the link angles and Figure 13 presents the thrust force on the x-axis.

In Figure 12, tail link behaviors are illustrated for different swimming motions. The robotic fish tail swings itself with sinusoidal undulations with pure sine waves.

The thrust force of the tail is determined as nearly 0.110 N and 0.177 N for forward and turning swimming modes.

3.2. Guidance and Trajectory Tracking

The proposed model of the robotic fish including guidance and trajectory tracking consists of a three-layered hierarchical architecture: inner loop, outer loop and trajectory generator. The inner loop is controlled with the simple and practical control structure inspired by Reaching Law Control (RLC) approach. The inner loop occurs in the body frame of the robotic fish. The outer loop deals with the guidance control structure. On the other hand, the trajectory generator builds up the 3D path tracking in earth-fixed frame. The block diagram of the closed loop system is shown in Figure 14.

In Figure 14, $T$, $\eta$, $G$ and $F$ indicate reference trajectory, current positions, reference of the guidance, and current orientations and speed, respectively. These parameters are defined as below:

$$\{\begin{array}{l}T={\left[\begin{array}{ccc}{X}_{ref}& Y_{ref}& Z_{ref}\end{array}\right]}^T\\ {\eta }^{\prime }={[\begin{array}{ccc}x& y& z\end{array}]}^T\\ {\eta }^{″}={[\begin{array}{ccc}\varphi & \theta & \psi \end{array}]}^T\\ G={\left[\begin{array}{ccc}{\theta }_{ref}& {\psi }_{ref}& v_{ref}\end{array}\right]}^T\\ F={\left[\begin{array}{ccc}\theta & \psi & v\end{array}\right]}^T\end{array}$$

(22)

It is noted that the reference matrix of the guidance is calculated by using a trigonometrical approach and that each element of the matrix can be given by Equation (23):

$$\{\begin{array}{l}{\theta }_{ref}=-{\tan }^{-1}\left(\left(z_{ref}-z\right)/\left(\sqrt{{\left(x_{ref}-x\right)}^2+{\left(y_{ref}-y\right)}^2}\right)\right)\\ {\psi }_{ref}={\tan }^{-1}\left(\left(y_{ref}-y\right)/\left(x_{ref}-x\right)\right)\\ v_{ref}=\zeta \sqrt{{\left(x_{ref}-x\right)}^2+{\left(y_{ref}-y\right)}^2+{\left(z_{ref}-z\right)}^2}\end{array}$$

(23)

here, $\zeta$ is a constant speed control gain. In the general RLC approach, the discrete time reaching law [36,37] can be given by:

$$\Delta S(k+1)=-\sigma \mathrm{sgn}\left(S(k)\right)-\alpha S(k)$$

(24)

where $\sigma$ and $\alpha$ are constant parameters. $\Delta S(k+1)$ is then expressed by:

$$\Delta S(k+1)=\frac{S(k+1)-S(k)}{T_s}$$

(25)

$T_s$ is the sampling time. With the assumption of $\{\Delta S(0)=0\}$, discrete time switching function is given as:

$$S(k)=\lambda e(k)+\Delta e(k)$$

(26)

here, $e(k)$ is the error of the guidance and $\Delta e(k)$ is the derivative of the error in the kth sampling step. $\lambda$ is the slope of the switching line and indicates a first-order dynamics of the error. The high-frequency oscillation problem known as chattering is occurred if Equation (24) is used directly in practice [37,38]. This chattering problem is presented in Figure 15 and it can be solved by using $sat\left(S(k)\right)$ function instead of $\mathrm{sgn}\left(S(k)\right)$. $sat(·)$ function is expressed by:

$$sat(\Omega )=\{\begin{array}{c}\Omega \\ \mathrm{sgn}(\Omega )\end{array}\begin{array}{c}\mathrm{if}|\Omega |\le 1\\ \mathrm{if}|\Omega |>1\end{array}$$

(27)

The new reaching law can be obtained as:

$$\Delta S(k+1)=-\sigma \mathrm{s}at(S(k))-\alpha S(k)$$

(28)

The reaching law is also rewritten as below with this $sat(·)$ function in the Boundary Layer (BL):

$$\Delta S(k+1)=-\beta S(k)$$

(29)

where $\beta =(\sigma /BL)+\alpha$. Thus, the control output becomes:

$$\Delta \upsilon (k)=KS(k)$$

(30)

Since $S$ is forced to reach zero with this RLC approach, exponentially reducing of the error to zero with time constant of $(1/\lambda )$ is provided.

Also, an anti-windup integrator is added to stop over-integration for the protection of the designed model as seen in Figure 16.

6-DoF motion responses completely include the generalized linear and angular positions and speeds of the robotic fish. Positions (X, Y, Z) are determined with the earth-fixed frame. Roll, pitch and yaw angles $(\varphi ,\theta ,\psi )$ are also rigid body orientations with respect to the earth-fixed frame. In order to validate the accuracy of the designed model, simulations are examined in this section.

Figure 17a shows the closed loop forward speed responses of the robotic fish from 0.1 m/s to 0.4 m/s. These speed values are chosen as observed speed ranges of the carp. Figure 17b shows the orientation control that includes yaw angles of the robotic fish model during 40 s simulation time. The control performance of the closed loop system is summarized in

Table 1. The closed loop control performance: Maximum overshoot (%) and settling time (s) of forward speed for (0.1 m/s)–(0.4 m/s) and yaw angles for (−40°)–(+40°).

	Interval	Maximum Overshoot (%)	Settling Time (s)
Forward Speed (m/s)	0.1	6.80	16.30
	0.2	4.40	11.80
	0.3	2.23	15.40
	0.4	1.65	12.20
Yaw Angles (°)	−40	37.50	18.30
	−30	16.66	8.80
	−20	15	16.65
	−10	10	8.20
	0	36	14.50
	10	27.50	10.80
	20	6.50	21.50
	30	15	17.60
	40	31.25	18.90

.

Another key point is to realize pitch and roll motions together. In Figure 18, these motions are performed for step and stairs references. The closed loop control performance of the roll and pitch motions are illustrated in

Table 2. The simultaneously closed loop control performance of the roll and pitch motions for step and stairs references.

Criteria	Step		Stairs
Value (°)	20	−15	−20		0		20
Motion	Roll	Pitch	Roll	Pitch	Roll	Pitch	Roll	Pitch
Maximum Overshoot (%)	-	-	-	-	27	-	27.50	-
Settling Time (s)	0.38	0.37	0.20	0.24	0.18	0.26	0.19	0.24

.

In order to validate the closed loop system, two different path trajectories are also generated in earth-fixed frame for 3D trajectory tracking. First path is a 3D sinusoidal path (P₁) reference and the second path is a 3D CCW circular path (P₂) reference. Trajectory tracking responses indicate the time-dependent performances of the robotic fish model. Figure 19 shows the closed loop P₁ tracking performance during a 60 s simulation time.

Amplitude and frequency of the P₁ reference are ±0.18 m and 0.2 Hz, respectively. The depth of P₁ is 3 m. The robotic fish tracks the P₁ reference with a speed of 0.42 m/s. Figure 20 presents the P₁ tracking errors of speed and orientation angles.

6-DoF motion responses including positions, yaw, roll, pitch and also link angles are given in Figure 21. The robotic fish swings the flexible tail to keep the P₁ reference during closed loop tracking.

Figure 22 shows the closed loop P₂ tracking performance during a 64 s simulation time. Diameter of the P₂ is 8 m and depth is 1.28 m. Initial positions of the robotic fish and P₂ reference begin the origin at X = 0 m, Y = 0 m, Z = 0 m.

The robotic fish turns around the CCW direction to track the P₂ reference with a speed of 0.4 m/s. Figure 23 also presents the P₂ tracking errors of speed and orientation angles.

6-DoF motion responses are given in Figure 24. When the robotic fish reaches the steady state, link angles become stable and generate small amplitudes.

As shown in

Table 3. The tracking errors for proposed closed loop control system: P₁ and P₂ indicate the generated 3D path trajectories.

Trajectory	ex (m)	ey (m)	ez (m)
P1	0.0968	0.0131	0.1058
P2	0.0827	0.0107	0.0920

, the tracking errors of P₁ and P₂ references on the x, y and z axis are evaluated by Root Mean Square Error (RMSE) performance index and given as:

$$RMSE=\sqrt{\frac{{\displaystyle \sum _{k=1}^{n}e{(k)}^2}}{n}}$$

(31)

In order to visualize the performance of the proposed closed loop control system in a virtual reality platform, the designed robotic fish model is also adopted to Virtual Reality Modeling Language (VRML). A simulation scenario is constructed consisting of a marine environment. The goal is to perform the P₂ tracking underwater. Figure 25 shows the top view above the sea, and the isometric view of the underwater environment is also shown in Figure 26.

All of the simulation results show that the bio-inspired non-linear dynamic model of the robotic fish is suitable to mimic biological features of Carangiform carp. In this way, biomimetic modeling can be easily implemented to the robotic fish prototype in the future work.

4. Conclusions

This study presents a complete 3D non-linear dynamic model of a biomimetic robotic fish based on real carp locomotion. The designed model is derived by using the Lagrange formulation approach and includes both kinematic and hydrodynamic effects acting on each link. In order to achieve six degrees of freedom motions in the earth fixed frame, an up-down mechanism is designed that controls the influence of pitch and roll. This mechanism gives the position and orientation of the robotic fish in 3D space. Moreover, linear and rotational behavior responses are obtained and controlled in the simulations. In order to mimic biological features of fish with the dynamic model, various swimming patterns of a real carp are observed in the experimental pool. Link length ratio, flapping frequency, speed range, displacement and Euler angles are analyzed to determine appropriate desired inputs for the designed model. According to optimized real carp values, link lengths of the anterior body, active links and passive caudal fins are determined as 0.044 m, [0.043 m 0.040 m 0.010 m] and 0.063 m, respectively. Flapping frequency of carp is observed from 1.818 Hz to 2.857 Hz and it is also chosen as 2.857 Hz in the simulations. Average speed of a carp ranges from 0 to 2.5 m/s in the ichthyology. In the experimental studies, the speed range of forward swimming is recorded from 0.105 m/s to 0.41 m/s and turning speed is obtained as 0.24 m/s, approximately. It is also seen that the simulation speed values remained in these ranges. The pitch angle of the descending pattern of the carp is measured as 35° during 1.05 s. According to simulation results, the descending motion of the robotic fish is performed with 0.24 s settling time for 20° pitch angle, which has very similar performance to carp.

Forward and turning speeds, descending-ascending motions and yaw, pitch, roll angles are also controlled to achieve real fish behaviors. Figure 11, Figure 17 and Figure 18 show the maneuverability of the fish-like motions. Similarly, two different reference trajectories (P₁ and P₂) are tracked to evaluate sufficient swimming performances in Figure 19 and Figure 22. The proposed closed loop control system is also adapted to VRML in order to validate the P₂ tracking in virtual marine environment. These simulation and experimental analyses show that the derived model achieves real carp motions to implement robotic fish prototypes for biomimetic design.

In future work, these experiences should be transferred to robotic fish prototypes in experimental platforms and path-tracking performance should be examined with various intelligent control algorithms.