*Article* **Method of Changing Running Direction of Cheetah-Inspired Quadruped Robot**

**Meng Ning 1, Jun Yang 1, Ziqiang Zhang 2,\* , Jun Li 2, Zhi Wang 2, Longxing Wei <sup>1</sup> and Pengjin Feng <sup>1</sup>**


**Abstract:** The rapid change of motion direction during running is beneficial to improving the movement flexibility of the quadruped robot, which is of great relevance to its research. How to make the robot change its motion direction during running and achieve good dynamic stability is a problem to be solved. In this paper, a method to change the running direction of the cheetah-inspired quadruped robot is proposed. Based on the analysis of the running of the cheetah, a dynamic model of the quadruped robot is established, and a two-level stability index system, including a minimum index system and a range index system, is proposed. On this basis, the objective function based on the stability index system and optimization variables, including leg landing points, trunk movement trajectory, and posture change rule, are determined. Through these constraints, the direction changes with good dynamic stability of the cheetah-inspired quadruped robot during running is realized by controlling the leg parameters. The robot will not roll over during high-speed movement. Finally, the correctness of the proposed method is proven by simulation. This paper provides a theoretical basis for the quadruped robot's rapid change of direction in running.

**Keywords:** quadruped robot; change of running direction; dynamic model; stability index system; simulation analysis

#### **1. Introduction**

Most quadrupeds have the ability to run fast. For example, the cheetah is the fastestrunning land animal in the world, and its speed can reach 104.4 km/h [1]. Antilocapra americana can run very fast, up to 100 km/h, and has good endurance [2]. In particular, to catch a fast-moving target or escape quickly, the running direction of the creature is not constant, thus its running is no longer a plane motion but a motion in a 3D space. Therefore, for the quadruped robot, how to achieve a rapid change of motion direction in running is a problem to be solved [3].

Many researchers have studied the movement mechanism of the quadruped during running [4]. For example, Kamimura et al. [5] hypothesized that the three characteristics of the small vertical movement of their center of mass, small whole-body pitching movement, and large spine bending movement enhance the running ability of the cheetah. The hypothesis was then verified by a model with a spine joint and a torsional spring. In addition, the running of bipedal creatures, such as birds [6,7] or humans [8,9], has also been studied. On this basis, many researchers have studied the running of bio-inspired quadruped robots. To make the quadruped robot have good dynamic performance in running, the research mainly focuses on structural design [10,11], a control algorithm based on the dynamic model [12–14], an energy transfer mechanism [15], and environmental adaptability [16–18]. In terms of prototype, the most representative quadruped robot with running ability is the Cheetah robot developed by the Massachusetts Institute of Technology (MIT). Based on the research on the design principles for highly efficient legged robots and

**Citation:** Ning, M.; Yang, J.; Zhang, Z.; Li, J.; Wang, Z.; Wei, L.; Feng, P. Method of Changing Running Direction of Cheetah-Inspired Quadruped Robot. *Sensors* **2022**, *22*, 9601. https://doi.org/10.3390/ s22249601

Academic Editors: Luige Vladareanu, Hongnian Yu, Hongbo Wang and Yongfei Feng

Received: 6 November 2022 Accepted: 4 December 2022 Published: 7 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

hierarchical controllers, the running speed of the Cheetah robot can reach 6 m/s and has good dynamic performance [19–22]. In addition to quadruped robots, many researchers have studied the running of biped robots [23,24] and hexapod robots [25,26], and have achieved good research results.

Creatures often do not run in only one direction, and they have the ability to move at high speed in 3D space [27,28]. The research on the motion abilities of bio-inspired robots has also expanded from plane motion to 3D space [29,30]. For examples, Di Carlo et al. [22] presented the implementation of model predictive control (MPC) to determine ground reaction forces for a torque-controlled quadruped robot, and the developed MIT Cheetah 3 can realize a full 3D gallop. Sullivan et al. [31] studied the effects of varying step width on the 3D running stability of a bipedal amputee-inspired robot. The research results showed that to obtain narrower step widths, as seen in human locomotion, a roll and yaw control would be needed. In addition to controlling the motion parameters of the robot itself, the robot can change motion direction during running by using auxiliary mechanisms. For example, Kim et al. [32] were inspired by a basilisk lizard's ability to run and steer on water surfaces for a hexapedal robot, which can steer on water by rotating its tail, and the controlled steering locomotion was stable. Kohut et al. [33] presented a running robot that used aerodynamic forces to turn. The research results showed that the robot is capable of stably turning in a 1.2 m radius at 1.6 ms, and the aerodynamic steering is superior for high-speed turns at high forward velocity. In particular, jumping is also a high-speed movement. Many researchers have studied the structure [34–36] and control algorithms [37,38] of the robot so that it can achieve fast steering when jumping.

For quadruped robots with running ability, the existing research mainly focuses on running in a plane. Research on the high-speed motion mechanism of the quadruped robot in 3D space is relatively rare. The difficulty of research on the 3D running of robots is mainly reflected in two aspects: the motion of the robot in 3D space involves many dynamic performance indices and variables to be optimized, and the coupling degree between them is high [39]; conversely, the change of direction in the high-speed motion of the robot can easily cause sudden changes in performance indices [40]. Guaranteeing the stable high-speed movement of the robot is difficult. To make the robot achieve good dynamic stability in high-speed steering, taking the steering running of the cheetah as a reference, a method of changing the running direction of a bio-inspired robot is proposed in this paper. A two-level stability index system, including minimum index system and range index system, is established based on the dynamic model of the robot, and the optimization variables, including leg landing points, trunk movement trajectory, and posture change rule, are determined. Then, the optimal leg input parameters can be obtained based on the improved bee colony algorithm. The analysis results show that the robot can turn quickly while running and has good dynamic stability by using the proposed method.

The remainder of the paper is structured as follows. Section 2 establishes the dynamic model and stability index system, and presents the optimization method of leg parameters. Section 3 shows examples to illustrate the feasibility of the method. Finally, Section 4 discusses the results. This paper provides a theoretical basis for the realization of rapid steering in the running of the quadruped robot.

#### **2. Methods**

#### *2.1. Research Objectives*

Cheetahs often need to change movement direction frequently during hunting. When the cheetah runs in a plane, the angle of leg adduction/abduction is almost zero. When the cheetah needs to change the motion direction during high-speed movement, the adduction/abduction angle is large. Figure 1a shows *θ*<sup>1</sup> and *θ*<sup>2</sup> are the angles between the leg and the vertical direction in the front view, which are 36.6◦ and 64.2◦, respectively, in the illustrated state. With the cooperation of muscle-driving forces, the cheetah can realize steering movement during running. In particular, for the movement gait, cheetahs use a

rotatory gallop with the footfall order of right fore, left fore, left hind, and right hind during curve running [41].

**Figure 1.** (**a**) Moment when the cheetah changes its movement direction during running. The legs of the cheetah have larger adduction/abduction angles than the cheetah moving in one plane. (**b**) 3D model of the cheetah-inspired quadruped robot. The hip joint has 2 DOFs, and the knee joint has 1 DOF.

The 3D model of the cheetah-inspired quadruped robot is shown in Figure 1b. The hip joints of each leg have 2 degrees of freedoms (DOF) for adduction/abduction and flexion/extension movements. The axes of the two hinges intersect. The knee joint has 1 DOF for flexion/extension movement. In addition, the leg is in point contact with the ground, which can be equivalent to a 3-DOF ball pair. At this time, each leg has 6 DOFs and no constraints on the trunk. By controlling the leg postures and driving forces of the cheetah-inspired robot, the robot can simulate the cheetah to change motion direction quickly during high-speed motion.

