Self-synchronization and Self-stabilization of 3D Bipedal Walking Gaits
Christine Chevallereau, Hamed Razavi, Damien Six, Yannick Aoustin, Jessy Grizzle
11 Self-synchronization and Self-stabilization of 3DBipedal Walking Gaits
Christine Chevallereau , Hamed Razavi , Damien Six , Yannick Aoustin , and Jessy Grizzle Abstract —This paper seeks insight into stabilization mecha-nisms for periodic walking gaits in 3D bipedal robots. Based onthis insight, a control strategy based on virtual constraints, whichimposes coordination between joints rather than a temporalevolution, will be proposed for achieving asymptotic convergencetoward a periodic motion. For planar bipeds with one degree ofunderactuation, it is known that a vertical displacement of thecenter of mass—with downward velocity at the step transition—induces stability of a walking gait. This paper concerns the qualitative extension of this type of property to 3D walking withtwo degrees of underactuation. It is shown that a condition onthe position of the center of mass in the horizontal plane atthe transition between steps induces synchronization between themotions in the sagittal and frontal planes. A combination of theconditions for self-synchronization and vertical oscillations leadsto stable gaits. The algorithm for self-stabilization of 3D walkinggaits is first developed for a simplified model of a walking robot(an inverted pendulum with variable length legs), and then itis extended to a complex model of the humanoid robot Romeousing the notion of Hybrid Zero Dynamics. Simulations of themodel of the robot illustrate the efficacy of the method and itsrobustness.
Index Terms —Robotics, Feedback Control, Self-stability,Legged Robots, Mechanical Systems, Hybrid Systems, PeriodicSolutions
I. I
NTRODUCTION
Despite a growing list of bipedal robots that are able towalk in a laboratory environment or even outdoors, stabilitymechanisms for 3D bipedal locomotion remain poorly under-stood. It would be very satisfying to be able to point at arobot and say, “it can execute an asymptotically stable walkinggait, and the stability is achieved through such and suchtheorem, principle, method, etc.” For planar (aka 2D) robotswith one degree of underactuation, virtual constraints andhybrid zero dynamics provide an integrated gait and controllerdesign method that comes with a formal closed-form stabilityguarantee [32][pp. 128-135]. Moreover, the stability conditionis physically meaningful: the velocity of the center of mass(CoM) at the end of the single support phase must be directeddownward [6].For fully actuated 3D bipeds and for 3D bipeds with aspecial form of underactuation, the 2D results extend nicely Christine Chevallereau, Damien Six and Yannick Aoustin are withLaboratoire des Sciences du Num´erique de Nantes (LS2N), CNRS,Ecole Centrale de Nantes, Universit´e de Nantes, Nantes, France, [email protected] Hamed Razavi is withBiorobotics Laboratory (Biorob) of the ´Ecole Polytechnique F´ed´erale deLausanne (EPFL), Lausanne, Switzerland, [email protected] , Jessy W. Grizzle is with Electrical Engineering and Computer ScienceDepartment of the University of Michigan, Ann Arbor, MI, USA, [email protected] [33], [30]. But for 3D robots with more than one degree ofunderactuation, not much is known.Basing the control law design of a fully actuated robot ona model with either a passive ankle or a point foot contactis an interesting intermediate view of the robot. Even in caseof a fully actuated robot, the ankle torques in the frontal andsagittal planes are limited by the size of the foot. A model withtwo degrees of underactuation is therefore useful in order toavoid the use of these torques for a nominal (unperturbed)gait. This allows the limited ankle torque to be saved foradapting the foot orientation in case of uneven terrain or forincreasing the robustness of the control strategy in the face ofperturbations.The possibility and interest of extending the method ofvirtual constraints to robot models with two degrees of un-deractuation have been shown in [7], [26], [17], [1], [12].While the implementation of the control strategy is quitestraightforward, the choice of the virtual constraints is notobvious. Their selection is often based on an optimizationprocess and choice of appropriate controlled outputs and/orthe introduction of a high-level event-based control strategy,neither of which provides insight into the stability mechanism.The objective of this paper is to provide some qualitativeresults on gait characteristics that when used to define virtualconstraints lead to asymptotically stable walking.In contrast to fully actuated bipedal walking, passive walk-ing [20], [21] can be seen as an emergent behavior, wherealternating leg impacts of a biped and the pull of gravitycombine to produce asymptotically stable motions on milddownward slopes. The obtained gait is very efficient froman energy consumption point of view and can be extendedto walking on flat ground [29], [8], [14]. One shortcomingof such an approach is limited robustness with respect toperturbations.As there is no actuation in passive robots, the obtainedwalking gaits may be called self-stable. The self-stabilityproperty is well understood in 2D walking gaits [5], [24],however, little is known on why a passive robot can demon-strate stable walking in 3D. Only a few studies have beendevoted to the investigation of analytical properties that leadto asymptotically stable gaits in 3D underactuated robots [16],[26]. In [7] for example, it has been shown that when acontroller design method that yielded asymptotically stablegaits in planar underactuated walking [5] is applied to 3Dwalking, unstable solutions are common.The notion of self-synchronization , introduced in [25] andgeneralized in [26], sheds some light on the stability mecha-nisms in 3D legged locomotion. a r X i v : . [ c s . R O ] J u l Fig. 1. The humanoid robot Romeo developed by Aldebaran robotics.
The aim of the paper is to provide insight towards under-standing the mechanisms of asymptotic stability in periodicwalking of 3D underactuated robots. Using a control lawbased on virtual constraints for 3D bipedal robots with twodegrees of underactuation, our objective is to propose virtualconstraints leading to self stabilization. In other words, withoutthe use of event-based control, the obtained gaits are asymp-totically stable.As opposed to the formal theorems proven in [32] for 2Drobots, and for the self-synchronization of the 3D LIP [25],[26] the results here for 3D will be numerical in nature. Bynormalizing the studied model, qualitative gait features leadingto asymptotically stable periodic orbits will be uncovered,nevertheless. The method will be first demonstrated on aninverted pendulum model and then will be extended to a fullmodel of the humanoid Romeo [27, Chapter 7], shown inFigure 1.The remainder of the paper is organized as follows. InSection II, as background, the control method based on virtualconstraints for a 3D humanoid robot with 2 degrees of under-actuation is recalled. The stability of the full-order model ofthe robot is discussed in relation to an induced reduced-ordermodel, called the Hybrid Zero Dynamics (HZD). In SectionIII, the difficulty of choosing appropriate virtual constraints ishighlighted. Moreover, an important difference between planarand 3D walking, which is due to an increase of the degreesof underactuation, is illustrated. In particular, it is shown thatunlike the planar case, where at the impact the set describingthe final configuration of the robot in the reduced modelconsists of a single point, in 3D walking with two degreesof underactuation, this set is a manifold of dimension one.In Section IV, where we study a very simplified model ofa humanoid robot, the role of this subset in achieving self-synchronization for the Linear Inverted Pendulum model (LIP)will be demonstrated. Moreover, a new geometric interpre- tation of the self-synchronization property will be providedin this section. Subsequently, in Section IV-D, we introducevertical oscillations of the CoM for the pendulum model.
Self-stabilization properties will be numerically studied in SectionIV-E. In Section V, the results obtained on the simplifiedmodel will be extended to a realistic model of a humanoidrobot. A numerical study shows how the stability propertiesobtained for the pendulum extend to a complete model of abipedal robot. Simulation results illustrate the efficiency of theapproach.II. B
ACKGROUND : A R
EDUCED - ORDER M ODEL A SSOCIATED WITH U NDERACTUATION
This section reviews how to create a reduced-order modelwithin the full-dimensional hybrid model of a bipedal robotwith the following properties: • the dimension of the reduced-order model is determinedby the number of degrees of underactuation of the robot; • the periodic orbits of the reduced-order model are peri-odic orbits of the robot; • locally exponentially stable periodic orbits in the reduced-order model can always be rendered locally exponentiallystable in the full-order model of the robot through the useof feedback control; and • no approximations are being made.Based on the above properties, it can be seen that the reduced-order model captures the effects of underactuation on thedesign of stabilizing feedback control laws for bipedal robots.Later sections of the paper develop properties of this reduced-order model that lead to existence of locally exponentiallystable walking gaits. A. Full-order hybrid model
Let q ∈ Q denote the generalized coordinates for ahumanoid robot in single support (one foot on the ground)and assume the stance ankle has two DoFs (pitch and roll),both assumed to be passive, that is, unactuated. Typically, q consists of body coordinates, such as joint angles of the stanceleg, swing leg, upper-body, and in our case, two additionalworld-frame coordinates capturing the pitch and roll degreesof freedom in the ankle. The robot is assumed to have right/leftsymmetry and the gait studied corresponds to walking along astraight line in the sagittal plane. Moreover, the gait is assumedto be composed of single support phases separated by impacts,which are modeled as instantaneous changes of support. Ateach transition, a relabeling of the joint variables is introduced.This relabeling allows us to work with a single model of therobot.The Lagrangian is assumed to have the form L ( q, ˙ q ) = 12 ˙ q T D ( q ) ˙ q − V ( q ) , and the single support model is given by the standard Lagrangeequations dd t ∂L ( q, ˙ q ) ∂ ˙ q − ∂L ( q, ˙ q ) ∂q = Bu, (1) where the N × ( N − torque distribution matrix B is constantand has rank N − , and u is an ( N − -dimensional vectorassociated to the torques of the joint coordinates.Equation (1) can also be written in the form D ( q )¨ q + H ( q, ˙ q ) = Bu, (2)where D ( q ) is the (positive definite) inertia matrix, and H ( q, ˙ q ) groups the centrifugal, Coriolis and gravity terms.Letting x = [ q (cid:62) , ˙ q (cid:62) ] (cid:62) , standard calculations lead to a state-variable model ˙ x = f ( x ) + g ( x ) u, (3)which is affine in the control torques.For simplicity, impacts are assumed to occur with the swingfoot parallel to the walking surface. Let S := { x | z s = 0 , ˙ z s < } , (4)where z s is the height of the swing foot above the walkingsurface. It is assumed that S is a smooth (2 N − -dimensionalmanifold in the state space of the robot. The widely usedimpact model of Hurmuzlu [15] leads to an algebraic rep-resentation of the jumps in the velocity coordinates when theswing foot contacts the ground, and hence to a hybrid modelof the robot: ˙ x = f ( x ) + g ( x ) u, x − / ∈ S ,x + = ∆ ( x − ) , x − ∈ S , (5)where x + := lim τ (cid:38) t x ( τ ) (resp. x − := lim τ (cid:37) t x ( τ ) ) is thestate value just after (resp. just before) impact. The control-affine ordinary differential equation (ODE) describes the dy-namics of the swing phase, whereas the algebraic equationdescribes the impacts of the swing foot with the ground. Thisalgebraic equation can be decomposed into as follows: q + = ∆ q ( q ) , ˙ q + = ∆ ˙ q ( q ) ˙ q − . (6)The first matrix equation of (6) describes the relabelling ofthe coordinates at impact, while the second define the jumpvelocity coordinates. B. Virtual constraints and an exact reduced-order model
It is well understood in mechanics that a set of (regular)holonomic constraints applied to a Lagrangian model of amechanical system leads to another Lagrangian model thathas dimension equal to the dimension of the original modelminus twice the number of independent constraints. In the caseof mechanical systems, the constraints are imposed by a vectorof generalized forces (aka Lagrange multipliers) that can becomputed through the principle of virtual work. When com-puting the resulting reduced-order model, no approximationsare involved, and solutions of the reduced-order model aresolutions of the original model along with the inputs arisingfrom the Lagrange multipliers.Virtual holonomic constraints are functional relations on theconfiguration variables of a robot’s model that are achievedthrough the action of the robot’s joint torques and feedbackcontrol instead of physical contact forces. They are called virtual because they can be re-programmed on the fly in the control software. Like physical constraints, under certain regu-larity conditions, virtual constraints induce a low-dimensionalinvariant model called the zero dynamics [3], [31], [32]. Adetailed comparison of “physical constraints” versus “virtualconstraints” is provided in [13]. In the following, a briefoverview is given.A total of N − virtual constraints can be generated by the N − actuators. The virtual constraints are first expressedas outputs applied to the model (3), and then a feedbackcontroller is designed that asymptotically drives them to zero.One such controller, based on the computed torque, is givenin the Appendix; others can be found in [2]. Here, the virtualconstraints will be expressed in the form y = h ( q ) := q c − h d ( q f ) , (7)where q c ∈ Q c ⊂ R ( N − , and q f ∈ Q f ⊂ R with ( q c , q f ) forming a set of generalized configuration variables for therobot, that is, such that a diffeomorphism F : Q c × Q f → Q exists. The variables q c represent physical quantities that onewishes to “control” or “regulate”, while the variables q f remain“free”. Later, a special case of h d ( q f ) will be employed basedon a gait phasing variable , which makes it easier to interpretthe virtual constraints in many instances.When the virtual constraints are satisfied, the relation q c = h d ( q f ) leads to the mapping F c : Q f → Q by F c ( q f ) := F ( h d ( q f ) , q f ) (8)being an embedding; moreover, its image defines a constraintmanifold in the configuration space, namely ˜ Q = { q ∈ Q | q = F c ( q f ) , q f ∈ Q f } . (9)The constraint surface for the virtual constraints is Z : = { ( q, ˙ q ) ∈ T Q | y = h ( q ) = 0 , ˙ y = ∂h ( q ) ∂q ˙ q = 0 } = { ( q, ˙ q ) ∈ T Q | q = F c ( q f ) , ˙ q = J c ( q f ) ˙ q f , ( q f , ˙ q f ) ∈ T Q f } , (10)where J c ( q f ) := ∂F c ( q f ) ∂q f . (11)The terminology “zero dynamics manifold” comes from [3]. Itis the state space for the internal dynamics compatible with theoutputs being identically zero. For the underactuated systemsstudied here, the dimension of the zero dynamics manifold is , due to the assumption of two degrees of underactuation.Let B ⊥ be the left annihilator of B , that is, a × N matrixof rank such that B ⊥ B = 0 . Lagrange’s equation (1) thengives dd t (cid:0) B ⊥ ∂L∂ ˙ q (cid:1) = B ⊥ ∂L∂q , because B ⊥ Bu = 0 . The term B ⊥ ∂L∂ ˙ q is a form of generalizedangular momentum .Restricting the full-system to the the zero dynamics mani-fold, gives ˙ σ f = κ ( q f , ˙ q f ) , where σ f := M ( q f ) ˙ q f κ ( q f , ˙ q f ) := B ⊥ ∂L ( q, ˙ q ) ∂q (cid:12)(cid:12)(cid:12)(cid:12) q = F c ( q f )˙ q = J c ( q f ) ˙ q f . The invertibility of M ( q f ) is established in [13] under theassumption that the decoupling matrix used in (63) is full rank. z = (cid:20) z z (cid:21) := (cid:20) q f σ f (cid:21) gives ˙ z = (cid:20) M − ( z ) z ¯ κ ( z , z ) (cid:21) (12) =: f zero ( z ) , where ¯ κ ( z , z ) = κ ( z , M − ( z ) z ) . The properties of thisfour-dimensional system are determined by the Lagrangaindynamics of the full-order model and the choice of the virtualconstraints through the mapping (8). The form of the modelis similar to the one used to analyze planar systems with onedegree of underactuation [32, Remark 5.2]. C. Reduced-order hybrid model (hybrid zero dynamics)
Analogous to the full-order hybrid model (5) is the so-calledhybrid zero dynamics in which the zero dynamics manifoldmust satisfy ∆ ( S ∩ Z ) ⊂ Z . (13)This condition means that when a solution evolving on Z meets the switching surface, S , the new initial conditionarising from the impact map is once again on Z .The invariance condition (13) is equivalent to h ◦ ∆ q ( q ) (14) ∂h (¯ q ) ∂q (cid:12)(cid:12)(cid:12)(cid:12) ¯ q = ∆ q ( q ) ∆ ˙ q ( q ) ˙ q (15)for all ( q ; ˙ q ) satisfying h ( q ) = 0 , ∂h ( q ) ∂q ˙ q =0 (16) z s ( q ) = 0 , ˙ z s ( q, ˙ q ) < . (17)At first glance, these conditions may appear to be difficult tomeet. In the case of models with one degree of underactuation,however, it is known that if a single non-trivial solution of thezero dynamics satisfies these conditions, then all solutions ofthe zero dynamics will satisfy them [32, Thm. 5.2]. In thecase of systems with more than one degree of underactuation,systematic methods have been developed which modify thevirtual constraints “at the boundary” and allow the conditionsto be met [22]. Very straightforward implementations of theresult are presented in a robotics context in [4] and [11], andwe refer the reader to these paper for the details. With the invariance condition of the impact map, the re-duced order hybrid system, called hybrid zero dynamics (HZD),can be written in the form Σ zero : (cid:40) ˙ z = f zero ( z ) , z − / ∈ S ∩ Z z + = ∆ zero ( z − ) , z − ∈ S ∩ Z , (18)with ∆ zero := ∆ | S∩Z . It is proven in [32] that solutions of(18) lift to solutions of (5); indeed, (8) results in q ( t ) = F c ( z ( t ))˙ q ( t ) = J c ( z ( t )) ˙ z ( t )= J c ( z ( t )) M − ( z ( t )) z ( t ) . It is also proven in [32] that locally exponentially stableperiodic solutions of (18) lift to solutions of (5), which canbe locally exponentially stabilized.
Remark 1:
Low-dimensional pendulum models are approx-imate representation of the the swing phase dynamics of arobot, and normally when they are used, the impact mapis ignored. The zero dynamics is an exact low-dimensionalmodel that captures the underactuated nature of the robot,and moreover, the hybrid zero dynamics is an exact low-dimensional “subsystem” of the full-order hybrid model thatcaptures the underactuated dynamics of the hybrid model.III. S
ELF - STABILITY AND THE
HZDThe use of the virtual constraints allows the stability ofthe full-order model to be studied on the HZD. For a systemof two degrees of underactuation, the HZD is reduced to asystem of dimension four for any number of actuated joints.The associated Poincar´e map has dimension three. This paperseeks conditions for the eigenvalues of the Jacobian of thePoincar´e map to lie within the unit circle, thereby assuringlocal exponential stability of a periodic orbit. More precisely,this paper seeks guidelines for the selection of the virtualconstraints so that eigenvalues of the Jacobian lie within theunit circle.Once the virtual constraints are defined, we know how towrite a control law that drives them to zero and the walkinggait will be stable or not depending on the stability of theHZD. If the chosen virtual constraints induce stable walking,we will call this self-stability , in analogy with passive walkers.
A. Choosing the virtual constraints by physical intuition
Choosing the virtual constraints (7) means not only selecting h d ( q f ) , but also q f and q c . One objective of this work is togive a physical interpretation of stability of walking; thus, wechoose a set of variables q c , q f that are physically meaningful.Several choices seem possible. Since the effect of gravityappears to be crucial and the stance ankle is supposed to beunactuated in the sagittal and frontal planes, the free variablesare chosen as the position of the CoM in the horizontal plane x M , y M . Moreover, to make our results independent of thestep width D and length S , we consider normalized variables,namely, q f = ( X, Y ) (cid:62) , X = x M S and Y = y M D .Concerning the controlled variables, in principle, we wishto use meaningful coordinates that appear in both models that we will study in the next sections, that is, a simplified invertedpendulum model and a realistic model of Romeo. Since thehorizontal position of the CoM is used in q f , it makes senseto consider the vertical position z M as one of the coordinates.Moreover, we also include the Cartesian position of the swingfoot tip as a part of the controlled variables because thesevariables have a significant contribution to the change ofsupport. Hence, for the simplified model of Section IV, the setof controlled variables will be q c = ( z M , X s , Y s , z s ) where z s is the vertical position of the swing foot and X s , Y s arenormalized horizontal positions of the swing foot. For thecomplete model of Romeo and other realistic bipeds, q c willconsist of the variables just listed for the simplified model,together with a minimal set of body coordinates to accountfor the remaining actuators on the robot. For Romeo, this willmean adding the orientation of the swing foot, and remainingjoints in the torso and stance leg as described in Section V. B. From one degree of underactuation to two degrees ofunderactuation
By definition, the condition of transition from one step to thenext is defined by (4). In this condition, the vertical velocityof the foot being negative at impact ensures a well-definedtransversal intersection of the solution with the transition set.In the following discussion, we compare the set of possibleconfigurations of the robot at the change of support in thecase of one degree of underactuation versus two degrees ofunderactuation.In the context of the virtual constraints, the evolution ofthe swing leg is imposed as a function of q f . Thus, theswitching condition can be viewed as a condition on thefree configuration q f . We will define a switching configurationmanifold S : S = (cid:8) q f | z ds ( q f ) = 0 (cid:9) . (19)where z ds is the virtual constraint for the height of the swingleg.In the case of one degree of underactuation, the switchingconfiguration manifold contains two points corresponding tothe final value q − f of the periodic motion and the initial one q +f , and a proper impact only requires ˙ z < at the end ofthe step. For planar walking in the sagittal plane [32], thestability condition is determined completely by the periodicmotion itself, in other words, it is independent of the choiceof virtual constraints that induce it. Moreover, the stabilitycondition is physically meaningful: the velocity of the CoMat the end of the single support phase (for the periodic motion)must be directed downward [6].In the case of two degrees of underactuation, the manifold(19) depends on at least two free variables X and Y , andtherefore under mild regularity conditions it includes an in-finite number of points. Hence, the switching configurationmanifold of the robot right before the impact is representedby the solutions of z ds ( q f ) = z ds ( X, Y ) = 0 . (20)Moreover, for two degrees of underactuation, it has beenshown that a given periodic gait can be asymptotically stable or unstable depending on the choice of the virtual constraintsused to induce the motion [7]. We want to understand howthis happens.As a first observation, the choice of the virtual constraint, z ds , clearly affects the set of free variables belonging to the switching manifold : S = (cid:8) ( X, Y ) | z ds ( X, Y ) = 0 (cid:9) . (21)We will show in this paper that the shape of the switchingmanifold (21) has a crucial effect on the stability of walking.In particular, we show that • in the case of the Linear Inverted Pendulum (LIP), avery common simplified model for a humanoid robot,a condition to obtain self-synchronization between themodel’s motions in the lateral and sagittal planes isobtained as a condition on the tangent of S at the pointwhere it is intersected by the periodic orbit; and • when the LIP model is extended to an inverted pendu-lum model with a stance leg of variable length, self-synchronization plus a downward motion of the CoMyields self-stability of a periodic motion.The results are then extended to a realistic model of thehumanoid Romeo and illustrated in simulation.IV. S IMPLE MODEL : V
ARIABLE L ENGTH I NVERTED P ENDULUM
Two of the main features involved in dynamic walking arethe role of gravity and the limited torque available at the stanceankle before the foot rolls about one of its edges. Emphasizingthese aspects, the inverted pendulum model has been usedin numerous humanoid walking control applications, such as[18], [9], and [19].
A. Modeling of a Variable Length Inverted Pendulum
In recognition of this, the proposed methodology is firstpresented on a two-legged inverted pendulum, which consistsof two telescoping massless legs and a concentrated mass attheir intersection, called the hips; see Figure 2. The stanceleg is free to rotate about axes s and n at the groundcontact and the length of each leg can be modified throughactuation, allowing a desired vertical motion of the pendulumto be obtained. We also assume that the actuation of the swingleg will allow a controlled displacement of the swing leg endvia hip actuators and control of the leg length.The configuration of the robot is defined via the positionof the concentrated mass ( x M , y M , z M ) with respect to areference frame attached to the stance foot and the posi-tion of the swing leg tip, denoted by ( x s , y s , z s ) . Angularmomenta along axes s and n are denoted by σ x and σ y . In order to explore simultaneously the existence andstability of periodic orbits as a function of step length andwidth, a dimensionless dynamic model of the pendulum willbe used [6]. The normalizing scaling factors applied alongaxes s and n depend on desired step length S , desiredstep width D , and the mass m of the robot. Thus, a newset of variables is defined: ( X, Y, z M , X s , Y s , z s , σ X , σ Y ) =( xS , yD , z M , x s S , y s D , z s , σ x mD , σ y mS ) . a s n xz y Fig. 2. A simplified model of 3D biped robot. Each leg is massless and hasvariable length. At the contact point, the stance leg rotates passively aroundaxes s and n , the rotation around axis a is not considered since thisrotation is inhibited by friction in normal biped locomotion. The swing leghas a fully actuated spherical joint with respect to the concentrated mass ofthe hip. The variables ( z M , X s , Y s , z s ) will be the controlled vari-ables q c , whose evolution will be determined by virtual con-straints as a function of q f = ( X, Y ) .The height of the CoM and its vertical velocity are thenexpressed as functions of its horizontal position and velocity: z dM = f ( X, Y )˙ z dM = ∂f ( X, Y ) ∂X ˙ X + ∂f ( X, Y ) ∂Y ˙ Y (22)Since the legs are massless, the position of the swing legtip and its time derivative do not affect the dynamic modeland thus have no effect on the evolution of the biped duringsingle support. The moment balance equation of the pendulumaround the rotation axis s and n directly yields the equationof the zero dynamics, (cid:40) ˙ σ x = − mgy M , ˙ σ y = mgx M . (23)For this simple model, the angular momentum is simply (cid:40) σ x = m ˙ z M y M − mz M ˙ y M ,σ y = − m ˙ z M x M + mz M ˙ x M , (24)Using the normalized coordinates, the zero dynamics model(12) with z = σ f = [ σ X , σ Y ] (cid:62) can be written as ˙ q f = M − XY σ f , ˙ σ X = − gY, ˙ σ Y = gX, (25)where M XY = (cid:34) ∂f ( X,Y ) ∂X Y ∂f ( X,Y ) ∂Y Y − f ( X, Y ) − ∂f ( X,Y ) ∂X X + f ( X, Y ) − ∂f ( X,Y ) ∂Y X (cid:35) . We assume that once the tip of the swing leg hits the ground,that is, when z ∈ S ∩Z , the legs immediately swap their roles;thus, we assume that the double support phase is instantaneous.The state right before (after) the swapping of roles of thelegs is denoted by the exponent − ( + ). During the swappingof the support leg, the configuration of the robot is fixed, but the reference frame is changed since it is attached to the newstance leg tip. From the geometry, for normalized variableswe have: (cid:40) X + = X − − X − s ,Y + = − Y − + Y − s , (26)where the change of sign in the second equation correspondsto the change of direction of the axis n . Since the velocityof the CoM of a biped with massless legs is conserved in thetransition, we also have (cid:40) ˙ X + = ˙ X − , ˙ Y + = − ˙ Y − . (27)The design of the switching manifold S defined via z ds ( X, Y ) , and the transition condition (26) affected by X − , Y − defined via X ds ( X, Y ) and Y ds ( X, Y ) , are discussed in thenext section. B. Virtual constraints for the swing leg
Assume that a desired periodic motion is characterized byan initial free configuration for the single support phase X , Y and a final free configuration X f , Y f . Due to the normalizationof the variables with respect to the gait length and width, andthe desired symmetry between left and right support, we have X f − X = 1 and Y + Y f = 1 . Our initial objective is todefine a virtual constraint for the height of the swing foot, z ds ,as a function of X and Y . This virtual constraint takes thevalue for X = X , Y = Y and X = X f , Y = Y f , and ispositive in between during the single support phase.Inspired by [25], we propose z ds = Z se ( X, Y ) = ν z S a ( X, Y ) , (28)where ν z is a negative parameter that adjusts the height of theswing foot during a step and S a ( X, Y ) is S a ( X, Y ) = ( X − X a ) + CY − (( X − X a ) + CY ) (29)with X a = ( X f + X ) + C ( Y f − Y )2 . (30)With this choice the ellipse-shaped switching manifold isdefined by S = { ( X, Y ) | S a ( X, Y ) = 0 } , (31)and is illustrated in Figure 3. Remark 2:
The shape of the switching manifold is some-what arbitrary. We will see in the stability analysis done inthe following section that only the tangent of the switchingmanifold at X f , Y f will be used. The orientation of thistangent is characterized by the parameter C , and its effectwill be studied in detail.The virtual constraint imposed for the height of the swing footdefines the switching manifold and is crucial for the stabilityof a gait. The virtual constraint associated with the horizontalposition of the swing foot will allow us to impose where theswing foot is placed on the ground.In the case studied, the horizontal position of the swing footis controlled to obtain an ( X , Y ) -invariant gait; that is, the a s n (a) The mass motion starts on theswitching manifold on one stanceleg a s n (b) When the mass crosses the switch-ing manifold, the support leg ischangedFig. 3. The switching manifold (red ellipse) defines the position of the CoMwhere the swing leg hits the ground, with an instantaneous change of supportleg. starting position of the CoM with respect to the stance footwill be ( X , Y ) independent of any error in the previous step.Based on (26), the horizontal motion of the swing leg must besuch that when S a ( X, Y ) = 0 (i.e., when the leg touches theground), the desired position X − s = − X + X , Y − s = Y + Y is reached.To achieve this objective, at the end of the motion of theswing foot, the virtual constraints are defined by: X ds = X se ( X, Y ) = (1 − ν X S a ( X, Y )) ( − X + X ) ,Y ds = Y se ( X, Y ) = (1 − ν Y S a ( X, Y )) ( Y + Y ) . (32)With these constraints, for any values of ν X and ν Y , thedesired initial position of the CoM at the beginning of thenext step, that is X and Y , will be obtained. Moreover, inorder to increase the robustness with respect to variations ofthe ground height, the parameters ν X and ν Y are chosen suchthat in a periodic gait at the end of the step the horizontalvelocity of the swing leg tip is zero.We note that (32) expresses the desired horizontal positionof the swing foot at the end of the step. At the beginning ofthe step, for X = X , Y = Y , equation (32) obviously doesnot provide the initial position of the swing leg tip, X + S = − , Y + S = 1 for the periodic gait. To solve this issue, the virtualconstraints concerning the motion of the swing leg are in factdecomposed into two parts: one concerning near the end of thestep, more specifically, for X > X l = 0 . X + X f to the endof the step and one concerning the beginning of the step forwhich a five order polynomial function of X is added to thevirtual constraints of the second part. Therefore, the motionalong the axis s is imposed by the following equations: X ds = X se ( X, Y ) + P x ( X ) for X < X l ,X ds = X se ( X, Y ) for X ≥ X l .Y ds = Y se ( X, Y ) + P y ( X ) for X < X l ,Y ds = Y se ( X, Y ) for X ≥ X l . (33)The coefficients of the polynomial P x are such that at X = X l ,continuity in positions, velocities and accelerations are ensured(that is, P x ( X l ) = 0 , ∂P x ( X l ) ∂X = 0 , ∂ P x ( X l ) ∂X = 0 , P y ( X l ) =0 , ∂P y ( X l ) ∂X = 0 , ∂ P y ( X l ) ∂X = 0 .), and at the beginning ofthe step, the desired position and velocity are obtained, i.e., X − s = − , Y − s = 1 , ˙ X − s = 0 , ˙ Y − s = 0 . An illustration ofthe use of a polynomial function to ensure continuity of thevirtual constraint is shown in Figure 8.In fact, for the height of the swing foot, (28) is used alsoonly at the end of the step. At the beginning of the step, a poly-nomial function of X is also added. This modification allowsus to take into account that since the leg was previously insupport, its initial velocity is zero. A fourth order polynomialis used to ensure the continuity in position and velocity at themiddle of the step: z ds = Z se ( X, Y ) + P z ( X ) for X < X + X f ,z ds = Z se ( X, Y ) for X ≥ X + X f . (34)This slight modification does not change the switching surface. Remark 3 : The selection of the parameters X f , Y f , ˙ X f , ˙ Y f , ν X , and ν Y depends on our expected duration of the periodicstep. These parameters are found by solving a boundary valueproblem describing the periodicity of the gait. A detailedillustration of the swing leg trajectory is given in Section V-A. C. Gait with constant height and self-synchronization
An inverted pendulum with a constant height z = z isknown as a 3D linear inverted pendulum (3D LIP) [18]. Sincewith the condition z = z the motions in the sagittal andfrontal planes are decoupled, from (25), the expression of the3D LIP hybrid zero dynamics reduces to: ˙ X = σ Y z , ˙ Y = − σ X z , ˙ σ X = − gY, ˙ σ Y = gX, (35)which can be analytically solved.The notion of self-synchronization for the 3D LIP model,which was defined and characterized in [25], is briefly recalled.For a 3D LIP model, it is well known that the orbital energies: E ( X ) = ˙ X − gz X E ( Y ) = ˙ Y − gz Y (36)are conserved quantities during a step [18]. Moreover, thequantity L ( X, Y ) = ˙ X ˙ Y − gz XY (37)is also conserved and is called the synchronization measure[25]. If this quantity is zero, the motion between the sagittaland frontal planes is synchronized (i.e., ˙ Y = 0 when X = 0 [25]). Any periodic motion of the LIP, with symmetric motionfor left and right support, is characterized by X = − , X f = Y = Y f = 12 . (38)In normalized variables, the periodic motion for differentvalues of step durations T are presented in Figure 4. Fig. 4. Periodic motions in normalized variables for several values of T ,orientation of the final velocity is defined by < − ˙ Y ˙ X < . In [6], it is shown that for a periodic motion of the LIP,self-stabilization is impossible. In this section, our objectiveis to show the influence of the virtual constraints chosen forthe swing leg, and especially the influence of the switchingmanifold on the self-synchronization of the walking gait.Indeed, since synchronization between the sagittal and frontalmotions implies coupling between these two motions and sincethese motions are decoupled during the single support phase,it is essential to introduce a coupling at the transition.For the 3D LIP model, due to the specific values of X , X f , Y , Y f , the switching configuration manifold corre-sponding to the virtual constraints defined in (31) becomes: S = (cid:8) ( X, Y ) | X + CY − ( X + CY ) = 0 (cid:9) (39)because X a = 0 using (30) with (38).The virtual constraints for the swing legs are also designedsuch that ( X + , Y + ) = ( X , Y ) at the beginning of theensuing step, so the 3D LIP is said to be ( X , Y ) -invariant.In scaled coordinates, we limit our study to the case X = − and Y = .For an ( X , Y ) -invariant 3D LIP, if the switching manifoldis defined by (39), reference [26] shows that the Jacobianof the Poincar´e return map at the fixed point, expressedin the coordinate system ( X − k , L k , K − k ) , where L k is thesynchronization measure at the end of step k and K − k is thekinetic energy at the end of the step k , is in the form J = ∗ − λ L ∗ , (40)where ∗ represents non-zero terms. Three eigenvalues areidentified in this Jacobian. One is zero due to the fact that theposition of the CoM at the end of one step does not affect thedynamics of the following step for an ( X , Y ) -invariant 3DLIP. One is a unit eigenvalue which corresponds to the kineticenergy of the pendulum, which indicates that a variation ofthe kinetic energy will be conserved, and the system will notconverge to its original fixed point. We will propose in this paper a new expression for λ L ,and we will give an original condition on the choice of theswitching manifold to induce synchronization and a physicalinterpretation of this condition.Suppose that the state of the robot is slightly per-turbed around the initial periodic state: [ − / , / , ˙ X , ˙ Y ] .After i steps the initial state of the robot is denotedby [ − / , / , ˙ X i , ˙ Y i ] , and its synchronization measure isdenoted L i . After the following step the state will be [ − / , / , ˙ X i +1 , ˙ Y i +1 ] , and its synchronization measure isdenoted by L i +1 . The calculation of L i +1 L i is developed inthe appendix. The derivation is based on the conservation oforbital energies (36), the synchronization measure (37) duringa single support phase, and the fact that the change of supportoccurs on the switching manifold (39). For small variationaround the periodic orbit, the final error for each step satisfies: δX − i = − CδY − i . (41)Direct calculation gives: λ L = lim i →∞ L i +1 L i = ( ˙ Y − ˙ X )( C ˙ Y + ˙ X )( ˙ X + ˙ Y )( − C ˙ Y + ˙ X ) (42)Two factors appear in (42): one varies as a function of C and can be modified by the design of the switching manifold,while the other depends on the periodic gait velocity, i.e onthe step duration of the gait.When C ˙ Y + ˙ X = 0 , the synchronization occurs in onestep. This case corresponds to a switching line co-linear withthe initial velocity of the periodic motion since in this casewe have: C = − δX − i δY − i = − ˙ X ˙ Y (43)Thus, the cross product between the final error in position ( δX − i , δY − i ) and initial periodic velocity ( ˙ X , ˙ Y ) is zero.The values of C that ensure synchronization of the motionmust be such that: | λ L | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( ˙ Y − ˙ X )( C ˙ Y + ˙ X )( ˙ X + ˙ Y )( − C ˙ Y + ˙ X ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) < (44)For a periodic velocity corresponding to Figure 4, namely,when ˙ X > , ˙ Y < , ˙ X > − ˙ Y , it can be shown that C has to satisfy: < C < (cid:32) ˙ X ˙ Y (cid:33) . (45)In Figure 5, for the same initial state, the behavior of the LIPis shown for several values of C. Three cases are illustratedfor C satisfying (45) or not. For C = 0 . , since condition(45) is not satisfied, a periodic motion is not obtained, and thedirection of walking is not along the axis s . For the two othercases tested, the condition (45) is satisfied, and convergenceto a periodic motion is observed, but in the two cases, theperiodic gait velocities differ. In the case C = 1 . , since thevalue of C is close to (43), the synchronization is faster than inthe case C = 1 . . It should be noted that the values of ˙ X , ˙ Y correspond to the periodic motion. In the case C = 1 . , the (a) (b)(c) (d)(e) (f)Fig. 5. Starting from the same initial state X = − . , Y = 0 . , ˙ X =2 . s − , ˙ Y = − . s − , the simulation of 10 steps shows severalbehaviors depending on the switching manifold defined by C . From top tobottom, we have C = 0 . ; C = 1 . , C = 1 . . The left-hand figuresshow the position of the stance foot, the evolution of the CoM, and a partof the switching manifold are shown. On the right-hand figures, the state ofthe HZD is shown in the Poincar´e section, just before the change of support:divergence or convergence is clearly illustrated. LIP motion converges to a motion such that ˙ X = 2 . s − , ˙ Y = − . s − , thus, − ˙ X ˙ Y = 1 . . In the case C =1 . , the LIP motion converges to a motion such that ˙ X =2 . s − , ˙ Y = − . s − , thus, − ˙ X ˙ Y = 1 . .In Figure 6 the level contours of | λ L | are drawn as a functionof the parameter C and the step duration T . The condition | λ L | < ensures convergence toward a synchronized motion,thus is a condition of self-synchronization. However, thewalking velocity of the gait is not controlled in the sense thata perturbation of such a gait will result in the convergenceof the gait back to a periodic motion, but with a differentgait velocity as illustrated in Figure 5. In the next section, wewill see how judiciously adding oscillations of the CoM canstabilize the walking velocity. D. VLIP model and self-stabilization
From study of planar walking robots, it is known thatvertical oscillations of the CoM can asymptotically stabilizea periodic walking gait [32]. Thus, here, in order to stabilize . . . . . . . . . . . . . . . . . . . . . . . . T C Fig. 6. λ L level contours expressed as function of C and step duration T for the LIP model and z = 0 . m. Xz Y
Fig. 7. Surface defined by the virtual constraint for the vertical motion themass motion z d ( X, Y ) . The red curve gives an example of the mass evolutionduring a single support phase. the walking velocity in an ( X , Y ) -invariant gait, oscillationsof the CoM will be introduced. These oscillations of the CoMare obtained by the following virtual constraint: z d ( X, Y ) = z − aS a ( X, Y ) , (46)with S a ( X, Y ) defined in (29) associated to the ellipse-shapedswitching manifold defined in (31)The expression of S a ensures that at the beginning and endof the step z d = z since S a ( X , Y ) = S a ( X f , Y f ) = 0 .Note that the case a = 0 with X a = 0 corresponds to the 3DLIP example described in Section IV-C.During the single-support phase, the horizontal positionof the mass is located inside the ellipse; thus, S a ( X, Y ) is negative and increasing when approaching the switchingmanifold ellipse. Choosing a > will ensure a negativevertical velocity of the mass at the transition (see Figure 7).This negative vertical velocity of the CoM right before thechange of support implies that the angular momenta σ X and σ Y decrease at the change of support. In order to obtain aperiodic motion, the angular momenta must therefore increase a Z X Fig. 8. Planar illustration of the modification of the virtual constraint atthe beginning of each step to avoid discontinuity in tracking error. At thebeginning of the step, the virtual constraint, in black, is continuous for theheight of the CoM, but the velocity is not continuous. A modification of thevirtual constraint z cor ( X ) is added to satisfy the conditions of continuity, theconstraint effectively used is shown in red. during the stance phase, and thus it is necessary to slightlyshift the relative position of the support leg and the CoM [5].The position of the CoM at the beginning of the step is thenwritten as X = − + D X , Y = − D Y and at the end of thestep X f = + D X , Y f = + D Y with D X = X + X f ≥ and D Y = Y f − Y ≥ .This choice of the switching manifold only ensures thecontinuity of the vertical position of the CoM. To ensure thecontinuity of the vertical velocity of the CoM as well, a thirdorder polynomial function of X , denoted by z cor ( X ) , is addedto expression (46) for z d . This function is null at the beginningof the step, where X = X + , because the continuity of theposition is already ensured, but its derivative should compen-sate the difference between the vertical velocity at the end ofthe previous step and the one corresponding to the reference ˙ z d ( X , Y ) in (46). This correction term only acts from thebeginning of the step where X = X + to the mid-step where X = D X . In the second half of the step, where X ≥ D X , z cor ( X ) ≡ , and smoothness in the mid-step is guaranteedby imposing the conditions z cor ( D X ) = ˙ z cor ( D X ) = 0 . Thistechnique, which has already been used in previous studies[22], [7], is illustrated in Figure 8.In the presence of vertical oscillations, periodic motionswill exist only for an appropriate set of virtual constraints.In particular, the values of D X and D Y cannot be chosenarbitrarily. In practice, the values of C , a , and T , the durationof the periodic step, are chosen, and D X , D Y , and theangular momentum at the end of the single support (or thefinal velocity of the CoM ( ˙ X + , ˙ Y + ) are deduced by solvinga boundary value problem describing the periodicity of theevolution of the CoM. The variations of D X and D Y asfunctions of a for fixed values of C and T are shown in Figure9(b). E. Stability of the gait obtained
The study of the gait stability is undertaken through numer-ical calculations. The evolution of the mass in single supportis integrated numerically from the dynamic model (25) takinginto account the conditions (46) and (31) developed in SectionIV. a e e e (a) a D X D Y (b)Fig. 9. Influence of the amplitude of vertical oscillation characterized by a on the stability via the magnitude of the three eigenvalues e , e , e (a) andon the support foot position shift (b). Illustrations are numerically obtainedfor the case z = 0 . m, C = 1 . and T = 0 . s. To study the effect of oscillations on stability, we first revisitthe 3D LIP case ( a = 0 ). The virtual constraints are completelydefined by the value of z and C . We fix z = 0 . m and will study the effect of the parameter C . Many periodicorbits exist for a = 0 corresponding to various step durations T , however, for a given value of T , a unique periodic gaitexists. As explained in Section IV-C, three eigenvalues areevaluated to examine the stability of the periodic motion. Inour numerical study, the Poincar´e return map is expressedusing our state variables χ k . Because the eigenvalues of theJacobian of the Poincar´e map are invariant under a changeof coordinates, the results that we observed are identical tothose of (40): among the three eigenvalues, one is equal tozero, one is equal to 1, and the last one is λ L . Selecting C and step duration T such that | λ L | < according to(45) will allow us to observe the effect of oscillations onthe eigenvalues. Figure 9(a) illustrates the evolution of thethree eigenvalues e , e , and e for C = 1 . and T = 0 . sas the oscillation amplitude parameter a increases. In ourexample, the magnitudes of the two non-null eigenvaluesdecrease as the amplitude of the oscillations increases. Thus,when a > , the absolute values of all eigenvalues aresmaller than 1. Consequently, we observe that introducingoscillations transforms the self-synchronization property of thewalking gait into self-stabilization. Figure 10 shows the effectof oscillations on stability; the stability criterion δ (defined asthe maximum magnitude of the eigenvalues) is shown via levelcontours for several values of the ellipse parameter C and agiven step duration T = 0 . s. The value δ decreases when a increases. We observe that when oscillations are introduced,the values of the ellipse parameter C for which stability existsare close to the ones where self-synchronization is observedfor the 3D LIP. Figure 11 shows the evolution of the stabilitycriterion δ for a = 0 . m. We can see that the stability effectobserved in our first example (Figure 9(a)) is extended to mostof the values of the parameter C and step duration T . Tocomplete this analysis, Figure 9(b) shows the evolution of theparameters D X and D Y , which introduce asymmetry in thegait, obtained for a periodic motion as the amplitude of theoscillations increases. As expected, the values of D X and D Y increase with oscillations to compensate the loss of angularmomentum at the support change.Figure 12 illustrates the evolution of the CoM in the hor- Fig. 10. Stability criterion level contours for the pendulum model as functionof C and amplitude parameter a for T = 0 . s and z = 0 . m.Fig. 11. Stability criterion level contours for the pendulum model as functionof C and step duration T for a = 0 . m and z = 0 . m. The red curvesrepresent the boundaries inside which | λ L | < (equation (45)) for the LIP. izontal, sagittal, and frontal planes for several step durationswith a = 0 . m. A figure-eight-shaped evolution of the CoMcan be observed in the frontal plane, which is qualitativelysimilar to the motion of the CoM of a human during walking[10].V. E XTENSION TO THE CONTROL OF A HUMANOID ROBOTMODEL
We saw in Section IV-E that by introducing oscillations ofthe CoM and with an appropriate change of support, it is pos-sible to generate asymptotically stable periodic walking gaitsfor the inverted pendulum model. In this section, we applythe same principles to a more complex model of a bipedalrobot. In particular, we consider a general n -DoF humanoidrobot with six joints per leg, where link mass parameters andactuator inertias are taken into account. Numerical simulationwill then allow us to discuss the similarities and differences X -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Y (a) X -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 z (b) Y z (c)Fig. 12. Motion of the mass in horizontal (a), sagittal (b), and frontal (c)planes ( a = 0 . , C = 1 . , z = 0 . ). Magenta, cyan, black, and blackhave periods of T = 0 . s , T = 0 . s , T = 0 . s, and T = 1 s, respectively.Green dots are initial positions and red dots final positions for one step. • the transition condition from step to step as defined by the C parameter of the switching manifold defined in (29);and • vertical oscillations of the CoM as defined by the param-eter a . A. A periodic walking gait for the humanoid robot Romeo
As we would like to reproduce the results obtained for thependulum, the same virtual constraints used previously andsummarized here will be used: • The height of the center of mass is given by: z d = z − aS a ( X, Y ) + P ( X ) for X < DX,z d = z − aS a ( X, Y ) for X ≥ DX (47) • The height of the swing foot is given by: z ds = ν z S a ( X, Y ) , + P z ( X ) for X < DX,z ds = ν z S a ( X, Y ) for X ≥ DX (48) • The horizontal position of the swing foot is given by: X ds = (1 − ν X S a ( X, Y )) (1 / − D X + X ) + P x ( X )for X < X l ,X ds = (1 − ν X S a ( X, Y )) (1 / − D X + X ) for X ≥ X l .Y ds = (1 − ν Y S a ( X, Y )) (1 / D Y + Y ) + P y ( X )for X < X l ,Y ds = (1 − ν Y S a ( X, Y )) (1 / D y + Y ) for X ≥ X l . (49)with S a ( X, Y ) = ( X − D x − CD y ) + CY − (( − / − CD y ) + C (1 / D Y ) .Here, for simplicity, the remaining DoFs are fixed to re-alistic constant values . The orientation of the swing foot iscontrolled to be constant (the foot remains flat with respect tothe ground). The upper-body joints are constrained to fixedpositions and the orientation of the torso is controlled tobe zero (in each of its Euler angles). The fixed upper-bodyconfiguration is visible in Figure 13. While the chosen valuesmay have an influence over the gait stability, examining theirpotential effects is out of the scope of this study.In contrast to the simple pendulum model, it is necessaryto take into account the step length and width in the designof the constraints due to limits in the robot’s workspace andbecause the motion of the swing leg affects the dynamics ofthe robot. Here, we study a fixed step length S = 0 . m andwidth D = 0 . m .The periodic motion of the robot is imposed by the param-eters defining the virtual constraints and the HZD model (18). In practice, we would use optimization to define them. Value Effecta several amplitude of oscillations of CoMC several shape of the switching manifoldT several duration of the single support z ν z X l . D X limit to adapt the horizontalmotion of the swing footTABLE IT HE FIXED PARAMETERS
Effect D X shift of the CoM position along axis x D Y shift of the CoM position along axis y ν x such that the horizontal velocities of ν y the swing leg tip are null at the end of SSTABLE IIT HE OPTIMIZED PARAMETERS
Part of the parameters are chosen arbitrarily, and some param-eters are defined by optimization to guarantee the existence ofa periodic motion compatible with the under-actuation of therobot satisfying some constraints such as duration of the step.Table I defines the chosen parameters and their effect, whileTable II presents the optimized parameters.An example of a walking gait obtained for the robot isillustrated in Figure 13. In this example, C = 1 , T = 0 . s ,and a = 0 . m . A periodic orbit of the system is found bycalculating a fixed point of the Poincar´e return map. Here wehave D X = 0 . , D Y = 0 . , σ − X = 20 . kg.m.s − ,and σ − Y = − . kg.m.s − .From the obtained periodic orbit we observe that the rangeof motion of the CoM in the frontal plane is relativelysmall. Moreover, the obtained walking gait is highly dynamicas the angular momentum in the sagittal plane is relativelyhigh, which on the other hand, makes a small initial velocityinsufficient for initiating the walking gait.We notice that when the walking velocity decreases, theamplitude of the lateral motion of the CoM increases, and theworkspace limitation for the leg prevents finding a periodicmotion for the chosen value of z . Z X -0.2 0 0.2 0.4 Z Y -0.1 0 0.1 0.2 0.3 Fig. 13. A periodic motion obtained with the proposed approach for C = 1 , T = 0 . s , a = 0 . m . . . . . . . . . . . . . . . . . . S = 0.3 m a=0.03 z0=0.65 T C Fig. 14. Stability criterion level contours for the complete model as a functionof C and step duration T for a = 0 . m , S = 0 . m. B. Stability of walking
Our main objective is to discover to what extent the effectsof C and T on the stability of walking for a pendulum modelcan be extrapolated to the stability of walking of a realistichumanoid. Figure 14 shows the evolution of the stabilitycriterion (maximum eigenvalue magnitude of the Jacobian ofthe Poincar´e return map) as a function of C and T for verticaloscillations with a = 0 . m . We observe that the shape ofthe stable area as a function of parameters C and T is quitesimilar to the one obtained for the pendulum model but witha shifting of the area of stability toward smaller values of C . C. Comparison between a realistic and a simplified model
To properly investigate the effect of having a distributedbody for the trunk and other links versus having a point massmodel, as in the inverted pendulum, we study a simplifiedmodel of the robot by concentrating all the mass of the robotto a point of the trunk. This point is placed in a position suchthat with straight legs, the CoM of the robot and this simplifiedmodel coincide. Thus, a model of the robot with similarkinematics but a different mass distribution is considered. Thestable area as a function of C and T is presented in Figure15. We observe that compared to the stability regions of therealistic model of Romeo, which are shown in Figure 14, thereis a shifting of the stability areas along the C -axis and aslight modification of the stability margin. Aside from thesedifferences, the stability regions for the two models appearquite similar. As a conclusion, the analysis of the ideal modelseems to be an appropriate tool to understand the key featuresleading to a stable gait, but in order to be able to implementthe control laws in practice, studying a more realistic model,like the one presented based on virtual constraints and HZD,is essential to obtain appropriate numerical values for the keyparameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simplified model of Roméo T C Fig. 15. Stability criterion level contours for the simplified model of Romeo,where a point mass replaces the distributed mass, as functions of C and stepduration T for a = 0 . m , S = 0 . m . D. Effect of the Impact
A main difference between the pendulum model and thecomplete humanoid model is the impact effect at supportchange. In the pendulum model, as the legs are massless, noimpact is transmitted to the CoM when changing the supportleg. Thus, the contact of the swing leg with ground does notaffect the pendulum motion. In the complete robot model,however, legs masses are considered, and a rigid impact ismodeled as the swing leg touches the ground. This impacthas a direct effect on the vertical velocity of the CoM ofthe humanoid robot when changing support and thus on thestability of the gait [26]. More specifically, loss of kineticenergy at the transition between two steps due to the impactor reduction of angular momentum at the transition [6] mayresult in stable gaits. To illustrate the contributions of thesetwo mechanisms on stability, the stable areas as functions of C and T are demonstrated in Figure 16. It can be clearly seenthat an increase of the vertical oscillations via parameter a leads to an expansion of the stable area and a reduction ofthe maximal eigenvalue. It can be noted that for a value of a = 0 . m , the real vertical oscillations of the CoM is onlyaround . m as it can be seen in Figure 20. In Figure 16, itis also clear that when the impact is reduced, via a reductionof the vertical velocity of the swing leg tip (induced by areduction of the height of the swing leg motion right beforeimpact), the area corresponding to the stable motion shrinksand the maximal eigenvalue increases. In the case studied, thehorizontal velocity of the swing leg tip at impact is chosen tobe null. E. Simulation of the walking gait
To illustrate the robustness of the proposed approach withrespect to the configuration and velocity of the robot, thissection presents some different results obtained with the com-plex model of the humanoid robot. Starting from rest, the firstdesired gait corresponds to a walking gait with S = 0 . m , . . . . . . . . . . . . . . . S = 0.3 m a=0.01, z0=0.65 T C (a) . . . S = 0.3 m a=0.01, impact reduced z0=0.65 T C (b) . . . . . . . . . . . . . . . . . . S = 0.3 m a=0.03 z0=0.65 T C (c) . . . . . . S = 0.3 m a=0.03, impact reduced z0=0.65 T C (d) . . . . . . . . . . . . . . . . . . . . . . S = 0.3 m a=0.05 z0=0.65 T C (e) . . . . . . . . . . S = 0.3 m a=0.05, impact reduced z0=0.65 T C (f)Fig. 16. Stability criterion level contours for the complete model as functionof C and step duration T for several amplitudes of oscillations and impacts.For top to down, the amplitude of oscillations increases a = 0 . m ; a =0 . m ; a = 0 . m . From left to right, the amplitude of impact is dividedby 2. D = 0 . m , and T = 0 . s for 30 steps, and then a transitionto a faster walking gait is imposed with T = 0 . s and thesame step sizes for 30 others steps.The robot starts its motion in double support with the twolegs in the same frontal plane. In double support, an initialforward motion of the CoM allows the humanoid to increaseits angular momentum around the axis n such that the barrierof potential energy can be overcome during the first step. Thisinitial motion places the CoM in front of the ankle axis andlaterally closer to the next stance foot. This choice of the CoMposition at the beginning of the single support is such thatduring the first single support phase the contribution of gravityto the dynamics will be coherent with the evolution expectedduring the periodic motion. Based on this initial state, thevirtual constraints are built as presented in Section IV-B. Theinitial configuration of the robot and a representation of thefirst steps in the sagittal plane are shown in Figure 17, wherethe schematic of the robot is presented at each transition withthe trace of the swing leg tip and the CoM motions.The convergence of the solution to the desired periodic orbitcan be illustrated in various ways. In Figure 18, the evolutionof the state error for the unactuated variables (defined by thedifference between the current state and the one corresponding Z X Fig. 17. Illustration of the first steps of the ROMEO gait starting from doublesupport. to the desired periodic motion) are shown; the convergence ofthe error to zero is clear. It can be observed that the conver-gence to the periodic motion is faster for the first motion (with T = 0 . s ) compared to the second one ( T = 0 . s ). Thisresult is coherent with the value of the maximal eigenvaluefor the two cases illustrated in Figure 16(e).We notice that a faster convergence can be obtained for theerror on position than for the error on angular momentum (orvelocity). This can be interpreted as a fast synchronizationbetween the motion in the frontal and sagittal planes followedby a stabilization to the desired motion in terms of velocity.This interpretation is reinforced by observing the evolution ofthe CoM in a frame attached to the stance foot as presented inFigure 19, where the evolutions corresponding to the first foursteps are numbered. It can be observed that since the swingleg touches the ground with an imposed distance with respectto the position of the CoM, the initial position of the CoMat the consecutive steps is constant (except for the first step).During the first steps a synchronization is observed due tothe choice of the transition between steps as shown in sectionIV-C. Then there is a slower convergence toward the exceptedperiodic motion while synchronization is conserved.We note that even though in the two periodic motionscorresponding to T = 0 . s and T = 0 . s the CoM travelsthe same distance in the sagittal plane with S = 0 . m , thedistances travelled by the CoM in the frontal plane by theCoM are different.The evolution of the CoM drawn in the frontal plane in afixed world frame is shown in Figure 20. It can be noticed thatthe vertical oscillation of the CoM is limited to approximately . m for a = 0 . m ; this amplitude is less than the valueobserved in human walking.We can also observe that while the walking gait starts witha foot placed on the point (0 , , in a fixed world frame,and the desired distance between feet in the frontal plane is D = 0 . m , the mean value of the lateral oscillations of theCoM for the last steps is different from D/ (see Figure 20).This can be explained by a step width and length that can be Step Number
10 20 30 40 50 60-3-2-1012345
State Error in Poincaré Section
100 X100 Y < X < Y Fig. 18. The convergence toward the desired periodic gait is illustrated viathe errors on the unactuated state
X, Y, σ X , σ Y in the Poincar´e section asfunction of the step number. Due to the scale used, the errors in position aremultiplied by on the graphic. X [m] -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 Y [ m ] COM evolution top view
Fig. 19. The evolution of the CoM of the robot in shown in a top view ina frame attached to the stance leg and lateral axis directed toward the CoM.Each step starts with the desired initial position of the CoM. The first stepsare numbered, and illustrate a self-synchronization of the walking gait asobserved in [26] then the convergence toward the desired motion is obtained.The green and red curves correspond respectively to the desired motions with T = 0 . s and T = 0 . s . slightly different from S or D during the transient phases. Acorrection of the path followed can be easily obtained basedon [28].The convergence toward the desired motion is also illus-trated in Figure 21 for the motion of one joint (the knee)described in its phase plane and for the evolution of theduration of the step versus the step number.VI. C ONCLUSIONS
With a control law based on virtual constraints and a planarrobot with one degree of underactuation, it is known that thestability of a walking gait depends only on the periodic orbititself: if the velocity of the CoM at the end of the singlesupport phase is directed downward, the gait is stable. For X [ m ] Y [ m ] COM Evolution in the Frontal Plane
Fig. 20. The convergence is shown for the evolution of CoM in a frontalview in a fixed frame. Starting for the state denoted by a cyan star, the motionconverge toward the green curve corresponding to the periodic motion with T = 0 . s and then to the red curve corresponding to periodic motion with T = 0 . s . a
10 20 30 40 50 60 S t e p D u r a t i o n [ s ] Step Duration (a)
Knee Angle [ rd ] -1 -0.9 -0.8 -0.7 -0.6 K n ee V e l o c i t y [ r d = s ] -3-2-10123 Phase Plane for One Knee (b)Fig. 21. The convergence is shown for the duration of the steps as functionof the number of steps in (a) and by the evolution of one knee motion drawnin its phase plane (a). Starting for the state denoted by a cyan star, the motionconverges toward the green curve corresponding to the periodic motion with T = 0 . s and then to the red curve corresponding to periodic motion with T = 0 . s . a 3D humanoid robot with two degrees of underactuation,the condition of stability of a walking gait no longer dependsonly on the periodic motion but also on the choice of virtualconstraints used to create the walking gait. Moreover, theswitching manifold plays a crucial role in the synchronizationbetween motion in the sagittal and frontal planes, and thus onstability.This paper studied the influence of the horizontal positionof the CoM at the transition from one step to the next and ofvertical oscillations of the CoM within a step on the inherentstability of a walking gait. Through internal constraints on asimplified two-legged pendulum model, we were able to obtainasymptotically stable periodic gaits. The stability propertiesemerge from the definition of appropriate virtual constraints.In particular, the walking gait does not require stabilizationthrough high-level control commands.We were able to transfer these properties to a completerobot model and observe qualitatively similar results. In doingthis, we did encounter new dynamic effects due to energy lossassociated with impacts at leg swapping in the complete model which was not considered in the simple pendulum model. Ithas to be noted that even if the principles of the approachtested on the simple model can be transferred to the complexmodel of the robot, the values of our key parameters C have tobe slightly modified to take into account the dynamics of therobot. As a consequence, a study of the Hybrid Zero Dynamicscorresponding to the specific robot considered is useful forobtaining stable walking gaits.Finally, a complete model of a humanoid has many degreesof freedom that were not exploited in this study. Here, theupper body and torso joints have been frozen to arbitraryvalues. In future work, we can exploit the upper body joints—and especially the orientation of the trunk—to increase theenergy efficiency of the gaits. Once this is done, our algorithmwill be tested experimentally on the humanoid robot Romeo.A PPENDIX
The control law in single support
From the definition of the virtual constraints, the output andits derivatives are y = q c − h d ( q f ) , (50) ˙ y = ˙ q c − ∂h d ( q f ) ∂q f ˙ q f , (51) ¨ y = ¨ q c − ∂h d ( q f ) ∂q f ¨ q f − ∂∂q f (cid:18) ∂h d ( q f ) ∂q f ˙ q f (cid:19) ˙ q f . (52)The control objective is to ensure that the components ofthe output vector y converge to zero sufficiently rapidly. Letdenote ν ( y, ˙ y ) such a control law, yielding ¨ y = ν ( y, ˙ y ) . (53)A special case could be ν ( y, ˙ y ) = − ( K d (cid:15) ˙ y + K p (cid:15) y ) for some K p and K d positive definite and (cid:15) > [22].In the following, we show how one can obtain the controllaw (i.e., the actuator torques) ensuring the satisfaction ofequation (53).The joint coordinates, velocities and accelerations of therobot can be expressed as function of controlled and freevariables: q = F ( q c , q f ) , (54) ˙ q = ∂F ( q c , q f ) ∂q c ˙ q c + ∂F ( q c , q f ) ∂q f ˙ q f , (55) ¨ q = ∂F ( q c , q f ) ∂q c ¨ q c + ∂F ( q c , q f ) ∂q f ¨ q f + Ψ( q c , q f , ˙ q c , ˙ q f ) , (56)where Ψ( q c , q f , ˙ q c , ˙ q f ) = ∂∂ ( q c , q f ) (cid:18) ∂F ( q c , q f ) ∂ ( q c , q f ) (cid:20) ˙ q c ˙ q f (cid:21)(cid:19) (cid:20) ˙ q c ˙ q f (cid:21) . Using equations (52) and (53), it follows that in the closed-loop system, the accelerations of the controlled variablessatisfy the equation ¨ q c = ∂h d ( q f ) ∂q f ¨ q f + ∂∂q f (cid:18) ∂h d ( q f ) ∂q f ˙ q f (cid:19) ˙ q f + ν ( y, ˙ y ) , (57)where y and ˙ y are computed in (50) and (51). The actuatortorques able to achieve this desired closed-loop behavior can be computed based on the dynamic model(2). The torquevector Γ allowing the n − controlled variables to followthe desired closed-loop behavior (57), using (54), (55) (56)and (2), satisfies the following equation: ¯ D ( q c , q f ) J r ( q c , q f )¨ q f + Ω( q c , q f , ˙ q c , ˙ q f , ν ( y, ˙ y )) = B Γ , (58)where ¯ D ( q c , q f ) = D ( q ) (cid:12)(cid:12)(cid:12)(cid:12) q = F ( q c ,q f ) ,J r ( q c , q f ) = ∂F ( q c , q f ) ∂q c ∂h d ( q f ) ∂q f + ∂F ( q c , q f ) ∂q f , Ω = ¯ D ( q c , q f (cid:20) ∂F ( q c , q f ) ∂q c (59) (cid:18) ∂∂q f (cid:18) ∂h d ( q f ) ∂q f ˙ q f (cid:19) ˙ q f + ν ( y, ˙ y ) (cid:19) + Ψ( q c , q f , ˙ q c , ˙ q f ) (cid:21) + ¯ H ( q c , q f , ˙ q c , ˙ q f ) , ¯ H ( q c , q f , ˙ q c , ˙ q f ) = H ( q, ˙ q ) (cid:12)(cid:12)(cid:12)(cid:12) q = F ( q c , q f )˙ q = ∂F ( q c ,q f ) ∂q c ˙ q c + ∂F ( q c ,q f ) ∂q f ˙ q f Multiplying (58) on the left by the full rank matrix (cid:20) B ⊥ B + (cid:21) , where B + is the pseudo inverse of B and B ⊥ is a × n matrixof rank such that B ⊥ B = 0 , results in B ⊥ ¯ D ( q c , q f ) J r ( q c , ˙ q f )¨ q f + B ⊥ Ω( q c , q f , ˙ q c , ˙ q f , y, ˙ y ) = 0 (60) B + ¯ D ( q c , q f ) J r ( q c , q f )¨ q f + B + Ω( q c , q f , ˙ q c , ˙ q f , y, ˙ y ) = Γ . (61)It follows that a feedback control law ensuring (53) can beobtained from the following equations Γ = B + ¯ D ( q c , q f ) J r ( q c , q f ) v + B + Ω( q c , q f , ˙ q c , ˙ q f , y, ˙ y ) (62) v = − (cid:0) B ⊥ ¯ D ( q c , q f J r ( q c , q f ) (cid:1) − B ⊥ Ω( q c , q f , ˙ q c , ˙ q f , y, ˙ y ) . (63)where v denoted the acceleration ¨ q f satisfying the equation(60). This form of the feedback law only requires the inver-sion of a matrix whose size corresponds to the number ofunactuated coordinates, here a 2 × Condition of synchronization for the LIP model
This section details of the condition of synchronization forthe LIP model and the proposed control strategy, with thevirtual constraints proposed in section IV-B.The initial state of the robot after step i is written as: X + i = − / , (64) Y + i = 1 / (65) ˙ X + i = ˙ X + δ ˙ X + i (66) ˙ Y + i = ˙ Y + δ ˙ Y + i (67) Due to the discrete invariance of the control law, there is noerror on position. At the end of the step the state of the robotwill be denoted X − i = 1 / δX − i (68) Y − i = 1 / δY − i (69) ˙ X − i = ˙ X + δ ˙ X − i (70) ˙ Y − i = − ˙ Y + δ ˙ Y − i (71)Using the fact that ( ˙ X and ˙ Y ) and L i = 0 define asynchronized motion and neglecting the second order terms,we obtain: L i ≈ ˙ X δ ˙ Y + i + ˙ Y δ ˙ X + i (72)Since the velocity of the CoM is conserved through change ofsupport for the LIP model: δ ˙ X + i +1 = δ ˙ X − i δ ˙ Y + i +1 = − δ ˙ Y − i and L i +1 can be expressed as L i +1 ≈ ˙ X δ ˙ Y + i +1 + ˙ Y δ ˙ X + i +1 (73)or L i +1 ≈ − ˙ X δ ˙ Y − i + ˙ Y δ ˙ X − i (74) L i +1 L i ≈ − ˙ X δ ˙ Y − i + ˙ Y δ ˙ X − i ˙ X δ ˙ Y + i + ˙ Y δ ˙ X + i (75)We will now express the final error in velocity as function ofthe initial error for the step i . The orbital energies, E x and E y and synchronization measure L are conserved quantities.Therefore, we have, (cid:16) ˙ X − i (cid:17) − ω (cid:0) X − i (cid:1) = (cid:16) ˙ X + i (cid:17) − ω (cid:0) X + i (cid:1) (76) (cid:16) ˙ Y − i (cid:17) − ω (cid:0) Y − i (cid:1) = (cid:16) ˙ Y + i (cid:17) − ω (cid:0) Y + i (cid:1) (77) ˙ X − i ˙ Y − i − ω X − i Y − i = L i (78)or (cid:16) ˙ X + δ ˙ X − i (cid:17) − ω (cid:0) / δX − i (cid:1) = (cid:16) ˙ X + δ ˙ X + i (cid:17) − ω / (cid:16) − ˙ Y + δ ˙ Y − i (cid:17) − ω (cid:0) / δY − i (cid:1) = (cid:16) ˙ Y + δ ˙ Y + i (cid:17) − ω / (cid:16) ˙ X + δ ˙ X − i (cid:17) (cid:16) − ˙ Y + δ ˙ Y − i (cid:17) − ω (cid:0) / δX − i (cid:1) (cid:0) / δY − i (cid:1) = L i . Using these equations and neglecting the second order termswe obtain: δ ˙ X − i ≈ ω X δX − i + δ ˙ X + i (79) δ ˙ Y − i ≈ − ω Y δY − i − δ ˙ Y + i (80) ω (cid:0) δX − i + δY − i (cid:1) ≈ ˙ X δ ˙ Y − i − ˙ Y δ ˙ X − i − L i . (81)The last equation (81) combines with the previous ones (79)(80) can give us a relation between δX − i , δY − i and L i . ω (cid:0) δX − i + δY − i (cid:1) ≈ − ω ˙ X Y δY − i − ω ˙ Y X δX − i − L i (82) (cid:32) Y ˙ X (cid:33) δX − i + (cid:32) X ˙ Y (cid:33) δY − i ≈ − L i ω . (83)Equation (83) confirms that synchronized motion occurs with δX − i = δY − i = 0 and in this case via (79) and (80), we seethat the error on velocity is conserved. Now we will write theratio L i +1 L i as function of C . The change of support occurs onthe switching manifold (31). For small variation around theperiodic motion, the final position error for each step satisfies: δX − i = − CδY − i . (84)Using (83), the error in position at the end of the step, candeduced as function of L i . We have: δY − i ≈ − L i ˙ X ˙ Y ω ( − ˙ Y C + ˙ X )( ˙ X + ˙ Y ) . (85)The others errors can also be deduced: δX − i ≈ C L i ˙ X ˙ Y ω ( − ˙ Y C + ˙ X )( ˙ X + ˙ Y ) (86) δ ˙ X − i ≈ C L i ˙ X ˙ Y ˙ X ( − ˙ Y C + ˙ X )( ˙ X + ˙ Y ) + δ ˙ X + i (87) δ ˙ Y − i ≈ L i ˙ X ˙ Y ˙ Y ( − ˙ Y C + ˙ X )( ˙ X + ˙ Y ) − δ ˙ Y + i (88)We can then calculate the change of synchronization valuestep by step, we obtain: L i +1 L i ≈ − ˙ X δ ˙ Y − i + ˙ Y δ ˙ X − i L i (89)Using (79) and (80) equation (89) becomes: L i +1 L i ≈ ˙ X δ ˙ Y + i + ˙ Y δ ˙ X + i L i +2( C ˙ Y ˙ X − ˙ X ˙ Y ) ˙ X ˙ Y ( − ˙ Y C + ˙ X )( ˙ X + ˙ Y ) (90)With (72) equation can be rewritten (90) as follows: L i +1 L i ≈ ( ˙ Y − ˙ X )( C ˙ Y + ˙ X )( ˙ X + ˙ Y )( − C ˙ Y + ˙ X ) (91) VII. A
CKNOWLEDGMENTS
We gratefully acknowledge partial support by ANR EquipexRobotex. The work of J. Grizzle was supported by NSFAwards NRI-1525006 and Inspire-1343720.R
EFERENCES[1] K. Akbari-Hamed and J.W. Grizzle. Reduced-order framework forexponential stabilization of periodic orbits on parameterized hybrid zerodynamics manifolds: Application to bipedal locomotion.
IFAC journalNonlinear Analysis: Hybrid Systems (NAHS) , 25:227–245, 2017.[2] A. D. Ames, K. Galloway, K. Sreenath, and J.W. Grizzle. Rapidlyexponentially stabilizing control Lyapunov functions and hybrid zerodynamics.
Automatic Control, IEEE Transactions on , 59(4):876–891,2014.[3] C. Byrnes and A. Isidori. Asymptotic stabilization of nonlinear minimumphase systems. 376:1122–37, 1991.[4] C. Chevallereau, J. W. Grizzle, and C.-L. Shih. Asymptotically stablewalking of a five-link underactuated 3D bipedal robot. 25(1):37–50,2009.[5] Christine Chevallereau, Gabriel Abba, Yannick Aoustin, Franck Plestan,E R Westervelt, Carlos Canudas-de wit, and J W Grizzle. RABBIT : ATestbed for advanced Control Theory.
IEEE Control Systems Magazine ,pages 57–79, 2003.[6] Christine Chevallereau and Yannick Aoustin. Self-stabilization of 3Dwalking via vertical oscillations of the hip. In
International Conferenceon Robotics and Automation , pages 5088–5093, 2015.[7] Christine Chevallereau, J. W. Grizzle, and Ching Long Shih. Asymp-totically stable walking of a five-link underactuated 3-D bipedal robot.
IEEE Transactions on Robotics , 25(1):37–50, 2009.[8] S. Collins, A. Ruina, R. Tedrake, and M. Wisse. Efficient bipedal robotsbased on passive-dynamic walkers.
Science , 307(1082–1085), 2005.[9] M. Doi, Y. Hasegawa, and T. Fukuda. 3D dynamic walking based onthe inverted pendulum model with two degree of underactuation. ,pages 4166–4171, 2005.[10] Saunders J. et al. The major determinants in normal and pathologicalgait.
The Journal of Bone and Joint Surgery , 35(3):543–558, 1953.[11] B. Griffin and J.W. Grizzle. Walking gait optimization for accom-modation of unknown terrain height variations. In
American ControlConference , 2015.[12] Brent Griffin and Jessy Grizzle. Nonholonomic virtual constraints andgait optimization for robust walking control.
The International Journalof Robotics Research , Accepted 2017.[13] Jessy.W Grizzle and Christine Chevallereau.
Virtual Constraints and Hy-brid Zero Dynamics for Realizing Underactuated Bipedal Locomotion .https://arxiv.org/pdf/1706.01127.pdf, 2017.[14] D. G. E. Hobbelen and M. Wisse. Swing-leg retraction for limit cyclewalkers improves disturbance rejection.
IEEE Transactions on Robotics ,24(2):377–389, April 2008.[15] Y Hurmuzlu and D B Marghitu. Rigid body collisions of planarkinematic chains with multiple contact points.
The International Journalof Robotics Research , 13(1):82–92, 1994.[16] Jianjuen J.Hu.
Stable Locomotion Control of Bipedal Walking Robots:Synchronization with Neural Oscillators and Switching Control . PhDthesis, Massachusetts Institute of Technology, 2000.[17] B.G. Buss K. Akbari-Hamed and J.W. Grizzle. Exponentially stabiliz-ing continuous-time controllers for periodic orbits of hybrid systems:Application to bipedal locomotion with ground height variations.
TheInternational Journal of Robotics Research , 35(8):977–999, 2015.[18] Shuuji Kajita, Fumio Kanehiro, Kenji Kaneko, Kazuhito Yokoi, andHirukawa Hirohisa. The 3D Linear Inverted Pendulum Mode : Asimple modeling for a biped walking pattern generation.
IEEE/RSJInternational Conference on Intelligent Robots and System , 4:239–246,2001.[19] T. Koolen, T. de Boer, J. Rebula, a. Goswami, and J. Pratt. Capturability-based analysis and control of legged locomotion, Part 1: Theory andapplication to three simple gait models.
The International Journal ofRobotics Research , 31:1094–1113, 2012.[20] T. McGeer. Passive dynamic walking.
The International Journal ofRobotics Research, 9(2), 62–82. , 1990.[21] T. McGeer. Passive walking with knees. In
Proc. 1990 IEEE int. conf.robot. autom. (ICRA) , pages 1640–1645, Cincinnati, May 1990. [22] Benjamin Morris and Jessy W. Grizzle. Hybrid invariant manifolds insystems with impulse effects with application to periodic locomotion inbipedal robots.
IEEE Transactions on Automatic Control , 54(8):1751–1764, 2009.[23] N. Pateromichelakis, A. Mazel, M. A. Hache, R. Gelin T. Koumpogian-nis, B. Maisonnier, and A. Berthoz. Head-eyes system and gaze analysisof the humanoid robot romeo. In
Proc of the 2014 IEEE/RSJ Int. Conf.on Intelligent Robots and System , pages 1374–1379, Chicago, USA,2014.[24] Matthew J. Powell, Wen-Loong Ma, Eric R. Ambrose, and Aaron D.Ames. Mechanics-based design of underactuated robotic walking gaits:Initial experimental realization. In , pages 981–986, Cancun,Mexico, Nov 15-17 2016.[25] Hamed Razavi, Anthony M. Bloch, Christine Chevallereau, and J. W.Grizzle. Restricted Discrete Invariance and Self-Synchronization ForStable Walking of Bipedal Robots. In
American Control Conference ,pages 4818–4824, July 2015.[26] Hamed Razavi, Anthony M. Bloch, Christine Chevallereau, and Jessy W.Grizzle. Symmetry in legged locomotion: a new method for designingstable periodic gaits.
Autonomous Robots , 41(5):1119–1142, 2017.[27] Seyed Hamed Razavi.
Symmetry Method for Limit Cycle Walking ofLegged Robots . PhD thesis, The University of Michigan, 2016.[28] Ching-Long Shih, J. W. Grizzle, and Christine Chevallereau. Fromstable walking to steering of a 3d bipedal robot with passive point feet.
Robotica , 30(7):119–1130, 2012.[29] M.W. Spong and F. Bullo. Controlled symmetries and passive walking.
IEEE Transactions on Automatic Control , 50(7):1025–1031, 2005.[30] P. van Zutven and H. Nijmeijer. Foot Placement Indicator for Balanceof Planar Bipeds with Point Feet.
International Journal of AdvancedRobotic Systems , 10:1–11, 2013.[31] E.R. Westervelt, J.W. Grizzle, and D.E. Koditschek. Hybrid zerodynamics of planar biped walkers.
IEEE Transactions on AutomaticControl , 48(1):42–56, jan 2003.[32] Eric R Westervelt, Jessy W Grizzle, Christine Chevallereau, Jun HoChoi, and Benjamin Morris.
Feedback Control of Dynamic BipedalRobot Locomotion . CRC Press, 2007.[33] Chen. Z., Lakbakbi Elyaaqoubi N., and Abba G. Optimized 3D stablewalking of a bipedal robot with line-shaped massless feet and sagittalunderactuation.