An Improved Lyapunov Stability Test of Equilibria Under Frictional Unilateral Contact by Sums of Squares Programming

Reliable quasi-static object manipulation and robotic locomotion require the verification of the stability of equilibria under unilateral contacts and friction. In a recent paper, Posa et al. (2016) demonstrated that sums-of-squares (SOS) programming can be used to verify the Lyapunov stability of planar systems via Lyapunov’s direct method if impacts are inelastic. However, this method appears to be too conservative to verify the stability of some remarkably simple examples. In this article, an extension of Lyapunov’s direct method is proposed, which makes use of a piecewise smooth Lyapunov function and allows a temporary increase of the Lyapunov function along a motion trajectory. In addition, a modified SOS formulation enables the investigation of planar systems with partially elastic contacts. The proposed method remains compatible with SOS programming techniques. The improved stability test is successfully applied to a point mass on a slope and to a rigid body with two contact points. For the latter, several mechanisms of instability have been demonstrated experimentally, but the exact conditions of Lyapunov stability are unknown.


I. INTRODUCTION
M ANY tasks in robotics involve unilateral contact with friction between hard objects.Small perturbations of an equilibrium state may induce stick-slip transitions, contact separation and impacts [1], [2].Within the framework of rigid-body dynamics, nonsmoothness and discontinuity of the response prevent one from using classical tools of stability analysis like linearization.At the same time, dry friction induces continuous sets of equilibrium states [3], which means that stable individual points within a set are usually not associated with local minima of potential energy.Moreover, the emerging hybrid dynamics has special features, such as Zeno The author is with the Department of Mechanics Materials and Structures, Budapest University of Technology and Economics, H-1111 Budapest, Hungary (e-mail: varkonyi.peter@epk.bme.hu).
Digital Object Identifier 10.1109/TAC.2023.3333738sequences of impact events [4], [5], and issues of nonexistence and nonuniqueness [6].For example, in the simple case of a planar rigid body with two contact points, conditions of stability are known for ideally inelastic impacts [7], [8], but only highly conservative conditions have been found in the case of partially elastic impacts [9], [10].Stability conditions of most systems with more than two contact points, extended areas of contact, as well as of 3-D models remain unknown.
There are numerous extensions of Lyapunov's direct method [11] tailored to deal with nonsmooth and hybrid dynamical systems [12], [13].Many of them involve nonsmooth [14], [15], [16] or multiple Lyapunov functions [17], [18], [19], [20].In this article, we test limitations of a recently proposed computational method [21] based on a single smooth Lyapunov function, which is tailored to deal with the stability analysis of planar perfectly inelastic contact systems.A modified version of that method using a nonsmooth Lyapunov function is proposed, and its wider applicability is demonstrated.
Sums-of-squares (SOS) optimization offers an efficient way to construct certificates of Lyapunov stability and invariance of sets in the case of hybrid systems [21], [22], [23].The algorithmic approach to the stability analysis of contact systems by Posa et al. [21] uses SOS polynomials and semidefinite programming to develop sufficient conditions of local Lyapunov stability.In addition, their method delivers estimations of the basin of attraction of a stable point, and it can also be used to verify the positive invariance of a set of states.This approach is quite general as it is formally applicable to any system for which the equations of motion and the constraints are expressed as polynomial equations or inequalities.At the same time, success is not guaranteed as the proposed conditions are conservative.Notably, most of the equilibria successfully tested by Posa et al. [21] correspond to a local minimum of potential energy.It remains an open question if more challenging problems can be handled by this method.
In this article, it is found that the method proposed by Posa et al. [21] fails to verify the stability of equilibria for two simple model problems, including a point mass on a slope and a planar rigid body with two contact points on a slope.This finding motivates the development of two improvements using ideas borrowed from the Lyapunov theory of nonsmooth and multiple Lyapunov functions.
First, the Lyapunov function is allowed to increase temporarily along trajectories, provided that it has a decreasing overall trend.Nonmonotonic Lyapunov candidates have been explored extensively in the context of stability analysis using multiple Lyapunov functions [17].
Second, the proposed extension makes use of a piecewise smooth Lyapunov function, constructed as the upper envelope of several polynomials.Our approach is somewhat unusual because the switching manifolds of the Lyapunov function need not correspond to switching manifolds of the hybrid dynamics in contrast to the standard approach using nonsmooth Lyapunov functions.Accordingly, each function is required to decrease in every subdomain of the hybrid dynamics.In turn, the individual smooth components of the Lyapunov function need not be positive definite as required by the standard theory involving multiple Lyapunov functions [17].
Finally, the present work makes use of a standard approximation of piecewise continuous systems: the so-called zero-order dynamics (ZOD).In each contact mode of motion, the equations of motion are approximated by their lowest order (constant) terms, which gives accurate local description of the system in a small neighborhood of an equilibrium state.This approach is highly similar to linearization techniques of smooth dynamical systems.
The rest of this article is organized as follows.In Section II, general notation and problem statement are introduced, and key results of [21] are reviewed.The existing stability theory is tested on two new model problems in Section III.After the extension of the existing stability theory in Section IV, the two test problems are revisited in Section V. Finally, Section VI concludes this article and points out related open problems.

