---
abstract: 'The paper studies the problem of making Getz’s bicycle model traverse a strictly convex Jordan curve with bounded roll angle and bounded speed. The approach to solving this problem is based on the virtual holonomic constraint (VHC) method. Specifically, a VHC is enforced making the roll angle of the bicycle become a function of the bicycle’s position along the curve. It is shown that the VHC can be automatically generated as a periodic solution of a scalar periodic differential equation, which we call virtual constraint generator. Finally, it is shown that if the curve is sufficiently long as compared to the height of the bicycle’s centre of mass and its wheel base, then the enforcement of a suitable VHC makes the bicycle traverse the curve with a steady-state speed profile which is periodic and independent of initial conditions. An outcome of this work is a proof that the constrained dynamics of a Lagrangian control system subject to a VHC are generally not Lagrangian.'
author:
- 'Luca Consolini, Manfredi Maggiore,[^1]'
bibliography:
- 'citazioni.bib'
title: Control of a Bicycle Using Virtual Holonomic Constraints
---

Introduction
============

This paper investigates the problem of maneuvering a bicycle along a closed Jordan curve $\C$ in the horizontal plane in such a way that the bicycle does not fall over and its velocity is bounded. The simplified bicycle model we use in this paper, developed by Neil Getz [@Get94; @GetMar95], views the bicycle as a point mass with a side slip velocity constraint, and models its roll dynamics as those of an inverted pendulum, see Figure \[fig:bicycle\]. The model neglects, among other things, the steering kinematics and the wheels dynamics with the associated gyroscopic effect. In [@HauSacFre04], Hauser-Saccon-Frezza investigate the maneuvering problem for Getz’s bicycle using a dynamic inversion approach to determine bounded roll trajectories. They constrain the bicycle on the curve and, given a desired velocity signal $v(t)$, they find a trajectory with the property that the velocity of the bicycle is $v(t)$ and its roll angle $\varphi$ is in the interval $(-\pi/2,\pi/2)$, i.e, the bicycle doesn’t fall over. In [@HauSac06], Hauser-Saccon develop an algorithm to compute the minimum-time speed profile for a point-mass motorcycle compatible with the constraint that the lateral and longitudinal accelerations do not make the tires slip, and apply their algorithm to Getz’s bicycle model.

The problem of maneuvering Getz’s bicycle along a closed curve is equivalent to moving the pivot point of an inverted pendulum around the curve without making the pendulum fall over. On the other hand, the seemingly different problem of maneuvering Hauser’s PVTOL aircraft [@HauSasMey92] along a closed curve in the vertical plane can be viewed as the problem of moving the pivot of an inverted pendulum around the curve without making the pendulum fall over. The two problems are, therefore, closely related, the main difference being the fact that in the former case the pendulum lies on a plane which is orthogonal to the plane of the curve, while in the latter case it lies on the same plane. In [@ConMagNieTos10], the path following problem for the PVTOL was solved by enforcing a virtual holonomic constraint (VHC) which specifies the roll angle of the PVTOL as a function of its position on the curve. In this paper we follow a similar approach for the bicycle model, and impose a VHC relating the bicycle’s roll angle to its position along the curve. However, rather than finding one VHC, as we did in [@ConMagNieTos10], we show how to automatically generate VHCs as periodic solutions of a scalar periodic differential equation which we call the [*virtual constraint generator*]{}. We show that if the path is sufficiently long compared to the height of the bicycle’s centre of mass and the wheel base, then the VHC can be chosen so that on the constraint manifold the bicycle traverses the entire curve with bounded speed, and its speed profile is periodic in steady-state. Finally, we design a controller that enforces the VHC, and recovers the asymptotic properties of the bicycle on the constraint manifold.

![Getz’s bicycle model.[]{data-label="fig:bicycle"}](FIGURES/bicycle){width=".5\textwidth"}

The concept of VHC is a promising paradigm for motion control. It is one of the central ideas in the work of Grizzle and collaborators on biped locomotion (e.g., [@PleGriWesAbb03] and [@WesGriKod03]), where VHCs are used to encode different walking gaits. The work of Shiriaev and collaborators in [@Can04; @ShiPerWit05; @FreRobShiJoh08] investigates VHCs for Lagrangian control systems, i.e., systems of the form [@Blo03] $$\label{eq:lagrangian}
\frac{d}{dt} \frac{\partial L}{\partial \dot q} - \frac{\partial
  L}{\partial q} = G u,$$ with control input $u$ and smooth Lagrangian $L(q,\dot q) = (1/2) \dot
q\trans D(q) \dot q - V(q)$, with $D = D\trans$ positive definite. In [@ShiPerWit05], the authors consider systems of the form  with degree of underactuation one. They find an integral of motion for the constrained dynamics, and use it to select a desired closed orbit on the constraint manifold. This orbit is then stabilized by linearizing the control system along it, and designing a time-varying controller for the linearization. In [@ShiFreGus10], these ideas are extended to systems with degree of underactuation greater than one, and in [@ShiFreRobJohSan07] they are applied to the stabilization of oscillations in the Furuta pendulum. In [@MagCon13], we investigated VHCs for Lagrangian control systems with degree of underactuation one. We introduced and characterized a notion of regularity of VHCs, and we presented sufficient conditions under which the reduced dynamics on the constraint manifold (described by a second-order unforced system) are Lagrangian (i.e., they satisfy the Euler-Lagrange equations, which have the form  with zero right-hand side). An outcome of this paper (Proposition \[prop:internal\_dynamics\_general\_properties\]) is a simple sufficient condition under which the reduced dynamics are [ *not*]{} Lagrangian. We refer the reader to Remarks \[rem:shiriaev\], \[rem:Lagrange\], and \[rem:not\_EL\] for more details.

This paper is organized as follows. In Section \[sec:problem\] we present Getz’s bicycle model and we formulate the maneuvering problem investigated in this paper. In Section \[sec:vcg\] we present the virtual constraint generator idea. The main result is Proposition \[prop:generator:periodic\_solutions\] which gives a constructive methodology to find VHCs for Getz’s bicycle that meet the requirements of the maneuvering problem. In Section \[sec:constrained\_motion\] we analyze the motion of Getz’s bicycle on the virtual constraint manifold. In Proposition \[prop:internal\_dynamics\_general\_properties\] we provide a general result with sufficient conditions under which an unforced second-order system of a certain form possesses an asymptotically stable closed orbit. In Proposition \[prop:internal\_dynamics\_bicycle\] we apply this result to Getz’s bicycle model. In Section \[sec:solution\] we bring these results together and solve the maneuvering problem. Finally, in Section \[sec:numerics\] we make remarks on numerical implementation of the proposed controller.

[**Notation.**]{} Throughout this paper we use the following notation. If $x$ is a real number and $T>0$, the number $x$ modulo $T$ is denoted by $[x]_T$. We let $[\Re]_T : = \{[x]_T: x \in \Re\}$. The set $[\Re]_T$ is diffeomorphic to the unit circle. We let $\pi: \Re
\to [\Re]_T$ be defined as $\pi(t) = [t]_T$. Then, $\pi$ is a smooth covering map (see [@Lee13 p.91]). If $M$ is a smooth manifold, and $h: [\Re]_T \to M$ is a smooth function, we define $\tilde h:=h
\circ \pi : \Re \to M$. This is a $T$-periodic function. Moreover, by [@Lee13 Theorem 4.29], $\tilde h$ is smooth if and only if $h$ is smooth. Finally, $\image(h)$ denotes the image of a function $h$.

