Decomposed optimal trajectory control of hydraulic actuated quadruped robot based on model prediction

Abstract

Hydraulic actuated quadruped robot has strong load capacity, good flexibility, and strong adaptability to complex environments. Pointing at the multi coupling and multi variable characteristics of hydraulic quadruped robot, this paper proposes a decomposed trajectory tracking control strategy based on model prediction for hydraulic quadruped robot. The tracking model is derived by the structure of the hydraulic quadruped robot. The inequality constraints of the quadruped robot motion are built. And the trajectory tracking control problem of the robot is transformed into a quadratic programing problem. The error model of the robot is built. Combining the minimization goals of robot motion speed and acceleration error, the stability and continuity constraints for robot motion are planned. The decomposed trajectory tracking controller with trajectory constraints is designed. The performance of proposed control method is verified by simulation and prototype experiment. The simulation and prototype comparative experiments of classical control strategy and the proposed control strategy are conducted on the hydraulic robot. The results show that compared with the traditional method, the motion trajectory error of the hydraulic robot is relatively small by the proposed optimal trajectory tracking control strategy. The rolling angle accuracy is ±0.066 rad, and the pitching angle accuracy is −0.073 to 0.044 rad. The effectiveness of the proposed control algorithm is verified by the simulation and experimental results. The accuracy of the hydraulic robot model in the control algorithm is improved, the impact during hydraulic robot motion is reduced, and the service life of the robot was increased, as well as the stability of robot motion is enhanced.

Introduction

With the development of artificial intelligence, human-computer interaction and other technologies, robots play an important role in aerospace, earthquake rescue, and other scenarios.^1–3 Hydraulic actuated quadruped robots have strong load capacity, good flexibility, and strong adaptability to complex environments, with broad application prospects.^4,5 The motion control of quadruped robot is the basic problem of robot system control. The trajectory tracking control problem is the focus of motion control research. And it is also a key technology for achieving intelligent control of robots. The trajectory tracking control of the robot can ensure that the actual trajectory track the desired trajectory, and improve the dynamic stability of the robot system.

It is unable to obtain the unknown environmental information for a quadruped robot in an unknown and complex environment. And it is difficult to describe the external interference of robot by an accurate mathematical model. In addition, due to the characteristics of underdrive, multivariable, strong coupling and nonlinearity, the trajectory tracking control performance of quadruped robot needs to be improved.

The trajectory tracking control methods of quadruped robot mainly include sliding mode control method,^6,7 fuzzy control method,^8–10 neural network control method,^11–13 and so on. The above control methods have high computational efficiency. However, the chattering phenomenon will appear in the high frequency switching by using the sliding mode control method. The uncertainty of system parameters can be overcome by the fuzzy control method, but the control accuracy of the system will be affected by the chattering of the controller input. Neural network control method needs online learning, and its real-time performance is poor.

Furthermore, Bellicoso et al. ¹⁴ proposed an online real-time motion planning method based on zero moment point of quadruped robot. According to the contact time and the real-time state of the robot, the reference trajectory was constantly updated. The reference foothold of each leg of the robot was obtained by solving the optimization problem. A variety of gait movements and smooth switching between each gait were realized on the ANYmal quadruped robot platform.

The above control method only considers the output tracking of the robot system, and the impact of control increments on the system is ignored. The kinematic and dynamic constraints of robots are often ignored or simplified in traditional control algorithms.^1,15,16 Such as constraints of motion position, constraints of motion speed, constraints of the joint speed, constraints of non-sliding of the foot, and so on. Nevertheless, these kinematic constraints typically have a significant impact on the stable motion of the robot.

The problems with multiple variables and constraints could be solved using model predictive control.^17–19 By optimizing the objective function, the model predictive control algorithm could combine the kinematic constraints of the quadruped robot, and real-time online calculation can be carried out. The future output can be predicted according to the historical value and the current value of the system by model predictive control. And the model predictive control method has good anti-interference ability and tracking performance.

Carron et al. ²⁰ proposed a trajectory tracking control method based on robot kinematics model, dynamics model and model predictive control. The linearized inverse dynamics feedback and error model of the manipulator were established. The residual interference was estimated online through the additional interference model, and the non-offset tracking of the manipulator was realized.

Grandia et al. ²¹ proposed a model predictive control method based on sequential linear quadratic programing for the optimal control of mobile robots. The robust solution was obtained according to the frequency cost function, which significantly improved the tracking performance of the kinematic trajectory, torque trajectory, and force trajectory of the ANYmal quadruped robot. Based on this, Grandia R et al. ²² proposed a model predictive feedback control method based on differential dynamic programing, which compensated the gap between the update rate of low model predictive control and the driving instruction. The sensitivity of the system to the high frequency model errors and actuator bandwidth constraints was reduced by the frequency dependent cost function. And the stable motion control of quadruped robot was realized by this method.

Minniti et al. ²³ used the model predictive control method to optimize the balance and interaction of the unstable manipulator. The optimal control problem was transformed into the end effector space. The position and direction tasks were processed in the model predictive control planner, and the contact force of the end effector was planned. The effectiveness of the attitude tracking controller was proved by experiments. Neunert et al. ²⁴ proposed a model predictive control method for the whole nonlinear contact rigid body system. A full dynamic system model including explicit contact dynamics was established. The contact position and order of quadruped robot system were calculated in real time by the model predictive controller. The validity and robustness of the proposed method were verified by the experiment of quadruped robot prototype.

As can be seen from the above research, the model predictive control method is an effective method to solve the trajectory tracking control problem of mobile robot. And the robustness of the robot control system can be improved.^25,26 The combination of model predictive control method and other methods was used by some scholars to complete the control of robots.^19,27,28 Villarreal et al. ²⁹ combined convolutional neural networks with model predictive control methods to determine the locations of the foothold of the hydraulic quadruped robot HyQReal. Qin et al. ³⁰ proposed a control method based on central pattern generator and model predictive control for quadruped robot, and the position tracking was achieved.

Furthermore, the MPC method can be used for short-term prediction and rolling optimization, and the model is re-linearized in each iteration to reduce the cumulative error caused by contact mutations. Through Lyapunov functions or constraint contraction theory, some studies have demonstrated that linearized MPC closed-loop systems can ensure the stability of the system under finite disturbances or the known contact sequence.^31–34 And the robot system could demonstrate good performance. Such as, the robot dynamics were simplified to formulate the problem as convex optimization by Carlo et al. ³⁵ Linearization error was alleviated by real-time state information feedback of the Cheetah3 robot.

However, there are currently few methods for trajectory control of hydraulic quadruped robots. The hydraulic quadruped robot system is complex, and the trajectory tracking accuracy of the foot end is affected by the position tracking accuracy of joint actuators.

In order to improve the accuracy and stability of the trajectory of quadruped robot, and restrain the influence of disturbance on the trajectory of robot, this paper proposed a decomposed trajectory tracking control strategy based on model predictive control. The tracking model and error model of quadruped robot are established. According to the inequality constraints in the motion process of quadruped robot, the trajectory tracking control problem is transformed into a quadratic programing problem. For the problem of disturbance rejection trajectory tracking control of quadruped robot, the trajectory tracking problem of the whole robot is divided into the trajectory tracking control problem of each foot end and the trajectory tracking control problem of the torso. The trajectory tracking target and constraint conditions of quadruped robot are established. The effectiveness of the algorithm proposed in this paper is verified through simulation and experiments.

Structural analysis and tracking model of hydraulic robot

Structure composition

The mechanical structure of a hydraulic quadruped robot is designed based on the structure of a quadruped mammal, consisting of a body and four bionic legs, as shown in Figure 1. The single leg has three active degrees of freedom, namely, yaw joint, hip joint roll, and knee joint. The motion pair between the thigh and the body is the yaw joint, the motion pair between the thigh and shank is the hip joint, and the motion pair between the ankle and shank is the knee joint. Linear bearings and shock-absorbing springs are installed between the knee joint and the foot end. And the shock-absorbing spring fixed on the shank leg of the robot is a passive degree of freedom. Each active joint is driven by the same electro-hydraulic servo actuator, which integrates a single rod hydraulic cylinder, electro-hydraulic servo valve, displacement sensor, and pressure sensor.

Figure 1.

Simulation prototype of the hydraulic quadruped robot.

