2.1. Equation of motion

Herein, we utilize a rimless wheel robot as a walking robot to simplify the analysis. For the rimless wheel robot, considering the next collision leg as the swing leg results in a walking system similar to that of a biped walking robot [Reference Asano and Kawamoto19–Reference Asano and Suguro24]. The model of the rimless wheel robot is shown in Fig. 1; the robot has a torso and can actuate it using a rotary actuator at the hip joint. The equation of the motion of the rimless wheel robot is expressed as:

(1) \begin{equation} \boldsymbol{M}_{o}(\boldsymbol{q})\ddot{\boldsymbol{q}}+\boldsymbol{H}_{cg}(\boldsymbol{q},\dot{\boldsymbol{q}})=\boldsymbol{S}u_1, \end{equation}

where $\boldsymbol{q}=[\theta _{1},\,\theta _{2}]^{\textrm{T}}$ is the generalized coordinate vector, the inertial matrix $\boldsymbol{M}_{o}(\boldsymbol{q})\in \mathbb{R}$ $^{2\times 2}$ is represented as

\begin{equation*} \boldsymbol{M}_{o}(\boldsymbol{q})=\left [ \begin {array}{c@{\quad}c} \displaystyle I_{1} + (m_{1} + m_{2})l_{1}^{2} & m_{2}dl_{1}\mathrm {cos}(\theta _1 - \theta _2) \\ \displaystyle m_{2}dl_{1}\mathrm {cos}(\theta _1 - \theta _2) & m_{2}d^{2} \\ \end {array} \right ], \end{equation*}

the Coriolis, centrifugal, and gravitational force vectors $\boldsymbol{H}_{cg}(\boldsymbol{q},\dot{\boldsymbol{q}})\in \mathbb{R}$ $^{2}$ are represented as

\begin{equation*} \boldsymbol{H}_{cg}(\boldsymbol{q},\dot {\boldsymbol{q}}) =\left [ \begin {array}{c} m_{2}dl_{1}\sin (\theta _{1}-\theta _{2})\dot {\theta }_{2}^{2} -(m_{1}+m_{2})l_{1}g \sin \theta _{1} \\ -m_{2}dl_{1}\sin (\theta _{1}-\theta _{2})\dot {\theta }_{1}^{2} -m_{2}dg \sin \theta _{2} \\ \end {array} \right ], \end{equation*}

$u_1$ is the input, and the driving matrix $\boldsymbol{S} \in \mathbb{R}$ $^{2}$ is represented as

Figure 1. Schematic of the rimless wheel robot with a torso.

\begin{equation*} \boldsymbol{S}=\left [ \begin {array}{c} -1 \\ 1 \\ \end {array} \right ]. \end{equation*}

2.2. Impact equation

In the 2D simulation on MATLAB, the collision between a leg and the ground is assumed instantaneous and inelastic. After each collision, the robot’s velocity is calculated using an impact equation [Reference Grizzle, Abba and Plestan8]. Further, the constraint equation $\boldsymbol{J}_I \in \mathbb{R}^{2 \times 4}$ of the robot when the leg tip touches the ground expressed as:

(2) \begin{align} \boldsymbol{J}_{I}\dot{\boldsymbol{q}}_a&= \begin{bmatrix} l_1 \cos \theta _1 - l_1 \cos (\alpha - \theta _1) &0&1&0\\ -l_1 \sin \theta _1 - l_1 \sin (\alpha - \theta _1) &0&0&1 \end{bmatrix}\dot{\boldsymbol{q}}_a=\textbf{0}_{2\times 1}. \end{align}

The impulse vector $\boldsymbol{\lambda }_I \in \mathbb{R}^{2}$ and the velocity vector immediately after impact $\dot{\boldsymbol{q}}_a^{+} \in \mathbb{R}^4$ are expressed as:

(3) \begin{align} \boldsymbol{\lambda }_{I} &=-\boldsymbol{X}_{I}(\boldsymbol{q}_a)^{-1}\boldsymbol{J}_{I}\dot{\boldsymbol{q}}_a^{-}, \end{align}

