Asymptotically Stable Walking of a FiveLink Underactuated 3D Bipedal Robot
Abstract
This paper presents three feedback controllers that achieve an asymptotically stable, periodic, and fast walking gait for a 3D bipedal robot consisting of a torso, revolute knees, and passive (unactuated) point feet. The walking surface is assumed to be rigid and flat; the contact between the robot and the walking surface is assumed to inhibit yaw rotation. The studied robot has 8 DOF in the single support phase and 6 actuators. In addition to the reduced number of actuators, the interest of studying robots with point feet is that the feedback control solution must explicitly account for the robot’s natural dynamics in order to achieve balance while walking. We use an extension of the method of virtual constraints and hybrid zero dynamics, a very successful method for planar bipeds, in order to simultaneously compute a periodic orbit and an autonomous feedback controller that realizes the orbit, for a 3D (spatial) bipedal walking robot. This method allows the computations for the controller design and the periodic orbit to be carried out on a 2DOF subsystem of the 8DOF robot model. The stability of the walking gait under closedloop control is evaluated with the linearization of the restricted Poincaré map of the hybrid zero dynamics. Most periodic walking gaits for this robot are unstable when the controlled outputs are selected to be the actuated coordinates. Three strategies are explored to produce stable walking. The first strategy consists of imposing a stability condition during the search of a periodic gait by optimization. The second strategy uses an eventbased controller to modify the eigenvalues of the (linearized) Poincaré map. In the third approach, the effect of output selection on the zero dynamics is discussed and a pertinent choice of outputs is proposed, leading to stabilization without the use of a supplemental eventbased controller.
I Introduction
The primary objective of this paper^{1}^{1}1A preliminary version of this paper was presented in [1]. is to contribute to the feedback control of 3D bipedal robots that do not rely on large feet and slow movement for achieving stability of a walking gait. We assume here an unactuated point contact at the leg end and, for a simple 5link robot, seek a timeinvariant feedback controller that creates an exponentially stable, periodic walking motion. Our approach is based on an extension of the method of virtual constraints, which was developed in [2, 3, 4, 5] for planar robots, and is extended here to the case of spatial robots. Virtual constraints are holonomic constraints on the robot’s configuration that are asymptotically achieved through the action of a feedback controller. Their function is to coordinate the evolution of the various links of the robot throughout a stride—which is another way of saying that they reduce the degrees of freedom. By using virtual constraints to achieve link coordination on a bipedal robot, different gaits can be more easily programmed than if the links were coordinated by hardware constraints.
The work most closely related to ours is [6], where the control of a 3D walker was decomposed into the study of its motion in the sagittal plane and the frontal plane; see also [7] for a related decomposition result on control in the frontal plane. The method of virtual constraints was applied in [6] to regulate the sagittal plane motion of the biped, while an inverted pendulum approximation of the dynamics was used to design a controller for the frontal plane. An eventbased controller was then introduced to synchronize the phasing of the independently designed sagittal and frontal plane controllers. The overall closedloop system was shown to be stable through simulation and subsequently through experimentation. In our approach, we do not decompose the model into sagittal and frontal plane motions, and coupling of the sagittal and frontal plane dynamics is introduced into the controller from the very beginning.
A very interesting study of the feedback control of underactuated spatial robots has been given in [8], where a controller for a fivelink 3D robot with unactuated point feet has been designed on the basis of linearizing the robot’s dynamic model along a periodic orbit. So that the controller would be timeinvariant, the orbit was parameterized with a configuration variable that is strictly monotonic throughout a normal gait, as in [2, 3, 4, 5], before linearization was applied. The (withinstride) control law is designed on the basis of a discretetime approximation of the linearized model, which makes stability of the closedloop system difficult to assess.
Other important work includes [9] and references therein, where the analysis of passive spatial bipeds is presented. The emphasis in their work is on energy efficiency and underactuation; the role of feedback control in achieving a wide range of behaviors is not emphasized. On the other hand, the work in [10, 11] seeks energy efficiency and a large basin of attraction under the assumption of full actuation; in particular, full actuation between the leg and ground is assumed (pitch, roll and yaw), as opposed to the unactuated assumption made here. Very careful stability analysis of the closedloop system is provided through geometric (Routhian) reduction. This work is taken one step further in [12], where, starting from a 2D (sagittal plane) passive limit cycle, the authors use geometric reduction to first achieve control of the frontal plane motion and then a second stage of geometric reduction to achieve steering within the walking surface.
To the best of our knowledge, other work on the control of spatial robots either assumes full actuation or does not provide significant analysis of the closedloop system. There are many control strategies based on the zero moment point ZMP [13], with one of the more famous users being the robot ASIMO [14]. In this approach, a desired trajectory of the ZMP is defined and successive inner control loops are closed on the basis of the ZMP. In the work of [15], predictive control is performed on the basis of the position of the center of mass and a simplified model of the robot in order to achieve a desired ZMP trajectory. Recently, online adjustment of the ZMP has been added [16]; this control method is implemented on the robot HRP2. The control of the ZMP ensures that the supporting foot will not rotate about its extremities, but this does not ensure stability in the sense of convergence toward a periodic motion, as proved in [17].
Ii Model
A simplified model of a spatial bipedal robot is given here. The model was chosen to be complex enough to capture interesting features of gait control that do not occur in planar robots, and simple enough that the presentation of the ideas will remain transparent. It is our expectation that the ideas presented here apply to a wider class of bipeds, but proving such a conjecture is not the objective of this paper.
Iia Description of the robot and the walking gait
The 3D bipedal robot discussed in this work is depicted in Fig. 1. It consists of five links: a torso and two legs with revolute one DOF knees that are independently actuated and terminated with “pointfeet”. Each hip consists of a revolute joint with two degrees of freedom and each degree of freedom (DOF) is independently actuated. The width of the hips is nonzero. The stance leg is assumed to act as a passive pivot in the sagittal and frontal planes, with no rotation about the zaxis (i.e., no yaw motion), so the leg end is modeled as a point contact with two DOF and no actuation. As discussed below, this model corresponds to the limiting case of robot with feet when the size of the feet decreases to zero. The unactuated DOF at the leg ends correspond to the classical DOF of an ankle. The DOF corresponding to the swingleg ankle are not modelled. In total, the biped in the single support phase has eight DOF, and there are two degrees of underactuation.
In a more complete model with feet, one would include at a minimum two degrees of freedom in the ankles, corresponding to motion in the sagittal and frontal planes. Moreover, it would be assumed that the friction between the foot and the ground is sufficient to prevent sliding, and hence in particular, rotation of the foot in the yaw direction [18, 19]. In most control studies, which assume flatfooted walking, the key limiting factor is the determination of ankle torques that respect the ZMP condition, namely the ground reaction forces must remain with the convex hull of the foot [20, 21, 13]. The lighter and smaller feet, the tighter are the constraints on the allowable ankle torques, increasing the difficulty of determining feasible walking trajectories and subsequently, stabilizing feedback controllers.
The objective of our study is twofold: to show that stable flatfooted walking is possible with zero ankle torques, thereby removing an important obstacle to previous studies on walking; and to prepare for a future study of walking gaits that allow phases with rotation about the heel and toe [22, 17]. In the limiting case of a point foot, the allowable ankle torque becomes zero, leading to an unactuated contact at the leg end. Given that the envisioned application of our results is to robots with feet, where yaw rotation is naturally inhibited by friction, it makes sense to assume in the pointfoot model that yaw rotation is not allowed.
In summary, the following assumptions are made in the present study:

Each link is rigid and has mass.

Walking consists of two alternating phases of motion: single support and double support.

The double support phase is instantaneous and occurs when the swing leg impacts the ground.

At impact, the swing leg neither slips nor rebounds.

The swing and stance legs exchange their roles at each impact.

The gait is symmetric in steady state.

Walking takes place on a flat surface.
A more detailed list of hypotheses is given in [23, Chap. 3] for planar robots, and with the obvious modifications for spatial robots, those hypotheses apply equally well here.
Since the gait is composed primarily of single support phases, the variables used to describe the robot are adapted to this phase of motion. The robot is represented as a tree structure. The stance foot, which is fixed on the ground, is the base of the tree structure. A set of generalized coordinates is shown in Fig. 1. Absolute angles are roll and pitch angles of the stance leg, respectively. Angles and are the relative joint angles of the stanceleg knee and swingleg knee, respectively. Angles and are the joint angles of the stance leg relative to the torso along the axis and the axis, respectively, and angles and are the joint angles of the swing leg relative to the torso along the axis and the axis, respectively. The coordinates are unactuated (due to the passive contact), while are independently actuated.
The position of the robot with respect to an inertial frame is defined by adding the four variables , where , and are the Cartesian coordinates of the stance foot^{2}^{2}2The leg ends are referred to as feet or point feet., and defines the rotation along the zaxis of the stance leg. These variables are constant during each single support phase.
We have chosen to define the generalized coordinates with respect to the contact point of the current stance foot. When leg2 is the supporting leg, the variables are defined as shown in Fig. 2 and the same notation is employed as when the supporting leg is leg1, viz. Fig. 1. Hence, at each leg exchange (i.e., impact), the variables undergo a jump due to the change in location of the reference frame.
During each single support phase, only one set of coordinates is used, depending on which leg is the supporting leg. In double support, either set of coordinates may be used. The transformation from one set of coordinates to the other is nonlinear [8], but it can be computed in closed form by standard means^{3}^{3}3The transformation from one set of coordinates at the end of a step, for example, to the other set of coordinates is done as follows. Compute the orientation and the angular velocity of the swing leg shin. From this, one deduces , and that are compatible with this orientation, and then one deduces , and yielding the angular velocity of the swing shin. The angles to exchange their roles viz ..
The legs exchange roles from one step to the next. If is the duration of a step, on a periodic walking cycle, due to the choice of coordinates in Figs. 1 and 2, we must have
(1) 
The last condition yields a motion along the axis.
Remark: Even though the model does not include rotation about the supporting foot, nor any other vertical axis, such as at the hip, coupling between the rotations in the sagittal and frontal planes can yield a net rotation about the vertical axis from one step to the next. Hence it is important to keep track of . As an example, Fig. 3 shows the projection onto the plane of a solution of the model yielding a circular motion. This was obtained by modifying the control law of VIII to have unequal step lengths on the right and left legs.
IiB Dynamic model
The dynamic models for single support and impact (i.e., double support) are derived here assuming support on leg . The models for support on leg can be written in a similar way by using a hip width of in place of (this places the swing leg on correct side of the stance leg). The EulerLagrange equations yield the dynamic model for the robot in the single support phase as
(2) 
where is the positivedefinite massinertia matrix, is the vector of Coriolis and gravity terms, is an fullrank, constant matrix indicating whether a joint is actuated or not, and is the vector of input torques. Following standard practice in the literature, the double support phase is assumed to be instantaneous. However, it actually consists of two distinct subphases: the impact, during which a rigid impact takes place between the swing foot and the ground, and coordinate relabeling. During the impact, the biped’s configuration variables do not change, but the generalized velocities undergo a jump. The derivation of the impact model in double support phase requires the use of the vector . Conservation of angular momentum of the robot about the end of the swing leg during the impact process, in combination with the swing leg neither slipping nor rebounding at impact, yields
(3) 
where and are the extended velocities before and after the impact, respectively, is the reaction force at the contact point, is the extend massinertia matrix, and is the Jacobian matrix for the position of the swing foot and its orientation in the plane. Analogously to [2], the overall impact model is written as
(4) 
and
(5) 
and is obtained from solving (3) and projecting down to the generalized coordinates for support on leg .
Define state variables as , and let and , where the subscript denotes the stance leg number. Then a complete walking motion of the robot can be expressed as a nonlinear system with impulse effects, as shown in Fig. 4 and written as
(6) 
where is the switching surface,
and
When leg2 is the support leg, the same derivation produces , , and .
Iii Virtual constraints
The method of virtual constraints, which has proven very successful in designing feedback controllers for stable walking in planar bipeds [2, 3, 4, 5], will be applied to the 3D biped of the previous section. In this method, one holonomic constraint per actuator is proposed in the form of an output that, when zeroed by a feedback controller, enforces the constraint. The most direct form of the constraint is
(7) 
where is the vector of actuated coordinates, is a quantity that is strictly monotonic (i.e., strictly increasing or decreasing) along a typical walking gait, and is the desired evolution of the actuated variables as a function of . Roughly speaking, is used to replace time in parameterizing a periodic motion of the biped. In a forward walking motion, the coordinate of the hip is monotonically increasing. Hence, if the virtual stance leg is defined by the line that connects the stance foot to the stance hip, the the angle of this leg in the sagittal plane is monotonic. When the shin and the thigh have the same length, the angle of the virtual leg in the sagittal plane can be selected as
(8) 
(the minus sign is used to make strictly increasing over a step).
The torque required to remain on the virtual constraint surface corresponding to can be computed as^{4}^{4}4As shown in [24], [23, pp. 60], an expression without inversion of the massinertia matrix is also possible. The expression given in (9) is more compact.
(9) 
This leads to an inputoutput linearizing controller to asymptotically drive the state of the robot to the constraint surface, assuming it does not initially start there [25] [23, Chap. 5],
(10) 
which results in
(11) 
In other words, determining the constraints is equivalent to the design of a feedback controller in the single support phase, up to the choice of the gains , , and such that (11) is exponentially stable and converges sufficiently rapidly with respect to the duration of a single support phase; see [23, Chap. 4].
The next objective is to determine the behavior of the robot under the virtual constraints. This task is simplified by noting that enforcing the virtual constraints, , results in and reduces the dimension of the dynamics.
Let denote the unactuated joints and denote the controlled joints, which are selected here to be the actuated joints. A linear relation exist between , and ,
(12) 
where is an invertible matrix. Then (2) can be rewritten as
(13) 
The first two lines of the RHS of this equation are zero, yielding
(14) 
where is the upper left submatrix of , is the upper right submatrix of and consists of the first two lines of . Substituting the expressions of , and corresponding to the virtual constraints, the dynamic model of the single support phase is now reduced to a lowdimensional, 2DOF, autonomous system,
(15) 
which is called the swing phase zero dynamics [26], [23, Chap. 5].
One can clearly see that the dynamic properties of the swing phase zero dynamics depend on the particular choice of the virtual constraint . How to determine a choice for that results in a periodic walking motion is summarized in the next section.
Iv Design for a Symmetric Periodic Gait
The objective of this section is to design virtual constraints that correspond to a periodic motion of the robot. The gait considered is composed of single support phases separated by impacts as described in Fig. 4. The legs exchange roles from one step to next, and due to symmetry, the study of a gait can be limited to a single step and the use of the symmetry relation (1).
Iva Virtual constraints and Bezier polynomials
The problem of designing the virtual constraints will be transformed into a parameter optimization problem as in [23, Chap. 6]. Here, our main goal is to obtain a periodic motion; optimality is not so crucial. To simplify the optimization process, the number of variables used in the optimization problem is first reduced. This is accomplished by exploiting boundary conditions that arise from periodicity. Bezier polynomials are parametric functions that allow one to easily take into account boundary conditions on the configuration and velocity at the beginning and end of a step.
The initial and final configuration and velocity of the robot for a single support phase are important for defining the passage between the single and double support phases. Because the terminal configuration of the robot is chosen to be the instant before the double support configuration, both legs are in contact with the ground and therefore only seven independent variables are needed to describe this configuration (a closed kinematic chain). These variables parameterize the final configuration of the first step denoted . The eight joint velocities are independent and are also added.
Knowing the final state of the single support phase, the impact model (4) and (5) determines the initial state of the ensuing single support phase. The symmetry condition (1) then gives the initial state of the first step: , . The initial orientation of the robot is calculated such that the orientation for the second step is symmetric to the orientation for the first step in order that no yaw rotation is observed during the nominal (periodic) gait.
To obtain a periodic gait, the single support must be such that the state of the robot evolves from , to , . For given desired initial and final state values, virtual constraints can be easily deduced to connect the desired values of the actuated variables. However, the evolution of the unactuated variables is known only by integration of the dynamics (15); a desirable dynamic behavior is imposed on these variables by the use of equality and inequality constraints in the optimization process.
IvB Specifics
Here, Bezier polynomials of degree are chosen to define the virtual constraints^{5}^{5}5A degree greater than can also be chosen, in which case the number of optimization variables increases [27].. The virtual constraints are expressed as functions of the variable ; see (8). From and , the initial and final values of , denoted and , can be calculated. Let
(16) 
where
(17) 
is the normalized independent variable. The coefficients of the Bezier polynomials, , are vectors of real numbers. They must be determined so as to join to and to , (the additional subscript “” denotes the actuated variables) when varies from to , yielding
(18) 
The evolution of the unactuated variables is calculated by integration of the dynamic subsystem (15), that is, the stance phase zero dynamics, starting from the initial state and terminating at , where denotes the initial value of for the first step.
When the evolution of the unactuated variables is calculated, because the evolution of the actuated variable is given by (7) and (16), the torque required to zero the constraints, i.e, in (9), can be calculated by lines 3 through 8 of (13), and the ground reaction force expressed in the inertial reference frame (see Fig. 1) can be calculated as well.
The search for a periodic walking motion can now be cast as a constrained nonlinear optimization problem: Find the 15 optimization parameters prescribing that minimize the integralsquared torque per step length^{6}^{6}6Torque being proportional to current in a DC motor, integralsquared torque is a rough approximation of energy dissipated in the motors.,
(19) 
where is the walking period
and is the step length, while satisfying symmetry (1), and subject to the following:
inequality constraints

