Gait optimization and energy-based stability for biped locomotion using large-scale programming

Abstract

This paper presents a gait optimization method to generate the locomotion pattern for biped and discuss its stability. The main contribution of this paper is a newly proposed energy-based stability criterion, which permits the dynamic stable walking and could be straight-forwardly generalized to different locomotion scenarios and biped robots. The gait optimization problem is formulated subject to the constraints of the whole-body dynamics and kinematics. The constraints are established based on the modelling of bipedal hybrid dynamical systems. Following the whole-body modelling, the system energy is acquired and then applied to create the stability criterion. The optimization objective is also established on the system energy. The gait optimization is solved by being converted to a large-scale programming problem, where the transcription accuracy is improved via the spectral method. To further reduce the dimensionality of the large-scale problem, the whole-body dynamics is re-constructed. The generalization of the optimized gait is improved by the design of feedback control. The optimization examples demonstrate that the stability criterion naturally leads to a cyclic biped locomotion, though the periodicity was not previously imposed. Two simulation cases, level ground walking and slope walking, verify the generalization of the stability criterion and feedback control. The stability analyses are carried out by investigating the motions of centre of gravity and centre of pressure. It is revealed that if the tracked speed is above 0.3 m/s or the biped accelerates/climbs the slope, the stability criterion accomplishes the dynamic stable walking, where zero moment point criterion is not strictly complied.

1. Introduction

Over the last decades, the legged robots have attracted more attention than other categories of mobile robots, since they can adapt to more complicated terrains and surroundings. Among legged robots, the bipedal robots have the least practical application, due to many research difficulties. The biped itself is a highly complex system with poor stability of locomotion. Its stability region is even small if the biped is designed to achieve high flexibility and mobility. To overcome this problem, many researchers are dedicated to study how to generate a reliable gait for bipedal walking [1–5]. The gait of bipedal robot describes how it moves in one or several walking steps using the robot states.

Gait planning based on one or several optimal criteria is one popular research methodology for walking pattern generation of the bipedal robot. Chow and Jacobson were the first to introduce the optimal control method to prescribe the hip and knee trajectory [6]. Energy expenditure were minimized in their study. Considering the constraints at the beginning and the end of the walking step, Bessonnet et al. applied Pontryagin’s Maximum Principle (PMP) to optimize the gait motion of a seven-link planar biped [7]. The optimal control-based methods guarantee the principle of optimality, but have difficulties in chandling time-dependent constraints. Therefore, finding an optimal gait is transformed into a numerical programming problem [8–11]. To solve the numerical programming problem, both gradient-based methods and heuristic-based methods have their pros and cons. The heuristic-based methods, such as Central Pattern Generators (CPG) [12], Genetic Algorithm (GA) [13], Particle Swarm Optimization (PSO) [14], JAYA [15], and Whale Optimization [16] etc., perform outstanding when optimizing the black-box type systems, due to no need of explicit expressions of problems.

The gradient-based method is preferred for bipedal robot since the system can be explicitly expressed, thereby the optimization being more efficient. The explicit dynamics of bipedal robots is a non-linear hybrid model which is tricky to acquire the solution. Most active research has converged to formulating the inverse dynamics as a Quadratic Programming (QP) problem [17–19]. Thus, the optimization of joint torque and contact force can be achieved simultaneously. Even managing all kinematics constraints, the inverse dynamics problem is still solvable in real time if formulated as a QP problem [17]. Currently, the QP-based optimization method is applied to solve the bipedal state at one instant time, yet cannot provide the complete gait of one or several footsteps ahead with one optimization. The gait generation can predict future locomotion if the Model Predictive Control (MPC) scheme is introduced to the trajectory optimization [20]. The MPC provides the capability to consider the step locations as design variables of the gait optimization, thereby being able to predict the whole foot trajectory along specified time horizon [21]. Normally, the MPC is mainly adopted to solve linear or convex problems, due to its expensive computational process. For bipedal gait generation, the whole-body hybrid dynamics would not be integrated into MPC-based optimization in most studies [20–22]. Guo et al. constructed a whole-body MPC problem via synthesizing the gait library that is optimized off-line [23].

One potential scheme to include the whole-body hybrid dynamics is to formulate the gait optimization as a non-linear programming (NLP) problem [24–26]. The NLP formulation can parameterize all system states within a period of time and solve the complete gait of many steps in one optimization. Furthermore, since all system states can be parametrized using the time variable, the jerky variations of optimal results could be avoided. Continuing on ref. [7], Bessonnet et al. converted the gait control problem to a NLP and adopted the fourth-order spline function to construct the joint trajectory [25]. Similar to MPC, the NLP method heavily suffers from high computational expense (even higher due to involving the whole-body dynamics). Under the shadow of expensive computing cost, only a few NLP-based studies are concern with stability criterion [27].

To guarantee the stable locomotion of the biped, many researchers integrated the stability criterion into the gait optimization. Similar to non-optimization-based locomotion control, the zero moment point (ZMP) is the most widely used criterion [13, 28, 29] for walking pattern optimization. Generally speaking, the ZMP can guarantee static stable locomotion, but its dynamic stability has not been fully discussed [30]. Another broadly adopted stability criterion is the capture point or capture region that are proposed to decide how to take steps to stabilize the bipeds [1]. The capture point is also referred as divergent component of motion (DCM) as it indicates the unstable dynamics of centre of mass (CoM). Here, the CoM dynamics is simplified using well-known linear inverted pendulum model (LIPM). Therefore, typically the full-order dynamics would not be involved when stabilizing the DCM [31, 32]. Mesesan et al. designed a whole-body torque controller to track the reference resulting from DCM criterion. A middle ground found between full-order model and a simplified model is centroidal dynamics which employs the linear and angular momentum of CoM as the stability criterion. Different from the simplified model like LIPM, the centroidal model respects the bipedal real dynamics, thereby be suitable to highly dynamic motion. Chignoli et al. accomplished acrobatic behaviours of a humanoid robot by regulating the centroidal momentum [33]. Many researches prove that this stability criterion endues the computational tractability and dynamic expressiveness [33–35].

