A theoretical–algorithmic framework for the construction of balance stability boundaries of biped robots with multiple contacts with the environment is proposed and implemented on a robotic platform. Comprehensive and univocal definitions of the states of balance of a generic legged system are introduced with respect to the system's contact configuration. Theoretical models of joint-space and center of mass (COM)-space dynamics under multiple contacts, distribution of contact wrenches, and robotic system parameters are established for their integration into a nonlinear programing (NLP) problem. In the proposed approach, the balance stability capabilities of a biped robot are quantified by a partition of the state space of COM position and velocity. The boundary of such a partition provides a threshold between balanced and falling states of the biped robot with respect to a specified contact configuration. For a COM state to be outside of the stability boundary represents the sufficient condition for falling, from which a change in the system's contact is inevitable. Through the calculated stability boundaries, the effects of different contact configurations (single support (SS) and double support (DS) with different step lengths) on the robot's balance stability capabilities can be quantitatively evaluated. In addition, the balance characteristics of the experimental walking trajectories of the robot at various speeds are analyzed in relation to their respective stability boundaries. The proposed framework provides a contact-dependent balance stability criterion for a given system, which can be used to improve the design and control of walking robots.

## Introduction

When biped robots act in constrained environments, physical interactions are instrumental. The modulation of contact interactions (stepping, grasping for support, etc.) is a vital strategy for robots that are required to navigate and manipulate in the operating environment. Current research is focused on the modeling and control of available contact forces to the advantage of the robot's performance [1]. On one hand, the redundancy among motion, actuator control inputs, and contacts' distribution in time and space existing in the dynamic model of a constrained biped system could be exploited in order to plan energy-efficient motion and control through contacts [2,3]. On the other hand, the systems' underactuated behavior that can arise in different contact scenarios plays a critical role in the ability of bipeds to maintain balance [4].

Given the complexity and redundancy of contact interactions, determining whether a biped robot is balanced or falling is not a trivial task. Common ground reference point criteria for dynamic balance of a biped system, such as zero-moment point (ZMP) and its extensions [5,6], are neither necessary nor sufficient conditions for balance [7,8], since they do not take into account all the relevant physical components (torque and joint limits, effects of large perturbations, complete contact constraints, etc.) and cannot be easily extended to more general systems and tasks. More general approaches have been attempted by employing the concepts of viability theory [9,10]. In the existing definitions, a biped system's current state is viable (balanced), if and only if it lies within the system's viability kernel, defined as the set of all initial states from which there exists at least one evolution that never results in a nonviable or failed (falling) state [11]. For typical biped systems of interest, the evaluation of the viability kernel is computationally intractable due to the high dimensionality, nonlinearity, and redundancy of the systems' dynamics and multicontact interactions [11,12].

A tractable approach is to evaluate the capturability of a biped robot's state, which relates to the ability of a system to avoid ever reaching a failed state through consecutive contacts. Specifically, the N-step capturability was introduced [12] as the ability of a biped system to come to a stop after taking at most N steps. Capturability analyses of biped motion address the difficulty of high dimensionality by relying on reduced-order models and the system's centroidal dynamics. The linear inverted pendulum model (LIPM) has been studied extensively to obtain analytic solutions for capturability [13]. The necessary conditions on initial center of mass (COM) positions and velocities for balance have been derived for the LIPM in planar motion with point feet [14]. For the LIPM with extended feet, a capture point exists such that applying a contact force at the point can bring the system to a statically stable configuration [15]. The capture point is extended to the three-dimensional LIPM for step recovery as the divergent component of motion [16]. Additionally, viable capture basins have been derived using sums-of-squares optimization for centroidal momentum based planar walking models, which explore the effects of COM height, angular momentum, and impact dynamics on capturability [17].

The aforementioned capturability analyses and the associated models are applicable to discrete stepping on level ground, but are of limited use for general multicontact whole-body balancing. The presence of multiple noncoplanar contacts in the general case introduces additional difficulties in enforcing proper no-slip and no-separation (i.e., unilaterality and no-roll) constraints for each contact surface at all times. For instance, a study that addressed multicontact balancing using capturability models for level walking [18] has evolved to an online model predictive control algorithm subject to static balance constraints and linearized centroidal dynamics [19]. In particular, its method for balance preservation consists in calculating joint control at any given state such that the COM velocity, acceleration, and angular momentum can be driven to zero within the model preview while satisfying contact constraints. The use of centroidal dynamics in multicontact configurations has been successful in the application of model predictive control to humanoid robots [20,21] and for the whole-body motion planning and control of robots in a wide range of tasks, such as climbing and quadruped locomotion [22]. A recent study [23] uses a motion planning approach that includes the contact constraints in a compact form, expressed as a dynamic balance matrix, to obtain a feasible condition on the total COM wrench, from which actuator forces can be solved. Similarly, another study [24] has shown that contact constraints over uneven terrain could be maintained as long as, for a given trajectory, the system's three-dimensional COM acceleration lies within the bounds of a conical volume.

These recent studies highlight the crucial role of contacts and COM dynamics in addressing methods for balance preservation during motion generation. On the other hand, while these approaches for multicontact balance stability provide methods for the online generation of feasible trajectories and control for a desired task, they do not indicate whether a system is balanced or falling at any given instant in time. In this regard, their conditions for balance stability in multicontact situations are subject to similar limitations intrinsic to the ground reference point criteria (e.g., ZMP) used in level walking. Establishing a univocal criterion to discriminate between a balanced and a falling state, irrespective of past and future (previewed) trajectories, remains a challenge, partially due to the lack of general and comprehensive mathematical definitions related to the states of balance of constrained legged systems. Such a general balance stability criterion would lead to the evaluation of a threshold in the state space between balanced and falling that is trajectory- and control-independent.

This study investigates the balance stability capabilities of a biped robot subject to various contact interactions with the environment. A theoretical balance stability framework recently introduced by the authors [25] is extended to the case of multicontact configurations and demonstrated on a biped robot platform. Given (a) the biped robot model, (b) its COM state at a given time instant, and (c) a specified (multi-)contact configuration, the proposed approach identifies a contact-specific partition of the COM state space (position and velocity) as an indicator of whether the current state is balanced or falling. An iterative optimization-based algorithm is established to identify a partition of the COM state space that includes all possible balanced states of the given system in the specified contact configuration. This partition is defined by its boundary, called the balance stability boundary, which is calculated from the solutions of a sequence of nonlinear programing (NLP) problems.

The following contributions are made in order to implement the proposed general balance stability criterion in multicontact configurations: (1) rigorous and univocal definitions of balanced and falling states of legged systems; (2) general constrained dynamic models for biped robots in single- and multicontact configurations, in which the indeterminacy between the system's motion, control, and contact reactions is resolved; (3) a numerically efficient and kinematically consistent method for evaluating the complete COM workspace for robots in multicontact configurations; and (4) an algorithm for the construction of contact-specific balance stability boundaries. The balance stability boundary results are presented for a biped robot in various contact configurations, along with several experimental walking trajectories measured at different walking speeds.

## Definitions and Problem Statement

The proposed method identifies a partition of a robot's COM state space as an indicator of whether the current state is balanced or falling. In this regard, it is essential to establish clear and univocal definitions of balanced and falling states within the context of legged systems, based on which, the general optimization-based framework to construct the balance stability boundary is introduced.

### Contact-Dependent States of Balance.

The stability analysis of legged systems during generic tasks has been addressed in several studies with respect to the conditions of balance versus fall. The evaluation of these conditions is commonly carried out in state space, where several notions of balanced (e.g., viable, captured) and falling (e.g., nonviable, failed) states have been introduced [9,12,2528], each in its own context and perspective of quantification. Therefore, the existing notions are often nonunivocal. In general, the concepts of balance and falling should be complementary. In other words, a system that maintains balance avoids falling, while a system that is falling has lost its balance. Based on the authors' previous preliminary notions [25] and other existing problem-specific concepts [8], in this study, comprehensive and univocal definitions of the states of balance of biped robots are introduced with respect to contact configurations with the environment (Fig. 1).