Tracking model

The position vector of quadruped robot in the inertial coordinate system is p  = {x, y, z}^T, the attitude vector of quadruped robot in inertial coordinate system is q _t  = {q₁, q₂, q₃}^T. Hence, the position and attitude vector of quadruped robot in inertial coordinate system is q  = {x, y, z, q₁, q₂, q₃}^T. The linear velocity vector of quadruped robot in the torso coordinate system is v _b = {v _x , v _y , v _z }^T, the angular velocity vector is ω _b = {ω_bx, ω_by, ω_bz}^T. Hence, the velocity vector of the quadruped robot in the torso coordinate system is u  = {v _x , v _y , v _z , ω_bx, ω_by, ω_bz}^T. It can be concluded that the kinematic transformation equation of quadruped robot is

$$\stackrel{·}{\mathit{q}}=\mathit{J}(q)\mathit{u}$$

(1)

In which, J (q) is the Jacobian transformation matrix between the inertial coordinate system and the torso coordinate system.

In this paper, the trajectory tracking problem of robot is studied. The forward motion, lateral motion and yaw motion of robot are considered. The position and attitude vector of quadruped robot in inertial coordinate system is q  = {x, y, q₃}^T. The velocity vector of quadruped robot in torso coordinate system is u  = {v _x , v _y , ω_bz}^T. According to formula (1), we can get

$$\stackrel{·}{q}=\left[\begin{array}{c}\stackrel{·}{x}\hfill \\ \multicolumn{1}{c}{\stackrel{·}{y}}\\ \multicolumn{1}{c}{\stackrel{·}{q}}\end{array}\right]=J(q)\left[\begin{array}{c}{v}_x\hfill \\ \multicolumn{1}{c}{v_y}\\ \multicolumn{1}{c}{{\omega }_{b_2}}\end{array}\right]=\left[\begin{array}{ccc}\cos q_3\hfill & -\sin q_3\hfill & 0\hfill \\ \multicolumn{1}{c}{\sin q_3}& \cos q_3\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right]\left[\begin{array}{c}{v}_x\hfill \\ \multicolumn{1}{c}{v_y}\\ \multicolumn{1}{c}{{\omega }_{b_2}}\end{array}\right]$$

(2)

The Jacobian matrix of a hydraulic quadruped robot is

$$\mathit{J}=\left[\begin{array}{ccc}0\hfill & l_2{c}_2+l_3{c}_{23}\hfill & l_3{c}_{23}\hfill \\ \multicolumn{1}{c}{c_1(l_1+l_2{c}_2+l_3{c}_{23})}& -s_1(l_2{s}_2+l_3{s}_{23})\hfill & -l_3{s}_1{s}_{23}\hfill \\ \multicolumn{1}{c}{-s_1(l_1+l_2{c}_2+l_3{c}_{23})}& c_1(l_2{s}_2+l_3{s}_{23})\hfill & l_3{c}_1{s}_{23}\hfill \end{array}\right]$$

(3)

In which, l₁,l₂, and l ₃ are the lengths of the hip link, thigh link and shank link, respectively. And θ_1i, θ_2i, and θ_3i are joint angles of roll hip joint, pitch hip joint, and pitch knee joint, respectively. s₁ = sinθ₁, c₁ = cosθ₁, s₂₃ = sin(θ₂+θ₃), c₂₃ = cos(θ₂+θ₃), and the other angles are similar.

The dynamic equation of the quadruped mobile robot is as follows

$$\mathit{M}(\mathit{R},\mathit{q})\mathit{q}+\mathit{N}(\mathit{R},\mathbf{\omega },\mathit{q},\stackrel{·}{\mathit{q}})=\mathbf{\tau }+{\mathit{J}}_c\mathbf{\lambda }$$

(4)

Where, M is the mass matrix of the quadruped robot. R is the rotation transformation matrix of the body coordinate system relative to the inertial coordinate system. q is the position and attitude vectors of the robot body relative to the inertial coordinate system are respectively. N is a non-contact force vector matrix, including Coriolis force, centrifugal force, and gravitational vector. ω is the angular velocity of the body of a hydraulic quadruped robot in an inertial coordinate system. τ is the joint torque vector. J _c is the Jacobian matrix. λ is the generalized contact force. The rotation transformation matrix R is

$$\mathit{R}=\left[\begin{array}{ccc}c{q}_{2}c{q}_3\hfill & -c{q}_{1}s{q}_3+s{q}_{1}s{q}_{2}c{q}_3\hfill & s{q}_{1}s{q}_3+c{q}_{1}s{q}_{2}c{q}_3\hfill \\ \multicolumn{1}{c}{c{q}_{2}s{q}_3}& c{q}_{1}c{q}_3+s{q}_{1}s{q}_{2}s{q}_3\hfill & -s{q}_{1}c{q}_3+c{q}_{1}s{q}_{2}s{q}_3\hfill \\ \multicolumn{1}{c}{-s{q}_2}& s{q}_{1}c{q}_2\hfill & c{q}_{1}c{q}_2\hfill \end{array}\right]$$

(5)

Where, sq₃ = sin(q₃), cq₃ = cos(q₃), and the other angles are similar.

The combined force and torque acting on the robot body are

$$\left[\begin{array}{c}{\mathit{F}}_b\hfill \\ \multicolumn{1}{c}{{\mathit{N}}_b}\end{array}\right]={\mathit{R}}^{-1}\left[\begin{array}{cc}{m}_b\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\mathit{I}}_{\mathrm{cb}}+{\mathit{I}}_{\mathrm{cb}}\mathbf{\omega }\hfill \end{array}\right]\left[\begin{array}{c}g\hfill \\ \multicolumn{1}{c}{\mathbf{\omega }}\end{array}\right]$$

(6)

In which, F _b is the combined force acting on the body of the quadruped robot, N _b is the combined torque acting on the body of the quadruped robot, m_b is the mass of the body, g is the gravitational acceleration, and I _b is the physical inertia of the robot body.

The dynamic equations of each link of the quadruped robot are as follows

$$\left[\begin{array}{c}{\mathit{F}}_{\mathrm{ik}}\hfill \\ \multicolumn{1}{c}{{\mathit{N}}_{\mathrm{ik}}}\end{array}\right]=\left[\begin{array}{cc}{m}_{\mathrm{ik}}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\mathit{I}}_{\mathrm{ick}}+{\mathit{I}}_{\mathrm{ick}}\mathbf{\omega }\hfill \end{array}\right]\left[\begin{array}{c}{\mathit{v}}_{\mathrm{ik}}\hfill \\ \multicolumn{1}{c}{{\mathbf{\omega }}_{\mathrm{ik}}}\end{array}\right]$$

(7)

Where, F _ik , N _ik , m _ik , and I _ik respectively represent the external force, torque, mass, and physical inertia of the k-th connecting rod on the i-th leg of the quadruped robot. Therein, the right front leg, right behind leg, left behind leg, and left front leg of the hydraulic robot are defined as the n-th leg, respectively, n = 1,2,3,4. The hip link, thigh link, and shank link are defined as the k-th link, respectively, k = 1,2,3.

State space model

Error model

The simplification of the model will produce model error, and the robot will also produce model error when walking in unstructured or in unknown environment. When the hydraulic quadruped robot uses the established model to perform dynamic tasks, the unmodeled part will make the robot run with deviation, which will lead to the instability of the robot movement or failure to complete the specified tasks. Therefore, in order to improve the stability of the robot motion, it is necessary to establish an error model for the hydraulic quadruped robot.

The trajectory tracking control of the hydraulic quadruped robot is realized by controlling the moving speed in the forward direction, the moving speed in the lateral direction and the yaw speed of the robot. In order to track the desired trajectory of the robot, the error between the actual trajectory and the desired trajectory should be zero. The relationship between the actual and the expected position and attitude of the robot is

$$e_q(t)=q_d(t)-q(t)$$

(8)

Where q(t) is the change of position and attitude with time of hydraulic quadruped robot, q _d (t) the desired position and attitude of the hydraulic quadruped robot, e _q (t) is the error between the expected and the actual position and attitude of the hydraulic quadruped robot.

When the motion time t of the hydraulic quadruped robot tends to infinity, the error between the desired pose and the actual pose of the robot should tend to zero.

$$\underset{t\to \infty }{lim}{e}_q(t)=0$$

