Optimal Dynamics Control in Trajectory Tracking of Industrial Robots Based on Adaptive Gaussian Pseudo-Spectral Algorithm

Abstract

A pseudo-spectral control algorithm based on adaptive Gauss collocation point reconstruction is proposed to efficiently solve the optimal dynamics control problem of industrial robots. A mathematical model for the kinematic relationship and dynamic optimization control of industrial robots has been established. On the basis of deriving the Legendre–Gauss collocation formula, a two-stage adaptive Gauss collocation strategy for industrial robot dynamics control variables was designed to improve the dynamics optimization control effect of industrial robot by improving the solution efficiency of constrained optimization problems. The results show that compared with the control variable parameterization method and the traditional Gaussian pseudo-spectral method, the proposed dynamic optimal control method based on an adaptive Gaussian point reconstruction algorithm can effectively improve the solving time and efficiency of constrained optimization problems, thereby further enhancing the dynamic optimization control and trajectory tracking effect of industrial robots.

1. Introduction

The research on the optimal control of industrial robot dynamics is a cutting-edge interdisciplinary direction in robotics technology, automatic control, industrial control, and other fields. Its research aims to achieve the efficient and stable operation of industrial robots in underwater environments through effective control of their motion states [1,2]. With the continuous development of industrial robot technology, its application in underwater operations, detection, rescue, and other fields is becoming increasingly widespread. In underwater environments, due to the complexity of fluid dynamics, the stability and maneuverability of industrial robots have become crucial factors in the design process. The dynamic performance and control accuracy of industrial robots directly affect the efficiency and safety of task execution. The research on high-performance dynamic optimal control for industrial robots is of great value [3,4].

In terms of optimal control of underwater robot dynamics, the establishment of dynamic models, calculation and verification of hydrodynamic coefficients, motion stability control, and performance analysis are the most important aspects related to operational control effectiveness. During the movement of underwater robots, when subjected to external interference, it is necessary to maintain a stable motion state while achieving rapid tracking of the predetermined trajectory. At this point, underwater robots should have good posture control and trajectory tracking capabilities [5,6]. On the basis of achieving precise modeling and dynamic analysis, high-performance optimized tracking control algorithms are the key to achieving the above objectives. The core goal of optimal control theory is to seek an optimal strategy that enables industrial robots to dynamically adjust their motion state based on predetermined goals and task requirements when facing complex and changing environments. In order to achieve this goal, optimal control theory usually adopts various optimization algorithms and strategies [7,8,9]. Considering multiple factors such as the dynamic characteristics of industrial robots, environmental constraints, and task requirements, the optimal control strategy is found to ensure that industrial robots can achieve their goals at the minimum cost [10,11].

In existing research, some optimization algorithms focus on minimizing the energy consumption of industrial robots. Dynamic programming and Pontryagin’s maximum principle are two commonly used optimization algorithms [12,13,14,15]. Dynamic programming algorithms transform complex multi-stage-decision-making problems into a series of simple subproblems by constructing state transition equations, in order to find the global optimal solution. The Pontryagin maximum principle utilizes the Lagrange multiplier method to transform the constraints into additional terms to the objective function and then obtains the optimal control strategy by solving the optimization equation [16,17]. These two algorithms have their own advantages and disadvantages in practical applications and need to be selected based on specific task requirements and system characteristics. In the design process of energy consumption minimization optimization algorithms, multiple factors need to be considered, such as the dynamic characteristics of the robot and the task constraints. By considering these factors comprehensively, optimization algorithms that are more in line with practical applications can be designed to minimize energy consumption [18,19,20].

In summary, existing research on optimization control algorithms for underwater robots mostly focuses on optimizing and solving a single control objective, while research on achieving multi-objective collaborative optimization [21,22] from a global perspective is still rare. The operating environment of underwater robots is diverse and the dynamic coupling characteristics are complex. Therefore, it is of great significance to conduct research on multi-objective optimization control algorithms while considering the improvement of their comprehensive performance [23,24]. The GPM has been extensively studied and applied in solving optimal control problems in industrial control and other fields due to its advantages of high solving accuracy, low computational complexity, simple structure, and high solving efficiency [25,26]. The idea of this method is to discretize the continuous optimal control problem at orthogonal collocation points and approximate the state variables and control variables through global interpolation polynomials, thereby transforming the optimal control problem into an NLP problem for solution [27,28]. In the field of industrial robot motion control, nonlinear optimization problems have become an important research direction due to the complexity of underwater environments and the dynamic nature of robot systems. Compared to traditional linear optimization problems, nonlinear optimization problems are more complex, so they have stricter requirements for the efficiency of solving optimization control problems [29,30].

However, traditional GPM methods use a globally unified point allocation approach to discretize state variables and control variables. In order to obtain accurate control results, a large number of point allocation points are usually required, which at the same time affects the optimization calculation efficiency and leads to a decrease in the control performance of industrial robots [31,32]. On the other hand, unified positioning can easily cause oscillations in the control curve, thereby affecting the stability of the industrial robot’s driving posture. Therefore, the coordination and smoothing of the problem between increasing discrete point allocation and computational efficiency in the GPM method is an important technological breakthrough direction for improving the dynamic optimization control problem of industrial robots [33,34].

Therefore, in order to achieve an efficient solution to the dynamic optimization control and optimization constraint problems for industrial robots, this study combines the Gauss pseudo-spectral control theory and proposes a multi-objective optimal control method for industrial robots based on an adaptive Gauss point reconstruction pseudo-spectral method, which solves the problem of increasing discrete point allocation and computational efficiency in GPM and improves the comprehensive motion control performance of underwater industrial robots. Firstly, based on the kinematic coordinate transformation relationship and force analysis results of industrial robots, a mathematical model for the multi-degree of freedom dynamic control of industrial robots is established. Then, the Legendre–Gauss collocation solution formula was derived based on the Legendre polynomial, and a two-stage adaptive Gauss collocation reconstruction strategy was proposed for the problem of discrete collocation of control variables. By adaptively optimizing the collocation number of control variables, the solving efficiency was improved. Finally, a comparative simulation test was conducted on the dynamic optimization control problem of industrial robots, and the proposed control method was compared with the CVP algorithm and the traditional GPM method to effectively reflect the effectiveness of the proposed method.

The rest of this paper is organized as follows. The dynamic modeling of industrial robots is presented in Section 2. The dynamic optimal control of industrial robots is designed in Section 3. The verification and analysis are provided in Section 4, and followed by the conclusion in Section 5.

2. Dynamic Modeling of Industrial Robot

2.1. Motion Coordinate System and Coordinate Transformation of Robot

To achieve dynamic modeling and optimal control of industrial robots, two orthogonal coordinate systems are established, including the geodetic coordinate system and the local coordinate system. As shown in Figure 1, in the geodetic coordinate system, with the starting point of motion as the origin, the Ox and Oy axes are in the horizontal plane, and the Oz axis is vertically downward along the vertical direction. To facilitate the description of the rigid body motion behavior of the robot, a local coordinate system for the industrial robot was established with the center of mass of the robot as the origin of the coordinate system. Among them, the Ox axis is along the longitudinal axis of the industrial robot, the Oy axis is perpendicular to the Ox axis and points towards the starboard direction, and the Oz axis is perpendicular to the plane of the Oxy axis, with its direction determined according to the right-hand rule. On the basis of establishing a coordinate system, the velocity and angular velocity transformation matrices of the industrial robot from the local coordinate system to the geodetic coordinate system can be obtained as follows:

$T_v=\left(\begin{array}{ccc}\cos \phi \cos \theta & -\cos \varphi \sin \phi +\cos \phi \sin \varphi \sin \theta & \sin \phi \sin \varphi +\cos \phi \cos \varphi \sin \theta \\ \sin \phi \cos \theta & \cos \phi \cos \varphi +\sin \phi \sin \varphi \sin \theta & -\cos \phi \sin \varphi +\cos \varphi \sin \phi \sin \theta \\ -\sin \theta & \cos \theta \sin \varphi & \cos \theta \cos \varphi \end{array}\right)$

$T_{\omega }=\left(\begin{array}{ccc}1& \tan \theta \sin \varphi & \tan \theta \cos \varphi \\ 0& \cos \varphi & -\sin \varphi \\ 0& \sin \varphi /\cos \theta & \cos \varphi /\cos \theta \end{array}\right)$

