Multi-Objective Multidisciplinary Design Optimization of a Robotic Fish System

Abstract

Biomimetic robotic fish systems have attracted huge attention due to the advantages of flexibility and adaptability. They are typically complex systems that involve many disciplines. The design of robotic fish is a multi-objective multidisciplinary design optimization problem. However, the research on the design optimization of robotic fish is rare. In this paper, by combining an efficient multidisciplinary design optimization approach and a novel multi-objective optimization algorithm, a multi-objective multidisciplinary design optimization (MMDO) strategy named IDF-DMOEOA is proposed for the conceptual design of a three-joint robotic fish system. In the proposed IDF-DMOEOA strategy, the individual discipline feasible (IDF) approach is adopted. A novel multi-objective optimization algorithm, disruption-based multi-objective equilibrium optimization algorithm (DMOEOA), is utilized as the optimizer. The proposed MMDO strategy is first applied to the design optimization of the robotic fish system, and the robotic fish system is decomposed into four disciplines: hydrodynamics, propulsion, weight and equilibrium, and energy. The computational fluid dynamics (CFD) method is employed to predict the robotic fish’s hydrodynamics characteristics, and the backpropagation neural network is adopted as the surrogate model to reduce the CFD method’s computational expense. The optimization results indicate that the optimized robotic fish shows better performance than the initial design, proving the proposed IDF-DMOEOA strategy’s effectiveness.

1. Introduction

Most conventional autonomous underwater vehicles (AUVs) use propellers for propulsion [1], but this propulsion mode has some disadvantages, such as poor concealment and low efficiency. Due to the evolution over millions of years, fish can perform fast and very efficient swimming motions [2]. Moreover, appealing morphological properties for moving through the water with prominent speed and maneuverability indicate that the biological features and locomotion abilities of fish can be used to design biomimetic robotic fish. As a new conceptual AUV, the robotic fish has the advantages of flexibility and adaptability, and it can be defined as an intelligent underwater propulsion system, relying on oscillatory and undulatory motions to move through the water [3,4]. A robotic fish is typically a complex system that involves many disciplines. Its design is a multi-objective and multidisciplinary problem.

Multidisciplinary design optimization (MDO) provides an effective way for solving conceptual design problems of complex systems [5]. Using MDO approaches, complex systems can be decomposed into several disciplines (or subsystems), and the interactions between disciplines are well considered [6]. Therefore, the optimal design scheme can be obtained, and the design cost can also be reduced to a great extent. Many well-known MDO approaches have been proposed in the past few decades. These MDO approaches can be classified into two categories: Monolithic approaches and distributed approaches. Monolithic approaches, including multidisciplinary feasible (MDF), all-at-once (AAO) and individual discipline feasible (IDF), only have one optimizer in the system layer for optimization, and there is only analysis in each discipline [7]. Distributed approaches such as collaborative optimization (CO) [8], concurrent subspace optimization (CSSO) [9] and bilevel integrated system synthesis (BLISS) [10] have an optimizer at the system layer and the discipline layer for optimization. MDO has been increasingly applied to solve the conceptual design problems of underwater robots in recent years. Luo and Lyu [11] applied CO with particle swarm optimization algorithm (PSO) to the hydrodynamics performance optimization of an autonomous underwater vehicle. Chen et al. [12] divided the optimal design of an AUV into three disciplines, including control discipline, hydrodynamics discipline, and power and mass distribution discipline, and the MDF method was adopted for the conceptual design of the AUV. By combining the PSO algorithm and the MDF approach, Bidoki et al. [13] presented a PSO-MDF strategy and applied it to the multidisciplinary design optimization of an AUV. Zhang et al. [14] applied a concurrent subspace design (CSD) approach to establish the MDO architecture of an intelligent ocean exploration underwater vehicle, and then they applied the surrogate model to reduce the computational cost of the MDO problem.

The above applications focused on solving single-objective optimization problems (SOPs). However, there are usually several conflicting optimization objectives in the design optimization process of underwater robots; these problems are multi-objective optimization problems (MOPs). As for MOPs, the performance improvement of one objective often leads to the performance reduction of other objectives, and it is difficult to obtain the optimal solution for multiple objectives simultaneously. Thus, there is a set of compromise solutions, which is called the Pareto optimal solution set [15]. In the past few decades, many multi-objective optimization algorithms have been proposed to solve MOPs, including non-dominated sorting genetic algorithm version 2 (NSGA-II) [16], region-based selection in evolutionary multi-objective optimization (PESA-II) [17], multi-objective particle swarm optimization (MOPSO) [18], multi-objective grey wolf optimization algorithm (MOGWO) [19], multi-objective whale optimization (MOWOA) [20], multi-objective hyper-heuristic algorithm based on adaptive epsilon-greedy selection (HH_EG) [21] and disruption-based multi-objective equilibrium optimization algorithm (DMOEOA) [22], etc. Wang et al. [23] developed a multi-objective multidisciplinary design optimization (MMDO) strategy for the shape design optimization of an AUV. This strategy utilized a CSD approach as the MDO architecture, and the unified-objective method was utilized to change the multi-objective optimization problem into a single-objective optimization problem. Liu et al. [24] carried out the research on the MMDO of a heavier-than-water underwater vehicle by using the computational fluid dynamics (CFD) method and the surrogate model, and NSGA-II was adopted as the optimizer. This paper presented a multi-objective multidisciplinary design optimization strategy called IDF-DMOEOA to solve the conceptual design problem of a three-joint robotic fish system. The robotic fish system is divided into four disciplines, including hydrodynamics discipline, propulsion discipline, weight and equilibrium discipline and energy discipline. The IDF approach requires fewer function calls and additional information than other MDO approaches in the formulation process [25]. Therefore, the IDF approach and a novel and efficient multi-objective optimization algorithm named DMOEOA in our previous work [22] are integrated into the proposed IDF-DMOEOA strategy. The hydrodynamics characteristics analysis of the robotic fish is established by using the CFD method. To reduce the computational cost of the CFD method in the optimization process, a backpropagation neural network (BPNN), which has a strong ability of nonlinear mapping, is adopted as the surrogate model in the hydrodynamics discipline [26].

