Optimal Trajectory Tracking of Nonholonomic Mechanical Systems: a geometric approach

TL;DR

This paper presents a geometric optimal control approach for trajectory tracking in nonholonomic systems, minimizing position and velocity errors while respecting system constraints, demonstrated on a fully actuated particle example.

Contribution

It introduces a coordinate-free geometric framework and applies Pontryagin's Minimum Principle to derive optimal trajectories for nonholonomic systems.

Findings

Effective trajectory tracking achieved with the proposed method

Framework applicable to a class of nonlinear nonholonomic systems

Validated on a fully actuated particle example

Abstract

We study the tracking of a trajectory for a nonholonomic system by recasting the problem as an optimal control problem. The cost function is chosen to minimize the error in positions and velocities between the trajectory of a nonholonomic system and the desired reference trajectory evolving on the distribution which defines the nonholonomic constraints. We prepose a geometric framework since it describes the class of nonlinear systems under study in a coordinate-free framework. Necessary conditions for the existence of extrema are determined by the Pontryagin Minimum Principle. A nonholonomic fully actuated particle is used as a benchmark example to show how the proposed method is applied.

Equations

$$\frac{d}{dt}\left(\frac{\partial L}{\partial {\dot{q}}^i}\right)-\frac{\partial L}{\partial q^i}=0,1\le i\le n.$$

$${\Gamma }_{ij}^k={\left({\mathcal{G}}^{-1}\right)}_{ks}\left(\frac{\partial {\mathcal{G}}_{sj}}{\partial q^i}+\frac{\partial {\mathcal{G}}_{si}}{\partial q^j}+\frac{{\mathcal{G}}_{i}j}{\partial q^s}\right)$$

$$L(v_q)=\frac{1}{2}\mathcal{G}(v_q,v_q)-V(q),$$

$${\mathcal{D}}_x=span\{e_A(x)\},x\in U\subset Q.$$

$${\Gamma }_{AB}^C=\frac{1}{2}({\mathcal{C}}_{CA}^B+{\mathcal{C}}_{CB}^A+{\mathcal{C}}_{AB}^C)$$

$${\dot{q}}^i(t)={\rho }_A^i(q(t))v^A(t).$$

$$\ell (v)=\frac{1}{2}{\mathcal{G}}^{\mathcal{D}}(v,v)-V({\tau }_D(v)),withv\in \mathcal{D}.$$

$${\nabla }_{\gamma (t)}^{{\mathcal{G}}^{\mathcal{D}}}\gamma (t)+gra{d}_{{\mathcal{G}}^{\mathcal{D}}}V({\tau }_{\mathcal{D}}(\gamma (t)))=0.$$

$${\mathcal{G}}^{\mathcal{D}}(gra{d}_{{\mathcal{G}}^{\mathcal{D}}}V,X)=X(V),\text{mbox}foreveryX\in \Gamma ({\tau }_{\mathcal{D}}).$$

$${\dot{q}}^i$$

$${\dot{v}}^C$$

$${\dot{v}}_x=$$

$${\dot{v}}_z=$$

$${\dot{v}}_x=$$

$$v_z=$$

$$Y_1(q)$$

$$Y_2(q)$$

$$e=v^1\frac{\partial }{\partial y}+v^2\left(\frac{\partial }{\partial z}-y\frac{\partial }{\partial x}\right)$$

$$[[Y_1,Y_2]]=\mathcal{P}[Y_1,Y_2]=\mathcal{P}\left(-\frac{\partial }{\partial x}\right)=\frac{y}{1+y^2}\left(\frac{\partial }{\partial z}-y\frac{\partial }{\partial x}\right).$$

$$\ell (x,y,z,y_1,y_2)=\frac{1}{2}\left((v^1{)}^2+(v^2{)}^2(y^2+1)\right).$$

$${\dot{v}}^1=0,{\dot{v}}^2=-\frac{y}{1+y^2}{v}^1{v}^2$$

$$x(t)=$$

$$z(t)=$$

$$-z_0{c}_2,$$

$$v^1(t)=$$

$${\nabla }_{\gamma (t)}^{{\mathcal{G}}^{\mathcal{D}}}\gamma (t)+gra{d}_{{\mathcal{G}}^{\mathcal{D}}}V({\tau }_{\mathcal{D}}(\gamma (t)))\in \Gamma ({\tau }_D),$$

$${\nabla }_{\gamma (t)}^{{\mathcal{G}}^{\mathcal{D}}}\gamma (t)+gra{d}_{{\mathcal{G}}^{\mathcal{D}}}V({\tau }_{\mathcal{D}}(\gamma (t)))=u^A(t)e_A({\tau }_{\mathcal{D}}(\gamma (t))),$$

