Proc. R. Soc. A (2008) 464, 2341–2363 doi:10.1098/rspa.2008.0040 Published online 24 April 2008

Optimal tracking control of nonlinear dynamical systems B Y F IRDAUS E. U DWADIA * ,1,2,3,4,5 1

Department of Aerospace and Mechanical Engineering, 2Department of Civil and Environmental Engineering, 3Department of Mathematics, 4 Systems Architecture Engineering, and 5Information and Operations Management, University of Southern California, 430K Olin Hall, Los Angeles, CA 90089-1453, USA

This paper presents a simple methodology for obtaining the entire set of continuous controllers that cause a nonlinear dynamical system to exactly track a given trajectory. The trajectory is provided as a set of algebraic and/or differential equations that may or may not be explicitly dependent on time. Closed-form results are also provided for the realtime optimal control of such systems when the control cost to be minimized is any given weighted norm of the control, and the minimization is done not just of the integral of this norm over a span of time but also at each instant of time. The method provided is inspired by results from analytical dynamics and the close connection between nonlinear control and analytical dynamics is explored. The paper progressively moves from mechanical systems that are described by the second-order differential equations of Newton and/or Lagrange to the ﬁrst-order equations of Poincare´, and then on to general ﬁrst-order nonlinear dynamical systems. A numerical example illustrates the methodology. Keywords: exact control; tracking control; nonlinear systems

1. Introduction The development of controllers for nonlinear mechanical systems has been an area of intense research over the last two decades or so. Many controllers that have been developed for trajectory tracking of complex nonlinear and multi-body systems rely on some approximations and/or linearizations (Slotine & Li 1991; Sastry 1999; Naidu 2003). Most control designs restrict controllers for nonlinear systems to be afﬁne in the control inputs (Brogliato et al. 2007). Often, the system equations are linearized about the system’s nominal trajectory and then the linearized equations are used along with various results from the well-developed theories of linear control. While this often works well in many situations, there are some situations in which better controllers may be needed. This is especially so when highly accurate trajectory tracking is required to be done in real time on systems that are highly nonlinear. Examples include the exact trajectory control of orbital, attitudinal and elastic motions of a multi-body spacecraft system that is required to perform precision tumbling. *Address for correspondence: Department of Aerospace and Mechanical Engineering, University of Southern California, 430K Olin Hall, Los Angeles, CA 90089-1453, USA ([email protected]). Received 28 January 2008 Accepted 28 March 2008

2341

This journal is q 2008 The Royal Society

2342

F. E. Udwadia

To place this paper within the context of the enormous literature that has been generated in the area of tracking control of nonlinear systems and to highlight what is new in it, we provide a brief review of the methods that have been developed so far and that have been applied to numerous areas of application ranging from chemical process control to robotics. Brogliato et al. (2007) provide an exhaustive review of the methods that have been developed to date for the tracking control of systems along with over 500 references. They point out that methods developed to date rely heavily on PID-type control and, most often, a linear feedback is provided to track a given trajectory. Chaou & Chang (2004) also provide recent developments on trajectory tracking and deal with the same basic theme (linear feedback) along with numerous applications. The optimality criterion considered in the literature to date is the minimization of the control cost integrated over a suitable span of time. In the robotics literature (Brogliato et al. 2007), trajectory tracking using inverse dynamics and model reference control has been used for some time now, and the methods developed therein can be seen as particular subclasses of the formulation discussed in the present work. Nonlinear control methods using controlled Lagrangian and Hamiltonians have also been explored along with passivity theory (Brockett 1977). These methods limit the structure of the control for a mechanical system to be a nonlinear function of its generalized displacement and they usually do not address the issue of control optimality. No such assumptions are made in this paper. Trajectory tracking in the adaptive control context (which is not the subject of this paper) has also been explored together with speciﬁc parametrizations to guarantee linearity in the unknown parameters of a system (Sadegh 1990). Thus, the methods used to date primarily rely on linearizations and/or PID-type control, and they posit assumptions on the structure of the control effort. By contrast, this paper takes a widely different approach that is based on recent results from analytical dynamics. Here the complete nonlinear problem is addressed with no assumptions on the type of controller that is to be used, except that it be continuous. In particular, we do not posit that the control is afﬁne in the inputs, we do not use linearized equations of motion about the desired trajectory, nor do we assume any ‘feedback structure’ to the nonlinear control effort. Furthermore, the optimality criterion used is the minimization of the control cost at each instant of time. As far as is known, the results provided here yield new and explicit methods for the control of highly nonlinear systems. Moreover, we provide for the ﬁrst time the entire class of nonlinear Lipschitzcontinuous (LC) controllers that would track the desired trajectory of the system, and we identify from among these controllers the one that minimizes, at each instant of time, a speciﬁed weighted norm of the control effort. Since we illustrate the power of the technique developed herein to ensure the tracking of a Rossler chaotic system by another Lorenz chaotic system, it is appropriate to indicate the state of the art in synchronization of nonlinear chaotic systems. Nearly all the work that has been done to date on the synchronization of chaotic systems deals with synchronization of identical chaotic oscillators that start from different initial conditions. This synchronization is performed by linear feedback between one or more of the phase states of the system (see Chen 2002; Lei et al. 2005). A recent monograph (with 350 references) points this out in detail (Boccaletti et al. 2002). It is only very Proc. R. Soc. A (2008)

Tracking control of nonlinear systems

2343

recently that the synchronization between two non-identical, low-dimensional chaotic systems has begun to be investigated (Boccaletti et al. 2002). One of the few papers on this is by Pyras (1996), which uses linear feedback. Rulkov et al. (1995) showed that this type of synchronization could exist. Such synchronization studies have resulted in the so-called imperfect phase synchronization (Zaks et al. 1999) in which the systems get ultimately de-synchronized by having intermittent phase slips, and lag synchronization (Rosenblum et al. 1997) which results in a time lag between the synchronized systems. Linear feedback has been used in all these studies, with the aim of synchronizing all the states of one chaotic system with those of the other. In our example we use non-identical chaotic systems. Our methodology does not assume any a priori structure on the control for the synchronization; no restrictions on the chaotic systems being of low dimensionality is required; as illustrated in the example, we can opt to synchronize (at will) one or more of the phase space variables; and, we can obtain the entire set of controllers in closed form that would do the job. Finally, among all these controllers, in the example we show that we can explicitly provide the controller that synchronizes and minimizes a weighted norm of the control effort at each instant of time. We begin by reformulating the trajectory control problem as a problem of constrained motion in the Lagrangian framework and use the underlying inspiration with which constrained motion in analytical dynamics is orchestrated by Nature (Udwadia & Kalaba 2002). We then expand and further develop this view by ﬁrst considering mechanical systems described by ﬁrst-order Poincare´ equations and then general nonlinear systems. Closed-form expressions for all the continuous controllers required for trajectory tracking for nonlinear systems that do not make approximations, either in describing the nonlinear system or in the nature of the nonlinear controller employed, are obtained. Such closed-form results appear to be new. Furthermore, no approximations or linearizations are made here with respect to the trajectory that is being tracked, which may be described in terms of nonlinear algebraic equations, nonlinear differential equations or a combination thereof; these descriptions could explicitly involve time also. Moreover, the approach arrives not just at one nonlinear controller for controlling a given nonlinear system, but also at the entire set of continuous controllers that would cause a given set of trajectory descriptions to be exactly satisﬁed. Furthermore, we show that when the cost function is the weighted norm of the generalized control input, its minimization can be done to yield the optimal controller that minimizes not just the integrated cost over a span of time, but also the cost at each instant of time. Explicit closed-form expressions for the optimal control are obtained. Section 2a of this paper begins with mechanical systems described by the second-order differential equations of the Newtonian and/or Lagrangian mechanics. Section 2b deals with Poincare´’s ﬁrst-order differential equations that describe the motion of mechanical systems; this takes us a step further towards general ﬁrst-order nonlinear dynamical systems. Section 3 deals with the close connection that the set of controllers developed in §2 have with the way Nature seems to orchestrate the control of mechanical systems subjected to trajectory requirements (constraints). The development of such controllers and their close correspondence with (i) Gauss’s principle of least constraint and (ii) the recently developed equations of motions for non-ideal constraints Proc. R. Soc. A (2008)

2344

F. E. Udwadia

are provided. Lessons from the way Nature seems to control mechanical systems are adduced. Section 4 deals with the exact control of general dynamical systems described by a set of nonlinear, non-autonomous ordinary differential equations. Section 5 deals with a numerical example that demonstrates the closed-form development of the control required to be applied to a nonlinear chaotic dynamical system so that it tracks the motions of another different nonlinear chaotic system. Section 6 concludes the paper with some remarks and observations.

2. Development of the entire set of controllers that cause a mechanical system to track a given trajectory (a ) Lagrangian and Newtonian descriptions Consider an unconstrained nonlinear mechanical system described by the secondorder differential equation of motion _ Z q_0 ; _ tÞ; ð2:1Þ M ðq; tÞ€ q Z Qðq; q; qð0Þ Z q0 ; qð0Þ where q(t) is the n-vector (n by 1 vector) of generalized coordinates; the dots indicate differentiation with respect to time; and the matrix M(q, t) is a positivedeﬁnite n by n matrix. Equation (2.1) can be obtained using either Newtonian or Lagrangian mechanics (Lagrange 1811; Hamel 1949; Goldstein 1976). The n-vector Q on the r.h.s. of equation (2.1) is a ‘known’ vector in the sense that it is a known function of its arguments. By ‘unconstrained’ we mean that the components of the initial velocity q_ 0 of the system can be independently assigned. We next require that this mechanical system be controlled so that it tracks a trajectory that is described by the following consistent set of m equations: fi ðq; tÞ Z 0; i Z 1; .; h ð2:2Þ and _ tÞ Z 0; ji ðq; q; i Z h C 1; .; m: ð2:3Þ We shall assume that the mechanical system’s initial conditions are such as to satisfy these relations at the initial time. The latter set of equations, which are non-integrable, is non-holonomic (Hamel 1949). In order to control the system so that it exactly tracks the required trajectory—i.e. satisﬁes equations (2.2) and (2.3)—we apply an appropriate _ tÞ, so that the equation of motion of the controlled control n-vector Q c ðq; q; system becomes _ tÞ C Q c ðq; q; tÞ; M ðq; tÞ€ q Z Qðq; q;

qð0Þ Z q0 ;

_ Z q_0 ; qð0Þ

ð2:4Þ

where now the components of the n-vectors q0 and q_0 satisfy equations (2.2) and (2.3) at the initial time, tZ0. Throughout this paper, we shall, for brevity, drop the arguments of the various quantities, unless needed for clarity. We begin by expressing equation (2.4) in terms of the weighted accelerations of the system. For any positive-deﬁnite n by n matrix N(q, t), we deﬁne the matrix Gðq; tÞ d½N 1=2 ðq; tÞM ðq; tÞ K1 Z M K1 ðq; tÞN K1=2 ðq; tÞ; Proc. R. Soc. A (2008)

ð2:5Þ

Tracking control of nonlinear systems

2345

and pre-multiplying equation (2.4) by N 1/2 (q, t), the ‘scaled’ equation, which we denote using the subscript ‘s’, is obtained as ð2:6Þ q€s Z as Cq€sc ; where € ð2:7Þ q€s dG K1q; a s dG K1 a Z ðN 1=2 M ÞðM K1 QÞ Z N 1=2 Q

ð2:8Þ

and ð2:9Þ q€sc dG K1q€c Z ðN 1=2 M ÞðM K1 Q c Þ Z N 1=2 Q c : In equation (2.8), we denote the acceleration of the uncontrolled system by _ tÞZ M K1 ðq; tÞQðq; q; _ tÞ: In equation (2.9), q€c ðq; q; _ tÞZ M K1 ðq; tÞQ c ðq; q; _ tÞ aðq; q; can be viewed as the deviation of the acceleration of the controlled system from that of the uncontrolled system. We now differentiate equation (2.2) twice with respect to time t, and equation (2.3) once with respect to time, giving the set of equations _ tÞ€ _ tÞ; Aðq; q; q Z bðq; q; ð2:10Þ where A is an m by n matrix of rank k and b is an m-vector. Noting equation (2.7), equation (2.10) can be further expressed as _ tÞ€ _ tÞ; Bs ðq; q; q s Z bðq; q; ð2:11Þ _ tÞ dAðq; q; _ tÞGðq; tÞ. We now express the n-vector q€s in terms of its where Bs ðq; q; orthogonal projections on the range space of BsT and the null space of Bs, so that q€s Z BsCBsq€s C ðI KBC q s: s Bs Þ€

ð2:12Þ

In equation (2.12), the matrix X denotes the Moore–Penrose (MP) generalized inverse of the matrix X (Moore 1910; Penrose 1955). It should be noted that equation (2.12) is a general identity that is always valid since it arises from the orthogonal partition of the identity matrix I Z BsCBs C ðI KBsCBs Þ. Using equation (2.11) in the ﬁrst member on the r.h.s. of equation (2.12), and equation (2.6) to replace q€s in the second member, we get ð2:13Þ q€s Z BsCb C ðI KBsCBs Þða s Cq€sc Þ; which, owing to equation (2.6), yields ð2:14Þ BsCBsq€sc Z BC s ðbK Bs a s Þ: But the general solution of the linear set of equations (2.14) is given by (Graybill 2001) C C C C C q€sc Z ðBC s Bs Þ Bs ðbK Bs a s Þ C ½I KðBs Bs Þ ðBs Bs Þz C

Z BsCðbK Bs a s Þ C ðI KBsCBs Þz;

ð2:15Þ

_ tÞ is any arbitrary n-vector. In the second equality where the n-vector zðq; q; above, we have used the property that ðBsCBs ÞCZ BsCBs in the two members on the r.h.s., along with the second MP-inverse property that BsCBs BsCZ BsC. Using equation (2.9), the explicit control force that exactly tracks the trajectory by exactly satisfying the given relations (2.2) and (2.3) is then explicitly given by K1=2 Q c Z N K1=2q€sc Z N K1=2 BC ðI KBsCBs Þz; s ðbK Bs a s Þ C N

Proc. R. Soc. A (2008)

ð2:16Þ

2346

F. E. Udwadia

_ tÞ dAðq; q; _ tÞGðq; tÞ. We may take zðq; q; _ tÞ to be C 1 (or, more where Bs ðq; q; generally, Lipschitz continuous (LC)) to ensure a unique solution to the system of equations (2.4). We hence have the following result. Result 2.1. Consider the mechanical system, which is described by the Lagrange (or Newtonian) equations of motion _ tÞ; M ðq; tÞ€ q Z Qðq; q;

ð2:17Þ

where M is an n by n positive-deﬁnite matrix and q is an n-vector. The system is required to exactly track the trajectory described by the equations fi ðq; tÞ Z 0;

i Z 1; .; h

ð2:18Þ

and _ tÞ Z 0; ji ðq; q;

i Z h C 1; .; m:

ð2:19Þ

The controlled system is described by the relation _ tÞ C Q c ðq; q; _ tÞ; M ðq; tÞ€ q Z Qðq; q;

ð2:20Þ

where Q c is the control. Assuming that the initial conditions of the mechanical system satisfy these _ tÞ (or controllers) trajectory requirements, the set of all possible controls Q c ðq; q; that causes the controlled system (2.20) to exactly track the required trajectory is explicitly given by Q c Z N K1=2 BsCðbKBs a s Þ C N K1=2 ðI KBC s Bs Þz;

ð2:21Þ