In this paper, the discipline analysis is presented in Section 2, and the multidisciplinary design optimization architecture of the robotic fish’s conceptual design is established in Section 3. The adopted multi-objective optimization algorithm, DMOEOA, is introduced in Section 4. Then, the optimization results of the proposed MMDO problem and analysis are given in Section 5. Finally, Section 6 provides the concluding remarks.

2. Discipline Analysis

The robotic fish system is decomposed into four disciplines: hydrodynamics discipline, propulsion discipline, weight and equilibrium discipline, and energy discipline. Discipline analysis models of the above four disciplines are established in this section.

2.1. Hydrodynamics Discipline

Fluid resistance plays an important role in underwater robots, and it is necessary to design the shape of the robotic fish with the least resistance. Thus, the hydrodynamics analysis model is presented for fluid resistance estimation.

2.1.1. Hydrodynamics Analysis Model

The hydrodynamics analysis model consists of three parts: parametric modeling, mesh generation and CFD numerical simulation. Three kinds of software are utilized to construct the analysis model, including Unigraphics NX 10.0, Ansys ICEM 20.0 CFD and Ansys FLUENT 20.0. The whole analysis process will be repeated with the change of input parameters. The flowchart of the hydrodynamics analysis model is shown in Figure 1.

2.1.2. Parametric Modeling of the Hull Shape

As shown in Figure 2, the hull of the robotic fish is shaped like a torpedo. It consists of the fore-body, parallel middle body and after-body. $L_a$ and $L_f$ represent the length of the after-body and fore-body, respectively, and R is the maximum radius of the hull. In addition, $L_{ct}$ refers to the length of cut off from the after-body section, and $L_p$ indicates the length of parallel middle body. The line shape of the hull is defined by the Glanville lines [27]. The formulas of the after-body section and fore-body section are described by Equations (1) and (2), respectively.

$$z_a=q{a}_{1}x{(x-1)}^4+1-{(x-1)}^4(4x+1)+\frac{q{a}_2{x}^2{(x-1)}^3}{6},(0\le x\le 1)$$

(1)

$$\begin{array}{cc}\hfill & z_f=|x^4-1|+q{f}_1\frac{(1-x)(5{(1-x)}^3-21{(1-x)}^2-35x)}{16}+q{f}_1\sqrt{(1-x)}-\frac{q{f}_2(1-x)x^3}{6},\hfill \\ \\ & (0\le x\le 1)\hfill \end{array}$$

(2)

where $q{a}_1$ and $q{a}_2$ are shape factors of the after-body section. $q{f}_1$ and $q{f}_2$ represent shape factors of the fore-body section. The Glanville lines are dimensionless geometry lines, and they must be expressed in physical forms. The physical line shape of the hull is described below.

$$Z_a=R{z}_a\left(\frac{x}{L_{ct}+L_a}\right),(0\le x\le L_{ct}+L_a)$$

(3)

$$Z_p=R,(L_a+L_{ct}\le x\le L_a+L_{ct}+L_p)$$

(4)

$$Z_f=R{z}_f\left(\frac{x-L_a-L_p-L_{ct}}{L_f}\right),(L_a+L_{ct}+L_p\le x\le L_a+L_{ct}+L_p+L_f)$$

(5)

Initial values of parameters that will influence the shape of the robotic fish are listed in

Table 1. Initial values of the shape parameters.

Shape Parameter	Value
$L_f$	20 (cm)
$L_p$	27 (cm)
$L_a$	6 (cm)
$L_{ct}$	0.2 (cm)
R	5 (cm)
$q{f}_1$	2
$q{f}_2$	1
$q{a}_1$	3
$q{a}_2$	12

.

UG 10.0 is applied to read the shape parameters of the hull and create the model, and then the generated parasolid file is exported to the mesh generation part. In this work, four shape factors, including $q{a}_1$, $q{a}_2$, $q{f}_1$ and $q{f}_2$, and the length of the parallel middle body $L_p$ are input variables in the hydrodynamics discipline. The objective of this discipline is to design the line shape of the hull with the least fluid resistance $F_{drag}$. The input and output variables of the hydrodynamics discipline are summarized in

Table 2. Summary of the hydrodynamics discipline.

Input Variable	Output Variable
$L_p$ (24 cm $\le L_p\le$ 30 cm)	$F_{drag}$
$q{f}_1$$(2\le q{f}_1\le 4)$
$q{f}_2$$(0\le q{f}_2\le 4)$
$q{a}_1$$(2\le q{a}_1\le 4)$
$q{a}_2$$(10\le q{a}_2\le 20)$

.

2.1.3. Mesh Generation

In this part, ICEM 20.0 is employed to generate meshes, to balance the accuracy and calculation efficiency of the CFD method. The computational domain shown in Figure 3 is selected, and it consists of a cylinder and a semi-sphere. The flank area of the cylinder and the semi-spherical are set as the velocity inlet (5 m/s), and the bottom surface of the cylinder is set as the pressure outlet with reference pressure 0. The hull of the robotic fish is defined as the wall. The distance between the vertex of the hull and the vertex of the hemispherical sphere is equal to the length of the hull. The distance between the hull’s end and the outlet boundary is twice the length of the hull. Because of the regular geometry of the hull, structured meshes are applied to generate elements and nodes in the computational domain. The mesh arrangement around the hull is shown in Figure 4.

2.1.4. CFD Numerical Simulation

Fluent 20.0 is utilized as the solver, the CFD numerical simulation adopts RANS equations as the control equations and the finite volume method is employed to discrete the control equations. Standard $k-\epsilon$ is adopted as the turbulence model, and the standard wall function method is applied to treat areas near the wall. The wall in the computational domain is set as a no-slip wall. The turbulent viscosity ratio and the intensity for the velocity inlet and pressure outlet are set to 2 and 2%, respectively, and the standard discretization scheme is used for pressure [11]. The SIMPLEC algorithm is utilized as the solution method. The second-order upwind scheme is applied for turbulence kinetic energy, momentum and turbulence dissipation rate, and under-relaxation factors are set to default values.

2.1.5. CFD Grid Convergence Study

