# Lyapunov stability analysis of rigid body systems with multiple contacts by
sums-of-squares programming

Péter L. Várkonyi *This work was supported by the National Research,
Innovation and Development Office of Hungary under grant $K124002$.P. L.
Várkonyi is with the Department of Mechanics Materials and Structures,
Budapest University of Technology and Economics, H-1111 Budapest, Hungary
<EMAIL_ADDRESS>
###### Abstract

Reliable quasi-static object manuipulation and robotic locomotion require
verification of the stability of equilibria under rigid contacts and friction.
In a recent paper, M. Posa, M. Tobenkin, and R. Tedrake demonstrated that
sums-of-squares (SOS) programming can be used to verify Lyapunov stability via
Lyapunov’s direct method. This test was successfully applied to several simple
problems with a single point contact. At the same time it has been found that
this method is too conservative for several multi-contact systems. In this
paper, an extension of Lyapunov’s direct method is proposed, which makes use
of several Lyapunov functions, and which allows _temporary_ increase of those
Lyapunov function along a motion trajectory. The proposed method remains
compatible with SOS programming techniques. The improved stability test is
successfully applied to a rigid body with 2 point contacts, for which the
exact conditions of Lyapunov stability are unknown.

##  I INTRODUCTION

Many tasks in robotics involve unilateral contact with friction between hard
objects. Small perturbations of an equilibrium state may induce stick-slip
transitions, as well as contact separation and impacts [1, 2]. Within the
framework of rigid body dynamics, non-smoothness, and discontinuity of the
response prevents 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 not associated with local minima of potential energy. Moreover, the
emerging hybrid dynamics has special features like Zeno sequences of impact
events [4, 5], and issues of non-existence and non-uniqueness [6]. As a
result, general methods to test the Lyapunov stability of equilibria under
contact and friction remain unavailable [7, 8].

Stability analysis of a planar model with one single point contact is
straightforward: an equilibrium is stable under any type of small state
perturbation if it is stable under rolling motion, and additionally the
Newtonian coefficient of restitution of light impacts in the neighborhood of
that state is below 1. In contrast, stability analysis of models with multiple
contact points becomes an open problem. In the case of 2 contact points and
ideally inelastic impacts, an almost exact semi-analytic condition has been
derived [9, 10]. Some highly conservative conditions have been found in the
case of partially elastic impacts [11, 12]. The last work used a custom-made
Lyapunov function for the verification of stability. Systems with more than 2
contact points, models with extended areas of contact, as well as 3D models
remain unexplored.

A novel algorithmic approach to stability analysis by Posa et al. [13] uses
sums-of-squares polynomials and semi-definite pogramming to develop sufficient
conditions of local Lyapunov stability of various systems. In addition, their
method delivers estimations of the basin of attraction of a stable point, and
it can also be used to verify 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 not guaranteed as the
proposed conditions are conservative. Notably, all of the models successfully
tested by [13] have only one point of contact, and it remains an open question
if more challenging problems can be addressed by this method.

In this paper we focus on Lyapunov stability, i.e. resilience against local
state perturbations. We find that the method proposed by [13] fails to verify
stability of an important class of model problems: a planar model of a rigid
body with 2 point contacts. This finding motivates the development of several
improvements of the stability theory. Most importantly, the semi-definite
programming techniques of [13] are combined with an extension of Lyapunov’s
direct method inspired by prior work of the author [12]. In its original form,
Lyapunov’s direct method uses Lyapunov functions defined over the state space
of the system, which are non-increasing along motion trajectories. The
proposed extension makes use of Lyapunov functions which may temporarily
increase along trajectories, provided that they have an decreasing overall
trend. In addition the theory is also extended here to multiple polynomial
Lyapunov functions, which improves the applicability of semi-definite
programming techniques. Finally, the present work makes use of a standard
approximation of piecewise continuous systems: the so-called zero-order
dynamics (ZOD). In each 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.

In Section II, general notation, and problem statement are introduced, and key
results of reference [13] are reviewed. This is followed by proposed
extensions of existing stability theory in Section III. Then, a family of test
problems is introduced and analyzed (Section IV). The paper is closed by
Conclusions and by pointing out related open problems.

##  II PROBLEM STATEMENT AND PREVIOUS RESULTS

###  II-A Problem statement

We consider a rigid body or a rigid multibody system in 2 dimensions with
$c\geq 1$ unilateral point contacts and $n$ degrees of freedom. The state $x$
of the system is given by the state vector $x=(q,v)$ composed of generalized
coordinates $q\in\mathbb{R}^{n}$ and generalized velocities $v=\dot{q}$ where
dot means derivation with respect to time. Hence we have a $2n$-dimensional
state space $x\in\mathcal{S}\equiv\mathbb{R}^{2n}$. We will consider an
equilibrium state given by $x_{0}=(q_{0},0)$. Without loss of generality
$q_{0}=0$ is assumed.

We use lower indices to denote elements of vectors, and columns of matrices.
For example $\gamma_{i}$ is the $i^{th}$ element of vector $\gamma$, and
$M_{i}$ is the $i^{th}$ column of matrix $M$. Furthermore $J_{\gamma}(q)$
means the Jacobian of $\gamma(q)$ if $\gamma$ is a vector or the gradient of
$\gamma(q)$ if $\gamma$ is a scalar.

The admissible set $\mathcal{A}$ is defined as a subset of the state space
given by $c$ scalar inequality constraints

$\displaystyle\gamma_{i}(q)\geq 0;\;i\in\\{1,2,...,c\\}$ (1)

where $\gamma_{i}(q)$ 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 $\gamma_{i}(0)=0$.

The unilateral contacts may give rise to non-negative normal contact forces
$\lambda_{iN}$ subject to the linear complementarity condition

$\displaystyle\lambda_{iN}\cdot\gamma_{i}=0\leq\lambda_{iN},\gamma_{i}$ (2)

furthermore if $\gamma_{i}=0$ then we also have an analogous condition at the
velocity level:

$\displaystyle\lambda_{iN}\cdot\dot{\gamma}_{i}=0\leq\lambda_{iN},\dot{\gamma}_{i}$
(3)

In the presence of friction, it is convenient to introduce tangental
displacement functions $\sigma(q)$ associated with each contact, since the
behaviors of contacts depend on the relative tangential velocities given by
$\dot{\sigma}(q,v)$. The sign of $\dot{\sigma}_{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

$\displaystyle\lambda_{iT}\in-\mathrm{sign}(\dot{\sigma}_{i})\mu_{i}\lambda_{iN}\;i\in\\{1,2,...,c\\}$
(4)

where $\mu_{i}\geq 0$ are friction coefficients associated with the contacts
and $\mathrm{sign}(\cdot)$ denotes the set-valued sign function.

The dynamics of the system is governed by the manipulator equations

$\displaystyle\dot{q}$ $\displaystyle=v$ (5) $\displaystyle
H(q)\dot{v}+C(q,v)$ $\displaystyle=J_{N}(q)\lambda_{N}+J_{T}(q)\lambda_{T}$
(6)

where the mass matrix $H$ depends on the current configuration $q$ of the
system. The term $C$ includes the effect of free forces as well as gyroscopic
terms. The two terms on the right-hand side represent the effects of normal
and tangential contact forces. If $H$ is invertible, one can express the
instanataneous acceleration $\dot{v}$ in terms of the (unknown) contact
forces.

As an example, consider a planar rigid body with two sharp vertices resting on
two straight surfaces as in Fig. 1. We fix a global coordinate frame such that
the origin conincides with the center of mass in the initial equilibrium
configuration and the $x$ 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 $\alpha$ acting at the center of mass as well as a resultant torque $T$.
Without loss of generality we assume 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}$, $l_{2}$, the angles $\phi_{1}$, $\phi_{2}$
and the friction coefficients $\mu_{1}$, $\mu_{2}$. We can use the global
coordinates $(x,z)$ of the center of mass, and the rotation angle $\theta$ of
the body as generalized coordinates, i.e. $q=(x,z,\theta)$. The parameters of
the manipulator equation and the gap functions and tangential displacement
functions are:

$\displaystyle H=\left[\begin{array}[]{ccc}1&0&0\\\ 0&1&0\\\
0&0&1\end{array}\right],J_{N}=\left[\begin{array}[]{cc}-\sin\phi_{1}&-\sin\phi_{2}\\\
\cos\phi_{1}&\cos\phi_{2}\\\ \zeta_{1}&\zeta_{2}\end{array}\right]$ (13)
$\displaystyle C=\left[\begin{array}[]{c}F\sin\alpha\\\ -F\cos\alpha\\\
T\end{array}\right],J_{T}=\left[\begin{array}[]{cc}\cos\phi_{1}&\cos\phi_{2}\\\
\sin\phi_{1}&\sin\phi_{2}\\\ \psi_{1}&\psi_{2}\end{array}\right]$ (20)
$\displaystyle\gamma_{i}(q)=-y\sin\phi_{i}+z\cos\phi_{i}-\xi_{i}+\psi_{i}$
(21)
$\displaystyle\sigma_{i}(q)=y\cos\phi_{i}+z\sin\phi_{i}-\eta_{i}+\zeta_{i}$
(22)

where the following notation has been used:

$\displaystyle\eta_{i}$ $\displaystyle=l_{i}\cos\phi_{i}-h\sin\phi_{i}$ (23)
$\displaystyle\xi_{i}$ $\displaystyle=-h\cos\phi_{i}-l_{i}\sin\phi_{i}$ (24)
$\displaystyle\zeta_{i}$
$\displaystyle=l_{i}\cos(\phi_{i}-\theta)-h\sin(\phi_{i}-\theta)$ (25)
$\displaystyle\psi_{i}$
$\displaystyle=-h\cos(\phi_{i}-\theta)-l_{i}\sin(\phi_{i}-\theta)$ (26)

Figure 1: A planar rigid body on two point contacts. Left: equilibrium
configuration with model parameters. Right: general configuration with state
variables, gap functions and tangential displacement functions.

###  II-B The 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 adresses all cases of the sign function
in a unified framework. In contrast, detailed analysis of the emerging motion
is usually done within the framework of hybrid dynamics (especially for
systems with moderate number of equilibria). 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_.

###  II-C Contact modes

