Similarity Evaluation Rule and Motion Posture Optimization for a Manta Ray Robot

Abstract

The current development of manta ray robots is usually based on functional bionics, and there is a lack of bionic research to enhance the similarity of motion posture. To better exploit the characteristics of bionic, a similarity evaluation rule is constructed herein by a Dynamic Time Warping (DTW) algorithm to guide the optimization of the control parameters of a manta ray robot. The Central Pattern Generator (CPG) network with time and space asymmetry oscillation characteristics is improved to generate coordinated motion control signals for the robot. To optimize similarity, the CPG network is optimized with the genetic algorithm and particle swarm optimization (GAPSO) to solve the problems of multiple parameters, high non-linearity, and uncertain parameter coupling in the CPG network. The experimental results indicate that the similarity between the forward motion pose of the optimized manta ray robot and the manta ray is improved to 88.53%.

1. Introduction

To further improve the stability, mobility, and bioaffinity of robots, the development of underwater vehicles using bionic means has become an international scientific research hotspot [1,2]. Through the long-term observation and study of the shape, structure, and movement characteristics of various organisms, researchers have found that natural selection and evolution have enabled fish to exhibit extraordinary swimming abilities in the water, with strong explosive power and flexible mobility [3,4]. The adaptability of organisms to their natural environment is the result of gradual evolution over millions of years. Among fish, manta rays exhibit high propulsion efficiency, stability, and maneuverability due to their streamlined, flattened bodies and large span-to-chord ratio pectoral fins [5,6]. Manta rays are excellent bionic objects for underwater robots.

A variety of robots that mimic the morphology and movement of manta rays have emerged [7], and a fundamental problem that needs to be solved is how to control the movement of the pectoral and caudal fins to achieve better locomotor capabilities. Neurobiologists believe that the rhythmic behavior of animals is controlled by the central pattern generator (CPG) network, which is a local oscillatory network composed of intermediate neurons [8]. It achieves self-excited oscillations through mutual inhibition between neurons, generating stable periodic signals that lead to the rhythmic movements of limbs and trunk positions. A few studies have examined the biologically inspired swimming control of robotic fish using CPG [9,10]. A CPG network based on the Hopf oscillator was used to control the robot fish Roman-ii, and the modes of forward, backward, turn, and switch between the modes were realized [11,12]. The CPG-based bionic motion control method avoids the establishment of a complex dynamics model of a robot, can realize the control of basic motion modes, and has good adaptability to the environment [13]. The determination of CPG model parameters is mainly based on experience. Each control parameter plays an important role in the realization of motion mode, and parameter optimization is crucial in bionic control based on CPGs. Cui, X [14] realized the motion planning of hexapod robots through a multi-objective genetic algorithm based on the CPG model. Wang, M [15,16] used a particle swarm optimization (PSO) algorithm to optimize CPG parameters, which enables the caudal swinging robot fish to achieve a higher swimming speed.

From the above studies, scholars on CPG optimization have focused on improving the speed and efficiency of robots, ignoring the bionic robot’s motion posture optimization. In contrast, the extraordinary locomotor ability shown by the locomotion and gestures formed by organisms in the evolutionary process is of great guidance for the optimization direction of the locomotion control of bionic robots.Zhu, J [17] designed and built a robotic fish that mimics yellowfin tuna (Thunnus albacares) and Atlantic mackerel (Scomber scombrus), demonstrating the ability of the organisms to swim at high frequencies by measuring body kinematics, speed, and power. Romano, D. [18] designed a female-mimicking robotic that mimics the behavior and appearance of female fish color patterns, successfully achieving a response from real male fish to the robotic fish, providing an advanced strategy for feature-based ecological studies. Liu, S [19] designed a fluid-driven soft robotic fish that mimics the red muscular system of fish. Wu, J [20] presented a new approach to the implementation of active flapping wing motion by imitating the movement of manta rays, and they verified that the trajectory of the bionic robot is consistent with the fluttering shape of the manta ray. According to the above research, the development of biomimetic robotic fish research has been facilitated by a large number of studies on fish behavior. However, there is no clear quantitative standard for the degree of similarity of mimicry.