The research group led by Grizzle creatively introduced the full-order hybrid zero dynamics (HZD) to study the dynamic stability of bipedal robots. The biped robot was modelled as a hybrid dynamical system. The system stability was achieved by the feedback linearization, which can enforce the system to reach its zero dynamics [36]. Owing to HZD based gait optimization, several well-known biped robots were established and accomplished the dynamic stability, e.g. ATRIAS [37] and Cassie [38]. Yet HZD-based gait design were mainly applied to underactuated robots or the robots were formulated as underactuated. Lyapunov theorem is another golden criterion for dynamic stable walking. Kakaei and Salarieh proposed a hybrid-dynamics-based-control for five-link underactuated-biped robot and the stability of generated locomotion was proven using Lyapunov function [39]. To achieve adaptive walking of semi-passive bipeds, Liu et al. stabilized the transition gait established on discrete Lyapunov control. Similar to HZD-based gait optimization, Lyapunov-based dynamic locomotion planning was mostly applied to underactuated bipeds.

This paper presents a gait optimization method to generate the locomotion pattern for bipeds. The bipedal locomotion pattern is described using the robot position, orientation, and joints. The gait optimization problem will be constrained to the whole-body dynamics and kinematics. The main contribution of this paper is a newly proposed stability criterion, which permits the dynamic stable walking. The proposed stability criterion is system energy-based, thereby involving the full-order or whole-body dynamics. By combining the stability criterion and whole-body dynamics, a large-scale programming problem that could predict the gait pattern of several steps ahead is formulated. To reduce the computational expense of the large-scale programming, the whole-body dynamics is re-constructed. Furthermore, there is no further assumption on establishing the whole-body dynamics so that the delivered gait optimization method could be effortlessly generalized to different-legged robots, either underactuated, or fully actuated ones.

The paper starts by introducing the developed biped robot – Cosmos, with its configurations and coordinates. In Section 2, the gait generation problem of Cosmos is formulated by firstly modelling the whole-body dynamics and kinematics. Then, the new energy-based stability criterion is proposed and extended to build the optimization objective. Section 3 also converts the original gait generation problem to a large-scale programming problem, where the whole-body dynamics and the optimal variables are re-formulated. The discretization methodology is presented afterwards. Section 4 describes the feedback control that combines the gait optimization and task/joint space-based control. Lastly, the optimization results and the simulation outcomes are provided and compared in Section 5.

2. Problem formulation

The Cosmos Robot has two configurations: bipedal and wheeled, as shown in Fig. 1. The bipedal configuration is the main concern of this paper due to its high adaptability to the terrain and its scientific difficulty. The detailed bipedal configurations of Cosmos Robot, including joints, links, and actuators, are presented in Appendix A.1.

Figure 1. Cosmos Robot that has two configurations: bipedal and wheeled.

2.1. Hybrid dynamics model

The gait of legged robot describes how it moves in one or several walking steps. The bipedal gait is mainly represented by the corresponding coordinates. The Cartesian coordinates are the most widely applied representation for studying complex and high Degrees-of-Freedom (DoF) systems. However, the generalized/joint coordinates are favoured in robotics, since it is computationally more efficient and accurate and can avoid the violation of joint constraints [40]. For Cosmos Robot, the minimum dimensionality of the configuration space is 10, and the corresponding coordinates are designed as

(1) \begin{equation}\boldsymbol{{q}}\,:\!=\, [p_{x}, p_{y}, p_{y},\varphi, \theta, \psi,q_{1L}, q_{2L}, q_{3L}, q_{4L}, q_{5L},q_{1L}, q_{2L}, q_{3L}, q_{4L}, q_{5L}]^{\mathrm{T}},\end{equation}

in which $p_{x}, p_{y}, p_{y},$ $\varphi, \theta, \psi$ are appended to directly express the robot position and orientation in world inertia frame. While $q_{1}$ , $q_{2}$ , $q_{3}$ are the hip roll, hip yaw, and hip pitch angles, respectively, $q_{4}$ is the knee pitch angle, and $q_{5}$ is the ankle pitch angle. The coordinates $q_{1}$ , $q_{2}$ , $q_{3}$ , $q_{4}$ , $q_{5}$ also indicates the relative angles between links of the robot, i.e., joint angles (as shown in Fig. 16).

Several assumptions are made before the modelling of locomotion dynamics of the bipedal Cosmos. It is first assumed that the stance leg stays pinned (not penetrate) to the ground and does not move when another leg is swinging [41]. Secondly, impact is assumed that happens instantaneously. Besides, only velocities vary at impact moment, while no instantaneous changes in position [42].

In the paper, the biped walking is divided into swing phase and impact phase, which are recognized as a continuous process and discrete process, respectively. The impact phase is defined when the swing leg touches the ground [43, 44]. The continuous dynamics of the floating-base model of biped robot can be addressed in the form of Lagrange as

(2) \begin{equation}\boldsymbol{{D}}\!\left(\boldsymbol{{q}}\right)\ddot{\boldsymbol{{q}}}+\boldsymbol{{C}}\!\left(\boldsymbol{{q}},\dot{\boldsymbol{{q}}}\right)+\boldsymbol{{G}}\!\left(\boldsymbol{{q}}\right)=\boldsymbol{{Bu}}+\boldsymbol{{J}}(\boldsymbol{{q}})^{\mathrm{T}}\boldsymbol{\lambda },\end{equation}

where $\boldsymbol{{D}}$ is the mass–inertia matrix, $\boldsymbol{{C}}$ is the Coriolis matrix, and $\boldsymbol{{G}}$ relates to the gravity. On the other side of dynamics equation, $\boldsymbol{{u}}$ is driver (motor) torque vector, and $\boldsymbol{{B}}$ denotes how the driver torques are mapped to joint torques. The mechanism design of Cosmos Robot determines that the mappings are independent between different motors. Meanwhile, the motor torques cannot directly affect the robot position and orientation. Thus, the matrix $\boldsymbol{{B}}$ can be represented by $\boldsymbol{{B}}=[\begin{array}{c@{\quad}c} \mathbf{0}_{6\times 10} & \boldsymbol{{R}}_{10\times 10} \end{array}]^{\mathrm{T}}$ , where $\boldsymbol{{R}}_{10\times 10}$ only depends on the gear ratio of motors. $\boldsymbol{\lambda }$ is contact wrench vector which contributes to variation of the system energy, while $\boldsymbol{{J}}$ is the Jacobian matrix related to contact and matches wrench to joint torques.

