Simulation and Controller Design for a Fish Robot with Control Fins

In this paper, a nonlinear simulation block for a fish robot was designed using MATLAB Simulink. The simulation block incorporated added masses, hydrodynamic damping forces, restoring forces, and forces and moments due to dorsal fins, pectoral fins, and caudal fins into six-degree-of-freedom equations of motion. To obtain a linearized model, we used three different nominal surge velocities (i.e., 0.2 m/s, 0.4 m/s, and 0.6 m/s). After obtaining output responses by applying pseudo-random binary signal inputs to a nonlinear model, an identification tool was used to obtain approximated linear models between inputs and outputs. Utilizing the obtained linearized models, two-degree-of-freedom proportional, integral, and derivative controllers were designed, and their characteristics were analyzed. For the 0.4 m/s nominal surge velocity models, the gain margins and phase margins of the surge, pitch, and yaw controllers were infinity and 69 degrees, 26.3 dB and 85 degrees, and infinity and 69 degrees, respectively. The bandwidths of surge, pitch, and yaw control loops were determined to be 2.3 rad/s, 0.17 rad/s, and 2.0 rad/s, respectively. Similar characteristics were observed when controllers designed for linear models were applied to the nonlinear model. When step inputs were applied to the nonlinear model, the maximum overshoot and steady-state errors were very small. It was also found that the nonlinear plant with three different nominal surge velocities could be controlled by a single controller designed for a linear model with a nominal surge velocity of 0.4 m/s. Therefore, controllers designed using linear approximation models are expected to work well with an actual nonlinear model.


Introduction
Autonomous underwater vehicles (AUVs) are designed to mimic the efficient and agile movements of fish.They are attracting attention because they can be effectively used for exploration, monitoring, and surveillance of underwater environments [1].These AUVs integrate advanced propulsion mechanisms such as Body Caudal Fin (BCF) and Median and Pectoral Fin (MPF) systems.The BCF system mimics the natural movements of fish, which propel themselves by flapping their bodies and tails.It is used by more than 85% of fish species [2].
Several studies have been conducted on the dynamics of fish robots.Humphreys [3] derived a detailed six-degree-of-freedom (DOF) equation of motion for surge, sway, heave, roll, pitch, and yaw motions, including moments of inertia and hydrodynamic forces.Nahon [4] proposed body lift and drag forces to model the hydrodynamic characteristics of vehicles.The dynamics of an underwater vehicle, Autolycus, have been well described by Tang [5].Vorticity control for propulsion and maneuvering was introduced by the Draper Laboratory, which allows fish to maintain stable swimming [6].Open Fish utilizes a simple and effective wire-driven active and passive compliant body to mimic highly efficient thunniform swimming [7,8].Tuna Bot aims to achieve high-performance fish and overcome the gap between robotic systems and fish swimming capabilities [9].
In UC-Ika 1 [10], inspired by the swimming motion of tuna fish, undulatory motions by the tail peduncle and caudal fin were used to generate propulsion force.This benefits from a tail mechanism that plays a crucial role in the dynamic behavior of the robot.A framework was developed to compute the steady-turning motion of a robotic fish undergoing periodic body or tail deformation [11].Flying fish performing gliding flights with tail-beating motion are being developed for future applications [12,13].Full-body-length swimming motion coordinates for anterior, midbody, and posterior displacements were proposed in iSplash-I to reduce large kinematic errors in existing free-swimming robotic fish [14].In iSplash-II, a new fabric technique was introduced to achieve high-speed propulsion at high frequencies [15].
In previous studies [16,17], pectoral fins were used for roll and pitch controls.Fish pectoral fins can produce high propulsive performance by driving active and passive fin deformation.A systematic approach was suggested to control pectoral fins by naturally accommodating fin constraints and automatically generating "intelligent" behavior (such as fin backing-up when required) for quick maneuvering [18].
Euler's angle speed control and 3D performance of a fish robot, including closed-loop control system responses, can be found in [19].A two-DOF barycenter mechanism was proposed to provide body stabilization and serve as an actuating device for active control design [20].Zeng et al. [21] proposed an underwater robot based on the hybrid propulsion of a quadrotor and undulating fin, which combines high-efficiency propulsion using a bionic undulating fin and stable control via the propeller of an underwater robot.The hybrid system enables independent heave motion, surge motion, in situ steering, and stable hovering.Kim et al. [22] proposed an integral sliding mode controller (ISMC) to stabilize an autonomous underwater vehicle (AUV) with modeling errors.
Aruna et al. [23] analyzed trajectory tracking control methods using a conventional proportional-integral-derivative (PID) control, H∞ control, and a feedforward (along with feedback) control.Xiang et al. [24] proposed a simple PID controller to drive individual AUVs following a 3D path by exploiting the 3D guidance law in the underactuated mode.Yang et al. [25] suggested the nonlinear formation keeping and mooring control of multiple autonomous underwater vehicles using the backstepping method.Qiao et al. [26] developed an effective control method to improve the trajectory tracking control in an underwater vehicle using an adaptive fast nonsingular integral terminal sliding mode control.Oktafianto et al. [27] presented a nonlinear model of an autonomous underwater vehicle (AUV) with six degrees of freedom, which was linearized using the Jacobian matrix.
In this paper, feedback control loops for a fish robot were developed using MAT-LAB Simulink [28].A nonlinear model was obtained by incorporating hydrodynamic effects (damping coefficients and added masses) into conventional six-DOF free body dynamic equations.
To design a linear feedback controller, a linear model was derived from the nonlinear model using system identification tools [29].Using the identified linear model and PID tuning, we obtained a robust PID controller with large gain and phase margins to overcome nonlinear uncertainties.The frequency responses of the linear model with a feedback controller were analyzed to check the gain margin and phase margin.After obtaining the linear model controller, we applied it to both a linearized model and a nonlinear model to compare its performances in both ideal and realistic situations.In the following sections, we will present the complete dynamic equation, the MATLAB Simulink model of the nonlinear model, linear models for surge, pitch, and yaw motions obtained from nonlinear model responses, closed-loop PID controllers for surge, pitch, and yaw motions using linear models, and the final simulation results of the controllers using nonlinear models.

