Constrained Predictive Tracking Control for Unmanned Hexapod Robot with Tripod Gait

Abstract

Since it is difficult to accurately track reference trajectories under the condition of stride constraints for an unmanned hexapod robot moving with rhythmic gait, an omnidirectional tracking strategy based on model predictive control and real-time replanning is proposed in this paper. Firstly, according to the characteristic that the stride dominates the rhythmic motion of an unmanned multi-legged robot, a body-level omnidirectional tracking model is established. Secondly, a quantification method of limb’s stretch and yaw constraints described by motion stride relying on a tripod gait is proposed, and then, a body-level accurate tracking controller based on constrained predictive control is designed. Then, in view of the low tracking efficiency of the robot under the guidance of common reference stride, a solution strategy of variable stride period and a real-time replanning scheme of reference stride are proposed based on the limb constraints and the integral mean, which effectively avoid the tracking deviation caused by the guidance of constant reference strides. Finally, the effectiveness and practicability of the proposed control strategy are demonstrated through the comparative analysis and simulation test of a hexapod robot WelCH with omnidirectional movement ability to continuously track the directed curve and the undirected polyline trajectory.

1. Introduction

As an important member of unmanned ground mobile devices, unmanned multi-legged robots have outstanding performance and broad application prospects in the fields of engineering operations, terrain adaptation, intelligent services, emergency rescue, and so on [1,2,3,4,5,6,7]. Thereupon, batches of well-manufactured, versatile and skilled multi-legged robots have been developed and gradually entered the public’s view, such as Spot [8], ANYmal [9], and Octopus [10]. However, as a series-parallel compound omnidirectional mobile device, unmanned multi-legged robots have a more complex structural layout and control system than the general mobile vehicles [11,12], and its theoretical research and technical system in all aspects are still immature. In particular, research on body trajectory tracking, which is the basis for safe movement and precise operation of unmanned robots [13,14,15], has not received enough attention.

Generally, the feet of an articulated multi-legged robot cannot rotate continuously like wheels, but instead propel the body by means of alternate supports and swings. Therefore, the motion trajectory of the legged robot is piecewise smooth rather than holistic smooth as that of the wheeled robot. This naturally leads to thinking about the segment size of the trajectory, that is, the consideration of the constraints of movement stride (i.e., the body displacement in a rhythmic period). Similar to the fact that humans need to adjust the appropriate step length according to their own motion limits when walking, the stride value sent to legged robots must also conform to the limitations of mechanical structure. However, the existing works have not reported how to connect the body stride with the robot’s control system, which is the primary problem addressed in this paper. In addition, to ensure the safety of motion control of hexapod robots, the quantitative modeling and evaluation of stride constraints should be taken as an important part of designing optimal control strategies, which are rarely mentioned in existing research.

Although some scholars have discussed the tracking control and constraint handling of legged robots, they are all very different from the core content of this paper. The authors of [16,17] deeply studied the smooth trajectory tracking of a single leg of hexapod robots, but this kind of local control methods cannot directly reflect the tracking ability of the robot’s body. Although Chen et al. [18,19] presented the accurate tracking of reference trajectories of the body position and posture of a hexapod robot based on sliding mode control method, it had not yet involved any constraints. In the recent past, model predictive control (MPC) has played an increasingly important role in optimal control and decision-making scenarios that need to consider various virtual or physical constraints [20,21,22,23]. In particular, MPC plays a pivotal role in dealing with torque limitations and friction cone constraints involved in the dynamic motion and gait planning of biped [24,25] and quadruped [26,27] walking robots. However, none of these works discussed the selection and restriction of stride length. When designing the push recovery and self-balancing algorithm for a quadruped robot, Dini et al. [28] used MPC to deal with various physical constraints including stability, self-collision avoidance, and step limitation. However, the modeling and solving process in this paper is essentially different from their approach. In addition, given that MPC has significant advantages in dealing with optimization control problems with multiple constraints [20,29], it has gradually emerged in research and the application of hexapod robots. For example, considering the mechanical restrictions, energy consumption, and physical constraints, Hu et al. [30] proposed a constrained MPC strategy to achieve acyclic and stable walking of a hexapod robot over irregular terrain. However, this work did not involve the constraints related to stride limitations. Aiming at the problem of stable walking and trajectory tracking of a six-wheel-legged robot under complex road conditions and heavy loads, an MPC scheme based on fuzzy approximation was developed in [31], which achieved satisfactory tracking performance under the constraints of velocity limit and input increment restriction. Unfortunately, the modeling analysis and controller design in this article are more like that of a wheeled mobile robot, neither reflecting the characteristics of the rhythmic motion control of legged robots, nor containing more constraints such as the swing limit of the limbs. In view of the weaknesses of current research works on hexapod robot’s trajectory tracking and the advantages of MPC in handling multi-constrained optimization, this paper establishes a trajectory tracking model based on the characteristics of the rhythmic motion of legged robots, and designs a tracking control scheme based on predictive optimization.

Compared with the existing research results on trajectory tracking of multi-legged robots, the main contributions of this paper are as follows:

(1)

A body-level trajectory tracking model based on the motion stride for hexapod robot is pioneeringly established, which quantitatively reveals the characteristics of the rhythmic motion of legged robots;

(2)

The quantitative relationship of stretch and yaw constraints of limb based on body stride and tripod gait is groundbreakingly modeled, which is used as an important consideration in designing the MPC-based optimal controller to ensure the structural integrity of the robot while achieving omnidirectional accurate trajectory tracking;

(3)

A method of defining the reference stride length based on the limb constraints and the integral mean is creatively proposed, which effectively takes into account the physical limitations and movement capabilities of the robot; meanwhile, a solution method of variable stride periods and the corresponding real-time replanning strategy are proposed, which effectively improve the tracking efficiency and real-time tracking ability of the robot.

The rest of this paper is organized as follows. Section 2 introduces the body-level kinematics of the multi-legged robot and the process of establishing the trajectory tracking error model. Section 3 constructs the stride-based model of limb constraints for a hexapod robot moving with tripod gait, and then designs an MPC-based optimal tracking controller. In Section 4, the strategies for reference trajectory segmenting, stride period solving and reference stride replanning are proposed based on the stride constraints. The demonstration results of the simulation analysis and tracking experiment of the robot are shown in Section 5. Section 6 summarizes the results of this work and discusses future research directions.

2. Body-Level Kinematics

The body-level kinematics of hexapod robot is mainly to establish the mathematical relationship between the motion stride and the body pose, which provides a model basis for the body-level trajectory tracking control. In view of the complexity of the structure and control of multi-legged robots and the fact that there are few direct results of the body-level trajectory tracking, this work, analogous to the omnidirectional mobility of a wheeled mobile robot on a plane, only considers the omnidirectional motion and tracking capabilities of the hexapod robot on a two-dimensional plane including the two-dimensional translation along the horizontal direction and the rotation around the plane normal.