When the biped is walking, it is subject to a certain of holonomic constraints [45]. For instance, feet and ground cannot interpenetrate. Also, the stance foot is assumed to be fixed to the ground. The above kinematic constraints can be addressed using

(3) \begin{equation}\boldsymbol{{J}}\!\left(\boldsymbol{{q}}\right)\dot{\boldsymbol{{q}}}=\mathbf{0},\end{equation}

(4) \begin{equation}\boldsymbol{{J}}\!\left(\boldsymbol{{q}}\right)\ddot{\boldsymbol{{q}}}+\dot{\boldsymbol{{J}}}\!\left(\boldsymbol{{q}},\dot{\boldsymbol{{q}}}\right)\dot{\boldsymbol{{q}}}=\mathbf{0},\end{equation}

where $\boldsymbol{{J}}$ is the Jacobian matrix for active contacts.

Refer to the discrete phase, the assumptions are (1) impact instantaneously happening and (2) the robot positions keeping while the velocities discontinuing. Define the robot state before and after impact are $\boldsymbol{{q}}^{-}$ and $\boldsymbol{{q}}^{+}$ , while $\boldsymbol{{q}}^{-}=\boldsymbol{{q}}^{+}$ as the position is invariant through impact. The velocity part is expressed as $\dot{\boldsymbol{{q}}}^{-}$ and $\dot{\boldsymbol{{q}}}^{+}$ . The dynamics for the discrete phase is represented by the impact map, whose inferences are provided in Appendix A.2.

Considering the continuous and discrete dynamics as a whole, the bipedal robot can be taken as a typical hybrid dynamical system. The overall hybrid model and its dynamics are governed by Eq. (5):

(5a) \begin{equation}\boldsymbol{{D}}\!\left(\boldsymbol{{q}}\right)\ddot{\boldsymbol{{q}}}+\boldsymbol{{C}}\!\left(\boldsymbol{{q}},\dot{\boldsymbol{{q}}}\right)+\boldsymbol{{G}}\!\left(\boldsymbol{{q}}\right)=\boldsymbol{{Bu}}+\boldsymbol{{J}}(\boldsymbol{{q}})^{\mathrm{T}}\boldsymbol{\lambda },\end{equation}

(5b) \begin{equation}\dot{\boldsymbol{{q}}}^{+}=\Delta _{\dot{\boldsymbol{{q}}}}(\boldsymbol{{q}}^{-})\dot{\boldsymbol{{q}}}^{-},\end{equation}

where the definition of $\Delta _{\dot{\boldsymbol{{q}}}}(\boldsymbol{{q}}^{-})$ can be found in Appendix A.2.

2.2. Kinematics constraints

The kinematic constraints of biped walking can be separated into holonomic and non-holonomic constraints. The holonomic constraints of stance foot are provided in the non-holonomic form, i.e., Eqs. (3) and (4). As for swing foot, the constraint of foot clearance at half step-time is given by Eq. (6), also in a non-holonomic form. The position of torso related to the stance foot is constrained in a range using Eq. (7). Equations (8) and (9) show that the pitch of swing foot is preferred to be slightly tilted, and the pitch of robot torso is allowed to be swayed in a small range.

(6) \begin{equation}p_{\mathrm{sw},\mathrm{clr}}\!\left(\frac{t_{s}}{2}\right)\geq c_{\mathrm{sw},\mathrm{clr}},\end{equation}

(7) \begin{equation}p_{b,\min }\leq p_{b}\leq p_{b,\max },\end{equation}

(8) \begin{equation}0\leq \theta _{\mathrm{sw}}\leq c_{\theta,\mathrm{sw}},\end{equation}

(9) \begin{equation}-c_{\theta,b}\leq \theta _{b}\leq c_{\theta,b}.\end{equation}

Other non-holonomic constraints at velocity level include the average sagittal/lateral velocity of robot torso and the impact velocity of swing foot, given as follows:

(10) \begin{equation}\overline{v}_{b,x}=c_{b,x}, \overline{v}_{b,y}=0,\end{equation}

(11) \begin{equation}v_{\mathrm{sw},x}(t_{s})=0, v_{\mathrm{sw},y}(t_{s})=0, v_{\mathrm{sw},z}\!\left(t_{s}\right)\leq 0.\end{equation}

On the other hand, the ground reaction force must respect the friction cone [46]. Following this, friction coefficient $\mu$ is introduced to formulate:

(12) \begin{equation}\frac{\left\| \boldsymbol{{F}}_{\mathrm{imp},\text{hori}}\right\| }{F_{\mathrm{imp},z}}\leq \mu,\end{equation}

where $F_{\mathrm{imp},z}$ is the magnitude of impact force in vertical direction and $\boldsymbol{{F}}_{\mathrm{imp},\text{hori}}$ denotes force vector in horizon plane. Note that most of constraints are expressed in the non-holonomic form.

Lastly, the boundary limits of optimal variables are imposed. The maximum and minimum values of actuators’ torques are limited to ensure that the optimized results can be practically implemented on the actual robot.

3. Gait optimization and stability

The gait of legged robot describes how it moves in one or several walking steps. The bipedal gait is mainly represented by defined coordinates, i.e. the system state. The robot state at each instant is a numerical solution of the equations of motion of biped dynamics and kinematics. Typically, the solution of system of non-linear differential equations that is subject to various equality and inequality constraints, can be found by numerical optimization methods. For time-varying systems, the scale of optimization could be largened significantly, since system states during a time period are generally requested to be solved simultaneously.

In the paper, the optimization of biped gait is transferred to a large-scale programming problem.

3.1. Optimal variables

The variables defining the biped walking gait include step time $t_{s}$ , driver torque $\boldsymbol{{u}}$ , system state variables $\boldsymbol{{q}}$ , and its derivatives $\dot{\boldsymbol{{q}}}$ , $\ddot{\boldsymbol{{q}}}$ . Besides, some other unsolved variables are needed to complete the optimization problem, like impact force $\boldsymbol{{F}}_{\textbf{imp}}$ .

