Modeling and Trajectory Tracking Model Predictive Control Novel Method of AUV Based on CFD Data

In this paper, a novel model predictive control (MPC) method based on the population normal probability division genetic algorithm and ant colony optimization (GA-ACO) method is proposed to optimally solve the problem of standard MPC with constraints that generally cannot yield global optimal solutions when using quadratic programming (QP). Combined with dynamic sliding mode control (SMC), this model is applied to the dynamic trajectory tracking control of autonomous underwater vehicles (AUVs). First, the computational fluid dynamics (CFD) simulation platform ANSYS Fluent is used to solve for the main hydrodynamic coefficients required to establish the AUV dynamic model. Then, the novel model predictive controller is used to obtain the desired velocity command of the AUV. To reduce the influence of external interference and realize accurate velocity tracking, dynamic SMC is used to obtain the control input command. In addition, stability analysis based on the Lyapunov method proves the asymptotic stability of the controller. Finally, the trajectory tracking performance of the AUV in an underwater, three-dimensional environment is verified by using the MATLAB/Simulink simulation platform. The results verify the effectiveness and robustness of the proposed control method.


Introduction
The autonomous underwater vehicle (AUV) is a type of underwater unmanned vehicle (UUV) that has become important to ocean exploration and as a development tool in which artificial intelligence and other advanced computing technologies are integrated [1,2]. Because AUVs operate without an umbilical cord, they are more flexible in underwater operations [3]. In the early stages of deployment, AUVs were mainly used for the development of offshore oil and gas facilities, etc., although later, this technology played an important role in the field of military salvage [4,5]. The AUV has increasingly become of interest to researchers in the marine, oil, industry and military fields of countries around the world [6,7]. So that the AUV can successfully complete its assigned underwater operations, the performance of its control system is critical, especially the design of its trajectory tracking control method [8,9]. Trajectory tracking is limited by time constraints [10]. This technology feature only needs to track the reference path within a certain error allowable range to guide the AUV from the initial state to the final state [11]. For the state error of the system to converge to zero, the selection and optimization of the control method is highly important [11]. At present, many control methods have been successfully applied in the field of AUV trajectory tracking control, including proportional-integral-derivative (PID) control, neural network (NN) control, sliding mode control (SMC) and model predictive control (MPC) [12,13].
PID control is the most widely used control method in engineering, as well as the most stable. Among the various methods, the proportional (P) aspect corrects the deviation so that the process can respond quickly [14,15]. The integral (I) represents information accumulated in the past, which eliminates static errors and improves the static characteristics of the system [16]. The derivative (D) plays a leading control role in signal changes, reduces overshoot, overcomes oscillations and improves system stability [17]. However, an AUV is a high-inertia system. When the expected control command changes suddenly, the output of the AUV cannot be adjusted in time, which results in large system errors [18]. Therefore, PID control is usually combined with adaptive control, fuzzy control or neural network control, among others [19].
NN control is a current frontier in the field of automatic control [20,21]. This method represents a new branch of intelligent control, offering a new way to solve the control problems of complex nonlinear, uncertain and uncertain systems [22,23]. The NN controller is actually a feedforward controller that establishes the inverse model of the controlled object and adjusts the weights between the NN by learning the input and output of the traditional controller as samples to allow the feedback control input to approach zero and replace the traditional controller [24]. The advantage of NN control is that a precise dynamic model of the control target is not needed; additionally, NN can be used to approximate the nonlinear system [25]. The disadvantage is that it is difficult to obtain learning samples. Because the AUV system is highly nonlinear and strongly coupled, these uncertainties challenge the ability of the offline neural network controller in regard to robustness; as a result, the development of an online NN controller should be emphasized [26].
SMC is a special kind of nonlinear control that can, in essence, overcome the uncertainty of the system and has strong robustness in the face of an unknown disturbance, thus it is suitable for the control of AUVs [27]. However, when the state trajectory of the SMC reaches the sliding mode surface, it is difficult to slide strictly along the sliding mode toward the equilibrium point and instead passes back and forth on both sides of the SMC to approach the equilibrium point, resulting in chattering [28]. Elmokadem et al. derived an SMC method for the trajectory tracking control of an AUV but did not solve the chattering problem exhibited by the controller [29]. Londhe and Patre proposed an adaptive fuzzy sliding mode controller by deriving fuzzy logic rules using the Lyapunov energy function to reduce the SMC chattering problem [30].
When using the MPC feedback control method, the current control action is obtained by solving a finite time-domain open-loop optimal control problem at each sampling moment [31,32]. The current state of the process is the initial state of the optimal control problem, and the optimal control sequence only implements the first control effect [33,34]. The advantage of MPC is that it is good at dealing with problems with constraints. In traditional MPC, the process of solving the optimal control sequence is usually transformed into a quadratic programming (QP) problem, which solves the minimum value of a positive definite quadratic function with linear equality constraints [35,36]. However, the QP process often cannot obtain the global optimal solution. To solve this problem, the improved particle swarm optimization (PSO) algorithm, genetic algorithm (GA) and ant colony optimization (ACO) algorithm are usually introduced into MPC to replace the QP solution process. Ref. [37] Gan et al. applied quantum particle swarm optimization (QPSO) to model a predictive controller to solve the optimal control quantity. The application of QPSO enhanced the global search ability of traditional PSO and applied the controller to the trajectory tracking control field of AUVs, achieving a good tracking effect. Ref. [38] Zhang et al. proposed an improved GA to optimize the Q and R matrices in the control objective function of the MPC method; unfortunately, the GA easily falls into a local optimal value due to the lack of population diversity. Additionally, it is difficult for the algorithm to converge quickly in the later stage, and it is difficult to obtain the optimal solution within the MPC sampling time. The ACO algorithm can compensate for the slow convergence rate of the GA in the later stage, although the initial convergence rate of the algorithm is slow due to its low initial pheromone concentration.
To solve the above problems, this paper proposes an MPC method optimized by the GA-ACO algorithm based on population normal probability division. The standard GA of the fitness value according to the normal probability interval is divided into the following steps: put forward the internal and external hybrid crossover operator and dynamic hybrid mutation operator to maintain the diversity of the population, and later, using the ACO algorithm to accelerate convergence, make MPC obtain the global optimal control sequence. Finally, combine the process with dynamic SMC, which is used in AUV trajectory tracking control. The controller recalculates the optimal input at the next time according to the systematic error at each sampling time. The hydrodynamic coefficients used in AUV modeling in this paper are computational fluid dynamics (CFD) data based on the actual solution of ANSYS Fluent.
The main content of this paper is as follows: Section 2 introduces the kinematic and dynamical model of the AUV. In Section 3, the solution methods and results for the main hydrodynamic coefficients of the AUV are introduced. Section 4 introduces the formulation of the optimization problem and the control method proposed in this article. Section 5 presents a stability analysis of the proposed control method based on terminal constraints. Section 6 gives the simulation results. Section 7 offers some concluding comments.