Unlike wheeled robots that rely on the continuous rotation and support of wheels to move, legged robots rely on the regular discrete support and coordinated swing of all limbs to drive the whole machine. This reflects that the accuracy of the global pose information of the legged robot in the world coordinate system directly depends on the reliability of the local segmental modeling of periodic rhythmic motion. Specifically, the pose of a multi-legged robot in the world coordinate system is related not only to its gait mode and motion state in the current stride period (that is, the time duration in which all limbs of the robot complete one swing), but also to the pose information at the end of the previous stride period. Therefore, only by accurately establishing the kinematic relationship and tracking model of the robot in a stride period, and then inheriting and analogizing this relationship and model, can the robot achieve satisfactory desired movement and tracking. Based on the above ideas, the specific quantitative analysis and modeling are carried out below.

The diagram of the kinematic modeling of a multi-legged robot in the i-th ($i\in {\mathbb{N}}_{+}$) stride period is shown in Figure 1. The body coordinate system at the beginning of this stride period, ${\Sigma }_{i-1}^{\mathrm{B}}\triangleq$$\left\{O_{0,i-1}{x}_{0,i-1}{y}_{0,i-1}{z}_{0,i-1}\right\}$, is regarded as the local reference coordinate system, and the body coordinate system at the end of this stride period is denoted as ${\Sigma }_i^{\mathrm{B}}\triangleq \left\{O_{0,i}{x}_{0,i}{y}_{0,i}{z}_{0,i}\right\}$. Assuming that in the current stride period, the total stride length of the body movement is $S_l$, the yaw angle of point $O_{0,i}$ away from $x_{0,i-1}$ axis is $\Psi$, and the total rotation angle of the body around $z_{0,i-1}$ axis is $S_z$. Then, at any moment in the current stride period, the instantaneous relative pose of the body center with respect to (w.r.t.) the frame ${\Sigma }_{i-1}^{\mathrm{B}}$ is