Contact modes are defined based on the signs of $\gamma_{i}$,
$\dot{\gamma}_{i}$, and $\dot{\sigma}_{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 contraints_ , 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,
point $i$ cannot be in slip or stick state if $\gamma_{i}>0$.

In addition, each point delivers exactly two equality constraints involving
contact forces ($\lambda_{N,i}$, $\lambda_{T,i}$) and/or accelerations
$\dot{v}$, which are combined with the manipulator equation in order to
determine the instantaneous values of the contact forces. Finally, a second
set of inequalitiy 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,0)$, in
which the contact mode with all contact points in stick state yields a
consistent solution. We also restrict our attention to those systems in which
kinematic constraints ensure immobility if all contacts are in S state. This
is true for the example of Figure 1, however not true if the same object has
only one point contact.

An important limitation of the hybrid dynamics approach is the non-uniqueness
and non-existence of consistent contact modes in some systems [14]. We will
distinguish between two types of this situation. The first one is given by .

###### Definition 1

_An equilibrium state $x_{0}$ is called ambiguous, if the system has (in
addition to the sustained equilibrium state) a non-static consistent contact
mode as well._

Ambiguity is quite common and it has been proven that ambiguous equilibria are
never Lyapunov-stable [15]. Conditions of ambiguity for the previous example
were also given by [15].

The second type involves solution non-existence as well as all forms of non-
uniqueness in non-static states. It is often referred to as Painlevé’s paradox
[6]. Non-existence is resolvable by considering impulsive contact forces in
non-impacting states. Non-uniqueness is in general not resolvable within the
framework of rigid models. Nevertheless, with the approach that we develop in
this paper 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$.

Letter | contact mode | kinematic | equalities | consistency
---|---|---|---|---
|  | admissibility |  | constraints
S | sticking | $\gamma_{i}=0$ | $\ddot{\gamma}_{i}=0$ | $|\lambda_{iT}|\leq\mu_{i}(\lambda_{iN})$
|  | $\dot{\gamma}_{i}\mbox{$\\!\,=\,\\!$}\dot{\sigma}_{i}\mbox{$\\!\,=\,\\!$}0$ | $\ddot{\sigma}_{i}=0$ |
F | free |  $\gamma_{i}\geq 0$ and | $\lambda_{iN}=\lambda_{iT}=0$ |  $\ddot{\gamma}_{i}>0$ if $\gamma_{i}\mbox{$\\!\,=\,\\!$}\dot{\gamma}_{i}\mbox{$\\!\,=\,\\!$}0$
|  |  $\dot{\gamma}_{i}\geq 0$ if $\gamma_{i}\mbox{$\\!\,=\,\\!$}0$ |  |
P | positive slip |  $\gamma_{i}\mbox{$\\!\,=\,\\!$}\dot{\gamma}_{i}\mbox{$\\!\,=\,\\!$}0$, |  $\ddot{\gamma}_{i}=0$, |  $\lambda_{iN}\geq 0$,
|  | $\dot{\sigma}_{i}\geq 0$ | $\lambda_{iT}=-\mu(\lambda_{iN})$ |  $\ddot{\sigma}_{i}>0$ if $\dot{\sigma}_{i}\mbox{$\\!\,=\,\\!$}0$
N | negative slip |  $\gamma_{i}\mbox{$\\!\,=\,\\!$}\dot{\gamma}_{i}\mbox{$\\!\,=\,\\!$}0$, |  $\ddot{\gamma}_{i}=0$, |  $\lambda_{iN}\geq 0$,
|  | $\dot{\sigma}_{i}\leq 0$ | $\lambda_{iT}=\mu(\lambda_{iN})$ |  $\ddot{\sigma}_{i}<0$ if $\dot{\sigma}_{i}\mbox{$\\!\,=\,\\!$}0$
TABLE I: Contact modes of a single contact point

###  II-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 violation of $\lambda_{iN}\geq
0$ triggers liftoff at point $i$ and violation of $\dot{\sigma}_{i}>0$
triggers slip-stick transition or slip reversal at point $i$. These
transitions are marked by continuous but non-smooth velocity functions.

In contrast, if any of the kinematic admissibility constraints related to
$\gamma_{i}$ and $\dot{\gamma}_{i}$ is violated (in other words, if the
boundary of the admissible set $\mathcal{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 non-identical left and right limits
$v(t^{-})$, $v(t^{+})$ at impact times.

Hence, an impact is modeled by a map assigning to each possible pre-impact
state $(q(t),v(t^{-}))$ with $\gamma_{i}=0$ and $\dot{\gamma}_{i}^{-}\leq 0$
for some $i$, an admissible post-impact state with updated velocity value
$v(t^{+})$ satisfying the condition $\dot{\gamma}_{i}\geq 0$.

For an _impact at a single point_ $i$, many algebraic impact models exist. We
will adopt the classical Whittaker-Kane-Levinson model [16], which assumes
that the Newtonian coefficient of restitution parameter $e$ is known in
advance. Hence the normal component of the post-impact normal velocity is
given by: $\dot{\gamma}_{i}(q(t),v(t^{+}))=-e\dot{\gamma}_{i}(q(t),v(t^{-}))$
Given that $\gamma$ is a function of $q$, the previous relation is
equivalently expressed as

$J_{\gamma,i}(q)^{T}v(t^{+})=-eJ_{\gamma i}(q)^{T}v(t^{-}).$ (27)

In addition, it is assumed that Coulomb’s law is satisfied in the following
form:

$\Lambda_{iT}\in-\mathrm{sign}(\dot{\sigma}_{i}(t^{+}))\Lambda_{iN}\ $ (28)

where $\Lambda_{Ni}$ and $\Lambda_{T,i}$ are signed normal and tangential
impulses at point $i$.

The constraints outlined above can be combined with a discrete version of the
manipulator equation

$H(q)(v(t^{+})-v(t^{+}))=J_{Ni}(q)\Lambda_{Ni}+J_{T,i}(q)\Lambda_{T,i}$ (29)

where $J_{N,i}$, $J_{T,i}$ are the $i^{th}$ column of $J_{N}$ and $J_{T}$,
respectively. By combining these equations the impact map can be expressed in
closed form [17]. In particular, an impact can be sticking
($\dot{\sigma}_{i}(t^{+})=0$) or slipping in either direction
($\pm\dot{\sigma}_{i}(t^{+})>0$) and it can be shown that exactly one of the 3
solution candidates satisfies the consistency conditions.

The impulse of a sticking impact is then given by

$\displaystyle\Lambda_{N,i}$
$\displaystyle=\frac{-(1+e)J_{\gamma,i}(q)^{T}v(t^{-})-J_{\gamma,i}(q)^{T}H^{-1}(q)J_{T,i}(q)A}{J_{\gamma,i}(q)^{T}H^{-1}(q)(J_{N,i}+J_{T,i}B)}$
(30) $\displaystyle\Lambda_{T,i}$ $\displaystyle=A+B\Lambda_{N,i}$ (31)

with

$A=\frac{J_{\sigma i}^{T}v(t^{-})}{J_{\sigma
i}^{T}H^{-1}(q)J_{T,i}},B=\frac{J_{\sigma i}^{T}H^{-1}(q)J_{N,i}}{J_{\sigma
i}^{T}H^{-1}(q)J_{T,i}}.$

The impulse of a slipping impact in the $\pm$ direction is given by

$\displaystyle\Lambda_{N,i}$
$\displaystyle=\frac{-(1+e)J_{\gamma,i}(q)^{T}v(t^{-})}{J_{\gamma,i}(q)^{T}H^{-1}(q)(J_{N,i}\mp\mu_{i}J_{T,i})}$
(32) $\displaystyle\Lambda_{T,i}$ $\displaystyle=\mp\mu_{i}\Lambda_{N,i}$ (33)

Simultaneous _impacts at several contact points_ can be modeled in several
different ways. For example, iterative solution techniques of the numerical
time-stepping methods of contact dynamics compute the final states of multi-
point impacts by considering sequences of impacts with single points, where
the subsequent impact point are chosen according to a priori rules [18]. In
this work, we adopt the idea that multi-point impacts are equivalent of 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.

###  II-E Zero-order dynamics

In a small neighborhood of the equilibrium state $x_{0}=0$, the dynamics can
be approximated by modified versions of equation (6),(29) 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)$:

$\displaystyle\gamma(q)$ $\displaystyle\approx
J_{\gamma}(0)q;\,\sigma(q)\approx J_{\sigma}(0)q$ (34) $\displaystyle H(q)$
$\displaystyle\approx H(0);\,C(q,v)\approx C(0,0)$ (35) $\displaystyle
J_{N}(q)$ $\displaystyle\approx J_{N}(0);\,J_{T}(q)\approx J_{T}(0)$ (36)

For the previously mentioned example of a planar biped, $H$, $C$ are state-
independent. $J_{N}$, and $J_{T}$ and the functions (21)-(22) are approximated
by:

$\displaystyle
J_{N}\approx\left[\begin{array}[]{cc}-\sin\phi_{1}&-\sin\phi_{2}\\\
\cos\phi_{1}&\cos\phi_{2}\\\ \eta_{1}&\eta_{2}\end{array}\right]$ (40)
$\displaystyle
J_{T}\approx\left[\begin{array}[]{cc}\cos\phi_{1}&\cos\phi_{2}\\\
\sin\phi_{1}&\sin\phi_{2}\\\ \xi_{1}&\xi_{2}\end{array}\right]$ (44)
$\displaystyle\gamma_{i}(q)\approx-y\sin\phi_{i}+z\cos\phi_{i}+\theta\eta_{i}$
(45) $\displaystyle\sigma_{i}(q)\approx
y\cos\phi_{i}+z\sin\phi_{i}-\theta\xi_{i}$ (46)

where we used the previously introduced shorthand notations (23), (24). This
approximation will be referred to as zero-order dynamics or ZOD. The ZOD
delivers linear admissibility conditions, and constant accelerations and
contact forces in each contact mode. On the one hand, the ZOD allows to solve
the equations of motion in all contact modes in closed from, 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$. Nevertheless this approximation brakes down
for systems with marginally consistent contact modes. For example, the ZOD
predicts $\ddot{\gamma}_{2}>\ddot{\gamma}_{1}=0$ for the planar biped in $FF$
mode if $\alpha=\pi/2$, $\phi_{1}=0$, and $-\pi<\phi_{2}<0$. Hence $FF$ is
marginally consistent in the state $q=v=0$. According to the exact equations,
the same system has $\ddot{\gamma}_{1}>0$ in each state with $q=0$ and
$\dot{\theta}\neq 0$ if the geometric parameter $h$ is strictly positive.
Hence the consistency of FF under the exact dynamics in a small neighborhood
of $x=0$ cannot be decided based on the ZOD.

A similar problem may occur in relation with marginally admissible contact
modes. For example,consider the planar biped with $\phi_{1}=\pm\pi/2$,
$l_{1}>l_{2}$, $h=0$ in a state such that point 2 is immobile,
$\dot{\theta}>0$, and $\theta\neq 0$. Contact mode $PS$ is marginally
kinematically admissible under the ZOD but the same mode under the exact
dynamics is not admissible as $\gamma_{1}>0$ .

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.

###  II-F Stability

We will use the standard concept of _Lyapunov stability_. For continuous-time
systems, Lyapunov stability is defined as

###### Definition 2

_Let $x_{0}=(q_{0},0)$ be an equilibrium state. This configuration is called
Lyapunov stable (LS) if for every arbitrarily small $\epsilon>0$ there exists
$\delta(\epsilon)>0$ such that for any kinematically admissible initial state
$x(0)$ that satisfies $|x(0)-x_{0}|<\delta$, all states $x$ along the emerging
motion trajectory satisfy $|x|<\epsilon$._

The trajectories of impacting (hybrid) systems are piecewise smooth, but
discontinuous due to impacts. The definition of Lyapunov stability stated
above is applicable to such systems.

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 the
state space such that $V$ has a local minimumpoint at the equilibrium state
under investigation, and $V$ is non-increasing along motion trajectories. The
adaptation of Lyapunov’s direct method to impacting systems [19] is summarized
below.

Let $\mathcal{B}_{\rho}$ denote a ball of radius $\rho$ in state space
centered around $x_{0}$. Furthermore, let cl denote the closure of a set. For
example, $\textrm{cl}(\mathcal{A})$ is the union of the admissible set and
possible pre-impact states where $\dot{\gamma}_{i}<\gamma_{i}=0$ for some $i$.

###### Theorem 1

_[[19], Thm. 6.23] If there exist_

  1. 1.

a positive value $h$

  2. 2.

a continuous, strictly increasing scalar function $\alpha$ with $\alpha(0)=0$

  3. 3.

a continuously differentiable function
$V:\mathbb{R}^{2n}\rightarrow\mathbb{R}$ such that for any initial state
$x\in\mathcal{B}_{h}\cap\textrm{cl}(\mathcal{A})$

    1. (a)

$V(x)\geq\alpha(|x-x_{0}|)$

    2. (b)

$V$ is non-increasing 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 below, as its logical steps will be
reused during the development of the extended stability theory.

###### Proof:

Consider an arbitrary $\epsilon>0$. Define $c$ as
$c=\alpha(\min(\epsilon,h))$. Then $V(x)\geq c$ along the boundary of
$\mathcal{B}_{\min(h,\epsilon)}$ hence $V(x)$ has a closed, connected sublevel
set $\Omega_{c}=\\{x:V(x)\leq c\\}$ such that $\Omega_{c}$ contains $x_{0}$ in
its interior, furthermore $\Omega_{c}\subset\mathcal{B}_{\min(h,\epsilon)}$.
Due to condition 3b), $\Omega_{c}$ is positively invariant set of the dynamics
(see [19], Proposition 6.5). Then we can find a scalar $\delta>0$ such that
$\mathcal{B}_{\delta}\subset(\Omega_{c}\cup\mathcal{A}^{c})$ where c means
complement. Now if the initial point of a trajectory is in
$\mathcal{B}_{\delta}$, then the positive invariance of $\Omega_{c}$ implies
that the trajectory remains in $\mathcal{B}_{\epsilon}$ at all times, which
proves stability. ∎