(4) \begin{align} \dot{\boldsymbol{q}}_a^{+} &=(\boldsymbol{I}_{4 \times 4}-\boldsymbol{M}_a(\boldsymbol{q}_a)^{-1}\boldsymbol{J}_{I}^{\textrm{T}}\boldsymbol{X}_{I}(\boldsymbol{q}_a)^{-1}\boldsymbol{J}_{I})\dot{\boldsymbol{q}_a}^{-}, \end{align}

where

(5) \begin{equation} \boldsymbol{X}_{I}(\boldsymbol{q}_a)\,:\!=\boldsymbol{J}_{I}\boldsymbol{M}_a(\boldsymbol{q}_a)^{-1}\boldsymbol{J}_{I}^{\textrm{T}}, \end{equation}

$\boldsymbol{q}_a =[\theta _1,\,\theta _2,\,x_1,\,z_1]^{\textrm{T}}$ is the augmented generalized variable, $x_1$ and $z_1$ denote the position of the stance leg, $\dot{\boldsymbol{q}}_a^{-} =[\dot{\theta }_1^-,\,\dot{\theta }_2^-,\,\dot{x}_1^-,\,\dot{z}_1^-]^{\textrm{T}}$ is the velocity immediately before impact, and $\boldsymbol{M}_a(\boldsymbol{q}_a)$ is the following inertia matrix in the augmented generalized variable.

\begin{align*} \boldsymbol{M}_{a}(\boldsymbol{q}_a)&= \begin{bmatrix} M_{11} & \, M_{12} & \, M_{13} & \, M_{14} \\ & \, M_{22} & \, M_{23} &\, M_{24} \\ & \, & \, M_{33} & \, M_{34} \\ \textit{Sym.} & \, & \, & \, M_{44} \end{bmatrix}, \end{align*}

where $M_{11} = I_1+m_2 l_1^2 m_h l_1^2$ , $M_{12} = m_2 l_1 a_2 \cos (\theta _1-\theta _2)$ , $M_{13} = (m_2 + m_h) l_1 \cos (\theta _1)$ , $M_{14} = -(m_2 + m_h) l_1 \sin (\theta _1)$ , $M_{22} =m_2 a_2^2$ , $M_{23} =-m_2 a_2 \cos (\theta _2)$ , $M_{24}= -m_2 a_2 \sin (\theta _2)$ , $M_{33}=m_2 + m_h$ , $M_{34}=0$ , $M_{44}= m_2 + m_h$ .

3.1. Generating limit cycle walking generation using NMPC

A robot with a torso, such as the rimless wheel robot, can generate propulsive force for walking by maintaining the forward-tilting torso posture. While the desired trajectory of the torso can be a continuous function, herein, we assumed a constant value ( $\theta _{2d} = \textit{const}$ ). When walking on an uneven terrain, the robot is susceptible to disturbances or state jumps that occur at unpredictable and irregular times intervals; for example, leg-ground collisions, gear backlash, obstacle contact, and leg-floor slippage. In other words, if we can realize a control that converges the torso posture to the desired trajectory with minimal energy consumption under these disturbances, periodically stable and energy-efficient walking can be achieved.

Therefore, NMPC is applied to optimize the convergence of the torso to the desired trajectory and its energy consumption within short prediction intervals because disturbances occur frequently within a gait cycle, making optimization over a long interval impractical. Additionally, short prediction intervals help reduce computational cost and facilitate implementation within a short control cycle of the walking robot. To solve this optimization problem within the control cycle, a fast algorithm is required; we use the C/GMRES method [Reference Ohtsuka12], which combines continuous deformation method and GMRES to solve nonlinear optimization problems at a high speed.

3.2. NMPC algorithm

The nonlinear system in a discrete time interval is obtained by

(6) \begin{align}\boldsymbol{x}(k+1) =\boldsymbol{f}(\boldsymbol{x}(k), \boldsymbol{u}(k)), \end{align}

