Analysis of Foot Slippage Effects on an Actuated Spring-Mass Model of Dynamic Legged Locomotion

The classical model of spring-loaded inverted pendulum (SLIP) and its extensions have been widely accepted as a simple description of dynamic legged locomotion at various scales in humans, legged robots and animals. Similar to the majority of models in the literature, the SLIP model assumes ideal sticking contact of the foot. However, there are practical scenarios of low ground friction that causes foot slippage, which can have a significant influence on dynamic behaviour. In this work, an extension of the SLIP model with two masses and torque actuation is considered, which accounts for possible slippage under Coulomb's friction law. The hybrid dynamics of this model is formulated and numerical simulations under representative parameter values reveal several types of stable periodic solutions with stick-slip transitions. Remarkably, it is found that slippage due to low friction can sometimes increase average speed and improve energetic efficiency by significantly reducing the mechanical cost of transport.


Introduction
Legged locomotion is ubiquitous in the motion of living creatures throughout a wide range of scales, from small insects to rodents, mammals and humans [16]. One type of legged locomotion that is typical to robotics is quasistatic motion, where the body and limbs slowly move while continuously maintaining static balance with external loads [38,46]. On the other hand, natural walking and running typically involve dynamic motion, where the body is constantly falling as an unstable inverted pendulum until a free foot reaches contact [8]. Inspired by biological locomotion capabilities in nature, many legged robots have been built and developed (cf. [49]), and a large body of theoretical work on the analysis of dynamic legged locomotion (DLL) [28] and on the non-linear control of legged robots [29,63] have been conducted.
Models of robotic legged locomotion often include a chain of rigid links, as in the well-known examples of passive dynamic walking [33,34] and compass biped [24]. Nevertheless, in biological legged locomotion, elastic energy storage often plays an important role in dynamic behaviour. A classic and simple model of DLL that captures this effect is the SLIP (spring loaded inverted pendulum) model, which is shown in Figure 1. This model, composed of a point mass and a massless linear spring, has been widely accepted as a low-dimensional description of the leg-body motion in the sagittal plane for running animals, humans and even some bio-inspired legged robots [2,7,8]. The dynamics of the SLIP model is composed of two different phases of motion: the stance phase, in which the spring's endpoint (foot) is in contact with the ground and the flight phase, in which the mass flies above the ground while the spring is uncompressed, as illustrated in Figure  1(b). Thus, this motion is described by a hybrid dynamic system that is governed by two different sets of differential equations, depending on the contact state, combined with conditions of discrete transitions between the states. In spite of its complex non-linear structure, the dynamics of the SLIP model can be simplified and reduced due to conservation of energy and invariance with respect to the x coordinate of forward motion, thus enabling analysis of periodic solutions by studying a scalar Poincaré map that encapsulates the system's discrete-time dynamics [3,27].
Many extensions of the SLIP model have been considered in the literature, such as incorporating the leg's mass [42], or a large body with a moment of inertia [53]. Some works introduced damping as a mechanism for dissipating energy [4,48], while others studied the influence of a suspended load on the dynamic behaviour of the SLIP model [1]. Several works studied the incorporation of mechanical actuation into the SLIP model by adding torque at the revolute joint [4,56], or by using linear actuation along the leg, which effectively changes the unstressed length of the linear spring [43,55]. An experimental realization of the latter concept has been demonstrated by the bow leg mechanism of a robotic hopper [11,15]. Most of the works on actuated versions of the SLIP models studied control laws that are needed for the stabilization of desired cyclic motion and for achieving robustness under variations in the terrain's height; cf. [44,52].
A key limitation of the theoretical DLL models noted mentioned is that they assume that the foot is always in sticking contact with the ground and acts as a stationary pivot point. This requires ground friction to be large enough, an assumption that is not always realistic on natural terrains. It has already been observed that slippage can have a prominent influence on biological legged locomotion across all macroscopic scales. At the millimetre length scale, weaver ants have evolved several adhesive mechanisms for overcoming slippage on the lubricated surfaces of pitcher plants [9]. Stick insects confronted by a slippery surface modulate their motor outputs to produce normal walking gaits, despite a drastic change in the loads that these limbs experience [25,26]. In the centimetre scale, the elephant shrew is remarkable at constructing known escape routes in loose, slippery dirt [57]. Clark et al. [13] studied slippage in running guinea fowl and showed that falling on slippery surfaces is a strong function of both speed and limb posture at touchdown. In the scale of metres, the largest body of work on slipping concerns human biomechanics. The conditions that cause slipping [36], its consequences [35] and dynamics [59] have been studied, with important applications in preventative care for the elderly, shoe design, safety and more. Finally, motion measurements of horses' feet during galloping revealed the existence of a significant phase of foot slippage during ground contact [58].
In legged robotics literature, few works have incorporated friction considerations in designing gaits that avoid slippage [12,30] or "reflexive" control laws for overcoming undesired slippage [10], while some works studied detection and the estimation of slippage using sensory data [31,61]. Nevertheless, observations show that sometimes, slippage is an inevitable and non-negligible part of motion on low-friction surfaces, and that it significantly alters gaits and overall dynamic behaviour.
Very few works have incorporated friction limitation and possible stick-slip transitions of foot contact into theoretical models of DLL. Recent works by Zarrouk and Fearing [64,65] incorporate friction and slippage into a theoretical model of hexapod robot motion in a horizontal plane. In a vertical plane, the work of Berkemeier and Fearing [6] studied control strategies for an underactuated two-link monopedal hopper with stick-slip contact transitions. Recent works by Gamus and Or [18,19] studied periodic solutions (gaits) with stick-slip transitions in simple rigidbody models of passive and actuated walking, and a systematic classification of such gaits has been presented in [60]. The work of Seipel and Holmes [51] studied a torqueactuated version of the SLIP model, and accounted also for Coulomb friction limit at the foot contact. However, when the friction limit was reached, it was assumed in [51] that the contact immediately detaches without any phase of foot slippage. The work [45] analysed SLIP dynamics under different contact models of frictional forces, which depend on slippage velocity. The parameters of the contact laws in [45] were extremely sensitive and had to be hand-tuned by running numerical simulations. On the other hand, Coulomb's law of dry friction [14,41] is a broadly accepted classic model of contact force, which has a very clear physical meaning and where the friction coefficient can be easily estimated empirically.
The goal of this work is to study the dynamics of a generalized version of the SLIP model with torque actuation, foot mass and damping, under Coulomb's friction model and possible stick-slip transitions of the foot. We analyse the model under representative values of physical parameters and a particular choice of actuation control policy. We show that under physically realistic values of friction coefficient, different types of stable periodic solutions (gaits) with stick-slip transitions evolve. We then numerically investigate the influence of large changes in the friction coefficient on three quantitative characteristics of the gaits: stability (in terms of the discrete system's linearization eigenvalues), average forward speed, and energetic efficiency in terms of the specific mechanical cost of transport [20,32,47]. It is found that decreasing the friction causes degradation in stability. On the other hand, slippage induced by realistic values of friction can significantly increase both speed and energetic efficiency.
The structure of the paper is as follows. The next section reviews known results pertaining to the classic SLIP model and demonstrates the concept of the Poincaré map for analysis of hybrid dynamical systems and the stability of periodic solutions. Section 3 incorporates Coulomb's law of friction for sticking and slipping contact into the classic SLIP model. It is then shown that one must add non-zero foot mass in order to enable the existence of periodic solutions that include foot slippage. This addition implies that the system's motion must undergo impact at ground collisions and following on, impact equations are introduced. Section 4 presents the generalized model TDMF-SLIP -torqued damped (foot-) mass frictional SLIP, and formulates its equations of motion. Open loop control is used for the actuation torque during stance phase, while a simple PD control law is proposed for the flight phase, which drives the leg angle to its desired value before landing. Section 5 presents simulation results and a numeric investigation of the influence of friction and slippage on periodic solutions, followed by a discussion of the findings. Finally, the closing section lists some limitations of the work and proposes directions for future extension of the research.

