# Equilibrium paths of mechanical systems with unilateral constraints I. Theory

M Schulz, S. Pellegrino

► To cite this version:

M Schulz, S. Pellegrino. Equilibrium paths of mechanical systems with unilateral constraints I. Theory. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2000, 456 (2001), pp.2223 - 2242. 10.1098/rspa.2000.0610. hal-01395595

HAL Id: hal-01395595

https://hal.science/hal-01395595

Submitted on 11 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
# Equilibrium paths of mechanical systems with unilateral constraints
## I. Theory

BY M. SCHULZ¹ AND S. PELLEGRINO²

¹Department FV/SLT2, Robert Bosch GmbH, PO Box 30 02 40,
70442 Stuttgart, Germany (Marcus.Schulz2@de.bosch.com)

²Department of Engineering, University of Cambridge, Trumpington Street,
Cambridge CB2 1PZ, UK (pellegrino@eng.cam.ac.uk)

The paper shows that the behaviour of mechanical systems subject to unilateral constraints differs from that of standard systems in subtle, and yet important, ways. Therefore, a proper theoretical formulation is required for simulating their behaviour. After showing that the equilibrium equations for a multibody system subject to unilateral constraints have the same form as the standard Kuhn-Tucker conditions in optimization theory, the first-order equilibrium equations are derived and their integration is discussed. At a general integration step, one has to distinguish between constraints that are strongly active, weakly active and inactive. Whereas strongly active constraints can be treated like bilateral constraints and inactive constraints can be neglected, weakly active constraints need to be constantly re-analysed to determine if they switch to a different state. The outcome is that, in addition to the well-known limit points and bifurcation points, a new type of limit point can exist, where the path is non-smooth and the first-order equilibrium equations—after elimination of any strongly active constraints—non-singular. Such points are called corner limit points. In analogy with common limit points, the degree of instability of the system changes by one at a corner limit point.

**Keywords:** multibody systems; unilateral constraint; equilibrium path; limit point; snap through; bifurcation

## 1. Introduction

A common problem in structural mechanics is the determination of the static equilibrium path of a structure for varying values of a certain parameter, which may be an external load, a displacement or the arc length of the path itself. The reason why it is useful to know the equilibrium path, and also which of its parts are stable or unstable, is that it characterizes the behaviour of the structure. Designers are, of course, interested in any special features of a structure that might cause it to suddenly collapse. Special points of the equilibrium path are very important. For example, in a shallow arch or dome, a catastrophic change of configuration occurs at the limit point, and the load-carrying capacity of a thin-walled cylindrical shell in compression suddenly decreases at the bifurcation point. Therefore, much work has
been done on the computation of the equilibrium paths of all kinds of structures, and
also on robust solution procedures, to handle special features that are encountered
in certain cases (see, for example, Crisfield (1991) for more details).

More recently, there has been a growing interest in mechanical systems that consist of rigid, or elastically deformable, elements, such as rods or shells, connected together by passive joints, such as hinges or sliders, driven by motors or other kinds of actuators. These are known as multibody systems and are covered by a now extensive literature (Garcia de Jalon & Bayo 1994; Shabana 1989). Multibody systems have a variety of applications, including robotic arms and automotive suspensions.

The particular type of system that we are interested in is a linkage of flexible
elements connected by mechanical joints that undergo large changes of configuration.
This is a model for a type of deployable structure that is packaged for storage and/or
transportation, and can be automatically deployed into its operational configuration
(Pellegrino & Guest 2000). Many such structures have a small mobility and hence can
be analysed by a purely kinematic approach (see, for example, Kumar & Pellegrino
1996). However, there is an increasing tendency to design systems with a mobility
of zero, whose flexible elements are able to store strain energy during folding and to
release it during deployment (Hayman et al. 1994).

To simulate such systems, it is useful to think in terms of an equilibrium path
in a configuration space that includes a set of independent coordinates augmented
with a parameter that controls the configuration of the system. Standard commercial
multibody packages are currently unable to handle the complex behaviour of these
structures. Recently, a general computational approach to systems with bilateral
constraints has been proposed by Crisfield (1997), and some simple examples have
been analysed by Cardona & Huespe (1998).

The research that is presented in this paper originated from the need to develop
a realistic simulation of a solid surface deployable antenna, recently developed in
the Deployable Structures Laboratory at Cambridge University (Guest & Pellegrino
1996). A physical model of the antenna had shown that the antenna follows different
paths during deployment and retraction; in Schulz & Pellegrino (1998) we described
this behaviour in detail. We showed, by analysing a simple multibody system, that
the observed behaviour can be explained by the existence of two limit points on the
equilibrium path; the first would cause a snap during deployment and the second a
snap during retraction. Thus it would be possible for a deployable structure to be
in a particular configuration at the beginning of a deployment-retraction cycle and
in a different configuration at the end of it. Further work on systems that model
more closely the real antenna has shown, however, that the two limit points that are
required to produce this behaviour exist only for a small range of system parame-
ters.

An important effect that was not included in our previous study and, as far as
we know, has not been previously considered in multibody systems, is the effect of
unilateral constraints on the equilibrium path. We have discovered that there are
some important differences between systems of this type and the better-known type
of systems subject to bilateral constraints. Thus the aim of this paper is to develop a
full understanding of the type of behaviour that can occur in systems with unilateral
constraints and, based on this understanding, to develop a proper formulation for
simulating such systems. Simple examples are used to illustrate key points, and a
full application is presented in the companion paper (Schulz & Pellegrino 2000).
The paper is laid out as follows. Section 2 derives the first-order equilibrium equations for a multibody system that is subject to bilateral constraints, by generalizing a formulation in terms of independent coordinates in Thompson & Hunt (1973). A way of switching the integration variable near a limit point is discussed.

Section 3 derives the first-order equilibrium equations for a conservative multibody system subject to unilateral constraints. In addition to requiring that the Lagrangian of the system be stationary in an equilibrium configuration, it is also required that the reactions of all unilateral constraints be either zero, if the constraint is not active, or non-negative, if the constraint is active. Mathematically, these conditions have the same form as the standard Kuhn-Tucker conditions in optimization theory.

It is important to distinguish between strongly active constraints, which, in a particular configuration, are active with non-zero reaction, and weakly active, whose reaction is zero. When integrating the first-order equilibrium equations, at a general step all strongly active constraints can be treated like bilateral constraints until they become weakly active or inactive, but weakly active constraints need to be re-analysed to determine if they switch to a different state.