(7) \begin{align} \boldsymbol{x}(0)=\boldsymbol{x}_0, \end{align}

where $\boldsymbol{x}(k) \in{\mathbb R}^{n}$ is the state vector and $\boldsymbol{u}(k) \in{\mathbb R}^{m}$ is the input vector. Then, we set the following objective function:

(8) \begin{align} J \,:\!= \varphi (\boldsymbol{x}_N(k)) + \Sigma ^{N-1}_{i=0}L(\boldsymbol{x}_i(k),\boldsymbol{u}_i(k)), \end{align}

where

(9) \begin{align} \varphi (\boldsymbol{x}_N(k))&=\frac{1}{2}(\boldsymbol{x}_{N}(k)-\boldsymbol{x}_{d}(k))^{\textrm{T}}\boldsymbol{S}_f(\boldsymbol{x}_{N}(k)-\boldsymbol{x}_{d}(k)), \end{align}

(10) \begin{align} L&=\frac{1}{2}(\boldsymbol{x}_{i}(k)-\boldsymbol{x}_{d}(k))^{\textrm{T}}\boldsymbol{Q}(\boldsymbol{x}_{i}(k)-\boldsymbol{x}_{d}(k))\nonumber \\ &\quad+ \boldsymbol{u}^{\textrm{T}}_{i}(k)\boldsymbol{R}\boldsymbol{u}_{i}(k), \end{align}

(11) \begin{align} \boldsymbol{x}_{i+1}(k)&=\boldsymbol{f}(\boldsymbol{x}_{i}(k),\boldsymbol{u}_{i}(k)), \end{align}

(12) \begin{align} \boldsymbol{x}_{0}(k)&=\boldsymbol{x}(k), \end{align}

$i=0, \cdots, N-1:$ $N$ is the number of the predictive steps, $\boldsymbol{x}_d$ is the desired state vector, $\boldsymbol{S}_f$ is the terminal weight, $\boldsymbol{Q}$ is the state weight, $\boldsymbol{R}$ is the input weight, and $\boldsymbol{x}_i(k)$ and $\boldsymbol{u}_i(k)$ are the state and input vectors, respectively, in the prediction interval $i$ -step. The Hamiltonian and optimal necessary conditions are obtained as follows:

(13) \begin{align} H(\boldsymbol{x}_i(k),\boldsymbol{u}_i(k),\boldsymbol{\lambda }_{i+1}(k)) & \,:\!= L + \boldsymbol{\lambda }^{\textrm{T}}_{i+1}(k)\boldsymbol{f}(\boldsymbol{x}_{i}(k),\boldsymbol{u}_{i}(k)), \end{align}

(14) \begin{align} \boldsymbol{\lambda }_N(k) & =\frac{\partial \varphi ^{\textrm{T}}}{\partial \boldsymbol{x}}(\boldsymbol{x}_N(k)), \end{align}

(15) \begin{align} \boldsymbol{\lambda }_i(k) & =\frac{\partial H^{\textrm{T}}}{\partial \boldsymbol{x}}(\boldsymbol{x}_i(k),\boldsymbol{u}_i(k),\boldsymbol{\lambda }_{i+1}(k)), \end{align}

(16) \begin{align} \frac{\partial H}{\partial \boldsymbol{u}}(\boldsymbol{x}_i(k),\boldsymbol{u}_i(k),\boldsymbol{\lambda }_{i+1}(k)) &= 0, \end{align}

where $\boldsymbol{\lambda }_i(k) \in{\mathbb R}^{n}$ is the companion variable vector at $i$ -step. The optimal control input is set as the actual input for the system, and then, the finite horizon is moved backward by a discrete interval. This sequence is repeated in the NMPC approach. In C/GMRES, $U(k) \in{\mathbb R}^{mN}$ and $F(\boldsymbol{U}(k),\boldsymbol{x}(k),k)\in{\mathbb R}^{mN}$ are defined by

