Descriptive parameter optimization method for energy-saving gait planning of biped robots

Abstract

Low energy efficiency is an important factor restricting the application and development of biped robots. In this article, an energy-saving walking control method using descriptive parameter optimization to obtain the gait type with minimum energy consumption is proposed. The algorithm evaluates the energy consumed during walking considering the load torque and angular velocity of all joint actuators of the five-mass simplified model of the robot. A gait database with the step length, gait type and stability margin as the input and gait parameters and energy consumption as the output is constructed, and the gait adapted to the zero-moment point region is dynamically adjusted during walking to realize a compromise between the robot walking stability and energy consumption. In the gait parameter optimization part of the algorithm, a mapping relationship between the descriptive parameters and body trajectory is established. Through parameter sampling and inverse kinematics calculations, seeds are selected from the sample set according to the stability margin, and the gradient descent method of directional acceleration is used to approximate the minimum energy consumption under the descriptive parameters in the neighborhood of the seeds. In the gait synthesis part of the algorithm, according to the given walking task and the energy consumption-related items in the gait database, the walking trajectory with minimum energy consumption is planned. In real-time walking, the database is queried according to the planned step sequence, the gait parameters are obtained, the robot joint movement is controlled, the feedback zero-moment point is calculated from the robot foot pressure, and the database input query is adjusted according to the trajectory deviation to simultaneously achieve walking stability and reduce energy consumption. To determine the effectiveness of the algorithm, dynamic simulation experiments and real robot walking experiments are carried out. The experimental results show that our algorithm has a significant energy-saving effect. The experimental videos are available at https://github.com/xkluzq/biped-robot.

Introduction

Currently, biped robots are being intensely researched because they are expected to be able to simulate humanoid walking, ¹ jumping, ² and even running ³ to replace humans in tasks such as accident rescue ⁴ and medical services. ⁵ After nearly 50 years of unremitting effort, researchers have made remarkable progress in the field of biped robots. ⁶ For example, the robot ASIMO designed by Kanazawa et al. ⁷ can climb vertical stairs, and the robot Atlas controlled by Koolen et al. ⁸ can traverse obstacles.

Biped walking is the most important and fundamental problem in research on biped robot motion. The proposed methods can be divided into three main approaches. One is the gait generation approach based on the external characteristics of human walking, which utilizes high-speed optical motion capture systems to obtain the motion trajectory of human body markers and applies the obtained data to robot walking. ⁹ Another approach is based on the underlying characteristics of human movement, which uses a central pattern generator to generate rhythm signals similar to the control of human walking and adjusts the mass center of the robot body and the online foot trajectory. ¹⁰ The third approach is based on the dynamic characteristics of a biped robot, in which the robot is simplified into a linear inverted pendulum model consisting of a point mass and a massless telescopic leg ¹¹ or a cart-table model consisting of a massless table and a cart driving on the table with all the mass concentrated within it. ¹²

There is still a large deficiency in the structural flexibility of biped robots compared with that of humans, which indicates that gait generation methods based directly on the external or internal characteristics of human walking have limitations. Most researchers prefer to simplify the biped robot model. However, in the model simplification process, some secondary information is ignored, and the accuracy of the simplification determines the upper bound of the control effect. ¹³ Shimmyo et al. ¹⁴ and Faraji et al. ¹⁵ decomposed the mass of a biped robot into the body and legs to form a three-mass inverted pendulum model and used preview control and model predictive control, respectively, for gait pattern generation. Luo and Chen ¹⁶ further proposed a three-mass angular momentum model using model predictive control to obtain better motion control accuracy.

In the gait control of biped robots, some constraint conditions for the robot body are often introduced to reduce the complexity of the algorithm, which causes unnatural walking motions, consumes large amounts of energy, and reduces the possible walking time of biped robots with limited energy storage. ¹⁷ Thus, some researchers are currently working toward reducing constraints and improving the walking efficiency. Hong et al. ¹⁸ relaxed the vertical constraint of the body center of mass (CoM) and effectively alleviated the degree of bending of the knee joint during walking. Shin and Kim ¹⁹ expanded the range of planned trajectories, which improved the walking speed and reduced the walking energy consumption. Li et al. ²⁰ studied a method for reducing energy consumption through trunk pitching motion during robot walking.

Fewer constraints to effectively reduce energy consumption will increase the parameter space in biped robot gait control. One promising way to find the optimal solution is to use artificial intelligence. Dau et al. ²¹ used polynomial interpolation to plan hip and foot trajectories and applied a genetic algorithm to optimize seven key parameters. Elhosseini et al. ²² designed a whale optimization algorithm with a random parameter “a” and a weight parameter “C” to find the optimal setting for the hip parameters. Wu and Li ²³ introduced fuzzy control between robot acceleration, pressure sensors, and gait mode generators, resulting in better responses to external force disturbances. Tao et al. ²⁴ proposed using the particle swarm optimizer of parallel comprehensive learning to select the key parameters affecting gait and realize rapid and stable walking of a robot. Wright and Jordanov ²⁵ summarized the applicability, advantages and disadvantages of common intelligent algorithms in robot motion control.

In summary, in this paper, a novel gait energy consumption optimization algorithm that utilizes descriptive parameters to characterize the three-dimensional motion space of a robot body is proposed. To enhance the optimization precision, a five-mass model is employed, which is located within the body, legs, and feet. The algorithm establishes a mapping relationship between inputs such as the walking step length, type, and stability margin and outputs such as the walking gait and energy consumption. Using this mapping, a robot gait database is constructed for both planning and real-time control of robot walking. The remainder of this paper is organized as follows. In Section 2, we describe the simplified biped robot model and formulate the gait control problem. Section 3 focuses on the gait control algorithm from two aspects: optimal minimization of the energy consumption and real-time motion control. According to the algorithm in this paper and various control methods proposed in the related literature, Section 4 presents the simulation and experimental results for performance analysis of the proposed algorithm. In Section 5, the algorithm is summarized, and possible follow-up work is presented.

Problem formulation

System structure

The biped robot used to verify the algorithm in this paper is shown in Figure 1. Each leg of the robot has five degrees of freedom (DoFs), including hip rolling angle $q_1/q_6$ , pitch angle $q_2/q_7$ , knee pitch angle $q_3/q_8$ , ankle pitch angle $q_4/q_9$ and roll angle $q_5/q_{10}$ , and the actuator vector can be defined as $\mathit{q}=[q_1{q}_2\cdots q_{10}{]}^T$ . The mass of the robot is divided into five CoMs: body b, left leg lk, right leg rk, left foot lf and right foot rf. The position $\left\{{\mathit{r}}_i|i\in P_c\right\}$ of each mass point is marked in Figure 1, where $P_c=\{b,\mathrm{lk},\mathrm{rk},\mathrm{lf},\mathrm{rf}\}$ . In the robot walking analysis, the origin in the global coordinate system is located at the midpoint between the two legs when the robot is upright. The forward, lateral and vertical directions of the robot are defined as the X-axis, Y-axis and Z-axis, respectively. Four FSR pressure sensors are installed on the sole of each foot to measure the force distribution on the feet. Table 1 lists the mass and length of the robot’s body and right leg, and the parameters of the robot left leg are the same as those of the right leg.

Figure 1.

Link structure of the experimental robot.

Table 1.

Basic parameters of the experimental robot.

