Flexible model predictive control based on multivariable online adjustment mechanism for robust gait generation

Abstract

The gait generation algorithm considering both step distance adjustment and step duration adjustment could improve the anti-disturbance ability of the humanoid robot, which is very important to the dynamic balance, but the step duration adjustment often brings non-convex optimization problems. In order to avoid this situation and improve the robustness of the gait generator, a gait generation mechanism based on flexible model predictive control is proposed in this article. Specifically, the step distance adjustment and step duration adjustment are set to be optimization objectives, while the change of pressure center is treated as the optimal input to minimize those objectives. With the current system state being used for online re-optimization, a feedback gait generator is formed to realize the strong stability of variable speed and variable step distance walking of the robot. The main contributions of this work are twofold. First, a gait generation mechanism based on flexible model predictive control is proposed, which avoids the problem of nonlinear optimization. Second, a variety of feasible optimization constraints were considered, they can be used on platforms with different computing resources. Simulations are conducted to verify the effectiveness of the proposed mechanism. Results show that as compared with those considering step adjustment only, the proposed method largely improves the compensation ability of disturbance and shortens the adjustment time.

Introduction

Leg movement affects the dynamic balance of humanoid robots, and it is the basic guarantee for robots to complete various advanced tasks. However, owing to the complexities of its hybrid dynamics, the unidirectional constraints on contact forces, as well as the high dimensionality and nonlinearity of robot general dynamics, the robot leg movement is widely regarded to be a difficult problem ¹ and various methods have been proposed. ^{1 –7} Specifically, to address the dynamic balancing problem of humanoid robots, a general method ^8,9 is to use a hierarchical structure shown in Figure 1 to control the robot. As illustrated, the bottom layer is a motion control layer, and it is responsible for tracking the planned robot trajectory by utilizing inverse kinematics or inverse dynamics. While the top level is the motion planning layer, which is utilized to generate a set of desired targets, for example, the Center of Mass (CoM), the Center of Pressure (CoP), the Center of Gravity (CoG), ¹⁰ limb movements as well as some other tracks. The top layer determines the walking gait of the robot and is a prerequisite for stable work. In the case of strong disturbance, it has to adjust the step of the robot quickly to the right position within a period of time like the human beings.

Figure 1.

Block diagram of the variable step duration method. At the beginning of each step, according to the actual state and expected speed, the Hessian matrix and the constraints in M-step prediction horizon are determined. Then, through a QP, the landing point of the swing foot and the next step duration are adjusted. The key point trajectories of CoM and CoP are produced. The swing foot trajectories are interpolated according to the landing point (not covered in this article). Only the first step in the generated M-walking trajectory is used for subsequent control. QP: quadratic programming; CoM: Center of Mass; CoP: Center of Pressure.

For top-level motion planning, many gait generation schemes perform zero moment point (ZMP, equivalent to CoP ² ) conditions by calculating the trajectory suitable for robot’s CoM. Most of them use simplified CoM dynamics models, such as the famous linear inverted pendulum (LIP) model ³ and Cart-Table (CT) model. ⁴ Using a simplified model to capture task-related dynamic processes to a set of linear equations is very useful for real-time generation of walk patterns. Specifically, Kajita et al. ⁴ put forward the concept of preview control and solved a linear quadratic (LQ) regulation problem for a dynamical extension of the CT. Wieber ⁵ further developed it into model predictive control (MPC) and selected the CoM jerk as control signal to achieve asymptotic tracking of the desired ZMP trajectory. However, both methods are not strong enough to disturbance. Recently, Englsberger et al. ⁶ and Takenaka et al. ⁷ split the LIP model into a stable first-order system and an unstable first-order system, and then introduced the intermediate variables of the divergence component of motion (DCM or capture point ¹¹ ) to solve the CoP and CoM trajectories online. Within these methods, the input of the unstable subsystem in LIP is planned online according to the predetermined footstep locations, while both the CoP and CoM trajectories were generated in real time by tracking the DCM according to the robot natural dynamics. Englsberger et al. ¹² extended the method to 3-D case. Yet, it is worth noting that although these methods have certain resistance to perturbations, their effects are still limited.

To enhance the robustness of the generated gait, Herdt et al. ¹³ and Diedam et al. ² used predictive control to include step distance to the optimization objective function, while made footsteps variable to generate both CoM and CoP trajectories under the CoP constraints. Results showed that the generated walking mode is robust to disturbances and achieves high level target tasks, such as desired step position or walking speed. However, the step distance adjustment contributes only to disturbance energy consumption, and they do not take into account the step timing to change the landing speed either. In practice, step duration adjustment could help adjust the robot foot movement speed for more robust gait generation. Since the CoM trajectory is a nonlinear function of time and of initial conditions, in order to deal with nonlinear optimization problems, many algorithms have been proposed in the literature. ^{14 –17} Kryczka et al. ¹⁸ and Aftab et al. ¹⁹ proposed to utilize a nonlinear optimization technique to modify both footstep positions and step-timing to maintain dynamic stability during the robot walking process. Specifically, Aftab et al. introduced a simple model of the mechanical cost in the objective function by penalizing the acceleration of the swing foot, while the acceleration can’t anyway exceed a given maximum value. Kryczka et al. rewrote the optimization objective to be a nonlinear function of the optimization variables. Although satisfactory results have been achieved, the nonlinear optimization process introduced high computational costs and cannot guarantee convergence to the minimum either.

Khadiv et al. ²⁰ adopted a DCM method by taking the constraints of both location and time of stepping into account in robot gait generation process and modeled the problem as a quadratic program that can be solved real time. The results showed that such the proposed strategy could help the robot recover from severe pushes. However, since the CoM motion constraint has not been taken into account in the objective function, the CoM speed tracking error cannot be minimized in the control process. Sun et al. ²¹ proposed a class of global and feasible projected Fletcher–Reeves conjugate gradient approach. This method guarantees the tracking accuracy of each physical quantity to the target value in the optimization process but requires the optimization constraints to be linear. Sun et al. ²² further proposed a superlinearly convergent trust region-sequential quadratic programming (QP) approach. The method incorporates a combination algorithm that allows both the trust region technique ^23,24 and the sequential QP method to be used. It avoids solving the QP subproblem for nonlinear constrained optimization problems, which gives the potential for fast convergence in the neighborhood of an optimal solution.