(17) \begin{align} \boldsymbol{U}(k) &\,:\!= [u_0^{\textrm{T}}(k) \cdots u_{N-1}^{\textrm{T}}(k)]^{\textrm{T}}, \end{align}

(18) \begin{align} \boldsymbol{F}(\boldsymbol{U}(k),\boldsymbol{x}(k),k) &\,:\!= \left [ \begin{array}{c} \boldsymbol{H}_u^{\textrm{T}}(\boldsymbol{x}_0(k),\boldsymbol{u}_0(k),\boldsymbol{\lambda _1}(k))\\ \vdots \\ \boldsymbol{H}_u^{\textrm{T}}(\boldsymbol{x}_{N-1}(k),\boldsymbol{u}_{N-1}(k),\boldsymbol{\lambda }_N(k)) \end{array} \right ], \end{align}

where $H_u = \frac{\partial H}{\partial u}$ . The state and companion variable vectors can be obtained using Eqs. (13)–(16). Then, $\Delta \boldsymbol{U}(k)$ is obtained when $\boldsymbol{F}=0$ is solved using the C/GMRES algorithm. Finally, $\boldsymbol{u}(k+1)=\boldsymbol{u}(k)+\Delta T\Delta \boldsymbol{u}(k)$ can be calculated.

Moreover, we introduce an integrator to suppress the oscillation of the control input as the input is unstable when there are disturbances or state jumps [Reference Onodera and Yamakita13–Reference Odachi, Oyama and Yamakita15]. To achieve this, we designed an augmented system that includes the input as a part of the state. The augmented system for the rimless wheel robot is expressed as follows:

(19) \begin{align} \boldsymbol{x}(k+1)= \left [ \begin{array}{c} y_1(k+1)\\ y_2(k+1)\\ y_3(k+1)\\ y_4(k+1)\\ u_1(k+1) \end{array} \right ]&= \left [ \begin{array}{c} y_1(k) + y_3(k)\Delta T\\ y_2(k) + y_4(k)\Delta T\\ y_3(k) + \ddot{\theta }_1 \Delta T\\ y_4(k) + \ddot{\theta }_2 \Delta T\\ u_1(k) \end{array} \right ]+ \left [ \begin{array}{c} 0 \\ 0 \\ 0 \\ 0 \\ \Delta T \end{array} \right ]\Delta u_1, \nonumber \\ &=\boldsymbol{f}(x(k)) + \boldsymbol{g}\Delta u_1, \end{align}

where

(20) \begin{align} \left [ \begin{array}{c} \ddot{\theta }_1\\ \ddot{\theta }_2\\ \end{array} \right ]= \boldsymbol{M}_{o}^{-1}(\boldsymbol{q}(k))(\boldsymbol{S}u_1(k)-\boldsymbol{H}_{cg}(\boldsymbol{q}(k),\dot{\boldsymbol{q}}(k))), \end{align}

$y_1(k)$ is the output of the angle of the stance leg, $y_2(k)$ is the output of the angle of the torso, $y_3(k)$ is the output of the angular velocity of the stance leg, $y_4(k)$ is the output of the angular velocity of the torso, $\Delta T$ is the discrete interval, and $\Delta u_1$ is the differential value of the input. The objective function is defined as follows:

(21) \begin{align} J \,:\!= \varphi (\boldsymbol{x}(N), \boldsymbol{y}_d(k)) + \Sigma ^{N-1}_{i=0}L(\boldsymbol{x}(k),\Delta \boldsymbol{u}(k)), \end{align}

where

(22) \begin{align} \varphi (\boldsymbol{x}(N), \boldsymbol{y}_d(k))&=\frac{1}{2}(\boldsymbol{y}(N)-\boldsymbol{y}_{d}(N))^{\textrm{T}} \boldsymbol{S}_f(\boldsymbol{y}(N)-\boldsymbol{y}_d(N)), \end{align}