#### *2.2. Establishment of Dynamic Model*

During running, the cheetah's two forelegs land first and its two hindlegs land later [41]. For simplicity, the two forelegs are assumed to land simultaneously. When the trunk moves to the lowest point, both forelegs leave the ground at the same time, and both hindlegs land. The trajectory of the trunk during the turning of the robot in the leg landing phase is shown in Figure 2a. When the forelegs of the robot land on the ground, its trunk moves along the trajectory *O*1*O*<sup>2</sup> (*O*<sup>1</sup> is the position of the center of mass of the trunk at the moment when the forelegs land, and *O*<sup>2</sup> is the lowest point of the trunk). When the center of mass of the trunk reaches *O*2, the movement direction of the trunk is changed, and the robot leaves the ground along the trajectory *O*2*O*<sup>3</sup> (*O*<sup>3</sup> is the position of center of mass of the trunk at the moment when the hindlegs leave the ground). In the following text, the "trunk descending phase" and "trunk ascending phase" refer to the above two processes. Trajectories *O*1*O*<sup>2</sup> and *O*2*O*<sup>3</sup> are not coplanar, and the trajectory is not necessarily a straight line. The mechanism diagram of the cheetah-inspired quadruped robot is shown in Figure 2b. The coordinate origin of the fixed coordinate system *O*0–*X*0*Y*0*Z*<sup>0</sup> coincides with the projection point of the lowest point *O*<sup>2</sup> of the trunk motion trajectory on the ground. The directions of the coordinate axes are shown in Figure 2b. The coordinate origin of the moving coordinate system *Ot–XtYtZt* coincides with the geometric center of the trunk, and the friction between the legs and the ground during the movement is ignored.

**Figure 2.** (**a**) Trajectory of the trunk of the cheetah-inspired quadruped robot in the process of changing the motion direction in the leg landing phase. The trajectory *O*1*O*<sup>2</sup> of the center of mass of the trunk in the descending phase and the trajectory *O*2*O*<sup>3</sup> in the ascending phase are not coplanar. (**b**) Mechanism diagram of the cheetah-inspired quadruped robot.

First, the motion law of the trunk must be determined. Trajectories *O*1*O*<sup>2</sup> and *O*2*O*<sup>3</sup> can be expressed as

$$\mathcal{S}\_t(\mathbf{x}, \mathbf{y}, \mathbf{z}) = \mathcal{S}\_t \left( \sum\_{q=0}^n d\_{qi} t^q \, \_\prime \sum\_{q=0}^n e\_{qi} t^q \, \_\prime \sum\_{q=0}^n f\_{qi} t^q \right) , \; i = 1, \; 2 \tag{1}$$

where *t* is the time, *aq*, *bq*, and *cq* are the polynomial coefficients, and (*x*1, *y*1, *z*1) and (*x*2, *y*2, *z*2) are the coordinates of points *O*<sup>1</sup> and *O*2, respectively. Equation (1) can also reflect the change of velocity and acceleration of the trunk by derivation. In particular, some boundary conditions are known. For example, at the moment when the forelegs of the robot land on the ground and the trunk reaches the lowest point *O*2, the position and the velocity of the center of mass of the trunk are known. The velocity of the trunk at point *O*<sup>1</sup> is also known according to the motion parameters of the previous cycle (a running cycle is defined from the moment of robot landing to the next landing moment). Therefore, a coupling relationship between the polynomial coefficients may exist.

During the movement of the trunk along trajectories *O*1*O*<sup>2</sup> and *O*2*O*3, the posture of the trunk can be represented by the *ZYX* Euler angle.

$$\mathbf{Q} = \begin{pmatrix} \csc\beta & \cos\beta s\gamma - \sec\gamma & \cos\beta c\gamma + \sin\gamma & \mathbf{x} \\ \sec\beta & \cos\beta s\gamma + \csc\gamma & \cos\beta c\gamma - \cos\gamma & y \\ -s\beta & c\beta s\gamma & c\beta c\gamma & z \\ 0 & 0 & 0 & 1 \end{pmatrix} \tag{2}$$

where "*s*" and "*c*" refer to "sin" and "cos", respectively, and *α*, *β*, and *γ* are the Euler angles of the motion coordinate system *Ot*–*XtYtZt* relative to the fixed coordinate system *O*0–*X*0*Y*0*Z*0. The change of trunk posture can be expressed as a polynomial function:

$$\Phi\_t(\boldsymbol{a}, \boldsymbol{\beta}, \gamma) = \Phi\_t\left(\sum\_{q=0}^n a\_{q\bar{i}} t^q, \sum\_{q=0}^n b\_{q\bar{i}} t^q, \sum\_{q=0}^n c\_{q\bar{i}} t^q\right) i = 1, 2 \tag{3}$$

Similarly, boundary conditions can be set according to actual requirements to reduce the number of polynomial coefficients to be optimized.

Φ Φ

If the position and the posture of the trunk are determined, the leg posture can be reflected through the contact points between the legs and the ground. Taking the foreleg landing as an example, the contact points can be determined by four parameters, *h*1, *h*2, *dk*, and *ϕk*, as shown in Figure 2b. *h*<sup>1</sup> and *h*<sup>2</sup> are the distances from the two landing points to the intersection *S* of *A*3*A*<sup>4</sup> (*A*<sup>3</sup> and *A*<sup>4</sup> are the two landing points of forelegs) and the *Z*<sup>0</sup> axis, respectively. *dk* is the distance from the coordinate origin of *O*0–*X*0*Y*0*Z*<sup>0</sup> to *S*. *ϕ<sup>k</sup>* is the angle between *A*3*A*<sup>4</sup> and the positive direction of the *Z*<sup>0</sup> axis. At this time, the coordinates of the two landing points *A*<sup>3</sup> and *A*<sup>4</sup> can be expressed as

$$A\_i = \begin{bmatrix} (-1)^i h\_i \sin \varphi\_k & 0 & d\_k + (-1)^i h\_i \cos \varphi\_k \end{bmatrix} i = 3, 4 \tag{4}$$

The calculation method of the landing points is the same for the case of hindleg landing. The kinematic equation of the *i*-th leg can be expressed as

$$\,^0T\_t = \,^0T\_{i1}\,^{i1}T\_{i2}\,^{i2}T\_{i3}\,^{i3}T\_{ti\prime} \tag{5}$$

where