where $u^{\prime }={\left(\begin{array}{ccc}u& v& w\end{array}\right)}^{\mathrm{T}}$ and $\mathrm{\Omega }={(\begin{array}{ccc}p& q& r\end{array})}^{\mathrm{T}}$.

2.2. Force Analysis of Robot

In order to achieve optimal dynamic control of an industrial robot in an underwater environment, it is necessary to first establish a dynamic model of the industrial robot. The robot system moving underwater is a nonlinear dynamic system that requires a large number of parameters to be determined. Due to limitations in technology and testing conditions, some parameters cannot be accurately measured or cannot be determined. In order to meet the needs of control, it is necessary to simplify the system and only consider the factors that have a major impact on system performance, such as gravity, buoyancy, thrust, and hydrodynamics. Firstly, the industrial robots are subject to the gravitational force of the Earth, and the resulting forces and moments need to be reflected in the local coordinate system, which can be expressed as

(1)${\left(\begin{array}{c}F\\ \tau \end{array}\right)}_G=G{\left(\begin{array}{cccccc}-\sin \theta & \cos \theta \sin \phi & \cos \theta \cos \phi & 0& 0& 0\end{array}\right)}^T,$

where due to the coincidence of the robot’s center of mass with the origin of the coordinate system, its gravitational moment is 0. At the same time, industrial robots are subjected to buoyancy in the water, and the resulting buoyancy and buoyancy moment can be represented in the local coordinate system as

(2)${\left(\begin{array}{c}F\\ \tau \end{array}\right)}_B=-B\left(\begin{array}{c}-\sin \theta \\ \cos \theta \sin \phi \\ \cos \theta \cos \phi \\ z_{\mathrm{b}}\cos \theta \sin \phi -y_{\mathrm{b}}\cos \theta \cos \phi \\ x_{\mathrm{b}}\cos \theta \cos \phi +z_{\mathrm{b}}\sin \theta \\ -y_{\mathrm{b}}\sin \theta -x_{\mathrm{b}}\cos \theta \sin \phi \end{array}\right),$

where ${\left(\begin{array}{ccc}{x}_{\mathrm{b}}& y_{\mathrm{b}}& z_{\mathrm{b}}\end{array}\right)}^{\mathrm{T}}$ is the coordinate of the center of buoyancy in the local coordinate system. In addition, it is necessary to consider the driving force and torque exerted on industrial robots. As Figure 1 shows, the direction of the thrust force is indicated by arrows, where the thrust forces F_T1 and F_T2 are located in the oxy plane and symmetrical with respect to the x-axis, achieving rotation around the z-axis and translation along the x-axis. The thrust forces F_T3 and F_T4 are located in the oyz plane and are symmetrical with respect to the y-axis, achieving rotation around the x-axis and translation along the y-axis. The thrust forces F_T5 and F_T6 are located in the oxz plane and are symmetrical with respect to the z-axis, achieving rotation around the y-axis and translation along the z-axis.

Assuming the resultant forces along the x, y, and z axes are T_x, T_y, and T_z, respectively, and the resultant moments acting on the x, y, and z axes are K_T, M_T, and N_T, respectively, the resultant forces and resultant moments on each axis can be represented as $T_x=F_{T1}+F_{T2}$, $T_x=F_{T1}+F_{T2}$, $T_z=F_{T5}+F_{T6}$, $K_T=R_{T3}{F}_{T3}+R_{T4}{F}_{T4}$, $M_T=R_{T5}{F}_{T5}+R_{T6}{F}_{T6}$, $N_T=R_{T1}{F}_{T1}+R_{T2}{F}_{T2}$.

Robots also experience hydrodynamic forces in water, and the resulting forces and moments need to be reflected in the local coordinate system. Let the velocity of water flow in the geodetic coordinate system be $u_f={\left(\begin{array}{ccc}{u}_f& v_f& w_f\end{array}\right)}^T$, and the transformation relationship between the velocity of water flow in the local coordinate system and the geodetic coordinate system is

(3)$u_f^{\prime }={\left(\begin{array}{ccc}{u}_1& v_{\mathrm{l}}& w_{\mathrm{l}}\end{array}\right)}^T=T_v^{-1}{\left(\begin{array}{ccc}{u}_f& v_f& w_f\end{array}\right)}^T.$

The velocity of the water flow relative to the robot can be expressed as

(4)$u_r^{\prime }={\left(\begin{array}{ccc}{u}_r& v_r& w_r\end{array}\right)}^T=u^{\prime }-u_1^{\prime }.$

The water flow resistance is directly proportional to the square of the relative water flow velocity of the industrial robot. Let the water flow resistance along the x, y, and z axes be F_x, F_y, and F_z, respectively, and they can be expressed as $F_x=0.5\rho C_{\mathrm{d}}{S}_{x^{\prime }}\left|u_{\mathrm{r}}\right|u_{\mathrm{r}}$, $F_y=0.5\rho C_{\mathrm{d}}{S}_{y^{\prime }}\left|v_{\mathrm{r}}\right|v_{\mathrm{r}}$,$F_z=0.5\rho C_{\mathrm{d}}{S}_{z^{\prime }}\left|w_{\mathrm{r}}\right|w_{\mathrm{r}}$. In the formula, $\rho$ is the water density, $C_{\mathrm{d}}$ is the dimensionless drag coefficient, and $S_{x^{\prime }}$,$S_{y^{\prime }}$, and $S_{z^{\prime }}$ are the cross-sectional areas of the industrial robot perpendicular to the x, y, and z axes, respectively. The resistance torque generated by water flow is proportional to the square of the robot’s angular velocity. Based on the resistance torque coefficient, $K_x,K_y,K_z$, the projections of the resistance torque generated by the fluid on the x, y, and z axes of the industrial robot can be obtained as $K_{x^{\prime }}\left|p\right|p,K_{y^{\prime }}\left|q\right|q,K_{z^{\prime }}\left|r\right|r$.

At the same time, when an industrial robot accelerates in a fluid, it also exerts a force on the surrounding fluid, causing acceleration. The fluid mass surrounding it that accelerates is called additional mass. It is approximately believed that industrial robots are symmetrical in terms of front, back, left, and right directions. Therefore, the additional mass of the robot along the three axes is only related to the acceleration and angular acceleration in their respective axis directions. Among them, the additional masses generated along the x, y, and z axis directions are $m_u$, $m_v$, and $m_w$, respectively, and the additional masses generated along the x, y, and z axis rotation directions are $m_p$, $m_q$, and $m_r$, respectively.

2.3. Assumptions and Description of Optimal Control Problem

Based on the establishment of the coordinate system and force analysis mentioned above, according to Newton’s second law and the moment of momentum theorem, the matrix representation of the dynamic equation of the industrial robot can be derived

(5)$M(h)\ddot{h}+N(h,\dot{h})=I(h)F,$

where $\dot{h}={\left(\begin{array}{cccccc}u& v& w& p& q& r\end{array}\right)}^{\mathrm{T}}$, $\ddot{h}={(\dot{u}\hspace{1em}\dot{v}\hspace{1em}\dot{w}\hspace{1em}\dot{p}\hspace{1em}\dot{q}\hspace{1em}\dot{r})}^{\mathrm{T}}$, $F={\left(\begin{array}{cccccc}{T}_x& T_y& T_z& K_T& M_T& N_T\end{array}\right)}^{\mathrm{T}}$,$N(h,\dot{h})={\left(\begin{array}{cccccc}N(1)& N(2)& N(3)& N(4)& N(5)& N(6)\end{array}\right)}^T$, $N(1)=F_x-(m+m_u)[(q{w}_1-v_{1}r)-(qw-vr)]$,$N(2)=F_y-(m+m_v)[(r{u}_1-w_{1}p)-(ru-wp)]$, $N(3)=F_z-(m+m_v)[(p{v}_1-u_{1}q)-(pv-uq)]$, $N(4)=K_{x^{\prime }}\left|p\right|p+(z_{\mathrm{b}}\cos \theta \sin \phi -y_{\mathrm{b}}\cos \theta \cos \phi )B-m{y}_{\mathrm{b}}(p{v}_1-u_{1}q)+m{z}_{\mathrm{b}}(u_{1}r-p{w}_1)$, $N(5)=K_y\left|q\right|q+(x_{\mathrm{b}}\cos \theta \cos \phi +z_{\mathrm{b}}\sin \theta )B+m{x}_{\mathrm{b}}(p{v}_1-u_{1}q)-m{z}_{\mathrm{b}}(w_{1}q-v_{1}r)$, $N(6)=K_z\left|r\right|r-(y_{\mathrm{b}}\sin \theta +x_{\mathrm{b}}\cos \theta \sin \phi )B-m{x}_{\mathrm{b}}(r{u}_1-w_{1}p)+m{y}_{\mathrm{b}}(w_{1}q-v_{1}r)$, $M(h)=diag\left(\begin{array}{cccccc}{m}_u+m& m_v+m& m_w+m& m_p+I_x& m_q+I_y& m_r+I_z\end{array}\right)$, $I(h)$ is the unit matrix, and m is the mass of industrial robot.