The accuracy of the CFD method has a significant impact on its practical application in engineering fields. In this work, the grid convergence analysis is carried out to ensure that the numerical simulation results are not affected by the grid size. Three different mesh sizes are adopted for the computational domain, including 942,208 cells, 1,492,608 cells and 2,090,538 cells. The boundary layer has a significant impact on the accuracy of the CFD method, so the $y^{+}$ of the first layer near the wall is the criteria for the mesh quality. The height of the first layer near the wall can be calculated by the estimation formula [11]:

$$y^{+}=\frac{0.172\Delta y{R}_e^{0.9}}{L^c}$$

(6)

where the desired $y^{+}$ value is set as 45. $R_e$ represents the Reynolds number, and $L^c$ is the characteristic length of the hull. In this simulation, the $R_e$ is $2.65\times 10^6$, and the density of the water $\rho$ is set as 998.2 kg/m³. In Figure 5, the x axis refers to the position of the wall surface of the robotic fish hull along the x axis direction, we can observe that $y^{+}$ values of the mesh with 942,208 cells are located in the range 30–60, so the mesh quality achieves the requirement [11].

Table 3. Comparison of $F_{drag}$ values using different numbers of mesh.

Number of Mesh	${\mathit{F}}_{\mathit{drag}}$$\left(\mathit{N}\right)$
942,208	12.95503
1,492,608	13.01365
2,090,538	13.03643

shows the CFD simulation results of fluid resistance $F_{drag}$ of the hull using different numbers of mesh. We can observe that the $F_{drag}$ obtained by the finest mesh differs by 0.624% from that of the coarsest mesh, which proves the CFD numerical simulation results are accurate and stable, and the mesh size with 942,208 cells is applied for the optimization to reduce computational costs. The surface pressure field of the hull with 942,208 cells is shown in Figure 6.

2.1.6. Surrogate Model for Hydrodynamics Analysis

As a high-fidelity analysis method, CFD simulation is time-consuming. If the CFD method is integrated into the hydrodynamics analysis model in the optimization process directly, it may lead to a large amount of calculation and low efficiency. Surrogate models provide an effective way to solve design optimization problems with a better balance between accuracy and efficiency. Many surrogate models have been applied successfully in engineering optimization fields, including the response surface method (RSM), Kriging model and artificial neural networks (ANNs). In this work, a backpropagation neural network (BPNN) [26] is utilized as the surrogate model to reduce the computational cost of the CFD simulaton. The process of constructing the surrogate model shown in Figure 7 consists of four steps.

• Design of experiment (DOE) can be applied to obtain sufficient input–output sample points for the construction of the surrogate model [5]. In this work, an optimal Latin hypercube design is adopted as the DOE method, and 200 sample points are collected.
• In this step, BPNN is employed as the surrogate model, and 70% of the sample points are selected randomly to train the surrogate model, and the number of hidden layer neurons is set to 20 in the BPNN model.
• In this step, 15% of the sample points are employed to validate the BPNN model, and 15% of the sample points are used to test the BPNN model.
• If the surrogate model meets the accuracy requirement, the BPNN model is exported to the hydrodynamics discipline model, or return to step 2 to retrain the BPNN model.

The accuracy of the surrogate model is important. The correlation coefficient R value measures the linear correlation between observed values and predicted values. An R value of 0 indicates no linear correlation relationship, and 1 indicates a full linear correlation. As shown in Figure 8, correlation coefficient values of the training set, validation set and test set are all greater than 0.9, and the coefficient of determination $R^2$ of the BPNN model is 0.9024. The mean squared error MSE of the BPNN model is listed in

Table 4. MSE values of the BPNN model.

Evaluation Index	Value
MSE for training set	0.0094
MSE for validation set	0.0247
MSE for test set	0.0234

. From Figure 8 and Table 4 we can observe that the BPNN model has sufficient precision for hydrodynamics discipline analysis.

2.2. Weight and Equilibrium Discipline

Weight distribution and static equilibrium are two factors that need to be considered in the conceptual design of robotic fish. As shown in Figure 9, the three-joint biomimetic robotic fish consists of four components, including the rigid fish head $fh$, the flexible fish body $fb$, the after-body $af$ and the caudal fin $sw$. The rigid fish head component includes the battery b, the control panel e and a junction plate $jp$. Three servo motors $s{m}_i(i=1,2,3)$ are utilized as joints in the flexible fish body to mimic fish-like swimming. There are also counterweights $P_i(i=1,2,3)$ inside the flexible fish body. The caudal fin and fish body are connected by the after-body. $l_1$ indicates the length of link $o_1{o}_2$, $l_2$ refers to the length of link $o_2{o}_3$ and $l_3$ represents the distance between joint $o_3$ and the after-body. The sum of lengths of $l_1$, $l_2$ and $l_3$ is equal to the length of the flexible fish body $L_p$. As for the weight and equilibrium discipline, the total weight and buoyancy should be balanced, and the trim angle ${\theta }_t$ also needs to be considered. ${\theta }_t$ can be calculated by the following equations:

$$r_{cm}=\frac{m_{fh}{r}_{fhm}+m_{fb}{r}_{fbm}+m_{af}{r}_{afm}+m_{sw}{r}_{swm}+m_h{r}_{hm}}{m_{fh}+m_{fb}+m_{af}+m_{sw}+m_h}$$

(7)

$$r_{cb}=\frac{V_{fh}{r}_{fhb}+V_{fb}{r}_{fbb}+V_{af}{r}_{afb}+V_{sw}{r}_{swb}+V_h{r}_{hb}}{V_{fh}+V_{fb}+V_{af}+V_{sw}+V_h}$$

(8)

$${\theta }_t=arctan\left(\frac{x_{cb}-x_{cm}}{z_{cb}-z_{cm}}\right)$$

(9)