<sup>0</sup>*Ti*<sup>1</sup> = ⎛ ⎜⎜⎝ *cθi*1*cθi*<sup>2</sup> *cθi*1*sθi*2*sθi*<sup>3</sup> − *sθi*1*cθi*<sup>3</sup> *cθi*1*sθi*2*cθi*<sup>3</sup> + *sθi*1*sθi*<sup>3</sup> *li*<sup>1</sup> + *li*3*cθi*1*cθi*<sup>2</sup> *sθi*1*cθi*<sup>2</sup> *sθi*1*sθi*2*sθi*<sup>3</sup> + *cθi*1*cθi*<sup>3</sup> *sθi*1*sθi*2*cθi*<sup>3</sup> − *cθi*1*sθi*<sup>3</sup> *li*3*sθi*1*cθi*<sup>2</sup> −*sθi*<sup>2</sup> *cθi*2*sθi*<sup>3</sup> *cθi*2*cθi*<sup>3</sup> *li*<sup>2</sup> − *li*3*sθi*<sup>2</sup> 00 0 1 ⎞ ⎟⎟⎠, *<sup>i</sup>*<sup>1</sup>*Ti*<sup>2</sup> = ⎛ ⎜⎜⎜⎜⎝ *sθi*<sup>4</sup> −*sθi*<sup>4</sup> 0 *li*4*cθi*<sup>4</sup> *sθi*<sup>4</sup> *cθi*<sup>4</sup> 0 *li*4*sθi*<sup>4</sup> 0 010 0 001 ⎞ ⎟⎟⎟⎟⎠ , *<sup>i</sup>*<sup>2</sup>*Ti*<sup>3</sup> = ⎛ ⎜⎜⎜⎜⎝ *cθi*<sup>5</sup> 0 *sθi*<sup>5</sup> 0 *sθi*<sup>5</sup> 0 −*cθi*<sup>5</sup> 0 01 0 0 00 0 1 ⎞ ⎟⎟⎟⎟⎠ , *<sup>i</sup>*<sup>3</sup>*Tti* = ⎛ ⎜⎜⎜⎜⎝ *cθi*<sup>6</sup> −*sθi*<sup>6</sup> 0 *li*5*cθi*<sup>6</sup> *sθi*<sup>6</sup> *cθi*<sup>6</sup> 0 *li*5*sθi*<sup>6</sup> 0 01 *li*<sup>6</sup> 0 001 ⎞ ⎟⎟⎟⎟⎠

*θi*1, *θi*2, and *θi*<sup>3</sup> are the rotation angles of the ball pair. *θi*<sup>4</sup> is the rotation angle of the knee joint. *θi*<sup>5</sup> and *θi*<sup>6</sup> are the rotation angles of hinges 2 and 1, respectively, as shown in Figure 1b. (*ai*1, *bi*1) is the position vector of point *Ai* in the fixed coordinate system *O*0–*X*0*Y*0*Z*0, *li*<sup>1</sup> is the length of link *AiBi*, *li*<sup>2</sup> is the length of link *BiCi*, and (*ai*2, *bi2*) is the position vector of point *Ci* in the moving coordinate system *Ot*–*XtYtZt*. For Equation (5), when the position and the posture of the trunk and the position of the landing points are determined, the joint angle *θij* (*j* = 1, 2, ... , 6) can be obtained by numerical solution. At this time, the position vector of any point on the link can be expressed as

$$r\_{ip} = r\_{ipx}\stackrel{\rightarrow}{i} + r\_{ipy}\stackrel{\rightarrow}{j} + r\_{ipz}\stackrel{\rightarrow}{k}\tag{6}$$

On the basis of solving the kinematics, the dynamic model of the robot should be established. By calculating the first and second derivatives of Equation (5), the angular velocities and angular accelerations of the joints can be expressed as

$$\begin{cases} \dot{\theta}\_{ij} = f\_1(\mathbf{V}\_t, \mathbf{W}\_t) | \mathbf{V}\_t = \left( v\_{t\mathbf{x}\prime} v\_{t\mathbf{y}\prime} v\_{t\mathbf{z}\prime} \right), \mathbf{W}\_t = \left( \dot{a}\_\prime \dot{\beta}\_\prime, \dot{\gamma} \right) \\\ \ddot{\theta}\_{ij} = f\_1(\mathbf{V}\_t, \mathbf{W}\_t, \mathbf{A}\_{t\prime} \mathbf{T}\_t) | \mathbf{A}\_t = \left( a\_{t\mathbf{x}\prime} a\_{t\mathbf{y}\prime} a\_{t\mathbf{z}\prime} \right), \mathbf{T}\_t = \left( \ddot{a}\_\prime \ddot{\beta}\_\prime \dot{\gamma} \right) \end{cases} \tag{7}$$

where *Vt* and *Wt* are the velocity and the angular velocity of the trunk, respectively; *At* and *Tt* are the acceleration and the angular acceleration of the trunk, respectively.

By calculating the first and second derivatives of Equation (6), the velocity and the acceleration of the joint points and centers of mass of the links can be obtained as

$$\mathbf{V}\_{i\bar{j}} = \sum\_{j=1}^{6} \left( a\_{i\bar{j}} \theta\_{i\bar{j}} + b\_{i\bar{j}} \dot{\theta}\_{i\bar{j}} \right) = \sum\_{j=1}^{6} \left( a\_{i\bar{j}} \mathbf{F}(\mathbf{a}\_{\prime} \boldsymbol{\beta}\_{\prime} \boldsymbol{\gamma}\_{\prime} \mathbf{x}, \mathbf{y}, z) + b\_{i\bar{j}} \dot{\mathbf{F}}(\mathbf{a}\_{\prime} \boldsymbol{\beta}\_{\prime} \boldsymbol{\gamma}\_{\prime} \mathbf{x}, \mathbf{y}, z) \right) \tag{8}$$

$$\mathbf{A}\_{i\bar{j}} = \sum\_{j=1}^{6} \left( c\_{i\bar{j}} \theta\_{i\bar{j}} + d\_{i\bar{j}} \dot{\theta}\_{i\bar{j}} + e\_{i\bar{j}} \ddot{\theta}\_{i\bar{j}} \right) = \sum\_{j=1}^{6} \left( c\_{i\bar{j}} F(a, \theta, \gamma\_r x, y, z) + d\_{i\bar{j}} \dot{F}(a, \theta, \gamma\_r x, y, z) + e\_{i\bar{j}} \ddot{F}(a, \theta, \gamma\_r x, y, z) \right) \tag{9}$$

On this basis, the angular velocity and the angular acceleration of each link can be obtained by

$$V\_{i,j+1} = V\_{i,j} + \omega\_i \times L\_i \tag{10}$$

$$
\dot{V}\_{i,j+1} = \dot{V}\_{i,j} + \dot{\omega}\_i \times L\_i + \omega\_i \times (\omega\_i \times L\_i) \tag{11}
$$

where *<sup>V</sup>i*,*<sup>j</sup>* and . *Vi*,*<sup>j</sup>* are the velocity and the acceleration of the *j*-th joint of the *i*-th link, respectively. *<sup>V</sup>i*,*j*+<sup>1</sup> and . *Vi*,*j*+<sup>1</sup> are the velocity and the acceleration of the (*j* + 1)-th joint of the *<sup>i</sup>*-th link, respectively. *<sup>ω</sup><sup>i</sup>* and . *ω<sup>i</sup>* are the angular velocity and the angular acceleration of the *i*-th link, respectively. *L<sup>i</sup>* is the direction vector of the *i*-th link.

The driving torques can be obtained by establishing the Lagrange dynamic equation. The total kinetic energy of the robot can be expressed as

$$E\_k = \sum\_{l=1}^2 \sum\_{j=1}^2 \left(\frac{1}{2} \,^0V\_{ij} \,^0m\_{ij} \,^0V\_{ij} + \frac{1}{2} \,^0\omega\_{ij} \,^T \,^0I\_{ij} \,^0\omega\_{ij}\right) + \left(\frac{1}{2} \,^0V\_t \,^T m\_t \,^0V\_t + \frac{1}{2} \,^0\omega\_t \,^T \,^0I\_t \,^0\omega\_t\right), \tag{12}$$

where *mij* and *mt* are the masses of the *j*-th link of the *i*-th leg and the trunk, respectively. <sup>0</sup>*Vij* and <sup>0</sup>*ωij* are the velocity and the angular velocity of the *j*-th link of the *i*-th leg in the fixed coordinate system, respectively. <sup>0</sup>*Vt* and <sup>0</sup>*ω<sup>t</sup>* are the velocity and the angular velocity of the trunk in the fixed coordinate system, respectively. <sup>0</sup>*Iij* and <sup>0</sup>*It* are the moment of inertia in the fixed coordinate system, which can be expressed as