From Eq. (5), the continuous dynamics of bipedal robots can be re-written to:

(13) \begin{equation}\ddot{\boldsymbol{{q}}}=\boldsymbol{{D}}\!\left(\boldsymbol{{q}}\right)^{-1}\left(\boldsymbol{{Bu}}+\boldsymbol{{J}}(\boldsymbol{{q}})^{\mathrm{T}}\boldsymbol{\lambda }-\boldsymbol{{C}}\!\left(\boldsymbol{{q}},\dot{\boldsymbol{{q}}}\right)-\boldsymbol{{G}}\!\left(\boldsymbol{{q}}\right)\right).\end{equation}

Defining new system state $\boldsymbol{{x}}=[\boldsymbol{{q}},\dot{\boldsymbol{{q}}}]^{\mathrm{T}}$ , Eq. (13) can be transferred to the state-space formulation:

(14) \begin{equation}\dot{\boldsymbol{{x}}}=f\!\left(\boldsymbol{{x}},\boldsymbol{{u}}\right)=\left[\begin{array}{c} \dot{\boldsymbol{{q}}}\\[4pt] \ddot{\boldsymbol{{q}}} \end{array}\right]=\left[\begin{array}{c} \dot{\boldsymbol{{q}}}\\[5pt] \boldsymbol{{D}}\!\left(\boldsymbol{{q}}\right)^{-1}\left(\boldsymbol{{Bu}}+\boldsymbol{{J}}(\boldsymbol{{q}})^{\mathrm{T}}\boldsymbol{\lambda }-\boldsymbol{{C}}\!\left(\boldsymbol{{q}},\dot{\boldsymbol{{q}}}\right)-\boldsymbol{{G}}\!\left(\boldsymbol{{q}}\right)\right) \end{array}\right].\end{equation}

The system states of pre- and post-impact are also re-defined as $\boldsymbol{{x}}^{-}$ and $\boldsymbol{{x}}^{+}$ . The impact map is then extended to

(15) \begin{equation}\boldsymbol{{x}}^{+}\,:\!=\,\Delta \!\left(\boldsymbol{{x}}^{-}\right)=\left[\begin{array}{c} \Delta _{\boldsymbol{{q}}}\boldsymbol{{q}}^{-}\\[4pt] \Delta _{\dot{\boldsymbol{{q}}}}\!\left(\boldsymbol{{q}}^{-}\right)\dot{\boldsymbol{{q}}}^{-} \end{array}\right],\end{equation}

where $\Delta _{\boldsymbol{{q}}}$ is an identity map since $\boldsymbol{{q}}^{-}=\boldsymbol{{q}}^{+}$ .

In sum, the state equations of biped dynamics can be transferred to

(16) \begin{align}\dot{\boldsymbol{{x}}} & =f\!\left(\boldsymbol{{x}},\boldsymbol{{u}}\right),\nonumber\\[5pt]\boldsymbol{{x}}^{+}& \,:\!=\,\Delta (\boldsymbol{{x}}^{-}),\end{align}

with new system state $\boldsymbol{{x}}$ . By introducing the new system state $\boldsymbol{{x}}$ , state variables in optimization are transferred from $[\boldsymbol{{q}},\dot{\boldsymbol{{q}}},\ddot{\boldsymbol{{q}}}]^{\mathrm{T}}$ to $[\boldsymbol{{q}},\dot{\boldsymbol{{q}}}]^{\mathrm{T}}$ . Therefore, the dimensionality of optimization space is reduced significantly, considering all time-varying variables would be discretized during the time duration of walking. In other words, the main benefit of reformulated system model is to reduce the computational expense of gait optimization.

The optimization of a physical system or its motion must be subject to its own system dynamics and kinematics. For biped, its locomotion meets the hybrid dynamics governed by Eq. (16) and kinematics constraints addressed by Eqs. (6)–(12).

3.2. Stability criterion

The dynamics and kinematics constraints can guarantee that the real bipedal robot could complete the optimized gait, but the stability of the biped performing the optimized gait is unknown.

As mentioned in Section 1, not only the stability of biped gait is difficult to define, but also the defined stability criterion is generally not easy to integrated into the gait optimization. In this paper, a straightforward principle is proposed to decide whether the optimized gait can be considered to be stable or not.

For a moving object, the system energy would not keep increasing as the time goes by, if the system is stable. Referred to the bipedal system, this generally acknowledged truth can be rephased. The system is stable if its kinetic energy of each cycle does not continue growing. Put it in a way to emphasize the sufficiency, if the biped is initially stable and the gait of each step contributes negatively to the system kinetic energy, the system cannot be diverging.

For a bipedal robot, the energy can only be changed by the work that is done by the drives and done by the impact forces, if the interpenetration of stance foot and ground is not permitted. In the swing phase, the variations of kinetic energy and potential energy are all powered from the drives (motors). While for the impact phase, only the kinetic energy is taken into consideration. The reason is that the potential energy keeps since the position of bipeds is invariant through impact.

The change of energy of bipedal system in one step is given as

(17) \begin{equation}\Delta E_{K,\mathrm{sw}}+\Delta E_{P,\mathrm{sw}}=\int _{0}^{t_{s}}\boldsymbol{{Bu}}\!\left(t\right)\cdot \dot{\boldsymbol{{q}}}\!\left(t\right)dt,\end{equation}

(18) \begin{equation}\Delta E_{K,\mathrm{imp}}=\frac{1}{2}m{\boldsymbol{{v}}_{c}}^{+}\cdot {\boldsymbol{{v}}_{c}}^{+}-\frac{1}{2}m{\boldsymbol{{v}}_{c}}^{+}\cdot {\boldsymbol{{v}}_{c}}^{+}+\frac{1}{2}{\boldsymbol{\omega }_{c}}^{+}\cdot \left(\boldsymbol{{I}}\cdot {\boldsymbol{\omega }_{c}}^{+}\right)-\frac{1}{2}{\boldsymbol{\omega }_{c}}^{-}\cdot \left(\boldsymbol{{I}}\cdot {\boldsymbol{\omega }_{c}}^{-}\right),\end{equation}

