# Mechanical Systems of Rigid Bodies Subject to Unilateral Constraints

Per Lötstedt

► To cite this version:

Per Lötstedt. Mechanical Systems of Rigid Bodies Subject to Unilateral Constraints. SIAM Journal on Applied Mathematics, 1982, 42 (2), pp.281-296. 10.1137/0142022. hal-01399632

HAL Id: hal-01399632

https://hal.science/hal-01399632

Submitted on 20 Nov 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# MECHANICAL SYSTEMS OF RIGID BODIES SUBJECT TO UNILATERAL CONSTRAINTS

PER LÖTSTEDT

**Abstract.** The properties of mechanical systems of rigid bodies subject to unilateral constraints are investigated. In particular, properties of interest for the digital simulation of the motion of such systems are studied. The constraints give rise to discontinuities in the solution. Under general assumptions on the system a unique solution is constructed using the linear complementarity theory of mathematical programming. A numerical method for solution of these problems and generalizations of the constraints studied in this paper are briefly discussed.

**1. Introduction.** The motion of a system of rigid bodies subject to constraints is of interest in many applications, e.g., biomechanics, King and Chou (1976), the performance of vehicles, Magnus (1978), analysis of mechanisms and machines, Paul (1975). It is of practical importance to know how to simulate such systems on digital computers. Sometimes the number of kinematical constraints satisfied as equalities varies in the time-interval of interest. A simple example will illustrate this. Let a rigid bar have one end hinged to the ceiling and the other end free. The hinge constrains the motion of the bar in that the end is not allowed to leave the ceiling, a bilateral constraint. The distance between the end and the ceiling is zero. The motion of the other end may be constrained by the floor. Either the distance between the floor and that end is positive or it is zero when the bar is in contact with the floor, a unilateral constraint. Another simple example is a ball rotating about a point by means of a string. If the string is slack, then no force from the string acts on the ball. The distance between the point and the ball is less than the length of the string. Otherwise, there is a force in the string preventing the distance between the ball and the point from being greater than the length of the string. More complicated systems of this kind are collapsing buildings, Lötstedt and Dahlquist (1978); and sliding rocks in a tunnel or on a slope, Cundall (1974).

Some mathematical properties of mechanical systems subject to unilateral con-
straints will be investigated here. The theory of complementarity in mathematical
programming has proven to be useful in this context. This theory and applications to
mechanics are surveyed in Cottle (1979). Other related areas where differential
equations and a number of constraints are to be satisfied are discretized problems in
variational inequalities, Duvaut and Lions (1972), and discretized contact problems,
Kalker (1977).

In § 2 the governing equations are stated and the discontinuities in the solution are discussed. The boundedness of the solution is studied in § 3. The equations in time-intervals where the number of constraints satisfied as equalities is constant are established in § 4. In § 5 the solution in the neighborhood of points at which this number is changed is constructed. The final section contains possible generalizations and some comments on the numerical solution of the problem.

**2. The governing equations.** The system of ordinary differential equations and inequalities describing the motion of a system of rigid bodies in two or three space
dimensions subject to scleronomic, holonomic and nonholonomic constraints are, for $t \ge 0$,

$$ (2.1a) \qquad M(q)\ddot{q} = f(q, \dot{q}, t) + G(q)\lambda(t), $$

$$ (2.1b) \qquad \phi(q) \ge 0, $$

$$ (2.1c) \qquad G_2^T(q) \dot{q} \ge 0. $$

The equations without constraints are found in Wittenburg (1977). Kilmister and Reeve (1966) and Pérès (1953) give conditions to be satisfied when unilateral constraints are included. $q \in \mathbf{R}^m$ is the vector of coordinates, and $f \in \mathbf{R}^m$ contains driving and inertia forces. We assume that $M \in \mathbf{R}^{m \times m}$ is symmetric, positive definite. However, it is possible that for certain values of $q$, $M$ is only semidefinite; see, e.g., Huston and Passerello (1976) and Wittenburg (1977). This difficulty can be avoided by redefining the coordinate system in the neighborhood of such a singularity. The components of $\phi \in \mathbf{R}^N$ are the holonomic, unilateral constraint functions. $G_2 \in \mathbf{R}^{m \times p}$ defines the nonholonomic constraints. The transposition of a matrix or vector $c$ is written $c^T$. Let the Jacobian of $\phi$ be denoted by

$$ G_1 = \frac{\partial \phi}{\partial q} \in \mathbf{R}^{m \times N} $$

and a column of $G_1$ and $G_2$ by $g_{1i}$ and $g_{2i}$ respectively. Furthermore, let the constraints be numbered consecutively

$$ G_1 = (g_{11}, g_{12}, \dots, g_{1N}) \quad \text{and} \quad G_2 = (g_{2,N+1}, g_{2,N+2}, \dots, g_{2,N+p}). $$

As an example, $\phi_i \ge 0$ may be the geometrical condition preventing a corner of one body from penetrating the edge or surface of another body. $\lambda_i$ is proportional to the normal force acting on the bodies at the contact point.

A constraint $i$ is called active if $\phi_i = 0$ or $g_{2i}^T \dot{q} = 0$, while for a passive constraint $j$, $\phi_j > 0$ or $g_{2j}^T \dot{q} > 0$. Let those constraints which were active and passive, respectively, in the beginning of the interval $T$, remain in the same state throughout the whole interval. Then

$$ JT = \{j | \phi_i = 0 \text{ or } g_{2j}^T \dot{q} = 0, t \in T\} $$

is called the active set in $T$. The Lagrange multiplier $\lambda_i$ corresponding to a constraint $i$, satisfies the following complementarity conditions

$$ (2.2a) \qquad \lambda_i \ge 0, \quad \phi_i \ge 0, \quad \lambda_i \phi_i = 0, $$

$$ (2.2b) \qquad \lambda_i \ge 0, \quad g_{2i}^T \dot{q} \ge 0, \quad \lambda_i g_{2i}^T \dot{q} = 0. $$

The physical interpretation of (2.2a) in the example above is that if the two bodies are in contact then $\phi_i = 0$ and $\lambda_i \ge 0$. If the corner is not in contact with the other body then $\phi_i > 0$ and $\lambda_i = 0$.

Let us investigate what happens when an active constraint $i \le N$ becomes passive at $t = t_0$. Assume that the first and second derivatives of $\phi$ exist. For $i \in J(0, t_0]$,

$$ \phi_i = 0, $$

$$ \frac{d\phi_i}{dt} = g_{1i}^T \dot{q} = 0, $$

