# Equations of motion for general constrained systems in Lagrangian mechanics

Firdaus Udwadia, Aaron Schutte

► To cite this version:

Firdaus Udwadia, Aaron Schutte. Equations of motion for general constrained systems in Lagrangian mechanics. Acta Mechanica, 2010, 213 (1-2), pp.111 - 129. 10.1007/s00707-009-0272-2. hal-01395969

HAL Id: hal-01395969

https://hal.science/hal-01395969

Submitted on 13 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
# Equations of motion for general constrained systems in Lagrangian mechanics

F. E. Udwadia
Department of Aerospace and Mechanical Engineering, Civil Engineering, Mathematics, Systems Architecture Engineering, and Information and Operations Management, University of Southern California, 430K Olin Hall, Los Angeles, CA 90089-1453, USA
E-mail: fudwadia@usc.edu

A. D. Schutte
Department of Aerospace and Mechanical Engineering, University of Southern California, Los Angeles, CA 90089-1453, USA
E-mail: schutte@usc.edu

**Abstract** This paper develops a new, simple, explicit equation of motion for general constrained mechanical systems that may have positive semi-definite mass matrices. This is done through the creation of an auxiliary mechanical system (derived from the actual system) that has a positive definite mass matrix and is subjected to the same set of constraints as the actual system. The acceleration of the actual system and the constraint force acting on it are then directly provided in closed form by the acceleration and the constraint force acting on the auxiliary system, which thus gives the equation of motion of the actual system. The results provide deeper insights into the fundamental character of constrained motion in general mechanical systems. The use of this new equation is illustrated through its application to the important and practical problem of finding the equation of motion for the rotational dynamics of a rigid body in terms of quaternions. This leads to a form for the equation describing rotational dynamics that has hereto been unavailable.

## 1 Introduction

Obtaining the explicit equation of motion for constrained mechanical systems with singular mass matrices has been a source of major difficulty in classical mechanics for some time now. The main difficulty in arriving at such an equation is the fact that the mass matrix does not have an inverse. Thus, the equations for the system's acceleration such as the explicit fundamental equation of motion for constrained systems developed by Udwadia and Kalaba [1] cannot be used. Recently, however, Udwadia and Phohomsiri [2] considered general mechanical systems described by Newtonian and/or Lagrangian mechanics, and derived an explicit equation of motion for such systems with singular mass matrices. Their equation differs in structure and form from the so-called fundamental equation (Udwadia and Kalaba [1]). They found that under certain restrictions on the structure of the mass matrix and the structure of the constraints, the equation of motion of the constrained system becomes unique.

In this paper, we show that under the conditions stipulated by Udwadia and Phohomsiri [2] for when the equation of motion of a constrained mechanical system becomes unique—a circumstance that must necessarily arise when one is modeling real-life physical systems in classical mechanics because the observed accelerations are always unique—a simpler set of equations can be obtained that again have the same form as the
fundamental equation. These equations are valid for systems whose mass matrices can be singular (positive semi-definite) and/or positive definite.

The basic idea is to consider an auxiliary mechanical system that has a positive definite mass matrix (instead of the actual system that has the singular mass matrix), which is subjected to the identical constraint conditions as the actual system. The auxiliary constrained system is then devised to directly provide at each instant of time the acceleration of the actual system and also the constraint force acting on the actual system.

While some may argue that systems with singular mass matrices are rare in classical mechanics, it should be pointed out that this is so only because of our predisposition, from years of scientific acculturation, in using the minimum number of coordinates to describe a constrained mechanical system. Often, for convenience, especially in modeling complex multibody systems, it is useful to use redundant coordinates, and in such situations singular mass matrices can and do arise when describing the unconstrained system. When more than the minimum number of coordinates is used in describing the motion of constrained mechanical systems, the coordinates are evidently related to one another, and systems like these can have singular mass matrices. For instance, this happens when obtaining the equations of rotational motion of a rigid body whose orientation is parameterized by unit quaternions so as to accommodate large angle rotations. Here, an additional coordinate to describe the orientation of a rigid body beyond the minimum number required provides the convenience of having no singularities, a malady that the Euler angles suffer from; but this, as we shall see, leads to singular mass matrices. An illustration of the use of our new equation is provided to this important application area of rotational dynamics, where we use a Lagrangian approach to directly and simply obtain a hereto unavailable form of the equations of rotational motion of a rigid body in terms of quaternions.

## 2 System description and general constraints

We shall envision the description of a constrained mechanical system S in a two-step process. First, we imagine the system to be unconstrained and prescribe the motion of this unconstrained system; then, we impose on this unconstrained system a set of constraints to yield the desired constrained system. In accordance with this idea, consider an unconstrained system whose configuration is described by the n-vector $q$ and whose equation of motion is given by

$$M(q, t)\ddot{q} = Q(q, \dot{q}, t), \quad q(t=0) = q_0, \quad \dot{q}(t=0) = \dot{q}_0, \qquad (1)$$

where the $n$ by $n$ matrix $M(q, t)$ is positive semi-definite ($M(q, t) \ge 0$) at each instant of time, and the $n$-vector $Q(q, \dot{q}, t)$, which is called the ‘given’ force, is a known function of its arguments. The dots in the above equation denote derivatives, and $t$ denotes time. By unconstrained we mean that the $n$ components of the initial velocity of the system, $\dot{q}_0$, can be independently assigned, and the acceleration $\ddot{q}$ in relation (1) refers to the acceleration of the unconstrained system.

We now impose on this unconstrained mechanical system a set of $m$ sufficiently smooth constraints given by

$$\varphi(q, t) = 0, \qquad (2)$$

and

$$\psi(q, \dot{q}, t) = 0, \qquad (3)$$

where $\varphi$ and $\psi$ are $h$- and $s$-vectors, respectively, and $m = h + s$. We shall assume that the initial conditions satisfy the constraint Eqs. (2) and (3) so that $\varphi(q_0, 0) = \psi(q_0, \dot{q}_0, 0) = 0$. Therefore, for the constrained system, the components of the $n$-vectors $q_0$ and $\dot{q}_0$ cannot all be independently assigned.

These $m$ constraint equations can be nonlinear functions of their arguments, and they are not required to be independent of one another. They can be differentiated with respect to time to obtain the relation

$$A(q, \dot{q}, t)\ddot{q} = b(q, \dot{q}, t), \qquad (4)$$

where $A(q, \dot{q}, t)$ is an $m$ by $n$ matrix of rank $r$. These are the constraints that the acceleration $\ddot{q}$ of the constrained system S is required to satisfy.

In the presence of these constraints, the equation of motion for the constrained system S becomes

$$M(q, t)\ddot{q} = Q(q, \dot{q}, t) + Q^c(q, \dot{q}, t), \qquad (5)$$
where the n-vector $Q^c(q, \dot{q}, t)$ is the additional force of constraint that is created by virtue of the presence of the constraints (2) and (3). We shall permit this force to be non-ideal. Denoting a virtual displacement at time $t$ as any n-vector $v(t) \neq 0$ that belongs to the null space of the matrix $A$, such a non-ideal force of constraint can do work under virtual displacements. The dynamical nature of this constraint force is prescribed for a given mechanical system by the specification of a sufficiently smooth n-vector, $C(q, \dot{q}, t)$, such that the work done by the force of constraint under virtual displacements is given by

$$v(t)^T Q^c(q, \dot{q}, t) = v(t)^T C(q, \dot{q}, t). \qquad (6)$$

Specification of the kinematical constraints (2) and (3) [or alternatively, Eq. (4)] along with the dynamical specification of the vector $C(q, \dot{q}, t)$ constitutes a description of the constraint conditions imposed on our system $S$. At those instants of time when $C(q, \dot{q}, t)$ is zero, the work done by the force of constraint under virtual displacements is zero, and the system satisfies D'Alembert's principle at those instants of time. When $C(q, \dot{q}, t) = 0$ for all time, we say that the system satisfies D'Alembert's principle, and the force of constraint is said to be ideal. Most of analytical dynamics deals with ideal constraints (see References [3,4]).

From here on, we shall suppress the arguments of the various quantities unless needed for clarity.

When the mass matrix $M$ in Eq. (1) is positive definite, the equation of motion of the constrained system $S$ is explicitly given by the so-called fundamental equation (Udwadia and Kalaba [3])

$$\ddot{q} = a + M^{-1/2} B^{+} (b - Aa) + M^{-1/2} (I - B^{+} B) M^{-1/2} C, \qquad (7)$$

where $a = M^{-1}Q$, $B = AM^{-1/2}$, and the superscript ‘+’ denotes the Moore-Penrose (MP) inverse of the matrix $B$. The above equation is valid for systems that may or may not satisfy D'Alembert's principle at each instant of time. However, when the mass matrix is positive semi-definite ($M \ge 0$), the matrix $M^{-1/2}$ may not exist, and Eq. 7 becomes inapplicable, in general.