Review of Classic SLIP Dynamics
In this section, the well-known SLIP model and its hybrid dynamics are briefly reviewed. The terminology of hybrid dynamics, Poincaré maps, periodic solutions and their stability is also briefly introduced by using the SLIP model as an example. The SLIP model, shown in Figure 1(a), consists of a point mass m connected to a massless springy leg with linear stiffness constant k and a free length of l 0 . The leg's orientation angle is denoted by θ and the spring's length is denoted by r. When in contact with the ground, the compressed spring applies a force of k(l 0 − r) > 0. When it is released at r = l 0 and ṙ > 0, the leg is detached from the ground and the mass begins ballistic flight under the influence of gravity acceleration g. It is assumed that during flight, the leg's angle can be set back to a prescribed value of θ = θ 0 before reaching the next touchdown event, where the uncompressed spring hits the ground. An illustration of motion snapshots and a typical trajectory of the mass are shown in Figure 1(b).
The hybrid dynamics of the SLIP is composed of flight phase and stance phase, which are formulated as follows.
During flight phase, the SLIP is governed by free ballistic motion: ÿ m = − g, ẍ m = 0, where (x m ,y m ) denote the position of the mass. The flight ends at the touchdown event, where y m = l 0 cos θ 0 , followed by transition to the stance phase. The dynamics of the stance phase is best described using polar coordinates (r, θ ) and is formulated as The stance phase ends at a lift-off event where r = l 0 .
Periodic solutions of the SLIP are analysed by using the notion of Poincaré maps (cf. [3,27]), summarized as follows. First, one has to choose a Poincaré section Σ, which is a defined by scalar constraint σ( x ) = 0 and where x is the state of the system, i.e., positions and velocities. The Poincaré map Π of the system then maps any initial state x that lies on section Σ to the state at which solution x (t) crosses Σ next. In classic analysis of the SLIP model (cf. [17]), the Poincaré section Σ is simply chosen as the state during flight phase where the mass reaches its apex point, i.e., ẏ m = 0. Thus, the Poincaré map Π induces a discrete-time dynamic system of states at apex times as x i+1 = Π( x i ), where i is the discrete-time index of the i th apex event. A periodic solution of the system corresponds to a fixed point x * of the Poincaré map, which satisfies Π ( x * ) = x * . Local stability of the periodic solutions is determined by eigenvalues of the matrix D Π ( x * ), which is the linearization of the Poincaré map Π( x ) about the fixed point x * . The periodic solution is locally asymptotically stable if and only if all eigenvalues λ i of D Π ( x * ) lie within the unit disc, i.e., they must all satisfy | λ i | < 1.
It has already been observed in the literature that SLIP dynamics are independent of horizontal position x m and that total mechanical energy E = 1 2 m( ẋ m 2 + ẏ m 2 ) + 1 2 k (r − l 0 ) 2 + mg y m is conserved during the entire motion. Therefore, the system's state x at an apex time can be fully described by the apex height ȳ for any given initial mechanical energy E 0 . The Poincaré map can thus be reduced to a scalar map of apex heights ȳ i+1 = Π( ȳ i ) and the stability condition becomes an inequality on its derivative | Π ′ ( ȳ = y * ) | < 1, where y * are fixed points of map Π. The Poincaré map Π( ȳ ) of the SLIP cannot be formulated in closed form, though some works have used additional simplifying assumptions in order to derive approximate analytical expressions; cf. [21,50]. Nevertheless, since Π( ȳ ) is a scalar map, it can be computed via numerical integra-tion and from plotting its graph, one can easily find fixed points and determine their stability.
As a simulation example, we chose the physical parameter values of k = 20,000N / m, l 0 = 1m, m = 80Kg and E 0 = 6000J, which are in the order of magnitude of a running human [17]. Graphs of the Poincaré map ȳ i+1 = Π( ȳ i ) for several values of touchdown angle θ 0 are plotted in Figure 2(a). The dash-dotted line of slope 45° is ȳ i+1 = ȳ i and the intersection points of this line with the graph curves are exactly fixed points, satisfying Π(y * ) = y * , which is associated with periodic solutions of the SLIP. It can be seen in the plot that each Poincaré map curve may have multiple fixed points. Moreover, the stability of each fixed point can be visually determined by the slope of the Poincaré map curve at the intersection points, which must lie within the range of ±45°t o achieve stability. The stable fixed points in Figure 2(a) are marked with '▲', while the unstable points are marked with '■'. It can be seen that there is a small range of touchdown angles θ 0 around 40°− 44° for which there exist stable fixed points. The points of maximum height on each Poincaré map curve correspond to vertical jumps ẋ m (t i+1 ) = 0, at which the entire mechanical energy is converted to potential energy. Continuation of the plots beyond this point (dotted curves) corresponds to a reversal of the horizontal velocity ẋ m , followed by penetration into the ground by crossing y=0 after the next touchdown; hence, they are excluded from the plot 1 . Figure 2(b) plots the curve of fixed points y * as a function of θ 0 , where stable fixed points are denoted by a solid curve while unstable points are denoted by a dashed curve.