$$\left\{\begin{array}{c}{x}_{\mathrm{B}}={\Gamma }_{\mathrm{B}}·S_l·cos\Psi ,\hfill \\ y_{\mathrm{B}}={\Gamma }_{\mathrm{B}}·S_l·sin\Psi ,\hfill \\ {\theta }_{\mathrm{B}}={\Gamma }_{\mathrm{B}}·S_z.\hfill \end{array}\right.$$

(1)

In the above equation, ${\Gamma }_{\mathrm{B}}=g(t,\delta ,{\Gamma }_{\mathrm{F}})$ is an interpolation function of time t related to the duty cycle $\delta$ and the foot swing interpolation ${\Gamma }_{\mathrm{F}}$. In addition, ${\Gamma }_{\mathrm{F}}$ is a linear or nonlinear function customized according to actual requirements. A typical case of it can be found in [32]. Without causing ambiguity, the character t representing time is omitted here and in the following texts to simplify writing and expression.

In the i-th stride period, if it is known that the pose of the body center at the initial time relative to the world coordinate system ${\Sigma }^{\mathrm{W}}\triangleq \left\{O_{\mathrm{w}}{X}_{\mathrm{w}}{Y}_{\mathrm{w}}{Z}_{\mathrm{w}}\right\}$ is $[X_{i-1},$$Y_{i-1},$${\Theta }_{i-1}{]}^{\mathrm{T}}$. Here ($X_{i-1},$$Y_{i-1}$) is the absolute position of the body center in the world frame ${\Sigma }^{\mathrm{W}}$; and ${\Theta }_{i-1}$ is the absolute orientation in the world frame ${\Sigma }^{\mathrm{W}}$, that is, the angle between the $x_{0,i-1}$ axis and the $X_{\mathrm{w}}$ axis. Then, at any moment in the current stride period, the instantaneous absolute pose of the center of body relative to the world frame ${\Sigma }^{\mathrm{W}}$ is

$$\left\{\begin{array}{c}X=x_{\mathrm{B}}cos{\Theta }_{i-1}-y_{\mathrm{B}}sin{\Theta }_{i-1}+X_{i-1},\hfill \\ Y=x_{\mathrm{B}}sin{\Theta }_{i-1}+y_{\mathrm{B}}cos{\Theta }_{i-1}+Y_{i-1},\hfill \\ \Theta ={\theta }_{\mathrm{B}}+{\Theta }_{i-1}.\hfill \end{array}\right.$$

(2)

Thus, in a stride period, the velocity kinematics model of the body relative to the world frame ${\Sigma }^{\mathrm{W}}$ is as follows:

$$\left\{\begin{array}{c}\dot{X}={\dot{\Gamma }}_{\mathrm{B}}{S}_{l}cos(\Theta -{\Gamma }_{\mathrm{B}}{S}_z+\Psi )\triangleq f_1(\mathit{z},\mathit{u}),\hfill \\ \dot{Y}={\dot{\Gamma }}_{\mathrm{B}}{S}_{l}sin(\Theta -{\Gamma }_{\mathrm{B}}{S}_z+\Psi )\triangleq f_2(\mathit{z},\mathit{u}),\hfill \\ \dot{\Theta }={\dot{\Gamma }}_{\mathrm{B}}{S}_z\triangleq f_3(\mathit{z},\mathit{u}),\hfill \end{array}\right.$$

(3)

where $\mathit{z}\triangleq$${[X,Y,\Theta ]}^{\mathrm{T}}$ is the state variable of the above-mentioned nonlinear control system, and $\mathit{u}\triangleq$${[S_l,\Psi ,S_z]}^{\mathrm{T}}$ is the control input composed of the stride triple. Define $\mathit{F}(\mathit{z},\mathit{u})\triangleq$$[f_1(\mathit{z},\mathit{u}),$$f_2(\mathit{z},\mathit{u}),$$f_3{(\mathit{z},\mathit{u})]}^{\mathrm{T}}$, then the nominal system model (3) can be abbreviated as $\dot{\mathit{z}}=\mathit{F}(\mathit{z},\mathit{u})$.

Taking the Taylor expansion of Equation (3) at the reference point $({\mathit{z}}^{\mathrm{r}},{\mathit{u}}^{\mathrm{r}})$ (where ${\mathit{z}}^{\mathrm{r}}=$${[X^{\mathrm{r}},Y^{\mathrm{r}},{\Theta }^{\mathrm{r}}]}^{\mathrm{T}}$, ${\mathit{u}}^{\mathrm{r}}=$${[S_l^{\mathrm{r}},{\Psi }^{\mathrm{r}},S_z^{\mathrm{r}}]}^{\mathrm{T}}$), i.e., the point on the reference trajectory, and ignoring the higher-order terms, the body-level trajectory tracking error model can be obtained in the following form:

$$\dot{\mathbf{\xi }}=\mathit{A}\mathbf{\xi }+\mathit{B}\mathbf{\upsilon },$$

(4)

where

$$\begin{array}{c}\begin{array}{cc}\hfill \mathbf{\xi }& \triangleq \mathit{z}-{\mathit{z}}^{\mathrm{r}}=\left[\begin{array}{c}X-X^{\mathrm{r}}\\ Y-Y^{\mathrm{r}}\\ \Theta -{\Theta }^{\mathrm{r}}\end{array}\right],\mathbf{\upsilon }\triangleq \mathit{u}-{\mathit{u}}^{\mathrm{r}}=\left[\begin{array}{c}{S}_l-S_l^{\mathrm{r}}\\ \Psi -{\Psi }^{\mathrm{r}}\\ S_z-S_z^{\mathrm{r}}\end{array}\right],\hfill \\ \hfill \mathit{A}& \triangleq \frac{\partial \mathit{F}}{\partial \mathit{z}}{|}_{\begin{array}{c}\mathit{z}={\mathit{z}}^{\mathrm{r}}\\ \mathit{u}={\mathit{u}}^{\mathrm{r}}\\ \phantom{{\mathit{u}}^{\mathrm{r}}}\end{array}}=\left[\begin{array}{ccc}0& 0& -sin\left(\varpi \right){\dot{\Gamma }}_{\mathrm{B}}{S}_l^{\mathrm{r}}\\ 0& 0& cos\left(\varpi \right){\dot{\Gamma }}_{\mathrm{B}}{S}_l^{\mathrm{r}}\\ 0& 0& 0\end{array}\right],\hfill \end{array}\end{array}$$

$$\begin{array}{c}\begin{array}{cc}\hfill \mathit{B}& \triangleq \frac{\partial \mathit{F}}{\partial \mathit{u}}{|}_{\begin{array}{c}\mathit{z}={\mathit{z}}^{\mathrm{r}}\\ \mathit{u}={\mathit{u}}^{\mathrm{r}}\\ \phantom{{\mathit{u}}^{\mathrm{r}}}\end{array}}=\left[\begin{array}{ccc}cos\left(\varpi \right){\dot{\Gamma }}_{\mathrm{B}}& -sin\left(\varpi \right){\dot{\Gamma }}_{\mathrm{B}}{S}_l^{\mathrm{r}}& sin\left(\varpi \right){\dot{\Gamma }}_{\mathrm{B}}{S}_l^{\mathrm{r}}{\Gamma }_{\mathrm{B}}\\ sin\left(\varpi \right){\dot{\Gamma }}_{\mathrm{B}}& cos\left(\varpi \right){\dot{\Gamma }}_{\mathrm{B}}{S}_l^{\mathrm{r}}& -cos\left(\varpi \right){\dot{\Gamma }}_{\mathrm{B}}{S}_l^{\mathrm{r}}{\Gamma }_{\mathrm{B}}\\ 0& 0& {\dot{\Gamma }}_{\mathrm{B}}\end{array}\right].\hfill \end{array}\end{array}$$

with $\varpi \triangleq {\Theta }^{\mathrm{r}}-{\Gamma }_{\mathrm{B}}{S}_z^{\mathrm{r}}+{\Psi }^{\mathrm{r}}$.

3. Constraints and Predictive Controller

It is well known that compared with periodic rhythm gaits such as quadruped and wave, the tripod gait of hexapod robots is the most efficient and, therefore, the most common gait. However, few works have explored the limb constraints of hexapod robots moving with rhythmic gait. Taking the tripod gait as a breakthrough, this paper focuses on the quantitative modeling method of the limb’s stretch and yaw constraints under this gait, and the predictive control strategy is designed to comply with these constraints. The proposed modeling idea may be able to provide theoretical guidance for further research on constrained tracking control of hexapod robots with rhythmic gaits or free gaits.

3.1. Stretch and Yaw Constraints of Limb

When the hexapod robot moves rhythmically with tripod gait, it is agreed that its body follows the motion rule of first translation and then rotation at every moment. Then, in a stride period, the robot has three key configurations—the starting configuration, the semi-periodic configuration, and the ending configuration. Their corresponding body coordinate systems are denoted as $\left\{O_{\mathrm{s}}\right\}$, $\left\{O_{\mathrm{m}}\right\}$, and $\left\{O_{\mathrm{e}}\right\}$, respectively, as shown in Figure 2. Assume that in each stride period, the robot is in the nominal pose (where all joint angles are 0) at the starting configuration and exactly recovers to the nominal pose at the ending configuration. Then, for a hexapod robot that moves periodically with a tripod gait, it is obvious that the amplitudes of stretch and yaw of each limb are the largest when the robot is in the semi-periodic configuration. Therefore, as long as the robot does not exceed its structural limits in this configuration, the movement at any time within a stride period is feasible and safe. Naturally, the following focuses on the quantitative analysis of the robot in the semi-periodic configuration.

Under the assumption of the above walking rules, when the robot is in the starting and ending configurations, the homogeneous position coordinates of the foot relative to the coordinate systems $\left\{O_{\mathrm{s}}\right\}$ and $\left\{O_{\mathrm{e}}\right\}$ are the same:

$${}^{O_{\mathrm{s}}}{\mathit{p}}_0={}^{O_{\mathrm{e}}}{\mathit{p}}_0=\left[\begin{array}{c}{p}_{x,0}\\ p_{y,0}\\ p_{z,0}\\ 1\end{array}\right],$$

(5)

where ${[p_{x,0},p_{y,0},p_{z,0}]}^{\mathrm{T}}$ represents the initial position coordinate of the foot relative to the body frame.

Then, in the semi-periodic configuration, the homogeneous position coordinate of the stance phase relative to the frame $\left\{O_{\mathrm{m}}\right\}$ is

$$\begin{array}{c}\hfill {}^{O_{\mathrm{m}}}{\mathit{p}}_{\mathrm{st}}={}_{O_{\mathrm{s}}}^{O_{\mathrm{m}}}\mathit{T}·{}^{O_{\mathrm{s}}}{\mathit{p}}_0={\left[\mathrm{Trans}(\frac{S_l}{2}cos\Psi ,\frac{S_l}{2}sin\Psi ,0)·\mathrm{Rot}(z,\frac{S_z}{2})\right]}^{-1}·{}^{O_{\mathrm{s}}}{\mathit{p}}_0,\end{array}$$

(6)

where $\mathit{T}$ represents the homogeneous transformation matrix synthesized by the homogeneous translation operator $\mathrm{Trans}(·)$ and the rotation operator $\mathrm{Rot}(·)$.

Meanwhile, the homogeneous position coordinate of the swing phase relative to the frame $\left\{O_{\mathrm{m}}\right\}$ is

$$\begin{array}{cc}\hfill {}^{O_{\mathrm{m}}}{\mathit{p}}_{\mathrm{sw}}& ={}_{O_{\mathrm{s}}}^{O_{\mathrm{m}}}\mathit{T}·{}_{O_{\mathrm{e}}}^{O_{\mathrm{s}}}\mathit{T}·{}^{O_{\mathrm{e}}}{\mathit{p}}_0\hfill \\ \hfill & ={\left[\mathrm{Trans}(\frac{S_l}{2}cos\Psi ,\frac{S_l}{2}sin\Psi ,0)·\mathrm{Rot}(z,\frac{S_z}{2})\right]}^{-1}·\left[\mathrm{Trans}(S_{l}cos\Psi ,S_{l}sin\Psi ,0)·\mathrm{Rot}(z,S_z)\right]·{}^{O_{\mathrm{e}}}{\mathit{p}}_0.\hfill \end{array}$$

(7)

Further, the homogeneous position coordinate of each foot relative to the frame $\left\{O_{\mathrm{m}}\right\}$ can be obtained by unifying Equations (6) and (7) as

$${}^{O_{\mathrm{m}}}{\mathit{p}}_{\mathrm{foot}}\triangleq \left[\begin{array}{c}{p}_{\mathrm{f},x}\\ p_{\mathrm{f},y}\\ p_{\mathrm{f},z}\\ 1\end{array}\right]=\left[\begin{array}{c}{p}_{x,0}cos\frac{S_z}{2}-\lambda ·p_{y,0}sin\frac{S_z}{2}+\frac{S_l}{2}cos(\Psi -\frac{S_z}{2})\\ \lambda ·p_{x,0}sin\frac{S_z}{2}+p_{y,0}cos\frac{S_z}{2}+\frac{S_l}{2}sin(\Psi -\frac{S_z}{2})\\ p_{z,0}\\ 1\end{array}\right],$$

(8)

where $\lambda =\left\{\begin{array}{cc}1,\hfill & \mathrm{for}\mathrm{swing}\mathrm{phase};\hfill \\ -1,\hfill & \mathrm{for}\mathrm{stance}\mathrm{phase}.\hfill \end{array}\right.$

Moreover, in the semi-periodic configuration, the homogeneous position coordinate of the hip joint w.r.t. the frame $\left\{O_{\mathrm{m}}\right\}$ is ${}^{O_{\mathrm{m}}}{\mathit{p}}_{\mathrm{hip}}={[p_{\mathrm{h},x},p_{\mathrm{h},y},0,1]}^{\mathrm{T}}$, which is usually only related to the structure dimensions of the body. Then, the horizontal coordinates from the foot end to the hip joint are as follows:

$$\left\{\begin{array}{cc}\hfill {}^{O_{\mathrm{m}}}{p}_{\mathrm{f}2\mathrm{h},x}& =p_{\mathrm{f},x}-p_{\mathrm{h},x}\hfill \\ \hfill & =p_{x,0}cos\frac{S_z}{2}-\lambda ·p_{y,0}sin\frac{S_z}{2}+\frac{S_l}{2}cos(\Psi -\frac{S_z}{2})-p_{\mathrm{h},x},\hfill \\ \hfill {}^{O_{\mathrm{m}}}{p}_{\mathrm{f}2\mathrm{h},y}& =p_{\mathrm{f},y}-p_{\mathrm{h},y}\hfill \\ \hfill & =\lambda ·p_{x,0}sin\frac{S_z}{2}+p_{y,0}cos\frac{S_z}{2}+\frac{S_l}{2}sin(\Psi -\frac{S_z}{2})-p_{\mathrm{h},y}.\hfill \end{array}\right.$$

(9)

Thus, the horizontal stretch length of a limb described by the motion stride is $L_s=\sqrt{{}^{O_{\mathrm{m}}}{p}_{\mathrm{f}2\mathrm{h},x}^2+{}^{O_{\mathrm{m}}}{p}_{\mathrm{f}2\mathrm{h},y}^2}$.

Following the above analysis, without considering the fluctuation of the robot’s body, the horizontal stretch length of a limb described by the stride should not exceed the structural constraints of the mechanism itself, that is,

$$0\le L_s\le S_{\mathrm{hor},\max }.$$

(10)

where $S_{\mathrm{hor},\max }$ represents the maximum horizontal stretch length from the foot end to the hip joint, which is an inherent quantity determined by the mechanical size and the joint limits.

In addition, the yaw angle solved by the motion stride is also subject to the inherent limitations of the robot mechanism, i.e.,

$$q_1=\mathrm{Atan}2({}^{O_{\mathrm{m}}}{p}_{\mathrm{f}2\mathrm{h},y},{}^{O_{\mathrm{m}}}{p}_{\mathrm{f}2\mathrm{h},x})-{\alpha }_j\in \left[q_{1,\min },q_{1,\max }\right],$$

(11)

where $\mathrm{Atan}2(·)$ is the four-quadrant inverse tangent function; $q_{1,\min }$ and $q_{1,\max }$ indicate the rotational range of Joint1; and ${\alpha }_j$ represents the azimuth angle of the j-th $(j=1,2,\cdots ,6)$ limb w.r.t. the body system.

Obviously, the constraints (10) and (11) are both functions of the stride parameter $(S_l,\Psi ,S_z)$. To ensure the structural integrity of the robot, the designed controller must satisfy the above constraints.

3.2. Model Predictive Controller

In view of the unique advantages of MPC in dealing with multi-constraint optimization problems, this section uses it to design an optimal control model with equality constraint, input saturation constraints and nonlinear inequality constraints. The detailed design process is described below.

Firstly, discretize the body-level tracking error model of the nominal linearization (4) to obtain

$${\mathbf{\xi }}_{k+1}={\mathit{G}}_k{\mathbf{\xi }}_k+{\mathit{H}}_k{\mathbf{\upsilon }}_k,$$

(12)

where ${\mathbf{\xi }}_k$ and ${\mathbf{\upsilon }}_k$ are the pose error and stride error at time k, respectively; ${\mathit{G}}_k\triangleq \mathit{I}+t_s{\mathit{A}}_k$, ${\mathit{H}}_k\triangleq t_s{\mathit{B}}_k$; $\mathit{I}$ represents an identity matrix with appropriate dimensions; $t_s$ is the sampling step size; ${\mathit{A}}_k,{\mathit{G}}_k\in {\mathbb{R}}^{N_x\times N_x}$; ${\mathit{B}}_k,{\mathit{H}}_k\in {\mathbb{R}}^{N_x\times N_u}$; $N_x$ and $N_u$ respectively represent the dimensions of the state variable and the control variable.

Then, define $\Delta {\mathbf{\upsilon }}_k\triangleq {\mathbf{\upsilon }}_k-{\mathbf{\upsilon }}_{k-1}$, and

$$\begin{array}{c}{\tilde{\mathbf{\xi }}}_k\triangleq \left[\begin{array}{c}{\mathbf{\xi }}_{k+1|k}\\ {\mathbf{\xi }}_{k+2|k}\\ ⋮\\ {\mathbf{\xi }}_{k+N_{\mathrm{p}}|k}\end{array}\right],{\tilde{\mathbf{\upsilon }}}_k\triangleq \left[\begin{array}{c}{\mathbf{\upsilon }}_{k|k}\\ {\mathbf{\upsilon }}_{k+1|k}\\ ⋮\\ {\mathbf{\upsilon }}_{k+N_{\mathrm{c}}-1|k}\end{array}\right],\Delta {\tilde{\mathbf{\upsilon }}}_k\triangleq \left[\begin{array}{c}\Delta {\mathbf{\upsilon }}_{k|k}\\ \Delta {\mathbf{\upsilon }}_{k+1|k}\\ ⋮\\ \Delta {\mathbf{\upsilon }}_{k+N_{\mathrm{c}}-1|k}\end{array}\right].\end{array}$$

In the above formulas, ${\mathbf{\xi }}_{k+m|k}$, ${\mathbf{\upsilon }}_{k+n|k}$ and $\Delta {\mathbf{\upsilon }}_{k+n|k}$ respectively represent the state, control variable and control increment at the time $k+m$ and $k+n$ predicted at time k ($m=1,2,\cdots ,N_{\mathrm{p}}$ and $n=0,1,\cdots ,N_{\mathrm{c}}-1$); $N_{\mathrm{p}}$ and $N_{\mathrm{c}}(\le N_{\mathrm{p}})$ represent the preview horizon and control horizon, respectively.

Suppose that for any $m\in \{1,2,\cdots ,N_{\mathrm{p}}\}$, there are ${\mathit{G}}_{k+m|k}={\mathit{G}}_k$ and ${\mathit{H}}_{k+m|k}={\mathit{H}}_k$. Then, the prediction model can be written in the following lumped form:

$${\tilde{\mathbf{\xi }}}_k={\tilde{\mathit{G}}}_k{\mathbf{\xi }}_k+{\tilde{\mathit{h}}}_k{\mathbf{\upsilon }}_{k-1}+{\tilde{\mathit{H}}}_k\Delta {\tilde{\mathbf{\upsilon }}}_k,$$

(13)

where both ${\tilde{\mathit{G}}}_k$, ${\tilde{\mathit{H}}}_k$, and ${\tilde{\mathit{h}}}_k$ are block matrices, and the block elements of these matrices are defined as follows:

$$\begin{array}{cc}\hfill {\tilde{\mathit{G}}}_k(m;1)& ={\mathit{G}}_k^m,\hfill \\ \hfill {\tilde{\mathit{H}}}_k(m;n)& =\left\{\begin{array}{cc}\mathbf{0}\hfill & ,\mathrm{if}n>m\hfill \\ {\displaystyle \sum _{\ell =n}^m}{\mathit{G}}_k^{m-\ell }{\mathit{H}}_k\hfill & ,\mathrm{if}n\le m\hfill \end{array}\right.\hfill \\ \hfill {\tilde{\mathit{h}}}_k(m;1)& ={\tilde{\mathit{H}}}_k(m;1),\hfill \\ \hfill m=1,2,& \cdots ,N_{\mathrm{p}};n=1,2,\cdots ,N_{\mathrm{c}}.\hfill \end{array}$$

Next, we define an optimization objective function as

$$\underset{\Delta {\tilde{\mathbf{\upsilon }}}_k}{min}{\mathit{J}}_k\triangleq {\tilde{\mathbf{\xi }}}_k^{\mathrm{T}}\tilde{\mathit{Q}}{\tilde{\mathbf{\xi }}}_k+\Delta {\tilde{\mathbf{\upsilon }}}_k^{\mathrm{T}}\tilde{\mathit{R}}\Delta {\tilde{\mathbf{\upsilon }}}_k,$$

(14)

where $\tilde{\mathit{Q}}={\mathit{I}}_{N_p}\otimes \mathit{Q}$ and $\tilde{\mathit{R}}={\mathit{I}}_{N_c}\otimes \mathit{R}$, with $\mathit{Q}\in {\mathbb{R}}^{N_x\times N_x}$ and $\mathit{R}\in {\mathbb{R}}^{N_u\times N_u}$ are weight matrices corresponding to the predictive state vector ${\tilde{\mathbf{\xi }}}_k$ and the predictive control increment $\Delta {\tilde{\mathbf{\upsilon }}}_k$, respectively; ⊗ represents the Kronecker product.

The corresponding constraints include at least the following five categories:

$$\mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}•{\tilde{\mathbf{\xi }}}_k={\tilde{\mathit{G}}}_k{\mathbf{\xi }}_k+{\tilde{\mathit{h}}}_k{\mathbf{\upsilon }}_{k-1}+{\tilde{\mathit{H}}}_k\Delta {\tilde{\mathbf{\upsilon }}}_k;\hfill \\ •{\tilde{\mathbf{\upsilon }}}_{\min }\le {\tilde{\mathbf{\upsilon }}}_k(=\tilde{\mathit{I}}·\Delta {\tilde{\mathbf{\upsilon }}}_k+\tilde{\mathbf{1}}·{\mathbf{\upsilon }}_{k-1})\le {\tilde{\mathbf{\upsilon }}}_{\max };\hfill \\ •\Delta {\tilde{\mathbf{\upsilon }}}_{\min }\le \Delta {\tilde{\mathbf{\upsilon }}}_k\le \Delta {\tilde{\mathbf{\upsilon }}}_{\max };\hfill \\ •L_{s,j}\left(\mathit{u}\right)\in [0,S_{\mathrm{hor},\max }],(j=1,2,\cdots ,6);\hfill \\ •q_{1,j}\left(\mathit{u}\right)\in [q_{1,\min },q_{1,\max }],(j=1,2,\cdots ,6);\hfill \end{array}\right.$$

(15)

where $L_{s,j}\left(\mathit{u}\right)$ and $q_{1,j}\left(\mathit{u}\right)$ are the stretch length and yaw angle of the j-th $(j=1,2,\cdots ,6)$ limb defined by (10) and (11), respectively, and they are both nonlinear functions about the stride input $\mathit{u}$; and

$$\begin{array}{c}\hfill \tilde{\mathit{I}}\triangleq {\left[\begin{array}{cccc}1& 0& \cdots & 0\\ 1& 1& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 1& 1& \cdots & 1\end{array}\right]}_{N_{\mathrm{c}}\times N_{\mathrm{c}}}\otimes {\mathit{I}}_{N_u},\tilde{\mathbf{1}}\triangleq {\left[\begin{array}{c}1\\ 1\\ ⋮\\ 1\end{array}\right]}_{N_{\mathrm{c}}\times 1}\otimes {\mathit{I}}_{N_u},\end{array}$$

$$\begin{array}{c}\hfill \mathit{u}\triangleq \left[\begin{array}{c}{S}_l\\ \Psi \\ S_z\end{array}\right]={\left[\begin{array}{cccccc}1& 0& 0& 0& \cdots & 0\\ 0& 1& 0& 0& \cdots & 0\\ 0& 0& 1& 0& \cdots & 0\end{array}\right]}_{N_u\times (N_u\times N_{\mathrm{c}})}·(\tilde{\mathit{I}}·\Delta {\tilde{\mathbf{\upsilon }}}_k+\tilde{\mathbf{1}}·{\mathbf{\upsilon }}_{k-1})+{\mathit{u}}^{\mathrm{r}}.\end{array}$$

In Equation (15), the first type of constraints corresponds to the predictive model constraint (13), which relates the predictive state vector to the predictive control increment to be solved. The second type of constraints is the boundary restriction of the predictive control vector, which indirectly reflects the limits of the body’s motion stride. The third type of constraints is the boundary limit of the predictive control increment to be solved, preventing the stride from changing too fast. The fourth and fifth types of constraints represent the stretch length restrictions and yaw constraints of all limbs, respectively, which are used to protect the physical structure of the robot itself.

It should be noted that the modeling process of the above-mentioned optimal controller (14)–(15) follows the design idea of MPC. Consequently, the proof of its stability is universal and trivial, which is described in detail in [33,34,35], and will not be repeated here.

At time k, we use the sequential quadratic programming (SQP) [36,37,38] to solve the above-mentioned optimization problem (14) with nonlinear constraints (15) to get the optimal control increment $\Delta {\mathbf{\upsilon }}_k^{*}$ (which is the first $N_u$ rows of the column vector $\Delta {\tilde{\mathbf{\upsilon }}}_k^{*}$), and then the optimal control input error ${\mathbf{\upsilon }}_k^{*}={\mathbf{\upsilon }}_{k-1}+\Delta {\mathbf{\upsilon }}_k^{*}$ can be obtained. At the next moment, by repeating the above process, a series of optimal control variables satisfying physical constraints can be obtained.

4. Stride Period and Reference Stride

The above has solved the problem of obtaining the optimal stride error. If the stride period and the reference stride are also known, the robot can be controlled to move at the corresponding speed and stride. A simple and straightforward idea is to take a predetermined constant stride period and constant reference stride. Unfortunately, in diverse scenarios, these constant reference values can seriously affect the tracking efficiency of the robot. For example, when the initial pose is significantly different from the reference starting point, the robot can only slowly approach the reference trajectory or suddenly send a completely unreasonable control signal that exceeds the motor’s capacity under the limitation of the constant stride period and the constant reference stride. Even if the robot has reached the reference trajectory at a certain time, it will appear tortuous and back-and-forth fluctuations along the reference trajectory due to the guidance of the constant reference values. To this end, this paper proposes a method to determine the variable stride period based on the reference stride length, and then proposes a strategy to solve the reference stride based on real-time replanning. The obtained stride period and reference stride depend on the reference trajectory and the actual pose information of the robot, which can effectively improve the robot’s tracking performance.

4.1. Determination of Stride Period

An “appropriate” reference stride length should take into account both the robot’s movement ability and locomotion efficiency. In other words, it cannot be too large, which easily exceeds the structural limitations of the robot; while it also cannot be too small, which easily causes the robot to move slowly and inefficiently. To this end, this paper proposes a method for determining the length of reference stride based on the stride constraints and integral mean.

Consider again the stretch and yaw constraints of limbs described by the motion stride of the body. Under these constraints, the stride length $S_l$ can be regarded as a function of the yaw angle $\Psi$ and the orientation angle $S_z$. Then, we take the integral mean of the spacial surface where the effective maximum stride length $S_{l,\max }$ is located (shown in Figure 3) as the reference stride length, that is,

$${\overline{S}}_{l,\max }=\frac{{\int \int }_D{S}_{l,\max }(\Psi ,S_z)\mathrm{d}\Psi \mathrm{d}{S}_z}{\left|D\right|},$$

(16)

where $\left|D\right|$ represents the projected area of the surface $S_{l,\max }(\Psi ,S_z)$ onto the plane $\Psi O{S}_z$. The practical significance of the reference stride length ${\overline{S}}_{l,\max }$ defined by Equation (16) is to give full play to the robot’s motion ability while taking into account the stride constraints.

Next, as shown in Figure 4, on the reference trajectory, starting from the start point (denoted as $P_{i-1}(t_{i-1},$$X_{i-1}^{\mathrm{r}},$$Y_{i-1}^{\mathrm{r}},$${\Theta }_{i-1}^{\mathrm{r}})$ ), we search for a trajectory point $P_i(t_i,$$X_i^{\mathrm{r}},$$Y_i^{\mathrm{r}},$${\Theta }_i^{\mathrm{r}})$ along the timeline such that the following equation holds:

$$\sqrt{{(X_i^{\mathrm{r}}-X_{i-1}^{\mathrm{r}})}^2+{(Y_i^{\mathrm{r}}-Y_{i-1}^{\mathrm{r}})}^2}={\overline{S}}_{l,\max }.$$

(17)

Then, the duration of the i-th stride period is

$$T_{\mathrm{sp},i}=t_i-t_{i-1}(i=1,2,\cdots ).$$

(18)

Then, taking $P_i$ as the start point of the search, we repeat the above process until the end of the trajectory. It should be noted that in the last search range, if the Euclidean distance from the track point $P_i$ to the end point $P_{\mathrm{end}}$ of the reference trajectory is less than the reference stride length ${\overline{S}}_{l,\max }$, then the duration corresponding to the last stride period is recorded as: $T_{\mathrm{sp},\mathrm{end}}=t_{\mathrm{end}}-t_i$. As a result, a series of key segmentation points, $\{P_0,P_1,\cdots ,P_{\mathrm{end}}\}$, on the reference trajectory and a variable stride period sequence, $\{T_{\mathrm{sp},1},T_{\mathrm{sp},2},\cdots ,T_{\mathrm{sp},\mathrm{end}}\}$, can be obtained.

4.2. Replanning of Reference Stride

As mentioned above, the tracking effect of a legged robot under the action of a certain constant reference stride ${\mathit{u}}^{\mathrm{r}}={[S_l^{\mathrm{r}},{\Psi }^{\mathrm{r}},S_z^{\mathrm{r}}]}^{\mathrm{T}}$ is not optimistic. For a time-varying reference trajectory, a slightly straightforward scheme is to determine the reference stride according to the key segmentation points on the trajectory, that is, the reference stride in the i-th stride period can be set as:

$$\left\{\begin{array}{c}{S}_{l,i}^{\mathrm{r}}=\sqrt{{(X_i^{\mathrm{r}}-X_{i-1}^{\mathrm{r}})}^2+{(Y_i^{\mathrm{r}}-Y_{i-1}^{\mathrm{r}})}^2},\hfill \\ {\Psi }_i^{\mathrm{r}}=\mathrm{Atan}2(Y_i^{\mathrm{r}}-Y_{i-1}^{\mathrm{r}},X_i^{\mathrm{r}}-X_{i-1}^{\mathrm{r}})-{\Theta }_{i-1}^{\mathrm{r}},\hfill \\ S_{z,i}^{\mathrm{r}}={\Theta }_i^{\mathrm{r}}-{\Theta }_{i-1}^{\mathrm{r}}.\hfill \end{array}\right.$$

(19)

However, experiments show that, when the actual pose of the robot is not much different from the initial reference pose, the tracking effect is good; but when the actual pose deviates greatly from the reference value, the tracking deviation under the guidance of the above reference stride is relatively large (which is visually analyzed in Section 5.2). To solve this problem, the actual pose of the robot is taken into account for real-time replanning and correction of the reference yaw angle and rotational orientation. That is, the improved reference stride is determined by

$$\left\{\begin{array}{c}{S}_{l,i}^{\mathrm{p}}={\overline{S}}_{l,\max },\hfill \\ {\Psi }_i^{\mathrm{p}}=\mathrm{Atan}2(Y_i^{\mathrm{r}}-Y_{i-1},X_i^{\mathrm{r}}-X_{i-1})-{\Theta }_{i-1},\hfill \\ S_{z,i}^{\mathrm{p}}={\Theta }_i^{\mathrm{r}}-{\Theta }_{i-1},\hfill \end{array}\right.$$

(20)

where ${\mathit{u}}_i^{\mathrm{p}}\triangleq {[S_{l,i}^{\mathrm{p}},{\Psi }_i^{\mathrm{p}},S_{z,i}^{\mathrm{p}}]}^{\mathrm{T}}$ represents the replanning reference stride in the i-th stride period. A graphical illustration of these symbols is shown in Figure 4.

In specific applications, the actual motion stride is solved according to the formula $\mathit{u}={\mathit{u}}^{\mathrm{p}}+{\mathbf{\upsilon }}^{*}$, which is then converted into the relative position of each foot by the gait mode generator [32]. Then, the rotation angle of each joint can be obtained by the inverse kinematics solver, and the robot is urged to perform desired movements through the drive/execution modules. Finally, the information such as actual pose of the robot is returned through sensors, so that the entire system forms a closed loop, which can significantly improve the accuracy of trajectory tracking. To sum up, the closed-loop structure of the trajectory tracking control of the hexapod robot body designed in this paper can be summarized by Figure 5.

5. Simulations and Demonstrations

5.1. Setups of Robot and Reference Trajectory

In this section, the effectiveness and practicability of the proposed predictive control strategy are illustrated by the experimental case on a walking-climbing hexapod robot WelCH, which is used to omnidirectionally track a composite reference trajectory.

The robot WelCH is a radial symmetrical hexapod robot used for inspection on glass curtain walls, and its six limbs have exactly the same structural design. Each limb contains four active rotating joints, of which the first three drive joints control the spatial position of the foot end, and the fourth drive joint mainly controls the pitch angle of the suction cup at the foot end. The structural sizes of this robot and the rotation range of each active joint are listed in

Table 1. Specifications of the robot WelCH.

Item	Symbol	Unit	Specification
Body radius	R	m	0.18
Length of Link1	$L_1$	m	0.09
Length of Link2	$L_2$	m	0.15
Length of Link3	$L_3$	m	0.16
Length of Link4	$L_4$	m	0.15
Range of Joint1	$q_1$	rad	$\left[-\pi /3,\pi /3\right]$
Range of Joint2	$q_2$	rad	$\left[-\pi /2,\pi /2\right]$
Range of Joint3	$q_3$	rad	$\left[-\pi /4,4\pi /9\right]$
Range of Joint4	$q_4$	rad	$\left[-\pi /2,\pi /2\right]$

. In addition, the other setups of hardware and software of the robot WelCH have been described in detail in [39,40,41] and will not be repeated here. The settings of the robot’s body coordinate system $\left\{O_0{x}_0{y}_0{z}_0\right\}$ and the j-th foot coordinate system $\left\{O_{\mathrm{f}j}{x}_{\mathrm{f}j}{y}_{\mathrm{f}j}{z}_{\mathrm{f}j}\right\}$ are shown in Figure 6. The specific expression of ${\Gamma }_{\mathrm{B}}$ in Equation (1) is consistent with that in [32].

As mentioned earlier, this paper focuses on the ability of multi-legged robots to track planar trajectories. To this end, it is set that the center of body of WelCH is always kept at the initial height, $H_{\mathrm{ver}}\equiv \left|p_{z,0}\right|=L_3+L_4$, when moving. Thus, according to the geometric diagram (shown in Figure 7) and simple algebraic calculation, it can be known that the maximum length of horizontal stretch from the foot end to the hip joint is $S_{\mathrm{hor},\max }=L_1+\sqrt{L_2^2+2L_2{L}_3{(sin{q}_3)}_{\max }}\approx 0.35$, which is jointly determined by the robot’s link length $L_i$ and the joint angle $q_i$.

To test the omnidirectional tracking ability of the multi-legged robot, the selected reference trajectory is composed of multi-segment time-varying curves, which are defined as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}0\le t\le 30\left[\mathrm{s}\right]\hfill \\ X^{\mathrm{r}}=0.15t\left[\mathrm{m}\right]\hfill \\ Y^{\mathrm{r}}=1.5cos\left(\frac{\pi }{15}t\right)\left[\mathrm{m}\right]\hfill \\ {\Theta }^{\mathrm{r}}=\mathrm{Atan}2({\dot{Y}}^{\mathrm{r}},{\dot{X}}^{\mathrm{r}})\left[\mathrm{rad}\right]\hfill \end{array}\right.\left\{\begin{array}{c}30<t\le 40\left[\mathrm{s}\right]\hfill \\ X^{\mathrm{r}}=X_{t=30}^{\mathrm{r}}+0.30(t-30)\left[\mathrm{m}\right]\hfill \\ Y^{\mathrm{r}}=Y_{t=30}^{\mathrm{r}}\left[\mathrm{m}\right]\hfill \\ {\Theta }^{\mathrm{r}}=0.0\left[\mathrm{rad}\right]\hfill \end{array}\right.\left\{\begin{array}{c}40<t\le 50\left[\mathrm{s}\right]\hfill \\ X^{\mathrm{r}}=X_{t=40}^{\mathrm{r}}\left[\mathrm{m}\right]\hfill \\ Y^{\mathrm{r}}=Y_{t=40}^{\mathrm{r}}-0.25(t-40)\left[\mathrm{m}\right]\hfill \\ {\Theta }^{\mathrm{r}}=0.0\left[\mathrm{rad}\right]\hfill \end{array}\right.\end{array}$$

(21)

where $(X_{t=30}^{\mathrm{r}},Y_{t=30}^{\mathrm{r}})$ and $(X_{t=40}^{\mathrm{r}},Y_{t=40}^{\mathrm{r}})$ represent the reference positions of the robot at 30s and 40s, respectively.

The cyan curve in Figure 8 shows the complete plot of this composite reference trajectory. Its first part ($0\le t\le 30\mathrm{s}$) is a directed curve and the reference orientations are its tangential angles (as shown by the green arrows on the left side in Figure 8); and the second part ($30\mathrm{s}<t\le 50\mathrm{s}$) is an undirected polyline trajectory for which the corresponding reference orientation always remains unchanged (as shown by the green arrows on the right side in Figure 8). In addition, the key segmentation points obtained according to the method proposed in Section 4.1 are also marked by the pink solid points in Figure 8.

5.2. Effect of the Constraints and the Replanning Strategy

Suppose the actual starting pose of the robot relative to the world coordinate system is set as ${\mathit{z}}_0=$${[X_0,Y_0,{\Theta }_0]}^{\mathrm{T}}=$${[0.0\mathrm{m},-1.0\mathrm{m},0.0\mathrm{rad}]}^{\mathrm{T}}$. The simulation process runs at a frequency of 100Hz (i.e., $t_s=0.01\mathrm{s}$). Based on the trial and error method [18], the parameters related to MPC are set as: $N_{\mathrm{p}}=30$, $N_{\mathrm{c}}=2$, $\mathit{Q}=10{\mathit{I}}_3$, and $\mathit{R}=500{\mathit{I}}_3$. In the following, the role and influence of the limb constraints established in Section 3.1 and the replanning strategy of stride proposed in Section 4 in the trajectory tracking of multi-legged robots are illustrated by simulations in MATLAB.

First, let us look at the role of limb constraints. Assuming that the same reference stride and the same controller parameters are used in the comparison cases, the simulation results with and without limb constraints are shown in Figure 9 and Figure 10. Figure 9 shows that, after a period of time, both the control methods with and without limb constraints can track the reference trajectory well. In addition, at the beginning stage of the tracking process, the motion trajectory without limb constraints can approach the reference position trajectory faster. However, this “seemingly better” tracking result ignores the robot’s locomotion ability. As shown in Figure 10, the angles of Joint3 calculated by the model without limb constraints exceed its rotational range. This means that when the stretch and yaw constraints of limbs are not considered, the model predictive controller outputs excessive strides in order to get closer to the reference trajectory faster, which exceeds the actual motion capability of the robot, and thus easily damages the robot’s structure. Figure 10 intuitively shows that all joint angles corresponding to the control method with limb constraints vary within their respective feasible ranges, which fully demonstrates that it is crucial to pay attention to limb constraints in controller design in this paper.

Next, let us analyze the effect of stride replanning. Using the same constrained predictive controller (14)–(15), the trajectory tracking results of the robot under the action of the common reference stride (19) and the replanned reference stride (20) are shown in Figure 11 and Figure 12. Figure 11 shows that when the initial error of the robot’s pose is large, the common reference stride defined by the reference trajectory easily leads to the robot sacrificing the position tracking accuracy in pursuit of smaller yaw error and smaller orientation error, resulting in low tracking efficiency. While the replanned reference stride can guide the robot to move towards the reference trajectory faster, because the robot can dynamically adjust its movement direction in time according to its real-time pose. Figure 12 demonstrates the above point more intuitively from the perspective of pose tracking error. That is, the tracking errors corresponding to the replanned reference stride are smaller than that of the common reference stride regardless of the translation or rotation on the plane. This means that the proposed replanning method of reference stride can effectively improve the tracking efficiency and tracking quality of the omnidirectional reference trajectory of the multi-legged robot.

5.3. Tracking Test

To reduce the loss of components and the cost of experiments, a virtual robot prototype with the same properties as the physical machine and as realistic as possible is built in the virtual robot experiment platform CoppeliaSim (also known as V-REP), and further used to test the practical application effect of the proposed constrained MPC algorithm.

We enable the Newton Dynamic engine in the CoppeliaSim environment, and set the robot WelCH to perform periodic rhythmic motion with tripod gait. Then, we directly transplant the tracking control algorithm proposed in this paper and the control parameters set in the numerical simulation to the robot prototype. Through the synchronous communication mechanism between MATLAB and CoppeliaSim to control the robot movement, the tracking results of WelCH in the simulation environment can be obtained, as shown in Figure 13. Figure 13a shows several snapshots of the robot WelCH during tracking of the reference trajectory. In addition, Figure 13b is a complete overview of the robot tracking the composite reference trajectory, which is consistent with the expected effect. Moreover, the actual position of the robot and the real-time replanning results recorded during the test are plotted in Figure 14. The experimental data again intuitively show that under the guidance of the MPC algorithm designed in Section 3 and the replanning strategy proposed in Section 4, the hexapod robot WelCH can effectively carry out omnidirectional motion and track any omnidirectional reference trajectory quickly and accurately. This proves that the tracking control algorithm proposed in this paper has strong portability and practicability.

6. Conclusions

In this paper, an MPC-based algorithm is proposed to solve the trajectory tracking problem of an unmanned hexapod robot with nonlinear stride constraints. For a hexapod robot that takes motion stride as control input and moves with tripod gait, the models of body-level-tracking control and the stretch and yaw constraints of limb about stride variables are studied and established. These constraints, together with input saturation constraints, serve as quantitative restrictive conditions for designing the optimal predictive controller to ensure that the robot moves without violating its own structural limits. In addition, in order to give full play to the movement capability of the robot in its feasible motion space, the reference stride length defined by the integral mean of surface is proposed; then, the variable stride periods and the key segmentation points of reference trajectory are determined, and then, the real-time replanned reference strides are calculated based on it, which effectively weakens the influence of the deviation of reference inputs on the tracking control. The hexapod robot WelCH adopting the designed tracking control strategy realizes the omnidirectional accurate tracking of composite trajectories, which demonstrates the rationality and effectiveness of the proposed tracking scheme.

Based on the research ideas and results presented in this paper, there are still some more challenging works worthy of further exploration. Firstly, this work only considers the stride constraints when the hexapod robot moves with the tripod gait. Does this quantitative result also apply to other rhythmic or free gaits? If not, how to make such a modification and generalization? Secondly, how to quantify the anti-collision constraints between multiple limbs and the physical constraints related to the self-collision between the robot’s components by using motion stride? In addition, in scenarios that focus on the stationarity and safety of robot motion, how to model dynamic constraints such as torque saturation and friction cone constraints of the support feet? Furthermore, if the degrees of freedom of the body’s fluctuation and pitching are also considered, how should the process of modeling and control in this paper be optimized and perfected?