$${\dot{q}}^i$$

$${\dot{v}}^A$$

$$\mathcal{J}(\gamma (t),u(t)):={\int }_0^T\mathcal{C}(\gamma (t),u(t))dt.$$

$$\mathcal{J}(\gamma (t))=$$

$$=\frac{1}{2}{\int }_0^T(||q^i(t)-q_r^i(t)|{|}^2+||v^A(t)-v_r^A(t)|{|}^2$$

$$+ϵ||u^A|{|}^2)dt+\omega \Phi (\gamma (T))$$

$$\mathcal{H}(q,v,\lambda ,\mu ,u)=$$

$$+{\mu }_A{\dot{v}}^A(q^i,v^A,u^A)$$

$${\mathcal{H}}_{opt}(q,v,\lambda ,\mu ,u^{\star })=\underset{u\in \mathcal{U}}{\min }\mathcal{H}(q,v,\lambda ,\mu ,u),$$

Full text

I Introduction

Nonholonomic optimal control problems arise in many engineering applications, for instance systems with wheels, such as cars and bicycles, and systems with blades or skates. There are thus multiple applications in the context of wheeled motion, space or mobile robotics and robotic manipulation. The earliest work on control of nonholonomic systems is by R. W. Brockett in [9]. A. M. Bloch [2], [3] has examined several control theoretic issues which pertain to both holonomic and nonholonomic systems in a very general form. The seminal works about stabilization in nonholonomic control systems were done by A. M. Bloch, N. H. McClamroch, and M. Reyhanoglu in [3], [5], [6], [7], and more recently by A. Zuryev [27].

A geometrical dynamical system of mechanical type is completely determined by a Riemannian manifold $Q$, a kinetic energy, which is defined through the Riemannian metric $\mathcal{G}$ on the manifold and the potential forces encoded into a potential (conservative) function $V:Q\to\mathbb{R}$. These objects, together with a non-integrable distribution $\mathcal{D}\subset TQ$ on the tangent bundle of the configuration space determines a nonholonomic mechanical system.

Stabilization of an equilibrium point of a mechanical system on a Riemannian manifold has been a problem well studied in the literature from a geometric framework along the last decades (see [2] and [11] for a review on the topic). Further extensions of these results to the problem of tracking a smooth and bounded trajectory can be found in [11] where a proportional and derivative plus feed forward (PD+FF) feedback control law is proposed for tracking a trajectory on a Riemannian manifold using error functions.

For trajectory tracking, the usual approach of stabilization of error dynamics [19], [22], [23], [25] cannot be utilized for nonholonomic systems. This is because there does not exist a $\mathcal{C}^{1}$ (even continuous) state feedback which can asymptotically stabilize the trajectory of a nonholonomic system about a desired equilibrium point. The closed loop trajectory violates Brockett’s condition [10], [7] which states that any system of the form $\dot{x}=f(x,u)$ must have a neighborhood of zero in the image of the map $x\to f(x,u)$ for some $u$ in the control set. This result appears in Theorem 4 in [7].

In this paper, we introduce a geometrical framework in nonholonomic mechanics to study tracking of trajectories for nonholonomic systems based on [12], [16], [17]. The application of modern tools from differential geometry in the fields of mechanics, control theory and numerical integration has led to significant progress in these research areas. For instance, the study in a geometrical formulation of the nonholonomic equations of motion has led to better understanding of different engineering problems such locomotion generation, controllability, motion planning, and trajectory tracking.

Combining the ideas of geometric methods in control theory, nonholonomic systems and optimization techniques, in this paper, we study the underlying geometry of a tracking problem for nonholonomic systems by understanding it as an optimal control problems for mechanical systems subject to nonholonomic constraints.

Given a reference trajectory $\gamma_{r}(t)=(q_{r}(t),v_{r}(t))$ on $\mathcal{D}$ the problem studied in this work consists on finding an admissible curve $\gamma(t)\in\mathcal{D}$, solving a dynamical control system, with prescribed boundary conditions on $\mathcal{D}$ and minimizing a cost functional which involves the error between the reference trajectory and the trajectory we want to find (in terms of both, positions and velocities), and the effort of the control inputs. This cost functional is accomplished with a weighted terminal cost (also known as Mayer term) which induces a constraint into the dynamics on $\mathcal{D}$. The interval length for the cost functional $T$ may either be fixed, or appear as a degree of freedom in the optimization problem, or be time horizon. In this work, we restrict to the case when $T$ is fixed.

To test the efficiency of the proposed method, we use a Runge Kutta integrator together with a shooting method in the solution of a trajectory optimization for a simple but challenging benchmark mechanical system: a fully actuated particle subject to a nonholonomic constraint into the dynamics.