###  II-G Semi-definite 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 non-
negative under certain inequality constraints. In general there is no
efficient computational method to test non-negativity of a function or even of
a polynomial function [20].

However for a polynomial function to be non-negative, it is sufficient to
prove that it can be written as the sum of squares of some polynomials, i.e.
it is SOS. Testing if a polynomial is SOS can be formulated as a semi-definite
programming task, for which efficient numerical implementations exist [21].
The theory of SOS polynomials has also been extended to optimization
algorithms over SOS polynomials [22]. 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, 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.

Condition to be satisfied |  Sufficient condition using SOS polynomials
---|---
Find $f(x):\mathbb{R}^{n}\rightarrow\mathbb{R}$ such that $\forall x\in\mathbb{R}^{n}:f(x)\geq 0$ |  Find $f(x)\in\mathcal{P}_{SOS}$
Given $g_{i}(x)\in\mathcal{P}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ ($i=1,2,...,m$), find $f(x):\mathbb{R}^{n}\rightarrow\mathbb{R}$ such that $\forall i:g_{i}(x)=0\implies f(x)\geq 0$ |  Find $h_{i}(x)\in\mathcal{P}$ ($i=1,2,...,m$) such that $f(x)-\sum_{i=1}^{m}h_{i}(x)g_{i}(x)\in\mathcal{P}_{SOS}$
Given $g_{i}(x)\in\mathcal{P}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ ($i=1,2,...,m$), find $f(x):\mathbb{R}^{n}\rightarrow\mathbb{R}$ such that $g(x)\geq 0\implies f(x)\geq 0$ |  Find $h_{i}(x)\in\mathcal{P}_{SOS}$ ($i=1,2,...,m$) such that $f(x)-\sum_{i=1}^{m}h_{i}(x)g_{i}(x)\in\mathcal{P}_{SOS}$
TABLE II: Sufficient conditions of (constrained) non-negativity using SOS
polynomials. $\mathcal{P}$, and $\mathcal{P}_{SOS}$ denote the set of
polynomials, and SOS polynomials.

In a recent paper [13], Posa et al., point 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. Even though the replacement of
non-negativity by the SOS property is conservative, the Lyapunov stability of
several systems with one single contact point was succesfully verified by
[13]. The same method could also provide conservative estimates of the basin
of attraction of a stable equilibrium as well as proof of positive invariance
of some subsets of $\mathcal{S}$, however we are not interested in these
extensions in the present paper. Attempts to test the stability of systems
with multiple contacts have not been reported.

##  III Extended stability theory

In order to test applicability to multi-contact systems, the stability test
proposed by [13] has been implemented using the Mosek solver [23] available
through the Yalmip toolbox of MatLab [24]. We tested the equilibria of the
previously introduced planar biped (Figure 1). Albeit this system has a highly
non-trivial behaviour, almost exact conditions of stability are available if
impacts are inelastic [9, 10]. Little is known about stability in the more
general case of $e>0$, however highly conservative stability conditions have
been reported in [11, 12]. These preliminary tests proved unsuccessful: no
certificates of stability were found by the algorithm, moreover the
corresponding semi-definite programming task was in most cases flagged by the
Mosek solver as provably infeasible. This negative result suggests that the
original stability test has limited ability to verify stability in multi-
contact systems. The failure of the test inspired the development of an
improved stability theory, which can also address these systems.

In a previous work of the author [12], a highly conservative condition of
stability was developed for the planar biped in the case of
$\phi_{1}=\phi_{2}=0$. That work used a Lyapunov-type function composed as the
envelope of 2 smooth functions: $V(x)=\max(V_{1}(x),V_{2}(x))$. 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 FF. For the FF 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 the present paper, these ideas are combined with the
results of [13], 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 semi-definite program. First, an improved stability
theory is developed following the ideas outlined above. Throughout the rest of
the paper the ZOD approximation of the dynamics is considered.

###  III-A The Lyapunov function $V(x)$ may increase temporarily along paths

Lyapunov’s direct method requires that $V(x)$ is a strictly non-increasing
function along solution trajectories and along jumps associated with impacts.
This condition can be replaced by a weaker condition allowing $V(x)$ to
increase temporarily along solution trajectories, provided that it has a non-
increasing overall trend. In particular, we will search for a countable set of
times $t_{i}$, $i=1,2,...$, such that the values of $V(x(t_{i}^{+}))$ form a
strictly non-increasing 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})$. In particular, the discrete set of times $t_{i}$ will be
chosen as times of impact. This choice is motivated by the observation, that
successful Lyapunov functions are often closely related to kinetic energy,
which is known to decrease during every impact.

In order to develop the idea outlined above into a formal statement, one needs
to establish upper bounds of time differences between contact mode
transitions. First we introduce the concept:

###### Definition 3

_A system is called impact-bounded if there exists a strictly increasing
function $\zeta(\xi):\mathbb{R}\rightarrow\mathbb{R}$ with $\zeta(0)=0$ such
that for any contact mode $\mathcal{M}$ and for any state $x\in\mathcal{S}$ in
which mode $\mathcal{M}$ is admissible and consistent, the next contact mode
transition or impact along the trajectory through $x$ occurs no later than at
time $t\leq\zeta(|x|)$._

Then, two useful results are formulated. The first one ensures that contact
mode transitions and impacts occur frequently

###### Lemma 1

_Every system is impact-bounded under the ZOD at an unambiguous equilibrium
$x_{0}=0$._

###### Proof:

Ambiguity means that every non-static contact mode $\mathcal{M}$ is
inconsistent in the equilibrium state $x_{0}$. It is straightforward to see
that for an unambiguous equilibrium, at least one of the following properties
must hold for each mode $\mathcal{M}$ under the ZOD:

  * •

at least one of the admissibility constraints, or one of the consistency
constraints referring to contact forces is violated. In this case,
$\mathcal{M}$ may not be realized at all.

  * •

at least for one contact $i$, the acceleration of the system in mode
$\mathcal{M}$ implies $\dot{\gamma_{i}}<0$.

  * •

at least for one contact $i$, the sign of $\sigma_{i}$ dictated by the
admissibility constraints and the sign of $\dot{\sigma_{i}}$ corresponding to
the acceleration of the system in mode $\mathcal{M}$ satisfy
$\sigma_{i}\dot{\sigma_{i}}<0$