As shown by Udwadia and Phohomsiri [2], when the matrix $M \ge 0$ in Eq. (1), the constrained equation of motion need not, in general, be unique. They prove that for it to be unique, a necessary and sufficient condition is that the $n$ by $n+m$ matrix $\hat{M}^T = [M|A^T]$ has (full) rank $n$. Since in classical mechanics we know from practical observation that the equation of motion is unique—that is, the acceleration of a given physical system that starts with a given set of initial conditions is unique at each subsequent instant of time—satisfaction of this requirement can then be taken to be an indicator of whether a given (constrained) physical system under consideration has been properly mathematically described to adequately model it.

Under this assumption of full rank of $\hat{M}$, Udwadia and Phohomsiri [2] provide the explicit (unique) equation of motion of the constrained system as

$$\ddot{q} = \begin{bmatrix} (I - A^+ A)M \\ A \end{bmatrix} + \begin{bmatrix} Q + C \\ b \end{bmatrix} := \bar{M} + \begin{bmatrix} Q + C \\ b \end{bmatrix}, \qquad (8)$$

which is valid when the $n$ by $n$ matrix $M \ge 0$. Here, the superscript ‘+’ again denotes the Moore-Penrose (MP) inverse of the matrix. The form of Eq. (8) when $M$ is positive semi-definite is therefore strikingly different from the form of the fundamental Eq. (7) obtained when $M$ is positive definite. A unified equation of motion that is applicable to both these situations is hereto unavailable, and this is the intent of this paper.

More specifically, in this paper we shall assume throughout that $M \ge 0$ and that the physical mechanical system is appropriately mathematically modeled so that the matrix $\hat{M}$ has full rank; therefore, the equation of motion of the constrained mechanical system is assumed to be unique, as demanded by physical observation. We shall obtain an alternative explicit equation of motion for such constrained systems that is valid when $M \ge 0$ and that has the same form as the aforementioned fundamental equation of motion given by relation (7). We provide a unified way of obtaining the explicit acceleration of a mechanical system irrespective of whether $M \ge 0$ or $M > 0$. More importantly, the new equation developed herein leads to new and somewhat unexpected physical insights into the nature of constrained motion.

An application of this equation is illustrated to the important field of rotational dynamics. This in turn leads to a new form of the equation of motion for the rotational dynamics of a rigid body.

### 3 Equation of motion for systems with positive semi-definite and positive definite mass matrices

We begin by proving the following lemmas.
**Lemma 1** Let $M \ge 0$, and denote the Moore-Penrose (MP) inverse of an $m$ by $n$ matrix $A$ by $A^+$. The $n$ by $n+m$ matrix $\hat{M}^T = [M|A^T]$ has full rank, if and only if the $n$ by $n$ augmented mass matrix $M_A = M + \alpha^2 A^+ A$ is positive definite for any real number $\alpha \ne 0$.

*Proof* (a) Assume first that $\hat{M}$ has full rank; we shall prove this implies that $M_A$ is a positive definite matrix. We first observe that the matrix $M_A$ is symmetric since $M$ is symmetric, and also because $(A^+ A)^T = A^+ A$ by the fourth MP condition [5].

Let the matrix $A$ have rank $r$, and consider the singular value decomposition of the $m$ by $n$ matrix $A = W \Lambda V^T$, where the $r$ by $r$ diagonal matrix $\Lambda$ contains the $r$ singular values of $A$. Then $A^T = V \Lambda W^T$, where the $n$ by $r$ matrix $V$ has rank $r$. Since $\text{rank}(A^T) = r = \text{rank}(V \Lambda W^T) = \text{rank}(V)$, we get

$$ \text{Col}(A^T) = \text{Col}(V \Lambda W^T) = \text{Col}(V) = \text{Col}(\alpha V), \qquad (9) $$

where $\text{Col}(X)$ denotes the column space of the matrix $X$. Hence,

$$ \text{Col}(\hat{M}^T) = \text{Col}([M|A^T]) = \text{Col}([M|\alpha V]). \qquad (10) $$

Furthermore, since $\text{Col}(M) = \text{Col}(M^{1/2})$, we have

$$ \text{Col}(\hat{M}^T) = \text{Col}([M^{1/2}|\alpha V]), \qquad (11) $$

and since we assume that $\text{rank}(\hat{M}) = n$, we find that

$$ n = \text{rank}(\hat{M}) = \text{rank}(\hat{M}^T) = \text{rank}([M^{1/2}|\alpha V]). \qquad (12) $$

Hence,

$$ M_A = M + \alpha^2 A^+ A = M + \alpha^2 V \Lambda^{-1} W^T W \Lambda V^T = M + \alpha^2 V V^T = [M^{1/2}|\alpha V] \begin{bmatrix} M^{1/2} \\ \alpha V^T \end{bmatrix}. \qquad (13) $$

From relation (12), we know that the rank of $[M^{1/2}|\alpha V]$ is $n$, so the rank of $M_A$ is $n$. Furthermore, the last equality in Eq. (13) points out that the $n$ by $n$ matrix $M_A$ must at least be positive semi-definite, and because its rank is $n$, our result is proved.

(b) Assume that $M_A$ is positive definite; we shall now prove that this implies that $\text{rank}(\hat{M}) = \text{rank}([M|A^T]) = n$. Here, we have

$$ \text{rank}([M|A^T]) = \text{rank}([M|A^T][M|A^T]^T) = \text{rank}(MM^T + A^T A). \qquad (14) $$

It suffices then to prove that $\text{rank}(MM^T + A^T A) = n$, or equivalently that

$$ \text{Kernel}(MM^T + A^T A) = \{0\}. \qquad (15) $$

Suppose now that $(MM^T + A^T A)x = 0$. Then,

$$ x^T(MM^T + A^T A)x = 0. \qquad (16) $$

We note that $MM^T$ and $A^T A$ are both positive semi-definite. Because of relation (16), we must then have

$$ x^T(MM^T x) = x^T A^T A x = 0, \qquad (17) $$

so that

$$ Mx = M^T x = Ax = 0, \qquad (18) $$

where the first equality follows because $M$ is symmetric. Since $Ax = 0 \Rightarrow A^+ Ax = 0$, we can use Eq. (18) to obtain

$$ (Mx + \alpha^2 A^+ Ax) = (M + \alpha^2 A^+ A)x = 0, \qquad (19) $$

which on pre-multiplication by $x^T$ gives

$$ x^T(M + \alpha^2 A^+ A)x = x^T M_A x = 0. \qquad (20) $$

From Eq. 20, it follows that $x = 0$ because $M_A$ is assumed to be positive definite. Hence, $\text{Kernel}(MM^T + A^T A) = \{0\}$, and so $\text{rank}([M|A^T]) = n$. □
**Lemma 2** Let A be an m by n matrix and denote $B = AM^{-1/2}$, where M is a positive definite matrix. Then,

$$ (I - B^{+}B)M^{-1/2}A^{+} = 0. \qquad (21) $$

*Proof* Taking the transpose of the left-hand side of Eq. (21) we get

$$ \begin{aligned} [(I - B^{+}B)M^{-1/2}A^{+}]^{\mathrm{T}} &= [(I - B^{+}B)M^{-1/2}A^{+}AA^{+}]^{\mathrm{T}} \\ &= (A^{+})^{\mathrm{T}}(A^{+}A)^{\mathrm{T}}M^{-1/2}(I - B^{+}B)^{\mathrm{T}} \\ &= (A^{+})^{\mathrm{T}}A^{+}(AM^{-1/2})(I - B^{+}B) \\ &= (A^{+})^{\mathrm{T}}A^{+}B(I - B^{+}B) = 0. \end{aligned} \qquad (22) $$

In the first equality above we have used the second MP condition [5], in the third equality we have used the fourth MP condition, and in the last equality we have used the first MP condition. Hence, the result. □

**Lemma 3** Let A be an m by n matrix and M a positive semi-definite matrix. Let $M_A = M + \alpha^2 A^+ A$. Then for any real number $\alpha$,

$$ (I - A^{+}A)M = (I - A^{+}A)M_A. \qquad (23) $$

*Proof* We have

$$ \begin{aligned} (I - A^{+}A)M_A &= (I - A^{+}A)(M + \alpha^2 A^{+}A) \\ &= \alpha^2(A^{+}A - [A^{+}AA^{+}]A) + (I - A^{+}A)M \\ &= \alpha^2(A^{+}A - A^{+}A) + (I - A^{+}A)M \\ &= (I - A^{+}A)M. \end{aligned} \qquad (24) $$

**Lemma 4** Using the notation established above, denote the $(m+n)$ by $n$ matrix

$$ \bar{M}_A = \begin{bmatrix} (I - A^{+}A)M_A \\ A \end{bmatrix}. \qquad (25) $$

Then,

$$ \bar{M}_A^+ \begin{bmatrix} (Q + A^+ z) + C \\ b \end{bmatrix} = \bar{M}_A^+ \begin{bmatrix} Q + C \\ b \end{bmatrix}, \qquad (26) $$

for any m-vector $z$.

*Proof*

