A grid gradient approximation method of energy-efficient gait planning for biped robots

Abstract

In this article, an energy-efficient gait planning algorithm that utilizes both 3D body motion and an allowable zero moment point region (AZR) is presented for biped robots based on a five-mass inverted pendulum model. The product of the load torque and angular velocity of all joint motors is used as an energy index function (EIF) to evaluate the energy consumption during walking. The algorithm takes the coefficients of the finite-order Fourier series to represent the motion space of the robot body centroid, and the motion space is gridded by discretizing these coefficients. Based on the geometric structure of the leg joints, an inverse kinematics method for calculating grid intersection points is designed. Of the points that satisfy the AZR constraints, the point with the lowest EIF value in each network line is selected as the seed. In the neighborhood of the seed, the point with the minimum EIF value in the motion space is successively approximated by the gradient descent method, and the corresponding joint angle sequence is stored in the database. Given a distance to be traveled, our algorithm plans a complete walking trajectory, including two starting steps, multiple cyclic steps, and two stopping steps, while minimizing the energy consumption. According to the preset AZR, the joint angle sequences of the robot are read from the database, and these sequences are adjusted for each step according to the zero-moment-point feedback during walking. To determine the effectiveness of the proposed algorithm, both dynamic simulation and walking experiment in the real environment were carried out. The experimental results show that compared with algorithms based on the fixed body height or vertical body motion, our gait algorithm has a significant energy-saving effect.

Keywords

Biped robot five-mass model gait planning zero moment point spatial gridding gradient approximation

Introduction

Bipedal robots have human-like structures and appearances, which can adapt to the human environment, and are ideal robots for replacing human work. In the past 50 years, research institutions and scholars at home and abroad have carried out much research on biped robots and have achieved remarkable progress.^1,2 For example, humanoid robot HRP-4 can drive vehicles in the middle of the road,³ and ATLAS can traverse obstacles and climb stairs.⁴

Bipedal walking is the most important and challenging problem in research on biped robot motion. According to different research ideas, the proposed biped walking methods can be divided into three approaches. One is to utilize high-speed optical motion capture systems to obtain data on human motion according to external characteristics of human walking and apply these features to the generation of robot motion modes.^5,6 The other approach is to use a central pattern generator to simulate the neural network control of human walking to generate rhythmic signals, thereby solving the problem of robot gait generation.^7
-9 The third approach is to simplify the biped robot into a linear inverted pendulum model consisting of a point mass and a massless telescopic leg^10

-13 or a cart-table model consisting of a table without mass and a cart driving on the table with all the mass concentrated in it.¹⁴

At present, the structural complexity of existing biped robots is much less than that of humans, which indicate that gait methods based on human characteristics and central pattern generators have limitations. Researchers prefer to adopt the method that involves simplifying the biped robot model. The accuracy of the simplified model determines the accuracy of the control effect. Based on the single-mass inverted pendulum model and considering the weight of the legs, Shimmyo et al. proposed biped walking pattern generation using preview control based on a three-mass model.¹⁵ Luo and Chen presented a three-mass angular momentum model¹⁶ and a five-mass momentum model¹⁷ using model predictive control to obtain better motion control accuracy. However, as model complexity increases, model nonlinearity increases and gait control becomes more difficult. Considering that the leg mass of the motor-driven biped robot is mainly distributed according to the position of the motor, the five-mass model, which includes five mass points representing the body, legs, and feet, is regarded as the best tradeoff.

To simplify the complexity of the algorithm used in gait control for biped robots, constraint conditions for the robot body are introduced, which cause unnatural walking motions, consume large amounts of energy, and limit the operation time of battery-powered robots.¹⁸ Hong et al. allowed vertical body motion and alleviated the problem of bending knee joints.¹⁹ Shin and Kim further removed constraints on the motion trajectory, which improved walking speed and reduced actuator energy consumption.²⁰ If the robot body is allowed to move with more degrees of freedom, a better control effect can be expected. In this article, the bipedal body center of mass is allowed to move in three dimensions.

Fewer constraint conditions result in an increase in parameter space in gait control. It would be optimal if a solution could be found by artificial intelligence. Dau et al. applied a genetic algorithm to optimize the seven key parameters defining the hip and foot trajectories.²¹ Elhosseinia et al. designed a whale optimization algorithm with a random parameter “a” and weight parameter “C” to find the optimal setting for the hip parameters.²² Wu and Li used fuzzy logic to control the dynamic gait pattern generator of a humanoid robot, resulting in better responses to external force disturbances.²³ Wright and Jordanov summarized the commonly used intelligent approaches to robot locomotion control.²⁴ Generally, the range and accuracy of parameters need to be included in the optimization algorithm. This article presents a hierarchical optimization algorithm, that is, the grid computing method, to find the best region in a large range of parameter value spaces. Then, the gradient approximation method is used to find the optimal gait parameters for the robot with high accuracy.

The rest of this article is organized as follows. In the second section, we describe the simplified biped robot model and formulate the problem of gait planning. The third section focuses on the gait planning algorithm from two aspects: optimal minimization of energy consumption and real-time motion planning. According to the algorithm in this article and various control methods proposed in the related literature, the fourth section presents the simulation and experimental results for the performance analysis of the proposed algorithm. In the fifth section, the algorithm is summarized, and possible follow-up work is presented.

Problem formulation

System structure

The humanoid robot used in this article has two arms and two legs, thereby imitating the walking motion of the human body. Each leg has five degrees of freedom (DoFs), including two DoFs at the hip, one DoF at the knee, and two DoFs at the ankle, and the joint vector can be expressed as $q = {[\begin{matrix} q_{1} & q_{2} & \dots & q_{10} \end{matrix}]}^{T}$ . The forward, lateral, and vertical directions of the robot are defined as the x-axis, y-axis, and z-axis, respectively. The origin in the global coordinate system is located at the midpoint of the two legs when the robot is upright. The structure of the experimental robot is shown in Figure 1. Four force-sensing resistor (FSR) sensors are installed under the sole of each foot to measure the distribution of force on the feet. The parameters for the body and right leg are summarized in Table 1, and the parameters for the left leg are the same as those for the right leg.