Problem formulation {#sec:problem}
===================

Consider the bicycle model depicted in Figure \[fig:bicycle\], with the following variable conventions (taken from [@HauSacFre04]):

-   $(x,y)$ - coordinates of the point of contact of the rear wheel

-   $\varphi$ - roll angle (a positive $\dot \varphi$ implies that the bicycle leans to the right)

-   $\psi$ - yaw angle (a positive $\dot \psi$ means that the bicycle turns right)

-   $\alpha$ - projected steering angle on the $(x,y)$ plane

-   $b$ - distance between the projection of the centre of mass and the point of contact of the rear wheel

-   $p$ - wheel base

-   $h$ - pendulum length

-   $v$ - forward linear velocity of the bicycle

-   $f$ - thrust force.

We denote $\bar \kappa = (\tan \alpha) / p = \dot \psi / v$. For a given velocity signal $v(t)$ and steering angle signal $\alpha(t)$, $\bar \kappa(t)$ represents the curvature of the path $(x(t),y(t))$ traced by the point of contact of the rear wheel. In [@GetMar95], the bicycle of Figure \[fig:bicycle\] was modelled by writing the Lagrangian of the unconstrained bicycle, incorporating the nonholonomic constraints that the wheels roll without slipping in the Lagrangian, and then extracting the model through the Lagrange-d’Alembert equations as in [@Blo03 Section 5.2]. The resulting model, which we’ll refer to as [*Getz’s bicycle model*]{}, reads as $$\begin{aligned}
&\dot {\bar \kappa} = \tau \nonumber\\
 &M \begin{bmatrix} \ddot \varphi \\ \dot v 
\end{bmatrix} = 
F + B  \begin{bmatrix}
  \tau \\ f
\end{bmatrix}, \label{eq:original_system}\end{aligned}$$ where $\tau$, the time derivative of the [curvature $\bar
  \kappa(t)$]{}, and $f$ are the control inputs and, denoting ${s_\varphi}=
\sin \varphi$, ${c_\varphi}=\cos \varphi$, $$\begin{aligned}
& M = \begin{bmatrix} h^2 & b h {c_\varphi}\bar \kappa \\
b h {c_\varphi}\bar \kappa & 1 + (b^2+h^2 {s_\varphi}^2) \bar \kappa^2 - 2 h {\bar \kappa} {s_\varphi}\end{bmatrix}, \\
& F= \begin{bmatrix} g h {s_\varphi}- (1 - h
   {\bar \kappa} {s_\varphi}) h {c_\varphi}{\bar \kappa} v^2 \\ (1 - h {\bar \kappa} {s_\varphi}) 2 h {c_\varphi}{\bar \kappa}
   v \dot \varphi + b h {\bar \kappa} {s_\varphi}\dot \varphi^2
 \end{bmatrix}, \\
& B =\begin{bmatrix} -b h {c_\varphi}v & 0 \\ -(b^2 {\bar \kappa} -
  h {s_\varphi}(1-h {\bar \kappa} {s_\varphi})) v & 1/m
  \end{bmatrix}.
\end{aligned}$$ In this model, $M$ is the inertia matrix, $F$ represents the sum of Coriolis, centrifugal and conservative forces, and $B$ is the input matrix. Now consider a $C^3$ closed Jordan curve $\C$ in the $(x,y)$ plane with regular parametrization $\sigma: [\Re]_T \to \Re^2$, not necessarily unit speed. Let $\kappa:[\Re]_T \to \Re$ denote the signed curvature of $\C$. Throughout this paper, we assume the following.

\[assump\_convex\] The curve $\C$ is strictly convex, i.e., $\kappa(s) > 0$ for all $s
\in [\Re]_T$.

In this paper we investigate the dynamics of the bicycle when the point $(x,y)$ is made to move along the curve $\C$ by an appropriate choice of the steering angle. In order to derive the constrained dynamics, suppose that $( x(0), y(0) ) \in \C$, i.e., $(x(0),y(0)) =
\sigma(s_0)$, for some $s_0 \in [\Re]_T$. A point $\sigma(s(t))$ moving on $\C$ has linear velocity $$\label{eqn_for_v}
v(t) = \|\sigma'(s(t))\| \dot s(t)$$ and acceleration $$\label{eqn_for_v_dot}
\dot v(t) = \|\sigma'(s(t))\| \ddot s(t) + \frac{\dot s^2(t)}{
 \|\sigma'(s(t))\|} \sigma'(s(t))\trans \sigma''(s(t)).$$

For an arbitrary input signal $f(t)$, $(x(t), y(t))$ traverses $\C$ with velocity $v(t)$ if and only if $(x(0),y(0)) \in \C$, $(\dot
x(0),\dot y(0))$ is tangent to $\C$, $\bar \kappa(0) = \kappa(s_0)$, and the input signal $\tau(t)$ is chosen to be $\tau(t) = \kappa'(s(t)) \dot s(t)$, where $\kappa'(s(t)) =
\frac{d\kappa} {ds} (s(t))$. With this choice, we obtain $$\label{eqn_for_bar_kappa}
{\bar \kappa}(t) = \kappa(s(t)),$$ where $s(t)$ and $\dot s(t)$ are solutions of a differential equation to be specified later. The motion of the bicycle on the curve $\C$ is now found by substituting (\[eqn\_for\_v\]), (\[eqn\_for\_v\_dot\]), (\[eqn\_for\_bar\_kappa\]) and $\tau = \kappa'(s) \dot s$ in (\[eq:original\_system\]): $$\label{eq:constrained_system}
\bar{M} \begin{bmatrix} \ddot \varphi \\ \|\sigma^\prime\| \ddot {s} +
\frac{(\sigma')\trans \sigma''}{\|\sigma'\|} \dot{s}^2 
\end{bmatrix} = 
\bar F + \begin{bmatrix} 0 \\1/m \end{bmatrix} f\;,$$ where $\bar M=M|_{\bar \kappa =\kappa(s)}$ and $$\bar F = \left(F+B \begin{bmatrix} 1 \\ 0
\end{bmatrix}\kappa'(s) \dot s
\right)\bigg|_{\bar \kappa =\kappa(s), v=\|\sigma'(s)\|\dot s}.$$ Note that we have used the control input $\tau$ to make the curve $\C$ invariant, so in  we are left with one control input, the thrust force $f$. One can check that system  is a Lagrangian control system with input $f$, i.e., it has the form  with $q=
(\varphi,s)\in \cQ =S^1 \times [\Re]_T$, $G = [0 \ \ 1/m]\trans$, and $$L(q,\dot q) = \frac 1 2  [\dot
  \varphi \ \ \|\sigma^\prime(s)\| \dot{s}] \bar{M}(q) \begin{bmatrix}
  \dot \varphi \\ \|\sigma^\prime(s)\| \dot{s}
\end{bmatrix}- g h \cos \varphi.$$ Since the control force $f$ enters nonsingularly in the $\ddot s$ equation, we can define a feedback transformation for $f$ in  such that $\ddot s=u$, where $u$ is the new control input. With this transformation, the motion of the bicycle when its rear wheel is made by feedback control to follow $\C$ reads as $$\label{eq:sys}
\begin{aligned}
& \ddot \varphi = h^{-1}g {s_\varphi}- h^{-1}\Big[ (1 - h \kappa(s)
    {s_\varphi})\kappa(s) \|\sigma'(s)\| + b \kappa'(s) \\
&  +\frac{b \kappa(s)}{\|\sigma'(s)\|^2} {\sigma'(s)}\trans
   \sigma''(s)\Big]{c_\varphi}\|\sigma'(s)\| \dot s^2 - a(s) {c_\varphi}u \\
& \ddot s = u,
\end{aligned}$$ where $a(s) = b h^{-1} \kappa(s) \|\sigma'(s)\|$. In the above equation, $u$ is the new control input, and it represents the acceleration of the curve parameter $s$. We will denote $\cX := \{
(q,\dot q) \in S^1 \times [\Re]_T \times \Re^2\}$ the state space of .

[If $\sigma(s)$ is a unit speed parametrization of $\C$ (i.e. it satisfies $\|\sigma'(s)\|=1$), then the first differential equation in  reduces to $h \ddot \varphi = g {s_\varphi}- \Big[ (1 - h \kappa(s) {s_\varphi})\kappa(s) + b
 \kappa'(s)\Big]{c_\varphi}\dot s^2 - b \kappa(s) {c_\varphi}u$. ]{}

[Maneuvering Problem]{} Find a feedback $u(q,\dot q)$ for system  such that there exists a set of initial conditions $\Omega$ with the property that, for all $( q(0),\dot q(0) ) \in \Omega$, the bicycle does not overturn, i.e., $|\varphi(t)| < \pi/2$ for all $t \geq 0$, and traverses the entire curve $\C$ in one direction, i.e., there exists $\bar t>0$ such that $| \dot s(t) | >0$ for all $ t \geq \bar t$. Moreover, the speed $\dot s(t)$ of the bicycle on $\C$ should remain bounded.

Our solution of this problem relies on the notion of VHC.

A [**virtual holonomic constraint (VHC)**]{} for system  is a relation $\varphi = \Phi(s)$, where $\Phi: [\Re]_T \to S^1$ is smooth and the set $\Gamma = \{(q,\dot q)\in \cX: \varphi = \Phi(s), \ \dot \varphi = \Phi'(s)
\dot s\}$ is controlled invariant. That is, there exist a smooth feedback $u(q,\dot q)$ such that $\Gamma$ is invariant for the closed-loop system. The set $\Gamma$ is called the [**constraint manifold**]{} associated with the VHC $\varphi = \Phi(s)$.

The definition above formalizes the notion of VHC used in [@WesGriKod03] in the context of biped locomotion. The constraint manifold $\Gamma$ is a two-dimensional submanifold of $\cX$ parametrized by $(s,\dot s)$, and therefore diffeomorphic to the cylinder $[\Re]_T \times \Re$. It is the collection of all those phase curves of  that can be made to satisfy the constraint $\varphi = \Phi(s)$ via feedback control. In order to solve the maneuvering problem, we look for VHCs $\varphi = \Phi(s)$ such that $| \Phi (s) | < \pi/2$ for all $s \in [\Re]_T$, and then stabilize the associated virtual constraint manifold. In [@HauSacFre04], Hauser-Saccon-Frezza find “bounded roll trajectories,” i.e., controlled trajectories of  along which the roll angle $\varphi$ is bounded in the interval $(-\pi/2,\pi/2)$. In our context, each VHC $\varphi = \Phi(s)$ provides a [*family*]{} of bounded roll trajectories. Once $\Gamma$ has been made invariant via feedback control, bounded roll trajectories can be obtained by picking arbitrary $(s_0,\dot s_0) \in
[\Re]_T \times \Re$, and picking as initial condition in , $(\varphi(0),s(0)) = (\Phi(s_0),s_0)$, $(\dot
\varphi(0),\dot s(0)) = (\Phi'(s_0)\dot s_0,\dot s_0)$. The resulting solution $(q(t),\dot q(t))$ will satisfy $\varphi(t) = \Phi(s(t))$, implying that the roll angle trajectory $\varphi(t)$ is bounded in the interval $(-\pi/2,\pi/2)$.

The virtual constraint generator {#sec:vcg}
================================

In this section we show that VHCs for  can be generated as solutions of a first-order differential equation, which we call the [*VHC generator.*]{} This idea was first presented in our previous work [@ConMag10_1]. We begin with a sufficient condition for a relation $\Phi$ to be a VHC for .

\[lem:feasibility\] A relation $\varphi = \Phi(s)$, where $\Phi:[\Re]_T \to S^1$ is smooth, is a VHC for system  if $$\label{eq:feasibility}
(\forall s \in [\Re]_T) \ \Phi'(s) + a(s) \cos \Phi(s) \neq 0.$$

Letting $H(\varphi,s):=\varphi - \Phi(s)$, condition  is the requirement that system  with output $H(\varphi,s)$ has a well-defined uniform relative degree 2 on the set $\{(q,\dot q)\in \cX: \varphi =
\Phi(s)\}$. The associated zero dynamics manifold is precisely $\Gamma$, and it is controlled invariant.

The foregoing lemma inspires the following observation. Instead of guessing a relation $\varphi = \Phi(s)$ and checking whether it is a VHC for , as in Lemma \[lem:feasibility\], we will assign a nonzero right-hand side to , view the resulting identity as an ODE, and [*generate*]{} VHCs by solving the ODE. More precisely, recall the notation $\tilde h:=h \circ \pi$, and consider the scalar $T$-periodic differential equation $$\label{eq:constraint_generator}
\dot x = - \tilde a(t)\cos x + \tilde \delta(t),$$ where $\tilde \delta:\Re \to \Re \backslash \{0\}$ is a $T$-periodic function to be assigned. Then, $T$-periodic solutions of  give rise to VHCs.

\[lem:vcg\] Suppose that $x(t)$ is a $T$-periodic solution of , where $\tilde \delta :\Re \to \Re
\backslash \{0\}$ is smooth and $T$-periodic. Then, the relation $\varphi = \Phi(s)$, with $\Phi:=x \circ \pi^{-1}:[\Re]_T \to \Re$, is a VHC for .

Let $x(t)$ be a $T$-periodic solution of . Then, $x(t)$ is smooth because the right-hand side of  is smooth. Since $x(t)$ is $T$-periodic, $\Phi:=x \circ \pi^{-1}$ is a well-defined function $[\Re]_T \to \Re$, and by [@Lee13 Theorem 4.29] it is smooth. Since $x = \Phi \circ \pi$, we have $\dot x(t) =
\Phi'(\pi(t)) \dot \pi(t) = \Phi'(\pi(t))$, so that $\Phi'(s) =\dot
x(t)|_{t = \pi^{-1}(s)}$. Using , we get $$\begin{aligned}
\Phi'(s) & = -\tilde a(\pi^{-1}(s)) \cos \Phi(s) + \tilde
\delta(\pi^{-1}(s)) \nonumber \\
&=- a(s) \cos \Phi(s) +
\delta(s). \label{eq:Phiprime}\end{aligned}$$ Since $\delta(s) \neq 0$ for all $s \in [\Re]_T$, by Lemma \[lem:feasibility\] the relation $\varphi = \Phi(s)$ is a VHC for system .

In light of Lemma \[lem:vcg\], we call the differential equation  a [*VHC generator*]{}, for which one is to pick a $T$-periodic $\tilde \delta:\Re \to \Re
\backslash \{0\}$ yielding a $T$-periodic solution. The next proposition shows how to pick $\tilde \delta$.

\[prop:generator:periodic\_solutions\] Consider the smooth and $T$-periodic differential equation , and set $\tilde \delta(t) =
\epsilon \tilde \mu(t)$, where $\tilde \mu: \Re \to (0,\infty)$ is smooth and $T$-periodic. Let $K^+ = \max_{t \in [0, T]}[ \tilde
  \mu(t)/\tilde a(t)]$, $K^- = \min_{t \in [0,T]} [\tilde
  \mu(t)/\tilde a(t)]$. Then, for any $x_0 \in (0,\pi/2)$ and $t_0
\in \Re$, the following two properties hold:

(i) there exists a unique $\epsilon \in [\epsilon^-,\epsilon^+] = [
      ( c_{x_0}) / K^+,$ $( c_{x_0}) / K^- ]$, such that the solution $x(t)$ of  with initial condition $x(t_0) = x_0$ is $T$-periodic, and setting $\Phi=x \circ \pi^{-1}$, the relation $\varphi = \Phi(s)$ is a VHC for .

(ii) If $\tilde\mu:\Re \to (0,\infty)$ is chosen so that $K^+ / K^- <
      (\cos x_0)^{-1}$, then the VHC $\varphi=\Phi(s)$ in part (i) satisfies $\image(\Phi) \subset (\Phi^-,\Phi^+) \subset (0,\pi/2)$, where $\Phi^-=\cos^{-1}\big(\frac{K^+}{K^-} c_{x_0}\big)$ and $\Phi^+=\cos^{-1}\big(\frac{K^-}{K^+} c_{x_0}\big)$.

\[rem:mu\] [By choosing $\tilde \mu(t) \equiv \tilde a(t)$, we have $K^+ = K^- =1$, $\epsilon^+ = \epsilon^- = c_{x_0}$. Proposition \[prop:generator:periodic\_solutions\] implies that, for all $x_0 \in (0,\pi/2)$, setting $\tilde \delta(t) = \cos (x_0)
  \tilde a(t)$, the VHC generator produces a $T$-periodic solution $x(t)$ with image contained in $(0,\pi/2)$. As a matter of fact, the solution in question is $x(t) \equiv x_0$, so that the relation $\varphi =x_0$ is a VHC for . The bicycle subject to this VHC has a constant roll angle as it travels around $\C$. The proposition provides great flexibility in finding VHCs with the property that the roll angle is confined within the interval $(0,\pi/2)$. All such constraints are compatible with the maneuvering problem.]{}

Pick an arbitrary $x_0 \in (0,\pi/2)$ and $t_0 \in \Re$, and consider the VHC generator $$\label{eq:constraint_generator:proof}
\dot x = - \tilde a(t) \cos x + \epsilon \tilde \mu(t),$$ where $\tilde \mu(\cdot)>0$. Let $\cS(\epsilon)$ denote the solution of  at time $t_0+T$ with initial condition $x(t_0) = x_0$. Since the right-hand side of  is $T$-periodic, the solution of  with initial condition $x(t_0) =x_0$ is $T$-periodic if and only if $\cS(\epsilon) = x_0$. For all $\epsilon > \epsilon^+$ we have $$\begin{aligned}
\dot x & > -\tilde a(t) \cos x + \tilde \mu(t) \cos(x_0) / K^-\\
 & > -\tilde a(t) \left[ \cos x - \cos x_0 (\tilde \mu(t) / \tilde
  a(t) ) / K^- \right]] \\
& > -\tilde a(t) (\cos x - \cos x_0).
\end{aligned}$$ The solution of $\dot x = -\tilde a(t) (\cos x - \cos x_0)$ with initial condition $x(t_0) = x_0$ is $x(t) \equiv x_0$. Therefore, by the comparison lemma ([@Kha02 Lemma 3.4]), for all $\epsilon>
\epsilon^+$ the solution of  with initial condition $x(t_0) = x_0$ satisfies $x(t)> x_0$ for all $t >
t_0$, so that $\cS(\epsilon)> x_0$ for all $\epsilon> \epsilon^+$. A similar argument can be used to show that for all $\epsilon<
\epsilon^-$ the solution of  satisfies $x(t) < x_0$ for all $t>t_0$, so that $\cS(\epsilon)< x_0$ for all $\epsilon<\epsilon^-$. By continuity of solutions with respect to parameters, the map $\cS$ is continuous, and therefore there exists $\epsilon \in [\epsilon^-,\epsilon^+]$ such that $\cS(\epsilon) =
x_0$. The corresponding solution $x(t)$ is $T$-periodic. The same comparison argument shows that if $0<\epsilon_1 < \epsilon_2$, and for $i=1,2$, $x_i(t)$ is the solution of  with $\epsilon = \epsilon_i$, then $x_1(t_0+T) < x_2(t_0+T)$, so that $\cS(\epsilon_1) <
\cS(\epsilon_2)$. In other words, $\cS(\epsilon)$ is a monotonically increasing function, and so the value of $\epsilon$ above is unique. This concludes the proof of part (i). Now suppose that $K^+ /
K^- < (\cos x_0)^{-1}$, and let $\epsilon \in [\epsilon^-,\epsilon^+]$ be such that the solution $x(t)$ of  with initial condition $x(t_0)=x_0$ is $T$-periodic, as in part (i). Let $x^+ =
\cos^{-1}\big(\frac{K^-}{K^+} c_{x_0}\big)$, and consider the subset of the extended phase space $S^+= \{(x,t): x > x^+\}$. This subset is positively invariant since $\tilde a>0$ and $$\begin{aligned}
\dot x\big|_{x=x^+} &= -\tilde a(t) \frac{K^-}{K^+} \cos x_0 +
\epsilon \tilde \mu(t) \\
&\hspace*{-5ex}\geq -\tilde a(t)
\frac{K^-}{K^+} \cos x_0 + (\cos x_0 / K^+) \tilde \mu(t)\\
&\hspace*{-5ex}\geq \frac{\tilde a(t)}{K^+} \cos x_0 \left( -K^- +
\tilde \mu(t)/\tilde a(t)\right)\\
&\hspace*{-5ex}\geq 0.
\end{aligned}$$ Similarly, letting $x^- = \cos^{-1}\big(\frac{K^+}{K^-} c_{x_0}\big)$, the subset $S^-= \{(x,t): x < x^-\}$ is positively invariant. Since $\{ (x,t): x= x_0\}$ has empty intersection with $S^+ \cup S^-$, and since $S^+ \cup S^-$ is positively invariant, it follows that the $T$-periodic solution $x(t)$ with initial condition $x(t_0)=x_0$ must be contained in the complement of $S^+ \cup S^-$, i.e., for all $t\in\Re$, $x^- \leq x(t) \leq x^+$. For, if for some $\bar t \in
\Re$, $x(\bar t) \in S^+ \cup S^-$, then the fact that $x(t) \in S^+
\cup S^-$ for all $t\geq \bar t$ would contradict the periodicity of $x(t)$.

\[ex:ellipse\]

Suppose $\C$ is an ellipse with major semiaxis $A$, minor semiaxis $B$, and $2 \pi$-periodic parametrization $\sigma(s) = (A
  \cos s, B \sin s)$, with $A=15$, $B=10$. The curvature is [$\kappa(s) = AB / ( A^2 \sin^2 s + B^2 \cos^2 s
    )^{3/2}$]{}. For the initial condition of the VHC generator , we pick $x(0) = \pi/8$. Following Proposition \[prop:generator:periodic\_solutions\], we need to choose a $2 \pi$-periodic function $\tilde \mu(t)>0$, set $\tilde \delta(t) = \epsilon \tilde \mu(t)$, and find the unique value of $\epsilon>0$ guaranteeing that the solution with initial condition $x(0) = \pi/8$ is $2\pi$-periodic. There is much freedom in the choice of $\tilde \mu(t)$. For instance, picking $\tilde
  \mu(t) = 1$, we numerically find [$\epsilon \approx
    0.6482$]{}. The corresponding VHC is depicted in Figure \[fig:ellipse:virtual\_constraint\]. The condition, in Proposition \[prop:generator:periodic\_solutions\](ii), that $K^+/
  K^- < (\cos x_0)^{-1}$ is conservative. Indeed, with our choice of $\tilde \mu$ we have $K^+ =3/2$, $K^-=2/3$, and thus the condition is violated. Yet, we see in Figure \[fig:ellipse:virtual\_constraint\] that $ \Phi(s) \in (0,
  \pi/2)$ for all $s \in [\Re]_T$.

![VHC for the ellipse in Example \[ex:ellipse\].[]{data-label="fig:ellipse:virtual_constraint"}](FIGURES/ellipse_constraint_pend "fig:"){width=".495\textwidth"} ![VHC for the ellipse in Example \[ex:ellipse\].[]{data-label="fig:ellipse:virtual_constraint"}](FIGURES/ellipse_constraint_fun "fig:"){width=".495\textwidth"}

Motion on the constraint manifold {#sec:constrained_motion}
=================================

Having chosen $\tilde \delta(t)= \epsilon \tilde \mu(t)$ as in Proposition \[prop:generator:periodic\_solutions\](ii), and having obtained a VHC $\varphi= \Phi(s)$, the next step is to analyse the [*reduced dynamics*]{} of  on the constraint manifold $\Gamma = \{ (q,\dot q) \in \cX : \varphi = \Phi(s), \ \dot \varphi =
\Phi'(s) \dot s\}$. These are the zero dynamics of  with output function $H(\varphi,s)=\varphi - \Phi(s)$. To this end, we impose that $[(d/dt) (\Phi'(s) \dot s)] |_\Gamma = \ddot
\varphi|_\Gamma$. Expanding both sides of the equation above, using the expression of $\ddot \varphi$ in , identity , and the fact that $\delta \neq 0$, we obtain the feedback making $\Gamma$ invariant $$\label{eqn_inv_feedback}
\begin{aligned}
& u = \frac{h^{-1}g \sin \Phi}{\delta} - \frac{ \dot s^2
}{\delta} \Big[ \Phi'' + \frac 1 h ( (1-h \kappa \sin \Phi) \kappa
  \|\sigma'\| \\
& + b \kappa' +b\kappa {\sigma'}\trans\sigma'' / \|\sigma'\|^2)
  \|\sigma'\| \cos \Phi \Big].
\end{aligned}$$ Substituting feedback  in the $s$ dynamics, we get the reduced dynamics on $\Gamma$ $$\label{eq:motion_on_Gamma}
\ddot s = \Psi_1(s) + \Psi_2(s) \dot s^2,$$ with $$\label{eq:Psi}
\begin{aligned}
&\Psi_1 = \frac{h^{-1}g s_{\Phi}}{\delta} \\
& \Psi_2= - \frac{1}{\delta} \Big[ \Phi'' + \frac 1 h ( (1-h
    \kappa s_{\Phi}) \kappa\|\sigma'\| \\ & + b \kappa'
    + b \kappa \sigma'{}\trans \sigma'' / \|\sigma'\|^2 )
    c_{\Phi} \|\sigma'\| \Big].
\end{aligned}$$ System  describes the motion of system  on the constraint manifold $\Gamma$ in the following sense. For a given initial condition $(s(0),\dot s(0)) =
(s_0,\dot s_0) \in [\Re]_T \times \Re$, let $(s(t),\dot s(t))$ be the corresponding solution of , and let $(\varphi(t),\dot \varphi(t), s_1(t),\dot s_1(t))$ be the solution of  with initial condition $(\varphi(0),\dot
\varphi(0),s(0),\dot s(0)) = ( \Phi(s_0), \Phi'(s_0) \dot s_0,s_0,\dot
s_0) \in \Gamma$. Then, for all $t \geq 0$ $(s_1(t),\dot s_1 (t)) =
(s(t),\dot s(t))$, and $(\varphi(t),\dot \varphi(t)) =$ $(\Phi(s(t)),$ $\Phi'(s(t)) \dot s(t))$. In order for the VHC $\varphi=\Phi(s)$ to be compatible with the maneuvering problem, we need to verify whether or not the rear wheel of the bicycle traverses the entire curve $\C$ with bounded speed, i.e., we need to show that for any solution $(s(t),\dot s(t))$ of , there exist $\bar
t>0$ and $\epsilon>0$ such that $\dot s(t)>\epsilon > 0$ for all $ t
\geq \bar t$, and $\lim\sup_{t \to \infty} \dot s(t) < \infty$. The next result explores general properties of systems of the form .

\[prop:internal\_dynamics\_general\_properties\] Consider a differential equation of the form  with state space $\{(s,\dot s) \in
[\Re]_T \times \Re\}$. Assume that $\Psi_1, \Psi_2:[\Re]_T \to \Re$ are smooth functions such that $\Psi_1(s) > 0$ for all $s$ and $\int_0^T \tilde \Psi_2(\tau) d\tau <0$. Then, there exists a smooth function $\nu:[\Re]_T \to (0,\infty)$ such that the closed orbit $\cR
= \{(s,\dot s) \in [\Re]_T \times \Re: \dot s = \nu(s)\}$ is exponentially stable for , with domain of attraction containing the set $\cD = \{(s,\dot s) \in [\Re]_T \times
\Re: \dot s \geq 0\}$. [Thus,]{} for all initial conditions in $\cD$, the solution $(s(t),\dot s(t))$ of  converges to the unique asymptotically stable limit cycle $\cR$.

\[rem:doa\] [It can be shown that the domain of attraction of the limit cycle $\cR$ is $\{(s,\dot s) : \dot s > -
  \nu(s)\}$.]{}

\[rem:shiriaev\] [Proposition \[prop:internal\_dynamics\_general\_properties\] has general implications. Consider a Lagrangian control system of the form  with configuration vector $q=
  (q_1,\ldots, q_n)$ and degree of underactuation one. Consider a VHC of the form $q_i = \Phi_i(q_n)$, $i=1,\ldots,n-1$, and suppose $q_n
  \in [\Re]_T$. For this system, it is shown in [@WesGriKod03; @ShiPerWit05; @MagCon13] that the reduced dynamics have the form , with $s$ replaced by $q_n$. In this context, Proposition \[prop:internal\_dynamics\_general\_properties\] gives sufficient conditions under which the reduced dynamics have an asymptotically stable closed orbit, implying that they are not Lagrangian. On the other hand, in Theorem 2 of [@Shiriaev2006900] (see also [@ShiPerWit05]) it is shown that for any $(s_0,\dot s_0)$, system  possesses an “integral of motion”[^2] $I(s,\dot s) = \dot s^2 +2 V(s)/ M(s) - \dot
  s_0^2/M(s)$, where $ M(s) = \exp \left\{ -2 \int_{s_0}^s
  \Psi_2(\tau) d \tau \right\}$, $V(s) = -\int_{s_0}^s \Psi_1(\mu)
  M(\mu) d\mu$. This integral of motion predicted by Theorem 2 of [@Shiriaev2006900] seemingly contradicts Proposition \[prop:internal\_dynamics\_general\_properties\]. Indeed, if $M(s)$ and $V(s)$ are well-defined functions, letting $L(s,\dot
  s) = (1/2) M(s) \dot s^2 - V(s)$ one has that $\frac{d}{dt}
  \frac{\partial L}{\partial \dot s} - \frac{\partial L}{\partial s} =
  0$, and hence the reduced dynamics  are Lagrangian. This fact rules out the existence of isolated closed orbits for  predicted by Proposition \[prop:internal\_dynamics\_general\_properties\]. This contradiction is due to the fact that Theorem 2 in [@Shiriaev2006900] implicitly assumes that the functions $M(s)$ and $V(s)$ are well-defined. This is always true if $s \in
  \Re$. However, when $s \in [\Re]_T$, the functions $M(s)$ and $V(s)$ may be multi-valued, in which case the Lagrangian $L(s,\dot s)$ is undefined. Indeed, the assumptions of Proposition \[prop:internal\_dynamics\_general\_properties\] imply that both $M(s)$ and $V(s)$ are multi-valued. We refer the reader to Section IV of [@MagCon13] for a detailed discussion on this subject, and to [@MohMagCon13] for necessary and sufficient conditions under which  is Lagrangian.]{}

\[rem:Lagrange\] [In Lagrangian mechanics, the Lagrange-d’Alembert principle implies that the enforcement of an ideal holonomic constraint (i.e., a constraint for which the constraint forces do not produce work) on a Lagrangian system gives rise to a reduced Lagrangian system. Proposition \[prop:internal\_dynamics\_general\_properties\] shows that this result is not true when enforcing VHCs on Lagrangian control systems. There is, therefore, a sharp difference between [*ideal*]{} and [*virtual*]{} holonomic constraints when it comes to the reduced motion they induce. This difference is due to the fact that VHCs are enforced through forces produced by feedback control. For underactuated systems such as the constrained Getz’s bicycle model in , the control forces produce work. ]{}

[In the context of walking robots, [@WesGriKod03] proved that under certain conditions, the hybrid zero dynamics subject to a VHC exhibit an exponentially stable hybrid limit cycle. The mechanism enabling this exponential stability property is the jump map representing the impact of the robot’s foot with the ground. Proposition \[prop:internal\_dynamics\_general\_properties\] shows that VHCs may induce stable limit cycles even when the control system has no jumps.]{}

Consider the differential equation $$\label{eq:constrained_system_lifted}
\ddot x = \tilde \Psi_1 (x) + \tilde \Psi_2(x) \dot x^2.$$ Let $\bar \pi: \Re \times \Re \to [\Re]_T \times \Re$ be defined as $\bar \pi(x,\dot x) = (\pi(x),\dot x)$. Then, it is readily seen that the vector fields in  and  are $\bar
\pi$-related [@Lee13]. It follows that if $(x(t),\dot x(t))$ is a solution of , then $(s(t),\dot
s(t)) = (\pi(x(t)),\dot x(t))$ is a solution of . We will show that there exists a smooth $T$-periodic function $\tilde \nu : \Re \to (0,+\infty)$ such that the set $\tilde \cR = \{(x,\dot x) \in \Re^2: \dot x = \tilde
\nu(x)\}$ is an exponentially stable orbit of , with domain of attraction containing $\tilde \cD = \{(x,\dot x)\in \Re^2: \dot x\geq 0\}$. Then, setting $\nu = \tilde \nu \circ \pi^{-1}$, by [@Lee13 Theorem 4.29] we will obtain a smooth function $[\Re]_T \to (0,+\infty)$, and the set $\cR= \{(s,\dot s) \in [\Re]_T
\times \Re: \dot s = \nu(s)\}$ will be an exponentially stable closed-orbit of  with domain of attraction containing $\cD$, proving the proposition. The set $\tilde \cD$ is positively invariant for  because, by assumption, $\ddot x|_{\dot x=0} = \tilde \Psi_1(x) >0$ for all $x\in \Re$. In the rest of the proof we will restrict initial conditions to lie in $\tilde \cD$. Letting[^3] $z = \dot
x^2$, we have $\dot z = 2\dot x (\tilde \Psi_1(x) + \tilde \Psi_2(x)
z)$. For all $(x,\dot x) \in \tilde \cD$, we have $\dot x>0$, so we can use $x$ as a time variable: $$\label{eq:core_subsystem}
\frac { d z } {dx}  = 2\tilde \Psi_1(x) + 2\tilde \Psi_2(x) z.$$ The above is a scalar linear $T$-periodic system. Letting $\phi(x) =
\exp(2\int_0^x \tilde \Psi_2(\tau) d \tau)$, the solution of  with initial condition $z(x_0)=z_0$ is $$z(x) = \phi(x) \phi(x_0)^{-1} z_0 + 2\int_{x_0}^x \phi(x) \phi^{-1}
(\tau) \tilde \Psi_1(\tau) d \tau.$$ Note that, for any integer $k$, $\phi(x+k T) = \phi(T)^k \phi(x)$. System  has a $T$-periodic solution if and only if there exists $\bar z_0$ such that the solution $z(x)$ of  with initial condition $z(x_0) = \bar
z_0$ satisfies, $z(x_0+T) = \bar z_0$, or $$\label{eq:z0}
\bar z_0 = \phi(T) \bar z_0 + 2 \int_{x_0}^{x_0+T}\phi(x_0+T)
\phi^{-1}(\tau) \tilde \Psi_1(\tau) d \tau.$$ By assumption, $0< \phi(T) <1$, so there is a unique $\bar z_0>0$ solving , implying that there is a unique $T$-periodic solution $\bar z:\Re \to \Re$ of . Since for all $x\in \Re$, $z=0 \implies (d z/ dx) = 2\tilde \Psi_1(x)>0$, the set $\{(z,x): z>0\}$ is positively invariant for , and therefore the $T$-periodic function $\bar z$ satisfies $\image(\bar z) \subset (0,\infty)$. Let $z:\Re
\to \Re$ be any other arbitrary solution of  (therefore, not necessarily $T$-periodic). Let $k$ be a nonnegative integer and denote $z_k:=z(x_0+k T)$, $\phi_k:=\phi(x_0+kT)$, $x_k:=x_0+kT$. Then, $$\begin{aligned}
&z_{k+1} = \phi_{k+1} \phi_k^{-1} z_k + 2\int_{x_k}^{x_{k+1}}
  \phi_{k+1} \phi^{-1}(\tau) \tilde \Psi_1(\tau) d \tau \\
& = \phi(T)  z_k +2\int_{x_0}^{x_1} \phi_1 \phi(T)^k
\phi^{-1} (\tau + kT) \tilde \Psi_1(\tau + kT) d \tau \\
& = \phi(T)  z_k +2\int_{x_0}^{x_1} \phi_1
\phi^{-1}(\tau) \tilde \Psi_1(\tau) d \tau \; \text{($\tilde
  \Psi_1$ is $T$-periodic)}\\
& = \phi(T)  ( z_k - \bar z_0) + \bar z_0. \quad
\text{(by~\eqref{eq:z0})} 
\end{aligned}$$ In light of the above, letting $\tilde z_k = z_k -\bar z_0$, we have $\tilde z_{k+1} = \phi(T) \tilde z_k$. Since $0<\phi(T)<1$, the origin of this discrete-time system is globally exponentially stable, proving that the $T$-periodic solution $\bar z:\Re\to(0,\infty)$ is globally exponentially stable for . Define $\tilde \nu:\Re \to (0,\infty)$ as $\nu:= \sqrt{\bar z}$, and return to system . For an arbitrary initial condition $(x(0),\dot x(0)) \in \tilde \cD$, the solution $(x(t),\dot
x(t))$ satisfies $\dot x(t) > 0$ for all $t >0$, and $\dot x(t)
=\sqrt{ z(x(t))}$, where $z:\Re\to\Re$ is the solution of  with initial condition $z(x(0)) =(\dot
x(0))^2$. Since the solution $\bar z:\Re \to(0,\infty)$ is globally exponentially stable for , the set $\tilde
\cR$ is exponentially stable for  with domain of attraction containing $\tilde \cD$. Going back to $(s,\dot s)$ coordinates, this implies that $\cR = \{(s,\dot s) \in
[\Re]_T \times \Re: \dot s = \nu(s)\}$ is exponentially stable for , and its domain of attraction contains $\cD$. Note that $\cR$ is a simple closed curve in $[\Re]_T \times
\Re$. We are left to show that $\cR$ is a closed orbit of . The set $\cR$ is closed, invariant, and no proper subset of it has these properties. For, if there existed a closed and invariant proper subset $\cR' \subset \cR$, then any solution $(s(t),\dot s(t))$ originating in $\cR'$ would not traverse the entire curve $\cR$, contradicting the fact that, for all $(s,\dot
s)\in\cR$, $\dot s >0$. $\cR$ is, therefore, a minimal set for . By [@Har02 Theorem 12.1], $\cR$ is a closed orbit.

We now show that if $\C$ is sufficiently long as compared to $b$ and $h$, the bicycle satisfies the hypotheses of Proposition \[prop:internal\_dynamics\_general\_properties\], and so the reduced motion  satisfies the requirements of the maneuvering problem.

\[prop:internal\_dynamics\_bicycle\] Assume that the curvature of $\C$ satisfies the inequality $$\label{eq:curvature_bound}
\frac 1 L \int_0^T \tilde \kappa(\tau) \|\tilde \sigma'(\tau)\| d\tau
< \frac{h}{b^2+h^2},$$ where $L=\int_0^T \|\tilde \sigma'(\tau)\| d \tau$ is the length of $\C$. Then, any VHC $\varphi = \Phi(s)$ with $\image(\Phi) \subset
(0,\pi/2)$ yields functions $\Psi_1,\Psi_2:[\Re]_T \to \Re$ in  that satisfy the hypotheses of Proposition \[prop:internal\_dynamics\_general\_properties\]. Hence, in particular, the reduced dynamics  have an exponentially stable closed orbit.

[The integral $(1/L)\int_0^T \tilde \kappa(\tau) \|\tilde
  \sigma'(\tau)\| d\tau$ is equal to ([*turning number*]{} of $\C)
  \times 2\pi/{L}$. [The turning number of $\C$ is the number of counterclockwise revolutions that the tangent vector to $\C$ makes as its base point is moved once around $\C$ in a way consistent with the orientation of $\C$. The turning number of a Jordan curve is $\pm 1$, and for curves satisfying Assumption \[assump\_convex\] it is always $1$. Thus, condition  can be written as $L > 2 \pi
    (b^2+h^2)/h$. It requires the curve to be sufficiently long as compared to the bicycle parameters $b$ and $h$.]{}]{}

\[rem:not\_EL\] [Proposition \[prop:internal\_dynamics\_bicycle\] implies that when a VHC $\varphi = \Phi(s)$ with $\image(\Phi) \subset (0,\pi/2)$ is enforced on the Lagrangian control system , the resulting reduced dynamics are not Lagrangian.]{}

The VHC $\varphi = \Phi(s)$ arising from Proposition \[prop:generator:periodic\_solutions\], part (ii), satisfies $\Phi(s) \in (0, \pi/2)$ for all $s\in [\Re]_T$, so that $\sin \Phi(s) >0$. Recall that $\Phi$ satisfies  with $\delta(s) =\epsilon \mu(s)>0$ for all $s \in [\Re]_T$. This implies that $\Psi_1(s)>0$. Moreover, $\Phi'' = \delta' -r \kappa'
\|\sigma'\| c_\Phi - (r \kappa\|\sigma'\|)^2 s_\Phi c_\Phi +
r \kappa \|\sigma'\|\delta s_\Phi -r \kappa {\sigma'}\trans
\sigma'' / \|\sigma'\|$. Substituting this expression in , we get $$\begin{aligned} 
\Psi_2 &= -\frac{\delta'}{\delta} - \frac {\kappa\|\sigma'\|} {\delta
  h} \Big[ b \delta s_\Phi - \|\sigma'\|c_\Phi \big( \kappa s_\Phi
  \frac{b^2+h^2} h -1 \big) \Big] \\
&\leq -\frac{\delta'}{\delta} + \frac{\kappa\|\sigma'\|^2} {\delta h}
c_\Phi \big( \kappa s_\Phi \frac{b^2+h^2} h -1 \big) .
\end{aligned}$$ Since $\int_0^T \tilde \delta'(\tau)/\tilde \delta(\tau) = \ln \tilde
\delta(T) - \ln \tilde \delta(0) =0$, using the fact that $\image(\tilde \Phi) \subset (\Phi^-,\Phi^+) \subset (0,\pi/2)$, we have $$\begin{aligned}
&\int_0^T \tilde \Psi_2(\tau) d\tau \leq \int_0^T \frac{\tilde
    \kappa\|\tilde \sigma'\|^2} {\tilde \delta h} c_{\tilde \Phi}
  \left( \tilde \kappa s_{\tilde \Phi} \frac{b^2+h^2} h -1 \right)
  d\tau \\
&\leq \max_t \left(\frac{\tilde \kappa\|\tilde \sigma'\|} {\tilde
    \delta h}\right) c_{\Phi^-} \int_0^T \left(\tilde \kappa
  \frac{b^2+h^2} h -1 \right) \|\tilde \sigma'\| d\tau \\
&\leq \max_t \left(\frac{\tilde \kappa\|\tilde \sigma'\|^2} {\tilde
    \delta h}\right) c_{\Phi^-} \left( \frac{b^2+h^2} h \int_0^T
  \tilde \kappa \|\tilde \sigma'\| d\tau -L \right)\\
&<0.
\end{aligned}$$

We return to the ellipse of Example \[ex:ellipse\] and the VHC displayed in Figure \[fig:ellipse:virtual\_constraint\]. For this example, $(1/2\pi)\int_0^{2\pi} \tilde \kappa(\tau) d\tau \approx 0.0792$, and $h/(b^2+h^2)=0.6711$ and thus  is satisfied. Indeed, one can numerically check that $\int_0^{2\pi} \tilde
  \Psi_2(\tau) d\tau \approx -105.1<0$, and Proposition \[prop:internal\_dynamics\_general\_properties\] applies. The phase portrait of the dynamics on the constraint manifold is displayed in Figure \[fig:ellipse:phase\_portrait\]. The figure illustrates the set $\cR$, corresponding to the steady-state velocity profile of the bicycle on $\Gamma$. The domain of attraction of $\cR$, shaded in the figure, is the set $\{(s,\dot s):
  \dot s > -\nu(s)\}$, as pointed out in Remark \[rem:doa\].

![Phase portrait of the dynamics on $\Gamma$ and set $\cR$ for the ellipse in Example \[ex:ellipse\]. The shaded region is the domain of attraction of $\cR$. Since $s \in [\Re]_{2\pi}$, the $s$ axis wraps around, and points on the lines $s = 0$ and $s = 2\pi$ are identified in the figure, from which it follows that $\cR$ is a closed orbit.[]{data-label="fig:ellipse:phase_portrait"}](FIGURES/ellipse_phase_portrait){width=".5\textwidth"}

[Let now $\C$ be a circle of radius $R$ with parametrization $\sigma(s) =(R \cos s,R \sin s)$. We have $T=2\pi$, $\|\sigma'(s)\|
  \equiv R$ and $\kappa(s) \equiv 1/R$, and so $a(s) \equiv
  bh^{-1}$. For any $x_0 \in (0,\pi/2)$, picking $\tilde \delta= bh^{-1}
  \cos x_0$, as in Remark \[rem:mu\], we obtain the constant VHC $\varphi=x_0$. Equation  becomes $$\frac{dz}{dx} =\frac{2g}{b} \tan x_0 - \frac{2}{bR} ( 1-
(h/R)\sin x_0) z.$$ The above is a linear time-invariant system with constant input which is stable if $R > h \sin x_0$. The periodic solution $\bar z(t)$ in this case is simply the equilibrium of the system above, $\bar z= g
R^2 \tan x_0 / (R - h \sin x_0) $, and thus the asymptotic velocity of the bicycle on $\Gamma$ is constant, and reads as $\nu = R \sqrt{g
  \tan x_0 / (R - h \sin \Phi_0)}$. It can be verified that $\nu$ is an increasing function of $x_0$. The conclusion is that the bicycle can travel around the circle with any constant roll angle in the interval $(0,\pi/2)$. In steady-state, its speed is constant. The larger is the roll angle $x_0$, the higher is the asymptotic speed of the bicycle. ]{}

Solution of the maneuvering problem {#sec:solution}
===================================

\[thm:solution\] Suppose that the curvature of $\C$ satisfies inequality . If $\varphi=\Phi(s)$ is a VHC such that $\Phi(s) \in
(0,\pi/2)$ for all $ s \in [\Re]_T$, then the feedback $$\label{eq:stabilizing_feedback}
\begin{aligned}
u& =\frac{1}{\Delta(q)}\Big( \frac 1 h g {s_\varphi}- \Big( \Phi'' + \frac 1
h\Big( (1-h \kappa {s_\varphi}) \kappa \|\sigma'\|\\
&+ b \kappa'+\frac{b \kappa {\sigma'}\trans \sigma''}{\|\sigma'\|^2}\Big) {c_\varphi}\|\sigma'\|\Big) \dot s^2 + K_1 e + K_2 \dot e\Big),
\end{aligned}$$ where $\Delta(q) = \Phi'(s) + a(s) {c_\varphi}$, $e= \varphi - \Phi(s)$, $\dot e =\dot \varphi - \Phi'(s)\dot s$, and $K_1$, $K_2$ are positive design parameters, solves the maneuvering problem and has the following properties:

(i) The constraint manifold $\Gamma$ is invariant and locally exponentially stable for the closed-loop system , .

(ii) There exists a $C^1$ function $\nu: [\Re]_T
      \to (0,\infty)$ such that the closed orbit $\bar \cR =
      \{(q,\dot q) \in \Gamma: \dot s=\nu(s)\}$ is asymptotically stable for the closed-loop system, and its domain of attraction is a neighbourhood of the set $\{(q,\dot q) \in \Gamma : \dot s>0\}$.

(iii) For initial conditions in the domain of attraction of $\bar
      \cR$, the rear wheel of the bicycle traverses the entire curve $\C$, and its speed is periodic in steady-state.

[Control law (\[eq:stabilizing\_feedback\]) is [an input-output linearizing feedback for system  with output $H(\varphi,s)=\varphi-\Phi(s)$. Letting $e= \varphi - \Phi(s)$, the closed-loop system satisfies $\ddot{e}=-K_1 e - K_2 \dot e$, and it exponentially stabilizes $\Gamma = \{(q,\dot q): e=\dot
    e=0\}$. Substituting $e=\dot e=0$ in (\[eq:stabilizing\_feedback\]) we recover the feedback in (\[eqn\_inv\_feedback\]) rendering $\Gamma$ invariant.]{}]{}

The map $(\varphi,\dot \varphi,s,\dot s) \mapsto (e,\dot e,s,\dot s)$ is a diffeomorphism, and the image of $\Gamma$ under this map is the set $\tilde \Gamma :=\{(e,\dot e,s,\dot s): e=\dot e=0\}$. The feedback  is smooth in a neighbourhood of $\Gamma$, and it yields $\ddot e + K_1 e + K_2 =0$, with $K_1,K_2>0$. Therefore, the set $\tilde \Gamma$ is exponentially stable, implying that $\Gamma$ is exponentially stable as well, and proving part (i). As for part (ii), by Propositions \[prop:internal\_dynamics\_general\_properties\] and \[prop:internal\_dynamics\_bicycle\], there exists a smooth function $\nu:[\Re]_T \to (0,\infty)$ such that the closed orbit $\cR
= \{(s,\dot s): \dot s = \nu(s)\}$ is exponentially stable for system . Therefore, $\bar \cR$ is asymptotically stable relative to $\Gamma$ (i.e., asymptotically stable when the initial conditions are restricted to lie on $\Gamma$). In order to prove that $\bar \cR$ is asymptotically stable for initial conditions in a neighbourhood of $\Gamma$, note that $\bar
\cR$ is a closed curve, and hence a compact set. Owing to the reduction principle for asymptotic stability of compact sets (see [@SeiFlo95] and [@ElHMag13 Theorem 10]), the asymptotic stability of $\bar \cR$ relative to $\Gamma$, together with the asymptotic stability of $\Gamma$, imply that $\bar \cR$ is asymptotically stable for . By Propositions \[prop:internal\_dynamics\_general\_properties\] and \[prop:internal\_dynamics\_bicycle\], its domain of attraction contains the set $\{(q,\dot q) \in \Gamma : \dot s>0\}$. Since the domain of attraction of a closed set is an open set, the domain of attraction of $\Gamma$ is a neighbourhood of $\{(q,\dot q) \in \Gamma
: \dot s>0\}$. Finally, concerning part (iii), for all $(q,\dot q)
\in\bar \cR$ we have $\dot s = \nu(s)>0$. Hence, for all initial conditions in the domain of attraction of $\bar \cR$ there exist $\bar
t>0$ and $\epsilon >0$ such that $\dot s(t)>\epsilon>0$ for all $ t
\geq \bar t$, and hence the bicycle traverses the entire curve $\C$. Since $\bar \cR$ is diffeomorphic to $S^1$, since it is asymptotically stable, and solutions originating on it are periodic, $\bar \cR$ is a stable limit cycle of the closed-loop system. Therefore, solutions in the domain of attraction of $\bar \cR$ converge asymptotically to a periodic orbit.

Numerical implementation {#sec:numerics}
========================

To implement the solution presented in Theorem \[thm:solution\] one needs an analytical expression of the function $\Phi: [\Re]_T \to
\Re$. This function is found through numerical integration of the virtual constraint generator  using a numerical procedure. We will now outline this procedure and make informal deductions about the impact of the approximation. We begin our design by picking a function $\tilde \mu:\Re \to (0,\infty)$ according to Proposition \[prop:generator:periodic\_solutions\], part (ii), and an initial condition $x_0$. Using a one-dimensional search, we find the unique value of $\epsilon$ such that the solution of  is $T$-periodic with a desired accuracy. The function $\Phi:[\Re]_T \to \Re$ is determined as a spline interpolation from the samples of the solution of  over a period. With a spline interpolation at hand, one can compute $\Phi', \Phi''$ analytically and implement . The numerical approximation process used to find $\Phi$ introduces a bounded error in the feedback controller  which can be made arbitrarily small. The effect of this error is to perturb the constraint manifold $\Gamma$ without affecting its exponential stability. The dynamics on the perturbed manifold are still governed by , where the functions $\Psi_1,\Psi_2$ are affected by small perturbations. Such perturbations have no effect on the hypotheses of Proposition \[prop:internal\_dynamics\_general\_properties\], because they involve strict inequalities, $\Psi_1 >0$, $\int_0^T \tilde \Psi_2
d \tau <0$. Thus, the conclusion of Theorem \[thm:solution\] still holds. The approximation of $\Phi$ perturbs the constraint manifold $\Gamma$ and the asymptotically stable limit cycle $\bar \cR$.

We return to example of the ellipse, with the VHC depicted in Figure \[fig:ellipse:virtual\_constraint\]. The simulation results for the closed-loop system with controller  and $K_1=5$, $K_2=2$ are shown in Figures \[fig:ellipse:closed\_loop2\], \[fig:ellipse:closed\_loop1\] for the initial condition $(\varphi(0),\dot \varphi(0),s(0),\dot
  s(0)) = (0.1,0,0,0.2)$. Figure \[fig:ellipse:closed\_loop2\] illustrates the exponential convergence of $\varphi(t)- \Phi(s(t))$ to zero. Figure \[fig:ellipse:closed\_loop1\] displays the projection of the phase curve on the $(s,\dot s)$ plane and its convergence to the submanifold $\cR$.

![Simulation of the closed-loop system for the ellipse example. The solution converges to the constraint manifold $\Gamma$.[]{data-label="fig:ellipse:closed_loop2"}](FIGURES/ellipse_closed_loop2bis){width=".5\textwidth"}

![Simulation of the closed-loop system for the ellipse example. The solution converges to the submanifold $\cR \subset \Gamma$.[]{data-label="fig:ellipse:closed_loop1"}](FIGURES/ellipse_closed_loop1){width=".5\textwidth"}

\[example\_non\_convex\]

In this second example, we consider the $2 \pi$-periodic curve shown in Figure \[fig:non\_convex\], parametrized as $(5+1.5 \cos
  (2 t)) (\cos t, \sin t)$, and then reparameterized with respect to the arc-length $s$. The curve has length $L=39.129$. Since $\sigma:[\Re]_L \to \Re^2$ is not convex, Assumption \[assump\_convex\] is not satisfied and the results presented in the paper cannot be applied directly. Nevertheless, the maneuvering problem can still be solved applying the VHC method. For the virtual constraint generator (\[eq:constraint\_generator\]), we pick $\tilde \delta(t)\equiv \epsilon$, and choose $x(0)=0.35$. Numerically, we find that setting $\epsilon=0.1194$ the solution of  is $L$-periodic, thus giving rise to a valid VHC.

![The non-convex curve considered in example \[example\_non\_convex\].[]{data-label="fig:non_convex"}](FIGURES/non_convex){width=".5\textwidth"}

Again, the simulation results for the closed-loop system with controller  and $K_1=5$, $K_2=2$ are shown in Figures \[fig:non\_convex:closed\_loop2\], \[fig:non\_convex:closed\_loop1\] for the initial condition $(\varphi(0),\dot \varphi(0),s(0),\dot s(0))
= (0,0,0,2)$.

![Simulation of the closed-loop system for the curve considered in example \[example\_non\_convex\]. The solution converges to the constraint manifold $\Gamma$.[]{data-label="fig:non_convex:closed_loop2"}](FIGURES/non_convex_closed_loop2bis){width=".5\textwidth"}

![Simulation of the closed-loop system for the curve considered in example \[example\_non\_convex\]. The solution converges to the submanifold $\cR \subset \Gamma$.[]{data-label="fig:non_convex:closed_loop1"}](FIGURES/non_convex_closed_loop1bis){width=".5\textwidth"}

Acknowledgements {#acknowledgements .unnumbered}
================

We are grateful to Alireza Mohammadi for his helpful comments on this paper.

[^1]: L. Consolini is with the Department of Information Engineering, University of Parma, Viale Usberti 181/A, Parma, 43124 Italy. E-mail: [lucac@ce.unipr.it]{}. M. Maggiore is with the Department of Electrical and Computer Engineering, University of Toronto, 10 King’s College Road, Toronto, Ontario, M5S 3G4, Canada. E-mail: [maggiore@control.utoronto.ca]{}.

[^2]: An integral of motion of , or first integral, is a real-valued function of the state whose value is constant along all solutions of  (see, e.g., [@Har02]). The function $I(s,\dot s)$ is only constant along the solution of  with initial condition $(s(0),\dot
        s(0)) = (s_0,\dot s_0)$. It is not constant along other solutions, because the Lie derivative of the function $I(s,\dot s)$ along the vector field of  is not zero. In our previous work [@MagCon13] we have shown that if $M$ and $V$ are well-defined, then the correct integral of motion of  is $E(s,\dot s) = (1/2) M(s) \dot
        s^2 + V(s)$.

[^3]: This substitution is standard. See, for instance, [@Pol03 Section 2.9.3-2, item 25]. It is also used in [@Shiriaev2006900].