(23) \begin{align} L(\boldsymbol{x}_k,\Delta u_1(k))&=\frac{1}{2}(\boldsymbol{y}(k)-\boldsymbol{y}_{d}(k))^{\textrm{T}}\boldsymbol{Q}(\boldsymbol{y}(k)-\boldsymbol{y}_{d}(k))\nonumber \\ &\quad+ \frac{1}{2} R \Delta u_1^2(k), \end{align}

$\boldsymbol{y}_d(k)$ is the desired output vector. Then, the Hamiltonian condition is defined as

(24) \begin{align} H(\boldsymbol{x}(k),u_1(k),\boldsymbol{\lambda }(k+1)) & \,:\!= L + \boldsymbol{\lambda }^{\textrm{T}}(k+1)(\boldsymbol{f}(\boldsymbol{x}(k))+\boldsymbol{g}\Delta u_1). \end{align}

From Eq. (24), the optimal necessary condition is obtained by

(25) \begin{align} \frac{\partial H}{\partial \Delta u}&=\Delta u_1^{\textrm{T}}R+\lambda ^{\textrm{T}}(k+1)\boldsymbol{g}. \end{align}

(26) \begin{align} &= 0. \end{align}

The optimal control input, $\Delta u_1$ , can be calculated, and $u_1$ is obtained by integrating $\Delta u_1$ .

3.3. Output-zeroing control

For comparison with the proposed method, a conventional method of output-zeroing control is introduced here [Reference Westervelt, Grizzle, Chevallereau, Choi and Morris10]. This method is input-output linearization and output-zeroing control for the torso’s attitude control. Although the actual robot is 3D, the mathematical model for the control is defined as 2D since the rimless wheel robot performs 2D motion constrained to the sagittal plane. The affine system is defined as:

(27) \begin{align} \frac{d}{dt} \begin{bmatrix} \boldsymbol{q} \\ \dot{\boldsymbol{q}} \end{bmatrix} \,:\!=f_a(\boldsymbol{q}, \dot{\boldsymbol{q}}) + g_a(\boldsymbol{q})u_1.\end{align}

\begin{align*}f_a(\boldsymbol{q},\dot{\boldsymbol{q}}) & = \begin{bmatrix} \dot{\boldsymbol{q}} \\ \boldsymbol{M}_{o}^{-1}(\boldsymbol{q})\,(-\boldsymbol{H}_{cg}(\boldsymbol{q},\dot{\boldsymbol{q}})) \end{bmatrix} \nonumber \\ g_a(\boldsymbol{q}) & = \begin{bmatrix} 0 \\ \boldsymbol{M}_{o}^{-1}(\boldsymbol{q}) \boldsymbol{S} \end{bmatrix} \end{align*}

To converge the torso angle $\theta _{2}$ to the desired torso angle $\theta _{2d}$ , the output function is defined as

(28) \begin{align} y \,:\!= h(q) = \theta _2 - \theta _{2d} = 0. \end{align}

Differentiating this function with time yields the following functions.

(29) \begin{align} \frac{dy}{dt} & = \begin{bmatrix} \frac{\partial h}{\partial q} & 0 \end{bmatrix} \begin{bmatrix} \dot{\boldsymbol{q}} \\ \boldsymbol{M}_{o}^{-1}(\boldsymbol{q})\,(-\boldsymbol{H}_{cg}(\boldsymbol{q},\dot{\boldsymbol{q}}))\end{bmatrix} \nonumber \\ &\quad +\begin{bmatrix} \frac{\partial h}{\partial q} & 0 \end{bmatrix} \begin{bmatrix} 0 \\ \boldsymbol{M}_{o}^{-1}(\boldsymbol{q}) \boldsymbol{S} \end{bmatrix} u_1 \nonumber \\ &= L_fh + L_gh, \end{align}

