Human-Like Balance Recovery Based on Numerical Model Predictive Control Strategy

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 (<inline-formula> <tex-math notation="LaTeX">$i$ </tex-math></inline-formula>) the ankle strategy, (<inline-formula> <tex-math notation="LaTeX">$ii$ </tex-math></inline-formula>) the hip-ankle strategy, (<inline-formula> <tex-math notation="LaTeX">$iii$ </tex-math></inline-formula>) 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.


I. INTRODUCTION
H UMAN 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 anteriorposterior (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  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 hipankle 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 variablestep 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 2 VOLUME 4, 2016 is proposed to implement a human-like balance strategy and autonomous switching between ankle strategy and hipankle 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 longtime 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 2. Illustration of the three-link, two-joint model. m0, m1, m2 represent the masses of foot, lower body and upper body respectively. L0 represents the length from the toe to the ankle. L1, L2 represent the lengths of lower body and upper body respectively. q0, q1, q2 represent the toe angle, the ankle angle and hip angle respectively. 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: where T is the total kinetic energy, V is the total potential energy, τ a is the ankle torque and τ 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 VOLUME 4, 2016 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 The total kinetic energy of the whole body is given by: Second, the potential energies of the foot V 0 , the lower body V 1 and the upper body V 2 are derived as follows, The total potential energy of the whole body is given by: 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 =q 0 =q 0 = 0. The vertical displacement, velocity, and acceleration of the toe are also set to zero, i.e., 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: Where, 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 X and the control input vector τ are defined as follows: Then, the dynamic equation of motion can be converted via linearization toẊ where A is the 4 × 4 state matrix and B is the 4 × 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: where x(k+1) represents the 4×1 vector of the angles and angular velocities of the ankle and hip joints at time k + 1.
x(k) represents the 4 × 1 vector of the angles and angular velocities of the ankle and hip joints at time k and τ (k) represents the 2 × 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 where Q > 0 and Q f > 0 are 4 × 4 real symmetric matrices and R > 0 is a 2 × 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 Riccatti 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 x(0), τ (0,N −1) subject to the following constraints: 1) The discrete time state space function: x(k + 1) = Ax(k) + Bτ (k).
2) For all i = 1, 2 and k = 0, 1, 2, ..., N − 1, the torques should satisfy 3) For all i = 1, ..., 4 and k = 0, ..., N , the system states satisfy:  N = 20). Then, the prediction-based optimal control sequence for the whole horizon can be obtained as The sequence of the predicted states is given by Then, the first sample of the obtained optimal control sequence τ opt (0) is applied to the system and produces the states x(1). Here, x(1) are the observed states, which can be identical or different from the predicted states x opt (1). In the next sampling time, 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 closedloop 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.
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 Then, V f = (1/2)x T Q f 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 interiorpoint (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. 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 po- 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 humanlike 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    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 , 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 similiar 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 1 , y 1 , 0, F 1 ), (x 2 , y 2 , 0, F 2 ), (x 3 , y 3 , 0, F 3 ), and (x 4 , y 4 , 0, F 4 ). 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: The CoM is calculated with the CoMs in the x-axis di-  rection 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:   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. 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, τ is the joint input torque andq is the angular velocity, t is time. 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

[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 is 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. 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.  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 1 and R 1 yielded the minimum energy consumption. Thus, we choose weight settings Q 1 and R 1 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

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.