Adding Coulomb Friction, Slippage and Impacts to the SLIP Model
In this section, the effect of Coulomb friction bounds on the dynamics of the SLIP model is considered. It is shown that non-zero foot mass must be added to the model in order to obtain slippage for a finite non-zero duration. Next, inelastic impact of the foot mass with the ground is considered and a condition for avoiding post-impact foot slippage is formulated.

Coulomb friction and slippage
Coulomb's friction law [14,41] imposes physical limitations on the contact force that the ground applies at the spring's endpoint. Denoting the tangential and normal components of the contact force as f = ( f x , f y ), the condition under which sticking contact with the ground can be maintained is given by where μ is called the coefficient of friction, which is assumed to be a property of the two surfaces in contact.
Since in the SLIP model, the spring is a massless link, it can apply forces only along its axis; the ground force is given by f = ( − f sin θ , f cos θ ) and the no-slip constraint (2) directly leads to an inequality in the orientation angle of the SLIP: Geometrically, (3) implies that the spring's orientation is restricted to lie within a friction cone with a half angle of tan −1 (μ) about the vertical direction, as shown in Figure 3(ab). Along periodic solutions of the SLIP, the no-slip condition reduces to the bound | θ 0 | ≤ tan −1 (μ), since θ (t) varies monotonically within the symmetric interval − θ 0 , θ 0 .
Next, we consider the behaviour of the SLIP model in the case where the inequality constraint (3) is violated, either due to low friction μ or due to a large orientation angle | θ | . In this case, the spring's endpoint begins to accelerate 1 There also exist lower bounds on apex heights ȳ i for which the Poincaré map gives non-physical results; see Figure 2 in [17] for details. and slips horizontally in velocity v x ≠ 0. The tangential contact force is given by f x = − sgn(v x )μ f y , that is, it attains its maximal allowed value and opposes the direction of slippage 2 . However, since the spring is massless, this model dictates that when the foot mass is free to slip horizontally, it will move with infinite acceleration, while the upper mass does not change its position and velocity. That is, the spring will extend instantaneously to its free length l 0 , leading to immediate contact detachment. Two different scenarios of exceeding the bound (3) along a typical trajectory of the SLIP are shown in Figures 3(a-b). In Figure 3(a), the SLIP lands with the leg pointing forward and its angle satisfying θ > tan −1 (μ), while the velocity of the mass points downward ẏ m < 0. In this case, the lower endpoint of the spring immediately slips forward until reaching the free length l 0 (dashed). Then, the contact detaches and the mass simply falls down on the ground, precluding the existence of a periodic solution. In the second scenario, shown in Figure  3(b), the friction cone bound (3) is reached at a negative orientation θ = − tan −1 (μ) during stance, while the spring is already in its decompression stage and the velocity of the mass points upward ẏ m > 0. In this case, the foot immediately slips backward until the spring reaches its free length l 0 (dashed). Then, the contact detaches and the mass flies up. That is, this scenario simply leads to an "early take-off", as assumed in the work of Seipel and Holmes [51], and the existence of a periodic solution is still possible. Note that this early take-off also leads to instantaneous loss of the remaining potential energy stored in the spring.
The obvious conclusion is that the SLIP model with a massless spring results in instantaneous slippage and cannot capture any finite non-zero duration of slippage during legged locomotion, which has been observed experimentally in various scales and scenarios, as discussed in section 1. This motivates generalization of the SLIP model by adding a small point mass at the foot, i.e., the lower end of the spring. This addition necessarily leads to the occurrence of impact upon collision with the ground, as discussed next.