Equation of Motion and Dynamics of Fins
In this section, a complete six-DOF equation of motion was derived to describe the motion of an underwater vehicle with pectoral, dorsal, and caudal fins.

Equation of Motion
The vehicle can move freely with six-DOF motions: surge, sway, heave, roll, pitch, and yaw motions.Figure 1 illustrates the body coordinate system of a fish robot.The origin of the body frame is at the geometrical center of the fish robot's body.The prototype structure is a modified version of Open Fish [7,8].The center of gravity (CG) is located below the origin.The center of buoyancy (CB) coincides with the origin of the body coordinate.In our fish robot, the buoyancy force (B) is equal to the weight (W) to maintain a neutral buoyant force condition.

Equation of Motion and Dynamics of Fins
In this section, a complete six-DOF equation of motion was derived to describe the motion of an underwater vehicle with pectoral, dorsal, and caudal fins.

Equation of Motion
The vehicle can move freely with six-DOF motions: surge, sway, heave, roll, pitch, and yaw motions.Figure 1 illustrates the body coordinate system of a fish robot.The origin of the body frame is at the geometrical center of the fish robot's body.The prototype structure is a modified version of Open Fish [7,8].The center of gravity (CG) is located below the origin.The center of buoyancy (CB) coincides with the origin of the body coordinate.In our fish robot, the buoyancy force (B) is equal to the weight (W) to maintain a neutral buoyant force condition.Assuming the  plane as a plane of symmetry, six-DOF equations of motion are given by Equations ( 1)-( 6).The right-hand side of each equation represents the sum of external forces and moments arising from added masses, hydrodynamic damping, buoyancy, weight, and control.The left-hand side represents rigid body dynamics.The control forces and moments generated by the control fins are described in Section 2.2.For a complete derivation, refer to [5].Assuming the xz plane as a plane of symmetry, six-DOF equations of motion are given by Equations ( 1)-( 6).The right-hand side of each equation represents the sum of external forces and moments arising from added masses, hydrodynamic damping, buoyancy, weight, and control.The left-hand side represents rigid body dynamics.The control forces and moments generated by the control fins are described in Section 2.2.For a complete derivation, refer to [5].
Biomimetics 2024, 9, 317 4 of 16 Variables u, v, w, p, q, r, ϕ, and θ indicate surge velocity, sway velocity, heave velocity, roll rate, pitch rate, yaw rate, roll angle, and pitch angle, respectively.Body mass and the moment of inertia along the x, y, and z axes are m, I xx , I yy , and I zz , respectively.The positions of the center of buoyancy and center of gravity are (x B , y B , z B ) and (x G , y G , z G ), respectively.
The added masses corresponding to surge, sway, heave, roll, pitch, and yaw motions are denoted as r , and Z .q , respectively.The ellipsoidal body is assumed to estimate the added mass and moments of inertia of the vehicle [30,31].The effects of the control fins were neglected.Assuming a symmetric body, we have q .The damping forces are represented by X u|u| u|u|, Y v|v| v|v|, Z w|w| w|w|, K p|p| p|p|, M w|w| w|w|, M q|q| q|q|, N r|r| r|r|, and N v|v| v|v|, illustrating the nonlinear and coupled damping characteristics of an underwater vehicle in six-DOF motions at a high speed .The forces and moments due to the control fins (i.e., X c , Y c , Z c , K c , M c , and N c ) will be detailed in Section 2.2.
The relationship between the Euler angle and the body axis rate is given as Equations ( 7)-( 9) [5].The yaw angle is denoted as ψ.
Tables 1 and 2 show parameters used in the control simulations.They were estimated based on the size and shape of the body.The damping coefficients were approximated assuming a flat ellipsoidal body.To be more realistic, a more accurate calculation of these parameters and verification via real experiments are necessary.
Table 2. Hydrodynamic parameters of the fish robot.

