Multi-Objective Multidisciplinary Design Optimization of a Robotic Fish System

: Biomimetic robotic ﬁsh systems have attracted huge attention due to the advantages of ﬂexibility and adaptability. They are typically complex systems that involve many disciplines. The design of robotic ﬁsh is a multi-objective multidisciplinary design optimization problem. How-ever, the research on the design optimization of robotic ﬁsh is rare. In this paper, by combining an efﬁcient 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 ﬁsh 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 ﬁrst applied to the design optimization of the robotic ﬁsh system, and the robotic ﬁsh system is decomposed into four disciplines: hydrodynamics, propulsion, weight and equilibrium, and energy. The computational ﬂuid dynamics (CFD) method is employed to predict the robotic ﬁsh’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 ﬁsh shows better performance than the initial design, proving the proposed IDF-DMOEOA strategy’s effectiveness.


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

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.

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.

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.

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.
where qa 1 and qa 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.
Initial values of parameters that will influence the shape of the robotic fish are listed in Table 1.

Shape Parameter Value
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 qa 1 , qa 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.

Input Variable
Output Variable 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.

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

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]: 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 × 10 6 , and the density of the water ρ is set as 998.2 kg/m 3 . 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 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.

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.

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

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

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

4.
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. From Figure 8 and Table 4 we can observe that the BPNN model has sufficient precision for hydrodynamics discipline analysis.

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 f h, the flexible fish body f b, the after-body a f 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 sm 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 θ t also needs to be considered. θ t can be calculated by the following equations: where r cm = [ x cm y cm z cm ] T and r cb = [ x cb y cb z cb ] T represent the barycenter and buoyant center of the robotic fish, respectively, m f h , m f b , m a f , m sw represent the masses of the four components of the robotic fish, respectively, and V f h , V f b , V a f , V sw are the volumes of the four components, respectively. r f hm , r f bm , r a f m , r swm indicate the barycenters of those four components, respectively, and r f hb , r f bb , r a f b , 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. Table 5. Summary of weight and equilibrium discipline. Figure 9. General overview of the robotic fish.

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.

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, . . . , 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, . . . , 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 θ i (i = 1, 2, . . . , nl − 1); in particular, θ 0 indicates the angle between the 0th link and the X axis of the world coordinate system and R 0 w can be expressed as follows: Figure 10. Coordinate systems and notations of the robotic fish.

Kinematic Analysis
Let the coordinate of the joint O 0 in WCS be (x, y, 0) T . The position vector r i oi (i = 0, 1, 2, . . . , 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: Therefore, the position vector r i ci (i = 0, 1, 2, . . . , 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: The angular velocity w i and the angular acceleration α i of the ith link relative to the origin of WCS can be expressed in O i − X i Y i Z i as: 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:

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 i ri and the relative acceleration a i ri are as follows: where v i f i and a i f 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: 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: The fluid drag F di of the ith link can be calculated by the following equation [31]: where J ci indicates the rotational inertia of the ith link. The motion rule of the ith joint is as follows: where A i and ω i represent the angle of amplitude and angular velocity of the ith joint, respectively. ϕ 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: 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. Table 6. Summary of the proplusion discipline.

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]: where m b and ρ 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 sm 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: where ω 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 g1 max 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: where ρ 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 PP 1 can be calculated as follows: 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: The main objective of energy discipline is to estimate the endurance of the robotic fish H en . The calculated torque of sm 1 should not exceed the maximum torque T dmax . The input variables and output variables of the energy discipline are summarized in Table 7.

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: 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 , ω 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.

Hydrodynmaics discipline
Weight & equilibrium discipline Propulsion discipline Energy discipline

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.

Individual Discipline Feasible Approach
The optimization model of the IDF approach is described as follows [32]: where X = [X 0 , X 1 , . . . , X nd ] T refers to the vector of design variables, nd is the number of disciplines. X 0 and X i (i = 1, 2, . . . , 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 t i is the copy of y i . y i represents the state variable of discipline i. c i and c c i 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. Objective functions Figure 14. The XDSM of the IDF architecture.

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. Using IDF architecture, the optimization process can be performed in parallel.

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.

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 (x) and min f i (x) 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: 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:
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: 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.
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.

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: 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. λ = [λ 1 , λ 2 , . . . , λ n ] T refers to a random vector in the interval of [0, 1] and n represents dimensions of the individual's C o .

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

Exponential Term F e
Exponential term F e that controls the position updating rule is calculated as follows:

Generation Rate G e
Generation rate is applied to enhance the exploitation ability of DMOEOA, and it can be calculated as follows: 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].

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: where γ 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 SN i smallest GCPD values are regarded as the mass M k , and individuals that will be disrupted have the total mass m k . SN i is calculated as follows: 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: The disruption equation is described as follows: where C oj represents the position of individual j, and Cau is a matrix that is composed of a set of Cauchy random numbers.

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

Pseudo Code of DMOEOA
The pseudo code of DMOEOA is described in Figure 17. The computational complexity of DMOEOA is o(kN 2 ), 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].

Parameter Setting
The parameters used in DMOEOA are summarized in Table 9. 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.

Optimization Results and Discussion
The optimization results of the robotic fish are shown in Figures 18-20 and  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.

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.