Dynamical Systems via

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.


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 ann-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 ofleast 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 3n-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) 2 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-ofconstraint 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 nonspecific 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 belive, is not available in the current literature (e.g., in Whittaker (1917), Synge (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 illus- trates 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.

Gauss's Principle
Consider a holonomic mechanical system with n-degrees-offreedom whose generalized coordinates are q 1 , q 2 , q 3 , • • , q 11 • The Lagrange equations describing the motion of the system may be written as where T denotes the kinetic energy and Q, is the generalized impressed force.The kinetic energy can be expressed as where, in general, the au and b; and c are functions of the generalized coordinates and time.
Assume now that the system is subjected to an additional where 01k, and f3kt are functions of the generalized coordinates and time.We note that the constraints may be scleronomic or rheonomic, catastatic or a catastatic (Rosenberg, 1972).These p constraints may be thought of as imposing additional constraint forces, Q/ , on our system, thereby altering the set of Eqs.
+ LJ -+ L J -q;. s ~l at J ~l aq; Expanding the second term we similarly get (5) Denoting q: = [q1 Q2 q 3 .. q 11 f, and substituting expressions (5) and (6) in relation (4), we obtain Lagrange's equations as where the vector function f is in general a nonlinear function of its arguments and Q The vector Q' is similarly defined.The n x 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, and q" let q; be any kinematically admissible ac- celeration which satisfies the p nonholonomic constraints given by equation set (3).Thus, the set (q;, q" q,) satisfies the differential constraint equations Furthermore, if q, are the actual generalized accelerations of the mechanical system satisfying both the Lagrange equa-tions and the constraints (3), then the set (ij, q, q,) satisfies the set of equations Qs s~l Qs f., aCikr .a{3kt .
Subtracting the corresponding equations from the sets (Sa) and (8b), we obtain n 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 -l/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.
~ Cikr(ij: -ij,) =0, k= 1, 2, ... , p. 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;, under these virtual accelerations equals zero.This requires where in the last expression we have denoted Q' Q; .. g;,f.Using Eq. ( 7), this in turn entails This condition is equivalent to +(q; -q,)TA(q: -q,).(13 Since A is positive definite, the second term on the righthand side is always positive.Hence, we obtain the condition that the generalized accelerations, q,(t) of the constrained mechanical system are such as to minimize IIJ(q(t) I q (t), q(t) )II~:= IIA -112 [Aq(t) +f(q, q, t)-Q]II~, (14) at each instant of timet, while satisfying the set of constraints (3).We note that both q,(t) and q(t) are known at timet.
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 (xi> Yi, Zi), i = 1, 2, .. , n, for the 3n degreesof-freedom of a discrete system of n masses mi, i = 1, 2, .. , n then expression ( 14) would reduce to where X, Y, Z, refer to the x, y, z-componems of the impressed forces on mass m,.This expression was first ennuciated in words by Gauss ( 1829) who called it the'' constraint,'' thereby ennuciating the "principle of minimum constraint" (Whittaker, 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 ij Using Gauss's Principle, we have thus reduced problems in mechanics to finding the accelerations q,(t) at each instant, t, given q (t) and q (t), so that we require to Minimize (IIA-112 [Aq(t)+f(q, q, t)-Q]II~), (17) ij while satisfying the constraints (8b).These constraints are again linear in the accelerations and can be written in matrix form as where we have denoted by D the p X n matrix [c;_]u 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) where the n X p matrix nis any generalized inverse (ginverse) of D which satisfies the relation The vector his arbitrary.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), where the matrix H,-sis the generalized ''least-squares'' inverse defined as satisfying the relations ( 25) and ( 26) The vector w is again an arbitrary vector.l)sing expression (24) in ( 19), we thus obtain the explicit solution to the constrained minimization problem given by Eqs. ( 17) and ( 18) as (27a) We now express the matrix H, as in Eq. ( 22 (29) Relations ( 29) follow from the definition of H in Eq. ( 22) and that of H 1 ~ in relations ( 25) and ( 26).Using relation ( 28) in (27 b), we observe that the last two terms on the right-hand side cancel out, yielding 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 1 ~, 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 ginverse of D (i.e., satisfying relation ( 20)), while the matrix H 1 ~ 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{; might be (HTH]-HT.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 1 ~ in Eq. (30a).
Even if the matrixD has rankp 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 PI < 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.

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 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 where we have denoted by the matrix X the quantity (I -D-D)H 1 ~A 112 .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~,, by Guass's Principle, must therefore minimize IIA -112 Q' II .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, q, and t, and therefore constitutes closed-loop control.We note that the elements of matrices D and A depend on the coordinates q; and time.Similarly, elements of the vector g depend on q;, q; 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;, and q; 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.

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, fr_ee 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, catastatic constraint y=zx.
(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 ~ 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" The value of m is taken to be unity, that of L was chosen to be 20 units, and that of g to be 47!'.Figures 1 (a-d) show the results obtained by using Eq. ( 31), for a=20 units, with J.t = 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 I (f) show the phase plots in the x-x and they-y planes.Figure I (g) shows the extent to which the constraint is satisfied during the numerical computations.
Here, the error is defined as the difference between the lefthand 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 J.t = -0.1.We again integrate Eqs.(31) using the fourthorder RK method.Figures 2 (e) and 2 (f) show the phase plots.A:s 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 , I, 2, and 4.
We aim to determine the tracking contr9l forces required so that the relative displacements of the masses are constrained to follow the exponentially decaying sinusoidal trajectory given by (35c) The parameters describing the system and its constraints are 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 IX 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 (j 1 and j 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 J ------------------: 0 Time  -h) show similar results when the value of IX 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 a~d the dashed lin~s show quantities relevant to the "2" coordmate.The error m 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 constraints.By using Gauss's Principle of Least Constraint, we have recast the Lagrange equations which describe the constrained 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 equations that describe the constrained motion of discrete, mechanical systems subjected to Pfaffian constraints.These equations can be numerically integrated using standard numerical techniques like the Runge-Kutta and Predictor-Corrector methods.We summarize our finding as follows: 1 We have used Gauss's Principle to obtain a new conceptualization 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 constraints.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 developed here can be used in situations where the derivation of the equations of motion (using Lagrange multipliers) for constrained 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 nand H 1 ;, 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 control (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 -iQ', must be minimized at each instant of time.This sheds light on the reason why leastsquares 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.

Fig. 1 Fig. 2
Fig.1Response of a pendulum of constant length to large amplitude motions using the new constrained equations of motion