We propose a geometric derivation of the equations of motion for tracking a trajectory of a nonholonomic system as an optimal control problem find we find necessary conditions via the Pontryagin Minimum Principle (PMP), where the optimal Hamiltonian is defined on the cotangent bundle of the constraint distribution. This approach allow for the reduction in the degrees of freedom of the equations for the optimal control problem, compared with typical methods describing the dynamics of a nonholonomic system, as the ones arising from the application of Lagrange-d’Alembert principle. The main advantages in this geometric framework consist in the use of a basis of vector fields on $\mathcal{D}$ allowing the reduction of some degrees of freedom in the dynamics for a nonholonomic mechanical system. The paper is structured as follows: we introduce mechanical systems on a manifold, connections on a Riemannian manifold and the geometry of nonholonomic dynamical systems on Section II, together with the example we used as benchmark the nonholonomic particle. Section III introduces the details of the problem under study motivated by the non-existence of a $\mathcal{C}^{1}$ feedback control to asymptotically stabilize the error dynamics in nonholonomic systems. Necessary conditions for the existence of extrema in the proposed optimal control problem are studied from the PMP in Section IV. We also show numerical results and analyze the results we obtain.

II Nonholonomic Mechanical Systems

II-A Preliminaries

Let $Q$ be a $n$-dimensional differentiable manifold with local coordinates $(q^{i})$, with $1\leq i\leq n$, the configuration space of a mechanical system. Denote by $TQ$ its tangent bundle with induced local coordinates $(q^{i},\dot{q}^{i})$. Given a Lagrangian function $L:TQ\rightarrow\mathbb{R}$, its Euler-Lagrange equations are

[TABLE]

These equations determine a system of implicit second-order differential equations in general. If we assume that the Lagrangian is regular, that is, the ${n\times n}$ matrix $\left(\frac{\partial^{2}L}{\partial\dot{q}^{i}\partial\dot{q}^{j}}\right)$ is non-degenerate, the local existence and uniqueness of solutions is guaranteed for any given initial condition.

Vector fields are used to calculate the directional derivative of a function defined on $Q$. In the realm of differential geometry a more general operator is defined to perform derivation of a wider range of geometric objects (tensors). This operator is called connection (linear, covariant, or affine connection). The definition of the connection is a wish list of properties which it is expected to have

Definition II.1

An (affine) connection on a smooth manifold $Q$ is a map which takes a pair consisting of a vector (or a vector field), and a $(p,q)$-tensor field, $T$, and returns a $(p,q)$-tensor field, such that it satisfies the following axioms

• •

$\nabla_{X}f=X(f)$, for $f\in\mathcal{C}^{\infty}(Q)$,

• •

$\nabla_{X}(T+S)=\nabla_{X}T+\nabla_{X}S$, for $T$ and $S$ tensors of the same type,

• •

$\nabla_{X}T(f,g)=(\nabla_{X}T)(f,g)+T(\nabla_{X}f,g)+T(f,\nabla_{X}g)$.

This definition of a connection is complete, i.e., this list of properties results in a uniquely defined geometric operator; however, an extra structure on the manifold is needed to define this object in a chart. To do so, we need to know how it acts on the basis of the tangent vector space. The result is a tangent vector field, and at each point it is spanned by the basis of the tangent space at that point $\displaystyle{\nabla_{\frac{\partial}{\partial q^{i}}}\left(\frac{\partial}{\partial q^{j}}\right)=\Gamma^{k}_{ij}\frac{\partial}{\partial q^{k}}.}$

Denote by $\mathfrak{X}(Q)$ the set of vector fields on $Q$. A metric $\mathcal{G}$ on a smooth manifold is a $(0,2)$-tensor field satisfying

• •

Symmetry: $\mathcal{G}(X,Y)=\mathcal{G}(Y,X)$ $X,Y\in\mathfrak{X}(Q)$,

• •

Non-degeneracy: $\mathcal{G}(X,Y)=0$ if and only if when $X=0$ then $Y=0$.

Locally, the metric is determined by the matrix $M=(\mathcal{G}_{ij})_{1\leq i,j\leq n}$ where $\mathcal{G}_{ij}=\mathcal{G}(\partial/\partial q^{i},\partial/\partial q^{j})$.

Using the metric $\mathcal{G}$ we may compute the Christoffel symbols associated with the metric as

[TABLE]

where $\mathcal{G}^{-1}$ is defined as the inverse of the metric with components determined by the inverse matrix of $M$.

II-B Nonholonomic mechanical systems