A. Notation and Problem Statement
Column vectors and matrices are denoted by boldface lowercase and uppercase letters, respectively.Lower indices are used to denote elements of vectors.For example, v i is the ith element of vector v. Furthermore, J γ (q) means the Jacobian of γ(q).The left and right limits of piecewise continuous functions are denoted by upper indices ±, i.e., v − (t) and v + (t), respectively.Transpose is denoted by an upper index T .Derivation with respect to time is denoted by a dot, i.e., q.
A rigid body or a rigid multibody system is considered in 2-D with c ≥ 1 unilateral contact points and n degrees of freedom.The state x of the system is given by the state vector x T = [q T , v T ] composed of generalized coordinates q ∈ R n and generalized velocities v = q.Hence, we have a 2n-dimensional state space x ∈ S ≡ R 2n .We will consider an equilibrium state given by x 0 = (q 0 , 0).Without loss of generality, q 0 = 0 is assumed.
The admissible set A is defined as a subset of the state space given by the vector-valued inequality constraint where γ ∈ R c and γ i (q) (i ∈ {1, 2, . .., c}) are gap functions associated with the unilateral contacts.As we perform local analysis of the equilibrium q 0 = 0, we are only interested in constraints with γ i (0) = 0.The unilateral contacts may give rise to nonnegative normal contact forces λ Ni subject to the linear complementarity condition ( Furthermore, if γ i = 0, then there is also an analogous condition at the velocity level The normal velocity γi at contact i is determined by using the chain rule as γi = J γ i (q)v(t).
In the presence of friction, it is convenient to introduce the vector-valued tangential displacement function σ(q), since the behaviors of contacts depend on the relative tangential velocities given by σ = J σ i (q)v(t).The sign of σi may be positive, zero, or negative, which correspond to slip in the positive direction, stick, and negative slip at contact i. Coulomb friction force is subject to the constraint where μ i ≥ 0 are friction coefficients associated with the contacts and sign(•) denotes the set-valued sign function defined as sign(s) = 1 for s > 0, sign(s) = −1 for s < 0, and sign(s) = [−1, 1] for s = 0. Throughout this article, static and kinetic friction coefficients are assumed to be equal.During episodes of continuous motion, the dynamics of the system is governed by the standard manipulator equations [24] q = v (5) which can be derived from Lagrange's equations.Here, λ N and λ T are vectors containing normal and tangential forces at all the contacts.The concepts introduced so far are illustrated by the following two examples.
Example 1 (Point mass on a slope): Consider a point mass resting on a slope of angle α subject to gravity and a coordinate system yz with its y-axis aligned with the slope (Fig. 1).Units of time and mass are chosen such that the mass of the object and the gravitational constant are equal to 1. Let (y, z) denote position coordinates of the point mass.Then, q = [y, z] T and v = q = [u, w] T .The parameters of the manipulator equation and displacement functions are Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Example 2 (Rigid biped on arbitrary terrain): Consider a planar rigid body with two sharp vertices resting on two straight surfaces, as shown in Fig. 2 .This example will be referred to in the following as "biped."We fix a global coordinate frame yz such that the origin coincides with the center of mass in the initial equilibrium configuration and the y-axis is parallel to the line through the initial positions of the contact points.The object is subject to constant external forces, which are lumped into a resultant force of size F and directional angle α acting at the center of mass as well as a resultant torque T .Units of length, time, and mass are chosen such that the mass, the radius of gyration, and the gravitational constant are equal to 1.The remaining model parameters include the lengths h, l 1 , and l 2 , the angles φ 1 and φ 2 , and the friction coefficients μ 1 and μ 2 .We can use the global coordinates (y, z) of the center of mass and the rotation angle θ of the body as generalized coordinates, i.e., q = [y, z, θ] T .The generalized velocities are denoted by v = [u, w, ω] T .The parameters of the manipulator equation and the gap functions and tangential displacement functions are where the following notation has been used: For these two examples, c and H are state independent.

B. Hybrid Dynamics Approach
In order to find contact forces and acceleration simultaneously, several standard methods can be used.Numerical methods may treat the problem as a linear complementarity problem, which addresses all the cases of the sign function in a unified framework.In contrast, a detailed analysis of the emerging motion is usually done within the framework of hybrid dynamics (especially for systems with moderate number of contact points), as, for example, in [7] and [8].The hybrid dynamics approach is used here, i.e., it is assumed that the system undergoes episodes of continuous motion in one of its contact modes.Such episodes are interrupted by contact mode transitions and impacts.

C. Contact Modes
Contact modes are defined based on the signs of γ i , γi , and σi .In particular, each individual contact point of a planar model is in one of the following states: free flight (F), stick (S), and slip in either one of two directions (P,N); see Table I.Thus, the system has 4 c contact modes.Each contact point and contact mode has a set of kinematic admissibility constraints, which are equalities or inequalities in the state variables.These constraints reduce the number of contact modes to be considered in a given state of the system.For example, Example 2 has 16 contact modes, which are two-letter words composed of the set of letters {P, N, F, S}.Among these, exactly ten modes are kinematically admissible in some parts of the state space.For example if cos φ 1 , cos φ 2 > 0, then modes P N, NP , P S, SP , SN , and NS are kinematically inadmissible for all the nonstatic states of the object.The remaining ten modes are admissible in some parts of the state space.
In addition to admissibility conditions, each contact mode delivers exactly two equality constraints involving contact forces (λ Ni , λ T i ) and/or accelerations v, which are combined with the manipulator equation in order to determine the instantaneous values of the contact forces.Finally, a second set of inequality constraints involves contact forces and accelerations.Those consistency constraints are tested after the contact forces and accelerations have been determined.
We consider an equilibrium state, i.e., a static state x 0 = 0, in which the contact mode with all the contact points in the stick state yields a consistent solution.For example, x 0 = 0 is the equilibrium of Example 1, if 0 ≤ α ≤ 90 • and μ ≥ tan α.For Example 2 with φ 1 = φ 2 = 0 (biped on a slope), it is a necessary condition for x 0 = 0 to be an equilibrium [25] that We restrict our attention to those systems in which kinematic constraints ensure immobility if all the contacts are in S state.This is true for the two examples introduced above; however, it is not true for example if a planar rigid body has only one contact point.
An important limitation of the hybrid dynamics approach is the nonuniqueness and nonexistence of consistent contact modes in some systems [26].We will distinguish between two types of this situation.The first one is given as follows.
Definition 1: An equilibrium state x 0 is called ambiguous, if the system has (in addition to the sustained equilibrium state) a nonstatic consistent contact mode as well.

TABLE I CONTACT MODES OF A SINGLE CONTACT POINT
Ambiguity is quite common among multicontact equilibria, and it has been proven that ambiguous equilibria are never Lyapunov stable (LS) [25].Conditions of ambiguity for Example 2 were also given by Or and Rimon [25].Example 1 is never ambiguous.
The second type involves solution nonexistence as well as all the forms of nonuniqueness in nonstatic states.It is often referred to as Painlevé paradox [6].Nonexistence is resolvable by considering impulsive contact forces in nonimpacting states.Nonuniqueness is, in general, not resolvable within the framework of rigid models.Nevertheless, with the approach that we develop in this article, it is possible to verify Lyapunov stability even when the hybrid dynamics approach fails to identify a unique solution.The second important limitation of the hybrid dynamics approach is complexity.The number of contact modes is exponential in c, which explains why this approach is limited to moderate values of c.

D. Contact Mode Transitions and Impacts
Continuous motion in a given contact mode is not possible unless the corresponding constraints are satisfied.If a consistency constraints or an admissibility constraint related to tangential velocities is violated, a contact mode transition occurs.For example, the violation of λ Ni ≥ 0 triggers liftoff at point i, and the violation of σi > 0 triggers slip-stick transition or slip reversal at point i.These transitions are marked by continuous but nonsmooth velocity functions.
In contrast, if any of the kinematic admissibility constraints related to γ i and γi is violated (in other words, if the boundary of the admissible set A is reached), then an impact occurs.The hybrid dynamics approach requires rigid impact models, which treat impacts as instantaneous velocity jumps generated by instantaneous impulses at the contact points.
Hence, velocity becomes a discontinuous piecewise smooth function of time, with nonidentical left and right limits v − (t), v + (t) at impact times.An impact is modeled by a map assigning to each possible preimpact state (q(t), v − (t)) with γ i (q(t)) = 0 and γ− i ≤ 0 for some i an admissible postimpact state with an updated velocity value v + (t) satisfying the condition γ+ i ≥ 0. For an impact at a single point i, many algebraic impact models exist.We will adopt the classical Whittaker-Kane-Levinson model [27], which assumes that the Newtonian coefficient of restitution parameter e is known in advance.Hence, the normal component of the postimpact normal velocity is given by In addition, it is assumed that Coulomb's law is satisfied in the following form: where Λ and Λ T are signed normal and tangential impulses transferred at point i, respectively.The assumptions outlined above are combined with a discrete-time version of the manipulator equation By combining these equations, the impact map can be expressed in closed form [28].The piecewise linearity of Coulomb's law gives rise to three types of impact (sticking and slipping in positive or negative direction).Exactly one of them is consistent with (21) for any potential preimpact state.Simultaneous impacts at several contact points can be modeled in various ways.For example, iterative solution techniques of the numerical time-stepping methods of contact dynamics compute the final states of multipoint impacts by considering sequences of impacts with single points, where the subsequent impact points are chosen according to a priori rules [29].In this work, we adopt the idea that multipoint impacts are equivalent to some (finite or infinite) sequence of single-point impacts, but we do not make any particular assumption about the order of impact points in those sequences.This assumption allows us to restrict our attention to single-point impacts.We note, however, that more general approaches to simultaneous impacts also exist [30], [31], and they are compatible with the optimization tools discussed in the following [21].

E. Zero-Order Dynamics
In a small neighborhood of the equilibrium state x 0 = 0, the dynamics can be approximated by modified versions of ( 6), (22), and of the constraints, in which the gap and tangential displacement functions are linearized around q = 0, and all other state-dependent terms are approximated by their nominal values at (q, v) = (0, 0) This approximation will be referred to as ZOD.
In the case of Example 1, equations of motion under the ZOD are identical to the exact equations.For Example 2, H and c are left unchanged by the ZOD approximation.J T γ , J T σ , and the functions ( 10) and ( 11) are approximated by where the previously introduced shorthand notations ( 14) and ( 15) have been used.The ZOD delivers linear admissibility conditions and constant accelerations and contact forces in each contact mode.On the one hand, the ZOD allows us to solve the equations of motion in all the contact modes in closed form, which will be exploited later.On the other hand, the consistency conditions of each contact mode related to accelerations are state independent: they can be verified for each contact mode in a single step.The ZOD approximation is also applied to impact maps.
The ZOD approximation yields a close approximation of the real dynamics in a small neighborhood of state x 0 = 0. Nevertheless, this approximation breaks down for systems with marginally admissible contact modes (see Fig. 3).Consider Example 2 with φ 1 = ±π/2, l 1 > l 2 , and h = 0 in a state such that point 2 is immobile, θ > 0, and θ = 0.The contact mode P S is marginally kinematically admissible under the ZOD, but the same mode under the exact dynamics is not admissible as γ 1 > 0. A similar problem may occur in relation with marginally consistent contact modes (see Fig. 3) as in the case of Example 2 with α = π/2, T = 0, φ 1 = 0, and −π < φ 2 < 0. The ZOD predicts γ2 > γ1 = 0 in the F F mode.Hence, F F is marginally consistent in the state q = v = 0.According to the exact equations, the same system exhibits γ1 > 0 in each state with q = 0 and θ = 0 if the geometric parameter h is strictly positive, but it exhibits γ1 < 0 if h < 0. Hence, consistency under the exact dynamics in a small neighborhood of x 0 = 0 cannot be predicted based on the ZOD.
In what follows, we exclude the marginal cases mentioned above, and we consider dynamics under the ZOD.It is not proven formally that local stability under ZOD implies local stability under the exact equations of motion.Nevertheless, the numerical techniques used for stability verification are robust against small numerical error; hence, it is plausible to assume that the (locally) small error introduced by the ZOD approximation does not invalidate the stability tests with the exception of the degenerate cases outlined above.

F. Stability
We will use the standard concept of Lyapunov stability.For continuous-time systems, Lyapunov stability is defined as follows.
Definition 2: Let x 0 be an equilibrium state.This configuration is called LS if for every arbitrarily small > 0, there exists δ( ) > 0 such that for any kinematically admissible initial state x(0) that satisfies |x(0) − x 0 | < δ, the emerging motion trajectory satisfies |x(t)| < for all t > 0.
The trajectories of impacting (hybrid) systems are piecewise smooth but discontinuous due to impacts.The definition of Lyapunov stability outlined above is applicable to such systems.
Contact systems with translational symmetry often have continuous sets of identical equilibrium states, which appear as single points in a reduced state obtained by removing cyclic position coordinates.In this case, one can also consider Lyapunov stability in the reduced state space.For example, both Examples 1 and 2 have translation symmetry with the equilibrium sets given by z = 0 and z = θ = 0, respectively.In both the cases, the equilibrium set shrinks to a single point in a reduced space lacking the y coordinate (but including u = ẏ).For these systems, the stability of individual points in the full state space implies stability in the reduced space.The opposite is not necessarily true, since slow "creep" motion along the y-direction while all other state variables remain close to 0 is allowable in the second case.The two notions of stability appear to be equivalent in many cases, but, in general, their exact relation is unknown according to the best of the author's knowledge.
Lyapunov stability is most often demonstrated with the help of Lyapunov's direct method, i.e., by constructing a scalar Lyapunov function V (x) over full or reduced state space such that V has a local minimum point at the equilibrium state under investigation, and V is nonincreasing along motion trajectories.The adaptation of Lyapunov's direct method to impacting systems [32] is summarized in the following.
Let B ρ denote a ball of radius ρ in the state space centered around x 0 .Furthermore, let cl denote the closure of a set.For example, cl(A) is the union of the admissible set and possible preimpact states, where γi < γ i = 0 for some i.
Theorem 1 (see [32,Th. 6.23]):If there exist: 1) a positive value h; 2) a continuous strictly increasing scalar function α with α(0) = 0; 3) a continuously differentiable function V : R 2n → R such that for any initial state V is nonincreasing along any continuous piece of trajectory through x and at jumps associated with impacts then x 0 is stable in the sense of Lyapunov.The proof of this statement is reviewed in the following, as its logical steps will be reused during the development of the extended stability theory.
Proof: Consider an arbitrary > 0. Let c = α(min( , h)).Then, V (x) ≥ c along the boundary of B min(h, ) ; hence, V (x) has a closed connected sublevel set Ω c = {x : V (x) ≤ c} such that Ω c contains x 0 in its interior; furthermore, Ω c ⊂ B min(h, ) .Owing to condition 3b), Ω c is positively invariant set of the dynamics (see [32,Proposition 6.5]).Then, we can find a scalar δ > 0 such that B δ ⊂ (Ω c ∪ A c ), where c means complement.Now, if the initial point of a trajectory is in B δ , then the positive invariance of Ω c implies that the trajectory remains in B at all times, which proves stability.