_ tÞ is any arbitrary n-vector whose components are continuously where zðq; q; differentiable—or LC—functions of its arguments; N(q, t) is any arbitrary n by n _ tÞZ Aðq; q; _ tÞ½N 1=2 ðq; tÞM ðq; tÞ K1 is positive-deﬁnite matrix a s Z N 1=2 Q; Bs ðq; q; _ tÞ is an m by n matrix of rank k, and bðq; q; _ tÞ is the an m by n matrix; Aðq; q; m-vector deﬁned in equation (2.10). _ tÞ given We can abbreviate the two components of the control vector Q c ðq; q; in equation (2.21) as _ tÞ Z Q1c ðq; q; _ tÞ C Q2c ðq; q; _ tÞ; Q c ðq; q;

ð2:22Þ

_ tÞ dN K1=2 BC Q1c ðq; q; s ðbK Bs a s Þ;

ð2:23Þ

_ tÞ dN K1=2 ðI KBsCBs Þz: Q2c ðq; q;

ð2:24Þ

where

and

Corollary 2.2. Relation (2.23) can also be given as Q1c Z N K1=2 ðAM K1 N K1=2 ÞCðbKAaÞ:

ð2:25Þ

Proof. Using the relation Bs Z AG Z AM N and equation (2.8), we have _ tÞZ M K1 ðq; q; _ tÞQðq; q; _ tÞ is the acceleration Bs a s Z AGG K1 aZ Aa, where aðq; q; of the uncontrolled system. The result then follows. & K1

Proc. R. Soc. A (2008)

K1=2

Tracking control of nonlinear systems

2347

Remark 2.3. Owing to the generality of the decomposition (2.13), relation (2.21) (alternatively, relations (2.22)–(2.24)) provides the entire set of continuous tracking controllers that cause the system to track the trajectory described by equations (2.18) and (2.19). Remark 2.4. The explicit closed-form tracking control obtained in equation (2.21), which causes the system to exactly track the given trajectory described by equations (2.18) and (2.19), does not contain any Lagrange multipliers, nor does our derivation invoke the notion of a Lagrange multiplier. We shall next show the following, somewhat remarkable, result. Result 2.5. We again consider the mechanical system described by the nonlinear Lagrange or Newtonian equation (2.17), which needs to be controlled _ tÞ, so that the trajectory through the addition of a control, n -vector Q c ðq; q; described by equations (2.18) and (2.19) is exactly tracked. Assuming that the system satisﬁes the trajectory requirements initially, the optimal controller that causes the system to (i) exactly track the required trajectory and (ii) minimize at each instant of time t, the cost _ tÞT N ðq; tÞQ c ðq; q; _ tÞ; J ðt Þ Z ½Q c ðq; q;

ð2:26Þ

for a given n by n positive-deﬁnite matrix N, is explicitly provided by K1=2 C _ tÞ Z N K1=2 BC Bs ðbKAaÞ; Q c ðq; q; s ðbKBs a s Þ Z N

ð2:27Þ

where Bs, b and as are as deﬁned before. Proof. Let us deﬁne the vector rðtÞ dN 1=2 Q c :

ð2:28Þ

We note that by equation (2.26), J ðtÞ Z ½Q c ðq; q; tÞT N ðq; tÞQ c ðq; q; tÞ Z krðt Þk2 :

ð2:29Þ

Then in view of equation (2.20), equation (2.28) can be rewritten as rðt Þ dN 1=2 ðM q€ KQÞ;

ð2:30Þ

so that we obtain, using equation (2.8), q€ Z ½N 1=2 M K1 ðr C N 1=2 QÞ Z Gðr C a s Þ:

ð2:31Þ

We assume that at the initial time the values of q0 and q_ 0 satisfy the trajectory requirements, and upon differentiating equations (2.18) and (2.19) we have, as before, Aq€ Z b:

ð2:32Þ

Using relation (2.31) and denoting BsZAG, equation (2.32) becomes Bs r Z bKBs a s : Proc. R. Soc. A (2008)

ð2:33Þ

2348

F. E. Udwadia

The solution of equation (2.33), subject to the condition that J ðtÞZ krðtÞk2 is a minimum, is then given by (Graybill 2001) rðt Þ Z BsCðbK Bs a s Þ;

ð2:34Þ

from which we obtain, by using relation (2.28), the explicit expression for the optimal control as _ tÞ Z ½N ðq; tÞ K1=2 BsCðq; q; _ tÞ½bðq; q; _ tÞK Bs ðq; q; _ tÞa s ðq; q; _ tÞ; Q c ðq; q;

ð2:35Þ

where we have written out the explicit result in extensio. The second equality in equation (2.27) follows from corollary 2.2. & Remark 2.6. For a given choice of the weighting matrix N, we have obtained the explicit closed-form expression for the exact full state controller that tracks the trajectory described by equations (2.18) and (2.19) and minimizes the cost function J ðtÞZ ½Q c T N ðq; tÞQ c , under the proviso that the initial conditions of the system satisfy the trajectory description. The minimum cost J(t) is given by 2 K1 K1=2 C J ðtÞ Z ½Q1c T NQ1c Z kBC N Þ ðbKAaÞk2 : s ðbKAaÞk Z kðAM

ð2:36Þ

Note, as before, that the explicit result given in equation (2.27) does not contain any Lagrange multipliers, nor does our derivation invoke anywhere the notion of a Lagrange multiplier. Remark 2.7. Result 2.5 and equation (2.23) show that for a given positivedeﬁnite matrix N(q, t) the optimal control that minimizes J(t) at each instant of _ tÞ h 0 in equation time is explicitly obtained in closed form by setting Q2c ðq; q; (2.22). One way (see corollary 2.8 below) of doing this would be by setting _ tÞ h 0 in equation (2.24) (or, in equation (2.21)). More precisely, we have zðq; q; the following result. _ tÞ Corollary 2.8. At any instant of time t, at which the arbitrary LC vector zðq; q; _ tÞ, the control given by (2.21) becomes optimal, belongs to the range space of BsT ðq; q; in the sense that it minimizes J(t) at that time. _ tÞ belongs to the range space of BsT ðq; q; _ tÞ at time t, we can Proof. When zðq; q; T _ tÞ becomes, give z Z Bs g for some m-vector g. Hence, by equation (2.24), Q2c ðq; q; Q2c dN K1=2 ðI KBsCBs Þz Z N K1=2 ½I KBsCBs BsT g T T K1=2 T T ½I KBsT ðBC Z N K1=2 ½I KðBC s Bs Þ Bs g Z N s Þ Bs g Z 0:

ð2:37Þ

The second equality follows from the fourth MP-inverse property, and the last equality follows from the ﬁrst MP-inverse property (Graybill 2001). & Corollary 2.9. (i) The two components Q1c and Q2c of the control vector Q given in equation (2.22) are N-orthogonal to one another and (ii) Q1c is (MN )-orthogonal to the null space of the matrix A. c

Proc. R. Soc. A (2008)

Tracking control of nonlinear systems

2349

Proof. (i) Since ½Q2c T NQ1c Z z T ðI KBsCBs ÞT N K1=2 NN K1=2 BC s ðbK Bs a s Þ Z z T ðI KBsCBs ÞBC s ðbK Bs as Þ Z 0:

ð2:38Þ

The result follows by using the fourth and then the second MP-inverse property (Graybill 2001). (ii) Substituting for Bs and using the identity X CZ X T ðXX T ÞC we get K1=2 ðAM K1 N K1=2 ÞC N K1=2 BC s ZN

Z N K1=2 ðAM K1 N K1=2 ÞT ½ðAM K1 N K1=2 ÞðAM K1 N K1=2 ÞT C Z ðMN ÞCAT ðAM K1 N K1 M K1 AT ÞC:

ð2:39Þ

Any vector belonging to the null space of the matrix A satisﬁes the relation AvZ0, whose explicit solution is v Z ðI KACAÞw;

ð2:40Þ

where w is any n-vector. Hence v T MNQ1c Z w T ðI KACAÞT ðMN ÞðMN Þ K1 AT ðAM K1 N K1 M K1 AT ÞCðbK Bs a s Þ Z w T ½I KAT ðAT ÞCAT ðAM K1 N K1 M K1 AT ÞCðbK Bs a s Þ Z 0; since ½I KAT ðAT ÞCAT Z 0.

ð2:41Þ &

Q1c

Corollary 2.10. The component as deﬁned in (2.23) belongs to the range K1=2 C space of the matrix N Bs , and the component Q2c as deﬁned in equation (2.24) is N-orthogonal to the range space of N K1=2 BsT . Proof. The ﬁrst result is obvious from the expression for Q1c in equation (2.23). The second follows because any vector a in the range space of N K1=2 BsT can be expressed as a Z N K1=2 BsT g;

ð2:42Þ

aT NQ2c Z gT Bs N K1=2 NN K1=2 ðI KBC s Bs Þz Z 0:

ð2:43Þ &

for some m-vector g, and so Hence the result.

We note that the range space of N K1=2 BsC is the same as the range space of N BsT . This follows because for any vector g, we can always ﬁnd a vector a such that N K1=2 BsT aZ N K1=2 BsCg, and vice-versa. K1=2

Remark 2.11. The addition of the second component Q2c to the control so that Q c Z Q1c C Q2c , whatever the LC n-vector z may be, still ensures that the description of the trajectory given by equations (2.18) and (2.19) is exactly satisﬁed. However, it contributes an additional amount given by kðI KBsCBs Þzk2 to the cost of the control, since the two components Q1c and Q2c are N-orthogonal Proc. R. Soc. A (2008)

2350

F. E. Udwadia

_ tÞ to each other (corollary 2.9). Furthermore, if at any time t the function zðq; q; belongs to the null space of Bs this additional cost at that time simply becomes kzk2 . It is sometimes advantageous to describe the motion in terms of a system of ﬁrst-order differential equations instead of a system of second-order equations. _ so One obvious way of doing this is to deﬁne a new variable, an n-vector v Zq, that the equation of motion takes the so-called state space form given by q_ Z v; qð0Þ Z q0 ; ð2:44aÞ M ðq; tÞv_ Z Qðq; v; tÞ;

vð0Þ Z v0 :

ð2:44bÞ

The ﬁrst equation of this set is simply a deﬁnition, whereas the second equation of the set is the one that contains the actual dynamics. And while this is the form often used when numerically integrating equation (2.1), the description of the motion of a mechanical system in terms of a system of ﬁrst-order differential equations goes far beyond just its use in numerical procedures. For, often these ﬁrst-order descriptions (i.e. descriptions using ﬁrst-order differential equations) arise when one wants to use descriptions of motion in terms of coordinates that may be more physically meaningful in the context of a particular problem. For example, the ﬁrst-order Hamilton’s equations describing the motion of mechanical systems are often useful when dealing with systems described by a Hamiltonian, and, more generally, the ﬁrst-order set of Poincare´ equations are often useful when dealing with rigid body and multi-body dynamics. We next obtain the entire set of explicit closed-form controllers for exactly tracking the trajectory of mechanical systems described by Poincare´’s equations of motion. (b ) Poincare´ descriptions The Poincare´ equations (Poincare´ 1901) are obtained by deﬁning a new _ where H~ ðq; tÞ is a non-singular matrix whose variable, an n-vector sZH~ ðq; tÞq, elements are known functions of q and t. Denoting the inverse of the matrix H~ by H(q, t), so that q_ Z H ðq; tÞs, the Lagrange equation (2.1) can be rewritten in ﬁrst-order form as (Udwadia & Phohomsiri 2007) q_ Z H ðq; tÞs; qð0Þ Z q0 ;

ð2:45Þ

Mp ðq; tÞ_s Z Sðq; s; tÞ; sð0Þ Z s0 ;

ð2:46Þ

where the state variables are now the n-vectors q and s. In rigid body dynamics, the 3-vector s is often chosen to constitute the three components of the angular velocity of the body. Here the matrix Mp is again positive deﬁnite. Though both ﬁrst-order equations are, mathematically speaking, on a par, the ﬁrst equation of the set (equation (2.45)) is a kinematic relation involving the deﬁnition of the new variable s—hence, an identity as in equation (2.44a)—while the second equation (equation (2.46)) is the dynamical equation of motion. The system is required to track the trajectory described by the consistent relations fpi ðq; tÞ Z 0; jpi ðq; s; tÞ Z 0; Proc. R. Soc. A (2008)

i Z 1; .; h;

ð2:47Þ

i Z h Z 1; .; m:

ð2:48Þ

Tracking control of nonlinear systems

2351

The subscripts and superscripts ‘p’ indicate that we are dealing with the Poincare´ description of the motion of the dynamical system. In order to control the system so that it exactly tracks this trajectory, we apply a control so that the equation of motion of the controlled system now becomes q_ Z H ðq; tÞs; qð0Þ Z q0 ; ð2:49Þ Mp ðq; tÞ_s Z Sðq; s; tÞ C S c ðq; s; tÞ;

sð0Þ Z s0 ;

ð2:50Þ

and the initial conditions q0 and s0 satisfy the trajectory requirements (2.47) and (2.48). The explicit expression for the control force, yielding the entire set of all controllers that will cause the dynamical system described by equations (2.45) and (2.46) to track the required trajectory described by equations (2.47) and (2.48) will now be obtained. As before, we deﬁne the matrix Gp ðq; tÞ d½N K1=2 ðq; tÞMp ðq; tÞ K1 Z MpK1 ðq; tÞN K1=2 ðq; tÞ; ð2:51Þ where the matrix N(q, t) is any positive-deﬁnite matrix, and the scaled variables 9 _ s_s Z GpK1s; > > > > = K1 1=2 K1 1=2 a s Z Gp a Z ðN Mp ÞðMp SÞ Z N S ð2:52Þ > and > > > ; s_cs Z N 1=2 S c ; so that relation (2.50) on pre-multiplication with N 1/2 becomes s_s Z a s Cs_sc ; ð2:53Þ which we note is of the same form as equation (2.6), except that we now have s instead of q, and ﬁrst derivatives instead of second derivatives with respect to time. On differentiating equation (2.47) twice with respect to time and differentiating equation (2.48) once with respect to time, and using relation (2.49), we then obtain the matrix equation Ap ðq; s; tÞ_s Z bp ðq; s; tÞ;

ð2:54Þ

which, upon using the ﬁrst equation of the set (2.52), can be rewritten as Bs ðq; s; tÞ_ss Z bp ðq; s; tÞ;

ð2:55Þ

where we deﬁne the m by n matrix Bs of rank k by the relation Bs ðq; s; tÞ Z Ap ðq; s; tÞGp ðq; tÞ:

ð2:56Þ

Again we note the similarity between equations (2.10) and (2.11) and equations (2.54) and (2.55). We next decompose the n-vector s_s as s_s Z BsCBss_s C ðI KBC s Bs Þs_s ;

ð2:57Þ

in a manner similar to equation (2.12), with the matrix Bs now deﬁned by equation (2.56). Proceeding along with similar lines as before (equations (2.13)–(2.16)), we then ﬁnd that K1=2 ðI KBsCBs Þz p dS c1 C S c2 ; ð2:58Þ Sc Z N K1=2s_sc Z N K1=2 BC s ðbp K Bs a s Þ C N

Proc. R. Soc. A (2008)

2352

F. E. Udwadia

where z p ðq; s; tÞ is any arbitrary n -vector. To ensure a unique solution of equations (2.49) and (2.50), we may then take the components of z p ðq; s; tÞ to be C 1 functions (or, more generally, LC functions) of q, s and t. Following the same lines as in §2a, for any Poincare´ system described by equation (2.45) that is required to exactly track the trajectory described by relations (2.47) and (2.48), we now obtain results that are analogous to those given in results 2.1 and 2.5. Using as deﬁned in relation (2.52), and Bs in (2.56), _ s, Q/ S, b/ bp , M / Mp , we simply make the following variable changes: q/ z / z p , A/ Ap , and G / Gp , in the expressions given in equations (2.21) and (2.27) to obtain the corresponding control forces. Remark 2.12. We can prove corollaries similar to corollaries 2.2, 2.8–2.10. The same goes for remarks 2.3, 2.4, 2.6, 2.7, 2.11.