Most nonholonomic systems have linear constraints, and these are the ones we will consider. Linear constraints on the velocities (or Pfaffian constraints) are locally given by equations of the form $\phi^{a}(q^{i},\dot{q}^{i})=\mu^{a}_{i}(q)\dot{q}^{i}=0,\quad 1\leq a\leq m$, depending, in general, on their configuration coordinates and their velocities. From an intrinsic point of view, the linear constraints are defined by a regular distribution ${\mathcal{D}}$ on $Q$ of constant rank $n-m$ such that the annihilator of ${\mathcal{D}}$ is locally given at each point of $Q$ by ${\mathcal{D}}^{o}_{q}=\operatorname{span}\left\{\mu^{a}(q)=\mu_{i}^{a}dq^{i}\;;1\leq a\leq m\right\}$ where the one-forms $\mu^{a}$ are independent at each point of $Q$.

Now we restrict ourselves to the case of nonholonomic mechanical systems where the Lagrangian is of mechanical type, that is, a Lagrangian systems $L:TQ\to\mathbb{R}$ defined by

[TABLE]

with $v_{q}\in T_{q}Q$, where $\mathcal{G}$ denotes a Riemannian metric on the configuration space $Q$ representing the kinetic energy of the systems and $V:Q\rightarrow\mathbb{R}$ is a potential function.

Next, assume that the system is subject to nonholonomic constraints, defined by a regular distribution $\mathcal{D}$ on $Q$ with corank$(\mathcal{D})=m$. Denote by $\tau_{\mathcal{D}}:\mathcal{D}\rightarrow Q$ the canonical projection of $\mathcal{D}$ onto $Q$ and by $\Gamma(\tau_{\mathcal{D}})$ the set of sections of $\tau_{D}$ which in this case is just the set of vector fields $\mathfrak{X}(Q)$ taking values on $\mathcal{D}.$ If $X,Y\in\mathfrak{X}(Q),$ then $[X,Y]$ denotes the standard Lie bracket of vector fields.

Definition II.2

A nonholonomic mechanical system on a smooth manifold $Q$ is given by the triple $(\mathcal{G},V,\mathcal{D})$, where $\mathcal{G}$ is a Riemannian metric on $Q,$ representing the kinetic energy of the system, $V:Q\rightarrow\mathbb{R}$ is a smooth function representing the potential energy and $\mathcal{D}$ a non-integrable regular distribution on $Q$ representing the nonholonomic constraints.

Given $X,Y\in\Gamma(\tau_{\mathcal{D}})$ that is, $X(x)\in\mathcal{D}_{x}$ and $Y(x)\in\mathcal{D}_{x}$ for all $x\in Q,$ then it may happen that $[X,Y]\notin\Gamma(\tau_{\mathcal{D}})$ since $\mathcal{D}$ is nonintegrable. We want to obtain a bracket definition for sections of $\mathcal{D}.$ Using the Riemannian metric $\mathcal{G}$ we can define two complementary orthogonal projectors ${\mathcal{P}}\colon TQ\to{\mathcal{D}}$ and ${\mathcal{Q}}\colon TQ\to{\mathcal{D}}^{\perp},$ with respect to the tangent bundle orthogonal decomposition $\mathcal{D}\oplus\mathcal{D}^{\perp}=TQ$. Therefore, given $X,Y\in\Gamma(\tau_{\mathcal{D}})$ we define the nonholonomic bracket $[\![\cdot,\cdot]\!]:\Gamma(\tau_{\mathcal{D}})\times\Gamma(\tau_{\mathcal{D}})\rightarrow\Gamma(\tau_{\mathcal{D}})$ as $[\![X_{A},X_{B}]\!]:=\mathcal{P}[X_{A},X_{B}]$. This Lie bracket verifies the usual properties of a Lie bracket except the Jacobi identity (see [4], [14] for example).

Definition II.3

Consider the restriction of the Riemannian metric $\mathcal{G}$ to the distribution $\mathcal{D}$, $\mathcal{G}^{\mathcal{D}}:\mathcal{D}\times_{Q}\mathcal{D}\rightarrow\mathbb{R}$ and define the Levi-Civita connection $\displaystyle{\nabla^{\mathcal{G}^{\mathcal{D}}}:\Gamma(\tau_{\mathcal{D}})\times\Gamma(\tau_{\mathcal{D}})\rightarrow\Gamma(\tau_{\mathcal{D}})}$ determined by the following two properties:

$[\![X,Y]\!]=\nabla_{X}^{\mathcal{G}^{\mathcal{D}}}Y-\nabla_{Y}^{\mathcal{G}^{\mathcal{D}}}X,$** 2. 2.

$X(\mathcal{G}^{\mathcal{D}}(Y,Z))=\mathcal{G}^{\mathcal{D}}(\nabla_{X}^{\mathcal{G}^{\mathcal{D}}}Y,Z)+\mathcal{G}^{\mathcal{D}}(Y,\nabla_{X}^{\mathcal{G}^{\mathcal{D}}}Z).$**

