# Equations of Motion for Nonholonomic, Constrained Dynamical Systems via Gauss's Principle

Robert Kalaba, Firdaus Udwadia

► To cite this version:

Robert Kalaba, Firdaus Udwadia. Equations of Motion for Nonholonomic, Constrained Dynamical Systems via Gauss's Principle. Journal of Applied Mechanics, 1993, 60, pp.662 - 662. 10.1115/1.2900855. hal-01406292

HAL Id: hal-01406292

https://hal.science/hal-01406292

Submitted on 1 Dec 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 Nonholonomic, Constrained Dynamical Systems via Gauss's Principle

In this paper we develop an analytical set of equations to describe the motion of discrete dynamical systems subjected to holonomic and/or nonholonomic Pfaffian equality constraints. These equations are obtained by using Gauss's Principle to recast the problem of the constrained motion of dynamical systems in the form of a quadratic programming problem. The closed-form solution to this programming problem then explicitly yields the equations that describe the time evolution of constrained linear and nonlinear mechanical systems. The direct approach used here does not require the use of any Lagrange multipliers, and the resulting equations are expressed in terms of two different classes of generalized inverses—the first class pertinent to the constraints, the second to the dynamics of the motion. These equations can be numerically solved using any of the standard numerical techniques for solving differential equations. A closed-form analytical expression for the constraint forces required for a given mechanical system to satisfy a specific set of nonholonomic constraints is also provided. An example dealing with the position tracking control of a nonlinear system shows the power of the analytical results and provides new insights into application areas such as robotics, and the control of structural and mechanical systems.

R. E. Kalaba¹

Professor of Economics,
Department of Eletrical Engineering
and Biomedical Engineering.

F. E. Udwadia¹

Professor of Mechanical Engineering,
Civil Engineering and Decision Systems,
Mem. ASME.

University of Southern California,
Los Angeles, CA 90089-1453

## 1 Introduction

D'Alembert's principle, which gives a complete conceptual solution to problems of classical mechanics, hinges upon the first-order virtual work done by the impressed (given) forces and that done by the forces of inertia (Lanczos, 1970). The former can often be expressed in terms of the variation of a potential energy function (Lanczos, 1970). By integrating with respect to time, the virtual work done by the forces of inertia can be transformed into a true variation (Rosenberg, 1972). Thus for holonomic systems, D'Alembert's principle can be reformulated as Hamilton's variational principle, which requires that a definite integral be stationary (Lanczos, 1970). The set of Lagrangian equations of motion that follow remain invariant under arbitrary, one-to-one point transformations.

It was in 1829 that Gauss (1829) gave an aesthetic and ingenious reinterpretation of D'Alembert's principle, changing it into a true minimum principle. This principle is applicable

to systems with general constraints, including configuration constraints (Rosenberg, 1972). Gauss argued that the determination of the motion of an *n*-degree-of-freedom system in which positions and velocities were known, hinged on our ability to determine the accelerations under the given applied forces. He formulated the principle of "least constraint" for describing the motion of mechanical systems. This principle is closely analogous to his celebrated "method of least squares," a method he developed and applied to the adjustment of errors in measurements. Unlike Hamilton's principle, the principle of least constraint has the additional advantage of not requiring any integration in time. Hertz gave a geometrical interpretation of Gauss's principle for the special case when the impressed forces vanish (Hertz, 1917). He showed that in this case Gauss's "constraint" can be interpreted as the geodesic curvature of the configuration point in 3*n*-dimensional space. Appell and Gibbs (see Pars, 1979) further extended the principle to apply to nonholonomic conditions and in cases where it may be advantageous to use kinematical variables (Lanczos, 1970). They used the idea of pseudo-coordinates (see, Pars 1979) which has, more recently, been again explored by Shan (1975)<fnref n="2" />. Synge (1926) has also provided an alternative set of equations of motion of nonholonomic systems in terms of the geometry

<footnote n="2">²The authors are thankful to an anonymous reviewer for pointing out this reference to them.</footnote>
of the resultant trajectories. As such, his formulation is difficult to directly apply to engineering problems.

From a practical standpoint, however, the computational difficulties of directly solving a minimization problems at each instant of time to describe the motion of a mechanical system made Gauss's principle unattractive at the time. This caused mechanicians of the late 18th, 19th, and 20th centuries to expound on, and mainly utilize the methods of Jacobi and Hamilton in the solution of problems in mechanics. Modern day texts in classical mechanics usually concentrate on these two latter approaches (e.g., Arnold, 1980), often relegating Gauss's principle to the position of a theoretically insightful approach, yet practically speaking, an unusable novelty.