Forces and Moments Due to Control Fins
Caudal fins were utilized to generate thrust, whereas pectoral and dorsal fins were used to control the pitch and yaw motions of the fish robot.The caudal fin produced the propulsion force necessary for propelling the fish forward.We set the tail beat angle to 50 degrees and estimated the average thrust produced by the tail when the flapping frequency varied from 1 Hz to 7 Hz.The surge force (X c ) can be expressed using Equation (10).The terms T cd , D b , D p f , and D d f represent thrust by the caudal fin, body drag, pectoral fin drag, and dorsal fin drag, respectively.
Figure 2 shows thrusts generated by the caudal fin.These thrusts were calculated using the added mass method [12,34,35], first proposed by Lighthill [35] and then modified for calculating the thrust generated by a forked-shaped caudal fin [12].The fin used in [12] was replaced by the current fin geometry.The driving frequency varied from 1 to 7 Hz for the thrust calculation, as shown in Figure 2.
Biomimetics 2024, 9, x FOR PEER REVIEW 5 of 17 Table 2. Hydrodynamic parameters of the fish robot.

Forces and Moments Due to Control Fins
Caudal fins were utilized to generate thrust, whereas pectoral and dorsal fins were used to control the pitch and yaw motions of the fish robot.The caudal fin produced the propulsion force necessary for propelling the fish forward.We set the tail beat angle to 50 degrees and estimated the average thrust produced by the tail when the flapping frequency varied from 1 Hz to 7 Hz.The surge force ( ) can be expressed using Equation (10).The terms  ,  ,  , and  represent thrust by the caudal fin, body drag, pectoral fin drag, and dorsal fin drag, respectively.
Figure 2 shows thrusts generated by the caudal fin.These thrusts were calculated using the added mass method [12,34,35], first proposed by Lighthill [35] and then modified for calculating the thrust generated by a forked-shaped caudal fin [12].The fin used in [12] was replaced by the current fin geometry.The driving frequency varied from 1 to 7 Hz for the thrust calculation, as shown in Figure 2. To calculate the drag and lift forces exerted by the fins, we used the general drag force ( ) equation in [4], given by Equations ( 11)-( 13): where  is the forward velocity,  is the relevant reference area of the fin, and  is the drag coefficient given in Equation (12).In Equation ( 12),  is the parasite drag,  is the aspect ratio of a fin,  is the Oswald efficiency factor, which typically has a value of 0.85 to 0.9 for any airfoil section, and  is the coefficient of the lift.To calculate the drag and lift forces exerted by the fins, we used the general drag force (F d ) equation in [4], given by Equations ( 11)-( 13): where U is the forward velocity, A is the relevant reference area of the fin, and C d is the drag coefficient given in Equation (12).In Equation ( 12), C D o is the parasite drag, A r is the aspect ratio of a fin, e is the Oswald efficiency factor, which typically has a value of 0.85 to 0.9 for any airfoil section, and C l is the coefficient of the lift.
To estimate the coefficient of the lift, we referred to Equation (13), where C lα is the slope of the curve with the angle of attack α.We assumed that fin shapes were roughly similar to the shape of the NACA0012 airfoil [36].C lα is obtained by assuming that the maximum α is 15 degrees.
Based on Equations ( 11)-( 13), we calculated the drag forces of the pectoral and dorsal fins.For the body drag, we used Equation (11), assuming a drag coefficient C d of 0.032, which corresponded to the wetted surface area of a fish robot [37].
To control the pitch moment in fish robots, we drove pectoral fins so that the angle of attack was between −15~+15 degrees.The heave force Z cp f and pitch moment M cp f of the pectoral fins are given by Equations ( 14) and ( 15): In Equations ( 14) and ( 15), L p f represents combined lift forces from both left and right pectoral fins, as shown in Figure 3a.The position of the pectoral fin was denoted as (x p , y p , z p ). D p f denotes the drag force exerted by the left and right pectoral fins.During the pitching motion, heave velocity is induced by the lift force, as shown in Equation ( 14).The lift force (F l ) of a fin is given by Equation ( 16): Biomimetics 2024, 9, x FOR PEER REVIEW 6 To estimate the coefficient of the lift, we referred to Equation (13), where  i slope of the curve with the angle of attack .We assumed that fin shapes were rou similar to the shape of the NACA0012 airfoil [36]. is obtained by assuming tha maximum  is 15 degrees.
Based on Equations ( 11)-( 13), we calculated the drag forces of the pectoral and d fins.For the body drag, we used Equation (11), assuming a drag coefficient  of 0 which corresponded to the wetted surface area of a fish robot [37].
To control the pitch moment in fish robots, we drove pectoral fins so that the a of attack was between −15~+15 degrees.The heave force  and pitch moment  the pectoral fins are given by Equations ( 14) and ( 15): In Equations ( 14) and ( 15),  represents combined lift forces from both left right pectoral fins, as shown in Figure 3a.The position of the pectoral fin was denote ( ,  ,  ). denotes the drag force exerted by the left and right pectoral fins.Du the pitching motion, heave velocity is induced by the lift force, as shown in Equation The lift force ( ) of a fin is given by Equation ( 16): The lift and drag forces of the dorsal fin are shown in Figure 3b.Force and mom by the dorsal fin are given by Equations ( 19) and ( 20): The lift and drag forces of the dorsal fin are shown in Figure 3b.Force and moments by the dorsal fin are given by Equations ( 19) and (20): In Equation ( 17), the lift force exerted by the dorsal fin, L d f , also works for the sway directional acceleration force Y cd f .K cd f in Equation ( 18) represents the roll moment by the dor-sal fin.The position of the dorsal fin is represented by (x d , y d , z d ).M cd f in Equation ( 19) denotes the pitch moment resulting from the drag force of the dorsal fin.N cd f in Equation (20) represents the yaw moment by the dorsal fin.Combining Equations ( 15) and ( 19), we obtained the pitch control moment for the pectoral and dorsal fins, as shown in Equation (21).