$$ \begin{aligned} \bar{M}_A^+ \begin{bmatrix} (Q + A^+ z) + C \\ b \end{bmatrix} &= [\bar{M}_A^T \bar{M}_A]^+ \bar{M}_A^+ \begin{bmatrix} (Q + A^+ z) + C \\ b \end{bmatrix} \\ &= [\bar{M}_A^T \bar{M}_A]^+ [M_A(I - A^+ A)|A^T] \begin{bmatrix} (Q + A^+ z) + C \\ b \end{bmatrix} \\ &= [\bar{M}_A^T \bar{M}_A]^+ [M_A(I - A^+ A)|A^T] \begin{bmatrix} Q + C \\ b \end{bmatrix} = \bar{M}_A^+ \begin{bmatrix} Q + C \\ b \end{bmatrix}. \end{aligned} \qquad (27) $$

The third equality above follows because of the second MP condition [5]. □

**Lemma 5** Let $M_A := M + \alpha^2 A^+ A > 0$, M be symmetric, and $H := AM_A^{-1}A^T$. Then we have the following identities:

(a)

$$ \frac{\partial M_A^{-1}}{\partial \alpha} = -2\alpha M_A^{-1} A^+ A M_A^{-1}, \qquad (28) $$
(b)

$$ (I - A^T H^+ A M_A^{-1}) A^+ = 0, \qquad (29) $$

(c)

$$ H(I - H H^+) = 0, \qquad (30) $$

(d)

$$ (I - H^+ H) A M_A^{-1} A^+ = 0, \qquad (31) $$

(e)

$$ \frac{\partial H}{\partial \alpha} = -2\alpha A M_A^{-1} A^+ H, \qquad (32) $$

(f)

$$ \left(\frac{\partial H}{\partial \alpha}\right)^T (I - H H^+) = (I - H^+ H) \left(\frac{\partial H}{\partial \alpha}\right)^T = 0, \qquad (33) $$

(g)

$$ \frac{\partial H^+}{\partial \alpha} = 2\alpha H^+ A M_A^{-1} A^+ H H^+, \qquad (34) $$

(h)

$$ \frac{\partial M_A^{-1} A^T H^+}{\partial \alpha} = 0. \qquad (35) $$

*Proof* (a) Since

$$ \frac{\partial M_A^{-1}}{\partial \alpha} = -M_A^{-1} \frac{\partial M_A}{\partial \alpha} M_A^{-1}, \qquad (36) $$

the result follows from the definition of $M_A$.

(b) By Lemma 2, since $M_A$ is positive definite, we have $0 = (I - B_A^+ B_A) M_A^{-1/2} A^+$, where $B_A = A M_A^{-1/2}$. Premultiplying by $M_A^{1/2}$ we get

$$ \begin{aligned} 0 &= M_A^{1/2} (I - B_A^+ B_A) M_A^{-1/2} A^+ \\ &= A^+ - M_A^{1/2} B_A^+ B_A M_A^{-1/2} A^+ \\ &= [I - A^T (A M_A^{-1} A^T)^+ A M_A^{-1}] A^+ = [I - A^T H^+ A M_A^{-1}] A^+, \end{aligned} \qquad (37) $$

where we have used the relation $B_A^+ = M_A^{-1/2} A^T (A M_A^{-1} A^T)^+ $ to get the third equality.

(c) Since $H$ is symmetric,

$$ \begin{aligned} H(I - H H^+) &= H^T [I - (H H^+)^T] \\ &= H^T [I - (H^+)^T H^T] = 0, \end{aligned} \qquad (38) $$

where we have used the third MP condition to get the first equality, and the first MP condition to get the last [5].
(d) Using the second MP condition, namely $A^+ = A^+AA^+$, we get

$$ 
\begin{aligned}
(I - H^+H)AM_A^{-1}A^+ &= (I - H^+H)AM_A^{-1}A^+AA^+ \\
&= (I - H^+H)AM_A^{-1}(A^+A)^T A^+ \\
&= (I - H^+H)(AM_A^{-1}A^T)(A^+)^T A^+ \\
&= (I - H^+H)H(A^+)^T A^+ \\
&= [I - (H^+H)^T]H^T(A^+)^T A^+ \\
&= [I - H^T(H^+)^T]H^T(A^+)^T A^+ = 0,
\end{aligned} 
\quad (39) $$

where we have used the fourth MP condition in the second and fifth equalities and the first MP condition to get the last equality [5].

(e) From part (a) above, we have

$$ 
\begin{aligned}
\frac{\partial H}{\partial \alpha} &= A \frac{\partial M_A^{-1}}{\partial \alpha} A^T \\
&= -2\alpha A M_A^{-1} A^+ A M_A^{-1} A^T \\
&= -2\alpha A M_A^{-1} A^+ H.
\end{aligned} 
\quad (40) $$

(f) Using part (e) above, and noting that $H$ is symmetric, the left-hand side of Eq. (33) becomes

$$ 
\begin{aligned}
\left(\frac{\partial H}{\partial \alpha}\right)^T (I - H H^+) &= \left(\frac{\partial H}{\partial \alpha}\right) (I - H H^+) \\
&= -2\alpha A M_A^{-1} A^+ H (I - H H^+) = 0,
\end{aligned} 
\quad (41) $$

where the last equality follows from part (c) above. Similarly, we have

$$ 
\begin{aligned}
(I - H^+H) \left(\frac{\partial H}{\partial \alpha}\right)^T &= (I - H^+H) \frac{\partial H}{\partial \alpha} \\
&= -2\alpha(I - H^+H)AM_A^{-1}A^+H = 0,
\end{aligned} 
\quad (42) $$

where the last equality follows from part (d) above.

(g) To find the derivative of the MP inverse of $H$ with respect to $\alpha$, we use the formula [6,7]

$$ \frac{\partial H^+}{\partial \alpha} = -H^+ \frac{\partial H}{\partial \alpha} H^+ + (I - H^+H) \left(\frac{\partial H}{\partial \alpha}\right)^T (H^+)^T H^+ + H^+ (H^+)^T \left(\frac{\partial H}{\partial \alpha}\right)^T (I - H H^+). \quad (43) $$

From part (f) above, the second and third members on the right-hand side of (43) are zero, and using the result of part (e) in the first member, the result follows.

(h) Using the result in parts (a) and (g) above, we have

$$ 
\begin{aligned}
\frac{\partial M_A^{-1} A^T H^+}{\partial \alpha} &= \frac{\partial M_A^{-1}}{\partial \alpha} A^T H^+ + M_A^{-1} A^T \frac{\partial H^+}{\partial \alpha} \\
&= -2\alpha M_A^{-1} A^+ (A M_A^{-1} A^T) H^+ + 2\alpha M_A^{-1} A^T H^+ A M_A^{-1} A^+ H H^+ \\
&= -2\alpha M_A^{-1} [I - A^T H^+ A M_A^{-1}] A^+ H H^+ = 0.
\end{aligned} 
\quad (44) $$

The last equality follows from part (b) above.

### 3.1 Explicit equation for constrained acceleration

We are now ready to begin by developing explicit equations for the constrained acceleration of the system S that has a positive semi-definite mass matrix M.
**Result 1:** For any mechanical system S

1. whose unconstrained equation of motion is described by $M(q, t)\ddot{q} = Q(q, \dot{q}, t)$, where the $n$ by $n$ mass matrix $M$ is positive semi-definite,

2. and that is further subjected to the $m$ consistent constraints described by Eqs. (2) and (3) [or alternatively, by Eq. (4)], which are satisfied by the initial conditions $q_0$ and $\dot{q}_0$,

3. such that the $n$ by $n+m$ matrix $\dot{M}^T = [M(q, t)|A^T(q, \dot{q}, t)]$ has full rank, $n$, at each instant of time,

4. and the dynamic nature of the non-ideal force of constraint is specified for the given mechanical system through a prescription of the $n$-vector $C(q, \dot{q}, t)$ as in Eq. (6),

the explicit acceleration of the constrained system $S$ is the same as the explicit acceleration of the constrained auxiliary mechanical system $S_A$ whose unconstrained equation of motion is described by the equation

$$M_A(q, t)\ddot{q} = Q(q, \dot{q}, t) + A^+(q, \dot{q}, t)z(q, \dot{q}, t) := Q_z(q, q, t), \quad (45)$$

which is subjected to the same constraints as system $S$ (specified in condition 2 above), and the same prescription of the $n$-vector $C$ as system $S$ (specified in condition 4 above). The mass matrix of the auxiliary system $S_A$ is given by

$$M_A = M(q, t) + \alpha^2 A^+(q, \dot{q}, t) A(q, \dot{q}, t) > 0. \quad (46)$$

The (sufficiently smooth, say, $C^2$) $m$-vector $z$ in Eq. (45) and the scalar $\alpha \neq 0$ in Eq. (46) are both arbitrary.

*Proof* Since the constrained acceleration of the constrained system is described by [2]

$$\ddot{q} = \begin{bmatrix} (I - A^+ A) M \\ A \end{bmatrix}^+ \begin{bmatrix} Q + C \\ b \end{bmatrix} := \bar{M}^+ \begin{bmatrix} Q + C \\ b \end{bmatrix}, \quad (47)$$