(9)

According to the kinematic model of the hydraulic quadruped robot, the input of the hydraulic quadruped robot system is u = [v_x v_yω_bz]^T, The state quantity is q = [x y q₃]^T, The state space expression of the hydraulic quadruped robot system is

$$\stackrel{·}{q}=f(q,u)$$

(10)

Therefore, every point (q_d, u_d) on the reference trajectory of the hydraulic quadruped robot satisfies the following kinematic equation

$${\stackrel{·}{q}}_d=f(q_d,u_d)$$

(11)

Where q _d  = [x_d y_d q_3d]^T, u _d  = [v_xd v_ydω_bzd]^T.

The first order Taylor expansion is performed at the reference point (q _d , u _d ) of the hydraulic quadruped robot

$$\stackrel{·}{q}\approx f(q_d,u_d)+\frac{\partial f(q,u)}{\partial q}{|}_{\begin{array}{c}q=q_d\hfill \\ \multicolumn{1}{c}{u=u_d}\end{array}}(q-q_d)+\frac{\partial f(q,u)}{\partial u}{|}_{\begin{array}{c}q=q_d\hfill \\ \multicolumn{1}{c}{u=u_d}\end{array}}(u-u_d)$$

(12)

The linear error model of hydraulic quadruped robot can be obtained by combining equation (6) and equation (12), as follows

$$\begin{array}{c}\stackrel{·}{\stackrel{~}{q}}=\left[\begin{array}{c}\stackrel{·}{x}-{\stackrel{·}{x}}_{\mathrm{d}}\hfill \\ \multicolumn{1}{c}{\stackrel{·}{y}-{\stackrel{·}{y}}_{\mathrm{d}}}\\ \multicolumn{1}{c}{{\stackrel{·}{q}}_3-{\stackrel{·}{q}}_{3\mathrm{d}}}\end{array}\right]=\left[\begin{array}{ccc}0\hfill & 0\hfill & -v_{\mathrm{x}}\sin q_3-v_{\mathrm{y}}\cos q_3\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & v_{\mathrm{x}}\cos q_3-v_{\mathrm{y}}\sin q_3\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \end{array}\right]\left[\begin{array}{c}x-x_{\mathrm{d}}\hfill \\ \multicolumn{1}{c}{y-y_{\mathrm{d}}}\\ \multicolumn{1}{c}{q_3-q_{3\mathrm{d}}}\end{array}\right]\hfill \\ +\left[\begin{array}{ccc}\cos q_3\hfill & -\sin q_3\hfill & 0\hfill \\ \multicolumn{1}{c}{\sin q_3}& \cos q_3\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right]\left[\begin{array}{c}{v}_x-v_{x\mathrm{d}}\hfill \\ \multicolumn{1}{c}{v_y-v_{y\mathrm{d}}}\\ \multicolumn{1}{c}{{\omega }_{\mathrm{bz}}-{\omega }_{\mathrm{bzd}}}\end{array}\right]\hfill \end{array}$$

(13)

Equation (13) is discretized, and an augmented model based on the state space model is adopted. The form is as follows

$$\{\begin{array}{c}\mathit{q}(k+1)={\mathit{A}}_{k,t}\mathit{q}(k)+{\mathit{B}}_{k,t}\mathit{u}(k)\hfill \\ \multicolumn{1}{c}{\mathit{y}(k)={\mathit{C}}_{k,t}\mathit{q}(k)}\end{array}$$

(14)

In which, the expression of A _k,t , B _k,t , and C _k,t are as follows

$${\mathit{A}}_{k,t}=\left[\begin{array}{ccc}1\hfill & 0\hfill & (-v_x\sin q_3-v_y\cos q_3)T\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & (v_x\cos q_3-v_y\sin q_3)T\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right]$$

(15)

$${\mathit{B}}_{k,t}=\left[\begin{array}{ccc}T\cos q_3\hfill & -T\sin q_3\hfill & 0\hfill \\ \multicolumn{1}{c}{T\sin q_3}& T\cos q_3\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right]$$

(16)

$${\mathit{C}}_{k,t}=\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right]$$

(17)

Where T is sampling time.

Performance index of controller

By setting a reasonable optimization target, the deviation between the motion track and the reference track of each joint of the hydraulic quadruped robot is minimized, and then the optimal control sequence is obtained. The objective function should ensure that the quadruped robot can track the desired trajectory quickly and stably. The form of the objective function is

$$\begin{array}{c}J(k)=\sum _{i=1}^{N_p}\Vert \mathit{q}(k+i|k)-{\mathit{q}}_{\mathrm{ref}}(k+i|k){\Vert }_{\mathit{P}(i)}^2\hfill \\ +\sum _{i=1}^{N_c-1}\Vert \mathrm{\Delta }\mathit{u}(k+i|k){\Vert }_{\mathit{S}(i)}^2\hfill \end{array}$$

(18)

In which, N _c is the control time domain, q (k+i|k) is the predictive output for the quadruped robot system under control, q _ref (k+i|k) is the desired output for the quadruped robot system, Δ u (k+i|k) is the control increment of the quadruped robot system at every control moment in the predictive time domain. P (i) and S (i) are the corresponding weight matrix, and P (i) >0, S (i)>0.

The definition matrices Q (k), Q _ref (k) and Δ U (k) are as follows

$$\mathit{Q}(k)=[\mathit{q}(k+1|k),\mathit{q}(k+2|k),\cdots ,\mathit{q}(k+N_p|k){]}^{\mathrm{T}}$$

(19)

$${\mathit{Q}}_{\mathrm{ref}}(k)=[{\mathit{q}}_{\mathrm{ref}}(k+1|k),{\mathit{q}}_{\mathrm{ref}}(k+2|k),\cdots ,{\mathit{q}}_{\mathrm{ref}}(k+N_p|k){]}^{\mathrm{T}}$$

(20)

$$\mathrm{\Delta }\mathit{U}(k)=[\mathrm{\Delta }\mathit{u}(k|k),\mathrm{\Delta }\mathit{u}(k+1|k),\cdots ,\mathrm{\Delta }\mathit{u}(k+N_c-1|k){]}^{\mathrm{T}}$$

(21)

Combining equation (18)∼equation (21), equation (18) can be written as follows

$$J(k)=\Vert \mathit{Q}(k)-{\mathit{Q}}_{\mathrm{ref}}(k){\Vert }_{\mathit{P}}^2+\Vert \mathrm{\Delta }\mathit{U}(k){\Vert }_{\mathit{S}}^2$$

(22)

Where P and S are as follows

$$\mathit{P}=\left[\begin{array}{cccc}\mathit{P}(1)\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \mathit{P}(1)\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{⋮}& \ddots \hfill & \ddots \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \cdots \hfill & 0\hfill & \mathit{P}(N_p)\hfill \end{array}\right]$$

(23)

$$\mathit{S}=\left[\begin{array}{cccc}\mathit{S}(1)\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \mathit{S}(1)\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{⋮}& \ddots \hfill & \ddots \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \cdots \hfill & 0\hfill & \mathit{S}(N_p)\hfill \end{array}\right]$$

(24)

According to the linear kinematic model of the hydraulic quadruped robot, the predicted output of the quadruped robot system is

$$\mathit{Q}(k)=\Upsilon \mathit{q}(k)+\mathbf{\Omega }\mathbf{u}(k-1)+\mathbf{\Phi }\mathrm{\Delta }\mathit{U}(k)$$

(25)

Where $\Upsilon$ , $\mathbf{\Omega }$ and Φ are as follows

$$\Upsilon ={\left[\begin{array}{cccccc}{\mathit{A}}_{k,t}{\mathit{C}}_{k,t}\hfill & {\mathit{A}}_{k,t}^2{\mathit{C}}_{k,t}\hfill & \cdots \hfill & {\mathit{A}}_{k,t}^{N_c}{\mathit{C}}_{k,t}\hfill & \cdots \hfill & {\mathit{A}}_{k,t}^{N_p}{\mathit{C}}_{k,t}\hfill \end{array}\right]}^{\mathrm{T}}$$

(26)