Impact of SLIP with a small foot mass
We now consider a SLIP model with a small non-zero mass m f added at the spring's lower endpoint, the position of which is denoted as (x,y). When the mass collides with the ground, an impact event occurs. This causes an instantaneous change in the velocity v of m f according to the Geometrically, this means that the direction of pre-collision velocity v − must lie within the downward-pointing friction cone with angle 2tan −1 (μ), as shown in Figure 3(c). Note that this condition is fundamentally different from the kinematic condition in (3). When condition (5) is violated, the impulse cannot impose sticking and the post-collision velocity satisfies ẏ + = 0, ẋ + ≠ 0, while the contact impulse satisfies F x = − sgn( ẋ − )μF y . Substituting the normal impulse 2 For the sake of simplicity, we do not distinguish between static and kinetic coefficients of friction throughout this work. Figure 3. (a) The SLIP exceeding the friction cone limit at touchdown. The foot instantaneously slips forward until the massless spring extends to its free length. Then, the mass falls down on the ground. (b) The SLIP reaching the edge of the friction cone at take-off. The foot instantaneously slips backward until the massless spring extends to its free length. Then, the mass flies up and contact is detached. (c) The SLIP with small foot mass approaching touchdown. In order to avoid slippage at the impact, the direction of pre-collision velocity v must lie within the downward-pointing friction cone.
When the foot mass is very small, i.e., m f → 0, the no-slip condition during stance approaches the geometric condition (3) ; see eq. (10) in the next section. However, the condition for sticking after impact in (5) is independent of m f and thus, it adds a fundamentally different constraint on SLIP trajectories that maintain sticking contact directly after impacts. As an example, we consider periodic trajectories of the SLIP model with the same physical parameters as in the simulation example of section 2, and compute the minimal value of friction coefficient μ, which is required in order to maintain sticking contact along a periodic solution as a function of the touchdown angle θ 0 .
The results are plotted in Figure 4(a), where the dashdotted line is tan −1 ( θ 0 ), which is the minimal value of μ for maintaining sticking contact during stance. The solid (stable) and dashed (unstable) lines are the minimal values of μ for maintaining sticking impact under an infinitesimally small foot mass m f according to condition (5), along the periodic solutions given in Figure 2(b). One can conclude from the plot that it is very likely that the SLIP model will suffer from foot slippage for realistic values of μ, particularly for stable periodic solutions. Moreover, it can be seen that the condition for sticking impact is much more restrictive than the condition for sticking contact during stance. For instance, the stable solution for θ 0 = 40°u nder an infinitesimally small foot mass requires a friction coefficient of μ > tan −1 (40°) = 0.84 for maintaining sticking contact during continuous-time motion, whereas an unrealistically large value of μ > 4.5 is required for obtaining sticking at impact. The findings of this section motivate the incorporation of friction and slippage into a generalized SLIP model with two masses, as presented next.