Section 4 shows that, in addition to the well-known type of limit points, there is a new type where the path is non-smooth and the first-order equilibrium equations non-singular. Such points are named corner limit points; they are associated with the change of state of a weakly active constraint.

In § 5 it is shown that, in analogy with common limit points, at a corner limit point the degree of instability of a system also changes by one. A simple example is presented in § 6. Section 7 analyses the different cases that can be encountered when the first-order equilibrium equations become singular.

A discussion, including a brief comparison of bilaterally and unilaterally constrained systems, concludes the paper.

## 2. Systems with bilateral constraints

We consider a conservative multibody system which is described by a set of $n$ generalized position coordinates

$$ \varphi = [\varphi_1 \ \varphi_2 \ \dots \ \varphi_n]^T $$

and $m$ scleronomic bilateral constraints

$$ g_j(\varphi) = 0, \quad j = 1, \dots, m. \qquad (2.1) $$

The potential energy $U$ of the system is assumed to be a function of the generalized coordinates and a single control parameter $p$. It is assumed that the constraint functions, as well as the energy function, are $C^2$-continuous. A standard procedure for enforcing $m$ constraints is to introduce the Lagrangian function

$$ \mathcal{L}(\varphi, p, \lambda) = U(\varphi, p) - g^T \lambda, \qquad (2.2) $$

where $\lambda$ is a vector of Lagrange multipliers and

$$ g = [g_1 \ g_2 \ \dots \ g_m]^T. \qquad (2.3) $$
The equations of equilibrium for this system can be obtained by differentiating (2.2) with respect to $\varphi$ and $\lambda$, which gives

$$ \nabla_{\varphi}\mathcal{L} = \nabla_{\varphi}U - E\lambda = 0, \qquad (2.4) $$

$$ \nabla_{\lambda}\mathcal{L} = -g = 0. \qquad (2.5) $$

Here, $E$, of size $n \times m$, is a constraint matrix, defined by

$$ E = [\nabla g_1 \ \nabla g_2 \ \dots \ \nabla g_m]. \qquad (2.6) $$

It will be assumed that the gradients of the constraint functions are independent in all configurations of interest, and hence $E$ has full rank (regularity assumption). The term $E\lambda$ represents the generalized constraint reactions.

Our aim is to determine equilibrium configurations of the system for varying values of the control parameter $p$. Excluding bifurcational behaviour, for the moment, the tangent to the equilibrium path in the $\varphi$-$p$ space is uniquely defined at every point of the path. The equilibrium path can be expressed with the aid of a parameter $s$,

$$ \varphi = \varphi(s), \quad \lambda = \lambda(s), \quad p = p(s). \qquad (2.7) $$

A suitable choice for $s$ would be the arc length of the path itself or—away from limit points—the control parameter $p$. Note that the symbol $\varphi$ is being used to denote both a general configuration of the system and an equilibrium configuration.

Starting from a given equilibrium configuration ($s=0$), the equilibrium path can be traced by applying a continuation method, i.e. the parameter $s$ is incremented in finite steps and the nonlinear algebraic equations (2.4) and (2.5) are solved iteratively for the current value of $s$. Efficient techniques for doing this exist (Riks 1979; Crisfield 1980).

A different approach will be followed in this study. The idea is to define the equilibrium path in terms of a set of differential equations in $\varphi$ and $\lambda$ that are numerically integrated with a standard solver. This is conceptually similar to using a continuation algorithm that evaluates derivatives in order to obtain improved initial estimates in each step. These differential equations are the *first-order* equilibrium equations for the system (Thompson & Hunt 1973), obtained by differentiating (2.4)–(2.5) with respect to $s$,