$$\Omega ={\left[\begin{array}{ccccc}{\mathit{B}}_{k,t}{\mathit{C}}_{k,t}\hfill & {\mathit{A}}_{k,t}{\mathit{B}}_{k,t}{\mathit{C}}_{k,t}\hfill & \cdots \sum _{i=0}^{N_c}{\mathit{A}}_{k,t}^i{\mathit{B}}_{k,t}{\mathit{C}}_{k,t}\hfill & \cdots \hfill & \sum _{i=0}^{N_p-1}{\mathit{A}}_{k,t}^i{\mathit{B}}_{k,t}{\mathit{C}}_{k,t}\hfill \end{array}\right]}^{\mathrm{T}}$$

(27)

$$\mathbf{\Phi }=\left[\begin{array}{ccc}{\mathit{B}}_{k,t}{\mathit{C}}_{k,t}\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{\sum _{i=0}^{N_c}{\mathit{A}}_{k,t}^i{\mathit{B}}_{k,t}{\mathit{C}}_{k,t}}& \cdots \hfill & {\mathit{B}}_{k,t}{\mathit{C}}_{k,t}+{\mathit{B}}_{k,t}{\mathit{C}}_{k,t}\hfill \\ \multicolumn{1}{c}{⋮}& \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{\sum _{i=0}^{N_p-1}{\mathit{A}}_{k,t}^i{\mathit{B}}_{k,t}{\mathit{C}}_{k,t}}& \cdots \hfill & \sum _{i=0}^{N_p-N_c}{\mathit{A}}_{k,t}^i{\mathit{B}}_{k,t}{\mathit{C}}_{k,t}\hfill \end{array}\right]$$

(28)

When the system model is linear and discrete, the optimal control problem in the predictive time domain of the model prediction control algorithm can be transformed into the optimal control problem in the finite time domain. The quadratic programing method is used to solve the optimal control sequence of a quadruped robot with multiple constraints, The form of quadratic programing problem is as follows

$$\begin{array}{c}\underset{\mathrm{\Delta }\mathit{U}(k)}{min}J(k)=\mathrm{\Delta }\mathit{U}(k{)}^{\mathrm{T}}\mathit{M}\mathrm{\Delta }\mathit{U}(k)+{\mathit{G}}^{\mathrm{T}}\mathrm{\Delta }\mathit{U}(k)\hfill \\ \mathrm{s}.\mathrm{t}.{\mathit{H}}_i\mathrm{\Delta }\mathit{U}(k)\le {\mathit{h}}_i\hfill \end{array}$$

(29)

In which, Δ U (k) is the input the increment matrix for each time in the predictive time domain of the MPC algorithm, M and G are matrices related to the input model of the control system, H _i is the constant coefficient matrix of inequality system, h _i is the constant term vector of inequality system, The subscript i represents the i-th row or the i-th element of the matrix.

According to the model of hydraulic quadruped robot, the following definitions are given

$$\mathbf{\eta }(k)\stackrel{\mathrm{\Delta }}{=}{\mathit{Q}}_{\mathrm{ref}}(k)-\Upsilon \mathit{q}(k)-\mathbf{\Omega }\mathbf{u}(k-1)$$

(30)

Combining equation (30), the objective function of hydraulic quadruped robot can be expressed as

$$\begin{array}{c}J(k)=\Vert \mathit{Q}(k)-{\mathit{Q}}_{\mathrm{ref}}(k){\Vert }_{\mathit{P}}^2+\Vert \mathrm{\Delta }\mathit{U}(k){\Vert }_{\mathit{S}}^2\hfill \\ =\Vert \mathbf{\Phi }\mathrm{\Delta }\mathit{U}(k)-\mathbf{\eta }(k){\Vert }_{\mathit{P}}^2+\Vert \mathrm{\Delta }\mathit{U}(k){\Vert }_{\mathit{S}}^2\hfill \end{array}$$

(31)

According to the form of quadratic programing problem, the following definitions are given

$$\mathit{M}={\mathbf{\Phi }}^{\mathrm{T}}\mathit{P}\Phi +\mathit{S}$$

(32)

$$\mathit{G}=-2{\mathbf{\Phi }}^{\mathrm{T}}\mathit{P}\eta (k)$$

(33)

$$\mathrm{\Gamma }=\mathbf{\eta }(k{)}^{\mathrm{T}}\mathit{P}\eta (k)$$

(34)

Where, Г is a constant.

Combining equation (31) to equation (34), the objective function of the hydraulic quadruped robot can be expressed as

$$J(k)=\mathrm{\Gamma }+\mathrm{\Delta }\mathit{U}(k{)}^{\mathrm{T}}\mathit{G}+\mathrm{\Delta }\mathit{U}(k{)}^{\mathrm{T}}\mathit{M}\mathrm{\Delta }\mathit{U}(k)$$

(35)

From equation (35), it can be seen that the objective function of the trajectory tracking of the hydraulic quadruped robot based on the model predictive control is transformed into the quadratic programing form of the input increment matrix Δ U (k) in the predictive time domain.

Inequality constraints

In the control process of the hydraulic quadruped robot, due to the complexity of the robot working environment, the input constraints of the robot leg joints, the joint angle range constraints, the motion speed constraints of the torso and the motion angle speed constraints of the torso should be considered.

Considering the motion safety of the hydraulic quadruped robot, the motion speed and angular velocity of the robot torso should be within a certain range, that is, the motion speed and angular velocity of the hydraulic quadruped robot should meet the following constraints

