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 [Reference Erez, Tassa and Todorov40]. 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 [Reference Griffin41]. Secondly, impact is assumed that happens instantaneously. Besides, only velocities vary at impact moment, while no instantaneous changes in position [Reference Westervelt42].

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 [Reference Bijalwan, Semwal, Singh and Mandal43, Reference Raj, Semwal and Nandi44]. 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 [Reference Hereid45]. 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 [Reference Grizzle, Chevallereau, Sinnet and Ames46]. 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.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 [Reference Kumar and Dutta11]. 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 [Reference Benson47]:

(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 [Reference Benson47]:

(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 [Reference Benson47]:

(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.