Symbol	Description	Value
$m_b$	Upper body mass	1.2 kg
$m_{\mathrm{rk}}$	Right leg mass	0.2 kg
$m_{\mathrm{rf}}$	Right foot mass	0.2 kg
$l_b$	Distance between ${\mathit{r}}_b$ and ${\mathit{r}}_h$	0.04 m
$l_h$	Hip width	0.08 m
$l_1$	Link length from $q_1$ to $q_2$	0.04 m
$l_2$	Link length from $q_2$ to $q_3$	0.085 m
$l_3$	Link length from $q_3$ to $q_4$	0.075 m
$l_4$	Link length from $q_4$ to the foot	0.04 m
$l_{\mathrm{fl}}$	Foot length	0.12 m
$l_{\mathrm{fw}}$	Foot width	0.06 m

There are many phases during biped robot walking, such as the single support phase (SSP), double support phase (DSP), and swinging foot landing collision phase. The mathematical model of a walking biped robot is a mixed system with continuous and discrete dynamics. ²⁶ To simplify the analysis of biped robot walking, our research is based on the following assumptions:

• (a) The motion process of the biped robot is continuous, there are no major mutations in the motion state, the force between the links of the robot is limited, and there are no large deformations of the links. Furthermore, the elastic effect of the robot is ignored, the joints between the links are regarded as rigid connections, and the internal clearance is ignored.

• (b) The robot body maintains an upright posture. Most relevant studies show that the pitch angle of the human trunk is generally within 3°, ²⁷ and with forward inclination of the body, the energy consumption in robot walking increases. ²⁸

• (c) The robot feet are always parallel to the ground. Most common biped robots do not have toes, so their feet cannot play the role of toes in improving the driving performance when lifting and landing. ²⁹ In the contact and collision between the foot and the ground, maintaining the foot parallel to the ground helps reduce sliding and bouncing. ³⁰

• (d) The duty ratio of the DSP is $\sigma =25\%$ . The state cycle of biped robot walking is defined as the gait cycle T _s , and the duty ratio $\sigma =T_{\mathrm{DSP}}/T_s$ is used to represent the proportion of the DSP time $T_{\mathrm{DSP}}$ within $T_s$ Through a study of human gait, the range of $\sigma$ was found to be approximately 15%–25%. ²⁷

Allowable zero-moment point (ZMP) region

The ZMP is the most widely used stable walking criterion for biped robots. The ZMP is the point on the ground at which the horizontal component of the total moment generated by gravity and inertia forces is zero.^{3–6,8,10,12–25} Ignoring the influence of the angular acceleration of robot actuators, the calculation of the ZMP position $r_{\mathrm{zmp}}=[x_{\mathrm{zmp}}{y}_{\mathrm{zmp}}0{]}^T$ is shown in equation (1). ³¹

