2.1. System analysis

In this section, a description of the hexapod robot ATHENA considered in this study is presented. The physical prototype of the ATHENA robot is depicted in Fig. 1. From the mechanical point of view, the robot is composed of 25 rigid bodies, which are interconnected by 24 kinematic joints, as can be observed in the corresponding multibody model illustrated in Fig. 2. In this model, body b is the torso, being the remaining bodies axisymmetrically distributed around it. Each limb includes three revolute joints that are actuated by servo motors with a stall torque of 1.89 N·m. The revolute joints aim at connecting the torso with the coxa, femur, and tibia segments – c _i , f _i , and t _i – which are, respectively, named Thorax-Coxa (TC), Coxa-Trochanterofemur (CTr), and Femur-Tibia (FTi). Besides that, the parameter p _i represents the hexapod feet, which are linked with t _i through prismatic joints. The passive connection between the two bodies aims at controlling the feet’s relative displacement caused by their interaction with the ground, using compression springs with a spring rate of 9.28 × 10³ N·m. The displacement associated with the compression spring does not provide significant position deviations in the foot trajectory. These two components are combined, and the relative motion between them is considered to be null. The inertia properties of the bodies of the ATHENA robot are listed in Table I. The presented principal and secondary moments of inertia are estimated using the local reference frame of each body, which is attached to the respective center of mass.

Figure 2. Schematic representation of the multibody model of the ATHENA robot.

2.2. Model design for the computational analysis

In order to have an efficient and accurate contact formulation between the feet and the ground, a simplified multibody model of the robot was developed and implemented in the Gazebo software. Bearing in mind that the goal of this work is to examine the presented posture control in inclined surfaces, the system geometry does not have a significant influence on the results, as long as the inertial properties are equivalent to the physical prototype. Therefore, the hexapod was modeled in Gazebo using simple geometries, such as spheres and planes (see Fig. 3). Most of the components were converted into prisms. On the contrary, the feet were modeled as spheres and therefore, during the collision with the ground, there is only one contact point, which is also beneficial for the computational point of view. All revolute joints were set to a maximum torque of 1.89 N·m, to simulate the selected servo motors and the connection between the tibia and the foot is assumed fixed. Table II lists the model dimensions and the respective inertial properties. The presented principal and secondary moments of inertia are estimated using the local reference frame of each body, which is attached to the respective center of mass.

3.1. Equations of motion

This section discusses the defined equations of motion for the hexapod robot that are used to assess the control behavior. The computational simulations performed in Gazebo resort to the external library Open Dynamics Engine (ODE) to estimate the motion and forces applied to a system. The general equations of motion involving contact-impact events are expressed as [Reference Flores, Leine and Glocker46]:

(1) \begin{equation} \mathbf{M}\ddot{\mathbf{q}}-\mathbf{N}_{\mathrm{n}}\mathbf{f}_{\mathrm{n}}-\mathbf{N}_{\mathrm{t}}\mathbf{f}_{\mathrm{t}}=\mathbf{0} \end{equation}

in which $\mathbf{M}^{(150,150)}$ is the system’s global mass matrix, $\ddot{\mathbf{q}}^{(150)}$ denotes the vector that contains the bodies’ linear and angular accelerations, the matrices $\mathbf{N}_{\mathrm{n}}, \mathbf{N}_{\mathrm{t}}\in \mathrm{\mathbb{R}}^{150\times n}$ contain respectively the normal and tangent vectors with respect to the contact surfaces, and $\mathbf{f}_{\mathrm{n}}, \mathbf{f}_{\mathrm{t}}\in \mathrm{\mathbb{R}}^{n}$ are respectively the n-dimensional vector of the normal and tangential contact forces. It must be mentioned that the dimension of $\mathbf{f}_{\mathrm{n}}$ and $\mathbf{f}_{\mathrm{t}}$ depends on the number of contacting points. The equations of motion are calculated at the impulse level. In each time step, contact-impact events are detected and the equations of motion of the model are evaluated with respect to the bodies’ velocities. The numeric simulation’s time discretization is based on the Moreau’s time-stepping methodology [Reference Moreau47].

Table I. Inertial properties of ATHENA robot.

Figure 3. Computational multibody model of the ATHENA robot used in the dynamic analysis.

Table II. Geometrical and inertial properties of the simplified model.

In turn, the contact resolution closely follows the approach proposed by Stewart-Trinkle [Reference Stewart and Trinkle48] and Anitescu-Potra [Reference Anitescu and Potra49], which is formulated as a Linear Complementarity Problem (LCP). This type of formulation can be expressed as:

(2) \begin{equation} \mathbf{y}=\mathbf{Ax}+\mathbf{b} \end{equation}

where $\mathbf{A}^{(n,n)}$ and $\mathbf{b}^{(n)}$ are constants with n-dimension,and $\mathbf{x}$ and $\mathbf{y}$ denote the system’s unknowns [Reference Flores, Leine and Glocker46]. In order to obtain the quantities of x and y, the following complementary conditions must be verified [Reference Flores, Leine and Glocker46]:

(3) \begin{equation} \mathbf{y}^{\mathrm{T}}\mathbf{x}=0\wedge \mathbf{y}\geq \mathbf{0}\wedge \mathbf{x}\geq \mathbf{0} \end{equation}

This methodology is implemented for the calculation of the normal contact forces. In this strategy, the normal contact forces are modeled using Newton’s impact law, which considers the existence of compression and restitution phases during the collision between two bodies. When the compression phase takes place, there is impulse conservation. The model’s equations of motion are expressed at level velocity and impulse level can be formulated as [Reference Stewart and Trinkle48, Reference Anitescu and Potra49]:

(4) \begin{equation} \mathbf{M}\left(\dot{\mathbf{q}}^{0}-\dot{\mathbf{q}}^{-}\right)-\left(\mathbf{N}_{\mathrm{n}}\hat{\mathbf{f}}_{\mathrm{n}}^{0}+\mathbf{N}_{\mathrm{t}}\hat{\mathbf{f}}_{\mathrm{t}}^{0}\right)=\mathbf{0} \end{equation}

where the superscript 0 denotes the end of the compression phase and $\dot{\mathbf{q}}^{-}$ is the pre-collision velocity. The variables $\hat{\mathbf{f}}_{\mathrm{n}}$ and $\hat{\mathbf{f}}_{\mathrm{t}}$ denote the impulses of the normal and tangential contact forces, respectively. During the restitution phase, the bodies move in opposite directions, and the equations of motion at the velocity and impulse level are expressed as:

(5) \begin{equation} \mathbf{M}\left(\dot{\mathbf{q}}^{+}-\dot{\mathbf{q}}^{0}\right)-\left(\mathbf{N}_{\mathrm{n}}\hat{\mathbf{f}}_{\mathrm{n}}^{0}+\mathbf{N}_{\mathrm{t}}\hat{\mathbf{f}}_{\mathrm{t}}^{0}\right)=\hat{\mathbf{f}}_{\mathrm{r}} \end{equation}

in which $\dot{\mathbf{q}}^{+}$ is the post-collision velocity.

The vector $\hat{\mathbf{f}}_{\mathrm{r}}$ , which represents the impulsive force, can be defined as:

(6) \begin{equation} \hat{\mathbf{f}}_{\mathrm{r}}=\xi \mathbf{N}_{\mathrm{n}}\hat{\mathbf{f}}_{\mathrm{n}}^{0} \end{equation}

where $\xi \in [0,1]$ is the coefficient of restitution. In this study, $\xi$ equals 0.5. In order to solve the Newton’s impact law as an LCP, the following complementary conditions are defined [Reference Stewart and Trinkle48, Reference Anitescu and Potra49]:

(7) \begin{equation} \mathbf{N}_{\mathrm{n}}\dot{\mathbf{q}}=\mathbf{0}\wedge \hat{\mathbf{f}}_{\mathrm{n}}\geq \mathbf{0} \end{equation}

(8) \begin{equation} \mathbf{N}_{\mathrm{n}}\dot{\mathbf{q}}\geq \mathbf{0}\wedge \hat{\mathbf{f}}_{\mathrm{n}}=\mathbf{0} \end{equation}

These conditions evaluate both the normal relative distance between bodies and the consequent existence of a normal contact force. If $\mathbf{N}_{\mathrm{n}}\dot{\mathbf{q}}$ is null, there is a collision, and a positive normal contact force exists. On the contrary, if the $\mathbf{N}_{\mathrm{n}}\dot{\mathbf{q}}$ is positive, then the bodies are separated, and the normal contact force is null.

A similar strategy is implemented for the tangential contact analysis. In this study, the tangential contact forces are modeled using Coulomb’s friction law. In order to increase the model accuracy, the isotropic friction model is approximated to a polyhedral [Reference Chakraborty, Berard, Akella and Trinkle50]. The complementary conditions that are set to solve the friction forces as an LCP are expressed as:

(9) \begin{equation} \boldsymbol{\lambda }+\mathbf{N}_{\mathrm{n}}\dot{\mathbf{q}}\geq \mathbf{0}\,\bot\, \hat{\mathbf{f}}_{\mathrm{t}}\geq \mathbf{0} \end{equation}

(10) \begin{equation} \mu \hat{\mathbf{f}}_{\mathrm{n}}-\hat{\mathbf{f}}_{\mathrm{t}}\geq \mathbf{0}\,\bot\, \boldsymbol{\lambda }\geq \mathbf{0} \end{equation}

where $\boldsymbol\lambda$ denotes the Lagrange multipliers, and μ is the friction coefficient, which is equal to 0.9. The presented inequalities state that if $\boldsymbol\lambda$ is greater than 0, there is sliding between the two surfaces. Otherwise, there is stick slippage between the surfaces [Reference Flores, Leine and Glocker46].

3.2. Kinematic model

This section presents the kinematic modeling process to control the robot. The implementation of the proposed method requires obtaining the relation between the joints’ angular positions and the robot’s feet coordinates. Considering the models presented in refs. [Reference Coelho, Dias, Lopes, Ribeiro and Flores37, Reference Coelho, Dias, Lopes, Ribeiro and Flores51], the torso’s center of mass is defined as the reference for all calculations. The kinematic model of one limb is depicted in Fig. 4. The angles $\theta _{2}^{\mathrm{r}}, \theta _{3}^{\mathrm{r}}$ , and $\theta _{4}^{\mathrm{r}}$ are the relative angles of the TC, CTr, and FTi joints, respectively. The hexapod motion is studied by setting the motion range of the TC joint to $\theta _{2}^{\mathrm{r}}\in [\!-\!0.612,0.612]$ rad, and the angular position of the remaining joints is constrained to $\theta _{3}^{\mathrm{r}},\theta _{4}^{\mathrm{r}}\in [\!-\!1.484,1.484]$ rad.

Figure 4. Schematic representation of the joints’ angular motion.

The locomotion control requires two different estimations from the kinematic model: (i) get the current feet coordinates based on the joints’ angular positions to verify if the desired motion is within the system’s workspace, and (ii) calculate the joints’ configuration according to the trajectory applied to the feet. In the first analysis, the Denavit–Hartenberg convention is used to obtain the relative position of the foot w.r.t. the reference frame B. Using the Denavit–Hartenberg parameters presented in Table III, the relative transformation between B and the hexapod feet can be established as:

(11) \begin{equation} T_{4}^{0}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} \cos\!\left(\theta _{1}^{r}+\theta _{2}^{r}\right)\cos\!\left(\theta _{3}^{r}+\theta _{4}^{r}\right) & -\!\cos\!\left(\theta _{1}^{r}+\theta _{2}^{r}\right)\sin\!\left(\theta _{3}^{r}+\theta _{4}^{r}\right) & \sin\!\left(\theta _{1}^{r}+\theta _{2}^{r}\right) & x_{\mathrm{p}}^{B}\\[5pt] \sin\!\left(\theta _{1}^{r}+\theta _{2}^{r}\right)\cos\!\left(\theta _{3}^{r}+\theta _{4}^{r}\right) & -\!\sin\!\left(\theta _{1}^{r}+\theta _{2}^{r}\right)\sin\!\left(\theta _{3}^{r}+\theta _{4}^{r}\right) & -\!\cos\!\left(\theta _{1}^{r}+\theta _{2}^{r}\right) & y_{\mathrm{p}}^{B}\\[5pt] \sin\!\left(\theta _{3}^{r}+\theta _{4}^{r}\right) & \cos\!\left(\theta _{3}^{r}+\theta _{4}^{r}\right) & 0 & z_{\mathrm{p}}^{B}\\[5pt] 0 & 0 & 0 & 1 \end{array}\right] \end{equation}

where the vector $\mathbf{p}_{i}^{B}=\{x_{\mathrm{p}}^{B},y_{\mathrm{p}}^{B},z_{\mathrm{p}}^{B}\}^{\mathrm{T}},i\in 1,\ldots,6$ contains the feet’ relative coordinates, which are expressed as:

(12) \begin{equation} \mathbf{p}_{\mathrm{i}}^{B}=\left\{\begin{array}{c} \cos\!\left(\theta _{1}^{\mathrm{r}}+\theta _{2}^{\mathrm{r}}\right)\left(l_{4}\cos\!\left(\theta _{3}^{\mathrm{r}}+\theta _{4}^{\mathrm{r}}\right)+l_{3}\cos \theta _{3}^{\mathrm{r}}+l_{2}\right)+l_{1}\cos \theta _{1}^{\mathrm{r}}\\[5pt] \sin\!\left(\theta _{1}^{\mathrm{r}}+\theta _{2}^{\mathrm{r}}\right)\left(l_{4}\cos\!\left(\theta _{3}^{\mathrm{r}}+\theta _{4}^{\mathrm{r}}\right)+l_{3}\cos \theta _{3}^{\mathrm{r}}+l_{2}\right)+l_{1}\sin \theta _{1}^{\mathrm{r}}\\[5pt] l_{4}\sin\!\left(\theta _{3}^{\mathrm{r}}+\theta _{4}^{\mathrm{r}}\right)+l_{3}\sin \theta _{3}^{\mathrm{r}} \end{array}\right\} \end{equation}

The limb’s inverse kinematics calculation is essential to control the driven joints’ positions. The inverse kinematics formulation uses the feet’s position to obtain the relative angular positions of the joints. For this matter, the value of $\theta _{1}^{\mathrm{r}}$ is constant and depends on the leg position as:

(13) \begin{equation} \theta _{1}^{\mathrm{r}}=\left\{\begin{array}{c} 0,i=1\\[5pt] \text{atan}\left(\dfrac{l_{0}}{w}\right), i=2\\[12pt] \pi -\text{atan}\left(\dfrac{l_{0}}{w}\right), i=3\\[8pt] \pi, i=4\\[8pt] \pi +\text{atan}\left(\dfrac{l_{0}}{w}\right), i=5\\[12pt] -\text{atan}\left(\dfrac{l_{0}}{w}\right), i=6 \end{array}\right. \end{equation}

The calculation of $\theta _{2}^{\mathrm{r}}$ depends on the fixed reference frame Q, which is placed on the TC joint (see Fig. 5). Thus, the angular position of this driven joint is defined as:

(14) \begin{equation} \theta _{2}^{\mathrm{r}}=\text{atan}\left(\dfrac{x_{\mathrm{p}}^{Q}}{y_{\mathrm{p}}^{Q}}\right) \end{equation}

in which $\mathbf{p}_{i}^{Q}=\{x_{\mathrm{p}}^{Q},y_{\mathrm{p}}^{Q},z_{\mathrm{p}}^{Q}\}^{\mathrm{T}}$ denotes the relative position of p _i w.r.t. Q, and is expressed as:

(15) \begin{equation} \mathbf{p}_{i}^{Q}=\left(\mathbf{T}_{1}^{0}\right)^{-1}\mathbf{p}_{i}^{B} \end{equation}

where $\mathbf{T}_{1}^{0}$ is the transformation matrix between the reference frames B and Q. The estimation of the remaining angular positions requires knowing the position of the reference frame U, which is placed on the CTr joint. Using the Denavit–Hartenberg convention, the position of U w.r.t. Q is defined as:

Table III. Denavit–Hartenberg parameters for all relative transformations between the reference frame B and the hexapod’s foot $\mathbf{p}_{i}$ .

Figure 5. Schematic representation of the limbs model.

(16) \begin{equation} \mathbf{u}_{i}^{Q}=\left(\mathbf{T}_{1}^{0}\right)^{-1}\mathbf{T}_{2}^{0} \end{equation}

Considering the vector s _p between the reference frame U and the foot, the value of the CTr and the FTi joints is as follows:

(17) \begin{equation} \theta _{3}^{\mathrm{r}}=\text{acos}\left(\frac{-l_{4}^{2}+l_{3}^{2}+\left\| \mathbf{s}_{\mathrm{p}}\right\| ^{2}}{2l_{3}\left\| \mathbf{s}_{\mathrm{p}}\right\| }\right)-\text{asin}\left(\frac{z^{U}-z^{\mathrm{p}}}{\left\| \mathbf{s}_{\mathrm{p}}\right\| }\right) \end{equation}

(18) \begin{equation} \theta _{4}^{\mathrm{r}}=\pi +\text{acos}\left(\frac{-\left\| \mathbf{s}_{\mathrm{p}}\right\| ^{2}+l_{3}^{2}+l_{4}^{2}}{2l_{3}l_{4}}\right) \end{equation}

where z ^U and z ^p represent, respectively, the vertical position of the reference frame U and the foot. The presented equations of motion and kinematic model of the ATHENA robot contain fundamental concepts for the design and assessment of the developed posture control strategy.

4.1. Gait generation

The locomotion of the ATHENA robot adopts the tripod gait. Among the commonly adopted gaits in hexapod robots, the tripod pattern provides the highest velocity of the main body, despite decreasing the system’s stability due to the number of supporting limbs. If the designed control for the posture adjustment provides good results for the tripod gait, then the methodology is considered assessed for other locomotion types. Considering the gait state, the limbs execute a different motion. In the swing phase, the feet reach a new foothold position by executing a cubic Bézier trajectory expressed as:

(19) \begin{equation} \mathbf{p}_{i}^{B}\left(t\right)=\sum _{j=0}^{n}\left(\mathrm{t}_{\mathrm{g}}-t\right)^{n-j}t^{j}\mathbf{p}_{i,j}^{B}, t\in \left[0,\mathrm{t}_{\mathrm{g}}\right]\wedge i\in 1,\ldots,6 \end{equation}

where n represents to the number of control points, t_g is the trajectory’s time interval, and $\mathbf{p}_{i,j}^{B}$ denote the control points of the trajectory. The first control point corresponds to the foot’s initial position $\mathbf{p}_{i,0}^{B}$ which is defined by Eq. (12). The remaining control points are determined as:

(20) \begin{equation} \mathbf{p}_{i,1}^{B}=\left\{0,\frac{s}{4},\frac{h}{3}\right\}^{\mathrm{T}}+\mathbf{p}_{i,0}^{B} \end{equation}

(21) \begin{equation} \mathbf{p}_{i,2}^{B}=\left\{0,\frac{3}{4}s,\frac{h}{3}\right\}^{\mathrm{T}}+\mathbf{p}_{i,0}^{B} \end{equation}

(22) \begin{equation} \mathbf{p}_{i,3}^{B}=\left\{0,s,0\right\}^{\mathrm{T}}+\mathbf{p}_{i,0}^{B} \end{equation}

where h is the maximum height and s denotes the stride.

In the stance phase, the feet are in contact with the ground while propelling the torso forward. Thus, the foot’s trajectory is defined as:

(23) \begin{equation} \mathbf{p}_{i}^{B}\left(t\right)=\left\{0,-\frac{s}{\mathrm{t}_{\mathrm{g}}}t,0\right\}^{\mathrm{T}}+\mathbf{p}_{i,0}^{B}, t\in \left[0,\mathrm{t}_{\mathrm{g}}\right] \end{equation}

In both phases, the result of the foot trajectory must be converted into motor commands, using the inverse kinematics model established in Section 3.2.

4.2. Posture adjustment

The aforementioned footpath planning of the limbs considers that the hexapod robot walks across a regular and smooth surface. Thus, the displacement of the torso’s height and orientation is neglected, and consequently the system is stable. Nonetheless, the variation of the terrain topology affects the hexapod’s achievable foothold positions and corresponding torso’s posture. The terrain complexity increases the displacement of the torso’s position, which constrains the hexapod’s ability to surmount inclined surfaces. Bearing that in mind, the proposed strategy aims at adjusting the hexapod’s height and orientation by estimating the ground inclination.

The terrain perception is estimated by the relative position of the feet w.r.t. the torso, which is given by Eq. (12). For this purpose, the control evaluates the relative coordinates of the feet, which are in contact with the ground. Using the relative coordinates of the supporting feet, the terrain can be approximated to a characteristic plane expressed as:

(24) \begin{equation} a_{\mathrm{p}}\mathbf{x}+b_{\mathrm{p}}\mathbf{y}+c_{\mathrm{p}}\mathbf{z}-d_{\mathrm{p}}=\mathbf{0} \end{equation}

in which a _p, b _p, c _p, and d _p are constants, and $\mathbf{x}, \mathbf{y}, \mathbf{z}\in {\mathbb{R}}^{f}$ are f-dimensional vectors that contain the coordinates of the supporting feet. With the values of a _p, b _p, c _p, and d _p, the norm vector of the plane g is obtained. The developed algorithm for the detection of the supporting limbs, and the calculation of the norm vector, is described by the following steps (see Fig. 6). In each gait cycle, the algorithm reads the normal contact forces using the force sensors placed on each foot. In a regular surface, there are exactly three legs with positive normal contact forces; hence, the ground plane equation can be estimated using Eq. (24). Nonetheless, the number of supporting limbs changes with the terrain irregularities. If the terrain topology forces the hexapod to have more than three supporting limbs, a random subset with three of the legs in contact with the ground is applied for the estimation of the ground plane. In case of having less than three legs in the stance phase, the algorithm searches among the remaining limbs which should be in contact with the ground, and their trajectory is extended until a positive normal contact force is read by the force sensors. After that, there are sufficient supporting limbs to estimate the terrain’s characteristic plane using Eq. (24). Figure 7 shows a generic configuration of the ATHENA hexapod walking across an inclined and irregular surface. Considering the torso’s local reference, and a two-dimensional representation of g in the $\eta$ $\zeta$ and $\xi$ $\zeta$ planes, the angular displacement of the pitch and roll angles is determined as:

(25) \begin{equation} \delta _{\alpha}=\arccos \left(\frac{\mathbf{g}'\mathbf{n}_{{\zeta} }}{\left\| \mathbf{g}'\right\| \left\| \mathbf{n}_{{\zeta} }\right\| }\right),\mathbf{g}'=\left\{x_{\mathrm{g}},z_{\mathrm{g}}\right\}^{\mathrm{T}}\wedge \mathbf{n}_{{\zeta} }=\left\{0,1\right\}^{\mathrm{T}} \end{equation}

(26) \begin{equation} \delta _{\alpha }=\arccos \left(\frac{\mathbf{g}'\mathbf{n}_{{\zeta} }}{\left\| \mathbf{g}'\right\| \left\| \mathbf{n}_{{\zeta} }\right\| }\right),\mathbf{g}'=\left\{y_{\mathrm{g}},z_{\mathrm{g}}\right\}^{\mathrm{T}} \end{equation}

where δ $_\alpha$ represents the pitch displacement, $\boldsymbol{n}_{\zeta}$ denotes a unitary vector, and δ _β is the angular deviation in the roll orientation.

Figure 6. Schematic representation of the contact detection algorithm.

Figure 7. Schematic representation of the plane formed by the supporting limbs.

Figure 8. Angular displacement between the torso and the ground with simple control methodology.

The terrain’s inclination affects the relative motion between the torso and the ground. Figure 8 shows the angular displacement performed by the ATHENA hexapod climbing a 10-degree ramp in Gazebo, without controlling the torso posture. From the analysis of this plot, it can be observed that in each gait cycle, the transition between gait phases causes an abrupt change of the roll and pitch angles due to the alternation between the supporting limbs and consequent inaccurate planning of the foothold positions. Thus, it is of paramount importance to ensure that during the gait phase transition, all limbs are in contact with the ground and their height is adjusted to the estimated surface plane. Subsequently, the results variation of the roll and pitch angles in each gait cycle assess the torso’s instability when the hexapod climbs the ramp. The posture variation is influenced by the limbs’ trajectory in the stance phase, which is responsible for adjusting the torso’s angular positions. Considering the outcome from the computational simulation with simple control, it can be concluded that the proposed control strategy is able to adjust the feet’ coordinates regarding the torso’s orientation and height. Since the swing phase is mainly influenced by the torso’s height, the final control point of the cubic Bézier trajectory is defined as:

(27) \begin{equation} \mathbf{p}_{i,3}^{O}=\left\{0,s,k_{\mathrm{z}}\delta _{\mathrm{z}}\right\}+\mathbf{p}_{i,0}^{O}, k_{\mathrm{z}}\in \left[0,1\right] \end{equation}

where k _z is a control parameter, and δ _z denotes the height displacement, being the difference between the relative height w.r.t. the ground measured by an infrared sensor placed on the torso, and the reference height, which denotes the torso’s height considering the hexapod’s standing posture. Regarding the stance phase, the control strategy requires an adjustment using Eq. (23) according to the torso orientation in inclined surfaces. Thus, the feet’s final position in the stance phase is expressed as follows:

(28) \begin{equation} \mathbf{p}_{i,\mathrm{f}}^{O}=\mathbf{t}_{\mathrm{r}}+\mathbf{p}_{i,0}^{O} \end{equation}

where t _r denotes the transformation vector, that are expressed by:

(29) \begin{equation} \mathbf{t}_{\mathrm{r}}=\mathbf{R}_{\mathrm{x}}^{\left(4,4\right)}\left(k_{\alpha}\delta _{\alpha}\right)\mathbf{R}_{\mathrm{y}}^{\left(4,4\right)}\left(k_{\beta }\delta _{\beta }\right)\mathbf{T}^{\left(4,4\right)}\left(k_{\mathrm{z}}\delta _{\mathrm{z}}\right), k_{\alpha}, k_{\beta }\in \left[0,1\right] \end{equation}

where R _x and R _y denote the rotation matrix along the axes x and y, respectively, k $_\alpha$ and k _β are the corresponding control parameters for the pitch and roll angles, and T is the transformation matrix defined as:

(30) \begin{equation} \mathbf{T}^{\left(4,4\right)}=\left[\begin{array}{cc} \mathbf{I}^{\left(3,3\right)} & \mathbf{t}\\[5pt] 0 & 1 \end{array}\right], \mathbf{t}=\left\{0,-s,\mathrm{k}_{\mathrm{z}}\delta _{\mathrm{z}}\right\}^{\mathrm{T}} \end{equation}

Bearing in mind that the transformation of $\mathbf{p}_{i,\mathrm{f}}^{O}$ , the stance phase is defined as:

(31) \begin{equation} \mathbf{p}_{i}^{O}\left(t\right)=\left\{\frac{s_{\mathrm{x}}}{\mathrm{t}_{\mathrm{g}}}t,\frac{s_{\mathrm{y}}}{\mathrm{t}_{\mathrm{g}}}t,\frac{s_{\mathrm{z}}}{\mathrm{t}_{\mathrm{g}}}t\right\}+\mathbf{p}_{i,0}^{O}, t\in \left[0,\mathrm{t}_{\mathrm{g}}\right]\wedge \mathbf{t}_{\mathrm{r}}\equiv \left\{s_{\mathrm{x}},s_{\mathrm{y}},s_{\mathrm{z}}\right\}^{\mathrm{T}} \end{equation}

5.1. Computational setup

The proposed computational procedure to control the ATHENA robot has been implemented and evaluated in the Gazebo software. The ODE library is considered in the computational simulations with a time step equal to 50 ms to solve the equations of motion and provide data regarding the torso posture. The controller implementation in the Gazebo environment uses a ROS-based client–server infrastructure with a rate of 20 Hz to simulate real-time data acquisition. An IMU and an infrared sensor are placed on the hexapod’s torso, to assess the body’s angular velocities and height, which are estimated by the ODE library.

The computational simulations are performed for a 10-degree ramp and an irregular terrain. The two scenarios are presented in Fig. 9. The surface’s inclination defined for the ramp scenario is representative and is appropriate to assess the control methodology (see Fig. 9(a)). A similar strategy is utilized in the testbed scenario portrayed in Fig. 9(b). The second situation contains blocks with different heights between 0.04 m and 0.13 m and aims at studying the proposed control in a more demanding scenario in terms of terrain topology variation. In both computational cases, the hexapod adopts a tripod gait with t_g = 1 s, and the kinetic friction coefficient is set to 0.9. For the height adjustment, the hexapod’s standing posture, which is used as reference for the infrared sensor readings, is set when $\theta _{2}^{\mathrm{r}}=0$ rad, $\theta _{3}^{\mathrm{r}}=0$ rad, and $\theta _{4}^{\mathrm{r}}=-\frac{\pi }{2}$ rad.

Figure 9. Environments tested in the computational simulations: (a) Ramp with 10-degree slope and (b) Irregular terrain.

Figure 10. Results obtained for the height control in the ramp scenario: (a) Height variation; (b) Pitch velocity variation; and (c) Roll velocity variation.

5.2. Discussion

In each environment, the influence of the parameters k _z, k $_\alpha$ , and k _β on the control system’s response is analyzed. It must be noticed that the value of the parameters can be defined by an infinite number of combinations. The main goal of this work is to assess the feasibility of the proposed control; hence, k _z, k $_\alpha$ , and k _β are restrained to a set of discrete values (k _z, k $_\alpha$ , k _β ∈ {0.0, 0.2, 0.4, 0.6, 0.8, 1.0}), which provides a general overview into the influence of the constants k _z, k $_\alpha$ , and k _β in the posture stability in the two defined scenarios. The hexapod’s posture is evaluated according to the torso’s height w.r.t. the ground and the angular velocities in the roll and pitch directions. In each computational simulation, the root-mean-square value of the aforementioned variables is estimated per gait cycle, and the resultant variation is examined.

5.2.1. Ramp scenario

The first computational simulations that were performed consider k $_\alpha$ = k _β = 0.0 and k _h ∈ {0.0, 0.2, 0.4, 0.6, 0.8, 1.0}. If k _z = 0.0 takes place, the control system does not adjust the torso’s height, and when k _z = 1.0, there is full implementation of the estimated δ _z in the limbs’ trajectories. The obtained results for the relative height and angular velocities are presented in Fig. 10. Using the simple control as reference (k _z = 0.0), it can be observed that the implementation of the height adjustment improves the hexapod robot’s stability due to the significant decrease of the height’s interquartile range when k _z > 0.0 (see Fig 10(a)). In comparison to the median value for k _z = 0.0, which corresponds to 0.068 m, there is, on average, an increase of 58.42% in the torso’s height when the proposed control strategy is implemented. Although the highest median value occurs when k _z = 1.0, which corresponds to 0.112 m, the height variation is smaller when k _z ∈ {0.6, 0.8}. The height adjustment does not affect the torso’s orientation and angular velocity. Figure 10(b) depicts the plots for the pitch velocity. The obtained results show that the pitch velocity’s median increases when k _z = 1.0. Thus, it can be concluded that, for the maximum height, the variation of the torso’s pitch angle increases, and the hexapod’s posture becomes unstable. For the remaining cases, k _z ∈ {0.2, 0.4, 0.6, 0.8}, the hexapod’s stability improves, since there is has an average reduction of 18.91% of the median value, and a decrease of the interquartile range. Among the obtained results for the pitch velocity, the robot achieves its stable configuration when k _z = 0.8, due to the smallest interquartile range and median value, which corresponds to 0.064 rad/s. Regarding the roll angle, the implementation of k _z = 1.0 increases the variation of the angular velocity, which penalizes the robot’s stability (see Fig. 10(c)). However, in the remaining cases, an improvement of the velocity variation is observed, corresponding to a decrease of 30.74% of the median value. In this case study, the smallest variation is observed when k _z ∈ {0.2, 0.8}.

A similar strategy is implemented for the analysis of the computational simulations performed considering k _z = k _β = 0.0 and k $_\alpha$ ∈ {0.0, 0.2, 0.4, 0.6, 0.8, 1.0} (see Fig. 11). Taking the simple control (k $_\alpha$ = 0.0) as a reference, the adjustment of the torso’s pitch orientation does not have a significant influence on the height’s variation (see Fig. 11(a)). However, for the results regarding the angular velocities in the pitch and roll directions, the configuration k $_\alpha$ = 0.6 provides a good improvement in the system’s stability, both in terms of the interquartile range and the median value, which decreases 39.56% and 35.7%, respectively, in the pitch and roll velocities (see Fig. 11(b) and (c)).

Figure 11. Results obtained for the pitch angle control in the ramp scenario: (a) Height variation; (b) Pitch velocity variation; and (c) Roll velocity variation.

Although control of the torso’s roll angle is not important for the tested environment, the computational simulations considering k _z = k $_\alpha$ = 0.0 and k _β ∈ {0.0, 0.2, 0.4, 0.6, 0.8, 1.0} were also performed. Figure 12 shows the outcome of the computational simulations performed. Similarly to the results presented for k $_\alpha$ , the control of the torso’s roll angle does not have a significant effect on the height variation (see Fig. 12(a)), despite the considerable decrease of the median value when k _β ∈ {0.2, 0.6}. Regarding pitch velocity, the interquartile range between configurations is considered identical, which means that there is no improvement in the system’s stability (see Fig. 12(b)). Nonetheless, the minimum median value takes place when k _β = 0.4 and corresponds to 0.077 rad/s. In terms of the results for the roll velocity, which is plotted in Fig. 12(c), the variable’s median improves for k _β ∈ {0.4, 0.6, 1.0}, which corresponds to a decrease respectively of 14.54%, 2.32%, and 15%. Nonetheless, by analyzing the maximum values and interquartile ranges, it is not possible to conclude that the implementation of k _β > 0.0 provides good results for the hexapod’s posture in the ramp environment.

Figure 12. Results obtained for the roll angle control in the ramp scenario: (a) Height variation; (b) Pitch velocity variation; and (c) Roll velocity variation.

Figure 13. Results obtained for the height control in the irregular terrain scenario: (a) Height variation; (b) Pitch velocity variation; and (c) Roll velocity variation.

5.2.2. Irregular terrain scenario

In comparison to the study presented to the previous section, a similar approach to the performance of the computational simulations is implemented in the testbed scenario. In this case, it is expected to observe a higher instability of the torso due to the terrain irregularities. The first analysis takes into account that k $_\alpha$ = k _β = 0.0 and k _z ∈ {0.0, 0.2, 0.4, 0.6, 0.8, 1.0}. Figure 13 depicts the outcome of the computational simulations performed. Considering the simple control (k _z = 0.0) as a reference for this analysis and observing the interquartile ranges, the implementation of the torso’s height control significantly improves the hexapod robot stability (see Fig. 13(a)). The minimum variation of the results occurs when k _z ∈ {0.6, 0.8, 1.0} and the maximum median value corresponds to 0.107 m, presented when k _z = 0.8. Nonetheless, the variation of torso’s angular velocity increases (see Fig. 13(b) and (c)). Without considering k _z = 0.2, the average increase of the pitch velocity’s median value corresponds to 66.66%. Regarding the roll velocity, the maximum torso’s instability takes place when k _z ∈ {0.2, 0.8}, and the most stable configuration occurs when k _z = 0.4. Thus, despite ensuring a more stable height of the torso, the correction of the feet’ relative height increases the variation of the hexapod’s orientation.

Similarly to the ramp, the influence of k $_\alpha$ on the torso posture does not provide significant enhancement in terms of height stability control (see Fig. 14(a)). Considering k $_\alpha$ = 0.0 as a reference and observing the height variation, the optimal configuration corresponds to k $_\alpha$ = 0.4. By observing the results depicted in Fig. 14(b), the control of the torso’s pitch angle decreases the hexapod’s stability. The most unstable configurations take place when k $_\alpha$ ∈ {0.8, 1.0}, which means that the torso cannot have the same inclination as the ground. Nonetheless, if k $_\alpha$ ∈ {0.4, 0.6}, the median roll velocity of the torso decreases, respectively, 8.43% and 16.00%, providing a possible improvement of the system’s stability (see Fig. 14(c)).

Figure 14. Results obtained for the pitch angle control in the irregular terrain scenario: (a) Height variation; (b) Pitch velocity variation; and (c) Roll velocity variation.

The influence of the roll angle control is studied by considering k _z = k $_\alpha$ = 0.0 and k _β ∈ {0, 0.2, 0.4, 0.6, 0.8, 1.0}. The implementation of this parameter does not have a significant effect on the height variation (see Fig. 15(a)), despite the decrease of 14.85% of the median value when k _β = 0.2. The roll control provides good results in terms of the variation of the pitch velocity (see Fig. 15(b)). When k _β ∈ {0.2, 0.4} is considered, the median value decreases respectively 20.39% and 38.80%. Additionally, there is a consequent decrease in the interquartile range, which shows evidence that the hexapod’s locomotion is adjusted to the terrain topology. Regarding the roll velocity, the interquartile range is improved when compared to k _β = 0.0, except for k _β = 0.2. Thus, in comparison to the ramp environment, it is concluded that the roll angle control plays a key role in the hexapod’s stability for more complex scenarios.

Figure 15. Results obtained for the roll angle control in the irregular terrain scenario: (a) Height variation; (b) Pitch velocity variation; and (c) Roll velocity variation.

Figure 16. Results obtained for the ramp scenario: (a) Height variation; (b) Pitch velocity variation; and (c) Roll velocity variation.

5.2.3. Enhanced posture control for the ramp and irregular terrain scenarios

Considering the results obtained in the previous sections, in that follows computational simulations for two demonstrative applications of the designed posture control. In each analysis, the values of k _z, k $_\alpha$ , and k _β are selected and compared against the outcome of a simple control (k _z = k $_\alpha$ = k _β = 0.0). For the ramp environment, the control of the torso’s height and pitch orientation has a higher impact on the system’s stability than the roll orientation control. Thus, k _β is considered null. Considering the height stability, the value of k _z which presents the most stable results for the height variation also has a negative influence on the torso’s angular velocities. Thus, the ideal value for k _z is 0.8. Regarding the pitch orientation, k $_\alpha$ does not influence the height displacement. This parameter is selected according to its impact on the torso’s orientation. Among the collected data, k $_\alpha$ = 0.6 provides good results for both the pitch and roll velocities. Therefore, the final configuration of the control strategy for the ramp scenario is defined as k _z = 0.8, k $_\alpha$ = 0.6, and k _β = 0. Figure 16 presents the performed computational simulations. Overall, the torso’s posture and stability improved when compared to the locomotion with simple control. The height’s median increased by 60.77% and the interquartile range significantly reduced, which means that the hexapod had a similar height throughout the entire simulation (see Fig. 16(a)). The implementation of the proposed control strategy caused a reduction of 11.92% and 35.98%, respectively, of the median of the pitch and roll angular velocities (see Fig. 16(b) and (c)). Therefore, the angular displacement of the torso decreased, and the gait can be considered optimized.

In the irregular terrain environment, the selection of the control parameters takes into account that the hexapod robot must have a stable height to overcome the terrain. Bearing that in mind, the value of k _z is selected regardless of its influence on the torso’s angular velocities. Thus, the adequate value for the height adjustment’s control parameter is 0.6. The value of k $_\alpha$ does not affect the variation of the height and pitch velocity. Nonetheless, k $_\alpha$ is set to 0.4. In opposite to the ramp scenario, roll control is relevant for ensuring the correct adaptation of the hexapod’s gait to the terrain. The variable k _β is set to 0.4 due to the stability improvement for all variables. The obtained results are presented in Fig. 17. The proposed control shows a different behavior in terms of the control of the torso’s orientation due to the terrain irregularities. When compared against a simple control, the height variation drastically reduced, and the median raised from 0.03 to 0.01 m (see Fig. 17(a)). Nonetheless, the results for the angular velocities were a problem (see Fig. 17(b) and (c)). For the pitch velocity, there is an increase of 54.33% in the median value, and the roll velocity’s median is 82.16% higher. Despite the adjustments of the control strategy, the terrain irregularity causes a higher variation of the torso’s position due to the variation of the surface’s position. With the simple control, the torso collided with the ground due to an insufficient height. On the contrary, with the developed posture the torso had a collision-free navigation. Bearing that in mind, it can be concluded that with posture control, the hexapod is less prone to collide with the terrain, and thus the posture becomes stable.

Figure 17. Results obtained for the irregular terrain scenario: (a) Height variation; (b) Pitch velocity variation; and (c) Roll velocity variation.

Considering that the physical prototype of the hexapod robot is currently under development, and it is not possible to examine the designed control method in real-world conditions, the posture control validation follows the Zero-Moment-Point (ZMP) approach, that is commonly found in the literature [Reference Deng, Xin, Zahng and Mistry52–Reference Asif and Iqbal54]. The presented methodology evaluates the locomotion stability by measuring the distance between the ZMP projection and the edges of the support polygon in each time step. If the ZMP coordinates are close to the polygon’s center, the robot is considered stable. Since the equations of motion of the Gazebo software are modeled at the velocity level, the ZMP cannot be calculated using the linear acceleration of the torso. Thus, for the computational simulations, the ZMP is estimated using the following expression [Reference Deng, Xin, Zahng and Mistry52, Reference You, Fan and Liu53]:

(32) \begin{equation} \mathbf{p}_{\mathrm{ZMP}}=\left\{\frac{\sum\limits_{i=1}^{n}f_{\mathrm{n},i}x_{i}^{\mathrm{p}}}{\sum\limits_{i=1}^{n}f_{\mathrm{n},i}},\frac{\sum\limits_{i=1}^{n}f_{\mathrm{n},i}y_{i}^{\mathrm{p}}}{\sum\limits_{i=1}^{n}f_{\mathrm{n},i}}\right\} \end{equation}

where n represents the number of limbs, f _n,i denotes the normal contact force of the ith limb, and $\{x_{i}^{\mathrm{p}},\mathrm{y}_{i}^{\mathrm{p}}\}$ are the relative of each foot. The stability margin is given by [Reference Deng, Xin, Zahng and Mistry52]:

(33) \begin{equation} s_{\mathrm{ZMP}}=\min\!\left(d_{1},\ldots,\,d_{\mathrm{e}}\right) \end{equation}

where d_i represents the distance of the ZMP toward the edges of the support polygon and e denotes the number of edges of the support polygon. The hexapod robot stability improves when the quantity of s _ZMP increases. For the posture control validation, the deviation of the ZMP and the support polygon’s center and the stability margin are compared against a simple locomotion control.

In Fig. 18(a), the variance of the displacement between the ZMP and the support polygon’s center is presented for the ramp environment. With the implementation of the enhanced posture control, there is a significant decrease in the interquartile range of the ZMP distance to the polygon’s center, which supports the stability improvement of the locomotion. In comparison to the simple control, the variation of the median quantity is considered negligible. The adjustment of the torso’s pitch and roll angles affects the final coordinates of the feet, which causes the amplitude growth observed in the posture control. In terms of stability margin, the parameter’s median improves, and the data oscillation decreases (see Fig. 18(b)). Subsequently, the ZMP is within the support polygon, and it is considered that the hexapod robot improves its stability with the developed posture control.

Figure 18. Results obtained for the stability analysis in the ramp scenario: (a) Displacement of the ZMP with respect to the geometric center of the support polygon and (b) Stability margin.

In the testbed environment, the terrain topology causes higher locomotion instability when compared to the ramp scenario. The distance between the ZMP and the polygon’s center in each time step is portrayed in Fig. 19(a). Despite the augmentation of torso’s instability regarding the roll orientation, the implementation of the enhanced posture control reduces drastically the variance of the distance between the ZMP and the polygon’s center. The fact that the torso’s position is adjusted to the terrain topology affects the force distribution of each foot, which supports the results observed. Similarly to the ramp environment, the posture control decreases the interquartile range of the stability margin and the median quantity increases. From the two parameters examined, there is a strong improvement in the hexapod robot’s stability.

Figure 19. Results obtained for the stability analysis in the testbed scenario: (a) Displacement of the ZMP with respect to the geometric center of the support polygon and (b) Stability margin.