In the first case, there is no need to investgate impact times in mode
$\mathcal{M}$. In the second and third cases, either the normal velocity
$\dot{\gamma_{i}}$ or the tangential velocity $\sigma_{i}$ approaches 0 at a
constant rate. Hence one of the admissibility constraints referring to the
signs $\gamma_{i}$, $\dot{\gamma_{i}}$ or $\sigma_{i}$ is violated after
bounded time. This finding implies the statement of the theorem. The
interested reader may find explicit expression of time bounds for a system
with two contact points in [9]. ∎

The second statement establishes bounds for trajectories completed in short
time intervals:

###### Lemma 2

_There exists a strictly increasing function
$\eta(\xi):\mathbb{R}\rightarrow\mathbb{R}$ with $\eta(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:
$|x(t_{2})|\leq|x(t_{1})|+\eta(t_{2}-t_{1})$ under the ZOD._

###### Proof:

In each contact mode $\mathcal{M}$ the system has a constant acceleration
$\dot{v}_{\mathcal{M}}$ under the ZOD. Let
$a=\max_{\mathcal{M}}|\dot{v}_{\mathcal{M}}|$ denote the largest absolute
value among those accelerations.

By using kinematics of constantly accelerating motion, the following bounds
are obtained:

$\displaystyle|v(t_{2})-v(t_{1})|$ $\displaystyle\leq a(t_{2}-t_{1})$ (47)
$\displaystyle|q(t_{2})-q(t_{1})|$
$\displaystyle\leq|v(t_{1})|(t_{2}-t_{1})+\frac{1}{2}a(t_{2}-t_{1})^{2}$ (48)
$\displaystyle\leq 1\cdot(t_{2}-t_{1})+\frac{1}{2}a(t_{2}-t_{1})^{2}$ (49)

which implies

$|x(t_{2})-x_{0}|\leq|x(t_{1})-x_{0}|+(a+1)(t_{2}-t_{1})+\frac{1}{2}a(t_{2}-t_{1})^{2}.$

proving the statement. ∎

As we will see in Sec. III.C, these results together enable us to prove
stability by constructing a function $V(x)$, which is not monotonic, but whose
values evaluated at times of impacts and contact mode transitions form a
decreasing sequence.

###  III-B Stability certificates via multiple Lyapunov functions

The application of the original stability condition relies on finding a
Lyapunov function $V(x)$ with the properties specified in Theorem 1. Stability
can also be verified if we find a finite set of functions $V_{i}(x)$
($i=1,2,...,n$), none of which satisfies condition 3a) of Theorem 1 provided
that a related condition

$\max_{i}V_{i}(x)\geq\alpha(|x-x_{0}|)$ (50)

is satisfied, furthermore each one of the functions satisfies condition 3b).

This extension may appear somewhat superfluous, because the existence of a set
of function with these properties implies that the non-smooth function defined
by $V(x)=\max_{i}V_{i}(x)$ satisfies the original conditions 3a-b) of the
Theorem. Nevertheless 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. In particular 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 (50).

###  III-C An improved stability condition

Now we are ready to state the main result of the paper:

###### Theorem 2

_Let $x_{0}=(0,0)$ be an unambiguous equilibrium state. If there exists_

  1. 1.

a positive value $h$

  2. 2.

a continuous, strictly increasing scalar functions $\alpha$ with $\alpha(0)=0$

  3. 3.

a finite set of continuously differentiable functions
$V_{i}:\mathbb{R}^{2n}\rightarrow\mathbb{R}$ ($i=1,2,...,\upsilon$) such that
for any initial state $x\in\mathcal{B}_{h}\cap\textrm{cl}(\mathcal{A})$,

    1. (a)

$\max_{i}V_{i}(x)\geq\alpha(|x-x_{0}|)$

    2. (b)

$V(x)\geq V(x^{*})$ where $x^{*}$ denotes the next point of contact mode
transition or the next post-impact state along the 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 $\upsilon=1$, and condition 3) of Theorem 2 is considered in
the limit $t_{2}\searrow t_{1}$.

The proof of Theorem 2 follows the logical steps of Theorem 1. As it has been
pointed out, the relaxed condition 3a) in the statement does not require any
significant change in the proof, as the non-smooth function $\max_{i}V_{i}(x)$
satisfies the requirements of the original theorem. The relaxed condition 3b)
requires some adaptation of the proof.

###### Proof:

Consider an arbitrary $\epsilon>0$. Define $c$ as
$c=\alpha(\min(\epsilon,h))$. Then $V(x)\geq c$ along the boundary of
$\mathcal{B}_{\min(h,\epsilon)}$ hence $V(x)$ has a closed, connected sublevel
set $\Omega_{c}:V(x)\leq c$ such that $\Omega_{c}$ contains $x_{0}$ in its
interior, furthermore $\Omega_{c}\subset\mathcal{B}_{\min(h,\epsilon)}$. Next
we can find a scalar $\beta>0$ such that
$\mathcal{B}_{\beta}\subset(\Omega_{c}\cup\mathcal{A}^{c})$ where c means
complement.

Since $\eta$ and $\zeta$ are strictly monotonic functions, there is a unique
positive value of $\delta$ satisfying the equation

$\delta+\eta(\zeta(\delta)))=\beta$ (51)

If the initial point of a trajectory at time $t_{1}$ is in
$\mathcal{B}_{\delta}$, then by Lemma 1 the next contact mode switch or impact
occurs no later than at time $t_{2}=t_{1}+\zeta(\delta))$. According to Lemma
2, and equation (51), the state of the object remains within the bound given
by $|x(t)-x_{0}|\leq\beta$ during the time interval $t\in(t_{1},t_{2})$, which
is a subset of $\mathcal{B}_{\epsilon}$. By using Condition 3b), we also know
that $|x(t_{2}^{+})|\leq|x(t_{1})|$, i.e. the state $|x(t_{2}^{+})$ is also
inside $\mathcal{B}_{\delta}$ 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 final conclusion that the state of the systems remains within
$\mathcal{B}_{\epsilon}$ at all times, implying Lyapunov stability. ∎

##  IV STABILITY OF A PLANAR BIPED

In this section we apply the extended stability test to a planar model of a
rigid body with 2 contact points, which we will refer to in the sequel as
’biped’.

###  IV-A Implementation of the stability condition

The planar biped has 16 contact modes, which are two-letter words composed of
the set of letters $\\{P,N,F,S\\}$. Among these, exactly 10 are kinematically
admissible in some parts of state space. For example if
$\cos\phi_{1},\cos\phi_{2}>0$. then modes $PN$, $NP$, $PS$, $SP$, $SN$, $NS$
are kinematically inadmissible for all states of the object.The remaining 10
modes 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 the $FF$ contact mode there are no contact forces and the
acceleration is given by the value of $C$, see equation (20).

After the preliminary steps outlined above, the conditions of Theorem 2 are
formulated as an SOS program and tested. Such an implementation requires

  * •

specification of $h$ in condition 1). We use $h=1$.

  * •

Lyapunov function candidates. We choose $V_{1}(x)$ to be total mechanical
energy that is

$V_{1}(x)=\frac{1}{2}(u^{2}+v^{2}+\omega^{2}+z\cos\alpha-x\sin\alpha+T\theta)$
(52)

$V_{2}(x)$ is chosen as

$V_{2}(x)=x+P_{1}^{D}(\phi,z,u,v,\omega)$ (53)

where $P_{1}^{D}$ denotes a polynomial containing all 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=3$ as larger values of $D$ result
in too many unknown coefficients, and high computational cost. In addition,
the nonlinear terms in the variable $x$ are omitted in order to reduce the
number of unknown coefficients. This choice is motivated by the fact that $x$
is a cyclic coordinate of the ZOD dynamics. With the restrictions above
$V_{2}(x)$ has altogether 55 unspecified coefficients to be determined.

  * •

existence of an appropriate function $\alpha(\xi)$ satisfying condition 3a).
In the Appendix, we prove that such a function exists if

$\displaystyle V_{2}(x)$ $\displaystyle\geq x+u^{4}+v^{4}+\omega^{4}$ (54)
$\displaystyle 0<\alpha$ $\displaystyle<\pi/2$ (55) $\displaystyle L_{1}L_{2}$
$\displaystyle<0$ (56)

In what follows, we restrict our attention to systems satisfying (55), (56).
The polynomial inequality (54) is added to the non-negativity conditions to be
fulfilled.

  * •

Implementation of condition 3b) for all possible forms of motion. The function
$V_{1}$ is non-increasing, hence all we need is to prove monotonicity of
$V_{2}$. If $e=0$, then the motion does not have free flight phase. In this
case, we use the following conditions:

    * –

$\frac{d}{dt}V_{2}(x(t))\leq 0$ for all of the non-static contact modes
involving an active contacts $\\{PF,SF,NF,FP,FS,FN,PP,NN\\}$, and for all
states $x\in\mathcal{B}_{h}\cap\textrm{cl}(\mathcal{A})$ where that contact
mode is admissible and consistent

    * –

$V_{2}(M(x(t^{-})))-V_{2}(x(t^{-}))\leq 0$ for all possible pre-impact states
$x(t^{-})\in\mathcal{B}_{h}\cap\textrm{cl}(\mathcal{A})$ where $M$ is the
impact map.

If $e>0$ then free-flight is also possible. Free-flight is always followed by
an impact, and we use the following condition:

    * –

$V(M(x(t^{-})))-V(t-\Delta)\leq 0$ for all possible pre-impact states
$x(t^{-})$ and for all possible durations $\Delta$ of a free-flight phase
immediately before that impact. As before, $M$ means the impact map.

All of these conditions as well as the corresponding conditions of
admissibility and consistency are given in closed form in the Appendix. 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 $\Delta$ in the case of free flight). Hence the whole set of
conditions can be recast as an SOS programming task as demonstrated in
Subsection II-G.

The extended stability test has been implemented in Yalmip optimization
toolbox. In the following section, some test results are presented.

###  IV-B Regions of stability in parameter space