Let $(q^{i})$ be local coordinates on $Q$ and $\{e_{A}\}$ be independent vector fields on $\Gamma(\tau_{D})$ (that is, $e_{A}(x)\in{\mathcal{D}}_{x}$) such that

[TABLE]

Then, we can determine the Christoffel symbols $\Gamma_{BC}^{A}$ of the connection $\nabla^{\mathcal{G}^{\mathcal{D}}}$ by $\displaystyle{\nabla_{e_{B}}^{\mathcal{G}^{\mathcal{D}}}e_{C}=\Gamma_{BC}^{A}(q)e_{A}.}$

As when we work in tangent bundles, it is possible to determine the Christoffel symbols associated with the connection $\nabla^{{\mathcal{G}}^{\mathcal{D}}}$ by $\nabla^{{\mathcal{G}}^{\mathcal{D}}}_{e_{B}}{e_{C}}=\Gamma^{A}_{BC}e_{A}$. Note that the coefficients $\Gamma_{AB}^{C}$ of the connection $\nabla^{{\mathcal{G}}^{\mathcal{D}}}$ are (see [1] for details)

[TABLE]

where the constant structures $\mathcal{C}_{AB}^{C}$ are defined as $[\![X_{A},X_{B}]\!]=\mathcal{C}_{AB}^{C}X_{C}$.

Definition II.4

A curve $\gamma:I\subset\mathbb{R}\rightarrow\mathcal{D}$ is admissible if $\displaystyle{\gamma(t)=\frac{d\sigma}{dt}(t)}$, where $\tau_{\mathcal{D}}\circ\gamma=\sigma$.

Given local coordinates on $Q,$ $(q^{i})$ with $i=1,\ldots,n;$ and $\{e_{A}\}$ sections on $\Gamma(\tau_{\mathcal{D}})$, with $A=1,\ldots,n-m$, such that $\displaystyle{e_{A}=\rho_{A}^{i}(q)\frac{\partial}{\partial q^{i}}}$ we introduce induced coordinates $(q^{i},v^{A})$ on $\mathcal{D}$, where, if $e\in\mathcal{D}_{x}$ then $e=v^{A}e_{A}(x).$ Therefore, $\gamma(t)=(q^{i}(t),v^{A}(t))$ is admissible if

[TABLE]

Consider the restricted Lagrangian function $\ell:\mathcal{D}\rightarrow\mathbb{R},$

[TABLE]

Definition II.5

A solution of the nonholonomic problem is an admissible curve $\gamma:I\rightarrow\mathcal{D}$ such that

[TABLE]

Here the section $grad_{{\mathcal{G}}^{\mathcal{D}}}V\in\Gamma(\tau_{\mathcal{D}})$ is characterized by

[TABLE]

These equations are equivalent to the nonholonomic equations. Locally, these equations are given by

[TABLE]

where $(\mathcal{G}^{\mathcal{D}})^{AB}$ denotes the coefficients of the inverse matrix of $(\mathcal{G}^{\mathcal{D}})_{AB}$ where $\mathcal{G}^{\mathcal{D}}(e_{A},e_{B})=(\mathcal{G}^{\mathcal{D}})_{AB}.$

Remark II.6

The nonholonomic equations only depend on the coordinates $(q^{i},v^{A})$ on $\mathcal{D}$. Therefore the nonholonomic equations are free of Lagrange multipliers. These equations are equivalent to the nonholonomic Hamel equations (see [8], for example, and references therein).

II-C Example: The nonholonomic particle

Consider a particle of unit mass evolving in $Q=\mathbb{R}^{3}$ with Lagrangian $\displaystyle{L(x,y,z,\dot{x},\dot{y},\dot{z})=\frac{1}{2}(\dot{x}^{2}+\dot{y}^{2}+\dot{z}^{2})}$, and subject to the constraint $\dot{x}+y\,\dot{z}=0$.

This nonholonomic system is defined by the annhilation of the one-form $\mu(x,y,z)=(1,0,y)$. The nonholonomic equations, derived from the Lagrange-d’Alembert principle, are given by

[TABLE]

which, after substituting the Lagrange multiplier $\displaystyle{\lambda=-\frac{v_{z}v_{y}}{1+y^{2}}}$, lead to

[TABLE]

such that $v_{x}+y\,v_{z}=0$. Let $\mathcal{D}\subset T\mathbb{R}^{3}\simeq\mathbb{R}^{3}\times\mathbb{R}^{3}$ denote the nonholonomic distribution corresponding to this system. Then these equations define a time-continuous flow $F_{t}:\mathcal{D}\rightarrow\mathcal{D}$, i.e. $F_{t}((q(0),v(0)))=(q(t),v(t))$, where $q(t)=(x(t),y(t),z(t))^{T}$ and $v(t)=(v_{x}(t),v_{y}(t),v_{z}(t))^{T}$, $(q(0),v(0))\in\mathcal{D}$.