use of Lemmas 3 and 4 show that this constrained acceleration $\ddot{q}$ can also be expressed as

$$\ddot{q} = \begin{bmatrix} (I - A^+ A)M \\ A \end{bmatrix}^+ \begin{bmatrix} Q + C \\ b \end{bmatrix} = \begin{bmatrix} (I - A^+ A)M_A \\ A \end{bmatrix}^+ \begin{bmatrix} Q + C \\ b \end{bmatrix} \\ = \begin{bmatrix} (I - A^+ A)M_A \\ A \end{bmatrix}^+ \begin{bmatrix} (Q + A^+ z) + C \\ b \end{bmatrix}. \quad (48)$$

But, as seen in the previous section, this is exactly the constrained acceleration that would arise for the auxiliary mechanical system $S_A$ whose mass matrix is $M_A$ and that is subjected to the 'given' force $Q_z = Q + A^+z$ (instead of $Q$), and that has (1) the same constraints imposed on it as the original system $S$ (given in condition 2 above) and (2) the same specification regarding the nature of the non-ideal constraint force as for the system $S$ (given in condition 4 above). Last, because of item 3 above and Lemma 1, the mass matrix of the auxiliary system $M_A > 0$. $\square$

**Result 2:** The explicit acceleration of the constrained system $S$ described in Result 1 is explicitly given by

$$\ddot{q} = a_{A,z} + M_A^{-1/2} B_A^+ (b - A a_{A,z}) + M_A^{-1/2} (I - B_A^+ B_A) M_A^{-1/2} C, \quad (49)$$

where

$$M_A = M + \alpha^2 A^+ A > 0, \quad (50)$$

$$a_{A,z} = M_A^{-1}(Q + A^+z) = M_A^{-1}Q_z, \quad (51)$$

and

$$B_A = A M_A^{-1/2}. \quad (52)$$

The sufficiently smooth $m$-vector $z$ and the scalar $\alpha \neq 0$ are both arbitrary, as before.
Proof From Result 1, we know that the acceleration of the constrained system S is the same as that of the constrained system S<sub>A</sub>. But the system S<sub>A</sub> has a positive definite mass matrix, and so applying the fundamental equation to the unconstrained system described by Eq. (45) we get [3]

$$
M_A \ddot{q} = Q + A^+ z + M_A^{1/2} B_A^+ (b - A a_{A,z}) + M_A^{1/2} (I - B_A^+ B_A) M_A^{-1/2} C, \quad (53)
$$

from which Eq. (49) follows upon pre-multiplication on both sides by $M_A^{-1}$.
□

We have already noted that when ˆM has full rank, the acceleration of the constrained system is unique, and hence cannot depend on either the value chosen for α ≠ 0 and/or on the arbitrary vector z; yet the right hand side of Eq. (49) appears to have the arbitrary vector z in it. We shall now show that the right-hand side of Eq. (49) is indeed independent of this arbitrary vector z and also the parameter α ≠ 0.

**Corollary 1** The acceleration of the constrained system S does not depend on

(a) the choice of the sufficiently smooth vector z, and

(b) the parameter $\alpha \neq 0$.

It can also be expressed as

$$
\ddot{q} = a_A + M_A^{-1/2} B_A^+ (b - A a_A) + M_A^{-1/2} (I - B_A^+ B_A) M_A^{-1/2} C, \quad (54)
$$

where $M_A = M + \alpha^2 A^+ A$, $a_A = M_A^{-1} Q$, and $B_A = A M_A^{-1/2}$.

Proof (a) We only need to consider the first two members on the right-hand side of Eq. (49). Expanding them, we have

$$
\begin{align*}
& M_A^{-1}(Q + A^+z) + M_A^{-1/2} B_A^+[b - A M_A^{-1}(Q + A^+z)] \\
&= a_A + M_A^{-1/2} B_A^+(b - A a_A) + M_A^{-1} A^+ z - M_A^{-1/2} B_A^+ A M_A^{-1} A^+ z \\
&= a_A + M_A^{-1/2} B_A^+(b - A a_A) + M_A^{-1/2} [I - B_A^+ B_A] M_A^{-1/2} A^+ z \\
&= a_A + M_A^{-1/2} B_A^+(b - A a_A). \tag{55}
\end{align*}
$$

In the first equality above we have used the definition of the *n*-vector $a_A$ given by Corollary 1, and in going from the second to the third equality above we have used Lemma 2. Hence, the result. Since the right-hand side of Eq. (49) does not depend on the choice of the vector $z$, Eq. (54) is obtained from Eq. (49) by simply setting $z \equiv 0$.

(b) Noting that $B_A^+ = M_A^{-1/2} A^T (A M_A^{-1} A^T)^+ = M_A^{-1/2} A^T H^+$ and relation (51), Eq. (49) can be rewritten as

$$
\ddot{q} = M_A^{-1} (I - A^T H^+ A M_A^{-1}) (Q + A^+ z + C) + M_A^{-1} A^T H^+ b. \quad (56)
$$

We show that the right-hand side of Eq. (49) does not depend on the parameter $\alpha \neq 0$. To prove this we show that

$$
\frac{\partial[M_A^{-1}(I - A^T H^+ A M_A^{-1})]}{\partial \alpha} = 0 \quad \text{and} \quad \frac{\partial[M_A^{-1} A^T H^+]}{\partial \alpha} = 0. \qquad (57)
$$

In Lemma 5, part (h), we have shown that the second equality in Eq. (57) is true, since $M_A > 0$. To prove the first equality we note that the left-hand side of the equality can be written as

$$
\begin{align*}
& \frac{\partial M_A^{-1}}{\partial \alpha} - \frac{\partial [M_A^{-1} A^T H^+]}{\partial \alpha} A M_A^{-1} - M_A^{-1} A^T H^+ A \frac{\partial M_A^{-1}}{\partial \alpha} \\
&= \frac{\partial M_A^{-1}}{\partial \alpha} - M_A^{-1} A^T H^+ A \frac{\partial M_A^{-1}}{\partial \alpha} \\
&= -2\alpha M_A^{-1} A^+ A M_A^{-1} + 2\alpha M_A^{-1} A^T H^+ A M_A^{-1} A^+ A M_A^{-1} \\
&= -2\alpha M_A^{-1} [I - A^T H^+ A M_A^{-1}] A^+ A M_A^{-1} = 0. \tag{58}
\end{align*}
$$

In Eq. (58), the first equality follows again from Lemma 5, part (h), the second from Lemma 5, part (a), and the last equality follows from Lemma 5, part (b). Hence, the acceleration given in Eq. (49) does not depend on the value of the parameter $\alpha \neq 0$.
□
**Corollary 2** Results 1 and 2, as well as Corollary 1 are applicable when the mass matrix M of the system S is positive definite.

*Proof* When the mechanical system S has a mass matrix M that is positive definite, then condition 3 in Result 1 is automatically satisfied, and hence the results obtained so far are all applicable. □

Using Result 2 and Corollaries 1 and 2, we have the following useful result:

*Result 3:* When the mass matrix M of the unconstrained system S is

(i) positive definite, or

(ii) positive semi-definite and $\hat{M}$ has rank n,

one can obtain the acceleration of the constrained system S by simply replacing its mass matrix M (in the unconstrained equation of motion) by the positive definite mass matrix $M_A = M + A^+A$, and then using the fundamental equation [3]. The explicit expression for the acceleration of the constrained system S is then given by

$$ \ddot{q} = a_A + M_A^{-1/2} B_A^+(b - Aa_A) + M_A^{-1/2} (I - B_A^+ B_A) M_A^{-1/2} C, \quad (59) $$

where $M_A = M + A^+A$, $a_A = M_A^{-1}Q$, and $B_A = AM_A^{-1/2}$.

*Proof* Set $z = 0$ in relation (49), and $\alpha = 1$ in relation (50). Now the auxiliary system $S_A$ is identical to S except that its mass matrix is now $M_A = M + A^+A > 0$ instead of M. Equation (59) follows from Eq. (49) on making these substitutions. □

Since in much (nearly all) of analytical dynamics, D'Alembert's principle is assumed to hold at each instant of time, we provide, for easy reference, the following narrower result:

*Result 4:* For a mechanical system D:

1. that satisfies D'Alembert's principle at each instant of time, and

2. whose unconstrained equation of motion is described by $M(q, t)\ddot{q} = Q(q, \dot{q}, t)$, where the n by n mass matrix M is positive semi-definite,

3. and that is further subjected to the m consistent constraints described by Eqs. (2) and (3) [or alternatively, by Eq. (4)], which are satisfied by the initial conditions $q_0$ and $\dot{q}_0$,

4. such that the n by n+m matrix $\hat{M}^T = [M(q, t)|A^T(q, \dot{q}, t)]$ has full rank n at each instant of time,

the acceleration of the constrained system D is the same as the acceleration of the auxiliary constrained mechanical system $D_A$ whose

1. unconstrained equation of motion is described by the equation

$$ M_A(q, t)\ddot{q} = Q(q, \dot{q}, t) + A^+(q, \dot{q}, t)z(q, \dot{q}, t) := Q_z(q, q, t), \quad (60) $$