3. The close connection between nonlinear control and analytical dynamics Since almost all mechanical systems are nonlinear in their behaviour, including even simple ones like a pendulum, we will be mainly addressing nonlinear systems here. The problem of control can be placed within the context of analytical mechanics by reinterpreting constrained motion in mechanical systems. Consider a mechanical system described by equation (2.17). When the system is further subjected to the trajectory requirements (2.18) and (2.19)—i.e. subjected to further constraints—additional control (constraint) forces are brought into play by Nature so that the controlled (constrained) system moves in such a manner that it satisﬁes these trajectory requirements (constraints). Thus, the additional _ tÞ that Nature provides may be thought of as the control (constraint) Q c ðq; q; control it generates in order for the system to satisfy the trajectory requirements given by equations (2.18) and (2.19). One might imagine Nature as a control engineer, attempting to control the mechanical system so that it satisﬁes the given trajectory requirements described by equations (2.18) and (2.19). However, to entertain such an interpretation, one is led to ask the following three questions. (i) To what extent might Nature be perceived as acting like a control engineer? (ii) Does Nature appear to be performing the control in any kind of an optimal manner? (iii) If so, what is the cost function (or functional) it appears to minimize? We shall start answering these questions by ﬁrst asserting that Nature appears to choose the weighting matrix N ðq; tÞZ ½M ðq; tÞ K1 dNNature ðq; tÞ in equation (2.21). Here M is the so-called ‘mass matrix’ and appears on the l.h.s. of equation (2.17). Result 3.1. Nature seems to control a mechanical system described by equation (2.17), so it exactly satisﬁes the trajectory described by the requirements (2.18) and (2.19) by choosing the weighting matrix to be N ðq; tÞ Z ½M ðq; tÞ K1 : Proc. R. Soc. A (2008)

ð3:1Þ

Tracking control of nonlinear systems

2353

Proof. Setting N ðq; tÞZ ½M ðq; tÞ K1 equation (2.21) now yields _ tÞ Z M 1=2 BCðbKAaÞ C M 1=2 ðI KBCBÞz; Q c ðq; q;

ð3:2Þ

where B dAG Z AM K1=2 and Bs a s Z AGG K1 aZ Aa. But equation (3.2) is exactly the equation of motion of a general constrained mechanical system, as given by Udwadia (2000) and Udwadia & Kalaba (2002). Hence the result. & Result 3.2. If we assume that Nature observes d’Alembert’s principle, then it appears to be minimizing the cost _ tÞT ½M ðq; tÞ K1 Q c ðq; q; _ tÞ; JNature ðt Þ Z ½Q c ðq; q; ð3:3Þ at each instant of time while controlling the system deﬁned by equation (2.17) so it exactly satisﬁes the trajectory described by the requirements (2.18) and (2.19). Proof. Setting N ðq; tÞZ ½M ðq; tÞ K1 in result 2.5, equation (2.27), which gives the optimal control while minimizing J(t), yields _ tÞ Z M 1=2 BCðbKAaÞ Z Q1c : Q c ðq; q; ð3:4Þ c Using this expression for Q in equation (2.20), we ﬁnd that we obtain the correct equation of motion of a constrained mechanical system that obeys d’Alembert’s principle, as given by Udwadia & Kalaba (1992, 1996). & Remark 3.3. Result 3.2 connects directly with the basic principles of analytical dynamics. In fact, d’Alembert’s principle, which is a principle that leads to a mathematical description of motion, which is in close conformity with observations, and which is one of the pivotal assumptions of analytical dynamics, is equivalent to Gauss’s principle (Gauss 1829; Udwadia & Kalaba 1996). And Gauss’s principle states that: of the entire set of constraint (control) forces that cause a constrained (controlled) mechanical system (2.17) to exactly satisfy the constraints (trajectory) described by requirements (2.18) and (2.19), Nature seems to choose that constraint force that minimizes JNature (t) at each instant of time. Thus we ﬁnd that (i) Nature seems to choose the weighting matrix N ðq; tÞZ ½M ðq; tÞ K1 and (ii) if we assume that d’Alembert’s principle is true, then Nature seems to pick the one controller given in result 3.2 that minimizes the cost JNature (t). Nature appears to go well beyond what most modern control engineers would try to do, by minimizing the cost JNature given in relation (3.3) at each instant of time, rather than minimizing the integral of this cost over any given span of time, as is the common practice in the ﬁeld of controls. Remark 3.4. When d’Alembert’s principle is assumed to be satisﬁed, Nature appears to be performing acceleration feedback control, where the control force is given by Q c Z Q1c Z M 1=2 BCðbKAaÞ ZKKNature ðAaKbÞ: ð3:5Þ _ tÞ dM 1=2 ðq; tÞBCðq; q; _ tÞ and the It is made up of the gain matrix KNature ðq; q; feedback error e(t)d(AaKb). The quantity e(t) is simply the extent to which _ tÞ does not satisfy the the acceleration of the uncontrolled system aðq; q; trajectory requirements imposed on it by relation (2.10). Hence, Nature appears to be behaving like a control engineer, using feedback control. The minus sign in equation (3.5) is taken so as to be compatible with the control concept of negative feedback. While most modern control engineers usually use integral, Proc. R. Soc. A (2008)

2354

F. E. Udwadia

velocity and proportional feedback in mechanical systems, few use acceleration feedback as Nature appears to be using. This acceleration feedback does not _ Þ, require the measurement of acceleration because a is a function of qðt Þ and qðt _ tÞZM K1 ðq; tÞQðq; q; _ tÞ. and can be obtained from their measurement, since aðq; q; Also, the gain matrix KNature used by Nature is complex and its elements are, in general, highly nonlinear functions of q, q_ and t. Remark 3.5. The component Q1c of the constraint force used by Nature is always orthogonal to the null space of the matrix A given in equation (2.10). By corollary 2.9(ii), Q1c is always MN-orthogonal to the null space of A. Since Nature picks N ðq; tÞZ ½M ðq; tÞ K1 , the result follows. The null space of A in analytical dynamics is called the space of virtual displacements and d’Alembert’s principle simply posits the assumption that the control (constraint) force n-vector Q1c is orthogonal to the null space of A, or in more analytical dynamics terms: the ‘work done’ by the constraint force Q c under virtual displacements at each instant of time is zero. Remark 3.6. Many mechanical systems, however, do not satisfy d’Alembert’s principle (Goldstein 1976). In such systems, the control (constraint force) does work under virtual displacements. To obtain the equation of motion for such systems one requires additional information about the work done by the control forces under virtual displacements. One then needs to prescribe, for a speciﬁc _ tÞ such that mechanical system, the C 1 vector C ðq; q; _ tÞ Z w T ðtÞC ðq; q; _ tÞ; w T ðtÞQ c ðq; q;

ð3:6Þ

where w(t) is any virtual displacement, i.e. any n-vector in the null space of A. In that case, the equation of motion of the constrained mechanical system is known to be described by the relation (Udwadia & Kalaba 2002) _ tÞ C M 1=2 BCðbKAaÞ C M 1=2 ðI KBCBÞM K1=2 C: Mq€ Z Qðq; q;

ð3:7Þ

This same result also follows directly from equations (2.20) and (2.21) by setting N ðq; tÞZ ½M ðq; tÞ K1 and zðt ÞZ M K1=2 C . Thus, if a mechanical system is non_ tÞ is ideal, the non-idealness being described by relation (3.6) in which C ðq; q; speciﬁed at each instant of time, then Nature chooses the n-vector z(t) in relation (2.21) to be zðt ÞZ M K1=2 C , so that condition (3.6) is satisﬁed along with relations (2.18) and (2.19) at each instant of time for the speciﬁc system at hand!

4. General systems described by ﬁrst-order differential equations In this section, we proceed to general dynamical systems that are described by n ﬁrst-order, non-autonomous, differential equations given by Mg ðx; tÞx_ Z f ðx; tÞ; xð0Þ Z x 0 ; ð4:1Þ where we shall again take Mg to be a positive-deﬁnite n by n matrix. We require that this dynamical system track a trajectory described by the equations fi ðx; tÞ Z 0; i Z 1; .; m; ð4:2Þ where we assume that the trajectory described is feasible and the system of m equations is consistent. In order to track this trajectory, we apply a control f c so Proc. R. Soc. A (2008)

Tracking control of nonlinear systems

2355

that the equation describing the time evolution of the controlled dynamical system becomes ð4:3Þ Mg ðx; tÞx_ Z f ðx; tÞ C f c ðx; tÞ; xð0Þ Z x 0 ; where we assume that the initial conditions satisfy the trajectory requirements (4.2). As before, we differentiate these m equations (4.2) with respect to time to give the relation Agx_ Z bg ðx; tÞ; ð4:4Þ where Ag(x, t) is an m by n matrix of rank k, and set Gg Z ½N 1=2 ðx; tÞMg ðx; tÞ K1 : Pre-multiplying equation (4.3) by N 1/2 we get x_s Z a s Cx_cs ; where _ x_s Z GgK1x;

a s Z GgK1 a Z ðN 1=2 Mg ÞðMgK1 f Þ Z N 1=2 f

ð4:5Þ ð4:6Þ and x_cs Z N 1=2 f c :

ð4:7Þ Here, aZ MgK1 f : We note that equations (4.6) and (4.7) have the same form as equations (2.53) and (2.52), respectively. Furthermore, equation (4.4) can be expressed as Bs ðx; tÞx_ s Z bg ðx; tÞ; ð4:8Þ where Bs ðx; tÞ Z Ag ðx; tÞGg ðx; tÞ: ð4:9Þ Expressing x_s in terms of orthogonal components, we get x_s Z BsCBsx_ s C ðI KBsCBs Þx_ s ; ð4:10Þ c and replacing x_ s on the r.h.s. of (4.10) by as Cx_ s yields ð4:11Þ BsCB sx_ s Z BsCðbg K Bs a s Þ; from which we get the following result by following the same lines as in equations (2.13)–(2.16). Result 4.1. Consider the nonlinear, non-autonomous dynamical system (4.1). The system is required to track the trajectory described by equation (4.2). Assuming that the initial conditions satisfy the trajectory described by equation (4.2), all the possible LC continuous controls that exactly track this trajectory are explicitly given by K1=2 f c Z N K1=2 BC ðI KBC ð4:12Þ s ðbg K Bs a s Þ C N s Bs Þz g ; where zg(x, t) is any arbitrary n-vector whose components are continuously differentiable (or more generally are LC) functions of their arguments; N(x, t) is any arbitrary n by n positive-deﬁnite matrix, a s Z N 1=2 ðx; tÞf ðx; tÞ; Bs ðx; tÞZ Ag ðx; tÞ½N 1=2 ðx; tÞMg ðx; tÞ K1 is an m by n matrix; Ag is the m by n matrix of rank k, and bg(x, t) is the m-vector deﬁned in equation (4.4).

Result 4.2. Consider the nonlinear dynamical system described by equation (4.1) that needs to be controlled through the addition of a control f c ðx; tÞ, so that the trajectory described by equation (4.2) is exactly tracked. Assuming that the system satisﬁes the trajectory requirements initially, the optimal controller that Proc. R. Soc. A (2008)

2356

F. E. Udwadia

causes the system to (i) exactly track the required trajectory and (ii) minimize at each instant of time t, the cost J ðt Þ Z ½ f c ðx; tÞT N ðx; tÞf ðx; tÞ; ð4:13Þ for a given n by n positive-deﬁnite matrix N, is explicitly provided by ð4:14Þ f c ðx; tÞ Z N K1=2 BC s ðbg K Bs a s Þ; where Bs, bg and a s are deﬁned in relations (4.9), (4.8) and (4.7), respectively. Proof. Similar to, and along with the same lines as, the proof of result 2.5.

&

Denoting f c Z f1c C f2c ;

ð4:15Þ

where f1c dN K1=2 BsCðbg K Bs a s Þ Z N K1=2 ½Ag MgK1 N K1=2 Cðbg K Ag aÞ

ð4:16Þ

and f2c dN K1=2 ðI KBsCBs Þz g ; we have the following remarks.

ð4:17Þ

Remark 4.3. We can prove corollaries similar to corollaries 2.2, 2.8–2.10. Four important results that emerge from this are (i) the two components f1c and f2c of the control vector f c are N-orthogonal to one another, (ii) f1c is (MN )-orthogonal to the null space of the matrix A, (iii) the minimum cost is J ðt ÞZ kBsCðbg K Ag aÞk2 and it occurs at those times when f2c Z0, and (iv) the control cost contributed by f2c is given by kðI KBsCBs Þz g k2 . Remark 4.4. So far, it has been assumed that the equations that describe the trajectory are consistent. This may not happen in practical situations; errors due to numerical computations, for example, could make these equations inconsistent. Hence, instead of equation (4.8) we would have the equation Bs x_s Z bg C 3ðx; tÞ;

ð4:18Þ

where 3ðx; tÞ is the error caused by the inconsistency of the trajectory equation (4.2) (as also, analogously, for the equation sets (2.18)–(2.19), and (2.47)–(2.48)). We would then need to replace equation (4.11) with BsCBsx_cs Z BC ð4:19Þ s ðbK Bs a s Þ C dðx; tÞ; C where dZ Bs 3 is the error. The least-squares solution to this inconsistent equation remains x_ cs Z B sCðbK Bs a s Þ C ðI KBsCBÞz; ð4:20Þ as before, pointing out that results 4.1 and 4.2 (and similarly results 2.1 and 2.5) provide the control in this least-squares sense, even when the trajectory description is inconsistent. 5. Example In this section, we present an application of the results obtained by considering the trajectory tracking of a chaotic Rossler system (Rossler 1976) by a Lorenz system (Strogatz 1994). We begin by considering the scaled Lorenz oscillator as a Proc. R. Soc. A (2008)

2357

Tracking control of nonlinear systems

system that we would like to control. Its description is given by the three ﬁrstorder differential equations, 2 3 2 3 x_ 1 sðx 2 K x 1 Þ 6 7 6 7 x_ d 4 x_ 2 5 Z d4 rx 1 K x 2 K x 1 x 3 5 Z f ðxÞ; x 1 ð0Þ Z x 3 ð0Þ Z 5; x 2 ð0Þ Z 10; x 1 x 2 Kbx 3

x_ 3

ð5:1Þ where we take sZ10, rZ28, bZ8/3 and dZ10. This system exhibits chaotic motion for the chosen values of the parameters (Strogatz 1994). We want the motion of the ﬁrst two components, x 1(t) and x 2(t), of this chaotic Lorenz system to track the ﬁrst two components, y1(t) and y2(t), respectively, of a very different chaotic system—a Rossler system—which is described by the equations (Rossler 1976), 2 3 2 3 Kðy2 C y3 Þ y_1 6 7 6 7 y_ d 4 y_2 5 Z 4 y1 C ay2 5 Z hðyÞ; y1 ð0Þ Z 3; y2 ð0Þ Z 12; y 3 ð0Þ Z 6; y_3

d C y3 ðy1 KcÞ

ð5:2Þ with aZ0.1, cZ18 and dZ0.3. For these parameter values, the Rossler system is also known to be chaotic (Strogatz 1994). Thus, the aim is to ﬁnd the control inputs needed to be applied to one chaotic system (the Lorenz system here) so that it tracks two of the components of the motion of another different chaotic system (the Rossler system). The control input 3-vector f c that needs to be applied to the Lorenz system is required to be found so that the cost J ðt ÞZ ½f c T NL f c is minimized at each instant of time. The weighting matrix NL is taken to be the diagonal matrix, NL Z Diag½x 3 ðt Þ2 C 1; x 2 ðt Þ2 C 1; x 1 ðt Þ2 C 1. This weighting matrix is guaranteed to be positive deﬁnite. A simple way to formulate this nonlinear trajectory tracking problem is to consider the augmented dynamical system, # " # " f ðxÞ x_ : ð5:3Þ Z y_ hðyÞ Our task would be to ﬁnd a control input 6-vector F c Z ½ðf c ÞT ; ðh c ÞT T which causes the system # # " c " # " f ðx; yÞ f ðxÞ x_ ð5:4Þ C c Z y_ h ðx; yÞ hðyÞ to satisfy the following trajectory tracking requirements: ðiÞ f1 ðx; yÞ dx 1 ðt ÞK y1 ðt Þ Z 0;