Simulation and Controller Design Results
We initially conducted an open-loop simulation for our robotic fish, where control inputs were applied by caudal, pectoral, and dorsal fins using Equations ( 10)-(21).Figure 4 shows a simplified simulation block for a fish robot.The outputs measured were surge, sway, and heave velocities, as well as roll, pitch, and yaw rates, in addition to roll, pitch, and yaw angles.
In Equation ( 17), the lift force exerted by the dorsal fin,  , also works for the sway directional acceleration force  . in Equation ( 18) represents the roll moment by the dorsal fin.The position of the dorsal fin is represented by ( ,  ,  ). in Equation (19) denotes the pitch moment resulting from the drag force of the dorsal fin. in Equation ( 20) represents the yaw moment by the dorsal fin.Combining Equations ( 15) and (19), we obtained the pitch control moment for the pectoral and dorsal fins, as shown in Equation (21).

Simulation and Controller Design Results
We initially conducted an open-loop simulation for our robotic fish, where control inputs were applied by caudal, pectoral, and dorsal fins using Equations ( 10)-(21).Figure 4 shows a simplified simulation block for a fish robot.The outputs measured were surge, sway, and heave velocities, as well as roll, pitch, and yaw rates, in addition to roll, pitch, and yaw angles.As shown in Equation (22), two-DOF PID controllers were used for surge, pitch, and yaw controllers.
where , , , and  are the proportional gain, integral gain, derivative gain, and derivative filter time, respectively. and  are the set point weights on the proportional term and derivative term. is the reference input, and  is the output.() is the controller transfer function, and s is the complex frequency in the Laplace transform [38].For practical use, the control loop [39][40][41][42] was designed to satisfy the following preferred conditions: (1) the gain margin should be higher than 10 dB to cover the model uncertainties [43]; and (2) the phase margin should be larger than 45° to cover time delays and model uncertainties [44,45].As shown in Equation ( 22), two-DOF PID controllers were used for surge, pitch, and yaw controllers.
where P, I, D, and N are the proportional gain, integral gain, derivative gain, and derivative filter time, respectively.b and c are the set point weights on the proportional term and derivative term.r is the reference input, and y is the output.G(s) is the controller transfer function, and s is the complex frequency in the Laplace transform [38].For practical use, the control loop [39][40][41][42] was designed to satisfy the following preferred conditions: (1) the gain margin should be higher than 10 dB to cover the model uncertainties [43]; and (2) the phase margin should be larger than 45 • to cover time delays and model uncertainties [44,45].

