Dynamic walking control of humanoid robots combining linear inverted pendulum mode with parameter optimization

Abstract

To improve the robustness of biped walking, a model parameters optimization method based on policy gradient decent learning is presented. For the linear inverted pendulum mode-based model parameters optimization, firstly, select the input parameters of the inverted pendulum model and the torso attitude parameters of the robot as the correction variables and establish the correction equation. Then, using the tracking errors of center of mass (CoM) of the robot and the errors of the robot posture relative to the upright state of the body to establish the fitness function. According to the fitness function, the gain coefficients in the model parameters correction equation are optimized by using the strategy gradient learning method, and the modified gain parameters are substituted into the model parameters correction equation to obtain the correction amount. By applying the model parameters optimization strategy, the robot can quickly and in real time adjust the body posture and walking patterns under unknown disturbances, hence, the walking robustness can be enhanced. Simulation and experiments on a full-body humanoid robot NAO validate the effectiveness of the proposed method. The experiments show that the optimized model yields a more controlled, robust walk on NAO robot and on various surfaces without additional manual parameters tuning.

Introduction

For humanoid robots, stable and robust biped walking gait generation is very important. Many studies for dynamic biped walking are based on a simple humanoid robot as simplified dynamics models, ^{1 –6} like three-dimensional linear inverted pendulum mode (3D-LIPM), cart-table model, and so on. Using such an abstract model, there are various ways to generate appropriate gait patterns. Humanoid walking is commonly realized by planning the center of mass (CoM) trajectories, so that the resultant zero moment point (ZMP) trajectory follows a desired ZMP trajectory, which is normally determined by predefined foot positioning. These models based on control methods can realize biped walking under ideal environment. However, in complex condition, the balanced walking of a humanoid robot requires compensation in real time to guarantee the dynamic stability. To realize robust and dynamic walking, many strategies have been proposed to compensate for the nonzero variations in momentum ^7,8 or regulate CoM through sensory feedback. ^{9 –11}

Biped dynamic walking under disturbances is still an open question. To achieve external environment disturbance rebalancing, some well-known methods to stabilize walking using sensor feedback are proposed, ^{7 –9,12,13} including regulation of the states of CoM or by modification of the predefined foot positions online through sensory feedback. ^{10,11,14 –17} Considering humanoid momentum compensation, Englsberger and Ott ⁷ proposed a walking method based on a capture point, which integrates CoM vertical motion and angular momentum change. Ugurlu et al. ¹² created a strategy to rotate a humanoid’s upper body to compensate for undesired yaw movement when the robot walks. Chang et al. ⁸ estimated the robot state and regulated the angular momentum of CoM and steps to recover from a push. Castano et al. ¹⁰ proposed an online walking control that replanned the gait pattern based on foot placement control using the actual CoM state feedback. Wieber ^11,15 proposed an online walking gait generation approach, with both adjustable foot positioning and step duration based on the linear model predictive control scheme. Based on the preview control scheme, Nishiwaki and Kagami ¹⁶ proposed strategies to adjust future foot positioning and ZMP trajectory. Stephens ¹⁷ utilized an integral control method to maintain standing balance and viewed the disturbance as impulsive. Liu and Atkeson ¹⁸ studied the balance of a robot standing based on a trajectory library, which is time-consuming.

Open loop walk can easily be disturbed by external influences like hitting obstacles or more irregular terrains. Closed loop online modulation of the controller is the important strategy to solve the disturbance affect question during biped walking. Graf and Röfer ¹⁹ presented a closed loop 3D-LIPM-based gait generation method for the RoboCup standard platform. By adding the sensor feedback, robot can reach faster walking speeds and be robust to disturbances. On one hand, this allows reacting to unexpected disruptions such as unevenness of the ground or the forces exerted on the robot. On the other hand, it compensates the forces and hardware characteristics that were not considered within the model used. Using sensor feedback to observe the state of the robot and to adjust the inverted pendulum model to the observed state in order to improve the stability of the walk. Their approach not only allow compensating the inaccuracies in the model but also allows for reacting to external perturbations effectively, like with robustness against perturbations such as forces exerted on the robot. Urbann and Hofmann ⁵ presented a reactive stepping algorithm based on the generation of walking motions with sensor feedback. The closed form calculation of foot placement modifications, such as disturbances of the CoM position, could be balanced with negligible ZMP deviations.

In order to improve the robustness of the previous work on the biped dynamic walking based on the LIPM model, the main contribution of this article is presenting an improved closed loop dynamic gait generation method based on the previous work of Graf et al. for humanoid robots combining LIPM with parameter optimization. For the presented LIPM-based model parameters optimization, firstly, select the input parameters of the inverted pendulum model and the torso attitude parameters of the robot as the correction variables and establish the correction equation. Then, using the tracking errors of CoM of the robot and the errors of the robot posture relative to the upright state of the body as the fitness index of the robot to the current environment, and establish the fitness function. According to the fitness function, the gain coefficients in the model parameters correction equation are optimized by using the strategy gradient learning method, and the modified gain parameters are substituted into the model parameters correction equation to obtain the correction amount. By applying the model parameters optimization strategy, the robot can quickly and in real time adjust the body posture and walking pattern under unknown disturbance, hence, the walking robustness can be enhanced. Simulation and experiments on a full-body humanoid robot NAO validate the effectiveness of the proposed method. The experiments show that the optimized model yields a more controlled, robust walk from NAO robots and on various surfaces without additional manual parameters tuning.