$$\{\begin{array}{c}{v}_{xmin}\le v_x\le v_{xmax}\hfill \\ v_{ymin}\le v_y\le v_{ymax}\hfill \\ {\omega }_{\mathrm{bz}min}\le {\omega }_{\mathrm{bz}}\le {\omega }_{\mathrm{bz}max}\hfill \end{array}$$

(36)

In which, v _{x min} and v _{x max} are the minimum and maximum values of the forward velocity of the hydraulic quadruped robot, respectively. v _{y min} and v _{y max} are the minimum and maximum values of the lateral velocity of the hydraulic quadruped robot, respectively. ω_bzmin and ω_bzmax are the minimum and maximum values of the yaw angular velocity of the torso of the hydraulic quadruped robot, respectively.

Combining the motion stability of the hydraulic robot, the acceleration and angular acceleration of the robot torso should be within a certain range, that is, the acceleration a _vx , a _vxy , and angular acceleration a _{ω bz} of the hydraulic quadruped robot should meet the following constraints

$$\{\begin{array}{c}{a}_{\mathrm{vx}min}\le a_{\mathrm{vx}}\le a_{\mathrm{vx}max}\hfill \\ a_{\mathrm{vy}min}\le a_{\mathrm{vy}}\le a_{\mathrm{vy}max}\hfill \\ a_{\omega \mathrm{bz}min}\le a_{\omega \mathrm{bz}}\le a_{\omega \mathrm{bz}max}\hfill \end{array}$$

(37)

In which, a _{vx min} and a _{vx max} are the minimum and maximum values of the forward acceleration of the hydraulic quadruped robot, respectively. a _{vy min} and a _{vy max} are the minimum and maximum values of the lateral acceleration of the hydraulic quadruped robot, respectively. a _{ω bzmin} and a _{ω bzmax} are the minimum and maximum values of the yaw angular acceleration of the torso of the hydraulic quadruped robot, respectively.

In each control cycle of the hydraulic quadruped robot, the control quantity constraint can be expressed as

$$\begin{array}{c}{\mathit{u}}_{min}(t+k)\le \mathit{u}(t+k)\le {\mathit{u}}_{max}(t+k)\hfill \\ \multicolumn{1}{c}{(k=0,1,\dots ,N_c-1)}\end{array}$$

(38)

Where, u _min (t + k) and u _max (t + k) are the minimum and the maximum of the control quantity of the hydraulic quadruped robot

The acceleration incremental constraint of the hydraulic quadruped robot can be transformed into velocity incremental constraint and angular velocity incremental constraint, hence the control incremental constraint can be expressed as

$$\begin{array}{c}\mathrm{\Delta }\begin{array}{c}{\mathit{u}}_{min}(t+k)\le \mathrm{\Delta }\mathit{u}(t+k)\le \mathrm{\Delta }{\mathit{u}}_{max}(t+k)\hfill \end{array}\hfill \\ (k=0,1,\dots ,N_c-1)\hfill \end{array}$$

(39)

Where, Δ u _min (t + k) and Δ u _max (t + k) are the minimum control increment and the maximum control increment of the hydraulic quadruped robot.

Optimal control based on constraint priority

When considering the multi-constraints such as the velocity and angular velocity of the robot at the same time, the optimal control problem may have no solution. The slack variable ε is introduced to ensure the existence of optimal control solutions under multiple constraints. The trajectory tracking problem under multiple constraint conditions of quadruped robot can be represented as

$$\begin{array}{c}\underset{\mathrm{\Delta }\mathit{U}(k)}{min}J(k)=\mathrm{\Delta }\mathit{U}(k{)}^T{\mathit{M}}^{*}\mathrm{\Delta }\mathit{U}(k)+{{\mathit{G}}^{*}}^T\mathrm{\Delta }\mathit{U}(k)+\rho {\epsilon }^2\hfill \\ \mathrm{s}.\mathrm{t}.\hfill \\ \begin{array}{cc}{\mathit{u}}_{min}\le \mathit{u}\le {\mathit{u}}_{max}\hfill & (k=0,1,\dots ,N_c-1)\hfill \\ \multicolumn{1}{c}{\mathrm{\Delta }{\mathit{u}}_{min}\le \mathrm{\Delta }\mathit{u}\le \mathrm{\Delta }{\mathit{u}}_{max}}& (k=0,1,\dots ,N_c-1)\hfill \end{array}\hfill \end{array}$$

(40)

Where ρ is the weight coefficients, M * and G * are matrices related to the input model of the control system.

The objective function for the quadruped robot system in the predicting time domain as follows

$$\begin{array}{c}J(k)={[\mathrm{\Delta }{\mathit{U}}^T(k),\epsilon ]}^T{\mathit{M}}^{*}[\mathrm{\Delta }{\mathit{U}}^T(k),\epsilon ]+{\mathit{G}}^{*}[\mathrm{\Delta }{\mathit{U}}^T(k),\epsilon ]\hfill \\ \mathrm{s}.\mathrm{t}.\hfill \\ {\mathit{U}}_{min}\le \mathrm{\Delta }\mathit{U}(k)\mathrm{\Delta }{\mathit{U}}_{max}\hfill \\ \mathrm{\Delta }{\mathit{U}}_{min}\le {\mathit{H}}_i\mathrm{\Delta }\mathit{U}(k)+\mathit{U}(k)\mathrm{\Delta }{\mathit{U}}_{max}\hfill \end{array}$$

(41)

The expressions for M * and G * are respectively

$${\mathit{M}}^{*}=\left[\begin{array}{cc}{\mathbf{\Phi }}^{\mathrm{T}}{\mathbf{\Gamma }}_p\mathbf{\Phi }+\mathit{S}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \rho \hfill \end{array}\right]$$

(42)

$${\mathit{G}}^{*}=\left[\begin{array}{cc}2e_k^{\mathrm{T}}{\mathbf{\Gamma }}_u^{\mathrm{T}}\mathbf{\Phi }\hfill & 0\hfill \end{array}\right]$$

(43)

Where, ${\mathrm{\Gamma }}_u$ is the input error weighting matrix of the quadruped robot system, ${\mathrm{\Gamma }}_p$ is the output error weighting matrix of the quadruped robot system, e _k is the tracking error of the robot in the predicted time domain.

The corresponding control increment of the robot in the control time domain can be obtained according to the objective function. The input control increment of the system is the first control increment. And the control process is repeated until the trajectory tracking control of the hydraulic robot is completed.

$$\mathit{u}(t)=\mathit{u}(t-1)+\mathrm{\Delta }\mathit{u}(t)$$

(44)

Stability analysis

The stability of the system is proved using the Lyapunov stability theorem. Combining equation (18) and equation (40), the Lyapunov function is constructed as follows:

$$\begin{array}{c}{V}_N(k)=min[\sum _{i=1}^{N_p}\left(\Vert \mathit{q}(k)-{\mathit{q}}_{\mathrm{ref}}(k){\Vert }_{\mathit{P}(i)}^2\right)\hfill \\ +\sum _{i=0}^{N_c-1}\left(\Vert \mathrm{\Delta }\mathit{u}(k){\Vert }_{\mathit{S}(i)}^2+\Vert \mathit{\epsilon }(k){\Vert }_{\rho }^2\right)]\hfill \end{array}$$

(45)

The function (45) satisfies the following conditions

$$\{\begin{array}{c}{V}_N(k)=0,k=0\hfill \\ \multicolumn{1}{c}{V_N(k)>0,k\ne 0}\end{array}$$

(46)

From equation (19) and equation (21)

$$\begin{array}{c}\tilde{\mathit{u}}(k+1+i|k+1)=[\tilde{\mathit{u}}(k+1|k+1),\tilde{\mathit{u}}(k+2|k+1),\hfill \\ \dots ,\tilde{\mathit{u}}(k+N_c|k+1)]\hfill \end{array}$$

(47)

$$\begin{array}{c}△\tilde{\mathit{u}}(k+1+i|k+1)=[△\tilde{\mathit{u}}(k+1|k+1),\hfill \\ △\tilde{\mathit{u}}(k+2|k+1),\dots ,△\tilde{\mathit{u}}(k+N_c|k+1)]\hfill \end{array}$$

(48)

$$\mathit{\epsilon }(k+1)={\mathit{\epsilon }}_{N}k$$

(49)

In which, $\tilde{\mathit{u}}$ is the optimal control quantity of the system, $\mathrm{\Delta }\tilde{u}$ is the optimal control increment of the system, and $\mathit{\epsilon }(k+1)$ is slack variable.

From equation (45) to equation (48)

$$\begin{array}{c}\tilde{\mathit{u}}(k+1+i|k+1)=[{\tilde{\mathit{u}}}_N(k+1|k),{\tilde{\mathit{u}}}_N(k+2|k),\hfill \\ \cdots ,{\tilde{\mathit{u}}}_N(k+N_c|k)]\hfill \end{array}$$

(50)

$$\begin{array}{c}△\tilde{\mathit{u}}(k+1+i|k+1)=[△{\tilde{\mathit{u}}}_N(k+1|k),△{\tilde{\mathit{u}}}_N(k+2|k),\hfill \\ \cdots ,△{\tilde{\mathit{u}}}_N(k+N_c|k)]\hfill \end{array}$$

(51)

Based on the above analysis, it can be obtained that

$$\begin{array}{c}J(k+1)=\sum _{i=1}^{N_p}\Vert \mathit{q}(k+1+i|k+1)-{\mathit{q}}_{\mathrm{ref}}(k+1+i|k+1){\Vert }_{\mathit{P}(i)}^2\hfill \\ +\sum _{i=0}^{N_c-1}(\Vert \mathrm{\Delta }\mathit{u}(k+1+i|k+1){\Vert }_{\mathit{S}(i)}^2+\Vert \mathit{\epsilon }(k+1+i|k+1){\Vert }_{\rho }^2)\hfill \end{array}$$

(52)

Combining equation (45), equation (50), and equation (51), equation (52) can be transformed into

$$\begin{array}{c}J(k+1)=\sum _{i=1}^{N_p}\Vert {\mathit{q}}_N(k+i|k+1)-{\mathit{q}}_{\mathrm{ref}}(k+i|k){\Vert }_{\mathit{P}(i)}^2\hfill \\ +\sum _{i=0}^{N_c-1}(\Vert \mathrm{\Delta }{\mathit{u}}_N(k+i|k){\Vert }_{\mathit{S}(i)}^2+\Vert {\mathit{\epsilon }}_N(k+i|k){\Vert }_{\rho }^2)-\mathit{D}\hfill \\ =V_N\left(k\right)-\mathit{D}\hfill \end{array}$$

(53)

In which,

$$\mathit{D}=\Vert \mathit{q}(k|k)-{\mathit{q}}_{\mathrm{ref}}(k|k){\Vert }_{\mathit{P}(i)}^2+\Vert \mathrm{\Delta }\mathit{u}(k|k){\Vert }_{\mathit{S}(i)}^2+\Vert \mathit{\epsilon }(k|k){\Vert }_{\rho }^2$$

(54)

Due to V _N (k+1) is the optimal solution of the function, it can be obtained that

$$V_N(k+1)\le J(k+1)$$

(55)

From equation (53),

$$V(k+1)\le J(k+1)\le V(k)-\mathit{D}$$

(56)

According to equation (53), $\mathit{D}\ge 0$ , so it can be inferred that

$$V(k+1)\le V(k)$$

(57)

From the above analysis, it can be concluded that the system is stable.

Design of trajectory tracking controller

Parameters of trajectory tracking

The trajectory of the foot in the corresponding leg transition coordinate system is

$$\begin{array}{c}\underset{\xi }{\mathrm{minimize}}\frac{1}{2}{\xi }^{\mathrm{T}}\mathit{K}\xi +{\mathit{c}}^{\mathrm{T}}\xi \hfill \\ \begin{array}{cc}\mathrm{subject}\hfill & \begin{array}{cc}\mathrm{to}\hfill & \mathit{A}\hfill \end{array}\hfill \end{array}\xi =\mathit{b}\begin{array}{cc},\hfill & \mathit{D}\xi \hfill \end{array}\le \mathit{f}\hfill \end{array}$$

(58)

In which, ξ is the vector composed of optimize variables, K and c are matrices of appropriate dimension. when Kξ = b and Kξ ≤ f , the quadratic objective ξ is the Minimum.

The foot trajectory of swinging leg could be represented as

$$\mathit{r}(t)=\left[\begin{array}{ccc}{\zeta }^{\mathrm{T}}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\zeta }^{\mathrm{T}}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & {\zeta }^{\mathrm{T}}\hfill \end{array}\right]\left[\begin{array}{c}{\alpha }_x\hfill \\ \multicolumn{1}{c}{{\alpha }_y}\\ \multicolumn{1}{c}{{\alpha }_z}\end{array}\right]=\mathit{T}(t)\xi$$

(59)

where ζ T is the basic form of trajectory, α* is the coefficient of the trajectory. The trajectory of the robot along the x, y, and z motion directions is optimized to ensure that the coefficient is included in the vector ξ =[ α _x ^T α _y ^T α _z ^T] ^T.

The velocity of the supporting end in the x-direction of corresponding leg transition coordinate system is a polynomial, and the velocity in the other directions is zero, the velocity vector of the foot is

$$\stackrel{·}{\mathit{r}}(t)=\left[\begin{array}{ccc}1\hfill & t\hfill & t\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \end{array}\right]\left[\begin{array}{c}{\alpha }_{\mathrm{s}x}\hfill \\ \multicolumn{1}{c}{{\alpha }_{\mathrm{s}y}}\\ \multicolumn{1}{c}{{\alpha }_{\mathrm{s}z}}\end{array}\right]$$

(60)

Based on the initial position of the robot support phase, the motion position vector of the foot of the robot support phase could be obtained

$$\mathit{r}(t)=\left[\begin{array}{c}x(T_m)\hfill \\ \multicolumn{1}{c}{y(T_m)}\\ \multicolumn{1}{c}{0}\end{array}\right]+{\int }_{T_{\mathrm{m}}}^{T_{\mathrm{m}}+T_{\mathrm{s}}}\mathit{R}(t{\mathit{q}}_{3\mathrm{ref}})\stackrel{·}{\mathit{r}}(t)\mathrm{d}t$$

(61)

Where, q _3ref = [0 0 q_3ref] ^T, R (t q _3ref) is a rotation matrix, which describes the direction change of the foot caused by the expected yaw rate change.

$$\mathit{r}(t)=\mathit{T}(q_{3\mathrm{ref}},t)\xi$$

(62)

Where, the coefficient of the optimization variable is

$$\xi =\left[\begin{array}{ccccc}{\alpha }_{\mathrm{s}x}\hfill & {\alpha }_{\mathrm{s}y}\hfill & {\alpha }_{\mathrm{s}z}\hfill & x(T_m)\hfill & y(T_m)\hfill \end{array}\right]$$

(63)

Otherwise, the speed and acceleration of the supporting foot are

$$\stackrel{·}{\mathit{r}}(t)=\stackrel{·}{\mathit{T}}(q_{3\mathrm{ref}},t)\xi$$

(64)

$$\stackrel{··}{\mathit{r}}(t)=\stackrel{··}{\mathit{T}}(q_{3\mathrm{ref}},t)\xi$$

(65)

Target of trajectory tracking

The optimization objectives for robot trajectory tracking are as follows

Minimize of the acceleration

The smooth motion of the robot could be ensured by the minimization of acceleration. The objective function of the acceleration for supporting foot during movement time is

$$\frac{1}{2}{\xi }^{\mathrm{T}}(2{\int }_{T_{\mathrm{m}}}^{T_{\mathrm{m}}+T_{\mathrm{s}}}{\stackrel{··}{\mathit{T}}}^{\mathrm{T}}(t){\mathit{W}}_{\mathrm{acc}}\stackrel{··}{\mathit{T}}(t)\mathrm{d}t)\xi$$

(66)

If

$${\mathit{Q}}_{\mathrm{acc}}=2{\int }_{T_{\mathrm{m}}}^{T_{\mathrm{m}}+T_{\mathrm{s}}}{\stackrel{··}{\mathit{T}}}^{\mathrm{T}}(t){\mathit{W}}_{\mathrm{acc}}\stackrel{··}{\mathit{T}}(t)\mathrm{d}t$$

(67)

In which, Q _acc is Hessian matrix. Combining equations (61) and (62), it can be seen that Q _acc is determined by the expected yaw rate of the robot, and W acc is the corresponding weight matrix.

The linear term of the acceleration quadratic programing problem for quadruped robot foot is zero.

$${\mathit{c}}^{\mathrm{T}}\xi =0$$

(68)

Similarly, the objective function of the acceleration of the swing foot could be obtained.

Minimize of tracking errors of position, velocity and acceleration

The actual motion trajectory of a quadruped robot deviates significantly from the expected motion trajectory, which can result in the robot being unable to complete the expected task.

To ensure the stable motion of the quadruped robot, the optimization variable for the position error of the robot’s foot end in the control time domain N _p is

$$\sum _{k=1}^{N_p}\Vert \mathit{r}(t_k)-{\mathit{r}}_{\mathrm{ref}}(t_k){\Vert }_{{\mathit{W}}_{\mathrm{ref}}}^2\mathrm{\Delta }t$$

(69)

In which, r _ref (t _k ) is the expected position of the robot foot at sampling time, t _k is the sampling time, and W _ref is the corresponding weight matrix. Where, Δt = t _k −t _{k −1}.

The optimization variables for the velocity error and acceleration of the robot are

$$\sum _{k=1}^{N_p}\Vert \stackrel{·}{\mathit{r}}(t_k)-{\stackrel{·}{\mathit{r}}}_{\mathrm{ref}}(t_k){\Vert }_{{\mathit{W}}_{\mathrm{ref}}}^2\mathrm{\Delta }t$$

(70)

$$\sum _{k=1}^{N_p}\Vert \stackrel{··}{\mathit{r}}(t_k)-{\stackrel{··}{\mathit{r}}}_{\mathrm{ref}}(t_k){\Vert }_{{\mathit{W}}_{\mathrm{ref}}}^2\mathrm{\Delta }t$$

(71)

Minimize of the expected speed of the supporting phase

The robot should not slip during movement, and the foot speed should track the expected speed. The objective function for minimizing the velocity of the foot is

$$\frac{1}{2}{\xi }^{\mathrm{T}}(2q_{3\mathrm{ref}}{\mathit{\Psi }}^{\mathrm{T}}\mathit{\Psi })\xi +(-2q_{3\mathrm{ref}}{v}_{\mathrm{xref}}\mathit{\Psi })\xi$$

(72)

Combining with the standard form of secondary planning, it can be concluded that

$${\mathit{K}}_{\mathrm{ref}}=2q_{3\mathrm{ref}}{\mathit{\Psi }}^{\mathrm{T}}\mathit{\Psi }$$

(73)

$${\mathit{c}}_{\mathrm{ref}}^{\mathrm{T}}=-2q_{3\mathrm{ref}}{v}_{\mathrm{xref}}\mathit{\Psi }$$

(74)

Where $\mathit{\Psi }=\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \end{array}\right]\stackrel{·}{\mathit{T}}(q_{3\mathrm{ref}},0)$ .

Minimization of tracking error of the leg height

When the robot is in the swing phase, in order to ensure that the robot’s leg lift height tracks the desired height, the objective function for minimizing the leg lift height at the foot end is

$$\frac{1}{2}{\xi }^{\mathrm{T}}(2{\omega }_H{\mathit{\Psi }}^{\mathrm{T}}\mathit{\Psi })\xi +(-2{\omega }_H{H}_{\mathrm{ref}}\mathit{\Psi })\xi$$

(75)

In which,

$$K_{\mathrm{H}}=2{\omega }_{\mathrm{H}}{\Psi }^{\mathrm{T}}\Psi$$

(76)

$${\mathit{c}}_{\mathrm{H}}^{\mathrm{T}}=-2{\omega }_H{H}_{\mathrm{ref}}\Psi$$

(77)

Where, $\mathit{\Psi }=\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \end{array}\right]\stackrel{·}{\mathit{T}}(q_{3\mathrm{ref}},0)$ . And t_H is the time when the foot at its highest point.