G. Analytic Conditions of Stability
The two examples introduced above are among the few ones, for which exact or sufficient stability conditions are available.
Consider first Example 1. Owing to its simple geometry and single unilateral contact point, the equations of motion can be solved, and Lyapunov stability can be assessed by hand.The point mass has four contact modes (see Table I).During free flight, the normal acceleration is constant.The point mass accelerates toward the slope if α < π/2, and the impact time can be expressed in the closed form.In the case of ideally inelastic impacts, the first impact is followed by sustained contact (left or right slip or stick).In the case of partially elastic impacts (0 < e < 1), a sequence of exponentially decaying impacts terminating in a Zeno point also leads to sustained contact.The acceleration during slip is also constant, and slip velocity decays if μ > tan α.It is easy to show that individual equilibrium points in the full state space are LS if and only if Recall that Example 1 has translation symmetry (see Section II-F); hence, stability can also be investigated in a 3-D reduced state space.In the reduced state space, the equilibrium set shrinks to a single point, which is LS under (30), and additionally in the marginal cases of μ = tan α and e = 1.Note that in the marginal cases, a small perturbation may induce slow sustained slip (if μ = tan α) or bouncing (if e = 1) motion in the downhill direction; hence, the equilibria are not stable in the full state space.While Example 1 may seem somewhat artificial, it is important to note that under the ZOD, it is equivalent to a wide class of nonlinear planar systems with a single contact point, including, for example, a model of a perching glider (see Fig. 4) and of the woodpecker toy [33].For all those systems, a linear transformation of state variables under the ZOD gives rise to equations of motion analogous to Example 1, with the only difference being that the mass matrix H is, in general, an arbitrary positive-definite matrix rather than the identity matrix.
Next, consider Example 2, a planar rigid body with arbitrary geometry on two contact points.It covers many well-known model problems in contact dynamics, including, for example, the Painlevé-Klein problem, an early demonstration of Painlevé's paradox [34]; a rigid falling ladder [35]; planar models of peg-in-hole insertion [36]; and the rimless wheel (see Fig. 4).This system has 16 possible contact modes, yet the kinematic constraints associated with rigidity leave only ten of them admissible.Example 2 shows highly nontrivial behavior.Several mechanisms of instability have been uncovered, and almost exact conditions of stability have been found if impacts are inelastic [7], [8].Little is known about stability in the more general case of e > 0; however, some highly conservative stability conditions have been reported in [9] and [10].