In this paper we show that with our improved understanding of generalized inverses of rank-deficient matrices, Gauss's principle may offer a new, direct and oftentimes simpler approach to handling complex problems in mechanical systems. This is true in particular where nonholonomic and rheonomic constraints may be present. The key idea is that Gauss's Principle allows us to reformulate the equations of motion of constrained mechanical systems as a quadratic programming problem. In this paper we solve this quadratic programming problem, and thereby obtain a new set of explicit equations governing the motion of constrained, discrete dynamical systems. In contrast with the hereto used standard approach, which requires the use of Lagrange multipliers (e.g., see Rosenberg, 1972) or an expanded set of coordinates (Appell, 1925), the new approach developed here does away with the need for Lagrange multipliers. Furthermore, these equations are valid for both holonomic and nonholonomic constraints thereby treating both these types of constraints with equal consideration, and ease. The paper thus presents a unified approach to the handling of equality constraints in the analytical mechanics of discrete systems. In addition, an explicit expression is provided for the determination of the forces-of-constraint required so that a discrete mechanical system satisfies a given set of nonholonomic constraints.

Wang and Huston (1989) have looked at the representation of the equations of motion for nonholonomic systems, more from a matrix algebra standpoint. They also obtain equations of motion which do not involve any Lagrange multipliers. The equations obtained in this paper are, in a sense, generalizations of their results because we present the results in terms of non-specific generalized inverses which belong to certain classes. With the flexibility of choosing any generalized inverse from a given class of inverses, specific generalized inverses suitable for specific problem situations can often be found quickly and efficiently.

In Section 2 we present a simple, short derivation of Gauss's principle for nonholonomic systems. The constraints are taken to be in Pfaffian form. The exposition in this section, we believe, is not available in the current literature (e.g., in Whittaker (1917), Syngve (1926) and Pars (1979)), and provides some new insights. In Section 3 we use the results obtained in Section 2 to provide an exact solution to the constrained quadratic minimization problem governing the motion of constrained, discrete mechanical systems. In Section 4 we obtain explicit expressions for the constraint forces needed to satisfy the imposed constraints. Explicit equations for systems subjected to nonholonomic constraints are also provided. Section 5 illustrates our results using three numerical examples. The first deals with nonholonomic constraints, the second with the nonlinear oscillations of a pendulum subjected to nonlinear constraints. The third deals with the determination of the forces of constraint that need to be imposed on an oscillatory system described by a coupled Duffing's oscillator so that a specified time-dependent trajectory (constraint) is followed in configuration space. This latter example shows the power of our new formulation to possible applications in the field of robotics and position-tracking control of mechanical systems.

## 2 Gauss's Principle

Consider a holonomic mechanical system with n-degrees-of-freedom whose generalized coordinates are $q_1, q_2, q_3, \dots, q_n$. The Lagrange equations describing the motion of the system may be written as

$$ \frac{d}{dt} \left( \frac{\partial T}{\partial \dot{q}_r} \right) - \frac{\partial T}{\partial q_r} = Q_r, \quad r=1, 2, \dots, n, \qquad (1) $$

where T denotes the kinetic energy and $Q_r$ is the generalized impressed force. The kinetic energy can be expressed as

$$ T = \frac{1}{2} \sum_{i,j=1}^{n} a_{ij} \dot{q}_i \dot{q}_j + \sum_{i=1}^{n} b_i \dot{q}_i + c, \qquad (2) $$

where, in general, the $a_{ij}$ and $b_i$ and $c$ are functions of the generalized coordinates and time.

Assume now that the system is subjected to an additional $p$ ($p < n$) independent nonholonomic, Pfaffian constraints of the form

$$ \sum_{r=1}^{n} \alpha_{kr} dq_r + \beta_{kt} dt = 0, \quad k=1, 2, \dots, p, \qquad (3) $$

where $\alpha_{kr}$ and $\beta_{kt}$ are functions of the generalized coordinates and time. We note that the constraints may be scleronomic or rheonomic, catastrophic or a catastrophic (Rosenberg, 1972). These $p$ constraints may be thought of as imposing additional constraint forces, $Q_i$, on our system, thereby altering the set of Eqs. (1) to

$$ \frac{d}{dt} \left( \frac{\partial T}{\partial \dot{q}_r} \right) - \frac{\partial T}{\partial q_r} = Q_r + Q'_r, \quad r=1, 2, \dots, n. \qquad (4) $$

Expanding the first term in Eq. (4) we get

$$ \begin{aligned} \frac{d}{dt} \left( \frac{\partial T}{\partial \dot{q}_r} \right) &= \frac{d}{dt} \left\{ \sum_{s=1}^{n} (a_{rs} \dot{q}_s + b_r) \right\} \\ &= \sum_{s=1}^{n} a_{rs} \ddot{q}_s + \sum_{j,s=1}^{n} \frac{\partial a_{rs}}{\partial q_j} \dot{q}_j \dot{q}_s + \sum_{s=1}^{n} \frac{\partial a_{rs}}{\partial t} \dot{q}_s \\ &\quad + \sum_{s=1}^{n} \frac{\partial b_r}{\partial t} + \sum_{j=1}^{n} \frac{\partial b_r}{\partial q_j} \dot{q}_j, \end{aligned} \qquad (5) $$

Expanding the second term we similarly get

$$ \frac{\partial T}{\partial q_r} = \frac{1}{2} \sum_{i,j=1}^{n} \frac{\partial a_{ij}}{\partial q_r} \dot{q}_i \dot{q}_j + \sum_{i=1}^{n} \frac{\partial b_i}{\partial q_r} \dot{q}_i + \frac{\partial c}{\partial q_r}. \qquad (6) $$