and

2. that is subjected to the same constraints (specified in 3 above) as system D.

The mass matrix of the auxiliary system $D_A$ is given by

$$ M_A = M(q, t) + \alpha^2 A^+(q, \dot{q}, t) A(q, \dot{q}, t) > 0. \quad (61) $$

The sufficiently smooth m-vector z in Eq. (60) and the scalar $\alpha \neq 0$ in Eq. (61) are both arbitrary. The explicit expression for the acceleration of the constrained system D is then given by

$$ \ddot{q} = a_{A,z} + M_A^{-1/2} B_A^+(b - Aa_{A,z}), \quad (62) $$

where, as before, $a_{A,z} = M_A^{-1}(Q + A^+z) = M_A^{-1}Q_z$, and $B_A = AM_A^{-1/2}$.

*Proof* When the system satisfies D'Alembert's principle at each instant of time the n-vector $C(q, \dot{q}, t) = 0$. The result follows by setting $C = 0$ into Results 1 and 2. □
**Corollary 3** For the mechanical system D described in Result 4, which satisfies D'Alembert's principle, as long as the positive semi-definite (or positive definite) matrix M is such that $\hat{M}$ has rank n at each instant of time, the acceleration of the constrained system can be obtained explicitly by imagining that the matrix M is replaced by the matrix $M_A = M + A^+A$, and then using the fundamental equation. Indeed, the acceleration of the constrained system D is given by

$$ \ddot{q} = a_A + M_A^{-1/2} B_A^+ (b - Aa_A), \qquad (63) $$

where $M_A = M + A^+A$, $a_A = M_A^{-1}Q$, and $B_A = AM_A^{-1/2}$.

*Proof* Set $z = 0$ and $\alpha = 1$ in Eqs. (60) and (61), and use Eq. (62) to obtain Eq. (63). $\square$

## 3.2 Explicit equation for constraint force $Q^c$

By using the fundamental equation on the auxiliary system $S_A$, we have so far obtained the acceleration of a constrained mechanical system $S$ whose mass matrix $M$ could be positive semi-definite, and for which the matrix $\hat{M}$ has full rank $n$ at each instant of time. We now ask if we can also determine the constraint force $Q^c$ acting on $S$ directly from the equation of motion of the constrained auxiliary system $S_A$.

Once the acceleration $\ddot{q}$ of the constrained system is known from Eq. (49), the constraint force acting on the system $S$ at each instant of time can indeed be determined by using Eq. (5) as

$$ Q^c = M\ddot{q} - Q. \qquad (64) $$

Recall that the acceleration of the constrained system $S$ described in Result 1 was obtained from Eq. (53). We created the auxiliary system $S_A$ with a mass matrix $M_A > 0$ (instead of $M \ge 0$ for the system $S$) and with the 'given' force $Q_z$ (instead of $Q$). Application of the fundamental equation to this auxiliary system $S_A$ (that is subjected to same constraints as $S$ and the same non-ideal constraint forces prescribed by specifying the vector $C(q, \dot{q}, t)$ as $S$) then yields the correct acceleration of the constrained system $S$.

We can now, however, think of the auxiliary system $S_A$ as a constrained system in its own right. Its unconstrained equation of motion is $M_A\ddot{q} = Q + A^+z$. The matrix $M_A$ is positive definite and the 'given' force acting on this auxiliary unconstrained system is $Q_z = Q + A^+z$, where $z$ is a sufficiently smooth arbitrary $m$-vector. This unconstrained system is subjected to the same constraint specifications as the constrained system $S$. Since $M_A > 0$, using the fundamental equation, the equation of motion for the acceleration of this constrained system $S_A$ is then obtained from the relation

$$ M_A \ddot{q} = Q_z + M_A^{1/2} B_A^+ (b - Aa_{A,z}) + M_A^{1/2} (I - B_A^+ B_A) M_A^{-1/2} C := Q_z + Q_{A,z}^c. \qquad (65) $$

The presence of the constraints (which are identical for both systems $S_A$ and $S$) generates a constraint force $Q_{A,z}^c$ on system $S_A$ as shown by the last equality in Eq. (65). But from Result 2, upon pre-multiplication on both sides of Eq. (65) by $M_A^{-1}$ we obtain the correct acceleration of the constrained system $S$, since the accelerations of both the constrained systems $S$ and $S_A$ are identical at each instant of time, as shown in Result 1. Thus, the constrained acceleration of both systems $S$ and $S_A$ can be written compactly as

$$ \ddot{q} = M_A^{-1} (Q_z + Q_{A,z}^c). \qquad (66) $$

We next investigate how the forces of constraint for the systems $S$ and $S_A$ are related. More specifically, we ask how one might directly obtain the force of constraint $Q^c$ acting on system $S$ as given in Eq. (64) from the force of constraint $Q_{A,z}^c$ acting on system $S_A$, which we have explicitly obtained in Eq. (65). We now show that this can indeed be done by appropriately choosing the vector $z$.

**Result 5:** For the unconstrained auxiliary system $S_A$ described in Result 1, setting $z = \alpha^2 b$ (see Eq. (4) for definition of the vector $b$), where $\alpha \neq 0$ is arbitrary, will cause the constraint force on the constrained auxiliary system $S_A$ to equal the constraint force on the constrained system $S$ at each instant of time, so that
$$Q_{A,\alpha^2 b}^c = Q^c, \qquad (67)$$

where $Q_{A,\alpha^2 b}^c$ is obtained from relation (65) with $M_A = M + \alpha^2 A^+ A$, $B_A = A M_A^{-1/2}$, and $z = \alpha^2 b$, where $\alpha \neq 0$ is any real number. Alternatively stated, with the choice $z = \alpha^2 b$, the equations of motion of the constrained systems $S$ and $S_A$ become

$$M \ddot{q} = Q + Q^c \quad \text{and} \quad M_A \ddot{q} = Q_z + Q_{A,\alpha^2 b}^c = Q_z + Q^c, \qquad (68)$$

respectively. Both equations in Eq. (68) are equivalent to one another.

*Proof* At each instant of time, the constraint force $Q^c$ on the mechanical system $S$ due to the constraints imposed on it are given by relation (64), in which $\ddot{q}$ is the acceleration of the constrained system $S$ (and is also the acceleration of the auxiliary constrained system $S_A$). Now the constraint (4) implies, upon pre-multiplication on both sides by $\alpha^2 A^+$, that

$$\alpha^2 A^+ A \ddot{q} - \alpha^2 A^+ b = 0, \qquad (69)$$

where $\ddot{q}$ is again the acceleration at each instant of time of the constrained system $S$. Adding the left-hand side of Eq. (69) to the right-hand side of Eq. (64) we get

$$Q^c = M \ddot{q} + \alpha^2 A^+ A \ddot{q} - (Q + \alpha^2 A^+ b) = M_A \ddot{q} - (Q + \alpha^2 A^+ b) = M_A \ddot{q} - Q_{\alpha^2 b}. \qquad (70)$$

But the acceleration $\ddot{q}$ of the constrained system $S$ at each instant of time is also given by Eq. (66)! Substituting $\ddot{q}$ from Eq. (66) into the last equality in Eq. (70) then yields

$$Q^c = Q_z + Q_{A,z}^c - Q_{\alpha^2 b}. \qquad (71)$$

When $z = \alpha^2 b$, the result then follows.

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

For simplicity, we particularize our previous result by choosing $\alpha = 1$.

*Result 6:* Consider a mechanical system $S$ whose unconstrained equation of motion is $M \ddot{q} = Q$, with $M \ge 0$, and that is further subjected to the constraints given by Eqs. (2) and (3) [or alternately, by Eq. 4]. Given the $n$-vector $C$ in Eq. (6), the equation of motion of the constrained system $S$, under the proviso that $\hat{M}$ has full rank at each instant of time (so that the system has a unique equation of motion), is given by the equation of motion of the constrained auxiliary system $S_A$,

$$M_A \ddot{q} = (Q + A^+ b) + M_A^{1/2} B_A^+ (b - A a_{A,b}) + M_A^{1/2} (I - B_A^+ B_A) M_A^{-1/2} C := Q_b + Q^c, \qquad (72)$$

where,

$$M_A = M + A^+ A, \quad Q_b = Q + A^+ b, \quad B_A = A M_A^{-1/2}. \qquad (73)$$

The acceleration of the constrained system $S$ is explicitly given by

$$\ddot{q} = a_{A,b} + M_A^{-1/2} B_A^+ (b - A a_{A,b}) + M_A^{-1/2} (I - B_A^+ B_A) M_A^{-1/2} C, \qquad (74)$$

where $a_{A,b} = M_A^{-1} (Q + A^+ b) := M_A^{-1} Q_b$. The constraint force on $S$ due to the presence of the constraints is explicitly given by

$$Q^c = M_A^{1/2} B_A^+ (b - A a_{A,b}) + M_A^{1/2} (I - B_A^+ B_A) M_A^{-1/2} C. \qquad (75)$$