Surge Response and Velocity Controller Design
To obtain surge responses, the nominal velocity commands were given as 0.2 m/s (Condition 1), 0.4 m/s (Condition 2), and 0.6 m/s (Condition 3).At 10 s, after a steadystate surge velocity was reached, a pseudo-random binary signal (PRBS) control input of ±0.05 m/s was added to the nominal input to check the plant response variations, as shown in Figure 5a.In Figure 5a, the nominal velocities were subtracted to show the velocity variations only.Input and output data were used to obtain transfer functions between velocity command variations and output velocity variations.After obtaining surge open-loop responses from the nonlinear model, we applied the system identification tool in MATLAB to obtain a linearized transfer function.Figure 5b-d show verifications of the identified linear model using different PRBS input signals.The output response of the identified linear model and that of the original nonlinear model were compared.The fits to estimations with nominal velocities of 0.2 m/s, 0.4 m/s, and 0.6 m/s were 61.9%, 91.3%, and 92.8%, respectively.Table 3 shows the linearly approximated transfer functions obtained when the nominal surge velocities were 0.2 m/s, 0.4 m/s, and 0.6 m/s. Figure 5e shows the frequency responses of linearly identified plants.The low-frequency gains decreased as the nominal surge velocity decreased.As the nominal surge velocity increased, the loop gains and bandwidths also increased.Thus, the responses to the control inputs were faster.
the loop gains and bandwidths also increased.Thus, the responses to the control inputs were faster.Two-DOF PID controllers, as shown in Equation ( 22), were designed for the model transfer functions in Table 3. Table 4 shows the surge controller parameters obtained using two-DOF PID tuning.Two-DOF PID controllers, as shown in Equation ( 22), were designed for the model transfer functions in Table 3. Table 4 shows the surge controller parameters obtained using two-DOF PID tuning.Figure 5f shows the frequency responses of plants and PID controllers to check system stability margins.As shown in Table 5, the bandwidths of the three models were 0.90 rad/s, 2.32 rad/s, and 2.19 rad/s, and the phase margins were 69 degrees, 69 degrees, and 70 degrees, respectively.The gain margin was infinite.Figure 5g shows the step responses of a closed-loop linearized plant with the obtained PID controllers.The rise times were about 4, 2, and 2 s.Both the maximum overshoot and steady-state error were small enough to be ignored.Figure 5h shows the step responses of closed-loop nonlinear plant responses when linearly obtained PID controllers are applied.The results showed that PID controllers work well even when they are applied to nonlinear models.Figure 6 shows the step responses of the closed-loop nonlinear plant when the linear controller designed using a nominal velocity of 0.4 m/s is applied to nonlinear plants with different nominal velocities.The results showed that the PID controller, originally optimized for 0.4 m/s, performs effectively with different nominal velocities.This robust performance is due to the fact that the controller has large gain and phase margins, which ensure stability and responsiveness under varying conditions.

Nominal Caudal Frequency (Hz)
Biomimetics 2024, 9, x FOR PEER REVIEW 10 of 17 Figure 6 shows the step responses of the closed-loop nonlinear plant when the linear controller designed using a nominal velocity of 0.4 m/s is applied to nonlinear plants with different nominal velocities.The results showed that the PID controller, originally optimized for 0.4 m/s, performs effectively with different nominal velocities.This robust performance is due to the fact that the controller has large gain and phase margins, which ensure stability and responsiveness under varying conditions.

Pitch Response and Pitch Controller Design
To obtain the pitch response, a PRBS signal of ±0.5 radians was used as the pitch input.Surge velocities were set to 0.2 m/s, 0.4 m/s, and 0.6 m/s.Initially, an open-loop simulation for the pitch was conducted, and input-output data were collected to identify