$$\prescript{0}{}{I}\_{\vec{ij}(t)} = \prescript{0}{}{T}\_{\vec{ij}(t)} \prescript{ij}{}{I}\_{\vec{ij}(t)} \prescript{0}{}{T}\_{\vec{ij}(t)} \prescript{T}{}{\prime} \tag{13}$$

where <sup>0</sup>*Tij*(*t*) is the transformation matrix of the *j*-th link of the *i*-th leg (or the trunk) in the fixed coordinate system. The total potential energy of the robot can be expressed as

$$E\_p = \sum\_{i=1}^{2} \sum\_{j=1}^{2} \left( m\_{ij} \text{g} h\_{ij} \right) + m\_t \text{g} h\_{t\prime} \tag{14}$$

where *hij* and *ht* are the distances from the center of mass of the link and the trunk to the ground, respectively, which can be obtained by kinematic analysis.

The trunk of the robot has 6 DOFs when two legs land at the same time. The drives are installed at the hip and knee joints of the legs because the ball pair formed by the contact between the legs and the ground is passive motion. At this time, the Lagrange dynamic equation can be written as

$$
\pi\_{ij} = \frac{d}{dt} \frac{\partial E\_k}{\partial \dot{q}\_{ij}} - \frac{\partial E\_k}{\partial q\_{ij}} + \frac{\partial E\_p}{\partial q\_{ij}},\tag{15}
$$

where *qij* <sup>=</sup> (*θi*1, *<sup>θ</sup>i*2, *<sup>θ</sup>i*3, *<sup>θ</sup>i*4, *<sup>θ</sup>i*5, *<sup>θ</sup>i*6), and . *qij* = . *θi*1, . *θi*2, . *θi*3, . *θi*4, . *θi*5, . *θi*6 .

Through the method detailed above, the dynamic model of the cheetah-inspired quadruped robot moving at a high speed in 3D space can be established.

#### *2.3. Establishment of Stability Index System*

On the basis of establishing the dynamic model, the stability index system must be set up so that the robot has good dynamic performance by optimizing the leg parameters. A two-level stability index system is proposed, including a minimum index system and range index system.

The indices contained in the minimum index system should be as small as possible during robot movement. It includes the total inertia moment, the angular velocity of the trunk, the zero moment point (ZMP), and the energy consumption of the robot in a motion cycle.

(1) Total inertia moment. During the high-speed movement of the robot, it should maintain good stability without overturning and rolling over. During the descending and ascending phases of the trunk, the mean and the variance of the total inertia moment and the total inertia moment of the robot at the moment of leaving the ground should be as small as possible. The total inertia moment at the *k*-th time can be expressed as

$$\mathbf{M}\_{lk} = \sum\_{i=1}^{2} \sum\_{j=1}^{2} \left( \mathbf{r}\_{ijk} \times \mathbf{F}\_{ijk} + \mathbf{M}\_{ijk} \right) + \mathbf{r}\_{tk} \times \mathbf{F}\_{tk} + \mathbf{M}\_{tk} \tag{16}$$

where *Fij* and *F<sup>t</sup>* are the inertia forces of the *j*-th link of the *i*-th leg and the trunk, respectively. *rij* and *r<sup>t</sup>* are the vectors of the center of mass of the *j*-th link of the *i*-th leg and the trunk in the fixed coordinate system, respectively. *Mij* and *M<sup>t</sup>* are the inertia moments. The above indices can be expressed as

$$\begin{cases} \quad D = \left| \forall D \left( \mathbf{M}\_{ij(t)} \right) - \forall D \left( \mathbf{M}\_{kp} \right) \right| \\ \quad E = \left| \forall E \left( \mathbf{M}\_{ij(t)} \right) \right| \\ \quad V = End \left( \mathbf{M}\_{ij(t)} \right) \end{cases} \tag{17}$$

where *D*, *E*, and *V* represent the mean, the variance, and the end value, respectively.

(2) Angular velocity of the trunk. A small total inertia moment can make the robot have a small angular acceleration, but further ensuring that the trunk has a small angular velocity at the moment of leaving the ground is still necessary to prevent the robot from turning during a long flight time. The angular velocity of the robot can be obtained according to Equations (10)–(11).

$$g(\omega\_r, \mathfrak{a}\_r) = g\left(\theta\_{ij}, \dot{\Phi}, \ddot{\Phi}\right) \tag{18}$$

(3) ZMP. ZMP can be expressed as [42]

