<figure><img src="image_1.png" /></figure>

HAL
open science

# Human-like Balance Recovery Based on Numerical Model Predictive Control Strategy

Keli Shen, Ahmed Chemori, Mitsuhiro Hayashibe

► To cite this version:

Keli Shen, Ahmed Chemori, Mitsuhiro Hayashibe. Human-like Balance Recovery Based on Numerical Model Predictive Control Strategy. IEEE Access, 2020, 8, pp.92050-92060. 10.1109/ACCESS.2020.2995104. lirmm-02610501

HAL Id: lirmm-02610501

https://hal-lirmm.ccsd.cnrs.fr/lirmm-02610501

Submitted on 17 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Date of publication xxxx 00, 0000, date of current version xxxx 00, 0000.

Digital Object Identifier 10.1109/ACCESS.2017.DOI

# Human-like Balance Recovery Based on Numerical Model Predictive Control Strategy

KELI SHEN<fnref n="1" />, (Member, IEEE), AHMED CHEMORI<fnref n="2" />, (Senior Member, IEEE), and MITSUHIRO HAYASHIBE<fnref n="1" />, (Senior Member, IEEE)

Corresponding author: Keli Shen (e-mail: shen.keli.s3@dc.tohoku.ac.jp).

This work was supported by the GP-mech Program of Tohoku University, Japan and the JSPS Grant-in-Aid for Scientific Research (B), no. 18H01399.

**ABSTRACT** The purpose of this study is to implement a human-like balance recovery controller and analyze its robustness and energy consumption. Three main techniques to maintain balance can be distinguished in humans, namely (i) the ankle strategy, (ii) the hip-ankle strategy, (iii) the stepping strategy. Because we only consider quiet standing balance, then stepping is not included in our balance recovery study. Numerical model predictive control (N-MPC) is proposed to predict the best way to maintain balance against various disturbance forces. To simulate balance recovery, we build a three-link model including a foot with unilateral constraints, the lower body, and the upper body. Subsequently, we derive the dynamical equations of the model and linearize them. Based on human balance capabilities, we set bound constraints on our model, including angles and balance torques of the ankle and hip. Unilateral constraints are set on the foot, which makes our model more similar to the human quiet standing case. Finally, we implemented a simulation of the proposed ankle and hip-ankle strategy in simulation and analyzed the obtained results from kinematic and dynamic indices as well as from an energy consumption perspective. The robustness of the proposed controller was verified through the obtained simulation results. Thus, this study provides a better understanding of human quiet standing balance that could be useful for rehabilitation.

**INDEX TERMS** Balance recovery, hip-ankle strategy, Numerical Model Predictive Control, energy consumption.

## I. INTRODUCTION

HUMAN balance recovery is an important topic in the human rehabilitation field. Human balance has been extensively studied for many years. Three main balance strategies have been observed in human push-recovery experiments for different extents of perturbation in the anterior-posterior (A/P) direction. For a small pushing force exerted on the back, the human tries to move the ankles ("ankle strategy") to maintain balance while keeping the knees, hips, and neck straight. For a larger pushing force exerted on the back, the human tends to rotate the hips ("hip strategy") when the ankle movement is not large enough to maintain balance. Finally, in cases where the human cannot maintain quiet standing, they need to step forward ("stepping strategy") for balance recovery. These three balance strategies are illustrated in Fig. 1. However, in certain environments or under

special circumstances, the human may not be able to step forward. In quiet standing cases, the ankle and hip strategies are the two possible choices for human balance control. After deriving general rules from human balance experiments, many researchers tried to model this kind of human posture behavior for further motor learning and control.

In the literature, one can find many studies about human balance control, involving clinical human experiments and numerical simulations.

Vukobratovic et al. stressed the importance of artificial locomotion systems for rehabilitation equipment design [1] and proposed the concept of zero moment point (ZMP) as a part of biped locomotion stability criteria [2]. Hemami et al. used nonlinear feedback to linearize a compound inverted pendulum system for postural stability analysis [3]. Goddard et al. studied the single-support postural stabilization of a

<footnote n="1">¹Department of Robotics, Tohoku University, Sendai, Japan (e-mail: shen.keli.s3@dc.tohoku.ac.jp, hayashibe@tohoku.ac.jp)</footnote>

<footnote n="2">²LIRMM, University of Montpellier, CNRS, Montpellier, France (e-mail: Ahmed.Chemori@lirmm.fr)</footnote>
<figure><img src="image_3.png" /><figcaption>FIGURE 1. The three balance strategies with different external disturbing forces in A/P direction: (a) ankle strategy, (b) hip strategy and (c) stepping strategy.</figcaption></figure>

three-link planar model of a biped while considering the system state and input constraints in the frontal plane [4].

Gatev et al. proposed to use a feed-forward approach to assess strategies for maintaining balance during quiet standing and evaluated the effects of narrow stance width and absence of sight. A feedforward controller predicts an external input or behaves using higher-order processing instead of the simple negative feedback of a variable [5]. The postural responses to unexpected small and slow or large and faster disturbances during quiet standing are defined as 'ankle strategy' and 'hip strategy' [6], respectively. The selection of a balance strategy against a disturbance is based on the available suitable sensory data [7]. However, Kuo et al. proposed a mixed hip-ankle strategy in the anteroposterior (A/P) direction, which was implemented in biomechanical optimization models, instead of a pure ankle strategy to recover postural balance versus different disturbances. Considering the large moment of inertia of the whole body and the difficulty of separating the control of the ankle and hip joints, the objective was to minimize the neural effort, and the predictions were based on the limited effort of the ankle joint torque to recover the balance from unstable postures [8], [9]. Studies of quiet standing [10], [11], [13] stated that the strategies for postural balance recovery should be divided into two categories: anteroposterior (A/P) and mediolateral (M/L). A/P balance studies include the ankle strategy, the hip strategy, and the stepping strategy. M/L balance studies focus solely on sway motion control [12]–[15].

Pai et al. [16] predicted the center-of-mass (CoM) velocity and position of an inverted pendulum with a foot segment within the limits of a base of support (BOS) for balance recovery from forward and backward falls. Kajita et al. [17] used a preview control of the ZMP to generate a bipedal walking pattern. A ZMP tracking servo controller based on the future-reference preview control was designed to compensate for the eventual ZMP error caused by the difference

between an amplified model and a precise multibody model. Azevedo et al. [18] studied the human walking and balance strategies and proposed the 'trajectory-free nonlinear model predictive control' to simulate various walking, and also the stable standing situations. Hofmann [19] studied humanoid walking and balance control in his thesis and highlighted the importance of horizontal motion control of the CoM for balance recovery. Stephens et al. [20]–[23] studied humanoid push balance recovery during walking and quiet stances from multiple perspectives, including balance indices and control methods. However, they did not consider the balance recovery against long-time disturbances and under-actuated humanoid feet in their simulations. Liu et al. [24] proposed a balance controller based on a trajectory library used for nonlinear systems control with constraints, such as a humanoid standing balance control. Kiemel [25] demonstrated that a humanoid robot can maintain balance against larger disturbances when using the hip-ankle strategy instead of the ankle strategy from preliminary experimental evidence and used a bang-bang controller to implement the proposed hip-ankle strategy. Nenchev [26] studied the deciding between the ankle and hip strategies for balance recovery depending on acceleration data measured during the impact. Aftab et al. [27]–[30] proposed a multi-step balance recovery scheme based on linear model predictive control (LMPC) by minimizing the horizontal CoM velocity and angular velocity of a flywheel, including the use of a hip strategy and a variable-step duration to correct large perturbations on an inverted pendulum or an inverted pendulum plus a flywheel. Thus, the typical kinematics of the human hip strategy could not be observed. Choi et al. [31] studied a trajectory-free reactive stepping controller using momentum control. The proposed controller was able to make a humanoid model move passively without following a planned trajectory in the direction of a disturbance and achieved natural stepping for different pushes. Ashtiani et al. [32] developed a control scheme based on model predictive control (MPC) and capture point (CP) for balance recovery against pushes. The proposed MPC was used to guide the CP to the desired position by regulating the ZMP and the centroidal moment pivot (CMP). Penco et al. [33] developed a retargeting framework to make an iCub robot mimic a human operator's motions for maintaining whole-body balance. Yamamoto [34] studied the maximal output admissible (MOA) set of a CP feedback controller for adaptive humanoid balance with external disturbances in both the M/L and A/P directions.