In this article, we propose to flexibly choose the analytical bounds of position, velocity, and acceleration of CoM in MPC frame to predict the stability and manually set the target footstep locations and step duration in the prediction horizon. The optimal outputs are allowed to deviate from the target, that is, to adjust the item weight coefficient in the optimization objective, and the variation of CoP are selected as the input to minimize the objective. With the current system state being used to recalculate the optimization online, a feedback controller as shown in Figure 1 is formed, which outputs the required CoM, CoP, footstep locations, and the next step timing. Finally, according to the practical robot structure, the maximum and minimum stride distances of the lateral plane and the sagittal plane are also included as the optimization constraint of QP to achieve variable speed and variable step.

The main contributions of this study could be summarized as follows. First, we improved the MPC objective function with adjustment of step duration taken into account, and therefore, such a step timing adjustment strategy could largely suppress the interference in forward visual field and stabilize CoM velocity mutation with the non-convex optimization problem avoided, as compared with those in the existing work. ^2,4,25 Second, the proposed method flexibly handles the stability optimization problem for biped robots and also presents a variety of gait optimization constraints, which can be easily implemented on control platforms with different computing resources. With the foot CoP stabilized without any mutations in the forward motion, the proposed mechanism mimics human walking, and thus is conducive to walking smoothly with small errors.

The article is organized as follows. In the “The MPC method for walking machine” section, we recall the LIP and the CoM trajectory characterization of standard MPC. The “System optimization model” section describes the motion model adopted, the “Walking planning with adjustable step duration” section achieves flexible walking optimization through carefully selected objective function and constraints, and the “Push recovery planning” section describes push recovery. To illustrate its benefits, all the simulation results are shown in “Results and discussion” section, the conclusion points out the next research work. For simplicity purpose, the abbreviations utilized in this article are summarized in Table 1 as below.

Table 1.

The full name of acronyms.

The acronyms	The full name
MPC	Model predictive control
CoM	Center of Mass
CoP	Center of Pressure
CoG	Center of Gravity
DCM	Divergence component of motion
LIP	Linear inverted pendulum
LQ	Linear quadratic
QP	Quadratic programming
ZMP	Zero moment point
CT	Cart-Table

The MPC method for walking machine

When a robot is walking on the ground, it comes down to the fact that the CoP can lie only within the convex hull of the contact points between the robot’s feet and the ground. For the LIP model used by most scholars, it describes the motion of the CoM of a robot when its height is constant while the rotation effect is not considered. In addition to its linear properties, the LIP has the advantage that dynamics along the x- and y-directions are decoupled and can be represented by the same differential equations. Many studies ^{4,26 –30} have proved that, despite its simplicity, the use of LIP (see Figure 2) for gait generation is effective. Simplified LIP model for real robots is shown in Figure 3. Establish LIP dynamic equation in the x-axis (the y-axis identical)

Figure 2.

The LIP model of biped. LIP: linear inverted pendulum.

Figure 3.

The biped robot.

$\ddot{x}={\omega }^2(x-x_{\text{cop}})$

1

where $\omega =\sqrt{g/z_c}$ , with g and z_c being the acceleration of gravity and the constant altitude of CoM, respectively. x is the horizontal position of the CoM, $\ddot{x}$ is horizontal acceleration of the CoM, and $x_{\text{cop}}$ is the position of the CoP.

It could be found in equation (1) that LIP has a right half plane pole, and there is always a divergence component in the system, which leads to the divergence of CoM. In addition, the change of CoP will directly lead to the change of the divergence velocity of CoM. However, CoP, the input signal has some constraints and must be within the range of stable polygons. The scheme proposed by Kajita et al. ⁴ generates a trajectory of the CoM under the constraint that the footsteps are fixed and impossible to change. CoM and CoP are programmed as a series of discrete points, with varied jerks $\stackrel{⃛}{x}$ over time intervals of constant lengths T. Therefore, the dynamics at $t_{k+1}=(k+1)T$ could be calculated as

$\begin{array}{l}{x}_{k+1} = A{x}_k+B\stackrel{⃛}{\mathit{x}}\\ y_{k+1} = C{x}_k\end{array}$

2

where $x=[x,\dot{x},\ddot{x}{]}^T$ , $y=x_{\text{cop}}$ , and

$\begin{array}{l}\begin{array}{l}A=\left[\begin{array}{ccc}1& T& T^2/2\\ 0& 1& \mathrm{T}\\ 0& 0& 1\end{array}\right],B=\left[\begin{array}{c}{T}^3/6\\ T^2/2\\ T\end{array}\right]\hfill \end{array}\\ \begin{array}{l}C=\left[\begin{array}{ccc}1& 0& -1/{\omega }^2\end{array}\right]\hfill \end{array}\end{array}$

This scheme is to minimize this jerk while maintaining a position of the CoP as close as possible to some prescribed reference positions, but in the process, there is no consideration of the robot’s foothold and step timing. Therefore, to address the problems of anti-disturbance and variable speed walking for the robot, we propose to add the velocity deviation and acceleration deviation of CoM, as well as landing point deviation and stepping time deviation to the optimization function, and select the variation of CoP as the optimal input. After the first QP solution is finished, the feedback value of the system state is recalculated online to form a feedback controller, and the optimal trajectory is constantly updated.

Gait generation strategy

System optimization model

For all of the implementations presented, the model is the LIP with 2-D dynamics given by equation (1). However, the math is generally independent of the model. It is straightforward to implement any other linear model. As shown in Figure 4, a complete step includes double support phase $t_{\text{double}}$ and single support phase $t_{\text{single}}$ , and NT time span may contain M steps. With the state space equation of the discrete linear system shown in equation (3) taken into account, the control input u is the variation of $x_{\text{cop}}$ . Given a sequence of control inputs $\overline{\mathit{u}}=[0,u_k,u_{k+1}\mathrm{,...},u_{k+N-1}{]}^T$ , the $X=[x_k{}^T,x_{k+1}{}^T,x_{k+2}{}^T,\mathrm{...},x_{k+N}{}^T{]}^T$ and $Y=[y_{k+1}{}^T,{\mathbf{y}}_{\mathit{k}+2}{}^T,\mathrm{...},{\mathit{y}}_{\mathit{k}+\mathit{N}+1}{}^{\mathit{T}}{]}^{\mathit{T}}$ sequences within NT time intervals can be calculated

Figure 4.

Predicted horizon of M-steps.

$\begin{array}{l}{x}_{k+1}=A{x}_k+B{u}_k\\ y_{k+1}=C{x}_k+D{u}_k\end{array}$

3