Pitch Response and Pitch Controller Design
To obtain the pitch response, a PRBS signal of ±0.5 radians was used as the pitch input.Surge velocities were set to 0.2 m/s, 0.4 m/s, and 0.6 m/s.Initially, an open-loop simulation for the pitch was conducted, and input-output data were collected to identify the plant model, as shown in Figure 7a.From the resulting input-output responses, we derived the linearized transfer function using the MATLAB Simulink System Identification Tool.As shown in Figure 7b-d, the nonlinear and linear model responses match well with another PRBS input.The fits to estimations for 0.2 m/s, 0.4 m/s, and 0.6 m/s were 86.14%, 84.3%, and 86.14%, respectively.simulation for the pitch was conducted, and input-output data were collected to identify the plant model, as shown in Figure 7a.From the resulting input-output responses, we derived the linearized transfer function using the MATLAB Simulink System Identification Tool.As shown in Figure 7b-d, the nonlinear and linear model responses match well with another PRBS input.The fits to estimations for 0.2 m/s, 0.4 m/s, and 0.6 m/s were 86.14%, 84.3%, and 86.14%, respectively.
The linear transfer functions obtained are shown in Table 6. Figure 7e shows the frequency responses of the obtained linear models.A two-DOF PID controller was designed using PID tuning for the linear transfer function, as shown in Table 6.The parameters of the two-DOF PID controller obtained are shown in Table 7.The linear transfer functions obtained are shown in Table 6. Figure 7e shows the frequency responses of the obtained linear models.A two-DOF PID controller was designed using PID tuning for the linear transfer function, as shown in Table 6.The parameters of the two-DOF PID controller obtained are shown in Table 7. Figure 7f shows the frequency responses of the linear pitch plant models with feedback controllers to check the system's performances.Table 8 shows the gain and phase margins and bandwidths of the closed-loop control systems.The gain margins are large enough, even though they are not infinite, as in the surge control.The phase margins are larger than 76 degrees, which is more than enough.Figure 7g shows the step responses of closed-loop pitch controllers for linear plant models with different nominal surge velocity conditions.The maximum rise time was about 20 s.Both maximum overshoots and steady-state errors were negligibly small.Figure 7h shows step responses when the PID controllers derived from linear models are applied to the nonlinear model.From the start to 10 s, pitch reference inputs remained at zero.Due to the surge velocities being changed from zero to nominal velocities, there were some oscillations in the pitch angle.At 10 s, 0.01-rad pitch commands were applied to the pitch controllers.It could be seen that all pitch-angle outputs were following the commands.The rise time was about 5 s, demonstrating that the obtained PID controller also performs well with the nonlinear plant.Figure 8 shows the step response of a nonlinear plant when the linearly derived pitch controller for the 0.4 m/s model is applied with different nominal surge velocities.When the nominal surge velocities are 0.2 m/s, 0.4 m/s, and 0.6 m/s, the rise time is about 15 s, 2 s, and 1 s, respectively.When the surge velocity was low, it took more time to maneuver the pitch because the pitching force is proportional to the square of the surge velocity.
Figure 8 shows the step response of a nonlinear plant when the linearly derived pitch controller for the 0.4 m/s model is applied with different nominal surge velocities.When the nominal surge velocities are 0.2 m/s, 0.4 m/s, and 0.6 m/s, the rise time is about 15 s, 2 s, and 1 s, respectively.When the surge velocity was low, it took more time to maneuver the pitch because the pitching force is proportional to the square of the surge velocity.

Yaw Response and Yaw Controller Design
In the yaw control loop simulation, the surge nominal velocities were given as 0.2 m/s, 0.4 m/s, and 0.6 m/s.The PRBS input with an amplitude of 0.5 rad was given as a yaw control input after 10 s, as shown in Figure 9a.After applying the system identification tool in MATLAB Simulink to the yaw command input and output data, we obtained a linear transfer function describing characteristics between the yaw command input and output.The obtained linear transfer functions are shown in Table 9. Figure 9b, 9c, and 9d show the responses of the identified linear and nonlinear models.For the 0.2 m/s, 0.4 m/s, and 0.6 m/s models, the fits to estimations were 72.59%, 64.09%, and 64.37%, respectively.Figure 9e shows the frequency responses of the obtained linear models.

Yaw Response and Yaw Controller Design
In the yaw control loop simulation, the surge nominal velocities were given as 0.2 m/s, 0.4 m/s, and 0.6 m/s.The PRBS input with an amplitude of 0.5 rad was given as a yaw control input after 10 s, as shown in Figure 9a.After applying the system identification tool in MATLAB Simulink to the yaw command input and output data, we obtained a linear transfer function describing characteristics between the yaw command input and output.The obtained linear transfer functions are shown in Table 9. Figure 9b, Figure 9c, and Figure 9d show the responses of the identified linear and nonlinear models.For the 0.2 m/s, 0.4 m/s, and 0.6 m/s models, the fits to estimations were 72.59%, 64.09%, and 64.37%, respectively.Figure 9e shows the frequency responses of the obtained linear models.Two-DOF PID controllers were designed and tuned in Simulink using transfer functions obtained from the system identification, as shown in Table 9.The parameters of these two-DOF PID controllers are shown in Table 10. Figure 9f shows the frequency responses of the linearized yaw plant with feedback controllers.Table 11 shows the gain and phase margins for the obtained controller.Figure 9g shows the closed-loop responses of the linearized models with the obtained controller.The results showed that for yaw control systems with surge velocities of 0.2 m/s, 0.4 m/s, and 0.6 m/s, the rise times were 2 s, 2.5 s, and 3 s, respectively.The maximum overshot and steady-state error were small when the nominal surge velocities were 0.4 m/s and 0.6 m/s.When the nominal surge velocity was 0.2 m/s, the output had some overshoot and slowly reached a steady state.Figure 9h shows the step responses of a nonlinear plant with linearly obtained PID controllers.At surge velocities of 0.2 m/s, 0.4 m/s, and 0.6 m/s, the rise times were 2 s, 2.7 s, and 5 s, respectively.Two-DOF PID controllers were designed and tuned in Simulink using transfer functions obtained from the system identification, as shown in Table 9.The parameters of these two-DOF PID controllers are shown in Table 10.  Figure 10 shows the step responses of the linearly derived yaw controller with 0.4 m/s nominal surge velocity applied to the nonlinear model with different nominal surge velocities.When the nominal surge velocities were 0.2 m/s, 0.4 m/s, and 0.6 m/s, the rise time was about 4.5 s, 2.5 s, and 2.5 s, respectively.It shows that the applied controller worked well even with different conditions.