After reviewing previous works, we found that they did not cover long-time disturbing forces, the robustness of MPC, the energy consumption of the hip and ankle joints, and the evolution of the ground reaction force and feet with unilateral constraints. The contributions of this paper are as follows.

1) We built a three-link simplified human model in which the foot is unilaterally constrained to remain in contact with the ground. This makes our model maintain a more human-like balance.

2) Numerical MPC with system states and control constraints
is proposed to implement a human-like balance strategy and autonomous switching between ankle strategy and hip-ankle strategies during quiet standing balance for the different disturbing forces. This model illustrates that N-MPC is similar to the behavior elicited by the human brain and nervous system from a neuroscience viewpoint. N-MPC is also endowed with a predictive aspect that enables it to predict future behavior and select a control balance strategy through the minimization of the energy consumption of the whole body. N-MPC can handle state and input constraints simultaneously. This is crucial to fulfilling realistic requirements because body limitations, including joint ranges and input torques, can be considered. N-MPC is also a robust controller able to optimize the balance strategy and to deal with different types of external disturbances, such as small and large disturbing forces and short-time and long-time disturbing forces.

3) The CoM and the center of pressure (CoP) are used as evaluation indices and constraints to maintain an upright orientation. Different disturbing forces are considered, namely small and large disturbing forces and short-time and long-time disturbing forces.

4) From the obtained simulation results, we analyzed the mixed hip-ankle strategy regarding two aspects: kinematic and dynamical indices and energy consumption. In addition, we tested the robustness of the proposed controller and verified that our model and control approach can implement a much more human-like balance behavior. Thus, the proposed controller could shed light on human motor control on the ankle and hip may become an efficient guide to understand the elderly's quiet standing balance.

The rest of this paper is organized as follows. In Section 2, the three-link model is described and its dynamic equation is derived and linearized. The proposed N-MPC is described particularly in Section 3. In Section 4, simulation settings and results are presented and discussed. Our conclusions are presented in Section 5.

## II. DYNAMIC EQUATION OF THE THREE-LINK MODEL

To implement quiet standing balance recovery, we consider the human body as a three-link simplified model comprising a unilaterally constrained foot, an ankle joint, a lower body, a hip joint, and an upper body, as illustrated in Fig. 2. The physical parameters of our model are summarized in Tables 1. Based on an existing anthropometric database [35], the total body height is 1.6 [m] and the total body mass is 66.3 [kg]. $m_0$, $m_1$, and $m_2$ represent the masses of the foot, the lower body, and the upper body respectively. $L_0$, $L_1$, and $L_2$ represent the lengths of the foot, the lower body and the upper body respectively. $q_0$, $q_1$, and $q_2$ represent toe, ankle and hip angles, respectively. It is worth noting that we ignore the body segments between the ankle joint and the hip joint, and between the hip joint and the head. This is consistent with the case of human quiet standing balance because humans maintain their knee joint angle within a certain range of disturbing forces acting on their body. However, if these

<figure><img src="image_4.png" /><figcaption>FIGURE 2. Illustration of the three-link, two-joint model. $m_0$, $m_1$, $m_2$ represent the masses of foot, lower body and upper body respectively. $L_0$ represents the length from the toe to the ankle. $L_1$, $L_2$ represent the lengths of lower body and upper body respectively. $q_0$, $q_1$, $q_2$ represent the toe angle, the ankle angle and hip angle respectively.</figcaption></figure>

<table><caption>TABLE 1. Summary of the physical parameters of the three-link, two-joint model inspired by human.</caption><thead><tr><td>Link</td><td>Mass [kg]</td><td>Length [m]</td></tr></thead><tbody><tr><td>Foot</td><td>1.3</td><td>0.3</td></tr><tr><td>Lower body</td><td>30</td><td>1.0</td></tr><tr><td>Upper body</td><td>35</td><td>0.6</td></tr></tbody></table>

disturbing forces become too large, they need to bend their knee and step forward to avoid falling down.

First, we use Lagrange formalism [36]–[38] to derive the dynamic equation of motion for this two-joint, three-link model controlled by the ankle and hip torques. The Lagrange equations are as follows:

$$ \frac{d}{dt} \left( \frac{\partial L}{\partial \dot{q}_1} \right) - \frac{\partial L}{\partial q_1} = \tau_a, \quad (1) $$

$$ \frac{d}{dt} \left( \frac{\partial L}{\partial \dot{q}_2} \right) - \frac{\partial L}{\partial q_2} = \tau_h, \quad (2) $$

$$ L = T - V. \quad (3) $$

where $T$ is the total kinetic energy, $V$ is the total potential energy, $\tau_a$ is the ankle torque and $\tau_h$ is the hip torque. In this mechanical system, we assume that the toe rotates with $q_0$ and $h$ is the vertical displacement of the toe with respect to its original point on the ground. Then, we set unilateral constraints on the foot to make the toe and the heel of the foot remain in contact with the ground. There are two degrees of freedom (DOFs) for the foot constraints: the toe angle $q_0$ and the displacement of the toe with respect to the ground $h$. The
constraint forces can be solved via the Lagrange method for $q_0$ and $h$.

Let us now express the complete dynamic equation of motion. First, the kinetic energies of the foot $T_0$, of the lower body $T_1$, and of the upper body $T_2$ are computed separately as follows

$$T_0 = \frac{1}{2} m_0 (\frac{1}{3} L_0^2 \dot{q}_0^2 + L_0 \dot{q}_0 \dot{h} + \dot{h}^2), \quad (4)$$

$$T_1 = \frac{1}{2} m_1 [L_0^2 \dot{q}_0^2 + \frac{1}{3} L_1^2 \dot{q}_1^2 + \dot{h}^2 + L_0 L_1 \dot{q}_0 \dot{q}_1 \\ \sin(q_0 - q_1) + 2\dot{h}L_0 \dot{q}_0 \cos q_0 - \dot{h}L_1 \dot{q}_1 \sin q_1], \quad (5)$$

$$T_2 = \frac{1}{2} m_2 [L_0^2 \dot{q}_0^2 + L_1^2 \dot{q}_1^2 + \frac{1}{3} L_2^2 \dot{q}_2^2 + \dot{h}^2 + L_0 L_2 \dot{q}_0 \dot{q}_2 \\ \sin(q_0 - q_2) + L_1 L_2 \dot{q}_1 \dot{q}_2 \cos(q_2 - q_1) - \dot{h} L_2 \dot{q}_2 \\ \sin q_2 + 2\dot{h} L_0 \dot{q}_0 \cos q_0 - 2\dot{h} L_1 \dot{q}_1 \sin q_1]. \quad (6)$$