where $r_{cm}={\left[\begin{array}{ccc}{x}_{cm}& y_{cm}& z_{cm}\end{array}\right]}^T$ and $r_{cb}={\left[\begin{array}{ccc}{x}_{cb}& y_{cb}& z_{cb}\end{array}\right]}^T$ represent the barycenter and buoyant center of the robotic fish, respectively, $m_{fh},m_{fb},m_{af},m_{sw}$ represent the masses of the four components of the robotic fish, respectively, and $V_{fh},V_{fb},V_{af},V_{sw}$ are the volumes of the four components, respectively. $r_{fhm},r_{fbm},r_{afm},r_{swm}$ indicate the barycenters of those four components, respectively, and $r_{fhb},r_{fbb},r_{afb},r_{swb}$ are the buoyant centers of the four components, respectively. $m_h,V_h,r_{hm},r_{hb}$ represent the mass, volume, barycenter and buoyant center of the hull, respectively. This discipline is constructed to compute the mass of each component and to balance the buoyancy and weight of the robotic fish. Three parameters, including $l_1$, $l_2$ and $l_3$ are input variables in the weight and equilibrium discipline. The input and output variables of this discipline are summarized in

Table 5. Summary of weight and equilibrium discipline.

Input Variable	Output Variable
$l_1$ (10 cm $\le l_1\le$ 12 cm)	$m_{fb}$
$l_2$ (10 cm $\le l_2\le$ 12 cm)	$V_{fb}$
$l_3$ (4 cm $\le l_3\le$ 8 cm)	$m_h$
	$V_h$
	${\theta }_t$

.

2.3. Propulsion Discipline

Most conventional AUVs use propellers as the propulsion mode, which may lead to low efficiency and poor concealment. In contrast, the multi-joint robotic fish use multiple joints and links to mimic fish-like swimming, in which rotating joints are actuated by servo motors. The propulsion velocity of robotic fish is a key indicator of swimming performance, and the dynamic model is employed to analyze the velocity of the three-joint robotic fish.

2.3.1. Determination of Coordinate Frames

Several methods have been employed to establish the dynamic model of multi-joint robotic fish successfully, such as the Lagrange method [28], the Kane method [29] and the Newton–Euler method [30]. Due to the advantages of simplicity and easy implementation, the Newton–Euler method is applied to establish the dynamic model of the multi-joint robotic fish in this section.

To analyze the swimming dynamics of the multi-joint robotic fish, the coordinate frames are shown in Figure 10. Assume that the multi-joint robotic fish consists of $nl$ links, $L_i(i=0,1,2,\dots ,nl-1)$ represents the length of the ith link, the center of mass of each link that coincides with its geometrical center is denoted as $C_i$. $m_i(i=0,1,2,\dots ,nl-1)$ represents the mass of the ith link. Links are connected by joints (i.e., output shafts of servo motors). The moving coordinate system $O_i-X_i{Y}_i{Z}_i$ is attached to the ith link, and the direction of the $O_i{X}_i$ is parallel to the axis of the ith link. $O_i$ is fixed in the starting point of the ith link. The angle between the $(i-1)$th link and the ith link is ${\theta }_i(i=1,2,\dots ,nl-1)$; in particular, ${\theta }_0$ indicates the angle between the 0th link and the X axis of the world coordinate system $O_w-X_w{Y}_w{Z}_w$ (WCS). $R_i^{i+1}$ is the transform matrix from $O_i-X_i{Y}_i{Z}_i$ to $O_{i+1}-X_{i+1}{Y}_{i+1}{Z}_{i+1}$. $R_w^0$ indicates the transform matrix from $O_w-X_w{Y}_w{Z}_w$ to $O_0-X_0{Y}_0{Z}_0$. $R_i^{i+1}$ and $R_w^0$ can be expressed as follows:

$$R_i^{i+1}=\left[\begin{array}{ccc}cos{\theta }_{i+1}& sin{\theta }_{i+1}& 0\\ -\sin {\theta }_{i+1}& \cos {\theta }_{i+1}& 0\\ 0& 0& 1\end{array}\right](i=0,1,2,\dots ,nl-1)$$

(10)

$$R_w^0=\left[\begin{array}{ccc}cos{\theta }_0& sin{\theta }_0& 0\\ -\sin {\theta }_0& \cos {\theta }_0& 0\\ 0& 0& 1\end{array}\right]$$

(11)

2.3.2. Kinematic Analysis

Let the coordinate of the joint $O_0$ in WCS be ${(x,y,0)}^T$. The position vector $r_{oi}^i(i=0,1,2,\dots ,nl-1)$ of joint $O_i$ relative to the origin of the world coordinate system $O_w$ can be expressed in its moving coordinate system as:

$$r_{o0}^0=R_w^0{\left[\begin{array}{ccc}x& y& 0\end{array}\right]}^T$$

(12)

$$r_{oi}^i=R_{i-1}^i(r_{oi-1}^{i-1}+{\left[\begin{array}{ccc}{L}_{i-1}& 0& 0\end{array}\right]}^T)$$

(13)

Therefore, the position vector $r_{ci}^i(i=0,1,2,\dots ,nl-1)$ of the center of mass of link i relative to the origin of WCS can be described in its moving coordinate system as:

$$r_{ci}^i=r_{oi}^i+{\left[\begin{array}{ccc}\frac{1}{2}{L}_i& 0& 0\end{array}\right]}^T$$

(14)

The angular velocity $w_i$ and the angular acceleration ${\alpha }_i$ of the ith link relative to the origin of WCS can be expressed in $O_i-X_i{Y}_i{Z}_i$ as:

$$w_i^i={\left[\begin{array}{ccc}0& 0& {\sum }_0^{i}d{\theta }_i\end{array}\right]}^T$$

(15)

$${\alpha }_i^i=\frac{d{w}_i^i}{dt}$$

(16)

The velocity $v_i$ and the acceleration $a_i$ of the center of mass of link i relative to the origin of WCS in $O_i-X_i{Y}_i{Z}_i$ can be expressed as:

$$v_{ci}^i=\frac{d{r}_{ci}^i}{dt}+w_i^i\times r_{ci}^i$$

(17)

$$a_{ci}^i=\frac{d{v}_{ci}^i}{dt}+w_i^i\times v_{ci}^i$$

(18)

2.3.3. Dynamic Analysis

The fluid resistance on a link depends on the relative velocity of the link with respect to the fluid. The relative velocity $v_{ri}^i$ and the relative acceleration $a_{ri}^i$ are as follows:

$$v_{ri}^i=v_{ci}^i-v_{fi}^i$$

(19)