$$\begin{cases} \begin{aligned} \label{10} X\_{ZMP} &= \frac{\displaystyle \sum\_{ij}^{4} m\_{ij} \Big( \bar{\mathbf{y}}\_{ij} + \mathbf{g} \Big) x\_{ij} + m\_{t} \Big( \bar{y}\_{t} + \mathbf{g} \Big) x\_{t} - \sum\_{i=1}^{4} m\_{ij} \bar{x}\_{ij} y\_{i} - m\_{t} \bar{x}\_{t} y\_{t} \\ & \sum\_{i=1}^{4} m\_{i} \Big( \bar{y}\_{ij} + \mathbf{g} \Big) + m\_{t} \Big( \bar{y}\_{t} + \mathbf{g} \Big) \end{aligned} }{\begin{aligned} \label{10} Y\_{ZMP} &= 0 \\ Z\_{ZMP} &= \frac{\displaystyle \sum\_{i=1}^{4} m\_{ij} \Big( \bar{y}\_{ij} + \mathbf{g} \Big) z\_{ij} + m\_{t} \Big( \bar{y}\_{t} + \mathbf{g} \Big) z\_{t} - \sum\_{i=1}^{4} m\_{ij} \bar{x}\_{ij} y\_{ij} - m\_{t} \bar{x}\_{t} y\_{t} \\ & \sum\_{i=1}^{4} m\_{i} \Big( \bar{y}\_{ij} + \mathbf{g} \Big) + m\_{t} \Big( \bar{y}\_{t} + \mathbf{g} \Big) \end{aligned} \tag{19}$$

where (*xij*, *yij*, *zij*) and (*xt*, *yt*, *zt*) are the position coordinates of the center of mass of the *j*-th link of the *i*-th leg and the trunk in the fixed coordinate system, respectively.

(4) Energy consumption. The energy consumption of the robot in high-speed motion should be as small as possible. The total energy consumption can be expressed as

$$C = \int\_0^T \sum\_{i=1}^2 \sum\_{j=1}^3 \left| P\_{ij}(t) \right| dt = \int\_0^T \sum\_{i=1}^2 \sum\_{j=1}^3 \left| \pi\_{ij}(t) \omega\_{ij}(t) \right| dt,\tag{20}$$

where *Pij* is the instantaneous power of the *j*-th joint of the *i*-th leg, *τij* is the joint torque, and *ωij* is the joint angular velocity.

The indices contained in the range index system are considered to meet the requirements within given ranges. This includes the driving torques of the legs and the leg swing angle.

(1) Driving torques. The mean value of the discrete points of the driving torques for different joints should be in a small range. In this way, motors with the same model can be selected, reducing the difficulty of robot prototype development and control. The variances of the joint torque should also be in a small range to ensure the smoothness of torque changes, and prevent excessive torque changes from affecting the service life of the motor. The above indices can be expressed as

$$\begin{cases} \quad D = \left| \forall D \left( \pi\_{ij} \right) - \forall D \left( \pi\_{kp} \right) \right| \\\ E = \left| \forall E \left( \pi\_{ij} \right) \right| \end{cases} \tag{21}$$

where *D* and *E* represent the mean and the variance, respectively.

(2) Leg swing angle. The leg swing angle refers to the angle between the line between the hip joint and the landing point and the vertical direction. If the leg swing angle is too large, the robot easily loses stability due to small friction. Therefore, the leg swing angle should be smaller than the given values. The leg swing angle can be expressed as

$$\mathbf{Y} = \left( \arccos \frac{z(R\_{i6}) - z(R\_{i1})}{y(R\_{i6}) - y(R\_{i1})} \right),$$

where *Rij* is the position vector of the *j*-th joint of the *i*-th leg.

Through the above analysis, a two-level stability index system including the minimum index system and the range index system is established. Among them, the constraints for total inertia moment, angular velocity of the trunk, ZMP, and leg swing angle determine the feasibility of robot motion, and the constraints for energy consumption and driving torques determine the performance advantages of robot long-term movement. The stability index system provides a basis for the subsequent optimization of the motion parameters of the robot.

#### *2.4. Leg parameter Optimization Method*

According to the analysis results of the biological mechanism in Section 2.1, the leg postures and the driving torques of the robot during high-speed movement must be determined, which can be obtained through optimization.

The optimization variables include leg posture parameters and trunk motion parameters. The former includes *h*1*i*, *h*2*i*, *dki*, and *ϕki* (*I* = 1, 2), as shown in Figure 2b. The latter includes the polynomial coefficients shown in Equations (1) and (3). In particular, the polynomial coefficients are different in the descending and ascending phases of the trunk.

The optimization objective function can be expressed as

$$\begin{aligned} Z &= \underset{i=1}{\text{Min}} \sum\_{i=1}^{n} w\_i f\_i(\mathbf{x}\_1, \dots, \mathbf{x}\_n) \\ \text{s.t.} &\Gamma \end{aligned} \tag{23}$$

where

$$\begin{split} f\_1 &= \left( |D|\mathbf{M}\_{\mathcal{R}}|-\min(|D|\mathbf{M}\_{\mathcal{R}}|) \right) / \left( \max(D|\mathbf{M}\_{\mathcal{R}}|) - \min(D|\mathbf{M}\_{\mathcal{R}}|) \right) \\ f\_2 &= \left( |\mathbb{E}|\mathbf{M}\_{\mathcal{R}}| - \min(\mathbb{E}|\mathbf{M}\_{\mathcal{R}}|) \right) / \left( \max(\mathbb{E}|\mathbf{M}\_{\mathcal{R}}|) - \min(\mathbb{E}|\mathbf{M}\_{\mathcal{R}}|) \right) \\ f\_3 &= \left( \mathbb{E}nd(|\mathbf{M}\_{\mathcal{R}}|) - \min(|\mathbf{M}\_{\mathcal{R}}|) \right) / \left( \max(|\mathbf{M}\_{\mathcal{R}}|) - \min(|\mathbf{M}\_{\mathcal{R}}|) \right) \\ f\_4 &= \left( |D(\mathbb{Z}MP, A\_i A\_{i+1})| - \min(\mathbb{D}(\mathbb{Z}MP, A\_i A\_{i+1})) \right) / \max(|D(\mathbb{Z}MP, A\_i A\_{i+1})| - \min(\mathbb{D}(\mathbb{Z}MP, A\_i A\_{i+1})) \\ f\_5 &= \left( |\mathbb{E}(\mathbb{Z}MP, A\_i A\_{i+1})| - \min(\mathbb{E}(\mathbb{Z}MP, A\_i A\_{i+1})) \right) / \max(|\mathbb{E}(\mathbb{Z}MP, A\_i A\_{i+1})| - \min(\mathbb{E}(\mathbb{Z}MP, A\_i A\_{i+1})) \\ f\_6 &= \left( \mathbb{E} - \min \mathbb{C} \right) / (\max \mathbb{C} - \min \mathbb{C}) \end{split}$$

*f* <sup>1</sup> and *f* <sup>2</sup> represent the mean and variance of the discrete points of the total inertia moment of the robot in the descending and ascending phases, respectively. *f* <sup>3</sup> represents the total inertia moment of the robot at the moment when it leaves the ground. *f* <sup>4</sup> and *f* <sup>5</sup> represent the mean and the variance of the distance from the ZMP to the line between the two landing points, respectively. *f* <sup>6</sup> represents the energy consumption. In particular, Equation (3) is derived, and the zero angular velocity of the trunk at the end of the descending and ascending phases can be taken as the boundary condition instead of being listed as the objective function. The relationship between polynomial coefficients and time can be obtained and used as the constraints for optimization. This can ensure that the angular velocity of the trunk of the robot is zero when it leaves the ground, and that the robot has good stability. *fi* (*i* = 1, 2, ... , 6) should be as small as possible, which corresponds to the minimum index system. *wi* is the weight coefficient, which is determined by analytic hierarchy process (AHP). The weight coefficient can be expressed as

$$w\_i = \frac{w\_j^0}{\sum\_{j=1}^n w\_j^0} \,\tag{24}$$

where *w*<sup>0</sup> *<sup>i</sup>* are the values obtained by adding rows after the standardization of the judgment matrix. In particular, the consistency of the judgment matrix must be checked to ensure that the scoring of experts is logical and does not appear contradictory. Consistency index can be expressed as

$$CR = \frac{\lambda\_{\text{max}} - n/n - 1}{RI},\tag{25}$$

where *λ*max is the maximum eigenvalue of the judgement matrix, and *RI* is an average random consistency index, which can be obtained by looking up the table.

For Equation (23), **Γ** is the constraints, which can be expressed as

$$\left\{ \begin{array}{l} \left| \forall D \left( \mathfrak{r}\_{ij} \right) - \forall D \left( \mathfrak{r}\_{kp} \right) \right| \leq Z\_{D} and \left| \forall E \left( \mathfrak{r}\_{ij} \right) \right| \leq Z\_{E} \\ \Psi \leq \Psi\_{o} \end{array} \right. \tag{26}$$

where the first formula indicates that the mean and the variance of the driving torques should meet the range requirements, and *ZD* and *ZE* are the given reference values. The second formula indicates that the angle between the leg and the ground should be less than the given value **Ψ**0. Equation (26) is consistent with the range index system.

For the above optimization variables and objective functions, an improved bee colony algorithm is applied in this paper, and the optimization is shown in Algorithm 1. First, the initial ranges of optimization variables *Q*, the kinematic feasible region *O* (make sure the trunk is in the workspace), and the maximum value *Z*(*i*)*max* and minimum value *Z*(*i*)*min* of the objective function are given. A set of optimization variables *W*(*j*) is taken from the initial range *Qrange*, and each item in the optimization objective function *Z*(*j*) and the motion parameters *O*(*j*) are calculated based on bee colony algorithm *H<sup>G</sup> rule* and range constraints *M*. By judging the value of the objective function, *Qrange*, *Z*(*i*) *min*, and *Z*(*i*) *max* are updated and assigned to *Qrange new* ,*Z*((i)) *min new*, and *Z*(*i*) *max new* , respectively, to obtain the approximate accurate values. While carrying out the accurate dimensionless processing of the objective function, the design efficiency is improved through the accurate constraints of the ranges. On this basis, values are taken from *Qrange new* , the objective function *Z*(*j*) is calculated, and *Qrange new* ,*Z*(*i*) *min new*, and *Z*(*i*) *max new* are simultaneously updated to improve the constraint accuracy continuously. In particular, when the ratio of the total number of cycles to the current number of cycles is a positive natural number, the result obtained by the previous generation calculation is used as a reference to reduce the ranges by multiplying the scale coefficient *k* (*k* < 1) to improve the calculation efficiency. Finally, the minimum value of the objective function is obtained, and the optimization variables are output.


In particular, for the above optimization process, parallel calculation is used in the process of employed bees, on-looker bees, and scout bees to find honey sources, and the extreme value of the objective function is dynamically updated after the calculation for each kind of bee is completed. At the same time, the dynamic parameters for the scout bees are added, and the working threshold of the scout bees is adjusted dynamically according to the convergence of the objective function. The above process can improve the convergence speed and enhances the ability of global optimal search.

#### **3. Results**

#### *3.1. Examples*

To prove the feasibility of the method proposed in this paper, two examples are given. The structural parameters of the robot are shown in Table 1. The variable (*a*2, *b*2) represents the coordinates of the hip joint in the moving coordinate system *Ot*–*XtYtZt*, as shown in Figure 2b. The thigh and the calf legs are cylinders, and the section radius *r* and length *h* are given in Table 1. For example 1, the known parameters are shown in Table 2. *v*<sup>1</sup> and *v*<sup>2</sup> are the velocities of the trunk at points *O*<sup>1</sup> and *O*3, respectively; (*x*1, *y*1, *z*1) and (*x*2, *y*2, *z*2) are the coordinates of points *O*<sup>1</sup> and *O*2, respectively. **Φ**<sup>0</sup> is the trunk posture angle at the moment of landing. **Ψ**0, *ZD*, and *ZE* are the given values shown in Equation (26).

**Table 1.** Structural parameters of the cheetah-inspired quadruped robot.


**Table 2.** Known values during optimization for example 1.


To ensure the motion stability of the robot, the trunk movement trajectory and posture change rule are assumed cubic functions, and the polynomial coefficients need to be determined. In the descending phase of the trunk, the trajectory equation and the posture equation of the trunk have 12 undetermined polynomial coefficients. The position and the velocity of the trunk at *O*<sup>1</sup> and *O*<sup>2</sup> are known, the angle and angular velocity of the trunk at *O*<sup>1</sup> and *O*<sup>2</sup> are known, and the trunk does not rotate around the *Zt* axis. By substituting the boundary conditions into Equations (1) and (3), all polynomial coefficients can be expressed as functions of time. In the ascending phase of the trunk, the position and the velocity of the trunk at *O*2, the velocity direction of the trunk at *O*3, the trunk angle and the angular velocity at *O*2, and the trunk angular velocity at *O*<sup>3</sup> are known. The relationship between the undetermined coefficients and the movement time can also be obtained by substituting the boundary conditions into Equations (1) and (3). However, not all polynomial coefficients can be expressed in time. Two coefficients in Equation (1) and three coefficients in Equation (3) still need to be determined. To sum up, the optimization variables involved in the trunk motion to be determined include *t*<sup>1</sup> (movement time of trunk in the descending phase), *t*<sup>2</sup> (movement time of trunk in the ascending phase), *a*12, *a*22, *d*12, *d*22, and *d*32. The optimization variables also include leg landing point parameters *h*1*i*, *h*2*i*, *dki*, and *ϕki*. The initial ranges of the optimization variables are listed randomly in Table 3. The weight coefficients are *w*<sup>1</sup> = 0.0755, *w*<sup>2</sup> = 0.0464, *w*<sup>3</sup> = 0.5984, *w*<sup>4</sup> = 0.1305, *w*<sup>5</sup> = 0.0623, and *w*<sup>6</sup> = 0.0869. For the hierarchical bee colony algorithm, the number of honey sources is 100, the number of leading bees is 100, and the maximum number of iterations is 100. The optimization results obtained by the method proposed in this paper are shown in Table 3.


**Table 3.** Initial parameter ranges and optimization results for example 1.

Figure 3a,b show the motion sequence when the robot changes motion direction during running. In the descending phase of the trunk, the two forelegs of the robot are in contact with the ground, and the center of mass of the trunk moves 243.72 mm in 0.8 s. In the ascending phase of the trunk, the two hindlegs of the robot are in contact with the ground, and the center of mass of the trunk moves 492.58 mm in 0.3 s. Figure 3c shows the trajectory of the center of mass of the trunk. When the forelegs of the robot touch the ground, the motion direction vector of the trunk is (0, 0, 1). At the moment when the hindlegs of the robot leave the ground, the motion direction vector of the trunk is (1.74, 0.5, 1.84). From the top view, the included angle of the direction vector is 28.08◦, and the running direction of the robot changes clearly. Figure 3d shows the change of trunk posture. The proper change of body posture is conducive to keeping the good dynamic stability of the robot. The angles of the robot around the three axes at the moment of leaving the ground are −13.01◦, 32.69◦, and 19.81◦. The changes of angular velocities obtained by the method described in Section 2.4 are shown in Figure 3e. The angular velocity of the robot at the moment of leaving the ground is zero, and the trunk of the robot will not rotate significantly in the flight phase. Figure 3f shows the positions of the landing points of the legs. In the descending phase of the trunk, the coordinates of the landing point of the two forelegs are (200, 0, and 243.81 mm) and (−62.01, 0, and −107.17 mm). In the ascending phase of the trunk, the landing point coordinates of the two hindlegs are (−151.93, 0, and 0.04 mm) and (−32.88, 0, and 99.98 mm). The leg landing points are no longer symmetrical along the *Z*<sup>0</sup> axis, and the legs have evident adduction/abduction angles. The maximum leg swing angles of the forelegs are 41.04◦ and 41.09◦, and the maximum leg swing angles of the hindlegs are 45.09◦ and 41.00◦. This outcome is consistent with the analysis results of the movement mechanism of the cheetah when it turns during running, as shown in Figure 1a.

The dynamic performance of the cheetah-inspired quadruped robot during steering is shown in Figure 4. The change of the total inertia moment of the robot is shown in Figure 4a. In the descending phase of the trunk, the amplitude of the total inertia moment of the robot is small and changes gently. When the trunk reaches the lowest point *O*2, the total inertia moments of the robot along the three axes are 0.99, −0.12, and −0.73 N·m. In the ascending phase of the trunk, the total inertia moment increases substantially but then decreases rapidly because the robot needs to obtain a large acceleration in a short time. When the trunk reaches point *O*3, the total inertia moments of the robot along the three axes are −4.64, 0.94, and −0.31 N·m. Figure 4b shows the changes of the total inertia moment before and after optimization. "*B*" and "*A*" refer to before and after optimization, respectively. *E*(*MtI*), *m*(*MtI*), and *V*(*MtI*) represent the end value, the mean, and the variance of the total inertia moment, respectively. The total inertia moment before optimization is calculated by substituting the initial parameters. Figure 4b shows that the total inertia moment decreases considerably after optimization. In the descending phase of the trunk, the maximum reductions of *E*(*MtI*), *m*(*MtI*), and *V*(*MtI*) after optimization are 54.22%, 47.79%, and 78.96%, respectively. In the ascending phase of the trunk, the maximum reductions of *E*(*MtI*), *m*(*MtI*), and *V*(*MtI*) after optimization are 99.6%, 97.7%, and 99.9%, respectively. The dynamic stability of the robot is remarkably improved. Figure 4c shows the distance from ZMP to the connecting line between the two landing points during the descending phase of the trunk. The average value of the distance is 9.08 mm. ZMP is near the connecting line of two points. The deviation is small compared with the size of the robot, and the robot has good dynamic stability. Figure 4d shows the mean and variance of driving torques. "J-1," "J-2," and "J-3" represent hinges 1, 2, and 3, respectively, as shown in Figure 1b. The maximum difference between the mean values of driving torques of the different joints is only 2.8 N/m, and the driving torques change smoothly with a slight difference in amplitude. Moreover, the energy consumption of the robot during movement is 27.96 J. The above analysis results reveal that the optimized indices that are contained in the minimum index system are very small, and the indices contained in the range index system are within reasonable ranges. The robot has good dynamic stability by using the parameters of the leg postures and the driving torques obtained by the method proposed in this paper.

For example 1, the robot turns left during running from the top view, thus the projection of the motion direction vector on the ground of the robot at the moment of leaving the ground is counterclockwise relative to that of the robot at the moment of landing. To prove the feasibility of the method proposed in this paper further, an example of the robot turning to the right is given. The known values remain unchanged, as shown in Table 2. The initial range of optimization variables and optimization results are shown in Tables 4 and 5.

**Figure 3.** (**a**) Motion sequence of the robot in the descending phase of the trunk for example 1. (**b**) Motion sequence of the robot in the ascending phase of the trunk for example 1. (**c**) Movement trajectory of the trunk for example 1. In the descending and ascending phases of the trunk, the trajectories are cubic functions. (**d**) Changes of trunk posture for example 1. In the descending phase of the trunk, the angle changes of the trunk around the three axes are 25◦, 10◦, and 0◦. In the ascending phase of the trunk, the angle changes of the trunk around the three axes are 38.01◦, −22.69◦, and −20.00◦. (**e**) Changes of trunk angular velocity for example 1. The angular velocity of the robot at the moment of leaving the ground is zero. (**f**) Landing points of the forelegs and the hindlegs for example 1.

**Figure 4.** (**a**) Change of the total inertia moment of the robot in the descending and the ascending phases after optimization for example 1. (**b**) Changes of the total inertia moment of the robot before and after optimization for example 1. The maximum reduction of the end value, the mean, and the variance of the total inertia moment of the robot after optimization are 99.6%, 97.7%, and 99.9%, respectively, and the decreases are evident. (**c**) Change of ZMP during the descending phase of the trunk for example 1. ZMP changes near the connecting line of the landing points of two legs, showing good stability. (**d**) Mean and variance of the driving torques for example 1. The maximum difference of the mean value is 2.8 N·m, and the maximum variance is 7.16 N·m. The driving torques are within reasonable ranges.

**Table 4.** Known values during optimization for example 2.


**Table 5.** Initial parameter ranges and optimization results for example 2.


The optimization results for example 2 are shown in Figure 5. Figure 5a,b show the motion sequence of the robot for example 2. In the descending phase of the trunk, the center of mass of the trunk moves 243.72 mm in 0.8 s. In the ascending phase of the trunk, the center of mass of the trunk moves 337.3 mm in 0.22 s. Figure 5c shows the change of trunk posture. Similarly, the angular velocity of the trunk at the moment of leaving the ground is zero. Figure 5d shows the positions of the landing points of the legs. The maximum leg swing angles of the forelegs and hindlegs are (53.30◦, 48.10◦) and (33.75◦, 27.71◦), respectively. The dynamic performance of the cheetah-inspired quadruped robot during steering is shown in Figure 5e,f. Figure 5e shows the changes of the total inertia moment before and after optimization. Similarly, the total inertia moment before optimization is calculated by substituting the initial parameters. Compared with those before optimization, in the descending phase of the trunk, the maximum reductions of *E*(*MtI*), *m*(*MtI*), and *V*(*MtI*) after optimization are 45.8%, 46.9%, and 68.8%, respectively. In the ascending phase of the trunk, the maximum reductions of *E*(*MtI*), *m*(*MtI*), and *V*(*MtI*) after optimization are 98.89%, 90.2%, and 97.57%, respectively. Moreover, the total inertia moments of the robot around the three axes at the moment of leaving the ground are −4.7069, −0.5166, and −3.9464 N·m within small ranges. This finding shows that the robot has good dynamic stability. Figure 5f shows the mean and the variance of the driving torques. The maximum difference between the mean values of driving torques of the different joints is only 5.88 N/m, and the variance is not too large.

#### *3.2. Simulation*

Two examples are simulated with Webots to verify that the robot can turn quickly while running and has good dynamic stability, and the simulation videos are shown in Supplementary Materials. Each example contains two continuous running cycles. The structural parameters of the robot are consistent with theoretical calculation. Since the robot has a long flight time during running, the stability of the robot is directly reflected by the rotation angle of its trunk. For simulation example 1, the robot turns left continuously while running. Taking the joint angles and driving torques obtained by theoretical calculation as input for the first running cycle, the input parameters of the robot in the second running cycle can be calculated by the same method. The top view of the motion sequences of the robot in two cycles is shown in Figure 6a. *O*1, *O*<sup>2</sup> and *O*<sup>3</sup> refer to the position of the center of mass of the trunk shown in Figure 2a, and *O*<sup>4</sup> refers to the highest point of the robot in the flight phase. The robot rotates 28.08◦ around the vertical axis in both cycles, and the motion direction changes substantially. The trunk rotation angles corresponding to Figure 6a is shown in Figure 6c. In the descending and ascending phases of the trunk of the robot in the first running cycle, the trunk rotation angles are exactly the same as those shown in Figure 3d. From one perspective, it proves the correctness of theoretical calculation; conversely, it can also show that the robot moves according to the predetermined rules in the descending and ascending phases of trunk, without movement failure, such as rollover. In the flight phase, the trunk posture of the robot almost remains unchanged. The maximum rotation angles of the trunk in two cycles around the three axes are −2.95◦, 3.75◦, and 3.32◦. This finding shows that the robot has small angular velocity at the moment of leaving the ground, which can be seen in Figure 6e, and it also proves the correctness of the proposed two-level stability index system. For simulation example 2, the top views of the motion sequences of the robot and the trunk rotation angle are shown in Figure 6b,d, respectively. The robot rotates −28.08◦ around the vertical axis in both cycles. In the descending and ascending phases of the trunk of the robot in the first cycle, the trunk rotation angles are exactly the same as those shown in Figure 4c. In the flight phase, the maximum rotation angles of the trunk in two cycles around the three axes are 1.88◦, 0.92◦, and −3.72◦. The change of the angular velocity of the trunk corresponding to simulation example 2 is shown in Figure 6f, and the angular velocity of the robot at the moment of leaving the ground in two cycles is approximately zero. The rotation angles of the trunk in the flight phases

and the angular velocity of the trunk at the moment of leaving the ground are within small ranges, and the robot shows good dynamic stability.

**Figure 5.** (**a**) Motion sequence of the robot in the descending phase of the trunk for example 2. (**b**) Motion sequence of the robot in the ascending phase of the trunk for example 2. (**c**) Changes of trunk posture for example 2. In the descending phase of the trunk, the angle changes of the trunk around the three axes are 25◦, −10◦, and 0◦. In the ascending phase of the trunk, the angle changes of the trunk around the three axes are −40.01◦, −12◦, and 12.88◦. (**d**) The landing points of the forelegs and hindlegs for example 2. (**e**) Changes of the total inertia moment of the robot before and after optimization for example 2. The maximum reductions of the end value, the mean, and the variance of the total inertia moment of the robot after optimization are 98.9%, 90.2%, and 97.6%, respectively, and the decreases are evident. (**f**) Mean and variance of the driving torques for example 2. The maximum difference of the mean value is 5.88 N·m, and the maximum variance is 6.12 N·m. The driving torques are within reasonable ranges.

**Figure 6.** *Cont*.

**Figure 6.** (**a**) Top view of the motion sequences of the robot for simulation example 1. The robot rotates 28.08◦ around the vertical axis in both cycles, and the direction of motion changes considerably. (**b**) Top view of the motion sequences of the robot for simulation example 2. The robot rotates −28.28◦ around the vertical axis in both cycles. (**c**) Change of trunk posture for simulation example 1. The changes of the maximum rotation angles of the trunk in the flight phase are −2.95◦, 3.75◦, and 3.32◦. (**d**) The change of trunk posture for simulation example 2. The changes of the maximum rotation angles of the trunk in the flight phase are 1.88◦, 0.92◦, and −3.72◦. (**e**) Change of trunk angular velocity of the trunk for simulation example 1. (**f**) Change of trunk angular velocity of the trunk for simulation example 2.

The above simulation results show that the proposed method in this paper can make the trunk posture of the robot stable and achieve good dynamic stability in high-speed steering motion by controlling the leg postures and the driving torques. The cheetahinspired quadruped robot does not overturn or roll over due to excessive velocity and change of movement direction, so that the movement fails.

#### **4. Discussion**

In this paper, the main research objective is to propose a method to maintain the dynamic stability of the robot during steering running. Therefore, a two-level stability index system, including a minimum index system and a range index system, is proposed based on the dynamic model, and optimization objective functions are established based on the index system. The optimization variables include not only leg posture parameters, but also the trunk movement trajectory and posture parameters. Through the coordination of leg postures and driving torques obtained by the improved bee colony algorithm, the legged robot can achieve good dynamic performance [42–44]. The method proposed in this paper can make the quadruped robot achieve fast steering in running, and the following factors need to be considered.

(1) Changes in trunk posture. Figures 3 and 5 show that the posture of the trunk changes during the descending and ascending phases. If the trunk is forced to remain horizontal without a posture change, although the trunk looks more stable, keeping the total inertia moment within a small range at the moment of the robot leaving the ground is difficult. The robot turns over evidently in the flight phase, which leads to motion failure. Figure 7a shows the change of the total inertia moment of the robot when the trunk is forced horizontally. The total inertia moments of the robot at the moment of leaving the ground are −26.43, 5.41, and 7.40 N·m. They are substantially larger than those shown in Figures 4 and 5. Moreover, the change of the trunk angle should be reasonable. For example, Figure 7b shows a set of calculated results for trunk angle changes. Although the total inertia moment of the robot corresponding to Figure 7b is within a reasonable range and the robot is stable, the pitch angle of the robot at the moment of leaving the ground is 41.28◦, which is not conducive to the stability of the robot in the next cycle. Figures 3d and 5c show that the pitch angle of the trunk is relatively small at the moment of leaving the ground, and this is conducive to the robot maintaining good dynamic stability.

**Figure 7.** (**a**) Change of the total inertia moment when the trunk is forced horizontally. The maximum value of the total inertia moment is 27.97 N·m, which is much greater than the value shown in Figure 4a. (**b**) A possible trunk posture change rule. Although the corresponding total inertia moment is within a reasonable range, the maximum pitching angle of the trunk is 41.28◦, which is not conducive to the stability of the robot.

(2) Coupling of multiple parameters. High-speed motion in 3D space has many dynamic stability indices, and the coupling degree between indices is high. For example, a large coupling relationship exists between the steering angle, the velocity, the total inertia moment of the robot at the moment of leaving the ground, and ZMP. When the velocity of the robot at the moment of leaving the ground increases—that is, when the forward distance of the robot in a cycle increases—or the steering angle is too large, always staying near the connecting line the landing points of the two legs is difficult for the ZMP of the robot and the total inertia moment increases, thus the robot has difficulty maintaining good stability. For example, when the velocity of the trunk increases from 3.5 m/s to 5 m/s at the moment of leaving the ground, the total inertia moment after optimization remarkably increases from 6.16 N·m to 14.21 N·m. Although it can increase the movement time of the trunk in the ascending phase to reduce the total inertia moment, the difficulty of obtaining the optimal solution increases. Therefore, the motion parameters of the robot must be reasonably determined to achieve continuous, stable motion.

(3) Determination of weight coefficients. For the optimization objective function shown in Equation (23), the weight coefficients influence the results. For the examples shown in Section 3.1, the weight coefficients are determined by AHP. *f* <sup>3</sup> and *f* <sup>4</sup> have a great influence on dynamic stability, and their weight coefficients are large; *f* 1, *f* 2, *f* 5, and *f* <sup>6</sup> have minimal influence on dynamic stability, and their weight coefficients are relatively small. If the weight coefficients are changed to *wi*=1/*i*, the optimization results show that the total inertia moments of the robot after optimization at the moment of leaving the ground are −21.31, 4.35, and −5.29 N·m. The dynamic stability clearly deteriorates. Therefore, using experts' experience to determine the importance of the indices to determine the weight coefficients is reasonable.

However, due to the complexity of the optimization objective functions and the large number of optimization variables, the current optimization efficiency cannot meet the real-time requirements. The robot needs to complete the motion planning in advance under the known terrain to achieve complex high-speed movement. In the future, on the basis of the method proposed in this paper, the data set can be established and the method to maintain the dynamic stability of the robot during steering running based on deep neural network can be further proposed. The efficiency of the algorithm will be further improved, making it possible for the robot to complete high-speed steering movement in real time.

#### **5. Conclusions**

The steering of the quadruped robot during high-speed running is of great importance for improving its movement flexibility. However, too many optimization variables, high coupling of multiple performance indices, and high velocity make the research difficult. Therefore, taking the cheetah-inspired quadruped robot as the research object, the method of changing the running direction of the robot was proposed to make the robot turn quickly in the process of high-speed movement and have good dynamic stability. (1) On the basis of establishing the dynamic model of the cheetah-inspired quadruped robot running, a twolevel dynamic stability index system was proposed, including a minimum index system and a range index system, which cover almost all of the indices that affect the dynamic stability of the robot. (2) The optimization objective function based on the dynamic stability index system and optimization variables are determined. Then, the optimal values were obtained based on the improved bee colony algorithm. By controlling the leg posture parameters and the corresponding driving torques, the robot can change the motion direction during high-speed movement. (3) According to the method proposed in this paper, two examples were given: The robot turned 28.08◦ to the left and −28.08◦ to the right during forward running when viewed from the top view. The calculation results showed that the total inertia moment of the robot was in a small reasonable range, and the angular velocity of the robot at the moment of leaving the ground was approximately zero, which proved that the robot had good dynamic stability. The simulation results show that there is no obvious change in the posture of the trunk of the robot during the flight phase, and the robot can land stably, which also proved the correctness of the method. The method proposed in this paper can provide a theoretical basis for the realization of high-speed movement of the robot in 3D space and had good applicability.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/s22249601/s1, Video S1: Simulation results of rapid steering in the running of the quadruped robot.

**Author Contributions:** Methodology, M.N. and Z.Z.; software, J.Y. and Z.W.; validation, J.L. and L.W.; writing, original draft preparation, Z.Z. and P.F.; writing—review and editing, Z.Z. and J.Y.; project administration, M.N.; funding acquisition, Z.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant Number 51805010).

**Data Availability Statement:** The original data contributions presented in the paper are included in the article. Further inquiries can be directed to the corresponding authors.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