referred to swing phase and impact phase, respectively. In above energy equations, $m$ and $I$ are the lumped mass and inertia matrix to CoG, while $\boldsymbol{{v}}_{c}$ and $\boldsymbol{\omega }_{c}$ are the linear and angular velocity of CoG. Note that dot product $\cdot$ in this paper follows the rule: $\boldsymbol{{a}}\cdot \boldsymbol{{b}}=\boldsymbol{{a}}^{\mathrm{T}}\boldsymbol{{b}}$ . Meanwhile, $\Delta E_{K,\mathrm{sw}}$ and $\Delta E_{K,\mathrm{imp}}$ , respectively, denote the kinetic energy change for the swing phase and impact phase, while $\Delta E_{P,\mathrm{sw}}$ represents the potential energy change.

For the commonest scenario, such as level ground walking, the variation of potential energy can be neglected so that the change of energy only includes the kinetic energy. On the other hand, for the cases like climbing a stair or walking up/down a slope surface, the potential energy does not keep constant. The kinetic energy change summed by the swing phase and impact phase can be easily computed, if only the change of potential energy can be pre-acquired or estimated with the aid of the navigation system.

Lemma 1. If the change of system kinetic energy of biped during one step is not positive, the biped is stable in this step or the gait of this step is “stable”, i.e.:

(19) \begin{equation}\Delta E_{K,\mathrm{sw}}+\Delta E_{K,\mathrm{imp}}\leq 0.\end{equation}

3.3. Objective

The above-proposed stability criterion only guarantees the non-positive energy variation between two continuous steps. However, the energy variation during the whole walking scenario has not been taken into consideration. The variation of the total energy expenditure is minimized in this study to achieve the gradual energy change [11]. The cost function of gait optimization can be expressed by

(20) \begin{equation}J=\int _{t_{0}}^{t_{f}}\left\| \Delta E\!\left(t\right)\right\| ^{2}dt,\end{equation}

where $\Delta E(t)$ is the change of energy at each instant time, calculated using Eqs. (17) and (18). The time $t\in [t_{0}, t_{f}]$ , where $t_{f}$ is the final time and the initial time $t_{0}$ generally starts at zero.

3.4. Discretization

Direct transcription method is adopted in the paper to discretize the continuous dynamics. As mentioned in Section 3, optimization of the biped gait would form a large-scale programming. Therefore, the spectral method is employed in this study, instead of the basic Runge–Kutta transcription.

The spectral method is a non-finite difference technique, in theoretical which can reach the global accuracy. The Legendre pseudospectral transcription is one of the most classic spectral methods. Under the scheme of Legendre pseudospectral method, general problems on the interval $t\in [t_{0}, t_{f}]$ can be scaled to time interval $\tau \in [{-}1,1]$ using the relation [47]:

(21) \begin{equation}t=\frac{\left(t_{f}-t_{0}\right)\cdot \tau +\left(t_{f}+t_{0}\right)}{2}.\end{equation}

Under the time scaling, the continuous dynamics of biped walking is re-scaled to:

(22) \begin{equation}\frac{d\boldsymbol{{x}}}{d\tau }=\frac{\left(t_{f}-t_{0}\right)}{2}\cdot f\!\left(\boldsymbol{{x}}\!\left(\tau \right),\boldsymbol{{u}}\!\left(\tau \right)\right).\end{equation}

On the other hand, the states and controls are approximated using a set of Lagrange interpolating polynomials, given by [47]:

(23) \begin{align}\boldsymbol{{x}}\!\left(\tau _{k}\right) & \approx \sum _{i=1}^{N_{c}}\boldsymbol{{x}}\!\left(\tau _{i}\right)\cdot L_{\boldsymbol{{i}}}\!\left(\tau _{k}\right),\nonumber\\\boldsymbol{{u}}(\tau _{k}) & \approx \sum _{i=1}^{N_{c}}\boldsymbol{{u}}(\tau _{i})\cdot L_{\boldsymbol{{i}}}(\tau _{k}),\end{align}

where $\tau _{i}$ , $i=1,\ldots,k$ , are Lagrange–Gauss–Lobatto points and $N_{c}$ is number of collocation points. $\boldsymbol{{L}}_{\boldsymbol{{i}}}(\tau _{k})$ represents the Lagrange polynomials that satisfy:

(24) \begin{equation}L_{\boldsymbol{{i}}}\!\left(\tau _{k}\right)=\delta _{ki},\end{equation}

where $\delta _{ik}$ denotes Kronecker delta.

Following this, the derivative of state can also be approximated using the interpolating polynomial, shown as [47]:

(25) \begin{equation}\frac{d\boldsymbol{{x}}}{d\tau }\left(\tau _{k}\right)\approx \sum _{i=1}^{N_{c}}\boldsymbol{{x}}\!\left(\tau _{i}\right)\cdot \boldsymbol{{D}}_{ki},\end{equation}

where $\boldsymbol{{D}}_{ki}=\dfrac{dL_{\boldsymbol{{i}}}}{d\tau }\left(\tau _{k}\right)$ , named as derivative matrix. Considering this, the dynamics constraint can be re-written to:

(26) \begin{equation}\sum _{i=1}^{N_{c}}\boldsymbol{{x}}\!\left(\tau _{i}\right)\cdot \boldsymbol{{D}}_{ki}=\frac{\left(t_{f}-t_{0}\right)}{2}\cdot f\!\left(\boldsymbol{{x}}\!\left(\tau _{k}\right),\boldsymbol{{u}}\!\left(\tau _{k}\right)\right).\end{equation}

Similarly, the object function is represented using discretization form, such as:

(27) \begin{equation}J=\frac{\left(t_{f}-t_{0}\right)}{2}\sum _{k=1}^{N_{c}}\left\| \Delta E\!\left(\tau _{k}\right)\right\| ^{2}\cdot w_{k},\end{equation}

where $w_{k}$ stands for weight for each collocation node.

4. Feedback control