Denoting $\dot{q} := [q_1 \ q_2 \ q_3 \ \dots \ q_n]^T$, and substituting expressions (5) and (6) in relation (4), we obtain Lagrange's equations as

$$ A \ddot{q} + f(q, \dot{q}, t) = Q + Q', \quad q(0) = q_0, \quad \dot{q}(0) = \dot{q}_0 \qquad (7) $$

where the vector function $f$ is in general a nonlinear function of its arguments and $Q := [Q_1 \ Q_2 \ \dots \ Q_n]$. The vector $Q'$ is similarly defined. The $n \times n$ matrix $A$ is positive definite and symmetric, and is related to the inertial properties of the system (Rosenberg, 1972, pp. 202).

Given the generalized coordinates and the generalized velocities $q_r$ and $\dot{q}_r$, let $\ddot{q}_r$ be any kinematically admissible acceleration which satisfies the $p$ nonholonomic constraints given by equation set (3). Thus, the set $(\ddot{q}_r, \dot{q}_r, q_r)$ satisfies the differential constraint equations

$$ \begin{aligned} \sum_{r=1}^{n} \alpha_{kr} \ddot{q}_r + \sum_{r=1}^{n} \sum_{s=1}^{n} \frac{\partial \alpha_{kr}}{\partial q_s} \dot{q}_s \dot{q}_r + \sum_{s=1}^{n} \frac{\partial \beta_{kt}}{\partial q_s} \dot{q}_s \\ + \sum_{r=1}^{n} \frac{\partial \alpha_{kr}}{\partial t} \dot{q}_r + \frac{\partial \beta_{kt}}{\partial t} = 0, \quad k=1, 2, \dots, p. \end{aligned} \qquad (8a) $$

Furthermore, if $\ddot{q}_r$ are the actual generalized accelerations of the mechanical system satisfying both the Lagrange equa-
tions and the constraints (3), then the set $(\ddot{q}_r, \dot{q}_r, q_r)$ satisfies the set of equations

$$ \sum_{r=1}^{n} \alpha_{kr} \ddot{q}_r + \sum_{r=1}^{n} \sum_{s=1}^{n} \frac{\partial \alpha_{kr}}{\partial q_s} \dot{q}_s \dot{q}_r + \sum_{s=1}^{n} \frac{\partial \beta_{kt}}{\partial q_s} \dot{q}_s + \sum_{r=1}^{n} \frac{\partial \alpha_{kr}}{\partial t} \dot{q}_r + \frac{\partial \beta_{kt}}{\partial t} = 0, k=1, 2, \dots, p. \quad (8b) $$

Subtracting the corresponding equations from the sets (8a) and (8b), we obtain