Constraints for trajectory tracking

To achieve trajectory tracking control of the foot end, the constraint conditions are as follows

Initial state of the supporting phase

Quadruped robots should effectively cope with their own disturbances and external disturbances. The constraints of initial state for the support phase of each motion cycle of the robot as follows

$$\mathit{T}(0)\xi ={\left[\begin{array}{ccc}{x}_{\mathrm{in}}\hfill & y_{\mathrm{in}}\hfill & 0\hfill \end{array}\right]}^{\mathrm{T}}$$

(78)

In the formula, x _in and y _in represent the actual positions of the foot, respectively.

Continuity of the trajectory

The position and speed of the foot of the quadruped robot should be continuous. At the end of the swing phase and the beginning of the support phase, the position and velocity constraints of the foot are

$$\left[\begin{array}{cc}\mathit{T}(T_{\mathrm{m}})\hfill & -\mathit{T}(T-T_{\mathrm{s}})\hfill \\ \multicolumn{1}{c}{\stackrel{·}{\mathit{T}}(T_{\mathrm{m}})}& -\stackrel{·}{\mathit{T}}(T-T_{\mathrm{s}})\hfill \end{array}\right]\left[\begin{array}{c}\xi \hfill \\ \multicolumn{1}{c}{\xi }\end{array}\right]=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]$$