First, the stability test is tested and verified using the example of a biped
on a slope of angle $30^{\circ}$ (that is $\phi_{1}=\phi_{2}=0$,
$\alpha=30^{\circ}$) subject to gravitational load ($T=0$, $F=1$) with
inelastic impacts ($e_{1}=e_{2}=0$). The friction coefficient $\mu_{2}$ is
kept fixed while $\mu_{1}$, and the distance $2d$ between the two contact
points are varied systematically. The positions of the contact points are
chosen as $h=1$, $L_{1}=h\tan\alpha+d$, $L_{2}=h\tan\alpha-d$. Note that this
choice of $L_{1}$, and $L_{2}$ ensures that the object does not topple in $SF$
or $FS$ mode in a static state.

This system has been investigated previously by several works. If
$\mu_{2}>\tan\alpha$ then the object is in frictional equilibrium for all
values of $\mu_{1}$ and $d$. Nevertheless the equilibrium is trivially
unstable due to ambiguity [15] if the following two conditions are satisfied:

$\displaystyle d$ $\displaystyle<h\tan\alpha$ (57) $\displaystyle\mu_{1}$
$\displaystyle<\tan\alpha-d/h$ (58)

In this case, the PF mode is consistent in addition to SS, i.e. a microscopic
perturbation of the equilibrium state may initiate downhill slip motion at
point 1 accompanied by separation of point 2.

Even when $\mu_{1}$ is above the limit (58), the object may be unstable due to
self-excited inverse chatter motion as first pointed out by [12]. The exact
range of parameter values corresponding to instability has been determined by
semi-analytical investigation of an appropriately defined Poincaré map [9,
10]. The newly proposed Lyapunov stability test is able to verify stability in
a significant portion of the stable region of model parameters (Fig. 2). This
test indicates that the new test is conservative, but it is successful in some
highly non-trivial cases.

The exact conditions of Lyapunov stability are available for the previous
example, which possesses two point contacts with ideally inelastic impacts
unknown. In fact, to the author’s best knowledge this is the only known case
of rigid (multi-)body systems with multiple, unilateral frictional contacts,
where conditions of Lyapunov stability are known.

As second example we choose the same model with slope angle
$\alpha=20^{\circ}$, partially elastic impacts $e_{1}=e_{2}=0.1$ and equal
coefficients of friction $\mu_{1}=\mu_{2}$ at the two contact points.

This system is in frictional equilibrium if $\mu_{1},\mu_{2}\geq\tan
20^{\circ}\approx 0.36$ and the equilibrium is unambiguous in all cases. The
exact parameter range implying Lyapunov stability is unknown, however this
example has been investigated analytically by [12] with the aid of a manually
constructed Lyapunov-type function chosen based on physical intuition. That
function could verify stability in a range of model parameters illustrated by
Fig. 3.

We test the stability using the newly proposed SOS program for various
combinations of the parameters $d$ and $\mu_{1}=\mu_{2}$. The test verified
stability in a significantly larger range than the region of stability,
verified by [12].

Figure 2: Equilibrium and stability charts of a planar biped with
$e_{1}=e_{2}=0$, $T=0$, $F=1$, $\alpha=30^{\circ}$, $h=1$,
$\phi_{1}=\phi_{2}=0$, and $L_{1},L_{2}=h\tan\alpha\pm d$, $\mu_{2}=2$ as a
function of $\mu_{1}$ and $d$. Background colors represent analytical results:
unstable equilibrium due to ambiguity (dark), unstable equilibrium due to
inverse chatter (middle) and stable equilibrium (white). Circles indicate
stability region verified by the SOS program. Figure 3: Equilibrium and
stability charts of a planar biped with $e_{1}=e_{2}=0.1$, $T=0$, $F=1$,
$\alpha=20^{\circ}$, $h=1$, $\phi_{1}=\phi_{2}=0$, and
$L_{1},L_{2}=h\tan\alpha\pm d$, as a function of $\mu_{1}=\mu_{2}$ and $d$.
Solid line indicates minimum value of friction coefficient for frictional
equilibrium. White background corresponds to region of stability succesfully
verified by [12]. Circles indicate larger stability region verified by the SOS
program.

##  V CONCLUSIONS

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 initial attempts to predict
stability of some simple model systems, general stability tests 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 [13], 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 Lyapunov stability of simple multi-contact test problems. In
the present work, several improvements of the method have been proposed, which
enabled successful verification of Lyapunov stability in those cases. It has
been found for a particular example that the algorithmic test outcompeted a
previously proposed manual stability test based on physical intuition. At the
same time, the stability test appears to be conservative in those cases where
the exact condition of stability are known.

The author foresees that future developments will radically extend the
applicability of algorithmic stability tests to contact-induced dynamics. A
fundamental challenge to be solved is reduction of high 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.

## APPENDIX

###  V-A Verification of condition 3a) in Theorem 2

The non-penetration constraints under the ZOD (45) now yields
$\gamma_{i}=z_{i}+L_{i}\theta\geq 0$ for $i=1,2$ since $\phi_{i}=0$. This
relation is combined with (56) to obtain the bound

$\displaystyle z\geq\min_{i}|L_{i}|\cdot|\theta|$ (59)

which will be used in the sequel.

Now we separate state space to two regions along the plane
$y=\frac{1}{2}\cot\alpha z$ and the following bounds are established.

If $y\geq\frac{1}{2}\cot\alpha z$ then this inequality along with (54), (59)
implies

$\displaystyle V_{2}(x)$ $\displaystyle\geq y=|y|$ (60)
$\displaystyle\geq\frac{1}{2}\cot\alpha\cdot|z|$ (61)
$\displaystyle\geq\frac{1}{2}\cot\alpha\min_{i}|L_{i}|\cdot|\theta|$ (62)
$\displaystyle V_{2}(x)$ $\displaystyle\geq u^{4},v^{4},\omega^{4}$ (63)

On the other hand, if If $y\leq\frac{1}{2}\cot\alpha z$ then this inequality
and (52), (59) imply

$\displaystyle V_{1}(x)\geq z\cos\alpha-y\sin\alpha$
$\displaystyle\geq\frac{1}{2}\cos\alpha\cdot z$ (64)
$\displaystyle\geq\frac{1}{2}\cos\alpha\min_{i}|L_{i}|\cdot|\theta|$ (65)
$\displaystyle V_{1}(x)\geq z\cos\alpha-y\sin\alpha$
$\displaystyle\geq\sin\alpha\cdot|y|$ (66) $\displaystyle
V_{1}(x)\geq\frac{1}{2}u^{2},\frac{1}{2}v^{2},\frac{1}{2}\omega^{2}$ (67)

In both cases, we have found for each one of the state variables
$\sigma\in\\{{y,z,\theta,u,v,\omega\\}}$ an associated strictly increasing
function $\alpha_{\sigma}$ satisfying $\alpha_{\sigma}(0)=0$ such that
$\max_{i}V_{i}(x)\geq\alpha_{\sigma}(|\sigma|)$. It is straightforward to show
that his property implies condition 3a.

###  V-B State-dependent conditions associated with contact modes

We present state-dependent admissibility and consistency condition of each
contact mode. All of these conditions are formulated as polynomial
inequalities in an appropriate set of variables.

PF, NF, FP, FN modes: Let $i$ denote the active contact, and $j=3-i$ the
inactive one. The normal contact force is determined from the consistency
constraint $\ddot{\gamma}_{i}=0$ using the equation of motion (6) and the
kinematic relations given in Sec. II-E:

$\lambda_{\pm,i}=\frac{\cos\alpha}{1+L_{i}^{2}\mp\mu_{i}L_{i}H}$

where the $\pm$ signs should be understood as $+$ for slip in the positive
direction and $-$ for negative slip. If $\lambda_{\pm,i}<0$ then the
corresponding mode is not consistent in any state, hence the condition given
below can be disregarded. In the opposite case, we express the time derivative
of the vector $x$ of state veriables as

$\displaystyle\dot{x}_{\pm,i}$
$\displaystyle=(u,v,\omega,\dot{u},\dot{v},\dot{\omega})$ (68)
$\displaystyle=(u,v,\omega,\sin\alpha\mp\mu_{i}\lambda_{\pm,i},-\cos\alpha+\lambda_{\pm,i},(L_{i}\mp\mu_{i}H)\lambda_{\pm,i})$
(69)

Finally, condition 3b) is given by

$\displaystyle-J_{V}(x)\dot{x}_{\pm,i}\geq 0$ whenever
$\displaystyle\gamma_{i}(q)=0$ (70) $\displaystyle\dot{\gamma}_{i}(q)=0$ (71)
$\displaystyle\gamma_{j}(q)\geq 0$ (72)
$\displaystyle\dot{\sigma}_{i}(q,v)\begin{array}[]{c}\geq\\\ \leq\end{array}0$
(75) $\displaystyle 1-x^{T}x\geq 0$ (76)

where $J_{V}(x)$ is the Jacobian of the function $V(x)$, which can be
expressed in terms of the 55 unknown coefficients of $V(x)$.

SF, FS modes: Let $\lambda$ and $\lambda_{t}$ denote normal and tangential
contact forces at the sticking contact $i$, and let $j$ denote the other
contact. Then, the contact forces can be determined using the kinematic
constraints $\ddot{\gamma}_{i}=\dot{\sigma}_{i}=0$:

$[\lambda,\lambda_{t}]^{T}=-P[-\cos\alpha,\sin\alpha]^{T}$

with

$P=\left[\begin{array}[]{cc}1+L_{i}^{2}&L_{i}H\\\
L_{i}H&1+H^{2}\end{array}\right]$

If $\lambda<0$ or $|\lambda_{t}|>\mu_{i}\lambda$, then the contact mode under
investigation can be omited otherwise the time derivative of $x$ is