$$a_{ri}^i=a_{ci}^i-a_{fi}^i$$

(20)

where $v_{fi}^i$ and $a_{fi}^i$ represent the velocity and acceleration of the fluid in $O_i-X_i{Y}_i{Z}_i$, respectively.

Then, the corresponding moment of fluid resistance can be expressed as follows:

$$T_{di}={\left[\begin{array}{ccc}\frac{1}{2}{L}_i& 0& 0\end{array}\right]}^T\times F_{di}$$

(21)

Joints of the robotic fish are employed to connect links, and there are forces and moments generated by the interaction between links. Suppose that the interaction force and moment of the ith joint are ${}^j{F}_i$ and ${}^j{T}_i$, respectively. Based on Newton–Euler method, the dynamic equations of the ith link can be formulated as follows:

$${}^j{F}_i+F_{di}=m_i{a}_{ci}^i$$

(22)

The fluid drag $F_{di}$ of the ith link can be calculated by the following equation [31]:

$$F_{di}=\left[\begin{array}{c}-0.5\rho c_x|v_{rix}^i|v_{rix}^i{s}_{ix}\\ -0.5\rho c_y|v_{riy}^i|v_{riy}^i{s}_{iy}\\ 0\end{array}\right]$$

(23)

$${}^j{T}_i+T_{di}=J_{ci}·{\dot{w}}_i^i+w_i^i\times J_{ci}·w_i^i$$

(24)

where $J_{ci}$ indicates the rotational inertia of the ith link. The motion rule of the ith joint is as follows:

$${\theta }_i=A_{i}sin({\omega }_{i}t+\phi )$$

(25)

where $A_i$ and ${\omega }_i$ represent the angle of amplitude and angular velocity of the ith joint, respectively. $\phi$ indicates the phase of joint motion. By solving the above equations, we are able to analyze the velocity of the rigid head of the robotic fish. As shown in Figure 9, the length of the ith $(i=0,1,2,3)$ link can be described as follows:

$$L_0=L_f$$

(26)

$$L_1=l_1$$

(27)

$$L_2=l_2$$

(28)

$$L_3=l_3+L_a+L_{sw}$$

(29)

Through the dynamic model, we can establish the propulsion discipline analysis framework of the robotic fish. The main objective of the propulsion discipline is to analyze the forward velocity of the three-joint robotic fish v. The input and output variables in this discipline are listed in

Table 6. Summary of the proplusion discipline.

Input Variable	Output Variable
${\omega }_i(i=1,2,3)(\pi \le {\omega }_i\le 2\pi )$	v
$A_i(i=1,2,3)(\frac{\pi }{6}\le A_i\le \frac{\pi }{3})$
$\phi (\frac{\pi }{6}\le \phi \le \frac{\pi }{2})$
$m_i(i=0,1,2,3)$
$l_1$ (10 cm $\le l_1\le$ 12 cm)
$l_2$ (10 cm $\le l_2\le$ 12 cm)
$l_3$ (4 cm $\le l_3\le$ 8 cm)

.

2.4. Energy Discipline

A battery is utilized to provide energy for the robotic fish. Among various kinds of batteries, in this work, a lithium battery is adopted due to its long service life and large capacity density. The task of the energy discipline is to estimate the endurance of the robotic fish based on the total energy provided by the battery. The total energy of the battery can be calculated as follows [24]:

$$E_b=\frac{m_b·{\rho }_{mb}·(1-SMF)}{1000}$$

(30)

where $m_b$ and ${\rho }_{mb}$ represent the mass and mass energy density of the lithium battery. $SMF$ refers to the spare margin factor, which is set to 5% according to [24]. To calculate the power of each motor in the robotic fish, we need to estimate the torque of each motor, as shown in Figure 11, and each joint of the robotic fish is driven by a motor. As for the motor $s{m}_1$, assume that the barycenter of its load is $g_1$, the distance between $g_1$ and joint $o_1$ is $r_1$, then, the motion rule of $g_1$ can be calculated as follows:

$$y_{g1}=r_1{A}_{1}sin\left({\omega }_{1}t\right)$$

(31)

$$v_{g1}={\omega }_1{r}_1{A}_{1}cos\left({\omega }_{1}t\right)$$

(32)

$$v_{g1max}={\omega }_1{r}_1{A}_1$$

(33)

where ${\omega }_1$ is the angular velocity and $A_1$ is the amplitude of the motion of $g_1$. $y_{g1}$ and $v_{g1}$ represent the location and velocity of $g_1$, respectively. $v_{g1max}$ represents the maximum velocity of $g_1$. Then the drag $D_{g1}$ and drag moment $T_{d1}$ applied to the load 1 can be estimated as follows:

$$D_{g1}=\frac{\rho C_{d1}{S}_{c1}{v}_{{}_{g1max}}^2}{2}$$

(34)

$$T_{d1}=D_{g1}{r}_1$$

(35)

where $\rho$ is the density of the water, $C_{d1}$ and $S_{c1}$ represent the drag coefficient and wetted surface of load 1, respectively. Therefore, the power of motor 1 $P{P}_1$ can be calculated as follows:

$$P{P}_1=T_{d1}{\omega }_1$$

(36)

In the same way, the power of other motors can be estimated through the above equations. The endurance of the robotic fish can be calculated as follows:

$$H_{en}=\frac{E_b}{{\displaystyle \sum _{i=1}^3}P{P}_i}$$

(37)

The main objective of energy discipline is to estimate the endurance of the robotic fish $H_{en}$. The calculated torque of $s{m}_1$ should not exceed the maximum torque $T_{dmax}$. The input variables and output variables of the energy discipline are summarized in

Table 7. Summary of the energy discipline.

Input Variable	Output Variable
${\omega }_i(i=1,2,3)(\pi \le {\omega }_i\le 2\pi )$	$H_{en}$
$A_i(i=1,2,3)(\frac{\pi }{6}\le A_i\le \frac{\pi }{3})$	$T_{d1}$
$l_1$ (10 cm $\le l_1\le$ 12 cm)
$l_2$ (10 cm $\le l_2\le$ 12 cm)
$l_3$ (4 cm $\le l_3\le$ 8 cm)

.

2.5. MMDO Model of the Robotic Fish