f2 ðx; yÞ dx 2 ðt ÞK y2 ðt Þ Z 0

and

ð5:5Þ

ðiiÞ y_ Z hðyÞ; implying thereby that the control force component h c ðxðtÞ; yðt ÞÞ h 0:

ð5:6Þ

Differentiating f1 ðx; yÞ and f2 ðx; yÞ with respect to time we get the relations f_ 1 ðx; yÞ dx_ 1 ðtÞKy_ 1 ðtÞZ 0 and f_ 2 ðx; yÞ dx_ 2 ðtÞKy_ 2 ðtÞZ 0. Since our theory requires that the given initial conditions satisfy the above-mentioned two Proc. R. Soc. A (2008)

2358 first component of Lorenz and Rossler system

F. E. Udwadia 30 20 10 0 −10 −20 −30

0

10

20

30

40

50

time Figure 1. The dynamics of the ﬁrst component of the Lorenz system x 1(t) and the Rossler system y1(t) are shown by the solid line and the thick dashed line, respectively. Fifty seconds of response are shown.

trajectory requirements, and they do not, we shall modify the trajectory requirements to f_ 1 ZKaf1 ;

f_ 2 ZKaf2 ;

ð5:7Þ

where aO0 is chosen to be a suitable parameter. We note that asymptotic solutions of equation (5.7) as t/N are fi (t)Z0, iZ1, 2 as required by (5.5). The parameter a in the numerical example is chosen to be 0.5. The trajectory description given by (5.7) and (5.6) then leads to (see equation (4.4)) 3 2 1 0 0 K1 0 0 7 6 3 2 Kaðx 1 K y1 Þ 6 0 1 0 0 K1 0 7 7 6 7 7 6 6 ð5:8Þ Ag Z 6 0 0 0 1 0 0 7 and bg Z 4 Kaðx 2 K y2 Þ 5: 7 6 60 0 0 0 1 07 hðyÞ 5 4 0 0 0 0 0 1 The weighting matrix N for the augmented six-dimensional dynamical system given in equation (5.4) will also need to be augmented so it is a diagonal 6 by 6 matrix. Owing to the trajectory requirement hc(t)h0, we can choose its last three diagonal entries to be each equal to unity; the ﬁrst three diagonal entries remain the same as those of the matrix NL. The explicit control input that causes the trajectory described by (5.6) and (5.7) to be tracked is then given by equation (4.14) with Bs Z Ag N K1=2 , since the 6 by 6 matrix Mg is an identity matrix. All the computations are performed using MATLAB and the integration is carried out using a relative error tolerance of 10K10 and an absolute error tolerance of 10K13. Figure 1 shows the dynamics of the ﬁrst component of the two separate dynamical systems described by equations (5.1) and (5.2). Proc. R. Soc. A (2008)

2359

Tracking control of nonlinear systems (a) 2.0

(b) tracking error (×10 –12)

1.5 tracking error

1.0 0.5 0 − 0.5 −1.0 −1.5 −2.0

0

10

20

30

40

50

6 4 2 0 −2 −4 −6 80

85

90 time

time

95

100

(a) 3 2 1 0 −1 −2 0

10

20

30 time

40

50

control inputs to Lorenz system (× 10 4)

control inputs to Rossler system (×10 −11)

Figure 2. (a) The solid line shows e1(t) over 50 s of integration. The dashed line shows e2(t). (b) Tracking errors e1(t) and e2(t) over the time interval [80–100] seconds using the same line conventions as in (a).

(b) 3 2 1 0 −1 −2 −3 −4 −5

0

10

20

30

40

50

time

Figure 3. (a) Components of the control inputs acting on the Rossler system. Note that the vertical scale ranges from approximately K3!10K11 to 3!10K11, making hcz0. The ﬁrst component of the control input is shown by a solid line, the second by a dashed line and the third by a thick solid line. (b) Components of the control inputs acting on the Lorenz system, using the same line notation.

On application of the optimal control input, the tracking errors e 1(t)d x 1(t)Ky1(t) and e 2(t)dx 1(t)Ky1(t) are shown in ﬁgure 2. These errors, as time increases, go down to the same order of magnitude as the tolerance used in the numerical integration of the differential equations as shown in ﬁgure 2b. Figure 3a shows the control inputs acting on the Rossler system. They are theoretically supposed to be zero, as required by the second trajectory requirement (5.6). They are seen to be very small, and of the same order of magnitude as the tolerances with which the integration is carried out. Figure 3b shows the control inputs required to be given to the Lorenz system so that it tracks the ﬁrst two components of the Rossler system. The third component of the control input to the Lorenz system is seen to be zero, since we require only the ﬁrst two components to be tracked. The optimal control cost J(t) is shown in ﬁgure 4. Proc. R. Soc. A (2008)

2360

F. E. Udwadia 10

8

6

4

2

0

10

20

30

40

50

time Figure 4. The optimal control cost JðtÞZ ½F c T N ðxÞF c as a function of time.

30 20 10 0 −10 −20 −30 −30

−20

−10

0

10

20

30

Figure 5. Projection of phase space trajectories over the time interval [0 –100] seconds of the controlled Lorenz system (dashed line) on the (x 1, x 2) plane and those of the Rossler system (solid line) on the (y1, y2) plane.

The projection of the phase trajectory of the controlled Lorenz system on to the (x 1, x 2) plane (solid line) and the projection of the corresponding phase trajectory of the Rossler system on to the (y1, y 2) plane (dashed line) are shown in ﬁgure 5. The ﬁgure shows the manner of convergence of these projected trajectories, which start with different initial conditions. Because the third component of the Rossler system is not tracked, the three-dimensional phase portrait looks very different for the controlled Lorenz system. The phase space portraits of these two systems are shown in ﬁgure 6. Proc. R. Soc. A (2008)

2361

Tracking control of nonlinear systems

(a) 50

(b) 50 40

0 −50

3

−100

30 20 10 0

20

20

0 −20

−20

0

20

0 −20

−20

0

20

Figure 6. The phase portraits of (a) the controlled Lorenz system and (b) the chaotic Rossler system over the time interval [0–100] seconds. The projections of these phase portraits on the horizontal plane are shown in ﬁgure 5.

6. Conclusions and remarks The methodology for the tracking control of nonlinear systems proposed herein has been inspired by results in analytical dynamics. This paper begins by developing this methodology for systems described by second-order differential equations, as commonly found in the Lagrangian and Newtonian mechanics, as well as ﬁrst-order differential equations, as found in the Hamiltonian and Poincare´ formulations of mechanics. It then extends the methodology to full state control of general nonlinear dynamical systems. The main contributions of the paper are the following: (i) The development of an explicit closed-form expression that provides the entire set of continuous tracking controllers that can exactly track a given trajectory description, assuming that the system’s initial conditions satisfy the description of the trajectory. We obtain explicit closed-form expressions for the controllers, which can be computed in real time. (ii) The development of a simple formula that explicitly gives the tracking control that minimizes the control cost J(t) at each instant of time. An explicit expression for the minimal cost is also obtained. (iii) For a general, ﬁrst-order, nonlinear system Mg ðx; tÞx_ Z f ðx; tÞ, x(0)Zx 0, the entire set of controllers needed to satisfy the trajectory described by the consistent equations fi(x, t)Z0, iZ1, ., m, is explicitly given by ð6:1Þ f c df1c C f 2c Z N K1=2 BsCðbg K Bs a s Þ C N K1=2 ðI KBsCB s Þz g ; where z g(x, t) is any LC function; N(x, t) is any positive-deﬁnite weighting matrix; and a s, bg and Bs are as deﬁned in equations (4.7)–(4.9). (iv) For a given weighting matrix N, the total control input can be split into two parts: a part that solves the optimal control problem that minimizes J(t) while exactly tracking the trajectory, and a second additive part that is N-orthogonal to the ﬁrst. While the addition of the second part to the optimal control effort allows the trajectory requirements to be still exactly satisﬁed, the norm of the control cost increases, in general, when it is added. The minimum cost J ðtÞZ kBsCðbg K Ag aÞk2 , and the corresponding controller that yields this minimum cost is f c df1c Z N K1=2 BsCðbg K Bs a s Þ. Proc. R. Soc. A (2008)

2362

F. E. Udwadia

The addition of the second part f2c Z N K1=2 ðI KBsCBs Þzg , so that f cd f1c C f2c , increases the cost beyond the optimal by kðI KBsCBs Þzg k2 . At each instant of time t, when z g(x, t) belongs to the range space of BsT , f2c ðt ÞZ 0, so that f c ðt Þ df1c ðtÞ, hence making the control at that instant of time optimal. (v) The close connection between nonlinear control and analytical mechanics is pointed out. Here we see that Nature seems to control a mechanical system so that it satisﬁes a given set of trajectory requirements—or alternatively stated, tracks a given trajectory—by using feedback control, much like a control engineer, except that instead of proportional, integral or derivative feedback control, which is commonly used by the control engineer, it uses acceleration feedback. While considerable work has been done on PID control, acceleration feedback seems to be far less studied by modern-day control engineers. Following Nature’s cue, the results developed in this paper point perhaps towards the need for more work in the area of acceleration feedback in mechanical systems. Furthermore, Nature seems to minimize the control cost J ðtÞZ ½Q c T N Q c at each instant of time, and it appears to use M K1 ðq; tÞ for the weighting matrix N. Nature’s choice of this weighting matrix can be understood when thought of in terms of a multibody mechanical system that is required to satisfy a given trajectory description (a set of constraints), and so track a given trajectory. Since it takes a larger control effort to move a body belonging to the multi-body system that has a larger inertia, Nature, in its effort to make the entire system satisfy the given trajectory description, appears to prefer applying control forces to bodies with the smaller inertias. Again, following Nature’s cue, the use of M K1 ðq; tÞ as a weighting matrix for deﬁning the control cost may be useful in many other dynamical systems, especially mechanical ones. (vi) While we have demonstrated the methodology by illustrating its use in determining the control required to be applied to a chaotic Lorenz system so that it tracks some components of the motion of another chaotic Rossler system, the general methodology can be used for more complex nonlinear mechanical systems dealing with, for example, orbital mechanics (Lam 2006). (vii) Finally, we note that the methodology presented herein does not include any magnitude constraints on the control, and work on this topic is currently being pursued. Furthermore, the results provided herein deal with full state control and further developments along the lines pursued herein to underactuated robotic systems would be useful. The effects of model errors and disturbance control also need to be investigated.

References Boccaletti, S., Kurths, J., Osipov, G., Valladares, D. & Zhou, C. 2002 The synchronization of chaotic oscillators. Phys. Rep. 366, 1–101. (doi:10.1016/S0370-1573(02)00137-0) Brockett, R. W. 1977 Control theory and analytical mechanics. In Geometric control theory (eds C. Martin & R. Herman), pp. 1–66. Brookline, MA: Math. Sci. Press. Brogliato, B., Lozano, R. & Maschke, B. 2007 Dissipative systems analysis and control. Berlin, Germany: Springer. Proc. R. Soc. A (2008)

Tracking control of nonlinear systems

2363

Chaou, Y. & Chang, W. 2004 Trajectory control. Lecture Notes on Control and Information Sciences. Berlin, Germany: Springer. Chen, H. K. 2002 Chaos and synchronization of a symmetric gryo with linear plus cubic damping. J. Sound Vib. 255, 719 –740. (doi:10.1006/jsvi.2001.4186) Gauss, C. F. 1829 Uber ein neues algemeines Grundgesetz der Mechanik. Zeitschrift fur die Reine und Angewandte Mathematik 4, 232–235. Goldstein, H. 1976 Classical mechanics. New York, NY: Addison-Wesley. Graybill, F. 2001 Matrices with application to statistics Duxbury classics series. Paciﬁc Grove, CA: Duxbury Publishing Company. Hamel, G. 1949 Theoretische mechanik. Berlin, Germany: Springer. Lagrange, J. L. 1811 Mechanique analytique. Paris, France: Mme Ve Coureier. Lam, T. 2006 New approach to mission design based on the fundamental equation of motion. J. Aerosp. Eng. 19, 59– 67. (doi:10.1061/(ASCE)0893-1321(2006)19:2(59)) Lei, Y., Xu, W. & Zheng, H. 2005 Synchronization of two chaotic gyros using active control. Phys. Lett. A 343, 153 –158. (doi:10.1016/j.physleta.2005.06.020) Moore, E. H. 1910 On the reciprocal of the general algebraic matrix. Abstr. Bull. Am. Soc. 26, 394– 395. Naidu, D. 2003 Optimal control systems. London, UK: CRC Press. Penrose, R. 1955 A generalized inverse of matrices. Proc. Camb. Philos. Soc. 51, 406 – 413. Poincare´, H. 1901 Sur une forme nouvelle des equations de la mechanique. CR Acad. Sci. Paris 132, 369 – 371. Pyras, K. 1996 Weak and strong synchronization of chaos. Phys. Rev. E 54, 4508–4511. (doi:10. 1103/PhysRevE.54.R4508) Rosenblum, M., Pilkovsky, A. & Kruths, J. 1997 From phase to lag synchronization in coupled chaotic oscillators. Phys. Rev. Lett. 78, 4193–4196. (doi:10.1103/PhysRevLett.78.4193) Rossler, O. 1976 An equation for continuous chaos. Phys. Lett. A 57, 397–398. (doi:10.1016/03759601(76)90101-8) Rulkov, N., Sushchik, M., Tsimring, L. & Abarbanel, H. 1995 Generalized synchronization of chaos in directionally coupled chaotic systems. Phys. Rev. E 51, 980–994. (doi:10.1103/PhysRevE.51.980) Sadegh, N. 1990 Stability analysis of a class of adaptive controllers for robotic manipulators. Int. J. Robot. 9, 74–92. (doi:10.1177/027836499000900305) Sastry, S. 1999 Nonlinear systems analysis, stability, and control. Berlin, Germany: Springer. Slotine, J. & Li, W. 1991 Applied nonlinear control. New York, NY: Englewood Cliffs. Strogatz, S. 1994 Nonlinear dynamics and chaos. Boulder, CO: Westview Press. Udwadia, F. E. 2000 Fundamental principles of lagrangian dynamics: mechanical systems with non-ideal, holonomic, and non-holonomic constraints. J. Math. Anal. Appl. 252, 341– 355. (doi:10.1006/jmaa.2000.7050) Udwadia, F. E. & Kalaba, R. E. 1992 A new perspective on constrained motion. Proc. R. Soc. A 439, 407– 410. (doi:10.1098/rspa.1992.0158) Udwadia, F. E. & Kalaba, R. E. 1996 Analytical dynamics: a new approach. Cambridge, UK: Cambridge University Press. Udwadia, F. E. & Kalaba, R. E. 2002 On the foundations of analytical dynamics. Int. J. Nonlin. Mech. 37, 1079 –1090. (doi:10.1016/S0020-7462(01)00033-6) Udwadia, F. E. & Phohomsiri, P. 2007 Explicit Poincare equations of motion for general constrained systems. Part I. Analytical results. Proc. R. Soc. A 1825, 1–14. (doi:10.1098/rspa. 2007.1825) Zaks, A., Park, E.-H., Rosenblum, M. & Kruths, J. 1999 Alternating locking ratios in imperfect phase synchronization. Phys. Rev. Lett. 82, 4228–4231. (doi:10.1103/PhysRevLett.82.4228)