The Actuated Two-mass SLIP Model -TDMF-SLIP
In this section, we present a generalized SLIP model where a small foot mass m f > 0 is added at the lower endpoint of the spring. The addition of the mass results in two fundamental differences from the original SLIP with massless spring, as follows. First, the impact dissipates energy and the mechanical system is no longer conservative as the original SLIP model. Second, the assumption that during the flight phase the leg angle θ can be kinematically prescribed at no cost is no longer plausible, since the leg is not massless. Therefore, it is necessary to incorporate some mechanical actuation into the model. Inspired by the work of [4,56], we have chosen to add actuation of a hip torque τ(t) at the joint connecting the spring to the upper mass. The roles of this torque are to inject energy in order to compensate for impact losses, and to control the θ angle during flight. Another addition (which is not essential) is a viscous damper c connected in parallel with the spring in order to attenuate its oscillations during flight. Inspired by the names M-SLIP [42] and TD-SLIP [4], our model is named TDMF-SLIP -torqued damped (foot-) mass frictional SLIP. The model is shown in Figure 4 Using standard Euler-Lagrange formulation, the dynamics of the system is given in matrix form as The M ( q ) matrix is referred to as the system's matrix of inertia, B ( q , q ) contains velocity-dependent inertial terms, G ( q ) contains potential forces and F q contains nonconservative generalized forces. Expressions for all of these terms are given as: = .
x y q f f In order to avoid further complication of the model, it is assumed here that the upper mass m does not rotate. This assumption is quite reasonable for multi-limbed animals such as insects, dogs, horses, etc., where the influence of limbs' actuation torques on the rotation of a massive body is negligible. This implies that τ is effectively an external torque rather than an internal one and thus, the total angular momentum of the two-mass system is not conserved during its application (the same assumption is implicitly used in previous works on hip-actuated versions of the SLIP; cf. [51,56]). The contact forces f x , f y in (9) depend on the contact state, as follows. In case of sticking contact, the foot mass is stationary y = ẏ = ẋ = 0. Substituting the constraints ẍ = ÿ = 0 into (8), (9) and eliminating the accelerations r , θ , the contact forces are obtained as: hand, as f y vanishes, tangential force f x is generically nonzero and thus, it necessarily exceeds inequality (2). Therefore, it is concluded that for any m f > 0 and a finite coefficient of friction μ, onset of slippage must occur before the contact is detached.
In the state of contact slippage ẋ , ẍ ≠ 0, ÿ = 0, the normal contact force f y is still given by (10), while the tangential force satisfies f x = − μsgn( ẋ ) f y . Three different transitions can occur from this contact state: when the slippage velocity vanishes ẋ = 0, the contact either reverses its slippage direction or switches back to sticking, depending on the consistency of contact forces 3 . Finally, when the normal contact force f y in (10) vanishes, the contact detaches and the system switches to the flight state, as discussed next.
In the state of flight, the equations of motion (8)  velocities ( ṙ , θ ), which represent the relative motion between the two masses, undergo a discontinuous jump and their post-impact values can be obtained from ẋ + , ẏ + , and differentiation of the kinematic relations (7). The dynamics of this system thus constitutes a hybrid system [5,23] with smooth transitions between three contact states during continuous-time dynamics (stick, slip, flight), and one non-smooth discrete-time transition due to impact.
Next, we specify a control law for the torque actuation τ in this model. In the stance phase during ground contact (either sticking or slipping), we adopt the concept of halfstance actuation introduced in [55] and the torque is given as That is, a constant torque τ 0 is applied only until reaching the upright position. The role of this simple control is to inject energy during stance. During flight, the role of the control is to enforce change in the angle θ towards its desired value of θ 0 at the next impact via a simple PD feedback law. The chosen control law during flight is given by: The control law in (12) is divided into two parts, before and after the centre of mass reaches its apex point ( ẏ c = 0). In both parts, the control law has a form that is similar to a proportional-derivative (PD) feedback, with control gains of k p and k d , and different choices for desired set point values for the leg angle θ . In the first part, the desired angle is η θ LO , where η > 1 and θ LO is the orientation angle at liftoff. The reason for this choice is that under zero hip torque, the leg angle θ tends to keep rotating monotonically after lift-off ( θ < 0) due to the conservation of total angular momentum; thus, the positive torque is necessary for preventing this rotation. Choosing η > 1 turns out to be necessary for avoiding additional collisions with the ground directly after lift-off. After the centre-of-mass reaches its apex, the torque actuation simply moves the leg angle to its desired value of θ 0 , preparing for the next touchdown. Under this control law, the two-mass model dynamics can have solution trajectories that are quite similar to those of the original SLIP model, except for the addition of impacts and the possibility of slippage, as shown next.

