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

Gonca Ozmen Koca 1,*, Cafer Bal 1, Deniz Korkmaz 2 ID , Mustafa Can Bingol 1, Mustafa Ay 1 ID , Zuhtu Hakan Akpolat 1 and Seda Yetkin 3 1 Department of Mechatronics Engineering, University of Firat, Elazig 23119, Turkey; caferbal@gmail.com (C.B.); mustafacanbingol@gmail.com (M.C.B.); mustafaay023@gmail.com (M.A.); z.h.akpolat@gmail.com (Z.H.A.) 2 Department of Electrical and Electronics Engineering, University of Firat, Elazig 23119, Turkey; denizkorkmaz17@gmail.com 3 Department of Electronics and Automation, University of Bitlis Eren, Bitlis 13000, Turkey; syetkin@beu.edu.tr * Correspondence: gonca.ozmen@gmail.com; Tel.: +90-424-237-0000


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].
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.

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. 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 nonlinear 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.

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 1 [ ,..., ] ( 1, 2,..., ) i l l i N   as one body link, three active tail links and one passive oscillating caudal fin link. 5 N  is the number of links. 1  is the angle of head during swimming, is the each active link angle of the propulsion mechanism and 5  is angle of the caudal fin. 1 m is the body mass and each link mass of The designed robotic fish consists of five links. We will define [l 1 , . . . , l i ] → (i = 1, 2, . . . , N) as one body link, three active tail links and one passive oscillating caudal fin link. N = 5 is the number of links. θ 1 is the angle of head during swimming, θ k (k = 2, 3, 4) is the each active link angle of the propulsion mechanism and θ 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 g i , y g i ) and (x i , y i ) are link coordinates of the center of mass and end points of the link coordinates, respectively. Also, ϕ i = θ i + θ 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 τ 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. i i x y are link coordinates of the center of mass and end points of the link coordinates, respectively. Also, 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. i a is the distance from joint to center of the link mass and 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.

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. 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.

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: 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:  Angles of the active links are also shown in Figure 5 for forward and turning swimming modes during one period.
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:  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: is the damping coefficient matrix, K(θ i ) is the spring coefficient matrix and τ i is the force vector that includes hydrodynamic forces.

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 .

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 V lateral force is C F and drag force is D F . It is considered that tail fin of the robotic fish flaps in constant flow representing by m U 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, V F 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); here, the chord length is 2C, the span of the caudal fin is L, density of the water is  , 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: 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;  Figure 6. Hydrodynamic force distributions acting on the tail.
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): here, the chord length is 2C, the span of the caudal fin is L, density of the water is ρ, 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: 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: Relative speed on the center of caudal fin U is determined by 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): 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 T which affect to the ith link, and can be expressed by: 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 u i = A max sin(2π f T t ± ∆ i ). A max is the amplitude of input torque, f T is the frequency of input torque and ±∆ 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): M ij is the moment element of the inertia matrix and each value of θ i can be obtained from Equation (9).

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. Relative speed on the center of caudal fin U is determined by , because m U and f u are perpendicular to each other. The thrust force F F pushes the robotic fish in the forward direction. The drag force D F acts on the fish body in the flow direction by the friction. The definition of the drag force is presented in Equation (7); where, D C is the drag coefficient, f V 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 which affect to the ith link, and can be expressed by: 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 max sin(2 ) is the amplitude of input torque, T f is the frequency of input torque and 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) ij M is the moment element of the inertia matrix and each value of i  can be obtained from Equation (9).

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.
Appl. Sci. 2018, 7, x FOR PEER REVIEW 8 of 22 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 ( , ) x y m m 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.
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.

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.
Appl. Sci. 2018, 7, x FOR PEER REVIEW 9 of 22 (( cos( ) sin( )) ( sin( )) ) Here, g is the acceleration of the gravity, ( , ) x y l l are the distance from the origin and ( , ) x y h h are the heights from the axis.

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 kinetic translation between the body-fixed and the earth-fixed frames for the robotic fish motion is expressed as;  The generalized speed and positions of the robotic fish are completely given by the state variable vector s = u v w p q r x y z φ θ ψ 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 (φ, θ, ψ) 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: where, η = x y z φ θ ψ T and v = u v w p q r T are the generalized positions and speeds, respectively. η is relative with the body-fixed frame and speed vector depends on the earth-fixed frame. R Θ (η) is the rotation matrix and T Θ (η) is the transformation matrix. In this way, the 3D equation of the rigid body motion is given by: 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. τ f is the vector that includes the hydrodynamic forces. Rigid body mass matrix and added mass matrix are defined as: In these equations, m is the total body mass. The position of the center of gravity is equal to The moments of inertia are defined as I 0 . Asymmetric matrix S is used for easier computation. X . u and so forth are the zero-frequency added mass coefficients. C RB (v) matrix indicates the Coriolis matrix p q r T × u v w T and it is expressed by: C A (v) matrix contains Munk moments and it is given by: 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 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 = ρg∆ is the buoyancy force. Here, ∆ is the submerged volume. Thus, hydrostatic effect g(η) is calculated as: Note that (s(), c(), t()) are equal to (sin(), cos(), tan()). r B = x B y B z B T is the position vector from the body origin to center of buoyancy. The force vector τ f includes the hydrodynamic forces and moments acting on the fish body and it is expressed as: where, (F x , F y , F z ) represent the surge, sway and heave forces. (τ x , τ y , τ z ) are the yaw, pitch and roll moments, respectively. It is considered that F x = F F cos(θ 2 + θ 3 + θ 4 ). F y and F z are chosen as zero, and τ z = ((l 1 + l 2 ) cos(θ 2 ) + l 3 cos(θ 2 + θ 3 ) + l 4 cos(θ 2 + θ 3 + θ 4 ))F F sin(θ 2 + θ 3 + θ). 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.
is the position vector from the body origin to center of buoyancy. The force vector f  includes the hydrodynamic forces and moments acting on the fish body and it is expressed as; where, ( , , ) 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.

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: 1. Proportional optimum link lengths according to actual size, 2. Flapping frequency, 3. Swimming speed performance, 4. Proportional physical parameters according to the carp, 5. Trajectory tracking. By solving these equations, whole 3D non-linear dynamic model of the robotic fish is obtained.

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: 1.
Proportional optimum link lengths according to actual size, 2.
Proportional physical parameters according to the carp, 5.
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.

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.
Appl. Sci. 2018, 7, x FOR PEER REVIEW 12 of 22 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.

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 12 shows the link angles and Figure 13 presents the thrust force on the x-axis.   Figure 12 shows the link angles and Figure 13 presents the thrust force on the x-axis.
Appl. Sci. 2018, 7, x FOR PEER REVIEW 12 of 22 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.

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 12 shows the link angles and Figure 13 presents the thrust force on the x-axis. The thrust force of the tail is determined as nearly 0.110 N and 0.177 N for forward and turning swimming modes.

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 ,  , G and F indicate reference trajectory, current positions, reference of the guidance, and current orientations and speed, respectively. These parameters are defined as below: 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) 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.

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. The thrust force of the tail is determined as nearly 0.110 N and 0.177 N for forward and turning swimming modes.

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 ,  , G and F indicate reference trajectory, current positions, reference of the guidance, and current orientations and speed, respectively. These parameters are defined as below: 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); In Figure 14, T, η, G and F indicate reference trajectory, current positions, reference of the guidance, and current orientations and speed, respectively. These parameters are defined as below: 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): here, ζ is a constant speed control gain. In the general RLC approach, the discrete time reaching law [36,37] can be given by: where σ and α are constant parameters. ∆S(k + 1) is then expressed by: T s is the sampling time. With the assumption of {∆S(0) = 0}, discrete time switching function is given as: here, e(k) is the error of the guidance and ∆e(k) is the derivative of the error in the kth sampling step.
λ 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(S(k)) function instead of sgn(S(k)). sat(·) function is expressed by: (27) here,  is a constant speed control gain. In the general RLC approach, the discrete time reaching law [36,37] can be given by; where  and  are constant parameters.
( 1) S k   is then expressed by; s T is the sampling time. With the assumption of   (0) 0 S   , discrete time switching function is given as; here, ( ) e k is the error of the guidance and ( ) e k  is the derivative of the error in the kth sampling step.  is the slope of the switching line and indicates a first-order dynamics of the error. The highfrequency 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 The new reaching law can be obtained as: The reaching law is also rewritten as below with this   sat  function in the Boundary Layer (BL); where Thus, the control output becomes: Since is forced to reach zero with this RLC approach, exponentially reducing of the error to zero with time constant of   1/  is provided. The new reaching law can be obtained as: The reaching law is also rewritten as below with this sat(·) function in the Boundary Layer (BL): where β = (σ/BL) + α. Thus, the control output becomes: Since S is forced to reach zero with this RLC approach, exponentially reducing of the error to zero with time constant of (1/λ) 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. Also, an anti-windup integrator is added to stop over-integration for the protection of the designed model as seen in Figure 16.  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. 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.
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 (P1) reference and the second path is a 3D CCW circular path (P2) reference. Trajectory tracking responses indicate the timedependent performances of the robotic fish model. Figure 19 shows the closed loop P1 tracking performance during a 60 s simulation time.
Amplitude and frequency of the P1 reference are ±0.18 m and 0.2 Hz, respectively. The depth of P1 is 3 m. The robotic fish tracks the P1 reference with a speed of 0.42 m/s. Figure 20 presents the P1 tracking errors of speed and orientation angles.  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. Also, an anti-windup integrator is added to stop over-integration for the protection of the designed model as seen in Figure 16.  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. 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.
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 (P1) reference and the second path is a 3D CCW circular path (P2) reference. Trajectory tracking responses indicate the timedependent performances of the robotic fish model. Figure 19 shows the closed loop P1 tracking performance during a 60 s simulation time.
Amplitude and frequency of the P1 reference are ±0.18 m and 0.2 Hz, respectively. The depth of P1 is 3 m. The robotic fish tracks the P1 reference with a speed of 0.42 m/s. Figure 20 presents the P1 tracking errors of speed and orientation angles. 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.
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 1 ) reference and the second path is a 3D CCW circular path (P 2 ) reference. Trajectory tracking responses indicate the time-dependent performances of the robotic fish model. Figure 19 shows the closed loop P 1 tracking performance during a 60 s simulation time.
Amplitude and frequency of the P 1 reference are ±0.18 m and 0.2 Hz, respectively. The depth of P 1 is 3 m. The robotic fish tracks the P 1 reference with a speed of 0.42 m/s. Figure 20 presents the P 1 tracking errors of speed and orientation angles.   Table 2. The simultaneously closed loop control performance of the roll and pitch motions for step and stairs references.