The total kinetic energy of the whole body is given by:

$$T = T_0 + T_1 + T_2. \quad (7)$$

Second, the potential energies of the foot $V_0$, the lower body $V_1$ and the upper body $V_2$ are derived as follows,

$$V_0 = m_0 g (\frac{1}{2} L_0 \sin q_0 + h), \quad (8)$$

$$V_1 = m_1 g (\frac{1}{2} L_1 \cos q_1 + L_0 \sin q_0 + h), \quad (9)$$

$$V_2 = m_2 g (\frac{1}{2} L_2 \cos q_2 + L_1 \cos q_1 + L_0 \sin q_0 + h). \quad (10)$$

The total potential energy of the whole body is given by:

$$V = V_0 + V_1 + V_2. \quad (11)$$

Here, $g$ represents the gravity coefficient. Because our balance implementation is for quiet standing balance, our model foot cannot leave the ground, which means that stepping and rotation of the foot are not allowed. We set the angle, angular velocity, and acceleration of the toe equal to zero, i.e., $q_0 = \dot{q}_0 = \ddot{q}_0 = 0$. The vertical displacement, velocity, and acceleration of the toe are also set to zero, i.e., $h = \dot{h} = \ddot{h} = 0$. These unilateral constraints are required to meet with the condition of quiet standing. Based on the above equations of unilateral constraints, the complete dynamic equation of motion can be expressed as follows:

$$\begin{bmatrix} M_{11} & M_{12} \\ M_{21} & M_{22} \end{bmatrix} \begin{bmatrix} \ddot{q}_1 \\ \ddot{q}_2 \end{bmatrix} + \begin{bmatrix} C_1 \\ C_2 \end{bmatrix} = \begin{bmatrix} \tau_a \\ \tau_h \end{bmatrix}. \quad (12)$$

Where,

$$M_{11} = \frac{1}{3} m_1 L_1^2 + m_2 L_1^2,$$

$$M_{12} = M_{21} = \frac{1}{2} m_2 L_1 L_2 \cos(q_2 - q_1),$$

$$M_{22} = \frac{1}{3} m_2 L_2^2,$$

$$C_1 = -\frac{1}{2} m_2 L_1 L_2 \dot{q}_2^2 \sin(q_1 - q_2) - \frac{1}{2} m_1 g L_1 \sin q_1 \\ - m_2 g L_1 \sin q_1,$$

$$C_2 = -\frac{1}{2} m_2 L_1 L_2 \dot{q}_1^2 \sin(q_1 - q_2) - \frac{1}{2} m_2 g L_2 \sin q_2.$$

In (12), $M_{11}$ and $M_{22}$ are the effective inertia terms, $M_{12}$ and $M_{21}$ are the coupling inertia terms. $C_1$ and $C_2$ are the total of centrifugal, Coriolis, and gravity forces.

### III. NUMERICAL MODEL PREDICTIVE CONTROL (N-MPC)

In this section, we propose a linear MPC scheme for our balance recovery problem. First, we propose to linearize the dynamic model around the unstable vertical equilibrium point, which is $q_1 = q_2 = 0$ in our system. Because the disturbance of the upright standing body is considered to be small enough, it is possible to linearize the dynamic equation of motion. The state vector $\mathbf{X}$ and the control input vector $\boldsymbol{\tau}$ are defined as follows:

$$\mathbf{X} = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{bmatrix} = \begin{bmatrix} q_1 \\ q_2 \\ \dot{q}_1 \\ \dot{q}_2 \end{bmatrix}, \quad \boldsymbol{\tau} = \begin{bmatrix} \tau_1 \\ \tau_2 \end{bmatrix} = \begin{bmatrix} \tau_a \\ \tau_h \end{bmatrix}.$$

Then, the dynamic equation of motion can be converted via linearization to

$$\dot{\mathbf{X}} = \mathbf{A}\mathbf{X} + \mathbf{B}\boldsymbol{\tau}, \quad (13)$$

where $\mathbf{A}$ is the $4 \times 4$ state matrix and $\mathbf{B}$ is the $4 \times 2$ control matrix. After obtaining the state-space representation of our model, we introduce the concept of MPC and propose our N-MPC approach with boundary conditions.

MPC is also referred to as receding horizon predictive control [39]. It should be noted that we use discrete-time MPC because the proposed MPC is implemented in discrete time through a combination of discrete-time state space functions. Here, the discrete-time state space equation is given by:

$$\mathbf{x}(k+1) = \mathbf{A}\mathbf{x}(k) + \mathbf{B}\boldsymbol{\tau}(k), \quad (14)$$

where $\mathbf{x}(k+1)$ represents the $4 \times 1$ vector of the angles and angular velocities of the ankle and hip joints at time $k+1$. $\mathbf{x}(k)$ represents the $4 \times 1$ vector of the angles and angular velocities of the ankle and hip joints at time $k$ and $\boldsymbol{\tau}(k)$ represents the $2 \times 1$ vector of the ankle and hip torques at time $k$. To implement N-MPC with boundary conditions, the cost function and constraints need to be defined as follows in a finite time $N$ as follows.

The cost function is
$$J(\mathbf{x}(0), \tau_{(0,N-1)}) = \sum_{0}^{N-1} l(\mathbf{x}, k, \tau) + V_f \quad (15)$$

$$l(\mathbf{x}, k, \tau) = \frac{1}{2} (\mathbf{x}^T(k) \mathbf{Q} \mathbf{x}(k) + \tau^T(k) \mathbf{R} \tau(k)),$$

$$V_f = \frac{1}{2} \mathbf{x}^T(k) \mathbf{Q}_f \mathbf{x}(k),$$

where $Q > 0$ and $Q_f > 0$ are $4 \times 4$ real symmetric matrices and $R > 0$ is a $2 \times 2$ real symmetric matrix. $Q$ and $R$ can be used as tuning parameters to penalize the states and the control inputs. The terminal weighting $Q_f$ is defined to be equal to the solution of the algebraic Ricatti equation (ARE) [40]. This makes $V_f$ become a Lyapunov function to achieve stable MPC performance. Then, by tuning $Q$ and $R$ into suitable values, the MPC controller can be improved.

The objective is to minimize $J(\mathbf{x}(0), \tau_{(0,N-1)})$ subject to the following constraints:

1) The discrete time state space function:

$$\mathbf{x}(k+1) = \mathbf{A}\mathbf{x}(k) + \mathbf{B}\tau(k).$$

2) For all $i = 1, 2$ and $k = 0, 1, 2, ..., N-1$, the torques should satisfy

$$\tau_{min}(i) \le \tau_{i}(k) \le \tau_{max}(i),$$

where $\tau_{min}(1) = -20$ [$Nm$], $\tau_{min}(2) = -100$ [$Nm$], $\tau_{max}(1) = 20$ [$Nm$], and $\tau_{max}(2) = 100$ [$Nm$].

3) For all $i = 1, ..., 4$ and $k = 0, ..., N$, the system states satisfy:

$$x_{min}(i) \le x_{i}(k) \le x_{max}(i),$$