The distribution $\mathcal{D}$ is determined by \mathcal{D}=\hbox{span}\{Y_{1},Y_{2}\}=\hbox{span}\Big{\{}\frac{\partial}{\partial y},\frac{\partial}{\partial z}-y\frac{\partial}{\partial x}\Big{\}}. Then, $\mathcal{D}^{\perp}=\{\frac{\partial}{\partial x}+y\frac{\partial}{\partial z}\}$. Let $(x,y,z,v^{1},v^{2})$ be induced coordinates on $\mathcal{D}$.

Given the vector fields $Y_{1}$ and $Y_{2}$ generating the distribution we obtain the relations for $q\in\mathbb{R}^{3}$

[TABLE]

Then, $\rho_{1}^{1}=\rho_{1}^{3}=\rho_{2}^{2}=0,\quad\rho_{1}^{2}=\rho_{2}^{3}=1,\quad\rho_{2}^{1}=-y$.

Each element $e\in{\mathcal{D}}_{q}$ is expressed as a linear combination of these vector fields: $e=v^{1}Y_{1}(q)+v^{2}Y_{2}(q),\quad q\in\mathbb{R}^{3}$. Therefore, the vector subbundle $\tau_{\mathcal{D}}:{\mathcal{D}}\rightarrow\mathbb{R}^{3}$ is locally described by the coordinates $(x,y,\theta;v^{1},v^{2})$; the first three for the base and the last two, for the fibers. Observe that

[TABLE]

and, in consequence, ${\mathcal{D}}$ is described by the conditions (admissibility conditions): $\dot{x}=-yv^{2},\quad\dot{y}=v^{1},\quad\dot{z}=v^{2}$ as a vector subbundle of $TQ$ where $v^{1}$ and $v^{2}$ are the adapted velocities relative to the basis of $\mathcal{D}$ defined before.

The nonholonomic bracket given by $[\![\cdot,\cdot]\!]=\mathcal{P}([\cdot,\cdot])$ satisfies

[TABLE]

Therefore, by using (2) all the Christoffel symbols for the connection $\displaystyle{\nabla^{\mathcal{G}^{\mathcal{D}}}}$ vanish except $\Gamma_{12}^{2}$ which is given by $\displaystyle{\Gamma_{12}^{2}=\frac{y}{1+y^{2}}}$. The restriction of the Lagrangian function $L$ on $\mathcal{D}$ in the adapted coordinates $(v^{1},v^{2})$ is given by

[TABLE]

Then, the nonholonomic equations for the constrained particle are given by

[TABLE]

together with the admissibility conditions $\dot{x}=-yv^{2}$, $\dot{y}=v^{1}$ and $\dot{z}=v^{2}.$ Then these equations define a time-continuous flow $F_{t}^{\mathcal{D}}:\mathcal{D}\rightarrow\mathcal{D}$, i.e. $F_{t}((q(0),v(0)))=(q(t),v(t))$, where $q(t)=(x(t),y(t),z(t))^{T}$ and $v(t)=(v_{1}(t),v_{2}(t))^{T}$, $(q(0),v(0))\in\mathcal{D}$.

The previous systems can be integrated explicitly, and solutions are given by:

[TABLE]

for $x_{0},y_{0},z_{0},c_{1},c_{2}$ constants to be determined by the initial conditions.

Remark II.7

Note that previous equations have a singularity at $c_{1}=0$. The constant $c_{1}$ arrises from the equation for $\dot{v}^{1}$. If $c_{1}=0$, and therefore ${v}^{1}(t)=0$, then the solution for the system of equations is given by $x(t)=-y_{0}v_{0}^{2}t+x_{0}$, $y(t)=y_{0}$, $z(t)=v_{0}^{2}t+z_{0}$, $v^{2}(t)=v_{0}^{2}$, where $x_{0},y_{0},z_{0},v_{0}^{2}$ are constants. $\diamond$

III Optimal trajectory tracking problem

The purpose of this section is to present the tracking problem for nonholonomic systems as an optimal control problem. The objective is the tracking of a suitable reference trajectory $\Upsilon(t)$ for a mechanical system with nonholonomic velocity constraints as described in the previous section. It is assumed that $\Upsilon(t)\in\mathcal{D}$.

We will analyze the case when the dimension of the inputs set or control distribution is equal to the rank of $\mathcal{D}$. If the rank of $\mathcal{D}$ is equal to the dimension of the control distribution, the system will be called a fully actuated nonholonomic system.