Fig. 1
Fig. 1

Definition: Falling or Unbalanced. If all (controlled) trajectories starting from a given state of a biped robot lead to an inevitable change in the robot's contacts, the state is defined as falling or unbalanced with respect to the specified contact configuration.

Definition: Balanced. If there exists a (controlled) trajectory starting from a given state of a biped robot such that the robot can come to a complete stop without ever altering its contacts, the state is defined as balanced with respect to the specified contact configuration.

It can be noted that the above definitions are formulated in a complementary manner, and are valid for any general legged system with multiple contacts with the environment. From a balanced state, the biped robot can reach a static equilibrium while preserving its contacts as enabled by its balanced initial conditions and available actuation. On the other hand, the falling or unbalanced state represents a set of initial conditions for which there does not exist a control input that can preserve the biped robot's contact configuration (Fig. 1). In the unbalanced case, the inevitable change of contacts could manifest either as the impending appearance of additional contacts (e.g., step recovery, reaching for support, etc.) or as the unavoidable detachment of one or more contact points from the original contact configuration (e.g., lifting of the foot or slipping).

The inevitable change of contacts resulting from an unbalanced (falling) state can be (i) intentional (i.e., losing balance is part of the planned task of the biped) or (ii) unintentional (i.e., losing balance is not part of the original plan for the biped's task). Case (i) can be observed, for example, during the late stance phase of normal human walking, in which the human body is intentionally falling onto the next step [29] to achieve efficient dynamic walking [30]. Case (ii) is closest to the intuitive notion of an accidental fall. For instance, this case of unintended change in contacts can be observed whenever a biped system standing on two feet unexpectedly reaches for support by leaning on a wall, as a result of receiving a large unanticipated disturbance.

### Contact-Specific Balance Stability Boundary.

The robot's COM state $(r¯,r¯̇)∈ℝ6$ is used to discriminate between balanced and unbalanced (or falling) states, where $r¯=[x¯y¯z¯]T$ and $r¯̇=[x¯̇y¯̇z¯̇]T$ are the system's COM Cartesian position and velocity vectors, respectively. Based on an earlier fundamental study [25], a method is proposed to identify a partition of the COM state space that includes all possible balanced states of a given biped robot with multiple contacts. The boundary of such a partition, called the balance stability boundary$S$, is constructed from an optimization-based algorithm as the set of COM velocity extrema components (along any given direction) calculated at every COM position, that allow the robot to reach a final static equilibrium while preserving its contacts and satisfying all necessary physical constraints and system requirements. In other words, $S$ represents a threshold between balanced and unbalanced COM initial conditions and quantifies, at any given COM position $r¯$, the maximum allowable velocity perturbation along a given direction $û$ such that the robot can come to a static equilibrium, while avoiding falling. This threshold represents the contact-specific balance stability boundary $S$ of a given biped robot, irrespective of a specific trajectory, task, or motion controller.

For a biped robot, if the components of the velocity perturbations at the COM are within a certain threshold such that balance can be maintained (e.g., the system can come to a stop) on a single foot (Fig. 1; top-left), then the robot is said to be balanced with respect to the single support (SS) contact configuration. Conversely, if the components of velocity perturbations at the COM surpass a certain threshold, then the robot will not be able to stop unless the SS contact is altered, for example, by placing the nonstance foot on the ground (i.e., stepping) or a hand on the wall (Fig. 1; top-right). In this case, the initial COM state is said to be a falling (or unbalanced) state with respect to the specified SS contact configuration, and will end up in a change of contacts from SS. When evaluated for SS, the stability boundary identifies the region of all possible 0-step capturable states [12,31] for the robot. The analysis can be extended and generalized to any other specified contact configuration, including multiple contacts, such as double support (DS) with both feet in contact with the ground (Fig. 1).

## Models

The formulation and analysis of the proposed balance stability criterion are based on the complete dynamic model of the system under consideration and its mechanical interaction with the environment. A general model for constrained robotic systems with multiple contacts with the environment is established based on joint- and COM-space dynamics, including the distribution of contact forces and moments. The model includes kinematic constraints, which are used to construct the COM workspace, and system-specific joint actuation bounds.

### Joint-Space Dynamics of Biped Robots With Multiple Contacts.

Under rigid-body assumptions, the joint-space dynamics of a generic n-degree-of-freedom (DOF) biped robot under contact constraints can be formulated as follows:
$ATτ=M(q)q̈+h(q,q̇)+c(q,⋅)$
(1)

where $q=[qrTqbT]T∈ℝn+6$ is the generalized coordinates vector (dots indicate time derivatives) that includes the vector of joint angle variables $qr∈ℝn$ and the position and orientation $qb∈ℝ6$ of frame {b} attached to the system's underactuated floating base expressed with respect to the global (inertial) frame {X, Y, Z} (Fig. 2). The mass-inertia matrix is $M(q)∈ℝ(n+6)×(n+6)$, while $h(q,q̇)∈ℝn+6$ is the vector that includes Coriolis–centrifugal forces and generalized gravitational torques, A = [Inxn0nx6] is the actuated joint selection matrix, and $τ∈ℝn$ is the vector of joint actuator torques.

Fig. 2
Fig. 2
The term $c(q,⋅)$ represents the generalized torques associated with the contact interactions between the system and the environment. In general, this term describes the contact wrenches cS, due to the specified support contacts, and cF, due to falling if a change in the specified contacts becomes inevitable:
$c(q,⋅)=cS(q,⋅)+cF(q,⋅)$
(2)
When N specified support contacts are active, the effect of all contact wrenches on the robot is formulated as:
$cS(q,⋅)=−∑p=1NΦpT(q)[FpMp]$
(3)
where the index p = 1, …, N identifies the specified contacts and $Φp(q)∈ℝ6×(n+6)$ is the Jacobian matrix mapping each contact wrench into generalized torques. The contact wrenches include resultant reaction forces Fp and moments Mp exerted on the robot by the environment due to the pressure distribution at the contact interface. These reactions are expressed with respect to the global reference frame:
$Fp=[Fp,xFp,yFp,z];Mp=[Mp,xMp,yMp,z]$
(4)

Assuming planar contacts between the robot's parts and the contact surfaces under no-slip and no-tip-over conditions, the resultant reactions Fp and Mp can be applied at a selected point rp(q) of the system that is in contact with the pth contact surface $CSp$.

A fixed local reference frame {p} is associated with each contact surface, with its origin Pp belonging to the contact surface and its axes oriented according to the unit vectors $t̂p$, $n̂p$, and $ŝp$ (Fig. 2). The vectors Pp and $Θp$ define the specified position and orientation of {p} relative to the global frame {X, Y, Z}, where the elements of $Θp$ are the Euler angles of $t̂p$, $n̂p$, and $ŝp$ with respect to the global frame. In order to maintain the specified contact configuration (and thus all kinematic contact constraints are active), the biped robot must be subject to:
$C1. φp(q)=[rp(q)θp(q)]=[PpΘp] for all p=1, …,N$
(5)

where $φp(q)$ describes the position and orientation of the rigid link in contact with the environment as a function of the generalized coordinates. The vector rp(q) = [rp,x rp,y rp,z]T is the global position of a point on the link and $θp(q)$ is the link's global orientation. As a result, each matrix $Φp=∂φp/∂q$ corresponds to the Jacobian of the six linearly independent constraints associated with the planar contact p.

Due to the unilateral nature of contacts and the presence of Coulomb friction, the reaction force Fp is subject to the following physical constraints:
$C2. Fp•n̂p≥0$
(6)
$C3. ||Fp−(Fp•n̂p)n̂p||≤μp(Fp•n̂p)$
(7)
where $μp$ is the coefficient of static friction at the contact surface p = 1, …, N.
For each planar contact where $Fp•n̂p>0$, the center of pressure (COP) position with respect to the local reference frame {p} is:
$C4. [Mp•ŝpFp•n̂p0−Mp•t̂pFp•n̂p]T∈CSp$
(8)
which must be contained within its corresponding finite contact surface $CSp$ (usually approximated as a convex polygon). For a rectangular contact surface defined by $dp,t̂$ and $dp,ŝ$ as the lower (LB) and upper (UB) bounds along $t̂p$ and $ŝp$, respectively, relative to frame {p}, the above constraint becomes:
$dp,t̂LB≤Mp•ŝpFp•n̂p≤dp,t̂UB and dp,ŝLB≤−Mp•t̂pFp•n̂p≤dp,ŝUB$
(9)

In this method, the COP position constraint is imposed locally, i.e., within each corresponding contact surface, as opposed to the commonly used constraint that limits the resultant global COP (or an equivalent ZMP) position within the overall base of support (BOS), which is the convex hull formed by all contacting areas in a plane [15,16]. This approach allows handling of the general case of noncoplanar contact surfaces (for which a global COP or ZMP position cannot always be defined) and ensures that all the conditions for physically consistent COP are satisfied at each contact.

In order to prevent rotational slip, i.e., relative rotation of the contacting surfaces about $n̂p$, the torsional friction constraint should be imposed at each contact:
$C5. |Mp•n̂p|≤(Mp•n̂p)UB$
(10)

The normal component $Mp•n̂p$ of the reaction moment is the integral of incremental moments about $n̂p$ due to the infinitesimal tangential forces acting on the corresponding surface elements and subject to friction cone constraint. Therefore, the upper bound of the normal component $(Mp•n̂p)UB$ of the reaction moment before rotational slip occurs should be identified, in general terms, as a function of the friction coefficient, the dimensions of the contact surface, and the distribution of the infinitesimal normal and tangential forces.

### Distribution of Biped Contact Wrenches Using Normalized Nonlinear Programing Variables.

In addition to joint-space dynamics, a biped robotic system is subject to COM-space dynamics, which, under N contact interactions, is:
$∑p=1NFp=m(r¯̈−g)$
(11)
$∑p=1NMp+∑p=1N(rp−r¯)×Fp=H¯̇$
(12)

where g = [0, −g, 0]T is the gravitational acceleration, m is the system's total mass, $r¯̈=[x¯̈y¯̈z¯̈]T$ is the system's COM acceleration vector, and $H¯̇$ is the rate of change of angular momentum about the COM. The COM kinematic variables are expressed with respect to the global frame and are functions of joint kinematics.

Each planar contact introduces six unknowns (reaction forces and moments) into the robot's dynamic model. If a given robot's motion at a time instant has a single contact (N = 1), as in SS configuration, the resultant contact force and moment can be uniquely determined using, for instance, the above COM dynamics. In case of multiple contact supports (N > 1), the calculation of contact wrenches under rigid-body assumption constitutes a mechanically indeterminate problem in the COM space (where the system is underdetermined with more unknowns than equations), in the joint space (where the co-dependence between actuator torques and contact wrenches for a given motion may cause kinetic redundancy), and even in their combined spaces. As a result, in general multicontact configurations, the distribution of contact wrenches cannot be solved uniquely.

In order to resolve this mechanical indeterminacy, a general approach consists in formulating a nonlinear optimization problem, in which contact wrenches [2,32,33] (or their normalized form [34]), joint kinematics, and actuator torques required for a given task are solved for as the optimal variables. While these approaches are general and comprehensive for solving the contact-related redundancies, the calculations are often affected by numerical difficulties due to the nonlinearity and nonconvexity of the problem, broad scales in the NLP variables, and high sensitivity with respect to the initial solution of the NLP. Here, an efficient method for resolving the mechanical indeterminacy of contact wrenches is proposed for a two-dimensional (2D) biped robot model in its typical contact configurations in the sagittal plane (X, Y): SS with one foot in full contact with the ground and DS with both feet in full contact with the ground. While the method is illustrated below for the case of coplanar contacts (Fig. 3), it can be generalized to noncoplanar contacts as well. In two dimensions, each contact introduces three unknown reactions, while the COM-space dynamics provides a system of three equations.

Fig. 3
Fig. 3
For a given motion of the robot in the SS contact configuration (N = 1), the reactions at the stance foot are directly calculated from the above COM dynamics, as follows:
$F1,x=mx¯̈ F1,y=m(y¯̈+g) M1,z+(r1,x−x¯)F1,y−(r1,y−y¯)F1,x=H¯̇z$
(13)
while the reactions at the swing foot are null.

When the robot with a given motion is in DS (N = 2), the distribution of contact wrenches between the two feet is indeterminate with six unknowns (Fp,x, Fp,y, Mp,z, for p = 1, 2). This indeterminacy is solved in this method through the introduction of three normalized ($∈[0,1]$) time-dependent variables $α(t)$, $β(t)$, and $γ(t)$ whose values at each time step are included in the vector of NLP variables and are determined through optimization.

In DS, the COM equations for translation yield:
$F1,y+F2,y=m(y¯̈+g)$
(14)
$F1,x+F2,x=mx¯̈$
(15)
The time-varying variable $0≤α(t)≤1$, along with Eq. (14), describes the normal contact forces distribution between both feet at any given time instant, as follows:
$F1,y=αm(y¯̈+g); F2,y=(1−α)m(y¯̈+g)$
(16)
With this distribution, constraints C2 (for p = 1, 2) can be rewritten more compactly as one constraint equation:
$y¯̈+g≥0$
(17)
For the distribution of tangential reaction forces in the X direction, subject to Eq. (15), the time-varying variable $0≤β(t)≤1$ is used to model the friction force acting at every time instant at one of the contact points, point 1 for instance, within its allowable range $±μ1F1,y$:
$F1,x=(2β−1)μ1F1,y; F2,x=mx¯̈−F1,x$
(18)
Since the range of $β$ implicitly restricts F1,x within its friction cone, the constraints C3 (for p = 1, 2) are reduced to the following static friction constraint for F2,x:
$−μ2F2,y≤F2,x≤μ2F2,y$
(19)
The COM equation for rotational dynamics about the Z direction is used to model the distribution of tangential reaction moments:
$M1,z+M2,z+(r1,x−x¯)F1,y−(r1,y−y¯)F1,x+(r2,x−x¯)F2,y−(r2,y−y¯)F2,x=H¯̇z$
(20)
where (r1,x, r1,y) and (r2,x, r2,y) are the global position coordinates of the link contact points r1 = P1 and r2 = P2 that are contained in the corresponding contact surfaces $CS1$ and $CS2$, respectively. According to the above forces distribution, this can be rewritten as a function of COM kinematics, contact positions, and the variables $α$ and $β$, as follows:
$M1,z+M2,z+m(y¯̈+g)[−x¯+r1,xα+r2,x(1−α)+(2β−1)μ1α(r2,y−r1,y)]+mx¯̈(y¯−r2,y)=H¯̇z$
(21)
In the two-dimensional system, the COP constraint C4 in the tangential direction Z for contact point 1 is $d1,t̂LB≤M1,z/F1,y≤d1,t̂UB$ (Fig. 3). Based on this, the ground reaction moment M1,z acting about P1 at every time instant can be reformulated by means of the time-varying variable $0≤γ(t)≤1$, as follows:
$M1,z=[γ(d1,t̂UB−d1,t̂LB)+d1,t̂LB]F1,y$
(22)
where $γ=0$ and $γ=1$ represent the cases of minimum ($d1,t̂LB$) and maximum ($d1,t̂UB$) admissible COP position, respectively, within the contact surface $CS1$ and for F1,y > 0 (for F1,y = 0, the COP is not defined). Then, the tangential ground reaction moment acting on $CS2$ about P2 is uniquely determined from Eq. (21). Since the range of $γ$ implicitly bounds the COP of contact 1 within its contact surface $CS1$, the COP constraints C4 (for p = 1, 2) are reduced to the following inequality for M2,z:
$d2,t̂LBF2,y≤M2,z≤d2,t̂UBF2,y$
(23)

The contact wrenches that are distributed in this COM-dynamics-based approach are incorporated into the above joint-space dynamics and the associated models and constraints. The proposed strategy for contact wrenches distribution during DS reduces the number of the contact-related time-varying unknowns from six to three in a two-dimensional system. As a result, the reduced number of the NLP variables, their normalized form, and the reduced number of constraint equations can improve the numerical performance of the associated NLP problem.

### Biped Robot System Models and Constraints.

A biped robot DARwIn-OP [35] during SS and DS contact configurations is modeled in the (X, Y) sagittal plane with a 7-DOF kinematic chain (Fig. 4) where the link lengths and inertia parameters are adapted from the literature [36]. The lower-body DOFs are represented by the ankle, knee, and hip joints. The upper body (head, torso, and arms) is modeled with a rigid link attached between the hips (connected through a link of 0.001 m [37]) and with an equivalent mass.

Fig. 4
Fig. 4

The global reference frame {X, Y} is attached at the center of a stance foot (with length fl = 10.4 cm), which maintains full contact with the ground in both SS and DS contact configurations. A 3-DOF system of fictitious joints connects the global frame to the base frame {b}, and represents the biped robot's in-plane translation and rotation (Fig. 4). The fictitious joints system has zero link lengths and masses, and is unactuated. The fictitious joints are connected to the biped system through joint 1, which has zero range of motion and velocity $q1(t)=q̇1(t)=0$ at all times with its torque $τ1(t)=M1,z(t)$ as the internal reaction of the rigid foot structure. In this case, the base frame is attached to the system such that {b} = {X, Y}; hence, qb = 0 ensures that the position of the center of the stance foot is fixed at the global origin at all times.

The kinematic contact constraint C1 for the given system in planar SS and DS contacts on a level ground can be represented as follows:
${φ1(q)=[r1(q)θ1(q)]=[P10]; r2(q)•n̂1>0 for SSφ1(q)=[r1(q)θ1(q)]=[P10]; φ2(q)=[r2(q)θ2(q)]=[P20] for DS$
(24)

where P1 = [0 0]T and P2 = [sl 0]T (step length sl is the distance between the center of each foot; Fig. 3). In addition, ground penetration avoidance for points of the system other than the specified contacts is ensured. For SS, the BOS dimension is equal to fl, while for DS, it is sl + fl. Assuming same biped feet dimensions, the COP limits are $−dp,t̂LB=dp,t̂UB=0.052$ m, for p = 1, 2, where $dp,t̂UB−dp,t̂LB=fl$. Constraint C5 does not apply in the planar model.

#### Joint Limits and Constrained Center of Mass Workspace.

The ranges of motion for the system joints 2–7 are obtained from experimental measurements, and are implemented in the biped robot model as follows:
$C6. qiLB≤qi(t)≤qiUB for i=2–7$
(25)

The contact-specific COM workspace is the kinematic set of all points in space reachable by the robot's COM subject to joint limits and specified contact configuration. The COM position $r¯$ is a function of joint variables and system parameters, such as link lengths and mass distributions. In this two-dimensional case, the set of reachable COM positions $r¯=[x¯y¯]T$ forms a region in the (X, Y) sagittal plane, subject to the joint limits C6, constraint C1 for p = 1, 2 (Eq. (24)), and ground penetration avoidance for points of the system other than the specified contacts (Fig. 5).

Fig. 5
Fig. 5

For such a high-DOF system, where analytical solutions are not readily available, the contact-specific COM workspace can be numerically evaluated by generating a set of feasible system joint configurations (i.e., joint angles) and calculating the corresponding COM positions from the constrained forward kinematics, which results in a dense cloud of workspace points. For the SS configuration, which forms a single open kinematic chain, the set of feasible joint configurations is obtained by uniformly sampling points within the hyper-rectangle in joint space bounded by the lower and upper bounds of each joint angle ($qiLB$ and $qiUB$, respectively), for i = 2–7. The joint configurations that result in the ground penetration of the points (such as joints, link COMs, and swing foot) other than those of the stance foot are excluded.

For the DS configuration, the position and the orientation of the front foot are fixed with r2(q) =P2 and $θ2(q)=0$. The joint angles q2, q3, and q4 are uniformly sampled in the joint-space hyper-rectangle bounded by their respective limits, in order to calculate the X, Y position of joint 5. The remaining joint angles q5, q6, and q7 are determined through inverse kinematics such that links 5 and 6 connect joint 5 to joint 7. The solutions that result in the ground intersection with points (such as joints and link COMs) other than those of the feet are excluded from the feasible joint configurations. The resulting feasible joint configurations are used to calculate the COM workspace points in DS with forward kinematics.

#### Joint Actuation Limits.

The torque limits of joint actuators are a sensitive factor in determining the stability boundary of a legged system [26]. Compared to the commonly used constant bounds for joint torques, here a more refined formulation of the actuation limits of the biped robot is considered. DARwIn-OP is actuated by Dynamixel MX-28 servomotors, which include a DC motor and a controller. The DC motor is generally required to work in motoring mode and braking mode, and its operation in these modes is limited by the motor characteristics. In order to properly constrain the multiquadrant operation of the DC motor, both the motoring and the braking modes are considered.

The motoring mode is constrained according to the torque-speed characteristic curve, which defines its safe (continuous and intermittent) operational ranges. The slope K of the linear torque-speed curve in the servomotor datasheet [38] is 35.71 RPM/Nm (or 3.74 rad/Nm·s), resulting in:
$q̇=−Kτ+q̇NL$
(26)
where $q̇$ and $q̇NL=5.76$ rad/s are the motor speed and no-load speed at 12V, respectively. It should be noted that in this equation the torque value $q̇NL/K$ is smaller than the stall torque (2.5 Nm), since, in general, the linear relationship fails near stall torque. Considering the permissible steady-state operation, any (safely) realizable joint torque and speed satisfy the following relationships:
$C7−I. If q̇≥0 and τ≥0: q̇≤−Kτ+q̇NL$
(27)
during forward motoring mode (quadrant I, positive mechanical power). A similar linear relationship exists for reverse motoring mode (quadrant III, positive mechanical power):
$C7−III. If q̇≤0 and τ≤0: q̇≥−Kτ−q̇NL$
(28)
While the motoring mode characteristics are provided by the experimental motor data sheet, the information required to characterize the braking mode is not readily available. Hence, constant torque and joint speed limits delimited by the quasi stall torque and no-load speed, respectively, are assumed, and the braking mode operation is constrained by:
$C7−II. 0≤q̇≤q̇NL and−q̇NL/K≤τ≤0$
(29)
$C7−IV.−q̇NL≤q̇≤0 and 0≤τ≤q̇NL/K$
(30)

The resulting joint torque-velocity constraint C7 (Fig. 6) is imposed to the robot's joints 2–7.

Fig. 6
Fig. 6

## Optimization

The balance stability boundary of a biped robot in various contact configurations is constructed through an optimization-based algorithm. The method consists of iteratively solving a series of NLP problems, into which the aforementioned dynamic and constraint models in joint- and COM-spaces are implemented. For compactness of notations and without loss of generality, the two-dimensional version of the algorithm is described here for the SS and DS configurations in (X, Y) sagittal plane, which can be easily extended to general multicontact three-dimensional cases.

### Algorithm for Balance Stability Boundary Construction.

For each COM position $r¯=[x¯y¯]T$ at a given time t0 and any given directional vector $û$, the following constrained optimization problem OP is formulated over the time interval [t0, tf] to find the corresponding COM velocity extrema components:
$OP. Maximizeq(t) for SS; {q(t),α(t),β(t),γ(t)} for DS û•r¯̇(t0)$
(31)
subject to constraints C2–C7, which ensure that each solution satisfies physical and system constraints at all time instants, and constraints C1 and others (C0, C8, and C9 as described below), which identify the contact-specific stability threshold.
The COM position at a given time t0 is iteratively assigned to points of interest within the robot's contact-specific COM workspace
$C0. r¯(t0)=[xkyj]$
(32)

where k and j are the indices for the X and Y coordinates, respectively (Fig. 7). In this algorithm, a rectangular grid with a uniform spacing $Δx$ and $Δy$ is used to discretize the entire workspace area, while any general discretization strategy can be used.

Fig. 7
Fig. 7
According to the above definition of balanced state, the robotic system at current state ($r¯(t0),r¯̇(t0)$) should have the capability of reaching a final static equilibrium at the terminal time tf:
$C8. r¯(tf)=r¯home; r¯̇(tf)=r¯̈(tf)=0$
(33)

where $r¯home=[x¯homey¯home]T$ is a statically stable (i.e., $x¯∈BOS$ [39]) COM position corresponding to a physically feasible home joint configuration qhome selected for each specified contact configuration (Fig. 5). For the SS contact configuration, in which the robot's kinematic chain is highly redundant, additional terminal constraints P2(tf) = 0 and $x¯(tf)=xCOP(tf)$ are imposed in order to further enhance the final static equilibrium condition. The terminal time tf is found through numerical experiments to allow sufficient time for the robot to reach the final static equilibrium and to ensure that no time restriction for balancing is imposed [25].

Finally, maintaining the specified contact configuration is imposed through the kinematic constraint C1 along with the following kinetic constraint:
$C9. [FiMi]T=0$
(34)
for all i that are not associated with the specified set of contacts, resulting in cF = 0 in Eq. (2).
The constraints C0, C1, C8, and C9 ensure that every trajectory starting from the calculated COM velocity extrema at the given COM position reaches the final static equilibrium without entering the falling domain (or altering the specified contacts). By iteratively solving (Algorithm 1) the constrained optimization problem OP for all COM positions of interest along the grid and for any direction of interest $û$, the complete contact-specific balance stability boundary is obtained:
$S={(r¯*,r¯̇*)|(r¯*,r¯̇*)=(r¯(t0),r¯̇(t0)) from solutions to OP, for all selected û and (xk,yj)∈COM workspace}$
(35)

where a point $(r¯*,r¯̇*)$ of the balance stability boundary $S$ is the COM's extreme initial conditions resulting from the solution to each OP. It should be noted that each constrained OP is solved for the optimal time profiles in [t0, tf] of joint kinematics (angle variables, velocities, and accelerations), normalized variables ($α(t)$, $β(t)$, and $γ(t)$, for DS only, resolving the mechanical indeterminacy), and actuator torques (through Eq. (1)). From the optimal joint angle variables q*(t), the respective COM balanced state at time t0 is calculated and stored as a point of $S$ (lines 23–26 in Algorithm 1).

Table 1
Algorithm 1 Contact-Dependent construction of $S$ via grid points
Require: Specified contact configuration ▷ Cl, C9
Require: Statically stable joint pose qhome
1: Compute $r¯home$ (qhome)
2: $(x0,y0)≅(x¯home,y¯home)$$▷$ center of workspace grid
3: $S=∅$$▷$initialize balance stability boundary set
4: repeat
5:  û ← desired direction
6:  j = 0$▷$initialize index j
7:  repeat
8:   k = 0$▷$initialize index k
9:   repeat
10:    $(x¯(t0),y¯(t0))←(xk,yj)$$▷$ C0
11:    ifk = 0 then$▷$ lines 11-22: assign NLP initial solution
12:     ifj = 0 then
13:      $qis(t)←qhome$
14:     else
15:      find statically stable q close to qhome
16:      $qis(t)←q$
17:     end if
18:     $αis(t)←1−xk/sl; βis(t)←0.5; γis(t)←0.5$
19:    else
20:     $qis(t)←qk−1,j*(t)$
21:     $αis(t)←αk−1,j*(t); βis(t)←βk−1,j*(t); γis(t)←γk−1,j*(t)$
22:    end if
23:    $qk,j*(t),αk,j*(t), βk,j*(t), γk,j*(t)←solution to OP from initial NLP solution$
24:    calculate $r¯k,j(t) and r¯̇k,j(t)from qk,j*(t)$
25:    $(r¯k,j*,r¯̇k,j*)←(r¯k,j(t0),r¯̇k,j(t0))$
26:    $S=(r¯k,j*,r¯̇k,j*)∪S$$▷$store value of boundary point
27:    k = k + 1
28:    $xk=xk−1+Δx$$▷$ move to next grid point with xk > x0
29:   until OP results infeasible
30:   k = 0$▷$ reset index
31:   repeat$▷$ move to next grid point with xk < x0
32:    k = k – l
33:    $xk=xk+1−Δx$
34:    $(x¯(t0),y¯(t0))←(xk,yj)$ ▷ C0
35:    $qis(t)←qk+1,j*(t)$
36:    $αis(t)←αk+1,j*(t); βis(t)←βk+1,j*(t); γis(t)←γk+1,j*(t)$
37:    Lines 23 - 26
38:   until OP results infeasible
39:   j = j + 1$▷$ increase index j
40:   $yj=yj−1−Δy$
41:  until OP results infeasible or all desired grid rows are evaluated
42: until All desired û are evaluated
43: return$S$$▷$ balance stability boundary set
Algorithm 1 Contact-Dependent construction of $S$ via grid points
Require: Specified contact configuration ▷ Cl, C9
Require: Statically stable joint pose qhome
1: Compute $r¯home$ (qhome)
2: $(x0,y0)≅(x¯home,y¯home)$$▷$ center of workspace grid
3: $S=∅$$▷$initialize balance stability boundary set
4: repeat
5:  û ← desired direction
6:  j = 0$▷$initialize index j
7:  repeat
8:   k = 0$▷$initialize index k
9:   repeat
10:    $(x¯(t0),y¯(t0))←(xk,yj)$$▷$ C0
11:    ifk = 0 then$▷$ lines 11-22: assign NLP initial solution
12:     ifj = 0 then
13:      $qis(t)←qhome$
14:     else
15:      find statically stable q close to qhome
16:      $qis(t)←q$
17:     end if
18:     $αis(t)←1−xk/sl; βis(t)←0.5; γis(t)←0.5$
19:    else
20:     $qis(t)←qk−1,j*(t)$
21:     $αis(t)←αk−1,j*(t); βis(t)←βk−1,j*(t); γis(t)←γk−1,j*(t)$
22:    end if
23:    $qk,j*(t),αk,j*(t), βk,j*(t), γk,j*(t)←solution to OP from initial NLP solution$
24:    calculate $r¯k,j(t) and r¯̇k,j(t)from qk,j*(t)$
25:    $(r¯k,j*,r¯̇k,j*)←(r¯k,j(t0),r¯̇k,j(t0))$
26:    $S=(r¯k,j*,r¯̇k,j*)∪S$$▷$store value of boundary point
27:    k = k + 1
28:    $xk=xk−1+Δx$$▷$ move to next grid point with xk > x0
29:   until OP results infeasible
30:   k = 0$▷$ reset index
31:   repeat$▷$ move to next grid point with xk < x0
32:    k = k – l
33:    $xk=xk+1−Δx$
34:    $(x¯(t0),y¯(t0))←(xk,yj)$ ▷ C0
35:    $qis(t)←qk+1,j*(t)$
36:    $αis(t)←αk+1,j*(t); βis(t)←βk+1,j*(t); γis(t)←γk+1,j*(t)$
37:    Lines 23 - 26
38:   until OP results infeasible
39:   j = j + 1$▷$ increase index j
40:   $yj=yj−1−Δy$
41:  until OP results infeasible or all desired grid rows are evaluated
42: until All desired û are evaluated
43: return$S$$▷$ balance stability boundary set

### Nonlinear Programing.

The optimization problem OP posed at each iteration of the above algorithm is solved as a finite-dimensional NLP problem using a direct method. The collocation method is used to parameterize each joint variable trajectory with a recursive third degree B-spline curve and v control vertices [26]. The control vertices constitute the vector $x∈ℝn⋅v$ of NLP variables and their calculated values from optimization can uniquely determine joint kinematics, from which COM kinematics is obtained. The time interval [t0, tf] is discretized with h time nodes, at which the constraints and objective functions are evaluated. Additional NLP variable vectors $α,β,γ∈ℝh$ represent the time discretization of the normalized variables $α(t)$, $β(t)$, and $γ(t)$ used for the distribution of contact wrenches in the DS configuration. Therefore, for SS configuration, the complete vector of NLP variables is $x∈ℝn⋅v$, while for DS, it is $[xTαTβTγT]T∈ℝ(n⋅v+3h)$. The actuator torques are evaluated at each time step as functions of the NLP variables (Eq. (1)).

The resulting NLP problem is solved numerically using sequential quadratic programming (SQP) method. The SQP uses gradients and quasi-Newton approximations to the Hessian of the augmented Lagrangian for curvature information [40] and, starting from a given initial solution for the NLP variables, obtains search directions from a sequence of quadratic programming subproblems. The analytical forms of the gradients for the objective and constraint functions with respect to the NLP variables are specified for reliable performance of the gradient-based algorithm.

Due to the nonlinearity and nonconvexity of the problem, each OP result is highly sensitive to the choice of initial solutions for joint variables qis(t) (which are then mapped into initial solutions for control vertices x) and normalized variables $αis(t)$, $βis(t)$, and $γis(t)$ in the DS case (which are then mapped into initial solutions for $α$, $β$, and $γ$, respectively). In order to resolve this issue, a systematic iterative method for assigning initial solution to each NLP problem is established (Algorithm 1) to improve its computational performance.

## Results and Applications

As a demonstration of the proposed contact-dependent balance stability approach, the stability boundaries for DARwIn-OP biped robot in SS and DS contact configurations are evaluated in the sagittal plane. The contact-dependent boundary results are analyzed, and their uses as a stability criterion are discussed along with the experimental walking data of the robot.

### Robotic Walking Experiments.

The robot's COM trajectory data in the (X, Y) plane, as well as in the COM state space, were obtained from experimental walking trials. Using the controller walk_tuner [41] (ROBOTIS Inc.), DARwIn-OP performed five walking trials at increasing speeds (Fig. 8). Each walking trial consists of a minimum of five consecutive periodic steps, from which the data of one full step cycle—from SS to DS phase—is extracted for analyses. Transitions between different phases of the step cycle were detected by embedded force sensing resistors located at the four corners of each rectangular foot sole. For each walking trial, the measured joint angle data were mapped into COM kinematics using a custom forward kinematics model of the robot. Five step cycles were analyzed, one at each speed, to obtain the COM (X, Y) paths with step lengths sli, for i = 1–5, and the corresponding walking parameters. At each step, the origin of the global reference frame is positioned at the center of the stance foot in the SS phase, which is the rear foot in DS in these examples.

Fig. 8
Fig. 8

### Center of Mass Workspace Construction and Discretization.

For six contact configurations (SS and DS with sli, i = 1–5), the corresponding constrained COM workspaces are evaluated in the (X, Y) sagittal plane (Fig. 9). The workspace area is largest when the system is in SS, since in this case, the swing foot is only subject to the inequality constraint $r2,y≥0$. For the DS configurations, the swing foot is constrained to the fixed position P2 = [sli 0]T for i = 1–5. As the step length sli increases, the DS workspace becomes smaller and is shifted forward, which is due to the reduced ranges of reachable joint configurations. Since every joint configuration that satisfies the DS kinematic constraints also satisfies those of SS, all the DS workspace areas are entirely included in the SS workspace area.

Fig. 9
Fig. 9

For the construction of the balance stability boundary $S$, the discretized grid points are evaluated for each workspace area using a uniform spacing $Δx=Δy=1$ cm (Fig. 9). The COM velocity extrema of the robot are evaluated at each grid point, which serves as an initial COM position for each problem OP. The resulting initial conditions are stored as the points of $S$ for the specified contact configuration. Then, the grid points nearest to the COM paths in the (X, Y) plane for the five walking motions (Fig. 9) are identified to be approximately at $y¯=0.2081$ m, and the corresponding ranges of X within the workspace limits are [−0.135, 0.130] for SS and [−0.068, 0.085], [−0.056, 0.086], [−0.035, 0.089], [−0.013, 0.092], and [0.006, 0.097], all in meters, for each DS step length. The velocity extrema at these COM positions of interest are retrieved from the entire set $S$.

### Contact-Specific Balance Stability Boundaries.

The balance stability boundaries for SS and DS contact configurations are constructed for the directions of interest $û=[1 0]T$ and $û=[−1 0]T$ to evaluate the stability characteristics of the robot in the sagittal plane against positive and negative perturbations along the X-axis. The results corresponding to the selected grid points of interest are projected onto the X-state space (Figs. 10 and 11). The coefficients of static friction $μ1=μ2=1.0$ and a time duration (tft0) of 3 s are used in the examples. Due to the high kinematic redundancy in the open loop chain, the SS stability boundary is evaluated in two stages, in which the solutions from the proposed algorithm that are close to the workspace limits are reprocessed for improved initial solutions.

Fig. 10
Fig. 10
Fig. 11
Fig. 11

The SS stability boundary extends to the entire workspace dimension and well beyond the size of the BOS (equal to fl) (Fig. 10). This illustrates that, for the given biped robot, there exists a significant region of the COM X-state space in which the states are balanced with respect to SS even when they are not statically stable. Furthermore, at the rear and front edges of the BOS, the robot's COM can sustain negative and positive X-velocity extrema, respectively, due to the available angular momentum (or multisegmental) regulation, as opposed to 1-DOF legged systems, for which $H¯̇=0$ [25,42].

The DS stability boundaries cover the entire ranges of their respective workspace limits (Fig. 11) and are highly dependent on the step lengths, while their shapes are consistent across various step lengths and are comparable with that for SS of the biped robot (Fig. 10) and simple lower-DOF legged systems [25]. In particular, for larger BOS, the extension of $S$ along the X-axis decreases as a result of the reduced workspace size. Only for a very small step length (sl1), the boundary $S$ includes both statically ($x¯∈BOS$) and dynamically ($x¯∉BOS$) balanced COM states ($x¯,x¯̇$). On the other hand, for $sl≥sl2$, all the balanced states are statically balanced. Depending on their velocity perturbations, statically stable COM positions may also result in falling states. These results indicate that, even for DS configurations, static stability (i.e., $x¯∈BOS$) of a COM state ($x¯,x¯̇$) is neither a necessary nor a sufficient condition for it to be balanced.

Analyses of the relationships among the stability boundaries for different contact configurations provide meaningful insights (Fig. 11; bottom-right). The region of the SS stability boundary is in general larger than those for DS, which is consistent with the respective workspace dimensions. This is because the robot in SS is less constrained than in DS, and thus can more freely employ the angular momentum of the nonstance leg in balancing. In DS, the allowable velocity extrema in the positive X direction increase as the BOS dimension increases, while they are not significantly affected by the BOS in the negative X direction. This asymmetry in the DS balancing capabilities in the anterior–posterior directions may be due to the kinematic and kinetic asymmetry in the biped robot's sagittal-plane model.

The biped robot in SS can recover from larger COM velocity perturbations in the negative X direction than that in DS, while simultaneously ensuring unaltered contact with the environment. For instance, when the robot is in DS and its COM velocity is perturbed in the negative direction such that its COM state is outside of the DS stability boundary but inside of that for SS, the system can still maintain balance using only foot/feet contacts by lifting the front stance foot. Conversely, during DS, the maximum allowable COM velocity perturbation along the positive X direction that the system can sustain without altering its contacts is smaller than in SS for relatively small BOS dimensions (sl1 and sl2), while it increases and surpasses the SS threshold as the BOS dimension increases ($sl≥sl3$). In other words, the maximum allowable positive velocity perturbation of the DS configuration can be either within or beyond that of the SS configuration depending on the step size. This implies that, for small step lengths ($sl≤sl3$), it is always possible for the robot to instantaneously transition from DS to SS (on the rear stance foot in this example) without ever resulting into an unbalanced (or falling) state. On the other hand, this is not always possible for larger step lengths, for which the transition from a balanced state in DS to a balanced state in SS might not be instantaneously achievable if the current DS balance state is outside of the SS stability boundary.

### Balance Stability Analysis of Robotic Walking Motions.

The balance stability of the biped robot during various walking motions can be characterized using their respective stability boundaries (Fig. 12). As the step size (and thus the walking speed) increases, the range of the robot's walking trajectory in the COM state space increases, but never reaches the stability boundaries or the respective SS and DS BOS limits. In all five walking motions, the transitions from the SS to DS, and then to the next SS phase occur well within the region inside the DS stability boundary, and the COM X-positions are maintained within their respective BOS during SS and DS phases. This illustrates that the DARwIn-OP robot's default gait generator results in walking motions that are not only balanced with respect to the SS (i.e., 0-step-captured) and DS contacts at all times, but are also statically stable [37] during the entire gait cycle with its foot/feet in full contact (i.e., no relative motion such as tip-over) with the ground, which is common for many (but not all) current biped robots. In contrast, normal human walking motions are highly dynamic [30,37] and mostly consist of states that are not 0-step capturable [26,43], with a foot contact surface that typically rocks from heel to toe such that advancing forward from an SS to the next DS configuration is inevitable [29].

Fig. 12
Fig. 12

## Concluding Remarks: Broader Implications in Balance Control

The proposed approach to contact-dependent balance stability analysis of a biped robot has been demonstrated for two types of contact configurations (SS and DS) in the sagittal plane, which are most relevant for bipedal walking. However, the framework could be extended to any multicontact configuration for stability analysis of three-dimensional generic tasks, as long as reasonably accurate reaction distributions are provided (either from prediction, as in this study, or from direct measurements). All possible contact-specific stability boundaries can be precomputed offline to become a property of a robot, and can be implemented as a control module in real-time applications (with low computational cost) as a readily available criterion for online fall prediction. In this regard, the proposed state-space partition serves as a “map” of balanced states for control of biped robots. In practice, the robot's actual capability of balancing depends not only on the calculated stability boundaries but also on the specific controller that is implemented in the system, which usually further limits its balancing capability. This map could be a useful reference for the design of balance controllers, whose performance region can be analyzed in the state space with respect to the stability boundaries evaluated beforehand. For instance, in robotic gait applications, this method could provide guidelines for generating less conservative walking phases. By exploiting the effects of limbs' passive dynamics, the robotic gait motion could be controlled such that it extends into the region of statically unstable states or even into the falling states, hence mimicking human locomotion.

## Acknowledgment

The authors would like to thank Wenjie Chen, Sunglyoung Kim, and Yunjiageng Chen for their help in the experiments and data processing. A preliminary portion of this material was presented at the 2017 ASME IDETC/CIE 41st Mechanisms and Robotics Conference (Cleveland, OH).

## Funding Data

• National Science Foundation, Directorate for Computer and Information Science and Engineering (Grant No. IIS-1427193).

• National Science Foundation, Directorate for Engineering (Grant No. CMMI-1436636).

## References

1.
Righetti
,
L.
,
Buchli
,
J.
,
Mistry
,
M.
,
Kalakrishnan
,
M.
, and
Schaal
,
S.
,
2013
, “
Optimal Distribution of Contact Forces With Inverse-Dynamics Control
,”
Int. J. Rob. Res.
,
32
(
3
), pp.
280
298
.
2.
Mummolo
,
C.
,
Mangialardi
,
L.
, and., and
Kim
,
J. H.
,
2015
, “Concurrent Contact Planning and Trajectory Optimization in One Step Walking Motion,”
ASME
Paper No. DETC2015-47745.
3.
Xi
,
W.
,
Yesilevskiy
,
Y.
, and
Remy
,
C. D.
,
2016
, “
Selecting Gaits for Economical Locomotion of Legged Robots
,”
Int. J. Rob. Res.
,
35
(
9
), pp.
1140
1154
.
4.
,
V.
,
Ivaldi
,
S.
,
Babič
,
J.
,
Mistry
,
M.
,
Peters
,
J.
, and
Nori
,
F.
,
2017
, “
Whole-Body Multi-Contact Motion in Humans and Humanoids: Advances of the CoDyCo European Project
,”
Rob. Auton. Syst.
,
90
, pp.
97
117
.
5.
Vukobratović
,
M.
, and
Borovac
,
B.
,
2004
, “
Zero-Moment Point - Thirty Five Years of Its Life
,”
Int. J. Humanoid Rob.
,
1
(
1
), pp.
157
173
.
6.
Hirukawa
,
H.
,
Hattori
,
S.
,
,
K.
,
Kajita
,
S.
,
Kaneko
,
K.
,
Kanehiro
,
F.
,
Fujiwara
,
K.
, and
Morisawa
,
M.
,
2006
, “
A Universal Stability Criterion of the Foot Contact of Legged Robots-Adios ZMP
,” IEEE International Conference on Robotics and Automation (
ICRA
), Orlando, FL, May 15–19, pp. 1976–1983.
7.
Firmani
,
F.
, and
Park
,
E. J.
,
2013
, “
Theoretical Analysis of the States of Balance in Bipedal Walking
,”
ASME J. Biomech. Eng.
,
135
(
4
), p.
041003
.
8.
Pratt
,
J.
, and
Tedrake
,
R.
,
2006
, “
Chapter: Velocity-Based Stability Margins for Fast Bipedal Walking
,”
Fast Motions in Biomechanics and Robotics Lecture Notes in Control and Information Sciences
, Vol.
340
,
Springer
,
Berlin
, pp.
299
324
.
9.
Aubin
,
J. P.
,
2009
,
Viability Theory
,
, New York.
10.
Wieber
,
P. B.
,
2008
, “
Viability and Predictive Control for Safe Locomotion
,” IEEE/RSJ International Conference on Intelligent Robots and Systems (
IROS
), Nice, France, Sept. 22–26, pp. 1103–1108.
11.
Wieber
,
P. B.
,
2002
, “
On the Stability of Walking Systems
,”
International Workshop on Humanoid and Human Friendly Robotics
, Tsukuba, Japan, pp. 53–59.
12.
Koolen
,
T.
,
De Boer
,
T.
,
Rebula
,
J.
,
Goswami
,
A.
, and
Pratt
,
J.
,
2012
, “
Capturability-Based Analysis and Control of Legged Locomotion—Part 1: Theory and Application to Three Simple Gait Models
,”
Int. J. Rob. Res.
,
31
(
9
), pp.
1094
1113
.
13.
Lanari
,
L.
,
Hutchinson
,
S.
, and
Marchionni
,
L.
,
2014
, “
Boundedness Issues in Planning of Locomotion Trajectories for Biped Robots
,”
14th IEEE-RAS International Conference on Humanoid Robots
(
Humanoids
), Madrid, Spain, Nov. 18–20, pp. 951–958.
14.
Koolen
,
T.
,
Posa
,
M.
, and
Tedrake
,
R.
,
2016
, “
Balance Control Using Center of Mass Height Variation: Limitations Imposed by Unilateral Contact
,”
16th IEEE-RAS International Conference on Humanoid Robots
(
Humanoids
), Cancun, Mexico, Nov. 15–17, pp. 8–15.
15.
Pratt
,
J.
,
Carff
,
J.
,
Drakunov
,
S.
, and
Goswami
,
A.
,
2006
, “
Capture Point: A Step Toward Humanoid Push Recovery
,” Sixth
IEEE-RAS
International Conference on Humanoid Robots, Genova, Italy, Dec. 4–6, pp. 200–207.
16.
Englsberger
,
J.
,
Ott
,
C.
, and
Albu-Schäffer
,
A.
,
2015
, “
Three-Dimensional Bipedal Walking Control Based on Divergent Component of Motion
,”
IEEE Trans. Rob.
, 31(2), pp.
355
368
.
17.
Posa
,
M.
,
Koolen
,
T.
, and
Tedrake
,
R.
,
2017
, “
Balancing and Step Recovery Capturability Via Sums-of-Squares Optimization
,”
Robotics: Science and Systems
,
Cambridge, MA
, July 12–16.
18.
Sherikov
,
A.
,
Dimitrov
,
D.
, and
Wieber
,
P. B.
,
2014
, “
Whole Body Motion Controller With Long-Term Balance Constraints
,”
14th IEEE-RAS International Conference on Humanoid Robots
(
Humanoids
), Madrid, Spain, Nov. 18–20, pp. 444–450.
19.
Sherikov
,
A.
,
Dimitrov
,
D.
, and
Wieber
,
P. B.
,
2015
, “
Balancing a Humanoid Robot With a Prioritized Contact Force Distribution
,”
15th IEEE-RAS International Conference on Humanoid Robots
(
Humanoids
), Seoul, South Korea, Nov. 3–5, pp. 223–228.
20.
Nagasaka
,
K.
,
Fukushima
,
T.
, and
Shimomura
,
H.
,
2012
, “
Whole-Body Control of a Humanoid Robot Based on Generalized Inverse Dynamics and Multi-Contact Stabilizer That Can Take Account of Contact Constraints
,” 17th
Robotics Symposium
(in Japanese), pp. 134–141.
21.
Audren
,
H.
,
Vaillant
,
J.
,
Kheddar
,
A.
,
Escande
,
A.
,
Kaneko
,
K.
, and
Yoshida
,
E.
,
2014
, “
Model Preview Control in Multi-Contact Motion-Application to a Humanoid Robot
,” IEEE/RSJ International Conference on Intelligent Robots and Systems (
IROS
), Chicago, IL, Sept. 14–18, pp. 4030–4035.
22.
Dai
,
H.
,
Valenzuela
,
A.
, and
Tedrake
,
R.
,
2014
, “
Whole-Body Motion Planning With Centroidal Dynamics and Full Kinematics
,” 14th IEEE-RAS International Conference on Humanoid Robots (
Humanoids
), Madrid, Spain, Nov. 18–20, pp.
295
302
.
23.
Nikolić
,
M.
,
Borovac
,
B.
, and
Raković
,
M.
,
2017
, “
Dynamic Balance Preservation and Prevention of Sliding for Humanoid Robots in the Presence of Multiple Spatial Contacts
,”
Multibody Syst. Dyn.
, 42(2), pp.
1
22
.
24.
Caron
,
S.
, and
Kheddar
,
A.
,
2016
, “
Multi-Contact Walking Pattern Generation Based on Model Preview Control of 3D COM Accelerations
,”
16th IEEE-RAS International Conference on Humanoid Robots
(
Humanoids
), Cancun, Mexico, Nov. 15–17, pp. 550–557.
25.
Mummolo
,
C.
,
Mangialardi
,
L.
, and
Kim
,
J. H.
,
2017
, “
Numerical Estimation of Balanced and Falling States for Constrained Legged Systems
,”
J. Nonlinear Sci.
,
27
(
4
), pp.
1291
1323
.
26.
Kim
,
J. H.
, and
Joo
,
C. B.
,
2014
, “
Numerical Construction of Balanced State Manifold for Single-Support Legged Mechanism in Sagittal Plane
,”
Multibody Syst. Dyn.
,
31
(
3
), pp.
257
281
.
27.
Kalyanakrishnan
,
S.
, and
Goswami
,
A.
,
2011
, “
Learning to Predict Humanoid Fall
,”
Int. J. Humanoid Rob.
,
8
(
2
), pp.
245
273
.
28.
Patton
,
J. L.
,
Pai
,
Y.
, and
Lee
,
W. A.
,
1999
, “
Evaluation of a Model That Determines the Stability Limits of Dynamic Balance
,”
Gait Posture
,
9
(
1
), pp.
38
49
.
29.
Perry
,
J.
,
1992
,
Gait Analysis. Normal and Pathological Function
,
Slack Incorporated
,
Thorofare, NJ
.
30.
Mummolo
,
C.
,
Mangialardi
,
L.
, and
Kim
,
J. H.
,
2013
, “
Quantifying Dynamic Characteristics of Human Walking for Comprehensive Gait Cycle
,”
ASME J. Biomech. Eng.
,
135
(
9
), p.
091006
.
31.
Pratt
,
J.
,
Koolen
,
T.
,
De Boer
,
T.
,
Rebula
,
J.
,
Cotton
,
S.
,
Carff
,
J.
,
Johnson
,
M.
, and
Neuhaus
,
P.
,
2012
, “
Capturability-Based Analysis and Control of Legged Locomotion—Part 2: Application to M2V2, a Lower Body Humanoid
,”
Int. J. Rob. Res.
,
31
(
10
), pp.
1117
1133
.
32.
Posa
,
M.
,
Cantu
,
C.
, and
Tedrake
,
R.
,
2014
, “
A Direct Method for Trajectory Optimization of Rigid Bodies Through Contact
,”
Int. J. Rob. Res.
,
33
(
1
), pp.
69
81
.
33.
Mummolo
,
C.
,
Mangialardi
,
L.
, and
Kim
,
J. H.
,
2014
, “Contact Status Optimization of Multibody Dynamic Systems Using Dual Variable Transformation,”
ASME
Paper No. DETC2014-34193.
34.
Kim
,
J. H.
,
Abdel‐Malek
,
K.
,
Xiang
,
Y.
,
Yang
,
J. J.
, and
Arora
,
J. S.
,
2011
, “
Concurrent Motion Planning and Reaction Load Distribution for Redundant Dynamic Systems Under External Holonomic Constraints
,”
Int. J. Numer. Methods Eng.
,
88
(
1
), pp.
47
65
.
35.
Ha
,
I.
,
Tamura
,
Y.
,
Asama
,
H.
,
Han
,
J.
, and
Hong
,
D.
,
2011
, “
Development of Open Humanoid Platform DARwIn-OP
,” SICE Annual Conference (
SICE
), Tokyo, Japan, Sept. 13–18, pp.
2178
2181.
36.
Ha
,
I.
,
Tamura
,
Y.
, and
Asama
,
H.
,
2011
, “
Gait Pattern Generation and Stabilization for Humanoid Robot Based on Coupled Oscillators
,” IEEE/RSJ International Conference on Intelligent Robots and Systems (
IROS
), San Francisco, CA, Sept. 25–30, pp. 3207–3212.
37.
Mummolo
,
C.
, and
Kim
,
J. H.
,
2013
, “
Passive and Dynamic Gait Measures for Biped Mechanism: Formulation and Simulation Analysis
,”
Robotica
,
31
(
4
), pp.
555
572
.
38.
ROBOTIS
, 2010, “ROBOTIS e-Manual v1.31.01,” ROBOTIS, Seoul, South Korea, accessed Mar. 9, 2017, http://support.robotis.com/en/product/actuator/dynamixel/mx_series/mx-28at_ar.htm
39.
Mummolo
,
C.
,
Park
,
S.
,
Mangialardi
,
L.
, and
Kim
,
J. H.
,
2016
, “
Computational Evaluation of Load Carriage Effects on Gait Balance Stability
,”
Comput. Methods Biomech. Biomed. Eng.
,
19
(
11
), pp.
1127
1136
.
40.
Gill
,
P. E.
,
Murray
,
W.
, and
Saunders
,
M. A.
,
2005
, “
SNOPT: An SQP Algorithm for Large-Scale Constrained Optimization
,”
SIAM Rev.
,
47
(
1
), pp.
99
131
.
41.
ROBOTIS, 2010, “Robotis e-Manual,” ROBOTIS, Seoul, South Korea, accessed Jan. 30, 2018, http://darwinop.sourceforge.net/
42.
Mummolo
,
C.
,
Mangialardi
,
L.
, and
Kim
,
J. H.
,
2015
, “
Identification of Balanced States for Multi-Segmental Legged Robots Using Reduced-Order Model
,” 15th IEEE-RAS International Conference on Humanoid Robots (
Humanoids
), Seoul, South Korea, Nov. 3–5, pp. 914–919.
43.
Mummolo
,
C.
,
Cursi
,
F.
, and
Kim
,
J. H.
,
2016
, “
Balanced and Falling States for Biped Systems: Applications to Robotic versus Human Walking Stability
,” 16th IEEE-RAS International Conference on Humanoid Robots (
Humanoids
), Cancun, Mexico, Nov. 15–17, pp. 1155–1160.