$\begin{array}{l}\begin{array}{l}A=\left[\begin{array}{c}\begin{array}{ccccc}1+\frac{1}{2}{T}^2{\omega }^2& T& -\frac{1}{2}{T}^2{\omega }^2& 0& 0\\ T{\omega }^2& 1& T{\omega }^2& 0& 0\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\end{array}\right],\hspace{0.17em}B=\left[\begin{array}{c}\begin{array}{c}-\frac{1}{2}{T}^2{\omega }^2\\ -T/{\omega }^2\\ 1\\ 0\\ 0\end{array}\end{array}\right]\hfill \end{array}\\ \begin{array}{l}C=\left[\begin{array}{c}\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ {\omega }^2& 0& -{\omega }^2& 0& 0\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\end{array}\right],\hspace{0.17em}D=0\hfill \end{array}\end{array}$

where $x=\hspace{0.17em}{[x,\dot{x},x_{\text{cop}},x_{\text{fl}},x_{\text{fr}}]}^T$ , $y=\hspace{0.17em}{[x,\dot{x},\ddot{x},x_{\text{cop}},x_{\text{fl}},x_{\text{fr}}]}^T$ , $x_{\text{fl}},x_{\text{fr}}$ are the left and right foot position

$\begin{array}{l}X=\hspace{0.17em}\overline{\mathit{A}}{x}_k+\overline{\mathit{B}}\overline{\mathit{u}}\\ Y\hspace{0.17em}=\hspace{0.17em}\overline{\mathit{C}}{x}_k+\overline{\mathit{D}}\overline{\mathit{u}}\end{array}$

4

$\begin{array}{l}\begin{array}{l}\overline{A}=\left[\begin{array}{c}\begin{array}{c}{I}_{5*5}\\ A\\ A^2\\ A^3\\ ⋮\\ A^N\end{array}\end{array}\right],\overline{B}=\left[\begin{array}{c}\begin{array}{ccccc}0& 0& 0& \cdots & 0\\ B& 0& 0& 0& 0\\ AB& B& 0& 0& ⋮\\ ⋮& ⋮& \ddots & ⋮& 0\\ A^{N-2}B& \cdots & AB& B& 0\\ A^{N-1}B& A^{N-2}B& \cdots & AB& B\end{array}\end{array}\right]\hfill \end{array}\\ \begin{array}{l}\overline{C}=\left[\begin{array}{c}\begin{array}{c}C\\ CA\\ C{A}^2\\ C{A}^3\\ ⋮\\ C{A}^{N-1}\end{array}\end{array}\right],\overline{D}=\left[\begin{array}{c}\begin{array}{ccccc}0& 0& 0& \cdots & 0\\ CB& 0& 0& 0& 0\\ CAB& CB& 0& 0& 0\\ ⋮& ⋮& \ddots & ⋮& 0\\ C{A}^{N-2}B& \cdots & CAB& CB& 0\\ C{A}^{N-1}B& C{A}^{N-2}B& \cdots & CAB& CB\end{array}\end{array}\right]\hfill \end{array}\end{array}$

In order to determine the best trajectory, a cost function is constructed. The function is a scalar function that scores the trajectory, by defining a quadratic cost on the elements of X and $\overline{\mathit{u}}$ such that

$\begin{array}{l}J=\frac{1}{2}{X}^{T}QX+\frac{1}{2}{\overline{\mathit{u}}}^{T}R\overline{\mathit{u}}\hfill \end{array}$

5

Introducing it into equation (4), then

$\begin{array}{l}J=\frac{1}{2}{\overline{\mathit{u}}}^T(R+{\overline{\mathit{B}}}^{T}Q\overline{\mathit{B}})\overline{\mathit{u}}+x_k^T{\overline{\mathit{A}}}^{T}Q\overline{\mathit{B}}\overline{\mathit{u}}+\frac{1}{2}{x}_k^T{\overline{\mathit{A}}}^{T}Q\overline{\mathit{A}}{x}_k\end{array}$

6

It has a simple quadratic form

$\begin{array}{l}J=\frac{1}{2}{\overline{\mathit{u}}}^{T}H\overline{\mathit{u}}+f\overline{\mathit{u}}+J_0\hfill \end{array}$

7

where

$\begin{array}{l}H=R+{\overline{\mathit{B}}}^{T}Q\overline{\mathit{B}},\\ f=x_k^T{\overline{\mathit{A}}}^{T}Q\overline{\mathit{B}},\\ J_0=\frac{1}{2}{x}_k^T{\overline{\mathit{A}}}^{T}Q\overline{\mathit{A}}{x}_k\end{array}$

Note that the J ₀ term is constant with respect to $\overline{u}$ and has no effect on the location of the minimum of J.

Adding equality and inequality dynamic equilibrium constraints

$\begin{array}{l}{c}_{eq}\overline{\mathit{u}}=b_{eq}\hfill \end{array}$

8

$\begin{array}{l}{c}_{in}\overline{\mathit{u}}\le b_{in}\hfill \end{array}$

9

For some constrained states, ${c\prime }_{in}X\le {b\prime }_{in}$ can be transformed into

$\begin{array}{l}{c\prime }_{in}\overline{\mathit{B}}\overline{\mathit{u}}\le {b\prime }_{in}-{c\prime }_{in}\overline{\mathit{A}}{x}_k\hfill \end{array}$

10

Finally, the problem of finding CoM trajectory in a period of time NT is converted into a QP. There are many methods can be used to solve this standard QP, such as interior points, active sets, conjugant gradient, or simplex methods ³¹