Proc. R. Soc. A (2008)

Optimal tracking control of nonlinear dynamical systems B Y F IRDAUS E. U DWADIA * ,1,2,3,4,5 1

Department of Aerospace and Mechanical Engineering, 2Department of Civil and Environmental Engineering, 3Department of Mathematics, 4 Systems Architecture Engineering, and 5Information and Operations Management, University of Southern California, 430K Olin Hall, Los Angeles, CA 90089-1453, USA

This paper presents a simple methodology for obtaining the entire set of continuous controllers that cause a nonlinear dynamical system to exactly track a given trajectory. The trajectory is provided as a set of algebraic and/or differential equations that may or may not be explicitly dependent on time. Closed-form results are also provided for the realtime optimal control of such systems when the control cost to be minimized is any given weighted norm of the control, and the minimization is done not just of the integral of this norm over a span of time but also at each instant of time. The method provided is inspired by results from analytical dynamics and the close connection between nonlinear control and analytical dynamics is explored. The paper progressively moves from mechanical systems that are described by the second-order differential equations of Newton and/or Lagrange to the ﬁrst-order equations of Poincare´, and then on to general ﬁrst-order nonlinear dynamical systems. A numerical example illustrates the methodology. Keywords: exact control; tracking control; nonlinear systems

1. Introduction The development of controllers for nonlinear mechanical systems has been an area of intense research over the last two decades or so. Many controllers that have been developed for trajectory tracking of complex nonlinear and multi-body systems rely on some approximations and/or linearizations (Slotine & Li 1991; Sastry 1999; Naidu 2003). Most control designs restrict controllers for nonlinear systems to be afﬁne in the control inputs (Brogliato et al. 2007). Often, the system equations are linearized about the system’s nominal trajectory and then the linearized equations are used along with various results from the well-developed theories of linear control. While this often works well in many situations, there are some situations in which better controllers may be needed. This is especially so when highly accurate trajectory tracking is required to be done in real time on systems that are highly nonlinear. Examples include the exact trajectory control of orbital, attitudinal and elastic motions of a multi-body spacecraft system that is required to perform precision tumbling. *Address for correspondence: Department of Aerospace and Mechanical Engineering, University of Southern California, 430K Olin Hall, Los Angeles, CA 90089-1453, USA ([email protected]). Received 28 January 2008 Accepted 28 March 2008

2341

This journal is q 2008 The Royal Society

2342

F. E. Udwadia

To place this paper within the context of the enormous literature that has been generated in the area of tracking control of nonlinear systems and to highlight what is new in it, we provide a brief review of the methods that have been developed so far and that have been applied to numerous areas of application ranging from chemical process control to robotics. Brogliato et al. (2007) provide an exhaustive review of the methods that have been developed to date for the tracking control of systems along with over 500 references. They point out that methods developed to date rely heavily on PID-type control and, most often, a linear feedback is provided to track a given trajectory. Chaou & Chang (2004) also provide recent developments on trajectory tracking and deal with the same basic theme (linear feedback) along with numerous applications. The optimality criterion considered in the literature to date is the minimization of the control cost integrated over a suitable span of time. In the robotics literature (Brogliato et al. 2007), trajectory tracking using inverse dynamics and model reference control has been used for some time now, and the methods developed therein can be seen as particular subclasses of the formulation discussed in the present work. Nonlinear control methods using controlled Lagrangian and Hamiltonians have also been explored along with passivity theory (Brockett 1977). These methods limit the structure of the control for a mechanical system to be a nonlinear function of its generalized displacement and they usually do not address the issue of control optimality. No such assumptions are made in this paper. Trajectory tracking in the adaptive control context (which is not the subject of this paper) has also been explored together with speciﬁc parametrizations to guarantee linearity in the unknown parameters of a system (Sadegh 1990). Thus, the methods used to date primarily rely on linearizations and/or PID-type control, and they posit assumptions on the structure of the control effort. By contrast, this paper takes a widely different approach that is based on recent results from analytical dynamics. Here the complete nonlinear problem is addressed with no assumptions on the type of controller that is to be used, except that it be continuous. In particular, we do not posit that the control is afﬁne in the inputs, we do not use linearized equations of motion about the desired trajectory, nor do we assume any ‘feedback structure’ to the nonlinear control effort. Furthermore, the optimality criterion used is the minimization of the control cost at each instant of time. As far as is known, the results provided here yield new and explicit methods for the control of highly nonlinear systems. Moreover, we provide for the ﬁrst time the entire class of nonlinear Lipschitzcontinuous (LC) controllers that would track the desired trajectory of the system, and we identify from among these controllers the one that minimizes, at each instant of time, a speciﬁed weighted norm of the control effort. Since we illustrate the power of the technique developed herein to ensure the tracking of a Rossler chaotic system by another Lorenz chaotic system, it is appropriate to indicate the state of the art in synchronization of nonlinear chaotic systems. Nearly all the work that has been done to date on the synchronization of chaotic systems deals with synchronization of identical chaotic oscillators that start from different initial conditions. This synchronization is performed by linear feedback between one or more of the phase states of the system (see Chen 2002; Lei et al. 2005). A recent monograph (with 350 references) points this out in detail (Boccaletti et al. 2002). It is only very Proc. R. Soc. A (2008)

Tracking control of nonlinear systems

2343

recently that the synchronization between two non-identical, low-dimensional chaotic systems has begun to be investigated (Boccaletti et al. 2002). One of the few papers on this is by Pyras (1996), which uses linear feedback. Rulkov et al. (1995) showed that this type of synchronization could exist. Such synchronization studies have resulted in the so-called imperfect phase synchronization (Zaks et al. 1999) in which the systems get ultimately de-synchronized by having intermittent phase slips, and lag synchronization (Rosenblum et al. 1997) which results in a time lag between the synchronized systems. Linear feedback has been used in all these studies, with the aim of synchronizing all the states of one chaotic system with those of the other. In our example we use non-identical chaotic systems. Our methodology does not assume any a priori structure on the control for the synchronization; no restrictions on the chaotic systems being of low dimensionality is required; as illustrated in the example, we can opt to synchronize (at will) one or more of the phase space variables; and, we can obtain the entire set of controllers in closed form that would do the job. Finally, among all these controllers, in the example we show that we can explicitly provide the controller that synchronizes and minimizes a weighted norm of the control effort at each instant of time. We begin by reformulating the trajectory control problem as a problem of constrained motion in the Lagrangian framework and use the underlying inspiration with which constrained motion in analytical dynamics is orchestrated by Nature (Udwadia & Kalaba 2002). We then expand and further develop this view by ﬁrst considering mechanical systems described by ﬁrst-order Poincare´ equations and then general nonlinear systems. Closed-form expressions for all the continuous controllers required for trajectory tracking for nonlinear systems that do not make approximations, either in describing the nonlinear system or in the nature of the nonlinear controller employed, are obtained. Such closed-form results appear to be new. Furthermore, no approximations or linearizations are made here with respect to the trajectory that is being tracked, which may be described in terms of nonlinear algebraic equations, nonlinear differential equations or a combination thereof; these descriptions could explicitly involve time also. Moreover, the approach arrives not just at one nonlinear controller for controlling a given nonlinear system, but also at the entire set of continuous controllers that would cause a given set of trajectory descriptions to be exactly satisﬁed. Furthermore, we show that when the cost function is the weighted norm of the generalized control input, its minimization can be done to yield the optimal controller that minimizes not just the integrated cost over a span of time, but also the cost at each instant of time. Explicit closed-form expressions for the optimal control are obtained. Section 2a of this paper begins with mechanical systems described by the second-order differential equations of the Newtonian and/or Lagrangian mechanics. Section 2b deals with Poincare´’s ﬁrst-order differential equations that describe the motion of mechanical systems; this takes us a step further towards general ﬁrst-order nonlinear dynamical systems. Section 3 deals with the close connection that the set of controllers developed in §2 have with the way Nature seems to orchestrate the control of mechanical systems subjected to trajectory requirements (constraints). The development of such controllers and their close correspondence with (i) Gauss’s principle of least constraint and (ii) the recently developed equations of motions for non-ideal constraints Proc. R. Soc. A (2008)

2344

F. E. Udwadia

are provided. Lessons from the way Nature seems to control mechanical systems are adduced. Section 4 deals with the exact control of general dynamical systems described by a set of nonlinear, non-autonomous ordinary differential equations. Section 5 deals with a numerical example that demonstrates the closed-form development of the control required to be applied to a nonlinear chaotic dynamical system so that it tracks the motions of another different nonlinear chaotic system. Section 6 concludes the paper with some remarks and observations.

2. Development of the entire set of controllers that cause a mechanical system to track a given trajectory (a ) Lagrangian and Newtonian descriptions Consider an unconstrained nonlinear mechanical system described by the secondorder differential equation of motion _ Z q_0 ; _ tÞ; ð2:1Þ M ðq; tÞ€ q Z Qðq; q; qð0Þ Z q0 ; qð0Þ where q(t) is the n-vector (n by 1 vector) of generalized coordinates; the dots indicate differentiation with respect to time; and the matrix M(q, t) is a positivedeﬁnite n by n matrix. Equation (2.1) can be obtained using either Newtonian or Lagrangian mechanics (Lagrange 1811; Hamel 1949; Goldstein 1976). The n-vector Q on the r.h.s. of equation (2.1) is a ‘known’ vector in the sense that it is a known function of its arguments. By ‘unconstrained’ we mean that the components of the initial velocity q_ 0 of the system can be independently assigned. We next require that this mechanical system be controlled so that it tracks a trajectory that is described by the following consistent set of m equations: fi ðq; tÞ Z 0; i Z 1; .; h ð2:2Þ and _ tÞ Z 0; ji ðq; q; i Z h C 1; .; m: ð2:3Þ We shall assume that the mechanical system’s initial conditions are such as to satisfy these relations at the initial time. The latter set of equations, which are non-integrable, is non-holonomic (Hamel 1949). In order to control the system so that it exactly tracks the required trajectory—i.e. satisﬁes equations (2.2) and (2.3)—we apply an appropriate _ tÞ, so that the equation of motion of the controlled control n-vector Q c ðq; q; system becomes _ tÞ C Q c ðq; q; tÞ; M ðq; tÞ€ q Z Qðq; q;

qð0Þ Z q0 ;

_ Z q_0 ; qð0Þ

ð2:4Þ

where now the components of the n-vectors q0 and q_0 satisfy equations (2.2) and (2.3) at the initial time, tZ0. Throughout this paper, we shall, for brevity, drop the arguments of the various quantities, unless needed for clarity. We begin by expressing equation (2.4) in terms of the weighted accelerations of the system. For any positive-deﬁnite n by n matrix N(q, t), we deﬁne the matrix Gðq; tÞ d½N 1=2 ðq; tÞM ðq; tÞ K1 Z M K1 ðq; tÞN K1=2 ðq; tÞ; Proc. R. Soc. A (2008)

ð2:5Þ

Tracking control of nonlinear systems

2345

and pre-multiplying equation (2.4) by N 1/2 (q, t), the ‘scaled’ equation, which we denote using the subscript ‘s’, is obtained as ð2:6Þ q€s Z as Cq€sc ; where € ð2:7Þ q€s dG K1q; a s dG K1 a Z ðN 1=2 M ÞðM K1 QÞ Z N 1=2 Q

ð2:8Þ

and ð2:9Þ q€sc dG K1q€c Z ðN 1=2 M ÞðM K1 Q c Þ Z N 1=2 Q c : In equation (2.8), we denote the acceleration of the uncontrolled system by _ tÞZ M K1 ðq; tÞQðq; q; _ tÞ: In equation (2.9), q€c ðq; q; _ tÞZ M K1 ðq; tÞQ c ðq; q; _ tÞ aðq; q; can be viewed as the deviation of the acceleration of the controlled system from that of the uncontrolled system. We now differentiate equation (2.2) twice with respect to time t, and equation (2.3) once with respect to time, giving the set of equations _ tÞ€ _ tÞ; Aðq; q; q Z bðq; q; ð2:10Þ where A is an m by n matrix of rank k and b is an m-vector. Noting equation (2.7), equation (2.10) can be further expressed as _ tÞ€ _ tÞ; Bs ðq; q; q s Z bðq; q; ð2:11Þ _ tÞ dAðq; q; _ tÞGðq; tÞ. We now express the n-vector q€s in terms of its where Bs ðq; q; orthogonal projections on the range space of BsT and the null space of Bs, so that q€s Z BsCBsq€s C ðI KBC q s: s Bs Þ€

ð2:12Þ

In equation (2.12), the matrix X denotes the Moore–Penrose (MP) generalized inverse of the matrix X (Moore 1910; Penrose 1955). It should be noted that equation (2.12) is a general identity that is always valid since it arises from the orthogonal partition of the identity matrix I Z BsCBs C ðI KBsCBs Þ. Using equation (2.11) in the ﬁrst member on the r.h.s. of equation (2.12), and equation (2.6) to replace q€s in the second member, we get ð2:13Þ q€s Z BsCb C ðI KBsCBs Þða s Cq€sc Þ; which, owing to equation (2.6), yields ð2:14Þ BsCBsq€sc Z BC s ðbK Bs a s Þ: But the general solution of the linear set of equations (2.14) is given by (Graybill 2001) C C C C C q€sc Z ðBC s Bs Þ Bs ðbK Bs a s Þ C ½I KðBs Bs Þ ðBs Bs Þz C

Z BsCðbK Bs a s Þ C ðI KBsCBs Þz;

ð2:15Þ

_ tÞ is any arbitrary n-vector. In the second equality where the n-vector zðq; q; above, we have used the property that ðBsCBs ÞCZ BsCBs in the two members on the r.h.s., along with the second MP-inverse property that BsCBs BsCZ BsC. Using equation (2.9), the explicit control force that exactly tracks the trajectory by exactly satisfying the given relations (2.2) and (2.3) is then explicitly given by K1=2 Q c Z N K1=2q€sc Z N K1=2 BC ðI KBsCBs Þz; s ðbK Bs a s Þ C N

Proc. R. Soc. A (2008)

ð2:16Þ

2346

F. E. Udwadia

_ tÞ dAðq; q; _ tÞGðq; tÞ. We may take zðq; q; _ tÞ to be C 1 (or, more where Bs ðq; q; generally, Lipschitz continuous (LC)) to ensure a unique solution to the system of equations (2.4). We hence have the following result. Result 2.1. Consider the mechanical system, which is described by the Lagrange (or Newtonian) equations of motion _ tÞ; M ðq; tÞ€ q Z Qðq; q;

ð2:17Þ

where M is an n by n positive-deﬁnite matrix and q is an n-vector. The system is required to exactly track the trajectory described by the equations fi ðq; tÞ Z 0;

i Z 1; .; h

ð2:18Þ

and _ tÞ Z 0; ji ðq; q;

i Z h C 1; .; m:

ð2:19Þ

The controlled system is described by the relation _ tÞ C Q c ðq; q; _ tÞ; M ðq; tÞ€ q Z Qðq; q;

ð2:20Þ

where Q c is the control. Assuming that the initial conditions of the mechanical system satisfy these _ tÞ (or controllers) trajectory requirements, the set of all possible controls Q c ðq; q; that causes the controlled system (2.20) to exactly track the required trajectory is explicitly given by Q c Z N K1=2 BsCðbKBs a s Þ C N K1=2 ðI KBC s Bs Þz;

ð2:21Þ

_ tÞ is any arbitrary n-vector whose components are continuously where zðq; q; differentiable—or LC—functions of its arguments; N(q, t) is any arbitrary n by n _ tÞZ Aðq; q; _ tÞ½N 1=2 ðq; tÞM ðq; tÞ K1 is positive-deﬁnite matrix a s Z N 1=2 Q; Bs ðq; q; _ tÞ is an m by n matrix of rank k, and bðq; q; _ tÞ is the an m by n matrix; Aðq; q; m-vector deﬁned in equation (2.10). _ tÞ given We can abbreviate the two components of the control vector Q c ðq; q; in equation (2.21) as _ tÞ Z Q1c ðq; q; _ tÞ C Q2c ðq; q; _ tÞ; Q c ðq; q;