There are three optimization objectives in the conceptual design of the robotic fish; the first objective is to design an optimal shape of the robotic fish with the least fluid resistance, and the other two objectives are to maximize the forward velocity and endurance of the robotic fish. Considering the main optimization objectives of the robotic fish, the constraints on mechanical limits and size requirements, the values of design variables are set in a reasonable range. The MMDO model of the robotic fish can be specified as follows:

$$\begin{array}{c}MinF=\{F_{drag},-v,-H_{en}\}\hfill \\ \\ w.r.t\hfill \\ 10\le x_1=l_1\le 12,10\le x_2=l_2\le 12\hfill \\ \\ 4\le x_3=l_3\le 8,2\le x_4=q{f}_1\le 4\hfill \\ \\ 0\le x_5=q{f}_2\le 4,2\le x_6=q{a}_1\le 4\hfill \\ \\ 10\le x_7=q{a}_2\le 20,\pi \le x_8={\omega }_1\le 2\pi \hfill \\ \\ \pi \le x_9={\omega }_2\le 2\pi ,\pi \le x_{10}={\omega }_3\le 2\pi \hfill \\ \\ \frac{\pi }{6}\le x_{11}=\phi \le \frac{\pi }{2},\frac{\pi }{6}\le x_{12}=A_1\le \frac{\pi }{3}\hfill \\ \\ \frac{\pi }{6}\le x_{13}=A_2\le \frac{\pi }{3},\frac{\pi }{6}\le x_{14}=A_3\le \frac{\pi }{3}\hfill \\ \\ s.t.\hfill \\ g_1=|{\theta }_t|-10\le 0\hfill \\ \\ g_2=T_{d1}-3.5\le 0\hfill \end{array}$$

(38)

Based on the above MMDO model and analysis of the four disciplines, coupling relationships between these four disciplines are shown in Figure 12. We can observe that $l_i,{\omega }_i,A_i(i=1,2,3)$ are design variables shared by different disciplines. The objective v is obtained not only by the analysis of propulsion discipline but also by decoupling the interactions between the weight and equilibrium discipline and the propulsion discipline. Traditional optimization algorithms are not suitable for solving this coupled MMDO problem because those methods do not take the coupling relationships into consideration. Therefore, the MMDO strategy is needed to solve the above optimization problem. The flowchart of the global proposed approach for the MMDO of the robotic fish is shown in Figure 13.

3. Multidisciplinary Design Optimization

Because the IDF approach requires less additional information and fewer function calls than other multidisciplinary design optimization approaches, it is adopted in this work.

3.1. Individual Discipline Feasible Approach

The optimization model of the IDF approach is described as follows [32]:

$$\begin{array}{c}MinF(X,y(X,y^t))\hfill \\ \\ w.r.t\hfill \\ X,y^t\hfill \\ \\ s.t.\hfill \\ c_i(X_0,X_i,y_i(X_0,X_i,y_{j\ne i}^t))\le 0\hfill \\ \\ c_i^c=y_i^t-y_i(X_0,X_i,y_{j\ne i}^t)=0\hfill \end{array}$$

(39)

where $X={[X_0,X_1,\dots ,X_{nd}]}^T$ refers to the vector of design variables, $nd$ is the number of disciplines. $X_0$ and $X_i(i=1,2,\dots ,nd)$ represent the shared variables and design variables local to the ith discipline, respectively. The coupling design variable supplied by the ith discipline is denoted as $y_i$, and $y_i^t$ is the copy of $y_i$. ${\overline{y}}_i$ represents the state variable of discipline i. $c_i$ and $c_i^c$ represent the constraint and the consistency constraint of discipline i, respectively. Coupling variable copies and consistency constraints are introduced into IDF architecture to allow the discipline analysis to run in parallel. The extended design structure matrix (XDSM) of the IDF architecture [33] is shown in Figure 14. The IDF architecture enables the discipline analysis to run independently, which can reduce the computation time, accelerate the computation speed and improve the computation efficiency.

3.2. IDF Architecture of the MMDO Problem

In this work, the IDF approach is applied to solve the MMDO of the robotic fish. The MMDO framework of the robotic fish based on IDF architecture is shown in Figure 15. The expressions of variables and constraints contained in the proposed MMDO problem in IDF architecture are listed in

Table 8. Expressions of notations contained in the proposed MMDO problem in IDF architecture.

Variable and Constraint	Notation in IDF Architecture
$l_i,{\omega }_i,A_i(i=1,2,3)$	$X_0$
$q{f}_1,q{f}_2,q{a}_1,q{a}_2$	$X_1$
$\phi$	$X_3$
$m_i(i=0,1,2,3)$	$y_2$
$F_{drag}$	${\overline{y}}_1$
${\theta }_t$	${\overline{y}}_2$
v	${\overline{y}}_3$
$H_{en},T_{d1}$	${\overline{y}}_4$
$g_i(i=1,2)$	c

. Using IDF architecture, the optimization process can be performed in parallel.

4. Multi-Objective Optimization Algorithm

In this work, a novel and efficient multi-objective optimization algorithm named disruption-based multi-objective equilibrium optimization algorithm (DMOEOA) is adopted as the optimizer in the proposed MMDO strategy. DMOEOA was proposed in our previous work [22]. It is an effective optimization algorithm that shows superiority in solving multi-objective optimization problems with a better balance between distribution and convergence of obtained Pareto solutions.

4.1. Grid Mechanism

In DMOEOA, a grid mechanism is applied to record the distribution and convergence of obtained solutions. Some definitions used in this work are introduced as follows [34]:

Definition 1.

(Grid boundary): $max{f}_i\left(x\right)$ and $min{f}_i\left(x\right)$ represent maximum and minimum values of the ith objective, respectively, lower limit $l_{wi}$ and upper limit $u_{pi}$ of grids in the ith objective space can be calculated as follows:

$$l_{wi}=min{f}_i\left(x\right)-\frac{max{f}_i\left(x\right)-min{f}_i\left(x\right)}{2Gd}$$

(40)

$$u_{pi}=max{f}_i\left(x\right)+\frac{max{f}_i\left(x\right)-min{f}_i\left(x\right)}{2Gd}$$