H. Semidefinite Programming and SOS Polynomials
According to Theorem 1, proving Lyapunov stability with the aid of Lyapunov's direct method requires the construction of functions, which are provably nonnegative under certain inequality constraints.In general, there is no efficient computational method to test the nonnegativity of a function or even of a polynomial function [37].
However, for a polynomial function to be nonnegative, it is sufficient to prove that it can be written as the SOS of some polynomials, i.e., it is SOS.In addition, all the nonnegative polynomials can be approximated by SOS polynomials [38].A polynomial is SOS if and only if it can be recast as a generalized quadratic form involving a positive-semidefinite matrix.Owing to this connection between SOS polynomials and positive-semidefinite matrices, testing if a polynomial is SOS can be formulated as a convex optimization task, for which efficient numerical implementations exist [39].The theory of SOS polynomials has also been extended to optimization algorithms over SOS polynomials, all of which have dual problems amenable to semidefinite programming [40].Hence, SOS polynomials are highly useful tools for the verification of Lyapunov stability.Polynomials, which are positive under polynomial equality or inequality constraints, can also be constructed by searching for unconstrained SOS polynomials by using the so-called S-procedure (see Table II).A combination of equality and inequality constraints can also be treated in a similar fashion, but this straightforward extension is not shown in the table.

I. Existing Lyapunov Stability Test Using SOS Programming
In a recent paper, Posa et al. [21] pointed out that the admissibility and consistency constraints of impacting systems are often polynomial functions of appropriately chosen variables.Hence, Lyapunov's direct method can be implemented as an SOS optimization problem.
In addition, it was also shown that for inelastic impacts, condition 3b) need not be prescribed for each contact mode separately, but instead a reduced set of inequality constraints can be verified (see [21,Th. 4]).This finding was termed separability of contacts.
Even though the replacement of nonnegativity by the SOS property is conservative, the Lyapunov stability of several systems was successfully verified by Posa et al. [21].
The set of examples investigated by Posa et al. [21] (see Fig. 4) includes a point mass (Example 1) in the special case of α = e = 0 and the "rimless wheel," which is a special case of Example 2 (see Fig. 2) with φ 1 = φ 2 = α = 0, l 1 = −l 2 , τ = 0, and e = 0.In addition, a modified version of the rimless wheel termed a balancing robot as well as a simple model of a perching glider is also investigated.
All of these examples belong to special subclasses of problems, where unstable equilibria in the sense of Lyapunov are unlikely to occur.In particular, the point mass and the rimless wheel possess translational symmetry.Their equilibria correspond to a local minimum of potential energy in the reduced state space (see Section II-F).The dissipative nature of friction and impacts guarantee Lyapunov stability in the reduced state space.The balancing robot has an added actuator, which may be dissipative or nondissipative depending on the control law.In the first case, Lyapunov stability remains obvious, whereas in the second case, a badly designed controller could, in principle, destabilize the equilibrium.The optimal control laws found in [21] by an SOS program are not strictly dissipative, yet Lyapunov stability is preserved.Hence, this example represents the verification of stability in a nontrivial case.Finally, the perching glider has a continuous set of nonidentical equilibria.One of them corresponds to a local minimum of potential energy, i.e., it is obviously stable.The stability of this point was confirmed by an SOS program, but the Lyapunov stability of other equilibrium points was not tested systematically in [21].
It should also be emphasized that the novel contributions of [21] extend far beyond Lyapunov stability verification, as the same method also provides conservative estimates of the basin of attraction of a stable equilibrium as well as proof of positive invariance of some subsets of S. All of these questions are highly challenging, and they are beyond the scope of this work.Instead, we will focus on the verification of Lyapunov stability in more challenging situations.