(79)

Trajectory tracking controller

The purpose of trajectory tracking control for hydraulic robots is to achieve speedy and stable motion along the desired trajectory. The principle of trajectory tracking control is shown in Figure 2. The motion parameters of the robot are given, and the position and attitude information of the actuators is obtained through the tracking model. Based on the error model of the system, the joint control quantities are calculated using the trajectory tracking controller, and the quadruped robot moved along the desired trajectory.

Figure 2.

The trajectory tracking controller of quadruped robot.

Simulation

To verify the effectiveness of the trajectory tracking algorithm, by establishing the simulation model, the hydraulic quadruped robot underwent the trajectory tracking control experiments in trotting gait, and pacing gait, respectively. And the comparative analysis of the control effects between not using trajectory tracking control method and using the proposed trajectory tracking control method were completed.

Trotting gait simulation experiment

The lifting height of the foot is 70 mm, the step length is 260 mm, the running cycle is 0.8 s, and the running speed of the robot is 0.325 m/s. The process of the hydraulic quadruped robot in trotting gait is shown in Figure 3. The horizontal direction displacement curve of the torse of the hydraulic robot in trotting gait is shown in Figure 4, and the vertical displacement curve is shown in Figure 5. The figure includes three kinds of displacement data for the hydraulic robot, the target displacement, the original displacement without trajectory tracking control method, and the optimal displacement with trajectory tracking control method, respectively.

Figure 3.

The simulation process of the hydraulic quadruped robot in trotting gait.

Figure 4.

The horizontal direction displacement of the hydraulic quadruped robot torso in trotting gait.

Figure 5.

The vertical direction displacement of the hydraulic quadruped robot torso in trotting gait.

The simulation experiment results in Figures 4 and 5 could thus be obtained. As shown in Figure 4, the target displacement of the hydraulic robot in the forward direction within 20 s is 6.5 m. The original displacement without trajectory tracking control method of the robot is 6.446 m, and the trajectory tracking error is 0.83%. The optimal displacement with the proposed trajectory tracking control method of the robot is 6.477 m, and the trajectory tracking error is 0.35%. The trajectory tracking error with the proposed method of the hydraulic robot in the forward direction is significantly reduced.

According to Figure 5, the target displacement in the vertical direction is 0.78 m. The original displacement variation in the vertical direction without trajectory tracking control method is 0.735–0.821 m, and the variation range of displacement is 87 mm. The optimal displacement variation with trajectory tracking control method is 0.753–0.805 m, and the variation range of displacement is only 52 mm. The vertical displacement variation of the hydraulic robot is reduced, and the trajectory tracking error is optimized. The motion trajectory of hydraulic robot can better track the expected trajectory by the proposed method. The variation range of the center of mass of the hydraulic robot torso in the vertical direction is smaller, and the stability of the robot torso is improved.

Pacing gait simulation experiment

The lifting height of the foot is 70 mm, the step length is 260 mm, the running cycle is 0.8 s, and the running speed of the robot is 0.325 m/s. The process of the hydraulic quadruped robot in pacing gait is shown in Figure 6. The horizontal direction displacement curve of the torse of the hydraulic robot in pacing gait is shown in Figure 7, and the vertical displacement curve is shown in Figure 8.

Figure 6.

The simulation process of the hydraulic quadruped robot in pacing gait.

Figure 7.

The horizontal direction displacement of the hydraulic quadruped robot torso in pacing gait.

Figure 8.

The vertical direction displacement of the hydraulic quadruped robot torso in pacing gait.

As shown in Figure 7, the target displacement of the hydraulic robot in the forward direction within 18 s is 4.5 m. The original displacement without trajectory tracking control method is 4.442 m, and the trajectory tracking error is 1.29%. The optimal displacement with the proposed trajectory tracking control method of the robot is 4.469 m, and the trajectory tracking error is 0.68%. The trajectory tracking error with the proposed method of the hydraulic robot in the forward direction in pacing gait is optimized.

As shown in Figure 8, the target displacement in the vertical direction is 0.78 m. The original displacement in the vertical direction in pacing gait without trajectory tracking control method is 0.767–0.792 m, and the variation range of displacement is 25 mm. The optimal displacement variation with trajectory tracking control method is 0.773–0.787 m, and the variation range of displacement is just 14 mm. The vertical displacement variation of the hydraulic robot in pacing gait is reduced, and the trajectory tracking error is optimized.

Based on the above simulation experiment results, it can be concluded that the error between the actual motion trajectory and the expected trajectory of the hydraulic quadruped robot by the proposed method is very small. The motion of hydraulic robots is stable. The effectiveness of the proposed trajectory tracking control algorithm based on model prediction is demonstrated.

Experiments

Experimental parameters

The trajectory tracking error of the quadruped robot during motion is analyzed. The lengths of the hip link, thigh link and shank link of the robot are 60, 340, and 370 mm, respectively. Experimental parameters are set based on the structural parameters and workspace of the hydraulic quadruped robot.

In order to verify the adaptability of the control strategy, a relatively large range of constraint parameters is set. The constraints of the robot system are shown in Tables 1 and 2. The control parameters for trajectory tracking are shown in Table 3.

Table 1.

The bound of velocity of hydraulic quadruped robot.

Speeds	u min	u max
v x	1 m/s	−1 m/s
v y	0.5 m/s	−0.5 m/s
ω bz	0.5 rad/s	−0.5 rad/s

Table 2.

The bound of velocity change of hydraulic quadruped robot.

Speed increments	Δumin	Δumax
Δv x	0.3 m/s	−0.3 m/s
Δv y	0.3 m/s	−0.3 m/s
Δω bz	0.3 rad/s	−0.3 rad/s

Table 3.

The control parameters based on model predictive control of the robot.

Parameters	Values
Prediction step	80
Control steps	20
Control cycle	0.8 s
Weight coefficient	ρ = 10
Slack variable	ε = 10

To verify the effectiveness and robustness of the trajectory tracking strategy for the hydraulic actuated quadruped robot, the experiments under different speeds of the quadruped robot are conducted.

Experimental analysis under 0.225 m/s