Modeling of the AUV
This paper focuses on AUV modeling and trajectory tracking control. As shown in Figure 1, the size and mechanical structure of the AUV proposed in this paper are designed using the SOLIDWORKS design platform, with the relevant parameter values shown in Table 1. The arrangement of thrusters supports the movement of the AUV with five degrees of freedom. To avoid the singularity of the Euler angle, the pitch angle θ of the AUV is bounded and satisfies −π/2 ≤ θ ≤ π/2. The AUV roll angle φ is also self-stable due to the recovery torque caused by buoyancy and other factors. Therefore, the default expected angle and velocity values of φ and θ in this paper are 0. In the forward direction, one main propeller is used to generate the thrust of the surge, two vertical propellers complete the thrust of the heave and the pitch moment, and two lateral propellers complete the thrust of the sway and the yaw moment.

Kinematic Modeling
To establish the kinematic model of the AUV to study its trajectory tracking control, the body coordinate frame B and the inertial coordinate frame I are introduced, both of which are defined in right-handed Cartesian coordinate systems. The origin of B is the center of gravity of the AUV, which is used to describe the velocity information of the AUV. The velocity vector of the AUV is defined as v = [u, v, w, q, r] T , where u, v and w are the velocities in the surge, sway and heave directions, respectively, and q and r are the angular velocities in the pitch and yaw directions. The origin of I may be any point on Earth, which is used to define the position and attitude information of the AUV in B with respect to I. The position and attitude vector of the AUV is defined as η = [x, y, z, θ, ψ] T , where the positive direction of the x, y and z axes points in the horizontal north direction, the horizontal east direction and the geocentric direction, respectively, and θ and ψ are the pitch and yaw angle attitudes of the AUV.
The kinematic model can be defined in the following form: Here, J(η) is the transition matrix between the two coordinate frames:

Dynamic Modeling
The dynamic model of the AUV, as proposed by Fossen, is expressed [35]: Here, M = M RB + M A is defined as the inertia matrix. M RB = diag(m, m, m, I y , I z ) is the rigid-body inertial mass matrix, which is a property of the AUV itself and does not change with time, where m is the mass of the AUV and I y and I z are the inertia tensors of the y B and z B axes. M A = diag(−Xu, −Yv, −Zẇ, −Mq, −Nṙ) is the additional mass matrix, which represents the hydrodynamic force due to the acceleration of the AUV underwater, where −Xu, −Yv, −Zẇ, −Mq and −Nṙ are the inertial hydrodynamic coefficients. C(v) = C RB + C A is the Coriolis and centripetal force matrix. C RB is the Coriolis and centripetal force matrix due to the rigid-body inertial mass, and C A is the hydrodynamic Coriolis and centripetal force matrix due to the additional mass. Specifically: is the fluid damping matrix, including the linear damping matrix D L (v) and the quadratic damping matrix D N (v), where D L (v) = diag(X u , Y v , Z w , M q , N r ) and D N (v) = diag(X |u|u |u|, Y |v|v |v|, Z |w|w |w|, M |q|q |q|), N |r|r |r|). Here, X u , Y v , Z w , N r , X |u|u , Y |v|v , Z |w|w , M |q|q |q|) and N |r|r are measurable viscous hydrodynamic coefficients. g(η) = [0, 0, 0, bZ g sin θ, 0] T is a vector representing the restoring forces and moments due to the forces of gravity and buoyancy, b. In this paper, the force of gravity of the AUV is set equal to its buoyancy. Z g is the vertical coordinate of the center of gravity. τ = [τ u , τ v , τ w , τ q , τ r ] T is the vector of the input forces and moments for each DOF.

Solution Method and Fluid Area Division
To realize the accurate modeling of the AUV as proposed in this paper, it is necessary to analyze its hydrodynamic performance and solve the main viscous hydrodynamic parameters and inertial hydrodynamic parameters. This part is based on the solution method of [39]. The solution process of this method is shown in Figure 2. Using the CFD solution ANSYS Fluent platform, the AUV simulates the corresponding translational motion and rotational motion in the seawater environment. Figure 3 shows the division of the solution domain for the translational motion, and Figure 4 shows the solution domain division of the rotational motion. The control method and parameter settings of the mesh divisions are shown in Table 2.