where $x_{min}(1) = -0.26$ [$rad$], $x_{min}(2) = -0.35$ [$rad$], $x_{min}(3) = -\infty$ [$rad/s$], $x_{min}(4) = -\infty$ [$rad/s$], $x_{max}(1) = 0.5$ [$rad$], $x_{max}(2) = 1.4$ [$rad$], $x_{max}(3) = \infty$ [$rad/s$], and $x_{max}(4) = \infty$ [$rad/s$].

4) For $k = 0, ..., N$, the CoP should satisfy:

$$CoP_{min} \le CoP_{i}(k) \le CoP_{max}(i),$$

where $CoP_{min} = -0.15$ [$m$] and $CoP_{max} = 0.15$ [$m$].

5) For $k = 0, ..., N$, the CoM should satisfy:

$$CoM_{min} \le CoM_{i}(k) \le CoM_{max}(i).$$

where $CoM_{min} = -0.15$ [$m$] and $CoM_{max} = 0.15$ [$m$].

The N-MPC problem described above can be solved as an iterative open-loop optimal control problem with a finite horizon and an observable initial state for each sampling time. For instance, let N-MPC starts at $k = 0$ with the observed initial states $\mathbf{x}(0) = \mathbf{x}$ and a prediction horizon $k = N$ (here $N = 20$). Then, the prediction-based optimal control sequence for the whole horizon can be obtained as

$$\tau_{opt} = [\tau_{opt}(0), \tau_{opt}(1), \tau_{opt}(2), \dots, \tau_{opt}(N-1)] \quad (16)$$

The sequence of the predicted states is given by

$$\mathbf{x}_{opt} = [\mathbf{x}_{opt}(1), \mathbf{x}_{opt}(2), \dots, \mathbf{x}_{opt}(N)] \quad (17)$$

Then, the first sample of the obtained optimal control sequence $\tau_{opt}(0)$ is applied to the system and produces the states $\mathbf{x}(1)$. Here, $\mathbf{x}(1)$ are the observed states, which can be identical or different from the predicted states $\mathbf{x}_{opt}(1)$. In the next sampling time, $\mathbf{x}(1)$ becomes the new initial variables for the new optimal control problem at the sampling time $k = 1$. Then, the N-MPC repeats the above described optimal process and obtains the new optimal control inputs for the current system. Afterward, the new initial state variables can be observed for the forthcoming optimal process. Thus, N-MPC is an iterative optimal control algorithm.

Stability analyses of MPC have been discussed from different perspectives in the literature [43]–[46]. Here, the stability of the proposed N-MPC is analyzed in a concise form. Here, the sufficient conditions here that ensure closed-loop asymptotic stability are obtained from a previous work [43]:

A1: state constraint satisfied in the terminal constraint set.

A2: control constraint satisfied in the terminal constraint set.

A3: the terminal constraint set is positively invariant under the control law.

A4: $V_f(x(k+1)) - V_f(x(k)) + l(\mathbf{x}, k, \tau) \le 0$, where $V_f(\cdot)$ is a local Lyapunov function.

With the constraints set for the proposed in N-MPC, conditions A1 through A3 are satisfied. Let $Q_f > 0$ satisfy the Lyapunov equation

$$A^T Q_f A + Q = 0.$$

Then, $V_f = (1/2)\mathbf{x}^T Q_f \mathbf{x}$ satisfies A4 with equality. Thus, the closed-loop system with N-MPC is asymptotically stable, which means that all the states converge to the origin. This is also verified by the obtained results, which are presented in the next section.

To implement the control scheme N-MPC scheme in the hip-ankle balance recovery simulations, we used CasADi 3.4.5 to solve the numerical optimization problems [41]. CasADi is an open-source tool that implements algorithmic differentiation (AD) on user-defined symbolic expressions. CasADi also provides standardized interfaces to a variety of numerical routines, such as the simulation, optimization, and solution of linear and nonlinear equations. The IPOPT solver, which is based on the primal-dual nonlinear interior-point (IP) method, was used in the proposed N-MPC scheme. IPOPT can solve optimization problems with boundary constraints for all variables. A multiple-shooting technique was applied for faster numerical integration and optimization. A fourth-order Runge-Kutta method (RK4) was used for the numerical integration of ordinary differential equations (ODE) [42].

The control method described in [47] is effective for three kinds of constraints, namely actuated state constraints, underactuated state constraints, and constraints on some specific composite variables. However, in this study, we need to simultaneously consider the state and input constraints, to meet the requirements of a human-like balance behavior. This is also one of the advantages of the proposed controller.
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Evolution of the ankle angle versus time for different disturbing forces: 0 [N], 20 [N], 50 [N], 80 [N], 100 [N], 120 [N].</figcaption></figure>

To the best of our knowledge, the method proposed in the aforementioned study is not straightforwardly applicable to solve our problem.

## IV. NUMERICAL SIMULATION OF THE PROPOSED RECOVERY STRATEGY

In this section, we describe our implementation of the proposed human-like hip-ankle strategy for balance recovery and the control scheme based on N-MPC in the Mujoco simulation environment [48]. The sampling period used was 0.01 [s]. The obtained simulation results are described and analyzed from kinematic, dynamic, energy consumption perspectives.

### A. KINEMATIC AND DYNAMIC ANALYSIS

In this section, we analyze the kinematic and dynamic aspects of the proposed hip-ankle strategy, such as joint angles and velocities, CoM, CoP, and control inputs. We pushed the position of the CoM of the upper body with different disturbing forces along the same direction within 0.1 [s]. The disturbing forces were set as follows: 0 [N], 20 [N], 50 [N], 80 [N], 100 [N], and 120 [N]. The model can recover balance after a perturbation within a recovery time of 12 [s]. The state weight $Q$ and the input weight $R$ are unchangeable, $I$ is a 4 × 4 identity matrix.

$$Q = 10^3 * I, \\ R = 10^{-4} * \begin{bmatrix} 1000 & 0 \\ 0 & 1 \end{bmatrix}.$$

The evolution of the ankle angle for different disturbing forces is shown in Fig. 3. Based on these results, we compare the influence of different external forces on the amplitude of the ankle angle and balance recovery time. For small disturbing forces, the ankle joint amplitude changes slightly, and the recovery time is also short. This is similar to human-like balance because, for a small pushing force, our body sways a little and maintains balance easily. For larger disturbing forces, the ankle joint changes considerably, and thus balance recovery takes a longer time. Figure 4 illustrates the relationship between the ankle joint position and velocity in a phase portrait representation. It can be seen that no matter

<figure><img src="image_6.png" /><figcaption>FIGURE 4. Evolution of the ankle phase portrait versus time for different disturbing forces: 0 [N], 20 [N], 50 [N], 80 [N], 100 [N], 120 [N].</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 5. Evolution of the hip angle versus time for different disturbing forces: 0 [N], 20 [N], 50 [N], 80 [N], 100 [N], 120 [N].</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 6. Evolution of the hip phase portrait versus time for different disturbing forces: 0 [N], 20 [N], 50 [N], 80 [N], 100 [N], 120 [N].</figcaption></figure>

how large the disturbance force is, the cycle finally converges to the origin (0, 0).

Figure. 5 shows the evolution of the hip joint angle for different disturbing forces. From the changes in amplitude of the hip joint angle, it is worth noting that for a disturbing force of 20 [N], the hip starts to react to achieve balance recovery, which indicates that the ankle strategy contains a small amount of hip rotation. Similar results have been reported in human postural balance experiments by Nashner et al. [49] and Horak et al. [50]. Moreover, for a disturbing force of 50 [N], hip rotation plays an important role, as evidenced by comparing the ankle joint amplitudes shown
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Screenshot of simulation animation for the disturbing force: 120 [N].</figcaption></figure>