*Proof* This is just a combined form of Results 2 and 5, specialized for $\alpha = 1$ and $z = b$.

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

We have thus proved the following somewhat remarkable result. When a mechanical system $S$ (1) has a positive semi-definite or a positive definite mass matrix $M$, (2) is subjected to the constraint $A \ddot{q} = b$, and the non-ideal constraint force is specified by the vector $C(q, \dot{q}, t)$, (3) has a matrix $\hat{M}$ being full rank, all that is needed to obtain the explicit acceleration of the constrained system and the explicit constraint force acting on it is to:
(i) imagine that the unconstrained system is replaced by an auxiliary system $S_A$ with a mass matrix $M_A := M + A^+A$,

(ii) add $A^+b$ to the 'given' force $Q$ applied to $S$ to obtain the 'given' force $Q + A^+b$ acting on the auxiliary unconstrained system $S_A$, and

(iii) apply the fundamental equation to this auxiliary system subjected to the given constraint conditions to which the system $S$ is subjected. This will give Eq. (72).

The acceleration of the constrained auxiliary system $S_A$ is the same, at each instant of time, as that of the constrained system $S$; and, the constraint force $Q^c$ at each instant of time acting on $S$ is simply the constraint force calculated on the auxiliary system to which the 'given' force $Q + A^+b$ is applied. In short, the motion of the auxiliary system $S_A$ will exactly mimic the actual system $S$. Indeed, the Lagrange equation describing the constrained system $S$ is exactly the equation for the constrained auxiliary system $S_A$.

We note that for systems for which D'Alembert's principle is assumed to be true, we simply set the vector $C = 0$ in Eqs. (72), (74) and (75). When the matrix $M$ is positive definite, Results 5 and 6 are still valid since then $\dot{M}$ has full rank.

## 4 Application to rotational dynamics

To place the analytical results obtained hereto in perspective and to show their practical applicability, we consider the rotational dynamics of a rigid body. In this section, we obtain in a simple and direct manner a new form for the equation for the rotational motion of a rigid body in terms of quaternions. We shall assume that the system satisfies D'Alembert's principle. As we shall see, the rigid body system $D$ has a singular mass matrix, and we shall use the auxiliary system $D_A$ to directly obtain the requisite equations of motion of $D$, and also the required constraint torque $Q^c$.

Consider a rigid body that has an absolute angular velocity, $\omega \in \mathbb{R}^3$, with respect to an inertial coordinate frame. The components of this angular velocity with respect to its body-fixed coordinate frame whose origin is located at the body's center of mass are denoted $\omega_1, \omega_2$, and $\omega_3$. Let us assume, without loss of generality, that the body-fixed coordinate axes attached to the rigid body are aligned along its principal axes of inertia, where the principal moments of inertia are given by $J_1, J_2$, and $J_3$. The rotational kinetic energy of the rigid body is then simply

$$T = \frac{1}{2}\omega^T J \omega, \qquad (76)$$

where $\omega = [\omega_1, \omega_2, \omega_3]^T$ and $J = \text{diag}(J_1, J_2, J_3)$.

We shall describe the orientation of the rigid body through the use of quaternions, also known as the Euler parameters. They permit the description of large angle rotations without any geometric singularities. The unit quaternion $u$ written as a $4 \times 1$ column vector by

$$u = [u_0, \mathbf{u}^T] = \left[ \cos \frac{\theta}{2}, \mathbf{e}^T \sin \frac{\theta}{2} \right]^T \qquad (77)$$

describes the rotation of the body. By a unit quaternion, we mean that the square of its norm is unity so that

$$\mathcal{N}(u) := u^T u = 1. \qquad (78)$$

The 3-vector $\mathbf{e}$ in Eq. (77), which is described by any unit vector relative to an inertial coordinate frame of reference, is the rotation axis. The rotation of the body about this rotation axis is defined by the rotation angle $\theta$.

The body-fixed angular velocity $\omega$ and the unit quaternion $u$ of the rigid body are related by

$$\omega = 2E\dot{u} = -2\dot{E}u \quad \text{and} \quad \dot{u} = \frac{1}{2}E^T\omega, \qquad (79)$$

where the $3 \times 4$ matrix $E$ is

$$E = \begin{bmatrix} -u_1 & u_0 & u_3 & -u_2 \\ -u_2 & -u_3 & u_0 & u_1 \\ -u_3 & u_2 & -u_1 & u_0 \end{bmatrix}. \qquad (80)$$
The rotational kinetic energy of the rigid body given by Eq. (76) then becomes

$$T = 2\dot{u}^T E^T J E \dot{u} = 2u^T \dot{E}^T J \dot{E} u. \qquad (81)$$

In Eq. (81), we note that the symmetric matrices $E^T J E$ and $\dot{E}^T J \dot{E}$ are not positive definite since we can have

$$u^T E^T J E u = \dot{u}^T \dot{E}^T J \dot{E} \dot{u} = 0, \quad u, \dot{u} \neq 0. \qquad (82)$$

Using Eq. (78) and the definition of the matrix $E$ given in Eq. (80), it is easy to show that the identities

$$E u = 0, \qquad (83)$$

$$\dot{E} \dot{u} = 0, \qquad (84)$$

$$E^T E + u u^T = I_4, \qquad (85)$$

$$E E^T = I_3, \qquad (86)$$

are true, where $I_n$ denotes the $n \times n$ identity matrix.

To obtain Lagrange's rotational equations of motion, we begin by assuming that the components of the generalized coordinate vector $u$ are independent of one another. Lagrange's equations for this unconstrained system $D$—the system is unconstrained since the components of $u$ are all assumed to be independent of one another—then becomes

$$\frac{d}{dt} \left( \frac{\partial T}{\partial \dot{u}} \right) - \frac{\partial T}{\partial u} = \Gamma_u, \qquad (87)$$

which yields four second-order nonlinear differential equations. The 4-vector $\Gamma_u$ in Eq. (87) represents the generalized impressed quaternion torque that would exist if all the components of $u$ were actually independent. Its relation to a physically applied body torque 3-vector $\Gamma_B$ will be presented later on in this section. By appropriately carrying out the differentiations in Eq. (87), we obtain

$$\frac{d}{dt} \left( \frac{\partial T}{\partial \dot{u}} \right) = 4E^T J E \ddot{u} + 4 \dot{E}^T J E \dot{u} + 4E^T J \dot{E} \dot{u} \qquad (88)$$

and

$$\frac{\partial T}{\partial u} = 4 \dot{E}^T J \dot{E} u = -4 \dot{E}^T J E \dot{u}. \qquad (89)$$

Lagrange's equation, assuming the components of the unit quaternion $u$ are all independent, is then

$$4E^T J E \ddot{u} + 8 \dot{E}^T J E \dot{u} + 4E^T J \dot{E} \dot{u} = \Gamma_u. \qquad (90)$$

Using Eq. (84), Eq. (90) reduces to

$$M \ddot{u} := 4 E^T J E \ddot{u} = -8 \dot{E}^T J E \dot{u} + \Gamma_u := Q, \qquad (91)$$

where the 4 × 4 matrix $M = 4E^T J E$ multiplying the 4-vector $\ddot{u}$ is symmetric and positive semi-definite (singular) since its rank is three.

Recall that the equations of motion in Eq. (91), which are obtained directly from Lagrange's equation, presume that the four components of $u$ are independent of one another. Since $u$ represents a physical rotation, it must in actuality satisfy Eq. (78) because it must be a unit quaternion. Thus, the unconstrained system $D$ whose motion we have described by Eq. (91) must be further subjected to the constraint

$$\varphi := \mathcal{N}(u) - 1 = u^T u - 1 = 0. \qquad (92)$$

Differentiating Eq. (92) twice with respect to time, we obtain the constraint matrix equation

$$A \ddot{u} := u^T \ddot{u} = -\mathcal{N}(\dot{u}) := b. \qquad (93)$$

Using Result 4 from the previous section and noting that $A^+ = u$, we now obtain the auxiliary system $D_A$ whose unconstrained equation of motion is given by

$$M_A \ddot{u} := (4E^T J E + \alpha^2 u u^T) \ddot{u} = -8 \dot{E}^T J E \dot{u} + u z + \Gamma_u := Q_z, \qquad (94)$$
where $\alpha \neq 0$. The acceleration of this unconstrained system is given by $a_{A,z} = M_A^{-1} Q_z$. We shall, in what follows, subject this system $D_A$, whose unconstrained equation of motion is given by Eq.(94), to the same constraints that our actual system $D$ is subjected to, namely, to the constraint specified in Eq. (92) [or alternatively in Eq. (93)].

We first show that the 4 × 4 matrix $M_A = 4E^T J E + \alpha^2 uu^T$ is positive definite. To show this, we use Lemma 1 and prove that the rank of the matrix

$$ \hat{M} = \begin{bmatrix} M \\ A \end{bmatrix} = \begin{bmatrix} 4E^T J E \\ u^T \end{bmatrix}, \qquad (95) $$

is four. Consider the matrix product