The rest of this article is organized as follows. Second section presents the LIPM-based biped gait generation, including the LIPM model, the CoM trajectory and foot trajectory planning, and observation of the CoM by Kalman filter. Third section introduces the design method of model parameters corrector based on policy gradient decent learning method. Fourth section verifies the real-time performance and validity of the presented compensation algorithm by experiments. Fifth section provides the conclusions and discussion.

LIPM-based gait generation

The dynamics of a humanoid robot can be approximated as a linear inverted pendulum with its mass concentrated at CoM and supported at ZMP (Figure 1). In this model, the robot’s leg is assumed as a weightless scalable limb and its kinematics constraints are not considered. When not considering the constraint of the limited area of supporting polygon and neglecting the vertical moment of the CoM, the dynamics of LIPM in the x–z plane can be formulated as a strictly proper dynamic system as follows

${\dot{x}}_c=A{x}_c+B{u}_x$

1

$x_z=C{x}_c$

2

Figure 1.

Linear inverted pendulum model.

where $x_c={[x_c,{\dot{x}}_c,{\ddot{x}}_c]}^{\text{T}}$ contains the position, velocity, acceleration of CoM along the x-axis in ∑ _w , u_x is the jerk of x_c, $p_c={[x_c,y_c,z_c]}^{\text{T}}$ is the position of CoM in ∑ _w and $p_z={[x_z,y_z,z_z]}^T$ is the position of ZMP in ∑ _w , let g be the gravitational acceleration, and $A=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right],B=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right],C=\left[\begin{array}{ccc}1& 0& -\frac{z_c-z_z}{g}\end{array}\right]$ .

Given a set of future foot positioning, a reference ZMP trajectory can be generated on fulfillment of the dynamic stability requirement, and the desired CoM trajectory can be generated online using preview control (or the build in PID control) based on the LIPM. Using the inverse kinematics to obtain joint commands from the CoM trajectory, an open-loop humanoid walking control method can be realized. However, this just can realize biped walking under ideal environment, like flat terrain walking without external disturbance.

The LIPM-based trajectory planning

The position and velocity of the CoM relative to the origin of the inverted pendulum are given by Graf and Röfer ¹⁹ as follows

$x(t)=x_0\cdot \text{cosh}(k\cdot t)+{\dot{x}}_0\cdot 1/k\cdot \text{sinh}(k\cdot t)$

3

$\dot{x}(t)=x_0\cdot k\cdot \text{sinh}(k\cdot t)+{\dot{x}}_0\cdot \text{cosh}(k\cdot t)$

4

where $k=\sqrt{g/h}$ , g is the gravitational acceleration, h is the height of the CoM above the ground, x₀ and ${\dot{x}}_0$ are the position and the velocity of the CoM relative to the origin of the inverted pendulum at t = 0, respectively.

The support leg is constantly switching during walking, so the origin of the pendulum need to be constantly switching because the origin of the pendulum is on the support foot. As shown in Figures 1(a) and 2, Q is the coordinate system of the robot, and its origin is located between the feet. To eliminate the double support phase, we need to calculate the switching time point of the support leg. t = 0 is defined as a starting time point for altering the support leg, when t = 0, the y-component of the velocity is 0 ( ${({\dot{x}}_0)}_y=0$ ). The position of the CoM at this point ${(x_0)}_y$ is an arbitrary parameter and has a value of greater or lower than 0 depending on the active support leg. Hence, the function of CoM’s position in y-direction during the range $t_b\le t\le t_e$ (which means a single support phase starts at t = t_b and ends at t = t_e) is as follows

$x_y(t)={(x_0)}_y\cdot \text{cosh}(k\cdot t)$

5

Figure 2.

xy-cross section showing the step size $\overline{s}$ and the inverted pendulum origins r and $\overline{r}$ .

The inverted pendulum motion diagram in the y–z plane is shown in Figure 3. The swing foot should have crossed a distance of ${\overline{r}}_y+{\overline{s}}_y-r_y$ with respect to the support leg (Figure 2) at the end of the single support phase. The switching time point of the support leg can be determined by the ending time of the single support phase t_e as follows

${(x(t_e))}_y-{(\overline{x}({\overline{t}}_b))}_y={\overline{r}}_y+{\overline{s}}_y-r_y$

6

${(\dot{x}(t_e))}_y\text{=}{(\overline{\dot{x}}({\overline{t}}_b))}_y$

7

Figure 3.

The inverted pendulum motion in y–z plane.

where $\overline{x}(\overline{t})\text{ and }\overline{\dot{x}}(\overline{t})$ are position and velocity of the CoM relative to the next pendulum origin. t_e and ${\overline{t}}_b$ can be calculated using an iterative method. ²⁰

The support leg switching in x-direction have the same regularity with the y-direction as follows

${(x(t_e))}_x-{(\overline{x}({\overline{t}}_b))}_x={\overline{r}}_x+{\overline{s}}_x-r_x$

8

${(\dot{x}(t_e))}_x={(\overline{\dot{x}}({\overline{t}}_b))}_x$

9