$$\{\begin{array}{c}\multicolumn{1}{c}{x_{\mathrm{zmp}}=\frac{\sum _{i\in {\mathit{P}}_c}{m}_i({\stackrel{··}{z}}_i{x}_i+g{x}_i-z_i{\stackrel{··}{x}}_i)}{\sum _{i\in {\mathit{P}}_c}{m}_i({\stackrel{··}{z}}_i+g)}}\\ \multicolumn{1}{c}{y_{\mathrm{zmp}}=\frac{\sum _{i\in {\mathit{P}}_c}{m}_i({\stackrel{··}{z}}_i{y}_i+g{y}_i-z_i{\stackrel{··}{y}}_i)}{\sum _{i\in {\mathit{P}}_c}{m}_i({\stackrel{··}{z}}_i+g)}}\end{array}$$

(1)

Here $g=9.80m/s^2$ , and ${\mathit{r}}_i=[x_i{y}_i{z}_i{]}^T$ is the position of mass point i.

Robot modeling errors and disturbances in walking will affect the stability of motion. Setting the ZMP trajectory at the centerline of the supporting foot during gait planning can achieve the most stable walking, ¹⁶ but this is not an energy-efficient method. The method of allowable ZMP region (AZR) draws out some edge regions in the region of the support foot to compensate for the disturbance, so that the planned ZMP trajectory is located in the subregion of the support foot to obtain a better trade-off between walking stability and energy consumption. ¹⁹

For a biped robot with foot length $l_{\mathrm{fl}}$ and foot width $l_{\mathrm{fw}}$ , the step length is $s$ , and the Y-axis distance between the two feet is $w$ . During the first DSP in which the left foot is supporting and the right foot is swinging (as an example), AZR is a hexagon with ( ${\mathit{r}}_{a0}$ , ${\mathit{r}}_{a1}$ , …, ${\mathit{r}}_{a5}$ ) which, during SSP, is a square with ( ${\mathit{r}}_{a3}$ , ${\mathit{r}}_{a4}$ , ${\mathit{r}}_{a5}$ , ${\mathit{r}}_{a6}$ ) and, during the second DSP, is a hexagon with ( ${\mathit{r}}_{a5}$ , ${\mathit{r}}_{a6}$ , …, ${\mathit{r}}_{a9}$ ). AZR is shown in Figure 2, where $l_{\mathrm{al}}$ and $l_{\mathrm{aw}}$ are used to represent the length and width of the AZR, respectively. In our gait control algorithm, $l_{\mathrm{al}}$ takes a constant value, dynamically adjusts $l_{\mathrm{aw}}$ to balance the relationship between robot walking stability and energy consumption, and defines the unit value variable $\eta =l_{\mathrm{aw}}/l_{\mathrm{fw}}$ describing the change in dynamic AZR.

Figure 2.

AZR of biped robots.

Energy consumption index

The main component of energy consumption for biped robot motion is the energy consumption $E_m$ of the joint actuators, which is equal to the integral of the instantaneous power vector $p_m(t)$ of all the joint actuators over time t, that is,

$$E_m={\int }_0^t{\mathit{p}}_m\left(\xi \right)d\xi$$

(2)

Through gait trajectory optimization of robot motion, $E_m$ can be significantly reduced. ³² However, for a biped robot with many actuators, $p_m(t)$ is difficult to accurately measure.

${\mathit{p}}_m(t)$ is equal to the sum of the actuator output power ${\mathit{p}}_2(t)$ and additional loss ${\mathit{p}}_{\mathrm{loss}}(t)$ , where ${\mathit{p}}_2(t)$ is the product of the actuator output torque ${\mathit{\tau }}_2(t)$ and angular velocity $\mathit{\omega }(t)$ . ³³ ${\mathit{\tau }}_2(t)$ is mainly used to overcome the load torque ${\mathit{\tau }}_g(t)$ formed by all the CoMs on the rotating shaft of the actuator, where ${\mathit{\tau }}_g\left(t\right)={[{\mathit{\tau }}_1\left(t\right),{\mathit{\tau }}_2\left(t\right),\cdots ,{\mathit{\tau }}_{10}(t)]}^T$ and ${\mathit{\tau }}_i(t)$ can be expressed as:

$${\mathit{\tau }}_i(t)={\mathit{l}}_i^T(t){\mathit{C}}_i(t){\mathit{g}}_m,i=1,2,\cdots ,10$$

(3)

Here ${\mathit{g}}_m={\left[\begin{array}{c}{m}_{b}g{m}_{\mathrm{lk}}g{m}_{\mathrm{rk}}g{m}_{\mathrm{lf}}g{m}_{\mathrm{rf}}g\hfill \end{array}\right]}^T$ is the gravity vector of the CoMs, $l_i(t)$ is the effective distance vector between ${\mathit{g}}_m$ and the i-th rotating shaft, and $C_i(t)$ is the $5\times 5$ diagonal matrix describing the fraction of $g_m$ acting on the rotating shaft.

The sampling period $t_s$ of the gait control system is assumed to be small enough to meet $\mathit{\tau }(n)=\mathit{\tau }(n·t_s)\approx \mathit{\tau }((n-1)·t_s)$ and $\mathit{\omega }(n)\approx [\mathit{q}(n)-\mathit{q}(n-1)]/t_s$ . When the biped robot moves in discrete time N, the energy consumption index E for evaluating the gait algorithm is defined as:

$$E=\frac{1}{t_s}\sum _{n=1}^N{\mathit{\tau }}^T(n)\mathit{\omega }(n)$$

(4)

Problem definition

Based on the above analysis, the energy-saving walking control problem of the biped robot studied in this article is expressed as follows:

• (a) Given the walking step length, gait type and AZR variable of the biped robot, the gait parameters needed to minimize E are determined.

• (b) Given the target distance, the biped walking gait trajectory with minimum energy consumption is planned, and the AZR of the gait is controlled in real time to overcome interference and realize energy-saving walking of the biped robot.

Gait control algorithm

Overview

The biped robot walking control algorithm proposed in this paper takes the gait database as the core, in which the database is constructed through offline gait parameter optimization (GPO), and queries the database according to online gait synthesis (GSYN) to realize stable and energy-saving walking of the robot, as shown in Figure 3.

• (a) The step length s defined in the gait database can be divided into three types: starting step ${\alpha }_0$ , cyclic step ${\alpha }_c$ and stopping step ${\alpha }_e$ . In the GPO algorithm, the body trajectory described by the parameters is determined according to the gait type $\alpha$ , and the low-energy seed set $P_s^{\eta }$ satisfying the $\eta$ condition is selected from the sample subset $R_s$ that satisfies the stability through sampling of the parameter space U and inverse kinematics calculation. Then, in the neighborhood of seed $p\in P_s^{\eta }$ , the minimum E point in the parameter space is iteratively approached according to the gradient descent method, and the minimum $E_s^{\eta }$ and the corresponding actuator angle sequence ${\mathit{g}}_s^{\eta }$ in the result set are stored in the database, that is, the ID_G in the database corresponds to the quintuple $I_G=(s,\alpha ,\eta ,E_s^{\eta },\{g_s^{\eta })$ .

• (b) According to the $(s,\alpha ,\eta ,E_s^{\eta })$ stored in the gait database, the trajectory planning part of the GSYN algorithm calculates the minimum E that meets the walking distance d and stability $\eta$ requirements and then outputs $(s_i,{\alpha }_i)$ in the corresponding step sequence S. The gait controller queries the database according to input $(s_i,{\alpha }_i,{\eta }_i)$ , obtains the actuator angle sequence ${\mathit{g}}_i$ , and successively inputs this sequence into the joint actuator according to discrete time n to drive the biped robot to walk. The ground pressure $F_i$ is collected during walking to calculate the ZMP trajectory ${\mathit{r}}_{\mathrm{zmp}}$ . The preset $\eta$ corresponding to ${\mathit{r}}_{\mathrm{azr}}$ is taken as the target curve of the AZR controller, and the deviation between ${\mathit{r}}_{\mathrm{zmp}}$ and ${\mathit{r}}_{\mathrm{azr}}$ is corrected according to the PD law during robot walking to realize a tradeoff between low energy consumption and robustness.

Figure 3.

Overall gait planning algorithm.

GPO algorithm

The GPO algorithm performs mapping calculation of the data items $(s,\alpha ,\eta )\to (E_s^{\eta },{\mathit{g}}_s^{\eta })$ in the construction of the gait database, which is the basis of our algorithm. The distances between the swinging foot and the supporting foot in the starting and stopping walking processes are not equal, which is different from the cyclic step with the same swinging distances; the gait type α defines these situations. The range of the AZR variable $\eta$ is $[0,1]$ , where a smaller $\eta$ corresponds to higher walking stability and energy consumption. ¹⁹ During robot walking control, s and $\eta$ are discretized. After GPO parameter space sampling, seed neighborhood optimization and other operations, the $E_s^{\eta }$ for trajectory planning and ${\mathit{g}}_s^{\eta }$ for gait control are obtained.

Parametric description of the motion trajectory

When the biped robot walks with step length s, taking the left foot ${\mathit{r}}_{\mathrm{lf}}(n)$ as the supporting foot as an example, the distance between the left hip ${\mathit{r}}_{\mathrm{lh}}(n)$ and ${\mathit{r}}_{\mathrm{lf}}(n)$ is less than the leg length, that is,

$$\max \{\Vert {\mathit{r}}_{\mathrm{lh}}(n)-{\mathit{r}}_{\mathrm{lf}}\Vert {(n)}_2|n=1,2,\cdots ,N\}\le \sum _{j=1}^4{l}_j$$

(5)

If the minimum length $l_{\min }$ of the leg after bending is defined, then the hip joint will be limited to a bounded space ${\Phi }_{\mathrm{lh}}$ during walking, as shown in Figure 4.

Figure 4.

Limited space of hip movement.

In the robot walking process, if the body ${\mathit{r}}_b(n)={\left[\begin{array}{c}{x}_b(n)y_b(n)z_b(n)\hfill \end{array}\right]}^T$ meets constraint (b), then there is a fixed relationship between ${\mathit{r}}_b(n)$ and ${\mathit{r}}_{\mathrm{lh}}(n)$ , and the limited space ${\Phi }_b$ of ${\mathit{r}}_b(n)$ can be obtained by translating ${\Phi }_{\mathrm{lh}}$ . The same is true when the right foot is the supporting foot. In ${\Phi }_b$ , the ${\mathit{r}}_b(n)$ of periodic motion can be expressed by a finite order Fourier series, that is:

$$\{\begin{array}{c}{x}_b\left(n\right)=l_{s}n+\underset{k=1}{{\displaystyle \sum }^{0.5K_x}}{a}_k\sin \left(k{\omega }_{0}n\right)+\underset{k=0.5K_x+1}{{\displaystyle \sum }^{K_x}}{a}_k\cos (\left(k-0.5K_x\right){\omega }_{0}n)\\ y_b\left(n\right)=\underset{k=K_x+1}{{\displaystyle \sum }^{K_x+K_y}}{a}_k\sin ((k-K_x){\omega }_{0}n)\\ z_b\left(n\right)=a_{z0}+\underset{k=K_{z1}+1}{{\displaystyle \sum }^{K_{z2}}}{a}_k\sin \left((k-K_{z1}){\omega }_{0}n\right)+\underset{k=K_{z2}+1}{{\displaystyle \sum }^K}{a}_k\cos (\left(k-K_{z2}\right){\omega }_{0}n)\end{array}$$

(6)

where α = α _c , and $l_s=s/N$ ; when α = α₀ or α _e , $l_s=0.5s/N$ . $N$ is the gait period, $n=1,2,\cdots ,N$ , ${\omega }_0=\pi /N$ , $K_{z1}=K_x+K_y+1$ , $K_{z2}=K_{z1}+0.5(K_z+1)$ , and $K=K_x+K_y+K_z$ . The different trajectories of ${\mathit{r}}_b(n)$ in ${\Phi }_b$ can be described by the parameter vector $\mathit{a}=[a_1{a}_2\cdots a_K]$ . Therefore, the problem of meeting the $(s,\alpha ,\eta )$ condition and extracting the optimal ${\mathit{r}}_b(n)$ trajectory to obtain the minimum E and corresponding $\mathit{g}$ can be transformed into an optimization problem in a -range space U. That is, for $\forall a\in U$ , there is $\mathit{a}\to (E,\mathit{g})$ , and there is an optimal ${\mathit{a}}^{\ast }\to (E^{\ast },{\mathit{g}}^{\ast })$ meeting $E^{*}\le E$ ; $(E^{\ast },{\mathit{g}}^{\ast })$ is $(E_s^{\eta },{\mathit{g}}_s^{\eta })$ in the data item.

Sample calculation

Without loss of generality, the algorithm in this article takes the left foot ${\mathit{r}}_{\mathrm{lf}}(n)=[x_{\mathrm{lf}}(n)y_{\mathrm{lf}}(n)z_{\mathrm{lf}}(n){]}^T$ as the supporting foot and the right foot ${\mathit{r}}_{\mathrm{rf}}(n)=[x_{\mathrm{rf}}(n)y_{\mathrm{rf}}(n)z_{\mathrm{rf}}(n){]}^T$ as the swinging foot during robot walking as an example. If the times of the DSP of the biped robot walking are 1 ∼N₁ and N₂∼N, then the motion trajectories of ${\mathit{r}}_{\mathrm{rf}}(n)$ and ${\mathit{r}}_{\mathrm{lf}}(n)$ are as shown in equations (7) and (8), in which the constants c _{x 1}∼c _{x 5} depend on the gait type α, and their values are defined in Table 2.

$$\{\begin{array}{c}{x}_{rf}(n)=\{\begin{array}{cc}{c}_{x1}s,& 1\le n\le N_1\\ \frac{s}{c_{x3}}\left(c_{x2}\cos \frac{\left(n-N_1\right)\pi }{N_2-N_1}\right),& N_1<n<N_2\\ c_{x4}s,& N_2\le n\le N\end{array}\\ \begin{array}{cc}{y}_{rf}(n)=-\frac{l_h}{2}& \begin{array}{c}1\le n\le N\\ 1\le n\le N_1\end{array}\end{array}\\ \begin{array}{cc}{z}_{rf}(n)=\{\begin{array}{c}0,\\ h_f\sin \frac{\left(n-N_1\right)\pi }{N_2-N_1},\\ 0,\end{array}& \begin{array}{c}{N}_1<n<N_2\\ N_2\le n\le N\end{array}\end{array}\end{array}$$

(7)

$$\{\begin{array}{c}\multicolumn{1}{c}{x_{\mathrm{lf}}(n)=c_{x5}s,1\le n\le N}\\ \multicolumn{1}{c}{y_{\mathrm{lf}}(n)=0.5l_h,1\le n\le N}\\ \multicolumn{1}{c}{z_{\mathrm{lf}}(n)=0,1\le n\le N}\end{array}$$

(8)

where h _f is the maximum height of the robot foot from the ground during walking.

Table 2.

Values of constants c _{x 1}∼c _{x 5} in the planned foot trajectory of the biped robot.

α	c x 1	c x 2	c x 3	c x 4	c x 5
αc	0	1	1	1	0
α0	−0.5	0.5	1	1.5	0.5
αe	−0.5	0.5	2	0.5	0.5

If $a_k\in [d_k,u_k]$ in the parameter set a describing the trajectory of ${\mathit{r}}_b(n)$ is sampled with interval ${\mu }_k$ to form a set $A_k=\{d_k,d_k+{\mu }_k,d_k+2{\mu }_k,\cdots ,u_k\}$ , then the sample space $U_s$ is the Cartesian product of $A_k$ , that is,

$$U_s=A_1\times A_2\times \cdots \times A_K$$

(9)

For $\forall a\in U_s$ , from equation (6), there is $\mathit{a}\leftrightarrow {\mathit{r}}_b(n)$ , and according to constraint (b), the right hip ${\mathit{r}}_{\mathrm{rh}}(n)={\mathit{r}}_b(n)-[00.5l_h{l}_b{]}^T$ . According to constraint (c), there is a geometric relationship between ${\mathit{r}}_{\mathrm{rh}}(n)$ and ${\mathit{r}}_{\mathrm{rf}}(n)$ as shown in Figure 5, based on which the motion angle of the right leg actuator can be directly solved. The inverse kinematics solution process is shown in algorithm 1. Similarly, the algorithm between the left hip ${\mathit{r}}_{\mathrm{lh}}(n)$ and ${\mathit{r}}_{\mathrm{lf}}(n)$ is the same.

Figure 5.

Geometric relations between $r_{\mathrm{rh}}$ and $r_{\mathrm{rf}}$ .

Algorithm 1.

Inverse kinematics algorithm of the right leg.

given ${\mathit{r}}_{\mathrm{rh}}(n)$ , ${\mathit{r}}_{\mathrm{rf}}(n)$

1  ${\mathit{r}}_a(n)\leftarrow {\mathit{r}}_{\mathrm{rf}}(n)+{\left[\begin{array}{c}00l_4\hfill \end{array}\right]}^T$

2  ${\left[\begin{array}{c}{l}_x(n)l_y(n)l_z(n)\hfill \end{array}\right]}^T\leftarrow r_{\mathrm{rh}}(n)-r_a(n)$

3  $q_1(n)\leftarrow \arctan ({\mathrm{l}}_{\mathrm{y}}(\mathrm{n})/l_z(n))$

4  ${\mathit{r}}_c(n)\leftarrow {\mathit{r}}_{\mathrm{rh}}(n)-{\left[\begin{array}{c}0l_1\sin q_1(n)/l_1\cos q_1(n)\hfill \end{array}\right]}^T$

5  $l_0(n)\leftarrow \Vert {\mathit{r}}_c(n)-{\mathit{r}}_a(n){\Vert }_2$

6  ${\theta }_x(n)\leftarrow \arcsin ((x_c(n)-x_a(n))/l_0(n))$

7  ${\theta }_b(n)\leftarrow \arccos ((l_2^2+l_3^2-l_0^2(n))/\left(2l_2{l}_3\right))$

8  ${\theta }_c(n)\leftarrow \arccos ((l_0^2(n)+l_2^2-l_3^2)/(2l_2{l}_0(n)))$

9  $\left[\begin{array}{c}{q}_2(n)q_3(n)\hfill \end{array}\right]\leftarrow \left[\begin{array}{c}{\theta }_c(n)-{\theta }_x(n)\pi -{\theta }_b(n)\hfill \end{array}\right]$

10  $[\begin{array}{c}{q}_4(n)q_5\hfill \end{array}(n)]\leftarrow [\begin{array}{c}{q}_3(n)-q_2(n)-q_1\hfill \end{array}(n)]$

Seed extraction

From algorithm 1 and equations (3) and (4), $\mathit{a}\to (E,\mathit{g})$ and the corresponding ZMP trajectory can be calculated, where $\mathit{g}=[\mathit{q}(1)\mathit{q}(2)\cdots \mathit{q}(N)]$ . Taking the robot shown in Figure 1 as an example, $\mathit{q}=[q_1{q}_2\cdots q_{10}{]}^T$ . Let $R_s\subseteq U_s$ ; then, ${\mathit{a}}_s\in R_s$ meets the following requirements:

• (a) Physical constraint: corresponding to ${\mathit{a}}_s\to (E_s,{\mathit{g}}_s)$ , the maximum angle difference of the actuator between the sampling points in ${\mathit{g}}_s$ is less than the maximum angular velocity ${\omega }_{\max }$ of the actuator times the sampling time $t_s$ , that is:

$$\max \left\{{\Vert {\mathit{q}}_s(n)-{\mathit{q}}_s(n-1)\Vert }_{\infty }|n=1,2,\cdots ,N\right\}\le t_s{\omega }_{\max }$$

(10)

• (b) ZMP condition: $r_{\mathrm{zmp}}(n)$ calculated based on ${\mathit{a}}_s$ is located in the foot support region of the biped robot.

According to the $\eta$ requirements of ${\mathit{r}}_{\mathrm{zmp}}(n)$ , $R_s$ can be divided the subsets $R_s^{\eta }$ , in which ${\mathit{r}}_{\mathrm{zmp}}(n)$ corresponding to ${\mathit{a}}_s^{\eta }\in R_s^{\eta }$ is located in the AZR corresponding to $\eta$ . If the calculated E range contains finite local minimum points, then the seed set $P_s^{\eta }=\{{\mathit{p}}_1,{\mathit{p}}_2,\dots ,{\mathit{p}}_{\Lambda }\}$ is taken, where ${\mathit{p}}_{\lambda }\in P_s^{\eta }$ satisfies the following:

• (a) The number of seeds is less than or equal to the sum of the numbers of values of all the parameters. That is,

$$\Lambda \le n_1+n_2+\cdots +n_K$$

(11)

where $n_k$ is the number of values of $a_k$ in $\forall {\mathit{a}}_s^{\eta }\in R_s^{\eta }$ , with $k=1,2,\cdots ,K$ .

• (b) Seed ${\mathit{p}}_{\lambda }$ has the smallest E among all the elements with the same $a_k$ value. Namely, ${\mathit{p}}_{\lambda }={\mathit{a}}_k^i$ if and only if $E(a_k^i)=\min [E(R_k^i)]$ . $R_k^i\subseteq R_s^{\eta }$ , and the element contains the i-th value of $a_k$ , where $1\le i\le n_k$ .

Parameter optimization

If E satisfies convexity in the neighborhood of ${\mathit{p}}_{\lambda }\in P_s^{\eta }$ , then an extended space $L_{\lambda }^0=B_1^0\times B_2^0\times \cdots \times B_K^0$ is formed with ${\mathit{p}}_{\lambda }=[p_1,p_2,\cdots ,p_k]$ as the center and ${\mathit{\delta }}_{\mathit{\lambda }}=[{\delta }_1,{\delta }_2,\cdots ,{\delta }_k]$ as the interval, where $B_k^0=\{p_k-{\delta }_k,p_k,p_k+{\delta }_k\}$ and $k=1,2,\cdots ,K$ . $\forall \mathit{\xi }\in L_{\mathit{\lambda }}^0$ , $\mathit{\xi }\to (E(\mathit{\xi }),g(\mathit{\xi }))$ is calculated.

If $\exists E_{\min }({\mathit{\xi }}^0)=\min \{E(\mathit{\xi })|\forall \mathit{\xi }\in L_{\mathit{\lambda }}^0\}$ and $E_{\min }({\mathit{\xi }}^0)<E(p_{\lambda })$ in $L_{\lambda }^0$ , then $B_k^1$ is constructed in the direction of ${\xi }_k$ , and the composition space $L_{\lambda }^1=B_1^1\times B_2^1\times \cdots \times B_K^1$ is used to accelerate the optimization calculation. Let

$$\mathrm{\Delta }E({\xi }_k)=\frac{\partial E(\mathit{\xi })}{\partial {\xi }_k}$$

(12)

When $\mathrm{\Delta }E({\xi }_k)<0$ , ${\beta }_k={\xi }_k-p_k$ is taken, and $B_k^1=\{{\xi }_k,{\xi }_k+{\beta }_k,\cdots ,{\xi }_k+\varphi {\beta }_k\}$ . When $\mathrm{\Delta }E({\xi }_k)\ge 0$ , $B_k^1=\left\{{\xi }_k\right\}$ . $\varphi$ is the adjustment factor controlling the calculation scale, which can be adaptively adjusted in the algorithm according to the number of ${\xi }_k$ satisfying the acceleration condition. If ${\mathit{\xi }}^1$ satisfies $E_{\min }({\mathit{\xi }}^1)=\min \{E(\mathit{\xi })|\forall \mathit{\xi }\in L_{\mathit{\lambda }}^1\}$ in $L_{\lambda }^1$ , then a new ${\mathit{p}}_{\lambda }$ is taken for ${\mathit{\xi }}^1$ and the calculation is redone. The specific calculation process is shown in algorithm 2.

Algorithm 2.

Parameter optimization algorithm.

given $P_s^{\eta }$ , ${\mathit{\delta }}_{\mathit{\lambda }}$

1 for $\lambda =1,2,\cdots ,\Lambda$

2  compute ${\mathit{p}}_{\lambda }\to \left(E({\mathit{p}}_{\lambda })g({\mathit{p}}_{\lambda })\right)$ .

3  repeat

4   expand ${\mathit{p}}_{\lambda }$ to form $L_{\lambda }^0$ .

5   for $\forall \mathit{\xi }\in L_{\mathit{\lambda }}^0$ , compute $E(\mathit{\xi })$ .

6   form $L_{\lambda }^1$ according to $\mathrm{\Delta }E({\xi }_k)$ and $\varphi$ . 7   for $\forall \mathit{\xi }\in L_{\mathit{\lambda }}^1$ , compute $E(\mathit{\xi })$ .

8   update ${\mathit{p}}_{\lambda }\leftarrow {\mathit{\xi }}^1$ .

9 until $|E({\mathit{\xi }}^1)-E(p_{\lambda })|\le \epsilon$ .

10  $E(\mathit{p}*)\leftarrow \min \{E({\mathit{p}}_1),E({\mathit{p}}_2),\cdots ,E({\mathit{p}}_{\Lambda })\}$ .

11  $(E_s^{\eta },{\mathit{g}}_s^{\eta })\leftarrow (E(\mathit{p}*),\mathit{g}(\mathit{p}*))$

Given ${\mathit{p}}_{\lambda }$ , in the neighborhood with convexity, the parameter optimization algorithm judges the descending direction of E according to step ${\mathit{\delta }}_{\mathit{\lambda }}$ , accelerates the progression, and must quickly converge. In algorithm 2, ε is the accuracy of the E approximation. After $\forall {\mathit{p}}_{\lambda }\in P_s^{\eta }$ is calculated, the seed $\mathit{p}*$ with minimum E can be taken to obtain the optimized gait ${\mathit{g}}_s^{\eta }$ corresponding to s and η.

GSYN algorithm

Using the gait database, according to the biped robot walking distance d and stability η requirements, the GSYN algorithm realizes trajectory planning of moving targets and real-time gait control during walking.

Trajectory planning

The format of the data stored in the gait database is $I_G$   $=(s,\alpha ,\eta ,E_s^{\eta },{\mathit{g}}_s^{\eta })$ , and the trajectory planning requirements are as follows: satisfy η and minimize E; thus, the main associated data subset in the calculation is $(s,\alpha ,\eta ,E_s^{\eta })$ . If the robot must keep two feet standing side by side at the starting position ${\mathit{r}}_0$ and the stopping position ${\mathit{r}}_e$ , then the planned walking step sequence S is as shown in Figure 6.

Figure 6.

Step sequence of trajectory planning.

In Figure 6, $s_0$ is the starting step, and the energy consumption is $E_0^{\eta }(s_0)$ ; $s_e$ is the stopping step and the energy consumption is $E_e^{\eta }(s_e)$ . $s_1$ to $s_c$ belong to the cyclic step set $S_c$ in the database, where $s_1$ and $s_{c-1}$ are the transition step set to smooth the walking process; $s_i=s_m$ is the intermediate step, and $i=2,\cdots ,c-2$ . If $s_0=s_c=s_e$ , $s_1=s_{c-1}$ , and the length of $d_0=2s_0+2s_1$ is controlled within $[s_m,2s_m]$ , then the calculation method of the step length is as shown in equation (13).

$$\{\begin{array}{c}\multicolumn{1}{c}{c=d/s_m+1}\\ \multicolumn{1}{c}{d_0=d-s_m(c-3)}\\ \multicolumn{1}{c}{s_0=(d_0-s_m)/3}\\ \multicolumn{1}{c}{s_1=(d_0+2s_m)/6}\end{array}$$

(13)

Knowing $s_m$ , if its energy consumption is $E_c^{\eta }(s_m)$ , then the total energy consumption $J_s^{\eta }(s_m)$ in robot walking is:

$$J_s^{\eta }(s_m)=(c-3)E_c^{\eta }(s_m)+E_a^{\eta }(s_m)$$

(14)

where $E_a^{\eta }(s_m)=E_0^{\eta }(s_0)+E_c^{\eta }(s_0)+2E_c^{\eta }(s_1)+E_e^{\eta }(s_e)$ .

Additionally, $\exists s_m^{*}\in S_c$ , such that,

$$J_s^{\eta }(s_m^{*})=\min \{J_s^{\eta }(s_j)|\forall s_j\in S_c\}$$

(15)

Gait control

The gait controller of GSYN queries the gait database according to the input $(s_i,{\alpha }_i,{\eta }_i)$ in step i and sends the obtained actuator angle sequence ${\mathit{g}}_i$ to the joint actuator according to the sampling interval $t_s$ to drive the biped robot to walk. During walking, the movement during each step $T_S=N{t}_s$ is taken as the control cycle; $N_1<n<N_2$ corresponds to the SSP, and $1\le n\le N_1$ and $N_2\le n\le N$ correspond to the DSP, for which different algorithms are applicable. ³⁴ The FSR data $F_i(n)$ generated by plantar pressure can be used to calculate feedback ZMP data ${\mathit{r}}_z(n)={\left[\begin{array}{c}{x}_z(n)y_z(n)0\hfill \end{array}\right]}^T$ , ³⁵ as shown in equation (16).

$$\{\begin{array}{c}{x}_z(n)=\sum _{j=1}^{c_n}{x}_i^j(\mathrm{n})f_i^j(n)/\sum _{j=1}^{c_n}{f}_i^j(n)\hfill \\ \multicolumn{1}{c}{y_z(n)=\sum _{j=1}^{c_n}{y}_i^j(n)f_i^j(n)/\sum _{j=1}^{c_n}{f}_i^j(n)}\end{array}$$

(16)

where for the SSP, $c_n=4$ , and for the DSP, $c_n=6$ .

Generally, during the DSP, the robot is in the overdrive state, and ${\mathit{g}}_i(n)$ easily meets the stability requirements. In AZR control, ${\mathit{r}}_z(n)$ must be located in the AZR. During the SSP, if $x_{\min }$ and $x_{\max }$ represent the minimum and maximum values, respectively, in the X direction in the AZR of the planning for η, then $x_z(N_1+1)\ge x_{\min }$ and $x_z(N_2-1)\le x_{\max }$ must be satisfied. The AZR controller corrects ${\eta }_{i+1}$ according to the deviation $e_i$ between $\{y_z(n)|N_1<n<N_2\}$ and the inner boundary $y_{\min }$ in the Y direction in the AZR of the planning, that is:

$$e_i=\frac{2}{l_{\mathrm{fw}}}\min \{c_o(y_z(n)-y_{\min })|N_1<n<N_2\}$$

(17)

where for $y_{\min }>0$ , $c_o=1$ , whereas for $y_{\min }<0$ , $c_o=-1$ .

To avoid large fluctuations, the AZR controller adopts the incremental method to correct ${\eta }_i$ , that is:

$${\eta }_i={\eta }_{i-1}+\mathrm{\Delta }{\eta }_i$$

(18)

The PD adjustment method of incomplete differentiation is adopted for each step increment $\mathrm{\Delta }{\eta }_i$ and $e_i$ , that is:

$$\mathrm{\Delta }{\eta }_i=K_{\alpha }\mathrm{\Delta }{\eta }_{i-1}+K_1{e}_i+K_0{e}_{i-1}$$

(19)

where, $K_0=-k_P{T}_D/(T_L+T_S)$ and $K_1=k_P[T_D/(T_L+T_S)+1]$ , with $k_P$ being the proportional coefficient, $T_D$ being the differential time constant, and $T_L$ being the inertial time constant.

Experimental evaluation

According to the biped robot structure described in section 2.1, 3D simulation models are constructed in the MATLAB and Webots environments to evaluate our proposed algorithm. The following is defined in the gait database: starting step set $S_0=\{s|s=1,2,\cdots ,6\}$ , cyclic step set $S_c=\{s|s=1,2,\cdots ,13\}$ , stopping step set $S_e=\{s|s=1,2,\cdots ,6\}$ , and AZR set $H=\{\eta |\eta =0.1,0.2,\cdots ,1\}$ . When the GPO algorithm is used to calculate $(E_s^{\eta },{\mathit{g}}_s^{\eta })$ with $(s,\alpha ,\eta )$ as input, the length K of the parameter description is a key parameter. Although in theory, a sufficiently long series can represent the arbitrary periodic motion of the robot body, considering nonlinear factors such as the physical constraints of the motion process and the dead zone of the actuator, a large K value may not be required in practice.

Acquisition of data items

${\mathit{r}}_b(n)$ described by equation (6) and ${\mathit{r}}_{\mathrm{rf}}(n)$ described by equation (8) are used to calculate the data item $(E_s^{\eta },{\mathit{g}}_s^{\eta })$ with cyclic step α _c and $s=10$ , where $K=24$ ( $K_x=8$ , $K_y=7$ , $K_z=9$ ), $N=16$ , and $h_f=1\mathrm{cm}$ . When sampling at the interval $\{{\mu }_k=0.2\mathrm{cm}|k=1,2,\cdots ,K\}$ , considering the monotonic change in ${\mathit{r}}_b(n)$ during walking, the effective number of samples in $U_s$ is 314,510,350,926. After the inverse kinematics calculation described in section 3.2.3 and screening according to physical constraints and ZMP conditions, 44,140,884 elements of $R_s$ are obtained, from which seeds corresponding to different η are selected.

When $\eta =0.8$ , the number of seeds selected is $\Lambda =123$ , as shown in Table 3, where $k$ is the label of $a_k$ , and $n_k^{*}$ is the actual number of seeds obtained by removing the repeatedly selected seeds. $P_k\subset P_{10}^{0.8}$ and contains $n_k^{*}$ seeds and $E_{\min }(P_k)$ is the minimum E of the seeds in $P_k$ . To obtain a gait with lower energy consumption, parameter optimization with algorithm 2 is performed twice with ${\delta }_k^1=0.05\mathrm{cm}$ and ${\delta }_k^2=0.01\mathrm{cm}$ . $E_{\min }^1(P_k)$ and $D_1\%$ , are the minimum value of E and the percentage reduction of E in the first calculation, respectively, while $E_{\min }^2(P_k)$ and $D_2\%$ are those in the second calculations. In the first parameter optimization calculation, $\forall {\mathit{p}}_{\lambda }\in P_{10}^{0.8}$ , $E_{\min }^1(P_{10}^{0.8})=182.0\mathrm{mJ}$ . In the second calculation, $E_{\min }^2(P_{10}^{0.8})=177.0\mathrm{mJ}$ . The $E_{\min }$ change curves of seeds ${\mathit{p}}_1$ ∼ ${\mathit{p}}_5$ and ${\mathit{p}}_{11}$ ∼ ${\mathit{p}}_{15}$ in the parameter optimization process are shown in Figure 7. The blue “o” indicates the minimum value of 182.8 mJ in the first optimization calculation of ${\mathit{p}}_{14}$ , and the red “o” indicates the minimum value of 177.0 mJ in the second calculation.

Table 3.

Analysis of $E_{\min }$ in gradient optimization.

$k$	$n_k$	$n_k^{*}$	$E_{\min }(P_k)$	$E_{\min }^1(P_k)$	$D_1\%$	$E_{\min }^2(P_k)$	$D_2\%$
1	5	5	197.2	182.2	7.64	177.4	2.65
2	2	0	/	/	/	/	/
3	1	0	/	/	/	/	/
4	1	0	/	/	/	/	/
5	12	10	197.6	182.8	7.47	177.0	3.21
6	3	1	210.7	183.9	12.71	179.1	2.66
7	1	0	/	/	/	/	/
8	1	0	/	/	/	/	/
9	36	35	201.1	182.2	9.42	177.0	2.85
10	1	0	/	/	/	/	/
11	10	9	202.8	182.2	10.16	177.7	2.47
12	1	0	/	/	/	/	/
13	10	9	212.1	182.2	14.13	177.4	2.65
14	1	0	/	/	/	/	/
15	7	4	228.7	182.5	20.19	177.7	2.64
16	12	10	198.3	182.2	8.13	177.4	2.65
17	13	12	200.2	182.2	9.00	177.2	2.71
18	5	3	219.1	182.9	16.52	177.4	3.02
19	3	2	247.4	182.9	26.09	177.4	3.02
20	3	2	280.0	182.2	34.93	177.8	2.42
21	14	13	201.7	182.0	9.77	177.0	2.75
22	5	4	213.3	182.2	14.58	177.2	2.71
23	3	2	246.3	182.3	26.01	177.3	2.74
24	3	2	278.1	182.2	34.51	177.8	2.40

Figure 7.

Convergence curve of $E_{\min }$ in parameter optimization.

Effect of K

In the parameter description of biped robot body motion, K is an important factor affecting the optimization results. Taking the optimization calculation of data item $(10,c,0.8)$ as an example, the component values of K are $K_x=\{2,4,6,8\}$ , $K_y=\{1,3,5,7\}$ and $K_z=\{3,5,7,9\}$ . Figure 8 shows the solved $E_{10}^{0.8}$ , in which $K_x$ and $K_y$ are taken as coordinate axes, and different graphs correspond to different values of $K_z$ . For comparison, some typical energy consumption values are marked in the figure, such as at $K_x=2$ and $K_y=7$ where $K_z$ is 3, 5, 7, and 9, and the energy consumption values are 1.856, 1.829, 1.820 and 1.814 $\times 100\mathrm{mJ}$ , respectively. When $K_x=2$ , $K_y=1$ , and $K_z=5$ , the energy consumption is 2.589 $\times 100\mathrm{mJ}$ ; when $K_x=8$ , $K_y=7$ and $K_z=5$ , the energy consumption is 1.779 $\times 100\mathrm{mJ}$ . With increasing K, the energy consumption obtained by the optimization calculation decrease, and compared with $K_x$ and $K_z$ , the change caused by the increased in $K_y$ is greater.

Figure 8.

$E_{10}^{0.8}$ for different K values.

In related research, Luo and Chen ¹⁶ adopted the fixed hip height (FHH) algorithm,16 whereas Shin and Kim ¹⁹ allowed the hip height to fluctuate according to a cosine function (CHH).19 The FHH and CHH algorithms are compared with the optimization algorithm with $K=24$ in this paper. Taking the $E_{10}^{0.8}$ data of $(10,c,0.8)$ as an example, the $E_{\min }^{*}$ values for the FHH and CHH algorithms are 222.8 and 198.7 mJ, respectively, and the energy consumption is reduced by 20.56% and 10.92% when using our algorithm, respectively. The $x_{\mathrm{zmp}}$ trajectory of the sampling point and the energy consumption in the sampling period during a single cycle step period are shown in Figure 9. Compared with the FHH and CHH algorithms, the energy consumption in biped robot walking is more stable with our algorithm, and the total energy consumption is lower.

Figure 9.

Comparison of $E_{10}^{0.8}(n)$ and $x_{\mathrm{zmp}}(n)$ for the three algorithms.

Gait planning

In the gait planning calculation, the associated database subset is $(s,\alpha ,\eta ,E_s^{\eta })$ ; because $s_0$ and $s_e$ only act in the starting and stopping stages, except for the movement over very short distances, the stable cyclic step stage contributes the main part of the walking energy consumption. Figure 10 shows the results of the optimization calculation of $s\in S_c$ and $\eta \in \mathrm{H}$ for parameter length $K=18$ (where $K_x=6$ , $K_y=5$ , and $K_z=7$ ). To facilitate the comparison of energy consumption values, the representative parameter of the relative energy consumption $E/s$ is adopted in the figure. When $\eta =0.8$ and s is {7, 8, 9, 10, 11}, the calculation result E is {128.9, 143.7, 159.7, 178.4, 197.9}, then, $E/s$ is {18.41, 17.96, 17.74, 17.84, 17.99}. The $E/s$ of $E_9^{0.8}$ is the lowest. In Figure 10, the position of the minimum value of $E/s$ along the same η line is marked with “*.” The calculation results show that with decreasing η, the s corresponding to the minimum value of $E/s$ gradually increases.

Figure 10.

$E/s$ analysis when $s\in S_c$ and $\eta \in \mathrm{H}$ .

According to equation (12), the planning step sequence S with walking distance d and $s_m\in S_c$ can be calculated. Then, using equation (13), the predicted energy consumption of S can be calculated, and $s_m^{*}$ corresponding to the minimum walking energy consumption can be obtained. Figure 11 plots the calculated $s_m^{*}$ change curve when $d=3,4,\cdots ,{10}^4$ cm and $\eta =0.2,0.4,\cdots ,1$ . To expand the representation range, a logarithmic scale is adopted for d in the figure, and the difference between adjacent coordinates $d_{i-1}$ and $d_i$ is greater than 1. The value of $(d_i,s_m^{*})$ is the average value of all $s_m^{*}$ between $d_{i-1}$ and $d_i$ . The figure shows that when d is small, s _m is mainly selected to reduce the number of steps. With increasing d, the relative energy consumption in the starting and stopping stages decreases, and $s_m^{*}$ depends on the minimum $E/s$ .

Figure 11.

$s_m^{*}$ relation curve when $d=3,4,\cdots ,{10}^4\mathrm{cm}$ .

According to the GSYN, if $d=80\mathrm{cm}$ , then $\eta =0.8$ , and when $s_m=10\mathrm{cm}$ , the total energy consumption $J_s^{0.8}(10)=1.540J$ is the smallest. Here, $c=9$ , $s_0=3\mathrm{cm}$ , $s_1=7\mathrm{cm}$ , and $S=\{3,7,10,10,10,10,10,10,7,3\}$ in the planned walking, and the starting steps $(3,0,0.8,{\mathit{g}}_3^{0.8})$ , $(3,c,0.8,{\mathit{g}}_3^{0.8})$ , $(7,c,0.8,{\mathit{g}}_7^{0.8})$ , $(10,c,0.8,{\mathit{g}}_{10}^{0.8})$ and $(3,e,0.8,{\mathit{g}}_{10}^{0.8})$ are queried from the database to control the movement of the robot actuator. The corresponding walking dynamic simulation of the biped robot is shown in Figure 12, in which the black spot at the trunk of the robot is ${\mathit{r}}_b(n)$ .

Figure 12.

Dynamic simulation of robot walking.

Robot walking experiment

The robot simulation model is designed in the virtual reality simulation software Webots according to the physical robot structure, and the dynamic performance of the above gait optimization results is evaluated. In the biped robot walking experiment, $d=80\mathrm{cm}$ , the initial $\eta =0.5$ , the step set $S=\{3,7,10,10,10,10,10,10,7,3\}$ is planned, the database is queried, and the starting step $(3,0,0.5,{\mathit{g}}_3^{0.5})$ is obtained. If the walking stability requires $\eta \le 0.8$ , then the AZR controller will adjust the gait during walking to meet the stability requirements. Some snapshots of the robot walking are shown in Figure 13.

Figure 13.

Walking experiment of the biped robot in a virtual scene.

In the walking control process, $T_S=1.6s$ , $k_P=0.6$ , $T_D=0.4s$ , $T_L=1.6s$ , and the initial conditions $e_0=e_{-1}=0$ and ${\eta }_0=0.5$ . The value detected by the foot pressure sensor during robot walking are shown in Figure 14. In the figure, the blue marks represent the right foot force, the red marks represent the left foot force, and the values obtained from the four pressure measurement points on the sole are represented by the symbols *, o, ×, and ◊.

Figure 14.

Plantar force during robot walking.

At the beginning of robot walking, to ensure stability, η is set low, and the planned ZMP trajectory is close to the center of the sole. As shown in Figure 14, the initial left foot support pressure is relatively balanced. Then, according to the detected ZMP correction, in step s₁, η is corrected to 0.7, corresponding to sampling points 16–32 in Figure 15; the pressure begins to be dispersed, and the feedback ZMP is close to the edge of the sole of the foot. After step s₂, η is stable at 0.8. Although among the feedback ZMP values, there is a smaller η corresponding to the AZR edge line, the small short-term deviation is not enough to change η. The planned ZMP trajectory and feedback ZMP trajectory of robot walking are shown in Figure 15.

Figure 15.

Planned ZMP trajectory and feedback ZMP trajectory during robot walking.

Based on the virtual reality simulation, a performance test of physical robot walking is carried out. The calculation data for the robot walking process are shown in Table 4.

Table 4.

GSYN data of robot walking $d=80\mathrm{cm}$ .

$i$	$s_i/\mathrm{cm}$	$e_i$	$\mathrm{\Delta }{\eta }_i$	${\eta }_i$	$E_s^{\eta }/\mathrm{mJ}$	$E_s^{0.5}/\mathrm{mJ}$
0	3(0 + 3)	0	0	0.5	90.46	90.46
1	7(3 + 7)	0.27	0.18	0.7	110.99	127.72
2	10(7 + 10)	0.09	0.13	0.8	153.68	176.04
3	10(10 + 10)	–0.02	0.046	0.8	178.41	200.32
4	10(10 + 10)	–0.04	–0.003	0.8	178.41	200.32
5	10(10 + 10)	0.03	0.022	0.8	178.41	200.32
6	10(10 + 10)	0.01	0.015	0.8	178.41	200.32
7	10(10 + 10)	0.03	0.027	0.8	178.41	200.32
8	7(10 + 7)	–0.02	–0.002	0.8	153.68	176.04
9	3(7 + 3)	0.03	0.021	0.8	107.34	127.72
e	3(3 + 0)	0.02	0.022	0.8	66.79	94.85

In the table, $i=0$ is the starting step, $i=e$ is the stopping step, and the indices in between correspond to the cyclic step. The data in brackets in column $s_i$ correspond to the moving distance of the swinging leg behind and in front of the body, the data in column $E_s^{\eta }$ correspond to the energy consumption according to the calculated ${\eta }_i$ , and the data in column $E_s^{0.5}$ refer to the energy consumption calculated according to the initial setting ${\eta }_i=0.5$ . $J_s^{\eta }(10)=1.575J$ and $J_s^{0.5}(10)=1.794J$ are calculated; that is, the energy consumption is reduced by 12.21% through real-time adjustment of ${\eta }_i$ by the AZR controller. The video frames of the biped robot walking experiment are shown in Figure 16.

Figure 16.

Video frames of physical robot walking.

Conclusion

A biped robot is a robot system designed with reference to the structure of the human body. It is expected to replace or accompany human beings in some environments and carry out rescue in case of danger. Biped walking is one of the basic problems in biped robot research. The gait planning algorithm of a biped robot proposed in this article comprehensively considers the stability and energy efficiency of robot walking. The algorithm has universal applicability for different forms of biped robots. The following conclusions can be drawn from the theoretical analysis and experimental results:

• (a) To obtain the gait with the minimum energy consumption that meets the stability requirements in the movement process of a biped robot, the concept of a limited space of robot body motion is put forward for the GPO algorithm, which establishes a mapping relationship between the descriptive parameters and robot body motion trajectory. Through sampling and inverse kinematics calculation of the range of the descriptive parameter values and seed selection based on the sample energy consumption, combined with the seed domain space, the directional gradient optimization algorithm is applied. The optimization problem of obtaining the minimum energy consumption in walking within the large space of robot body motion is solved.

• (b) The length K of the parameter description has an important influence on the obtained gait parameters and energy consumption. Compared with the X and Z directions, the optimization calculation results are more sensitive to changes in the parameter length in the Y direction. For sampling of the descriptive parameters, there cannot be multiple peaks within the sampling gap. When the value space is large, a multistep gradient approximation method can be selected to accelerate the convergence process. The GPO algorithm requires many calculations in the optimization process. By using the geometric relationship of leg joints in the walking process of biped robots, an acceleration algorithm suitable for GPU parallel operation is established, which can quickly determine the optimal gait values for robot walking. With the progress of computer science, the sampling approximation method for descriptive parameters will show superior performance compared with the random optimization algorithm.

• (c) For the GSYN algorithm, a gait database is established with the walking step size, gait type and AZR variables as input items and the minimum energy consumption and gait parameters as output. Given the walking distance of a biped robot, the GSYN algorithm proposes a robot walking trajectory with the minimum energy consumption as the optimization index according to the gait database. For real-time walking control of a robot, a method to modify the AZR variable in the database query according to the feedback ZMP data to adjust the gait parameters of the robot is proposed. The GSYN algorithm dynamically adjusts the gait to meet different AZRs, effectively overcomes the interference of robot modeling and walking environment errors, realizes a compromise between low energy consumption and high robustness, and better solves the energy-saving walking problem of biped robots with highly nonlinear characteristics.