$\displaystyle\dot{x}_{S,i}$
$\displaystyle=(u,v,\omega,\sin\alpha+\lambda_{t},-\cos\alpha+\lambda,L_{i}\lambda+H\lambda_{t}$
(77)

and the stability conditions are formulated as

$\displaystyle-J_{V}(x)\dot{x}_{S,i}\geq 0$ whenever
$\displaystyle\gamma_{i}(q)=0$ (78) $\displaystyle\dot{\gamma}_{i}(q)=0$ (79)
$\displaystyle\gamma_{j}(q)\geq 0$ (80) $\displaystyle 1-x^{T}x\geq 0$ (81)

We note that $\dot{x}_{S,i}$ is a convex combination of $\dot{x}_{+,i}$ and
$\dot{x}_{-,i}$, i.e.
$\dot{x}_{S,i}=\kappa\dot{x}_{+,i}+(1-\kappa)\dot{x}_{-,i}$ for some
$0\leq\kappa\leq 1$. Hence the conditions (69) for positive and negative slip
together imply (81).

PP, NN modes: In this case, the two normal contact forces are determined by
the constraints $\ddot{\gamma}_{1}=\ddot{\gamma}_{2}=0$:

$[\lambda_{1},\lambda_{2}]^{T}=-P^{-1}[\cos\alpha,\cos\alpha]^{T}$

where

$P=\left[\begin{array}[]{cc}1+L_{1}^{2}\mp\mu_{1}HL_{1}&1+L_{1}L_{2}\mp\mu_{2}HL_{1}\\\
1+L_{1}L_{2}\mp\mu_{1}HL_{2}&1+L_{2}^{2}\mp\mu_{2}HL_{2}\end{array}\right]$

If any of the two contact force values is negative, than the conditions
associated with the contact mode under investigation can be disregarded. In
the opposite case, we have

$\displaystyle\dot{x}_{\pm}$
$\displaystyle=(u,v,\omega,\sin\alpha\mp\sum_{i}\mu_{i}\lambda_{i},-\cos\alpha+\sum_{i}\lambda_{i},\sum_{i}(L_{i}\mp\mu_{i}H)\lambda_{i})$
(82)

and the stability conditions are formulated as

$\displaystyle-J_{V}(x)\dot{x}_{\pm}\geq 0$ whenever
$\displaystyle\gamma_{1}(q)=\gamma_{2}(q)=0$ (83)
$\displaystyle\dot{\gamma}_{i}(q)=\dot{\gamma}_{2}=0$ (84)
$\displaystyle\sigma_{1}(q,v)\begin{array}[]{c}\geq\\\ \leq\end{array}0$ (87)
$\displaystyle 1-x^{T}x\geq 0$ (88)

Impacts: The impact maps of sticking and slipping impacts have been given in
closed form in the main text. Instead of treating the three types of impacts
separately, slipping and sticking impacts are addressed now in a common
framework as follows. Let $i$ be the contact point undergoing an impact. Let
$\Lambda_{N,i}$ and $\Lambda_{T,i}$ denote net impulses in the normal and
tangential directions during that impact.

Then $\Lambda_{N,i}$ can be expressed in terms of the pre-impact state and of
$\Lambda_{T,i}$ by using the condition
$\dot{\gamma}_{i}(x(t^{+}))=-e\dot{\gamma}_{i}(x(t^{-}))$ and (29):

$\Lambda_{N,i}=\frac{-(1+e)(v_{z}(t^{-})+l_{i}\omega(t^{-}))+\Lambda_{T,i}hl_{i}}{1+l_{i}^{2}}$

Then, the post-impact state becomes

$\displaystyle x^{+}$
$\displaystyle=x^{-}+(0,0,0,\Lambda_{T,i},\Lambda_{N,i},\Lambda_{T,i}H+\Lambda_{N,i}L_{i})$
(89)

and the stability conditions are formulated as

$\displaystyle V(x^{-})-V(x^{+})\geq 0$ whenever
$\displaystyle\gamma_{i}(q)=0$ (90) $\displaystyle\dot{\gamma}_{i}(x^{-})\leq
0$ (91) $\displaystyle\gamma_{j}(q)\geq 0$ (92)
$\displaystyle(u^{-}+\omega^{+}H)\lambda_{t}\leq 0$ (93)
$\displaystyle\mu_{i}^{2}\lambda^{2}-\lambda_{t}^{2}\geq 0$ (94)
$\displaystyle(\mu_{i}^{2}\lambda^{2}-\lambda_{t}^{2})$
$\displaystyle\cdot(u^{-}+\omega^{+}H)=0$ (95)

The last three conditions express Coulomb’s law for the impact. The conditions
listed above are polynomial in the set of variables
$\\{y,z,\theta,u,v,\omega,\lambda_{t}\\}$

FF mode followed by an impact:

## References

  * [1] B. Brogliato and B. Brogliato, Nonsmooth mechanics, vol. 3.  Springer, 1999.
  * [2] R. I. Leine, “The historical development of classical stability concepts: Lagrange, poisson and lyapunov stability,” Nonlinear Dynamics, vol. 59, no. 1, pp. 173–182, 2010.
  * [3] R. Leine and N. van de Wouw, “Stability properties of equilibrium sets of non-linear mechanical systems with dry friction and impact,” Nonlinear Dynamics, vol. 51, no. 4, pp. 551–583, 2008.
  * [4] J. Zhang, K. H. Johansson, J. Lygeros, and S. Sastry, “Zeno hybrid systems,” International Journal of Robust and Nonlinear Control: IFAC-Affiliated Journal, vol. 11, no. 5, pp. 435–451, 2001.
  * [5] Y. Or and A. D. Ames, “Stability and completion of zeno equilibria in lagrangian hybrid systems,” IEEE Transactions on Automatic Control, vol. 56, no. 6, pp. 1322–1336, 2010.
  * [6] A. R. Champneys and P. L. Várkonyi, “The Painlevé paradox in contact mechanics,” IMA Journal of Applied Mathematics, vol. 81, no. 3, pp. 538–588, 2016.
  * [7] A. M. Liapunov, Stability of motion.  Elsevier, 2016.
  * [8] Y. V. Orlov, Discontinuous systems: Lyapunov analysis and robust synthesis under uncertainty conditions.  Springer Science & Business Media, 2008.
  * [9] P. L. Várkonyi and Y. Or, “Lyapunov stability of a rigid body with two frictional contacts,” Nonlinear dynamics, vol. 88, no. 1, pp. 363–393, 2017.
  * [10] Y. Or and P. L. Várkonyi, “Experimental verification of stability theory for a planar rigid body with two unilateral frictional contacts,” IEEE Transactions on Robotics, vol. 37, no. 5, pp. 1634–1648, 2021.
  * [11] Y. Or and E. Rimon, “On the hybrid dynamics of planar mechanisms supported by frictional contacts. ii: Stability of two-contact rigid body postures,” in 2008 IEEE International Conference on Robotics and Automation, pp. 1219–1224, IEEE, 2008.
  * [12] P. L. Várkonyi, D. Gontier, and J. W. Burdick, “On the Lyapunov stability of quasistatic planar biped robots,” in 2012 IEEE International Conference on Robotics and Automation, pp. 63–70, IEEE, 2012.
  * [13] M. Posa, M. Tobenkin, and R. Tedrake, “Stability analysis and control of rigid-body systems with impacts and friction,” IEEE Transactions on Automatic Control, vol. 61, no. 6, pp. 1423–1437, 2015.
  * [14] D. E. Stewart, “Rigid-body dynamics with friction and impact,” SIAM review, vol. 42, no. 1, pp. 3–39, 2000.
  * [15] Y. Or and E. Rimon, “On the hybrid dynamics of planar mechanisms supported by frictional contacts. i: Necessary conditions for stability,” in 2008 IEEE International Conference on Robotics and Automation, pp. 1213–1218, IEEE, 2008.
  * [16] T. R. Kane and D. A. Levinson, Dynamics, theory and applications.  McGraw Hill, 1985.
  * [17] T. Kane and D. Levinson, “An explicit solution of the general two-body collision problem,” Computational mechanics, vol. 2, no. 1, pp. 75–87, 1987.
  * [18] V. Acary and B. Brogliato, Numerical methods for nonsmooth dynamical systems: applications in mechanics and electronics.  Springer Science & Business Media, 2008.
  * [19] R. I. Leine and N. Van de Wouw, Stability and convergence of mechanical systems with unilateral constraints, vol. 36.  Springer Science & Business Media, 2007.
  * [20] K. G. Murty and S. N. Kabadi, “Some np-complete problems in quadratic and nonlinear programming,” tech. rep., 1985.
  * [21] M. Marshall, Positive polynomials and sums of squares.  No. 146, American Mathematical Soc., 2008.
  * [22] P. A. Parrilo, “Semidefinite programming relaxations for semialgebraic problems,” Mathematical programming, vol. 96, no. 2, pp. 293–320, 2003\.
  * [23] E. D. Andersen and K. D. Andersen, “The mosek interior point optimizer for linear programming: an implementation of the homogeneous algorithm,” in High performance optimization, pp. 197–232, Springer, 2000.
  * [24] J. Lofberg, “Yalmip: A toolbox for modeling and optimization in matlab,” in 2004 IEEE international conference on robotics and automation (IEEE Cat. No. 04CH37508), pp. 284–289, IEEE, 2004.


$FH^{1}(Y,\mathcal{O}_{Y})$.

###  10.2. Computing Frobenius on $H^{1}(X,\mathcal{O}_{X})$ (review)

Fix an effective, non-special divisor $D$ on $X$. (The most “efficient” choice
is to take $D$ of degree $g_{X}$, but this may not be possible over a small
ground field.) Thus we have $H^{1}(X,\mathcal{O}_{X}(D))=0$, and the exact
sequence

$0\to\mathcal{O}_{X}\to\mathcal{O}_{X}(D)\to\mathcal{O}_{X}(D)/\mathcal{O}_{X}\to
0$

is an acyclic resolution of $\mathcal{O}_{X}$. Taking cohomology yields an
isomorphism

$\frac{H^{0}(X,\mathcal{O}_{X}(D)/\mathcal{O}_{X})}{H^{0}(X,\mathcal{O}_{X}(D))}\cong
H^{1}(X,\mathcal{O}_{X}).$

The divisor $pD$ is also non-special, and we find an isomorphism

$\frac{H^{0}(X,\mathcal{O}_{X}(pD)/\mathcal{O}_{X})}{H^{0}(X,\mathcal{O}_{X}(pD))}\cong
H^{1}(X,\mathcal{O}_{X}).$

The composed isomorphism

$\frac{H^{0}(X,\mathcal{O}_{X}(pD)/\mathcal{O}_{X})}{H^{0}(X,\mathcal{O}_{X}(pD))}\cong\frac{H^{0}(X,\mathcal{O}_{X}(D)/\mathcal{O}_{X})}{H^{0}(X,\mathcal{O}_{X}(D))}$

can be computed explicitly by taking a meromorphic function $t$ on a
neighborhood of $D$ representing an element of
$H^{0}(X,\mathcal{O}_{X}(pD)/\mathcal{O}_{X})$ and “correcting” it by the
principal part of an element of $H^{0}(X,\mathcal{O}_{X}(pD))$ so that it has
poles no worse than $D$, and therefore defines an element of
$H^{0}(X,\mathcal{O}_{X}(D)/\mathcal{O}_{X})$. The result is well-defined up
to an element of $H^{0}(X,\mathcal{O}_{X}(D))$.

The Frobenius endomorphism

$F:H^{1}(X,\mathcal{O}_{X})\to H^{1}(X,\mathcal{O}_{X})$

can then be computed as the composition

$H^{1}(X,\mathcal{O}_{X})\cong\frac{H^{0}(X,\mathcal{O}_{X}(D)/\mathcal{O}_{X})}{H^{0}(X,\mathcal{O}_{X}(D))}\to\frac{H^{0}(X,\mathcal{O}_{X}(pD)/\mathcal{O}_{X})}{H^{0}(X,\mathcal{O}_{X}(pD))}\\\
\cong H^{0}(X,\mathcal{O}_{X}(D)/\mathcal{O}_{X})\cong
H^{1}(X,\mathcal{O}_{X})$

where the second homomorphism is induced by $t\mapsto t^{p}$.

Summarizing, to compute $F$ on $H^{1}(X,\mathcal{O}_{X})$, it suffices to know
the principal parts of elements of $H^{0}(X,\mathcal{O}_{X}(pD))$. It will
transpire later in this section that to compute Frobenius on
$H^{1}(Y,\mathcal{O}_{Y})$, it will suffice to know the principal parts of
elements of $H^{0}(X,\mathcal{O}_{X}((2p-1)D))$ plus some simple linear
algebra.

###  10.3. Making $\pi_{*}\mathcal{O}_{Y}$ explicit

We make the description of $\mathcal{F}=\pi_{*}\mathcal{O}_{Y}$ given in the
previous section more explicit.

Let $D$ be an effective, non-special divisor on $X$. Let $U=X\setminus D$ and
choose an affine open neighborhood $V$ of $D$. Then $\\{U,V\\}$ is a cover of
$X$ by affine opens. Choose elements $f,g\in k(X)$ with $f$ regular on $U$,
$g$ regular on $V$, such that

$k(Y)\cong\frac{k(X)[y]}{(\wp(y)-f)}\cong\frac{k(X)[z]}{(\wp(z)-g)}$

and such that $f-g=\wp(h)$ where at each closed point of $D$, $h$ has poles no
worse than $D$.

Over $U$, the sections $1,y,\dots,y^{p-1}$ are a basis of $\mathcal{F}$ as
$\mathcal{O}_{X}$-module. Similarly, over $V$, the sections
$1,z,\dots,z^{p-1}$ are a basis. Over $U\cap V$, we compute a transition
matrix $H$ as follows:

$\left(\alpha_{0},\dots,\alpha_{p-1}\right)\begin{pmatrix}1\\\ y\\\ \vdots\\\
y^{p-1}\end{pmatrix}=\left(\alpha_{0},\dots,\alpha_{p-1}\right)\begin{pmatrix}1\\\
z+h\\\ \vdots\\\
(z+h)^{p-1}\end{pmatrix}=\left(\alpha_{0},\dots,\alpha_{p-1}\right)H\begin{pmatrix}1\\\
z\\\ \vdots\\\ z^{p-1}\end{pmatrix}$

where

$H=\begin{pmatrix}1&0&0&\dots&0\\\ h&1&0&\dots&0\\\ h^{2}&2h&1&\dots&0\\\
\vdots&\vdots&\vdots&\ddots&0\\\
h^{p-1}&(p-1)h^{p-2}&\binom{p-1}{2}h^{p-3}&\dots&1\end{pmatrix}.$ (10.1)

Numbering the rows and columns of $H$ from $0$ to $p-1$, the $(i,j)$ entry of
$H$ is $\binom{i}{j}h^{i-j}$.

Thus, if $W\subset X$ is an open subset, a section of $\mathcal{F}$ over $W$
is determined by a tuple of functions
$\left(\alpha_{0},\dots,\alpha_{p-1}\right)$ such that each $\alpha_{i}$ is
regular on $U\cap W$ and the tuple

$\left(\beta_{0},\dots,\beta_{p-1}\right):=\left(\alpha_{0},\dots,\alpha_{p-1}\right)H$

has each $\beta_{i}$ regular on $V\cap W$.

The $\mathcal{O}_{X}$-module $\mathcal{F}$ is self-dual. Indeed, letting $A$
be the $p\times p$ anti-diagonal matrix with non-zero entries equal to $\pm
1$:

$A=\begin{pmatrix}0&0&\dots&0&1\\\ 0&0&\dots&-1&0\\\
\vdots&\vdots&\reflectbox{$\ddots$}&\vdots&\vdots\\\ 0&-1&\dots&0&0\\\
1&0&\dots&0&0\end{pmatrix},$

one computes that ${}^{t}H^{-1}=AHA^{-1}=AHA$. The induced pairing

$\mathcal{F}\otimes_{\mathcal{O}_{X}}\mathcal{F}\cong\mathcal{F}\otimes_{\mathcal{O}_{X}}\mathcal{F}^{\vee}\to\mathcal{O}_{X}$

is determined by

$\langle y^{i},y^{j}\rangle=(-1)^{i}\delta_{p-1,i+j}$

and satisfies $\langle z^{i},z^{j}\rangle=(-1)^{i}\delta_{p-1,i+j}$ as well.

The $\mathcal{O}_{X}$-module $\mathcal{F}$ also carries an action of
$\operatorname{Gal}(Y/X)$. Let

$\gamma=\begin{pmatrix}1&0&0&\dots&0\\\ 1&1&0&\dots&0\\\ 1&2&1&\dots&0\\\
\vdots&\vdots&\vdots&\ddots&0\\\ 1&(p-1)&\binom{p-1}{2}1&\dots&1\end{pmatrix}$

be the matrix whose $(i,j)$ entry is $\binom{i}{j}$ (numbering the rows and
columns from 0 to $p-1$). Then the element of $\operatorname{Gal}(Y/X)$
corresponding to $1\in\mathbb{Z}/p\mathbb{Z}$ acts as

$\left(\alpha_{0},\dots,\alpha_{p-1}\right)\mapsto\left(\alpha_{0},\dots,\alpha_{p-1}\right)\gamma\left(\beta_{0},\dots,\beta_{p-1}\right)\mapsto\left(\beta_{0},\dots,\beta_{p-1}\right)\gamma.$

(Note that $g$ and $H$ commute.)

###  10.4. Computing $H^{1}(Y,\mathcal{O}_{Y})$

We keep the notation of the preceding subsections. Since $\pi$ is finite, we
have

$H^{1}(Y,\mathcal{O}_{Y})\cong
H^{1}(X,\pi_{*}\mathcal{O}_{Y})=H^{1}(X,\mathcal{F}).$

We compute the latter using an acyclic resolution.

Since the $\mathcal{O}_{X}$-module $\mathcal{F}$ is a repeated extension of
copies of $\mathcal{O}_{X}$, an easy inductive argument shows that
$H^{1}(X,\mathcal{F}(D))=0$ and
$h^{0}(X,\mathcal{F}(D))=ph^{0}(X,\mathcal{O}_{X}(D))$. Thus, the exact
sequence

$0\to\mathcal{F}\to\mathcal{F}(D)\to\mathcal{F}(D)/\mathcal{F}\to 0$

is an acyclic resolution of $\mathcal{F}$, and we have an isomorphism

$\frac{H^{0}(X,\mathcal{F}(D)/\mathcal{F})}{H^{0}(X,\mathcal{F}(D))}\cong
H^{1}(X,\mathcal{F}).$

We make the left hand side more explicit. The description of $\mathcal{F}$
over the neighborhood $V$ of $D$ in the last section yields an identification

$H^{0}(X,\mathcal{F}(D)/\mathcal{F})\cong
H^{0}(X,\mathcal{O}_{X}(D)/\mathcal{O}_{X})^{p}$

which sends the class of $\left(\beta_{0},\dots,\beta_{p-1}\right)$ on the
right to the section

$\left(\beta_{0},\dots,\beta_{p-1}\right)\begin{pmatrix}1\\\ z\\\ \vdots\\\
z^{p-1}\end{pmatrix}$

of $\mathcal{F}(D)/\mathcal{F}$. To obtain $H^{1}(X,\mathcal{F})$, we need to
take the quotient by $H^{0}(X,\mathcal{F}(D))$. The section $s_{0}:=1\in
H^{0}(X,\mathcal{F}(D))$ maps to zero in
$H^{0}(X,\mathcal{O}_{X}(D)/\mathcal{O}_{X})^{p}$, and the section $s_{1}:=y$
maps to

$\left(\beta_{0},\dots,\beta_{p-1}\right)=(h,1,0,\dots,0)=(h,0,\dots,0).$

Other sections of $H^{0}(X,\mathcal{F}(D))$ are somewhat less explicit, but
can be constructed as follows. For $\ell=2,\dots,p-1$, consider a section
$s_{\ell}$ of $\mathcal{F}(D)$ corresponding to tuple
$\left(\alpha_{0},\dots,\alpha_{p-1}\right)$ with $\alpha_{j}=0$ for $j>\ell$,
$\alpha_{\ell}=1$, $\alpha_{\ell-1}=0$, and $\alpha_{j}$ with $j\leq\ell-2$ to
be chosen. In order to obtain a section of $\mathcal{F}(D)$ the tuple

$\left(\beta_{0},\dots,\beta_{p-1}\right)=\left(\alpha_{0},\dots,\alpha_{p-1}\right)H$

should have entries with poles no worse than $D$ on $V$. This is automatic for
$\beta_{j}$ with $j>\ell$ (since $\beta_{j}=0$) as well as for $j=\ell$ (since
$\beta_{\ell}=1$) and $j=\ell-1$ (since $\beta_{\ell-1}=\ell h$). There is a
function $\alpha_{\ell-2}\in H^{0}(X,\mathcal{O}_{X}(2D))$ (unique up to
addition of a scalar) such that

$\beta_{\ell-2}=\alpha_{\ell-2}+\binom{\ell}{2}h^{2}\alpha_{\ell}$

has poles no worse that $D$ on $V$. We continue to choose $\alpha_{j}\in
H^{0}(X,\mathcal{O}_{X}((\ell-j)D))$ with $j<\ell-2$ in descending order to
satisfy the condition that $\beta_{j}$ have poles no worse than $D$ on $V$,
thus obtaining a section $s_{\ell}$ in $H^{0}(X,\mathcal{F}(D))$.

It remains to consider the images of the $s_{j}$ in
$H^{0}(X,\mathcal{F}(D)/\mathcal{F})$. Inspection of the process for choosing
the $\alpha_{j}$ in the last paragraph then shows that

$\displaystyle s_{0}$
$\displaystyle\mapsto\left(1,0,\dots,0\right)=\left(0,0,\dots,0\right)$
$\displaystyle s_{1}$
$\displaystyle\mapsto\left(h,1,0,\dots,0\right)=\left(h,0,\dots,0\right)$
$\displaystyle s_{2}$
$\displaystyle\mapsto\left(*,2h,1,0\dots,0\right)=\left(*,2h,0,\dots,0\right)$
$\displaystyle\vdots$ $\displaystyle s_{\ell}$
$\displaystyle\mapsto\left(*,\dots,*,\ell h,0,\dots,0\right)$

where the last non-zero entry in the image of $s_{\ell}$ occurs in column
$\ell-1$ (numbering from 0).333We know nothing about the entries $*$ other
than that they have poles no worse than $D$ on $V$. They can of course be
computed for a given $Y$, but it is not clear how to say anything explicit
about them in general, except for $s_{2}$.

###  10.5. Computing Frobenius on $H^{1}(Y,\mathcal{O}_{Y})$

With notation as before, we have isomorphisms

$H^{1}(Y,\mathcal{O}_{Y})\cong
H^{1}(X,\mathcal{F})\cong\frac{H^{0}(X,\mathcal{F}(D)/\mathcal{F})}{H^{0}(X,\mathcal{F}(D))},$

and the Frobenius endomorphism of $H^{1}(Y,\mathcal{O}_{Y})$ transported to
the right hand group above is given by

$\frac{H^{0}(X,\mathcal{F}(D)/\mathcal{F})}{H^{0}(X,\mathcal{F}(D))}\to\frac{H^{0}(X,\mathcal{F}(pD)/\mathcal{F})}{H^{0}(X,\mathcal{F}(pD))}\cong\frac{H^{0}(X,\mathcal{F}(D)/\mathcal{F})}{H^{0}(X,\mathcal{F}(D))}$

where the first arrow is $s\mapsto s^{p}$ and the second is the isomorphism
obtained by noting that

$0\to\mathcal{F}\to\mathcal{F}(pD)\to\mathcal{F}(pD)/\mathcal{F}\to 0$

is another acyclic resolution of $\mathcal{F}$. Note that making this explicit
requires understanding the principal parts along $D$ of sections of
$\mathcal{F}(pD)$.

We want to further transport Frobenius via the isomorphism

$H^{0}(X,\mathcal{F}(D)/\mathcal{F})\cong
H^{0}(X,\mathcal{O}_{X}(D)/\mathcal{O}_{X})^{p}$

as in the previous section. Note that if $s=\sum\beta_{i}z^{i}$, then

$s^{p}=\sum\beta_{i}^{p}z^{pi}=\sum\beta_{i}^{p}(z+g)^{i}=\left(\beta_{0}^{p},\dots,\beta_{p-1}^{p}\right)G\begin{pmatrix}1\\\
z\\\ z^{2}\\\ \vdots\\\ z^{p-1}\end{pmatrix}$

where

$G=\begin{pmatrix}1&0&0&\dots&0\\\ g&1&0&\dots&0\\\ g^{2}&2g&1&\dots&0\\\
\vdots&\vdots&\vdots&\ddots&0\\\
g^{p-1}&(p-1)g^{p-2}&\binom{p-1}{2}g^{p-3}&\dots&1\end{pmatrix}$

(Numbering the rows and columns of $G$ from $0$ to $p-1$, the $(i,j)$ entry of
$G$ is $\binom{i}{j}g^{i-j}$.)

Thus to compute the image of the class of
$\left(\beta_{0},\dots,\beta_{p-1}\right)$, we should form

$\left(\beta_{0}^{p},\dots,\beta_{p-1}^{p}\right)G\in
H^{0}(V,\mathcal{O}_{X}(pD))^{p}$

and then use sections of $H^{0}(X,\mathcal{F}(pD))$ to “reduce” this quantity
so that it lies in $H^{0}(V,\mathcal{O}_{X}(D))$.

We end by examining what needs to be known to write down global sections of
$\mathcal{F}(pD)$: They are given by tuples

$\left(\alpha_{0},\dots,\alpha_{p-1}\right)\in\mathcal{O}_{X}(U)^{p}$

such that

$\left(\beta_{0},\dots,\beta_{p-1}\right)=\left(\alpha_{0},\dots,\alpha_{p-1}\right)H$

has components lying in $\mathcal{O}_{X}(pD)(V)$. To find them, choose
$\alpha_{p-1}\in H^{0}(X,\mathcal{O}_{X}(pD))$ arbitrarily. Then choose
$\alpha_{p-2}\in H^{0}(X,\mathcal{O}_{X}((p+1)D)$ such that
$\alpha_{p-2}+(p-1)h\alpha_{p-1}$ has poles no worse than $pD$ in $V$. (The
set of such choices is a homogeneous space for
$H^{0}(X,\mathcal{O}_{X}(pD))$.) Iterating, one sees that $\alpha_{i}$ lies in
$H^{0}(X,\mathcal{O}_{X}((2p-1-i)D)$ and is uniquely determined up to addition
of an element of $H^{0}(X,\mathcal{O}_{X}(pD))$. Thus, if we have good control
on the principal parts of elements of $H^{0}(X,\mathcal{O}_{X}((2p-1)D))$, we
can compute Frobenius on $H^{1}(Y,\mathcal{O}_{Y})$.

###  10.6. Completing the Hasse-Witt triple of $Y$

Since $\pi$ is finite and étale, we have
$\pi^{*}\Omega^{1}_{X}\cong\Omega^{1}_{Y}$ and

$\pi_{*}\Omega^{1}_{Y}\cong\Omega^{1}_{X}\otimes_{\mathcal{O}_{X}}\pi_{*}\mathcal{O}_{X}=\Omega^{1}_{X}\otimes_{\mathcal{O}_{X}}\mathcal{F}.$

We write $\mathcal{F}^{1}$ for
$\Omega^{1}_{X}\otimes_{\mathcal{O}_{X}}\mathcal{F}$. The auto-duality of
$\mathcal{F}$ induces a bilinear map of $\mathcal{O}_{X}$-modules

$\mathcal{F}\otimes_{\mathcal{O}_{X}}\mathcal{F}^{1}\to\Omega^{1}_{X}.$

Let $Q=H^{1}(Y,\mathcal{O}_{Y})\cong H^{1}(X,\mathcal{F})$ and let $\Phi$ be
the Frobenius endomorphism as computed in the preceding section. Then Serre
duality says that $Q^{\vee}\cong H^{0}(Y,\Omega^{1}_{Y})$ and since $\pi$ is
finite, $H^{0}(Y,\Omega^{1}_{Y})\cong H^{0}(X,\mathcal{F}^{1})$.

Recall that we have fixed an effective, non-special divisor $D$. We have

$H^{0}(X,\Omega^{1}_{X}(-D))\cong H^{1}(X,\mathcal{O}_{X}(D))^{\vee}=0,$

and an easy argument by induction shows that $H^{0}(X,\mathcal{F}^{1}(-D))=0$
as well. We find an injection444Note that the right hand side here is an
explicit $pg$-dimensional vector space, and we can compute in it, rather than
worrying about writing down an explicit basis of $H^{0}(X,\mathcal{F}^{1})$

$H^{0}(\mathcal{F}^{1})\to H^{0}(X,\mathcal{F}^{1}/\mathcal{F}^{1}(-D)).$

Its image is easily seen to be orthogonal to

$\operatorname{Im}\left(H^{0}(X,\mathcal{F}(D))\to
H^{0}(X,\mathcal{F}(D)/\mathcal{F})\right)=\operatorname{Ker}\left(H^{0}(X,\mathcal{F}(D)/\mathcal{F})\twoheadrightarrow
H^{1}(X,\mathcal{F})\right)$

under the pairing

$H^{0}(X,\mathcal{F}(D)/\mathcal{F})\times
H^{0}(X,\mathcal{F}^{1}/\mathcal{F}^{1}(-D))\to
H^{0}(X,\Omega^{1}_{X}(D)/\Omega^{1}_{X})\overset{\operatorname{Res}}{\longrightarrow}k$

where the first arrow is induced by the bilinear map mentioned above, and the
second arrow is the sum of residues along $D$.

Now suppose $s\in H^{0}(X,\mathcal{F}(D)/\mathcal{F})$ maps to an element in
$H^{1}(X,\mathcal{F})$ which is killed by Frobenius. This means that there is
a global section $t\in H^{0}(X,\mathcal{F}(pD))\cong
H^{0}(Y,\mathcal{O}_{Y}(\pi^{*}D))$ whose principal parts along $\pi^{*}D$ are
given by $s^{p}$. It is then immediate that $\omega=dt$ is a _regular_ 1-form
on $Y$, i.e., an element of $H^{0}(Y,\Omega^{1}_{Y})\cong
H^{0}(X,\mathcal{F}^{1})$. The map $\Psi$ is then given by

$\Psi:\operatorname{Ker}(\Phi)\to Q^{\vee}\qquad[s]\mapsto\omega.$ (10.2)

Since $\omega$ is exact, it is orthogonal to $\operatorname{Im}(\Phi)$, so our
$\Psi$ has the required properties.

Summing up, we have proven:

######  Proposition 10.7.

The Hasse-Witt triple associated to $H^{1}_{dR}(Y)$ is $(Q,\Phi,\Psi)$ with
$Q=H^{1}(X,\mathcal{F})$ computed explicitly in Section 10.4, with Frobenius
$\Phi$ defined in Section 10.5, and with $\Psi$ defined in the paragraph
before equation (10.2).

## References