ð2:22Þ

_ tÞ dN K1=2 BC Q1c ðq; q; s ðbK Bs a s Þ;

ð2:23Þ

_ tÞ dN K1=2 ðI KBsCBs Þz: Q2c ðq; q;

ð2:24Þ

where

and

Corollary 2.2. Relation (2.23) can also be given as Q1c Z N K1=2 ðAM K1 N K1=2 ÞCðbKAaÞ:

ð2:25Þ

Proof. Using the relation Bs Z AG Z AM N and equation (2.8), we have _ tÞZ M K1 ðq; q; _ tÞQðq; q; _ tÞ is the acceleration Bs a s Z AGG K1 aZ Aa, where aðq; q; of the uncontrolled system. The result then follows. & K1

Proc. R. Soc. A (2008)

K1=2

Tracking control of nonlinear systems

2347

Remark 2.3. Owing to the generality of the decomposition (2.13), relation (2.21) (alternatively, relations (2.22)–(2.24)) provides the entire set of continuous tracking controllers that cause the system to track the trajectory described by equations (2.18) and (2.19). Remark 2.4. The explicit closed-form tracking control obtained in equation (2.21), which causes the system to exactly track the given trajectory described by equations (2.18) and (2.19), does not contain any Lagrange multipliers, nor does our derivation invoke the notion of a Lagrange multiplier. We shall next show the following, somewhat remarkable, result. Result 2.5. We again consider the mechanical system described by the nonlinear Lagrange or Newtonian equation (2.17), which needs to be controlled _ tÞ, so that the trajectory through the addition of a control, n -vector Q c ðq; q; described by equations (2.18) and (2.19) is exactly tracked. Assuming that the system satisﬁes the trajectory requirements initially, the optimal controller that causes the system to (i) exactly track the required trajectory and (ii) minimize at each instant of time t, the cost _ tÞT N ðq; tÞQ c ðq; q; _ tÞ; J ðt Þ Z ½Q c ðq; q;

ð2:26Þ

for a given n by n positive-deﬁnite matrix N, is explicitly provided by K1=2 C _ tÞ Z N K1=2 BC Bs ðbKAaÞ; Q c ðq; q; s ðbKBs a s Þ Z N

ð2:27Þ

where Bs, b and as are as deﬁned before. Proof. Let us deﬁne the vector rðtÞ dN 1=2 Q c :

ð2:28Þ

We note that by equation (2.26), J ðtÞ Z ½Q c ðq; q; tÞT N ðq; tÞQ c ðq; q; tÞ Z krðt Þk2 :

ð2:29Þ

Then in view of equation (2.20), equation (2.28) can be rewritten as rðt Þ dN 1=2 ðM q€ KQÞ;

ð2:30Þ

so that we obtain, using equation (2.8), q€ Z ½N 1=2 M K1 ðr C N 1=2 QÞ Z Gðr C a s Þ:

ð2:31Þ

We assume that at the initial time the values of q0 and q_ 0 satisfy the trajectory requirements, and upon differentiating equations (2.18) and (2.19) we have, as before, Aq€ Z b:

ð2:32Þ

Using relation (2.31) and denoting BsZAG, equation (2.32) becomes Bs r Z bKBs a s : Proc. R. Soc. A (2008)

ð2:33Þ

2348

F. E. Udwadia

The solution of equation (2.33), subject to the condition that J ðtÞZ krðtÞk2 is a minimum, is then given by (Graybill 2001) rðt Þ Z BsCðbK Bs a s Þ;

ð2:34Þ

from which we obtain, by using relation (2.28), the explicit expression for the optimal control as _ tÞ Z ½N ðq; tÞ K1=2 BsCðq; q; _ tÞ½bðq; q; _ tÞK Bs ðq; q; _ tÞa s ðq; q; _ tÞ; Q c ðq; q;

ð2:35Þ

where we have written out the explicit result in extensio. The second equality in equation (2.27) follows from corollary 2.2. & Remark 2.6. For a given choice of the weighting matrix N, we have obtained the explicit closed-form expression for the exact full state controller that tracks the trajectory described by equations (2.18) and (2.19) and minimizes the cost function J ðtÞZ ½Q c T N ðq; tÞQ c , under the proviso that the initial conditions of the system satisfy the trajectory description. The minimum cost J(t) is given by 2 K1 K1=2 C J ðtÞ Z ½Q1c T NQ1c Z kBC N Þ ðbKAaÞk2 : s ðbKAaÞk Z kðAM

ð2:36Þ

Note, as before, that the explicit result given in equation (2.27) does not contain any Lagrange multipliers, nor does our derivation invoke anywhere the notion of a Lagrange multiplier. Remark 2.7. Result 2.5 and equation (2.23) show that for a given positivedeﬁnite matrix N(q, t) the optimal control that minimizes J(t) at each instant of _ tÞ h 0 in equation time is explicitly obtained in closed form by setting Q2c ðq; q; (2.22). One way (see corollary 2.8 below) of doing this would be by setting _ tÞ h 0 in equation (2.24) (or, in equation (2.21)). More precisely, we have zðq; q; the following result. _ tÞ Corollary 2.8. At any instant of time t, at which the arbitrary LC vector zðq; q; _ tÞ, the control given by (2.21) becomes optimal, belongs to the range space of BsT ðq; q; in the sense that it minimizes J(t) at that time. _ tÞ belongs to the range space of BsT ðq; q; _ tÞ at time t, we can Proof. When zðq; q; T _ tÞ becomes, give z Z Bs g for some m-vector g. Hence, by equation (2.24), Q2c ðq; q; Q2c dN K1=2 ðI KBsCBs Þz Z N K1=2 ½I KBsCBs BsT g T T K1=2 T T ½I KBsT ðBC Z N K1=2 ½I KðBC s Bs Þ Bs g Z N s Þ Bs g Z 0:

ð2:37Þ

The second equality follows from the fourth MP-inverse property, and the last equality follows from the ﬁrst MP-inverse property (Graybill 2001). & Corollary 2.9. (i) The two components Q1c and Q2c of the control vector Q given in equation (2.22) are N-orthogonal to one another and (ii) Q1c is (MN )-orthogonal to the null space of the matrix A. c

Proc. R. Soc. A (2008)

Tracking control of nonlinear systems

2349

Proof. (i) Since ½Q2c T NQ1c Z z T ðI KBsCBs ÞT N K1=2 NN K1=2 BC s ðbK Bs a s Þ Z z T ðI KBsCBs ÞBC s ðbK Bs as Þ Z 0:

ð2:38Þ

The result follows by using the fourth and then the second MP-inverse property (Graybill 2001). (ii) Substituting for Bs and using the identity X CZ X T ðXX T ÞC we get K1=2 ðAM K1 N K1=2 ÞC N K1=2 BC s ZN

Z N K1=2 ðAM K1 N K1=2 ÞT ½ðAM K1 N K1=2 ÞðAM K1 N K1=2 ÞT C Z ðMN ÞCAT ðAM K1 N K1 M K1 AT ÞC:

ð2:39Þ

Any vector belonging to the null space of the matrix A satisﬁes the relation AvZ0, whose explicit solution is v Z ðI KACAÞw;

ð2:40Þ

where w is any n-vector. Hence v T MNQ1c Z w T ðI KACAÞT ðMN ÞðMN Þ K1 AT ðAM K1 N K1 M K1 AT ÞCðbK Bs a s Þ Z w T ½I KAT ðAT ÞCAT ðAM K1 N K1 M K1 AT ÞCðbK Bs a s Þ Z 0; since ½I KAT ðAT ÞCAT Z 0.

ð2:41Þ &

Q1c

Corollary 2.10. The component as deﬁned in (2.23) belongs to the range K1=2 C space of the matrix N Bs , and the component Q2c as deﬁned in equation (2.24) is N-orthogonal to the range space of N K1=2 BsT . Proof. The ﬁrst result is obvious from the expression for Q1c in equation (2.23). The second follows because any vector a in the range space of N K1=2 BsT can be expressed as a Z N K1=2 BsT g;

ð2:42Þ

aT NQ2c Z gT Bs N K1=2 NN K1=2 ðI KBC s Bs Þz Z 0:

ð2:43Þ &

for some m-vector g, and so Hence the result.

We note that the range space of N K1=2 BsC is the same as the range space of N BsT . This follows because for any vector g, we can always ﬁnd a vector a such that N K1=2 BsT aZ N K1=2 BsCg, and vice-versa. K1=2

Remark 2.11. The addition of the second component Q2c to the control so that Q c Z Q1c C Q2c , whatever the LC n-vector z may be, still ensures that the description of the trajectory given by equations (2.18) and (2.19) is exactly satisﬁed. However, it contributes an additional amount given by kðI KBsCBs Þzk2 to the cost of the control, since the two components Q1c and Q2c are N-orthogonal Proc. R. Soc. A (2008)

2350

F. E. Udwadia

_ tÞ to each other (corollary 2.9). Furthermore, if at any time t the function zðq; q; belongs to the null space of Bs this additional cost at that time simply becomes kzk2 . It is sometimes advantageous to describe the motion in terms of a system of ﬁrst-order differential equations instead of a system of second-order equations. _ so One obvious way of doing this is to deﬁne a new variable, an n-vector v Zq, that the equation of motion takes the so-called state space form given by q_ Z v; qð0Þ Z q0 ; ð2:44aÞ M ðq; tÞv_ Z Qðq; v; tÞ;

vð0Þ Z v0 :

ð2:44bÞ

The ﬁrst equation of this set is simply a deﬁnition, whereas the second equation of the set is the one that contains the actual dynamics. And while this is the form often used when numerically integrating equation (2.1), the description of the motion of a mechanical system in terms of a system of ﬁrst-order differential equations goes far beyond just its use in numerical procedures. For, often these ﬁrst-order descriptions (i.e. descriptions using ﬁrst-order differential equations) arise when one wants to use descriptions of motion in terms of coordinates that may be more physically meaningful in the context of a particular problem. For example, the ﬁrst-order Hamilton’s equations describing the motion of mechanical systems are often useful when dealing with systems described by a Hamiltonian, and, more generally, the ﬁrst-order set of Poincare´ equations are often useful when dealing with rigid body and multi-body dynamics. We next obtain the entire set of explicit closed-form controllers for exactly tracking the trajectory of mechanical systems described by Poincare´’s equations of motion. (b ) Poincare´ descriptions The Poincare´ equations (Poincare´ 1901) are obtained by deﬁning a new _ where H~ ðq; tÞ is a non-singular matrix whose variable, an n-vector sZH~ ðq; tÞq, elements are known functions of q and t. Denoting the inverse of the matrix H~ by H(q, t), so that q_ Z H ðq; tÞs, the Lagrange equation (2.1) can be rewritten in ﬁrst-order form as (Udwadia & Phohomsiri 2007) q_ Z H ðq; tÞs; qð0Þ Z q0 ;

ð2:45Þ

Mp ðq; tÞ_s Z Sðq; s; tÞ; sð0Þ Z s0 ;

ð2:46Þ

where the state variables are now the n-vectors q and s. In rigid body dynamics, the 3-vector s is often chosen to constitute the three components of the angular velocity of the body. Here the matrix Mp is again positive deﬁnite. Though both ﬁrst-order equations are, mathematically speaking, on a par, the ﬁrst equation of the set (equation (2.45)) is a kinematic relation involving the deﬁnition of the new variable s—hence, an identity as in equation (2.44a)—while the second equation (equation (2.46)) is the dynamical equation of motion. The system is required to track the trajectory described by the consistent relations fpi ðq; tÞ Z 0; jpi ðq; s; tÞ Z 0; Proc. R. Soc. A (2008)

i Z 1; .; h;

ð2:47Þ

i Z h Z 1; .; m:

ð2:48Þ

Tracking control of nonlinear systems

2351

The subscripts and superscripts ‘p’ indicate that we are dealing with the Poincare´ description of the motion of the dynamical system. In order to control the system so that it exactly tracks this trajectory, we apply a control so that the equation of motion of the controlled system now becomes q_ Z H ðq; tÞs; qð0Þ Z q0 ; ð2:49Þ Mp ðq; tÞ_s Z Sðq; s; tÞ C S c ðq; s; tÞ;

sð0Þ Z s0 ;

ð2:50Þ

and the initial conditions q0 and s0 satisfy the trajectory requirements (2.47) and (2.48). The explicit expression for the control force, yielding the entire set of all controllers that will cause the dynamical system described by equations (2.45) and (2.46) to track the required trajectory described by equations (2.47) and (2.48) will now be obtained. As before, we deﬁne the matrix Gp ðq; tÞ d½N K1=2 ðq; tÞMp ðq; tÞ K1 Z MpK1 ðq; tÞN K1=2 ðq; tÞ; ð2:51Þ where the matrix N(q, t) is any positive-deﬁnite matrix, and the scaled variables 9 _ s_s Z GpK1s; > > > > = K1 1=2 K1 1=2 a s Z Gp a Z ðN Mp ÞðMp SÞ Z N S ð2:52Þ > and > > > ; s_cs Z N 1=2 S c ; so that relation (2.50) on pre-multiplication with N 1/2 becomes s_s Z a s Cs_sc ; ð2:53Þ which we note is of the same form as equation (2.6), except that we now have s instead of q, and ﬁrst derivatives instead of second derivatives with respect to time. On differentiating equation (2.47) twice with respect to time and differentiating equation (2.48) once with respect to time, and using relation (2.49), we then obtain the matrix equation Ap ðq; s; tÞ_s Z bp ðq; s; tÞ;

ð2:54Þ

which, upon using the ﬁrst equation of the set (2.52), can be rewritten as Bs ðq; s; tÞ_ss Z bp ðq; s; tÞ;

ð2:55Þ

where we deﬁne the m by n matrix Bs of rank k by the relation Bs ðq; s; tÞ Z Ap ðq; s; tÞGp ðq; tÞ:

ð2:56Þ

Again we note the similarity between equations (2.10) and (2.11) and equations (2.54) and (2.55). We next decompose the n-vector s_s as s_s Z BsCBss_s C ðI KBC s Bs Þs_s ;

ð2:57Þ

in a manner similar to equation (2.12), with the matrix Bs now deﬁned by equation (2.56). Proceeding along with similar lines as before (equations (2.13)–(2.16)), we then ﬁnd that K1=2 ðI KBsCBs Þz p dS c1 C S c2 ; ð2:58Þ Sc Z N K1=2s_sc Z N K1=2 BC s ðbp K Bs a s Þ C N

Proc. R. Soc. A (2008)

2352

F. E. Udwadia

where z p ðq; s; tÞ is any arbitrary n -vector. To ensure a unique solution of equations (2.49) and (2.50), we may then take the components of z p ðq; s; tÞ to be C 1 functions (or, more generally, LC functions) of q, s and t. Following the same lines as in §2a, for any Poincare´ system described by equation (2.45) that is required to exactly track the trajectory described by relations (2.47) and (2.48), we now obtain results that are analogous to those given in results 2.1 and 2.5. Using as deﬁned in relation (2.52), and Bs in (2.56), _ s, Q/ S, b/ bp , M / Mp , we simply make the following variable changes: q/ z / z p , A/ Ap , and G / Gp , in the expressions given in equations (2.21) and (2.27) to obtain the corresponding control forces. Remark 2.12. We can prove corollaries similar to corollaries 2.2, 2.8–2.10. The same goes for remarks 2.3, 2.4, 2.6, 2.7, 2.11.