is strictly increasing (i.e, along the solution);

the swing foot is positioned above the ground (;

a notakeoff constraint, ;

a friction constraint, ;
equality constraints
and a set of conditions imposing periodicity,
where and result from the integration of the zero dynamics and the walking period is such that .
The above procedure can be performed in MATLAB with the FMINCON function of the optimization toolbox. A fixedpoint solution minimizing defines a desired periodic walking cycle (or nominal orbit). The criterion being optimized (19) has many local minima and the optimization technique used is local. Thus, the obtained optimal periodic motion depends on the initial set of optimization parameters.
IvC An example periodic motion minimizing integralsquared torque
The physical parameters of the 3D biped studied here are given in Table I. For these parameters, a periodic orbit was computed following the technique presented in the previous subsection. We obtained a periodic motion defined by ), where
g  W  L1  L2  L3  m1  m2  m3 
9.81  0.15  0.275  0.275  0.05  0.875  0.875  5.5 
A stickfigure diagram for the first step of the periodic walking gait is presented in Fig. 5. The walking gait has a period of seconds, a step size of m, and an average walking speed of m/sec, or body lengths per second. The step width is m, close to the hip width. The nominal gait’s joint profiles and angular velocities over two consecutive steps are shown in Figs. 6 and 7, respectively. The unactuated and actuated variables are presented; note that is monotonic over each step. Fig. 8 shows the torque required to produce the periodic motion, which is less than Nm for each joint. Fig. 9 shows the profile of the ground reaction force on the stance foot and the profile of the swing leg tip; this figure shows that the inequality constraints are satisfied on the nominal motion.
V Evaluating the Stability of a Walking Cycle
The stability of a fixedpoint can be tested numerically by linearizing the Poincaré map about the fixedpoint as presented in [1]. This numerical stability test has a high computational cost, however, because it requires the estimation of the Jacobian of the Poincaré map, in a space of dimension 1, where is the number of independent joint coordinates; here . We propose a slight modification of the control law in order to be able to study the stability of the closedloop system in a reduceddimensional state space.
Va Hybrid zero dynamics (HZD) and a stability test in a reduced space
The control law (10) is such that, on the periodic orbit, the virtual constraints (7) are identically satisfied. However, off the periodic orbit, even if the virtual constraints are satisfied at the end of given step, they will not in general be satisfied^{7}^{7}7This may be true for several reasons, one of which is that the virtual constraints may not have been chosen to be compatible with the impact map. at the beginning of the next step. Consequently, the behavior of the robot cannot be deduced from the behavior of the uncontrolled variables and the simulation of the complete model is required to predict the behavior of the robot. In the language of [27], [23, Chap. 5], while the feedback control law (10) has created a zero dynamics of the stance phase dynamics, it has not created a hybrid zero dynamics, that is, a zero dynamics of the full hybrid model (6).
If the control law could be modified so as to create a hybrid zero dynamics, then the study of the swing phase zero dynamics (15) and the impact model would be sufficient to determine the stability of the complete closedloop behavior of the robot, thereby leading to a reduceddimension stability test. A modification of the control law to achieve a hybrid zero dynamics was first proposed in [28]; a second more easily implementable method has been given in [29], along with a complete stability analysis.
Following [29], the virtual constraints are modified stride to stride so that they are compatible with the initial state of the robot at the beginning of each step. The new output for the feedback control design is
(20) 
This output consists of the previous output (7), and a correction term that depends on (7) evaluated at the beginning of the step, specifically, and , where the subscript “” denotes the initial value for the current step. The values of are updated at the beginning of each step and held constant throughout the step. The function is taken to be a threetimes continuously differentiable function of such that^{8}^{8}8In our specific application, we used a fifth order polynomial for ; continuity of position, velocity and acceleration is ensured at .
(21) 
With designed in this way, the initial errors of the output and its derivative are smoothly joined to the original virtual constraint at the middle of the step. In particular, for any initial error, the initial virtual constraint is exactly satisfied by the end of the step.
The robot’s behavior with the new control law is very close to its behavior with the fixed virtual constraint, the difference being that the initial output error is zeroed by the choice of the function instead of being approximately zeroed by equation (11). Our choice of zeroes smoothly the initial error, and the initial virtual constraints are joined at the middle of the step. Consequently, the torque required to zero the error is less than with a highgain controller, the variation of ground reaction forces is reduced, and sliding or takeoff is avoided.
Under the new control law defined by (20), the behavior of the robot is completely defined by the impact map and the swing phase zero dynamics (15), where is replaced by . The stability of a fixedpoint can now be tested numerically using a restricted Poincaré map defined from to , where and is the switching surface. The key point is that in , the state of the robot can be represented using only three independent variables, .
The restricted Poincaré map induces a discretetime system . From [29], for sufficiently small in (10), the linearization of about a fixedpoint determines exponential stability of the fullorder closedloop robot model. Define . The Poincaré map linearized about a fixedpoint gives rise to a linearized system,
(22) 
where the () square matrix is the Jacobian of the Poincaré map and is computed as follows
where
(23) 
and The quantities , , are small perturbations introduced to calculate the linearized model, and the denominator in (23) must of course be interpreted as the scalar perturbation used in computing . In general, the calculation of the Jacobian is sensitive to the amplitude of the perturbation . The calculation of matrix requires six evaluations of the function . Each evaluation of this function is composed of the reconstruction of the vector from with , the calculation of the impact map (4)(5), the calculation of , the calculation of and the integration of the swing phase zero dynamics. A fixedpoint of the restricted Poincaré map is locally exponentially stable, if, and only if, the eigenvalues of have magnitude strictly less than one [23, Chap. 4].
VB Example of the periodic motion minimizing integralsquared torque
We consider the virtual constraints corresponding to the optimal periodic motion obtained in Section IVC, and the control law defined by (20) is used.
To study the stability of this control law around the periodic motion, the eigenvalues of the linearized restricted Poincaré map are computed (for , ), yielding
(24) 
The three eigenvalues are:
One eigenvalue has magnitude greater than one and hence the gait is unstable under this controller.
We have found that for most periodic motions optimized with respect to integralsquared torque per step length, (19), the obtained gait is unstable under the control law defined by (20). In the following sections, three strategies to obtain stable walking will be considered. In the first strategy, closedloop stability is directly considered in the design of the periodic motion. In the second strategy, a stridetostride controller is introduced to stabilize the gait. In the third strategy, freedom in the selection of the controlled outputs^{9}^{9}9The controlled outputs are no longer the actuated variables as in (10), but a judiciously chosen linear combination of . A convenient choice of outputs is given. is used to obtain a stable walking cycle using only withinstride control.
Vi Periodic motion optimized with respect to a stability criterion
Because we are using optimization to compute a periodic solution of the model, it is possible to consider the stability condition (i.e., magnitude of eigenvalues less than one) either as a criterion in the optimization process, or as a constraint.
Using the stability test above, which has a reasonable calculation cost, the selection of a periodic motion that minimizes the maximum of the magnitudes of the eigenvalue of was performed,
where is the set of eigenvalues of . The optimization process is similar to the optimization process described in Section IV, only the criterion is changed. Starting from the same initialization for the optimization process, a new periodic trajectory is obtained. This trajectory is defined by
The stickdiagram for a step is presented in Fig. 10. The walking cycle has a period of seconds, a step length of m, and an average walking speed of m/sec. The step width is m, close to the hip width. The nominal gait’s joint profiles and angular velocities over two steps are shown in Figs. 11 and 12, respectively. The unactuated and actuated variables are presented, and it is shown that is monotonic. The duration of a step is shorter than in the case of torque optimization (i.e., the motion is faster). The movement in the frontal plane has smaller amplitude; indeed, the maximal sidetoside variation of the center of mass in the frontal direction is less than cm (see Fig. 10), while this deviation was more than cm in the case of torque optimization (see Fig. 5).
The eigenvalues of the matrix are
All of the eigenvalues have magnitude less than 1.0, indicating that the obtained nominal orbit is locally exponentially stable. To illustrate the orbit’s local stability, the 3D biped’s full model in closedloop is simulated with an initial state perturbed from the fixedpoint . An initial error of is introduced on each joint.
Fig. 13 shows the evolution of the final values of the uncontrolled variables from one step to the next. These variables converge toward the periodic motion. With the modification of the virtual constraints by the introduction of the polynomial , the output and its derivative are zero throughout each step.
Fig. 14 shows the phaseplane evolution of the first four variables as an illustration of the behavior of the state of the robot. The discontinuity at impact manifests itself as a straight line on the plots. The desired motion of the robot is periodic, with exchange of the legs at impact. This behavior corresponds to a motion with a period of one step for the variables , and and a period of two steps for ; see (1). The convergence towards a periodic motion is clear for both the controlled and uncontrolled variables.
Vii StridetoStride controller
If a desired periodic gait is not stable, or if the corresponding rate of convergence is not sufficiently rapid, then eventbased control may be designed and integrated with the continuous, stance phase controller [30].
Viia General method
Let be a vector of parameters that are held constant during the stance phase and updated at each impact. The parameters could be a subset of the parameters used to specify the virtual constraints, , or an auxiliary set of parameters. Moreover, the parameter updates can be done on the basis of the fullstate of the robot as in [1], or on the basis of the state of the hybrid zero dynamics, which is done here.
The output in (20) is augmented with an additional term depending on a vector of parameters :
(25) 
with^{10}^{10}10In our specific application, we used a sixth order polynomial for . Continuity of position, velocity and acceleration is ensured at .
(26) 
This choice is convenient because it allows the value of the virtual constraints at midstance to be updated, without requiring a redesign of the eventbased controller that created the hybrid zero dynamics. As long as the impact occurs for , which happens when the motion is close to the periodic orbit, the final state of the robot is on the original hybrid zero dynamics independently of any modification for the previous step.
The Poincaré map can now be viewed as a nonlinear control system on with inputs , where is the value of during the step , namely . Linearizing this nonlinear system about the fixed point and the nominal parameter value leads to
(27) 
where and is the Jacobian of with respect to . Designing a feedback matrix
(28) 
such that the eigenvalues of have magnitude strictly less than one will exponentially stabilize the fixed point .
ViiB Example of the periodic motion minimizing the integralsquared torque
The periodic motion described in Section IVC can be stabilized using a stridetostride controller. The virtual constraints are modified according to equation (25); for the th step the virtual constraints are
(29) 
The matrix was computed numerically, analogously to (22), yielding
(30) 
The matrix is given in (24).
The gain matrix was calculated via DLQR so that the statefeedback law minimized the cost function , subject to the state dynamics (27). The coefficient allows a tradeoff to be made between the convergence rate and the amplitude of the change of the virtual constraint, which affects the region of convergence. For we obtain
(31) 
The eigenvalues of the linearized Poincaré map in closed loop indicate that the gait is stable with the addition of the stridetostride controller^{11}^{11}11The controller of Sec. VB is used within stride and the eventbased controller is applied at each impact.. The eigenvalues become
To illustrate the orbit’s local stability, the 3D biped’s complete model in closedloop is simulated with the initial state perturbed away from the fixedpoint . An initial error of is introduced on each joint and a velocity error of is introduced on each joint velocity.
Fig. 15 shows the evolution of the values of the uncontrolled variables at the end of each step, when eventbased DLQR controller is used. These variables clearly converge toward the periodic motion.
Fig. 16 shows phaseplane plots of the first four variables. The convergence towards a periodic motion is clear for the controlled and uncontrolled variables.
Viii Improved Selection of the Outputs to be Controlled
In the previous two sections, the controlled variables driven by the virtual constraints are simply the actuated variables, ; see (7). The choice of the controlled variables directly affects the zero dynamics in (15). It is shown here that for the same desired periodic motion, the stability of the closedloop system can be dramatically improved through a judicious choice of the controlled variables.
Viiia Effect on the swing phase zero dynamics
For simplicity, we limit our analysis to the case of controlled variables that are linear with respect to the configuration variables. Thus the controlled variables are
(32) 
where is a constant matrix with invertible. A known periodic motion can be reparameterized^{12}^{12}12This assumes that is monotonic. as function of the variable , yielding . The virtual constraint for the new controlled variables then yields the output
(33) 
When the constraint is satisfied, , equation (33) allows us to solve for , giving
(34) 
Substituting this equation into (14), we obtain for the swing phase zero dynamics
(35) 
The nominal periodic motion satisfies both equations (35) and (15), but the two equations produce different solutions away from the periodic motion. When the principle of virtual constraints is applied to a system with only one degree of underactuation, namely , which is common for example in planar bipeds, the swing phase zero dynamic is not affected by the choice of the output, and therefore the stability of a periodic orbit (i.e., walking motion) is not modified; only the transient motion can be different. In other words, for planar robots with one degree of underactuation, the stability depends only on the trajectory of the periodic orbit and not on the choice of virtual constraints used to achieve it [31], [23, pp. 160].
In the case of a system with two degrees of underactuation, the choice of the controlled output can affect the stability of the gait via the choice of . In order to illustrate this property, a new choice of output is proposed. This choice is based on the following physical reasoning: The motion in the frontal direction is difficult to stabilize. The position of the center of mass in the frontal direction is important. If, at touchdown, the center of mass is not between the feet, but outside the position of the next supporting foot, the robot will topple sideways. Thus, the control of the variable (which regulates step width on the swing leg) is replaced by the control of the distance between the swing leg end and the center of mass along the frontal direction. To obtain a linear output, this function is linearized around the touchdown configuration to define in (33).
ViiiB Example of the periodic motion minimizing integralsquared torque
The periodic motion described in Section IVC can be stabilized using the new controlled output. As mentioned in the previous subsection, the actuated joints , , , , and are controlled via virtual constraints just as in the original control law. A new output is considered, with this output no longer based on but instead on distance between the swing leg end and the center of mass along the frontal direction.
For this trajectory, for support on leg 1, the linearization around of the distance between the swing leg end and the center of mass along the frontal direction yields
On the periodic orbit, this distance is evaluated and approximated by a function of , denoted . The new controlled output is then
(36) 
When the control law is defined using this new output, the walking gait is stable, as can be shown via the calculation of the eigenvalues of the linearization of the restricted Poincaré map, . The eigenvalues are:
To illustrate the orbit’s local exponential stability, the 3D biped’s model in closedloop is simulated with the initial state perturbed from the fixedpoint . An initial error of is introduced on each joint and a velocity error of is introduced on each joint velocity.
Fig. 17 shows the evolution of values of the uncontrolled variables at the end of each step when the new output is used. These variables clearly converge toward the periodic motion.
Fig. 18 shows phaseplane plots of the first four variables. The convergence towards a periodic motion is clear for the controlled and uncontrolled variables.
Ix Conclusions
A simple 3D bipedal model has been studied, with the objective of developing a timeinvariant feedback control law that induces asymptotically stable walking, without relying on the use of large feet. For this reason, a biped consisting of five links, connected to form two legs with knees and a torso, was assumed to have point feet with no actuation between the feet and ground. Inspired by its success in solving similar problems for planar robots, the method of virtual constraints was applied to the 3D robot, with the virtual constraints chosen via optimization as suggested in [5].
The main contributions of the paper are:

The development of an efficient optimization process to obtain periodic motions in 3D for a robot with point contact feet.

The computation of humanlike periodic walking motions that can be stable or unstable, depending on the choice of actuated variables and corresponding virtual constraints.

The numerical study of stability on the basis of a lowdimensional subsystem corresponding to the hybrid zero dynamics. The Poincaré return map was computed in a space of dimension three for a robot with two degrees of underactuation.

The use of a stridetostride controller to stabilize a walking motion that was not naturally stable for a given choice of controlled outputs.

The discovery of the importance of the selection of the controlled outputs on the stability of a given periodic motion.
Points (1), (3) and (4) can be viewed as direct extensions of the work previously done on planar walking robots. However, the points (2) and (5) are specific to the study of robots in 3D.
The control strategy developed here for the control of robots with point feet can be extended to the case of robots with actuated, nontrivial feet, as was done in [17], [32], [22] for planar robots. In this case, foot rotation can be included as part of a natural gait and the explicit control of the ZMP position can be addressed. Another possible extension would be to continue with the unactuated pointfoot hypothesis, this time including yaw rotation about the foot. It is unclear at this time whether this extension is straightforward or not.
Acknowledgment: The work of C. Chevallereau is supported by ANR grants for the PHEMA project. The work of J.W. Grizzle is supported by NSF grant ECS0600869. The work of C.L. Shih is supported by Taiwan NSC grant NSC962221E011126.
References
 [1] C. Shih, J. Grizzle, and C. Chevallereau, “Asymptotically stable walking of a simple underactuated 3D bipedal robot,” in The 33rd Annual Conference of the IEEE Industrial Electronics Society (IECON), Taipei, Taiwan, Novembre 2007, pp. 2766–2771.
 [2] J. Grizzle, G. Abba, and F. Plestan, “Asymptotically stable walking for biped robots: Analysis via systems with impulse effects,” IEEE Transactions on Automatic Control, vol. 46, pp. 51–64, January 2001.
 [3] F. Plestan, J. Grizzle, E. Westervelt, and G. Abba, “Stable walking of a 7DOF biped robot,” IEEE Transactions on Robotics and Automation, vol. 19, no. 4, pp. 653–668, August 2003.
 [4] C. Chevallereau, G. Abba, Y. Aoustin, E. Plestan, F. Westervelt, C. Canduasde Wit, and J. Grizzle, “RABBIT: A testbed for advanced control theory,” IEEE Control Systems Magazine, vol. 23, no. 5, pp. 57–79, October 2003.
 [5] E. Westervelt, J. Grizzle, and D. Koditschek, “Hybrid zero dynamics of planar biped walkers,” IEEE Transactions on Automatic Control, vol. 48, no. 1, pp. 42–56, January 2003.
 [6] T. Fukuda, M. Doi, Y. Hasegawa, and H. Kajima, Fast motions in Biomechanics and Robotics, ser. Lecture Notes in Control and Information Sciences. Heidelberg, Allemagne: Springer, 2006, ch. MultiLocomotion Control of Biped Locomotion and Brachiation Robot, pp. 121–145.
 [7] A. Kuo, “Stabilization of lateral motion in passive dynamic walking,” International Journal of Robotics Research, vol. 18, no. 9, pp. 917–930, 1999.
 [8] G. Song and M. Zefran, “Underactuated dynamic threedimensional bipedal walking,” in Proceedings of the IEEE International Conference On Robotics and Automation. Orlando, Florida: IEEE Press, May 2006, pp. 854–859.
 [9] S. Collins, A. Ruina, R. Tedrake, and M. Wisse, “Efficient bipedal robots based on passivedynamic walkers,” Science, no. 307, pp. 1082–1085, 2005.
 [10] M. Spong and F. Bullo, “Controlled symmetries and passive walking,,” IEEE Transactions on Automatic Control, vol. 50, no. 7, pp. 1025–1031, July 2005.
 [11] A. Ames and R. Gregg, “Stably extending twodimensional bipedal walking to three,” in Proc. of the 2007 American Control Conferenc, New York, NY, 2007, pp. 177–182.
 [12] R. Gregg and M. Spong, “Reductionbased control with application to threedimensional bipedal walking robots,” in Preprint, October 2007.
 [13] M. Vukobratović, B. Borovac, and V. Potkonjak, “ZMP: A review of some basic misunderstandings,” International Journal of Humanoid Robotics, vol. 3, no. 2, pp. 153–175, June 2006.
 [14] A. Takanishi, Y. Egusa, M. Tochizawa, T. Takeya, and I. Kato, “Realisation of dynamic walking stabilized with trunk motion,” in Proceedings of ROMANSY 7, 1988, pp. 68–79.
 [15] S. Kajita, M. Morisawa, K. Harada, K. Kaneko, F. Kanehiro, K. Fujiwara, and H. Hirukawa, “Biped walking pattern generation by using preview control of zeromoment point,” in Proceedings of the ICRA ’03 IEEE International Conference on Robotics and Automation, vol. 2, 2003, pp. 1620–1626.
 [16] S. Kajita, F. Kanehiro, K. Kaneko, K. Fujiwara, K. Harada, K. Yokoi, and H. Hirukawa, “Biped walking pattern generator allowing auxiliary zmp control,” in Proceedings of the 2006 IEEE/RSJ International Conference on Intelligent Robots and Systems, 2006, pp. 2993–2999.
 [17] J. H. Choi and J. W. Grizzle, “Planar bipedal walking with foot rotation,” in Proc. of the American Control Conference, 2005, pp. 4909–4916.
 [18] Y. Fujimoto and A. Kawamura, “Simulation of an autonomous biped walking robot including environmental force interaction,” IEEE Robotics and Automation Magazine, pp. 33–42, June 1998.
 [19] J. Furusho and A. Sano, “Sensorbased control of a ninelink biped,” International Journal of Robotics Research, vol. 9, no. 2, pp. 83–98, 1990.
 [20] A. Goswami, “Postural stability of biped robots and the footrotation indicator (FRI) point,” International Journal of Robotics Research, vol. 18, no. 6, pp. 523–533, June 1999.
 [21] M. Vukobratovic, B. Borovac, D. Surla, and D. Stokic, Biped Locomotion. Berlin: SpringerVerlag, 1990.
 [22] C. Chevallereau, D. Djoudi, and J. Grizzle, “Stable bipedal walking with foot rotation through direct regulation of the zero moment point,” IEEE Transactions on Robotics, vol. 25, no. 2, pp. 390–401, April 2008.
 [23] E. Westervelt, J. Grizzle, C. Chevallereau, J. Choi, and B. Morris, Feedback Control of Dynamic Bipedal Robot Locomotion, ser. Control and Automation. Boca Raton: CRC Press, June 2007.
 [24] M. Spong, “Energy based control of a class of underactuated mechanical systems,” in Proc. of IFAC World Congress, San Francisco, CA, 1996, pp. 431–435.
 [25] B. Morris and J. Grizzle, “A restricted Poincaré map for determining exponentially stable periodic orbits in systems with impulse effects: Application to bipedal robots,” in IEEE Conf. on Decision and Control. Seville, Spain: IEEE Press, December 2005.
 [26] A. Isidori, Nonlinear Control Systems, 3rd ed. Berlin: SpringerVerlag, 1995.
 [27] E. Westervelt, J. Grizzle, and D. Koditschek, “Zero dynamics of planar biped walkers with one degree of under actuation,” in IFAC 2002, Bareclona, Spain, July 2002.
 [28] B. Morris, E. Westervelt, C. Chevallereau, G. Buche, and J. Grizzle, Fast Motions Symposium on Biomechanics and Robotics, ser. Lecture Notes in Control and Information Sciences. Heidelberg, Germany: SpringerVerlag, 2006, ch. Achieving Bipedal Running with RABBIT: Six Steps toward Infinity, pp. 277–297.
 [29] B. Morris and J. Grizzle, “Hybrid invariant manifolds in systems with impulse effects with application to periodic locomotion in bipedal robots,” IEEE Transactions on Automatic Control, vol. 54, no. 8, pp. 1751–1764, 2009.
 [30] J. W. Grizzle, “Remarks on eventbased stabilization of periodic orbits in systems with impulse effects,” in Second International Symposium on Communications, Control and Signal Processing, 2006.
 [31] C. Chevallereau, A. Formal’sky, and D. Djoudi, “Tracking of a joint path for the walking of an underactuated biped,” Robotica, vol. 22, pp. 15–28, 2004.
 [32] D. Djoudi, C. Chevallereau, and J. Grizzle, “A pathfollowing approach to stable bipedal walking and zero moment point regulation,” in IEEE International Conference On Robotics And Automation, 2007.