Definition III.1

A solution of a fully actuated nonholonomic problem is an admissible curve $\gamma:I\rightarrow\mathcal{D}$ such that

[TABLE]

or, equivalently,

[TABLE]

where $u^{A}$ are the control inputs.

Locally, the above equations are given by

[TABLE]

As we mentioned in the Introduction, For trajectory tracking, the usual approach of stabilization of error dynamics [19], [22], [23], [25] cannot be utilized for nonholonomic systems because the closed loop trajectory violates Brockett’s condition. A common approach to trajectory tracking for nonholonomic systems found in the literature is the backstepping procedure [15], [18]. This approach is done on a per example basis, in particular, mobile robots or unycicle models. In [15], [18] the error dynamics of the unicycle model is shown to be in strict feedback form. Thereafter, integrator backstepping is employed to choose an appropriate Lyapunov function for stabilization of the error dynamics. This error dynamics does not evolve on the constrained manifold (unlike our approach). Therefore, Brockett’s condition is not violated. However, since $\rho^{i}_{A}(q)$ is unknown in a general framework, the approach can not be generalized to solve the tracking problem for a general nonholonomic system with our method and then backstepping needs to be studied for each system.

So we propose a new approach to consider tracking problem as an optimal control problem and we call this optimal tracking.

In the following, we shall assume that all the control systems under consideration are controllable in the configuration space, that is, for any two points $q_{0}$ and $q_{f}$ in the configuration space $Q$, there exists an admissible control $u(t)$ defined on the control manifold $\mathcal{U}\subseteq\mathbb{R}^{n}$ such that the system with initial condition $q_{0}$ reaches the point $q_{f}$ at time $T$ (see [2] for more details). Given a cost function $\mathcal{C}:\mathcal{D}\times\mathcal{U}\rightarrow\mathbb{R}$ the optimal control problem consists of finding an admissible curve $\gamma:I\rightarrow\mathcal{D}$ which is a solution of the fully actuated nonholonomic problem given initial and final boundary conditions on $\mathcal{D}$ and minimizing the cost functional

[TABLE]

For trajectory tracking of a nonholonomic system we consider the following problem

Problem (optimal trajectory tracking): Given a reference trajectory $\gamma_{r}(t)=(q_{r}(t),v_{r}(t))$ on $\mathcal{D}$, find an admissible curve $\gamma(t)\in\mathcal{D}$, solving (10)-(11), with prescribed boundary conditions on $\mathcal{D}$ and minimizing the cost functional

[TABLE]

where $\epsilon>0$ is a regularization parameter, $\Phi:TQ\to\mathbb{R}$ is a terminal cost (Mayer term), $\omega>0$ is a weight for the terminal cost. $\mathcal{C}$ and $\Phi$ are assumed to be continuously differentiable functions, and the final state $\gamma(T)$ is required to fulfill a constraint $r(\gamma(T),\gamma_{r})=0$ with $r:\mathcal{D}\times\mathcal{D}\to\mathbb{R}^{d}$ and $\gamma_{r}\in\mathcal{D}$ given. The interval length $T$ may either be fixed, or appear as degree of freedom in the optimization problem. In this work we restrict to the case when $T$ is fixed.

Remark III.2

Note that if $\epsilon=0$ then the optimal control problem turns into a singular optimal control problem (see [21] Section $3.2$) $\diamond$.

IV Necessary conditions for optimality

In this section we apply Pontryagin’s minimization principle to the optimal tracking problem. The Hamiltonian $\mathcal{H}:T^{*}\mathcal{D}\times\mathcal{U}\to\mathbb{R}$ for the problem is given by

[TABLE]

where $\dot{v}^{A}$ comes from equation (11). Note that $\lambda_{i}$ and $\mu_{A}$ are the costate variables or Lagrange multipliers. The last two terms in (12) corresponds with the nonholonomic dynamics given in equations (3) and (4) paired with the costate variables, which represents the standard construction of the Hamiltonian for the PMP. Also note that $\mathcal{H}$ is defined on a subset of $T^{*}(TQ)$.

Denote by $t\mapsto u^{\star}(t)$ a curve that satisfies along a trajectory $t\mapsto(q(t),v(t),\lambda(t),\mu(t))\in T^{*}\mathcal{D}$,

[TABLE]

then $u^{*}$ may be determined implicitly as a function of $(q(t),v(t),\lambda(t),\mu(t))\in T^{*}\mathcal{D}$ using the previous equation and then we may define the optimal Hamiltonian $\mathcal{H}_{opt}:T^{*}\mathcal{D}\to\mathbb{R}$ by prescribing the control $u$ as $u^{\star}$.