Numerical Simulation Results and Analysis
In  Table 1. While the choice of values of α , β , γ , δ ,ζ and θ 0 has some physical and biological relevance, the parameters k p ,k d ,η -related to the control law -were chosen by trial and error in order to achieve a reasonably fast convergence of the leg's motion during flight to the desired trajectory, which is similar to that of the massless leg model.  Plots of the normalized heights of the foot y(t) and centreof-mass y c (t) along a full period from apex to apex are shown in Figure 5(a), where t is the normalized time. The two vertical dashed lines enclose the time duration of the stance phase. Figure 5(b) plots the normalized leg length r(t) and the leg angle θ (t) in degrees along a full period from apex to apex. One can see that r(t) converges to its free length of 1 (normalized by l 0 ) under the influence of damping and that it is compressed to values smaller than 1 during stance phase. Due to the action of the controlled hip torque τ in (12), the angle θ (t) converges to the desired value θ 0 after apex time and to the constant value of η θ LO after lift-off. Figure 5(c) plots the foot's tangential velocity ẋ (t). One can see that the impact results in immediate sticking ẋ + = 0, which persists for almost the entire duration of the stance phase. It can also be seen that ẋ (t) changes continuously but not smoothly at the beginning and end of the period, due to the discontinuous change in the control law (12) at apex time. Figure 5(d) plots the ground force ratio f x / f y during the stance phase.
One can see that when the normal force f y tends to zero directly before lift-off, this ratio reaches its lower bound of − μ and a very short period of positive slippage evolves, followed by immediate lift-off. The ' × ' at the left end of the plot denotes the ratio of impulses F x / F y at the impact.
Since this ratio is above − μ, condition (4) is satisfied and the impact results in sticking, ẋ + = 0. Thus, this periodic solution can be classified as type ' 0 + '. The vertical dashdotted line denotes the mid-stance time at which θ = 0, where the switch in the control law in (11) imposes a discontinuous jump in the ground reaction forces.    When the friction coefficient μ is decreased, condition (4) is no longer satisfied and the impact results in positive slip.
As an example, Figure 6 shows the simulation results of a stable periodic solution under friction coefficient μ = 1.2. Figure 6(a) plots the tangential foot velocity ẋ (t). The plot indicates that ẋ does not vanish due to impact; hence, the stance phase begins with a short episode of positive slippage followed by sticking and another brief slippage right before lift-off. Thus, this periodic solution can be classified as type ' +0 + '. The ratio of tangential to normal contact forces during stance phase f x / f y is shown in Figure   6(b), where the dashed line for which f x = − μ f y denotes slip, while the solid line denote sticking. One can observe that the transition from slip to stick involves discontinuity of the contact forces, while the converse transition (directly before liftoff) does not.
When μ is further decreased, another interesting type of stable periodic solution evolves, which is classified as ' + 0 − 0 + '. As an example, Figure  consists of three phases of slippage separated by two phases of sticking. Upon further decrease of μ, the two episodes of sticking gradually shrink until they disappear. For μ = 0.43, the tangential velocity ẋ (t) during the stance phase is plotted in Figure 8(a). This periodic solution is classified as type ' +0 − ', where the latter phase of sticking disappears and liftoff occurs after negative slip. For μ = 0.2, the tangential velocity ẋ (t) during the stance phase is plotted in Figure 8(b). This periodic solution is classified as type ' + − ', where the first phase of sticking also disappears and the foot slips during the entire stance phase. A reversal of the slippage direction occurs during stance, which implies a discontinuous jump in the contact forces and a nonsmooth change of velocities. When μ is decreased further, no qualitative change in the periodic solution is observed; however, numerical simulations for μ < 0.1 become extremely sensitive and no periodic solutions can be numerically obtained.

