## Abstract

A safety-critical measure of legged locomotion performance is a robot's ability to track its desired time-varying position trajectory in an environment, which is herein termed as “global-position tracking.” This paper introduces a nonlinear control approach that achieves asymptotic global-position tracking for three-dimensional (3D) bipedal robots. Designing a global-position tracking controller presents a challenging problem due to the complex hybrid robot model and the time-varying desired global-position trajectory. Toward tackling this problem, the first main contribution is the construction of impact invariance to ensure all desired trajectories respect the foot-landing impact dynamics, which is a necessary condition for realizing asymptotic tracking of hybrid walking systems. Thanks to their independence of the desired global position, these conditions can be exploited to decouple the higher-level planning of the global position and the lower-level planning of the remaining trajectories, thereby greatly alleviating the computational burden of motion planning. The second main contribution is the Lyapunov-based stability analysis of the hybrid closed-loop system, which produces sufficient conditions to guide the controller design for achieving asymptotic global-position tracking during fully actuated walking. Simulations and experiments on a 3D bipedal robot with twenty revolute joints confirm the validity of the proposed control approach in guaranteeing accurate tracking.

## 1 Introduction

A robot's global position represents its absolute position in an environment. Poor global-position tracking can potentially put the safety of both humans and robots at risk, for example, by causing robots' failure to avoid pedestrians in human-populated environments. To achieve accurate global-position tracking, the zero-moment-point control approach has been introduced based on the zero-moment-point balance criterion and the continuous-time dynamic model of bipedal walking [13]. Yet, bipedal walking is inherently a hybrid process involving both continuous motions (e.g., foot swinging) and discrete impact dynamics (e.g., sudden joint-velocity jumps upon a foot landing) [47]. Achieving reliable global-position tracking by explicitly addressing the hybrid robot dynamics presents substantial challenges.

This study focuses on addressing the challenges associated with: (a) lower-level trajectory generation (i.e., level 2 in Fig. 1) and (b) controller design (i.e., level 3 in Fig. 1). It is assumed that the desired global path and the desired time-varying position trajectory along the path have both been provided by a higher-level planner (i.e., level 1 in Fig. 1) without impact dynamics considered.

Fig. 1
Fig. 1
Close modal

One challenge in the lower-level trajectory generation (i.e., level 2 in Fig. 1) for a hybrid robot model is to respect both continuous dynamics and the discrete impact dynamics, which is computational heavier than just respecting the continuous dynamics [8]. For the controller to achieve asymptotic tracking based on a hybrid robot model, the desired trajectories need to agree with the impact dynamics; i.e., their pre- and postimpact values should satisfy the impact map. This is because the impact dynamics cannot be directly controlled due to their infinitesimally short duration [912]. Yet, the additional computational load caused by respecting the nonlinear impact map could be amplified when frequent replanning is needed during real-world tasks such as dynamic obstacle avoidance in cluttered environments.

Another challenge is the closed-loop stability analysis of the hybrid dynamical system that produces sufficient conditions to guide the controller derivation (i.e., level 3 in Fig. 1). Such a stability analysis is complex because a closed-loop system capable of stabilizing a time-varying global-position trajectory is hybrid, nonlinear, and time-varying with uncontrolled, state-triggered impact dynamics.

### 1.1 Related Work on Orbitally Stabilizing Control.

The most widely studied control approach that explicitly addresses the hybrid walking dynamics is the hybrid zero dynamics (HZD) method [1318]. The HZD method provably stabilizes dynamic walking motions through orbital stabilization of the hybrid closed-loop control system. It has realized remarkable performance for various gait types such as periodic underactuated [19,20], fully actuated [21], and multidomain walking [22].

The HZD framework introduces virtual constraints to represent the evolution of a robot's desired configuration with respect to a phase variable that indicates how far a step has progressed. To enforce the impact dynamics on the desired gait, the HZD approach introduces a method termed “impact invariance construction” to produce an equality constraint under which the desired gait respects the impact dynamics and incorporates the constraint in the optimization-based generation of virtual constraints. For systems with a single degree of underactuation [14], the HZD approach ensures the impact invariance of the zero dynamics manifold by imposing the impact invariance of a single point on it. Yet, because the encoding of the global-position trajectory is inherently different from that of virtual constraints, the previous impact invariance construction cannot be directly applied or extended to ensure the agreement with impact dynamics for the desired global-position trajectory. Specifically, the virtual constraints are encoded by a local phase variable that is reset at the beginning of a walking step, while the desired global-position trajectory is usually encoded by a global phase variable that evolves continuously and monotonically across all walking steps.

To analyze the closed-loop stability for guiding controller designs, the HZD approach exploits the Poincaré section method to examine the asymptotic convergence of a robot's state to the desired periodic orbit representing the desired gait in the state space. Recently, the HZD framework has been extended to achieve asymptotic tracking of the desired global path during 3D underactuated bipedal walking [23]. Yet, an orbitally stabilizing controller cannot stabilize a prespecified time-varying trajectory [24] such as the desired global-position trajectory.

### 1.2 Related Work on Trajectory Tracking Control.

Our previous trajectory tracking controller designs either focus on individual joint trajectory tracking [25,26] or only consider 2D walking [2729]. In particular, our previous work on 2D walking, including the impact invariance construction and stability analysis, is not valid for 3D robots. Specifically, the walking dynamics of 3D robots are nonlinearly coupled in the heading and lateral directions of the robot's global path, but 2D walking does not exhibit lateral motion, and accordingly, the coupling is trivial. This nonlinear coupling significantly increases the complexity of controller derivation in addressing 3D robots compared with 2D robots. Furthermore, experimental validation of these previous controllers has been missing.

Beyond the scope of global-position tracking control for bipedal walking robots, trajectory tracking control of general hybrid systems with state-triggered jumps is an active research topic [3035]. Lyapunov-based controller design methodologies have been introduced to provably achieve asymptotic trajectory tracking for linear hybrid systems [31,32]. In this study, to guide the needed controller design, we will extend the previous Lyapunov-based stability analysis to nonlinear hybrid systems that include 3D bipedal robots during fully actuated walking.

### 1.3 Contributions.

This study aims to derive and experimentally validate a nonlinear walking control approach for 3D bipedal robots that achieves asymptotic global-position tracking by explicitly addressing the hybrid robot dynamics. The main contributions of this study are summarized as follows:

• Constructing impact invariance conditions that are independent of the desired global-position trajectory and yet ensure all desired trajectories respect the impact dynamics. They can be used to decouple the planning of virtual constraints and global position, thus improving trajectory generation efficiency.

• Establishing sufficient conditions based on the multiple Lyapunov stability analysis [36] of the hybrid system for guiding the design of a continuous state-feedback control law to achieve asymptotic global-position tracking.

• Demonstrating the global-position tracking accuracy of the proposed control approach both through simulations and experimentally on a 3D bipedal walking robot.

• Experimentally validating the inherent robustness of the proposed control design in addressing irregular walking surfaces such as moderately slippery floors.

Some of the results presented in this paper were initially reported in Refs. [37] and [38]. This paper includes substantial, new contributions in the following aspects: (a) the proof of the main theorem (i.e., Theorem 1) is updated with a new choice of Lyapunov function to properly analyze the convergence of the robot's lateral foot placement, and Proposition 4 is added along with its full proof to support the proof of the main theorem; (b) fully developed proofs of all theorems and propositions are presented, which were missing in Refs. [37] and [38]; (c) comparative experiment results are added to show the reliable global-position tracking performance of the proposed control approach; and (d) robustness evaluation is newly included to illustrate the capability of the proposed method in handling relatively slippery grounds.

This paper is structured as follows. Section 2 describes the problem formulation. Section 3 explains the proposed continuous-phase tracking control law. Section 4 presents the proposed construction of impact invariance conditions for designing virtual constraints. Section 5 introduces the closed-loop stability analysis based on multiple Lyapunov functions. Section 6 reports the simulation and experiment results. Section 7 discusses the proposed approach and potential directions of future work. Proofs of all theorems and propositions are given in the Appendix.

## 2 Problem Formulation

This section presents the problem formulation of global-position tracking control. The formulation includes dynamics modeling and tracking error design.

### 2.1 Full-Order Robot Model.

This section describes a full-order model that accurately captures the dynamic behaviors of all degrees-of-freedom (DOFs) involved in bipedal walking. Thanks to the model's accuracy, a controller that is effective for the model would also be valid for the physical robot. Hence, we use the full-order model as a basis of the proposed control approach.

The full-order model is naturally hybrid and nonlinear, because walking dynamics are inherently hybrid, involving both nonlinear continuous behaviors (e.g., leg-swinging motions) and state-triggered discrete behaviors (e.g., the joint-velocity jumps caused by foot-landing impact).

In this study, we assume that the swing and the stance legs immediately switch roles upon a foot landing, with the new swing leg beginning to move in the air and the new stance leg remaining in a full, static contact with the ground until the next landing occurs [14]. The assumption is valid when the double-support phase is sufficiently short and when the stance foot does not notably slip on the ground.

Under this assumption, if all of the robot's (revolute or prismatic) joints are directly actuated, then the robot is fully actuated; i.e., its full DOFs can be directly commanded within continuous phases.

This study focuses on the relatively simple gait, fully actuated gait, for two main reasons. First, asymptotic tracking of time-varying global-position trajectories for the 3D hybrid robot model is still an open control problem for this simple gait. Second, using a simple gait allows us to focus on addressing the complexity of the controller design problem induced by the hybrid, nonlinear robot dynamics, and the time-varying global-position trajectory.

Continuous-phase dynamics. As illustrated in Fig. 2, a complete walking cycle comprises: (a) a fully actuated continuous phase during which one foot contacts the ground and the other swings in the air and (b) a landing impact.

Fig. 2
Fig. 2
Close modal
Walking dynamics during continuous phases can be described by usual ordinary differential equations. Lagrange's method is used to obtain the following nonlinear full-order model during continuous phases [13]
$M(q)q¨+c(q,q˙)=Buu$
(1)

where $q∈Q$ is the joint-position vector, $M:Q→ℝn×n$ is the symmetric, positive-definite inertia matrix, $c:TQ→ℝn$ is the sum of Coriolis, centrifugal, and gravitational terms, $Bu∈ℝn×m$ is the joint-torque projection matrix with full column rank, and $u∈U$ is the joint-torque vector. Here, $Q⊂ℝn$ is the configuration space of the robot, TQ is the tangent bundle of Q, and $U⊂ℝm$ is the admissible joint-torque set. Note that m = n when a robot is fully actuated.

Impact dynamics. When the swing foot lands on the ground, the swing and stance legs immediately switch their roles. Here, we model the swing-foot landing impact as the contact between rigid bodies [14]. This assumption is valid for dynamic walking on relatively stiff surfaces (e.g., concrete and ceramic floors) during which the swing foot strikes the surface at a relatively significant downward velocity.