Figure 1.

(a, b) 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_{r u}$	Right leg mass	0.2 kg
$m_{r d}$	Right foot mass	0.2 kg
l_h	Hip width	0.08 m
l ₁	Link length from q₁ to q₂	0.04 m
l ₂	Link length from q₂ to q₃	0.085 m
l ₃	Link length from q₃ to q₄	0.075 m
l ₄	Link length from q₄ to foot	0.04 m
$l_{f l}$	Foot length	0.12 m
$l_{f w}$	Foot width	0.06 m

In general, some constraints are needed to simplify the scope of an analysis of robot walking. The proposed walking gait assumes the following conditions:

The upper body remains upright at all times. The pitch rotation of the human trunk is small, within 3°,²⁵ and energy consumption increases as the trunk is leaned forward.²⁶ Therefore, most relevant research studies^{4,5,10

-23} have shown that this assumption is acceptable.

Both feet are always parallel to the ground. Most common humanoid robots have no toes and cannot utilize the toes to improve driving the lifting and planting of the feet.^27,28

One step of duration T contains a double support phase (DSP) $T_{DSP}$ and a single support phase (SSP) $T_{SSP}$ and defines the duty ratio of the DSP as $σ \frac{T_{DSP}}{T}$ . In the process of human walking, $σ$ is approximately 15–25%²⁵; $σ = 25 %$ was selected for the algorithm proposed in this study.

Zero moment point equations

Among various dynamic stability standards for biped robots, the most widely used technique is the zero moment point (ZMP).²⁹ The ZMP is the point on the ground where the horizontal component of the total moment generated by gravity and inertia forces is zero. The five mass points of the biped robot are the body b, the left leg $l u$ , the left foot $l d$ , the right leg $r u$ , and the right foot $r d$ , and they constitute set $P_{c} = {b, l u, l d, r u, r d}$ . The positions of the six joints form a position set $R_{l} = {r_{l h}, r_{l k}, r_{l f}, r_{r h}, r_{r k}, r_{r f}}$ . These positions are shown in Figure 2.

Figure 2.

Five-mass biped robot model: (a) Sagittal plane and (b) frontal plane.

At ZMP position $r_{ZMP} = {[x_{ZMP}, y_{ZMP}, 0]}^{T}$ , the resultant moment $M_{ZMP}$ is generated by the resultant ground reaction force $F$ . The resultant moment $M_{s}$ is exerted on the robot body due to gravity and inertia of all centers of mass (CoMs). The moment balance equation for $r_{ZMP}$ is expressed as

n \times M_{ZMP} = n \times (M_{s} - r_{ZMP} \times F) = 0

where

F = \sum_{i \in P_{c}} m_{i} (g - {\ddot{r}}_{i})

M_{s} = \sum_{i \in P_{c}} [r_{i} \times m_{i} (g - {\ddot{r}}_{i}) - I_{i} \times {\ddot{θ}}_{i}]

Here, n = [0,0,1]^T is the ground-surface unit normal vector, g = [0,0,-g]^T is the gravitational acceleration vector, g = 9.80 m/s², and I _i and ${\ddot{θ}}_{i}$ are the rotational inertia and angular acceleration at position i, respectively. Kajita et al. proved that the influence of $\ddot{θ}$ is small and can be ignored.¹⁰ By substituting equations (2) and (3) into equation (1), the simplified ZMP equations are