By solving the dynamic equations, the motion laws described by the robot in the local coordinate system (i.e., u, v, w, p, q, and r) can be obtained, and then converted to the geodetic coordinate system to obtain the motion laws described by the robot in the geodetic coordinate system

(6)${\dot{x}}_h=T{\left(\begin{array}{cccccc}\dot{x}& \dot{y}& \dot{z}& \dot{\varphi }& \dot{\theta }& \dot{\phi }\end{array}\right)}^T=T{\left(\begin{array}{cccccc}u& v& w& p& q& r\end{array}\right)}^T,$

where $T=\left(\begin{array}{cc}{T}_v& 0\\ 0& T_{\omega }\end{array}\right)$. $\dot{x}$,$\dot{y}$,$\dot{z}$ are the velocity of the robot along three axes in the geodetic coordinate system.

In order to ensure that the robot can track the preset trajectory as soon as possible and smoothly and avoid the phenomenon of sudden changes in the control variables of the control system that affect the continuity of the control variables, the deviation of the system state variables and the control increment are added to the objective function. In actual robot tracking and control systems, in order to ensure the smooth operation of the robot, it is necessary to limit the speed and speed increment of the robot to a certain range of variation. Therefore, it is necessary to increase the constraints on the control quantity and control increment, which are $\left\{\begin{array}{l}{u}_{\min }^{\prime }\le u^{\prime }\le u_{\max }^{\prime }\\ {\mathrm{\Omega }}_{\min }\le \mathrm{\Omega }\le {\mathrm{\Omega }}_{\max }\end{array}\right.$ and $\left\{\begin{array}{l}\mathrm{\Delta }{F}_{\min }\le \mathrm{\Delta }F\le \mathrm{\Delta }{F}_{\max }\\ \mathrm{\Delta }{\tau }_{\min }\le \mathrm{\Delta }\tau \le \mathrm{\Delta }{\tau }_{\max }\end{array}\right.$, respectively.

In the actual tasks of industrial robots, in order to shorten the attitude maneuvering time, it is inevitable to consume more energy. Therefore, it is necessary to consider both the performance improvement requirements in terms of maneuvering time and energy, and design a control system that balances efficiency and energy comprehensive indicators. Therefore, the optimal control problem for trajectory tracking of industrial robots can be expressed as: finding the control input, under certain constraint conditions, moving from the given initial state x(t₀) = x₀, t₀ = 0 to the desired end state x(t_f) = x_f, and minimizing energy consumption. It is defined that $x_J={\left(\begin{array}{cc}h& x_h\end{array}\right)}^T$ and $u_J=F$, so trajectory tracking optimal control is essentially an energy optimal functional extremum problem. Therefore, the performance indicators that balance optimal energy and efficiency can be expressed as

(7)$J={\displaystyle {\int }_{t_0}^{t_f}F}dt+{\displaystyle {\int }_{t_0}^{t_f}dt}.$

At the same time, the initial value constraint that needs to be satisfied is x (t₀) = x₀, the state constraint is $x_{\min }\le x_J\le x_{\max }$, the control constraint is $u_{\min }\le u_J\le u_{\max }$, and the boundary constraint is $t\le t_f$ and $x\left(t_f\right)\le x_f$. Nonlinear optimal control problem refers to the optimal control problem with nonlinear function terms in performance indicators, state equations, or constraint conditions. The trajectory tracking control problem of industrial robots can be equivalent to the optimal control solution problem of the following optimization objective function

(8)$\min J=\varphi \left(x_0,t_0,x_f,t_f\right)+{\displaystyle {\int }_{t_0}^{t_f}g\left(x\left(t\right),u\left(t\right),t\right)dt}.$

3. Dynamic Optimal Control of Industrial Robot in Trajectory Tracking

3.1. Gauss Pseudo-Spectral Method

The Gauss pseudo-spectral method is used to solve the trajectory tracking optimization control problem of industrial robots. The main process is to discretize the state and control variables at Legendre–Gauss points and then approximate the state and control variables through Lagrange interpolation polynomials, thereby converting the system state equation and performance index function integral terms into algebraic operations. Finally, the optimization manipulation problem is transformed into an NLP problem with node state and control variables as the parameters to be optimized.

3.1.1. Lagrange Interpolation Polynomial Approximation

In this work, the discrete approximation of state variables and control variables is completed using the Lagrange interpolation method. Given N + 1 sets of discrete points and their corresponding function values are $\{({\tau }_1,y_1),({\tau }_2,y_2),\cdots ,({\tau }_{N+1},y_{N+1})\}$, the variable can be obtained through N-degree polynomial interpolation

(9)$y(\tau )\approx Y(\tau )=\sum _{i=1}^{N+1}{\widehat{L}}_i(\tau )y_i,$

where $Y(\tau )$ is the N-degree polynomial approximation, ${\widehat{L}}_i(\tau )$ is the N-degree Lagrange interpolation polynomial and which can be calculated using the following formula: ${\widehat{L}}_i(\tau )={\displaystyle \prod _{j=1,j\ne i}^{N+1}}{\displaystyle \frac{\tau -{\tau }_j}{{\tau }_i-{\tau }_j}}$. It is defined that ${\widehat{L}}_i({\tau }_j)={\delta }_{ij}$, and analysis show that the Lagrange interpolation polynomials have the following characteristics: ${\widehat{L}}_i({\tau }_j)={\delta }_{ij}=\left\{\begin{array}{c}1,\hspace{0.33em}\hspace{0.33em}i=j\\ 0,\hspace{0.33em}\hspace{0.33em}i\ne j\end{array}\right.$.

Based on the above equation, it can be seen that the interpolation polynomial satisfies $Y({\tau }_i)=y_i,i=1,2,\cdots ,N+1$ at discrete points, so it is necessary to select appropriate discrete points to obtain the values of $y_i$. However, traditional methods generally divide the time interval into N equal segments and approximate them using a linear combination of Lagrange interpolation polynomials, but this equidistant separation can easily cause the Runge phenomenon. Therefore, the Legendre polynomial with a better approximation effect is selected for discrete point selection in the [–1, 1] interval, and its Nth Legendre polynomial is

(10)$\begin{array}{l}{P}_{n+1}\left(\tau \right)=(\tau -{\alpha }_n)P_n\left(\tau \right)-{\beta }_n^2{P}_{n-1}\left(\tau \right)\\ P_0\left(\tau \right)=1,P_{-1}\left(\tau \right)=0,n=0,1,\cdots ,N\end{array}.$

Furthermore, by calculating the solution of Equation (10), N discrete points located within the open interval (–1,1) can be obtained, denoted as ${\tau }_i$(i = 1, 2, …, N), which are called Legendre–Gauss (LG) collocation points. Analysis shows that it is the root of Equation (8). Given the Legendre polynomial of degree N is $\begin{array}{l}{P}_{n+1}\left(\tau \right)=(\tau -{\alpha }_n)P_n\left(\tau \right)-{\beta }_n^2{P}_{n-1}\left(\tau \right) \\ P_0\left(\tau \right)=1,P_{-1}\left(\tau \right)=0,n=0,1,\cdots ,N\end{array}$, its roots can be solved by the eigenvalues of matrix H, where $H=\left[\begin{array}{ccccc}{\alpha }_0& {\beta }_1& & & \\ {\beta }_1& {\alpha }_1& {\beta }_2& & \\ & & \ddots & & \\ & & & {\alpha }_{n-2}& {\beta }_{n-1}\\ & & & {\beta }_{n-1}& {\alpha }_{n-1}\end{array}\right]$.

Proof. 

According to the definition, the roots of an Nth Legendre polynomial are the solutions of the following equation

(11)$P_{n+1}\left(\tau \right)=\left(\tau -{\alpha }_n\right)P_n\left(\tau \right)-{\beta }_n^2{P}_{n-1}\left(\tau \right)=0,$

when n = 0, there is $P_1\left(\tau \right)=(\tau -{\alpha }_0)P_0\left(\tau \right)-{\beta }_0^2{P}_{-1}\left(\tau \right)$. And according to $P_{-1}\left(\tau \right)=0$ and $P_0\left(\tau \right)=1$, so it can be obtained that $P_1(\tau )=(\tau -{\alpha }_0)$. When n = 1, there is

(12)$P_2\left(\tau \right)=\left(\tau -{\alpha }_1\right)P_1\left(\tau \right)-{\beta }_1^2{P}_0\left(\tau \right).$

According to $P_0\left(\tau \right)=1$, we have $P_2\left(\tau \right)=(\tau -{\alpha }_1)P_1\left(\tau \right)-{\beta }_1^2$, thus the root of $P_2\left(\tau \right)=0$ is the solution of $(\tau -{\alpha }_1)(\tau -{\alpha }_0)-{\beta }_1^2=0$. It is defined that $P_0\left(\tau \right)=1$ and $\left.A_1=\left[\begin{array}{cc}{A}_0& {\beta }_1\\ {\beta }_1& {\alpha }_1\end{array}\right.\right]$, thus, the eigenvalue of A₁ is $|\lambda I-A_1|=(\lambda -A_0)(\lambda -{\alpha }_1)-{\beta }_1^2=0$.

Analysis shows that the root of $P_2\left(\tau \right)=0$ is the eigenvalue of matrix A₁. Similarly, the root of $P_3\left(\tau \right)=0$ is the solution of $(\tau -\alpha 2)|\tau I-A_1|-{\beta }_2^2|\tau I-A_0|=0$. So, according to $\left.A_2=\left[\begin{array}{cc}{A}_1& {\beta }_2\\ {\beta }_2& {\alpha }_2\end{array}\right.\right]$, the eigenvalue of A₂ is

(13)$|\lambda I-A_2|=\left|\begin{array}{cc}\lambda I-A_1& -{\beta }_2\\ -{\beta }_2& \lambda -{\alpha }_2\end{array}\right|=\left|\begin{array}{ccc}\lambda -{\alpha }_0& -{\beta }_1& \\ -{\beta }_1& \lambda -{\alpha }_1& -{\beta }_2\\ & -{\beta }_2& \lambda -{\alpha }_2\end{array}\right|=\lambda I-A_0|=0.$

Therefore, the root of $P_3\left(\tau \right)=0$ is the eigenvalue of matrix A₂. By recursion, it can be concluded that the roots of an Nth Legendre polynomial are the eigenvalues of matrix H. □

3.1.2. Discretization of State and Control Variables

Combining the Lagrange interpolation method from the previous section, Legendre polynomials are used to discretize and approximate the state variables and control variables. Firstly, a time variable is introduced as $\tau$, and the time interval is scaled in Equation (8) using the formula $t={\displaystyle \frac{t_f-t_0}{2}}\tau +{\displaystyle \frac{t_f+t_0}{2}}$. In this way, the time interval [t₀, t_f] is transformed into the [–1, 1] interval. By substituting the scale transformation result into Equation (8), the optimized control problem of the transformed industrial robot is obtained as follows:

(14)$\begin{array}{l}\min J=\varphi (x_0,t_0,x_f,t_f)+{\displaystyle \frac{t_f-t_0}{2}}{\displaystyle {\int }_{-1}^{1}g(x(\tau ),u(\tau ),\tau )d\tau }\\ \mathrm{s}.\mathrm{t}.{\displaystyle \frac{dx(\tau )}{d\tau }}={\displaystyle \frac{t_f-t_0}{2}}f(x(\tau ),u(\tau ),\tau )\\ E(x(-1),t_0,x(1),t_f)=0\\ C(x(\tau ),u(\tau ),\tau )\le 0,\tau \in [-1,1]\end{array}.$

Furthermore, by performing Lagrange interpolation approximation on the state variables at the LG point, we can obtain

(15)$\begin{array}{c}x(\tau )\approx X(\tau )={\displaystyle \sum _{i=1}^{N+1}}{L}_i(\tau )X({\tau }_i)={\displaystyle \sum _{i=1}^{N+1}}{L}_i(\tau )X_i\\ L_i(\tau )={\displaystyle \prod _{j=1,j\ne i}^{N+1}}{\displaystyle \frac{\tau -{\tau }_j}{{\tau }_i-{\tau }_j}}={\displaystyle \frac{b(\tau )}{(\tau -{\tau }_i)\dot{b}({\tau }_i)}}\end{array},$

where $b(\tau )={\displaystyle \prod _{i=1}^{N+1}}(\tau -{\tau }_i)$. Meanwhile, by discretizing the derivative of the state variable, thus it can be obtained that $\dot{x}(\tau )\approx \dot{X}(\tau )={\displaystyle \sum _{i=1}^{N+1}}{\dot{L}}_i(\tau )X_i$. Then, only considering the derivative values of the state variables at the LG distribution point and defining that $\tau ={\tau }_k$, and substituting it into equation (16), it can be obtained that

(16)$\dot{x}({\tau }_k)\approx \dot{X}({\tau }_k)=\sum _{i=1}^{N+1}{\dot{L}}_i({\tau }_k)X_i=\sum _{i=1}^{N+1}{D}_{k,i}{X}_i,$

(17) $D_{k,i}={\dot{L}}_i({\tau }_k)=\left\{\begin{array}{c}{\displaystyle \frac{\dot{b}({\tau }_k)}{({\tau }_k-{\tau }_i)\dot{b}({\tau }_i)}},\hspace{0.33em}\hspace{0.33em}k\ne i\\ {\displaystyle \frac{\dot{b}({\tau }_k)}{2\dot{b}({\tau }_k)}},\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}k=i\end{array}\right..$

Similarly, by approximating the control variable with LG point Lagrange interpolation over the time interval [–1, 1], we can obtain

(18)$u(\tau )\approx U(\tau )=\sum _{i=1}^{N+1}{L}_i(\tau )U({\tau }_i)=\sum _{i=1}^{N+1}{L}_i(\tau )U_i.$

Furthermore, the state equation can be replaced by the following constraint conditions

(19)$\sum _{i=1}^{N+1}{D}_{k,i}{X}_i={\displaystyle \frac{t_f-t_0}{2}}f(X_k,U_k,{\tau }_k).$

Meanwhile, the integral term of the Lagrange term can be approximated as

(20)${\int }_{-1}^{1}g(x(\tau ),u(\tau ),\tau )\mathrm{d}\tau \approx \sum _{k=1}^{N+1}{w}_{k}g(x({\tau }_k),u({\tau }_k),{\tau }_k),$

where $w_k$ is the integral weight in the Gauss integral formula. Ultimately, problem (14) is transformed into the following form of NLP solution problem

(21)$\begin{array}{l}\min J=\varphi (X_0,t_0,X_f,t_f)+{\displaystyle \frac{t_f-t_0}{2}}{\displaystyle \sum _{k=1}^{N+1}}{w}_{k}g(X_k,U_k,{\tau }_k)\\ \mathrm{s}.\mathrm{t}.{\displaystyle \sum _{i=1}^{N+1}}{D}_{k,i}{X}_i-{\displaystyle \frac{t_f-t_0}{2}}f(X_k,U_k,{\tau }_k)=0\\ E(X_1,t_0,X_f,t_f)=0\\ C(X_k,U_k,{\tau }_k)\le 0\\ X_f-X_1-{\displaystyle \sum _{i=1}^{N+1}}{X}_i{\displaystyle \sum _{i=2}^{N+1}}{w}_i{D}_{k,i}=0\end{array}.$

Analysis shows that after LG point discretization, problem (24) is a finite dimensional NLP problem. Therefore, gradient based nonlinear optimization methods such as SQP and interior point method are used to solve this problem.

3.2. Adaptive Gauss Collocation Reconstruction Method

In the traditional Gaussian pseudo-spectral method, the state variables and control variables are both interpolated using the same LG coordinates for polynomial interpolation. Therefore, the number of parameters that need to be optimized in NLP problem (21) are

(22)$NUM=(r_u+r_x)\times N,$

where $r_x$ and $r_u$, respectively, represent the dimensions of state variables and control variables. To ensure the accuracy of solving optimization problems, it is usually necessary to increase the number of LG allocation points. However, significantly increasing the number of allocation points will exponentially increase the number of optimization parameters NUM, reduce computational efficiency, and thus affect the application of GPM method in locomotive optimization operation. In addition, unified point allocation can easily cause oscillations in the control curve, especially for optimization problems involving jump nodes. Therefore, in order to solve the above problems, a two-stage adaptive Gauss point reconstruction strategy is proposed. In the first stage, the control variables are divided into LG distribution points to reduce the number of optimization parameters and perform global optimization calculations. In the second stage, introduced slope change adaptive reconstruction control variable allocation. The two-stage adaptive Gauss point reconstruction process is shown in Figure 2.

In the first stage, the state variables are discretized using Equation (17) for LG allocation points, and the control variables are discretized using pairwise Gaussian allocation points over the time interval [–1, 1], that is

(23)$U_{2l-1}=U_{2l},l=1,2,\cdots ,{\displaystyle \frac{N+1}{2}}.$

Among them, N is an odd number. Thus, the number of variable allocation points are reduced to (N + 1)/2, and the approximation formula is

(24)$u(\tau )\approx U(\tau )=\sum _{i=1}^{N+1}{L}_i(\tau )U({\tau }_i)=\sum _{l=1}^{{\displaystyle \frac{N+1}{2}}}\sum _{i=1}^{2l}{L}_i(\tau )U_l.$

After substituting the control variables obtained from the allocation points into Equation (21), it can be analyzed that the number of optimized parameters is reduced to $(r_u/2+r_x)\times N$.

In the second stage, based on the solution in the first stage, the control variable allocation is reconstructed. Firstly, define the slope change of the $\tau \mathrm{th}$ control variable at the collocation point ${\tau }_{2l+1}$ as

(25)$S_{r,l+1}={\displaystyle \frac{U_{l+1}^r-U_l^r}{{\tau }_{2l+1}-{\tau }_{2l-1}}},l=1,2,\cdots ,{\displaystyle \frac{N+1}{2}}-1.$

Among them $U_{l+1}^r$ and $U_l^r$, respectively, represent the values of the rth control variable at the collocation points ${\tau }_{2l+1}$ and ${\tau }_{2l-1}$. Set the upper and lower threshold parameters for slope changes as ${\epsilon }_{\mathrm{U}}$ and ${\epsilon }_L$, respectively, and then perform segmentation refinement and point merging judgments separately. If $S_{r,l+1}\ge {\epsilon }_{\mathrm{U}}$, the control variables $U_{l+1}^r$ and $U_l^r$ will form a mutation, which indicate that the number of control variable allocation points are insufficient to support the control variable’s requirements for control continuity and smoothness. At this point, it is necessary to refine the number of control variable allocation points in order to increase the number of allocation points between interval $[{\tau }_{2l-1},{\tau }_{2l+1}]$. By increasing the number of control points, the slope of the incremental change in the control variable can be adjusted to achieve smoothness in the motion control of underwater robots and suppress control fluctuations. Similarly, if $S_{r,l+1}\le {\epsilon }_{\mathrm{L}}$, then the control variables $U_{l+1}^r$ and $U_l^r$ are smooth in the interval $[{\tau }_{2l-1},{\tau }_{2l+1}]$. At this time, by setting $U_{l+1}^r=U_l^r$, it is possible to reduce the number of control points in smooth regions appropriately to reduce computational complexity and improve computational efficiency. Finally, the number of control variable allocation points after refinement and merging is obtained through adaptive iteration. Furthermore, after the two-stage adaptive Gauss collocation reconstruction, the LG collocation points of the state variables and control variables in the time interval [–1, 1] can be obtained separately, thereby achieving adaptive adjustment of the collocation points. The pseudocode of the adaptive Gaussian pseudo-spectral algorithm is shown in Table 1.

4. Verification and Analysis

The underwater robot in the simulation example [35] has a body length of 1.33 m, an inner diameter of 0.36 m, a mass of 31.41 kg, and a rotational inertia of I_z = 1.77 × 10⁻¹ kg·m² around the body coordinate system, I_y = I_z = 3.45 kg·m². The length of the y-axis and z-axis from the driver to the connected coordinate system is 0.075 m. In addition, assuming that the center of gravity and center of gravity of the robot coincide, the gravity and buoyancy of the robot are equal, and the density of water is 1000 kg/m³. Using a control moment gyroscope for attitude control of underwater vehicles, in order to avoid the loss of steering effect and reduced maneuverability caused by zero or low speed of underwater vehicles. The parameters of the control torque gyroscope are as follows: the radius of the control torque gyroscope rotor is 0.048 m, the width is 0.05 m, and the installation angle is 54.7°. The constant speed of a single control torque gyroscope rotor is about 15,000 r/min, with a moment of inertia of 3.25 × 10⁻³ kg·m², and the angular momentum of a single control torque gyroscope rotor is 5 kg·m²/s. The structural parameters of the pyramid shaped gyroscope group are as follows: the bottom layout width is 0.24 m, the frame height is 0.15 m, and it can be placed inside the underwater robot to achieve attitude control of the underwater robot. Due to the fact that the driving method used in this design relies on the interaction between the driver and water, the greater the driving force, the greater the impact on the surrounding flow field. Therefore, it is necessary to consider constraining the maximum values of multiple sets of driving forces and based on this, calculate the optimal control of the underwater robot. According to the driving performance of the driver, the control constraints for the resultant force and resultant torque applied to the robot are set to 20 N and 15 Nm, respectively.

Using Simulink software to build a simulation model of the underwater robots, corresponding robot dynamics models and control system models were established based on the above model parameters and setting conditions, and comparative simulation tests and analysis of dynamic optimal control were conducted. In the simulation, the velocity of water flow in the x-axis, y-axis, and z-axis directions is set to 0.1 m/s. The industrial robot starts from the coordinate origin and ends at (20, 10, 5). The CVP method utilizes the approximate control effect of piecewise functions to transform the optimal control problem into a class of NLP problems for a solution. The CVP method is widely used and has the advantage of a relatively simple solution, but it has the disadvantage of inaccurate optimization time nodes and control parameters due to uniform time grid partitioning. In order to intuitively reflect the effectiveness of the proposed control method, optimization constraint problem solving methods based on control vector parameterization (CVP), the Gaussian pseudo-spectral method, and the adaptive Gaussian pseudo-spectral method were carried out, and simulation results were compared and analyzed.

The driving trajectories of industrial robots obtained under different control methods are shown in Figure 3. As shown in the figure, the overall trend of the trajectory curves of industrial robots under the three algorithms is relatively consistent. Among them, the trajectory curve under the CVP algorithm has relatively larger curvature changes and driving mileage, while the trajectory curve under the GPM algorithm has significantly improved in curvature changes and driving mileage compared to the CVP algorithm. In addition, the relevant indicators of the trajectory curve under the AGPM algorithm have been further improved compared to the GPM algorithm. In order to further demonstrate the comparison of the driving trajectories of industrial robots, the projected trajectories of industrial robot trajectories on different coordinate planes were compared, as shown in Figure 4. Combining the projected trajectories on the OXY plane, the OXZ plane, and the OYZ plane, it can be seen that the trajectory projection curve under the AGPM algorithm tends to be closer to the diagonal line from the starting point to the ending point, followed by the trajectory projection curve under the GPM algorithm, while the curvature change of the trajectory projection curve under the CVP algorithm is the most significant, and its overall driving mileage is relatively longer. Based on the comparison of trajectory projection curves, it can be concluded that the APGM algorithm proposed in this study can achieve optimal control of industrial robot dynamics. The optimization results show that the overall mileage of the driving trajectory is the smallest, and the curvature change is smoother, thus effectively achieving energy and efficiency optimization of overall control.

The control results of the industrial robot’s driving posture angle under different control methods are shown in Figure 5. In the comparison of yaw angles, the magnitude of the yaw angle gradually increases during the initial start-up of the robot, reaches its maximum value in the middle of the driving trajectory, and then gradually decreases. Under the three control methods, the overall yaw angle obtained by the AGPM algorithm is the smallest, while the yaw angle obtained by the CVP algorithm is relatively the largest. The change in the size of the yaw angle is also consistent with the trajectory trend of the robot. According to the comparison of the roll angles of industrial robots, it can be seen that the roll angle obtained by the AGPM algorithm is also smaller compared to the other two control methods. In the comparison of pitch angle control results, there is a significant fluctuation in the pitch angle under the CVP algorithm during the initial start-up stage. In contrast, the pitch angle under the GPM algorithm and the AGPM algorithm is generally stable in terms of variation amplitude, with the pitch angle obtained by the AGPM algorithm being the smallest in comparison. Based on the trajectory comparison results in Figure 2 and Figure 3, as well as the attitude angle comparison results in Figure 4, it can be concluded that the AGPM algorithm can achieve trajectory optimization control while reducing the yaw angle during the industrial robot’s movement, ensuring optimal driving efficiency from the starting point to the ending point.

The control results of the industrial robot’s driving attitude and angular velocity under different control methods are shown in Figure 6. In the control results of lateral angular velocity, the AGPM algorithm has the smallest overall lateral angular velocity and the fastest convergence speed, proving that the application of the AGPM algorithm effectively improves the lateral stability of industrial robots. According to the comparison results of roll angle velocity, it can be seen that there are severe fluctuations of different amplitudes and frequencies in the roll angle velocity under both the CVP algorithm and the GPM algorithm. Among them, the roll angle velocity fluctuation under the CVP algorithm is the most obvious, with the highest amplitude and frequency of the fluctuation. In contrast, the overall change in roll angle velocity under the AGPM algorithm is the smoothest, with almost no rapid amplitude fluctuations or shaking. The comparison results of pitch angular velocity obtained under the three control methods also showed similar trends of change. The pitch angular velocity under the CVP algorithm exhibits significant and severe fluctuations in amplitude during the initial start-up stage. Although the fluctuation amplitude decreases in the subsequent process, it remains the largest among the three control algorithms. In contrast, the pitch angular velocity under the GPM and the AGPM algorithms is significantly smaller overall, with the AGPM algorithm having the smallest pitch angular velocity fluctuation and the smoothest variation curve. The above comparison results reflect that the AGPM algorithm can significantly improve the attitude angular velocity control results of industrial robots, making them more stable overall and effectively ensuring the driving stability of industrial robots. The resultant control force results of an underwater industrial robot under different control methods are shown in Figure 7. As shown in the figure, the trend of the comparison results with other state variables is consistent. The control force under the CVP algorithm is correspondingly the largest, while the control force of the underwater robot under the AGPM algorithm is significantly the smallest, which helps to achieve the energy-saving goal of the underwater industrial robot.

In order to further demonstrate the optimal control performance of the underwater industrial robot dynamics under different control algorithms, a comparative analysis of optimization solution capability performance indicators was conducted, and the energy consumption performance of an underwater industrial robot under different algorithms was also compared, as shown in Table 2. All three control algorithms use the interior point method for solving optimization constraint problems. The CVP algorithm uses the unified grid method to discretize the state and control variables, while the GPM algorithm and the AGPM algorithm use the LG collocation method and adaptive LG collocation method for discretization, respectively. The three control algorithms are optimized and solved under different allocation points of control variables. Among them, the number of allocation points under the CVP algorithm is 30, the number of allocation points under the GPM algorithm is 38, and the number of allocation points under the AGPM algorithm is 17. Meanwhile, based on the optimization solution times under the three control algorithms, it can be concluded that the optimization solution times for CVP, GPM, and AGPM are 19.83 s, 11.65 s, and 5.33 s, respectively. By comparison, it can be seen that the CVP algorithm has the longest optimization solution time. Compared with the CVP algorithm, the GPM algorithm can effectively improve optimization solving time, but its number of allocation points is relatively large. In contrast, the AGPM algorithm achieves a significant reduction in optimization solution time with fewer allocation points. The comparative results show that after adopting the two-stage adaptive Gauss collocation algorithm, the number of control variable collocation points are significantly reduced, thereby effectively reducing its optimization solution time. In addition, the energy consumption in maneuvering cruise of underwater industrial robot under the CVP, the GPM, and the AGPM algorithms is 206.79 J, 186.34 J, and 173.61 J, respectively, indicating that the proposed optimization control algorithm can achieve the control goal of energy-saving for underwater industrial robots. Based on the above data, it can be concluded that compared to the CVP and the GPM algorithms, the AGPM algorithm has improved energy consumption performance by 16.04% and 6.83%, respectively.

On the basis of verifying the efficiency performance indicators of different optimization allocation algorithms and the energy consumption of underwater robots, further statistical analysis of data is used to reflect the performance of underwater robot optimization control in terms of dynamic performance. Trajectory tracking error and attitude control error are selected as indicators to characterize the dynamic performance of underwater robots. At the same time, in the process of calculating dynamic performance errors, the error suppression effect is comprehensively reflected by obtaining the average error value and root mean square error of the corresponding variables. The calculation formula [36] can be expressed as

(26)$\begin{array}{c}{e}_{AVE}={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N\left({\widehat{x}}_i-x_i\right)}\\ e_{RMSE}=\sqrt{{\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N{\left({\widehat{x}}_i-x_i\right)}^2}}\end{array},$

where e_AVE is the error average value, e_RMSE is the error root mean square, N is the sampling amount, $x_i$ is the referenced value, and ${\widehat{x}}_i$ is the estimated value.

The e_AVE and e_RMSE results of the calculated trajectory tracking error and attitude control error are shown in Table 3 and Table 4, respectively. According to Table 3, in terms of trajectory tracking error, under the CVP algorithm, the average tracking error obtained in the projected trajectories of OXY, OXZ, and OYZ planes is relatively larger. On the contrary, the average trajectory tracking error under the AGPM algorithm is relatively small. Similarly, in terms of attitude control error of underwater robots, the average error of attitude tracking control obtained by using the AGPM algorithm is significantly smaller. The results in Table 3 further demonstrate through quantitative data comparison that the proposed AGPM algorithm can effectively improve the trajectory tracking and attitude control performance of underwater robots. Then, by comparing the root mean square values of trajectory error and attitude control error in Table 4, it can be further found that under the AGPM algorithm, the e_RMSE of trajectory tracking control and attitude control errors are relatively smaller. This indicates that the proposed AGPM algorithm can achieve smoother and more stable motion control effects, reducing the attitude fluctuations of underwater robots during maneuvering, which also helps to achieve comprehensive optimization of navigation stability and energy-saving effects. Based on the above data, it can be concluded that compared to the CVP and the GPM algorithms, the trajectory tracking accuracy of the AGPM algorithm has improved by 15.76% and 8.13%, respectively, and the attitude control accuracy of the AGPM algorithm has improved by 17.29% and 8.62%, respectively.

5. Conclusions

In order to achieve efficient solution goals for the dynamic optimization control and optimization constraint problems of underwater robots, this paper proposes a dynamic optimization control method for underwater robots based on an improved Gaussian pseudo-spectral method, which solves the contradiction between the increase in discrete points and the decrease in computational efficiency in the GPM method. In the dynamic modeling of industrial robots, the motion coordinate system of industrial robots and the corresponding kinematic coordinate transformation relationship in the coordinate system were established. Based on the force analysis of the robot, a six degree of freedom dynamic model of industrial robots was summarized. We designed performance indicators that balance energy and efficiency for the dynamic control of industrial robots. Then, within the Gauss pseudo-spectral control framework, we provided the Legendre–Gauss collocation solution formula and proposed a two-stage adaptive Gauss collocation reconstruction strategy for control variables. The comparative simulation results show that the proposed AGPM algorithm can effectively improve the dynamic optimization control effect of industrial robots. Compared with the CVP algorithm and the GPM algorithm, the driving range, driving attitude control stability, and optimization solution time under the AGPM algorithm can be optimized. Compared to the CVP and the GPM, the AGPM algorithm has improved energy consumption performance by 16.04% and 6.83%, trajectory tracking accuracy by 15.76% and 8.13%, and attitude control accuracy by 17.29% and 8.62%, respectively.

In future research, it is expected to achieve new technological breakthroughs in model accuracy and environmental adaptability, as well as robustness and intelligence of control algorithms. The underwater robot dynamics model established in this study, although considering hydrodynamic effects, may still have limitations in simulating complex marine environments. Future research can further refine the model, improve its adaptability to extreme water quality conditions, and explore the integration of new sensor technologies and data processing algorithms, such as using machine learning and deep learning techniques to enhance perception capabilities. With the development of artificial intelligence technology, future research can explore more intelligent control strategies based on deep learning and reinforcement learning to improve the autonomous learning and adaptability of underwater robot control systems. In addition, in practical applications, research on the safety and robustness of robot control systems is also crucial.

Author Contributions

Conceptualization, J.Z., X.Z., and T.C.; methodology, J.Z. and T.C.; software, J.Z. and X.Z.; validation, J.Z., X.Z., and G.D.; writing—original draft preparation, J.Z. and T.C.; writing—review and editing, J.Z. and G.D.; visualization, J.Z. and X.Z. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Nomenclature

GPM	Gaussian pseudo-spectral method
NLP	Nonlinear programming
CVP	Control vector parameterization
LG	Legendre–Gauss
SQP	Sequential quadratic programming
AGPM	Adaptive Gaussian pseudo-spectral method
φ	Roll angle of the robot
θ	Pitch angle of the robot
ϕ	Yaw angle of the robot
$u^{\prime }$	Velocity at the origin in the local coordinate system
u	Velocity along the x-axis in the local coordinate system
v	Velocity along the y-axis in the local coordinate system
w	Velocity along the z-axis in the local coordinate system
$\mathrm{\Omega }$	Angular velocity at the origin in the local coordinate system
p	Angular velocity along the x-axis in the local coordinate system
q	Angular velocity along the y-axis in the local coordinate system
r	Angular velocity along the z-axis in the local coordinate system
Tv	Velocity transformation matrix
$T_{\omega }$	Angular velocity transformation matrix
F	The resultant force exerted on the robot
$\tau$	The resultant moment exerted on the robot
B	Buoyancy
G	Gravity of the industrial robot
$u_f^{\prime }$	Velocity of water flow in the local coordinate system
$u_r^{\prime }$	Velocity of water flow relative to the robot
ax/ay	Longitudinal/lateral accelerations of vehicle

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Dynamic model of industrial robot.

Figure 2. Adaptive Gauss collocation reconstruction process.

Figure 3. Trajectory tracking results of industrial robot.

Figure 4. The projection of the trajectory of industrial robot on different planes. (a) OXY plane, (b) OXZ plane, and (c) OYZ plane.

Figure 5. The control results of the driving posture angle of industrial robot. (a) Yaw angle, (b) roll angle, and (c) pitch angle.

Figure 6. The control results of the driving attitude angular rate of industrial robot. (a) Yaw rate, (b) roll rate, and (c) pitch rate.

Figure 7. The resultant control force of industrial robot.

Table 1

Pseudocode of the adaptive Gaussian pseudo-spectral algorithm.

Adaptive Gaussian Pseudo-Spectral Algorithm

(1) Step 1: Import the mathematical models for robot dynamics optimization control.(2) Step 2: Set algorithm optimization parameters.(3) Step 3: The state variable LG allocation and control variable pair LG allocation in the first stage are used to obtain a finite dimensional NLP optimization solution problem.(4) Step 4: Solve the NLP problems. If the solution is successful, proceed to step 5. Otherwise, determine whether LG point reselection is needed. If so, proceed to step 3. Otherwise, the solution ends.(5) Step 5: In the second stage of point allocation reconstruction, the slope change is calculated based on the solution results of the first stage, and the control variable LG is reconstructed to obtain a new finite dimensional NLP problem. Then, the step 6 is entered.(6) Step 6: Solve the finite dimensional NLP problems. If the solution is successful, end the process. Otherwise, determine whether LG point reselection is needed. If it is necessary, proceed to step 3. Otherwise, end the process.

Table 2

Comparison of performance indicators of different control algorithms.

Method	Collocation Points	NLP Solution Method	Discrete Method	Solution Time	Energy Consumption
CVP	30	Interior point method	Unified grid	19.83 s	206.79 J
GPM	38	Interior point method	LG collocation	11.65 s	186.34 J
AGPM	17	Interior point method	Adaptive LG collocation	5.33 s	173.61 J

Table 3

e_AVE of performance indicators under different control algorithms.

Method	Trajectory Tracking Error			Attitude Control Error
	OXY (m)	OXZ (m)	OYZ (m)	Yaw Rate (deg/s)	Roll Rate (deg/s)	Pitch Rate (deg/s)
CVP	0.1912	0.0839	0.2177	0.3992	0.1217	0.0691
GPM	0.0865	0.0498	0.0587	0.2385	0.0932	0.0567
AGPM	0.0462	0.0305	0.0336	0.1679	0.0768	0.0409

Table 4

e_RMSE of performance indicators under different control algorithms.

Method	Trajectory Tracking Error			Attitude Control Error
	OXY (m)	OXZ (m)	OYZ (m)	Yaw Rate (deg/s)	Roll Rate (deg/s)	Pitch Rate (deg/s)
CVP	0.0286	0.0187	0.0388	0.1996	0.2901	0.3967
GPM	0.0237	0.0160	0.0216	0.1533	0.1649	0.2941
AGPM	0.0195	0.0135	0.0177	0.1171	0.1208	0.1536

References

1. Phung, M.D.; Ha, Q.P. Safety-enhanced UAV Path Planning with Spherical Vector-Based Particle Swarm Optimization. Appl. Soft Comput.; 2021; 107, 107376. [DOI: https://dx.doi.org/10.1016/j.asoc.2021.107376]

2. Li, Y.; He, J.; Shen, H.; Zhang, W.; Li, Y. Adaptive practical prescribed-time fault-tolerant control for autonomous underwater vehicles trajectory tracking. Ocean Eng.; 2023; 277, 114263. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2023.114263]

3. Li, J.; Du, J.; Chang, W.J. Robust time-varying formation control for underactuated autonomous underwater vehicles with disturbances under input saturation. Ocean Eng.; 2019; 179, pp. 180-188. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2019.03.017]

4. Liu, H.; Lyu, Y.; Lewis, F.L.; Wan, Y. Robust time-varying formation control for multiple underwater vehicles subject to nonlinearities and uncertainties. Int. J. Robust Nonlinear Control; 2019; 29, pp. 2712-2724. [DOI: https://dx.doi.org/10.1002/rnc.4517]

5. Liu, H.; Yan, S.C.; Shen, Y.; Li, C.H.; Zhang, Y.F.; Hussain, F. Model predictive control system based on direct yaw moment control for 4WID self-steering agriculture vehicle. Int. J. Agric. Biol. Eng.; 2021; 14, pp. 175-181. [DOI: https://dx.doi.org/10.25165/j.ijabe.20211402.5283]

6. Jin, Y.C.; Liu, J.Z.; Xu, Z.J.; Yuan, S.Q.; Li, P.P.; Wang, J.Z. Development status and trend of agricultural robot technology. Int. J. Agric. Biol. Eng.; 2021; 14, pp. 1-19. [DOI: https://dx.doi.org/10.25165/j.ijabe.20211404.6821]

7. Saravanakumar, T.; Lee, T. Hybrid-driven-based resilient control for networked T-S fuzzy systems with time-delay and cyber-attacks. Int. J. Robust Nonlinear Control; 2023; 33, pp. 7869-7891. [DOI: https://dx.doi.org/10.1002/rnc.6801]

8. Ghrab, S.; Ahmed Ali, S.; Benamor, A.; Langlois, N.; Messaoud, H. A novel robust discrete-time integral sliding mode tracking control design for time-varying delay MIMO systems with unknown uncertainties. ISA Trans.; 2024; 145, pp. 1-18. [DOI: https://dx.doi.org/10.1016/j.isatra.2023.11.024] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/38016883]

9. Liu, Y.; Cui, D. Optimal Control of Vehicle Path Tracking Problem. World Electr. Veh. J.; 2024; 15, 429. [DOI: https://dx.doi.org/10.3390/wevj15090429]

10. Cheng, S.; Li, L.; Chen, X.; Wu, J. Model-predictive-control-based path tracking controller of autonomous vehicle considering parametric uncertainties and velocity-varying. IEEE Trans. Ind. Electron.; 2020; 68, pp. 8698-8707. [DOI: https://dx.doi.org/10.1109/TIE.2020.3009585]

11. Vimal Kumar, S.; Kim, J. Robust H_∞ Control for Autonomous Underwater Vehicle’s Time-Varying Delay Systems under Unknown Random Parameter Uncertainties and Cyber-Attacks. Appl. Sci.; 2024; 14, 8827. [DOI: https://dx.doi.org/10.3390/app14198827]

12. Campos-Martínez, S.-N.; Hernández-González, O.; Guerrero-Sánchez, M.-E.; Valencia-Palomo, G.; Targui, B.; López-Estrada, F.-R. Consensus Tracking Control of Multiple Unmanned Aerial Vehicles Subject to Distinct Unknown Delays. Machines; 2024; 12, 337. [DOI: https://dx.doi.org/10.3390/machines12050337]

13. Zhang, X.; Xu, X.; Li, J.; Luo, Y.; Wang, G.; Brunauer, G.; Dustdar, S. Observer-based H_∞ fuzzy fault-tolerant switching control for ship course tracking with steering machine fault detection. ISA Trans.; 2023; 140, pp. 32-45. [DOI: https://dx.doi.org/10.1016/j.isatra.2023.05.021] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37295998]

14. Tian, Y.; Sun, J.; Zhou, X.; Wu, X.; Lu, B.; Dai, C. Research on apple origin classification based on variable iterative space shrinkage approach with stepwise regression-support vector machine algorithm and visible-near infrared hyperspectral imaging. J. Food Process Eng.; 2020; 43, e13432. [DOI: https://dx.doi.org/10.1111/jfpe.13432]

15. Pan, W.; Wang, Y.; Song, F.; Peng, L.; Zhang, X. UUV-Assisted Icebreaking Application in Polar Environments Using GA-SPSO. J. Mar. Sci. Eng.; 2024; 12, 1845. [DOI: https://dx.doi.org/10.3390/jmse12101845]

16. Liu, B.; Ouyang, W.; Yan, X.; Vanierschot, M. Numerical Investigation of Hydrodynamic Characteristics of a Rim-Driven Thruster Coupled with an Underwater Vehicle. J. Mar. Sci. Eng.; 2024; 12, 1838. [DOI: https://dx.doi.org/10.3390/jmse12101838]

17. Feng, J.; Xu, W.; Dong, J.; Yao, Y.; Hu, Z. A UUV Cluster Route-Planning Method for Dynamic Target Search. Electronics; 2024; 13, 4033. [DOI: https://dx.doi.org/10.3390/electronics13204033]

18. Xu, Z.; Shen, Y.; Xie, Z.; Liu, Y. Research on Autonomous Underwater Vehicle Path Optimization Using a Field Theory-Guided A* Algorithm. J. Mar. Sci. Eng.; 2024; 12, 1815. [DOI: https://dx.doi.org/10.3390/jmse12101815]

19. Zheng, X.; Wang, Y.; Tian, Q.; Zhang, Q.; Wang, X.; Xu, W.; Wang, G.; Yang, X.; Sun, Y. An Adaptive Controller with Disturbance Observer for Underwater Vehicle Manipulator Systems. Electronics; 2024; 13, 3938. [DOI: https://dx.doi.org/10.3390/electronics13193938]

20. Rex, P.T.; Abbott, K.J.; Prezgay, R.E.; Lowe, C.G. The Effects of Depth and Altitude on Image-Based Shark Size Measurements Using UAV Surveillance. Drones; 2024; 8, 547. [DOI: https://dx.doi.org/10.3390/drones8100547]

21. Chen, T.; Cai, Y.; Chen, L.; Xu, X. Trajectory and Velocity Planning Method of Emergency Rescue Vehicle Based on Segmented Three-Dimensional Quartic Bezier Curve. IEEE Trans. Intell. Transp. Syst.; 2023; 24, pp. 3461-3475. [DOI: https://dx.doi.org/10.1109/TITS.2022.3224785]

22. Chen, T.; Xu, X.; Cai, Y.; Chen, L.; Li, K. QPSOMPC-based chassis coordination control of 6WIDAGV for vehicle stability and trajectory tracking. J. Frankl. Inst.; 2024; 362, 107458. [DOI: https://dx.doi.org/10.1016/j.jfranklin.2024.107458]

23. Lu, E.; Ma, Z.; Li, Y.M.; Xu, L.Z.; Tang, Z. Adaptive backstepping control of tracked robot running trajectory based on real-time slip parameter estimation. Int. J. Agric. Biol. Eng.; 2020; 13, pp. 178-187. [DOI: https://dx.doi.org/10.25165/j.ijabe.20201304.5739]

24. Xu, G.; Fang, H.; Song, Y.; Du, W. Optimal Design and Analysis of Cavitating Law for Well-Cellar Cavitating Mechanism Based on MBD-DEM Bidirectional Coupling Model. Agriculture; 2023; 13, 142. [DOI: https://dx.doi.org/10.3390/agriculture13010142]

25. Liu, Y.J.; Cui, D.W.; Peng, W. Optimum Control for Path Tracking Problem of Vehicle Handling Inverse Dynamics. Sensors; 2023; 23, 6673. [DOI: https://dx.doi.org/10.3390/s23156673] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37571456]

26. Wang, B.; Zhang, Y.M.; Zhang, W. Integrated path planning and trajectory tracking control for quadrotor UAVs with obstacle avoidance in the presence of environmental and systematic uncertainties: Theory and experiment. Aerosp. Sci. Technol.; 2022; 120, pp. 107277-107294. [DOI: https://dx.doi.org/10.1016/j.ast.2021.107277]

27. Zhu, C.; Chen, J.; She, J.; Zhang, H. Actuator Fault Reconstruction for Quadrotors Using Deep Learning-Based Proportional Multiple-Integral Observer. IEEE Trans. Ind. Electron.; 2023; 71, pp. 7986-7995. [DOI: https://dx.doi.org/10.1109/TIE.2023.3294631]

28. Zuo, Y.; Chen, J.; Zhu, X.; Lee, C.H.T. Different Active Disturbance Rejection Controllers Based on the Same Order GPI Observer. IEEE Trans. Ind. Electron.; 2022; 69, pp. 10969-10983. [DOI: https://dx.doi.org/10.1109/TIE.2021.3118378]

29. Duan, G.; Yao, L.; Zhan, Z.; Kang, T.; Guo, C. Multi-Objective Optimal Control Method for the 6-DOF Robotic Crusher. Appl. Sci.; 2024; 14, 9397. [DOI: https://dx.doi.org/10.3390/app14209397]

30. Thuyen, N.A.; Thanh, P.N.N.; Anh, H.P.H. Adaptive finite-time leader-follower formation control for multiple AUVs regarding uncertain dynamics and disturbances. Ocean Eng.; 2023; 269, 113503. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2022.113503]

31. Zeng, X.; Yang, L.; Zhu, Y.; Yang, F. Comparison of Two Optimal Guidance Methods for the Long-Distance Orbital Pursuit-Evasion Game. IEEE Trans. Aerosp. Electron. Syst.; 2021; 57, pp. 521-539. [DOI: https://dx.doi.org/10.1109/TAES.2020.3024423]

32. Pérez-Alvarado, M.E.; Gómez-Espinosa, A.; González-García, J.; García-Valdovinos, L.G.; Salgado-Jiménez, T. Convolutional Long Short-Term Memory Predictor for Collaborative Remotely Operated Vehicle Trajectory Tracking in a Leader–Follower Formation Subject to Communication and Sensor Latency in the Presence of External Disturbances. Machines; 2024; 12, 691. [DOI: https://dx.doi.org/10.3390/machines12100691]

33. Li, Z.; Wu, X.; Zhang, W.; Yu, W. Obstacle Avoidance Trajectory Planning for Autonomous Vehicles on Urban Roads Based on Gaussian Pseudo-Spectral Method. World Electr. Veh. J.; 2024; 15, 7. [DOI: https://dx.doi.org/10.3390/wevj15010007]

34. Shojaei, K.; Chatraei, A. Robust platoon control of underactuated autonomous underwater vehicles subjected to nonlinearities, uncertainties and range and angle constraints. Appl. Ocean Res.; 2021; 110, 102594. [DOI: https://dx.doi.org/10.1016/j.apor.2021.102594]

35. Prestero, T. Verification of a Six Degree of Freedom Simulation Model for the REMUS Autonomous Underwater Vehicle. Ph.D. Thesis; Massachusetts Institute of Technology/Woods Hole Oceanographic Institution: Cambridge, MA, USA, 2001.

36. Chen, T.; Cai, Y.F.; Chen, L.; Xu, X. Sideslip angle fusion estimation method of three-axis autonomous vehicle based on composite model and adaptive cubature Kalman filter. IEEE Trans. Transp. Electrif.; 2024; 10, pp. 316-330. [DOI: https://dx.doi.org/10.1109/TTE.2023.3263592]