Due to the coordinate swap of the swing and stance legs as well as the impulsive rigid-body impact, both joint position and velocity vectors experience a sudden jump at a landing event. This state-triggered jump is described by the nonlinear reset map $Δq,q˙:TQ→TQ$ [13] given by
$[q+q˙+]=Δq,q˙(q−,q˙−):=[Δq(q−)Δq˙(q−)q˙−]$
(2)

where $⋆−$ and $⋆+$ represent the values of $⋆$ just before and just after the impact, respectively.

Switching surface. A swing-foot landing event is triggered when the robot's state reaches the switching surface Sq, which is expressed as
$Sq:={(q,q˙)∈TQ:zsw(q)=0, z˙sw(q,q˙)<0}$
(3)

where $zsw:Q→ℝ$ is the swing-foot height above the ground.

Combining the above equations yields the following hybrid full-order model:
${M(q)q¨+c(q,q˙)=Buu,if (q−,q˙−)∉Sq[q+q˙+]=Δq,q˙(q−,q˙−),if (q−,q˙−)∈Sq$
(4)

### 2.2 Global-Position Tracking Error.

A fully actuated, n-DOF bipedal robot can track n independent desired position trajectories, including the reference global-position trajectories.

In this study, we choose to use the position of a biped's base (e.g., trunk), $(xb,yb,zb)$, to represent its global position in an environment. The horizontal components of the base position are related to the stance-foot position as
$xb=xst+x¯b(q) and yb=yst+y¯b(q)$
(5)

where $(xst,yst,0)$ denotes the stance-foot position with $xst,yst∈ℝ$. The scalar variables $x¯b:Q→ℝ$ and $y¯b:Q→ℝ$ represent the x- and y-coordinates of the base position relative to the stance foot, respectively.

In real-world locomotion tasks, a higher-level planner typically specifies the desired global motions as:

• The centerline Γd of the desired global path.

• The desired smooth position trajectory $sd(t)$ along Γd.

As an arbitrary curved path can be approximated as a nonsmooth curve pieced together by straight lines, this study focuses on the tracking control of straight-line paths, which could be extended to the tracking of a curved path as briefly explained in Sec. 8.

Without loss of generality, suppose that the centerline Γd coincides with the Xw-axis of the world frame; that is
$Γd={(xb,yb)∈ℝ2:yb=0}$

Then, the global-position tracking error $h1(t,q)$ along Γd is defined as $h1(t,q):=x¯b(q)−(sd(t)−xst)$.

While $sd(t)$ and Γd are often provided by a higher-level path planner, the desired base motion in the direction lateral to Γd remains to be designed, which is explained next.

### 2.3 Virtual-Constraint Tracking Error.

Besides the desired global-position trajectory $sd(t)$, a legged robot typically has multiple directly actuated DOFs that can track additional desired motions. We choose to use virtual constraints to define the desired trajectories for the lateral base position yb and the remaining control variables $ϕc:Q→ℝn−2$.

Analogous to the HZD framework, we use virtual constraints to represent the desired configuration relative to a phase variable $θ:Q→Qf⊂ℝ$. The phase variable represents how far a step has progressed. Without loss of generality, the phase variable θ is chosen as the forward position of the base relative to the support foot, $x¯b(q)$; that is,
$θ=x¯b(q)$
(6)
The virtual constraints can be encoded by θ as
$[y¯b(q)ϕc(q)]−[yd(θ(q))−ystϕd(θ(q))]=0$
(7)

where the scalar function $yd:Qf→ℝ$ and the vector $ϕd:Qf→ℝn−2$ are the desired trajectories of yb and $ϕc$, respectively. Suppose that yd, $ϕd, ϕc$, and θ are all continuously differentiable in their respective arguments.

Thus, the tracking error $h2(q)$ corresponding to the virtual constraints is defined as $h2(q):=[y¯b(q)ϕc(q)]−[yd(θ(q))−ystϕd(θ(q))]$.

### 2.4 Control Objective.

The tracking errors can be compactly expressed as
$h(t,q)=[h1(t,q), h2T(q)]T=hc(q)−hd(t,θ(q))$
(8)

where the control variables hc and their desired trajectories hd are, respectively, defined as $hc:=[x¯b,y¯b,ϕcT]T$ and $hd:=[sd−xst,yd−yst,ϕdT]T$.

The control objective is to asymptotically drive the tracking error $h$ to zero for achieving asymptotic tracking of the desired motions, which are the desired global-position trajectory $sd(t)$ and the desired functions yd and $ϕd$ that define the virtual constraints, for 3D bipedal robots walking along the given straight line.

To achieve this objective, the proposed control approach (Fig. 1) comprises three main components: (a) continuous-phase controller design (for stabilizing the desired trajectories within continuous phases); (b) impact invariance construction (for satisfying a necessary condition of asymptotic tracking for the hybrid model); and (c) closed-loop stability analysis (for providing sufficient stability conditions that guide the controller design).

## 3 Continuous-Phase Control

This section presents a continuous state-feedback control law that asymptotically stabilizes the desired trajectories within continuous phases.

We choose to design a controller that directly regulates the continuous-phase walking dynamics instead of the impact dynamics because the impact dynamics cannot be directly commanded due to its infinitesimally short duration. We will show in Sec. 5 how the proposed continuous control law could be tuned to indirectly stabilize the desired trajectories for the overall hybrid system.

The proposed control law (Fig. 3) is synthesized based on the full-order model of bipedal walking dynamics. Analogous to the HZD framework, we utilize the input–output linearization technique [24] to linearize the nonlinear continuous-phase dynamics in Eq. (1) into a linear map, which allows us to exploit the well-studied linear system theory to design the needed controller for the continuous phase.

Fig. 3
Fig. 3
Close modal
We define the output function as the tracking error $h$
$y:=[y1, y2T]T:=[h1, h2T]T=h$
(9)
Then a continuous-phase control law synthesized via input–output linearization is given by
$u=(JhM−1Bu)−1(v+∂h∂qM−1c−∂2h∂t2−∂∂q(∂h∂qq˙)q˙)$
(10)

with $Jh(q):=∂h∂q(t,q),$ which yields the linearized dynamics $y¨=v$. Note that the variables $ϕc$, yd, and $ϕd$ can be chosen such that there exists an open subset $Q̃$ of the configuration space Q on which the Jacobian matrix $Jh(q)$ is invertible. Then, the matrix $JhM−1Bu$ is invertible on $q∈Q̃$.

Choosing $v$ as a proportional-derivative (PD) term
$v=−KPy−KDy˙$
(11)

where the proportional gain matrix $KP∈ℝn×n$ and the derivative gain matrix $KD∈ℝn×n$ are both positive-definite diagonal matrices, the linear closed-loop dynamics of the output function becomes $y¨+KDy˙+KPy=0$ during continuous phases.

Define the state of the output function dynamics as $x:=[yT,y˙T]T∈ℝ2n$. Then, the closed-loop error equation can be compactly expressed as
${x˙=Ax, if (t−,x−)∉Sx+=Δ(t−,x−), if (t−,x−)∈S$
(12)

Here, $A:=[0I−KP−KD]$ with $0$ a zero matrix and $I$ an identity matrix with appropriate dimensions. S and $Δ$ are the switching surface and impact map associated with the closed-loop dynamics, respectively. Note that $Δ$ is explicitly time-dependent because of the explicit time dependence of $h$. The expressions of S and $Δ$ can be obtained from their counterparts in the open-loop dynamics (Eqs. (3) and (2)) as well as the output function definition (Eq. (8)).

The origin (i.e., $x=0$) of the continuous-phase closed-loop dynamics (i.e., $x˙=Ax$) will be asymptotically stable if the PD gains are chosen such that $A$ is Hurwitz [24]. Then, there exist positive numbers c1, c2, and c3 and a Lyapunov function candidate $V(x)$ such that
$c1||x||2≤V(x)≤c2||x||2 and V˙(x)≤−c3||x||2$
(13)

hold for all $x$ within continuous phases. These inequalities indicate that $V(x)$ exponentially converges at the rate of $c3c2$ within a continuous phase.

While the proposed control law with properly chosen PD gains guarantees the asymptotic tracking of the desired trajectories within continuous phases, the impact dynamics (i.e., $x+=Δ(t,x−)$) remain uncontrolled, and thus, the stability of the hybrid closed-loop system is not yet ensured. To satisfy a necessary condition of asymptotic trajectory stabilization in the presence of uncontrolled impact dynamics, we introduce impact invariance construction next.

## 4 Impact Invariance Construction for Virtual Constraint Design

This section derives impact invariance conditions that can be incorporated in the trajectory generation of the desired functions yd and $ϕd$, which define the virtual constraints, for ensuring all desired trajectories (i.e., yd and $ϕd$ as well as sd) respect the impact dynamics.

For the proposed feedback controller to achieve asymptotic tracking for hybrid dynamical systems, the output function state variables $y$ and $y˙$ need to satisfy the following condition across an impact event at the steady-state: $y+=0$ and $y˙+=0$ hold just after the impact if $y−=0$ and $y˙−=0$ hold just before the impact. Suppose this condition is not met at the steady-state. Then, because the robot's impact dynamics cannot be directly regulated, the output function state may become nonzero just after an impact even if it is zero just before the impact, which means asymptotic tracking cannot be achieved.

To mathematically describe this condition, we introduce the manifold Z given by
$Z:={(t,q,q˙)∈ℝ×TQ:h(t,q)=0,h˙(t,q,q˙)=0}$
(14)

Here, Z is a one-dimensional embedded submanifold of $ℝ×TQ$. The corresponding guard Sa is defined by rewriting Sq as: $Sa:={(t,q,q˙)∈ℝ×TQ:zsw(q)=0,z˙sw(q,q˙)<0}.$

Definition 1. (Impact invariance) The manifold Z is called impact invariant if$(t+,Δq,q˙(q−,q˙−))∈Z$holds for any$(t−,q−,q˙−)∈Z∩Sa$; that is, the manifold is invariant across the impact event.

Remark 1. (Differentiation from the HZD approach) The concept of impact invariance was first introduced within the HZD framework, along with a systematic method of impact invariance construction (see Theorem 4 in Ref. [14]). The concept was later on termed as “impact invariance” [15].

The equations defining Z are time-dependent in this study whereas in the original HZD framework [14,15], the hybrid zero dynamics manifold is time-independent. The difference is essentially due to their different control objectives. In the HZD framework, the controller aims to track the desired walking pattern encoded by a configuration-based phase variable, resulting in a time-invariant definition of output functions and accordingly a time-invariant hybrid zero dynamics manifold. In contrast, the control objective here is to track time-varying global-position trajectories, and thus, the output function (specifically, $h1(t,q)$) is time-varying, inducing the time dependence of the submanifold Z.

Another difference is that, in the HZD framework [14,15], the system of interest is underactuated, and thus, the dimension of the impact invariant manifold, when restricted to the guard, can be higher than zero. Yet, in our case of fully actuated systems, the dimension of $Z∩Sa$ is zero because $Z∩Sa$ is a single point. Although ensuring that a single point respects the impact map is typically easier for trajectory optimization, the proposed impact invariance construction for Z has an attractive property for efficient planning, as revealed by Remark 3 in Sec. 4.4.

We choose to construct the manifold Z to be impact invariant by properly planning the desired function hd. As the desired global position sd is often supplied by a higher-level path planner without impact dynamics considered, the generation of the remaining desired functions yd and $ϕd$, which define the virtual constraints, needs to ensure the impact agreement for all trajectories (i.e., sd, yd, and $ϕd$). To this end, the proposed impact invariance construction boils down to the derivation of conditions that the virtual constraints should satisfy in order to ensure Z is impact invariant. We call these conditions “impact invariance conditions.”

The proposed impact invariance construction consists of two steps corresponding to two sets of impact invariance conditions. We first extend the existing HZD method (i.e., Theorem 4 in Ref. [14]) to derive conditions that ensure the impact invariance of the three-dimensional embedded submanifold of $ℝ×TQ$ associated with the virtual constraint, that is
$Z̃:={(t,q,q˙)∈ℝ×TQ:h2(q)=0,h˙2(q,q˙)=0}$

Then, we introduce a new, additional condition that, in combination with the first set of conditions, guarantees the impact invariance of Z. Note that both conditions are placed on the virtual constraints alone.

These two steps of impact invariance construction are introduced in Secs. 4.3 and 4.4, respectively. Before presenting them, we first explain the timing and the unique robot configuration associated with a landing impact in Secs. 4.1 and 4.2, which are needed for the derivation of the proposed impact invariance conditions.

### 4.1 Impact Timings.

Because the desired global-position trajectory sd is explicitly time-varying, we need to consider the impact timings in the proposed impact invariance construction. As the actual and desired impact timings generally do not coincide due to the state-triggered nature of a foot-landing event [10], they are individually defined as follows.

Definition 2. (Actual and desired impact timings) Let Tk be the timing of the kth ($k∈ℤ+$) actual landing impact, which is defined as the timing of the first intersection between the state$x$and the switching surface S on$t>Tk−1+$. Without loss of generality, define$T0=0$. Let τk denotes the kth desired impact timing, which is defined as the timing of the first intersection between$x$and S on$t>Tk−1+$assuming$x=0∀t>Tk−1+$.

The precise definition of Tk is given in Ref. [14]. Figure 4 shows an illustration of Tk. The variables $⋆(Tk−1−)$ and $⋆(Tk−1+)$ are, respectively, denoted as $⋆|k−1−$ and $⋆|k−1+$ in the rest of the paper where brevity is preferred.

Fig. 4
Fig. 4
Close modal

### 4.2 Unique Configuration.

The proposed impact invariance construction utilizes the uniqueness of the robot's joint position $q*$ just before an impact event when the virtual constraints in Eq. (7) are exactly satisfied. Note that although the proposed construction relies on the unique configuration, it does not require the joint velocity should be unique just before an impact.

The joint position $q*$ is mathematically defined as the solution to the following equations:
$F(q):=[y¯b(q)−(yd(θ(q))−yst)ϕc(q)−ϕd(θ(q))zsw(q)]=0$
(15)

on $S∩Q̃$. Note that the last equation in Eq. (15) holds because the swing-foot height $zsw(q)$ reaches zero at a touchdown.

Due to the nonlinearity of the function $F(q)$, Eq. (15) may have multiple solutions on $S∩Q̃$. Suppose that the output function is designed such that $∂F∂q(q*)$ is invertible on $S∩Q̃$. Then by the implicit function theorem, there exits $Q¯⊂Q̃$ such that $q*$ is a unique solution to $F(q)=0$ on $S∩Q¯$.

### 4.3 Impact Invariance Construction for Virtual Constraints.

We are now ready to introduce the conditions that ensure the impact invariance of $Z̃$. The impact invariance conditions are built upon the uniqueness of the joint position $q*$ on $S∩Q¯$. From Eq. (15), we know the value of $q*$ depends on the lateral foot placement yst. The following proposed impact invariance conditions use the value of $q*$ associated with the desired lateral foot placement ystd.

Proposition 1. (Impact invariance conditions for$Z̃$) Suppose that the desired functions yd and$ϕd$are planned to meet the following conditions:

1. $y¯b(q0)=yd(θ0)−ystd and ϕc(q0)=ϕd(θ0)$

2. $([∂y¯b∂q(q0)∂ϕc∂q(q0)]−[∂yd∂θ(θ0)∂ϕd∂θ(θ0)]∂x¯b∂q(q0))v0=0$

Here,$q0:=Δq(q*)$, $θ0:=x¯b(q0)$, $θ*:=θ(q*)$, $Jhc(q*):=∂hc∂q(q*)$, and$v0:=Δq˙(q*)Jhc−1(q*)[1∂yd∂θ(θ*)∂ϕd∂θ(θ*)]$. Then, under the lateral foot-placement condition yst = ystd, the impact invariance of$Z̃$holds.

From Sec. 4.2, we know that Eq. (15) has a unique solution on $S∩Q¯$ when yst = ystd; that is, the robot has a unique configuration $q*$ just before an impact if $y2(τk−)=0$ and yst = ystd. Given the uniqueness of $q*$, the equations in condition (A1), which are imposed on the virtual constraints, ensure that $y2(τk+)=0$ if $y2(τk−)=0$ and yst = ystd. Similarly, the equation in condition (A2) guarantees that $y˙2(τk+)=0$ if $y2(τk−)=0, y˙2(τk−)=0$, and yst = ystd.

Remark 2. (Differentiation from the HZD approach) The proposed construction of the impact invariant manifold $Z̃$ is analogous to the original HZD method [14,15]. The first difference lies in that the output function in our case is explicitly a function of the lateral foot placement yst and that the proposed construction is for the case where the desired foot placement yst = ystd is realized. The second difference is that we define the submanifold $Z̃$ based on $ℝ×TQ$ instead of just the tangent bundle TQ. This is because the impact invariance construction for $Z̃$ is used as a basis for rendering the manifold Z impact invariant, and Z is associated with the time-varying global-position tracking error $h1(t,q)$.

### 4.4 Impact Invariance Construction for Global-Position Tracking Error.

As the desired global-position trajectory sd is often supplied by a high-level planner without impact dynamics considered, we construct an additional condition, which is placed on the virtual constraints, to further ensure the impact invariance associated with the global-position error state, i.e., $x¯b−(sd−xst)$ and its first derivative. Note that $x¯b−(sd−xst)≡xb−sd$. This additional condition, together with those introduced in Sec. 4.3, guarantees that the submanifold Z is impact invariant.

The key to the proposed construction is to exploit the property of sd that it is commonly planned as a smooth function for any $t>T0$. Thanks to this property, $xb−sd=0$ automatically holds just after an impact if it holds just before the impact. This is because both the forward base position $xb$ and its desired trajectory sd are continuous across an impact.

To ensure $x˙b−s˙d=0$ holds just after an impact if it holds just before the impact, we choose to enforce the continuity of the global velocity $x˙b$ across the planned impact event. The rationale of this design choice is threefold. First, given the continuity of $s˙d$ for any $t>T0$, the continuity of $x˙b$ across the planned impact event guarantees the continuity of $x˙b−s˙d$, which then ensures that $x˙b−s˙d=0$ holds just after the planned impact if it holds just before the impact. Second, the continuity of $x˙b$ is equivalent to that of $x¯˙b$ because the stance foot does not move (i.e., $x˙st=0$). Third, $x¯˙b$ is a function of the joint position $q$ and velocity $q˙$ only, and thus its continuity across the planned impact event can be satisfied through virtual constraint design alone without explicitly relying on the profile of sd.

The proposed conditions for the impact invariance of Z is summarized as follows.

Proposition 2. (Impact invariance condition for Z) Suppose that the desired functions yd and$ϕd$satisfy conditions (A1) and (A2) and the following condition:

1. (A3)$∂x¯b∂q(q0)Δq˙(q*)Jh−1(q*)[1dyddθ(θ*)dϕddθ(θ*)]=1$

Then, under the lateral foot-placement condition yst = ystd, the impact invariance of Z holds.

Condition (A3) ensures that the base velocity does not jump (i.e., $x˙b(τk+)=x˙b(τk−)$) at the impact event if $y2(τk−)=0, y˙2(τk−)=0$, and yst = ystd and if conditions (A1) and (A2) in Proposition 1 also hold. Furthermore, if the virtual constraints are generated to meet the conditions in Proposition 2, which contains the conditions from Proposition 1, then under the lateral foot-placement condition yst = ystd, the impact invariance of Z holds; that is, if $x(τk−)=0$ then $x(τk+)=Δ(τk−,0)=0$.

Remark 3. (Independence from desired global-position trajectory) Propositions 1 and 2 indicate that the satisfaction of the impact invariance conditions only relies on the design of the virtual constraints but not the arbitrary global-position trajectory sd provided by a higher-level planner. For this reason, the design of virtual constraints does not need to explicitly consider sd and thus can be performed offline even when the higher-level planner updates sd online. This could reduce the computational load for online planning especially for mobility tasks that could frequently demand the replanning of sd (e.g., dynamic obstacle avoidance).

Remark 4. (Ensuring the desired lateral foot placement through controller design) Note that the foot-placement condition yst = ystd underlying the proposed impact invariance construction is only assumed in the virtual constraint planning but not the controller design. Indeed, Sec. 5 introduces sufficient conditions under which the proposed controller guarantees this foot-placement condition holds at the actual steady-state.

Remark 5. (Planning virtual constraints offline for impact invariance) There is a relatively simple two-step procedure to plan virtual constraints offline that meet the impact invariance conditions in Propositions 1 and 2. Recall the virtual constraints are given by: $y¯b−yd(θ)+yst=0$ and $ϕc−ϕd(θ)=0$. The first step is to plan desired time trajectories for the control variables $y¯b$ and $ϕc$ within a complete hybrid walking cycle (i.e., a continuous phase and a landing impact), which respect the impact dynamics with a constant forward velocity imposed across the impact. Let $ỹd(t)−ystd$ and $ϕ̃d(t)$ denote these time trajectories, respectively. Let $θ̃d(t)$ be the desired fictitious time trajectory for the phase variable θ (i.e., $x¯b$). The function $θ̃d(t)$ is only used for the offline planning and can be prespecified as any function that monotonically increases in time t within the planned walking cycle. Let $θ̃d−1$ represent the inverse of $θ̃d(t)$, which exists within the cycle. The next step is to obtain the virtual constraints by defining the desired functions $yd(θ)$ and $ϕd(θ)$ via $yd(θ)=ỹd(θ̃d−1(θ))$ and $ϕd(θ)=ϕ̃d(θ̃d−1(θ))$. Considering Remark 3, we can prove that the resulting virtual constraints automatically satisfy Propositions 1 and 2 even when the desired global-position trajectory $sd(t)$ is different from the fictitious desired trajectory $θ̃d(t)+xst$.

## 5 Stability Analysis

This section introduces Lyapunov-based stability analysis of the hybrid, nonlinear, time-varying closed-loop error dynamics (Eq. (12)) under the proposed continuous-phase control law (Eqs. (10) and (11)). The outcome of this stability analysis is a set of sufficient conditions under which the proposed control law provably realizes asymptotic stabilization of the desired global-position trajectory sd and the desired functions yd and $ϕd$ for the overall hybrid system.

### 5.1 Boundedness of Foot Placement and Impact Timing.

Before presenting the main theorem on closed-loop stability, we first introduce the boundedness of the impact timing Tk and the lateral stance-foot position yst. The boundedness of the impact timing is needed in the stability analysis to derive how much a Lyapunov function converges within a continuous phase. The boundedness of yst also needs to be explicitly considered, because yst = ystd underlies the proposed impact invariance conditions and should hold at the actual steady-state for achieving asymptotic tracking.

Proposition 3. (Boundedness of impact timing error) Let$x̃(t;t0,λ0)$be a solution of a fictitious continuous-time system$x̃˙=Ax̃$with the initial condition$x̃(t0)=λ0, ∀t>t0$. There exists a positive number r1 and a Lipschitz constant$LTx$such that the difference between the actual and the planned impact timings is bounded above in norm as
$|Tk−τk|≤LTx||x̃(τk;Tk−1+,x|k−1+)||$
(16)

for any$x|0+∈Br1(0):={x∈ℝ2n:||x||≤r1}$and any$k∈ℤ+$.

Proposition 4. (Boundedness of lateral foot-placement error) Suppose that the lateral swing-foot position ysw is chosen as an element of$ϕc$and is thus directly controlled. Then, there exist positive numbers βst and d1 such that the foot-placement error just after the kth swing-foot landing is bounded above in norm as
$|yst|k+−ystd|≤||x|k−||+βst||x̃(τk;Tk−1+,x|k−1+)||$
(17)

for any$x|0+∈Bd1(0):={x∈ℝ2n:||x||≤d1}$and any$k∈ℤ+$.

Rationale of proofs. The full proofs of Propositions 3 and 4 are given in the Appendix. The proof of Proposition 3 utilizes the implicit dependence of the actual impact timing Tk on the error state $x$. The proof of Proposition 4 mainly relies on the fact that the stance-foot position within the current step is the end position of the swing foot within the previous step. Thus, by including ysw as a control variable, the lateral foot-placement error $yst−ystd$ is also contained in the state $x$.▪

### 5.2 Main Theorem.

If the virtual constraints are designed to satisfy the impact invariance conditions in Propositions 1 and 2 and if the continuous-phase convergence rate of $x$ is sufficiently fast, then the origin of the hybrid closed-loop error system is asymptotically stable, as summarized in the main theorem:

Theorem 1. (Closed-loop stability conditions) Suppose that the virtual constraints satisfy the impact invariance conditions (A1)(A3). Also, suppose that the PD gains in Eq. (11) are chosen such that$A$is Hurwitz and that the continuous-phase convergence rate of$x$is sufficiently fast. Then, there exists a positive number d2 such that for any$x|0+∈Bd2(0):={x∈ℝ2n:||x||≤d2}$, the origin of the closed-loop error system in Eq. (12) is locally asymptotically stable; that is,$x(t)→0 as t→∞.$

Furthermore, both the lateral foot placement and actual impact timing asymptotically converge to their desired values; that is, $Tk−τk→0 and yst−ystd→0 as k→∞.$

Rationale of proof. The full proof of Theorem 1 is given in the Appendix. The proof utilizes the stability theory of the multiple Lyapunov functions [36], which prescribes how a Lyapunov function candidate should evolve in order for the origin of a hybrid dynamical system to be stable.

The stability analysis begins with the construction of the Lyapunov function candidate. Since the lateral foot-placement error $yst−ystd$ directly affects the satisfaction of the impact invariance conditions and thus the system stability, we choose to construct the Lyapunov function Va by augmenting V with a positive-definite function of the foot-placement error
$Va(x,yst−ystd):=V(x)+σ(yst−ystd)2$
(18)

where σ is a positive number to be specified in the proof.

Next, we analyze the evolution of Va during a continuous phase as well as through a hybrid transition. The last step is to derive the sufficient closed-loop stability conditions that the continuous-phase convergence rate should meet such that the divergence of Va caused by the uncontrolled impact is compensated by the continuous-phase convergence.

The convergence of the foot placement yst and impact timing Tk is proved based on Propositions 3 and 4 and the asymptotic convergence of the error state $x$. By Propositions 3 and 4, the deviations of the lateral foot placement and impact timing are bounded above by the norms of the actual state $x$ and the fictitious state $x̃$. Note that by definition, $x̃$ overlaps with $x$ within the given actual continuous phase. Thus, driving $x$ to zero will indirectly make $x̃$ diminish, which then eliminates the deviations $yst−ystd$ and $Tk−τk$ at the actual steady-state.▪

Remark 6. (Tuning continuous-phase convergence rate) By Theorem 1, the continuous-phase convergence rate of $x$ (or equivalently, Va) needs to be sufficiently fast for guaranteeing asymptotic trajectory tracking of the hybrid closed-loop system. The continuous-phase convergence rate of Va solely depends on that of V, because the stance foot is static during a continuous phase and $|yst−ystd|$ remains constant. We can construct V as $V=xTPx$, where $P$ is the solution to the Lyapunov equation [24]
$PA+ATP=−Q$

Here, $Q$ is any symmetric, positive-definite matrix satisfying $0<λQI≤Q$ with a positive number λQ. For simplicity, we can choose $Q$ as an identity matrix, and then λQ can be any number satisfying $0<λQ≤1$. Then, the bounds of V and $V˙$ in Eq. (13) become $c1=λmin(P), c2=λmax(P)$, and $c3=λQ$, where $λmin(P)$ and $λmax(P)$ are the smallest and the largest eigenvalues of $P$, respectively. Thus, the exponential convergence rate of V becomes $c3c2=λQλmax(P)$. Note that the value of the matrix $P$ depends on the PD gains, and thus $λmax(P)$ can be adjusted by tuning those gains. The full proof (Sec. 9.5) provides greater details about PD gain tuning. It also explains how to compute the lower bound of the convergence rate $c3c2$ for guaranteeing asymptotic error convergence of the hybrid closed-loop system.

Remark 7. (Satisfying lateral foot-placement condition) Theorem 1 indicates that the lateral foot-placement condition underlying the proposed impact invariance conditions in Propositions 1 and 2 is exactly met at the steady-state. Thus, the impact invariance of Z, which is the necessary condition for asymptotic trajectory tracking, is indeed satisfied at the steady-state; that is, if $x(τk−)→0$ then $x(τk+)=Δ(τk,0)→0$ as $k→∞$.

## 6 Simulations and Experiments

This section reports simulation and experiment results that demonstrate the global-position tracking performance of the proposed control approach.

The hardware platform used for controller validation is the OP3 bipedal humanoid robot developed by ROBOTIS Co., Ltd. (Seoul, South Korea) (Fig. 1). OP3 weighs 3.5 kg with a height of 0.51 m. It has twenty revolute joints comprising eight upper-body and 12 leg joints. As these joints (including ankles) are all independently actuated, the robot is fully actuated during a continuous phase of flat-foot walking without slippage.

### 6.1 Virtual Constraint Generation.

This section explains the lower-level, optimization-based trajectory generation of virtual constraints based on the proposed impact invariance conditions.

With full actuation, OP3's 12 leg joints can be directly commanded to track 12 independent desired trajectories, which are: (1) the desired global-position trajectory sd and (2) the desired functions yd and $ϕd$. As a higher-level planner supplies the desired global path on the walking surface and the desired position trajectory along the path, the objective of the trajectory generation is to plan the desired lateral base position yd and desired functions $ϕd$ that both define the virtual constraints.

Trajectory parameterization. The desired lateral base position yd is chosen as the following simple sinusoidal function to enable an oscillatory global motion about the centerline Γd during walking
$yd(x¯b):=α1 sin(α2x¯b+α3)$
(19)

with $α:=[α1 α2 α3]T∈ℝ3$ an unknown vector to be optimized.

The desired functions $ϕd$ are chosen as the desired trajectories for the following ten control variables $ϕc$:

• Height (zb) and roll, pitch, and yaw angles ($ψbroll, ψbpitch, ψbyaw$) of the base.

• Position (xsw, ysw, zsw) and roll, pitch, and yaw angles ($ψswroll, ψswpitch, ψswyaw$) of the swing foot.

This choice of control variables allows direct regulation of the poses (i.e., positions and orientations) of the trunk and swing foot to avoid overstretched leg joints, enforce a relatively steady trunk posture, and maintain a sufficient clearance between the swing foot and the walking surface.

The desired functions $ϕd(θ)$ are parameterized using Bézier curves [39]
$ϕd(θ):=∑k=0MakM!k!(M−k)!s(θ)k(1−s(θ))M−k$
(20)

where $M∈ℤ+$ is the order of the Bézier curves, $s(θ):=θ−θ+θ−−θ+$, $ak∈ℝ10$ is the unknown vector to be optimized, and $θ+$ and $θ−$ are the planned values of θ at the beginning and the end of a step, respectively. Recall that θ is chosen as the relative forward position of the base (Eq. (6)) and represents how far a step has progressed within a step.

Optimization formulation. The optimization variables are chosen as parameters $α$ in Eq. (19) and ak in Eq. (20). The constraints are set as:

• (B1) The proposed impact invariance conditions (A1)(A3) in Propositions 1 and 2.

• (B2) Feasibility constraints (e.g., joint-position limits, joint-torque limits, and ground-contact constraints).

• (B3) Gait parameters (e.g., step length and duration).

This list of constraints is not intended to be exhaustive as this study focuses on impact invariance construction and controller design instead of trajectory generation. matlab command fmincon is used to solve the optimization.

Desired trajectories. In the simulations and experiments, the planned virtual constraints are illustrated in Fig. 5. The centerline Γd of the desired path is the Xw-axis of the world reference frame. To test the capability of the proposed control approach in tracking desired position trajectories with constant or time-varying velocities, the following two desired position trajectories $sd(t)$ along Γd are considered:

Fig. 5
Fig. 5
Close modal
• $sd(t)=4.4t−3$ cm.

• $sd(t)=3.1t−1.5+1.5 sin(0.3t)−sin(0.8t)$ cm.

In theory, the proposed control approach can locally stabilize any profiles of $sd(t)$ that are differentiable in time. In practice, $sd(t)$ also needs to respect the robot's hardware constraints (e.g., actuation and kinematic limits).

### 6.2 Controller Implementation Procedure.

This subsection explains the experiment procedure that we adopt to implement the proposed controller on the physical OP3 robot using the ros package ($op3_manager$) developed by OP3's manufacturer.

Since the ros package does not support direct access to the output torques of joint motors, the proposed control law in Eq. (10), which is a torque command, cannot be directly implemented on OP3 and needs to be adapted for its implementation on the robot.

Considering that OP3's ros package allows users to send desired joint-position trajectories to individual joints and specify the PD gains of OP3's default joint controller, we adopt the following controller implementation procedure [21]: (a) to generate the desired position trajectories of individual joints, $qd(t)$ and (b) to send the desired trajectories to the default joint-position controller. The main steps of this procedure are shown in Fig. 6.

Fig. 6
Fig. 6
Close modal

Although the adapted controller directly tracks the individual joint trajectories qd instead of the original Cartesian-space trajectories hd, the controller implementation procedure still allows satisfactory tracking of hd. This is because $qd$ preserve the dynamic feasibility and desired features of hd as specified in B1–B3.

### 6.3 Simulation and Experiment Setup.

This subsection reports the setup of MATLAB and WEBOTS simulations and hardware experiments for controller validation.

MATLAB. To validate the theoretical controller design, we utilize matlab to implement the control law based on the full-order model of OP3 (Eq. (4)). The control gains are set as $KP=225·I$ and $KD=30·I$ to ensure the matrix $A$ is Hurwitz. The simulation results are shown in Figs. 7 and 8.

Fig. 7
Fig. 7
Close modal
Fig. 8
Fig. 8
Close modal

WEBOTS. To gain preliminary insights into the effectiveness of the proposed controller implementation procedure as explained in Sec. 6.2, we use webots to simulate a 3D realistic biped model that closely emulates OP3's graphical, physical, and dynamical properties (including its limited actuator accessibility). The control gains that the emulated robot system allows users to tune are the effective PD gains, whose physical meaning is different from KP and KD in Eq. (11). These effective gains are tuned to be “10” and “0” such that the resulting tracking performance is comparable with the matlab results. webots simulation results of the adapted controller are displayed in Figs. 8 and 9.

Fig. 9
Fig. 9
Close modal

Experiments. The experiment setup is shown in Fig. 10. With this setup, the robot's joint angles can be directly measured by joint encoders, and its global pose can be determined by: (a) using the 4 K prowebcam and apriltag [40] to obtain the stance-foot pose in the world reference frame and (b) using the obtained stance-foot pose to solve for the robot's global pose via forward kinematics. By providing relatively accurate measurement, the use of the overhead camera and apriltag allows us to focus on controller validation. The experiment is guided by the controller adaptation procedure from Sec. 6.2. The initial tracking error of the desired position trajectory sd is 3 cm, which is approximately 1/3 of a nominal step length. The initial path tracking error is 5 cm. Similar to the gain tuning in webots, the effective PD gains are, respectively, tuned to be “800” and “0” to ensure a relatively fast error convergence without violating the actuator's torque limit. Experiment results of OP3 walking on a concrete and a relatively slippery ceramic floor are shown in Fig. 11. Videos of the experiments can be accessed at following link.2

Fig. 10
Fig. 10
Close modal
Fig. 11
Fig. 11
Close modal

### 6.4 Discussions on Validation Results.

This subsection provides discussions on the controller evaluation results.

Tracking accuracy in simulations. The virtual-constraint tracking results in Figs. 7 (matlab) and 9 (webots) show that the proposed control law is capable of accurately enforcing the virtual constraints for 3D bipeds during fully actuated walking. Note that the tracking errors observed in the webots simulations are larger than the matlab results in the swing foot's lateral position ysw, base height zb, base yaw angles $ψbyaw$, and the swing foot's orientation ($ψswroll, ψswpitch, ψswyaw$). This is partly caused by the foot slippage of the robot during webots simulations, which is not captured by the robot model used in matlab simulations. The global-position tracking results in Fig. 8 validate that the proposed control law drives the robot to asymptotically converge to the desired global-position trajectory sd while moving along the centerline Γd of the global path. In particular, the accurate tracking results obtained in webots indicate the effectiveness of the proposed controller implementation procedure in guaranteeing reliable trajectory tracking in the presence of hardware limitations.

Tracking accuracy in experiments. As illustrated in Fig. 11(a) (top), under the proposed global-position tracking (GPT) controller, the robot's actual global position xb (labeled as “xb (GPT)”) converges to a relatively small neighborhood about its desired trajectory sd within 3 s when the robot walks on a concrete floor. Also, Fig. 11(a) (bottom) illustrates that despite an initial path tracking error of 5 cm, the robot remains close to the centerline Γd of the desired global path, as indicated by the footstep trajectories labeled as “foot placement (GPT).” Due to uncertainties such as hardware limitations, modeling errors, and floor surface irregularity, achieving an exactly zero steady-state tracking error on a physical robot may not be feasible. Thanks to the inherent robustness of feedback control, the proposed control approach achieves a small steady-state tracking error, although uncertainties are not explicitly considered in the controller design.

Robustness. To further test the limit of the inherent robustness of the proposed control approach, experiments of OP3 walking on a ceramic tile floor were conducted (Fig. 11(b)). As the surface of the ceramic tiles is relatively more slippery than the concrete floor, the robot's stance foot slips more frequently on the tile floor, causing a stronger violation of the modeling assumption of static stance foot. Yet, a relatively small global-position tracking error is still realized when the initial foot placement error is small, as shown in Fig. 11(b).

Necessity of global-position tracking control. To illustrate the need to explicitly address global-position tracking in controller design, a global-velocity tracking (GVT) controller, which is analogous to the orbitally stabilizing controller for fully actuated walking [21], is implemented. Its global-position tracking performance is displayed in Figs. 11. The figure shows that the GVT controller achieves accurate tracking of the desired global velocity $s˙d$. Yet, the global-position tracking performance is not satisfactorily guaranteed, as indicated by the relatively large deviations of the global position (labeled as “xb (GVT)”) and the footstep trajectories (labeled as “foot placement (GVT)”). Thus, this experiment result indicates the necessity to explicitly handle global-position tracking through controller design.

## 7 Discussion

This study has extended the previous method of impact invariance construction from orbital stabilization [14] to the stabilization of time-varying global-position trajectory for 3D fully actuated robots. The proposed method produces impact invariance conditions that can be imposed in the trajectory generation of virtual constraints for ensuring their agreement with impact dynamics. Moreover, although the impact maps of the virtual constraints and global trajectory are generally nonlinearly coupled through the robot's kinematic chains, these conditions can automatically ensure any arbitrary smooth desired global-position trajectory respects the impact dynamics. Indeed, as shown in Fig. 8, the proposed controller achieves asymptotic tracking of two different global-position trajectories under the same virtual constraints, indicating that the virtual constraints ensure the impact agreement for different desired global-position trajectories. Thus, the proposed impact invariance conditions can allow the decoupling between the lower-level trajectory generation of virtual constraints and the higher-level planning of global-position trajectory. The decoupling could permit offline planning of virtual constraints, thus reducing the computational load for online planning.

This study has also introduced the Lyapunov-based stability conditions for the hybrid closed-loop error system associated with 3D bipedal robots during fully actuated walking. Controller designs satisfying these conditions can accurately track the time-varying desired global-position trajectory, as demonstrated in Figs. 8 and 11. The proposed control approach can also indirectly drive the lateral foot placement yst to the desired location ystd, which is predicted by the asymptotic convergence of the Lyapunov function Va that explicitly contains the lateral foot-placement error. Note that our previous controller for 2D walking cannot address the convergence of $yst−ystd$ as it does not consider the robot's lateral movement. The capability of accurate foot placement could potentially be exploited to handle locomotion on discrete terrains (e.g., stepping stones [41]).

In theory, for the proposed control law to achieve zero tracking error, the desired trajectory needs to respect the proposed impact invariance conditions. Yet, in practice, achieving exactly zero final tracking error may not be necessary. Rather, achieving a final error within a reasonably small bound could be sufficiently satisfactory for practical applications. To this end, the proposed impact invariance conditions can be theoretically relaxed to allow bounded violation of the condition for stabilizing the origin of the tracking error system in the sense of Lyapunov stability rather than asymptotic stability. Specifically, we can incorporate the impact invariance conditions as an inequality constraint in trajectory generation instead of an equality constraint.

The proposed control approach, including continuous input–output linearizing control design and the impact invariance construction, builds upon a hybrid robot model that holds under several modeling simplifying assumptions. These assumptions are reasonable for robot walking in relatively structured environments. Yet, they may not hold if the environments are more complex, which will in turn affect the controller performance. Indeed, as the experiment results in Fig. 11 illustrate, when the floor is relatively slippery, the assumption that the foot and surface have a secured contact no longer holds, which is not explicitly addressed by the proposed control law. For this reason, these experiment results show a slower convergence rate and larger tracking error compared with the simulation results in Fig. 8. To this end, adapting the proposed controller to more complex environments is necessary. For instance, we could cast the proposed control approach into a quadratic program [42] for ensuring the feasibility of ground contact forces in the presence of foot slippage induced by surface irregularity. Another potential approach is to integrate the proposed control law with adaptive and robust control action [26,43,44] for enabling online model estimation and better disturbance rejection.

The proposed controller is built upon a fully actuated robot model, and is thus effective when the robot is fully actuated. For a bipedal robot to be fully actuated, its motion needs to satisfy certain necessary constraints. For instance, a bipedal robot with finite size feet (e.g., OP3) is fully actuated when its support foot is in a static, full contact with the ground during walking; that is, its motion satisfies the ground contact constraints (e.g., friction cone, unilateral, and center of pressure constraints). Thus, the planned motion (defined by the virtual constraints and the desired global trajectory) should meet these constraints with a reasonable margin. Here the margin ensures that the actual walking motion also satisfy those constraints when it is near the planned motion, so that the controller will be effective in driving the actual motion to the desired one. Note that the proposed controller is not intended for robots with limited-size feet (e.g., point feet) or passive ankles, which commonly adopt underactuated gait due to the lack of control authority compared with the robot's degrees-of-freedom. Still, the proposed control method could be extended to multidomain walking gait, which comprises subphases of full actuation, underactuation, and over actuation, as our preliminary theoretical and simulation studies indicate [29]. The key to this extension is the adaptation of the proposed stability conditions from single-mode to multimode hybrid robot dynamics.

## 8 Conclusion and Future Work

This paper has introduced a control approach that explicitly addresses the hybrid dynamics of 3D bipedal robots for achieving asymptotic global-position tracking during fully actuated walking. With the output function designed as the tracking error of the desired global-position trajectory and virtual constraints, a continuous input–output linearizing control law was synthesized to asymptotically drive the output function to zero within continuous phases. Impact invariance conditions were derived to guide the generation of virtual constraints such that the robot's desired motions defined by the virtual constraints and the desired global-position trajectory all respect the discrete landing impact dynamics. Sufficient conditions were derived based on Lyapunov theory under which the proposed continuous control law provably guarantees the asymptotic tracking performance of the hybrid closed-loop system. Simulation and experiment results demonstrated the effectiveness of the proposed control approach in realizing satisfactory global-position tracking.

Our future work will apply and extend the proposed approach from straight-line to curved-path locomotion as real-world applications of legged robots commonly require walking in varying directions. To enable efficient planning, we will construct a library [20] of virtual constraints offline that corresponds to a common range of direction-varying gait parameters and interpolate the virtual constraints online to fit the varying walking directions along a curved path.

## Acknowledgment

The authors would like to thank the anonymous reviewers for their insightful and constructive comments.

Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

## Funding Data

• Directorate for Engineering, National Science Foundation (Grant No. CMMI-1934280; Funder ID: 10.13039/100000084).

### Appendix: Proofs of All Propositions and Theorems

#### A.1 Proof of Proposition 1

Proof 1. With the pre-impact joint position $q*$, the postimpact joint position and phase variable are $q(τk+)=Δq(q*)=q0$ and $θ(τk+)=x¯b(q(τk+))=x¯b(q0)=θ0$, respectively.

Under the foot-placement condition yst = ystd and condition (A1), the postimpact value of the output function $y¯b−(yd−yst)$ is $y¯b(q0)−yd(θ0)+yst(τk+)=y¯b(q0)−yd(θ0)+ystd=0.$ Similarly, given condition (A1), the postimpact value of $ϕc−ϕd$ is $ϕc(q0)−ϕd(θ0)=0$.

Since $y˙(τk−)=0$, we have
$y˙(τk−)=h˙c(q*,q˙(τk−))−h˙d(t,θ*,θ˙(τk−))=Jhc(q*)q˙(τk−)−[s˙d(τk−)∂yd∂θ(θ*)θ˙(τk−)∂ϕd∂θ(θ*)θ˙(τk−)]=0$
(A1)
and $θ˙(τk−)=∂x¯b∂q(q*)q˙(τk−)=s˙d(τk−)−x˙st(τk−)=s˙d(τk−),$ where $Jhc(q*)=∂hc∂q(q*)$. Thus, the pre-impact joint velocity is
$q˙(τk−)=Jhc−1(q*)[1∂yd∂θ(θ*)∂ϕd∂θ(θ*)]θ˙(τk−)$
(A2)
Note that $q˙(τk+)=Δq˙(q*)q˙(τk−)$ and $θ˙(τk+)=∂x¯b∂q(q0)q˙(τk+)$. Then, by condition (A2), the postimpact value of $[y¯˙b−(y˙d−y˙st)ϕ˙c−ϕ˙d]$ becomes
$[∂y¯b∂q(q0)∂ϕc∂q(q0)]q˙(τk+)−[∂yd∂θ(θ0)∂ϕd∂θ(θ0)]θ˙(τk+)=0$

Thus, the impact invariance of $Z̃$ is met under conditions (A1) and (A2) from Proposition 1.▪

#### A.2 Proof of Proposition 2

Proof 2. Because xb and $sd(t)$ are both continuous in t, we obtain
$x¯b(q0)−sd(τk+)+xst=xb(τk+)−sd(τk+)=xb(τk−)−sd(τk−)$

Thus, if $xb(τk−)−sd(τk−)=0$, then $x¯b(q0)−sd(τk+)+xst=0$ automatically holds.

Because the stance foot remains static just before and after the impact, $x˙st(τk+)=x˙st(τk−)=0$ holds. Also, as the desired global velocity $s˙d$ is continuous in t, we have $s˙d(τk+)=s˙d(τk−)$.

From the proof of Proposition 1, we have
$θ˙(τk+)=∂x¯b∂q(q0)Δq˙(q*)Jh−1(q*)[1∂yd∂θ(θ*)∂ϕd∂θ(θ*)]θ˙(τk−)$
(A3)

which yields $θ˙(τk+)=θ˙(τk−)$ under condition (A3). Then, the postimpact value of the first time derivative of the output function $xb−sd$ becomes $x¯˙b(q0,q˙(τk+))−(s˙d(τk+)−x˙st+)=θ˙(τk+)−s˙d(τk+)+x˙st+=θ˙(τk−)−s˙d(τk−)+x˙st−$, which is zero if $θ˙(τk−)−s˙d(τk−)+x˙st−=0$.

Thus, under conditions (A1)(A3) from Propositions 1 and 2, the impact invariance of Z holds.▪

#### A.3 Proof of Proposition 3

Proof 3. Because the output function state $y$ and $y˙$ and the swing-foot height zsw defining the switching surface S are both continuously differentiable in their respective arguments, the function defining the switching surface Sx is continuously differentiable in its argument [28]. Also, note that the continuous-phase vector field (i.e., $Ax$) of the error state $x$ is continuously differentiable in $x$.

Then, by Lemma 2.1 and Corollary 2.4 in Ref. [28], the impact timing Tk is an implicit function of the state $x$, and is Lipschitz continuous with respect to $x$. Thus, there exists a positive number r1 and a Lipschitz constant $LTx$ such that $|Tk−τk|≤LTx||x̃(τk;Tk−1+,x|k−1+)||$ for any $x|0+∈Br1(0)$ and any $k∈ℤ+$.

#### A.4 Proof of Proposition 4

Proof 4. Let $ϕsw,y(θ)$ denote the desired trajectory of the control variable $ysw$. Because the stance-foot position during the $(k+1)th$ step is the swing-foot position at the end of the kth step, one has $yst|k+=ysw|k−$ and $ϕsw,y(θ*)=ystd$.

Let $θ̃(t;Tk−1+,x|k−1+)$ be the phase variable associated with the fictitious state $x̃(t;Tk−1+,x|k−1+)$ that overlaps with the actual state $x$ on $t∈[Tk−1+,Tk−]$. Then, by the triangular inequality, we can approximate the upper bound of the absolute lateral foot-placement error as
$|yst|k+−ystd|=|ysw|k−−ϕsw,y(θ*)|≤ |ysw|k−−ϕsw,y(θ(Tk−))|+|ϕsw,y(θ(Tk−))−ϕsw,y(θ̃(τk;Tk−1+,x|k−1+))|+|ϕsw,y(θ̃(τk;Tk−1+,x|k−1+))−ϕsw,y(θ*)|$
(A4)

The upper bounds of the three terms on the right-hand side of this inequality are derived next.

As $ysw−ϕsw,y$ is an element of the full error state $x$, its norm satisfies
$|ysw|k−−ϕsw,y(θ(Tk−))|≤||x|k−||$
(A5)
Because $θ(Tk−)=θ̃(Tk−;Tk−1+,x|k−1+)$ and because $ϕsw,y(θ)$ and $θ̃(t;Tk−1+,x|k−1+)$ are continuously differentiable in θ and t, respectively, there exists a positive number r2 and Lipschitz constants $Lϕsw,y$ and $Lθt$ such that
$||ϕsw,y(θ(Tk−))−ϕsw,y(θ̃(τk;Tk−1+,x|k−1+))||≤ Lϕsw,y||θ(Tk−)−θ̃(τk;Tk−1+,x|k−1+)||= Lϕsw,y||θ̃(Tk−;Tk−1+,x|k−1+)−θ̃(τk;Tk−1+,x|k−1+)||≤ Lϕsw,yLθt|Tk−τk|$
(A6)

for any $x|0+∈Br2(0):={x∈ℝ2n:||x||≤r2}$.

With $θ*=θ(q*)=θ̃(τk−;Tk−1+,0)=sd(τk−)−xst(τk−),$ we have
$|θ̃(τk;Tk−1+,x|k−1+)−θ*|=||θ̃(τk;Tk−1+,x|k−1+)−(sd(τk−)−xst(τk−))||≤ ||x̃(τk;Tk−1+,x|k−1+)||$
(A7)
Then
$|ϕsw,y(θ̃(τk;Tk−1+,x|k−1+))−ϕsw,y(θ*)|≤Lϕsw,y||θ̃(τk;Tk−1+,x|k−1+)−θ*||≤Lϕsw,y||x̃(τk;Tk−1+,x|k−1+)||$
(A8)

Let $βst:=Lϕsw,y(LθtLTx+1)$ and $d1:=min(r1,r2)$. From Proposition 3 and Eqs. (A4)(A8), we obtain $|yst|k+−ystd|≤||x|k−||+βst||x̃(τk;Tk−1+,x|k−1+)||$ for any $x|0+∈Bd1(0)$ and $k∈ℤ+$.▪

#### A.5 Proof of Theorem 1

Proof 5. As explained in Sec. 5.2, the Lyapunov function candidate is defined as $Va(x,yst−ystd):=V(x)+σ(yst−ystd)2,$ where σ is a positive number to be specified in this proof.

By the stability theory based on the construction of multiple Lyapunov functions [36], the origin of the hybrid time-varying system in Eq. (12) is locally asymptotically stable if there exists a positive number d2 such that for any $x|0+∈Bd2(0)$, Va is monotonically decreasing within each continuous phase, and ${Va|1+,Va|2+,Va|3+…}$ is a strictly decreasing sequence with $Va|k+→0$ as $k→∞$.▪

Evolution of Va during continuous phases. With the PD gains chosen such that $A$ is Hurwitz, Eq. (13) gives $V|k−≤e−c3c2(Tk+1−Tk)V|k−1+$ within the kth ($k∈ℤ+$) continuous phase. Since $yst−ystd$ remains constant within the step due to the static stance foot, Va monotonically decreases within the kth phase.

Evolution of Va across nonlinear impact maps. Consider the foot-landing event at the end of the kth walking step (i.e., $t=Tk−$). The tracking error expansion across the landing event is analyzed as follows.

Because the desired functions yd and $ϕd(θ)$ satisfy the conditions (B1)–(B3), the impact invariance associated with the error state $x$ holds, which leads to $Δ(τk−,0,ystd)=0$. Then, the value of $x$ just after the landing can be approximated by applying the triangular inequality as
$||x|k+||=||Δ(Tk−,x|k−,yst|k−)||=||Δ(Tk−,x|k−,yst|k−)−Δ(τk−,0,ystd)||≤||Δ(Tk−,x|k−,yst|k−)−Δ(τk−,x|k−,yst|k−)||+||Δ(τk−,x|k−,yst|k−)−Δ(τk−,0,yst|k−)||+||Δ(τk−,0,yst|k−)−Δ(τk−,0,ystd)||+||Δ(τk−,0,ystd)||$
(A9)
As the reset map $Δ(t,x,yst)$ is continuously differentiable in t, $x$, and yst, there exists a positive number r3 and Lipschitz constants $LΔt, LΔx$, and $LΔst$ such that the following inequalities hold for any $x|0+∈Br3(0)$:
$||Δ(Tk−,x|k−,yst|k−)−Δ(τk−,x|k−,yst|k−)||≤LΔt|Tk−τk|||Δ(τk−,x|k−,yst|k−)−Δ(τk−,0,yst|k−)||≤LΔx||x|k−||||Δ(τk−,0,yst|k−)−Δ(τk−,0,ystd)||≤LΔst|yst|k−−ystd|$
(A10)
From Eqs. (16) and (A10), we obtain
$||Δ(Tk−,x|k−,yst|k−)−Δ(τk−,x|k−,yst|k−)||≤ LΔtLTx||x̃(τk;Tk−1+,x|k−1+)||$
(A11)
for any $x|0+∈Bd2(0)$, where $d2=min{d1,r3}$. From Eqs. (A9)(A11), we have for any $x|0+∈Bd2(0)$
$||x|k+||=||Δ(Tk−,x|k−,yst|k−)||≤LΔx||x|k−||+LΔtLTx||x̃(τk;Tk−1+,x|k−1+)||+LΔst||yst|k−−ystd||$
(A12)
The upper bounds of $||x|k−||$ and $||x̃(τk;Tk−1+,x|k−1+)||$ with respect to the tracking error norm $||x|k−1−||$ can be derived based on Eq. (13) as
$||x|k−||≤c2c1e−c32c2(Tk−Tk−1)||x|k−1+|| and ||x̃(τk;Tk−1+,x|k−1+)||≤c2c1e−c32c2(τk−Tk−1)||x|k−1+||$
(A13)
Then, from Eqs. (A12) and (A13), the postimpact error norm can be approximated as
$||x|k+||≤c2c1(LΔtLTx+LΔxe−c32c2(Tk−τk))e−c32c2(τk−Tk−1)||x|k−1+||+ LΔst|yst|k−−ystd|$
(A14)
For any $ε>0$ there exist PD gains corresponding to a sufficiently high convergence rate $c32c2$ such that $e−c32c2(Tk−τk)≤1+ε$. Then, the approximation of the postimpact error norm can be simplified into
$||x|k+||≤αx||x|k−1+|| +αst|yst|k−−ystd|$
(A15)

where $αx:=c2c1(LΔtLTx+LΔx(1+ε))e−c32c2Δτk, Δτk:=τk−Tk−1$, and $αst:=LΔst$.

Now, we derive the upper bound of $|yst|k−−ystd|$ with respect to the tracking error norm $||x|k−1−||$. Because the stance foot remains static within a step, we have $yst|k−=yst|k−1+.$ Then, from Eq. (17)
$|yst|k+−ystd|≤||x|k−||+βst||x̃(τk;Tk−1+,x|k−1+)||≤γx||x|k−1+||$
(A16)

holds, where $γx:=c2c1(βst+(1+ϵ))e−c32c2Δτk$.

Finally, combining Eqs. (13), (A15), and (A16) provides the following approximation of the postimpact value of the Lyapunov function Va:
$Va|k+=V|k++σ(yst|k+−ystd)2≤c2||x|k+||2+σ(yst|k+−ystd)2≤B(c1||x|k−1+||2+σ(yst|k−1+−ystd)2)≤B(V|k−1++σ(yst|k−1+−ystd)2)≤BVa|k−1+$

where $B:=max(2c2αx2+σγx2c1,2c2αstσ)$.

Evolution of Va for the hybrid model. If the PD gains and σ are chosen such that
$2c2αx2+σγx2c1<1 and 2c2αstσ<1$
(A17)

hold (i.e., B < 1), then for any $x|0+∈Bd2(0)$, the sequence ${Va|1+, Va|2+, Va|3+…}$ is strictly decreasing with $Va|k+→0$ as $k→∞$. Thus, the closed-loop hybrid system is locally asymptotically stable if the PD gains ensure that the matrix $A$ is Hurwitz and that Eq. (A17) holds for any $x|0+∈Bd2(0)$.

To meet the two inequality conditions in Eq. (A17), we can choose the function $V(x)$ to be $V(x)=xTPx$ as explained in Remark 6. This choice results in the continuous-phase convergence rate of $V(x)$ as $c3c2=λQλmax(P)$, which can be tuned with the PD gains. Specifically, to satisfy the second inequality in Eq. (A17), we can specify σ as any positive number such that $σ>2c2αst=2λmax(P)αst$, where αst can be estimated from system dynamics. For instance, we can choose σ to be $2kσλmax(P)αst$ with any constant $kσ>1$. Then, we can tune the PD gains to meet the first inequality in Eq. (A17), by allowing a sufficiently high continuous-phase convergence rate that leads to sufficiently small values of αx and γx for satisfying $αx2+γx2≤c12λmax(P)max(1,kσαst)$.

Convergence of impact timings. When the state $x$ reaches zero at the steady-state, from Eq. (A13), the fictitious state satisfies $||x̃(τk;Tk−1+,x|k−1+)||≤c2c1e−c32c2(τk−Tk−1)||x|k−1+||→0$ as $k→∞$. Then, by Eq. (16), $|Tk−τk|≤LTx||x̃(τk;Tk−1+,x|k−1+)||→0$ as $k→∞$.

Convergence of lateral foot placement. By the definition of Va in Eq. (18), $Va(x,yst−ystd)=V(x)+σ(yst−ystd)2$, where σ is positive and $V(x)$ and $(yst−ystd)2$ are all bounded and non-negative. Thus, if $Va→0$ as $t→∞$, then $(yst−ystd)2→0$ as $t→∞$; that is, $yst→ystd$ as $t→∞$.

## References

1.
Vukobratović
,
M.
,
Borovac
,
B.
, and
Šurdilović
,
D.
,
2001
, “
Zero Moment Point-Proper Interpretation and New Applications
,”
Proceedings of IEEE International Conference on Humanoid Robots
, Tokyo, Japan, Nov. 22–24, pp.
237
244
.
2.
Kajita
,
S.
,
Kanehiro
,
F.
,
Kaneko
,
K.
,
Fujiwara
,
K.
,
,
K.
,
Yokoi
,
K.
, and
Hirukawa
,
H.
,
2003
, “
Biped Walking Pattern Generation by Using Preview Control of Zero-Moment Point
,”
Proceedings of IEEE International Conference Robotics Automation
, Taipei, Taiwan, Sept. 14–19, pp.
1620
1626
.10.1109/ROBOT.2003.1241826
3.
Kim
,
J.-Y.
,
Park
,
I.-W.
, and
Oh
,
J.-H.
,
2006
, “
Experimental Realization of Dynamic Walking of the Biped Humanoid Robot KHR-2 Using Zero Moment Point Feedback and Inertial Measurement
,”
,
20
(
6
), pp.
707
736
.10.1163/156855306777361622
4.
Golliday
,
C.
, and
Hemami
,
H.
,
1977
, “
An Approach to Analyzing Biped Locomotion Dynamics and Designing Robot Locomotion Controls
,”
IEEE Trans. Automat. Contr.
,
22
(
6
), pp.
963
972
.10.1109/TAC.1977.1101650
5.
Hürmüzlü
,
Y.
, and
Moskowitz
,
G. D.
,
1986
, “
The Role of Impact in the Stability of Bipedal Locomotion
,”
Dyn. Stabil. Syst.
,
1
(
3
), pp.
217
234
.10.1080/02681118608806015
6.
Bhounsule
,
P. A.
,
Zamani
,
A.
, and
Pusey
,
J.
,
2018
, “
Switching Between Limit Cycles in a Model of Running Using Exponentially Stabilizing Discrete Control Lyapunov Function
,”
Proceedings of American Control Conference
, Milwaukee, WI, June 27–29, pp.
3714
3719
.10.23919/ACC.2018.8431123
7.
Yeatman
,
M.
,
Lv
,
G.
, and
Gregg
,
R. D.
,
2019
, “
Decentralized Passivity-Based Control With a Generalized Energy Storage Function for Robust Biped Locomotion
,”
ASME J. Dyn. Syst., Meas., Control
,
141
(
10
), p.
101007
.10.1115/1.4043801
8.
Gao
,
Y.
,
Da
,
X.
, and
Gu
,
Y.
,
2020
, “
Impact-Aware Online Motion Planning for Fully-Actuated Bipedal Robot Walking
,”
Proceedings of American Control Conference
, Denver, CO, July 1–3, pp.
2100
2105
.10.23919/ACC45564.2020.9147498
9.
Gu
,
Y.
,
Yao
,
B.
, and
Lee
,
C.
,
2017
, “
Time-Dependent Orbital Stabilization of Underactuated Bipedal Walking
,”
Proceedings of American Control Conference
, Seattle, WA, May 24–26, pp.
4858
4863
.10.23919/ACC.2017.7963707
10.
Rijnen
,
M.
,
Biemond
,
J. B.
,
van de Wouw
,
N.
,
Saccon
,
A.
, and
Nijmeijer
,
H.
,
2020
, “
Hybrid Systems With State-Triggered Jumps: Sensitivity-Based Stability Analysis With Application to Trajectory Tracking
,”
IEEE Trans. Automat. Contr.
,
65
(
11
), pp.
4568
4583
.10.1109/TAC.2019.2961996
11.
Rijnen
,
M.
,
Chen
,
H. L.
,
van de Wouw
,
N.
,
Saccon
,
A.
, and
Nijmeijer
,
H.
,
2019
, “
Sensitivity Analysis for Trajectories of Nonsmooth Mechanical Systems With Simultaneous Impacts: A Hybrid Systems Perspective
,”
Proceedings of American Control Conference
, Philadelphia, PA, July 10–12, pp.
3623
3629
.10.23919/ACC.2019.8814388
12.
Wang
,
Y.
,
Dehio
,
N.
,
Tanguy
,
A.
, and
Kheddar
,
A.
,
2020
, “
,” arXiv preprint arXiv:2006.01987.
13.
Grizzle
,
J.
,
Abba
,
G.
, and
Plestan
,
P.
,
2001
, “
Asymptotically Stable Walking for Biped Robots: Analysis Via Systems With Impulse Effects
,”
IEEE Trans. Automat. Contr.
,
46
(
1
), pp.
51
64
.10.1109/9.898695
14.
Westervelt
,
E. R.
,
Grizzle
,
J. W.
, and
Koditschek
,
D. E.
,
2003
, “
Hybrid Zero Dynamics of Planar Biped Walkers
,”
IEEE Trans. Automat. Contr.
,
48
(
1
), pp.
42
56
.10.1109/TAC.2002.806653
15.
Morris
,
B.
, and
Grizzle
,
J. W.
,
2009
, “
Hybrid Invariant Manifolds in Systems With Impulse Effects With Application to Periodic Locomotion in Bipedal Robots
,”
IEEE Trans. Automat. Contr.
,
54
(
8
), pp.
1751
1764
.10.1109/TAC.2009.2024563
16.
Martin
,
A. E.
, and
Gregg
,
R. D.
,
2015
, “
Hybrid Invariance and Stability of a Feedback Linearizing Controller for Powered Prostheses
,”
Proceedings of American Control Conference
, Chicago, IL, July 1–3, pp.
4670
4676
.10.1109/ACC.2015.7172065
17.
Gong
,
Y.
,
Hartley
,
R.
,
Da
,
X.
,
Hereid
,
A.
,
Harib
,
O.
,
Huang
,
J.-K.
, and
Grizzle
,
J.
,
2019
, “
Feedback Control of a Cassie Bipedal Robot: Walking, Standing, and Riding a Segway
,”
Proceedings of American Control Conference
, Philadelphia, PA, July 10–12, pp.
4559
4566
.10.23919/ACC.2019.8814833
18.
Fevre
,
M.
,
Goodwine
,
B.
, and
Schmiedeler
,
J. P.
,
2019
, “
Terrain-Blind Walking of Planar Underactuated Bipeds Via Velocity Decomposition-Enhanced Control
,”
Int. J. Robot. Res.
,
38
(
10–11
), pp.
1307
1323
.10.1177/0278364919870242
19.
Hamed
,
K. A.
, and
Grizzle
,
J. W.
,
2014
, “
Event-Based Stabilization of Periodic Orbits for Underactuated 3-D Bipedal Robots With Left-Right Symmetry
,”
IEEE Trans. Robot.
,
30
(
2
), pp.
365
381
.10.1109/TRO.2013.2287831
20.
Da
,
X.
,
Harib
,
O.
,
Hartley
,
R.
,
Griffin
,
B.
, and
Grizzle
,
J. W.
,
2016
, “
From 2D Design of Underactuated Bipedal Gaits to 3D Implementation: Walking With Speed Tracking
,”
IEEE Access
,
4
, pp.
3469
3478
.10.1109/ACCESS.2016.2582731
21.
Ames
,
A. D.
,
Cousineau
,
E. A.
, and
Powell
,
M. J.
,
2012
, “
Dynamically Stable Bipedal Robotic Walking With NAO Via Human-Inspired Hybrid Zero Dynamics
,”
Proceedings of ACM International Conference on Hybrid Systems: Computation and Control
, Beijing China, Apr. 17–19, pp.
135
144
.https://dl.acm.org/doi/abs/10.1145/2185632.2185655
22.
Hamed
,
K.
,
Safaee
,
B.
, and
Gregg
,
R. D.
,
2019
, “
Dynamic Output Controllers for Exponential Stabilization of Periodic Orbits for Multidomain Hybrid Models of Robotic Locomotion
,”
ASME J. Dyn. Syst. Meas. Control
,
141
(
12
), p.
121011
.10.1115/1.4044618
23.
Xiong
,
X.
,
Reher
,
J.
, and
Ames
,
A. D.
,
2021
, “
Global Position Control on Underactuated Bipedal Robots: Step-to-Step Dynamics Approximation for Step Planning
,”
Proceeding of IEEE International Conference Robotics Automation
, Xi'an, China, May 30–June 5, pp.
2825
2831
.https://www.researchgate.net/publication/345788665_Global_Position_Control_on_Underactuated_Bipedal_Robots_Stepto-step_Dynamics_Approximation_for_Step_Planning
24.
Khalil
,
H. K.
,
1996
,
Nonlinear Control
,
Prentice Hall
, Hoboken, NJ.
25.
Gu
,
Y.
, and
Yuan
,
C.
,
2020
, “
Adaptive Robust Trajectory Tracking Control of Fully Actuated Bipedal Robotic Walking
,”
Proceedings of IEEE/ASME International Conference on Advanced Intelligent Mechatronics
, Boston, MA, July 6–9, pp.
1310
1315
.10.1109/AIM43001.2020.9158814
26.
Gu
,
Y.
, and
Yuan
,
C.
,
2021
, “
Adaptive Robust Tracking Control for Hybrid Models of Three-Dimensional Bipedal Robotic Walking Under Uncertainties
,”
ASME J. Dyn. Syst., Meas., Control
,
143
(
8
), p.
081007
.10.1115/1.4050259
27.
Gu
,
Y.
,
Yao
,
B.
, and
Lee
,
C. S. G.
,
2016
, “
Bipedal Gait Recharacterization and Walking Encoding Generalization for Stable Dynamic Walking
,”
Proceedings of IEEE International Conference on Robotics Automation
, Stockholm, Sweden, May 16–21, pp.
1788
1793
.https://dl.acm.org/doi/abs/10.1109/ICRA.2016.7487324
28.
Gu
,
Y.
,
Yao
,
B.
, and
Lee
,
C. S. G.
,
2018
, “
Exponential Stabilization of Fully Actuated Planar Bipedal Robotic Walking With Global Position Tracking Capabilities
,”
ASME J. Dyn. Syst. Meas. Control
,
140
(
5
), p.
051008
.10.1115/1.4038268
29.
Gao
,
Y.
, and
Gu
,
Y.
,
2019
, “
Global-Position Tracking Control of Multi-Domain Planar Bipedal Robotic Walking
,”
ASME
Paper No. DSCC2019-9117.10.1115/DSCC2019-9117
30.
Menini
,
L.
, and
Tornambè
,
A.
,
2001
, “
Asymptotic Tracking of Periodic Trajectories for a Simple Mechanical System Subject to Nonsmooth Impacts
,”
IEEE Trans. Automat. Control
,
46
(
7
), pp.
1122
1126
.10.1109/9.935068
31.
Biemond
,
J. J. B.
,
van de Wouw
,
N.
,
Heemels
,
W.
, and
Nijmeijer
,
H.
,
2013
, “
Tracking Control for Hybrid Systems With State-Triggered Jumps
,”
IEEE Trans. Automat. Control
,
58
(
4
), pp.
876
890
.10.1109/TAC.2012.2223351
32.
Forni
,
F.
,
Teel
,
A. R.
, and
Zaccarian
,
L.
,
2013
, “
Follow the Bouncing Ball: Global Results on Tracking and State Estimation With Impacts
,”
IEEE Trans. Automat. Control
,
58
(
6
), pp.
1470
1485
.10.1109/TAC.2013.2237952
33.
Naldi
,
R.
, and
Sanfelice
,
R. G.
,
2013
, “
Passivity-Based Control for Hybrid Systems With Applications to Mechanical Systems Exhibiting Impacts
,”
Automatica
,
49
(
5
), pp.
1104
1116
.10.1016/j.automatica.2013.01.018
34.
Rijnen
,
M.
,
van Rijn
,
A. T.
,
Dallali
,
H.
,
Saccon
,
A.
, and
Nijmeijer
,
H.
,
2016
, “
Hybrid Trajectory Tracking for a Hopping Robotic Leg
,”
IFAC-PapersOnLine
,
49
(
14
), pp.
107
112
.10.1016/j.ifacol.2016.07.993
35.
Rijnen
,
M.
,
de Mooij
,
E.
,
Traversaro
,
S.
,
Nori
,
F.
,
van de Wouw
,
N.
,
Saccon
,
A.
, and
Nijmeijer
,
H.
,
2017
, “
Control of Humanoid Robot Motions With Impacts: Numerical Experiments With Reference Spreading Control
,”
Proceedings of IEEE International Conference on Robotics Automation
, Marina Bay Sands, Singapore, May 29–June 3, pp.
4102
4107
.10.1109/ICRA.2017.7989472
36.
Branicky
,
M. S.
,
1998
, “
Multiple Lyapunov Functions and Other Analysis Tools for Switched and Hybrid System
,”
IEEE Trans. Automat. Contr.
,
43
(
4
), pp.
475
482
.10.1109/9.664150
37.
Gu
,
Y.
,
Yao
,
B.
, and
Lee
,
C. S. G.
,
2018
, “
Straight-Line Contouring Control of Fully Actuated 3-D Bipedal Robotic Walking
,”
Proceedings of American Control Conference
, Milwaukee, WI, June 27–29, pp.
2108
2113
.10.23919/ACC.2018.8431067
38.
Gao
,
Y.
, and
Gu
,
Y.
,
2019
, “
Global-Position Tracking Control of a Fully Actuated NAO Bipedal Walking Robot
,”
Proceedings of American Control Conference
, Philadelphia, PA, July 10–12, pp.
4596
4601
.10.23919/ACC.2019.8815144
39.
Westervelt
,
E. R.
,
Grizzle
,
J. W.
,
Chevallereau
,
C.
,
Choi
,
J. H.
, and
Morris
,
B.
,
2007
,
Feedback Control of Dynamic Bipedal Robot Locomotion
, Vol.
28
,
CRC Press
, Boca Raton, FL.
40.
Olson
,
E.
,
2011
, “
AprilTag: A Robust and Flexible Visual Fiducial System
,”
Proceedings of IEEE International Conference on Robotics Automation
, Shanghai, China, May 9–13, pp.
3400
3407
.10.1109/ICRA.2011.5979561
41.
Nguyen
,
Q.
,
Da
,
X.
,
Grizzle
,
J. W.
, and
Sreenath
,
K.
,
2020
, “
Dynamic Walking On Stepping Stones With Gait Library and Control Barrier Functions
,”
Algorithmic Foundations of Robotics XII
, Springer, Berlin/Heidelberg, pp.
384
399
.
42.
Galloway
,
K.
,
Sreenath
,
K.
,
Ames
,
A. D.
, and
Grizzle
,
J. W.
,
2015
, “
Torque Saturation in Bipedal Robotic Walking Through Control Lyapunov Function-Based Quadratic Programs
,”
IEEE Access
,
3
, pp.
323
332
.10.1109/ACCESS.2015.2419630
43.
Liao
,
J.
,
Chen
,
Z.
, and
Yao
,
B.
,
2017
, “
High-Performance Adaptive Robust Control With Balanced Torque Allocation for the Over-Actuated Cutter-Head Driving System in Tunnel Boring Machine
,”
Mechatronics
,
46
, pp.
168
176
.10.1016/j.mechatronics.2017.08.007
44.
Yuan
,
M.
,
Chen
,
Z.
,
Yao
,
B.
, and
Liu
,
X.
,
2019
, “
Fast and Accurate Motion Tracking of a Linear Motor System Under Kinematic and Dynamic Constraints: An Integrated Planning and Control Approach
,”
IEEE Trans. Control Syst. Tech.
, 29(2), pp.
804
811
.10.1109/TCST.2019.2955658