in Fig. 3 and the hip joint angle amplitudes shown in Fig. 5. For disturbing forces of 100 [N] and 120 [N], the sway of the ankle is not enough to maintain balance and therefore the hip sways as well for balance control. These results indicate that the hip-ankle strategy (not a pure hip strategy) is used to maintain balance against a certain range of disturbing forces, which is similar to the published results of human movement experiments published by Runge et al. [51]. It is also worth noting that the hip joint angle amplitude is larger than the ankle joint amplitude, which is similar to the results of human experiments published by Colobert et al. [52]. The hip joint angle also converges to zero but takes a longer time than the ankle joint angle. The relationship between hip velocity and hip angle is illustrated through a phase portrait in Fig. 6. It can be seen that the cycle shown in Fig. 6 is larger than that shown in Fig. 4. This comparison also shows that for larger disturbances, our model favors using more hip movement because ankle movement is not sufficient to maintain balance. In our simulations, the ankle and hip movements observed were different from the results of Aftab et al. [27]–[30], where the upper body was not included. A screenshot of a simulation animation for a disturbing force of 120 [N] is shown in Fig. 7. We note that the deviation of the hip angle is within ‘± 0.15 [rad]’ owning to our linearization assumption. This is an important limitation of the proposed model.

The definitions of the CoP and the CoM are given as follows. A schematic diagram for calculating the location of the CoP is shown in Fig. 8. The origin of the world coordinate system is point O at the center of the foot’s bottom. The positions and ground reaction forces of four load cells under pressure are defined by (x₁, y₁, 0, F₁), (x₂, y₂, 0, F₂), (x₃, y₃, 0, F₃), and (x₄, y₄, 0, F₄). Because the model sways in the x-axis direction, the CoP and CoM in the y-axis direction are always zeros and are omitted in the study. The formula for calculating the location of the CoP in the x-axis direction is as follows:

$$CoP = \frac{F_1x_1 + F_2x_2 + F_3x_3 + F_4x_4}{F_1 + F_2 + F_3 + F_4}. \quad (18)$$

The CoM is calculated with the CoMs in the x-axis direction

<figure><img src="image_10.png" /><figcaption>FIGURE 8. The schematic diagram for calculating the CoP location. The origin of world coordinate system is point O at the foot bottom center. The positions and ground reaction forces of four load cells under pressure are defined by (x₁, y₁, 0, F₁), (x₂, y₂, 0, F₂), (x₃, y₃, 0, F₃), (x₄, y₄, 0, F₄).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 9. Evolution of the CoP versus time for different disturbing forces: 0 [N], 20 [N], 50 [N], 80 [N], 100 [N], 120 [N].</figcaption></figure>

and the masses of the foot, the lower body, and the upper body. The CoMs in each part of the model can be obtained online during the simulation using the Mujoco API and are represented by $x_f$, $x_l$, and $x_u$. In addition, $m_0$, $m_1$, and $m_2$ represent the masses of the foot, the lower body, and the upper body, respectively, as mentioned in Section 2. The formula for calculating the location of the CoM in the x-axis direction is as follows:

$$CoM = \frac{m_0x_f + m_1x_l + m_2x_u}{m_0 + m_1 + m_2}. \quad (19)$$
<figure><img src="image_12.png" /><figcaption>FIGURE 10. Evolution of the CoM versus time for different disturbing forces: 0 [N], 20 [N], 50 [N], 80 [N], 100 [N], 120 [N].</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 11. Evolution of the ankle balance recovery torque versus time for different disturbing forces: 0 [N], 20 [N], 50 [N], 80 [N], 100 [N], 120 [N].</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 12. Evolution of the hip balance recovery torque versus time for different disturbing forces: 0 [N], 20 [N], 50 [N], 80 [N], 100 [N], 120 [N].</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIGURE 13. Evolution of the ankle, hip and total energy respectively for different disturbing forces: 0 [N], 20 [N], 50 [N], 80 [N], 100 [N], 120 [N].</figcaption></figure>

The changes of the ankle and hip reaction torques for different disturbing forces are depicted in Figs. 11 and 12, respectively. Here, we consider boundary constraints of the input torques (ie. saturation) for the ankle and hip joints. This makes the behavior of our balance model more similar to the hip-ankle strategy. From Fig. 11, it can be seen that the ankle input torque starts to be saturated at -20 [N] for a disturbing force of 50 [N]. This indicates that when the pushing force becomes larger enough, the ankle joint torque produces a maximum torque of -20 [N] to maintain balance. However, the ankle input torque seems insufficient to maintain balance, and thus the hip input torque is produced to help the body maintain balance. By observing the ankle and hip input torques for pushing forces of 100 [N] and 120 [N], we find that larger hip input torques are used and the balance recovery duration becomes longer.

In this subsection, we analyzed the evolution of the angles, phase portraits, CoP, CoM, and input torques of the ankle and hip joints. As shown through our simulation results, we implemented a hip-ankle strategy with a unilaterally constrained foot and analyzed the resulting balance behavior against the different pushing forces.

## B. ENERGY CONSUMPTION VIEWPOINT

In this subsection, we analyze our implementation of the hip-ankle strategy and from a new analysis perspectives: energy consumption. The joint energy is calculated as follows, where, $\tau$ is the joint input torque and $\dot{q}$ is the angular velocity, $t$ is time.

$$W = \int \tau \dot{q} dt. \qquad (20)$$

Figure 13 clearly shows that as the disturbing force increases, the energy consumption increases to maintain balance. We compared the energy consumption of the ankle joint with that of the hip joint and note three main concluding observations. First, the energy consumption of the hip is larger than that of the ankle, for disturbing forces greater than 50 [N]. This indicates that when the pushing force is considerable, the hip joint needs to make a higher effort to maintain balance. Secondly, for pushing forces lower than

The evolution of the CoP and the CoM is represented in Figs. 9 and 10, respectively. One can observe that the CoP and CoM amplitudes become progressively higher as the disturbing force becomes greater. In addition, it is worth noting that the CoP amplitude is generally larger than the CoM amplitude. In our simulations, We make the foot model not leave the ground via unilateral constraints. We also take the CoP as a criterion to evaluate the dynamic stability of the body. The CoP remains permanently inside the footprint. If the CoP were outside the footprint, our optimization problem would become infeasible, and maintaining the balance behavior could not be guaranteed.
<figure><img src="image_16.png" /><figcaption>FIGURE 14. Evolution of the ankle energy, hip energy and total energy of both joints for different state weight *Q* settings and a constant control weight *R*: *Q*₁, *Q*₂, *Q*₃, *Q*₄, *Q*₅.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIGURE 15. Evolution of the ankle energy, hip energy and total energy of both joints for different control weight *R* settings and a constant state weight *Q*: *R*₁, *R*₂, *R*₃, *R*₄, *R*₅.</figcaption></figure>

50 [*N*], the energy consumption of the ankle is larger than that of the hip. This indicates that the ankle a higher effort for balance control than the hip joint for small disturbing forces. Thirdly, as pushing forces increase, the energy consumption of our model for balance control increases.