Given that $u^{\star}$ minimizes $\mathcal{H}$, then $u^{\star}$ is a critical point for $\mathcal{H}$ and may be uniquely determined by

[TABLE]

The PMP applied to our particular problem gives the following necessary conditions

• •

Stationary condition: from equation (13) $\mu_{A}=-\epsilon u^{A},$

• •

State equation: Equations (10) and (11),

• •

Adjoint equations (or costate equations):

[TABLE]

• •

Constraint induced by terminal cost: $r(\gamma(T),\gamma_{r})=0$,

• •

Boundary conditions: $\gamma(0):=(q(0),v(0))\in\mathcal{D}$, $\frac{\partial\Phi}{\partial q^{i}}(\gamma(T))=\lambda_{i}(T)$, $\frac{\partial\Phi}{\partial v^{A}}(\gamma(T))=\mu_{A}(T)$.

IV-A Optimal trajectory tracking for the nonholonomic particle

Consider the situation of Example II-C. Let $\gamma_{r}$ be the reference trajectory, $\gamma_{r}=(x_{r}(t),y_{r}(t),z_{r}(t),v_{1,r},v_{2,r})$ which follows the constraint $\dot{x}_{r}=y_{r}\dot{z}_{r}$ at all $t$ and the dynamical equations for the nonholonomic particle. We wish to control the velocity of the nonholonomic particle. We add then control inputs in the fiber coordinates $v^{1}$ and $v^{2}$. Therefore the control dynamical system to study is given by

[TABLE]

together with the admissibility conditions $\dot{x}=-yv^{2}$, $\dot{y}=v^{1}$ and $\dot{z}=v^{2}.$

The cost function $\mathcal{C}:\mathcal{D}\times\mathcal{U}\to\mathbb{R}$ for the optimal control problem is given by

[TABLE]

and the terminal cost function is given by

[TABLE]

with $T\in\mathbb{R}^{+}$ fixed.

The Hamiltonian for the PMP is given as

[TABLE]

In order for $u(t)$ to be the optimal control we employ the stationary condition. Therefore, $\displaystyle{u_{1}^{\star}=-\frac{\mu_{1}}{\epsilon}\;\hbox{ and }\;u_{2}^{\star}=-\frac{\mu_{2}}{\epsilon}}$. The final cost is given by $\Phi(\gamma(T))=||\gamma(T)-\gamma_{r}||^{2}$ which induces the constraint

[TABLE]

Finally, the optimal Hamiltonian $\mathcal{H}_{opt}$ is given by

[TABLE]

The adjoint equations are $\dot{\lambda}_{1}=-(x-x_{r})$, $\dot{\lambda}_{3}=-(z-z_{r})$,

[TABLE]

The state equations were given in Example II-C in equation (8) together with the admissibility conditions. Boundary conditions must satisfy the constraints in order for the trajectory to evolve on $\mathcal{D}$, that is $\dot{x}_{0}+y_{0}\,\dot{z}_{0}=0$ where $\dot{x}_{0},y_{0},\dot{z}_{0}$ denotes the boundary conditions for the variables $\dot{x}$, $y$ and $\dot{z}$ respectively.

IV-B Numerical results

We now test with numerical simulations how the proposed methods work.

Denote $F_{\mu}^{\lambda}:[0,T]\times T^{*}\mathcal{D}\to T^{*}\mathcal{D}$, the integral flow given by equations (15) on $T^{*}\mathcal{D}$ and $\gamma(0)\in\mathcal{D}$ the initial condition for the state dynamics. The initial guess for the initial condition of the costate variables is denoted by $\alpha=F_{\mu}^{\lambda}(0)$. We wish to find the initial condition of the costates for which $F_{\mu}^{\lambda}(T,\gamma(0),\alpha)=(0_{1\times 5})^{T}$ . The goal is to find the root of the polynomial

[TABLE]

where $T\in\mathbb{R}^{+}$ is the final time, $\omega\in\mathbb{R}^{+}$ is a weight for the terminal cost and $F_{\mu}^{\lambda}(\tau,\gamma(0),p_{0})$ is the flow of the adjoint equations (15) starting at $(\gamma(0),p_{0})$. The root finder used in both situations was the fsolve routine in MATLAB.

For the intial condition $\gamma(0)=$ $\begin{pmatrix}0.5&0.2&0.7;&0.5&0.4\end{pmatrix}$ and reference trajectory $\gamma_{r}(t)=$ $\left(1,0,t+1,0,1\right)$, $p_{0}=0_{1\times 5}$, $T=4$ and $\epsilon=7$ we exhibit the results in Figure 1.

Acknowledgments: The authors wish to thank Prof. Ravi Banavar for fruitful discussions along the development of this work .