$$ \frac{d^2\phi_i}{dt^2} = g_{1i}^T \ddot{q} + \dot{g}_{1i}^T \dot{q} = 0. $$
Suppose that for $t > t_0$ we have $\phi_i > 0$. Then, at least for a short interval, $(t_0, t_0 + \Delta t'_0)$, $\Delta t'_0 > 0$, $g_{1i}^T \dot{q} > 0$ since

$$ \phi_i(t) = g_{1i}^T(t) \dot{q}(t) (t - t_0) + O((t - t_0)^2) > 0. $$

Similarly, it follows that there is a $\Delta t_0 < \Delta t'_0$ such that in $(t_0, t_0 + \Delta t_0)$,

$$ g_{1i}^T \ddot{q} + \dot{g}_{1i}^T \dot{q} > 0. $$

$\lambda_i$ is nonnegative as long as $\phi_i = 0$, but at $t = t_0$, $\lambda_i(t_0) = 0$, which prevails until $\phi_i$ vanishes again. If $\partial g_{2i}/\partial q_i$ exist for all $j$, then the same arguments can be applied to the constraints (2.1c). Hence we have with $g_i = g_{ji}$, $i \in J(0, t_0]$, $j = 1$ or 2, in $(0, t_0 + \Delta t'_0]$,

$$ (2.3) \qquad \lambda_i g_i^T \dot{q} = 0, \quad \lambda_i \ge 0, \quad g_i^T \dot{q} \ge 0, $$

and with $\delta_i = g_i^T \ddot{q} + \dot{g}_i^T \dot{q}$ in $(0, t_0 + \Delta t_0]$,

$$ (2.4) \qquad \lambda_i \delta_i = 0, \quad \lambda_i \ge 0, \quad \delta_i \ge 0. $$

If two or more constraints become passive simultaneously each one of them satisfies the relations (2.3) and (2.4). Note that in $[0, \infty)$, $\lambda_i$, $g_i^T \dot{q}$ and $\delta_i$ fulfil the relations

$$ \lambda_i \ge 0, \quad \lambda_i g_i^T \dot{q} = 0, \quad \lambda_i \delta_i = 0, $$

without the nonnegativity conditions on $g_i^T \dot{q}$ for $i \le N$ and on $\delta_i$ for all $i$. When $\phi_i > 0$, both $g_i^T \dot{q}$ and $\delta_i$ may be negative.

Insert $\ddot{q}$ from (2.1a) in the definition of $\delta_i$. Let $G$ consist of the columns $g_i$, $i \in J(0, t_0]$ and $\lambda$, $\delta$ of those $\lambda_i$, $\delta_i$ for which $i \in J(0, t_0]$. Then it follows from the discussion that the relations satisfied in an interval $[0, t_0 + \Delta t_0]$ where previously active constraints become passive are

$$ (2.5a) \qquad M \ddot{q} = f + G \lambda, $$

$$ (2.5b) \qquad G^T M^{-1} G \lambda + G^T M^{-1} f + G^T \dot{q} = \delta \ge 0, $$

$$ \lambda \ge 0, \quad \lambda^T \delta = 0. $$

(2.5) is a system of second order ordinary differential equations (2.5a) and a linear complementarity problem (LCP) (2.5b). See Cottle and Dantzig (1968) and Karamar-dian (1972) for an introduction to LCP. Since $\lambda_i = 0$, $i \in J(0, t_0] \setminus J(t_0, t_0 + \Delta t_0]$, in $(t_0, t_0 + \Delta t_0]$ the corresponding column $g_i$ in $G$ can be removed without affecting the solution. The complementarity problem (2.2a), which is valid for all $t$, is in general nonlinear in $\lambda$ since $\phi_i = \phi_i(q)$, where $q$ is determined by (2.5a).

A constraint is added to the set of active constraints when a previously positive constraint function $\phi_i$ or $g_{2i}^T \dot{q}$ attains zero. In general, some derivative of the solution is discontinuous when a constraint is added. When the corner reaches the surface of the other body in the example above, a collision occurs and the velocity $\dot{q}$ and the relative velocity $g_{1i}^T \dot{q}$ have jumps. There is no problem in deciding when and which constraint to add to the active set since there is a constraint function $\phi_i$ or $g_{2i}^T \dot{q}$ to base the decision on. The problem of dropping constraints is more delicate and will be treated further in § 5.

In order to discuss the discontinuities in the solution due to the changes in the active set let us consider the system

$$ (2.6a) \qquad M \ddot{q} = f + G \lambda, $$

$$ (2.6b) \qquad G^T \dot{q} = \Delta, $$

$$ (2.6c) \qquad G^T \ddot{q} + \dot{G}^T \dot{q} = \dot{\Delta} = \delta. $$
Let a superscript 1 or 2 on a variable denote its value at $t_0-0$ or $t_0+0$, respectively, where $t_0$ is the point at which the active set is altered. Assume that $M, f$ and $G$ and their first derivatives and the second derivatives of $G$ are continuous functions of their arguments.

It is intuitively evident that $q(t)$ in the rigid body model is continuous, but discontinuities in $\dot{q}, \ddot{q}$ and $\dddot{q}$ are possible,

1. Let $\phi_i \to 0$ continuously so that for $t < t_0$, $\phi_i(t) > 0$ but $\phi_i(t_0) = 0$. Then by (2.6b) the corresponding component in $\Delta$, $\Delta_i$, in general, satisfies

$$ \Delta_i^1 < 0 \quad \text{and} \quad \Delta_i^2 \ge 0. $$

If the collision is inelastic, then $\Delta_i^2 = 0$. Since

$$ \Delta_i^2 - \Delta_i^1 = (G^2 - G^1)^T \dot{q}^2 + (G^1)^T (\dot{q}^2 - \dot{q}^1) > 0 $$

and $G^2 = G^1$, $\dot{q}$ is also discontinuous.

2. If no holonomic constraint becomes active in the neighborhood of $t_0$ then $\dot{q}(t)$ is continuous, Kilmister and Reeve (1966). Suppose that for a nonholonomic constraint $g_{2i}^T \dot{q} \to 0$ so that when $t < t_0$, $g_{2i}^T \dot{q}(t) > 0$, but $g_{2i}^T \dot{q}(t_0) = 0$. The corresponding component in $\delta, \delta_i$, is discontinuous in general and

$$ \delta_i^1 < 0 \quad \text{but} \quad \delta_i^2 = 0. $$

After inserting $\ddot{q}$ from (2.6a) into (2.6c), the equation for $\lambda$ and $\delta$ is as in (2.5),

$$ (2.7) \qquad G^T M^{-1} G \lambda + G^T M^{-1} f + \dot{G}^T \dot{q} = \delta. $$

$M^{-1}$ has the same continuity properties as $M$ by Lemma A (in the Appendix), $\dot{q}$ and $\dot{G}$ are continuous and thus,

$$ (2.8) \qquad b(q, \dot{q}, t) = G^T M^{-1} f + \dot{G}^T \dot{q} $$

is continuous as well. It follows from the above arguments that $\lambda^2 - \lambda^1 \ne 0$ and $\lambda$ is discontinuous. Furthermore, (2.6c) implies that $\ddot{q}$ has a jump.

3. Suppose that no constraint becomes active in the neighborhood of $t_0$ and that $\delta_i$ is continuous when a previously active constraint $j$ becomes passive at $t_0$. Then $\delta_j(t) = 0$, $t \le t_0$, and $\delta_j(t) > 0$, $t > t_0$. In order to satisfy the complementarity conditions on $\lambda_j$ (2.4), the following holds: $\lambda_j \ge 0$, $t \le t_0$ and $\lambda_j = 0$, $t > t_0$. By (2.7) and the continuity of $q$ and $\dot{q}$, $\lambda^2$ and $\lambda^1$ satisfy

$$ G^T M^{-1} G (\lambda^2 - \lambda^1) = \delta^2 - \delta^1 = 0. $$

Hence, $\lambda^2 - \lambda^1 \in N(G)$, the null space of $G$, and since $0 \in N(G)$, $\lambda$ can be chosen continuous at $t_0$. In particular, $\lambda_j(t_0) = 0$. However, $\dot{\lambda}$ is usually discontinuous; $\dot{\lambda}_j^1 < 0$, but $\dot{\lambda}_j^2 = 0$. $\ddot{q}$ is continuous by (2.6a). By the smoothness of $M, G$ and $f, \dot{b}(q, \dot{q}, t)$ and $d/dt$ ($G^T M^{-1} G$) are continuous. The time derivative of (2.7), after introducing (2.8), is

$$ (2.9) \qquad G^T M^{-1} G \dot{\lambda} + \frac{d}{dt} (G^T M^{-1} G) \lambda + \dot{b} = \dot{\delta}. $$

In general, $\dot{\delta}$ is discontinuous, and from (2.6a), we infer that $\ddot{q}$ also is discontinuous.

In the first case the jump condition for $\dot{q}$ is, Kilmister and Reeve (1966),

$$ (2.10a) \qquad M(q)(\dot{q}^2 - \dot{q}^1) = G(q)\Lambda + F(q, \dot{q}, t_0), $$

where $G$ consists of those columns $g_i$ associated with active constraints at $t_0$ and $F \in \mathbf{R}^m$ contains externally generated impulses at $t_0$. External impulses are the counterparts of external forces in (2.1a). The Lagrange multipliers $\Lambda$ correspond to impulses in the rigid
body system. In the inelastic case they are determined such that the constraints, after the discontinuity, are fulfilled;

$$ (2.10b) \qquad G^T \dot{q}^2 = \Delta^2 \ge 0, \quad \Lambda \ge 0, \quad \Lambda^T \Delta^2 = 0, $$

or, equivalently,

$$ (2.11) \qquad G^T M^{-1} G \Lambda + G^T M^{-1} F + G^T \dot{q}^1 = \Delta^2 \ge 0, \quad \Lambda \ge 0, \quad \Lambda^T \Delta^2 = 0. $$

It is shown by Ingleton (1966) and Cottle (1968) that this LCP always possesses a solution $\Lambda$ and a unique $\Delta^2$ (cf. Lemma 5.2). There are no particular jump conditions in the other two cases, when $\ddot{q}$ and $\ddot{\bar{q}}$ are discontinuous. Equation (2.5) is to be satisfied on both sides of the discontinuity.

In addition to these jumps caused by the changes in the active set, $G$ and $f$ may be sources of discontinuities in $\dot{q}, \ddot{q}$ and $\ddot{\bar{q}}$. A jump in $\dot{q}$ determined by (2.10) will sometimes be necessary to fulfil all the constraints after such a discontinuity. Consider, e.g., the case when $g_i^1 \neq g_i^2$, $i \in J(0, t_0)$. $(g_i^T \dot{q})^1 = 0$ and if $(g_i^2)^T \dot{q}^1 > 0$ then the constraint $i$ is no longer active for $t > t_0$. $\dot{q}^1$ is a solution of (2.10) with $\Lambda = 0$. Hence, with $F = 0$ in (2.10a), then $\dot{q}^2 = \dot{q}^1$ and $\dot{q}$ is continuous. Conversely, if $(g_i^2)^T \dot{q}^1 < 0$, then $\dot{q}$ is discontinuous satisfying a jump condition such as (2.10). Simple examples are particles sliding on plane surfaces which have $L$-shaped corners. Also externally generated impulses cause discontinuities in $\dot{q}$. It follows from (2.6a) that $\ddot{q}$ has a jump when, for some $i$, $\lambda_i^1 > 0$, but $\lambda_i^2 = 0$ because of the geometry of the system. The corner of a body may, e.g., slide off its support.

**3. Bounds on the velocity vector.** In this section simple bounds on $\dot{q}$ are derived which are independent of the changes of the active set. A norm is introduced which is related to the kinetic energy of the system. For $v \in \mathbb{R}^n$, $\|v\|$ denotes the Euclidean norm and for $A \in \mathbb{R}^{n \times n}$, $\|A\|$ denotes the spectral norm.

Let $\alpha, \beta$ and $\gamma$ be positive and $q_0 = q(0), \dot{q}_0 = \dot{q}(0)$. Assume that for

$$ (3.1) \qquad \|q\| \le \|q_0\| + \alpha, \quad \|\dot{q}\| \le \|\dot{q}_0\| + \beta, \quad t \in [0, \gamma], $$

the following estimates are valid

$$ (3.2) \qquad a^2 \|x\|^2 \le x^T M x \le b^2 \|x\|^2, \quad \|f\| \le c, \quad \|\dot{M}\| \le d \|\dot{q}\|. $$

Furthermore, assume that there are no external impulses, no discontinuities in $G$ and that no holonomic constraints are added to the active set. Then the discussion in the end of the previous section shows that $\dot{q}$ is continuous.

Since $M$ is symmetric and positive definite, the Cholesky factorization of $M$ exists;

$$ (3.3) \qquad M = CC^T. $$

It follows from the construction of $C$, Dahlquist and Björck (1974), that $C$ is nonsingular and has the same smoothness properties as $M$. Let

$$ (3.4) \qquad S = \|C^T \dot{q}\|. $$

A bound on $S$ will be computed and if $t$ is sufficiently small then $q$ and $\dot{q}$ will remain within the bounds (3.1).

By (3.2), $S$ satisfies

$$ (3.5) \qquad a \|\dot{q}\| \le S \le b \|\dot{q}\|, $$

and by (2.1a),

$$ \frac{d}{dt} S^2 = 2S\dot{S} = 2\dot{q}^T M \ddot{q} + \dot{q}^T \dot{M} \dot{q} = 2\dot{q}^T f + 2\dot{q}^T G \lambda + \dot{q}^T \dot{M} \dot{q}. $$
The term $\dot{q}^T G \lambda$ vanishes identically since $\lambda$ and $G^T \dot{q}$ are complementary (2.3) (or “constraint forces $G\lambda$ perform no work”). If $\|\dot{q}\|$ satisfies (3.1) then by (3.2) there is a bound $d'$ such that $\|\dot{M}\| \le d'$. Therefore,

$$ (3.6) \qquad 2SS = 2\dot{q}^T f + \dot{q}^T \dot{M} \dot{q} \le 2\|\dot{q}\|\|f\| + \|\dot{M}\|\|\dot{q}\|^2 \le (2c/a)S + (d'/a^2)S^2. $$

Suppose that $S > 0$, divide the first and the last part of the relation (3.6) by $S$ and use Gronwall's lemma to establish the estimate

$$ (3.7) \qquad S(t) \le (S(0) + 2ac/d') \exp(0.5d't/a^2) - 2ac/d'. $$

In the neighborhood of $S=0$, and thus $\dot{q}=0$, $S^2$ is majorized by $W$ satisfying

$$ \dot{W} = 2\left(k + \frac{c}{a}\right)\sqrt{W} $$

for some $k > 0$. Hence, when $S$ is small

$$ S(t) \le \left(k + \frac{c}{a}\right)t. $$

Combining (3.5) and (3.7) yields an upper bound on $\|\dot{q}(t)\|$:

$$ (3.8) \qquad \|\dot{q}(t)\| \le S(0) \frac{\exp(c_1 t)}{a} + c_2(\exp(c_1 t) - 1), $$

where $c_1$ and $c_2$ are constants. This bound is independent of the number of active constraints and would be the same if the description of the mechanical system (2.1) contained no algebraic constraints at all. Moreover, a change in the rank of $G$, which sometimes occurs without $\dot{q}$ being discontinuous, has no influence on the bound.

When holonomic constraints $\phi_i$ have become active, after discontinuities in $G$ and external impulses, $\Lambda$ in (2.10) can, e.g., be determined according to Newton's impulse law, Kilmister and Reeve (1966), such that

$$ (3.9) \qquad \text{if } (g_i^2)^T \dot{q}^1 < 0 \quad \text{then} \quad (g_i^2)^T \dot{q}^2 = -e_i (g_i^2)^T \dot{q}^1, \quad 0 \le e_i \le 1, $$

where $e_i$ is the coefficient of restitution and superindices are to be interpreted as in § 2. If $e_i = 1$ then we have the perfectly elastic case, and if $e_i = 0$ then the collision is inelastic. Let $J_0 = \{i | \text{ the impulse law is (3.9), } \phi_i(t_0) = 0 \text{ and } (g_{1i}^2)^T \dot{q}^1 < 0 \text{ or } (g_{2i}^2)^T \dot{q}^1 < 0 \}$, and assume that the impulses at the other active constraints at $t_0$ and the other constraints $i$ for which $(g_{2i}^2)^T \dot{q}^1 < 0$ satisfy the complementarity condition,

$$ (3.10) \qquad \Lambda_i \ge 0, \quad (g_i^2)^T \dot{q}^2 \ge 0, \quad \Lambda_i (g_i^2)^T \dot{q}^2 = 0. $$

Furthermore, let

$$ J_1 = \{i | \phi_i(t_0) = 0 \text{ or } (g_{2i}^2)^T \dot{q}^1 \le 0\}. $$

$J_1$ is the set of all constraints $i$ for which a $\Lambda_i$ must be computed in (2.10). Observe that $\Lambda_i (g_i^2)^T \dot{q}^1 \le 0$ and for $i \in J_1 \setminus J_0$, $\Lambda_i (g_i^2)^T \dot{q}^2 = 0$ and that $J_0$ or $J_1 \setminus J_0$ may be empty. Consider the difference;

$$ (3.11) \qquad \Delta S^2 = (S^2)^2 - (S^1)^2 = (\dot{q}^2)^T M \dot{q}^2 - (\dot{q}^1)^T M \dot{q}^1 = (\dot{q}^2 + \dot{q}^1)^T M (\dot{q}^2 - \dot{q}^1). $$

$G$ in (2.10) consists of those $g_i^2$ such that $i \in J_1$. Insert (2.10) into the last expression in
(3.11) and use (3.9) and (3.10);

$$
\begin{align*}
\Delta S^2 &= (\dot{q}^2 + \dot{\dot{q}}^1)^T G \Lambda + (\dot{q}^2 + \dot{\dot{q}}^1)^T F \\
&= \sum_{i \in J_1} (\Lambda_i (g_i^2)^T \dot{q}^2 + \Lambda_i (g_i^2)^T \dot{\dot{q}}^1) + (\dot{q}^2 + \dot{\dot{q}}^1)^T F \\
&= \sum_{i \in J_0} (1 - e_i) \Lambda_i (g_i^2)^T \dot{q}^1 + \sum_{i \in J_1 \setminus J_0} \Lambda_i (g_i^2)^T \dot{\dot{q}}^1 + (\dot{q}^2 + \dot{\dot{q}}^1)^T F \\
&\leq (\dot{q}^2 + \dot{\dot{q}}^1)^T F.
\end{align*}
$$

Hence, if $F = 0$ then $S^2 \leq S^1$ and $\|C^T \dot{q}^2\| \leq \|C^T \dot{q}^1\|$.

Discontinuities in $\dot{q}$ due to the constraints will not increase $\dot{q}$ in the weighted norm defined by $M$ and $C$. Therefore, the upper bounds in (3.7) and (3.8) are completely independent of the constraints.

**4. Between the changes of the active set.** The motion of many rigid body systems is such that there are finite intervals $(t_1, t_2)$, $t_2 > t_1$, where no change in the active set takes place. Then the active constraints can be treated as if they were bilateral. A simple example where these intervals become infinitesimally small is a ball bouncing on a plane with a coefficient of restitution $0 < e < 1$. The elevation of the ball over the plane decreases geometrically and so does the length of the intervals until the ball is at rest on the plane. A point where the number of alterations of the active set is not finite in a finite interval containing the point, is termed an accumulation point. The bounds obtained in the previous section are not affected by accumulation points.

The equations satisfied in $(t_1, t_2)$ are

$$
(4.1a) \qquad M\ddot{q} = f + G\lambda,
$$

$$
(4.1b) \qquad \phi_i = 0, \quad i \in J(t_1, t_2),
$$

$$
(4.1c) \qquad g_{2i}^T \dot{q} = 0, \quad i \in J(t_1, t_2).
$$

$G$ consists of those columns $g_i$ for which $i \in J(t_1, t_2)$. Since $g_{1i}^T \dot{q} = 0, i \in J(t_1, t_2)$, (4.1b) and (4.1c) are equivalent to

$$
(4.2) \qquad G^T \dot{q} = 0,
$$

or including (4.1a) and taking the time-derivative of (4.2),

$$
(4.3a) \qquad M\ddot{q} = f + G\lambda,
$$

$$
(4.3b) \qquad G^T \ddot{q} + \dot{G}^T \dot{q} = 0.
$$

Insert $\ddot{q}$ from (4.3a) into (4.3b)

$$
(4.4) \qquad G^T M^{-1} G \lambda + G^T M^{-1} f + \dot{G}^T \dot{q} = 0.
$$

If $G$ has full column rank $(G^T M^{-1} G)^{-1}$ exists and $\lambda$ can be solved from (4.4);

$$
(4.5) \qquad \lambda = -(G^T M^{-1} G)^{-1} (G^T M^{-1} f + \dot{G}^T \dot{q}).
$$

After substituting $\lambda$ from (4.5) into (4.3a), the governing system of differential equations in $(t_1, t_2)$ is

$$
(4.6) \qquad M\ddot{q} = f - G(G^T M^{-1} G)^{-1} G^T M^{-1} f - G(G^T M^{-1} G)^{-1} \dot{G}^T \dot{q}.
$$

Using (3.3) and the notation $B = C^{-1}G$, a more convenient form is

$$
(4.7) \qquad C^T \ddot{q} = (I - B(B^T B)^{-1} B^T) C^{-1} f - B(B^T B)^{-1} \dot{G}^T \dot{q}.
$$
Let

$$
(4.8) \qquad P = (I - B(B^T B)^{-1} B^T)
$$

and $R(A)$ and $N(A)$ denote the range and nullspace of $A$, respectively. It follows that
$P^2 = P$ and $P^T = P$. $P$ is an orthogonal projection matrix, Ben-Israel and Greville
(1974), where

$$
\begin{align*}
R(P) &= N(B^T) = N(G^T C^{-T}), \\
N(P) &= R(B) = R(C^{-1}G).
\end{align*}
$$

Hence, the first term on the right-hand side of (4.7) lies in $N(G^T C^{-T})$ and the second
term belongs to the orthogonal space $R(C^{-1}G)$.

The classical uniqueness theorem for ordinary differential equations can be applied
to the particular case (4.6) to obtain:

THEOREM 4.1. Suppose that $M(q)$, $G(q)$, $\partial G(q)/\partial q$ and $f(q, \dot{q}, t)$ are bounded and Lipschitz continuous in $q$ and $\dot{q}$, that $f(q, \dot{q}, t)$ is continuous in $t$ and that $G$ has full column rank in the neighborhood of $t_0 \in (t_1, t_2)$, $q_0$ and $\dot{q}_0$. Then there is a unique solution $q(t)$, $\dot{q}(t)$ of (4.6) in a neighborhood of $t_0$ satisfying $q(t_0) = q_0$ and $\dot{q}(t_0) = \dot{q}_0$.

*Proof.* By Lemma A and Theorem 2.3 in Coddington and Levinson (1955) the statement follows. □

The elements of *M* and *G* are often sums of products of constants, sin *q*ᵢ and cos *q*ᵢ,
where *q*ᵢ represents some angle. Then *M*, *G* and ∂*G*/∂*q* will all be bounded and
Lipschitz continuous. Thus, if the columns of *G* are linearly independent in [*t*₁, *t*₂] and
*f*(*q*, *dot;q*, *t*) is smooth and bounded, *ddot;q* will remain bounded as well as *dot;q* and *q*. λ(*t*) in (4.5) is
continuous, and if λᵢ(*t*₁) > 0 for all active constraints at *t*₁ then there is certainly an
interval (*t*₁, *t*₂) where λᵢ > 0 and the active set is unaltered.

5. The existence of λ and a complementary δ. In the previous section we assumed that G had linearly independent columns. A simple example where this restriction is not satisfied is a table with four or more legs on a floor. Therefore, that assumption will be removed when the complementarity problem (2.5) is studied further to show that with sufficiently smooth M, G and f a solution to (2.5) can always be constructed.

To begin with, three lemmas are stated and proved. They will be used in the sequel
to investigate the properties of the solution of (2.5).

LEMMA 5.1. Let $x \ge 0$ be the solution of

$$
(5.1) \qquad G^T M^{-1} G x = z.
$$

$G \in \mathbf{R}^{m \times n}$, rank $(G) = r \le \min (m, n)$ and $M \in \mathbf{R}^{m \times m}$ is symmetric, positive definite. There are $k \le r$ columns of $G$ which are linearly independent and an $x' \in \mathbf{R}^k$, $x_i' > 0$, such that for $G' \in \mathbf{R}^{m \times k}$, consisting of these $k$ columns, $Gx = G'x'$.

*Proof.* Apply Lemma 4 in Cottle (1968) to (5.1). Then there is an $x_* \ge 0$ such that

$$
G^T M^{-1} G x = G^T M^{-1} G x_*, \quad x - x_* \in N(G).
$$

There is also a nonsingular $G'^T M^{-1} G'$ with rank $(G'^T M^{-1} G') = \text{rank}(C^{-1} G') =$
rank $(G') = k$, the number of $x_{*i} > 0$. $Gx = Gx_* = G'x'$ and the lemma is proved. $\square$

In some mechanical applications the following problem is of interest (see Cottle
(1979)):
Let $M \in \mathbf{R}^{m \times m}$, $b \in \mathbf{R}^m$, $G = (G_1, G_2)$, $G_1 \in \mathbf{R}^{m \times p}$, $G_2 \in \mathbf{R}^{m \times r}$, $n = p+r$, $\lambda^T = (\lambda_1^T, \lambda_2^T)$ and $\delta^T = (\delta_1^T, \delta_2^T)$, $\lambda_1$, $\delta_1 \in \mathbf{R}^p$, $\lambda_2$, $\delta_2 \in \mathbf{R}^r$. Find a solution of

$$
(5.2a) \qquad G_1^T M^{-1} G_1 \lambda_1 + G_1^T M^{-1} G_2 \lambda_2 + G_1^T b = 0,
$$

$$
(5.2b) \qquad G_2^T M^{-1} G_1 \lambda_1 + G_2^T M^{-1} G_2 \lambda_2 + G_2^T b = \delta_2,
$$

$$
(5.2c) \qquad \lambda_2 \ge 0, \quad \delta_2 \ge 0, \quad \lambda_2^T \delta_2 = 0.
$$

LEMMA 5.2. Assume that rank ($G_1$) = $q \leq \min (m, p)$ and rank ($G_2$) = $s \leq \min (m, r)$. Then a solution to (5.2) always exists with the following properties:

1. $\delta$ and $G\lambda$ are unique. If $G$ has full column rank, then $\lambda$ is also unique.

2. Let $\lambda_i = \lambda_i^1 + \lambda_i^2$ and $(\lambda^i)^T = ((\lambda_i^1)^T, (\lambda_i^2)^T)$, $i=1,2$, where $\lambda^1 \in \mathbf{R}(G^T)$, $\lambda^2 \in N(G)$. Then $\lambda^1$ is unique and $\lambda^2$ can be chosen arbitrarily but such that $\lambda_2 = \lambda_2^1 + \lambda_2^2 \ge 0$.

3. There are subsets $J_1 \subset \{1, 2, \dots, p\}$ and $J_2 \subset \{1, 2, \dots, r\}$ of cardinality $k \le q$ and $l \le s$, respectively, such that $G'_1 = (g_{1i_1}, g_{1i_2}, \dots, g_{1i_k})$, $i_j \in J_1$, rank $(G'_1) = k$, and $G'_2 = (g_{2i_1}, g_{2i_2}, \dots, g_{2i_l})$, $i_j \in J_2$, rank $(G'_2) = l$. There is a $\lambda'_1 \in \mathbf{R}^k$ and a $\lambda'_2 \in \mathbf{R}^l$ such that $\lambda'_{2i} > 0$. $\lambda'_1$ and $\lambda'_2$ satisfy $G_1\lambda_1 = G'_1\lambda'_1$ and $G_2\lambda_2 = G'_2\lambda'_2$ for any $\lambda_1$ and $\lambda_2$ solving (5.2). The rank of $G' = (G'_1, G'_2)$ is $k+l$.

4. There is a unique solution $\lambda_0$ such that $\|\lambda_0\| \le \|\lambda\|$ for every solution $\lambda$.

*Proof.* 1. To simplify the notation, let $B_1 = C^{-1}G_1$ and $B_2 = C^{-1}G_2$. Any particular $\lambda_2$, $\lambda_1$ in (5.2a) can be regarded as the solution of a linear least squares problem

$$
(5.3) \qquad \min_{\lambda_1} \|B_1\lambda_1 + B_2\lambda_2 + C^T b\|.
$$

The solution of (5.3) satisfying $\min \|\lambda_1\|$ can be written, Ben-Israel and Greville (1974),

$$
(5.4) \qquad \lambda_1 = -B_1^+ (B_2 \lambda_2 + C^T b),
$$

where $A^+$ denotes the unique Moore-Penrose inverse of $A$. Furthermore, after eliminating $\lambda_1$ from (5.2b), the LCP to be solved for $\lambda_2$ and $\delta_2$ is

$$
(5.5) \qquad (PB_2)^T (PB_2) \lambda_2 + (PB_2)^T C^T b = \delta_2, \quad \lambda_2 \ge 0, \quad \delta_2 \ge 0, \quad \lambda_2^T \delta_2 = 0,
$$

where $P = (I - B_1 B_1^+)$, an orthogonal projection matrix. According to Theorem 1 in
Cottle (1968), a solution to (5.5) always exists and $PB_2 \lambda_2$, $\delta_2$ and $\delta$ are unique. Since

$$
G\lambda = C(B_1\lambda_1 + B_2\lambda_2) = C(PB_2\lambda_2 - B_1B_1^+C^Tb),
$$

$G\lambda$ is also unique.

If $G$ has full column rank, then $G^T M^{-1} G$ is nonsingular, and since $G_1^T b, G_2^T b$ and $\delta_2$ are unique in (5.2), $\lambda_1$ and $\lambda_2$ are also unique.

2. $G\lambda$ is unique. Hence, $z$ in

$$
(5.6) \qquad G^T M^{-1} G \lambda = z
$$

is also unique.

From linear algebra we know that $\lambda^1$ is unique. $\delta \in R(G^T)$ and for any $\lambda^2$,

$0 = \delta^T \lambda^2 = \delta_2^T \lambda_2^2.$

Since

$$
0 = \delta_2^T \lambda_2 = \delta_2^T \lambda_2^1 + \delta_2^T \lambda_2^2 = \delta_2^T \lambda_2^1,
$$

(5.2c) is satisfied for any $\lambda_2^2$ such that $\lambda_2 = \lambda_2^1 + \lambda_2^2 \ge 0$.

3. It is easily shown that for any solution $\lambda$ to (5.6) there is a solution $\lambda'$ and a matrix $G'_1$ with linearly independent columns and rank $(G') = k \le s$ such that $G_1\lambda_1 = G'_1\lambda'_1$.
From (5.5) and part 1 of the proof it follows that $z_3$ in

$$ (5.7) \qquad (PB_2)^T(PB_2)\lambda_2 = z_3 $$

is unique. Apply Lemma 5.1 with $M = I$, the identity matrix, to (5.7). Let $B'_2$ consist of those columns in $B_2$ such that $(PB_2)' = PB_2'$. rank $(PB_2')$ satisfies

$$ l = \operatorname{rank}(PB_2') \le \operatorname{rank}(PB_2) \le \operatorname{rank}(B_2) = s. $$

According to the definition of $P$, $R(P) = N(B_1^T)$ and $N(P) = R(B_1)$. A column $b_i$ of $B_2'$ can be split into two parts;

$$ b_i = b_i^1 + b_i^2, \quad b_i^1 \in R(B_1), \quad b_i^2 \in N(B_1^T). $$

Since the columns of $PB_2'$ are linearly independent, $b_i^2$, $i = 1, 2, \dots, l$, are linearly independent. Therefore, the columns of $B'_2$ and $G'_2$ are linearly independent. Furthermore, since the columns of $C^{-1}G'_1$ belong to $R(B_1)$, $G'$ has linearly independent columns and rank $(G') = k + l$.

4. Choose $R > 0$ such that $\|\lambda_2\| \le R$, where $\lambda_2$ is the lower part of a solution $\lambda$ of (5.6). The lower parts of the solutions of (5.6) satisfying $\|\lambda_2\| \le R$ and $\lambda_2 \ge 0$ form a closed, convex set. The convex function $\|\lambda_2\|^2$ has a unique minimum $\lambda_{20}$ on that set. $B_2\lambda_2$ in (5.4) is unique. Hence, $\lambda_{10}$ from (5.4) is unique and $\lambda_0 = (\lambda_{10}^T, \lambda_{20}^T)$. $\square$

The next lemma is used to show that the LCP in (2.5) has a solution at a point $t = t_0$.

**LEMMA 5.3.** *Assume that the rank of $G(t)$ is constant in the neighborhood of $t = t_0$ and that $G$ has a sufficient number of derivatives. Then the following holds at $t = t_0$: If*

$$ \frac{d^p}{dt^p} (G^T v) = 0, \quad p = 0, 1, \dots, k-1, \quad k \ge 1, $$

then

1. If $X(t)$ has at least one derivative, then $(d^p/dt^p)(X^T G^T v) = X^T (d^p/dt^p)(G^T v) = 0$. $p = 0, 1, 2, \dots, k-1$.

2. $(d^k/dt^k)(G^T v) \in R(G^T)$.

*Proof.* 1. For $p = 0$ the statement is trivially true. Suppose that

$$ \frac{d^p}{dt^p} (G^T v) = 0 \quad \text{and} \quad \frac{d^p}{dt^p} (X^T G^T v) = X^T \frac{d^p}{dt^p} (G^T v) = 0. $$

Then

$$ (5.8) \qquad \frac{d^{p+1}}{dt^{p+1}} (X^T G^T v) = \dot{X}^T \frac{d^p}{dt^p} (G^T v) + X^T \frac{d^{p+1}}{dt^{p+1}} (G^T v) = X^T \frac{d^{p+1}}{dt^{p+1}} (G^T v). $$

If $(d^{p+1}/dt^{p+1})(G^T v) = 0$, then $(d^{p+1}/dt^{p+1})(X^T G^T v) = 0$ and statement 1 follows by induction on $p$.

2. Any $G \in \mathbf{R}^{m \times n}$ with rank $(G) = r > 0$ in the neighborhood of $t_0$ can be written after permutation of the columns

$$ G = (G_r, G_r X), \quad G_r \in \mathbf{R}^{m \times r}, \quad X \in \mathbf{R}^{r \times (n-r)}, $$

where $G_r$ has full column rank. The $k$th derivative of $X^T G_r^T v$ is

$$ \frac{d^k}{dt^k} (X^T G_r^T v) = X^T \frac{d^k}{dt^k} (G_r^T v), $$
according to (5.8) with $p = k - 1$. Hence,

$$ (5.9) \qquad \frac{d^k}{dt^k} (G^T v) = \begin{pmatrix} \frac{d^k}{dt^k} (G_r^T v) \\ X^T \frac{d^k}{dt^k} (G_r^T v) \end{pmatrix} = \begin{pmatrix} y \\ X^T y \end{pmatrix}. $$

Since the columns of $G_r$ are not linearly dependent, there is a $z$ such that $y = G_r^T z$. The conclusion from (5.9) is that

$$ \frac{d^k}{dt^k} (G^T v) = \begin{pmatrix} G_r^T z \\ X^T G_r^T z \end{pmatrix} \in \mathbb{R}(G^T). \quad \square $$

The complementarity problem to be solved at a certain time point $t_0$ is (2.5b). We can assume that $G^T \dot{q} = 0$ at $t_0$, since if $g_i^T \dot{q} > 0$, then the $i$th constraint is no longer active for $t > t_0$ and if $g_i^T \dot{q} < 0$, then either $\phi_i > 0$ or an impulse must be computed, see § 2 and § 3. Applying Lemma 5.3 with $k = 1$ to $G^T \dot{q} = 0$, we find that $G^T \ddot{q} + \dot{G}^T \dot{q} \in \mathbb{R}(G^T)$ and consequently $\dot{G}^T \dot{q} \in \mathbb{R}(G^T)$. Then (2.5b) has the form of (5.2), with $p = 0$ and $G_1 = \text{nil}$. Its solution has the properties listed in Lemma 5.2. In particular, $G\lambda$ is unique, and therefore, $\ddot{q}$ in (2.5) is unique.

A natural question to pose is: If $G^T \dot{q} = 0$ at $t_0$, is it possible to construct a solution to (2.5) which is valid in an interval $[t_0, t_0 + \Delta t)$, $\Delta t > 0$? The following observation is important. It is not sufficient to solve (2.5b) at $t_0$ and let the active set for $t > t_0$ be $J_1 = \{i | \lambda_i(t_0) > 0\}$ or $J_2 = \{i | \delta_i(t_0) = 0\}$. If the first alternative is chosen, there may be an index $j$ for which $\lambda_j(t_0) = 0$, but $\dot{\lambda}_j(t_0) > 0$ and $\delta_j(t) \equiv 0$ which should belong to the active set. If $j \in J_2$, but $\dot{\delta}_j(t_0) > 0$ and $\lambda_j(t) \equiv 0$, then $j$ should not be in the active set of $J(t_0, t_0 + \Delta t)$ since

$$ \delta_j(t) = \dot{\delta}_j(t_0)(t - t_0) + O((t - t_0)^2), \quad t \ge t_0. $$

The theorem in this section proves that under certain smoothness assumptions, there is a solution to (2.5) in $[t_0, t_0 + \Delta t)$ for $\Delta t > 0$ sufficiently small. The initial conditions at $t_0$ are

$$ q(t_0) = q_0, \qquad \dot{q}(t_0) = \dot{q}_0. $$

We assume that $q_0$ and $\dot{q}_0$ are compatible with the constraints such that $\phi_i(q_0) \ge 0$ and $g_{2i}^T(q_0)\dot{q}_0 \ge 0$, and if $\phi_i(q_0) = 0$ then $g_{1i}^T(q_0)\dot{q}_0 \ge 0$ (possibly after determining an appropriate $\Lambda$ in (2.10)).

**THEOREM 5.4.** Let $J = \{i | \phi_i = 0 \text{ and } g_{1i}^T \dot{q} = 0 \text{ or } g_{2i}^T \dot{q} = 0 \text{ at } t = t_0\}$, and let $G$ consist of those $n$ columns $g_i$ for which $i \in J$. Assume that $M(q), G(q)$ and $f(q, \dot{q}, t)$ are analytic in $t$ and all $q_i, \dot{q}_i$ separately, and that $M^{-1}, G, f$ and $\partial G / \partial q$ are bounded in the neighborhood of $q_0, \dot{q}_0$ and $t_0$. Furthermore, assume that $R = \text{rank}(G)$ is constant in the vicinity of $q_0$.

Then the active set is constant in $(t_0, t_0 + \Delta t)$ for a sufficiently small $\Delta t$. There is a unique solution $q(t)$ of (2.5) which is analytic in $[t_0, t_0 + \Delta t)$. If $G(q)$ has full column rank, $\lambda(t)$ is also unique and analytic. There is a subset $J' \subset J$ of cardinality $r \le R$ such that $G' = (g_{i_1}, g_{i_2}, \dots, g_{i_r})$, $i_j \in J'$, $\text{rank}(G') = r$, and there is a $\lambda'(t) \in \mathbb{R}^r$, $\lambda'_i(t) > 0$ in $(t_0, t_0 + \Delta t)$, such that $G\lambda(t) = G'\lambda'(t)$.

**Proof.** All variables are evaluated here at $t_0$ unless otherwise stated. Let $x^{(j)}$ denote $d^j x / dt^j$ and define

$$ J_{\lambda k} = \{i | \exists p \ni 0 \le p \le k, \lambda_i^{(j)} = 0, j = 0, 1, \dots, p-1, \lambda_i^{(p)} > 0\}, $$

$$ J_{\delta k} = \{i | \exists p \ni 0 \le p \le k, \delta_i^{(j)} = 0, j = 0, 1, \dots, p-1, \delta_i^{(p)} > 0\}. $$
Solve the LCP (2.5b) at $t = t_0$ to determine $J_{\lambda 0}$ and $J_{\delta 0}$. According to Lemma 5.2, the solution $\lambda$ can be chosen so that the rank of $G' = (g_{i_1}, g_{i_2}, \cdots, g_{i_s})$, $i_j \in J_{\lambda 0}$, is $s$, the number of positive $\lambda_i$. Moreover, $G\lambda$ and $\ddot{q}$ are unique.

Let us make the following induction hypothesis:

For $i = 0, 1, 2, \cdots, k-1, k \ge 1$,

$q^{(i+2)}$ is known and unique;

$G_k = (g_{i_1}, g_{i_2}, \cdots, g_{i_q})$, $i_j \in J_k = J \setminus J_{\delta, k-1}$, and $(G_k(q))^{(i)}$ are known and unique;

$G_{k1} = (g_{i_1}, g_{i_2}, \cdots, g_{i_r})$, $i_j \in J_{\lambda, k-1}$, has full column rank;

$\lambda^{(i)}$ is known, but not necessarily unique.

Let $\lambda_k$ and $\delta_k$ consist of those components $\lambda_i$ and $\delta_i$ of $\lambda$ and $\delta$, respectively, for which $i \in J_k$. For $i \notin J_k$, $\lambda_i^{(j)} = 0, j = 0, 1, 2, \cdots$. The equation for the $(k+2)$-th time derivative of $q$ is by (2.5a) and the definition of $G_k$

$$ (5.10) \qquad q^{(k+2)} = (M^{-1}f)^{(k)} + (M^{-1}G_k\lambda_k)^{(k)}. $$

Since $(G_k^T \dot{q})_j^{(i)} = \delta_j^{(i-1)} = 0, j \in J_k$, for $i=1, 2, \cdots, k$, and $G_k^T \dot{q} = 0$, it follows from Lemma 5.3 that

$$ (5.11) \qquad \delta_k^{(k)} = (G_k^T \dot{q})^{(k+1)} \in \mathbb{R} (G_k^T). $$

$G$ is analytic and all derivatives $q^{(i)}$ in (5.11) are known except for $q^{(k+2)}$. Split $G_k$ into two parts $G_{k1}$ and $G_{k2}$ such that after permutation of the columns

$$ G_k = (G_{k1}, G_{k2}), $$

and partition $\lambda_k$ and $\delta_k$ accordingly. It follows from (5.10) and (5.11) that there exists a $b_k$ such that the system of equations to be solved for $\lambda_k^{(k)}$ and $\delta_k^{(k)}$ can be written in the partitioned form

$$ (5.12) \qquad \begin{pmatrix} G_{k1}^T M^{-1} G_{k1} & G_{k1}^T M^{-1} G_{k2} \\ G_{k2}^T M^{-1} G_{k1} & G_{k2}^T M^{-1} G_{k2} \end{pmatrix} \begin{pmatrix} \lambda_{k1}^{(k)} \\ \lambda_{k2}^{(k)} \end{pmatrix} + \begin{pmatrix} G_{k1}^T b_k \\ G_{k2}^T b_k \end{pmatrix} = \begin{pmatrix} 0 \\ \delta_{k2}^{(k)} \end{pmatrix}, \\ \lambda_{k2}^{(k)} \ge 0, \quad \delta_{k2}^{(k)} \ge 0, \quad (\lambda_{k2}^{(k)})^T \delta_{k2}^{(k)} = 0. $$

$\delta_{k1}^{(k)}$ must be equal to zero since $\lambda_{k1}(t) > 0$ and $\delta_{k1}(t) \equiv 0$ at least for a short time interval $T = (t_0, t_0 + \Delta t')$. The sign of $\lambda_{k1i}^{(k)}$ is of no importance here, since there is a $p \le k-1$ such that $\lambda_{k1i}^{(p)} > 0$ and hence in $T$,

$$ \lambda_{k1i}(t) = \lambda_{k1i}^{(p)}(t_0)(t-t_0)^p + O((t-t_0)^{p+1}) > 0. $$

Furthermore, $\lambda_{k2}^{(k)}$ and $\delta_{k2}^{(k)}$ must be complementary since $\lambda_{k2}^{(i)} = \delta_{k2}^{(i)} = 0, i = 0, 1, \cdots, k-1$.

By Lemma 5.2, (5.12) has a solution where $\delta_k^{(k)}$ is unique. There is a matrix $G'_{k2}$ such that

$$ (5.13) \qquad G_{k+1,1} = (G_{k1}, G'_{k2}) $$

has full column rank.

The principle of construction of $C$ in (3.3) and Lemma A imply that $C$ and $C^{-1}$ are bounded and analytic. Consider the equation

$$ C^T \ddot{q} = C^{-1} f + C^{-1} G \lambda $$

derived from (2.5a) and its time derivatives

$$ (5.14) \qquad (C^T \ddot{q})^{(k)} = (C^{-1} f)^{(k)} + (C^{-1} G \lambda)^{(k)}. $$

From (5.11) and the uniqueness of $\delta_k^{(k)}$, the conclusion is that $G_k^T q^{(k+2)} = c_k \in \mathbb{R} (G_k^T)$,
where $c_k$ is unique. Suppose that there are two different $q_1^{(k+2)}$, $q_2^{(k+2)}$ and $q_2^{(k+2)}$. The difference $C^T(q_1^{(k+2)} - q_2^{(k+2)})$ must be an element of $N(G_k^T C^{-T})$ since $c_k$ is unique. From the uniqueness of $q^{(i)}$, $i=0, 1, 2, \dots, k+1$, follows that

$$ (5.15) \qquad C^T(q_1^{(k+2)} - q_2^{(k+2)}) = (C^{-1}G_k(\lambda_1 - \lambda_2))^{(k)} $$

and that $(C^{-1}G_k(\lambda_1 - \lambda_2))^{(i)} = 0$ for $j < k$. Lemma 5.3 asserts that the right-hand side of (5.15) belongs to $R(C^{-1}G_k)$, which is the orthogonal space to $N(G_k^T C^{-T})$. There is a contradiction unless $q_1^{(k+2)} = q_2^{(k+2)}$ and $q^{(k+2)}$ is unique.

Determine $J_{\lambda k}$, $J_{\delta k}$ and $J_{k+1}$. By the definitions of these sets we have $J_{\lambda i} \subset J_{\lambda, i+1}$, $J_{\delta i} \subset J_{\delta, i+1}$ and $J_{i+1} \subset J_i$. $q^{(k+2)}$ is known and unique, and $\lambda^{(k)}$ is known. $G_{k+1,1}$ is specified in (5.13). The inductive step is completed.

Since the number of indices in $J$ is finite, eventually in a step the iterative procedure for determination of the active set can be terminated for three reasons:

1. $J_k \setminus J_{\lambda k} = \emptyset$.

2. All columns $g_{k2i}$ of $G_{k2}$ are linear combinations of the columns $g_{k1i}$ in $G_{k1}$, $g_{k2i} = G_{k1i}x_i$. Then by Lemma 5.3, for all $i$,

$$ (g_{k2i}^T \dot{q})^{(i)} = (x_i^T G_{k1}^T \dot{q})^{(i)} = 0, \quad j=0, 1, \dots $$

3. $\lambda_2^{(k)} = 0$.

Let $G' = G_{k1}$, which has full column rank. Since the rank of $G$ is constant in the neighborhood of $q_0$, so is the rank of $G'$. Similarly to (4.6) the equation for $\ddot{q}$ in $\bar{T} = [t_0, t_0 + \Delta t')$ is

$$ (5.16) \qquad M\ddot{q} = f - G'(G'^T M^{-1} G')^{-1} G'^T M^{-1} f - G'(G'^T M^{-1} G')^{-1} \dot{G}'^T \dot{q}. $$

$M^{-1}$, $G'$, $f$ and $\partial G'/\partial q$ are bounded and analytic in all $q_i$, and by Lemma A, $(G'^T M^{-1} G')^{-1}$ has the same properties. $f$ is analytic in $\dot{q}_i$ and $t$. Then there exists an analytic solution $q(t)$ satisfying (5.16) and the initial conditions in an open neighborhood of $t_0$, Coddington and Levinson (1955, Thm. 8.1).

$\lambda'(t)$ consists of those $\lambda_i$ for which $\lambda_i^{(0)} > 0$ or $\lambda_i^{(j)} > 0$, but $\lambda_i^{(k)} = 0$, $k = 0, 1, \dots, j-1$. $\lambda'(t)$ satisfies (cf., (4.5))

$$ \lambda'(t) = -(G'^T M^{-1} G')^{-1} (G'^T M^{-1} f + \dot{G}'^T \dot{q}), $$

and is analytic. After permutation of the columns in $G$,

$$ (5.17) \qquad \lambda^T = (\lambda'^T, 0, 0, \dots, 0) \in \mathbb{R}^n $$

is a feasible solution. The computed $\lambda_i^{(i)}$ and $\delta_i^{(i)}$ determine $\lambda^{(i)}$ and $\delta^{(i)}$.

If the iteration terminated due to the third criterion, $\delta_i(t)$ may become negative for some $i$ when using (5.16) for $q(t)$, viz., if $\delta_i^{(j)} = 0, j = 0, 1, \dots, p \ge k$ but $\delta_i^{(p+1)} < 0$. The iterative procedure must be continued. In that situation compute $\lambda_1^{(k)}$, $\lambda_1^{(k+1)}$, $\dots$ with $\lambda_2^{(k)} = \lambda_2^{(k+1)} = \dots = 0$ until $p+1$ is reached. Determine the unique $q^{(k+2)}$ and the sets $J_{\lambda k}$ and $J_{\delta k}$ and proceed as before. The iteration is eventually terminated, since the number of $\delta_i$ is finite, $\delta(t)$ is analytic and if $\delta_i(t) > 0$ in $T$ there must be an $i(j)$ such that $\delta_i^{(j)} > 0$ for each $j$.

When the iteration has terminated due to criterion 1, 2 or 3, then

$$ \lambda(t) = \sum_{i=0}^{\infty} \lambda^{(i)}(t_0)(t-t_0)^i \quad \text{and} \quad \delta(t) = \sum_{i=0}^{\infty} \delta^{(i)}(t_0)(t-t_0)^i $$

fulfil the complementarity condition $\lambda^T(t)\delta^T(t) = 0$, $\lambda(t) \ge 0$, $\delta(t) \ge 0$, in $\bar{T}$ for $\Delta t'$ sufficiently small.
The derivatives $q^{(i)}(t_0)$, $i=0, 1, 2, \dots$ are calculated uniquely using the iterative process. They will coincide with the derivatives at $t_0$ for the solution of (5.16). If no previously passive constraint becomes active in $T$, then the solution of (5.16) is also the unique solution of (2.5) in $\bar{T}$. However, since the number of passive constraints is finite, there is a $\Delta t > 0 < \Delta t \le \Delta t'$ such that the active set in $(t_0, t_0 + \Delta t)$ is $\{i | \delta_i^{(j)} = 0, j = 0, 1, 2, \dots\}$. A possible but not always unique $\lambda(t)$ is given by (5.17). If $G$ has full column rank, then $G'_{k2}$ in (5.13) is unique in each step and $G'$ as well as $\lambda$ are unique. $\square$

The proof can easily be modified to allow some constraints $i$ to be bilateral, i.e., there is no sign restriction on $\lambda_i$ and $\delta_i = 0, t \ge 0$.

The active set has been determined at a point $t_1 \in (t_0, t_0 + \Delta t)$ and in its neighborhood. If $G$ has not full column rank then we can find a $G'$, where $G'^T M^{-1} G'$ is nonsingular, yielding the same solution. As long as the active set and rank ($G$) are constant, the system of differential equations to be solved is (4.6) with $G = G'$ commencing from $t_1$ with $q(t_1)$ and $\dot{q}(t_1)$ as initial values. The uniqueness of the solution is then also given by Theorem 4.1.

It follows from the remark after Theorem 4.1 that $M$ and $G$ are often analytic in all their arguments $q_i$ separately. In many applications $f$ is at least piecewise analytic.

The rank of $G$ is not necessarily constant between the updates of the active set. If $G$ is slightly perturbed so that the rank is changed at $t_0$, the solution before $t_0$ may be completely different from that obtained after $t_0$ although the active set is constant. It is easy to construct simple two-dimensional examples which illustrate this fact. Fortunately, these problems do not seem to arise very frequently in practice. For most systems $G$ has the representation

$$G(q) = A(q)B(q), \quad A \in \mathbf{R}^{m \times R}, \quad B \in \mathbf{R}^{R \times n},$$

where the columns of $A$ are linearly independent and $rank(B) = R$. Such a $G$ satisfies the rank condition in the theorem.

In the neighborhood of an accumulation point $t_2$, $\Delta t$ is necessarily very small. In an arbitrarily small $\epsilon$-neighborhood of $t_2$, the solution can be constructed using the above procedure to the point $t = t_2 - \epsilon/2$. Then estimates for $\dot{q}(t)$ and $q(t)$ are given by (3.7)

$$ \begin{aligned} & \|C^T \dot{q}(t)\| \le \|C^T \dot{q}(t_2 - \epsilon/2)\| (1 + O(\epsilon)), \\ & \|C^T q(t)\| \le \|C^T q(t_2 - \epsilon/2)\| (1 + O(\epsilon)), \\ & |t - t_2| \le \epsilon. \end{aligned} $$

**6. Concluding remarks.** A solution to the linear complementarity problem (2.5b) is also a solution to the quadratic programming problem

$$(6.1) \qquad \min_{\lambda \ge 0} \frac{1}{2} \lambda^T G^T M^{-1} G \lambda + \lambda^T (G^T M^{-1} f + \dot{G}^T \dot{q})$$

and vice versa, see Cottle and Dantzig (1968). The dual problem of (6.1) is, Dorn (1960),

$$(6.2) \qquad \min_{\lambda} \frac{1}{2} \lambda^T G^T M^{-1} G \lambda, \quad G^T M^{-1} G \lambda + G^T M^{-1} f + \dot{G}^T \dot{q} \ge 0.$$

After inserting the expression for $G\lambda$ from (2.5a) into (6.2), the new minimization problem is

$$(6.3) \qquad \min_{\ddot{q}} \frac{1}{2}(\ddot{q} - M^{-1}f)^T M(\ddot{q} - M^{-1}f), \quad G^T \ddot{q} + \dot{G}^T \dot{q} \ge 0.$$
In (6.1–3) $q$ and $\dot{q}$ are fixed variables. Equation (6.3) is a generalization of the principle of least constraint by Gauss, Whittaker (1944), to rigid body problems with unilateral constraints.

The formulation best suited for numerical computation seems to be (6.1) or (6.3). In order to compute approximations to $q(t)$ and $\dot{q}(t)$ at discrete time points, the time-derivatives are replaced by difference approximations. Then there are reliable algorithms for solving the optimization problem. For further details see Lötstedt (1979b).

In some applications there are upper and lower bounds on the multipliers

$$ (6.4) \qquad \lambda^L \le \lambda_i \le \lambda^U. $$

The corresponding constraint function $g_i^T \dot{q}$ satisfies

$$ (6.5) \qquad \begin{array}{ll} \text{if } \lambda^L < \lambda_i < \lambda^U & \text{then } g_i^T \dot{q} = 0, \\ \text{if } \lambda_i = \lambda^L & \text{then } g_i^T \dot{q} \ge 0, \\ \text{if } \lambda_i = \lambda^U & \text{then } g_i^T \dot{q} \le 0. \end{array} $$

The problem treated in §§ 2–5 is a special case of (6.4) and (6.5), where $\lambda^L = 0$ and $\lambda^U = +\infty$. The Coulomb friction forces, Kilmister and Reeve (1966), fulfil the relations (6.4) and (6.5). $\lambda^L$ and $\lambda^U$ depend on the normal force $\lambda_N \ge 0$ associated to the friction force in the following way

$$ \lambda^L = -\mu\lambda_N, \qquad \lambda^U = \mu\lambda_N, $$

where $\mu > 0$ is the friction coefficient. The results concerning the discontinuities in § 2 are applicable also when (2.5b) is replaced by conditions like (6.4) and (6.5). When $\lambda^L$ and $\lambda^U$ are constants, theorems similar to those in § 4 and § 5 can be proved to cover also systems governed by (2.5a), (6.4) and (6.5). In the Coulomb friction case, however, there are systems where no solution $\lambda$ exists and if it exists $\ddot{q}$ is not necessarily unique. This fact was first pointed out by Painlevé (1895).

**Appendix.** The lemma is stated for reference purposes only. It can be proved by using Cramer's rule for inversion of matrices.

**LEMMA A.** Let $A(q) \in \mathbf{R}^{n \times n}$ be nonsingular for $\|q - q_*\| \le b$ and have one of the following properties:

I. A has $p \ge 0$ continuous derivatives in all arguments $q_i$, $i = 1, 2, \dots, m$.

II. A is bounded and Lipschitz continuous in $\|q - q_*\| \le b$.

III. A is bounded in $\|q - q_*\| \le b$ and analytic in all the arguments separately in $\|q - q_*\| < b$.

Then $A^{-1}$ also has the same property.

**Acknowledgment.** This paper is a generalization of a part of a Ph.D. thesis, Lötstedt (1979a). The author wishes to express his gratitude to his thesis adviser, Professor G. Dahlquist, for many constructive discussions on the subject.

## REFERENCES

A. BEN-ISRAEL AND T. N. E. GREVILLE (1974), *Generalized Inverses*, John Wiley, New York.

E. A. CODDINGTON AND N. LEVINSON (1955), *Theory of Ordinary Differential Equations*, McGraw-Hill, New York.

R. W. COTTLE (1968), *On a problem in linear inequalities*, J. London Math. Soc., 43, pp. 378–384.

— (1979), *Numerical methods for complementarity problems in engineering and applied science*, Computing Methods in Applied Sciences and Engineering, 1977, I., R. Glowinski and J. L. Lions, eds., Springer, Berlin, pp. 37–52.
R. W. COTTLE AND G. B. DANTZIG (1968), *Complementary pivot theory of mathematical programming*, Lin. Alg. and Appl., 1, pp. 103–125.

P. A. CUNDALL (1974), *Rational Design of Tunnel Support*, Dept. of Civil and Mining Eng., University of Minnesota, Minneapolis, MN.

G. DAHLQUIST AND Å. BJÖRCK (1974), *Numerical Methods*, Prentice-Hall, Englewood Cliffs, NJ.

W. S. DORN (1960), *Duality in quadratic programming*, Quart. Appl. Math., 18, pp. 155–162.

G. DUVAUT AND J. L. LIONS (1972), *Les inéquations en mécanique et en physique*, Dunod, Paris.

R. L. HUSTON AND C. E. PASSERELLO (1976), *Eliminating singularities in governing equations of mechanical systems*, Mech. Res. Comm., 3, pp. 361–365.

A. W. INGLETON (1966), *A problem in linear inequalities*, Proc. London Math. Soc. (3), 16, pp. 519–536.

J. J. KALKER (1977), *A survey of the mechanics of contact between solid bodies*, ZAMM, 57, pp. T3–T17.

S. KARAMARDIAN (1972), *The complementarity problem*, Math. Prog., 2, pp. 107–129.

C. W. KILMISTER AND J. E. REEVE (1966), *Rational Mechanics*, Longmans, London.

A. I. KING AND C. C. CHOU (1976), *Mathematical modelling, simulation and experimental testing of biomechanical systems crash response*, J. Biomech., 9, pp. 301–317.

P. LÖTSTEDT AND G. DAHLQUIST (1978), *Interactive simulation of the progressive collapse of a house*, Numerical Methods for Differential Equations and Simulation, A. W. Bennet and R. Vichnevetsky, eds., North-Holland, New York, pp. 135–138.

P. LÖTSTEDT (1979a), *Analysis of some difficulties encountered in the simulation of mechanical systems with constraints*, TRITA-NA-7914, Dept. Numerical Analysis, Royal Institute of Technology, Stockholm.

——— (1979b), *A numerical method for the simulation of mechanical systems with unilateral constraints*, TRITA-NA-7920, Dept. Numerical Analysis, Royal Institute of Technology, Stockholm.

K. MAGNUS (1978), ed., *Dynamics of Multibody Systems*, Springer, Berlin.

P. PAINLEVÉ (1895), *Sur les lois du frottement de glissement*, C.R. Acad. Sci., 121, pp. 112–115; ibid. (1905), 140, pp. 702–707; ibid. (1905), 141, pp. 401–405 and pp. 546–552.

B. PAUL (1975), *Analytical dynamics of mechanisms—A computer oriented overview*, Mech. Mach. Theory, 10, pp. 481–507.

J. PÉRÈS (1953), *Mécanique générale*, Masson et Cie, Paris.

E. T. WHITTAKER (1974), *A Treatise on the Analytical Dynamics of Particles and Rigid Bodies*, 4th edition, Dover, New York.

J. WITTENBURG (1977), *Dynamics of Systems of Rigid Bodies*, Teubner, Stuttgart.