Simulation Process
As shown in Equation (3), through the analysis of each component in the dynamic model of the AUV, it can be concluded that the hydrodynamic force of the AUV when moving in the fluid includes the inertial hydrodynamic force or moment τ A caused by the additional mass and the viscous drag or moment τ D caused by the viscosity of the fluid. The expression is as follows: D L is the linear damping matrix, and D N is the nonlinear damping matrix. From Equation (6), it can be concluded that the viscous drag and moment of the AUV are related to velocity via a quadratic function. By simulating the translational and rotational motion of each degree of freedom of the AUV under different velocities, the values of drag and moment under different velocities can be solved, and the viscous hydrodynamic coefficients can be identified by curve fitting with the least squares method.
In all simulations, the fluid is seawater with a density of 1025 kg/m 3 and viscosity of 8.8871 × 10 −4 Pa · s. The turbulence model is K-epsilon.
The uniform translational and rotational motions of the AUV are simulated. Each simulation only gives the velocity value of one degree of freedom to the AUV, while the velocity value of any other degrees of freedom is zero.
The velocity cloud diagrams of the AUV at a uniform velocity of 1 m/s in the surge, sway and heave directions are shown in Figure 5. Table 3 shows the calculated viscous drag on the AUV at different velocities. The least square method was used to fit the velocity values of translational motion and the corresponding viscous resistance data to a quadratic curve through the origin, and the fitting results are shown in Figure 6.   The velocity cloud diagrams of the AUV moving around the z axis with a uniform rotational velocity of 0.5 rad/s are shown in Figure 7. The calculated damping moments of the AUV at different angular velocities are given in Table 4. The least square method was used to fit the rotational velocity values and the corresponding viscous moment data to a quadratic curve through the origin, and the fitting results are shown in Figure 8.     The approximate viscous hydrodynamic coefficients obtained by fitting the results are shown in Table 5.
According to Equation (7), it can be clearly seen that the inertial hydrodynamic force on the AUV has a linear relationship with its acceleration. The next simulation considers the uniform acceleration of the AUV to 2 m/s at different linear accelerations and to 1 rad/s at different angular accelerations, because the results obtained by solving for the state of the AUV at different accelerations contain the drag force and moment caused by the viscous hydrodynamic force. Therefore, the inertial hydrodynamic force due to the additional mass should be the difference between the results for the uniform-acceleration and uniformvelocity motion solutions. By fitting this difference to a linear relationship, the approximate inertial hydrodynamic coefficients can be obtained. The user-defined function (udf) is used to write programs to achieve uniform acceleration of the AUV translational motion and uniform acceleration of its rotational motion.  Figure 9 shows the velocity cloud diagrams of the AUV with a uniform acceleration of 0.1 m/s 2 in the surge, sway and heave directions. Table 6 shows the calculated results for the drag force on the AUV at different accelerations. The least square method was used to fit the acceleration values and drag data to a straight line through the origin, and the fitting results are shown in Figure 10.     The velocity cloud diagrams of the AUV in uniform rotational motion around the z axis with an angular acceleration of 0.01 rad/s 2 are shown in Figure 11. The calculated results for the moments of the AUV at different angular accelerations are given in Table 7. The least square method was used to fit the angular acceleration values and moment data to a straight line through the origin, and the fitting results are shown in Figure 12.    The approximate inertial hydrodynamic coefficients obtained by fitting the results are given in Table 8.

Design of the Controller
In this part, the AUV controller used to perform trajectory tracking tasks is designed in detail. The controller is composed of two cascaded controllers. The outer loop controller uses the novel GA-ACO algorithm based on the population division of the normal probability interval, combined with the constraints of velocity and position and the system output feedback, to solve the optimal control sequence of the model predictive controller, that is, the expected velocity command of the AUV, which is transmitted to the inner loop controller. The inner loop controller is a sliding mode controller that can obtain the control input instructions of the AUV and ensure tracking of the AUV along its entire trajectory. The control strategy block diagram is shown in Figure 13.

MPC Based on Novel GA-ACO
It is assumed that the trajectory is given in advance and is smooth and bounded. The expected trajectory is described as follows: The expected velocity command is obtained by inputting the error between the actual position of the AUV and the expected position. T s is the sampling time, and the motion model of the system is obtained after discretization as follows: To ensure stable motion of the AUV, the velocity increment is given as an input, and the state-space model of the AUV is expressed as follows: The velocity increment is defined as follows: Y(k) is defined as the output. By defining matrices Q A (k) and Q B (k) in the form given below, the equation of state can be rewritten as follows: with where I 5 is a 5 × 5 identity matrix and 0 5 is a 5 × 5 zero matrix. Equations (10) and (11) are the model prediction equations for the AUV, from which the future state of the system can be inferred. The prediction horizon is denoted by N p , and the control horizon is denoted by N c ; these horizons must satisfy N c ≤ N p . The control value is held constant outside the control horizon. The controller output sequences are as follows: The system input sequence is as follows: Based on the principle of MPC, the model prediction equations can be further expressed asX whereQ Constraints: The upper and lower bounds on the AUV state and the velocity increment are set as follows: By adjusting the current and future inputs to the system, the predicted performance cost function is defined as follows: where y 2 Q y = y T Q y y, ∆vs. 2 Q y and Q ∆v are symmetric positive-definite weight matrices. Within the sampling time T s , the cost function E(t) is discretized as follows: According to Equation (13), the above equation can be transformed into where C † is the pseudoinverse of the matrix C. According to Equation (17), the simplified equation can be expressed as follows: It can be seen from Equations (27)-(29) that E(k) is a function of ∆v(k); thus, the optimization problem can be transformed into the following problem of finding the minimum value of the cost function with constraints (20) and (21):