The diagram of feedback control flow for the Cosmos locomotion is shown in Fig. 2. The biped gaits in this paper are optimized off-line and integrated into walking control as a gait library. The off-line optimization-based gait generation improves the capability of generalization to different locomotion scenarios [38]. Also, the off-line-based method is especially suitable for platforms with limited computational resources.

Figure 2. Architecture diagram of feedback control for biped locomotion.

The gait library is re-formulated as a Gait Generation module, as highlighted with red colour in Fig. 2. The Gait Generation block outputs target walking gait based on the target velocity and the robot velocity. The target velocity derives from the Tasks module, while the robot actual velocity is estimated using Velocity Estimator that is developed based on invariant extended Kalman filter.

For the sagittal plane, the target gait is decomposed to target joint position/velocity via Inverse Kinematics module. For the lateral direction, the lateral target position from the Tasks block provides the target joint position/velocity of roll and yaw joint via the Lateral Adjust module. Combining both the target and feedback information ( $q_{d}$ , $\dot{q}_{d}$ , $q$ , $\dot{q}$ ), the regulation control is implemented at the joint level to generate the torque command $\tau _{d}^{q}$ for the robot actuator. The torque command $\tau _{d}^{q}$ from the joint control is calculated using PD controller:

(28) \begin{equation}\tau _{d}^{q}=K_{p}^{q}\!\left(q_{d}-q\right)+K_{d}^{q}\!\left(\dot{q}_{d}-\dot{q}\right),\end{equation}

where $K_{p}^{q}$ and $K_{d}^{q}$ , respectively, denote P-gain and D-gain of the designed Joint Control module.

Note that in current study, only optimal results in sagittal plane are collected to constitute the gait library and guide the control. For the lateral direction, the cyclic movement is mainly achieved by stabilizing the robot attitude $\Psi =\{\varphi, \psi \}$ . The Posture Control module is also formulated using PD control algorithm:

(29) \begin{equation}\tau _{d}^{\Psi }=K_{p}^{\Psi }\!\left(\Psi _{d}-\Psi \right)+K_{d}^{\Psi }\!\left(\dot{\Psi }_{d}-\dot{\Psi }\right),\end{equation}

where $K_{p}^{\Psi }$ and $K_{d}^{\Psi }$ , respectively, denote P-gain and D-gain to output an attitude stabilization-based torque command $\tau _{d}^{\Psi }$ . Normally, the desired lateral attitude equals to zero for most locomotion cases. The additional torques $\tau _{d}^{G}$ are supplemented to torque commands, in order to compensate the gravitational term.

Lastly, the final torque command inputting to the Robot is the sum of outputs from the joint control, posture control, and gravity compensation blocks: $\tau _{d}=\tau _{d}^{q}+\tau _{d}^{\Psi }+\tau _{d}^{G}$ .

5. Results and analysis

This section provides the optimal results of walk pattern generation of Cosmos robot. The effects of kinematics constraints and stability criterion on optimal results are illustrated. Afterwards, the optimization results are collected and integrated into the feedback control of Cosmos biped. The simulation results and optimization results are compared in the end.

5.1. Optimization

As mentioned in Section 3.1, the optimal variables contain step time, impact force, driver torque, joint position and joint velocity. Their bounds are listed in Table I. The range of step time is estimated based on desired height of CoG of the Cosmos, set between (0.3, 0.5) s. The joint position must be positive since it reveals the relative angle between two connected links. Meanwhile, its value would not exceed $\pi$ due to mechanical limits. The joint velocity is also constrained considering the motor driver limits. The total weight of the biped would not be over 30 kg with the battery, so the impact force is limited to lower than 300 N.

Table I. Boundary of optimal variables.

Note that the ranges of some variables are ten/hundred times larger than other variables. It would raise difficulties in searching the variable space. Therefore, the normalization method is performed in the paper to scale the ranges of all variables to the same magnitude.

Other parameters used in kinematic constraint are also assigned. The foot clearance is set as 0.1 m, while the friction coefficient is 0.7. On the other hand, the pitch of swing foot and torso, and the body height are all restricted in a quite small range.

The optimization results are illustrated as follows. Fig. 3 to Fig. 6 show the gait optimization outcomes with four-step cycles. Each cycle can be categorized as the stance phase or swing phase for different foots, based on the foot position/velocity. Note that the separation of step cycles can be recognized using the step time of each cycle.

Figure 3. Position and velocity along $x,z$ -coordinates of left foot in optimization.

Figure 4. Position and velocity along $x,z$ -coordinates of centre of gravity in optimization.

Figure 5. Rotational velocities and torques of thigh, knee and ankle joints in optimization.

Figure 6. Energy in optimization with and without stability criterion.

5.1.1. Constraints

Figure 3 displays the position and velocity of the left foot with reference to the $x$ -axis and $z$ -axis, respectively. The gait starts with left-leg-stance state, then switch to left-leg-swing, and repeat switching with continuing two-step cycles. The figure demonstrates that the kinematic constraints of stance foot and swing foot are well-satisfied in the optimization. For instance, the zero positions and velocities of stance foot imply that the foot is fixed to the ground during its stance phase. With regards to the swing foot, the velocity along $x$ -axis rises from zero at the beginning of one step and declines to zero at the end of one step, which follows the requirement of constraint (11). Constraint (11) is imposed so that the velocity along $z$ -axis would not be positive at the end of one step, as proved in Fig. 3.

Figure 4 presents the position and velocity of the CoG with reference to the $x$ -axis and $z$ -axis, respectively. The average target speed of torso in gait optimization is set as 1 m/s and the height of torso only varies in a small range, as constrained in Section 2.2. As shown in Fig. 4, the resulted CoG speed is slightly lower than 1 m/s, about 0.9 m/s; while the height of CoG basically keeps constant. It is interesting that the CoG velocity trajectory from optimization is similar to the motion of simplified LIPM.

5.1.2. Stability

As illustrated in Section 3.2, the change of system energy in the continuous process mainly depends on rotational velocities and torques of robot joints at each instant. Here, the rotational velocities and torques of thigh, knee and ankle joints for left leg are plotted in Fig. 5. The rotational velocities for stance leg are much smaller than ones for swing leg; on the contrary, the torques for stance leg are much larger than ones for swing leg. The relative change of knee velocities and torques during two steps is the biggest among three joints. It is worth mentioning that the ankle joint is barely actuated if this leg swings in this phase.