The state weight *Q* and input weight *R* in the cost function (15) were adjusted to find the minimum energy consumption for balance recovery. The condition for adjusting the weights is to recover balance within 12 [*s*]. First, *Q* varies over five cases while *R* is not changed as follows. *I* is a 4 × 4 identity matrix. The energy consumption of the ankle and hip joints, as well as the total energy consumed, are depicted in Fig. 14.

$$Q_1 = 10^3 * I,$$

$$Q_2 = 10^4 * I,$$

$$Q_3 = 10^5 * I,$$

$$Q_4 = 10^6 * I,$$

$$Q_5 = 10^7 * I,$$

$$R = 10^{-4} * \begin{bmatrix} 1000 & 0 \\ 0 & 1 \end{bmatrix}.$$

Then, *Q* is kept unchanged and *R* is adjusted over five cases as follows. Figure 15 shows the energy consumption of the ankle and hip joints and the total energy consumed.

$$Q = 10^3 * I,$$

$$R_1 = 10^{-4} * \begin{bmatrix} 1000 & 0 \\ 0 & 1 \end{bmatrix},$$

$$R_2 = 10^{-5} * \begin{bmatrix} 1000 & 0 \\ 0 & 1 \end{bmatrix},$$

$$R_3 = 10^{-6} * \begin{bmatrix} 1000 & 0 \\ 0 & 1 \end{bmatrix},$$

$$R_4 = 10^{-7} * \begin{bmatrix} 1000 & 0 \\ 0 & 1 \end{bmatrix},$$

$$R_5 = 10^{-8} * \begin{bmatrix} 1000 & 0 \\ 0 & 1 \end{bmatrix}.$$

In these calculations, the disturbing force was set to 100 [*N*] and its acting time was 0.1 [*s*]. By comparing the

<table><caption>TABLE 2. The ankle joint energy, hip joint energy and total energy of both joints for different disturbing forces and their application time under consideration of balance recovery within 12 [*s*].</caption><thead><tr><th>Force [N]</th><th>Force time [s]</th><th>Ankle joint energy [J]</th><th>Hip joint energy [J]</th><th>Total energy [J]</th></tr></thead><tbody><tr><td>5</td><td>0.10</td><td>0.286</td><td>0.042</td><td>0.356</td></tr><tr><td>5</td><td>8.00</td><td>0.340</td><td>0.227</td><td>0.567</td></tr><tr><td>10</td><td>0.10</td><td>0.290</td><td>0.066</td><td>0.356</td></tr><tr><td>10</td><td>7.00</td><td>1.851</td><td>2.554</td><td>4.405</td></tr><tr><td>15</td><td>0.10</td><td>0.287</td><td>0.083</td><td>0.370</td></tr><tr><td>15</td><td>6.00</td><td>2.386</td><td>2.408</td><td>4.794</td></tr><tr><td>20</td><td>0.10</td><td>0.291</td><td>0.110</td><td>0.401</td></tr><tr><td>20</td><td>0.70</td><td>4.018</td><td>8.016</td><td>12.034</td></tr><tr><td>30</td><td>0.10</td><td>0.354</td><td>0.245</td><td>0.600</td></tr><tr><td>30</td><td>0.45</td><td>4.351</td><td>8.818</td><td>13.169</td></tr><tr><td>40</td><td>0.10</td><td>0.487</td><td>0.455</td><td>0.942</td></tr><tr><td>40</td><td>0.30</td><td>3.630</td><td>7.236</td><td>10.866</td></tr><tr><td>50</td><td>0.10</td><td>0.724</td><td>0.822</td><td>1.546</td></tr><tr><td>50</td><td>0.25</td><td>4.196</td><td>8.491</td><td>12.687</td></tr><tr><td>60</td><td>0.10</td><td>1.074</td><td>1.358</td><td>2.432</td></tr><tr><td>60</td><td>0.20</td><td>4.131</td><td>8.162</td><td>12.193</td></tr><tr><td>80</td><td>0.10</td><td>1.842</td><td>3.054</td><td>4.896</td></tr><tr><td>100</td><td>0.10</td><td>3.129</td><td>5.693</td><td>8.822</td></tr><tr><td>120</td><td>0.10</td><td>4.762</td><td>9.395</td><td>14.157</td></tr></tbody></table>

energy consumption values shown in Figs. 14 and 15, which were obtained by respectively changing the settings of *Q* and *R*, we found that the weight settings *Q*₁ and *R*₁ yielded the minimum energy consumption. Thus, we choose weight settings *Q*₁ and *R*₁ as the initial settings in our balance simulations, which is representative of human behavior. When a human is pushed, their body tries to predict a way to maintain balance with low energy consumption.

We then analyzed the robustness of the proposed controller. We tested the longer time acting for different disturbances and found the maximum acting time that allows for balance recovery within 12 [*s*]. The results obtained are shown in Table. 2. For instance, when the disturbing force is 5 [*N*], the maximum acting time is 8.00 [*s*]. When the disturbing force is 10 [*N*], the maximum acting time is 7.00 [*s*]. For a disturbing force of 15 [*N*], the maximum acting time is 6.00 [*s*]. In contrast, for a disturbing force of 20 [*N*], the maximum acting time is 0.70 [*s*], etc. These results indicate that the proposed controller has good robustness. For an acting time of 8.00 [*s*] and a force of 5 [*N*], the ankle
joint energy (0.340 [J]) is larger than the hip joint energy (0.227 [J]). Therefore, the ankle strategy plays a crucial role in this case. However, for an acting time of 7 [s] and a force of 10 [N], the ankle joint energy (1.851 [J]) is lower than the hip joint energy (2.554 [J]). The hip strategy thus plays a major role when the disturbing force and acting time are increased.

### C. DISCUSSION

The stability of the closed-loop system with the proposed control strategy, which was previously analyzed, was confirmed through the results presented in this section.

In the cost function (15), the states and input torques are adjusted to desired values. In our case, the desired values are zeros. This makes the model remain upright around the vertical equilibrium point. Figures 3, 11, 5, and 12 respectively show that the angles and torques of the ankle and hip joints converge to zeros after several seconds. The limit cycles indicate that the evolution of the angles and angular velocities of the ankle and hip versus time in Figs. 4 and 6 converge to the origin. The CoM and CoP are determined by the model postures and input torques. Thus, similar results are shown in Figs. 10 and 9, respectively. In the proposed N-MPC, the ankle input torque constraints are given, and Fig. 11 shows the input saturation of the ankle. However, the hip input torque is not saturated simultaneously, which indicates that the hip makes a much higher effort to maintain balance in our model. The energy consumption of the model is determined by the balance recovery time and the angular velocities and input torques of the ankle and hip. As shown in Figs. 14 and 15, changing *Q* and *R* in the cost function (15) can affect the resulting energy consumption. The robustness of the proposed N-MPC is evident from the results shown in Table. 2, where it can be seen that the proposed model can achieve balance recovery within 12 [s] for different disturbing forces and acting times.

### V. CONCLUSION

In this paper, we proposed a new model with a unilaterally constrained foot and derived its dynamic equation of motion. Subsequently, we proposed an N-MPC scheme for our model and provided a detailed explanation of our implementation of N-MPC. We implemented the hip-ankle strategy based on the proposed model and controller in a simulated physical environment. Finally, we analyzed the obtained simulation results, which were found to be similar to those of previous human balance experiments in two perspectives: kinematic and dynamic aspects and energy consumption. This helped us gain a better understanding of the hip-ankle strategy from new perspectives. Notably, this study may also be meaningful for the control of exoskeleton devices because N-MPC is very useful for bio-mechanical optimization control. In our future works, we will consider implementing the stepping strategy. Changing the height of the center of mass is also a very interesting topic, which could help us test the robustness of the proposed controller.