(30) \begin{align} \frac{d^2y}{dt^2} & = \begin{bmatrix} \frac{\partial }{\partial q}(\frac{\partial h}{\partial q} \dot{q}) & \frac{\partial h}{\partial q} \end{bmatrix} \begin{bmatrix} \dot{\boldsymbol{q}} \\ \boldsymbol{M}_{o}^{-1}(\boldsymbol{q})\,(-\boldsymbol{H}_{cg}(\boldsymbol{q},\dot{\boldsymbol{q}}))\end{bmatrix} \nonumber \\ &\quad+\frac{\partial h}{\partial q} \boldsymbol{M}_{o}^{-1}(\boldsymbol{q}) \boldsymbol{S}u_1 \nonumber \\ &= L^2_fh + L_gL_fhu_1. \end{align}

Since the relative order of the system is two, the partially linearized feedback input for input-output linearization is given by

(31) \begin{align} u_1 = -\frac{1}{L_{g}L_{f}h}(L_{f}^2h - v), \end{align}

where $L_gL_fh$ and $L_{f}^2h$ are given by

\begin{align*} L_{g}L_{f}h &= \boldsymbol{C}\boldsymbol{M}_{o}^{-1}(\boldsymbol{q})\boldsymbol{S}, \\ L_{f}^2h &= -\boldsymbol{C}\boldsymbol{M}_{o}^{-1}(\boldsymbol{q})\boldsymbol{H}_{cg}(\boldsymbol{q},\dot{\boldsymbol{q}}), \\ \boldsymbol{C} &= \begin{bmatrix} 0 & 1 \end{bmatrix}. \end{align*}

Finally, we set the new control input $v$ as

(32) \begin{align} v = -K_p(\theta _2 - \theta _{2d})-K_d\dot{\theta }_2. \end{align}

4.1. Setting for simulation

We conducted a 2D walking simulation using MATLAB to compare the performance of NMPC and output-zeroing control. This simulation assumed level ground, fully inelastic collisions, and no slippage, while disturbances were applied to each joint (the contact point of the support leg and the hip joint) as normally distributed random torques. Tables I, II, and III show physical and control parameters on MATLAB simulations. We set the initial state as $\boldsymbol{X}_0=[-0.20\,\,0.35\,\,1.50\,\,0.00]^{\textrm{T}}$ , and The desired values as $\theta _{2d} = 0.35$ rad and $\dot{\theta }_{2d} = 0$ rad/s. Furthermore, we analyzed the walking speed and the energy efficiency of walking (20 steps) using NMPC and output-zeroing control. The specific resistance (SR) metric was used to evaluate the energy efficiency, which is given by:

(33) \begin{equation} \textrm{SR}\,:\!= \frac{p}{Mgv_a}, \end{equation}

where $p$ [J/s] is the average input power, $M$ [kg] is the total weight of the robot, and $v_a$ [m/s] is the average walking speed. The average input power is defined as

\begin{align*} p\,:\!=\frac{1}{T}\int _0^T |u_1(\dot{\theta }_2-\dot{\theta }_1)|\mathrm{d}t, \end{align*}

where $T$ [s] is the total walking time.

Table I. Physical parameters in the MATLAB simulations.

Table II. Control parameters of output-zeroing control in the MATLAB simulations.

Table III. Control parameters of NMPC in the MATLAB simulations.

4.2. Parameter setting for NMPC

The key points for determining the parameters in Table III are $S_f$ , $Q$ , $R$ , $T_f$ , and $N$ . First, $T_f$ is the length of the prediction interval and is set to a short value as demonstrated in Section 3.1. Various disturbances can occur at unexpected times during the robot’s gait. When the robot deviates from the target trajectory because of such disturbances, we aimed to realize an entrainment that quickly and efficiently pulls the robot’s state back to the target trajectory using NMPC. Therefore, herein, the generation of optimal entrainment to target trajectory by NMPC is main idea. In this study, $T_f=0.05$ was set by trial and error and $N$ is the number of divisions of this prediction interval, $T_f$ . Herein, this number of divisions is set to $N=5$ in terms of computation time. If a shorter computation time is desired, we have set $N=3$ or 5 ( $N=3$ in Section 5). For $S_f$ , $Q$ , $R$ a trial-and-error approach is used to find values for which the robot state converges sufficiently to the target values. From there, the input weight $R$ is increased to evaluate the values that improve energy efficiency. Other parameters are determined with reference to previous studies [Reference Ohtsuka12].