3. The close connection between nonlinear control and analytical dynamics Since almost all mechanical systems are nonlinear in their behaviour, including even simple ones like a pendulum, we will be mainly addressing nonlinear systems here. The problem of control can be placed within the context of analytical mechanics by reinterpreting constrained motion in mechanical systems. Consider a mechanical system described by equation (2.17). When the system is further subjected to the trajectory requirements (2.18) and (2.19)—i.e. subjected to further constraints—additional control (constraint) forces are brought into play by Nature so that the controlled (constrained) system moves in such a manner that it satisﬁes these trajectory requirements (constraints). Thus, the additional _ tÞ that Nature provides may be thought of as the control (constraint) Q c ðq; q; control it generates in order for the system to satisfy the trajectory requirements given by equations (2.18) and (2.19). One might imagine Nature as a control engineer, attempting to control the mechanical system so that it satisﬁes the given trajectory requirements described by equations (2.18) and (2.19). However, to entertain such an interpretation, one is led to ask the following three questions. (i) To what extent might Nature be perceived as acting like a control engineer? (ii) Does Nature appear to be performing the control in any kind of an optimal manner? (iii) If so, what is the cost function (or functional) it appears to minimize? We shall start answering these questions by ﬁrst asserting that Nature appears to choose the weighting matrix N ðq; tÞZ ½M ðq; tÞ K1 dNNature ðq; tÞ in equation (2.21). Here M is the so-called ‘mass matrix’ and appears on the l.h.s. of equation (2.17). Result 3.1. Nature seems to control a mechanical system described by equation (2.17), so it exactly satisﬁes the trajectory described by the requirements (2.18) and (2.19) by choosing the weighting matrix to be N ðq; tÞ Z ½M ðq; tÞ K1 : Proc. R. Soc. A (2008)

ð3:1Þ

Tracking control of nonlinear systems

2353

Proof. Setting N ðq; tÞZ ½M ðq; tÞ K1 equation (2.21) now yields _ tÞ Z M 1=2 BCðbKAaÞ C M 1=2 ðI KBCBÞz; Q c ðq; q;

ð3:2Þ

where B dAG Z AM K1=2 and Bs a s Z AGG K1 aZ Aa. But equation (3.2) is exactly the equation of motion of a general constrained mechanical system, as given by Udwadia (2000) and Udwadia & Kalaba (2002). Hence the result. & Result 3.2. If we assume that Nature observes d’Alembert’s principle, then it appears to be minimizing the cost _ tÞT ½M ðq; tÞ K1 Q c ðq; q; _ tÞ; JNature ðt Þ Z ½Q c ðq; q; ð3:3Þ at each instant of time while controlling the system deﬁned by equation (2.17) so it exactly satisﬁes the trajectory described by the requirements (2.18) and (2.19). Proof. Setting N ðq; tÞZ ½M ðq; tÞ K1 in result 2.5, equation (2.27), which gives the optimal control while minimizing J(t), yields _ tÞ Z M 1=2 BCðbKAaÞ Z Q1c : Q c ðq; q; ð3:4Þ c Using this expression for Q in equation (2.20), we ﬁnd that we obtain the correct equation of motion of a constrained mechanical system that obeys d’Alembert’s principle, as given by Udwadia & Kalaba (1992, 1996). & Remark 3.3. Result 3.2 connects directly with the basic principles of analytical dynamics. In fact, d’Alembert’s principle, which is a principle that leads to a mathematical description of motion, which is in close conformity with observations, and which is one of the pivotal assumptions of analytical dynamics, is equivalent to Gauss’s principle (Gauss 1829; Udwadia & Kalaba 1996). And Gauss’s principle states that: of the entire set of constraint (control) forces that cause a constrained (controlled) mechanical system (2.17) to exactly satisfy the constraints (trajectory) described by requirements (2.18) and (2.19), Nature seems to choose that constraint force that minimizes JNature (t) at each instant of time. Thus we ﬁnd that (i) Nature seems to choose the weighting matrix N ðq; tÞZ ½M ðq; tÞ K1 and (ii) if we assume that d’Alembert’s principle is true, then Nature seems to pick the one controller given in result 3.2 that minimizes the cost JNature (t). Nature appears to go well beyond what most modern control engineers would try to do, by minimizing the cost JNature given in relation (3.3) at each instant of time, rather than minimizing the integral of this cost over any given span of time, as is the common practice in the ﬁeld of controls. Remark 3.4. When d’Alembert’s principle is assumed to be satisﬁed, Nature appears to be performing acceleration feedback control, where the control force is given by Q c Z Q1c Z M 1=2 BCðbKAaÞ ZKKNature ðAaKbÞ: ð3:5Þ _ tÞ dM 1=2 ðq; tÞBCðq; q; _ tÞ and the It is made up of the gain matrix KNature ðq; q; feedback error e(t)d(AaKb). The quantity e(t) is simply the extent to which _ tÞ does not satisfy the the acceleration of the uncontrolled system aðq; q; trajectory requirements imposed on it by relation (2.10). Hence, Nature appears to be behaving like a control engineer, using feedback control. The minus sign in equation (3.5) is taken so as to be compatible with the control concept of negative feedback. While most modern control engineers usually use integral, Proc. R. Soc. A (2008)

2354

F. E. Udwadia

velocity and proportional feedback in mechanical systems, few use acceleration feedback as Nature appears to be using. This acceleration feedback does not _ Þ, require the measurement of acceleration because a is a function of qðt Þ and qðt _ tÞZM K1 ðq; tÞQðq; q; _ tÞ. and can be obtained from their measurement, since aðq; q; Also, the gain matrix KNature used by Nature is complex and its elements are, in general, highly nonlinear functions of q, q_ and t. Remark 3.5. The component Q1c of the constraint force used by Nature is always orthogonal to the null space of the matrix A given in equation (2.10). By corollary 2.9(ii), Q1c is always MN-orthogonal to the null space of A. Since Nature picks N ðq; tÞZ ½M ðq; tÞ K1 , the result follows. The null space of A in analytical dynamics is called the space of virtual displacements and d’Alembert’s principle simply posits the assumption that the control (constraint) force n-vector Q1c is orthogonal to the null space of A, or in more analytical dynamics terms: the ‘work done’ by the constraint force Q c under virtual displacements at each instant of time is zero. Remark 3.6. Many mechanical systems, however, do not satisfy d’Alembert’s principle (Goldstein 1976). In such systems, the control (constraint force) does work under virtual displacements. To obtain the equation of motion for such systems one requires additional information about the work done by the control forces under virtual displacements. One then needs to prescribe, for a speciﬁc _ tÞ such that mechanical system, the C 1 vector C ðq; q; _ tÞ Z w T ðtÞC ðq; q; _ tÞ; w T ðtÞQ c ðq; q;

ð3:6Þ

where w(t) is any virtual displacement, i.e. any n-vector in the null space of A. In that case, the equation of motion of the constrained mechanical system is known to be described by the relation (Udwadia & Kalaba 2002) _ tÞ C M 1=2 BCðbKAaÞ C M 1=2 ðI KBCBÞM K1=2 C: Mq€ Z Qðq; q;

ð3:7Þ

This same result also follows directly from equations (2.20) and (2.21) by setting N ðq; tÞZ ½M ðq; tÞ K1 and zðt ÞZ M K1=2 C . Thus, if a mechanical system is non_ tÞ is ideal, the non-idealness being described by relation (3.6) in which C ðq; q; speciﬁed at each instant of time, then Nature chooses the n-vector z(t) in relation (2.21) to be zðt ÞZ M K1=2 C , so that condition (3.6) is satisﬁed along with relations (2.18) and (2.19) at each instant of time for the speciﬁc system at hand!

4. General systems described by ﬁrst-order differential equations In this section, we proceed to general dynamical systems that are described by n ﬁrst-order, non-autonomous, differential equations given by Mg ðx; tÞx_ Z f ðx; tÞ; xð0Þ Z x 0 ; ð4:1Þ where we shall again take Mg to be a positive-deﬁnite n by n matrix. We require that this dynamical system track a trajectory described by the equations fi ðx; tÞ Z 0; i Z 1; .; m; ð4:2Þ where we assume that the trajectory described is feasible and the system of m equations is consistent. In order to track this trajectory, we apply a control f c so Proc. R. Soc. A (2008)

Tracking control of nonlinear systems

2355

that the equation describing the time evolution of the controlled dynamical system becomes ð4:3Þ Mg ðx; tÞx_ Z f ðx; tÞ C f c ðx; tÞ; xð0Þ Z x 0 ; where we assume that the initial conditions satisfy the trajectory requirements (4.2). As before, we differentiate these m equations (4.2) with respect to time to give the relation Agx_ Z bg ðx; tÞ; ð4:4Þ where Ag(x, t) is an m by n matrix of rank k, and set Gg Z ½N 1=2 ðx; tÞMg ðx; tÞ K1 : Pre-multiplying equation (4.3) by N 1/2 we get x_s Z a s Cx_cs ; where _ x_s Z GgK1x;

a s Z GgK1 a Z ðN 1=2 Mg ÞðMgK1 f Þ Z N 1=2 f

ð4:5Þ ð4:6Þ and x_cs Z N 1=2 f c :

ð4:7Þ Here, aZ MgK1 f : We note that equations (4.6) and (4.7) have the same form as equations (2.53) and (2.52), respectively. Furthermore, equation (4.4) can be expressed as Bs ðx; tÞx_ s Z bg ðx; tÞ; ð4:8Þ where Bs ðx; tÞ Z Ag ðx; tÞGg ðx; tÞ: ð4:9Þ Expressing x_s in terms of orthogonal components, we get x_s Z BsCBsx_ s C ðI KBsCBs Þx_ s ; ð4:10Þ c and replacing x_ s on the r.h.s. of (4.10) by as Cx_ s yields ð4:11Þ BsCB sx_ s Z BsCðbg K Bs a s Þ; from which we get the following result by following the same lines as in equations (2.13)–(2.16). Result 4.1. Consider the nonlinear, non-autonomous dynamical system (4.1). The system is required to track the trajectory described by equation (4.2). Assuming that the initial conditions satisfy the trajectory described by equation (4.2), all the possible LC continuous controls that exactly track this trajectory are explicitly given by K1=2 f c Z N K1=2 BC ðI KBC ð4:12Þ s ðbg K Bs a s Þ C N s Bs Þz g ; where zg(x, t) is any arbitrary n-vector whose components are continuously differentiable (or more generally are LC) functions of their arguments; N(x, t) is any arbitrary n by n positive-deﬁnite matrix, a s Z N 1=2 ðx; tÞf ðx; tÞ; Bs ðx; tÞZ Ag ðx; tÞ½N 1=2 ðx; tÞMg ðx; tÞ K1 is an m by n matrix; Ag is the m by n matrix of rank k, and bg(x, t) is the m-vector deﬁned in equation (4.4).

Result 4.2. Consider the nonlinear dynamical system described by equation (4.1) that needs to be controlled through the addition of a control f c ðx; tÞ, so that the trajectory described by equation (4.2) is exactly tracked. Assuming that the system satisﬁes the trajectory requirements initially, the optimal controller that Proc. R. Soc. A (2008)

2356

F. E. Udwadia

causes the system to (i) exactly track the required trajectory and (ii) minimize at each instant of time t, the cost J ðt Þ Z ½ f c ðx; tÞT N ðx; tÞf ðx; tÞ; ð4:13Þ for a given n by n positive-deﬁnite matrix N, is explicitly provided by ð4:14Þ f c ðx; tÞ Z N K1=2 BC s ðbg K Bs a s Þ; where Bs, bg and a s are deﬁned in relations (4.9), (4.8) and (4.7), respectively. Proof. Similar to, and along with the same lines as, the proof of result 2.5.

&

Denoting f c Z f1c C f2c ;

ð4:15Þ

where f1c dN K1=2 BsCðbg K Bs a s Þ Z N K1=2 ½Ag MgK1 N K1=2 Cðbg K Ag aÞ

ð4:16Þ

and f2c dN K1=2 ðI KBsCBs Þz g ; we have the following remarks.

ð4:17Þ

Remark 4.3. We can prove corollaries similar to corollaries 2.2, 2.8–2.10. Four important results that emerge from this are (i) the two components f1c and f2c of the control vector f c are N-orthogonal to one another, (ii) f1c is (MN )-orthogonal to the null space of the matrix A, (iii) the minimum cost is J ðt ÞZ kBsCðbg K Ag aÞk2 and it occurs at those times when f2c Z0, and (iv) the control cost contributed by f2c is given by kðI KBsCBs Þz g k2 . Remark 4.4. So far, it has been assumed that the equations that describe the trajectory are consistent. This may not happen in practical situations; errors due to numerical computations, for example, could make these equations inconsistent. Hence, instead of equation (4.8) we would have the equation Bs x_s Z bg C 3ðx; tÞ;

ð4:18Þ

where 3ðx; tÞ is the error caused by the inconsistency of the trajectory equation (4.2) (as also, analogously, for the equation sets (2.18)–(2.19), and (2.47)–(2.48)). We would then need to replace equation (4.11) with BsCBsx_cs Z BC ð4:19Þ s ðbK Bs a s Þ C dðx; tÞ; C where dZ Bs 3 is the error. The least-squares solution to this inconsistent equation remains x_ cs Z B sCðbK Bs a s Þ C ðI KBsCBÞz; ð4:20Þ as before, pointing out that results 4.1 and 4.2 (and similarly results 2.1 and 2.5) provide the control in this least-squares sense, even when the trajectory description is inconsistent. 5. Example In this section, we present an application of the results obtained by considering the trajectory tracking of a chaotic Rossler system (Rossler 1976) by a Lorenz system (Strogatz 1994). We begin by considering the scaled Lorenz oscillator as a Proc. R. Soc. A (2008)

2357

Tracking control of nonlinear systems

system that we would like to control. Its description is given by the three ﬁrstorder differential equations, 2 3 2 3 x_ 1 sðx 2 K x 1 Þ 6 7 6 7 x_ d 4 x_ 2 5 Z d4 rx 1 K x 2 K x 1 x 3 5 Z f ðxÞ; x 1 ð0Þ Z x 3 ð0Þ Z 5; x 2 ð0Þ Z 10; x 1 x 2 Kbx 3

x_ 3

ð5:1Þ where we take sZ10, rZ28, bZ8/3 and dZ10. This system exhibits chaotic motion for the chosen values of the parameters (Strogatz 1994). We want the motion of the ﬁrst two components, x 1(t) and x 2(t), of this chaotic Lorenz system to track the ﬁrst two components, y1(t) and y2(t), respectively, of a very different chaotic system—a Rossler system—which is described by the equations (Rossler 1976), 2 3 2 3 Kðy2 C y3 Þ y_1 6 7 6 7 y_ d 4 y_2 5 Z 4 y1 C ay2 5 Z hðyÞ; y1 ð0Þ Z 3; y2 ð0Þ Z 12; y 3 ð0Þ Z 6; y_3

d C y3 ðy1 KcÞ

ð5:2Þ with aZ0.1, cZ18 and dZ0.3. For these parameter values, the Rossler system is also known to be chaotic (Strogatz 1994). Thus, the aim is to ﬁnd the control inputs needed to be applied to one chaotic system (the Lorenz system here) so that it tracks two of the components of the motion of another different chaotic system (the Rossler system). The control input 3-vector f c that needs to be applied to the Lorenz system is required to be found so that the cost J ðt ÞZ ½f c T NL f c is minimized at each instant of time. The weighting matrix NL is taken to be the diagonal matrix, NL Z Diag½x 3 ðt Þ2 C 1; x 2 ðt Þ2 C 1; x 1 ðt Þ2 C 1. This weighting matrix is guaranteed to be positive deﬁnite. A simple way to formulate this nonlinear trajectory tracking problem is to consider the augmented dynamical system, # " # " f ðxÞ x_ : ð5:3Þ Z y_ hðyÞ Our task would be to ﬁnd a control input 6-vector F c Z ½ðf c ÞT ; ðh c ÞT T which causes the system # # " c " # " f ðx; yÞ f ðxÞ x_ ð5:4Þ C c Z y_ h ðx; yÞ hðyÞ to satisfy the following trajectory tracking requirements: ðiÞ f1 ðx; yÞ dx 1 ðt ÞK y1 ðt Þ Z 0;