Novel GA-ACO for Solving the Cost Function
In this section, we introduce a novel GA-ACO method to replace the QP method in the standard MPC to solve the extreme value of the cost function and obtain the optimal control variable ∆v(k). The proposed population normal probability division method is helpful to maintain the diversity of the GA population, to avoid falling into the local optimal solution and to improve the global search ability of the algorithm. To improve the convergence speed in the later stage of the algorithm, the initial optimization result of the GA is used as the initial pheromone of the ACO algorithm, and the global optimal solution is obtained by the rapid convergence of the ACO algorithm.

1.
Chromosome coding and fitness function setting. Chromosomes are encoded in real numbers. As shown in Figure 14, the length of the chromosomes is consistent with the control horizon of MPC, where each locus corresponds to the predicted control amount in each sampling period of MPC. Taking Equation (27) as the fitness function of the algorithm, through continuous iteration, the optimal control sequence for each sampling period is obtained.
To reasonably improve the diversity of the population in the GA algorithm, a method of dividing the population normal probability interval is proposed. Equation (31) is a one-dimensional normal distribution function, with a probability density function given as Equation (32), where σ is the standard deviation and µ is the mean value. According to the 3σ criterion, the probability that the value is distributed in (µ − σ, µ + σ) is 68.26%, the probability of being distributed in (µ − 2σ, µ + 2σ) is 95.44%, and the probability of being distributed in (µ − 3σ, µ + 3σ) is 99.74%, which basically represents the entire probability space. The initial population is sorted according to the fitness value from greatest to least, and the population is divided into six families according to the above method. In this paper, the initial population number is set to 200, giving the division results shown in Table 9. 3. Crossover, mutation and convergence.
(1) Intra-family crossover: The optimal individual within the family is retained without a crossover operation. In the same family, individuals were randomly selected for intrafamily crossover, and the elite retention strategy was adopted.
The parents and offspring were sorted according to the fitness function value, and the optimal retention was selected. (2) Cross between families: Set the maximum and minimum number of iterations, when the number of iterations is less than the maximum number of iterations, and update the optimum population fitness values that are either not changed or for which the rate of change is small. For the crossover operation between families, randomly select the best individual in the families of two different crossover operations, and take the elite reserved strategy after the cross. The idea of interfamily crossover controls the hemming distance of crossover individuals, avoids inbreeding, helps to maintain population diversity and prevents the algorithm from falling into a local minimum.
Mixed mutation operator: Hybrid mutation integrates the advantages of Gaussian mutation, Cauchy mutation, Levy mutation and single point mutation to make all of the mutation types work together. In the initial state, the selection probability of the four mutation operators is the same, where all are 0.25. The elite retention strategy is implemented after the mutation operation. If the offspring fitness of the individual after mutation is better than that of the parent, then the probability of the individual type occurring with the corresponding mutation strategy is increased; if the offspring fitness value of the individual is lower than that of the parent after mutation, then the probability of the individual type occurring in the corresponding mutation strategy is reduced. Finally, the probability of updating the four mutation modes after each population mutation operation is determined. (4) Algorithm convergence: In the proposed algorithm, by combining the advantages of the GA and ACO algorithms, the ant colony algorithm adopts the maximum-minimum ant algorithm. The preliminary optimal solution obtained by the GA is used as the initial pheromone of the ACO algorithm, allowing the algorithm to quickly converge to the global optimal control sequence ∆v * (k), which makes up for the disadvantage of the slow convergence of the GA algorithm in the later stage. Take the first predicted value ∆v * (k + 1|k) of the optimal control sequence as the expected velocity control command.