{\begin{matrix} x_{zmp} = \frac{\sum_{i \in P_{c}} m_{i} ({\ddot{z}}_{i} x_{i} + g x_{i} - z_{i} {\ddot{x}}_{i})}{\sum_{i \in P_{c}} m_{i} ({\ddot{z}}_{i} + g)} \\ y_{zmp} = \frac{\sum_{i \in P_{c}} m_{i} ({\ddot{z}}_{i} y_{i} + g y_{i} - z_{i} {\ddot{y}}_{i})}{\sum_{i \in P_{c}} m_{i} ({\ddot{z}}_{i} + g)} \end{matrix}

Allowable zero moment point region

Modeling errors inevitably exist when modeling biped robots. Hong et al. set the ZMP trajectory near the centerline of the support foot to achieve the greatest stability during walking,¹¹ but this is not an energy-efficient method. In the region of the support foot, some edge regions are drawn out to compensate for the modeling errors,²⁰ and the ZMP trajectory is located in the allowable ZMP region (AZR) of the middle subregion of the support foot. This is a better method for balancing modeling errors and walking efficiency. For a biped robot with foot length $l_{f l}$ and foot width $l_{f w}$ , the step length is s, and the y-axis distance between the two feet is w. The AZR when the left foot (LF) is supporting and the right foot (RF) is swinging is shown in Figure 3. The AZR during the first DSP is a hexagon with ( r ₁, r ₂,… r ₆), that during the SSP is a square with ( r ₃, r ₄, r ₅, r ₆), and that during the second DSP is a hexagon with ( r ₅, r _6… r ₁₀). In our gait planning algorithm, $\begin{array}{l} η_{A Z R}^{i} = l_{a l} / l_{f l} \end{array}$ and $η_{A Z R}^{w} = l_{a w} / l_{f w}$ are used to represent the unit values of the AZR subregions on the x and y axes, respectively, $η_{A Z R}^{l}$ takes a constant value, so the default is $η = η_{A Z R}^{w}$ .

Figure 3.

AZR of biped robots. AZR: allowable zero moment point region.

Energy consumption index function

The energy consumption E of a biped robot can be divided into two parts: the energy consumption E_m for robot motion and the energy consumption E_a for nonmotion. E_m is the main component of E, which is the integral of the instantaneous power consumption vector $p_{m} (t)$ for all joint motors over time t. E_a is used for sensors, controllers, and so on. Their power P_a is relatively stable during operation and can be expressed as a linear function of time. E can be described by

E = E_{m} + E_{a} = \int_{0}^{t} p_{m} (ξ) d ξ + P_{a} t

In general, E_m can be significantly reduced through gait trajectory optimization,³⁰ which is suitable for evaluating the performance of the gait algorithm. For biped robots with many motors, it is difficult to accurately measure p_m(t). p_m(t) is equal to the sum of the output power consumption p₂(t) and subsidiary electrical loss p_loss(t); p₂(t), its main component, can be calculated by multiplying the motor output torque τ₂(t) and angular velocity $\dot{q} (t)$ .³¹ When the robot is moving, τ₂(t) is used to overcome the load moment τ_g(t) formed by all CoMs on the motor rotor, where τ_g(t) = [τ₁(t) τ₂(t)…τ₁₀(t)]^T and τ_i(t) can be written as follows

τ_{i} (t) = {l_{i}}^{T} (t) C_{i} (t) m, i = 1, 2, \dots, 10

where $m = {[\begin{matrix} m_{b} & m_{l u} & m_{l d} & m_{r u} & m_{r d} \end{matrix}]}^{T}$ is the CoM vector, $l_{i} (t)$ is the effective distance vector between $m$ and the i’th rotor, and $C_{i} (t)$ is the $5 \times 5$ diagonal matrix describing the proportion of $m$ acting on the rotor.

It is assumed that the sampling period t_s of the control system is small enough between the two sampling points with $τ_{i} (n) = τ_{i} (n \cdot t_{s}) \approx τ ((n - 1) \cdot t_{s})$ , and ${\dot{q}}_{i} (n \cdot t_{s}) \approx \frac{q_{i} (n \cdot t_{s}) - q_{i} ((n - 1) \cdot t_{s})}{t_{s}}$ , where $i = 1, 2, \dots, 10$ . In general, the joint motors of the robot carry out cyclical motion with the period N, and even if gravity and the direction of motion are the same, most robots are not able to regenerate energy. Therefore, we define the energy consumption index function E for the gait algorithm as follows

E = \frac{1}{t_{s}} \sum_{n = 1}^{N} \sum_{i = 1}^{10} τ_{i} (n) [q_{i} (n) - q_{i} (n - 1)]

Problem definition

Based on the previous analysis, our energy-efficient gait planning problem can be formulated as follows:

Given the step length and AZR of a biped robot, find the optimal gait parameters, that is, the parameters that minimize E, while allowing 3D body motion and satisfying the constraints in “System structure” section.

For a goal distance d of bipedal walking, control the value of $η$ for the AZR in real time, thereby minimizing the total energy consumption of the complete gait trajectory, including two starting steps, multiple cyclic steps, and two stopping steps.

Gait planning algorithm

Overview

In the gait planning optimization (GPO) algorithm, given the step length set S and the AZR set H for bipedal walking, 18 parameters are represented by the body trajectory $s \in S$ , and the range U_s of the parameters is gridded. After the inverse kinematics calculation, the parameter subset R_s satisfying the stability and physical constraints is obtained, and the seed set $P_{s}^{η}$ meeting the $η \in H$ requirement is selected from R_s. Then, in the neighborhood of $p \in P_{s}^{η}$ , the joint angle sequence $g_{s}^{η} = {q_{s}^{η} (n) | n = 1, 2, \dots, N}$ with minimum E is obtained by iterative calculation according to the gradient descent method and stored in an offline database. The GPO algorithm takes a long time and is suitable for offline operation. It completes the calculation of S and H and copies the offline database to the online database, which is called in real time in the gait synthesis (GSYN) algorithm.

Given the required walking distance d and the value of η for the AZR, the GSYN algorithm plans the step sequence S* to achieve the minimum E. Each step length $s_{i} \in S^{*}$ and initial value $\hat{η}$ of the AZR are removed in turn, and the online database is queried to obtain the motor angle sequence g _i for gait control. By inputting g _i into the joint controllers, the robot can walk. During walking, according to the foot pressure set F_i, the real-time ZMP trajectory is calculated. In the AZR control, the deviation between r _ZMP and r_AZR is calculated from r_AZR, which is required for stability, and the PI correction method is used to correct η_i, optimizing the tradeoff between low energy consumption and robustness of robot motion.

Gait planning optimization algorithm

The GPO algorithm is used to generate a low-energy joint angle sequence $g_{s}^{n}$ for gait control; this is the core of our proposed algorithm. According to the accuracy requirements of bipedal walking, the range in step length, from minimum to maximum, is discretized to constitute the set S. H selects several levels from the interval [0,1] to adapt to different robustness requirements. The GPO algorithm selects the elements s and η from S and H, respectively, to calculate the distribution gradient and obtain $g_{s}^{n}$ for the minimum E.

Gridding method

Our proposed gridding method is a process of discretizing multidimensional parameters by appropriate intervals, constructing parameter space grids, and then calculating the value of the function for each gridded line intersection point. If the walking step length is s, the robot body b performs a 3D cyclic motion satisfying the constraint (a) in “System structure” section, and the position $r_{b} (n) = {[\begin{matrix} x_{b} (n) & y_{b} (n) & z_{b} (n) \end{matrix}]}^{T}$ can be expressed by a Fourier series of finite terms as follows

{\begin{cases} x_{b} (n) = \frac{s}{N} n + \sum_{k = 1}^{3} a_{k} sin (k ω_{0} n) \\ + \sum_{k = 4}^{6} a_{k} cos ((k - 3) ω_{0} n) \\ y_{b} (n) = \sum_{k = 7}^{11} a_{k} sin ((k - 6) ω_{0} n) \\ z_{b} (n) = a_{12} + \sum_{k = 13}^{15} a_{k} sin ((k - 12) ω_{0} n) \\ + \sum_{k = 16}^{18} a_{k} cos ((k - 15) ω_{0} n) \end{cases}

where N is the gait period, $ω_{0} = π / N$ , $n = 1, 2, \dots, 2 N$ , and $r_{b} (n)$ are represented by the parameters ${a_{k} | k = 1, 2, \dots, 18}$ . If $a_{k} \in [d_{k}, u_{k}]$ is discretized according to the interval $μ_{k}$ and its range is the set $A_{k} = {d_{k}, d_{k} + μ_{k}, d_{k} + 2 μ_{k}, \dots, u_{k}}$ , then the intersection set U_s formed by the gridding method is the Cartesian product of A_k, that is, $U_{s} = A_{1} \times A_{2} \times \dots \times A_{18}$ .

Inverse kinematics calculation

To utilize the parallel acceleration of the GPU, the inverse kinematics algorithm is designed for the geometric structure of the robot leg. Without loss of generality, suppose that the robot starts walking during the DSP with the left leg in front and the right leg in the rear; then, the right leg swings forward. The starting position of the RF $r f$ is $r_{r f} (0) = {[\begin{matrix} 0 & 0 & 0 \end{matrix}]}^{T}$ , and the motion satisfies constraint (c) in “System structure” section; then, the motion trajectory of $r_{r f} (n) = {[\begin{matrix} x_{r f} (n) & y_{r f} (n) & z_{r f} (n) \end{matrix}]}^{T}$ is

{\begin{cases} x_{r f} (n) = {\begin{cases} 0, 1 \leq n \leq N_{1} \\ s (1 - cos \frac{(n - N_{1}) π}{K + 1}), N_{1} < n \leq N_{2} \\ 2 s, N_{2} < n \leq 2 N \end{cases} \\ y_{r f} (n) = - \frac{w_{s}}{2}, 1 \leq n \leq 2 N \\ z_{r f} (n) = {\begin{cases} 0, 1 \leq n \leq N_{1} \\ h_{s} sin \frac{(n - N_{1}) π}{K + 1}, N_{1} < n \leq N_{2} \\ 0, N_{2} < n \leq 2 N \end{cases} \end{cases}

where w_s is the step width, h_s is the step height, N is the gait cycle, $N_{1} = N / 8$ , and $N_{2} = 7 N / 8$ . Cyclic bipedal walking has the property of symmetry, and the LF trajectory is $r_{l f} (n) = {[\begin{matrix} x_{r f} (n + N) + s & - y_{r f} (n + N) & z_{r f} (n + N) \end{matrix}]}^{T}$ .

Given $r_{b} (n)$ , according to constraint (a) in “System structure” section, the right hip joint trajectory is $r_{r h} (n) = r_{b} (n) - {[\begin{matrix} 0 & 0.5 l_{h} & l_{b} \end{matrix}]}^{T}$ . Figure 4 shows the geometric relations between the joint angles of the robot leg.

Figure 4.

Geometric relations between $r_{r h}$ and $r_{r f}$ .

Combining $r_{h} = [r_{r h}, r_{l h}]$ and $r_{f} = [r_{r f}, r_{l f}]$ that need to be calculated, the process for obtaining the inverse kinematics solution presented in this article is shown in Algorithm 1.

Algorithm 1

Inverse kinematics algorithm

Generally, U_s has a large number of elements. Considering the locomotion characteristics of $r_{b} (n)$ , unreasonable elements in U_s can be excluded, thus reducing the number of inverse kinematics calculations needed. Common characteristics of $r_{b} (n)$ include the following

$x_{b} (n)$ increases predictably during bipedal walking.

$y_{b} (n)$ is convex during RF swing N and concave during LF swing N.

$z_{b} (n)$ increases predictably and then decreases predictably during N.

During walking, the given $r_{h} (n)$ and $r_{f} (n)$ distance are less than the leg length, that is, $max {‖ r_{h} (n) - r_{f} (n) ‖}_{2} \leq \sum_{j = 1}^{4} l_{j}, n = 1, 2, \dots N$ .

Seed extraction

For $\forall a \in U_{s}$ , the map $a \to {q (0), q (1), \dots, q (N)}$ is formed by inverse kinematics calculation, and then, the corresponding ZMP trajectory can be obtained. $R_{s} \subseteq U_{s}$ , where $a_{s} \in R_{s}$ , is defined, which satisfies the following:

Physical constraint: for $a_{s} \to q_{s} (n)$ , $max | q_{s} (n) - q_{s} (n - 1) | \leq t_{s} θ_{max}$ , where $n = 1, 2, \dots N$ , and $θ_{max}$ is the maximum angular velocity of the joint motors.

ZMP condition: the ZMP trajectory corresponding to $a_{s}$ satisfies the stability requirement.

The wider AZR consumes less energy but yields less walking stability.²⁰ In gait control of the biped robot, multiple $η$ values are taken to constitute the set $Η$ . The dynamic adjustment of $η$ during walking can achieve a better compromise between energy consumption and stability. The element $a_{s}^{η}$ in R_s that meets the requirement of $η$ is used to calculate $E (a_{s}^{η})$ according to equation (7). If the range of $E (a_{s}^{η})$ contains a finite number of convex sets, the seed set $P_{s}^{η} = {p_{1}, p_{2}, \dots, p_{Λ}}$ , where $p_{λ} \in P_{s}^{η}$ , satisfies the following:

There is a linear relationship between the number of seeds $Λ$ and the number of parameters. That is, $Λ \leq n_{1} + n_{2} + \dots + n_{18}$ , where n_k is the number of a_k values in $\forall a_{s}^{η}$ and $k = 1, 2, \dots, 18$ .

The energy consumption of the same column grid is the lowest. Namely, $p_{λ} = a_{k}^{j}$ if and only if $E (a_{k}^{j}) = min [E (A_{k}^{j})]$ , where $A_{k}^{j}$ is the set containing the j’th value of a_k in $\forall a_{s}^{η}$ , where $j = 1, 2, \dots, n_{k}$ .

Gradient optimization

Obviously, the seed $p_{λ} \in P_{s}^{η}$ is a gridded intersection in U_s. If $E (p_{λ})$ of $p_{λ} = [α_{1} α_{2} \dots α_{18}]$ achieves convexity in adjacent grids, then the search space $L_{λ} = B_{1} \times B_{2} \times \dots B_{18}$ is formed with $p_{λ}$ as the center and $^{δ_{λ} = [δ_{1} δ_{2} \dots δ_{18}]}$ as the search interval, where $B_{k} = {α_{k} - δ_{k}, α_{k} + δ_{k}}$ and $k = 1, 2, \dots, 18$ .

For $ξ \in L_{λ}$ , the inverse kinematics and E are calculated. If $ξ^{*}$ satisfies $E (ξ^{*}) = min {E (ξ) | ξ \in L_{λ}}$ and its ZMP meets the requirements of $η$ , then taking the gradient $\nabla E (ξ^{*}) = ξ^{*} - p_{λ} = [\begin{matrix} β_{1} & β_{2} & \dots & β_{18} \end{matrix}]$ , a new search space $L_{λ}^{*} = B_{1}^{*} \times B_{2}^{*} \times \dots \times B_{18}^{*}$ is obtained, where $B_{k}^{*} = B_{k} + β_{k}$ . The specific calculation process is shown in Algorithm 2.

Algorithm 2

Gradient optimization algorithm

Given $p_{λ}$ , Algorithm 2 is progressive in the finite grid neighborhood, according to $δ_{λ}$ , so it must converge. In the algorithm, $ε$ is the approximation accuracy of E. After all, $p_{λ}$ is calculated, the seed with the minimum E can be used to obtain the optimal gait $g_{s}^{η}$ corresponding to s and $η$ .

Gait synthesis algorithm

For a given walking distance d, there is a sequence $S^{*} = {s_{1}, s_{2}, \dots, s_{c}}$ with minimum E. To make the walking process smooth, both s₁ and s₂ are planned as starting steps, both $s_{c - 1}$ and s_c are stopping steps, and s_m is a cyclic step, where $m = 3, 4, \dots, c - 2$ . We let $s_{1} = s_{c}$ , $s_{2} = s_{c - 1}$ , $\forall s_{m}$ cover equal distances and $d_{b} = 2 s_{1} + 2 s_{2}$ be controlled in $1 \sim 2 s_{m}$ , that is, $s_{m} < d_{b} \leq 2 s_{m}$ . When s_m is known, the calculation of s₁ and s₂ is shown in equation (10)

{\begin{cases} c = ⌈ d / s_{m} ⌉ + 2 \\ d_{b} = d - s_{m} (c - 4) \\ s_{1} = (d_{b} - s_{m}) / 3 \\ s_{2} = (d_{b} +2 s_{m}) / 6 \end{cases}

If the robot’s walking starts and ends with two feet standing side by side, the AZR is required to be $η$ , $E_{s}^{η} (s_{m})$ is the energy consumption of cycle step s_m, $E_{b}^{η} (s_{1})$ is the energy consumption of starting step s₁, and $E_{e}^{η} (s_{c})$ is the energy consumption of the stopping step s_c. Then, the total energy consumption $J_{s}^{η} (s_{m})$ is

J_{s}^{η} (s_{m}) = (c - 4) E_{s}^{η} (s_{m}) + E_{a}^{η} (s_{m})

where $E_{a}^{η} (s_{m}) = E_{s}^{η} (s_{1}) + 2 E_{s}^{η} (s_{2}) + E_{b}^{η} (s_{1}) + E_{e}^{η} (s_{c})$ .

Obviously, $\exists s_{m}^{*} \in S$ , such that $J_{s}^{η} (s_{m}^{*}) = min [J_{s}^{η} (S)]$ .

As shown in Figure 5, when the robot moves initially, the AZR should use a smaller $\hat{η}$ to ensure robustness during walking, but such walking will consume more energy.³² The GSYN algorithm uses closed-loop control to query gait trajectory $g_{i}$ from the database according to the s_i and $η_{i}$ of each step. While the robot is walking, the ZMP data $r_{ZMP} (n)$ are calculated from the feedback FSR data $F_{i} (n)$ ³³ and sent to the AZR controller to calculate $η_{i + 1}$ .

Figure 5.

Overall gait planning algorithm.

In the AZR controller, the time T_s of each step is taken as the control period, that is, T_s = N t_s, which includes the underactuated SSP and the overactuated DSP, which are applicable to different algorithms.³⁴ As shown in equation (9), N₁ < n < N₂ is the SSP, and both 1 ≤ n ≤ N₁ and N₂ < n ≤ N are the DSP. We let the AZR position during the i’th step satisfy robot walking stability as r _AZR (n). During the DSP, the AZR control algorithm requires r _ZMP(n) to be located in the AZR. During the SSP, the algorithm requires x_ZMP(N₁+1) ≥ x_AZR(N₁+1), x_ZMP(N₂) ≤ x_AZR(N₂), and η_i₊₁ is calculated from the deviation e_i in the Y direction, that is

e_{i} = \frac{2}{l_{f w}} min {c_{o} (y_{ZMP} (n) - y_{AZR} (n)) | N_{1} < n \leq N_{2}}

where $y_{ZMP} (n) = \sum_{j = 1}^{c_{n}} y_{i}^{j} (n) f_{i}^{j} (n) / \sum_{j = 1}^{c_{n}} f_{i}^{j} (n)$ , $y_{i}^{j} (n)$ , and $f_{i}^{j} (n) \in F_{i} (n)$ correspond to the sensor position and pressure in the Y direction, respectively. For the SSP, $c_{n} = 4$ , and for the DSP, $c_{n} = 6$ . When $y_{AZR} (n) \geq 0$ , $c_{o} = 1$ ; when $y_{AZR} (n) < 0$ , $c_{o} = - 1$ .

The incremental PI regulator with transfer function $D (s) = k_{P} (1 + \frac{1}{T_{I} s})$ is used to calculate $η_{i + 1}$ , that is

η_{i + 1} = η_{i} + Δ η_{i + 1}

where $Δ η_{i + 1} = k_{P} (e_{i} - e_{i - 1}) + k_{P} \frac{T_{S}}{T_{I}} e_{i}$ . Because $(e_{i} - e_{i - 1})$ easily causes high frequency interference, the first-order inertia link $G (s) = \frac{1}{T_{L} s + 1}$ is introduced to smooth its output, so that

Δ η_{i + 1} = K_{α} Δ η_{i} + K_{1} e_{i} + K_{2} e_{i - 1}

where $K_{α} = \frac{T_{L}}{T_{L} + T_{S}}$ , $K_{1} = k_{P} (\frac{T_{S}}{T_{L} + T_{S}} + \frac{T_{S}}{T_{I}})$ , and $K_{2} = - k_{P} \frac{T_{S}}{T_{L} + T_{S}}$ .

Experimental evaluation

According to the structure of the biped robot described in “System structure” section, a 3D simulation model is constructed in Webots to evaluate our proposed algorithm. In the GPO algorithm, we take $S = {s | s = 1, 2, \dots 13}$ and $H = {η | η = 0.1, 0.2, \dots 1}$ and calculate the gait trajectory $g_{s}^{n}$ from their combination. Taking s = 10 as an example, r_b(n) is described by a_k, as shown in equation (8), where N = 16, $μ_{k} = 0.2 cm$ , and k = 1, 2,…18; r_rf(n) is shown in equation (9), where $w_{s} = 8 cm$ and $h_{s} = 1 cm$ . Using the features extracted from r_b(n), the number of effective grid elements is 1 763 449 134, and the inverse kinematics calculation is carried out with Algorithm 1. According to the physical constraints and ZMP conditions in “Seed extraction” section, 17 024 963 elements of R_s are obtained, and the seeds corresponding to different η are selected.

When $η = 0.8$ , the number of seeds selected is $Λ = 102$ , as presented in Table 2, 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 seeds in P_k. To obtain a gait with lower energy consumption, gradient optimization is performed twice with Algorithm 2 with $δ_{k}^{1} = 0.05 cm$ and $δ_{k}^{2} = 0.0 1 cm$ . $E_{min}^{1} (P_{k})$ and $D_{1} %$ and $E_{min}^{2} (P_{k})$ and $D_{2} %$ are the minimum values of E and the percent reduction in E in the two calculations, respectively. The $E_{min}$ curve of $p_{1} \sim p_{20}$ in the gradient optimization calculation is shown in Figure 6. In the first gradient optimization calculation, $\forall p_{λ} \in P_{10}^{0.8}$ , $E_{min}^{1} (P_{10}^{0.8}) = 182.3 mJ$ , and the results are marked with “ .” In the second calculation, $E_{min}^{2} (P_{10}^{0.8}) = 178.4 mJ$ , and the results are marked with “ .”

Table 2.

Analysis of $E_{min}$ during 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.8	7.30	178.4	2.42
2	9	7	197.2	182.3	7.75	178.4	2.14
3	2	0	—	—	—	—	—
4	3	1	210.7	182.3	13.48	179.8	1.36
5	1	0	—	—	—	—	—
6	1	0	—	—	—	—	—
7	34	33	210.2	182.8	13.01	178.4	2.42
8	1	0	—	—	—	—	—
9	9	8	202.8	182.8	9.83	178.4	2.42
10	1	0	—	—	—	—	—
11	8	5	212.2	182.8	13.81	178.4	2.42
12	11	9	198.3	182.3	8.06	178.4	2.14
13	11	10	200.2	182.3	8.93	179.8	1.40
14	14	13	201.7	182.3	9.62	178.4	2.14
15	5	3	219.1	183.0	16.48	178.9	2.20
16	5	4	213.3	182.8	14.27	179.5	1.81
17	3	2	247.4	183.0	26.05	178.9	2.20
18	3	2	246.3	183.0	25.72	178.9	2.20

Figure 6.

Convergence curve for $E_{min}$ during gradient optimization.

To compare the performance of different algorithms, the fixed hip height (FHH) is taken as the FHH algorithm,¹⁶ the vertical variation in the hip height (VHH) according to the cosine waveform is taken as the VHH algorithm,²⁰ and the other conditions are the same as the algorithm in this article. After optimization, $E_{min}^{*} = 222.8 m J$ for FHH algorithm and $198.7 mJ$ for VHH algorithm, which means that the energy consumption is reduced by $19.93 %$ and $10.22 %$ , respectively. In a single gait cycle, the $x_{ZMP} (n)$ trajectories of the three algorithms and the corresponding E of each sampling point are identified in Figure 7, where $σ = 0.25$ , $N = 16$ , and $t_{s} = 0.1 s$ . It can be seen from the figure that compared with that of both FHH algorithm and VHH algorithm, the E distribution of our algorithm is more balanced, and the total energy consumption value is lower.

Figure 7.

Comparison of the $E (n)$ and $x_{ZMP} (n)$ of the three algorithms.

According to the GPO algorithm, $\forall s_{m} \in S$ and $\forall η = Η$ , the calculated $E_{min}^{*}$ is shown in Figure 8. Obviously, as the value of $η$ decreases, the required $E_{min}^{*}$ gradually increases to achieve the same s.

Figure 8.

$E_{min}^{*}$ corresponding to $s \in S$ and $η \in Η$ .

In the GSYN algorithm, if $d = 100 cm$ and $\forall s_{m} \in S$ , according to equations (10) and (11), $J_{s}^{η} (S)$ can be calculated, as shown in Figure 9, where $η = 0.2, 0.3, \dots, 1$ , and $s_{m} = 10 cm$ , $J_{s}^{η} (S)$ is the minimum in the same column data. In this case, $c = 8$ , $s_{1} = 3 cm$ and $s_{2} = 7 cm$ , that is, $S^{*} = {3, 7, 10, 10, 10, 10, 10, 10, 10, 10, 7, 3}$ .

Figure 9.

$J_{s}^{η} (S)$ analysis when $d = 100 cm$ .

If we let η = 0.8 and set the starting and stopping positions of walking to be two feet standing side by side, the swing length of the RF is L_rf = {3,17,20,20,17,3}, and the swing length of the LF is L_lf = {10,20,20,20,20,10}. The two 3-cm elements in L_rf correspond to the starting step and stopping step, respectively, where $E_{b}^{0.8} (3) = 65.2$ and $E_{e}^{0 .8} (3) = 66.8 m J$ . There are two 17-cm elements in L_rf and two 10-cm elements in L_lf, which correspond to the case when the movement distance of the swing foot is unequal before and after it becomes the support foot. From equation (11), $J_{s}^{0.8} (10) = 1.89 67 J$ , where $E_{a}^{0.8} (10) = 515.9 mJ$ . The gait array for each step g _i = [ q ₁ q ₂… q ₁₆], which is obtained from the database, is used to control the joint motion of the robot, where i = 1, 2,…13. The dynamic simulation of the corresponding biped robot walking is shown in Figure 10, and the black dots in the figure are the trajectory of $r_{b} (n)$ .

Figure 10.

Dynamic simulation of robot walking.

According to the above simulation data, the actual robot walking experiment is shown in Figure 11.

Figure 11.

Video frames of robot walking.

In the experiment, we set $d = 100 cm$ , $\hat{η} = 0.5$ , $T_{S} = 1.6 s$ , $k_{P} = 0.875$ , $T_{I} = 3.2 s$ , $T_{L} = 1.6 s$ , and we set $y_{AZR} (n)$ along the $η = 0. 8$ sideline during the SSP. The initial conditions for walking are $e_{0} = e_{- 1} = 0$ and $η_{0} = \hat{η} = 0.5$ . The $η$ value of the AZR controller is stable at $0.8$ after the third step. The calculation data are presented in Table 3, where the data in parentheses in the column titled s_i correspond to the movement distance of the swinging leg in front of and behind the support foot, $E_{s}^{η}$ is the energy consumption calculated according to the actual $η_{i}$ , and $E_{s}^{0.5}$ is the energy consumption value when $\hat{η} = 0.5$ . It is calculated that $J_{s}^{η} (10) = 1.929$ and $J_{s}^{0.5} (10) = 2.195 J$ ; that is, the AZR controller adjusts $η_{i}$ to reduce the energy consumption by $12.13 %$ .

Table 3.

GSYN data for robot walking, $d = 100 cm$ .

i	$s_{i} (cm)$	e_i	$Δ η_{i}$	$η_{i}$	$E_{s}^{η} (mJ)$	$E_{s}^{0.5} (mJ)$
	3(0+3)	0.25	0	0.5	90.5	90.5
2	7(3+7)	0.09	0.219	0.7	111.0	127.7
3	10(7+10)	0.03	0.079	0.8	153.7	176.0
4–10	10(10+10)	0.02	0.02	0.8	178.4	200.3
11	7(10+3)	0.01	0.009	0.8	153.7	176.0
12	3(7+3)	0.03	0.029	0.8	104.3	127.7
13	3(3+0)	0.02	0.016	0.8	65.2	94.8

GSYN: gait synthesis.

Conclusion

In this article, we have proposed a gait planning algorithm for biped robots that improve the 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:

The motor-driven robot has a dispersed mass, and the established five-mass robot model describes the physical characteristics of the robot better than other models. To balance the great computational complexity caused by the multimass model, we designed an energy consumption function with load torque as the main parameter; it could calculate the energy consumption index of the biped robot quickly during the simulation.

In the process of gridding the parameter space with the GPO algorithm, it is required that no grid has multiple peaks. If the appropriate mesh gap is selected, this method is suitable for optimizing most systems. During gradient approximation with the GPO algorithm, the step size of the successive approximations is determined according to the accuracy requirements. When the mesh size is large, the mesh can be meshed again to accelerate the approximation process. GPO requires many calculations during the optimization process. Using a GPU to accelerate the calculation can help determine the optimal value quickly. With the advancement of computer science, compared with many random optimization algorithms, our algorithm will show increasingly superior performance.

Given the walking distance of the robot, the GSYN algorithm plans the whole motion process, including the starting step, the intermediate steps, and the stopping step, while minimizing the energy consumption. In the process of robot walking, the algorithm adjusts the gait parameters according to the feedback ZMP data, and the AZR is set to overcome perturbations due to modeling and environmental errors. The algorithm can optimize the tradeoff between low energy consumption and high robustness and solve the problem of walking in biped robots with highly nonlinear characteristics.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Lu Zhiqiang

References

Saeedvand

Jafari

Aghdasi

, et al. A comprehensive survey on humanoid robot development. Knowl Eng Rev 2019; 34(20): 1–18.

Robio

Valero

Llopis-albert

. A review of mobile robots: concepts, methods, theoretical framework, and applications. Int J Adv Robot Syst 2019; 16: 1–22.

Paolillo

Cherubini

Keith

, et al. Toward autonomous car driving by a humanoid robot: a sensor-based framework. In: 2014 IEEE-RAS international conference on humanoid robots, Madrid, Spain, 18–20 November 2014, pp. 451–456. IEEE. ISBN: 978-1-4799-7174-9.

Koolen

Bertrand

Thomas

, et al. Design of a momentum-based control framework and application to the humanoid robot Atlas. Int J Humanoid Robot 2016; 13(1): 1650007.

Boutin

Eon

Zeghloul

, et al. From human motion capture to humanoid locomotion imitation application to the robots HRP-2 and HOAP-3. Robotica 2011; 29(2): 325–334.

Esfanabadi

Rostami

Rahmati

, et al. The average speed of motion and optimal power consumption in biped robots. Knowl Eng Rev 2019; 34(e25): 1–13.

Juang

Yeh

. Multiobjective evolution of biped robot gaits using advanced continuous ant-colony optimized recurrent neural networks. IEEE Trans Cybern 2018; 48(6): 1910–1922.

Lin

Yang

, et al. Rhythmic-reflex hybrid adaptive walking control of biped robot. J Intell Robot Syst 2019; 94: 603–619.

Wang

, et al. Hybrid CPG–FRI dynamic walking algorithm balancing agility and stability control of biped robot. Autonom Robot 2019; 43: 1855–1865.

10.

Kajita

Kanehiro

Kaneko

, et al. The 3D linear inverted pendulum mode: a simple modeling for a biped walking pattern generation. In: 2001 IEEE/RSJ international conference on intelligent robots and systems, Maui, HI, USA, 29 October–3 November 2001, pp. 239–246. IEEE. ISBN: 0-7803-6612-2.

11.

Hong

Kim

, et al. Real-time walking pattern generation method for humanoid robots by combining feedback and feedforward controller. IEEE Trans Ind Electron 2014; 61(1): 355–364.

12.

Akhond

Herzig

Wegiriya

, et al. A method to guide local physical adaptations in a robot based on phase portraits. IEEE Access 2019; 7: 78830–78841.

13.

Kang

, et al. Adaptive time-delay balance control of biped robots. IEEE Trans Ind Electron 2020; 67(4): 2936–2944.

14.

Shafii

Abdolmaleki

Lau

, et al. Development of an omnidirectional walk engine for soccer humanoid robots. Int J Adv Robot Syst 2015; 12: 193.

15.

Shimmyo

Sato

Ohnishi

. Biped walking pattern generation by using preview control based on three-mass model. IEEE Trans Ind Electron 2013; 60(11): 5137–5147.

16.

Luo

Chen

. Biped walking trajectory generator based on three-mass with angular momentum model using model predictive control. IEEE Trans Ind Electron 2016; 63(1): 268–276.

17.

Luo

Chen

. Quasi-natural humanoid robot walking trajectory generator based on five-mass with angular momentum model. IEEE Trans Ind Electron 2018; 65(4): 3355–3364.

18.

Kuo

. Choosing your steps carefully. IEEE Robot Autom Mag 2007; 14(2): 18–29.

19.

Hong

Park

Kim

. Stable bipedal walking with a vertical center-of-mass motion by an evolutionary optimized central pattern generator. IEEE Trans Ind Electron 2014; 61(5): 2346–2355.

20.

Shin

Kim

. Energy-efficient gait planning and control for biped robots utilizing vertical body motion and allowable ZMP region. IEEE Trans Ind Electron 2015; 62(4): 2277–2286.

21.

Dau

Chew

Poo

. Achieving energy- efficient bipedal walking trajectory through GA-based optimization of key parameters. Int J Humanoid Robot 2009; 6(4): 609–629.

22.

Elhosseinia

Haikal

Badawy

, et al. Biped robot stability based on an A-C parametric whale optimization algorithm. J Comput Sci 2019; 31: 17–32.

23.

. Fuzzy dynamic gait pattern generation for real-time push recovery control of a teen-sized humanoid robot. IEEE Access 2020; 8: 36441–36453.

24.

Wright

Jordanov

. Intelligent approaches in locomotion: a review. J Intell Robot Syst 2015; 80(2): 255–277.

25.

Bertram

JEA

. Constrained optimization in human walking: cost minimization and gait plasticity. J Exp Biol 2005; 208(6): 979–991.

26.

Paparisabet

Dehghani

Ahmadi

. Knee and torso kinematics in generation of optimum gait pattern based on human-like motion for a seven-link biped robot. Multibody Syst Dynam 2019; 47: 117–136.

27.

Sellaouti

Stasse

Kajita

, et al. Faster and smoother walking of humanoid HRP-2 with passive toe Joints. In: 2006 IEEE/RSJ international conference on intelligent robots and systems, Beijing, China, 9–15 October 2006, pp. 4909–4914. IEEE. ISBN: 1-4244-0258-1.

28.

Mobinipour

. Gate planning and control of a nine-link biped robot with toe joints via impact effects. Int J Humanoid Robot 2019; 16(1): 1950001.

29.

Vukobratovic

Borovac

. Zero-moment point - thirty five years of its life. Int J Humanoid Robot 2004; 01(01): 157–173.

30.

Srinivasan

Ruina

. Computer optimization of a minimal biped model discovers walking and running. Nature 2006; 439(5): 72–75.

31.

Lee

Park

Lee

, et al. Minimizing energy consumption of parallel mechanisms via redundant actuation. IEEE/ASME Trans Mechatron 2015; 20(6): 2805–2812.

32.

Lee

Stonier

Kim

, et al. Modifiable walking pattern of a humanoid robot by using allowable ZMP variation. IEEE Trans Robot 2008; 24(4): 917–925.

33.

Kim

. Multi-axis force-torque sensors for measuring zero-moment point in humanoid robots: a review. IEEE Sensors J 2020; 20(3): 1126–1141.

34.

Dehghani

Fattah

Abedi

. Cyclic gait planning and control of a five-link biped robot with four actuators during single support and double support phases. Multibody Syst Dynam 2015; 33: 389–411.