As for how to evaluate the simulation degree of bionic robots, researchers have defined a similarity evaluation method for 7-linked bipedal walking robots. They are based on spatial-time control from forward and inverse kinematic solving and motion model simplification and reorientation design methods [21]. From the perspective of motion rhythm control, researchers have proposed a similarity calculation rule for humanoid robots based on base segmentation [22,23]. Motor performance aspects such as maximum pitch angle (MPA), maximum reach height (MRH), maximum bending angle (MBA), and minimum bending distance (MBD) are used to establish the motor similarity evaluation model of mouse-like robots [24]. For underwater robots with obvious active and passive deformation of the motion mechanism, no scholars have yet established a complete evaluation system, which limits optimization research in the direction of the bionic similarity of manta ray robots.

In order for a manta ray robot to attain better bionic characteristics that are able to simulate real manta rays, the motion posture of the manta ray robot should be similar to that of a manta ray. To solve the problem of how to directionally improve the bionic similarity of robot motion posture, the similarity evaluation rules of the manta ray robot are established in this paper. The similarity evaluation rules of motion posture were established by feature point trajectory extraction and a dynamic time warping (DTW) algorithm. By introducing the bias equations and time asymmetry coefficients, a phase oscillator model characterizing the time and space asymmetry of the pectoral fin motion of a manta ray was established. A CPG topology network based on a phase oscillator was constructed to give a multi-drive structure to the pectoral fin of the manta ray robot. Based on this, an adaptation function for the similarity of the motion posture was established to improve the similarity of the motion posture of the manta ray robot, and a genetic algorithm and particle swarm optimization (GAPSO) method for CPG network parameters is proposed. The similarity evaluation rule and the parameter optimization results were verified by the pool experiments.

2. Method

2.1. Forward Swimming Posture of Manta Ray

Forward swimming is one of the most important movement postures of manta rays. In the forward swimming posture, the pectoral fins on both sides of the manta ray are flapping symmetrically, generating fluctuation transmission in both the spreading and chordal directions, thus generating forward thrust and upward lift, as seen in Figure 1b. The motion of the posterior edge of the pectoral fin lags behind the anterior edge. The maximum amplitude of the upstroke is much larger than that of the downstroke, and the time spent on upstrokes is greater than that spent on downstrokes, showing obvious space and time asymmetry.

To further analyze the kinematics of the forward swimming of manta rays in the lateral view, three typical points on the pectoral fin of manta rays were selected for kinematic analysis. The coordinate system is established with the symmetrical center of the manta’s head as the origin of the coordinate system. The X-axis is the longitudinal direction of the pectoral fin, the Y-axis is the chordal direction of the trunk, and the Z-axis perpendicular to the trunk plane can be obtained according to the rule. When observed from the side view angle, the chordal motion of the feature point is mainly reflected in the Z-axis motion, and the displacement change of the X-axis is negligible. The motion trajectory of the feature point is represented in Figure 2.

As can be seen in Figure 2, the curves of the motion characteristics of points $P_1$, $P_2$, and $P_3$ over time are sinusoidal-like. $P_1$ differs in the amplitude and time taken to beat up and down, and $P_1$, $P_2$ and $P_3$ differ in the time taken to reach the extremes. These differences indicate that the typical characteristic point movements on the pectoral fins of manta rays when they flutter forward are sinusoidal-like with temporal asymmetry, spatial asymmetry, and phase difference properties. The quantitative analysis of $P_i$ is shown in Equation (1). In order to better guide the manta robot to realize a motion posture similar to that of a manta ray, we conduct the parameter design of the motion posture of the manta robot in this paper. The parameters of feature points in Equation (1) are quantitatively analyzed. The kinematic equation is summarized in Equation (2).

$$z_i=R_{xi}+R_i\sin (2\pi v_{i}t+\Delta {\phi }_{ij})$$

(1)