(41)

where $Gd$ refers to the number of grids in each dimension of the objective space.

Definition 2.

(Grid location): An individual’s grid location can be defined as:

$$G{L}_i\left(x\right)=\left[\frac{f_i\left(x\right)-l_{wi}}{d_i}\right]$$

(42)

where $d_i$ indicates the width of grids in the ith objective, $[·]$ refers to the function of rounding up. For instance, in Figure 16, grid locations of individuals B and C are $(1,4)$ and $(2,3)$, respectively.

Definition 3.

(Grid ranking): An individual’s grid ranking ($GR$) is defined as follows:

$$GR\left(x\right)=\sum _{i=1}^{k}G{L}_i\left(x\right)$$

(43)

where k refers to the number of objectives, as shown in Figure 16, solution A’s grid ranking is four and solution D’s grid ranking is six, which means that solution A is closer to the Pareto optimal front than solution D.

Definition 4.

(Grid coordinate point distance): Grid coordinate point distance $\left(GCPD\right)$ is defined as the normalized Euclidean distance between the individual and minimum boundary point in its grid:

$$GCPD\left(x\right)=\sqrt{\sum _{i=1}^k{\left[\frac{(f_i\left(x\right)-(l_{wi}+(G{L}_i\left(x\right)-1)·d_i))}{d_i}\right]}^2}$$

(44)

As for solutions that have the same $GR$ value, the one that has a smaller $GCPD$ value should be preferred. For instance, in Figure 16, individuals F and E have the same $GR$ value, and the $GCPD$ value of solution E is smaller than that of solution F, so solution E should be selected.

4.2. Equilibrium Optimizer

DMOEOA involves the effective properties of the single-objective equilibrium optimizer (EO), which is inspired by the control volume mass balance model [35]. In EO, each individual in the population with its concentration $C_o$ is a search agent, the individual is similar to a solution in PSO, and the $C_o$ of each individual is similar to the position of a particle. The position updating rule of EO is shown below:

$$C_o=C_{oe}+\frac{G_e}{\lambda V}.(1-F_e)+(C_o-C_{oe}).F_e$$

(45)