Criteria
Step    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 P1 reference during closed loop tracking.  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 P1 reference during closed loop tracking. 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 1 reference during closed loop tracking.  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 P1 reference during closed loop tracking.  Figure 22 shows the closed loop P 2 tracking performance during a 64 s simulation time. Diameter of the P 2 is 8 m and depth is 1.28 m. Initial positions of the robotic fish and P 2 reference begin the origin at X = 0 m, Y = 0 m, Z = 0 m. Figure 22 shows the closed loop P2 tracking performance during a 64 s simulation time. Diameter of the P2 is 8 m and depth is 1.28 m. Initial positions of the robotic fish and P2 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 P2 reference with a speed of 0.4 m/s. Figure 23 also presents the P2 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 of P1 and P2 references on the x, y and z axis are evaluated by Root Mean Square Error (RMSE) performance index and given as: The robotic fish turns around the CCW direction to track the P 2 reference with a speed of 0.4 m/s. Figure 23 also presents the P 2 tracking errors of speed and orientation angles. The robotic fish turns around the CCW direction to track the P2 reference with a speed of 0.4 m/s. Figure 23 also presents the P2 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 of P1 and P2 references on the x, y and z axis are evaluated by Root Mean Square Error (RMSE) performance index and given as: 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 of P 1 and P 2 references on the x, y and z axis are evaluated by Root Mean Square Error (RMSE) performance index and given as:  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 P2 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  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 2 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.  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 P2 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.

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. Figures 11, 17 and 18 show the maneuverability of the fishlike motions. Similarly, two different reference trajectories (P1 and P2) are tracked to evaluate sufficient swimming performances in Figures 19 and 22. The proposed closed loop control system is also adapted to VRML in order to validate the P2 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. 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.

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. Figures 11, 17 and 18 show the maneuverability of the fish-like motions. Similarly, two different reference trajectories (P 1 and P 2 ) are tracked to evaluate sufficient swimming performances in Figures 19 and 22. The proposed closed loop control system is also adapted to VRML in order to validate the P 2 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.