Combining the system energy and its change for both swing and impact phase, its periodic behaviour is graphed in Fig. 6 with the black solid line. The system energy after impact process is certainly not higher than the one at the beginning. Therefore, the stability criterion that requests the change of system energy for one step not to be positive, is satisfied.

Another two optimizations not integrated with stability criterion are also presented in Fig. 6, using blue and cyan dotted lines, respectively. The optimization that enforces the periodic constraint, graphed by cyan dotted lines, is taken as a benchmark case. Note that the step time is an optimal variable so it would be diverse for different optimizations. Therefore, the ranges of time from three optimizations are re-scaled to [0, 1]. On the other hand, the objectives of the “blue” and “cyan” optimizations are replaced by the minimization of mean-squared error between CoG speed and target speed.

From the “blue” one, it is obvious that even if all dynamics and kinematics constraints are satisfied and the instant CoG speed is carefully regulated to around target value using new (speed-change-minimized) cost function, the optimized result does not necessarily exhibit periodicity for each step. Also, the changes in system energy vary for all four steps. Further state constraints between continuous steps are enforced to acquire the periodic “cyan” result. Still, the system energy of the “cyan” optimization slowly increases compared to its initial value. It is concluded that the proposed criterion constraint not only guarantee the gait stability, but also naturally leads to a cyclic biped locomotion, even though it does not take periodicity into consideration.

5.2. Simulation

The simulation completes the feedback control of Cosmos Robot(see Fig. 7). The diagram of feedback control that integrates offline optimized gait with low-level stabilized control is shown in Fig. 2. For sagittal plane, the gait library is implemented to generate reference joint trajectory, in order to achieve stable forward walking. On the other hand, the stabilization of lateral direction is accomplished by regulating roll and yaw attitude.

Figure 7. Cosmos biped in simulation environment.

5.2.1. Level ground walking

The test case 1 simulates a complete 20 s scenario where the robot walks on the flat surface. The scenario includes walking in place, accelerating, constant speed, deaccelerating, and finally walking in place. Different walking conditions are integrated to verify the generalization of the proposed optimization and feedback control methodology. The walking velocity command begins from zero, rises to 1 m/s, remains at the top and declines back to zero, as shown in Fig. 8.

Figure 8. Velocity along $x$ -coordinate and position along $z$ -coordinate of centre of gravity in the simulation of level ground walking.

As illustrated in Section 5.1, the average velocity along $x$ -axis of optimized gait would be slightly slower than the given target. This problem is solved by introducing the velocity feedback at the joint level control. Fig. 8 demonstrates that this policy succeeds and both instantaneous velocity and average velocity of CoG well follow the reference one.

The variation of CoG position/velocity with reference to $z$ -axis increases with the growth of the velocity along $x$ -axis, still limited between (−0.03,0.02) m. This indicates that the biped has difficult to comply with the kinematics of Linear Inverted Pendulum (LIP) when the robot speed is large.

Figure 9 presents the position of CoG, Centre of Pressure (CoP) and the boundary of stance foot, with reference to $x$ -axis. The upper and lower limit of the stance foot are determined by the foot size, which is 8.5 cm in length. If any position along $x$ -axis is smaller than the upper limit and larger than the lower limit, this position can be considered as in the support polygon of bipedal walking. Since the reference CoG velocity is axial symmetry regarding the time of 10 s (see Fig. 8), the curves of CoG positions are roughly central symmetrical with respect to the point (10, 5). Without loss of any features, the simulation results regarding half of the total time (between [0, 10] s) are plotted.

Figure 9. Position along $x$ -coordinate of centre of gravity, centre of pressure and boundaries of foot in the simulation of level ground walking.

During the time between 0 s and 4 s, the robot walks in a relatively slow velocity (lower than 0.3 m/s as given in Fig. 8). Figure 9 demonstrates that CoP falls inside the support polygon between [0, 4] s, when the robot experiences stepping-in-place and forward walking under 0.3 m/s. The robot CoG can also be projected inside the support polygon for walking in place. On the contrary, CoG is slightly ahead of the upper limit of stance foot when accelerating from 0 to 0.3 m/s.

The robot holds the average velocity of about 1 m/s during [9, 10] s. Different from slow velocity walking, the robot CoP coincides with the upper bound of support polygon during most of time for fast walking. Meanwhile, the robot CoG exceeds the upper limit of stance foot after half of stance time when the velocity remains around 1 m/s. While during the time between 5 s and 8 s, the exceeding time comes earlier since the robot walks in an acceleration. All these facts imply that the dynamic stable walking, in which ZMP criterion is not strictly complied, can be achieved with our proposed stability criterion.

Figure 10 presents the position of left foot, right foot and the whole robot, with reference to $x$ -axis. Note that the positions of foot are represented by the positions of ankle projected to the ground, while the body position stands for the whole robot position. Similar to the robot CoG, the body position falls near the edge of support polygon during most of walking time, with reciprocating support switch between left and right foot. Unlike the robot CoG, the body position is slightly behind the left and right foot at the stage of stepping-in-place (between [0, 4] s). It is indicated that the body position keeps behind the robot CoG if the target velocity is zero.

Figure 10. Position along $x$ -coordinate of left foot, right foot and body in the simulation of level ground walking.

5.2.2. Slope walking

The test case 2 carries out a locomotion that walks on a slope with 0.2 m/s target velocity. To verify the generalization of the proposed stability criterion, the robot is simulated to climb a slope that has 5-degree inclination.

Figure 11 graphs the instantaneous velocity and average velocity of CoG, with reference to $x$ -axis. It is clear that both of them track the reference velocity in a good performance. This indicates that the velocity feedback control at the joint level can compensate the error between the optimized velocity and the target velocity, when the robot moves on an inclined surface. The position of CoG along $z$ -axis is also plotted in Fig. 11. After the time of 4 s, the CoG position increases linearly as the CoG velocity keeps around a constant value. On the other hand, at the early process of speed tracking (between [0, 4] s), the velocity of CoG varies more intensively around the command so that the position of CoG goes up without well linearity.