Conclusions
In this paper, a simulation block for a fish robot was constructed using MATLAB Simulink.Linearly approximated transfer function models were obtained from surge, pitch, and yaw open-loop responses.Using these linear models, we designed a simple two-DOF PID controller for surge, pitch, and yaw control.A frequency analysis using Bode diagrams for the surge, pitch, and yaw controllers showed that the gain margin and phase margin were sufficiently large for the surge, yaw, and pitch controllers.
When surge, pitch, and yaw controllers were applied to the linear model, the maximum overshoot and steady-state error were small.The bandwidths of the surge pitch and yaw controllers were 2.3 rad/s, 0.17 rad/s, and 2 rad/s, respectively.When these controllers were applied directly to a nonlinear plant, the response characteristics were similar to those when the controllers were applied to a linear model.Therefore, it is expected that

Conclusions
In this paper, a simulation block for a fish robot was constructed using MATLAB Simulink.Linearly approximated transfer function models were obtained from surge, pitch, and yaw open-loop responses.Using these linear models, we designed a simple two-DOF PID controller for surge, pitch, and yaw control.A frequency analysis using Bode diagrams for the surge, pitch, and yaw controllers showed that the gain margin and phase margin were sufficiently large for the surge, yaw, and pitch controllers.
When surge, pitch, and yaw controllers were applied to the linear model, the maximum overshoot and steady-state error were small.The bandwidths of the surge pitch and yaw controllers were 2.3 rad/s, 0.17 rad/s, and 2 rad/s, respectively.When these controllers were applied directly to a nonlinear plant, the response characteristics were similar to those when the controllers were applied to a linear model.Therefore, it is expected that the proposed controller can be applied to control real fish robots with nonlinear coupling terms and uncertainties.
In simulations, we assumed a symmetric ellipsoidal body and obtained parameters simplifying the structure.Thus, these parameters need to be updated using more accurate methods and actual experimental data.In this research, we focused on single-input and single-output responses.However, in future work, we will consider multi-input and multioutput responses to see more complex system behaviors.Additionally, the current control law does not fully address interactions between surges and other forces, which could lead to instability.
In the future, we plan to fabricate an actual fish robot.With a real fish robot, we will perform experiments to identify dynamic parameters and transfer characteristics and implement controllers to verify their performances.

Figure 1 .
Figure 1.Fish robot coordinate system: (a) side view of the fish robot; (b) top view of the fish robot.

Figure 1 .
Figure 1.Fish robot coordinate system: (a) side view of the fish robot; (b) top view of the fish robot.

Figure 3 .
Figure 3. Pectoral and dorsal moments and forces: (a) side view of the fish robot; (b) top view o fish robot.

Figure 3 .
Figure 3. Pectoral and dorsal moments and forces: (a) side view of the fish robot; (b) top view of the fish robot.

Figure 5 .
Figure 5. Surge velocity responses: (a) open-loop responses to speed commands with different nominal velocities; (b) responses of the identified linear model and nonlinear plant at a nominal speed of 0.2 m/s; (c) responses of the identified linear model and nonlinear plant at a nominal speed of 0.4 m/s; (d) responses of the identified linear model and nonlinear plant at a nominal speed of 0.6 m/s; (e) frequency response of identified linear models; (f) Bode plots of identified plants with feedback controllers; (g) step responses of the closed-loop linearized models; (h) step responses of the closedloop nonlinear model.

Figure 5 .
Figure 5. Surge velocity responses: (a) open-loop responses to speed commands with different nominal velocities; (b) responses of the identified linear model and nonlinear plant at a nominal speed of 0.2 m/s; (c) responses of the identified linear model and nonlinear plant at a nominal speed of 0.4 m/s; (d) responses of the identified linear model and nonlinear plant at a nominal speed of 0.6 m/s; (e) frequency response of identified linear models; (f) Bode plots of identified plants with feedback controllers; (g) step responses of the closed-loop linearized models; (h) step responses of the closed-loop nonlinear model.