III. VERIFICATION OF STABILITY USING THE EXISTING METHODS
In order to test applicability in more challenging situations, the stability test proposed by Posa et al. [21] has been implemented using the Mosek solver [41] available through the Yalmip toolbox of MATLAB [42], and the previously introduced two examples have been tested.

A. Example 1: Point Mass on an Inclined Slope
The stability of individual equilibrium points in the full state space has been tested for several combinations of α and μ in the case of e = 0.The SOS program was evaluated with the following settings.
1) The Lyapunov function V was parameterized as a quartic polynomial in the state variables y, z, u, and w. 2) An explicit lower bound V (x) ≥ V min (x) = 0.1(y 4 + z 4 + u 4 + v 4 ) was prescribed implying conditions 2) and 3a) of Theorem 1. 3) Separability of contacts (see [21,Th. 4]) was used to develop a reduced set of nonnegativity conditions implying condition 3b) of Theorem 1. 4) All the conditional inequalities were expressed as SOS conditions using the S-procedure formulation presented in Table II.Altogether, seven auxiliary polynomials (see II) in the state variables and an additional variable representing tangential to normal contact force ratio were used.All of these polynomials have also been parameterized as quartic polynomials.5) All the conditions were prescribed within a unit ball |x| ≤ 1.This version of the program could not verify the stability of equilibrium points except in the trivial case of α = 0.In a second attempt, stability in the reduced state space was also tested using a very similar program, in which y was removed from the set of variables of all the polynomials.The modified version was successful: stability was verified by the program in all the cases satisfying (30).
This simple example illustrates the potential of SOS-based stability tests to verify stability in energetically nontrivial cases.At the same time, the failure of the test when applied in the full state space also highlights the fragility of the method even in very simple cases.It is also noteworthy that the verification of instability is not possible using this approach.

B. Example 2: Biped on an Inclined Slope
Next, Example 2 with an underlying inclined slope (φ 1 = φ 2 = 0) was tested under the ZOD for many combinations of model parameters.As before, points of equilibrium in the full (6-D) state space as well as points in a (5-D) reduced state space were both tested.The implementation of the SOS program followed principles similar to the stability test of Example 1. Separability of contacts [21] was used to decrease the number of SOS conditions.The lower bound of V was given by V (x) ≥ Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
This time, the S-procedure formulation involved altogether 16 auxiliary polynomials in the variables of V as well as two additional variables.The auxiliary polynomials have also been parameterized as quartic polynomials.
In the reduced state space, the SOS program successfully verified stability in the trivial case when the slope angle was 0. However, for α > 0, all the tests proved unsuccessful: no certificates of stability were found by the algorithm; moreover, the corresponding semidefinite programming task was in most cases flagged by the Mosek solver as provably infeasible.In the full state space, all the tests (including the special case of α = 0) were unsuccessful.

IV. EXTENDED STABILITY THEORY
The negative results outlined in Section III suggest that the original stability test has limited ability to verify stability in nontrivial situations.This observation inspires the development of an improved stability test, which can address a wider class of systems.
In a previous work of Várkonyi et al. [10], a conservative condition of stability was developed for Example 2 in the case of φ 1 = φ 2 = 0.That work used a Lyapunov-type function composed as the envelope of two smooth functions: ).The first one was the total mechanical energy of the system, and V 2 (x) was found by trial and error.Analytical conditions were derived, under which V 2 was decreasing over time in each contact mode except F F .For the F F mode, it was proved that V 2 may increase; however, each episode of free flight is followed by an impact, at which V 2 drops to a lower value such that the net change of V 2 is negative.Those properties could be used to prove Lyapunov stability in some cases.In this article, these ideas are combined with the results of [21] in order to develop a less conservative algorithmic stability test, in which manual search by trial and error is replaced by algorithmic search using a semidefinite program.Next, the proposed extensions of Lyapunov's direct method are introduced informally.This is followed by the precise statement of the new stability condition.Throughout the rest of this article, the ZOD approximation is considered.

A. Lyapunov Function V (x) May Increase Temporarily Along Paths
Multibody systems with unilateral contacts often exhibit motion trajectories, which are characterized by accelerating motion in some contact modes combined with significant loss of energy during impact as well as in other contact modes.This feature limits the applicability of Lyapunov's direct method, relying on functions V (x), which are strictly nonincreasing function along trajectories in all the contact modes.There might be situations in which a proper Lyapunov function does not exist; nevertheless, the high amount of energy loss during impacts still ensures stability.The extension of Lyapunov's direct method proposed in the following is tailored to exploit this general feature of contact systems.A similar observation in the context of neural nets was proposed by Pavlidis [43], and similar ideas are proposed within the stability theory of switching systems using multiple Lyapunov functions [17].
We propose an improved condition of stability, which allows V (x) to increase temporarily along solution trajectories, provided that it has a nonincreasing overall trend.In particular, we will prescribe that the values of V (x + (t i )) at impact times t i , i = 1, 2, . .., form a strictly nonincreasing sequence; furthermore, the time differences t i+1 − t i are sufficiently small to prevent the escape of trajectories from the proximity of the equilibrium state during the time interval (t i , t i+1 ).