$$ \sum_{r=1}^{n} \alpha_{kr}(\ddot{q}'_r - \ddot{q}_r) = 0, \quad k=1, 2, \dots, p. \quad (9) $$

Thus,

$$ \delta \ddot{q}_r = \ddot{q}'_r - \ddot{q}_r, \quad r=1, 2, \dots, n, \quad (10) $$

are kinematically admissible instantaneous variations of the acceleration; i.e., virtual accelerations which are consistent with the nonholonomic constraints. Thus, by D'Alembert's principle, the total virtual work done by the forces of constraint, $Q'_r$, under these virtual accelerations equals zero. This requires

$$ \sum_{r=1}^{n} Q'_r \delta \ddot{q}_r = \sum_{r=1}^{n} Q'_r (\ddot{q}'_r - \ddot{q}_r) = [Q']^T (\ddot{q}' - \ddot{q}) = 0, \quad (11) $$

where in the last expression we have denoted $Q' = [Q'_1, Q'_2, \dots, Q'_n]^T$. Using Eq. (7), this in turn entails

$$ [A\ddot{q} + f(q, \dot{q}, t) - Q]^T (\ddot{q}' - \ddot{q}) = 0. \quad (12) $$

This condition is equivalent to

$$ \begin{aligned} [A\ddot{q}' + f(q, \dot{q}, t) - Q]^T A^{-1} [A\ddot{q}' + f(q, \dot{q}, t) - Q] \\ = [A\ddot{q} + f(q, \dot{q}, t) - Q]^T A^{-1} [A\ddot{q} + f(q, \dot{q}, t) - Q] \\ + (\ddot{q}'_r - \ddot{q}_r)^T A (\ddot{q}'_r - \ddot{q}_r). \end{aligned} \quad (13) $$

Since $A$ is positive definite, the second term on the right-hand side is always positive. Hence, we obtain the condition that the generalized accelerations, $\ddot{q}_r(t)$ of the constrained mechanical system are such as to minimize

$$ \|J\{\ddot{q}(t) | \dot{q}(t), q(t)\}\|_2^2 = \|A^{-1/2} [A\ddot{q}(t) + f(q, \dot{q}, t) - Q]\|_2^2, \quad (14) $$

at each instant of time $t$, while satisfying the set of constraints (3). We note that both $\ddot{q}_r(t)$ and $q(t)$ are known at time $t$.

We have thus reduced the problem of the determination of the evolution of a mechanical system subjected to given forces to that of solving a constrained quadratic minimization problem at each instant of time. We note that while the expression $J$ in (14) may be nonlinear in terms of the generalized coordinates and velocities, it is always linear in the unknown accelerations.

We point out here that were we to have chosen rectangular coordinates $(x_i, y_i, z_i)$, $i=1, 2, \dots, n$, for the $3n$ degrees-of-freedom of a discrete system of $n$ masses $m_i$, $i=1, 2, \dots, n$ then expression (14) would reduce to

$$ \sum_{r=1}^{n} \left[ m_r \left\{ \ddot{x}_r - \frac{X_r}{m_r} \right\}^2 + m_r \left\{ \ddot{y}_r - \frac{Y_r}{m_r} \right\}^2 + m_r \left\{ \ddot{z}_r - \frac{Z_r}{m_r} \right\}^2 \right], \quad (15) $$

where $X_r, Y_r, Z_r$ refer to the $x, y, z$-components of the impressed forces on mass $m_r$. This expression was first ennuciated in words by Gauss (1829) who called it the "constraint," thereby ennuciating the "principle of minimum constraint" (Whitaker, 1917). It was further elaborated on by Hertz (1917) and many others.

Since the forces of constraint $Q'$ satisfy Eq. (7), the minimization in (14) can also be expressed as

$$ \text{Minimize } \{ \|A^{-1/2} Q' \|_2^2 \}. \quad (16) $$

Thus, we see that for the mechanical system to satisfy a given set of constraints at each instant of time, Gauss's principle requires that at each instant the norm of the constraint forces $Q'$ weighted with respect to $A^{-1/2}$ be minimized—hence the name, the principle of minimum constraint. In addition, the $p$ constraint equations need to be satisfied. We show in the next section how this minimization can be carried out explicitly and also the explicit expression that can be written down for the constraint forces.

### 3 Solution of the Constrained Quadratic Minimization Problem

Using Gauss's Principle, we have thus reduced problems in mechanics to finding the accelerations $\ddot{q}_r(t)$ at each instant, $t$, given $q(t)$ and $\dot{q}(t)$, so that we require to

$$ \text{Minimize } \{ \|A^{-1/2} [A\ddot{q}(t) + f(q, \dot{q}, t) - Q]\|_2^2 \}, \quad (17) $$

while satisfying the constraints (8b). These constraints are again linear in the accelerations and can be written in matrix form as

$$ \delta \ddot{q}(t) = g\{\dot{q}(t), q(t), t\} \quad (18) $$

where we have denoted by $D$ the $p \times n$ matrix $[\alpha]_{ij}$ and by $g$ the vector containing the remainder of the terms in Eq. (8b). The right-hand side of Eq. (18) is known. For convenience, we shall now drop reference to the independent variable $t$, remembering that Eqs. (17) and (18) need to be satisfied at each instant of time. The solution of the consistent set of Eqs. (18) is obtained as (Rao and Mitra, 1972)

$$ \ddot{q} = D^{-g} + (I - D^{-D})h \quad (19) $$

where the $n \times p$ matrix $D^{-}$ is any generalized inverse (g-inverse) of $D$ which satisfies the relation

$$ DD^{-}D = D. \quad (20) $$

The vector $h$ is arbitrary. Substituting relation (19) into relation (17), we obtain

$$ \text{Min}_h \|Hh - z\|_2^2 \quad (21) $$

where

$$ H = A^{1/2} (I - D^{-D}) := A^{1/2} \hat{H}, \quad (22) $$

and

$$ z = - \{ A^{1/2} D^{-g} + A^{-1/2} (f - Q) \}. \quad (23) $$

For brevity, we have dropped the arguments of the vector function $f$. We next obtain the solution, $h$, of the least squares problem (21) as (Rao and Mitra, 1972),

$$ h = H_{ls}^{-} z + (I - H_{ls}^{-} H) w, \quad (24) $$

where the matrix $H_{ls}^{-}$ is the generalized "least-squares" inverse defined as satisfying the relations

$$ (HH_{ls}^{-} H) = H, \quad (25) $$

and

$$ [HH_{ls}^{-}]^T = HH_{ls}^{-}. \quad (26) $$

The vector $w$ is again an arbitrary vector.

Using expression (24) in (19), we thus obtain the explicit solution to the constrained minimization problem given by Eqs. (17) and (18) as

$$ \begin{aligned} \ddot{q} &= D^{-g} + (I - D^{-D}) \{H_{ls}^{-} z + (I - H_{ls}^{-} H) w\} \quad (27a) \\ &= D^{-g} + (I - D^{-D}) H_{ls}^{-} z + (I - D^{-D}) w - (I - D^{-D}) H_{ls}^{-} Hw. \end{aligned} \quad (27b) $$

We now express the matrix $H$, as in Eq. (22), by $A^{1/2} \hat{H}$. The least squares $g$-inverse, $H_{ls}^{-}$, can now be expressed as

$$ H_{ls}^{-} := \hat{H}_{ls}^{-} A^{-1/2} = (I - D^{-D})_{ls}^{-} A^{-1/2}, \quad (28) $$
where

$$ \hat{H} \hat{H}_{ls}^{-} \hat{H} = \hat{H}, \quad \text{and} \quad [\hat{H} \hat{H}_{ls}]^T A = A [\hat{H} \hat{H}_{ls}]. \qquad (29) $$

Relations (29) follow from the definition of $H$ in Eq. (22) and that of $H_{ls}$ in relations (25) and (26). Using relation (28) in (27b), we observe that the last two terms on the right-hand side cancel out, yielding

$$ \begin{aligned} \ddot{q} &= D^{-}g + (I - D^{-}D)H_{ls}z \\ &= D^{-}g - (I - D^{-}D)(I - D^{-}D)_{ls} \{D^{-}g + A^{-1}(f - Q)\}. \end{aligned} \qquad (30a) \qquad (30b) $$

We have thus obtained explicit expressions for the accelerations at time $t$ for the constrained motion, given the generalized coordinates and velocities.

Heuristically speaking, the inverse, $D^{-}$, comes about (and is related to) to the constraints on the mechanical system, while the inverse, $H_{ls}$, comes about (and is related to) the motion of the dynamical system. These two inverses employed in expression (30a) are in general different in nature from each other (Rao and Mitra, 1972). The matrix $D^{-}$ refers to any g-inverse of $D$ (i.e., satisfying relation (20)), while the matrix $H_{ls}^{-}$ refers to any “least-squares g-inverse” of $H$ (satisfying relations (25) and (26)). From a practical standpoint, the flexibility of choosing any g-inverse belonging to the requisite classes stated above is a useful advantage in obtaining the equations of motion of complex dynamical systems, for, depending on the situation at hand, certain g-inverses are easier to determine than others.

For example, a possible candidate for $H_{ls}^{-}$ might be $[H^T H]^{-} H^T$. Another might be the often-used Moore-Penrose (MP) inverse. However, the MP inverse is an element of both the set of g-inverses and the set of least squares g-inverses, and is a useful candidate because of the several computer codes available for its ready determination. The MP inverse can thus be used for both $D^{-}$ and $H_{ls}^{-}$ in Eq. (30a).

Even if the matrix $D$ has rank $p_1 < p$, we are assured (Lawson and Hanson, 1974) that the accelerations thus obtained are unique because the matrix $A$ is of rank $n$. However, when $p_1 < p$, the equation set (18) may not be consistent (Dahlquist and Bjorck, 1974) for all right-hand sides $g$. When the rank of $D$ is $p$, a unique solution to the constrained minimization problem exists for all $g$.

## 4 Explicit Form of Constrained Equations of Motion and the Constraint Forces

We have obtained in Eq. (30) an explicit expression for the generalized accelerations at time $t$, given the generalized coordinates and the velocities. Using this, we can therefore express the constrained equations of motion, valid at any time, $t$, for a general system in first-order form, as

$$ \begin{aligned} \frac{d}{dt} \begin{bmatrix} q \\ \dot{q} \end{bmatrix} &= \begin{bmatrix} \dot{q} \\ D^{-}g + (I - D^{-}D)H_{ls}^{-}z \end{bmatrix}, & \text{with } q(0) \\ &= q_0, \quad \dot{q}(0) = \dot{q}_0. \end{aligned} \qquad (31) $$

The quantity $z$ is defined in Eq. (23). We note this explicit set of equations for the system include the effects of the Pfaffian constraints. They can therefore be thought of as the new equivalent equations of motion; they constitute a generalization of the equations found in Wang and Huston (1989).

The equation set (31), which in general will be nonlinear, can now be numerically solved using any of the standard numerical integration schemes, such as the fourth-order Runge-Kutta method, or other methods like the predictor corrector methods. The right-hand side of Eq. (31) guarantees that the accelerations satisfy both the constraints and Gauss's Principle simultaneously at each instant of time.

Furthermore, comparing Eqs. (7) and (30a), the forces of constraint can also be explicitly written as

$$ Q' = (I - AXA^{-1})[f(q, \dot{q}, t) - Q] + A(I - X)D^{-}g, \qquad (32) $$

where we have denoted by the matrix $X$ the quantity $(I - D^{-}D)H_{ls}^{-}A^{-1/2}$. Often the constraints require that the system follows a given trajectory in configuration space. The constraint forces can then be thought of as the control forces necessary to cause the system to follow this particular trajectory.

The satisfaction of the constraint equations at each instant of time for the mechanical system entails the development of constraint forces which, at each instant of time, satisfy Eq. (16); the constraint forces $Q'$, by Gauss's Principle, must therefore minimize $\|A^{-1/2}Q'\|^2$. We note the the constraint forces acting at any instant, thus require for their determination nothing other than the displacement and velocity information at that instant, along with information about the constraints, at that specific instant. These forces of constrain can hence be determined at each time instant as the system's dynamics evolve. This makes the approach useful in real-time control, especially when the complete constrained trajectory is not known a priori. Thus, use of relation (32) may be made in the determination of real-time control required for tracking a given trajectory.

Equation (32) also shows that, in general, the control force vector, $Q'$, is dependent on $q$, $\dot{q}$, and $t$, and therefore constitutes closed-loop control. We note that the elements of matrices $D$ and $A$ depend on the coordinates $q_i$ and time. Similarly, elements of the vector $g$ depend on $q_i$, $\dot{q}_i$ and time. In certain special situations the elements of $D$, $A$, and $g$ may depend solely on time; then, the second term on the right-hand side of Eq. (32) is not dependent on $q_i$, and $\dot{q}_i$ and may be thought of as the feed-forward component of the total control force which is required to generate the constrained motion. Thus, in this special case, the control force may be thought of as being composed of a feedback control force (the first term of Eq. (32)) and a feed-forward control force.

## 5 Numerical Examples

In this section we consider three examples, the last two of which are numerical. The first deals with the constrained motion of a particle free of any “given” forces, where the constraint is nonholonomic. The second deals with the large amplitude motion of a planar pendulum. The pendulum bob is provided with two degrees-of-freedom and a constraint relation is provided on the length of the pendulum. We first constrain the length of the pendulum to be a constant. We use this example as a base line to check our results with the direct use of Runge-Kutta integration where the angle coordinate is used to preserve the length constraint. We next consider a more general, nonlinear constraint on the length of the pendulum. The third example is related to the problem of controlling a dynamical system (e.g., a machine tool) so that it follows a given trajectory. Here we show the ease with which the feedback tracking control force can be obtained using Eq. (31). In all the computations, a variable-step Runge-Kutta integration scheme is used with a local error tolerance of $10^{-10}$. The Moore-Penrose inverse is used for each of the generalized inverses in Eqs. (31) and (32).

(1) Consider the motion of a particle of unit mass, free of any “given” forces, moving in three-dimensional euclidean space ($q_1 = x, q_2 = y, q_3 = z$). Let the particle be subjected to the nonholonomic, catastrophic constraint

$$ \dot{y} = z \dot{x}. \qquad (36) $$

At time $t = 0$, the initial conditions of the particle are compatible with this constraint. We want to find the equations of motion for the particle for $t \ge 0$.

The system has two degrees-of-freedom; yet, the nonholonomic nature of the constraint requires three coordinates for a specification of the system's configuration. This example is taken from Rosenberg (1972, p. 204). Since there are no “given”
<figure><img src="image_2.png" /><figcaption>Fig. 1 Response of a pendulum of constant length to large amplitude motions using the new constrained equations of motion</figcaption></figure>

<caption>Fig. 2 Response of a pendulum with nonlinear constraints with $\mu = -0.1$ and $a = 4.0$</caption>

forces, $Q=f=0$. The matrix $D = [-z \ 1 \ 0]$ here, and the matrix $D^- = \{1/(1+z^2)\}[1-z \ 1 \ 0]^T$. The matrix $A$ and the vector $g$ are scalars equal to unity and $\dot{z}\dot{x}$, respectively. The vectors $f$ and $Q$ are zero, and the $MP$ inverse of $(I - D^-D)$, obtained by using full rank factorization, is given by

$$ (I - D^{-}D)^{+} = \begin{bmatrix} \frac{1}{(1+z^2)} & \frac{z}{(1+z^2)} & 0 \\ \frac{z}{(1+z^2)} & \frac{z^2}{(1+z^2)} & 0 \\ 0 & 0 & 1 \end{bmatrix}. \quad (37) $$

Using Eq. (29), we obtain explicitly the constrained equations of motion as

$$ \begin{bmatrix} \ddot{x} \\ \ddot{y} \\ \ddot{z} \end{bmatrix} = \frac{\dot{z}x}{(1+z^2)} \begin{bmatrix} -z \\ 1 \\ 0 \end{bmatrix} \quad (38) $$

which are, of course, the same as Rosenberg's result.

(2) Using rectangular axes in an inertial frame of reference, we can express the motion of the pendulum bob in the $x$ and $y$-directions ($y$ taken downwards) as

$$ m\ddot{x} = 0, \quad m\ddot{y} = g; \quad x(0) = \dot{y}(0) = 0, \quad \dot{x}(0) = a, \quad y(0) = L, \quad (33) $$

along with the constraint given by

$$ x^2 + y^2 = L^2 \{1 - \mu \sin(x^2)\}. \quad (34) $$

The value of $m$ is taken to be unity, that of $L$ was chosen to be 20 units, and that of $g$ to be $4\pi$. Figures 1(a–d) show the results obtained by using Eq. (31), for $a=20$ units, with $\mu = 0$. The fourth-order Runge-Kutta scheme (RK) was used. The oscillations are in the nonlinear range; the maximum angle made by the pendulum bob with the vertical during the oscillation being about 80 degrees. These results are the same (to within the error tolerance) as those obtained using the direct Runge-Kutta (RK) integration taking the angle of rotation (about the vertical) of the pendulum bob as the generalized coordinate. Figures 1(e) and 1(f) show the phase plots in the $x$-$\dot{x}$ and the $y$-$\dot{y}$ planes. Figure 1(g) shows the extent to which the constraint is satisfied during the numerical computations. Here, the error is defined as the difference between the left-hand side and the right-hand side of Eq. (34), a quantity which should theoretically be zero.

Figures 2(a–d) show the response of the system when $a = 4$ and $\mu = -0.1$. We again integrate Eqs. (31) using the fourth-order RK method. Figures 2(e) and 2(f) show the phase plots. As before, Fig. 2(g) shows the error in satisfying the constraint.

Figure 3 shows the phase trajectories of the same system starting with different initial velocities, $a = 0.5, 1, 2$, and $4$.

(3) The third example deals with a coupled, damped Duffing's oscillator described by the equations

$$ m_1\ddot{x}_1 + k_1(x_1 - x_2) + c_1(\dot{x}_1 - \dot{x}_2) + \bar{k}_1(x_1 - x_2)^3 = 0 \quad (35a) $$

$$ m_2\ddot{x}_2 + k_2x_2 - k_1(x_1 - x_2) + c_2\ddot{x}_2 - c_1(\dot{x}_1 - \dot{x}_2) + \bar{k}_2x_2^3 - \bar{k}_1(x_1 - x_2)^3 = 0, \quad (35b) $$
<figure><img src="image_3.png" /><figcaption>Fig. 3 The phase plots for $\mu = -0.1$. Figs. 3(a) and 3(b) show the response for $a = 0.5$; Figs. 3(c) and 3(d) show the response for $a = 1.0$; Figs. 3(e) and 3(f) show the response for $a = 2$; and, Figs. 3(g) and 3(h) show the response for $a = 4.0$.</figcaption></figure>

with

$x_1(0) = a, x_2(0) = b, \dot{x}_1(0) = c, \text{ and } \dot{x}_1(0) = d$.

We aim to determine the tracking control forces required so that the relative displacements of the masses are constrained to follow the exponentially decaying sinusoidal trajectory given by

$$x_1(t) - x_2(t) = Ae^{-\alpha t}\sin(\omega t). \quad (35c)$$

The parameters describing the system and its constraints are

$m_1 = 2, m_2 = 1, k_1 = 10, k_2 = 12, \bar{k}_1 = 1, \bar{k}_2 = 2, c_1 = 0.1, c_2 = 0.15, A = 1, \omega = 2\pi$.

The initial conditions are taken to be $a = 1, b = 1$, and $d = 2$. We note that the initial conditions must satisfy the constraints and hence the parameter $c$ is determined from $d$ and Eq. (35c). Figures 4(a) and 4(b) show the time histories of the displacement and velocity, obtained by integrating Eq. (31), when the parameter $\alpha$ equals 4 in Eq. (35c). The solid lines show quantities relevant to the “1” coordinate (i.e., to mass $m_1$), and the dashed lines show quantities relevant to the “2” coordinate (i.e., to mass $m_2$). Figure 4(c) shows the control forces ($f_1$ and $f_2$), calculated using Eq. (32), required to be applied to masses $m_1$ and $m_2$, respectively, to track this trajectory appropriately in configuration space. Figure 4(d) shows the computed value of $[x_1(t) - x_2(t)]$. We find (see, Fig. 5(a)) that the constraint is tracked to within an error (i.e.,

<figure><img src="image_4.png" /><figcaption>Fig. 4 Figures 4(a) and 4(b) shows the constrained response of the coupled Duffing's Oscillator, with $\alpha = 4.0$; Fig. 4(c) shows the two components of the control force needed to have the system follow this constrained trajectory; Fig. 4(d) shows $[x_1(t) - x_2(t)]$. Figures 4(e-h) show similar results when the constraint is enforced with $\alpha = 0.4$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 5 (a) position tracking error for $\alpha = 4$; (b) position tracking error for $\alpha = 0.4$</figcaption></figure>

LHS-RHS of Eq. 35(c)) of $10^{-9}$, a number consistent with the local error tolerance level of $10^{-10}$ used for the RK integration. Figures 4(e-h) show similar results when the value of $\alpha$ in Eq. (35c) is now taken to be 0.4. All other parameter values are left unchanged. As before, the solid lines show quantities relevant to the “1” coordinate and the dashed lines show quantities relevant to the “2” coordinate. The error in tracking this trajectory is shown in Fig. 5(b), and is again found to be of the order of $10^{-9}$.

# 6 Conclusions and Discussion

This paper deals with discrete, dynamical systems which are
subjected to Pfaffian, holonomic, and nonholonomic con-
straints. By using Gauss's Principle of Least Constraint, we
have recast the Lagrange equations which describe the con-
strained motion of such a system in the form of a quadratic
programming problem. In this paper we present an explicit
analytical solution to this constrained quadratic minimization
problem and thereby obtain an exact and explicit set of equa-
tions that describe the constrained motion of discrete, me-
chanical systems subjected to Pfaffian constraints. These
equations can be numerically integrated using standard nu-
merical techniques like the Runge-Kutta and Predictor-Cor-
rector methods.

We summarize our finding as follows:

1 We have used Gauss's Principle to obtain a new concep-
tualization of the equations of motion of a constrained, discrete
dynamic system. The equations of motion that we develop,
directly yield the time evolution of the system; we do away
with the need to use any Lagrange multipliers. Thus Gauss's
principle, though largely neglected by the mechanicians of this
century, is shown to yield significant insights into the dynamics
of constrained systems.

2 Besides its aesthetic appeal, the method proposed herein has
special advantages when working with nonintegrable con-
straints. In fact it does away with the somewhat unnecessary
categorization of Pfaffian constraints into: (a) holonomic and
nonholonomic constraints and (b) rheonomic and scleronomic
constraints—the approach being able to handle all these types
of equality constraints with equal ease. The equations devel-
oped here can be used in situations where the derivation of
the equations of motion (using Lagrange multipliers) for con-
strained systems may become cumbersome and/or difficult to
implement computationally. The approach thus provides a
conceptual and practical simplicity in the formulation of the
equations of motion of complex mechanical systems, because
the constraints can be explicitly handled as additional equations
whose effect can be directly incorporated in the equations of
motion. The flexibility that this formulation affords in the
specific choices of the generalized inverses $D^{-}$ and $H_{is}$, is an
added feature which is new, and particularly helpful from a
practical standpoint. Response sensitivity studies related to
altering the constraints can thus be easily carried out.

3 Even for systems, where it may be possible to eliminate
certain variables directly from the equations of motion, the
method provides a direct and more aesthetic approach by not
favoring any particular subset of coordinates over any other.

4 The explicit expressions obtained for the constraint forces
may be used to advantage when dealing with the determination
of control forces required to control a system so that it follows
a certain trajectory in configuration space, or more generally,
satisfies a given set of Pfaffian constraints. Such problems
arise in many areas of application, like position tracking of

machine tools (Tomizuka, 1987) and robotic manipulator con-
trol (Seraji, 1987).

Furthermore, we obtain the additional insight from Gauss's
Principle that for the system to satisfy the constraint equations
*at each instant* of time, a *specific* quadratic function of the
constraint forces, namely, $Q'^T A^{-1} Q'$, must be minimized at
each instant of time. This sheds light on the reason why least-
squares formulations of the tracking control problem have
often *not* led to proper trajectory tracking when minimizing
the integrals of 'general quadratic functions of the control
forces.

5 The three examples considered here illustrate that the approach may be useful in answering the two commonly occurring problems in particle mechanics (Rosenberg, 1972): (a) finding the response of mechanical systems subjected to general types of time-dependent, Pfaffian equality constraints and (b) finding the control forces required to be imposed on a system, in real-time, so that it satisfies a given set of holonomic or nonholonomic Pfaffian constraints. Our third example shows that by using the new set of dynamical equations obtained herein, the accuracy with which the system is led to follow a constrained trajectory can indeed be high.

References

Appell, P., 1925, *Sur une forme generale des equations de la dynamique*,
Memorial des Sciences Mathem. Paris, Gauthier-Villars.

Arnold, V. I., 1980, *Mathematical Methods of Classical Mechanics*, Springer-
Verlag, New York.

Dahlquist, G., and Bjorck, A., 1974, *Numerical Methods*, Prentice-Hall,
Englewood Cliffs, N.J.

Hertz, H., 1917, *Principles of Mechanics*, Macmillan, London.

Gauss, C., 1829, "Über ein neues allgemeines Grundgesetz der Mechanik," *Journal fur Reine und Angewandte Mathematik*, pp. 232-235.

Lanczos, C., 1970, *The Variational Principles of Mechanics*, Dover Publi-
cations, New York.

Lawson, C., and Hanson, R., 1974, *Solving Least Squares Problems*, Prentice-
Hall, Englewood Cliffs, N.J.

Pars, L. A., 1979, *A Treatise on Analytical Dynamics*, Ox Bow Press, Con-
necticut.

Rao, C., and Mitra, S., 1972, *Generalized Inverse of Matrices and Its Ap-
plications*, John Wiley and Sons, New York.

Rosenberg, R., 1972, *Analytical Dynamics of Discrete Systems*, Plenum Press.

Seraji, H., 1987, "Design of Feedforward Controllers for Multi-variable
Plants," *International Journal of Control*, Vol. 46, pp. 1633-1651.

Shan, D., 1975, "Canonical Equations for Mechanical Systems With Second-
Order Linear Nonholonomic Constraints," *Prikladnaya Mekhanika*, Vol. 11,
No. 2, pp. 58-67.

Synge, J. L., 1926, "On the Geometry of Dynamics," *Proceedings of the
Royal Society, London*, Vol. A226, pp. 31-106.

Tomizuka, M., 1987, "Zero Phase Error Tracking Algorithm for Digital
Control," ASME *Journal of Dynamic Systems, Measurement and Control*, Vol.
109, pp. 65-68.

Wang, J. T., and Huston, R. L., 1989, "A Comparison of Analysis Methods
of Redundant Multibody Systems," *Mechanics Research Communications*, Vol.
16, pp. 175-182.

Whittaker, E., 1917, *A Treatise on the Analytical Dynamics of Particles and
Rigid Bodies*, Cambridge, Univ. Press, Cambridge, U.K.