Figure 6 .
Figure 6.Step responses of the closed-loop nonlinear surge model with a 0.4 m/s PID controller.

Figure 6 .
Figure 6.Step responses of the closed-loop nonlinear surge model with a 0.4 m/s PID controller.

Figure 7 .
Figure 7. Pitch responses: (a) open-loop pitch responses at surge velocities of 0.2 m/s, 0.4 m/s, and 0.6 m/s; (b) responses of the identified linear model and nonlinear plant at a nominal speed of 0.2 m/s; (c) responses of the identified linear model and nonlinear plant at a nominal speed of 0.4 m/s; (d) responses of the identified linear model and nonlinear plant at a nominal speed of 0.6 m/s; (e) frequency response of the identified linear models; (f) Bode plots of identified plants with feedback controllers; (g) step responses of closed-loop controllers with linearized models; (h) step responses of closed-loop controllers with a nonlinear model.

Figure 7 .
Figure 7. Pitch responses: (a) open-loop pitch responses at surge velocities of 0.2 m/s, 0.4 m/s, and 0.6 m/s; (b) responses of the identified linear model and nonlinear plant at a nominal speed of 0.2 m/s; (c) responses of the identified linear model and nonlinear plant at a nominal speed of 0.4 m/s; (d) responses of the identified linear model and nonlinear plant at a nominal speed of 0.6 m/s; (e) frequency response of the identified linear models; (f) Bode plots of identified plants with feedback controllers; (g) step responses of closed-loop controllers with linearized models; (h) step responses of closed-loop controllers with a nonlinear model.

Figure 8 .
Figure 8. Step responses of the nonlinear plant applied with a linearly obtained pitch controller of 0.4 m/s.

Figure 8 .
Figure 8. Step responses of the nonlinear plant applied with a linearly obtained pitch controller of 0.4 m/s.

Figure 9 .
Figure 9. Yaw simulation responses: (a) open-loop yaw angle responses when PRBS yaw command input is applied; (b) responses of the identified linear model and nonlinear plant when the nominal velocity is 0.2 m/s; (c) responses of the identified linear model and nonlinear plant when the nominal velocity is 0.4 m/s; (d) responses of the identified linear model and nonlinear plant when the nominal velocity is 0.6 m/s; (e) frequency responses of the identified linear models; (f) Bode plots of the identified linear models with feedback controllers; (g) step responses of the closed-loop linearized model with PID controllers; (h) step responses of the closed-loop nonlinear model with linearized PID controllers.

Figure 9 .
Figure 9. Yaw simulation responses: (a) open-loop yaw angle responses when PRBS yaw command input is applied; (b) responses of the identified linear model and nonlinear plant when the nominal velocity is 0.2 m/s; (c) responses of the identified linear model and nonlinear plant when the nominal velocity is 0.4 m/s; (d) responses of the identified linear model and nonlinear plant when the nominal velocity is 0.6 m/s; (e) frequency responses of the identified linear models; (f) Bode plots of the identified linear models with feedback controllers; (g) step responses of the closed-loop linearized model with PID controllers; (h) step responses of the closed-loop nonlinear model with linearized PID controllers.

Figure 10
Figure10shows the step responses of the linearly derived yaw controller with 0.4 m/s nominal surge velocity applied to the nonlinear model with different nominal surge velocities.When the nominal surge velocities were 0.2 m/s, 0.4 m/s, and 0.6 m/s, the rise time was about 4.5 s, 2.5 s, and 2.5 s, respectively.It shows that the applied controller worked well even with different conditions.

Figure 10 .
Figure 10.Step responses of the closed-loop nonlinear model with a 0.4 m/s yaw controller.

Figure 10 .
Figure 10.Step responses of the closed-loop nonlinear model with a 0.4 m/s yaw controller.

Table 1 .
Parameters of the fish robot.

Table 3 .
Linearly obtained surge velocity transfer functions with different nominal surge velocities.

Table 4 .
Parameters of two-DOF PID surge controllers.

Table 4 .
Parameters of two-DOF PID surge controllers.

Table 5 .
Gain and phase margins for surge controllers.

Table 6 .
Linearly obtained pitch transfer functions.

Table 7 .
Parameters of two-DOF PID pitch controllers.

Table 7 .
Parameters of two-DOF PID pitch controllers.

Table 8 .
Gain and phase margins for pitch controllers.

Table 10 .
Parameters of two-DOF PID yaw controllers.

Table 11 .
Gain and phase margins for yaw controllers.