B. Stability Certificates via a Nonsmooth Lyapunov Function
There are several extensions of Lyapunov's direct method, which make use of piecewise smooth Lyapunov functions.Nonsmooth points of the Lyapunov function typically coincide with switching points and impacts of the original system [15].Here, we take a slightly different approach and look for the upper envelope of a finite set of functions V i (x) (i = 1, 2, . .., n), none of which satisfies condition 3a) of Theorem 1.Nevertheless, it is required that Furthermore, each one of the functions satisfies condition 3b).It is not assumed that the switching manifolds are the same as the switching points of contact modes.
Stability tests using SOS programming can benefit from such an extension, as the upper envelope of two polynomial functions is usually not a polynomial.Hence, the use of several Lyapunov functions makes a wider set of Lyapunov candidates available for the solver.
Notably, a second Lyapunov function can be used to inform the solver about the dissipative nature of frictional dynamics.One can choose V 1 (x) as the total mechanical energy of the system and search for a second function V 2 (x), which together satisfy (31).

C. New Stability Condition
Now, we are ready to state the main result of this article.Theorem 2: Let x 0 be an unambiguous equilibrium state.If there exists: 1) a positive value h; 2) a continuous strictly increasing scalar functions α with α(0) = 0; 3) a finite set of continuously differentiable functions V i : R 2n → R (i = 1, 2, . .., υ) such that for any initial state , where x * denotes the next point of contact mode transition or the next postimpact state along the trajectory through x, whichever occurs earlier then x 0 is stable in the sense of Lyapunov.
Clearly, Theorem 2 is a generalization of Theorem 1.The latter is recovered if one requires υ = 1, and if condition 3) of Theorem 2 is replaced by the stronger requirement of monotonicity of V i (x).
The proof of Theorem 2 follows the logical steps of Theorem 1.As has been pointed out, the relaxed condition 3a) in the statement does not require any significant change in the proof, as the nonsmooth function max i V i (x) satisfies the requirements of the original theorem.The relaxed condition 3b) requires some adaptation of the proof.
The proof relies on two crucial properties of contact systems under the ZOD.The first one is an upper bound of time differences between contact mode transitions and impacts.
Lemma 1: For every planar contact system such that x 0 = 0 corresponds to an unambiguous equilibrium, there exists a continuous strictly increasing function ζ(ξ) : R → R with ζ(0) = 0 such that for any contact mode M and for any state x ∈ S in which mode M is admissible and consistent, the next contact mode transition or impact along the trajectory through x occurs no later than at time t ≤ ζ(|x|).
The second one is a bound for trajectories completed in short time intervals.
Lemma 2: There exists a continuous strictly increasing function η(ξ) : R → R with η(0) = 0 such that for any trajectory x(t) under the ZOD, which includes no impacts within a time interval (t 1 , t 2 ) and |x(t 1 )| < 1, the following bound is satisfied under the ZOD: The proofs of the lemmas are presented in Appendixes A and B.
Proof of Theorem 2: Consider an arbitrary 0 < < 1. Define c as c = α(min( , h)).Then, V (x) ≥ c along the boundary of B min(h, ) ; hence, V (x) has a closed connected sublevel set Ω c : V (x) ≤ c such that Ω c contains x 0 in its interior; furthermore, Ω c ⊂ B min(h, ) .Next, we can find a scalar β > 0 such that B β ⊂ (Ω c ∪ A c ), where c means complement.
Since η and ζ (see Lemmas 1 and 2) are strictly monotonic continuous functions and η(0) = ζ(0) = 0, there is a unique positive value of δ satisfying the equation If the initial point of a trajectory at time t 1 is in B δ , then by Lemma 1, the next contact mode switch or impact occurs no later than at time According to Lemma 2 and (32), the state of the object remains within the bound given by |x(t) − x 0 | ≤ β during the time interval t ∈ (t 1 , t 2 ), which is a subset of B .By using condition 3b), we also know that |x + (t 2 )| ≤ |x(t 1 )|, i.e., the state |x + (t 2 ) is also inside B δ just like the initial state x(t 1 ).This observation allows the recursive application of the previous arguments for each phase of motion free of impacts and contact mode transitions.This way, we arrive to the conclusion that the state of the systems remains within B at all times, implying Lyapunov stability.

V. VERIFICATION OF STABILITY USING THE EXTENDED STABILITY TEST
In this section, we apply the extended stability test to Example 1 (point mass on a slope) and Example 2 (planar biped on a slope).

A. Implementation of the Stability Conditions
Four contact modes of Example 1 and ten contact modes of Example 2 may be admissible or inadmissible depending on state.For each of them, the values of contact forces and instantaneous accelerations of the system are state independent under the ZOD, and thus, they can be calculated in advance.For example in the F F contact mode of Example 2, there are no contact forces, and the acceleration is given by the value of c [see (13)].After the preliminary steps outlined above, the conditions of Theorem 2 are formulated as an SOS program and tested.Such an implementation requires the following.

1) Specification of h in condition 1):
We use h = 1.
2) Lyapunov function candidates: For both the examples, we choose V 1 (x) to be total mechanical energy, i.e., for Example 1 and for Example 2. V 2 (x) is chosen as for Example 1 and for Example 2. Here, P D 1 denotes a polynomial containing all the terms with degrees within the interval (1, D).V 2 does not have constant terms due to the required property V (x 0 ) = 0. We chose D = 2.In addition, the nonlinear terms in the variable y are omitted.This choice has no strict theoretical reasons; however, it is plausible, because y is a cyclic coordinate in both the cases, and the system has a continuous set of identical equilibria parameterized by y.The lack of nonlinear terms containing y ensures that the derivatives of V 2 are independent of y; thus, a Lyapunov function verifying stability of one equilibrium point also verifies the stability of all others.With the restrictions above, V 2 (x) has altogether nine unspecified coefficients for Example 1 and 25 for Example 2.
3) Existence of an appropriate function α(ξ) satisfying condition 3a): At this point, we can exploit the fact that uphill motion is bounded by increase of the potential energy term in V 1 (x).In turn, V 2 (x) need not be positive definite; it should, however, increase within the region of state space characterized by negative potential energy.In Appendix C, we prove that in the case of Example 1, the bound V 2 (x) ≥ y (37) and in the case of Example 2, the bound imply condition 3a) of Theorem 2. This inequality is added to the set of nonnegativity conditions to be fulfilled.4) Implementation of condition 3b) for all possible forms of motion: The function V 1 is nonincreasing; hence, all we need is to prove monotonicity of V 2 .If e = 0, then the motion does Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
not have free flight phase.In this case, we use the following conditions.a) d dt V 2 (x(t)) ≤ 0 for all of the nonstatic contact modes involving an active contacts (that is {P, N } for Example 1 and {P F, SF, NF, F P, F S, F N, P P, NN} for Example 2) and for all the states x ∈ B h ∩ cl(A), where that the contact mode is admissible and consistent.b) V 2 (I(x − (t))) − V 2 (x − (t)) ≤ 0 for all the possible preimpact states x − (t) ∈ B h ∩ cl(A), where I is the impact map.If e > 0, then free flight is also possible.Free flight is always followed by an impact, and we also use the following additional condition: a) V (I(x − (t))) − V (x − (t − Δ)) ≤ 0 for all possible preimpact states x − (t) and for all possible durations Δ of a free flight phase immediately before that impact.As before, I means the impact map.All of these conditions as well as the corresponding conditions of admissibility and consistency are given in closed form for Example 2 in Appendix D. All of them are polynomial equations or inequalities in an appropriate set of variables (which includes state variables as well as some additional variables such as Δ in the case of free flight).Hence, the whole set of conditions can be recast as an SOS programming task, as demonstrated in Section II-H.The SOS program involving all the unknown coefficients generated via the S-procedure has altogether 479 unknown coefficients for Example 1 and 916 for Example 2.