$$ V \hat{M} W := \begin{bmatrix} E^T E & u \\ u^T & 0 \end{bmatrix} \begin{bmatrix} 4E^T J E \\ u^T \end{bmatrix} \begin{bmatrix} \frac{1}{4} E^T J^{-1} E + uu^T \end{bmatrix}. \qquad (96) $$

Using Eqs. (85), (86) and (92), we find that the matrix $V$ is nonsingular since

$$ V V^T = \begin{bmatrix} E^T E + uu^T & 0 \\ 0 & u^T u \end{bmatrix} = I_5, \qquad (97) $$

and so also the matrix $W$ since by Eqs. (83) and (85) we have

$$ WW^{-1} := \left( \frac{1}{4} E^T J^{-1} E + uu^T \right) (4E^T J E + uu^T) = E^T E + uu^T = I_4. \qquad (98) $$

Thus, the matrix $\hat{M}$ has the same rank as the matrix $V\hat{M}W$. But the right-hand side of Eq. (96) yields

$$ V \hat{M} W = \begin{bmatrix} I_4 \\ 0 \end{bmatrix}, \qquad (99) $$

whose rank is four. Hence, the augmented mass matrix $M_A$ is positive definite.

The rotational equations of motion in terms of quaternions are now explicitly found by simply using the fundamental equation given in Result 4 as

$$ (4E^T J E + \alpha^2 uu^T) \ddot{u} = -8 \dot{E}^T J E \dot{u} + uz + \Gamma_u + M_A^{1/2} B_A^+ (b - Aa_{A,z}). \qquad (100) $$

Since

$$ M_A^k = (4E^T J E + \alpha^2 uu^T)^k = 4^k E^T J^k E + \alpha^{2k} uu^T, \quad k = \pm 1, \pm 1/2, \qquad (101) $$

the acceleration $a_{A,z}$ in Eq. (100) is obtained from Eq. (94) as

$$ \begin{aligned} a_{A,z} &= M_A^{-1} Q_z \\ &= \left( \frac{1}{4} E^T J^{-1} E + \frac{1}{\alpha^2} uu^T \right) (-8 \dot{E}^T J E \dot{u} + uz + \Gamma_u). \end{aligned} \qquad (102) $$

We now proceed by assembling the various quantities required to determine the last member on the right-hand side of Eq. (100) by computing $Aa_{A,z}$, $b - Aa_{A,z}$, and finally $M_A^{1/2} B_A^+$. We find that

$$ \begin{aligned} Aa_{A,z} &= \left( \frac{1}{4} u^T E^T J^{-1} E + \frac{1}{\alpha^2} u^T u u^T \right) (-8 \dot{E}^T J E \dot{u} + uz + \Gamma_u) \\ &= \frac{1}{\alpha^2} u^T (-8 \dot{E}^T J E \dot{u} + uz + \Gamma_u) \\ &= \frac{1}{\alpha^2} (2 \omega^T J \omega + z + u^T \Gamma_u), \end{aligned} \qquad (103) $$

where to get the second equality above we have used Eqs. (78) and (83), and in the last equality we have used Eqs. (78) and (79). This yields

$$ b - Aa_{A,z} = -\mathcal{N}(\dot{u}) - \frac{1}{\alpha^2} (2\omega^T J \omega + z + u^T \Gamma_u). \qquad (104) $$
Next, using Eqs. (101) and (83) the matrix $B_A^+ := (AM_A^{-1/2})^+$ is found as

$$ 
\begin{aligned}
B_A^+ &= \left[ u^T \left( \frac{1}{2} E^T J^{-1/2} E + \frac{1}{\alpha} u u^T \right) \right]^+ \\
&= \left( \frac{1}{\alpha} u^T \right)^+ = \alpha u.
\end{aligned} 
\quad (105) $$

Finally, we obtain

$$ M_A^{1/2} B_A^+ = (2E^T J^{1/2} E + \alpha u u^T) \alpha u = \alpha^2 u, \quad (106) $$

by once again using Eqs. (78) and (83). Using relations (104) and (106), Eq. (100) becomes

$$ (4E^T J E + \alpha^2 u u^T) \ddot{u} = -8 \dot{E}^T J E \dot{u} - \alpha^2 \mathcal{N}(\dot{u}) u - 2(\omega^T J \omega) u + (I_4 - u u^T) \Gamma_u, \quad (107) $$

where as usual the real number $\alpha \neq 0$ is arbitrary.

Equation (107) constitutes a new set of Lagrange's equations for rotational motion in terms of quaternions, which are different in form from those found in Refs. [8–10]. In contrast to the formulations provided in these references, we have arrived at Eq. (107) by working directly with the kinetic energy that is not positive definite [Eq. (81)], and without the use of any Lagrange multipliers. Attempts in the past at using a non-positive definite kinetic energy have apparently resulted in the development of incorrect equations of motion [11]. Thus, use of the general equation of motion obtained in Sections 3 and 4 allows us to obtain in a simple, straightforward way the equation for rotational dynamics without the need to artificially augment the mass matrix to make it nonsingular as done in Refs. [8] and [9], or the need to determine the Lagrange multipliers as done in Refs. [9] and [10].

The generalized rotational acceleration $\ddot{u}$ is now explicitly obtained as

$$ \ddot{u} = \left( \frac{1}{4} E^T J^{-1} E + \frac{1}{\alpha^2} u u^T \right) \left[ -8 \dot{E}^T J E \dot{u} - \alpha^2 \mathcal{N}(\dot{u}) u - 2(\omega^T J \omega) u + (I_4 - u u^T) \Gamma_u \right]. \quad (108) $$

Expanding Eq. (108), and using Eqs. (78) and (83), we obtain

$$ \ddot{u} = -2E^T J^{-1} E \dot{E}^T J E \dot{u} - \mathcal{N}(\dot{u}) u + \frac{1}{4} E^T J^{-1} E \Gamma_u - \frac{8uu^T \dot{E}^T J E \dot{u} + 2u(\omega^T J \omega)}{\alpha^2}. \quad (109) $$

Also, we find that

$$ E \dot{E}^T = \frac{1}{2} \tilde{\omega}, \quad (110) $$

where

$$ \tilde{\omega} = \begin{bmatrix} 0 & -\omega_3 & \omega_2 \\ \omega_3 & 0 & -\omega_1 \\ -\omega_2 & \omega_1 & 0 \end{bmatrix}. \quad (111) $$

The first member on the right-hand side of Eq. (109) can therefore be expressed as $-\frac{1}{2}E^T J^{-1} \tilde{\omega} J \omega$. Also, using Eq. (81), $8uu^T \dot{E}^T J E \dot{u} = -2u(\omega^T J \omega)$, so that the explicit equation for the generalized acceleration of a rotating rigid body in terms of quaternions is

$$ \ddot{u} = -\frac{1}{2} E^T J^{-1} \tilde{\omega} J \omega - \mathcal{N}(\dot{u}) u + \frac{1}{4} E^T J^{-1} E \Gamma_u. \quad (112) $$

The connection between the generalized torque 4-vector $\Gamma_u$ and the physically applied torque 3-vector $\Gamma_B = [\Gamma_1, \Gamma_2, \Gamma_3]^T$, whose components $\Gamma_i$, $i = 1, 2, 3$ are about the body-fixed axes of the rotating body, is determined by the relations [8]

$$ \Gamma_u = 2E^T \Gamma_B, \quad (113) $$

$$ \Gamma_B = \frac{1}{2} E \Gamma_u. \quad (114) $$
In addition, by Eqs. (79) and (86) we have

$$ \mathcal{N}(\dot{u}) = \dot{u}^{\mathrm{T}}\dot{u} = \frac{1}{4}\omega^{\mathrm{T}}E E^{\mathrm{T}}\omega = \frac{1}{4}\omega^{\mathrm{T}}\omega = \frac{1}{4}\mathcal{N}(\omega). \qquad (115) $$

Hence, Eq. (112) becomes

$$ \ddot{u} = -\frac{1}{2}E^{\mathrm{T}}J^{-1}\tilde{\omega}J\omega - \frac{1}{4}\mathcal{N}(\omega)u + \frac{1}{2}E^{\mathrm{T}}J^{-1}\Gamma_B. \qquad (116) $$

Using Eqs. (79) and (85), we further note that upon pre-multiplying Eq. (110) by $E^T$ we get

$$ (E^{\mathrm{T}}E)\dot{E}^{\mathrm{T}} = (I_4 - uu^{\mathrm{T}})\dot{E}^{\mathrm{T}} = \dot{E}^{\mathrm{T}} + \frac{1}{2}u\omega^{\mathrm{T}} = \frac{1}{2}E^{\mathrm{T}}\tilde{\omega}, \qquad (117) $$

so that the matrix $\dot{E}^T$ becomes

$$ \dot{E}^{\mathrm{T}} = \frac{1}{2}E^{\mathrm{T}}\tilde{\omega} - \frac{1}{2}u\omega^{\mathrm{T}}. \qquad (118) $$

The equation of motion given by Eq. (107) then becomes

