Vibration Suppression Trajectory Planning of Underwater Flexible Manipulators Based on Incremental Kriging-Assisted Optimization Algorithm

Abstract

It is of great significance to expand the functions of submarines by carrying underwater manipulators with a large working space. To suppress the flexible vibration of underwater manipulators, an improved sparrow search algorithm (ISSA) combining an elite strategy and a sine algorithm is proposed for the trajectory planning of underwater flexible manipulators. In this method, the vibration evaluation function is established based on the precise dynamic model of the underwater flexible manipulator and considering complex motion and vibration constraints. Simulation results show that the ISSA algorithm requires only 1/3.68 of the time of PSO. Compared to PSO, SSA and the opposition-based learning sparrow search algorithm (OBLSSA), the optimization performance is improved by 17.3%, 13.1% and 9.7%, respectively. However, because the complex dynamics model of the underwater flexible manipulator leads to large computational effort and a long optimization time, ISSA is difficult to apply directly in practice. To obtain a large number of optimization results in a shorter time, an incremental Kriging-assisted ISSA (IKA-ISSA) is proposed in this paper. Simulation results show that IKA-ISSA has good nonlinear approximation ability and the optimization time is only 3% of that of the ISSA.

1. Introduction

Nowadays, the extension of the functions of large submersibles using underwater manipulators for the dynamic deployment and recovery of unmanned underwater vehicles (UUV) is a major research direction [1,2,3]. UUVs have the characteristics of strong concealment, high intelligence, lower construction costs and no casualties. The integration of submarines and UUVs can combine each one’s advantages and improve the overall performance of submarines. As new collaborative models between UUVs and submarines are gradually enriched and perfected, the question of how to carry, deploy, dock and recover one or more UUVs on a submarine is a key technical issue in the new collaborative model. To realize the above operations, using a submarine to carry the manipulator with large-scale operation capability is an important solution. The recovery and docking process of UUVs with the underwater manipulator is shown in Figure 1, in which the underwater manipulator can be stored in a small space inside the submarine after the operation is completed. Considering the occupation of the limited collection space of the submarine and the impact on its load capacity, the manipulator should have a large working space, a small storage space and a small weight so that it can be stored in the small space of the submarine. To meet the above technical requirements for large-scale operation and small-space storage of the underwater manipulator, the manipulator is often designed to be very slender. Due to the low stiffness of slender manipulators, the end effector is prone to vibration, and previous control methods that ignore the flexible deformation of manipulators are not able to meet the high-precision requirements of underwater operations [4,5]. In addition, when the manipulator is working, the speed and acceleration at the end of the manipulator should be strictly limited to meet the corresponding constraints. It is also necessary to overcome the interference of hydrodynamic forces. These factors bring great challenges to the control of underwater manipulators. Although extensive research has been carried out on the control of flexible manipulators and underwater manipulators, the vibration suppression of underwater flexible manipulators is still a difficult problem to solve.

At present, the vibration suppression of flexible manipulators is mainly studied from the perspective of control methods and trajectory planning. Sahu and Patra proposed an observer-based backstepping control method for a two-degree-of-freedom flexible manipulator [6]. Yang et al. proposed a hybrid control scheme composed of a trajectory planning approach and adaptive variable structure control to suppress the elastic vibration [7]. Rahmani and Belkheiri designed a neural network adaptive controller of flexible multi-link robots. The adaptive controller has good robustness under uncertain disturbances [8]. Yavuz proposed an improved vibration control method that suppresses residual vibration by shaping the speed input [9]. Guo et al. proposed a residual vibration suppression method for Delta robots based on input shaping technology [10]. Xu et al. proposed a compensation method; this method is based on the analytical solution of a second-order vibration system and designs compensation torque to suppress vibration [11]. Sands designed a whiplash compensator for flexible space robots, which had zero residual vibration at the end of the maneuver [12]. Trajectory planning for vibration suppression was firstly proposed by Park and Park [13]. Their method uses the excitation force, elastic deformation or elastic potential energy as the objective function to find the joint trajectory with the smallest flexible deformation or elastic potential energy, so as to achieve the vibration suppression of the manipulator. The vibration suppression trajectory optimization mainly defines the motion trajectory of the robot joint as polynomial functions or B-spline functions, transforms the discrete dynamic variables into parameter optimization problems and suppresses the flexible vibration of the flexible robot by optimizing the parameters [14,15,16].

The most commonly used optimization algorithms include traditional optimization algorithms (simplex method and compound shape method) and intelligent optimization algorithms (particle swarm optimization algorithm, genetic algorithm and sparrow search algorithm). The compound shape method is a direct method to solve the constrained nonlinear optimization problem. Lin et al. proposed an optimization method for the shortest total travel time under physical constraints such as joint velocity and acceleration based on the compound shape method [17]. Yin et al. optimized the two joint trajectories based on the Nelder–Mead simplex method in the feedforward control of the flexible manipulator to reduce the residual vibration of the flexible manipulator [18]. Since the selection and replacement of each vertex of the complex shape algorithm needs to meet the requirement of decreasing the value of the objective function and constraints, the search efficiency is relatively low. In processing optimization problems with complex models and constraints, the intelligent optimization algorithms have more advantages than traditional optimization algorithms. Wu et al. used a uniform aperiodic fourth-order B-spline curve to describe the trajectory of the manipulator, and conducted vibration suppression optimization for B-spline control points based on PSO [14]. Li et al. used a polynomial function to describe the trajectory of the manipulator, and used PSO for the optimization of polynomial interpolation points to suppress vibration [19]. Yue et al. and Kazem et al. studied the trajectory planning of robots based on the genetic algorithm (GA) to optimize the time in motion [20,21]. Li et al. used the GA to control the coefficients of cubic spline interpolation to suppress the vibration of the manipulator [22]. Different from GA, PSO does not use hybridization, mutation or replication for individuals, but instead treats each individual as a particle without volume in a multi-dimensional search space. In the process of vibration suppression trajectory optimization, the PSO algorithm has lower computational complexity than the GA. PSO has no special requirements for the objective function, while the GA requires a constant positive fitness function. However, the standard PSO algorithm is an unconstrained optimization algorithm, while the vibration suppression optimization of underwater flexible manipulators is constrained by the joint angle and flexible vibration. For the application of the PSO algorithm in constrained optimization problems, researchers have proposed a series of schemes to solve constrained optimization problems, among which the penalty function method and the constraint-based individual sorting method are the most widely used [23,24,25]. Cao et al. proposed a PSO algorithm with a compression factor as a penalty function to correct the joint trajectory and suppress the vibration of flexible joints [26]. Wang et al. studied the application of the PSO strategy under constraints of trajectory planning for a free-floating dual-arm robot [27]. Due to the constraints, PSO is more prone to premature convergence in the constrained optimization process. Although various improved PSO algorithms to solve constrained optimization problems have improved the optimization performance to a certain extent, premature convergence still exists. SSA is a new heuristic algorithm for swarm intelligence, proposed by Xue and Shen [28]. Zhang et al. proposed a path planning method for bionic mobile robots based on the sparrow search algorithm [29]. Liu et al. proposed an improved sparrow search algorithm to solve the obstacle avoidance problem of UAV route planning, and achieved good results [30]. Compared with the traditional heuristic search method, the SSA has more diversified position updating strategies, faster convergence speeds and more extensive application scenarios, indicating its great potential in dealing with robot vibration suppression trajectory planning.

Although there are some studies on the vibration suppression of flexible manipulators, in these studies, the suppression optimization effect remains to be further improved due to the relatively simple constraints and largely simplified models of flexible manipulators, and the vibration suppression of underwater manipulators is not involved. In addition, most of the current optimization algorithms suffer from premature convergence and poor accuracy in the constrained optimization process, and the optimization time is too long and costly for the case of complex models, making it difficult to be applied directly to practice.