To ensure the optimum original of the inverted pendulum in the next single support phase, and the CoM is located just above the origin, we design the conditions of ${\overline{r}}_x=0\text{ and }{(\overline{x}(0))}_x=0$ . When the CoM has the position $x_{t_b}$ and velocity ${\dot{x}}_{t_b}$ relative to the origin of Q at the beginning of a single support phase, r_x, ${(x_0)}_x,\text{ and }{({\dot{x}}_0)}_x$ can be computed by the following equations, so that all the pendulum parameters can be determined.

$\begin{array}{l}{r}_x+{(x_0)}_x\cdot C(k\cdot t_b)+{({\dot{x}}_0)}_x\cdot \frac{S(k\cdot t_b)}{k}={(x_{t_b})}_x\\ {(x_0)}_{x}k\cdot S(k\cdot t_b)+{({\dot{x}}_0)}_x\cdot C(k\cdot t_b)={({\dot{x}}_{t_b})}_x\\ {(x_0)}_{x}k\cdot S(k\cdot t_e)+{({\dot{x}}_0)}_x\cdot C(k\cdot t_e)-{({\overline{x}}_0)}_x\cdot k\cdot S(\overline{k}\cdot {\overline{t}}_b)\\ -{({\overline{\dot{x}}}_0)}_x\cdot C(\overline{k}\cdot {\overline{t}}_b)=0\\ r_x+{(x_0)}_x\cdot C(k\cdot t_e)+{({\dot{x}}_0)}_x\cdot \frac{S(k\cdot t_e)}{k}-{({\overline{x}}_0)}_x\cdot C(\overline{k}\cdot {\overline{t}}_b)\\ -{({\overline{\dot{x}}}_0)}_x\cdot \frac{S(\overline{k}\cdot {\overline{t}}_b)}{k}={\overline{r}}_x+{\overline{s}}_x\end{array}$

10

where C( ) and S( ) are the abbreviation of cosh( ) and sinh( ).

The $x_z(t)$ in z-direction can be defined as follows

$x_z(t)=h+l_z\cdot l(t)$

11

$\phi (t)=\frac{t-t_b}{t_e-t_b}$

12