B. Test Results for Example 1
In the case of Example 1, all the stability tests with parameter values satisfying (30) were successful.Hence, unlike the original SOS-based test, proposed by Posa et al. [21], the extended stability theory enables the verification of the stability of individual equilibrium points in the 4-D full state space.

C. Test Results for Example 2
Here, we report on the results of the extended stability test and compare them with previously published results on Lyapunov stability.
First, the stability test was executed for slope of angle 30 • (i.e., φ 1 = φ 2 = 0 and α = 30 • ) and gravitational load (T = 0 and F = 1) with inelastic impacts (e 1 = e 2 = 0).The positions of the contact points were chosen as h = 1, l 1 = h tan α + d, and l 2 = h tan α − d.The friction coefficient μ 2 = 1 was kept fixed, while μ 1 and the distance 2d between the two contact points were varied systematically.This particular case of Example 2 has been investigated previously.If μ 2 > 2 tan α, then the object is in frictional equilibrium for all the values of μ 1 and d.Nevertheless, the equilibrium is unstable due to ambiguity [25] if Even when μ 1 is above the limit (39), the object may be unstable due to self-excited inverse chatter motion as first pointed out by Várkonyi et al. [10].The exact range of parameter values corresponding to instability has been determined by semianalytical investigation of an appropriately defined Poincaré map [7], [8].The newly proposed Lyapunov stability test is able to verify stability in a significant portion of the stable region of model parameters (see Fig. 5).As a second step, we examined parameter values The SOS stability test was not able to verify stability in the first case, but a large portion of the equilibrium region was verified to be stable in the second case.All the sets of results presented so far indicate that the new test is conservative, but it is successful in some nontrivial cases.
As we have seen, the exact conditions of Lyapunov stability are not available from the literature for partially elastic impacts.Now, we return to l 1 = h tan α + d and l 2 = h tan α − d, but the slope angle is chosen as α = 20 • .The coefficients of friction are equal μ 1 = μ 2 , and they are varied systematically along with the coefficients of restitution e 1 = e 2 as well as d.The system is in frictional equilibrium if μ i ≥ tan 20 • ≈ 0.36, and the equilibrium is unambiguous in all the cases.The exact parameter range implying Lyapunov stability is unknown.This example has been investigated by Várkonyi et al. [10], where stability was verified for low values of the coefficients of restitution in a range of model parameters illustrated by Fig. 7 (see also [10,Fig. 4]).We tested stability using the new SOS-based algorithm.For e i = 0, we test verified stability in the same region as the analytical condition of [10].However, the range of stability verified by the extended SOS-based stability test appears to be independent of the coefficients of restitution within the range 0 ≤ e i ≤ 0.8.This finding is remarkable because the analytical method of [7] and [8] is not applicable when e i > 0, and the stability conditions  of [10] fail to verify stability for e i ≥ 0.4.This feature of the new test shows that an appropriately chosen Lyapunov function (significantly different from total energy of the system) may successfully verify stability even when the energy absorption associated with impacts is moderate.
One possible way to improve the performance of SOS-based stability tests is to increase the maximum degree D in the polynomial parameterization of the Lyapunov candidates (36) and the auxiliary functions used in the S-procedure.Many of the tests outlined above have been repeated with D = 3 and D = 4 for all the polynomials.However, the performance of the stability test did not improve for this particular example.Moreover, weaker performance was observed in a few cases, which is probably due to the excessive computational complexity of the semidefinite program in these cases.