To solve the problems of the premature convergence and low accuracy of traditional optimization algorithms in the process of constrained optimization, this paper proposes ISSA, which combines the sine algorithm and elite opposition-based learning strategy to optimize the trajectory of constrained underwater flexible manipulators. To obtain a large number of optimization results in a shorter time, this paper proposes IKA-ISSA, based on ISSA, to quickly generate vibration suppression trajectories of underwater flexible manipulators. The main contributions of this paper are summarized as follows. (i) In the past, the research object of vibration suppression trajectory planning was land-based flexible manipulators, while the research object of this paper is underwater flexible manipulators, considering the influence of hydrodynamics. (ii) The previous vibration suppression trajectory planning of flexible manipulators has greatly simplified the model. This paper is based on a more accurate dynamic model, and also considers the velocity, acceleration, maximum flexible displacement and other complex constraints. (iii) Based on the combination of the sine algorithm and OBL strategy, an improved sparrow algorithm is proposed, which is used for the first time to solve the vibration suppression trajectory planning problem for underwater flexible manipulators. (iv) IKA-ISSA is proposed for the problems of large computation demands and long optimization times for complex models. The algorithm has good robustness and nonlinear approximation capability, and can obtain accurate vibration suppression trajectories in a short time.

The rest of this paper is organized as follows. Section 2 presents the dynamic modeling approach for an underwater flexible manipulator. Section 3 describes the motion and vibration constraints, and establishes the manipulator vibration evaluation function. Section 4 designs the vibration suppression trajectory strategy of the underwater flexible manipulator based on PSO and ISSA and verifies the optimization effect of ISSA by simulation. Section 5 proposes IKA-ISSA and verifies its high optimization accuracy and efficiency by simulation. Finally, the conclusions can be found in Section 6.

2. Materials and Methods

The modeling of underwater flexible manipulator systems involves the complex integration of multiple disciplines, such as fluid mechanics and multi-body dynamics. The underwater flexible manipulator is usually disturbed by hydrodynamics, as well as some uncertain factors. These influences and disturbances are highly nonlinear and time-varying, which increases the difficulty of flexible manipulator modeling [31,32]. The authors of this paper fully consider the coupling effect of hydrodynamic force and flexible vibration, and we have established the dynamic equation of rigid–flexible coupling in our previous work [33,34].

The mechanism diagram of a two-link rigid–flexible coupling manipulator can be simplified as in Figure 2, in which arm 1 is abstracted as a rigid link and arm 2 is abstracted as a flexible link. The flexible link adopts the Euler–Bernoulli beam model, and the transverse shear deformation is ignored. ${\theta }_i$ and $l_i$, respectively, represent the rotation angle and length of the i-th link. $w$ represents the flexible displacement. According to the assumption of small elastic deformation, only the first two-order modes are considered when $w$ is expanded based on the assumed mode method (AMM).

By combining the Morrison formula and Lagrange equation, the dynamic equations of the underwater flexible manipulator can be expressed in the following compact form [33,34]:

$$\left(\begin{array}{cc}{\mathit{M}}_{\theta \theta }& {\mathit{M}}_{\theta q}\\ {\mathit{M}}_{q\theta }& {\mathit{M}}_{qq}\end{array}\right)\left(\begin{array}{c}\ddot{\mathit{\theta }}\\ \ddot{\mathit{q}}\end{array}\right)+\left(\begin{array}{c}{\mathit{C}}_r\\ {\mathit{K}}_q\mathit{q}+{\mathit{C}}_f\end{array}\right)=\left(\begin{array}{c}\mathit{\tau }\\ \mathit{0}\end{array}\right)+\left(\begin{array}{c}{\mathit{F}}_{\theta }\\ {\mathit{F}}_q\end{array}\right)$$

(1)

where $\mathit{\theta }={\left[\begin{array}{cc}{\theta }_1& {\theta }_2\end{array}\right]}^T$ is the joint angle, $\mathit{q}={\left[\begin{array}{cc}{q}_1& q_2\end{array}\right]}^T$ is the first two-order modes, ${\mathit{M}}_{**}$ is the mass matrices, ${\mathit{K}}_q$ is the stiffness matrix, ${\mathit{C}}_r$ and ${\mathit{C}}_f$ is the terms of centrifugal forces, Coriolis forces, and gravity. $\mathbf{\tau }$ is the joint torque, and ${\mathit{F}}_{\theta }$ and ${\mathit{F}}_q$ are the generalized forces related to hydrodynamic force.

3. Vibration Suppression Trajectory Planning

3.1. Cubic Polynomial Trajectory Planning

According to (1), the flexible motion part of the dynamic equation of the underwater flexible manipulator can be expressed as

$${\mathit{M}}_{qq}\ddot{\mathit{q}}+{\mathit{K}}_q\mathit{q}+{\mathit{C}}_f-{\mathit{F}}_q=-{\mathit{M}}_{q\theta }\ddot{\mathit{\theta }}$$

(2)

From (2), it can be known that the inertia moment generated during the rotation of the flexible manipulator arouses the elastic vibration of the manipulator, which can be controlled by optimizing the trajectory of the manipulator during the working process to suppress the vibration of the manipulator. In order to improve the optimization efficiency and the search speed, the basic displacement value of the trajectory control point can be obtained by discretely using a reference trajectory curve. In this paper, a quintic polynomial is used as the reference trajectory curve; it is expressed as

$${\theta }_d(t)=\left[6{(\frac{t}{t_f})}^5-15{(\frac{t}{t_f})}^4+10{(\frac{t}{t_f})}^3\right]\times \left({\theta }_f-{\theta }_0\right)+{\theta }_0$$

(3)

where $t_f$ is the ending time; ${\theta }_0$ and ${\theta }_f$ are the positions of the start and end times of the manipulator, respectively.

After obtaining the reference curve, sufficient nodes are taken in the reference trajectory curve, and then some floating changes to the basic value are made. By adding the interpolation point increment to change the position of each interpolation point, a new set of interpolation point values are obtained. Then, different trajectory curves are obtained by curve fitting. As shown in Figure 3, the dotted line represents the initial trajectory, and the solid line is the optimization trajectory. Since the speed and acceleration of the underwater manipulator are limited by the rated power and working environment, the speed and acceleration constraints need to be satisfied in the optimization process, and the cubic polynomial function has second-order derivability at both the interpolation point and the interpolation interval; it has good stability and it is easy to calculate the velocity and acceleration. Therefore, after the trajectory control points are obtained, the cubic polynomial interpolation function is used to fit the obtained two adjacent control points.

The joint angle at the i-th trajectory control point can be expressed in the following form:

$${\theta }_{Bi}=\left\{\begin{array}{cc}\hfill {\theta }_0,& i=1\hfill \\ \hfill {\tilde{\theta }}_{Bi}+\Delta {\theta }_{i-1},& i=3,4...n-2\hfill \\ \hfill {\theta }_f,& i=n\hfill \end{array}\right.$$

(4)

where ${\tilde{\theta }}_{Bi}$ is the basic value of the trajectory control point and $\Delta {\theta }_{i-1}$ is its floating value.

Let $t_i$ be the corresponding time series at node i, and $Q_i(t)$ is a cubic polynomial over the time interval $[t_i,t_{i+1}]$. Since the second derivative of $Q_i(t)$ is a linear function, it can be expressed as

$$Q_i{}^{″}(t)=\frac{t_{i+1}-t}{h_i}{Q}_i{}^{″}(t_i)+\frac{t-t_i}{h_i}{Q}_i{}^{″}(t_{i+1}),i=1,2,\dots ,n-1$$

(5)

where $h_i=t_{i+1}-t_i$, and the interpolation function $Q_i(t)$ can be obtained by integrating $Q_i{}^{″}(t)$ twice and substituting the node conditions

$$\begin{array}{l}{Q}_i(t)=\frac{Q_i{}^{″}(t_i)}{6h_i}{\left(t_{i+1}-t\right)}^3+\frac{Q_i{}^{″}(t_{i+1})}{6h_i}{\left(t-t_i\right)}^3+\left[\frac{{\theta }_{B(i+1)}}{h_i}-\frac{h_i{Q}_i{}^{″}(t_{i+1})}{6}\right]\left(t-t_i\right)\\ +\left[\frac{{\theta }_{Bi}}{h_i}-\frac{h_i{Q}_i{}^{″}(t_i)}{6}\right]\left(t_{i+1}-t\right),i=1,2,\dots ,n-1\end{array}$$

(6)

In order to determine the undetermined coefficient of the cubic fitting function, each selected interpolation point should satisfy the joint trajectory function over the interpolation point, the joint trajectory function continuity, the joint trajectory first-order derivative function continuity and the joint trajectory second-order derivative function continuity. Coupled with the upper boundary conditions, the following equations in matrix form can be obtained [17]:

$$\mathit{D}\left[\begin{array}{c}{Q}_2{}^{″}(t_2)\\ Q_3{}^{″}(t_3)\\ ⋮\\ Q_{n-1}{}^{″}(t_{n-1})\end{array}\right]=\mathit{b}$$

(7)

In this paper, the trajectory planning adopts uniform interpolation and $h_i=h$. $\mathit{D}$ is the matrix describing the continuous relationship of the interpolation curve, which can be expressed as

$$\mathit{D}=\left[\begin{array}{ccccccc}6h& h& & & & & \\ 0& 4h& h& & & & \mathrm{O}\\ & h& 4h& h& & & \\ & & & ⋮& & & \\ & & & & h& 4h& 0\\ \mathrm{O}& & & & & h& 6h\end{array}\right]$$

(8)

The matrix $\mathit{b}$ can be expressed as

$$\mathit{b}=\left[\begin{array}{c}\frac{6}{h^2}(-{\theta }_0+{\theta }_{B3}-2h{\omega }_0-\frac{5h^2}{6}{a}_0)\\ \frac{6}{h^2}({\theta }_0+h{\omega }_0+\frac{h^2}{3}{a}_0+{\theta }_{B4}-2{\theta }_{B3})\\ \frac{6}{h^2}({\theta }_{B(i+1)}+{\theta }_{B(i-1)}-2{\theta }_{Bi})\\ ⋮\\ \frac{6}{h^2}({\theta }_f-h{\omega }_f+\frac{h^2}{3}{a}_f-2{\theta }_{B(n-2)}+{\theta }_{B(n-3)})\\ \frac{6}{h^2}(-{\theta }_f+h{\omega }_f-\frac{5h^2}{6}{a}_f+{\theta }_{B(n-2)})\end{array}\right]$$

(9)

where ${\omega }_0$ and ${\omega }_f$ are the angular velocities at the start time and the end time, and $a_0$ and $a_f$ are the angular acceleration at the start time and the end time.

In (8), the matrix is an upper triangular matrix, and each element on the diagonal is greater than zero; thus, this is a non-singular matrix. Therefore, the unique solution of $[{Q_2}^{″}(t_2), \cdots , {Q_i}^{″}(t_i), \cdots , {Q_{n-2}}^{″}(t_{n-2})]$ can be calculated. By substituting the solution into (21), the values of trajectory control point 2 and point n-2 can be obtained.

$${\theta }_{Bi}=\left\{\begin{array}{l}{\theta }_0+h\times {\omega }_0+\frac{h^2}{3}{a}_0+\frac{h^2}{6}{Q_2}^{″}(t_2), i=2\\ {\theta }_f-h\times {\omega }_f+\frac{h^2}{3}{a}_f+\frac{h^2}{6}{Q_{n-1}}^{″}(t_{n-1}), i=n-1\end{array}\right.$$

(10)

The final optimal trajectory is obtained by combining (4), (6) and (10).

3.2. Constraint Condition

When the underwater manipulator is working, the angular velocity and angular acceleration are limited by the rated power and working environment. The joint angular velocity and angular acceleration constraints of the i-th trajectory control point of joint j can be expressed as

$$\left|Q_{ji}^{\prime }(t)\right|\le \omega C_j, i=1,2,\dots ,n-1$$

(11)

$$\left|Q_{ji}^{″}(t)\right|\le a{C}_j, i=1,2,\dots ,n-1$$

(12)

where $\omega C_j$ is the velocity constraint for joint j, $a{C}_j$ is the acceleration constraint for joint j, and the joint angular velocity is a quadratic polynomial. The maximum of the absolute value of the angular velocity exists at the endpoint $t_i$, $t_{i+1}$ or extreme point $\overline{t_i}$. The velocity constraint can be expressed as

$$\max \left|Q_{ji}^{\prime }(t)\right|=\max \left\{\left|Q_{ji}^{\prime }(t_i)\right|\left|Q_{ji}^{\prime }(t_{i+1})\right|\left|Q_{ji}^{\prime }(\overline{t_i})\right|\right\}\le \omega C_j$$

(13)

where

$$\left|Q_{ji}^{\prime }(t_i)\right|=\left|\frac{a_{ji}}{2}h+\frac{({\theta }_{B(j,i+1)}-{\theta }_{B(j,i)})}{h}+\frac{(a_{ji}-a_{j,i+1})h}{6}\right|$$

(14)

$$\left|Q_{ji}^{\prime }(t_{i+1})\right|=\left|\frac{a_{j,i+1}}{2}h+\frac{({\theta }_{B(j,i+1)}-{\theta }_{B(j,i)})}{h}+\frac{(a_{ji}-a_{j,i+1})h}{6}\right|$$

(15)

From the condition of the extreme point that $\left|Q_{ji}^{″}(\overline{t_i})\right|=0$, we can derive that

$$|Q_{ji}^{\prime }(\overline{t_i})|=\left\{\begin{array}{cc}\hfill \left|\frac{a_{ji}{a}_{j,i+1}}{2\left(a_{ji}-a_{j,i+1}\right)}h+\frac{\left({\theta }_{B(j,i+1)}-{\theta }_{B(j,i)}\right)}{h}+\frac{\left(a_{ji}-a_{j,i+1}\right)h}{6}\right|,\hfill & \mathit{if} a_{ji}\ne a_{j,i+1} \mathit{and} Q_{ji}^{″}\left(t_i\right)Q_{ji}^{″}\left(t_{i+1}\right)\le 0\hfill \\ \hfill 0\hfill & \mathit{if} a_{ji}=a_{j,i+1} \mathit{or} Q_{ji}^{″}\left(t_i\right)Q_{ji}^{″}\left(t_{i+1}\right)>0\hfill \end{array}\right.$$

(16)

Since the joint angular acceleration is a first-order polynomial, and the maximum value of the instantaneous velocity can only be obtained at the endpoint, the joint angular acceleration constraint can be expressed as

$$\max \left|Q_{ji}^{″}(t)\right|=\max \left\{\left|a_{j1}\right|,\left|a_{j2}\right|,\dots ,\left|a_{jn}\right|\right\}\le a{C}_j$$

(17)

In addition, in order to improve the safety performance of the flexible manipulator, it is necessary to limit the amplitude of the manipulator, and the vibration constraint can be expressed as

$$w=\left|{\varphi }_1(L_2)q_1(t)+{\varphi }_2(L_2)q_2(t)\right|\le w_{\max }$$

(18)

3.3. Objective Function

In order to minimize the vibration of the flexible manipulator, the objective function of vibration suppression trajectory planning based on the accumulation of vibration displacement at the end of the flexible manipulator is proposed:

$$f(\theta )={\alpha }_1{\displaystyle {\int }_0^{t_f}}w(\theta ,t)dt+{\alpha }_2{\displaystyle {\int }_{t_f}^{3t_f}w(\theta ,t)dt}$$

(19)

The first half of the objective function is used to measure the vibration deviation of the flexible manipulator during the movement process. The latter part of the objective function is used to measure the residual vibration of the flexible manipulator in time interval $t_f$ to $3t_f$ after the motion stops. ${\alpha }_1$ and ${\alpha }_2$ are weight factors.

4. Optimization Algorithm

4.1. PSO Algorithm

For a d-dimensional optimization problem, suppose that the position and velocity of the particle are ${\mathit{X}}_i=(x_{i,1} x_{i,2} \dots x_{i,d})$, ${\mathit{V}}_i=(v_{i,1} v_{i,2} \dots v_{i,d})$. The particle individual optimal solution is recorded as ${\mathit{P}}_i=(p_{i,1} p_{i,2} \dots p_{i,d})$, and the global optimal solution currently found from the whole group is recorded as ${\mathit{P}}_g=(p_{g,1} p_{g,2} \dots p_{g,d})$. The particles update their position and velocity according to the following rules [35]:

$$\left\{\begin{array}{l}{v}_{i,j}(t+1)=\varpi v_{i,j}(t)+c_1{\varsigma }_1[p_{i,j}-x_{i,j}(t)]+c_2{\varsigma }_2[p_{g,j}-x_{i,j}(t)]\\ x_{i,j}(t+1)=x_{i,j}(t)+v_{i,j}(t+1)\end{array}\right.$$

(20)

where $\varpi$ is the inertia weight, $c_1$ and $c_2$ are positive learning factors, and ${\varsigma }_1$, ${\varsigma }_2$ are random numbers evenly distributed from 0 to 1.

The speed update formula of a PSO algorithm with a shrinkage factor is [36]

$$\left\{\begin{array}{l}{v}_{i,j}(t+1)=\mathsf{\phi }\{\varpi v_{i,j}(t)+c_1{\varsigma }_1[p_{i,j}-x_{i,j}(t)]+c_2{\varsigma }_2[p_{g,j}-x_{i,j}(t)]\}\\ \mathsf{\phi }=\frac{2}{\left|2-C-\sqrt{C^2-4C}\right|}\end{array}\right.$$

(21)

where $C=c_1+c_2$ to ensure the smooth solution of the algorithm.

PSO is an optimization method without variable constraints. For the application of the PSO algorithm in constrained optimization problems, researchers have proposed a series of schemes to solve constrained optimization problems. Among them, the penalty function method is the most widely used because of its simple design and stable calculation results [23,24]. In this paper, the penalty function method is introduced to expand the objective function to deal with the constraints. The new objective function is expressed as

$$F(\theta )=f(\theta )+K_p{\displaystyle \sum _{i=1}^n{\left[\min (0,g_i(x))\right]}^2}$$

(22)

where $K_p$ is the penalty factor, and $g_i(x)$ is the constraint discriminant. When the constraint is satisfied, $g_i(x)=1$; otherwise, $g_i(x)=-1$.

The steps of the particle swarm optimization algorithm for vibration suppression trajectory planning are as follows. Firstly, when the trajectory of the end actuator of the manipulator is known, the joint motion trajectory can be obtained through inverse kinematics, as shown in the dashed blue frame in Figure 4. The initial value of the joint interpolation point increment is generated by tent chaotic mapping, and the value of the discrete trajectory control point is obtained by (4), (6) and (10). Secondly, the discrete trajectory control points are fitted segmentally by using a cubic polynomial to obtain the angular displacement, angular velocity and angular acceleration of the trajectory to be optimized. Next, the dynamic equation of the system is calculated through (2) to obtain the vibration mode of the flexible manipulator. Then, the objective function is calculated according to (22) and the speed is updated according to (21) until the maximum iterative steps are reached. Finally, the optimal variables are output and the optimal trajectory is obtained by cubic polynomial fitting. The entire process is shown in Figure 4.

4.2. Sparrow Search Algorithm (SSA)

SSA is a new heuristic optimization algorithm, the position update strategy of the sparrow algorithm is richer than that of PSO and it has better optimization ability in benchmark functions and certain applications [28]. There are two types of sparrows in sparrow algorithms, namely the discoverer and follower, and the proportion of discoverers in the population is set as PD. Discoverers with better fitness values will give priority to obtaining food in the search process and providing the direction of foraging for the followers. Therefore, the foraging search range of the discoverer is larger than that of the followers. For a d-dimensional optimization problem, the position information of the i-th sparrow can be expressed as ${\mathit{X}}_i=\left[\begin{array}{cccc}{x}_{i1}& x_{i2}& \dots & x_{id}\end{array}\right]$. During each iteration, the discoverer updates the position according to the following rules:

$${\mathit{X}}_i^{(t+1)}=\left\{\begin{array}{c}{\mathit{X}}_i^{(t)}\cdot exp(-\frac{i}{\gamma \cdot itermax}),\hspace{1em}\mathit{if}\hspace{1em}{R}_e<ST\\ {\mathit{X}}_i^{(t)}+\Psi \cdot \mathbf{L},\hspace{1em}\mathit{if}\hspace{1em}{R}_e\ge ST\end{array}\right.$$

(23)

where $itermax$ is the maximum number of iterations, $\gamma \in \left[0,1\right]$ is a random number, $R_e\in \left[0,1\right]$ represents the early warning value, $ST\in \left[0.5,1\right]$ represents the safety value, $\Psi$ is a random number subject to a normal distribution, and $\mathbf{L}$ stands for a 1 × d-dimensional matrix in which each element is l. When $R_e<ST$, there is no predator in the foraging environment, and the discoverer can conduct an extensive search operation; when $R_e\ge ST$, some sparrows spot predators and send an early warning to others, and all sparrows need to fly to a safe area.

During foraging, some followers will always monitor the discoverer. Once they realize that the discoverer has found better food, they will leave to compete with the discoverer. If they win, they can immediately obtain the food of the discoverer and execute the rule as the discoverer (23); otherwise, they continue to execute the rule (24). The position update of participants is described as follows:

$${\mathit{X}}_i^{(t+1)}=\left\{\begin{array}{c}\Psi \cdot exp(\frac{{\mathit{X}}_w-{\mathit{X}}_i^{(t)}}{i^2}),\hspace{1em}if\hspace{1em}i>n/2\hfill \\ {\mathit{X}}_p+\left|{\mathit{X}}_i^{(t)}-{\mathit{X}}_p\right|\cdot {\mathit{E}}^{+}\cdot \mathbf{L},\hspace{1em}if\hspace{1em}i\le n/2\end{array}\right.$$

(24)

where ${\mathit{X}}_p$ is the best position of the current discoverer, and ${\mathit{X}}_w$ is the worst position in the current global situation. $\mathit{E}$ stands for a 1 × d-dimensional matrix in which each element is randomly assigned 1 or −1 and ${\mathit{E}}^{+}={\mathit{E}}^T{(\mathit{E}{\mathit{E}}^T)}^{-1}$. When $i>n/2$, it indicates that the i-th follower has a low fitness value, and it needs to fly to other places to obtain more energy.

In addition, the sparrow algorithm assumes that some sparrows in the population perceive danger. The proportion of these sparrows is set as SD, and their initial positions are randomly generated in the population. If the individual fitness value of these sparrows is greater than the current global optimal fitness value, it means that the sparrow is on the verge of being attacked by predators. When the individual fitness value of these sparrows is equal to the current global optimal fitness value, it means that the sparrow is aware of the danger and needs to be close to other sparrows to minimize the risk of predation. It can be expressed as follows:

$${\mathit{X}}_i^{(t+1)}\left\{\begin{array}{c}{\mathit{X}}_g^{\left(t\right)}+\beta ·\left|{\mathit{X}}_i^{\left(t\right)}-{\mathit{X}}_g^{\left(t\right)}\right|,\hspace{1em}if\hspace{1em}{f}_i>f_g\hfill \\ {\mathit{X}}_i^{\left(t\right)}+K_r·\left(\frac{\left|{\mathit{X}}_i^{\left(t\right)}-{\mathit{X}}_w\right|}{\left(f_i-f_w\right)+\epsilon }\right),\hspace{1em}if\hspace{1em}{f}_i=f_g\hfill \end{array}\right.$$

(25)

where ${\mathit{X}}_g$ is the current global optimal location, $\beta$ is the step control parameter, $K_r\in [-1, 1]$ is a random number indicating the moving direction and step size, and $\epsilon$ is a small constant. $f_i$ is the fitness value of the current sparrow; $f_g$ and $f_w$ are, respectively, the current global best and worst fitness values.

4.3. Improved Sparrow Search Algorithm (ISSA)

In order to improve the global search ability of the sparrow optimization algorithm and reduce the risk of falling into local optima, this paper integrates the sine cosine algorithm (SCA) and elite opposition learning strategy into the sparrow algorithm, and proposes a new, improved sparrow algorithm. The SCA algorithm has a simple structure, few parameter settings and easy implementation, offering certain advantages over PSO and GA in some optimization examples. SCA carries out global exploration and local development according to the oscillation change of the sine cosine model. The position update formula of the standard SCA is as follows [37]:

$${\mathit{X}}_i^{(t+1)}=\left\{\begin{array}{l}{\mathit{X}}_i^{(t)}+{\varsigma }_3\times \sin ({\varsigma }_4)\times \left|{\varsigma }_5{\mathit{X}}_{best}-{\mathit{X}}_i^{(t)}\right|\hspace{1em}{\varsigma }_6<P_r\\ {\mathit{X}}_i^{(t)}+{\varsigma }_3\times \cos ({\varsigma }_4)\times \left|{\varsigma }_5{\mathit{X}}_{best}-{\mathit{X}}_i^{(t)}\right|\hspace{1em}{\varsigma }_6\ge P_r\end{array}\right.$$

(26)

where ${\mathit{X}}_{best}$ is the best individual of the group after the t-th iteration; ${\varsigma }_4$, ${\varsigma }_5$ and ${\varsigma }_6$ are random numbers, and ${\varsigma }_4\in [0,2\pi ]$, ${\varsigma }_5\in [-2,2]$, ${\varsigma }_6\in [0,1]$. $P_r$ is a constant indicating the probability of switching between sine and cosine. ${\varsigma }_3$ is the conversion parameter.

$${\varsigma }_3=a_r{e}^{-\frac{t}{itermax}}$$

(27)

where $a_r$ is a preset constant, $t$ is the current iteration, and $itermax$ is the maximum number of iterations.

The improved sparrow algorithm introduces the sine algorithm into the position update formula of the discoverer of the standard sparrow search algorithm; it is described as follows:

$${\mathit{X}}_i^{(t+1)}=\left\{\begin{array}{c}{\mathit{X}}_i^{(t)}+{\varsigma }_3\times \sin ({\varsigma }_4)\times \left|{\varsigma }_5{\mathit{X}}_{best}-{\mathit{X}}_i^{(t)}\right|,\hspace{1em}if\hspace{1em}{R}_e<ST\\ {\mathit{X}}_i^{(t)}+\Psi \cdot \mathbf{L},\hspace{1em}if\hspace{1em}{R}_e\ge ST\end{array}\right.$$

(28)

The elite opposition-based learning strategy seeks to find the reverse solution for some elite individuals, and then liberate the reverse solution of elite individuals into the population to participate in competition [38]. It is proven that this method can effectively improve the search ability of intelligent optimization algorithms. The improved sparrow algorithm uses the inverse solution of elite points to further update the sparrow position after sparrow optimization. The elite group is represented by ${\mathit{X}}_e$. The proportion of elite groups is $ED$. ${\mathit{X}}_{ei}=(x_{ei,1},x_{ei,2},\dots ,x_{ei,d})\in [\mathit{\ell },\mathit{u}]$ is an elite point in d-dimensional space, and its inverse solution ${\mathit{X}}_{ei}^{\prime }=(x_{i,1}^{\prime },x_{i,2}^{\prime },\dots ,x_{i,d}^{\prime })$ can be expressed as

$${\mathit{X}}_{ei}^{\prime }=k\cdot (\mathit{\ell }+\mathit{u})-{\mathit{X}}_{ei}$$

(29)

where $k$ is a constant. When $f({\mathit{X}}_{ei})>f(X_{ei}^{\prime })$, we take ${\mathit{X}}_{ei}^{\prime }$ as the elite individual of the next iteration to replace ${\mathit{X}}_{ei}$. Then, we calculate the global optimal position of the group and rearrange the new sparrow population.

The steps of trajectory optimization of flexible manipulators by using ISSA are as follows. Firstly, the initial value is generated by tent chaotic mapping, and then the objective function is calculated according to (19). The positions of the discoverer, follower and early warning sparrow are updated, respectively, according to (28), (24) and (25), and then the position of the sparrow is further updated through the opposition-based learning of the elite solution, until the optimization variable that minimizes the objective function of the flexible manipulator is found. Finally, the optimal variables are output and the optimal trajectory is obtained by cubic polynomial fitting. The entire process is shown in Figure 5.

4.4. Simulation

To verify the effectiveness of the proposed ISSA in vibration suppression trajectory optimization, the flexible manipulator shown in Figure 2 is used as an example for optimization simulation, and compared with the traditional PSO algorithm, SSA and OBLSSA. The simulation is carried out in the Matlab2020 environment, using the 10-core parallel operation provided by the National Supercomputing Center in Jinan, China. The size of the manipulators in the simulation is $l_1=0.975 \mathrm{m}, l_2=1 \mathrm{m}$; the density is ${\rho }_1=13.507, {\rho }_2=12.363$; the bending stiffness of the flexible link is $EI=125 \mathrm{N}\cdot {\mathrm{m}}^2$. The parameters of the manipulator are shown in

Table 1. Parameters of underwater flexible manipulator.

Parameter	Symbol	Value
Length (m)	$l$	$l_1=0.975,l_2=1$
Density (kg·m−1)	$\rho$	${\rho }_1=13.507,{\rho }_2=12.363$
Bending stiffness (N·m2)	$EI$	125
Water resistance coefficient	$C_D$	1.1
Additional mass coefficient	$C_M$	1

. In the selection of the initial parameters, we refer to some existing methods [19,26] and study the sensitivity of some important parameters. The selected parameters can ensure that the 4 optimization algorithms have good performance.

4.4.1. Trajectory Planning with Different Gravity Conditions

Firstly, the trajectory optimization effects of the traditional PSO, SSA, OBLSSA and ISSA for the underwater flexible manipulator are studied under two conditions: considering the influence of gravity (the resultant force of gravity and buoyancy is 0.1 times that of gravity) and not considering the influence of gravity (the resultant force of gravity and buoyancy is zero). In this case, the accelerations at the initial and termination moments are known and the accelerations are zero. The starting position is set as ${\mathit{\theta }}_0=$ [0,0], and the target position is set as ${\mathit{\theta }}_f=[\pi /4,\pi /2]$; the weight coefficient of objective function ${\alpha }_2=0$, seven groups of time points are evenly selected, and the optimization parameters are set as $[\Delta {\theta }_{13}, \Delta {\theta }_{14}, \Delta {\theta }_{15}, \Delta {\theta }_{23}, \Delta {\theta }_{24}, \Delta {\theta }_{25}]$. The dynamic equation of the flexible manipulator is solved by the fourth-order Runge–Kutta method, and the number of data points is taken as Num = 500. The corresponding sampling step is taken as $2\pi /499$; the joint angle, velocity and acceleration constraints are ${\left|\theta \right|}_{\max }={\left|\dot{\theta }\right|}_{\max }={\left|\ddot{\theta }\right|}_{\max }=2$; the maximum flexible vibration displacement is 0.1 m. The size of the population is N, and the dimension of the optimization problem is recorded as d; the parameters of several optimization algorithms are shown in

Table 2. Optimization algorithm parameters.

Case	Algorithm	Optimization Algorithm Parameters
Case 1 (without gravity)	PSO	N = 30, d = 6, c1 = 2.05, c2 = 2.05, K = 150, $\varpi =0.9$, Num = 500
	SSA	N = 30, d = 6, ST = 0.7, PD = 0.7, SD = 0.2, K = 150, Num = 500
	OBLSSA	N = 30, d = 6, ST = 0.7, PD = 0.7, SD = 0.2, ED = 0.2, K = 150, Num = 500
	ISSA	N = 30, d = 6, ST = 0.7, PD = 0.7, SD = 0.2, ED = 0.2, $a_r$ = 2, K = 150, Num = 500
Case 1 (with gravity)	PSO	N = 30, d = 6, c1 = 2.05, c2 = 2.05, K = 150, $\varpi =0.9$, Num = 500
	SSA	N = 30, d = 6, ST = 0.7, PD = 0.7, SD = 0.2, K = 150, Num = 500
	OBLSSA	N = 30, d = 6, ST = 0.7, PD = 0.7, SD = 0.2, ED = 0.2, K = 150, Num = 500
	ISSA	N = 30, d = 6, ST = 0.7, PD = 0.7, SD = 0.2, ED = 0.2, $a_r$ = 2, K = 150, Num = 500
Case 2 (${\alpha }_2=1$)	PSO	N = 30, d = 8, c1 = 2.05, c2 = 2.05, K = 150, $\varpi =0.9$, Num = 500
	SSA	N = 30, d = 8, ST = 0.7, PD = 0.7, SD = 0.2, K = 150, Num = 500
	OBLSSA	N = 30, d = 8, ST = 0.7, PD = 0.7, SD = 0.2, ED = 0.2, K = 150, Num = 500
	ISSA	N = 30, d = 8, ST = 0.7, PD = 0.7, SD = 0.2, ED = 0.2, $a_r$ = 2, K = 150, Num = 500
Case 2 (${\alpha }_2=0$)	PSO	N = 30, d = 8, c1 = 2.05, c2 = 2.05, K = 150, $\varpi =0.9$, Num = 500
	SSA	N = 30, d = 8, ST = 0.7, PD = 0.7, SD = 0.2, K = 150, Num = 500
	OBLSSA	N = 30, d = 8, ST = 0.7, PD = 0.7, SD = 0.2, ED = 0.2, K = 150, Num = 500
	ISSA	N = 30, d = 8, ST = 0.7, PD = 0.7, SD = 0.2, ED = 0.2, $a_r$ = 2, K = 150, Num = 500

.

After iterative operation, an optimal trajectory consisting of three cubic polynomial trajectories can be obtained. Figure 6a shows the evolution process of the four optimization algorithms without considering the influence of gravity, and Figure 6b shows the evolution process of the four optimization algorithms considering the influence of gravity. In order to eliminate the influence of accidental factors, each optimization algorithm is implemented three times to take the average. The optimization results are shown in

Table 3. Optimization results.

Case	Algorithm	Optimal Solution			Running Time/s
		Best	Worst	Mean
Case 1 (without gravity)	PSO	0.0096	0.0111	0.0101	66,838
	SSA	0.0162	0.0166	0.0163	16,434
	OBLSSA	0.0108	0.0121	0.0110	19,978
	ISSA	0.0095	0.0095	0.0095	17,518
Case 1 (with gravity)	PSO	0.0622	0.0622	0.0622	65,619
	SSA	0.0617	0.0617	0.0617	19,058
	OBLSSA	0.0630	0.0632	0.0631	22,286
	ISSA	0.0608	0.0612	0.0609	18,673
Case 2 (${\alpha }_2=1$)	PSO	0.0191	0.0191	0.0191	67,459
	SSA	0.0174	0.0174	0.0174	16,893
	OBLSSA	0.0206	0.0208	0.0207	20,974
	ISSA	0.0148	0.0148	0.0148	18,146
Case 2 (${\alpha }_2=0$)	PSO	0.0217	0.0217	0.0217	82,930
	SSA	0.0137	0.0145	0.0142	15,743
	OBLSSA	0.0132	0.0135	0.0133	19,876
	ISSA	0.0112	0.0113	0.0113	22,495

. It can be seen from Figure 6a and Table 3 that, compared with the other three algorithms, SSA falls into the local optimum earlier; the OBLSSA optimization results are improved compared with SSA; the optimization result of PSO is better than that of SSA and OBLSSA, but the optimization time is much longer than that of the other three algorithms; ISSA improves the global search ability, and the search accuracy and convergence speed are better than for the other three optimization algorithms. Combined with Figure 6b and Table 3, it can be seen that the optimization effect of the SSA optimization algorithm is better than that of the PSO optimization algorithm. Under the set parameters, the best fitness values of the optimization trajectories when using the OBLSSA optimization algorithm are not improved compared with the use of SSA. This may be due to the increase in random factors in the process of generating elite solutions and controlling the boundaries of elite solutions, which reduces the accuracy of the optimization results. The ISSA optimization algorithm combines the advantages of the strong global search ability of the sine and cosine algorithm and opposition elite learning to avoid falling into local optima, and it overcomes the influence of random factors. Figure 7 shows the optimized trajectory without considering the influence of gravity. It can be seen from Figure 7 that the interpolation point increment of the ISSA-optimized trajectory is much larger than that of the SSA optimization algorithm. Figure 8 is the optimized trajectory considering the influence of gravity. Compared with Figure 7, it can be seen that the curvature of the optimized trajectory when gravity’s influence is considered is greater than that when gravity’s influence is not considered, which is caused by the inherent vibration characteristics of the manipulator. Figure 9a shows the flexible displacement at the end of the manipulator under each optimized trajectory when the influence of gravity is not considered. Compared with the non-optimized trajectory, the flexible vibration at the end of the manipulator based on the SSA- and ISSA-optimized trajectories is significantly suppressed in the first 4 s, and the total vibration is smaller than that of the non-optimized trajectory. The vibration accumulation and maximum vibration displacement of the manipulator optimized by ISSA are obviously lower than those of the one optimized by SSA. Figure 9b shows the flexible displacement at the end of the manipulator generated by the SSA- and ISSA-optimized trajectories when the influence of gravity is considered. Compared with the non-optimized trajectory, the flexible vibration is effectively suppressed at 1.5–4 s. Due to the large flexible displacement caused by gravity, the vibration displacement reduced by trajectory planning is not obvious compared with the total flexible displacement.

4.4.2. Trajectory Planning with Different Weight Factors

When the acceleration at the initial and terminal times is not zero and is unknown, the residual vibration will interfere with the control of the manipulator after the joint stops moving. In order to study the influence of residual vibration on the optimization of the vibration suppression trajectory, the paper discuss the trajectory optimization effect of PSO, SSA, OBLSSA and ISSA for underwater flexible manipulators in ${\alpha }_2=1$ and ${\alpha }_2=0$. The starting position is set as ${\mathit{\theta }}_0=[0,0]$, and the target position is set as ${\mathit{\theta }}_f=[\pi /4,\pi /2]$. The resultant force of gravity and buoyancy is zero. Seven groups of time points are evenly selected, and the joint angle increment at the middle time point and the initial acceleration of joint2 are taken as the parameters to be optimized. The optimization parameters can be expressed as $[\Delta {\theta }_{13}, \Delta {\theta }_{14}, \Delta {\theta }_{15}, \Delta {\theta }_{23}, \Delta {\theta }_{24}, \Delta {\theta }_{25}, {\ddot{\theta }}_{20}, {\ddot{\theta }}_{2f}]$; the joint angle, velocity and acceleration constraints are ${\left|\theta \right|}_{\max }={\left|\dot{\theta }\right|}_{\max }={\left|\ddot{\theta }\right|}_{\max }=2$; the maximum flexible vibration displacement is 0.1 m.

Figure 10a shows the optimization processes of four optimization algorithms at ${\alpha }_2=1$. Combining these with Table 3, it can be seen that the optimization result of SSA is better than that of the PSO optimization algorithm. The optimization point obtained by the PSO optimization algorithm in the initial stage of optimization exceeds the constraint boundary. Due to the effect of the penalty function, it overcomes the constraints and obtains a good optimization result. The fitness value of the OBLSSA optimization result is lower than that of SSA, which is related to the increase in random factors in the process of generating and controlling the boundary of the elite solution. The ISSA optimization algorithm combines the global search ability of the sine and cosine algorithm and opposition elite learning to avoid falling into the local optimum and overcome the influence of random factors, obtaining better results than other algorithms. Figure 10b shows the optimization processes of four optimization algorithms at ${\alpha }_2=0$, The optimization points obtained by the PSO optimization algorithm in the initial stage of optimization also exceed the constraint boundary. According to Table 3, it can be seen that the optimization results of the three sparrow optimization algorithms are better than those of the PSO. The OBLSSA optimization algorithm reverse elite solution works in overcoming the local optimum. The optimization results are better than those of the SSA optimization algorithm, and ISSA improves the global search ability and avoids falling into the local optimum, and it obtains optimization results that are better than those of other algorithms. By comparing the optimal fitness values obtained by the four optimization algorithms in Table 3, it can be seen that the fitness values of SSA, OBLSSA and ISSA when ${\alpha }_2=1$ are greater than those when ${\alpha }_2=0$. This is because we add the residual vibration after the joint motion stops. The fitness value of the PSO algorithm when ${\alpha }_2=1$ is less than that when ${\alpha }_2=0$; this may be caused by the PSO algorithm falling into the local optimum. Figure 11 and Figure 12, respectively, show the joint optimization trajectories when ${\alpha }_2=1$ and ${\alpha }_2=0$. Figure 13 shows the comparison of the optimization trajectories of ISSA when ${\alpha }_2=1$ and ${\alpha }_2=0$. Compared with the joint optimization trajectory when ${\alpha }_2=0$, the joint optimization trajectory has smaller curvature and a smaller interpolation point increment when ${\alpha }_2=1$. Figure 14 shows the vibration displacement of the manipulator endpoint. It can be seen from Figure 14a that the residual vibration of the trajectory optimized by SSA and ISSA is significantly suppressed at $t=[t_f,3t_f]$. It can be seen from Figure 14a,b that the vibration suppression is mainly within 1–5 s, and the vibration amount of the optimized trajectory by ISSA is less than that by SSA. Although the amplitude of residual vibration is less than that in the process of motion, the accumulation of residual vibration will also have a great impact on the optimization results.

Remark 1. 

Since the initial value of the joint interpolation point increment is automatically generated by tent chaotic mapping, different optimization algorithms will have different initial fitness values, as shown in Figure 6 and Figure 10. The fitness value of the PSO algorithm in Figure 10 in the early iteration is much higher than that of other optimization algorithms, which is caused by the action of the penalty function.

Remark 2. 

The optimized objective function is the accumulation of vibration displacement in the whole optimization process. In Figure 14b, although the vibration accumulation of SSA is smaller than that of ISSA in the first 4 s, the vibration accumulation of the SSA optimization algorithm increases sharply in the following seconds, resulting in a total vibration accumulation greater than that of ISSA.

To further analyze the optimization performance of the four optimization algorithms, two optimization performance indicators, the average optimization percentage $\lambda$ and average relative optimization time $\kappa$, are introduced.

$$\lambda =\frac{1}{4}{\displaystyle \sum _{i=1}^4\frac{F_{n_i}-F_i}{F_{n_i}}}$$

(30)

$$\kappa =\frac{1}{4}{\displaystyle \sum _{i=1}^4\frac{t_i}{t_{I_i}}}$$

(31)

where $F_n$ is the non-optimized objective function value, $F$ is the optimized objective function value, $t$ is the optimization time, and $t_I$ is the optimization time of the ISSA algorithm.

The optimization performance indicators are shown in

Table 4. Optimization performance indicators.

Optimization Algorithm	$\mathit{\lambda }$	$\mathit{\kappa }$
PSO	25.0%	3.68
SSA	29.2%	0.90
OBLSSA	32.6%	1.10
ISSA	42.3%	1.00

. According to Table 4, the ISSA algorithm takes only 1/3.68 of the time of PSO. Compared to PSO, SSA and OBLSSA, the optimization performance is improved by 17.3%, 13.1% and 9.7%, respectively. This is because the PSO algorithm and the SSA algorithm are more likely to fall into the local optimum compared with OBLSSA and ISSA. The introduction of the opposition learning strategy improves the optimization ability of SSA to a certain extent, but random factors are added in the process of generating and controlling the boundary of the elite solution, which will affect the optimization accuracy of OBLSSA. As the application of SA in ISSA overcomes the effects of random factors, compared with the other three optimization algorithms, ISSA has higher optimization accuracy and a faster convergence speed and obtains a good vibration reduction effect.

5. Incremental Kriging-Assisted ISSA

The numerical simulation in Section 4 verifies that the ISSA-based vibration suppression trajectory planning method can obtain better optimization results. Due to the long optimization time, ISSA is difficult to be directly applied in practice. In order to achieve a fast online search of vibration suppression trajectories, neural networks and intelligent algorithms are used as optimal path generators [19], which can generate vibration suppression trajectories in real time after the start and end positions of the joint angle are given. These methods require a large number of optimization results as sample data. To obtain a large amount of optimized trajectory data in a shorter time, this paper uses an incremental Kriging model based on an improved additive point criterion to construct a global surrogate model, and combines the ISSA optimization algorithm to train the model.

5.1. Incremental Learning Kriging Model

Data-driven technology has great potential in engineering applications [39]. The Kriging model is widely used in engineering optimization because of its good robustness and nonlinear approximation capability [40,41]. The Kriging model introduces statistical assumptions to represent the objective function as a stochastic process as follows.

$$Y(\mathit{x})=\mu +Z(\mathit{x})$$

(32)

The following Gaussian correlation function is often used for the Kriging model:

$$R({\mathit{x}}^{(i)},{\mathit{x}}^{(j)},\mathit{\vartheta })=\exp [-{\displaystyle \sum _{i=1}^d{\vartheta }_i{({\mathit{x}}^{(i)}-{\mathit{x}}^{(j)})}^2}]$$

(33)

The estimated value of the model is expressed as

$$\widehat{y}(\mathit{x})=\widehat{\mu }+{\mathit{r}}^T{\mathit{R}}^{-1}(\mathit{y}-\mathbf{1}\widehat{\mu })$$

(34)

The mean square deviation of the estimated value is expressed as

$${\widehat{s}}^2(\mathit{x})={\widehat{\sigma }}^2\left[1-{\mathit{r}}^T{\mathit{R}}^{-\mathbf{1}}\mathit{r}+\frac{{(1-{\mathbf{1}}^T{\mathit{R}}^{-1}\mathit{r})}^2}{{\mathbf{1}}^T{\mathit{R}}^{-\mathbf{1}}\mathbf{1}}\right]$$

(35)

where $\mathit{\vartheta }$, $\mu$, ${\sigma }^2$ are hyperparameters and $\widehat{\mathit{\vartheta }}$, $\widehat{\mu }$, ${\widehat{\sigma }}^2$ are the corresponding maximum likelihood estimations.

The training complexity of the Kriging model mainly comes from the optimization of the hyperparameters, which requires solving the correlation matrix $\mathit{R}$ and its inverse. The Cholesky decomposition can improve the accuracy of the numerical solution. The Cholesky decomposition of the correlation matrix $\mathit{R}$ can be expressed as

$$\mathit{R}=\mathit{L}{\mathit{L}}^T$$

(36)

The training and updating of the Kriging model have $O(n^3)$ complexity. Since each new piece of data is a small amount, it is not necessary to completely retrain the current Kriging model, but only to learn the effect of the new data on the current Kriging model. The incremental Kriging model approach proposed by Zhan and Xing can significantly reduce the computational complexity and greatly shorten the surrogate model construction time by chunking the incremental relationship matrix $\tilde{\mathit{R}}$ with respect to old data and new data [42].

$$\tilde{\mathit{R}}=\left[\begin{array}{cc}\mathit{R}& \mathit{A}\\ {\mathit{A}}^T& \mathit{B}\end{array}\right]$$

(37)

where $\mathit{A}$, $\mathit{B}$ is the incremental matrix of new data. For the new data with q incremental points, we can obtain

$$\tilde{\mathbf{1}}={[\mathbf{1},1_1,\cdots ,1_q]}^T$$

(38)

$$\tilde{\mathit{y}}={[\mathit{y},y^{(n+1)},\cdots ,y^{(n+q)}]}^T$$

(39)

$$\tilde{\mathit{r}}={[\mathit{r},R(\mathit{x},{\mathit{x}}^{(\mathrm{n}+1)};\widehat{\mathit{\vartheta }}),\cdots ,R(\mathit{x},{\mathit{x}}^{(\mathrm{n}+\mathrm{q})};\widehat{\mathit{\vartheta }})]}^T$$

(40)

$$\tilde{\mathit{R}}=\tilde{\mathit{L}}{\tilde{\mathit{L}}}^T$$

(41)

$$\tilde{\mathit{L}}=\left[\begin{array}{cc}\mathit{L}& \mathbf{0}\\ {\mathit{A}}^T{({\mathit{L}}^T)}^{-1}& \mathrm{chol}(\mathit{B}-{\mathit{A}}^T{({\mathit{L}}^T)}^{-1}{\mathit{L}}^{-1}\mathit{A})\end{array}\right]$$

(42)

After obtaining the incremental relationship matrix $\mathit{R}$ and the incremental matrices $\mathit{A}$, $\mathit{B}$, the trained Kriging model can be obtained, and then the predicted values and predicted standard deviations of the unknown points can be obtained by the new Kriging model.

5.2. Infill Criteria

At present, the most commonly used infill criteria include the MSP and EI criteria. The MSP criterion aims to find the optimal solution ${\mathit{x}}^{\ast }$ of the objective function directly on the surrogate model, and then obtain the exact solution of ${\mathit{x}}^{\ast }$ as new sample data. The EI criterion finds the point with the highest probability of improvement of the objective function as a new sample point. The probability value of EI improvement can be expressed as

$$\mathrm{E}\mathrm{I}(\mathit{x})=(y_{\min }-\widehat{y}(\mathit{x}))\Phi (\frac{y_{\min }-\widehat{y}(\mathit{x})}{\widehat{s}(\mathit{x})})+\widehat{s}(\mathit{x})\phi (\frac{y_{\min }-\widehat{y}(\mathit{x})}{\widehat{s}(\mathit{x})})$$

(43)

where $y_{\min }$ is the current optimal response function value, $\Phi$ is the standard normal cumulative distribution function, and $\phi$ is the standard normal probability density function.

Similar methods are used to establish the Kriging model of constraint $g_i(\mathit{x})$, $g_i(\mathit{x})\in N[{\widehat{g}}_i(\mathit{x}),{s_{gi}}^2(\mathit{x})]$. Then, the probability of satisfying the constraint at any position is

$$P[-\zeta \le g_i(\mathit{x})\le \zeta ]=2\Phi (\frac{\zeta -{\widehat{g}}_i(\mathit{x})}{s_{gi}(\mathit{x})})-1$$

(44)

The EI value with constraints is expressed as

$${\mathrm{E}}_{\mathrm{c}}\mathrm{I}(\mathit{x})=\mathrm{E}\mathrm{I}(\mathit{x})\cdot \underset{i=1}{\overset{N_c}{\Pi }}P[-\zeta \le g_i(\mathit{x})\le \zeta ]$$

(45)

Studies have shown that the combination of the MSP criterion and EI criterion has a better optimization effect than a single infill criterion [43]. The IKA-ISSA used in this paper constructs an incremental Kriging model based on the MSP + EI hybrid infill criterion and combines the ISSA optimization algorithm to train the model. The IKA-ISSA optimization process is shown in Algorithm 1.

Algorithm 1. Optimization process of IKA-ISSA.

Parameters: number of sample points Nt, maximum number of sample points Nm, learning increment step of Kriging model $generation$

1. Initial sample set generation: Latin superelevation method to generate N initial design points $sample\_x$. the response value of vibration accumulation $sample\_y1$ and the maximum flexible displacement constraint value $sample\_y2$ obtained by the dynamic model.

2. Generate objective function values: generate objective function y with penalty factors based on response and constraint values and rank $y$ to obtain the optimal value $y\_min$.

3. While Nt < Nm do

4. If $generation$ = = 1

Generating the initial Kriging model: the initial sample set is used to generate the responding Kriging model $kriging\_model1$ and the constrained Kriging model $kriging\_model2$, respectively.

5. Else

Update the Kriging model: generate new Kriging models  $kriging\_model1$, $kriging\_model2$ based on the Kriging model from the previous learning increment step $generation$ and the sample set increments $infill\_x$, $infill\_y1$, $in fill\_y2$.

End if

6. ISSA optimization: set the optimal population P, obtain the predicted value of the population sample based on the initial Kriging model and perform $K$ iterations of ISSA optimization to obtain the optimized population P1.

7. Generate candidate increment set: merge similar individuals in population P1 to generate population P2. Remove points in population P2 that are similar to sample set $sample\_x$ as candidate increment set $P\_candi$.

8. Generate increment sets:  select sample set increments $infill\_x$, $infill\_y1$, $in fill\_y2$ from the candidate increment set $P\_candi$ based on MSP + EI infill criterion and Kriging model.

9. Generate objective function values of infill points: obtain the response value $infill\_y1$ and the constraint value $in fill\_y2$ for the sample set increment through the dynamic model and calculate the objective function $in fill\_y$ for the increment sample set.

10. Update variables: update the sample sets $sample\_x$, $sample\_y1$, $sample\_y2$; update Nt and $generation$; update the optimal value $y\_min$ and the corresponding optimal point $x\_min$.

11. End while

12. Output optimal solution: Output $y\_min$ and $x\_min$.

5.3. Simulation

Taking the example of case 1 (without gravity), six groups of time points are evenly selected, $\mathrm{d}=4$ and the optimization parameters are set as $[\Delta {\theta }_{13}, \Delta {\theta }_{14}, \Delta {\theta }_{23}, \Delta {\theta }_{24}]$. The maximum flexible displacement is selected as the constraint, and five groups of different joint angle initial and termination positions $[{\theta }_{11}, {\theta }_{1t}, {\theta }_{21}, {\theta }_{2t}]$ are used to verify the learning ability of IKA-ISSA under the constraint. The five groups of joint angle initial and termination positions are set as $[0,\pi /4,0,\pi /2]$, $[0,\pi /2,0,\pi /4]$, $[0,\pi /4,\pi /6,\pi /2]$, $[\pi /6,\pi /2,0,\pi /4]$, $[\pi /6,\pi /2,\pi /6,\pi /2]$. In IKA-ISSA, the number of initial sample points $N_K=100$, the maximum number of samples in the sample set $N_m=180$, the size of the sparrow population $N=30$, and the number of sparrow optimization iteration steps $K=10$. The parameter of normal ISSA is the same as that in case 1.

From Figure 15, it can be seen that the vibration accumulation of the optimized trajectory is significantly smaller than that of the non-optimized trajectory, and the optimized result of IKA-ISSA is slightly different from that of ordinary ISSA after a short learning and training period, and even better than that of ordinary ISSA in some cases. Figure 16 shows the optimization times of the two methods. IKA-ISSA takes less than 3% of the time of ISSA, among which the time to generate the initial sample set accounts for a large proportion, which is caused by the complex dynamic calculation. In addition, the learning and training time of IKA-ISSA is much less than that of ISSA optimization training, which indicates that IKA-ISSA integrates the strong nonlinear approximation ability of the incremental Kriging model and the high search efficiency of ISSA.

6. Conclusions

In this paper, ISSA combining an elite strategy and the sine algorithm is proposed for the trajectory planning of underwater flexible manipulators. The simulation results show that ISSA has the advantages of the global search ability of the sine and cosine algorithm and the reverse elite learning ability to avoid falling into local optima. It is more computationally efficient than the PSO algorithm, and the average time consumption is only 1/3.68 of that of the PSO algorithm. In order to solve the problem of the long optimization time of ISSA, this paper further proposes IKA-ISSA, which integrates the strong nonlinear approximation capability of the incremental Kriging model and the high search efficiency of ISSA to greatly reduce the time required for optimization and ensure the accuracy of optimization. This method can obtain a large number of optimization results in a shorter time to construct real-time optimal path generators.