$l(t)=\{\begin{array}{c}\frac{1}{2}\cdot \left(1-\text{cos}\left(\frac{2\pi \cdot (\phi (t)-{\phi }_s)}{{\phi }_D}\right)\right)\\ 0\end{array},\begin{array}{c}{\phi }_s\le \phi (t)\le {\phi }_s+{\phi }_D\\ \text{else}\end{array}$

13

where l_z is the lift height in z-direction of CoM, φ_s and φ_D are the proportional constant coefficients, which describe the beginning time and the duration of movement with the relationship of $2{\phi }_s+{\phi }_D=1,0\le {\phi }_s,{\phi }_D\le 1$ .

When the movement of the CoM trajectory relative to the origin of the pendulum motion is determined, the movement of the support leg with respect to the CoM will also be determined, which can be expressed as

$p_{\text{sup}\_\text{foot}}\text{=}{\left[\begin{array}{ccc}-x_x(t)& -x_y(t)& -x_z(t)\end{array}\right]}^{\text{T}}$

14

The movement of the swing foot with respect to the robot coordinate system Q can be divided into the xy-direction movement and the z-direction movement as follows

$\begin{array}{l}{p}_{\text{sw}\_\text{foo}{\text{t}}_{xy}}=(s_{xy}+{\overline{s}}_{xy})\cdot m(t)-s_{xy}\text{+}{f}_y\\ m(t)=\{\begin{array}{c}0\\ \frac{1}{2}\cdot \left(1-\text{cos}\left(\frac{\pi \cdot (\phi (t)-{\phi }_s)}{{\phi }_D}\right)\right)\\ 1\end{array},\begin{array}{c}\phi (t)>{\phi }_s\\ {\phi }_s\le \phi (t)\le {\phi }_s+{\phi }_D\\ \phi (t)<{\phi }_s+{\phi }_D\end{array}\end{array}$

15

$p_{{\text{sw\_foot}}_{\text{z}}}=s_z\cdot l(t)$

16

where s_xy and ${\overline{s}}_{xy}$ are the current and next step’s size in xy-direction, s_z is the lift height in z-direction, and f_y is the foot position which only contains y-component (positive or negative is determined by the left or right foot).

Meanwhile, we need to set the angular velocity of the support foot and the swing foot as follows

$\begin{array}{l}{\omega }_{\text{sup\_foot}}\text{=}{\left[\begin{array}{ccc}0& 0& 0\end{array}\right]}^{\text{T}}\\ {\omega }_{\text{sw\_foot}}\text{=}{\left[\begin{array}{ccc}{\omega }_x\cdot l(t)& \omega y\cdot l(t)& 0\end{array}\right]}^{\text{T}}\end{array}$

17

Finally, combining target foot positions relative to the CoM and CoM position relative to torso by computing current robot’s stance, the joint angles can be calculated by inverse kinematics.

Observing the CoM by Kalman filter

Based on Craf’s previous work, ¹⁹ for determining an observed position of the CoM, two Kalman filters are used for the x- and y-components of the position and the velocity of the CoM. The expected CoM position x_ei with respect to Q is computed using x(t) and the error with the measured position x_mi is calculated to get update the Kalman as follows

$x_{ei}=r+x(t_i)$

18

$\Delta x_i=x_{mi}-x_{ei}$

19

The state quantities and covariance matrices of the Kalman prediction phase are as follows

$u_i\prime ={\left[\begin{array}{cccc}{(x_{ei})}_x& {(x_{ei})}_y& {(\dot{x}(t_i))}_x& {(\dot{x}(t_i))}_y\end{array}\right]}^{\text{T}}$

20

${\Sigma }_i\prime =A\cdot {\Sigma }_{i-1}{A}^{{\rm T}}+{\Sigma }_{\epsilon i}$

21

where $A=\left[\begin{array}{cccc}1& 0& \Delta t_i& 0\\ 0& 1& 0& \Delta t_i\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right],\Delta t_i=t_i-t_{i-1}$ , and Σ _εi is a process noise covariance matrix.

The state quantity and covariance matrix of the Kalman update phase are as follows

$K_i={\Sigma }_i\prime \cdot C^{\text{T}}\cdot {(C\cdot {\Sigma }_i\prime \cdot C^T+{\Sigma }_{mi})}^{-1}$

22

$u_i={[\begin{array}{cccc}{(x_{fi})}_x& {(x_{fi})}_y& {({\dot{x}}_{fi})}_x& {({\dot{x}}_{fi})}_y\end{array}]}^T=u_i\prime +K_i\cdot \Delta x{}_i{}^{}$

23

${\Sigma }_i={\Sigma }_i\prime -Ki\cdot C\cdot {\Sigma }_i\prime$

24

where K_i is the Kalman gain, $C=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right]$ , and Σ _mi is the measured noise covariance matrix. Smaller Σ _εi and larger Σ _mi with less correction effectiveness for the inverted pendulum model parameters.

The filtered position x_fi and the filtered velocity ${\dot{x}}_{fi}$ of the Kalman filter that estimates the true position and velocity of the CoM are used to redetermine the parameters of the inverted pendulum of the current single support phase. Since the prediction of the natural motion of the CoM is based on the pendulum parameters, the main purpose of the correction is to improve the prediction in the next iteration. So the pendulum parameters r, x₀, ${\dot{x}}_0$ , and t_i are redetermined to find a pendulum function that fits to the estimated position and the estimated velocity of the CoM

$r\prime +x\prime (t_i\prime )=x_{fi}$

25

$\dot{x\prime }(t_i\prime )={\dot{x}}_{fi}$

26

The motion equation in the y-direction is corrected as follows

${(x_f)}_y-r_y\prime ={(x_0)}_y\prime \cdot \text{cos\hspace{0.17em}}h(k\prime \cdot t_i\prime )$

27

${({\dot{x}}_f)}_y={({\dot{x}}_0)}_y\prime \cdot k\cdot \text{sin\hspace{0.17em}}h(k\prime \cdot t_i\prime )$

28

The values of ${(x_0)}_y\prime$ and t_i′ can be obtained by solving equations as follows, then the parameters of the inverted pendulum model are obtained by the previous inverted pendulum parameter calculation step

${(x_0)}_y\prime =\sqrt{{({(x_{fi})}_y-r_y\prime )}^2-\frac{{(x_f{}_i)}_y{}^2}{{k\prime }^2}}$

29

$t_i\prime =\frac{1}{k\prime }\cdot \text{arcsin\hspace{0.17em}}h\left(\frac{{({\dot{x}}_f)}_y}{k\prime \cdot {(x_0)}_y}\right)$

30

Learning-based model parameters correction

In this work, a model parameters optimization method based on policy gradient decent learning is proposed. Usually, there are two methods of using policy gradient learning in the inverted pendulum model. One approach is to update the input parameters of the inverted pendulum every time a policy set is generated during the process of parameter optimization. When a robot via the virtual linkage model takes one step using this policy set, the corresponding value of the objective function and the parameter set of the next iteration are obtained. Once the iteration number reaches the preset iteration number N_iter, the parameter set is taken as the input parameters of the inverted pendulum, and the output parameters of the inverted pendulum are passed on to the robot. This method can directly adjust the gait of a robot. However, the computations of the robot model and the inverted pendulum model are involved in each iteration, which means N_iter previews for each real step of the robot. The previews will take a long time, preventing the robot’s pace from being smooth, and therefore the approach is not suitable for real-time robot walking.

The other approach is that after a robot takes a step, the objective function is used to evaluate the state of this step, and then the input parameters for the inverted pendulum are compensated. At this time, the method of policy gradient learning does not directly optimize the input parameters of the inverted pendulum but optimizes the gain coefficient of the compensation amount, thereby indirectly adjusting the gait of the robot. The second approach separates the gradient descent method for policy search from the virtual linkage model for the robot, and the computation speed is greatly improved. It is thus suitable for the real-time gait adjustment of the robot. However, the relational function for the parameter compensation of the inverted pendulum and the gain coefficient need to be established to find the internal relationship between the compensation amount and the virtual linkage model for robots, with a certain degree of increase in difficulty.

The overall procedure of humanoid walking with learning-based parameters correction can be summarized by the diagram of Figure 4. For determining an observed position of the CoM, a Kalman filter will be used for the x- and y-components of the position and the velocity of the CoM in the inner feedback loop. The filtered position x_f and the filtered velocity ${\dot{x}}_f$ of the Kalman filter that estimates the true position and velocity of the CoM are used to redetermine the parameters of the inverted pendulum of the current single support phase. For the parameters correction outside loop of policy gradient learning (PGL) compensator, select the input parameters of the inverted pendulum model and the torso attitude parameters of the robot as the correction variables and establish the compensator.

Figure 4.

The diagram of our humanoid walking control with model parameters correction.

Model parameters correction equation

There have been many parameters in the LIPM model which influence the biped walking, such as, the step stride s_x and s_y, the step height s_z, and the height of the inverted pendulum h. In addition, the robot is a multidimensional model with high degrees of freedom (DOFs), it cannot be completely described by the inverted pendulum, such as the trunk posture θ_B, and so on. Select the variables s_x, s_y, ${\theta }_B{}_x,\text{ and }{\theta }_B{}_y$ as the key parameters to effect the walking pattern and construct the parameters correction equation as follows

$\Delta s_x=K_1\cdot \frac{1}{N}{\displaystyle \sum _{i=1}^N\left(x_{f,x,i}-x_{e,x,i}\right)}+K_3\cdot \frac{1}{N}{\displaystyle \sum _{i=1}^N\left({\theta }_{B,y,i}-{\theta }_{B,y,i}^{\text{ref}}\right)}$

31

$\Delta s_y=K_2\cdot \frac{1}{N}{\displaystyle \sum _{i=1}^N\left(x_{f,y,i}-x_{e,y,i}\right)}+K_4\cdot \frac{1}{N}{\displaystyle \sum _{i=1}^N\left({\theta }_{B,x,i}-{\theta }_{B,x,i}^{\text{ref}}\right)}$

32

$\Delta {\theta }_{B,x}=K_5\cdot \frac{1}{N}{\displaystyle \sum _{i=1}^N\left(p_{\text{LHip},z,i}-p_{\text{RHip},z,i}\right)}$

33

$\Delta {\theta }_{B,y}=K_6\cdot \frac{1}{N}{\displaystyle \sum _{i=1}^N\left(p_{\text{SuppFoot},x,i}-p_{\text{Head},x,i}\right)}$

34

where $N=\frac{t_e-t_b}{\Delta T}$ is the interpolation number of steps in the single support phase, x_fi is the estimate value of CoM with Kalman filter, x_ei is the ideal value of the CoM, s is the step stride, and Δs is the correction of step size. θ_B is the actual torso angle inclination and Δθ_B is the correction of torso angle. ${\theta }_{Bx}^{\text{ref}}$ is the torso’s tilt angle when upright. ${(p_{\text{RHip}})}_z\text{ and }{(p_{\text{LHip}})}_z$ is the z-direction’s displacement of the right leg and left leg of the robot, respectively. ${(p_{\text{Head}})}_x\text{ and }{(p_{\text{SuppFoot}})}_x$ is the x-direction’s displacement of the head and support leg of the robot, respectively. $\underset{\_}{K}=\left\{K_1,\mathrm{..}{K}_6\right\}$ are the gain coefficients.

This parameter modification function is designed based on the correction of the foothold of the robot in the x- and y-directions and the correction of the straightness of the torso. The correction of the step length along the x-direction, $\Delta s_x$ , is related to the error in the x-axis centroid and the error in the y-axis torso angle because the y-axis torso angle represents the forward and backward inclination of the body. The correction of the step length along the y-direction, $\Delta s_y$ , is related to the error in the y-axis centroid and the error in the x-axis torso angle because the x-axis torso angle represents the left and right inclination of the body. The correction of the torso angle on the x-direction, $\Delta {\theta }_{B_x}$ , is related to the height difference between the left and right hip joints, which is because, when a robot needs to keep both feet on the ground or perform the walking task, the left and right tilts of the body will cause a certain degree of height difference between the left and right hip joints, and for a robot without a driving waist joint, the inclination of the body is mainly caused by the height difference of the two hip joints. The correction of the torso angle in the y-direction, $\Delta {\theta }_{B_y}$ , is related to the difference in distance between the head joint and the support foot in the x-direction, which is because, for a robot without a driving waist joint, the forward and backward inclination of the body of a robot is approximately equivalent to the linkage rotation of the robot with the support foot as the fulcrum, with the support foot as the origin of the linkage, and with the head as the end of the linkage. Therefore, the distance between the support foot and the head in the x-direction represents the degree of inclination of the body.

The output of the compensator is used to correct the inverted pendulum model and the robot model as follows

${\overline{s}}_{xy}=s_{xy}+\Delta s$

35

${\theta }_{B,}{}_i={\theta }_{B,i-1}+\Delta {\theta }_B$

36

where ${\overline{s}}_{xy}$ is the input parameter of the inverted pendulum model which represents the step size of the next single foot support phase and every single foot support phase corrects only once. ${\theta }_{B,}{}_i$ is the torso angle for each frame during the single-leg support phase and each single foot support stage will be compensated times.

Policy gradient reinforcement learning

Policy gradient learning method is a powerful machine learning algorithm to find the local optimal policy when there is no analytic mapping from the performance index to the pending parameters or when the gradient of the cost function cannot be directly calculated. This method has already been adopted in many robotics applications, especially for learning the walking of legged robot. ^{21 –23} In this work, the parameters of inverted pendulum is optimized by using the method of gradient learning. As we know, when the robot is disturbed by unknown disturbance, if after adjusting, it can maintain an high accuracy of CoM tracking and good body upright posture, then the robot’s with high adaptability to the disturbance. In the correction equation, $K_1\cdots K_6$ are the gains to control the output of the compensator. The fitness function $F(\underset{\_}{K})$ is set as follows, the selection of their value is to minimize the following performance

$F(\underset{\_}{K})={\alpha }_x\left(\left|\Delta s_x\right|+\left|\Delta {\overline{x}}_x\right|\right)+{\alpha }_y\left(\left|\Delta s_y\right|+\left|\Delta {\overline{x}}_y\right|\right)+{\beta }_x\left(\left|\Delta {\theta }_{B,x}\right|+\left|\Delta {\overline{\theta }}_{B,x}\right|\right)+{\beta }_y\left(\left|\Delta {\theta }_{B,y}\right|+\left|\Delta {\overline{\theta }}_{B,y}\right|\right)$

37

$\Delta {\overline{x}}_x=\frac{1}{N}{\displaystyle \sum _{i=1}^N\left(x_{f,x,i}-x_{e,x,i}\right)}$

38

$\Delta {\overline{x}}_y=\frac{1}{N}{\displaystyle \sum _{i=1}^N\left(x_{f,y,i}-x_{e,y,i}\right)}$

39

$\Delta {\overline{\theta }}_{B,x}=\frac{1}{N}{\displaystyle \sum _{i=1}^N\left({\theta }_{B,x,i}-{\theta }_{B,x,i}^{\text{ref}}\right)}$

40

$\Delta {\overline{\theta }}_{B,y}=\frac{1}{N}{\displaystyle \sum _{i=1}^N\left({\theta }_{B,y,i}-{\theta }_{B,y,i}^{\text{ref}}\right)}$

41

where $\underset{\_}{K}=\left\{K_1,\cdots ,K_6\right\}$ represents the set of gain parameters, α_x, α_y, β_x, and β_y, are the weight factors which satisfy ${\alpha }_x+{\alpha }_y=1\text{ and }{\beta }_x+{\beta }_y=1$ . The smaller the value of the fitness evaluation function, the higher the fitness of the robot in the gain parameter set.

Each term of the fitness function contains an absolute value of the compensation amount with the mean error, which is obtained by the linear addition of the absolute value of the output of the compensator, the absolute value of the centroid error, and the absolute value of the body inclination angle error, with different weights. The first two terms represent the following effect of the centroid, and the latter two terms represent the erecting effect of the body. In the process of gradient learning, if the value of the fitness function decreases gradually, the errors in the centroid following and the body inclination are gradually reduced, and the adaptability of the robot under the current parameter set is gradually increased. The corrections of the input parameters and the body inclination of the inverted pendulum cannot be too large, to ensure a reasonable range for allocation of the robot joint space. At the same time, it should be noted that the objective function is not directly associated with the robot walking time but a better erecting effect of the body indicates a longer implementation of centroid following in the robot.

Optimization learning of the parameters in the model

Policy gradient learning is used to optimize the gain parameter set of $K\text{=}\left\{K_1,\cdots ,K_6\right\}$ , which refers to the assigned value of each parameter in the parameter set of the fitness function, and the value of the fitness function is calculated based on the gain parameter set and the current state of the robot. Assuming that the objective function $F(\underset{\_}{K})$ is derivable for each parameter in $\underset{\_}{K}$ , the local optimal solution ${\underset{\_}{K}}^{*}$ can be obtained by calculating the gradient of $F(\underset{\_}{K})$ .

The specific procedures of the method of policy gradient learning are as follows.

Step 1

In the k th iteration, for the gain parameter set ${\underset{\_}{K}}^{k-1}$ obtained in the previous iteration, the partial derivatives of $F(\underset{\_}{K})$ at each parameter value within ${\underset{\_}{K}}^{k-1}$ are calculated, and n policies are randomly generated around ${\underset{\_}{K}}^{k-1}$ . The obtained policy set is represented as ${}_m{}^{}K{}^{k-1}$ (m = 1,…, n), with the number of policies n proportional to the searching space, and the equation for the generation of the policy set is

${}_m{}^{}\underset{\_}{K}{}^{k-1}={\underset{\_}{K}}^{k-1}+{}_m{}^{}\underset{\_}{\rho }$

42

where ${}_m{}^{}\underset{\_}{\rho }$ (m=1,…, n) refers to the disturbance set, each disturbance ρ_m in the disturbance set is randomly selected from the set of $\left\{-e_m,0,+e_m\right\}$ , and e_m represents the gain parameter of disturbance of the corresponding ρ_m.

Step 2

According to the values of −e_m, 0, +e_m for the disturbance ρ_m, the ${}_m{}^{}\underset{\_}{K}{}^{k-1}$ are divided into three corresponding groups: $G_{-{\epsilon }_m}$ , G₀, and $G_{+{\epsilon }_m},\text{ and }{}_m{}^{}\underset{\_}{K}{}^{k-1}$ is substituted into the fitness evaluation function to obtain the corresponding average of each group: ${\overline{F}}_{-{\epsilon }_m},{\overline{F}}_0,\text{ and }{\overline{F}}_{\text{+}{\epsilon }_m}$ .

Step 3

To calculate the approximate gradient value $\nabla {\underset{\_}{K}}^{k-1},\text{ if }{\overline{F}}_{-{\epsilon }_m}>{\overline{F}}_0\text{ and }{\overline{F}}_{\text{+}{\epsilon }_m}>{\overline{F}}_0,\nabla {\underset{\_}{K}}^{k-1}\text{=}0$ ; otherwise, $\nabla {\underset{\_}{K}}^{k-1}\text{=}{\overline{F}}_{\text{+}{\epsilon }_m}-{\overline{F}}_{-{\epsilon }_m}$ .

Step 4

With the orthogonal processing for $\nabla {\underset{\_}{K}}^{k-1}$ , the gradient value $\nabla {\stackrel{⌢}{\underset{\_}{K}}}^{k-1}$ is obtained by multiplying a fixed step factor η. The gradient value $\nabla {\underset{\_}{\stackrel{⌢}{K}}}^{k-1}$ is subtracted from the policy set ${\underset{\_}{K}}^{k-1}$ to yield the policy set K^k of this iteration, and then K^k is used for the next iteration.

Step 5

When the number of iterations reaches the preset value N_iter, the iteration ends.

Experiments and results

Humanoid robot

To validate the proposed method, NAO robot produced by Aldebaran robotics is used for the experiments (Figure 5(a)). The robot has 25 DOFs, with 11 DOFs for the two legs. The height of this robot is about 54 cm during walking and its weight is 4.35 kg. The sensors data and the walking control module update at a rate of 50 frames per second. Besides, the robot is equipped with an ×86 AMD GEODE 500 MHz CPU. The difficulty to develop walking approaches on this robot is that each leg weighs almost the same with its body and the links of the robot are made of plastic, which makes the dynamics of the robot differs a lot from the simplified dynamics model, and the disturbances brought with the deformation of plastic links are significant.

Figure 5.

Humanoid robot platform. (a) NAO robot. (b) The ODE-based simulation environment. (c) Coordinate systems of the humanoid robot.

Two coordinate systems are used (Figure 5(c), the first one is the world coordinate system ∑ _w , with its origin on the ground and its x-, y-axis formed a plane parallel to the ground; the second one is called the supporting coordinate system ∑_sup, which is attached on the supporting foot. Its origin is the origin point of the supporting foot with its x-axis pointing in the forward direction of the supporting foot, and y-axis to the outside. In this work, experiments on both the simulated and real full-body humanoid robot are developed to show the performance of proposed method. Due to the potential damages and time involved in PGL process on a real humanoid robot, the learning process is currently carried out in Webots simulator.

Simulated experiments

Two experimental setups are developed to test the robot performance and in both the experimental setups, the robot has no prior knowledge of the terrain. The first experimental setup is to let the robot disturbed by uncertain external force during walking. The impose point of the external force on the robot and the force moment, the size of the force and the duration of the force are all uncertain, which can verify the robustness of the robot to unknown disturbances during walking. The second experimental setup is to let the robot walk on different walking terrains to verify the walking robustness when the walking condition changes.

Experiements with only Kalman filter

Omnidirectional walking with only Kalman filter

Figure 6 shows an experiment that the robot turning counterclockwise during forward walk at a speed of 0.05 m/s and with only Kalman filter but no LIPM model parameters online correction. The robot can follow the specified pace size, and with good CoM trajectory tracking result, the ZMP always remain in the feet supporting polygons.

Figure 6.

Plot recorded when the robot walks both forward and turning counterclockwise. The plot on the top is the x–y plane data and the plot on the bottom is the CoM and ZMP trajectories. CoM: Centre of mass; ZMP: Zero moment point.

Flat terrain walking with unknown disturbance

The external impulse about 6.5 N happens at about 9 s, the direction is mainly along the y-axis positive direction. Figure 7 shows an experiment that the robot walk straight forward and with only Graf’s method of correcting the inverted pendulum parameters but no model parameters optimization method based policy learning. The robot swings and falls down after being hit by the disturbance. Since there is no online compensation to modify the state of the robot, the robot still follows the predefined footsteps and thus results in the same reference ZMP trajectory. Although with the Kalman filter, the robot can generate the appropriate CoM trajectory to track the reference ZMP, the CoM tracking requirement is however unable to be fulfilled with large external disturbance. The robot manages to swing for several steps and finally falls down.

Figure 7.

Plot recorded when the robot walks forward with only Graf’s method but no learning-based model parameters correction.

Experiments with parameters online correction

Flat terrain walking with unknown disturbance

The parameters correction is added outside loop of PGL compensator as shown in Figure 4 in the same experiment. By using the input parameters of the inverted pendulum model and the torso attitude parameters of the robot as the correction variables and establish the compensator, and optimizing the parameters of inverted pendulum by using the method of gradient learning. Figure 8 shows the snapshots of the motions performed in disturbance resistance. The external impulse about 6.5 N happens at about 5 s in the y-direction, the duration of force of about 0.5 s. The initial values and the weight factors are set as shown in Table 1. The robot begins to walk along the x-direction, after disturbance by the external force in the left direction, and with the online compensation, the robot immediately adjust the next few steps to adjust the pace, steps to the left, and through the trunk angle compensation, the robot gradually restores upright state.

Figure 8.

Snapshots of the motions performed in disturbance resistance.

Table 1.

FPC parameters.

Gain	Value
αx	0.3
αy	0.7
βx	0.6
βy	0.4
KO	$[000011]$
η	0.001
εj	0.001
Niter	480

The curve of the torso of the robot is shown in Figure 9. It can be seen from the plot that the robot is tilted to the right after being subjected to external force. After about 7s′ adjustment, the x-axis torso curve returns to the periodic fluctuation during normal walking. Since the robot inertial unit sensor cannot measure the z-axis torso angle, the z-axis torso angle compensation is not added to the parametric compensator design, so the z-axis torso angle is always zero. The expected and measured CoM trajectories and the measured ZMP trajectory are shown in Figure 10. When the robot is subjected to the right external force, the ZMP is shifted to the right and has a large fluctuation but always remains in the support area. The centroid has a good tracking effect, and the external disturbance has little effect on the centroid tracking. The recorded footstep of the robot in the world coordinate system is shown in Figure 11. After using the compensator, the pace of the robot changes due to external disturbances, and when the robot is subjected to right external forces, it will move to the right and the travel distance in the x-direction is slightly reduced. The value of the objective function of the strategy gradient learning method is shown in Figure 12 when the parameter compensator is called in the process of robot walking. The objective function curve shows a decreasing trend as a whole, indicating that the error of the centroid tracking and the error of the torso angle are gradually reduced, which verify the effectiveness of the presented model compensation.

Figure 9.

Body attitude of robot under external disturbance.

Figure 10.

The ideal CoM, measured CoM and measured ZMP trajectory under external disturbance. CoM: Centre of mass; ZMP: Zero moment point.

Figure 11.

The footstep of the robot in x–y plane under external disturbance.

Figure 12.

Fitness function curve of the policy gradient learning.

Irregular terrain adaptive walking

The second experiment setup is to test the ability of the presented method’s walking on irregular terrain. In the simulated experiment, a flat plate with size 0.5 × 0.5 × 0.016 m³ is laid in walking ground. During walking on the plate, the robot will be affected by the reaction of the plate. The parameter set initial value of the strategy gradient learning method is set as: ${\underset{\_}{K}}^0=\left\{0,0,0,0,1,2\right\}$ , η = 0.001, ${\epsilon }_j=0.001,\text{ and }{N}_{\text{iter}}=720$ . The weighting factor of the objective equation is set as: ${\alpha }_x=0.6,{\alpha }_y=0.4,{\beta }_x=0.6,\text{ and }{\beta }_y=0.6$ . Figure 13 shows the experiment snapshots of the robot walking up and down the plate. The robot can adjust the walking step size and trunk posture through the compensator in real time to successfully walking through the plate.

Figure 13.

Snapshots of the walking up and down the plate. (a) Without online parameters correction and (b) with online parameters correction.

Figure 14 shows the body attitude in pitch plane and roll plane during walking. The robot slides backward with shaking when walking up the plate. When walking on the plate, the robot’s x-axis torso curve returns to normal waveform, while the y-axis torso curve remains negative and fluctuates slightly, with fluctuations and slides in the forward direction when the robot walking down the plate. The ideal CoM of the robot, the measured CoM and the measured ZMP trajectory are shown in Figure 15. When the robot walks up and down the plate, the x-direction centroid tracking appears as a more obvious error, the actual CoM position lags behind the ideal CoM position. The ZMP trajectory has a large fluctuation while walking up and down the plate but always remain in the supporting polygon. In the world coordinate system, the foothold diagram of the robot with the compensator is shown in Figure 16. The yellow area is the walking on plate. When the robot moves up and down the plate, the pace of the robot changes with the change in the ground conditions.

Figure 14.

Body attitude angles.

Figure 15.

The ideal CoM, measured CoM, and measured ZMP trajectory during adaptive walking. CoM: Centre of mass; ZMP: Zero moment point.

Figure 16.

Footsteps in x–y plane during walking through the plate.

Real robot experiments

Two experimental setups are developed to test the real robot performance and in both the experimental steups the robot has no prior knowledge of the disturbance and the walking terrain. The video clips of the experiments can be found in the video attachment. The first experiment setup is to let the robot disturbed by the unknown external force as on Figure 17. The result of the experiment is shown in Figure 18. The robot is disturbed by an external force from the back with about 5 N during the robot walking at a speed of 0.12 m/s. After the disturbance, by applying the model parameter compensator, the robot can adjust the pace and body posture according to the target equation in real time, thus reducing the error of the CoM tracking and the inclination of the body. Several experiments show that the robot can withstand the external force of about 10 N in the y-direction in the simulation environment, and can withstand the force of about 6 N in the x-direction. The real experiment can resist about 5 N force in the x-direction.

Figure 17.

Real experiment setup of external disturbance during biped walking.

Figure 18.

Real experiment setup of external disturbance during biped walking

Conclusions and discussion

For the control problem of humanoid walking, many recent researchers focus on generating reference CoM trajectories according to the constraints, predefined ZMP or other dynamics. To resist with small disturbances, online adaptation of those trajectories are also added in some recent approaches, with heuristics that the robot can exactly track these trajectories. However, due to the nature of the robot’s contact with the environment, the ability of the robot to track these trajectories is unavoidably limited, which consequently leads to the out of control and correspondingly the fall of the robot. Although this problem has been noticed in some recent propositions and some approaches are also presented, an applicable and efficient method is yet to be explored. In this article, a novel solution is proposed to the problem raised above. By the online model parameters optimization strategy, the robot can quickly and in real time adjust the body posture and walking patterns under unknown disturbances, hence, the walking robustness can be enhanced. The robot can successfully retrieve its balance even after being pushed by external forces significantly. The gains of the model parameters corrector are learnt through the policy gradient reinforcement learning method, and it shows satisfactory results when applied on both the physical simulated and real full-body humanoid robots.