Influence of friction coefficient on performance
We now analyse the influence of changes in the coefficient of friction μ on three key characteristics of periodic solutions: stability, average speed and energetic efficiency. Stability is determined by the magnitude of the eigenvalues λ i of the linearization of the Poincaré map Π, where the condition λ max = max{ | λ i | } < 1 has to be satisfied for stability.
The linearization matrix D Π( x * ) was obtained by numerically computing the Poincaré map under perturbations of the state x about x * , in small variations of each state variable. Figure 9(a) plots the maximal eigenvalue magnitude λ max as a function of μ. The dashed vertical lines denote critical values of μ where the periodic solution undergoes a transition in its type, and the different types are shown as labels near the curve. The range μ ∈ 2,7 is cut out from the plot in order to improve the visibility of lower values of μ, since the periodic solution type remains as ' +0 + ' within this range, while the value of λ max changes monotonically. It can be seen that the stability of the periodic solution degrades as μ is decreased and slippage evolves. For values of μ below 0.3, computation of D Π ( x * ) became numerically sensitive, since perturbations about x * led to solutions with a second impact immediately after lift-off. This sensitivity can be mitigated by better tuning of the control parameters η,k p ,k d of the actuation torque τ in (12) for the flight phase. This insignificant tuning is not discussed in this work. The average speed of a periodic solution is defined as V = Δx / T , where Δx is the step length and T is the step's period time. Figure 9(b) plots V as a function of μ, where the range μ ∈ 2,7 is cut out from the plot. Remarkably, it can be seen that decreasing μ results in an increase in the average speed V , due to the evolution of forward slippage, which is not stopped by the impact. Maximal speed for the given parameter values is obtained at μ ≈ 0.5, which provides roughly a 12% increase compared to the speed at no-slip solutions for large friction. When μ is decreased below 0.5 the speed drops again, since the periodic solutions of type ' +0 − ' involve negative slip prior to lift-off, which cancels the effect of the positive slip after impact. When μ is decreased below 0.2, where solution type ' + − ' evolves, the negative slip becomes dominant and the speed V drops sharply.
The energetic efficiency of legged locomotion is evaluated based on the mechanical cost of transport [20,32,47], which is the mechanical work per distance, defined as W is the mechanical work expended by the hip torque actuation during a full period and Mg = (m + m f )g is the total weight. This is analogous to measuring "gallons per mile" in cars, where lower values of c mt correspond to greater energetic efficiency. Figure 9(c) plots the value of the mechanical cost of transport c mt of periodic solutions as a function of μ, where the range μ ∈ 2,7 is cut out from the plot. Remarkably, it can be seen that decreasing the friction coefficient reduces c mt substantially for all values of μ, in spite of the added energy dissipation due to slippage. The reasons for this increase are twofold: first, the onset of slipping impact, which does not fully stop the forward velocity results in increasing the step length Δx. Second, energy losses due to slipping impact are significantly smaller than those of sticking impact, at which the tangential velocity is always stopped. Optimal energetic efficiency is obtained at μ ≈ 0.5, which provides roughly a 25% decrease in c mt , compared to the value at no-slip solutions for large friction. When μ is decreased below 0.5, the duration of slip phases becomes large and the resulting dissipation causes an increase in energy expenditure W . Moreover, the onset of negative slip decreases the step length Δx and hence, the overall effect is that of increasing c mt .