The trajectory tracking control experiments in trotting gait are conducted by hydraulic quadruped robot experimental platform, including experiments without trajectory tracking control method and experiments with trajectory tracking control method. The process of the hydraulic quadruped robot in trotting gait is shown in Figure 9.

Figure 9.

The process of the hydraulic quadruped robot in trotting gait.

The displacement data of the leg actuators for the right front leg, right behind leg, left behind leg, and left front leg of the hydraulic robot are shown in Figures 10 to 13, respectively. Figure (a) shows the displacement data of the hip joint actuator of the hydraulic robot, and Figure (b) shows the displacement data of the knee joint actuator. There are three sets of data in each figure, which are the command displacement data of the actuator, the original data of the joint actuator without trajectory tracking control method, and the optimal data of the joint actuator using the trajectory tracking control method proposed in this paper.

Figure 10.

Displacement tracking experiment of right front leg of the hydraulic robot under 0.225 m/s: (a) hip joint and (b) knee joint.

Figure 11.

Displacement tracking experiment of right hind leg of the hydraulic robot under 0.225 m/s: (a) hip joint and (b) knee joint.

Figure 12.

Displacement tracking experiment of left hind leg of the hydraulic robot under 0.225 m/s: (a) hip joint and (b) knee joint.

Figure 13.

Displacement tracking experiment of left front leg of the hydraulic robot under 0.225 m/s: (a) hip joint and (b) knee joint.

The absolute values of the original trajectory maximum tracking errors for the hip joint actuators are 1.62, 1.032, 0.809, and 0.808 mm, respectively. And the absolute values of the optimal trajectory maximum tracking errors are 0.719, 0.432, 0.31, and 0.309 mm, respectively. The absolute values of the original trajectory maximum tracking errors for the knee joint actuators are 0.978, 0.797, 1.037, and 1.796 mm, respectively. And the absolute values of the optimal trajectory maximum tracking errors by the proposed method are 0.377, 0.346, 0.536, and 0.396 mm, respectively. The trajectory tracking error of the hip joint actuator of the hydraulic robot is reduced by the proposed trajectory tracking method.

Comparing the original data and the optimal data with the command data, the root mean square tracking errors of each joint actuator of the hydraulic quadruped robot are obtained. Table 4 shows the mean square error of hip joint actuators, and Table 5 shows the mean square error of knee joint actuators.

Table 4.

Root mean square error of hip joints under 0.225 m/s.

Legs	Right front	Right hind	Left hind	Left front
Original	2.986	3.538	3.134	5.131
Optimal	0.939	0.980	1.871	2.515

Table 5.

Root mean square error of knee joints under 0.225 m/s.

Legs	Right front	Right hind	Left hind	Left front
Original	1.494	1.482	1.411	2.148
Optimal	0.650	0.615	0.650	0.589

From the tables, it can be seen that using the proposed model predictive control method, the hip joint actuator errors of the right front, right hind, left hind, and left front are reduced by 68.55%, 72.3%, 40.28%, and 50.99%, respectively. And the knee joint actuator errors of the right front, right hind, left hind, and left front are reduced by 56.48%, 58.5%, 53.93%, and 72.58%, respectively.

Figure 14 gives the values of attitude angle of the hydraulic quadruped robot. As shown in Figure 14, by using the proposed trajectory tracking control method, the range of rolling angle for the hydraulic robot is −0.066 to 0.066 rad, and the range of pitching angle is −0.073 to 0.044 rad.

Figure 14.

Curves of attitude angle of the hydraulic robot under 0.225m/s: (a) rolling angle and (b) pitching angle.

Experimental analysis under 0.315 m/s

The displacement data of the leg actuators the hydraulic robot under 0.315 m/s are shown in Figures 15 to 18, respectively. From the figures, it can be seen that the trajectory tracking of the joint actuators of the hydraulic robot are closer to the expected trajectory by the proposed trajectory tracking method.

Figure 15.

Displacement tracking experiment of right front leg of the hydraulic robot under 0.315 m/s: (a) hip joint and (b) knee joint.

Figure 16.

Displacement tracking experiment of right hind leg of the hydraulic robot under 0.315 m/s: (a) hip joint and (b) knee joint.

Figure 17.

Displacement tracking experiment of left hind leg of the hydraulic robot under 0.315 m/s: (a) hip joint and (b) knee joint.

Figure 18.

Displacement tracking experiment of left front leg of the hydraulic robot under 0.315 m/s: (a) hip joint and (b) knee joint.

The absolute values of the original trajectory maximum tracking errors for the hip joint actuators are 1.54, 1.23, 0.84, and 0.956 mm respectively. And the absolute values of the optimal trajectory maximum tracking errors are 0.692, 0.465, 0.323, and 0.454 mm, respectively. The absolute values of the original trajectory maximum tracking errors for the knee joint actuators are 0.945, 0.827, 1.206, and 1.656 mm, respectively. And the absolute values of the optimal trajectory maximum tracking errors by the proposed method are 0.415, 0.387, 0.517, and 0.618 mm, respectively.

Comparing the original data and the optimal data with the command data, the root mean square tracking errors of each joint actuator of the hydraulic quadruped robot are obtained. Table 6 shows the mean square error of hip joint actuators, and Table 7 shows the mean square error of knee joint actuators.

Table 6.

Root mean square error of hip joints under 0.315 m/s.

Legs	Right front	Right hind	Left hind	Left front
Original	3.133	3.794	3.131	5.222
Optimal	1.365	1.324	1.098	1.853

Table 7.

Root mean square error of knee joints under 0.315 m/s.

Legs	Right front	Right hind	Left hind	Left front
Original	1.240	1.494	1.807	1.999
Optimal	0.686	0.674	0.865	0.547

From the tables, it can be seen that using the proposed model predictive control method, the hip joint actuator errors of the right front, right hind, left hind, and left front are reduced by 56.45%, 65.11%, 64.95%, 64.52%, respectively. And the knee joint actuator errors of the right front, right hind, left hind, and left front are reduced by 44.63%, 54.87%, 52.13%, 72.64%, respectively.

Figure 19 gives the values of attitude angle of the hydraulic quadruped robot. As shown in Figure 19, by using the proposed trajectory tracking control method, the range of rolling angle for the hydraulic robot is −0.023 to 0.086 rad, and the range of pitching angle is −0.031 to 0.1 rad. From the comparison between the original data and optimal data in the figures, it can be seen that the stability of the robot can be improved by the proposed trajectory tracking control method.

Figure 19.

Curves of attitude angle of the hydraulic robot under 0.315m/s: (a) rolling angle and (b) pitching angle.

From the above experimental results, it can be concluded that the hydraulic robot can track the command trajectory better by the proposed trajectory tracking control method. And the error between the optimal displacement and the expected displacement of each actuator is smaller. The trajectory tracking accuracy and the stability of hydraulic robot is improved by the proposed method. And the effectiveness of the proposed control method is proved.

Conclusions

In order to meet the high precision trajectory tracking performance requirements of hydraulic robot. Aiming at the multi coupling characteristics of hydraulic robot, a decomposed trajectory tracking control strategy with trajectory constraints for hydraulic robot is proposed. The effectiveness of proposed strategy was verified through simulation and prototype experiments. This paper mainly completes the following four aspects:

• (1) The tracking model of the hydraulic robot was derived. Considering the inequality constraints during the motion of the quadruped robot, the trajectory tracking control problem of the robot was transformed into a quadratic programing problem. The trajectory tracking controller with trajectory constraints was designed, and the accuracy of the robot model in the control algorithm was improved.

• (2) Combining the minimization goals of robot motion speed and acceleration error, stability and continuity constraints for robot motion were designed, the impact during hydraulic robot motion was reduced, and the service life of the robot was increased, as well as the stability of robot motion was enhanced.

• (3) The impact of classical control strategies and the proposed control strategy on robot performance was compared through simulation. The motion trajectory error of the robot is relatively small by the proposed trajectory tracking control strategy, and the effectiveness of the control strategy in simulation was proved.

• (4) The prototype experiments on trajectory tracking control of hydraulic robots were conducted. And the effectiveness and robustness of the proposed control algorithm in practical quadruped robot system was verified by the experimental results.

The above results provide a basis for further research on motion control algorithms and hydraulic quadruped robot control. Future work will involve the detailed design of stability strategies for robot, as well as conducting application studies on actual hydraulic robots.