$$ (\nabla_{\varphi}^2\mathcal{L})\dot{\varphi} + (\nabla_{\varphi}U')\dot{p} - E\dot{\lambda} = 0, \qquad (2.8) $$

$$ E^T \dot{\varphi} = 0, \qquad (2.9) $$

where the dot denotes total differentiation with respect to $s$, the prime denotes partial differentiation with respect to $p$, and the Hessian of the Lagrangian function

$$ \nabla_{\varphi}^2\mathcal{L} = \nabla_{\varphi}^2 U - \sum_{j=1}^{m} \lambda_j \nabla^2 g_j = \nabla_{\varphi}^2 U - [E_{,\varphi_1}\lambda \ E_{,\varphi_2}\lambda \ \dots] \qquad (2.10) $$

is not necessarily regular.

Equation (2.9) is equivalent to stating that $\dot{\varphi}$ is in the nullspace of the constraint matrix transposed. Hence, considering a matrix $T$ whose columns span the nullspace,

$$ \dot{\varphi} = Tx, \qquad (2.11) $$
where $x \in \mathbb{R}^{n-m}$. Substituting (2.11) into (2.8) and premultiplying by $T^T$ yields

$$Hx + T^T(\nabla_{\varphi} U')\dot{p} = 0, \qquad (2.12)$$

where

$$H = T^T(\nabla_{\varphi}^2 \mathcal{L})T \qquad (2.13)$$

is the reduced Hessian.

If the reduced Hessian is regular we can solve (2.12) for $x$ and substitute the result into (2.11),

$$\dot{\varphi} = -TH^{-1}T^T(\nabla_{\varphi} U')\dot{p}. \qquad (2.14)$$

Then the path derivative of $\lambda$ is obtained by pre-multiplying (2.8) by $E^T$ and solving the resulting equation,

$$\dot{\lambda} = (E^T E)^{-1} E^T \{(\nabla_{\varphi}^2 \mathcal{L})\dot{\varphi} + (\nabla_{\varphi} U')\dot{p}\}. \qquad (2.15)$$

Note that the matrix $E^T E$ is positive definite due to the regularity assumption.

If we choose the control parameter as the path parameter ($s = p$), equations (2.14) and (2.15) become

$$\frac{d\varphi}{dp} = -TH^{-1}T^T\nabla_{\varphi}U', \qquad (2.16)$$

$$\frac{d\lambda}{dp} = (E^T E)^{-1} E^T \left\{ (\nabla_{\varphi}^2 \mathcal{L}) \frac{d\varphi}{dp} + \nabla_{\varphi} U' \right\}. \qquad (2.17)$$

### (a) Singularities

In some equilibrium configurations the reduced Hessian $H$ may become singular. Let

$$W = [w_1 \quad \cdots \quad w_m]$$

be the matrix of the eigenvectors of $H$; the diagonalized form of $H$ is

$$W^T H W = \begin{bmatrix} d_{11} & & & \\ & d_{22} & & \\ & & \ddots & \\ & & & d_{mm} \end{bmatrix}. \qquad (2.18)$$

Introducing the transformation

$$x = Wy, \qquad (2.19)$$

the diagonalized form of (2.12) is

$$W^T H W y + W^T T^T \nabla_{\varphi} U' \dot{p} = 0. \qquad (2.20)$$

Assuming that only the first eigenvalue is zero, the first equation in (2.20), which corresponds to this eigenvalue, reads

$$0y_1 + w_1^T T^T \nabla_{\varphi} U' \dot{p} = 0. \qquad (2.21)$$
If $w_1^T T^T \nabla_\varphi U' = 0$, it is likely that the system is at a bifurcation point, in which case a higher-order analysis is required (Thompson & Hunt 1973). A corresponding equation has been derived, for the case of unconstrained systems, by Riks (1979). The post-critical behaviour can be determined by applying a branch switching procedure (Stein et al. 1990).

If $w_1^T T^T \nabla_\varphi U' \neq 0$, then (2.21) requires

$$ \dot{p} = 0, \qquad (2.22) $$

and the system is likely to be at a limit point (Thompson & Hunt 1973), in which case the solution of (2.20) is

$$ y = [y_1 \ 0 \ \dots \ 0]^T, \qquad (2.23) $$

and hence, from (2.19) and (2.11),

$$ x = w_1 y_1, \qquad (2.24) $$

$$ \dot{\varphi} = T w_1 y_1. \qquad (2.25) $$

Because $\dot{p} = 0$, the control parameter $p$ cannot be chosen as path parameter and hence an alternative choice is required. Any generalized coordinate such that $\dot{\varphi}_i \neq 0$ can be chosen and, considering (2.25), this requires that the $i$th component of $Tw_1$ be not equal to zero.

In practice, it is best to switch parameters before getting too close to a limit point, and a reasonable strategy is to choose the particular generalized coordinate that corresponds to the largest magnitude component of $Tw_1$. Letting $\varphi_i$ be the new path parameter, equation (2.12) is supplemented by

$$ \dot{\varphi}_i = 1. \qquad (2.26) $$

Hence (2.11) and (2.12) yield

$$ \begin{bmatrix} H & T^T \nabla_\varphi U' \\ t_i & 0 \end{bmatrix} \begin{bmatrix} x \\ \dot{p} \end{bmatrix} = \begin{bmatrix} 0 \\ 1 \end{bmatrix}, \qquad (2.27) $$

where $t_i$ denotes the $i$th row of $T$. After solving (2.27), we can compute the path derivatives of the generalized coordinates according to (2.11).

### 3. Systems with unilateral constraints

This section extends the formulation presented in §2 to systems with unilateral constraints, representing frictionless contacts between different parts of the system. Consider a system with $m$ unilateral constraints,

$$ g_j(\varphi) \geq 0, \quad j \in I, \qquad (3.1) $$

where $I = \{1, 2, \dots, m\}$ identifies the full set of constraints. For simplicity, it is assumed that there are no bilateral constraints and that the coordinates $\varphi$ are a set of independent coordinates when all unilateral constraints are inactive.

As in the previous section, the potential energy $U$ of this conservative system depends on the generalized coordinates $\varphi$ plus a single control parameter $p$. The same
assumptions on the continuity of $U$ and $g_j$ are made, and an analogous Lagrangian function is defined,

$$ \mathcal{L}(\varphi, p) = U - \sum_{j \in I} \lambda_j g_j. \qquad (3.2) $$

The differentiation of this function with respect to $\varphi$ is unchanged from §2, and yields an equation equivalent to (2.4). Equation (2.5) is now replaced by inequality constraints and, additionally, the Lagrange multipliers must be non-negative, and cannot be non-zero together with the corresponding constraint function (comple- mentarity condition). That is, a constraint is either inactive ($g_j > 0$) or it is active ($g_j = 0$) and its reaction is positive ($\lambda_j \ge 0$) (see Pfeiffer & Glocker (1996) for further details). Mathematically,

$$ \nabla_{\varphi} \mathcal{L} = \nabla_{\varphi} U - \sum_{j \in I} \lambda_j \nabla g_j = 0, \qquad (3.3) $$

$$ \frac{\partial \mathcal{L}}{\partial \lambda_j} = -g_j \le 0, \quad \lambda_j \ge 0, \quad \lambda_j g_j = 0, \quad j \in I. \qquad (3.4) $$

In optimization theory, equations (3.3), (3.4) are known as Kuhn-Tucker condi- tions (Fletcher 1987). It is well known that their validity is subject to a regularity condition (see below).

We distinguish between *active constraints*, for which $g_j = 0$, whose index set is

$$ I^* = \{j \in I \mid g_j(\varphi) = 0\}, \qquad (3.5) $$

and *inactive constraints*, for which $g_j > 0$. Active constraints can be further subdivided into *strongly active* ($\lambda_j > 0$), whose index set is

$$ I_+ = \{j \in I \mid g_j = 0, \lambda_j > 0\}, \qquad (3.6) $$

and *weakly active constraints* ($\lambda_j = 0$), whose set is simply $I^* \setminus I_+$. The regularity condition requires that the gradient vectors $\nabla g_j$ of all active constraints be linearly independent.

We are interested in tracing the equilibrium path of this system following an approach similar to §2. However, the path in the $\varphi$-$p$ space is, in general, non-smooth, due to the change of state of unilateral constraints from active to inactive, or vice versa. Therefore, whereas for systems with bilateral constraints the equilibrium path is smooth and the path derivatives are continuous, now we will have to consider left- and right-hand side derivatives, which may not coincide.

In analogy with §2, we have to differentiate (3.3) with respect to $s$. Due to the non-smoothness of the path, we replace ordinary derivatives with right-hand side derivatives with respect to $s$, denoted with a dot,

$$ (\nabla_{\varphi}^2 \mathcal{L}) \dot{\varphi} + (\nabla_{\varphi} U') \dot{p} - E_+ \dot{\lambda}_+ - E_0 \dot{\lambda}_0 = 0, \qquad (3.7) $$

$$ \nabla_{\varphi}^2 \mathcal{L} = \nabla_{\varphi}^2 U - \sum_{j \in I_+} \lambda_j \nabla^2 g_j, \qquad (3.8) $$

where the constraint matrices for the strongly and weakly active constraints, $E_+$ and $E_0$ respectively, and the corresponding column vectors of path derivatives of
Lagrange multipliers, are

$$E_+ = [\nabla g_j], \quad \dot{\lambda}_+ = \{\dot{\lambda}_j\}, \quad j \in I_+, \qquad (3.9)$$

$$E_0 = [\nabla g_j], \quad \dot{\lambda}_0 = \{\dot{\lambda}_j\}, \quad j \in I^* \setminus I_+. \qquad (3.10)$$

The reason why the last two terms in (3.7) include only the active constraints is because, for every point $\tilde{\varphi}$ of the equilibrium path, excluding end points, there exists a neighbourhood on the path where the strict inequality $g_j > 0$ holds for all constraints which are inactive at $\tilde{\varphi}$. In this neighbourhood, the corresponding Lagrange multipliers $\lambda_j$ are equal to zero, due to the complementarity condition $g_j\lambda_j = 0$, and thus their path derivatives $\dot{\lambda}_j$ are also zero at $\tilde{\varphi}$.

Next, we turn to (3.4). Since $\nabla_{\varphi}U$ and $\nabla g_j$ are continuous on the equilibrium path in the $\varphi$-p space, the regularity assumption and (3.3) imply that $\lambda_j$ is also continuous. For this reason, at every equilibrium point $\tilde{\varphi}$, there exists a neighbourhood on the equilibrium path—excluding end points—where the strict inequality $\lambda_j > 0$ holds for all constraints which are strongly active at $\tilde{\varphi}$. The corresponding constraint functions are equal to zero in this neighbourhood because of the complementarity condition, and hence also their path derivatives at $\tilde{\varphi}$. In conclusion,

$$\dot{g}_j = 0, \quad j \in I_+. \qquad (3.11)$$

The corresponding equation for systems with bilateral constraints would be identical. In analogy with §2, we replace (3.11) with

$$\dot{\varphi} = T_+x, \qquad (3.12)$$

where the columns of $T_+$ span the nullspace of $E_+^T$.

For weakly active constraints, the path derivatives of the constraint functions and of the corresponding Lagrange multipliers are greater than or equal to zero. Moreover, complementarity must hold for these path derivatives,

$$\dot{g}_j \ge 0, \quad \dot{\lambda}_j \ge 0, \quad \dot{\lambda}_j \dot{g}_j = 0, \quad j \in I^* \setminus I_+. \qquad (3.13)$$

Otherwise, the constraint functions and the Lagrange multipliers could become negative and there would exist points where the complementarity condition (3.4) would be violated.

At this point, we eliminate the strongly active constraints by introducing a reduced Hessian, as in §2,

$$H_+ = T_+^T (\nabla_{\varphi}^2 \mathcal{L}) T_+. \qquad (3.14)$$

Defining

$$\dot{g}_0 = \{\dot{g}_j\}, \quad j \in I^* \setminus I_+, \qquad (3.15)$$

we can write (3.7) and (3.13) as

$$H_+x + T_+^T (\nabla_{\varphi} U') \dot{p} - T_+^T E_0 \dot{\lambda}_0 = 0, \qquad (3.16)$$

$$\dot{g}_0 \ge 0, \quad \dot{\lambda}_0 \ge 0, \quad \dot{\lambda}_0^T \dot{g}_0 = 0, \qquad (3.17)$$

where $x$ and $\dot{g}_0$ are linked by

$$\dot{g}_0 = E_0^T T_+ x. \qquad (3.18)$$
The inequalities in (3.17) are to be interpreted in the sense that all elements of a vector are greater than or equal to zero.

Equations (3.16) and (3.17) are the first-order equilibrium equations of a multi-body system with unilateral constraints. Similar equilibrium equations were derived by Björkman (1992) for structures that come into contact with a single rigid obstacle. Björkman made use of the mathematical concept of B-differentiability. Our derivations are more straightforward, as no advanced mathematical concepts are required.

Consider the case of systems with bilateral constraints. The standard stability condition for a solution $\varphi$, $\lambda$ of (2.4), (2.5) is

$$z^T(\nabla_{\varphi}^2 \mathcal{L})z > 0 \qquad (3.19)$$

for any $z \neq 0$ in the nullspace of the constraint matrix transposed. This is a sufficient condition for uniqueness of the solution and regularity of the reduced Hessian. In the case of systems with unilateral constraints, the stability condition would have to be satisfied only for those elements $z \neq 0$ of the nullspace of the matrix of strongly active constraints transposed that also satisfy the condition

$$z^T \nabla g_j \geq 0 \quad \forall j \in I^* \setminus I_+. \qquad (3.20)$$

However, this will not guarantee the reduced Hessian $H_+$ to be positive definite, or, indeed, equations (3.16) and (3.17) to have a unique solution.

The situation can be illustrated with the aid of the example depicted in figure 1. The two generalized coordinates are subjected to two linear unilateral constraints. The potential energy $U$ is a quadratic positive semi-definite function of the generalized coordinates, whose contours are shown in the figure. For $s = 0$ (figure 1a), the intersection of the constraint functions corresponds to a stable equilibrium point and (3.19) holds for $z \neq 0$ in the nullspace of $E_+^T$ and satisfying (3.20). There is no other feasible equilibrium point. Both constraints are weakly active at the point of equilibrium.

Let us now consider a continuous transformation of the potential energy from $s = 0$ (figure 1a) to $s = \epsilon > 0$ (figure 1b). Obviously, all points on the segment AB correspond to equilibrium solutions and, thus, the path derivative of $\varphi$ at the equilibrium point of figure 1a is not unique. There are path derivatives with

(i) $\dot{g}_1 > 0, \dot{g}_2 > 0;$

(ii) $\dot{g}_1 > 0, \dot{g}_2 = \dot{\lambda}_2 = 0; \text{ and}$

(iii) $\dot{g}_1 = \dot{\lambda}_1 = 0, \dot{g}_2 > 0.$

If, instead, the continuous transformation from $s = 0$ (figure 1a) to $s = \epsilon > 0$ produces the result shown in figure 1c, the path derivative is unique and $\dot{\lambda}_0 > 0$. In the next sections we discuss the solution of the first-order equilibrium equations depending on the features of the reduced Hessian.

## 4. Regular $H_+$

If the reduced Hessian is non-singular, equation (3.16) can be solved for $x$ and the solution substituted into (3.18). Together with (3.17), this yields

$$\dot{g}_0 = G\dot{\lambda}_0 + \dot{p}\hat{g}, \qquad (4.1)$$

$$\dot{g}_0 \geq 0, \quad \dot{\lambda}_0 \geq 0, \quad \dot{\lambda}_0^T \dot{g}_0 = 0, \qquad (4.2)$$
<figure><img src="image_2.png" /></figure>

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

<figure><img src="image_4.png" /><figcaption>Figure 1. Solution of first-order equilibrium equations. (a) $s = 0$. (b) $s = \epsilon > 0$, the solution is non-unique. (c) $s = \epsilon > 0$, but the solution is unique.</figcaption></figure>

where

$$G = E_0^T T_+ H_+^{-1} T_+^T E_0, \quad (4.3)$$

$$\hat{g} = -E_0^T T_+ H_+^{-1} T_+^T \nabla_{\varphi} U'. \quad (4.4)$$
<figure><img src="image_5.png" /><figcaption>Figure 2. Solution of the LCP for positive definite $H_+$.</figcaption></figure>

This is a linear complementarity problem (LCP) in standard form, which has been dealt with by, for example, Murty (1988), Lemke (1970) and Cottle (1977). Contact problems, both of a static and dynamic nature, invariably lead to complementarity problems, in some cases of the linear type. Examples of LCPs in dynamic mechanical systems can be found in Pfeiffer & Glocker (1996) and Klarbring (1988); structural mechanics examples can be found in Björkman (1991).

In tracing the equilibrium path of a system with unilateral constraints, when one or more of these constraints are weakly active we have to determine in each step which of the weakly active constraints become inactive and which remain active, in order to update the index sets $I^*$ and $I_+$. In general, $\dot{\lambda}_j$ is discontinuous at a point where an inactive constraint becomes active, or vice versa. Such an event can induce the transformation of other weakly active constraints into strongly active or inactive.

For this reason, treating unilateral constraints as either bilateral if they are active, or altogether neglecting them when their Lagrange multiplier becomes negative, without solving the LCP can lead in some cases to non-convergence and errors.

After solving the LCP, the right-hand side path derivatives of $\varphi$ can be computed by substituting $\dot{\lambda}_0$ into (3.16), solving for $x$, and substituting into (3.12),

$$ \dot{\varphi} = T_+ H_+^{-1} T_+^T (E_0 \dot{\lambda}_0 - (\nabla_\varphi U') \dot{p}). \quad (4.5) $$

The right-hand side path derivatives of $\lambda_+$ are obtained by solving (3.7),

$$ \dot{\lambda}_+ = (E_+^T E_+)^{-1} E_+^T ((\nabla_\varphi^2 \mathcal{L}) \dot{\varphi} + (\nabla_\varphi U') \dot{p} - E_0 \dot{\lambda}_0). \quad (4.6) $$

Equations (4.5) and (4.6) are integrated. At each step, the constraint function of the inactive constraints and the Lagrange multipliers of the strongly active constraints are checked for zeros. In this way, the index sets $I^*$ and $I_+$ can be updated.

### (a) Positive definite $H_+$

Due to the regularity assumption, the matrix $G$ is also positive definite and, therefore, the LCP in (4.1) and (4.2) has a unique solution both for $\dot{p} < 0$ and for $\dot{p} > 0$ (Murty 1988). The two solutions that are obtained correspond to the two directions in which we can move on the path (see figure 2).
<figure><img src="image_6.png" /><figcaption>Figure 3. Solution of the LCP for one weakly active constraint and negative G.</figcaption></figure>

It can be shown that (4.1) and (4.2) represent the Kuhn–Tucker conditions for a solution of the associated quadratic programming problem

$$ \min_{\dot{\lambda}_0} (\frac{1}{2} \dot{\lambda}_0^T G \dot{\lambda}_0 + \hat{g}^T \dot{\lambda}_0), \quad \dot{\lambda}_0 \ge 0, \qquad (4.7) $$

which can be solved instead of the LCP.

### (b) Non-positive definite $H_+$

In this case, the uniqueness of the solution of the LCP is not guaranteed and hence bifurcational behaviour can occur, despite the fact that the reduced Hessian is not singular. This is an important difference from systems with bilateral constraints.

If $H_+$ is, for example, negative definite, G is also negative definite. In fact, G can be negative definite even in the case of an indefinite $H_+$.

In practical situations, often only one constraint is weakly active. In this case, if G is negative, the possible solutions of the LCP are illustrated in figure 3. The three lines plotted have the same slope G (less than zero). For $\hat{g} = 0$, the LCP has the same unique solution for $\dot{p} > 0$ and $\dot{p} < 0$, as $\dot{p}\hat{g} = 0$. The weakly active constraint remains weakly active and the left- and right-hand side path derivatives of $\varphi$ and $\lambda$ coincide (see (4.5) and (4.6)).

If $\hat{g} \neq 0$, the LCP has two solutions for $\dot{p}\hat{g} > 0$ and no solution for $\dot{p}\hat{g} < 0$. Thus, if there are two solutions for $\dot{p} < 0$ ($\dot{p} > 0$), there is no solution for $\dot{p} > 0$ ($\dot{p} < 0$). This means that the equilibrium path has reached a maximum (minimum) in the $\varphi$-$p$ space, which will be called a corner limit point because of the non-smoothness of the path at this point.

At a common limit point, the condition $\dot{p} = 0$ holds, but at a corner point it does not necessarily hold (see figure 4b). The reduced Hessian $H_+$ is singular at a common limit point, whereas at a corner point it can be regular but not positive definite.

Note that the equilibrium path may also have reached a minimum or maximum in the $\varphi$-$p$ space if $G = 0$ ($H_+$ regular but indefinite) and $\hat{g} \neq 0$. In this case, the LCP has a unique solution for $\dot{p} = 0$, corresponding to a horizontal tangent, and a unique solution for $\dot{p}\hat{g} > 0$, corresponding to a non-horizontal tangent. In order to
<figure><img src="image_7.png" /><figcaption>Figure 4. (a) End point (number of weakly active constraints greater than one). (b) Corner limit point. (c) Common limit point.</figcaption></figure>

find out if the path has indeed reached a maximum or minimum, we have to derive the second-order, or possibly even higher-order, equilibrium equations.

If the number of weakly active constraints is greater than one, it is still possible to find corner limit points. In addition, one can find mathematical examples where there is a unique solution of the LCP for $\dot{p} < 0$ ($\dot{p} > 0$) and no solution for $\dot{p} > 0$ ($\dot{p} < 0$). In such cases, the equilibrium path has reached an *end point* (figure 4a). This is not possible with a single active constraint, as in this case the plot in figure 3 shows that there are either two solutions or none.

## 5. Stability at corner limit points

The stability of a mechanical system that is subject only to bilateral or strongly active constraints can be analysed by considering the reduced Hessian $H_+$. The number of negative eigenvalues of this matrix is equal to the *degree of instability of the system* (Thompson & Hunt 1973).

At a common limit point, the degree of instability changes by one if a single stability coefficient changes sign, and an obvious question is whether the behaviour at a corner limit point is the same. Let us consider the case of a corner limit point $s = s_P$ where a single constraint is weakly active. It is assumed that, in advancing on the equilibrium path with increasing $s$, this constraint is transformed from strongly active into inactive. Therefore, the dimension of $H_+$ increases by one. Next, it will be shown that in this case the number of negative stability coefficients also increases by one, and hence the degree of instability increases.

At such a corner limit point, the basis of the nullspace of the matrix of strongly active constraints transposed is obtained by augmenting with a vector $t_0$ the basis of the nullspace of $[E_+ \ E_0]^T$,

$$ T_+ = [T \ t_0]. \qquad (5.1) $$

Note that $E_0$ is a column vector and that $t_0$ has the property

$$ t_0^T E_0 \neq 0, \qquad (5.2) $$

because of the regularity assumption. Next, we partition the reduced Hessian

$$ H_+ = \begin{bmatrix} H & T^T(\nabla_{\varphi}^2 \mathcal{L})t_0 \\ t_0^T(\nabla_{\varphi}^2 \mathcal{L})T & t_0^T(\nabla_{\varphi}^2 \mathcal{L})t_0 \end{bmatrix}, \qquad (5.3) $$
where

$$H = T^T(\nabla_{\varphi}^2 \mathcal{L})T \qquad (5.4)$$

can be interpreted as the reduced Hessian to the left of $s_P$. Since $H_+$ is regular, assuming that $H$ and $t_0^T(\nabla_{\varphi}^2 \mathcal{L})t_0$ are non-singular, the inverse matrix has the form

$$(H_+)^{-1} = \begin{bmatrix} \cdots & \cdots \\ \vdots & (H_+/H)^{-1} \end{bmatrix}, \qquad (5.5)$$

where

$$H_+/H = t_0^T(\nabla_{\varphi}^2 \mathcal{L})t_0 - t_0^T(\nabla_{\varphi}^2 \mathcal{L})TH^{-1}T^T(\nabla_{\varphi}^2 \mathcal{L})t_0 \qquad (5.6)$$

is the Schur complement of $H$ in $H_+$ (Cottle 1974). Substituting (5.5) into (4.3) and evaluating

$$T_+^T E_0 = \begin{bmatrix} 0 \\ t_0^T E_0 \end{bmatrix}, \qquad (5.7)$$

we arrive at

$$G = (t_0^T E_0)^2 (H_+/H)^{-1}. \qquad (5.8)$$

Since $G$ is negative at the corner limit point, the Schur complement must also be negative. Using the inertia formula (Cottle 1974)

$$\mathrm{In}(H_+) = \mathrm{In}(H) + \mathrm{In}(H_+/H), \qquad (5.9)$$

where the inertia In of a real symmetric matrix is the triple

[no. eigenvalues > 0, no. eigenvalues < 0, no. eigenvalues = 0],

we can conclude that the number of negative eigenvalues of $H_+$ increases by one at the corner limit point. Thus, if the equilibrium path is stable for $s = s_P - \epsilon$, it is unstable for $s = s_P + \epsilon$, where $\epsilon$ is a small positive number. This result agrees completely with the behaviour at a common limit point.

## 6. Example of corner limit point

A simple system exhibiting a corner limit point is shown in figure 5. This system consists of a rigid bar AB of length $L$, connected by a frictionless hinge A to a rigid foundation, and an extensional spring of stiffness $k$, whose undeformed length is zero. The spring is attached to the bar at B and to a vertical slider at D. A unilateral constraint on point B requires it to lie above a horizontal line through C.

The position coordinate is the angle $\varphi$ between the bar and the horizontal, and the control parameter $p$ is the displacement of the end of the spring, measured from the line AC and positive upwards.

In the initial configuration, shown in the figure, $\varphi = 0$ and $p = -L$. The bar is in contact with C and the unilateral constraint $g(\varphi) = \sin \varphi \ge 0$ is strongly active.

The equilibrium path is shown in figure 6 for initially increasing $p$. The initial configuration corresponds to $(0, -L)$. As $p$ increases from $-L$ to 0, there is no change
<figure><img src="image_8.png" /><figcaption>Figure 5. Example of system with corner limit point; in the initial configuration $p = -L$, $\varphi = 0$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 6. Equilibrium path of simple example.</figcaption></figure>

in $\varphi$. At $(0, 0)$, the constraint becomes weakly active, and the path derivatives of the gap, $\dot{g}_0$, and of the reaction, $\dot{\lambda}_0$, are found by solving the corresponding LCP,

$$ \dot{g}_0 = -(kL^2)^{-1} \dot{\lambda}_0 - \dot{p}L^{-1}, \qquad (6.1) $$

$$ \dot{g}_0 \ge 0, \quad \dot{\lambda}_0 \ge 0, \quad \dot{\lambda}_0 \dot{g}_0 = 0. \qquad (6.2) $$

Because the complementarity condition requires at least one of $\dot{\lambda}_0$, $\dot{g}_0$ to be zero, this LCP can be solved by trial and error. Thus, first we set $\dot{g}_0 = 0$ and find $\dot{\lambda}_0 = -kL\dot{p} \ge 0$, which requires $\dot{p} \le 0$. Second, we set $\dot{\lambda}_0 = 0$ and obtain $\dot{g}_0 = -\dot{p}L^{-1} \ge 0$, which also requires $\dot{p} \le 0$. So, there are two non-trivial solutions for $\dot{p} < 0$ and no solution for $\dot{p} > 0$, as shown in figure 6.†

So, from the corner point, the system can move in two different directions on the equilibrium path. If the system goes back towards the initial configuration, bar AB remains horizontal and in contact with the constraint at B; this part of the path is

<footnote>† Note that $\dot{p}$ is the right-hand side derivative, which is equal to the infinitesimal change of $p$ when *moving away* from the corner point.</footnote>
stable. Alternately, AB can rotate by an amount such that AB and BD are always parallel; this part of the path is unstable.

In practice, from O, the system will jump to a stable configuration with $\varphi = 180^{\circ}$, from which $\varphi$ can be increased to follow a stable path that is not shown in the figure.

## 7. Singular $H_+$

We introduce a transformation analogous to (2.19), where W contains the eigenvectors of $H_+$ at the singularity $s = s_0$, and diagonalize (3.16),

$$W^T H_+ W y + W^T T_+^T (\nabla_\varphi U') \dot{p} - W^T T_+^T E_0 \dot{\lambda}_0 = 0. \quad (7.1)$$

If the set of weakly active constraints is empty, the last term disappears and the problem becomes identical to one involving only bilateral constraints (see § 2 a).

If the set of weakly active constraints is not empty, the situation becomes more complex. Assuming that only one eigenvalue of $H_+$ is zero, the solution of the diagonalized equations is

$$c_1 \dot{p} - c_2^T \dot{\lambda}_0 = 0, \quad (7.2)$$

$$y_i = \frac{1}{d_{ii}} w_i^T T_+^T (E_0 \dot{\lambda}_0 - \dot{p} \nabla_\varphi U'), \quad i = 2, \dots, m_+, \quad (7.3)$$

where

$$c_1 = w_1^T T_+^T \nabla_\varphi U'|_{s=s_0}, \quad (7.4)$$

$$c_2 = E_0^T T_+ w_1|_{s=s_0}, \quad (7.5)$$

and $m_+$ is the number of constraints that are strongly active in the configuration $s = s_0$.

We have analysed the case in which the number of weakly active constraints is one, in which case $E_0$ is a column vector and $c_2$, $\dot{\lambda}_0$ and $\dot{g}_0$ are scalar quantities, and have found that bifurcational behaviour can occur only for

$$c_1 = 0. \quad (7.6)$$

This condition corresponds to (2.22) for systems with bilateral constraints, but recall that systems with unilateral constraints can show bifurcational behaviour also when $H_+$ is regular (see § 4 b).

Next we summarize the results of a detailed study of the various cases that can occur.

If $c_1 \neq 0$ and $c_2 \neq 0$, it is likely that we are at a corner limit point. The first-order equilibrium equations have two solutions, corresponding to two directions on the equilibrium path. For one of these directions, $\dot{p} = 0$ and a more detailed classification of the different cases that are possible (see figure 7) requires the introduction of

$$c_3 = W_+^T T_+^T (E_0 - (c_2/c_1) \nabla_\varphi U'), \quad (7.7)$$

$$c_4 = E_0^T W_+ c_3, \quad (7.8)$$

where

$$W_+ = \left[ \frac{w_2}{\sqrt{d_{22}}} \quad \frac{w_3}{\sqrt{d_{33}}} \quad \dots \quad \frac{w_{m_+}}{\sqrt{d_{m_+ m_+}}} \right]. \quad (7.9)$$

If $c_1 \neq 0$ and $c_2 = 0$, it is likely that we are at a common limit point, where the equilibrium path is smooth and the weakly active constraint does not change state.
<figure><img src="image_10.png" /></figure>

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

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

<figure><img src="image_13.png" /><figcaption>Figure 7. Different shapes of equilibrium path at a corner point, for $c_1 \neq 0$ and $c_2 \neq 0$. (a) $c_3 \neq 0, c_4 > 0$. (b) $c_3 \neq 0, c_4 < 0$. (c) $c_3 \neq 0, c_4 = 0$. (d) $c_3 = 0$.</figcaption></figure>

## 8. Discussion

The first-order equilibrium equations of a conservative multibody system with unilateral constraints that are functions of the coordinates, but not of the control parameter of the system, have been derived in a straightforward manner, without recourse to advanced mathematical concepts. If the reduced Hessian matrix of the system is regular, a standard linear complementarity problem (LCP) is obtained, whose dimension is equal to the number of weakly active constraints.

During the integration of the equilibrium path equations, when the index set of weakly active constraints is non-empty, we have to determine at every step of the integration which of the weakly active constraints become inactive and which remain active. This requires the LCP to be solved. Strongly active unilateral constraints can be treated in a straightforward way, like bilateral constraints.

For systems with bilateral constraints, a sufficient condition for stability can be formulated in terms of second-order derivatives. If this condition holds, the path derivatives are uniquely defined. A corresponding second-order condition exists for systems with unilateral constraints, but it does not ensure that the path derivatives are unique.

In contrast to bilaterally constrained systems, systems subject to unilateral constraints can be at a limit or bifurcation point, even if the reduced Hessian matrix is regular. Therefore, a new type of limit point, called a corner limit point, exists for such systems. At corner limit points the equilibrium path is non-smooth.

At an end point, the LCP has a unique solution when the control parameter changes in one direction, but no solution in the other direction. It has been shown that, if the constraint functions are independent of the control parameter, end points can exist only if at least two unilateral constraints are weakly active.
A simple example system exhibiting a corner limit point has been discussed. However, finding a good example of a system with constraints independent of the control parameter and which has an end point is more difficult. This difficulty can be understood by recalling that, under the assumptions made above, an end point requires at least two constraints to be weakly active, whereas a corner limit point requires only one.

Finally, note that a corner limit point may look exactly like an end point in a plot of the control parameter versus a particular coordinate that happens not to vary in the vicinity of the limit point. Such plots can be misleading, and care should be taken before reaching general conclusions. We believe that it is important that behind any simulation software for the analysis of multibody systems of the type described in this paper there should be a solid theoretical framework that can address the many subtle issues discussed in this paper.

Support from the European Commission, in the form of a TMR Marie Curie Research Training Grant for M.S., and from the Royal Academy of Engineering, in the form of a Foresight Award for S.P., is gratefully acknowledged.

## Nomenclature

<table><tr><td><i>E</i>, <i>E</i><sub>0</sub>, <i>E</i><sub>+</sub></td><td>constraint matrices of bilateral, weakly and strongly active unilateral constraints</td></tr><tr><td><i>G</i></td><td>matrix defined in equation (4.3)</td></tr><tr><td><i>H</i>, <i>H</i><sub>+</sub></td><td>reduced Hessian matrices</td></tr><tr><td><i>H</i><sub>+</sub>/<i>H</i></td><td>Schur complement of <i>H</i> in <i>H</i><sub>+</sub></td></tr><tr><td><i>I</i>, <i>I</i><sup>*</sup>, <i>I</i><sub>+</sub><sup>*</sup></td><td>index sets of unilateral constraints, and of active and strongly active constraints</td></tr><tr><td><i>L</i></td><td>Lagrangian function</td></tr><tr><td><i>T</i>, <i>T</i><sub>+</sub></td><td>matrices of vectors that span the nullspace of <i>E</i><sup>T</sup> or <i>E</i><sub>+</sub><sup>T</sup></td></tr><tr><td><i>U</i></td><td>potential energy</td></tr><tr><td><i>W</i></td><td>matrix of eigenvectors of reduced Hessian</td></tr><tr><td><i>W</i><sub>+</sub></td><td>matrix defined in equation (7.9)</td></tr><tr><td><i>c</i><sub>1</sub>, <i>c</i><sub>2</sub>, <i>c</i><sub>3</sub>, <i>c</i><sub>4</sub></td><td>coefficients</td></tr><tr><td><i>d</i><sub>ii</sub></td><td>elements of diagonalized reduced Hessian</td></tr><tr><td><i>g</i><sub>j</sub>, <i>g</i>, <i>g</i><sub>0</sub></td><td>constraint function, column vector of constraint functions, corresponding column vector for weakly active constraints</td></tr><tr><td><i>ĝ</i></td><td>vector defined in equation (4.4)</td></tr><tr><td><i>m</i></td><td>number of bilateral or unilateral constraints</td></tr><tr><td><i>m</i><sub>+</sub></td><td>number of unilateral constraints that are strongly active</td></tr><tr><td><i>n</i></td><td>number of generalized coordinates</td></tr><tr><td><i>p</i></td><td>control parameter</td></tr><tr><td><i>s</i></td><td>path parameter</td></tr><tr><td><i>t</i><sub>i</sub></td><td><i>i</i>th row of <i>T</i></td></tr><tr><td><i>t</i><sub>0</sub></td><td>column vector defined in equation (5.1)</td></tr><tr><td><i>w</i><sub><i>i</i></sub></td><td><i>i</i>th column of <i>W</i></td></tr></table>
$x, y, z$ vectors

$\varphi_j, \varphi$ generalized coordinates, column vector of generalized coordinates

$\lambda_j, \lambda, \lambda_0, \lambda_+$ Lagrange multipliers, column vector of Lagrange multipliers, corresponding column vectors for weakly (subscript 0) and strongly active (subscript +) constraints

## References

Björkman, G. 1991 The solution of large displacement frictionless contact problems using a sequence of linear complementarity problems. *Int. J. Numer. Meth. Engng* **31**, 1553–1566.

Björkman, G. 1992 Path following and critical points for contact problems. *Computational Mech.* **10**, 231–246.

Cardona, A. & Huespe, A. 1998 Continuation methods for tracing the equilibrium path in flexible mechanism analysis. *Engng Computations* **15**, 190–220.

Cottle, R. W. 1974 Manifestations of the Schur complement. *Linear Algebra Applications* **8**, 189–211.

Cottle, R. W. 1977 Numerical methods for complementarity problems in engineering and applied science. Technical report SOL 77-24, Department of Operations Research, Stanford University, October.

Crisfield, M. A. 1980 A fast incremental/iterative solution procedure that handles snap-through. *Comp. Struct.* **13**, 55–62.

Crisfield, M. A. 1991 *Non-linear finite element analysis of solids and structures*, vol. 1, *Essentials*. Wiley.

Crisfield, M. A. 1997 *Non-linear finite element analysis of solids and structures*, vol. 2, *Advanced topics*. Wiley.

Fletcher, R. 1987 *Practical methods of optimization*, 2nd edn. Wiley.

Garcia de Jalon, J. & Bayo, E. 1994 *Kinematic and dynamic simulation of multibody systems*. Springer.

Guest, S. D. & Pellegrino, S. 1996 A new concept for solid surface deployable antennas. *Acta Astronautica* **38**, 103–113.

Hayman, G. J., Hedgepeth, J. M. & Park, K. C. 1994 Design freedoms of articulating astromast and their optimization for improved performance. In *Proc. 35th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conf.*, Hilton Head, SC, 18–20 April 1994, pp. 811–817.

Klarbring, A. 1988 A mathematical programming approach to three-dimensional contact problems with friction. *Computer Meth. Appl. Mech. Engng* **58**, pp. 349–384.

Kumar, P. & Pellegrino, S. 1996 Deployment and retraction of a cable-driven rigid panel solar array. *J. Spacecraft Rockets* **33**, 836–842.

Lemke, C. E. 1970 Recent results in complementarity problems. In *Nonlinear programming* (ed. J. B. Rosen and others), pp. 349–384. Academic.

Murty, K. G. 1988 *Linear complementarity, linear and nonlinear programming*. Berlin: Heldermann.

Pellegrino, S. & Guest, S. D. 2000 Deployable structures: theory and applications. In *Proc. IUTAM-IASS Symp., Cambridge, September 1998*. Dordrecht: Kluwer.

Pfeiffer, F. & Glocker, Ch. 1996 *Multibody dynamics with unilateral contact*. Wiley.

Riks, E. 1979 An incremental approach to the solution of snapping and buckling problems. *Int. J. Solids Struct.* **15**, 529–551.
Schulz, M. & Pellegrino, S. 1998 'Multi-path' motion of deployable structures. In *Proc. 4th World Cong. on Computational Mechanics, Buenos Aires, 29 June-2 July* (ed. S. R. Idelsohn, E. Onate & E. N. Dvorkin).

Schulz, M. & Pellegrino, S. 2000 Equilibrium paths of mechanical systems with unilateral constraints. II. Deployable reflector. *Proc. R. Soc. Lond. A 456*, 2243–2262. (Following paper.)

Shabana, A. A. 1989 *Dynamics of multibody systems*. Wiley.

Stein, E., Wagner, W. & Wriggers, P. 1990 Nonlinear stability-analysis of shell and contact-problems including branch-switching. *Computational Mech. 5*, 428–446.

Thompson, J. M. T. & Hunt, G. W. 1973 *A general theory of elastic stability*. Wiley.