$$ (4E^{\mathrm{T}}JE + \alpha^2uu^{\mathrm{T}})\ddot{u} = -2E^{\mathrm{T}}\tilde{\omega}J\omega - \frac{\alpha^2}{4}\mathcal{N}(\omega)u + 2E^{\mathrm{T}}\Gamma_B, \qquad (119) $$

in terms of the body-fixed angular velocity 3-vector $\omega$, the unit quaternion 4-vector $u$, and the physically applied torque 3-vector $\Gamma_B$.

The requisite generalized torque of constraint needed to arrive at the correct equations of motion is found by utilizing Eqs. (104) and (106) so that

$$ \begin{aligned} Q_{A,z}^{c} &= M_{A}^{1/2} B_{A}^{+} (b - A a_{A,z}) \\ &= \alpha^2 u \left[ -\mathcal{N}(\dot{u}) - \frac{1}{\alpha^2} (2 \omega^{\mathrm{T}} J \omega + z + u^{\mathrm{T}} \Gamma_u) \right]. \end{aligned} \qquad (120) $$

Utilizing Result 5 and Eq. (113), the constraint torque is explicitly obtained by setting $z = \alpha^2 b = -\alpha^2 \mathcal{N}(\dot{u})$ so that

$$ Q^c = Q_{A,\alpha^2 b}^c = -2(\omega^T J \omega)u - uu^T \Gamma_u. \qquad (121) $$

Finally, we state that the equation of motion, which has (1) a positive definite mass matrix and (2) exactly the same acceleration—and also the same constraint torque acting on it—as the constrained system $D$ with $M \ge 0$, is explicitly found by considering the unconstrained system

$$ (4E^{\mathrm{T}}JE + \alpha^2uu^{\mathrm{T}})\ddot{u} = -8\dot{E}^{\mathrm{T}}JE\dot{u} - \alpha^2\mathcal{N}(\dot{u})u + \Gamma_u := Q_{\alpha^2 b} \qquad (122) $$

and further subjecting it to the constraint in Eq. (78) to obtain the explicit equation of motion of the constrained system $D_A$ given by

$$ \begin{aligned} (4E^{\mathrm{T}}JE + \alpha^2uu^{\mathrm{T}})\ddot{u} &= Q_{\alpha^2 b} + Q_{A,\alpha^2 b}^c, \\ &= \{-8\dot{E}^{\mathrm{T}}JE\dot{u} - \alpha^2\mathcal{N}(\dot{u})u + \Gamma_u\} + \{-2(\omega^{\mathrm{T}}J\omega)u - uu^{\mathrm{T}}\Gamma_u\}, \end{aligned} \qquad (123) $$

where $Q_{A,\alpha^2 b}^c$ and $Q_{\alpha^2 b}$ are obtained from Eqs. (121) and (122), respectively, and are shown by the two bracketed quantities.

It should be noted that the constrained equation of motion $D_A$ given by Eq. (123) explicitly provides the requisite equation of motion of the constrained system $D$—a system that has a singular mass matrix.

It also provides the requisite constraint torque the unconstrained system $D$ [Eq. (91)] is subjected to due to the enforcement of the unit quaternion constraint [Eq. (78)]. The closed form expression of this constraint torque is given by the second bracketed quantity on the right-hand side of Eq. (123).

In short, the equation of motion that has a positive definite mass matrix [Eq. (123)] is the equation of motion for the constrained system $D$ that has a singular, positive semi-definite mass matrix. Equation (119) [or (123)] is a new form for the equation of motion for rotational dynamics of a rigid body that has a positive definite mass matrix, and that is obtained directly through the use of a semi-definite kinetic energy quadratic form.
# 5 Conclusions and remarks

This paper develops a new and simple description for the constrained motion of general mechanical systems that may be described using mass matrices that may or may not be singular. In Lagrangian mechanics, these systems are often seen to arise when the number of coordinates chosen to describe the configuration of a system exceeds the minimum number required to describe it. The use of such additional coordinates is very useful when modeling complex mechanical systems since they provide considerable convenience and ease in the modeling process.

To date, the use of a general explicit equation for constrained motion that is broad enough to encompass systems in which $M \ge 0$ has been unavailable. The results obtained to date (see, for example, Refs. [1] and [2]) provide very different explicit equations for describing the constrained motion of general mechanical systems; these equations depend critically on whether the mass matrix $M$ is positive definite, or semi-definite (and singular).

The general equation of motion for constrained systems obtained herein yields both the unique generalized acceleration and the unique force of constraint when $M \ge 0$ and when the matrix $\hat{M}^T = [M|A^T]$ has full rank. It includes systems that may or may not satisfy D'Alembert's Principle at each instant of time, subjected to general holonomic and/or nonholonomic constraints. The nonholonomic constraints can be nonlinear in the generalized velocities, and the constraints need not be functionally independent of each other. The absence of the notion of Lagrange multipliers throughout this paper is noteworthy.

The main contributions of this paper are as follows:

1. The development of a unified equation of motion for general constrained mechanical systems that (1) are described with mass matrices that can be either positive semi-definite ($M \ge 0$) or positive definite ($M > 0$), (2) are constrained by the general set of $m$ constraint equations given in the form $A(q, \dot{q}, t)\ddot{q} = b(q, \dot{q}, t)$, (3) may or may not satisfy D'Alembert's Principle at each instant of time, and (4) yield a full rank matrix $\hat{M}^T = [M|A^T]$. This rank condition is a necessary and sufficient condition for the acceleration of the constrained system to be uniquely defined at each instant of time [2].

2. For a system $S$ that is described with an unconstrained equation of motion $M\ddot{q} = Q$, which may have a singular mass matrix $M$, its corresponding constrained generalized acceleration is found by the (virtual) creation of an auxiliary unconstrained system $S_A$ that always has a positive definite mass matrix and that is subjected to the same constraints as system $S$. The auxiliary unconstrained system $S_A$ is given by the equation $(M + \alpha^2 A^+ A)\ddot{q} = Q + A^+ z$, where $z$ is an arbitrary $C^2$ vector, and $\alpha \ne 0$. By subjecting this unconstrained system to the same constraints as those imposed on the system $S$, use of the so-called fundamental equation of constrained motion (to this auxiliary system) explicitly yields, in closed form, the constrained generalized acceleration of the system $S$.

3. It is somewhat surprising that the simplest way to get the correct acceleration of a constrained mechanical system no matter whether its mass matrix is positive definite or positive semi-definite is to simply add the matrix $A^+A$ to the mass matrix $M$ of the unconstrained system.

4. Under any set of general holonomic and/or nonholonomic constraints, the force of constraint acting on the system $S$—a system that may be described with a singular mass matrix—is the same force of constraint acting on the auxiliary system $S_A$ when the vector $z$ in item 2 is chosen so that $z = \alpha^2 b$. Hence, the auxiliary system $S_A$ when subjected to the same constraints as system $S$ directly furnishes the force of constraint acting on the system $S$.

5. The use of the general equation of motion developed in this paper for deriving the Lagrangian equations of motion describing the rotational dynamics of a rigid body in terms of quaternions is shown to have the advantages of directness and simplicity when compared to other derivations obtained to date [8–11]. The form of the equation of motion so obtained [Eq. (119)] appears to have been unavailable to date.

# References

1. Udwadia, F.E., Kalaba, R.E.: A new perspective on constrained motion. Proc. R. Soc. Lond. Ser. A **439**, 407–410 (1992)

2. Udwadia, F.E., Phohomsiri, P.: Explicit equations of motion for constrained mechanical systems with singular mass matrices and applications to multi-body dynamics. Proc. R. Soc. Lond. Ser. A **462**, 2097–2117 (2006)

3. Udwadia, F.E., Kalaba, R.E.: On the foundations of analytical dynamics. Int. J. Non-Linear Mech. **37**, 1079–1090 (2002)

4. Udwadia, F.E., Kalaba, R.E.: What is the general form of the explicit equations of motion for constrained mechanical systems? J. Appl. Mech. **69**, 335–339 (2002)

5. Graybill, F.: Matrices with Applications in Statistics. Wadsworth, USA (1983)
6. Wedin, P.A.: Perturbation theory for pseudo-inverses. BIT **13**, 217–232 (1973)

7. Golub, G.H., Pareyra, V.: The differentiation of pseudoinverses and nonlinear least squares problems whose variables separate. SIAM J. Numer. Anal. **10**, 413–432 (1973)

8. Udwadia F.E., Schutte A.D.: An alternative derivation of the quaternion equations of motion for rigid-body rotational dynamics. J. Appl. Mech. **77**(3), May (2010)

9. Morton, H.: Hamiltonian and Lagrangian formulations of rigid-body rotational dynamics based on the Euler parameters. J. Astronaut. Sci. **59**, 569–591 (1993)

10. Nikravesh, P.E., Wehage, R.A., Kwon, O.K.: Euler parameters in computational kinematics and dynamics. ASME J. Mech. Transm. Autom. Des. **107**, 358–369 (1985)

11. Vadali, S.: On the Euler parameter constraint. J. Astronaut. Sci. **36**, 259–265 (1988)