$$\{\begin{array}{l}{R}_1:R_2:R_3=1.6:0.6:0.3\\ R_{x1}:R_{x2}:R_{x3}=0.3:0.1:-0.5\\ v_1=v_2=v_3=0.4Hz\\ \Delta {\phi }_{12}=9.5°,\Delta {\phi }_{13}=52.3°\end{array}$$

(2)

Here, $R_{xi}$ is the bias value caused by the different amplitudes of up and down and can be calculated according to $(R_{i\_up}+R_{i\_down})/2-R_{i\_down}$. $R_i$ is the desired amplitude, which represents the maximum value that flaps with the bias value as the center and can be calculated according to $(R_{i\_up}+R_{i\_down})/2$. $v_i$ is the flapping frequency, which is calculated based on the time it takes to beat up and down, $v_i=1/(T_{i\_up}+T_{i\_down})$. $\Delta {\phi }_{ij}$ is the phase difference between $P_i$ and $P_j$, calculated from the time difference between $P_i$ and $P_j$ when they reach their extreme points.

2.2. Manta Ray Robot Design and Control

2.2.1. Robot Design and Its Pectoral Fin Flexible Deformation Analysis

To achieve a manta ray robot with a similar posture to the manta ray, the main morphological features of the manta ray robot should be the same as those of a manta ray. In this paper, the design of the manta ray robot is modeled on the biological appearance of the manta ray and the structural features of the pectoral and caudal fins. In order to resemble the manta ray as closely as possible, a 2D front, side, and top view of the manta ray was taken and used as a prototype to design and produce an accurate 3D model of the manta ray robot. The manta ray robot is mainly composed of the main structure, pectoral fin structure, and tail fin structure. The overall dimensions are 934 $\times$ 570 $\times$100 mm and the total weight is 10 kg, as shown in Figure 3.

The body of the manta ray robot has a flattened base and a raised back made of nylon; inside the body, there are basic electronics such as the main control board, attitude sensor, battery, and communication equipment. The pectoral fins are made of a three-stage flutter fin structure that mimics the pectoral fin skeleton of a manta ray robot, with each stage driven by a waterproof servo. The pectoral fins are connected to each other using silicone to mimic the musculature and structure of the pectoral fins. The caudal fin section is driven by a waterproof servo with a caudal fin plate, the outer side of which is silicone to mimic the shape of the biological caudal fin. Both the pectoral and caudal fins are made of silicone with a gradual change in thickness to mimic the flexibility of a biological caudal fin.

The forward swimming posture of the manta ray is mainly achieved via pectoral fin flapping, and its motion performance determines the swimming performance of the robot. The pectoral fins on both sides of the manta ray robot are completely symmetrical; therefore, the analysis of the motion performance of the pectoral fins is only performed on one side of the pectoral fins. If the flexible deformation of the fins is ignored, the displacement of the endpoint of the fins along the Z-axis direction $h_i$ is related to the rudder output angle ${\theta }_i$ as follows:

$$h_i=l_i\sin ({\theta }_i)$$

(3)

However, under the condition that the current and the flexibility of the fins themselves exist, the fins exhibit not only active deformation controlled by the rudder output but also passive flexibility along the spreading direction [25]. The actual flap of the fins underwater is smaller than that under active control. Taking the upstroke flapping of pectoral fins as an example, the gray part represents the actual motion deformation of strings due to fluid load, as shown in Figure 4a.

The flexible deformation of the fins varies approximately according to the parabolic law along the wingspan direction [26,27]. The magnitude of flexural deformation varies continuously with torque. Take the torque at the tip as an example, when the flap starts to flap down from the highest point, the torque increases gradually. Then, the torque is at its maximum when it passes the horizontal position, following which it starts to decrease, and it reaches the minimum value when the flap flaps up and passes the horizontal position again. Therefore, the expression of flexible deformation is shown in Equation (4).

$${\theta }_i(t)={\theta }_i+{\theta }_{tip}(t){(\frac{x}{b})}^2{(\frac{y}{c})}^2$$

(4)

Here, $b$ is half the span length, $c$ is the chord length of the fin root, ${\theta }_i$ is the angle of active deformation, and ${\theta }_{tip}$ is the wing tip torsion angle, both of which are shown in Equations (5) and (6):

$${\theta }_i=\sin (2\pi v_{i}t)$$

(5)

$${\theta }_{tip}(t)={\theta }_{tip\_\max }\sin (2\pi v_{i}t-\frac{\pi }{2})$$

(6)

Substituting Equations (5) and (6) into Equation (4), the expression for the displacement $h_i$ of the endpoint of the root fin along the Z-axis direction is shown in Equation (7). The simulation results for $h_i$ are shown in Figure 4b,c.

$$h_i=l_i\sin ({\theta }_i+{\theta }_{tip\_\max }\sin (2\pi v_{i}t-\frac{\pi }{2}){(\frac{x}{b})}^2{(\frac{y}{c})}^2)$$

(7)

2.2.2. Construction of CPG Phase Oscillator

The manta ray’s swimming posture is controlled by the CPG. The manta ray’s fin rays are driven by the corresponding muscle groups, which are controlled by CPG neural units. All units form CPG neural networks through complex connections to control the phase difference, swing frequency, swing amplitude, and other parameters between the fin rays, thus forming the manta ray’s complex kinematic characteristics. If an artificial CPG neural control network can be established to simulate the real CPG neural network, the manta ray robot will achieve bionics from morphology to motion posture. The artificial CPG control method needs to meet the kinematic characteristics of manta rays so that the pectoral fin can have the sinusoidal motion characteristics of time and space asymmetry, and the swing amplitude, frequency, and phase difference between fins can be controlled. Commonly used CPG non-linear oscillator models include the phase oscillator model, the recursive oscillator model, and the Hopf oscillator model. The phase oscillator model has meaningful parameters such as amplitude, phase difference, and frequency, and is more suitable for the kinematic bionic control of manta ray robots [8]. Therefore, the phase oscillator model is chosen to construct the CPG neural network in this paper. The traditional phase oscillator model is shown in Equation (8):

$$\{\begin{array}{l}{\ddot{r}}_i=a_i\left[\frac{a_i}{4}\left(R_i-r_i\right)-{\dot{r}}_i\right]\\ {\dot{\varphi }}_i=2\pi v_i+{\displaystyle \sum _j{\omega }_{ij}\sin \left({\varphi }_j-{\varphi }_i-\Delta {\phi }_{ij}\right)}\\ {\theta }_i=r_i\left[1+\cos \left({\varphi }_i\right)\right]\end{array}$$

(8)

Here, the first equation is the amplitude equation. $r_i$ is the amplitude, $a_i$ is the amplitude convergence coefficient, which can control the amplitude convergence rate, and $R_i$ is the expected amplitude. The second equation is the phase equation. ${\varphi }_i$ is the phase in the current state; ${\omega }_{ij}$ is the coupling weight, which represents the coupling weight of unit $j$ to unit $i$; and $\mathsf{\Delta }{\phi }_{ij}$ is the expected phase difference between oscillators $i$ and $j$. The third equation is the output equation, and ${\theta }_i$ is the output of oscillator i, which is determined by both ${\varphi }_i$ and $r_i$.

The forward motion of manta rays has the characteristics of space asymmetry and time asymmetry, while the phase oscillator equation can only output a cosine signal greater than 0. To better achieve the purpose of imitating real manta ray motion, this paper improves the initial phase oscillator model by introducing an amplitude bias to realize the space asymmetric flutter characteristics and setting different up and down strike frequencies to realize the time asymmetric characteristics.

The bias equation for the amplitude is constructed from the amplitude equation shown in Equation (9):

$${\ddot{r}}_{xi}=m_i\left[\frac{m_i}{4}\left(R_{xi}-r_{xi}\right)-{\dot{r}}_{xi}\right]$$

(9)

where $r_{xi}$ is the bias, $m_i$ is the bias convergence coefficient, and $R_{xi}$ is the desired bias.

At the same time, it is defined as the pectoral fin flapping mode when the space asymmetry coefficient ${\alpha }_i$ line $R_i$ is not equal to 0, and ${\lambda }_i$ represents the ratio of the maximum Z-axis displacement of the fin in the whole down (or up) stroke. The space asymmetry coefficient ${\lambda }_i$ is shown in Equation (10).

$${\lambda }_i=\frac{R_i+R_{xi}}{2R_i},R_i\ne 0$$

(10)

By combining Equations (5) and (6), the amplitude bias equation is shown in Equation (11):

$${\ddot{r}}_{xi}=m_i\left[\frac{m_i}{4}((2{\lambda }_i-1)R_i-r_{xi})-{\dot{r}}_{xi}\right],R_i\ne 0$$

(11)

The difference between the upstroke travel time and the downstroke travel time can be expressed in terms of frequency $v_i$. When the pectoral fin of the manta ray is flapping up $v_i=v_{i\_up}$, and when the pectoral fin of the manta ray is flapping down $v_i=v_{i\_down}$. To determine whether the pectoral fin is in an upstroke or a downstroke, the positive or negative value of $\stackrel{·}{{\theta }_i}$ can be used. At the same time, this paper introduces a time asymmetric flapping coefficient ${\eta }_i$, defined as the ratio of the fin up-tempo travel time to the flapping period, so that the frequency equation can be expressed as Equation (12):

$$\{\begin{array}{cc}{v}_i=v_{i\_up}=\frac{1}{2{\eta }_i{T}_i}, {\dot{\theta }}_i>0& \\ v_i=v_{i\_down}=\frac{1}{2(1-{\eta }_i)T_i}, {\dot{\theta }}_i\le 0& \end{array}$$

(12)

Substituting Equations (11) and (12) into the original phase oscillator model, the improved phase oscillator model is shown in Equation (13), the meanings of model parameters are shown in

Table 1. Parameters of the phase oscillator model.

Parameter	Meaning	Parameter	Meaning
$v_i$	flapping frequency	${\eta }_i$	time asymmetric coefficient
$T_i$	flapping period	$r_i$	amplitude
$a_i$	amplitude convergence coefficient	$R_i$	expected amplitude
$r_{xi}$	bias	$m_i$	bias convergence coefficient
${\lambda }_i$	pace asymmetry coefficient	${\varphi }_i$	phase
${\omega }_{ij}$	coupling weight	$\Delta {\phi }_{ij}$	expected phase

.

In the manta ray robot, the CPG unit corresponding to the first level of fins on one side plays a more prominent role, and the pectoral fin units on one side are closely linked, while the pectoral fin units on the left and right sides are relatively weakly linked, so this paper chooses the simplest form of connection to build the CPG network, as shown in Figure 5a.

$$\{\begin{array}{c}\{\begin{array}{cc}{v}_i=v_{i\_up}=\frac{1}{2{\eta }_i{T}_i}, {\dot{\theta }}_i>0& \\ v_i=v_{i\_down}=\frac{1}{2(1-{\eta }_i)T_i}, {\dot{\theta }}_i\le 0& \end{array}\\ {\ddot{r}}_i=a_i\left[\frac{a_i}{4}\left(R_i-r_i\right)-{\dot{r}}_i\right]\\ {\ddot{r}}_{xi}=m_i\left[\frac{m_i}{4}\left((2{\lambda }_i-1)R_i-r_{xi}\right)-{\dot{r}}_{xi}\right], R_i\ne 0\\ {\dot{\varphi }}_i=2\pi v_i+{\displaystyle \sum _j{\omega }_{ij}\sin \left({\varphi }_j-{\varphi }_i-\Delta {\phi }_{ij}\right)}\\ {\theta }_i=r_{xi}+r_i\sin ({\varphi }_i)\end{array}$$

(13)

It can be seen from the simulation result that the oscillator successfully realizes the signal output of space asymmetry and time asymmetry, as shown in Figure 5b. The characteristics of the output waveform of the improved phase oscillator model have been preliminarily consistent with that of a real manta ray. To maintain similarity with the motion trajectory of the pectoral fin of a manta ray when it swims forward, the parameters in the oscillator need to be further optimized.

2.3. Bionic Similarity Evaluation Rule Establishment

The manta ray robot can visually imitate the movements of a manta ray, but the visual similarity is not a scientific method of determining whether or not the imitation is accurate. The lack of a quantitative method to evaluate this similarity makes it difficult to improve the imitation of motion and limits the application of optimal control strategies for manta ray robots. The trajectory in space can often be directly expressed as the posture of the robot in the process of motion. Therefore, this paper establishes a similarity evaluation rule for quantitative motion states.

Due to poor controllability and reproducibility in animal experiments, it is difficult to obtain good reproducibility of the observed manta ray motion behavior. However, in order to reduce the variability between manta ray behaviors, we observed a lot of manta ray behavioral activities. For the quantitative analysis of the similarity of the motion posture, this paper will use the trajectory relationship between the angles of the pectoral fins to evaluate the similarity of the motion posture of the manta ray robot and the manta ray.

For the selection of feature points, this paper takes the endpoints of the pectoral fins as the feature points from the perspective of the manta ray robot, and the corresponding feature points of the pectoral fins of the manta ray are selected according to their proportional relationship, as shown in Figure 6. Then, the trajectory coordinate points of the manta ray feature points can be represented by a vector $M_i=\left\{q_{i1},q_{i2},q_{i3},\dots ,q_{in}\right\}(i=1,2,3)$. The trajectory coordinate points of the manta ray robot feature points can be represented by a vector $Q_i=\left\{w_{i1},w_{i2},w_{i3},\dots ,w_{in}\right\}(i=1,2,3)$. The difficulty of photographing underwater moving objects leads to differences in synchronization when sampling feature points of the manta ray mimicking robot and manta rays. The time series of feature points exist roughly similarly in terms of time but are not in one-to-one correspondence in the time series. Calculation with a traditional distance algorithm without considering the dynamic changes in time will cause big errors.

DTW is a dynamic programming algorithm for calculating the similarity of two sequences, especially those of different lengths, and is mainly applied to time-series data [28,29]. DTW matches the data points of a time-series by bending the time domain of the time-series, not only to obtain a better morphological measure but also to measure two sequences of unequal length. This paper assumes that two-time series A and B are represented as $A:\left\{a_1,a_2,\cdots ,a_i,\cdots ,a_n\right\}$ and $B:\left\{b_1,b_2,\cdots ,b_i,\cdots ,b_m\right\}$. DTW can be expressed as the following objective function:

$$\{\begin{array}{l}g(i,j)=d\left(a_i,b_j\right)+\min \{g(i-1,j-1),g(i-1,j),g(i,j-1)\}\\ DTW(A,B)=g(m,n)\\ i=2,3,\cdots n;j=2,3,\cdots m;\end{array}$$

(14)

where $d\left(i_k,j_k\right)=|a(i)-b(j)|$ denotes the absolute error between $a(i)$ and $b(j)$, and $g(i,j)$ is the cumulative distance of local distances on the path from $(a(1),b(1))$ to $(a(i),b(i))$.

The analysis of the DTW algorithm above shows that it can calculate the similarity of two time series of different lengths. In this section, the DTW algorithm is used to evaluate the similarity $S_{ma}(M_i,Q_i)$ between the motion pose of the manta ray robot and the manta ray.

$$S_{ma}(M_i,Q_i)=\frac{1}{1+\frac{DTW(M,Q)}{k}}$$

(15)

Here, $k$ is the number of feature points.

2.4. Optimization Study on the Similarity of Forward Swimming Posture

Forward swimming is the main motion posture of manta rays, which is the basic mode for manta rays to complete multimodal motion and the prerequisite for the manta ray robot to be more like manta rays. The improved CPG control model can mimic the multimodal motion control of real manta rays with good stability and adjustability. However, there are many parameters in the CPG network, and it is very difficult to adjust them by manual trial or experimental method. To improve the similarity of the forward swimming posture of the manta ray robot faster and more accurately, this paper establishes the adaptation function of the posture similarity and combines the optimization algorithm for the CPG network parameters. The forward swimming posture of the manta ray only requires both pectoral fins. The caudal fin is hardly involved. Pectoral fins flutter symmetrically when the manta ray swims forward, so this paper only optimizes the CPG model parameters for one side of the pectoral fins to complete the similarity optimization of the forward swimming posture.

The standard PSO algorithm accomplishes the optimal search by tracking the individual optimum and the population optimum, which is simple and converges quickly [30,31]. However, as the number of iterations increases, each particle becomes more and more similar, and it is easy to fall into a local optimum that cannot be jumped out of. GAPSO is used in this paper. In this algorithm, GA is introduced into particle swarm optimization [32,33] and then the crossover and mutation steps in the genetic algorithm are used to obtain a better next-generation population after the PSO algorithm finds the optimal individual and population. The GAPSO algorithm not only retains the original location transfer of the PSO algorithm but also integrates the powerful global search ability of the genetic algorithm, making the optimization process more efficient. It makes the search process more efficient and less likely to fall into a local optimum, and the obtained solution is thus more accurate.

In this paper, the trajectory relationship between the corresponding feature points is used to evaluate the similarity between the motion posture of the manta ray robot and the real manta ray. The trajectory coordinate points of manta ray feature points over time that can be represented by a vector $S_i=\left\{q_{i1},q_{i2},q_{i3},\dots ,q_{in}\right\}$, where $q_{ij}$ is the coordinate position of the manta ray’s pectoral fin feature points along the Z-axis over time.

The coordinate points of the trajectory of the pectoral fins of the manta ray robot can be represented by the vector $M_i=\left\{w_{i1},w_{i2},w_{i3},\dots ,w_{in}\right\}$, and $w_{ij}$ is the coordinate position of the pectoral fins of the manta ray robot along the Z-axis over time. Considering the flexible deformation of the pectoral fins of the manta ray robot when it flaps in the water, $w_{ij}$ can be calculated by Equation (16):

$$w_{ij}=l_i\sin ({\theta }_i+{\theta }_{tip}\sin (2\pi vt-\frac{\pi }{2}){(\frac{x}{b})}^2{(\frac{y}{c})}^2)$$

(16)

Since a simulation does not exist to compare two data points with different lengths of time series, the Euclidean distance algorithm is chosen to calculate the similarity value of the two models during the motion pose similarity simulation [34], and then the similarity function can be expressed by Equation (17). The fitness function established in this way is Equation (18):

$$d(S_i,M_i)=\sqrt{{\displaystyle \sum _{i=1}^n{(q_i-w_i)}^2}}$$

(17)

$$fitness=\min (d(S_i,M_i))$$

(18)

Due to the mechanical constraints of the manta ray robot, the maximum amplitude of pectoral fin flapping cannot exceed its mechanical limit. The stretching of the silicone pectoral fins cannot be too large, and to avoid damage to the silicone pectoral fins the phase difference between the two fins should be guaranteed to be within the safe stretching range, so the constraint is shown in Equations (19) and (20):

$$-w_{i\max }\le w_i\le w_{i\max }, w_{i\max }=60°$$

(19)

$$\Delta {\phi }_{ij}\le {35}^{\circ }$$

(20)

3. Results

3.1. Forward Swimming Posture Test of Manta Ray Robot

The forward swimming posture of a manta ray is often characterized by space asymmetric flapping and slip-flapping. The above multimodal motions involve the symmetric flapping of the pectoral fins on both sides. Therefore, this paper only needs to simulate one side of the pectoral fin flap to simulate the manta flap. Based on the forward travel parameters obtained in Equation (2) and the relationship between the displacement of the prototype fins along the Z-axis and the rudder output angle, the parameter initialization setting of the forward travel attitude is carried out in this paper.

$$\{\begin{array}{l}{R}_1=1.9R_l,R_2=R_3=R_l\\ R_4=1.9R_r,R_5=R_6=R_r,R_r=R_l=10°\\ v_i=v=0.4Hz\\ \Delta {\phi }_{12}=\Delta {\phi }_{45}=9.5°,\Delta {\phi }_{23}=\Delta {\phi }_{56}=42.8°,\Delta {\phi }_{14}=0°\\ \omega =2\end{array}$$

(21)

Firstly, the initial state of the CPG network is set as the space asymmetric coefficients ${\lambda }_1={\lambda }_2=0.67$, ${\lambda }_3=0.4$, and the time asymmetric coefficients ${\eta }_i=0.4$. The time taken from the starting state to the steady state is about 1 s, realizing the asymmetric cruise state forward in time and space. Six seconds later, the above bias state is realized as the gliding state, and 10 s later the acceleration frequency is 0.67, the amplitude increases to $1.3R_l$, and the space asymmetric coefficients ${\lambda }_1={\lambda }_2=0.5$. The whole process simulates the common manta ray forward cruise. The effectiveness of the CPG bionic control established in this paper is verified, as shown in Figure 7.

To verify the validity of the simulation, this paper conducts a forward swimming experiment with a manta ray robot. The manta ray robot was flapped with a space asymmetry factor of 0.7, the amplitude was set to 40°, and the frequency and phase difference were kept the same as the simulation settings. The motion frame diagram is shown in Figure 8 and the basic performance is shown in Figure 9.

To verify the effectiveness of the bionic motion control method constructed in this paper, the similarity calculation rule constructed in Section 2.3 is used to calculate the similarity of the fluttering forward swimming posture. In this paper, the similarity value between the manta ray robot and real manta ray forward swimming posture is 19.49% through the similarity evaluation of motion posture. The similarity value is still small, but the motion trend is consistent with that of the manta ray, which verifies the effectiveness of the bionic control.

3.2. Forward Swimming Posture Optimization Based on Similarity

The trajectory of the manta ray robot was extracted and the similarity to the manta ray motion was calculated to be only 19.49%. According to the manually calculated CPG control parameters that still need to be improved, this paper uses the GAPSO algorithm to optimize them while comparing the GA and PAO algorithms. Many parameters need to be optimized. To avoid the optimization algorithm falling into local optimization in the process of solving the objective function, this paper firstly optimized a single CPG unit to obtain the optimal solution range of each unit ${\eta }_i$, ${\lambda }_i$, and $R_i$. Then, coupling parameters are added to optimize the CPG network. The iteration speed and optimal values of the simulation are shown in Figure 10a. The optimized parameters are shown in

Table 2. Parameters after optimization.

${\lambda }_1$	${\eta }_1$	$R_1$	${\lambda }_2$	${\eta }_2$	$R_2$	${\lambda }_3$	${\eta }_3$	$R_3$	$\Delta {\phi }_{12}$	$\Delta {\phi }_{23}$
0.73	0.47	0.56	0.61	0.38	0.29	0.77	0.49	0.17	−0.77	−2

.

Comparing the three methods by the graph of iterative results, the GA has a slower convergence speed due to cross-variance and other operations, but it has a strong global search capability and a longer running period for the whole iteration. PSO has a strong search capability at the beginning of the iteration, converges faster, and reaches the optimum quickly but is prone to fall into a local optimum. The GAPSO algorithm shows a strong search capability and shows a global search capability at the end of the iteration. The search process is more efficient, less likely to fall into a local optimum, and the running speed is also improved over that of the GA algorithm. The three methods were run 10 times each, and the parameters with the best iterative results were selected and brought into the CPG network for simulation. The output waveforms of the optimized CPG network were compared with the corresponding trajectories of the real manta ray pectoral fin feature points, as shown in Figure 10b. The comparison shows that the waveform output from the CPG network is similar to the trajectory of the manta ray feature point, and the error occurs mainly at the transition between the upper and lower states (1.5 s).

In order to verify the effectiveness of the optimization simulation, the optimized CPG parameters were applied to the manta ray robot to perform forward swimming experiments and obtain the posture and performance data of the robot. The forward swimming motion sequence frame diagram is shown in Figure 11. The motion trajectory of the fin tip point of the manta ray robot is similar to the real manta ray trajectory, which verifies the effectiveness of CPG asymmetric characteristic network construction and GAPSO optimization, as shown in Figure 12a. The error is mainly generated at the point $P_1$. This point is the end of the first stage fin tip, and its flexible deformation is the largest at the beginning of the flutter. Thus, its error is the largest in the starting state. The specific motion pose similarity value is calculated; the similarity value of the improved motion pose at point $P_1$ is 80.96%, the similarity value at point $P_2$ is 85.33%, the similarity value at point $P_3$ is 99.30%, and the similarity value of the overall motion pose is 88.53%. Compared with the method of using manual settings of CPG parameters, the posture has been greatly improved, which verifies the effectiveness of this optimization, as shown in

Table 3. Comparison of similarity values before and after optimization.

	$P_1$	$P_2$	$P_3$	Average
Before optimization	19.57%	28.47%	10.43%	19.49%
After optimization	80.96%	85.33%	99.30%	88.53%
Improved	61.39%	56.86%	88.87%	69.04%

. Yet, the stability and maneuverability of the optimized robot were not significantly affected during forward swimming.

4. Conclusions

The manta ray robot’s kinematic posture is closer to that of the manta ray, which allows the manta ray robot to better express the bionic characteristics of a bionic robot, such as higher biophilic affinity and better stealth. To improve the bionic similarity of the manta ray robot, this paper establishes a bionic similarity evaluation rule for the manta ray robot. The traditional phase oscillator model is improved with space and time asymmetry, and a CPG network is constructed to realize the motion bionic of the manta ray robot. This paper proposes a GAPSO-based optimization method for CPG network parameters to improve the forward swimming posture imitation. Compared with the 19.49% similarity of the forward swimming posture of the manta ray robot under manual control parameter settings, the optimized similarity is improved to 88.53%, which verifies the effectiveness of the optimization method. This paper constructs a new and systematic high similarity bionic control method, which can effectively improve the bionic similarity of the manta ray robot and is of great significance for the future application of the manta ray robot in marine museums’ stealthy approach reconnaissance and other fields.

In the future, we will consider the characteristics of underwater robots to establish a more comprehensive similarity evaluation rule for bionic underwater robots in order to facilitate the creation of robots with the same or better performance as living things.