Figure 11. Velocity along $x$ -coordinate and position along $z$ -coordinate of centre of gravity in the simulation of slope walking.

Figure 12 compares the results of a specific time segment between slope walking (test case 2) and level ground walking (test case 1), including the CoG velocity and the CoG position along $x$ -axis. In addition to the CoG position, the position of CoP and the boundary of stance foot are also exhibited. The segmented time period of slope walking is from 9 s to 12 s, while the time segment is between 3 s and 5 s for level ground walking. For those two segments, the robot strides in a relatively slow speed, about 0.2 m/s on the slope and 0.3 m/s on the flat terrain.

Figure 12. Comparison in velocity and position along $x$ -coordinate of centre of gravity, centre of pressure and boundaries of foot, between the simulation of slope walking and level ground walking.

From test case 1, it is already concluded that the robot CoP falls inside the support polygon for most of one step cycle if the robot walks under 0.3 m/s. The conclusion is validated again in the test case 2 with the climbing velocity of about 0.2 m/s. On the other hand, the CoP positions of test case 2 during [9, 12] s are close to the edge of the support polygon, compared to ones of test case 1 during [3, 5] s. The reason might be that going up a slope requires larger friction than moving on a level floor. The further forward placement of CoP can provide the required force.

The difference between [9, 12] s of test case 2 and [3, 5] s of test case 1 is that the robot is regulated to a constant speed for test case 2, while for test case 1 the robot walks in acceleration. Consequently, considering the time when the robot CoG exceeds the upper bounds of stance foot, it is apparent that the time for test case 2 is behind the time for test case 1. The deduction is also demonstrated by the CoG position in Fig. 12.

5.3. Comparative study

The optimization results and simulation of walking on the level ground are compared in Figs. 13, 14, and 15. Four walking cycles between 7.8 s and 9.2 s are presented. Different from the optimization, the separation of each cycle does not rely on the optimized step time, but based on the event of foot touchdown, i.e. the detection of the moment of foot switch.

Figure 13. Comparison of velocity along $x$ -coordinate and position along $z$ -coordinate of centre of gravity between optimization and simulation.

Figure 14. Comparison of rotational velocity and torque of knee joint between optimization and simulation.

Figure 15. Comparison of total energy between optimization and simulation.

Fig. 13 shows CoG velocity along $x$ -coordinate and CoG position along $z$ -coordinate, respectively. The average simulation velocity is slightly higher than the optimization one, as the velocity feedback control is performed to reduce the static error. Similar to the optimization one, the CoG trajectory of forward direction roughly complies with the dynamics of simplified LIPM. On the other hand, the variation of position along $z$ -coordinate from the simulation is more intense than one from the optimization. It is implied that the constant height feature of LIPM is hard to hold for real-biped robots.

Fig. 14 selects the left knee joint information as comparison. The rotational velocity and torque of the joint are plotted. The simulation ones basically follow the same variation trend as the optimization results, but experience the “impulse response” when the foot switch from swing phase to stance phase. It can be seen that the influence of foot touchdown is underestimated in the optimization. For real-biped robots, the foot impact usually does not happen instantaneously and the joint positions may discontinue, as opposed to the assumptions made in Section 2.1. As a consequence, the total energy after impact might largely differ from the one before impact.

The energy between optimization and simulation are compared in Fig. 15. Different from the optimization result, the simulated energy of first walking step is not completely same as the second step, but basically same as the third step. Similarly, the fourth step is a repeat of the second step. This indicates that the period of simulated walking is two steps instead of one walking step. In fact, the simulated CoG velocity at the end of the first step would not recover to its initial value, which is different from LIP trajectory. As a result, the second step does not repeat the cycle of the first one, but the third step repeat with the aid of feedback control. Finally, the simulated biped walking exhibit periodicity at intervals of two steps.

6. Conclusion

This paper studied the large-scale programming-based gait optimization methodology. A system energy-based stability criterion and objectives were proposed to achieve dynamic stable walking. The bipedal gait stability was verified using numerical optimizations and forward simulations.

The optimization analysis demonstrated that the stability criterion constraint naturally leads to a cyclic biped locomotion, though the periodicity was not previously imposed. The optimized CoG velocity was slightly lower than the target one, but this error can be easily compensated by our feedback control scheme.

The simulation studies comprised two test cases: level ground walking and slope walking. The test case 1, on the flat terrain, demonstrated that the motions of CoP and CoG are quite different for slow walking and fast walking. To track a slow constant speed, CoP and CoG fall inside the support polygon. By contrast, CoP and CoG easily exceeds the bounds of support polygon. The test case 1 also verified the generalization capability of the proposed stability criterion to both acceleration and deacceleration locomotion. On the accelerating and deaccelerating stages, the CoP and CoG would not be strictly inside the support polygon even if the robot speed is slow. All these phenomena indicate that the proposed stability criterion can accomplish the dynamic stable walking for various locomotion (like tracking a fast constant speed, accelerating or deaccelerating). On the other hand, the test case 2 verified the generalization capability of the proposed stability criterion to move on a slope surface. The behaviours of CoP in the test case 2 differ from ones in test case 1, even though two cases have a close CoG velocity. The CoP positions of test case 2 are closer to the edge of the support polygon, compared to the ones of test case 1.

The optimal and simulated gait was compared with LIP motion trajectory while several interesting facts were revealed. For a relatively slow speed, LIPM could well simulate the CoG motion. While for large speed, LIPM could not represent the dynamics and kinematics of bipeds. In the comparisons between optimization and simulation results, it is concluded that the model of impact phase requires further study since its dynamics is not completely disclosed in the optimization.

Currently, the gait optimization methodology proposed in this paper only considers the stability of periodic-biped locomotion. Though the generalization of proposed method has been verified by the simulations under different scenarios (like acceleration and deacceleration walking), the corresponding theoretical stability is still not clear. It is worth to carry out the study on the stability criterion of aperiodic-biped locomotion, using the system energy information.

For further work, the stabled gait and feedback control will be implemented to our Cosmos prototype. The robustness to system uncertainties and external disturbances will also be evaluated.