f2 ðx; yÞ dx 2 ðt ÞK y2 ðt Þ Z 0

and

ð5:5Þ

ðiiÞ y_ Z hðyÞ; implying thereby that the control force component h c ðxðtÞ; yðt ÞÞ h 0:

ð5:6Þ

Differentiating f1 ðx; yÞ and f2 ðx; yÞ with respect to time we get the relations f_ 1 ðx; yÞ dx_ 1 ðtÞKy_ 1 ðtÞZ 0 and f_ 2 ðx; yÞ dx_ 2 ðtÞKy_ 2 ðtÞZ 0. Since our theory requires that the given initial conditions satisfy the above-mentioned two Proc. R. Soc. A (2008)

2358 first component of Lorenz and Rossler system

F. E. Udwadia 30 20 10 0 −10 −20 −30

0

10

20

30

40

50

time Figure 1. The dynamics of the ﬁrst component of the Lorenz system x 1(t) and the Rossler system y1(t) are shown by the solid line and the thick dashed line, respectively. Fifty seconds of response are shown.

trajectory requirements, and they do not, we shall modify the trajectory requirements to f_ 1 ZKaf1 ;

f_ 2 ZKaf2 ;

ð5:7Þ

where aO0 is chosen to be a suitable parameter. We note that asymptotic solutions of equation (5.7) as t/N are fi (t)Z0, iZ1, 2 as required by (5.5). The parameter a in the numerical example is chosen to be 0.5. The trajectory description given by (5.7) and (5.6) then leads to (see equation (4.4)) 3 2 1 0 0 K1 0 0 7 6 3 2 Kaðx 1 K y1 Þ 6 0 1 0 0 K1 0 7 7 6 7 7 6 6 ð5:8Þ Ag Z 6 0 0 0 1 0 0 7 and bg Z 4 Kaðx 2 K y2 Þ 5: 7 6 60 0 0 0 1 07 hðyÞ 5 4 0 0 0 0 0 1 The weighting matrix N for the augmented six-dimensional dynamical system given in equation (5.4) will also need to be augmented so it is a diagonal 6 by 6 matrix. Owing to the trajectory requirement hc(t)h0, we can choose its last three diagonal entries to be each equal to unity; the ﬁrst three diagonal entries remain the same as those of the matrix NL. The explicit control input that causes the trajectory described by (5.6) and (5.7) to be tracked is then given by equation (4.14) with Bs Z Ag N K1=2 , since the 6 by 6 matrix Mg is an identity matrix. All the computations are performed using MATLAB and the integration is carried out using a relative error tolerance of 10K10 and an absolute error tolerance of 10K13. Figure 1 shows the dynamics of the ﬁrst component of the two separate dynamical systems described by equations (5.1) and (5.2). Proc. R. Soc. A (2008)

2359

Tracking control of nonlinear systems (a) 2.0

(b) tracking error (×10 –12)

1.5 tracking error

1.0 0.5 0 − 0.5 −1.0 −1.5 −2.0

0

10

20

30

40

50

6 4 2 0 −2 −4 −6 80

85

90 time

time

95

100

(a) 3 2 1 0 −1 −2 0

10

20

30 time

40

50

control inputs to Lorenz system (× 10 4)

control inputs to Rossler system (×10 −11)

Figure 2. (a) The solid line shows e1(t) over 50 s of integration. The dashed line shows e2(t). (b) Tracking errors e1(t) and e2(t) over the time interval [80–100] seconds using the same line conventions as in (a).

(b) 3 2 1 0 −1 −2 −3 −4 −5

0

10

20

30

40

50

time

Figure 3. (a) Components of the control inputs acting on the Rossler system. Note that the vertical scale ranges from approximately K3!10K11 to 3!10K11, making hcz0. The ﬁrst component of the control input is shown by a solid line, the second by a dashed line and the third by a thick solid line. (b) Components of the control inputs acting on the Lorenz system, using the same line notation.

On application of the optimal control input, the tracking errors e 1(t)d x 1(t)Ky1(t) and e 2(t)dx 1(t)Ky1(t) are shown in ﬁgure 2. These errors, as time increases, go down to the same order of magnitude as the tolerance used in the numerical integration of the differential equations as shown in ﬁgure 2b. Figure 3a shows the control inputs acting on the Rossler system. They are theoretically supposed to be zero, as required by the second trajectory requirement (5.6). They are seen to be very small, and of the same order of magnitude as the tolerances with which the integration is carried out. Figure 3b shows the control inputs required to be given to the Lorenz system so that it tracks the ﬁrst two components of the Rossler system. The third component of the control input to the Lorenz system is seen to be zero, since we require only the ﬁrst two components to be tracked. The optimal control cost J(t) is shown in ﬁgure 4. Proc. R. Soc. A (2008)

2360

F. E. Udwadia 10

8

6

4

2

0

10

20

30

40

50

time Figure 4. The optimal control cost JðtÞZ ½F c T N ðxÞF c as a function of time.

30 20 10 0 −10 −20 −30 −30

−20

−10

0

10

20

30

Figure 5. Projection of phase space trajectories over the time interval [0 –100] seconds of the controlled Lorenz system (dashed line) on the (x 1, x 2) plane and those of the Rossler system (solid line) on the (y1, y2) plane.

The projection of the phase trajectory of the controlled Lorenz system on to the (x 1, x 2) plane (solid line) and the projection of the corresponding phase trajectory of the Rossler system on to the (y1, y 2) plane (dashed line) are shown in ﬁgure 5. The ﬁgure shows the manner of convergence of these projected trajectories, which start with different initial conditions. Because the third component of the Rossler system is not tracked, the three-dimensional phase portrait looks very different for the controlled Lorenz system. The phase space portraits of these two systems are shown in ﬁgure 6. Proc. R. Soc. A (2008)

2361

Tracking control of nonlinear systems

(a) 50

(b) 50 40

0 −50

3

−100

30 20 10 0

20

20

0 −20

−20

0

20

0 −20

−20

0

20

Figure 6. The phase portraits of (a) the controlled Lorenz system and (b) the chaotic Rossler system over the time interval [0–100] seconds. The projections of these phase portraits on the horizontal plane are shown in ﬁgure 5.

6. Conclusions and remarks The methodology for the tracking control of nonlinear systems proposed herein has been inspired by results in analytical dynamics. This paper begins by developing this methodology for systems described by second-order differential equations, as commonly found in the Lagrangian and Newtonian mechanics, as well as ﬁrst-order differential equations, as found in the Hamiltonian and Poincare´ formulations of mechanics. It then extends the methodology to full state control of general nonlinear dynamical systems. The main contributions of the paper are the following: (i) The development of an explicit closed-form expression that provides the entire set of continuous tracking controllers that can exactly track a given trajectory description, assuming that the system’s initial conditions satisfy the description of the trajectory. We obtain explicit closed-form expressions for the controllers, which can be computed in real time. (ii) The development of a simple formula that explicitly gives the tracking control that minimizes the control cost J(t) at each instant of time. An explicit expression for the minimal cost is also obtained. (iii) For a general, ﬁrst-order, nonlinear system Mg ðx; tÞx_ Z f ðx; tÞ, x(0)Zx 0, the entire set of controllers needed to satisfy the trajectory described by the consistent equations fi(x, t)Z0, iZ1, ., m, is explicitly given by ð6:1Þ f c df1c C f 2c Z N K1=2 BsCðbg K Bs a s Þ C N K1=2 ðI KBsCB s Þz g ; where z g(x, t) is any LC function; N(x, t) is any positive-deﬁnite weighting matrix; and a s, bg and Bs are as deﬁned in equations (4.7)–(4.9). (iv) For a given weighting matrix N, the total control input can be split into two parts: a part that solves the optimal control problem that minimizes J(t) while exactly tracking the trajectory, and a second additive part that is N-orthogonal to the ﬁrst. While the addition of the second part to the optimal control effort allows the trajectory requirements to be still exactly satisﬁed, the norm of the control cost increases, in general, when it is added. The minimum cost J ðtÞZ kBsCðbg K Ag aÞk2 , and the corresponding controller that yields this minimum cost is f c df1c Z N K1=2 BsCðbg K Bs a s Þ. Proc. R. Soc. A (2008)

2362

F. E. Udwadia

The addition of the second part f2c Z N K1=2 ðI KBsCBs Þzg , so that f cd f1c C f2c , increases the cost beyond the optimal by kðI KBsCBs Þzg k2 . At each instant of time t, when z g(x, t) belongs to the range space of BsT , f2c ðt ÞZ 0, so that f c ðt Þ df1c ðtÞ, hence making the control at that instant of time optimal. (v) The close connection between nonlinear control and analytical mechanics is pointed out. Here we see that Nature seems to control a mechanical system so that it satisﬁes a given set of trajectory requirements—or alternatively stated, tracks a given trajectory—by using feedback control, much like a control engineer, except that instead of proportional, integral or derivative feedback control, which is commonly used by the control engineer, it uses acceleration feedback. While considerable work has been done on PID control, acceleration feedback seems to be far less studied by modern-day control engineers. Following Nature’s cue, the results developed in this paper point perhaps towards the need for more work in the area of acceleration feedback in mechanical systems. Furthermore, Nature seems to minimize the control cost J ðtÞZ ½Q c T N Q c at each instant of time, and it appears to use M K1 ðq; tÞ for the weighting matrix N. Nature’s choice of this weighting matrix can be understood when thought of in terms of a multibody mechanical system that is required to satisfy a given trajectory description (a set of constraints), and so track a given trajectory. Since it takes a larger control effort to move a body belonging to the multi-body system that has a larger inertia, Nature, in its effort to make the entire system satisfy the given trajectory description, appears to prefer applying control forces to bodies with the smaller inertias. Again, following Nature’s cue, the use of M K1 ðq; tÞ as a weighting matrix for deﬁning the control cost may be useful in many other dynamical systems, especially mechanical ones. (vi) While we have demonstrated the methodology by illustrating its use in determining the control required to be applied to a chaotic Lorenz system so that it tracks some components of the motion of another chaotic Rossler system, the general methodology can be used for more complex nonlinear mechanical systems dealing with, for example, orbital mechanics (Lam 2006). (vii) Finally, we note that the methodology presented herein does not include any magnitude constraints on the control, and work on this topic is currently being pursued. Furthermore, the results provided herein deal with full state control and further developments along the lines pursued herein to underactuated robotic systems would be useful. The effects of model errors and disturbance control also need to be investigated.

References Boccaletti, S., Kurths, J., Osipov, G., Valladares, D. & Zhou, C. 2002 The synchronization of chaotic oscillators. Phys. Rep. 366, 1–101. (doi:10.1016/S0370-1573(02)00137-0) Brockett, R. W. 1977 Control theory and analytical mechanics. In Geometric control theory (eds C. Martin & R. Herman), pp. 1–66. Brookline, MA: Math. Sci. Press. Brogliato, B., Lozano, R. & Maschke, B. 2007 Dissipative systems analysis and control. Berlin, Germany: Springer. Proc. R. Soc. A (2008)

Tracking control of nonlinear systems

2363

Chaou, Y. & Chang, W. 2004 Trajectory control. Lecture Notes on Control and Information Sciences. Berlin, Germany: Springer. Chen, H. K. 2002 Chaos and synchronization of a symmetric gryo with linear plus cubic damping. J. Sound Vib. 255, 719 –740. (doi:10.1006/jsvi.2001.4186) Gauss, C. F. 1829 Uber ein neues algemeines Grundgesetz der Mechanik. Zeitschrift fur die Reine und Angewandte Mathematik 4, 232–235. Goldstein, H. 1976 Classical mechanics. New York, NY: Addison-Wesley. Graybill, F. 2001 Matrices with application to statistics Duxbury classics series. Paciﬁc Grove, CA: Duxbury Publishing Company. Hamel, G. 1949 Theoretische mechanik. Berlin, Germany: Springer. Lagrange, J. L. 1811 Mechanique analytique. Paris, France: Mme Ve Coureier. Lam, T. 2006 New approach to mission design based on the fundamental equation of motion. J. Aerosp. Eng. 19, 59– 67. (doi:10.1061/(ASCE)0893-1321(2006)19:2(59)) Lei, Y., Xu, W. & Zheng, H. 2005 Synchronization of two chaotic gyros using active control. Phys. Lett. A 343, 153 –158. (doi:10.1016/j.physleta.2005.06.020) Moore, E. H. 1910 On the reciprocal of the general algebraic matrix. Abstr. Bull. Am. Soc. 26, 394– 395. Naidu, D. 2003 Optimal control systems. London, UK: CRC Press. Penrose, R. 1955 A generalized inverse of matrices. Proc. Camb. Philos. Soc. 51, 406 – 413. Poincare´, H. 1901 Sur une forme nouvelle des equations de la mechanique. CR Acad. Sci. Paris 132, 369 – 371. Pyras, K. 1996 Weak and strong synchronization of chaos. Phys. Rev. E 54, 4508–4511. (doi:10. 1103/PhysRevE.54.R4508) Rosenblum, M., Pilkovsky, A. & Kruths, J. 1997 From phase to lag synchronization in coupled chaotic oscillators. Phys. Rev. Lett. 78, 4193–4196. (doi:10.1103/PhysRevLett.78.4193) Rossler, O. 1976 An equation for continuous chaos. Phys. Lett. A 57, 397–398. (doi:10.1016/03759601(76)90101-8) Rulkov, N., Sushchik, M., Tsimring, L. & Abarbanel, H. 1995 Generalized synchronization of chaos in directionally coupled chaotic systems. Phys. Rev. E 51, 980–994. (doi:10.1103/PhysRevE.51.980) Sadegh, N. 1990 Stability analysis of a class of adaptive controllers for robotic manipulators. Int. J. Robot. 9, 74–92. (doi:10.1177/027836499000900305) Sastry, S. 1999 Nonlinear systems analysis, stability, and control. Berlin, Germany: Springer. Slotine, J. & Li, W. 1991 Applied nonlinear control. New York, NY: Englewood Cliffs. Strogatz, S. 1994 Nonlinear dynamics and chaos. Boulder, CO: Westview Press. Udwadia, F. E. 2000 Fundamental principles of lagrangian dynamics: mechanical systems with non-ideal, holonomic, and non-holonomic constraints. J. Math. Anal. Appl. 252, 341– 355. (doi:10.1006/jmaa.2000.7050) Udwadia, F. E. & Kalaba, R. E. 1992 A new perspective on constrained motion. Proc. R. Soc. A 439, 407– 410. (doi:10.1098/rspa.1992.0158) Udwadia, F. E. & Kalaba, R. E. 1996 Analytical dynamics: a new approach. Cambridge, UK: Cambridge University Press. Udwadia, F. E. & Kalaba, R. E. 2002 On the foundations of analytical dynamics. Int. J. Nonlin. Mech. 37, 1079 –1090. (doi:10.1016/S0020-7462(01)00033-6) Udwadia, F. E. & Phohomsiri, P. 2007 Explicit Poincare equations of motion for general constrained systems. Part I. Analytical results. Proc. R. Soc. A 1825, 1–14. (doi:10.1098/rspa. 2007.1825) Zaks, A., Park, E.-H., Rosenblum, M. & Kruths, J. 1999 Alternating locking ratios in imperfect phase synchronization. Phys. Rev. Lett. 82, 4228–4231. (doi:10.1103/PhysRevLett.82.4228)

Proc. R. Soc. A (2008)