4.
Population normal probability judgment. Due to crossover and mutation within the population, it is necessary to monitor the normality of the population fitness value during the iteration process. It is proposed to use the absolute value difference of the median, mode and mean of the population fitness value of each iteration as the threshold for monitoring normality. According to the observations of several simulation results in the later period, the threshold value is set as 0.05. When the data normality quality is low and the minimum value of the difference value is greater than 0.05, the normal probability interval family is terminated. If convergence to the optimal solution has been achieved, then the optimal path is output and executed. Otherwise, the random crossover and mutation algorithm is adopted, and the elite retention strategy is implemented. As shown in Figure 15, using the simulation tool Quantile-Quantile Plots, the initial population number is 200, the maximum number of iterations is 50, and the 5th to 50th and the normality monitoring results of the population fitness value at the 50th generation. Excluding the special values, the results show that the sample points all fall near the red standard normal distribution line, which conforms to the normal probability distribution.

Dynamic Sliding Mode Velocity Controller
In this part, we design a dynamic sliding mode velocity controller. The conventional SMC control law is usually designed as a discontinuous function, which leads to chattering in the whole system. To overcome the chattering problem of SMC, we put the discontinuous function into the first derivative of the control input, so that the derivative of the control input is discontinuous, but the control input becomes a continuous function after integration, which can effectively suppress or eliminate the chattering problem of SMC.
The dynamic sliding mode controller is designed to overcome the external interference of THE AUV in the offshore environment and improve the robustness of the whole AUV control system. The force and torque required to drive the AUV are designed according to the input velocity error.
The difference between the expected velocity of the system according to Equation (33) and the actual velocity of the AUV is used as the input to the velocity controller. This error is defined as follows: where According to the velocity error e v , the dynamic sliding surface is designed as follows: where δ is the sliding surface coefficient. The following form can be obtained by deriving the sliding surface: We take the sliding mode surface s as the system state, design the sliding mode surface again as follows: The derivative of sliding mode surface in Equation (36) is obtained and is equal to the exponential asymptotic law: Here, sgn is a discontinuous symbolic function, denotes the velocity at which the moving point of the system approaches the switching surface σ = 0.
It can be concluded from Equations (36) and (37) that the discontinuous function causing SMC chattering only appears in sliding mode surface σ, that is, the discontinuous function only appears in the derivative of control input. After integrating the sliding mode surface σ, the control input is a continuous function, which eliminates the chattering problem of dynamic SMC. According to Equation (3), we obtaiṅ According to Equations (3)-(5), by combining Equations (34) and (35), the amount of control for each degree of freedom of the AUV can be obtained: At each sampling time k, ∆v * (k) is recalculated. Then, the optimal control forces and moments are repeatedly computed and executed by the AUV to achieve rolling optimization. The optimization process iterates until the AUV completes the trajectory tracking task.

Stability Analysis
The stability of the controller proposed in this paper is analyzed in this section. Lemma 1. The stability of the controller is proven using the traditional method of proving Lyapunov stability in control theory, i.e., by finding a Lyapunov function of the system that is positive definite and whose its derivative is negative definite. Theorem 1. Suppose that the following statements are true: (1) When k = 0, the constrained optimization problem (30) has a solution ∆v(k|k). (2) The zero state of the system output is measurable. For the nominal system (without considering external disturbances or system uncertainty), the following can be concluded: (a) For any k > 0, the constrained optimization problem (30) updated with the state measurement X(k) has a solution. (b) The closed-loop system described by Equations (12) and (13) are nominally stable.