## ACKNOWLEDGMENTS

This work was supported by the GP-mech Program of To-hoku University, Japan and the JSPS Grant-in-Aid for Scientific Research (B), no. 18H01399.

## REFERENCES

[1] M. Vukobratovic, A. A. Frank, and D. Juricic, "On the stability of biped locomotion," *IEEE Trans. Biomedical Engineering*, vol. bme-17, no. 1, pp. 25–36, Jan. 1970.

[2] M. Vukobratovic, B. A. Borovac, "Zero-Moment Point - thirty five years of its life," *Int. J. Humanoid Robotics*, vol. 1, no. 1, pp. 157-173, Mar. 2004.

[3] H. Hemami and P. Camana, "Nonlinear feedback in simple locomotion systems," *IEEE Trans. Automatic Control*, vol. 21, no. 6, pp. 855–860, Dec. 1976.

[4] R. Goddard, H. Hemami, and F. Weimer, "Biped side step in the frontal plane," *IEEE Trans. Automatic Control*, vol. 28, no. 2, pp. 179–187, Feb. 1983.

[5] P. Gatev, S. Thomas, T. Kepple, and M. Hallett, "Feedforward ankle strategy of balance during quiet stance in adults," *J. Physiology*, vol. 514, no. 3, pp. 915-928, Feb. 1999.

[6] L. M. Nashner, G. McCollum, "The organization of human postural movements: A formal basis and experimental synthesis," *Behav. Brain Sciences*, vol. 8, pp. 135-172, 1985.

[7] L. M. Nashner, C. L. Shupert, F.B. Horak, F. O. Black, "Organization of posture controls: an analysis of sensory and mechanical constraints," *Progr. Brain Research*, vol. 80, pp. 411-418, 1989.

[8] A. D. Kuo, F. E. Zajac, "Human standing posture: multi-joint movement strategies based on biomechanical constraints," *Progr. Brain Research*, vol. 97, pp. 349-358, 1993.

[9] A. Kuo, "An optimal control model for analyzing human postural balance," *IEEE Trans. Biomedical Engineering*, vol. 42, pp. 87–101, 1995.

[10] D. Winter, "Human balance and posture control during standing and walking," *Gait and Posture*, vol. 3, no. 4, pp. 193–214, Dec. 1995.

[11] D. A. Winter, F. Prince, J. S. Frank, C. Powell, K. F. Zabjek, "Unified theory regarding A/P and M/L balance in quiet stance," *J. Neurophysiology*, vol. 75, no. 6, pp. 2334-2343, Jun. 1996.

[12] A. Ishida, S. Imai, " Responses of the posture-control system to pseudorandom acceleration disturbances," *Medic. Biologic. Engineering and Computing*, vol. 18, no. 4, pp. 433-438, Jul. 1980.

[13] B. L. Day, M. J. Steiger, P. D. Thomson, C. D. Marsden, "Effect of vision and stance width on human body motion when standing: implications for afferent control of lateral sway," *J. Physiology*, vol. 469, pp. 479-499, Sep. 1993.

[14] R. C. Fitzpatrick, D. Burke, S. C. Gandevia, " Loop gain of reflexes controlling human standing measured with the use of postural and vestibular disturbances," *J. Neurophysiology*, vol. 76, no. 6, pp. 3994-4008, Dec. 1996.

[15] Y. Li, W. S. Levine and G. E. Loeb, "A Two-Joint Human Posture Control Model," *IEEE Trans. Neural Systems and Rehabil. Engineering*, vol. 20, no. 5, pp. 738-748, Sep. 2012.

[16] Y. C. Pai and J. Patton, "Center of mass velocity-position predictions for balance control," *J. Biomechanics*, vol. 33, no. 4, pp. 347-354, Apr. 1997.

[17] S. Kajita, F. Kanehiro, K. Kaneko, K. Fujiwara, K. Harada, K. Yokoi and H. Hirukawa, "Biped walking pattern generation by using preview control of zero-moment point," in *Proc. IEEE Int. Conf. Robot. Automat*, Taipei, Taiwan, 2003, pp. 1620-1626.

[18] C. Azevedo, P. Poignet and B. Espiau, "Artificial locomotion control: from human to robots," *Robot. Auton. Systems*, vol. 47, no. 4, pp. 203-223, Jul. 2004.

[19] A. Hofmann, "Robust execution of bipedal walking tasks from biomechanical principles," *PhD. Thesis*, Massachusetts Institute of Technology, 2006.

[20] C. G. Atkeson and B. Stephens, "Multiple balance strategies from one optimization criterion," in *Proc. IEEE-RAS Int. Conf. Humanoid Robots*, Pittsburgh, PA, USA, 2007, pp. 57-64.

[21] B. Stephens, "Humanoid push recovery," in *Proc. IEEE-RAS Int. Conf. Humanoid Robots*, Pittsburgh, PA, USA, 2007, pp. 589-595.

[22] B. Stephens and C. G. Atkeson, "Push Recovery by stepping for humanoid robots with force controlled joints," in *Proc. IEEE-RAS Int. Conf. Humanoid Robots*, Nashville, TN, USA, 2010, pp. 52-59.
[23] B. Stephens, "Push Recovery Control for Force-Controlled Humanoid Robots," *PhD. thesis*, Carnegie Mellon University, 2011.

[24] C. G. Liu and C. G. Atkeson, "Standing Balance Control Using a Trajectory Library," in *Proc. IEEE/RSJ Int. Conf. Intel. Robot. Systems*, St. Louis, MO, USA, 2009, pp. 3031-3036.

[25] S. Kiemel, "Balance maintenance of a humanoid robot using the hip-ankle strategy," *Master of Science Thesis*, Delft University of Technology, 2012.

[26] D. N. Nenchev, "Ankle and hip strategies for balance recovery of a biped subjected to an impact," *Robotica*, vol. 26, no. 5, pp. 643-653, Sep. 2008.

[27] Z. Aftab, T. Robert and P. B. Wieber, "Predicting multiple step placements for human balance recovery tasks," *J. Biomechanics*, vol. 45, no. 16, pp. 2804-2809, Sep. 2012.

[28] Z. Aftab, T. Robert and P. B. Wieber, "Simulating the effect of upper-body inertia on human balance recovery," *Comput. Method. Biomechanics. Biomedical Engineering*, vol. 15, no. S1, pp. 148-150, Sep. 2012.

[29] Z. Aftab, T. Robert and P. B. Wieber, "Ankle, hip and stepping strategies for humanoid balance recovery with a single Model Predictive Control scheme," in *Proc. IEEE/RAS Int. Conf. Humanoid Robots*, Osaka, Japan, Dec. 2012, pp. 159-164.

[30] Z. Aftab, T. Robert and P. B. Wieber, "Balance Recovery Prediction with Multiple Strategies for Standing Humans.," *PLoS One*, vol. 11, no. 3, e0151166, Mar. 2016.

[31] H. Choi, S. Lee, T. Jin and S. H. Lee "Trajectory-Free Reactive Stepping of Humanoid Robots Using Momentum Control," in *Proc. IEEE/RAS Int. Conf. Humanoid Robots*, Seoul, South Korea, Nov. 2015, pp. 1173-1178.