VI. CONCLUSION
The Lyapunov stability analysis of rigid multibody systems with unilateral frictional contacts is a fundamental question of robotics and object manipulation.Recent theoretical and experimental works demonstrated the existence of various types of instability in response to small perturbations initiating liftoff at the contacts.Despite successful attempts to verify the stability of some simple model systems, general stability tests, which provide an answer in all the cases, remain unavailable.The complexity of contact-induced rigid-body dynamics makes analytical investigation infeasible even for seemingly simple systems.Algorithmic stability tests have great potential in more complex cases.An important step in this respect was made by Posa et al. [21], where the compatibility of SOS programming with the special version of Lyapunov's direct method applicable to hybrid dynamical systems induced by rigid contact was recognized.Despite the success of this method in several related problems, it fails to verify the Lyapunov stability of the two model problems investigated here.In this article, several improvements of the method were proposed, which enabled the successful verification of Lyapunov stability in those cases.It was found that the algorithmic test outcompeted a previously proposed manual stability test based on physical intuition.At the same time, the stability test appeared to be conservative.
It is foreseen that future developments will radically extend the applicability of algorithmic stability tests to contact-induced dynamics.A fundamental challenge to be solved is the reduction of high computational complexity of the algorithm (currently exponential in the number of contact points), which emanates from the complexity of the hybrid dynamics approach to contact problems.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
where J V (x) is the Jacobian of the function V (x), which can be expressed in terms of the unknown coefficients of V (x).
Hence, the conditions (50) for positive and negative slips together imply (51).P P and NN modes: In this case, the two normal contact forces are determined by the constraints γ1 = γ2 = 0 [λ 1 , λ 2 ] T = −P −1 [cos α, cos α] T where If any of the two contact force values is negative, then the conditions associated with the contact mode under investigation can be disregarded.In the opposite case, we have The stability conditions are formulated as −J V (x) ẋ± ≥ 0 whenever γ 1 (q) = γ 2 (q) = 0 γi (q) = γ2 (q) = 0 Impacts: Slipping and sticking impacts can be addressed in a common framework as follows.Let i be the contact point undergoing an impact.Let Λ and Λ T denote net impulses in the normal and tangential directions during that impact, respectively.
Λ can be expressed in terms of the preimpact state and of Λ T by using the kinematic condition (20) and ( 22) Then, the postimpact state becomes The stability conditions are formulated as V (x − ) − V (x + ) ≥ 0 whenever γ i (q) = 0, γi (x − ) ≤ 0, γ j (q) ≥ 0 The last three conditions express Coulomb's law for the impact.The conditions listed above are polynomial in the set of variables {y, z, θ, u, w, ω, Λ t }.F F mode followed by an impact: Consider a motion trajectory of a system with partially elastic impacts (e > 0), which includes an episode of free flight of duration Δ followed by an impact at time t.
During free flight, contact forces are inactive, and the (constant) acceleration of the object under the ZOD is vFF = [sin α, − cos α, 0] T .Let x − = (q, v − ) denote the preimpact state.The state of the object at the beginning of the free flight phase under the ZOD is expressed in terms of x − by using the kinematics of uniformly accelerating motion whereas the postimpact state is given explicitly by (52).Condition 3b) of Theorem 2 is expressed as V (x + (t − Δ)) − V (x + ) ≥ 0 whenever γ i (q) = 0, γi (x − ) ≤ 0, γ j (q) ≥ 0 The last set of conditions expresses that the duration of free flight must be nonnegative; at the same time, it cannot be longer than −2 γi (q)/ cos α; otherwise, contact point i would be under the ground at the beginning of the free flight phase.The conditions listed above are polynomial in the set of variables {y, z, θ, u, w, ω, λ t , Δ}.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

An
Improved Lyapunov Stability Test of Equilibria Under Frictional Unilateral Contact by Sums of Squares Programming Péter L. Várkonyi Abstract-Reliable quasi-static object manipulation and robotic locomotion require the verification of the stability of equilibria under unilateral contacts and friction.In a recent paper, Posa et al. (2016) demonstrated that sumsof-squares (SOS) programming can be used to verify the Lyapunov stability of planar systems via Lyapunov's direct method if impacts are inelastic.However, this method appears to be too conservative to verify the stability of some remarkably simple examples.In this article, an extension of Lyapunov's direct method is proposed, which makes use of a piecewise smooth Lyapunov function and allows a temporary increase of the Lyapunov function along a motion trajectory.In addition, a modified SOS formulation enables the investigation of planar systems with partially elastic contacts.The proposed method remains compatible with SOS programming techniques.The improved stability test is successfully applied to a point mass on a slope and to a rigid body with two contact points.For the latter, several mechanisms of instability have been demonstrated experimentally, but the exact conditions of Lyapunov stability are unknown.Index Terms-Hybrid dynamics, Lyapunov stability, semidefinite programming, sums-of-squares (SOS) polynomials, unilateral contact.

Manuscript received 9
June 2023; revised 13 June 2023 and 11 October 2023; accepted 11 November 2023.Date of publication 16 November 2023; date of current version 30 May 2024.This work was supported by the National Research, Innovation and Development Office of Hungary under Grant K124002.Recommended by Associate Editor D. Efimov.

Fig. 2 .
Fig. 2. Planar rigid body on two contact points.Left: equilibrium configuration with model parameters.Right: general configuration with state variables, gap functions, and tangential displacement functions.

Fig. 3 .
Fig. 3. Left: an example of marginally kinematically admissible contact mode P S under the ZOD.Dashed contour illustrates the instantaneous velocity of the object.Right: an example of marginally consistent F F mode under the ZOD.The velocity of the object is 0. Solid arrow marks the external force acting upon it, and dashed contour illustrates the instantaneous acceleration.

Fig. 4 .
Fig. 4. Examples investigated by Posa et al. [21].(a) "bean bag": an inelastic point mass resting on a horizontal surface.(b) Rimless wheel: a planar rigid body with two symmetrically arranged contact points resting on a horizontal surface.(c) Balancing robot: a modified version of the rimless wheel equipped with a second rigid component attached to the wheel by an actuated pin joint driven by a controller.(d) Perching glider: a two-degree-of freedom model with two rigid components, a passive compliant joint, a pin joint support, and a unilateral contact point with a vertical surface.

25 ,
and l 2 = l 0 + h tan α + 0.25 for two different constant values μ 1 = {0.3,0.5}.The parameters h and l 0 were varied systematically, which corresponds to shifting of the center of mass relative to the contact points.This case has been investigated experimentally and analytically by Or and Várkonyi [8], and remarkably complex patterns of instability have been found for low values of μ 1 , such as μ 1 = 0.3 [see Fig. 6(a)].For μ 1 = 0.5, all the equilibria become stable [see Fig. 6(b)].

Fig. 6 .
Fig. 6.Equilibrium and stability charts of a planar biped with α = 25 • , μ 2 = 1, μ 1 = 0.3 (a) or 0.5 (b), e = 0, l 1 = l 0 + h tan α − 0.25, and l 2 = L 0 + h tan α + 0.25 as a function of l 0 and h.The diagram represents the variation of the position of the center of mass if the contact points are fixed at the two points marked by black circles.Labels indicate the characterization of individual points as follows: no equilibrium (0), unstable equilibrium due to ambiguity (1) or due to the inverse chatter phenomenon (2), and stable equilibrium (3).The characterization was done using the semianalytical stability test of [8].Empty circle markers in (b) indicate stability region verified by the SOS program.

Fig. 7 .
Fig. 7. Equilibrium and stability charts of a planar biped with T = 0, F = 1, α = 20 • , h = 1, φ 1 = φ 2 = 0, and l 1 , l 2 = h tan α ± d, as a function of μ 1 = μ 2 and d for various values of e 1 = e 2 .Dashed line indicates minimum value of friction coefficient for frictional equilibrium.The regions on the left of the solid curves correspond to stability successfully verified by Várkonyi et al. [10] for various values of e i .Those regions disappear for e i ≥ 0.4, but they are not monotonic in e i , and there is a sudden jump at e i = 0. Circles indicate stability region verified by the extended SOS stability test, which is independent of e i .