Lemma 2.
If the continuous differentiable function f(t), when t approaches infinity, has a finite limit value, and the derivative of f(t) is uniformly continuous, then lim x→+∞ (ḟ (t)) = 0 Proof of Lemmas 1 and 2 and Theorem 1. Suppose that the solution to the optimization problem at time k is expressed as follows: . . .
This approach satisfies the control constraints (20) and (21). The terminal constraint of the model predictive controller is set as follows: The predicted optimal state sequence and output sequence are as follows: These sequences satisfy the output constraints (20) and (21) and the terminal constraint (45). The optimal value of the cost function is: For the model predictive controller based on MPC-GA-ACO, the optimal cost function E * (k) is selected as the Lyapunov function L * m (k): Clearly, Equation (48) satisfies L * m (0) = 0 with k = 0 and L * m (k) > 0 with arbitrary k = 0.
For the dynamic sliding mode velocity controller, Since the expected trajectories of AUV designed in this paper are all functions of time t, the system is a nonautonomous system, and its stability can be proved according to Lemma 2. the Lyapunov function L * s (k) is defined as follows: From Equations (55) and (56), It can be readily shown that L * s (k) ≥ 0 anḋ L * s (k) < 0. Since σ andσ are bounded,L * s (k) is bounded andL * s (k) is uniformly continuous. From Lemma 2, we can get that when t approaches infinity,L * s (k), σandS approaches zero. It can be obtained from Equations (34), (35) and (36) thatė v = −e v . So we can get that as follows: It can be obtained from Equations (57) that when time t approaches infinity, e v approaches zero, that is, the output of the system can track the expected signal, so the system is asymptotically stable. Lyapunov function L * s (k) still satisfies theorem 1, that is, L * s (k) is positive definite and its derivative is negative definite. It can also prove the asymptotic stability of the dynamic sliding mode controller designed in this paper. The asymptotic stability of the inner loop controller and the outer loop controller ensures the overall stability of the AUV control system.

Simulation Results
In this section, the effectiveness and robustness of the controller proposed in this paper (MPC-GA-ACO) are verified using the MATLAB/Simulink simulation platform, and the simulation results for the proposed controller are compared with those for a GA model predictive controller (MPC-GA) [35] and a standard model predictive controller (MPC) [38]. The above three controllers all form a cascade control with dynamic SMC. The main hydrodynamic parameters of the AUV used in the simulations are shown in Tables 5 and 8.
Random disturbances are also introduced in the simulation as follows: where rand(1) is a normally distributed noise signal with mean 0 and variance 1.
The initial position and pose vector of the AUV is η(0) = [0, 5, 0, 0, 0] T . Expected trajectory 1 of the AUV in the simulation is as follows: The expected trajectory 2 of the AUV in the simulation is as follows: (60) Figure 16 shows the trajectory 1 tracking performance of the AUV without disturbance. The simulation results show that the three controllers can give the AUV good trajectory tracking performance. Obviously, the controller designed in this paper is the best of the three controllers in the initial speed of tracking the expected trajectory and the final tracking performance (red curve). Figure 17 shows a comparison of the tracking error results of each degree of freedom of the AUV and the tracking errors of the three control methods are bounded. Although the tracking error of the controller designed in this paper is lower than that of the standard model predictive controller in the initial stage, its overall tracking error result is better than that of the other two controllers (red curve). Table 10 shows the average position error of the three control methods in tracking the desired trajectory. It can be clearly seen from the table that the average error of the controller designed in this paper in five degrees of freedom is the smallest. The average position error data of surge, sway, heave, pitch and yaw are, respectively, 0.1431 m, 0.1654 m, 0.0999 m, 0.1178 m, 0.0987 m.
As shown in Figure 18, the control input stability of the controller proposed in this paper is the best in the comparison of the control input of the three controllers without disturbance (red curve).   Figure 19 shows the simulation results of trajectory 1 tracking performance with the random disturbance. Compared with the simulation results without disturbance, the random disturbance interferes with the motion of AUV in the three degrees of freedom of surge, sway and heave. Obviously, the tracking performance of the three control methods is slightly reduced. However, compared with the other two controllers, the controller proposed in this paper still has the best tracking performance (red curve). Figure 20 shows the comparison of the tracking error results of each degree of freedom of AUV, and the tracking error of the three control methods is bounded. The error fluctuates with the random disturbance, and the error fluctuation of standard MPC controller is the largest (green curve). However, the controller proposed in this paper has the minimum tracking error and the error fluctuation is stable (red curve). Table 11 shows the average position error of the three control methods when tracking the expected trajectory. It can be clearly seen from the table that the average position error of the controller designed in this paper is the smallest under five degrees of freedom. The average position error data of surge, sway, heave, pitch and yaw are 0.1893 m, 0.1867 m, 0.1445 m, 0.1498 m and 0.1612 m, respectively.  Figure 21 shows the control input fluctuation of the three controllers in the case of random disturbance. The control input stability of the proposed controller is still the best among the controllers and the fluctuation is the smallest (red curve). To further verify the performance of the controller, the trajectory tracking the performance of AUV is simulated again in aperiodic function trajectory 2. Figure 22 shows the trajectory 2 tracking performance of AUV without disturbance. The overall tracking effect of the standard MPC controller is poor (green curve). Although the overall tracking effect of the MPC-GA controller is improved compared with the standard MPC controller, the tracking speed is the slowest in the initial stage of trajectory tracking (blue curve). Because the new GA-ACO algorithm can reasonably maintain the diversity of the population, and ensure the local search and global search ability of the algorithm at the same time, so that the optimal control sequence can be obtained in each sampling period of MPC, the simulation results show that the controller proposed in this paper not only ensures the tracking speed of AUV in the initial stage of trajectory tracking but also the overall tracking effect is the best of the three (red curve).   Figure 23 shows the comparison of the tracking error results of each degree of freedom of AUV, and the tracking errors of the three control methods are bounded. Among them, the standard MPC controller error fluctuates greatly when the AUV bow direction change rate is large (green curve), but the controller proposed in this paper has the smallest tracking error (red curve). Table 12 shows the average position error of the three control methods when tracking the desired trajectory. It can be seen from the table that the average error of the controller designed in this paper is the smallest under five degrees of freedom.
The average position error data of surge, sway, heave, pitch and yaw are 0.1881 m, 0.1825 m, 0.1110 m, 0.0915 m, and 0.0873 m, respectively. Figure 24 shows the control input fluctuation of the three controllers without disturbance. The control input stability of the proposed controller is still the best and the fluctuation is the smallest (red curve).     Figures 25-27 show the trajectory 2 tracking performance of the AUV with random disturbance. Compared with the other two controllers, the controller achieves the optimal trajectory tracking performance under the condition of minimum error and smooth control intput under random disturbance. The random disturbance leads to a high-frequency variation of the desired speed command, but the proposed controller can still track the desired speed perfectly, which proves its effectiveness and robustness (red curve). Table 13 shows the average position error of the three control methods in tracking the desired trajectory. It can be clearly seen from the table that the average error of the controller designed in this paper in five degrees of freedom is the smallest. The average position error data of surge, sway, heave, pitch and yaw are, respectively, 0.1932 m, 0.1894 m, 0.1777 m, 0.0998 m, 0.1221 m.