[32] M. S. Ashtiani, A. Y. Koma, M. S. Panahi, M. Khadiv, "Push recovery of a humanoid robot based on model predictive control and capture point," in *Proc. Int. Conf. Robot. Mechatronics*, Oct. 2016, pp. 433-438.

[33] L. Penco, B. Clement, V. Modugno, E. Mingo Hoffman, G. Nava, D. Pucci, N. G. Tsagarakis, J. B. Mouret, S. Ivaldi, "Robust Real-Time Whole-Body Motion Retargeting from Human to Humanoid," in *Proc. IEEE/RAS Int. Conf. Humanoid Robots*, Beijing, China, Nov. 2018, pp. 425-431.

[34] K. Yamamoto, "Time-variant feedback controller based on capture point and maximal output admissible set of a humanoid," *Adv. Robotics*, vol. 33, no. 18, pp. 944-955, Jun. 2019.

[35] M. Kouchi, M. Mochimaru, H. Iwasawa and S. Mitani, "Anthropometric database for Japanese Population 1997-98," *Japanese Industrial Standards Center (AIST, MIT)*, 2000.

[36] Y. G. Sun, H. Y. Qiang, J. Q. Xu, D. S. Dong, "The nonlinear dynamics and anti-sway tracking control for offshore container crane on a mobile harbor," *J. Mar. Sci. Tech*, vol. 25, no. 6, pp. 656-665, 2017.

[37] B. Siciliano, L. Sciavicco, L. Villani, G. Oriolo, *Robotics: Modelling, Planning and Control*, Springer, Verlag London, 2009.

[38] M. W. Spong, S. Hutchinson, M. Vidyasagar, *Robot Dynamics and Control*, Second Edition, Wiley, 2004.

[39] S. V. Raković and W. S. Levine, *Handbook of Model Predictive Control*, Cham Switzerland, 2018.

[40] P. Lancaster and L. Rodman, *Algeb. Riccati Equa.*, Oxford University Press, 1995.

[41] J. A. E. Andersson, J. Gillis, G. Horn, J. B. Rawlings, M. Diehl, "CasADi: a software framework for nonlinear optimization and optimal control," *Math. Program. Comput.*, vol. 11, no. 1, pp. 1-36, Jul. 2018.

[42] J. B. Rawlings, D. Q. Mayne, M. M. Diehl, *Model Predictive Control: Theory, Computation, and Design*, 2nd ed, Nob Hill Publishing, Santa Barbara, 2019.

[43] D. Q. Mayne, J. B. Rawlings, C. V. Rao, P. O. M. Scokaert, "Constrained model predictive control: Stability and optimality," *Automatica*, vol. 36, no. 6, pp. 789-814, 2000.

[44] H. Chen, F. Allgöwer, "Nonlinear model predictive control schemes with guaranteed stability," *NATO ASI on nonlinear model based process control*, vol 353. Springer, pp. 465-494, 1998.

[45] H. Chen, F. Allgöwer, "A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability," *Automatica*, vol. 14, no. 10, pp. 1205-1217, 1998.

[46] S. Valluri, V. Kapila, "Stability Analysis for Linear/Nonlinear Model Predictive Control of Constrained Processes," in *Proc. ACC*, Philadelphia, Pennsylvania. 1998, pp. 1679-1683.

[47] H. Chen, N. Sun, "Nonlinear control of underactuated systems subject to both actuated and unactuated state constraints with experimental verification," *IEEE. Trans. Ind. Elec*, DOI: 10.1109/TIE.2019.2946541

[48] E. Todorov, T. Erez, Y. Tassa, "MuJoCo: A physics engine for model-based control," in *Proc. IEEE/RSJ Int. Conf. Intel. Robot. Systems*, Vilamoura, Portugal, 2012, pp. 5026-5033.

[49] L. M. Nashner, G. McCollum, "The organization of human postural movements: a formal basis and experimental synthesis," *Behav. Brain Sci.*, vol. 8, no. 1, pp. 135-150, 1985.

[50] F. B. Horak, L. M. Nashner, "Central programming of postural movements: adaptation to altered support-surface configurations," *J. Neurophysiol.*, vol. 55, no. 6, pp. 1369-1381, 1986.

[51] C. F. Runge, C. L. Shupert, F. B. Horak, F. E. Zajac, "Ankle and hip postural strategies defined by joint torques," *G. Posture*, vol. 10, no. 2, pp. 161-170, 1999.

[52] B. Colobert, A. Crétual, P. Allard, P. Delamarche, "Force-plate based computation of the ankle and hip strategies from a double-inverted pendulum model" *Clin. Biomech*, vol. 21, no. 4, pp. 427-434, 2006.

<figure><img src="image_18.png" /></figure>

KELI SHEN received his B.S. degree from China University of Mining and Technology, Xuzhou, China, in 2015, and M.S degree from Graduate School of Natural Science and Technology, Okayama University, in 2018. He was a visiting researcher at LIRMM-CNRS from Apr. to Jul. in 2019. He is currently a third-year Ph.D. student and a member of the GP-Mech program at the Department of Robotics, Tohoku University. His research interests include numerical optimization, predictive control, humanoid robotics, human motor control, and biomechanics.

<figure><img src="image_19.png" /></figure>

AHMED CHEMORI received his M.Sc. and Ph.D. degrees, both in automatic control from Polytechnic Institute of Grenoble, France, in 2001 and 2005 respectively. During the year 2004/2005 he has been a Research and Teaching Assistant at Laboratory of Signals and Systems (LSS - Centrale Supelec) and University Paris 11. Then he joined Gipsa-Lab (Former LAG) as a CNRS post-doctoral researcher. He is currently a tenured research scientist in Automatic control and Robotics for the French National Center for Scientific Research (CNRS), at the Montpellier Laboratory of Computer Science, Robotics and Microelectronics (LIRMM). His research interests include nonlinear (adaptive, robust and predictive) control and their real-time applications in robotics.

<figure><img src="image_20.png" /></figure>

MITSUHIRO HAYASHIBE (M'04—SM'14) received his B.S. degree from Tokyo Institute of Technology, in 1999, and the M.S. and Ph.D. degrees from the Graduate School of Engineering, University of Tokyo in 2001 and 2005, respectively. He was an Assistant Professor with the Department of Medicine, Jikei University School of Medicine from 2001 to 2006. He was a Post-doctoral fellow at Institute National de Recherche en Informatique et en Automatique (INRIA), and Laboratoire d'Informatique, de Robotique, de Microelectronique de Montpellier (LIRMM) of CNRS/University of Montpellier, France, in 2007. Since 2008, he has been a tenured Research Scientist with INRIA and University of Montpellier. He has been a visiting researcher at RIKEN Brain Science Institute and TOYOTA Collaboration Center since 2012. He is now a Professor at Department of Robotics, Tohoku University, Sendai, Japan since 2017. He is recipient of the 15th Annual Delsys Prize 2017 for Innovation in Electromyography from De Luca Foundation, USA. He is a Senior Member of IEEE Engineering in Medicine and Biology Society. He serves as the Editorial Board Member for ROBOMECH Journal, SpringerOpen and co-chair of IEEE Robotics and Automation Society Technical Committee on Human Movement Understanding.

<figure><img src="image_21.png" /></figure>