4.3. Walking using output-zeroing control

First, we analyzed the walking using the output-zeroing control. Figures 2 and 3 show the torso angle in walking without and with disturbance, respectively. The robot achieved periodic walking using output-zeroing control without and with disturbance, the walking speed was 0.757 m/s, the SR was 0.071 without disturbance, the walking speed was 0.806 m/s, and the SR was 0.158 with disturbance. The nature of the given disturbance caused a slight increase in the walking speed. Further, we confirmed that the disturbance greatly reduced the energy efficiency of walking.

Figure 2. Torso angle with respect to time.

Figure 3. Torso angle in disturbance with respect to time.

Figure 4. Torso angle with respect to time.

Figure 5. Torso angle with respect to time in a disturbance.

4.4. Walking using NMPC

We analyzed the walking performance using NMPC under the same conditions as the output-zeroing control by including the same disturbances. Figures 4 and 5 show the torso angle in walking without and with disturbance, respectively. The robot achieved periodic walking using NMPC without and with disturbance, the walking speed was 0.772 m/s, the SR was 0.076 without disturbance, the walking speed was 0.824 m/s, and the SR was 0.084 with disturbance. The walking speed using NMPC was almost the same as observed in output-zeroing control. However, the energy efficiency was improved by nearly double. This can be attributed NMPC optimizing the energy consumption under the disturbances.

Figure 6. SR with respect to the input weight of NMPC.

Figure 7. Walking speed with respect to the input weight of NMPC.

Figure 8. SR with respect to the input weight of NMPC in disturbance.

Figure 9. Walking speed with respect to the input weight of NMPC in disturbance.

Then, we analyzed the SR values and walking speeds for the input weight of NMPC. Figures 6 and 7 show SR values and walking speeds with respect to the input weights of NMPC, respectively. We observed that SR values and walking speeds increase slightly with respect to the input weight. Further, Figs. 8 and 9 show SR values and walking speeds with respect to the input weight of NMPC in disturbance, respectively. These figures indicate that SR monotonically decreases with increasing input weight, while the walking speed is monotonically increases with the increasing input weight. The increase in walking speed is attributed to the nature of the disturbance; by increasing the input weights, the motion was designed to take advantage of the disturbance, which increased the walking speed and further optimized energy consumption. Table IV and Fig. 10 summarize the results of the above analysis. We observed no significant difference in the walking speed of output-zeroing control and NMPC under disturbances. However, there was a significant difference in their energy consumption in the presence of a disturbance. The NMPC produced energy-efficient walking even in the presence of disturbance.

Table IV. Comparison of the SR values and average walking speeds.

Figure 10. Comparison of the SR values and average walking speeds (R = 0.045).

5.1. Simulation setting

To perform the walking analysis in an environment that closely resembles an actual environment, we carried out simulations in a 3D environment using Unity. We demonstrate the practicality of walking with NMPC through an uneven surface as shown in Fig. 11. This uneven surface was created using Terrain Toolkit 2017, with the environment set to “PATH OF THE FLESH” and the maximum ground altitude was set to $0.6$ m. The simulation communicated between Unity and MATLAB and the control input was computed by the same MATLAB program, as shown in Fig. 12. Table V shows the physical parameters of the robot in Unity. These parameters are those of our previous robots [Reference Hanazawa, Nishinami and Sagara22], and the robot inertia, $I_2$ , around the center of gravity of the body since the body is rectangular. Therefore, the inertia matrix, $M_o$ , is modified as follows:

\begin{equation*} \boldsymbol{M}_{o}(\boldsymbol{q})=\left [ \begin {array}{c@{\quad}c} I_{1} + (m_{1} + m_{2})l_{1}^{2} & m_{2}dl_{1}\mathrm {cos}(\theta _1 - \theta _2) \\ m_{2}dl_{1}\mathrm {cos}(\theta _1 - \theta _2) & I_2 + m_{2}d^{2} \\ \end {array} \right ]. \end{equation*}

Tables VI and VII show the control parameters of the output-zeroing control and NMPC in Unity, respectively. The desired values are $\theta _{2d}=0.65$ rad and $\dot{\theta }_{2d}=0$ rad/s, and the initial state is $\boldsymbol{x}_0=[0.00\,\,0.65\,\,8.0\,\,0.00]^{\textrm{T}}$ .

Table V. Physical parameters of the robot in the Unity simulations.

Figure 11. Screenshot of the walking environment in Unity.

Figure 12. Overview of the connection between MATLAB and Unity.

Table VI. Control parameters of output-zeroing control in the Unity simulations.

Table VII. Control parameters of NMPC in the Unity simulations.

Figures 13 and 14 show the robot’s walking on a flat surface and on the uneven terrain for 20 s, respectively. Here, three types of uneven terrains (no. 1, 2, and 3) were created with random shapes, and a walking analysis was performed.

Figure 13. Walking on the flat surface.

Figure 14. Walking on the uneven terrain.

Figure 15. Torso angle in the walking using the output-zeroing control on the flat surface.

Figure 16. Torso angle in the walking using the output-zeroing control on the uneven terrain no. 1.

5.2. Walking by output-zeroing control in Unity

First, we analyze the walking using the output-zeroing control on the flat surface and the uneven terrain. Figures 15 and 16 show the torso angle during walking on a flat surface and on an uneven terrain no. 1, respectively. These figures show the robot performs periodic walking on the flat surface and on the uneven terrain. Table VIII and Fig. 17 show the SR values and the average speeds of walking. It was observed that walking on a flat surface is the most energy-efficient and the fastest. However, the walking speed and energy efficiency were reduced in all uneven terrains.

Table VIII. Comparison of the SR values and average walking speeds using output-zeroing control in Unity.

Figure 17. Comparison of the SR values and average walking speeds using output-zeroing control in Unity.

5.3. Walking using NMPC in Unity

Next, we analyzed walking using NMPC in Unity in the same conditions. Figures 18 and 19 show the torso angle in walking using NMPC on the flat surface and the uneven terrain no. 1, respectively. As with output-zeroing control, we observed that the robot can periodically walk on a flat surface and uneven terrain. Table IX and Fig. 20 show the SR values and average speeds for walking using NMPC. Unlike the output-zeroing control, the energy efficiency did not almost decrease while walking on uneven terrains.

Further, we analyzed the walking with varying input weights of NMPC. Figures 21 and 22 show the SR values of walking using the NMPC on the flat surface and on the uneven terrain no. 1 with respect to the input weight, respectively. Similarly, Figs. 23 and 24 show the walking speeds of the NMPC walking on the flat surface and on the uneven terrain no. 1 with respect to the input weight, respectively. We observed that, as the input weight increased, the walking performance improved. In this case, setting the input weight to a value >1 enabled energy-efficient walking. Finally, we compared the SR values and walking speeds obtained for the output-zeroing control and NMPC in Table X and Fig. 25. Comparing the output-zeroing control and NMPC on a flat surface and on the uneven terrain, it is found that NMPC outperforms the output-zeroing control in energy efficiency and walking speed. Therefore, it can be concluded that NMPC is upward compatible with output-zeroing control in this simulation case.

Table IX. Comparison of the SR values and average walking speeds using NMPC in Unity.

Figure 18. Torso angle in the walking using NMPC on a flat surface.

Figure 19. Torso angle in the walking using NMPC on the uneven terrain no. 1.

Figure 20. Comparison of the SR values and average walking speeds using NMPC in Unity.

Figure 21. SR of the walking using NMPC on the flat surface with respect to the input weight ( $R$ ).