Conclusions
In this paper, a novel AUV double closed-loop trajectory tracking method based on model prediction and sliding mode cascade control is designed. First, a kinematics analysis of the fully driven AUV is carried out, and the kinematics and dynamics equations of five degrees of freedom are established. Second, the CFD software ANSYS Fluent is used to simulate the uniform velocity and uniform acceleration of AUV in surge, sway, heave, pitch and yaw degrees of freedom, and solve the main hydrodynamic coefficients required to establish the dynamic model. Then, a novel GA-ACO algorithm based on normal probability interval population division is proposed is proposed. Combined with the constraints of actual control input and state, the algorithm is applied to the outer loop MPC to solve the optimal velocity control sequence, generate the expected velocity command, which is passed to the inner loop. Compared with standard MPC and MPC-GA, this method can reasonably increase the diversity of GA population, improve the local and global search ability of the algorithm, and avoid the problem of falling into local optimal solution caused by QP or GA. Then combined with ACO algorithm, the convergence speed of the algorithm can be accelerated and the global optimal speed control sequence can be obtained for each sampling period. The inner loop uses dynamic SMC to generate the available control input of AUV to ensure the closed-loop tracking along the whole trajectory. A rolling time-domain implementation can be used to recalculate the optimal control input at each sampling instant, compensating for system uncertainties caused by modeling and external disturbances. Finally, two different types of reference desired trajectories are simulated to verify the effectiveness and robustness of the proposed controller under uncertain disturbances in the simulated marine environment. It is improved in reducing the position error of trajectory tracking and increasing the input stability of the controller. In the future research, we will further solve the adverse impact of SMC chattering on the controller, and strengthen the analysis of the impact of the actual marine environment on the control disturbance of AUV and the research on practical problems such as thruster output saturation.