$\begin{array}{l}\underset{\overline{u}}{\text{min}}J=\frac{1}{2}{\overline{\mathit{u}}}^{T}H\overline{\mathit{u}}+f\overline{\mathit{u}},\\ \text{st}\mathrm{.}\left\{\begin{array}{c}{c}_{eq}\overline{\mathit{u}}=b_{eq}\\ \left[\begin{array}{c}{c}_{in}\\ {c\prime }_{in}\overline{\mathit{B}}\end{array}\right]\overline{\mathit{u}}\le \left[\begin{array}{c}{b}_{in}\\ {b\prime }_{in}-{c\prime }_{in}\overline{\mathit{A}}{x}_k\end{array}\right]\end{array}\right.\end{array}$

11

Walking planning with adjustable step duration

The MPC scheme adopted by Kajita et al. ⁴ is basic in generating stable CoM trajectories. Its only mandatory feature is to adjust the magnitude of the motion derivative of CoM in the prediction horizon. However, any control variables that help to suppress CoM mutations contribute to more stable walking.

Objective function

For the dynamic balance of biped robot, the form of function will be different under different circumstances. In case of steady walking, this article takes the following form of objective function

$\begin{array}{ll}J=\hfill & \frac{a_1}{2}\parallel C_x^{\text{target}}-C_x{\parallel }^2+\frac{a_2}{2}\parallel {\dot{C}}_x^{\text{target}}-{\dot{C}}_x{\parallel }^2\hfill \\ \hfill & +\frac{a_3}{2}\parallel {\ddot{C}}_x^{\text{target}}-{\ddot{C}}_x{\parallel }^2+\frac{a_4}{2}\parallel {\overline{u}}_x{\parallel }^2\hfill \\ \hfill & +\frac{a_5}{2}\parallel x_f^{\text{target}}-x_f{\parallel }^2+\frac{a_6}{2}\parallel T_{x\text{step}}^{\text{target}}-T_{x\text{step}}{\parallel }^2\hfill \\ \hfill & +\frac{a_1}{2}\parallel C_y^{\text{target}}-C_y{\parallel }^2+\frac{a_2}{2}\parallel {\dot{C}}_y^{\text{target}}-{\dot{C}}_y{\parallel }^2\hfill \\ \hfill & +\frac{a_3}{2}\parallel {\ddot{C}}_y^{\text{target}}-{\ddot{C}}_y{\parallel }^2+\frac{a_4}{2}\parallel {\overline{u}}_y{\parallel }^2\hfill \\ \hfill & +\frac{a_5}{2}\parallel y_f^{\text{target}}-y_f{\parallel }^2+\frac{a_6}{2}\parallel T_{y\text{step}}^{\text{target}}-T_{y\text{step}}{\parallel }^2\hfill \end{array}$

12

where $(C_x,C_y)$ , $({\dot{C}}_x,{\dot{C}}_y)$ , and $({\ddot{C}}_x,{\ddot{C}}_y)$ are the vectors of CoM position, velocity, and acceleration in 2-D over the next N time steps, respectively. $(C_x^{\text{target}},C_y^{\text{target}})$ , $({\dot{C}}_x^{\text{target}},{\dot{C}}_y^{\text{target}})$ , and $({\ddot{C}}_x^{\text{target}},{\ddot{C}}_x^{\text{target}})$ are the target position, target velocity, and target acceleration of CoM, respectively. $({\overline{\mathit{u}}}_x,{\overline{\mathit{u}}}_y)$ in this case is the variation of the $x_{\text{cop}}$ in 2-D. $x_f^{\text{target}},y_f^{\text{target}},x_f$ , and y_f are vectors of the next M reference and actual foothold, respectively. $T_{x\text{step}}^{\text{target}},T_{y\text{step}}^{\text{target}}$ are the target step timing in the x-direction and in the y-direction, respectively. $T_{x\text{step}}$ and $T_{y\text{step}}$ are the optimal outputs in their respective directions, $a_1,a_2,a_3,a_4,a_5$ , and a ₆ are the weight coefficients of corresponding variables in J.

As shown in Figure 5, in the calculation at a time, the robot can walk one step or walk $M>1$ steps. So here, we chose to give it three steps, $M=3$ , where $x_f=[x_{f_1},x_{f_2},x_{f_3}{]}^T$ , $y_f=[y_{f_1},y_{f_2},y_{f_3}{]}^T$ . The reference value is selected as

Figure 5.

M-steps quadratic programming process.

$\{\begin{array}{l}{C}_x^{\text{target}}=\frac{1}{2}(x_{s_0}+x_{f_0})+v_{x\text{ref}}NT\hfill \\ C_y^{\text{target}}=\frac{1}{2}(y_{s_0}+y_{f_0})+v_{y\text{ref}}NT\hfill \end{array},$

$\begin{array}{l}\{\begin{array}{l}{\dot{C}}_x^{\text{target}}=v_{x\text{ref}}\hfill \\ {\dot{C}}_y^{\text{target}}=v_{y\text{ref}}\hfill \end{array},\{\begin{array}{l}{\ddot{C}}_x^{\text{target}}=0\hfill \\ {\ddot{C}}_y^{\text{target}}=0\hfill \end{array},\{\begin{array}{l}{T}_{x\text{step}}^{\text{target}}=t_{\text{step}}\hfill \\ T_{y\text{step}}^{\text{target}}=t_{\text{step}}\hfill \end{array}\end{array}$

$\begin{array}{l}\{\begin{array}{l}{x}_f^{\text{target}}=[x_{f_0}+\frac{1}{2}{v}_{x\text{ref}},x_{f_1}+\frac{1}{2}{v}_{x\text{ref}},x_{f_2}+\frac{1}{2}{v}_{x\text{ref}}{]}^T\hfill \\ y_f^{\text{target}}=[y_{f_0}+d,y_{f_1}-d,y_{f_2}+d]\hfill \end{array}\end{array}$

$d=0.2\text{sign}(y_{s_0}-y_{f_0})$ , $(x_{s_0},y_{s_0})$ is the position of initial swing foot, $(x_{f_0},y_{f_0})$ is the position of initial support foot, $v_{x\text{ref}}$ and $v_{y\text{ref}}$ are the reference velocity of CoM, $t_{\text{step}}$ is a known deterministic value. Note: The $(x_f,y_f)$ , $T_{x\text{step}}$ , and $T_{y\text{step}}$ are unknown in advance, the objective function J can be written as a function of u

$u=\hspace{0.17em}{[{\overline{\mathit{u}}}_x^T,x_f^T,T_{x\text{step}},{\overline{\mathit{u}}}_y^T,y_f^T,T_{y\text{step}}]}^T\in R^{2(N+M+1)}$

13

After solving the QP, the value of u is obtained. The CoM and CoP trajectories can be got from equation (4) and initial status x ₀. Take the minimum values in $T_{x\text{step}}$ and $T_{y\text{step}}$ for the next QP, that is, $T=\text{min}(T_{x\text{step}},T_{y\text{step}})$ . Although we worked out the CoM and CoP trajectories in the M-walk process, only the first foothold is taken as the current destination, the QP is recomputed in next NT time.

Constraints

Constraints play a major role in the computation of the optimal trajectory. In the robot walking process, CoP directly determines the balance state of the robot. It should also be noted that since the locations of the footsteps are variables in optimization process, the constraints on the CoP after the first step are unknown, and thus, it is necessary to dynamically set CoP constraints according to the optimized planning footstep locations. In addition, in order to maintain a simple linear form and use a fast QP solver, simple conservative constraints are chosen to approximate the real constraints. The conservative constraints represent smaller regions than the real constraints. Another advantage of this approach is that it introduces a safety margin and enhances the stability of optimization.

• CoP constraints

Since robot is in a different supporting phase at different times, the CoP constraints could be divided into a set $Z_{\text{cop}}^S$ consisting of multiple single support constraint blocks and a set $Z_{\text{cop}}^D$ consisting of multiple double support constraint blocks throughout the NT time (see Figure 5).

Single support phase (during $t_{\text{single}}$ )

When in single support, the CoP can only be within the range of the support foot, so the feasible range of the CoP is

$\begin{array}{l}{\displaystyle \sum _{i=0}^2}(x_{f_i-{\Delta }_1})U_i^S\le x_{cop}\le {\displaystyle \sum _{i=0}^2}(x_{f_i+{\Delta }_2})U_i^S\\ {\displaystyle \sum _{i=0}^2}(y_{f_i-{\Delta }_3})U_i^S\le y_{cop}\le {\displaystyle \sum _{i=0}^2}(y_{f_i+{\Delta }_4})U_i^S\end{array}$

14

where $U_i^S\in R^{N\times 1}$ be defined as a vector of zeros and ones with the ones corresponding to the time steps of the ith single phase, with the first phase being the initial double support phase. ${\Delta }_1,{\Delta }_2,{\Delta }_3$ , and ${\Delta }_4$ are the stable ranges, which must be less than half of the foot size. From equation (4), we get

$\begin{array}{l}{x}_{cop}={\overline{\mathit{A}}}_{x_{cop}}{x}_k+{\overline{\mathit{B}}}_{x_{cop}}\overline{\mathit{u}}\hfill \\ y_{cop}={\overline{\mathit{A}}}_{y_{cop}}{x}_k+{\overline{\mathit{B}}}_{y_{cop}}\overline{\mathit{u}}\hfill \end{array}$

15

${\overline{\mathit{A}}}_{x_{cop}},{\overline{\mathit{B}}}_{x_{cop}}$ and ${\overline{\mathit{A}}}_{y_{cop}},{\overline{\mathit{B}}}_{y_{cop}}$ are the corresponding line in $\overline{\mathit{A}},\overline{\mathit{B}}$ . Finally, we have

$\begin{array}{l}{\displaystyle \sum _{i=0}^2}(x_{f_i-{\Delta }_1})U_i^S\le {\overline{\mathit{A}}}_{x_{cop}}{x}_k+{\overline{\mathit{B}}}_{x_{cop}}\overline{u}\le {\displaystyle \sum _{i=0}^2}(x_{f_i+{\Delta }_2})U_i^S\hfill \\ {\displaystyle \sum _{i=0}^2}(y_{f_i-{\Delta }_3})U_i^S\le {\overline{\mathit{A}}}_{y_{cop}}{x}_k+{\overline{\mathit{B}}}_{y_{cop}}\overline{\mathit{u}}\le {\displaystyle \sum _{i=0}^2}(y_{f_i+{\Delta }_4})U_i^S\hfill \end{array}$

16

Double support phase (during $t_{\text{double}}$ )

When in double support, establish the coordinate system $o^{\prime }$ as shown in Figure 6(a) in the middle of the feet. $o^{\prime }$ rotates θ around the z-axis relative to o. Hence, for any CoP $(x_{\text{cop}},y_{\text{cop}})$ in planning trajectory within coordinate $o^{\prime }$ , it is $({x\prime }_{\text{cop}},{y\prime }_{\text{cop}})$ , and the constraints of double support are as follows

Figure 6.

CoP constraints in the double support. (a) The first form, which is closest to human beings, (b) the second form, which is conservative and simplified, and (c) the third form, which is similar to fast walking. CoP: Center of Pressure.

$\begin{array}{l}\Vert \overleftarrow{{x\prime }_{\text{cop}}}\Vert \le {\Delta \prime }_x\hfill \\ \Vert \overleftarrow{{y\prime }_{\text{cop}}}\Vert \le {\Delta \prime }_y\hfill \end{array}$

17

where ${\Delta \prime }_x,{\Delta \prime }_y$ are the stability margin in the direction of $x^{\prime }$ and $y^{\prime }$ as illustrated in Figure 6(a), and

$\begin{array}{ll}{x\prime }_{\text{cop}}\hspace{0.17em}=\hspace{0.17em}(x_{\text{cop}}\hspace{0.17em}-\hspace{0.17em}{x}_r)\text{cos}\theta \hspace{0.17em}-\hspace{0.17em}(y_{\text{cop}}\hspace{0.17em}-\hspace{0.17em}{y}_r)\text{sin}\theta \hfill & \hfill \\ {y\prime }_{\text{cop}}\hspace{0.17em}=\hspace{0.17em}(y_{\text{cop}}\hspace{0.17em}-\hspace{0.17em}{y}_r)\text{cos}\theta \hspace{0.17em}-\hspace{0.17em}(x_{\text{cop}}\hspace{0.17em}-\hspace{0.17em}{x}_r)\text{sin}\theta \hspace{0.17em}-\hspace{0.17em}\frac{l_{\text{pace}}}{2}\hfill & \hfill \end{array}$

18

Because $l_{\text{pace}}=\sqrt{{(x_l-x_r)}^2+{(y_l-y_r)}^2}$ , let $l^{\prime }=1/l_{\text{pace}}$ , the equation (18) can be converted to

$\begin{array}{ll}{x\prime }_{\text{cop}}\hspace{0.17em}=\hspace{0.17em}{l}^{\prime }[(x_{\text{cop}}\hspace{0.17em}-\hspace{0.17em}{x}_r)(y_l\hspace{0.17em}-\hspace{0.17em}{y}_r)\hspace{0.17em}-\hspace{0.17em}(y_{\text{cop}}\hspace{0.17em}-\hspace{0.17em}{y}_r)(x_l\hspace{0.17em}-\hspace{0.17em}{x}_r)]\hfill & \hfill \\ {y\prime }_{\text{cop}}\hspace{0.17em}=\hspace{0.17em}{l}^{\prime }[(y_{\text{cop}}\hspace{0.17em}-\hspace{0.17em}{y}_r)(y_l\hspace{0.17em}-\hspace{0.17em}{y}_r)\hspace{0.17em}-\hspace{0.17em}(x_{\text{cop}}\hspace{0.17em}-\hspace{0.17em}{x}_r)(x_l\hspace{0.17em}-\hspace{0.17em}{x}_r)]\hspace{0.17em}-\hspace{0.17em}\frac{1}{2l^{\prime }}\hfill & \hfill \end{array}$

19

Then combine equation (15) and (19), equation (17) are transformed into

$\begin{array}{l}{\displaystyle \sum _{i=0}^3}-l_{\text{pace}}{\Delta \prime }_x{U}_i^D\le [({\overline{\mathit{A}}}_{x_{cop}}{x}_k+{\overline{\mathit{B}}}_{x_{cop}}\overline{\mathit{u}}-x_r{U}_i{}^D)(y_l-y_r)\hfill \\ -({\overline{\mathit{A}}}_{y_{cop}}{x}_k+{\overline{\mathit{B}}}_{y_{cop}}\overline{\mathit{u}}-y_r{U}_i{}^D)(x_l-x_r)]\le {\displaystyle \sum _{i=0}^3}{l}_{\text{pace}}{{\Delta }^{\prime }}_x{U}_i^D\hfill \\ {\displaystyle \sum _{i=0}^3}-l_{\text{pace}}{\Delta \prime }_y{U}_i^D\le [({\overline{\mathit{A}}}_{y_{cop}}{x}_k+{\overline{\mathit{B}}}_{y_{cop}}\overline{\mathit{u}}-y_r{U}_i{}^D)(y_l-y_r)\hfill \\ +({\overline{\mathit{A}}}_{x_{cop}}{x}_k+{\overline{\mathit{B}}}_{x_{cop}}\overline{u}-x_r{U}_i{}^D)(x_l-x_r)]-\frac{l_{\text{pace}}^2}{2}{U}_i^D\hfill \\ \le {\displaystyle \sum _{i=0}^3}{l}_{\text{pace}}{{\Delta }^{\prime }}_y{U}_i^D\hfill \end{array}$

20

where $U_i^D\in R^{N\times 1}$ be defined as a vector of zeros and ones with the ones corresponding to the time steps of the ith double phase, Hence

$\begin{array}{l}{\displaystyle \sum _{i=0}^2}(U_i^S+U_i^D)+U_3^D=I_{N\times 1}\hfill \end{array}$

21

$I_{N\times 1}$ is a N-dimensional vector with all elements being 1.

The constraints of CoP in double support can also be selected in other forms. The optimal results obtained in different forms are slightly different in the macroscopical attitude of the robot. For example (as shown in Figure 6(b))

$\begin{array}{ll}-{\Delta }_x\le \frac{x_l+x_r}{2}-x_{\text{cop}}\le {\Delta }_x\hfill & \hfill \\ -{\Delta }_y\le \frac{y_l+y_r}{2}-y_{\text{cop}}\le {\Delta }_y\hfill & \hfill \end{array}$

22

Alternatively, the CoP constraints of the double support and the single support can be directly set up to coincide with each other as shown in Figure 6(c).

• Step distance constraints

Because of the limitation of mechanical structure, the step distance needs to be constrained

$\begin{array}{ll}-{\Delta }_{dx}\le x_{f_i}-x_{f_{i-1}}\le {\Delta }_{dx}\hfill & \hfill \\ -{\Delta }_{dy}\le y_{f_i}-y_{f_{i-1}}\le {\Delta }_{dy}\hfill & \hfill \end{array}$

23

where ${\Delta }_{dx},{\Delta }_{dy}$ are the maximum step distances of the robot.

• Foot placement constraints

When walking steadily, the landing point swings back and forth in the direction of the y-axis, yet too small swing amplitude may cause interference between the legs. Hence, in the y-axis, we set the distance between the two landing points is greater than a threshold.

$\begin{array}{ll}\hfill & y_{f_i}-y_{f_{i-1}}\ge {\Delta }_{\text{min}y}\hspace{0.17em}i=1,2,3\hfill \\ \text{or }\hfill & y_{f_i}-y_{f_{i-1}}\ge {\Delta }_{\text{min}y}\hspace{0.17em}i=1,2,3\hfill \end{array}$

24

where ${\Delta }_{\text{min}y}$ is the minimum step distance of the robot in the y-axis. In the next few steps, the landing point $y_{f_i}$ should be determined according to the size of initial conditions $y_{f_0}$ and $y_{s_0}$ . No crossing of feet allowed.

• Step duration constraints

In the case of large disturbance, humanoid robot is limited by mechanical structure. The maximum distance of step is equal to ${\Delta }_{dx},{\Delta }_{dy}$ under the condition of fixed step timing. If it is not enough to overcome disturbance, it will lead to body forward and CoM divergence. If step duration is changeable, the robot can choose the maximum step distance to land in a shorter step time to quickly enter the next step. Through several rapid steps, it has more advantages in dealing with high-speed and large disturbances, which is similar to human beings. So when the initial velocity is obtained, the distance of one step should not exceed ${\Delta }_{dx},{\Delta }_{dy}$

$\begin{array}{l}0\le x_0(2)T_{x\text{step}}\le {\Delta }_{dx}\hfill \\ 0\le y_0(2)T_{y\text{step}}\le {\Delta }_{dy}\hfill \end{array}$

25

The above method finds the real-time optimal trajectory by rolling the QP within the constrained framework, in NT time and minimizes CoM position error, velocity error, acceleration error, landing point error, and step timing error. In addition, the variation of CoP could also be minimized. Therefore, the robot can change its speed and walk with variable step size, and automatically adjust the step size and time according to the required speed or disturbance. As weight coefficients of the objective function, the parameters $[a_1,a_2,a_3,a_4,a_5,a_6]$ determines the constraint strength of the corresponding physical quantity. The larger the coefficient, the smaller the error between the trajectory output and the target. While if the robot is disturbed, different initial state x ₀ appears in each QP solution. Because of the online operation characteristics of the algorithm, the optimal output results are also changed to generate the optimal target trajectory under the disturbance.

Push recovery planning

In the same way, if the cost function and all kinds of constraints can be found manually, the rolling optimization method can be used to carry out the real-time gait planning. For push recovery, the large disturbance would grant the robot an initial velocity, and thus, the robot needs to be stabilized by step, which can be stabilized at one step or several times.

Adopt the following form of objective function

$\begin{array}{ll}J=\hfill & \frac{a_1}{2}\parallel C_x^{\text{target}}-C_x{\parallel }^2+\frac{a_2}{2}\parallel {\dot{C}}_x{\parallel }^2+\frac{a_3}{2}\parallel {\overline{\mathit{u}}}_x{\parallel }^2\hfill \\ \hfill & +\frac{a_4}{2}\parallel x_f^{\text{target}}-x_f{\parallel }^2\hfill \\ \hfill & +\frac{a_1}{2}\parallel C_y^{\text{target}}-C_y{\parallel }^2+\frac{a_2}{2}\parallel {\dot{C}}_y{\parallel }^2+\frac{a_3}{2}\parallel {\overline{\mathit{u}}}_y{\parallel }^2\hfill \\ \hfill & +\frac{a_4}{2}\parallel y_f^{\text{target}}-y_f{\parallel }^2\hfill \end{array}$

26

The variable is in line with the previous section

$\begin{array}{l}\left\{\begin{array}{l}{C}_x^{\text{target}}=x_{f_1}\hfill \\ C_y^{\text{target}}=y_{f_1}\hfill \end{array}\right.,\left\{\begin{array}{l}{x}_f^{\text{target}}=x_{f_0}\hfill \\ y_f^{\text{target}}=y_{f_0}+0.2\text{sign}(y_{s_0}-y_{f_0)}\hfill \end{array}\right.\hfill \end{array}$

(as shown in Figure 7(a)), or

Figure 7.

Push recovery planning process. (a) One-step stabilization mode and (b) multiple step stabilization mode.

$\begin{array}{l}\left\{\begin{array}{l}{C}_x^{\text{target}}=\frac{x_{f_0}+x_{f_1}}{2}\hfill \\ C_y^{\text{target}}=\frac{y_{f_0}+y_{f_1}}{2}\hfill \end{array}\right.,\left\{\begin{array}{l}{x}_f^{\text{target}}=x_{f_0}\hfill \\ y_f^{\text{target}}=y_{f_0}+0.2\text{sign}(y_{s_0}-y_{f_0)}\hfill \end{array}\right.\hfill \end{array}$

In the case of one-step stability, the robot reaches a stable state through one foot support, and the double support constraints can be selected in the form of Figure 7(a). While for multiple stride, the CoP constraint given by equation (22) in the double support, as shown in Figure 7(b), could be selected, and the optimization method of push recovery is set manually. Other structural constraints are identical to those in the “Walking planning with adjustable step duration” section.

Results and discussions

In this section, we present simulation results using the proposed walking pattern generation method. In the first scenario, the results show that the model can recover from the larger pushes when the step time adaptive controller is used. In the second scenario, we compare the maximize anti-disturbance ability of the step duration adjustable method with that of the fixed step duration method in Diedam et al., ² Kajita et al., ⁴ and Herdt et al. ^13,25

Walking planning simulation

In this scenario, we compare our controller with the one that uses fixed step durations. Walking simulation was carried out with Matlab (version R2017b) on a personal computer with x64 Win10 platform (Intel(R) Core(TM) i5-7200U CPU, 8G RAM). At the beginning, the robot is in a static state in the flat ground. Using LIP model, the height of CoM remains unchanged. The state space equation (3) is used, the structural parameters of the robot (see Figure 3) are shown in Table 2. The total weight of the robot is 25.5 kg and it has 12 degrees of freedom. Each leg has three degrees of freedom (pitch, roll, yam) on the hip joint, one pitch on the knee, the ankle joint has pitch and roll degrees of freedom. Each joint is equipped with an angular displacement and torque sensors to achieve position or force closed-loop. Inertial measurement unit (IMU) units are installed on floating base coordinates.

Table 2.

System parameters.

Parameter	Value
M	3
$T\hspace{0.17em}(\text{s})$	0.01
${\omega }^2=g/h$	9.81
$t_{x\text{step}}\hspace{0.17em}(\text{s})$	0.5
${\Delta }_{dx}={\Delta }_{dy}\hspace{0.17em}(\text{m})$	0.3
$[a_1,\hspace{0.17em}{a}_2,\hspace{0.17em}{a}_3,\hspace{0.17em}{a}_4,\hspace{0.17em}{a}_5,\hspace{0.17em}{a}_6]$	${[10}^{-2},\hspace{0.17em}{10}^{-2},\hspace{0.17em}0.3,\hspace{0.17em}{10}^{-3},\hspace{0.17em}{10}^{-2},\hspace{0.17em}{10}^{-4}]$
$[{\Delta }_1,\hspace{0.17em}{\Delta }_2,\hspace{0.17em}{\Delta }_3,\hspace{0.17em}{\Delta }_4]\hspace{0.17em}(\text{m})$	[0.06, 0.06, 0.04, 0.04]
$[{\Delta }_{x^{\prime }},\hspace{0.17em}{\Delta }_{y^{\prime }}]\hspace{0.17em}(\text{m})$	[0.4, 0.2]
$[{\Delta }_x,\hspace{0.17em}{\Delta }_y]\hspace{0.17em}(\text{m})$	[0.06, 0.04]
${\Delta }_{\text{min}y}\hspace{0.17em}(\text{m})$	0.2

Figures 8 to 11 show the simulation results of fixed step timing and adaptive step timing of the humanoid under disturbance. In the dual-support phase, the CoP constraint selection equation (20). $x_0=[\mathrm{0,0.2,0,0,0,0,0,0,0.10},-0.10]$ , $v_{x\text{ref}}=0.5$ m/s, $v_{y\text{ref}}=0.1$ m/s. At the initial time of 5th step, the x-direction is added a disturbance of 0.3 m/s, the y-direction is added a disturbance of 0.1 m/s. For the case when step timing is adjustable, as shown in Figures 8 and 10, the system shortens the step timing, has a recovery motion that starts with a large recovery step, and then converges to the reference motion in 1 or 2 steps. However, for the case when the step timing cannot be adjusted, as shown in Figures 9 and 11, because of the excessive disturbance, the humanoid robots uses the maximum step distance at each step, but it is still not enough to suppress the CoM divergence, which eventually leads to the robot dumping.

Figure 8.

The CoM and CoP trajectories in the presence of perturbations under variable step timing. CoM: Center of Mass; CoP: Center of Pressure.

Figure 9.

The CoM and CoP trajectories in the presence of perturbations under fixed step timing. CoM: Center of Mass; CoP: Center of Pressure.

Figure 10.

Variable step timing, the position and speed of CoM. CoM: Center of Mass.

Figure 11.

Fixed step timing, the position and speed of CoM. CoM: Center of Mass.

Comparisons with the existing methods

In the second scenario, we compare the robustness of the proposed approach with that of Herdt et al. ²⁵ The proposed walking pattern generation approach by Herdt et al. ²⁵ and variations ^2,4,13 of this approach are standard walking pattern generators. Herdt et al. ²⁵ calculates CoM trajectory and real-time landing position in a fixed prediction horizon and realizes the travel of predetermined speed. Here, we apply the same parameters to both methods and calculate the maximum perturbation velocities that each method can recover from different directions (see Figure 12). The measured data are based on the system parameters in the previous section. The robot has a ground contact as shown in Figure 12(a) at the beginning of the disturbance action of the 5th step, when the humanoid is in double support. At the beginning of the walking, both methods set the step timing to 0.5 s, $v_{x\text{ref}}=0.4$ m/s, and the disturbance in the forward field of vision (0, 180°) was applied when the humanoid reached the 5th step. That is to say, the forward is 90° direction of forward vision (positive direction of the x-axis), the 0° of forward vision is negative direction of the y-axis, and 180° is positive direction of the y-axis.

Figure 12.

Contrast of maximum restorable disturbance velocities and recovery time in the front horizon. (a) Definition of frontal horizon, (b) maximum recoverable speed in the forward horizon, and (c) recovery time in the forward horizon.

As can be observed in Figure 12, our approach with time adjustment is able to recover from much more severe pushes compared with Herdt et al. ²⁵ When the step time is fixed, because the robot is in double support, the resistance to disturbance in the vertical direction of the two-foot line is the weakest, and that in the direction of the two-foot line is the strongest, as shown in the green line in the Figure 12(b). With the method proposed in this article, the step timing can be adjusted, and the speed increases when the x-direction is disturbed. Under the restriction of equation (25), the step duration decreases, and the robot moves frequently, which enhances the anti-disturbance ability in the x-direction. However, in the y-direction, the minimum distance between the feet is within the range from 0.2 m to 0.3 m, because cross feet is not allowed. Even if the step time is reduced to a very small amount, it still can not produce a large CoP movement in space, as shown in Figure 13. So in the y-axis direction, the feasible area is more limited than the other direction.

Figure 13.

The instability of the humanoid lateral plane. (a) Definition of frontal horizon, (b) maximum recoverable speed in the forward horizon, and (c) recovery time in the forward horizon.

Figure 12(b) compares the advantages of this method presented in this article in terms of spatial. Figure 12(c) shows the advantages of this method in terms of time. In the forward 0–180° field of view, the common recoverable disturbance velocity is as shown in the shaded area in Figure 12(b). The recovery time corresponding to the two methods is shown in Figure 12(c). It can be seen that the recovery time used in this method is less than Herdt et al. ²⁵ when the two methods executed in the same disturbance speed in all directions. This is because the method reduces the stride time after the disturbance, but the number of steps used for recovery does not change. This in turn leads to a reduction in the recovery time used. Overall, the variable step duration method has an improved robustness as compared with the fixed step timing method.

Push recovery simulation

For push recovery, similar to human, the robot could take a big step to restore stability, or step by step to stability. As shown in Figure 14, the CoM has a initial velocity of 0.5 m/s in the x-direction and a initial velocity of 0.2 m/s in the y-direction, ${\Delta }_{dx}=0.5,{\Delta }_{dy}=0.3$ . By utilizing the scheme shown in Figure 7(a) under the disturbance, it can be seen that the CoM and CoP position appear overshoot, then gradually tend to stabilize. At the same initial condition, using the scheme shown by Figure 7(b), the CoM gradually stabilizes after many steps, as shown in Figure 15, without overshoot. Hence, it could be known that multi-step recovery is helpful to overcome the large disturbance and is more robust than one-step method.

Figure 14.

One-step stabilization mode for push recovery.

Figure 15.

Multi-steps stabilization mode for push recovery.

Discussion

Through the above simulation and comparative analysis, the following advantages could be achieved with the proposed method, as compared to the other existing methods,

Flexibility: The proposed method provides a variety of schemes to deal with the QP problem in the dynamic balance for humanoid robots. For the same task, there are different optimization models. While for different tasks, only minor modifications are needed to use similar optimization models with strong flexibilities.

Robustness: The simulation results show that the variable step duration method is more robust than the MPC approach with several preview steps without timing adjustment. In our method, a nominal step duration is set to allow the next prediction horizon to deviate from it in a linear constraint after disturbances, and then solve a convex optimization problem by looking at the next step location and timing.

Similarity with people: The gait generated by the method, in the single support phase, the CoP moves stably forward in the forward direction of the foot, which is quite similar to human walking. A similar approach to human is also used for the perturbation process, which results in larger CoP movements in a shorter period of time, and suppresses CoM divergence, while minimizes the speed and acceleration tracking errors.

It is worth noting that, since the main focus of this article is to test the feasibility of the proposed MPC algorithm onto the robot, we verified its robustness on flat ground only. When considering some other much more complicated forms of disturbances, like the uneven plane or walking on the slope, however, there may exist various other factors influencing on the performances of the robot. Those factors could impose a huge burden on the performance analysis of biped robots. Therefore, one aspect of our next-step work is to explore and improve the online anti-disturbance performances of the proposed algorithm under complex terrain.

Conclusion

In summary, we proposed an MPC method, which uses LIP as the motion model and the change of CoP as the input to minimize the variation of the CoP in the objective function, and thus guarantees the stable CoM and CoP trajectories, realizes the speed change and step duration change under the action of disturbance for humanoid robots. Compared with the fixed step duration method, the variable step distance and step timing method could further improve the compensation ability of the perturbations. Meanwhile, the flexible task of the robot is realized owing to the diversity of constraint settings. Moreover, the objective function can choose a different form by modifying some of the constraints that match it.

This study can be regarded as preliminary work for biped robots to enter human life. Our next-step work focuses mainly on the following aspects. First, we are to verify such proposed method on our lab-customized humanoid robot as shown in Figure 3, and then we would like to further extend our method to a three-dimensional space by considering the optimization in z-axis direction. Last but not least, the method could also be extended to uneven ground, within the slope or step environment to mimic the real human living environment, to further improve its adaptability.

Supplemental Material

Supplemental Material, files - Flexible model predictive control based on multivariable online adjustment mechanism for robust gait generation

Supplemental Material, files for Flexible model predictive control based on multivariable online adjustment mechanism for robust gait generation by Sheng Dong, Zhaohui Yuan, Xiaojun Yu, Muhammad Tariq Sadiq, Jianrui Zhang, Fuli Zhang and Cheng Wang in International Journal of Advanced Robotic Systems