Discussion
We now discuss the results of the numerical simulations presented above. First, the observation that periodic solutions with slippage due to low friction are less stable (larger λ max ) compared to those with large friction is in good agreement with the observed results in [18,19] for models of passive dynamic walking, except for the fact that in these works, stability was lost ( λ max > 1) below a critical value of μ. Note that lower linearization eigenvalues only indicate faster convergence to the periodic solution under small perturbations. Perhaps a more useful measure of stability will be in the region of attraction, which was not studied here, due to its high dimensional complexity in systems with many degrees of freedom. An important observation from the numerical simulations is that some of the numerically computed linearization eigenvalues of D Π ( x * ) were very close to zero, a fact that indicates that the effective dimension of the Poincaré map can be reduced. While the dimension of the state vector x is 8 (four degrees of freedom, coordinates x m ,y m ,x,y or x c ,y c ,r, θ and their velocities), the Poincaré section at apex is defined by the constraint ẏ c = 0 and hence, its dimension is 7. However, in all types of periodic solutions under different values of μ, four out of the seven numerically computed eigenvalues of D Π ( x * ) were very small, in the order of 10 −3 − 10 −5 . The reason for this is the presence of hidden constraints on the hybrid dynamics, which are more evident if one chooses different Poincaré sections along the periodic solutions.
According to the explanations in [62], the rank of the Poincaré map is determined by the Poincaré section with the smallest dimensions. If one chooses a Poincaré section at post-impact states, where y = 0 and ẏ > 0, the assumption of inelastic impact imposes an additional constraint on the state space, ẏ + = 0. Moreover, under sticking impact, which occurs in periodic solutions of type ' 0 + ', the constraint ẋ + = 0 is also imposed. That is, the Poincaré section at postimpact states has an effective dimension of 6 for slipping impact and 5 for sticking impact. On the other hand, if one chooses a Poincaré section at pre-impact states, where y = 0 and ẏ < 0, other constraints are imposed by the chosen control law at flight phase (12), which forces convergence of the foot angle θ (t) towards its desired value of θ 0 prior to touchdown, as seen in Figure 5(b). Moreover, the leg's damping attenuates the two-mass oscillations and the distance r(t) also converges towards its free length l 0 (normalized value of 1) as also seen in Figure 5(b). This implies that the state vector x at the pre-impact Poincaré section is very close to satisfy the constraints θ − = θ 0 ,r − = l 0 and θ ̣ − = ṙ − = 0. Therefore, according to the theory in [62], the effective rank of the Poincaré map drops from 7 to 3, as corroborated by the numerical results of the present study.
Finally, another interesting observation is that for some intermediate ranges of μ, the average forward speed V increases due to slippage, which contrasts the results obtained in [19] for the actuated compass biped walking under slip-stick transitions. A possible explanation for this difference is the specific choice made for the control law in [19], which used a reference trajectory with small step length that reduces walking speed. Nevertheless, the increase in energetic efficiency due to the presence of slipping impact in the periodic solution is in good qualitative agreement with the results in [19].

Conclusion
In this paper, we analysed a generalized version of the SLIP model that incorporates the effects of Coulomb friction and foot slippage by adding a foot mass, linear damping and controlled torque actuation. The influence of the friction coefficient on different types of periodic solutions with different contact transitions was analysed. It was found that periodic solutions under low friction and foot slippage are less stable in terms of linearization eigenvalue magnitude. On the other hand, for some intermediate ranges of μ, average speed and energetic efficiency were significantly increased in the presence of slippage, mainly due to slipping impact which does not completely stop the foot's tangential velocity and thus results in smaller losses of energy and forward speed, compared to a no-slip periodic solution. For the chosen parameter values, an optimal value of friction coefficient that maximizes the average speed and energetic efficiency is μ ≈ 0.5, which is a practically achievable value.
We now briefly discuss the limitations of this work and propose possible directions for future extension of the research. First, it is obvious that the proposed complex model of TDMF-SLIP lacks much of the elegance and simplicity of the original SLIP model, mainly due to the increase in dimensionality and the necessity of arbitrarily specifying a particular control law for the hip torque during flight. Nevertheless, the reduction in dimensionality of the Poincaré map due to four nearly-zero eigenvalues suggests that this model can be simplified into a lower-dimensional model. One may consider a reduced model, in which the details of two masses and torque actuation are taken into account only during the stance phase, while during flight phase, the two masses are approximated as a single rigid body, whose angle θ is prescribed kinematically upon touchdown. This model has a three-dimensional Poincaré section at the pre-impact state, which can be further reduced to two dimensions by ignoring the x coordinate, due to invariance of the dynamics with respect to it. This proposed reduction, which is currently under our investigation, results in a much simpler model that is amenable to more thorough analysis of the influence of the system's parameters on dynamic behaviour and performance.
Second, it is planned to systematically fit the parameters of the TMDF-SLIP model in order to match measured results from the legged motion of animals where significant foot slippage is observed, in order to test the biological relevance of the model.
Third, it is important to incorporate different methods of actuation into the model and test their influence on the modification of gaits in order to minimize foot slippage, or for harnessing slippage to improve performance. In addition, incorporation of control strategies developed for various generalizations of the SLIP model into the TMDF-SLIP model in order to test their interaction with possible foot slippage on varying terrains is a challenging open problem. In particular, it is worth revisiting the concept of swing leg retraction (SLR) [52]. Inspired by legged animal locomotion, the SLR concept involves retraction of the leg angle θ during flight according to a linear time profile, θ (t) = θ a − ω 0 (t − t a ), where t a is the time at apex and θ a ,ω 0 are given constants. The SLR method is known to significantly enhance the stability of the periodic solutions of the SLIP model and to also add robustness with respect to uncertainties in ground height [54]. Since this control concept affects the foot's velocity upon touchdown, it will be interesting to analyse its interaction with friction constraints and possible foot slippage. Finally, our model can be extended by considering more detailed models for contact, friction and slippage; cf. [22,45].