where $C_{oe}$ indicates the equilibrium candidate, $G_e$ and $F_e$ represent the generation rate and exponential term, respectively, and V is set to the unit. $\lambda =[{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n{]}^T$ refers to a random vector in the interval of [0, 1] and n represents dimensions of the individual’s $C_o$.

4.2.1. Equilibrium Pool $C_{oe,pool}$

The equilibrium state refers to the final convergence state of obtained solutions. In the initial stage of the search process, there is no knowledge about the convergence state, so the equilibrium candidate $C_{oe}$ is employed to provide the search guide for other individuals. In DMOEOA, an external archive $Arc$ in MOPSO is adopted as the equilibrium pool. The $Arc$ is used to store non-dominated solutions in the search process. Solutions in the $Arc$ are defined as equilibrium candidates. The equilibrium pool $C_{oe,pool}$ can be described as follows:

$$C_{oe,pool}=\left\{Arc\right\}$$

(46)

Individuals in the population update positions with a roulette wheel selection from equilibrium candidates. The more equilibrium candidates $C_{oe}$ with the same GR value in the $C_{oe,pool}$, the less likely they are to be chosen to guide the individuals in the population. This selection method can maintain diversity of the obtained non-dominated solutions in the optimization process.

4.2.2. Exponential Term $F_e$

Exponential term $F_e$ that controls the position updating rule is calculated as follows:

$$F_e=a_{1}sign(h_0-0.5)(e^{-\lambda t}-1)$$

(47)

$$t={(1-\frac{cit}{MT})}^{\left(a_2\frac{cit}{MT}\right)}$$

(48)

where $a_1$ and $a_2$ are constant values that control the exploration and exploitation abilities of DMOEOA, respectively, $sign(h_0-0.5)$ is employed to control the direction of exploitation and exploration and $h_0$ is a random number in the interval of [0,1]. t refers to the function of iterations, which declines with the number of iterations, and $MT$ and $cit$ represent the maximum iteration and the current iteration, respectively. The values of $a_1$ and $a_2$ are set to 2 and 1, respectively; in this work, the choice of these two values is consistent with original DMOEOA.

4.2.3. Generation Rate $G_e$

Generation rate is applied to enhance the exploitation ability of DMOEOA, and it can be calculated as follows:

$$G_e=GCP(C_{oe}-\lambda C_o)$$

(49)

$$GCP=\left\{\begin{array}{c}0.5h_1{h}_2\ge GP\hfill \\ 0h_2\le GP\hfill \end{array}\right.$$

(50)

where $GCP$ refers to generation rate control probability, $GP$ indicates the generation probability, which is set to 0.5 in this work. $h_1$ and $h_2$ represent two random numbers in the interval of [0,1].

4.3. Layered Disruption Method

In DMOEOA, a novel mutation method named layered disruption method (LDM) is proposed. LDM is inspired by the disruption phenomenon originated from astrophysics.

“When a swarm of gravitationally bound particles having a total mass, $m_k$, approaches too close to a massive object, $M_k$, the swarm tends to be torn apart. The same thing can happen to a solid body held together by gravitational forces when it approaches a much more massive object” [36]; this is the statement of disruption phenomena.

In LDM, individuals that have the same $GR$ value are regarded as one group, and there are different disruption conditions in different groups. The disruption coefficient $Q_i$ is calculated as follows:

$$Q_i=exp[-{\left(\frac{i}{\gamma }\right)}^2]$$

(51)

where $\gamma$ represents the number of groups, and i is the ith index after sorting different groups by increasing order based on the $GR$ values. In the ith group, individuals with $S{N}_i$ smallest $GCPD$ values are regarded as the mass $M_k$, and individuals that will be disrupted have the total mass $m_k$. $S{N}_i$ is calculated as follows:

$$S{N}_i=\left[Q_i{U}_i\right]$$

(52)

where $U_i$ represents the number of individuals in the ith group, $[·]$ indicates the function of rounding up. When the individual satisfies the disruption condition, the random number that obeys the Cauchy distribution is used to disrupt the individual. The $cauchyrnd$ is determined as follows:

$$f\left(x\right)=\frac{1}{\pi (1+x^2)},x\in [-\infty ,+\infty ]$$

(53)

The disruption equation is described as follows:

$$C_{oj}\left(cit\right)=\left(\frac{cit}{MT}\right)C_{oj}\left(cit\right)+Cau.C_{oj}\left(cit\right)(1-\frac{cit}{MT})$$

(54)

where $C_{oj}$ represents the position of individual j, and $Cau$ is a matrix that is composed of a set of Cauchy random numbers.

4.4. Constraint Handling

As there are constraints in the MMDO of the robotic fish, the constrained-dominate principle [16] is adopted as the constraint handling technique in DMOEOA.

Definition 5.

(Constrained-domination): Given two solutions i, j. i constrained-dominates j (denoted as $i\prec j$) $iff$:

1. 

Solution i is feasible and solution j is non-feasible.

2 

Both solution i and solution j are non-feasible, and the constraint violation of solution i is smaller than that of solution j.

3. 

Both solution i and solution j are feasible, but solution i dominates solution j.

4.5. Pseudo Code of DMOEOA

The pseudo code of DMOEOA is described in Figure 17. The computational complexity of DMOEOA is $o\left(k{N}^2\right)$, k indicates the number of objectives and N represents the number of individuals in the population. The computational complexity of DMOEOA is the same as some famous multi-objective optimization algorithms, including NSGA-II, MOPSO, MOWOA, MOGWO. More information about DMOEOA may refer to [22].

4.6. Parameter Setting

The parameters used in DMOEOA are summarized in

Table 9. Parameters of DMOEOA.

Parameter	Value
$Gd$	10
$a1$	2
$a2$	1
$GP$	0.5
$MT$	300
$PopsizeN$	30
$ArcsizeAr$	50

. The selection of these parameters is consistent with the original DMOEOA. The optimization process is performed on a desktop computer with an Intel Core i5-8400 CPU clocking at 2.80 GHz, and the RAM of the computer is 16 GB. The average running time of each generation in the optimization process is 15.52 s.

5. Optimization Results and Discussion

The optimization results of the robotic fish are shown in Figure 18, Figure 19 and Figure 20 and

Table 10. Comparison of inputs and outputs between the initial design and four selected cases.

	Initial Design	Case1	Case2	Case3	Case4
$l_1$	12	10	10.007	10.77	10.055
$l_2$	11	10	10.005	10	10.021
$l_3$	4	4	4.01	4	4.04
$w_1$	5.589	6.283	5.679	4.955	4.773
$w_2$	5.1	4.28	4.003	3.142	3.142
$w_3$	3.142	3.142	3.494	3.142	3.178
$A_1$	0.533	0.524	0.849	0.743	0.695
$A_2$	0.524	0.524	0.524	0.524	0.534
$A_3$	1.047	1.032	1.046	1.047	1.042
$\phi$	0.524	0.751	0.769	0.579	0.64
$q_{f1}$	2	2	2	2.078	2.313
$q_{f2}$	1	1.901	2.593	0.0213	0
$q_{a1}$	3	3.446	2.929	3.273	3.574
$q_{a2}$	12	16.326	14.405	14.105	14.71
$F_{drag}$ (N)	12.858	11.941	12.499	12.462	12.577
v (m/s)	0.104	0.2	0.4117	0.2827	0.2148
$H_{en}$ (h)	10.792	11.779	6.345	11.443	15.516

. The obtained Pareto solutions of different intermediate generations (50, 100, 200 and 300) in Figure 18 indicate the convergence process, and it can be observed that more and more non-dominated solutions are found with the number of iterations, and the obtained Pareto solutions are uniformly distributed in the objective space in the 300th generation. As shown in Figure 19, as the forward velocity of the robotic fish increases, the endurance of the robotic fish decreases. To balance the performances of different design objectives of the robotic fish, four cases in the middle of the Pareto front shown in Figure 19 are selected for analysis.

The comparison between optimization results of those four selected cases and the corresponding initial design is shown in Table 10. The values of optimization objectives are shown in bold. The obtained results listed in Table 10 suggest that all the four selected cases have less fluid resistance than the initial design, and the comparison between the optimal lines of those four cases and the initial line is shown in Figure 20. In addition, the four selected cases show better performance in the forward velocity than that of the initial design.

Among the four selected cases, case1 has the least fluid resistance, case2 has the fast forward velocity and case4 shows better performance in endurance than that of the other three cases. Compared to case2, case1 shows longer endurance and smaller fluid resistance, and the forward velocity of case2 is more than twice as fast as that of case1. Case3 and case4 have similar fluid resistance values, and the forward velocity of case3 is faster than that of case4; however, case4 shows better performance in endurance than case3. Although case2 has a 41% reduction in endurance compared to the initial design, the forward velocity of case2 is four times as fast as that of the initial design. All the four selected cases show better comprehensive performance than the initial design, which validates the effectiveness of the IDF-DMOEOA strategy.

6. Summary and Conclusions

In this paper, a novel multi-objective multidisciplinary design optimization strategy (IDF-DMOEOA) combining the IDF approach and the DMOEOA algorithm is proposed for the conceptual design of a three-joint robotic fish system. The robotic fish system is divided into four disciplines: hydrodynamics discipline, propulsion discipline, weight and equilibrium discipline, and energy discipline. The fluid resistance of the hull, forward velocity and endurance of the three-joint robotic fish are adopted as the optimization objectives. The CFD method is applied to predict the hull’s fluid resistance, and the grid convergence study is conducted to prove the stability of the mesh for simulation. The backpropagation neural network is employed as the surrogate model in the hydrodynamics discipline analysis to reduce the computational expense. The Pareto front is obtained after 300 generations. Four cases are chosen to compare with the initial design. The optimization results suggest that the optimized robotic fish shows better comprehensive performance than the initial design, which validates the proposed IDF-DMOEOA strategy’s effectiveness. In this work, only four disciplines are considered. In future work, more aspects of the multi-joint robotic fish’s conceptual design, such as maneuverability and structure, should be taken into account, and the proposed IDF-DMOEOA strategy is expected to have a wide application in other MMDO problems of complex systems.