# Feedback Equivalence and Time Minimal Geodesics for Controlled Euler Equation with a Single Actuator and Attitude Control

Bernard Bonnard, Cristian Mendico, Jérémy Rouot

► To cite this version:

Bernard Bonnard, Cristian Mendico, Jérémy Rouot. Feedback Equivalence and Time Minimal Geodesics for Controlled Euler Equation with a Single Actuator and Attitude Control. 2025. hal-05382034

HAL Id: hal-05382034

https://inria.hal.science/hal-05382034v1

Preprint submitted on 25 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Feedback Equivalence and Time Minimal Geodesics for Controlled Euler Equation with a Single Actuator and Attitude Control

Bernard Bonnard<fnref n="1" />, Cristian Mendico<fnref n="2" />, Jérémy Rouot<fnref n="3" />*

*Corresponding author(s). E-mail(s): jeremy.rouot@univ-brest.fr;
Contributing authors: bernard.bonnard@u-bourgogne.fr;
cristian.mendico@u-bourgogne.fr;

## Abstract

In this article we study the attitude transfer of a rigid spacecraft controlled by a single torque using time minimal control. The geodesics are parametrized by Pontryagin's maximum principle. We concentrate mainly our study to analyze the geodesics of the angular speed dynamics which are classified depending upon the position of the torque on the satellite. The important role in the study is related to the so-called abnormal geodesics which are branches of a ternary cubic which is analyzed in details in the frame of a real classification of cubic using symbolic and numerical computations. It is applied to describe attitude control techniques based on the construction of mechanical gates which in the future will be tested on low cost CubeSat. They can be interpreted as perturbations of the periodic or quasi-periodic solutions of the Euler-Poisot (uncontrolled) dynamics in relation with KAM theory.

**Keywords:** Attitude control, Pontryagin's maximum principle, Abnormal geodesics, Hamiltonian dynamics, KAM theorem

**MSC Classification:** 49K15, 70E60, 14H52

<footnote n="1">¹Institut de Mathématique de Bourgogne, UMR 5584 CNRS, Université Bourgogne Franche-Comté; INRIA, McTAO, Sophia Antipolis, France.</footnote>

<footnote n="2">²Institut de Mathématique de Bourgogne, UMR 5584 CNRS, Université Bourgogne Franche-Comté, France.</footnote>

<footnote n="3">³*Univ Brest, UMR CNRS 6205, Laboratoire de Mathématiques de Bretagne Atlantique, Brest, F-29200, France.</footnote>
# 1 Introduction

At the end of the 1970s a scientific project supported by E.S.T.E.C was launched with two main objectives: first to investigate the applications of the emerging geometric control theory in the attitude control problem using the geometric interpretation of Euler-Poinsot dynamics and, secondly, apply this analysis to understand the controllability properties in the case of under-actuated control systems due to failure of one or more gas jets or momentum wheel actuators. In addition to the confidential report delivered to the European agency [17], the results were mainly published later in two independent references [8, 16]. In particular it was proved that in general a single torque can be used to control the attitude and based on this a pertinent project is to analyze the time minimal control attitude problem using an under-actuated control system-even in the limit case of a single torque. Some revival of this study is now motivated by attitude control of low cost recent CubeSat projects started in 1999 and developed in different public or private space agencies since complex control laws can be tested in this new framework, and the destruction of the satellite causes no environmental damage.

The origin of this geometric study came from many contributors, see in particular, the seminal work by V. Arnold [1, 2], who gave a neat geometric presentation of the analysis of the inertial motion of a rigid body, revisiting the Euler analysis which has produced the so-called Euler-Poinsot geometric representation of the rigid body motion in classical mechanics [29]. The motion can be visualized making rotating without slipping on a plane the inertia ellipsoid. In particular the Euler-Poinsot dynamics is almost everywhere Poisson stable. This property allows to deduce general controllability properties thanks to [7]. More precisely, Euler-Poinsot computations amount to fully integrating the dynamics defining either periodic motion or quasi-periodic motions dense on two-dimensional tori which can be computed using elliptic functions. The Euler dynamics by itself led to polhode periodic motions which in fine represent $\wp$-Weierstrass functions real branches.

Such Poisson stability property of the uncontrolled dynamics can be generalized to any compact Lie group and in the non compact case to any semi-simple Lie group using the Kirillov-Konstant-Souriau theorem on action of the co-adjoint representation, which is also useful to compute first integrals of the Euler dynamics. This geometric framework comes from the general setting of left-invariant dynamics on Lie groups and this can be explained in the case of matrices Lie groups.

## Presentation of the system

Next, we use [13] for a short representation. Euler-Poinsot rigid body motions can be interpreted as the Euler-Lagrange equation for left-invariant Riemannian problems on the group of rotations $SO(3)$, that is,

$$ \min_{u(\cdot) \in L^2[0,T]} \int_0^T \sum_{i=1}^3 I_i u_i(t)^2 dt \qquad (1) $$
$$ \text{subject to} \quad \frac{dR}{dt} = \sum_{i=1}^{3} u_i R_i A_i, \qquad (2) $$

where $u(\cdot)$ is an accessory control parameter and

$$ A_1 = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & -1 \\ 0 & 1 & 0 \end{pmatrix}, \quad A_2 = \begin{pmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ -1 & 0 & 0 \end{pmatrix}, \quad A_3 = \begin{pmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} $$

with $RA_i$, for $i = 1, 2, 3$, forming a basis of left-invariant vector fields on $SO(3)$. Moreover, we assume the principal moments of inertia of the body to be ordered following the usual convention $I_1 \ge I_2 \ge I_3 > 0$.

Using [1, 2] the geodesics curves are solutions of the left-invariant Hamiltonian

$$ H = \frac{1}{2} \left( \frac{M_1^2}{I_1} + \frac{M_2^2}{I_2} + \frac{M_3^2}{I_3} \right) $$

where $M = (M_1, M_2, M_3)^T$ is the angular momentum of the body measured in a body-moving frame, denoted by $K(t)$, with principal axis with basis $(E_1, E_2, E_3)$, describing the geometric shape of the body. The attitude of the body is given by the matrix $R$ associated with the change of coordinates from a fixed Galilean frame $k := (e_1, e_2, e_3)$ to the moving frame.

This is the starting point of many generalizations concerning controlled dynamics on Lie groups which are detailed in the classical reference textbook [20]. Assume for simplicity that the control devices are opposing gas jets—the case of momentum wheels actuators can be treated similarly. One obtains a control system on an arbitrarily $n$-dimensional linear Lie group $G$ of the form

$$ \frac{dR}{dt} = \sum_{i=1}^{n} \Omega_i R A_i \qquad (3) $$

where $R \in G$ and $X_i = RA_i$ are a basis of left-invariant vector fields and $\Omega = (\Omega_1, ..., \Omega_n)^T$ represents the angular velocity measured in the moving frame satisfying the controlled dynamics defined by

$$ \frac{d\Omega}{dt} = Q(\Omega) + \sum_{i=1}^{p} u_i(t)b^i $$

with $b^i$ constants and $Q$ is the quadratic (homogeneous) vector field defined for $i = 1, ..., n$ by the $n$-equations

$$ I_i \frac{d\Omega_i}{dt} = \sum_{j,k=1}^{n} I_j c_{ij}^k \Omega_j \Omega_k. \qquad (4) $$

where $c_{ij}^k$ are given by

$$ [X_i, X_j] = \sum_{k=1}^{n} c_{ij}^{k} X_{k}. \qquad (5) $$
The left-invariant metric associated to the free motion is the kinetic energy

$$L \left( \sum_{i=1}^{n} \Omega X_i \right) = \frac{1}{2} \sum_{i=1}^{n} I_i \Omega_i^2. \qquad (6)$$

with $I_i > 0$.

## Some properties

To understand the Euler dynamics some deep properties come from Lie group representation theory, we refer to [1].

It is convenient to restrict on semi-simple Lie group for which the adjoint and co-adjoint representations are equivalent. The Euler dynamics is transported from the tangent space $T_gG$ to the cotangent space $T_g^*G$ and the Euler equation (4) is an Hamiltonian dynamics restricting to the orbits of the co-adjoint representation of $G$. This generalizes the $SO(3)$ case where orbits are constant momentum norm spheres.

Additional symmetries come from specific properties of the metric. The case $I_i = I_j$ for $i \neq j$ is the called the case of revolution symmetry while the case $I_i = I_j$ for all $i, j = 1, \dots, n$ is called the round or bi-invariant situation since it corresponds to a bi-invariant metric.

This leads in particular to a complete description of the geodesic dynamics on $SO(3)$ given by the next proposition.

### Proposition 1

1. On $SO(3)$, with the restriction $\| m \| = \| M \| = 1$, where $m, M$ are the momenta measured respectively in the fixed frame $k$ and the moving frame $K(t)$—hence satisfying $m = RM$— the dynamics is Hamiltonian on the orbits of the co-adjoint representation, with Hamiltonian given by the kinetic energy

$$L = \frac{1}{2} \sum_{i=1}^{3} \frac{M_i^2}{I_i}.$$

Equilibria from which the angular velocity $\Omega(t)$ is constant correspond to stationary rotations along the axes of the reference frame. If moreover $I_1 > I_2 > I_3$ the two rotations around $E_1$ and $E_3$ are stable while the rotation around the intermediate axis $E_2$ is unstable.

2. The trajectories are quasi-periodic except for the separatrices and are defined by two frequencies. Moreover, the motion is either periodic or dense on a two dimensional torus.

Such properties were used in the earliest contributions [16], [8] to design attitude transfer with two torques in the case of failure of one actuator, based for instance on
Euler angle representation [29], where the attitude $R \in SO(3)$ is decomposed as

$$R = \exp(\varphi A_1) \exp(\psi A_3) \exp(\varphi \theta A_1),$$

this product is called a mechanical gate, where we choose the two stable stationary rotations $A_1, A_3$ in the decomposition. Designs of attitude transfer are based on local accessibility properties of the angular velocity dynamics and this led to piecewise constant feedback [16].

In the case of a single torque the control problem is much more complex since the angular velocity is no longer locally controllable at zero and Poisson stability of the dynamics has to be used to design control laws and led to formulate the problem as a time-minimizing problem, where the time transfer depends on the quasi-periodicity of the uncontrolled motion. Time minimizing trajectories can be parameterized using Pontryagin's maximum principle and the main contribution of this article is to show that we can construct mechanical gate using the one associated to algebraic periodic geodesics which are precisely the so-called abnormal geodesics.

Finally this will lead to a KAM and weak KAM type theorem where the periodic and two-dimensional tori of the uncontrolled Euler-Poinsot motion will be preserved for a family of abnormal geodesics which are quasi-periodic and depending upon two frequencies. They will be interpreted as Euler-Poinsot motion where the metric corresponds to a reshaping of the geometry of the satellite defined by the $I_i$ by a rational abnormal singular feedback.

## Organization of the paper.

1. The second section studies the problem of classification of real ternary cubics, this in order to analyze the abnormal geodesics for the controlled Euler dynamics. In particular an algorithm called the chord and tangent method based on [18] is described to calculate the abnormal geodesics. It is applied to geodesic dynamics in the case of quadratic control systems.

2. The third section presents a detailed analysis of the geodesics dynamics in the case of controlled Euler dynamics. A bifurcation diagram is outlined where the bifurcating parameters are the position of the torque on the satellite. This study is briefly applied to control the full attitude dynamics.

3. The fourth section describes the $\alpha$ and $\omega$ limit sets for geodesics of single-input homogeneous systems in dimension 3, distinguishing the even from the odd case.

The article is complemented by two appendices. The second appendix presents symbolic computations and numerical simulations in the case of the controlled Euler dynamics.
# 2 Classification of quadratic systems and the controlled Euler equation

## 2.1 Quadratic differential equation and linear classification of commutative and non-associative algebras

Consider a quadratic differential equation in $\mathbb{R}^n$:

$$ \frac{dx_i}{dt}(t) = \sum_{1 \le j, k \le n} a_{jk}^i x_j(t) x_k(t) \qquad (7) $$

where the $n^3$ coefficients satisfy the conditions $a_{jk}^i = a_{kj}^i$. Such an equation can be identified to a (1, 2) tensor and a linear change of coordinates $x = Py$ induces the following action

$$ Q \to P * Q := P^{-1}Q \circ P. $$

It is a specific case of homogeneous polynomial differential equations

$$ \frac{d}{dt}x(t) = H(x(t)), \quad \text{with} \quad H(\lambda x) = \lambda^k H(x), \quad k \in \mathbb{N}, $$

extending the linear case. Those equations admit in general rays solutions, that is, solutions which are straight-lines. They are crucial in the classification and it comes from the following geometric property.

Since the dynamics is homogeneous the distribution $x \to \mathbb{R}H(x)$ is invariant with respect to homotheties and it can be projected on the sphere $S^{n-1}$. More precisely, if we set

$$ r = \sqrt{\sum_{i=1}^{n} x_i^2}, \quad \text{and} \quad x' = x/r, $$

then the dynamics has the form

$$ \frac{dr}{dt} = \langle x', H(x') \rangle r \qquad (8) $$

$$ \frac{dx'}{d\sigma} = H(x') - \langle x', H(x') \rangle x' \qquad (9) $$

where $\sigma$ corresponds to a time re-parameterization. Here, the equation (9) defines the projected vector field for which rays correspond to equilibria. This gives the general framework to classify quadratic dynamics. For instance, in the 2d-case rays delimit node (fan), saddle or elliptic sectors and this comes from the classification of planar degenerate singular points [26].

Additionally, in the quadratic case, L. Markus [27] has introduced the following algebra $\mathcal{A}$ in relation with the classification.

**Definition 1** Let $(e_i)$ be the $\mathbb{R}^n$ canonical basis and define the multiplication rule by
$$e_j \cdot e_k = \sum_{i=1}^{n} a_{jk}^{i} e_i. \qquad (10)$$

We immediately have that $\mathbf{A}$ is commutative and possibly non-associative. Moreover, the linear classification of quadratic differential equations is equivalent to the linear classification of such an algebra and the change of coordinates represented by the transition matrix $P$ induces the isomorphism.

Indeed, consider more generally a system defined by a quadratic mapping $E : \mathbb{R}^n \to \mathbb{R}^n$ of the form

$$E(x) = C + Tx + Q(x),$$

with $C$ constant, $T$ linear and $Q$ quadratic homogeneous map. From $Q$ we obtain an unique symmetric bilinear mapping

$$\beta : \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}^n$$

defined by

$$\beta(x, y) = \frac{1}{2}[Q(x + y) - Q(x) - Q(y)]$$

and the product defined by

$$x \cdot y = \beta(x, y).$$

The pair $\mathbf{A} = (\mathbb{R}^n, \beta)$ defines a commutative and non-associative algebra. Such an algebra is equivalent to L. Markus' definition of the algebra, since it depends only on the homogeneous quadratic form. Geometric definitions translate into algebraic concepts, as briefly presented next.

**Definition 2** A subalgebra $\mathbf{B}$ is a subspace such that $\mathbb{B}^2 \subseteq \mathbf{B}$. An ideal $\mathbf{l}$ is a subspace such that $\mathbf{l} \cdot \mathbf{A} \subseteq \mathbf{A}$. The algebra is simple if $\mathbf{A}^2 \neq 0$ and has no proper ideals. An algebra $\mathbf{A}$ is nilpotent if there exists an integer $N$ such that all products with $N$ factors are zero. The algebra $\mathbf{A}$ is solvable if $\mathbf{A} \supseteq \mathbf{A}^{(2)} \dots \supseteq \mathbf{A}^{(N)} = 0$ where $\mathbf{A}^{(k+1)}$ denotes $\mathbf{A}^{(k)} \cdot \mathbf{A}^{(k)}$.

Algebraic definitions are translated next for the dynamics.

**Proposition 2**

1. A subalgebra $\mathbf{E}$ is an invariant subspace for the solutions of the dynamics on which it can be restricted.

2. An ideal $\mathbf{I} \subset \mathbf{A}$ allows to project the dynamics on the quotient space $\mathbf{A}/\mathbf{I}$.

3. Idempotents $e^2 = e$ correspond to one dimensional subalgebra defined by $Q(\lambda e) = \lambda^2 Q(e)$ forming rays.

4. Lines of non-isolated equilibria correspond to nilpotent elements of order 2.

*Example 1* Consider the Euler field
$$Q(x) = [a_{23}^{1}x_{2}x_{3}, a_{13}^{2}x_{1}x_{3}, a_{12}^{3}x_{1}x_{2}]^{T}.$$

Then, we have that:

1. it is simple provided $a_{23}^1 a_{13}^2 a_{12}^3 \neq 0$ and the axes are the nilpotent elements of order 2;

2. the Euler algebra is solvable if $a_{23}^1 = 0$.

More results about such algebras will be given later since the controlled Euler dynamics is the core of this article.

In the general quadratic dynamics

$$\frac{dx}{dt} = E(x)$$

the above construction can be extended as a homogeneous quadratic mapping $Q'$ in the $n+1$-dimensional space with coordinates $(x, z)$ setting

$$Q'(x, z) = \begin{pmatrix} Cz^2 + zTx + Q(x) \\ 0 \end{pmatrix}.$$

### Example 2

(i) The standard Lotka-Volterra model for population dynamics, written in a specific form

$$\dot{x} = E(x) = \operatorname{diag}(x)(r + Ax)$$

where $x \in \mathbb{R}^n$ represents the population densities of $n$ different species, $r \in \mathbb{R}^n$ is the intrinsic growth rate vector, $A$ is the interaction matrix and its element $A_{ij}$ represents the effect of species $j$ on the growth rate of species $i$. Let us write it component-wise to see why it's a quadratic vector field:

$$\dot{x}_i = x_i \left( r_i + \sum_{j=1}^{n} A_{ij} x_j \right) = r_i x_i + \sum_{j=1}^{n} A_{ij} x_i x_j$$

This is clearly a quadratic function in $x$. The term $r_i x_i$ is linear, and the terms $A_{ij} x_i x_j$ are quadratic. This system corresponds to an algebra where the multiplication is defined to make $x \cdot x = E(x)$. The structure constants $a_{ij}^k$ would be derived from the coefficients of $x_i x_j$ in the expansion above.

(ii) Consider a quadratic system

$$\frac{dx}{dt} = Q(x) + bu$$

where $Q(x)$ is a homogeneous quadratic vector field, $b$ is the control input vector (it determines how the control $u$ influences the system), $u$ is the scalar control signal. We consider a control in the feedback form

$$u = \alpha(x) + \beta v$$
where $\alpha(x)$ is a quadratic function of the state $x$, $\beta$ is a constant scalar, $v$ is a new external control input. Substituting the feedback $u$ into the original system, we get

$$ \frac{dx}{dt} = Q(x) + b(\alpha(x) + \beta v) = Q(x) + \alpha(x)b + \beta vb. $$

By augmenting the state space to include $z$ such that

$$ \frac{dz}{dt} = 0 $$

we get the quadratic term $z^2\betavb$. This allows the entire controlled system to be written in the homogeneous quadratic form $\dot{x} = \tilde{Q}(\tilde{x})$, where $\tilde{x} = (x, z)$.

We refer to [22] for the computations of algebraic cycles for quadratic dynamics.

**Proposition 3** Non isolated periodic trajectories can occur for quadratic dynamics if $n \ge 3$ and they are isolated periodic trajectories for the projection of the dynamics onto sphere defined by (9). Moreover, they form in dimension 3 invariant conic sectors.

## 2.2 Lie algebraic properties of quadratic control systems [4, 8]

The following lemma relates properties of $\mathcal{A}$ to Lie brackets computations.

**Lemma 1** For every $v_1, v_2 \in \mathbb{R}^n$, we have

1. $[v_1, Q](v_2) = 2v_1 \cdot v_2$,

2. $[v_2, [v_1, Q]] = 2v_1 \cdot v_2$.

This leads to the following result from [11] and [4], showing the complexity of the study of quadratic systems.

Consider a quadratic control system

$$ \frac{dx}{dt}(t) = Q(x(t)) + Bu(t), \qquad (11) $$

let $\{b_1, ..., b_m\}$ be the columns of $B$. We denote by $\{H^i\}_{i \in \mathbb{N}}$ the set of homogeneous polynomial vector fields of $\mathbb{R}^n$ with respective degree $i$. Then, the following holds.

**Proposition 4** Denote by $\mathcal{L}$ the Lie algebra generated by the vector fields $\{Q, b_1, ..., b_m\}$ and we define:

(i) $L^0$ is the smallest subspace of $H^0$ such that: $\mathrm{span}\{b_i\} \subset L^0$. If $v_1, v_2 \in L^0$, then $[[Q, v_1], v_2] \in L^0$;

(ii) $L^1$ is the smallest Lie subalgebra of $H^1$ generated by the vector fields $\{[Q, v], v \in L^0\}$;
(iii) $L^2$ is the Lie subalgebra of $\oplus_{k \ge 2} H_k$ generated by the vector fields

$$ \{ \mathrm{ad} L_1 \circ \cdots \circ \mathrm{ad} L_p(q) \text{ such that } p \in \mathbb{N}, L_1 \in L^1 \}. $$

Then $\mathbf{L} = L^0 \oplus L^1 \oplus L^2$, $\mathbf{L}(0) = L^0$ and, moreover, $\dim \mathbf{L}(x) \ge \dim L^0$.

## 2.3 Some consequence for the controllability properties of generalized Euler equation

**Definition 3** Let $X$ be a complete vector field on a manifold $M$ and denote $\exp(tX)$ the one parameter group associated to $X$. A point $x_0$ is called Poisson stable if for every neighborhood $V$ of $x_0$ and every time $T > 0$ there exists $t_1, t_2 > T$ such that $\exp(t_1X(x_0))$, $\exp(-t_2X(x_0)) \in V$. A point $x_0$ is departing if the trajectory starting at $x_0$ has no $\omega$ and $\alpha$ limit points. The vector field $X$ is Poisson stable if almost every point is Poisson stable.

**Theorem 5** ([33, p. 454]) Let $X$ be a complete conservative vector field on a manifold $(M, \omega)$, where $\omega$ denotes the volume form. Then, almost all points $x_0$ are either Poisson stable or departing.

From the above section one can prove the following result.

**Proposition 6** Consider a quadratic control system on $\mathbb{R}^n$,

$$ \frac{dx}{dt} = Q(x) + Bu, $$

with $B = \operatorname{span}\{b_1, \dots, b_p\}$. Then, we have the following.

1. The Lie bracket $[b, [b, Q]]$ is collinear to $Q(b)$.

2. Assume $m = n - 1$, and let $H_+$ and $H_-$ be the positive and negative domains delimited by the hyperplane $B$. If $Q$ restricted to $B$ is pointing towards $H_+$ and $H_-$, then for all $t > 0$ the accessibility set $A(0, t)$ is a neighborhood of 0.

3. The Lie algebra $L$ generated by $\{Q, B\}$ is of full rank if and only $L^0$ has rank $n$. This is the case if and only $L^0$ doesn't belong to a proper subalgebra of $\mathbf{A}$.

4. If $Q$ is Poisson stable, then the system is controllable if and only $L^0$ has rank $n$.

5. If $Q$ is induced by a $G$-left invariant metric, then it is Poisson stable provided $G$ is compact or $G$ is a semi-simple.

*Proof* The first assertion follows from the ad-formulae:

$$ b * Q(0) = \sum_{i=1}^{n} \frac{1}{n!} ad^n b.Q(0) = [b, [b, Q]](0). $$

By construction the vector field $b * Q$ is induced by the translation of vector $x \to x + b$ at 0 is collinear to $Q(b)$. If $B$ is a hyperplane, then for all $t > 0$, the accessibility set in time $t$, denoted $A(0, t)$ is a neighborhood of 0 if there exists two vectors $b_1, b_2$ in $B$ so that $Q(b_1), Q(b_2)$ point respectively towards $H_+, H_-$.
From Poincaré-Hopf Theorem 5 taking a conservative vector field, then for almost every point the initial point of the dynamics is either Poisson stable or departing from every compact.

For left-invariant Riemannian metrics on Lie group the geodesic dynamics is Poisson stable if G is compact. If G is semi-simple, the Euler dynamics is Poisson stable restricting to each orbit of the co-adjoint representation, see [8]. □

## 2.4 Feedback classification for quadratic systems

Next we introduce the feedback classification for quadratic systems restricting to the finite dimensional action induced by linear change of coordinates. The appendix A gives properties for the general action on affine control systems. See [11] for the details.

Denote by $\mathcal{C}$ the set of quadratic systems on $\mathbb{R}^n$ defined as

$$ \frac{dx}{dt} = Q(x) + \sum_{i=1}^{m} u_i b_i = Q(x) + Bu, $$

with $m \le n$ and $B$ of full rank. It will be usually written in short $\{(Q, B)\}$. We consider the action of the feedback group denoted $G'_f$ restricted to triples $\{(P, \alpha(x), \beta)\}$ where $P$ is a linear mapping, $\alpha$ is quadratic and $\beta$ is constant. Recall that the linear action $x = Py$ induces the following on a quadratic vector field

$$ Q \to P^{-1} \circ Q \circ P $$

and feedback action $(\alpha, \beta)$ gives

$$ Q \to Q + B\alpha, \quad B \to B\beta. $$

The class $\mathcal{C}$ is invariant for this (finite) dimensional action and normal forms can be obtained by considering the action of the feedback group for the reduced action of linear transformations preserving the distribution $B$.

To simplify the notations one can consider the class of single-input quadratic systems

$$ \frac{dx}{dt} = Q(x) + ub. $$

where the vector $b$ can be identified to the vector $e_n = (0, 0, ..., 0, 1)^T$. The subgroup of automorphisms leaving the line $e_n$ invariant corresponds to matrices $P = (p_{ij})$ where we have $p_{ij} = 0$ for $j = n, i = 1, ..., n-1$. If $Q = (Q_1, ..., Q_n)^T$, then the subgroup acts by the action induced on the set of $n-1$ quadratic forms $(Q_1, ..., Q_{n-1})^T$.

### 2.4.1 Euler Equation and Elliptic Functions

We use [25, p. 131] using Jacobi elliptic functions, or [25, p. 137] using Weierstrass functions.

We set

$$ I_1\omega_1^2 + I_2\omega_2^2 + I_3\omega_3^2 = T $$
and

$$I_1^2 \omega_1^2 + I_2^2 \omega_2^2 + I_3^2 \omega_3^2 = H^2.$$

Assume

$$I_1 > I_2 > H^2/T > I_3$$

and define

$$P^2 = \frac{H^2 - I_3 T}{I_1(I_1 - I_3)},$$

$$Q^2 = \frac{H^2 - I_3 T}{I_2(I_2 - I_3)},$$

$$R^2 = \frac{I_1 T - H^2}{I_3(I_1 - I_3)},$$

$$p^2 = \frac{1}{I_1 I_2 I_3} (I_2 - I_3)(I_1 T - H^2),$$

$$k^2 = \frac{(I_1 - I_2)(H^2 - I_3 T)}{(I_2 - I_3)(I_1 T - H^2)}.$$

### On the Polhode motion

Using the previous notation, the solution writes

$$\omega_1 = P \operatorname{cn}[p(t - t_0)],$$

$$\omega_2 = -Q \operatorname{sn}[p(t - t_0)],$$

$$\omega_3 = R \operatorname{dn}[p(t - t_0)].$$

where cn, sn and dn stand respectively for the Jacobi cosine elliptic function, Jacobi sine elliptic function and Jacobi delta elliptic function with parameter $p$. Moreover, the limit separatrix case is when $H^2/T = B$, where the modulus $k = 1$ and the solution writes

$$\omega_1 = P \operatorname{sech}[p(t - t_0)],$$

$$\omega_2 = -Q \operatorname{tanh}[p(t - t_0)],$$

$$\omega_3 = R \operatorname{sech}[p(t - t_0)],$$

with

$$P^2 = \frac{(I_2 - I_3)T}{I_1(I_1 - I_3)},$$

$$Q^2 = \frac{T}{B},$$

$$R^2 = \frac{(I_1 - I_2)T}{I_3(I_1 - I_3)},$$

$$p^2 = \frac{1}{I_1 I_2 I_3} (I_1 - I_3)(I_2 - I_3)T.$$

One can treat similarly the other cases see [25] for the details.

(i) If $H^2/T = C$, then $P = Q$ and $k = 0$ and the solution is the stationary solution of the form $\omega_1 = \omega_2 = 0, \omega_3 = T/C$.

(ii) In case $I_1 > H^2/T > I_2 > I_3$, we have a solution of the form

$$\omega_1 = P \operatorname{dn}[p(t - t_0)],$$

$$\omega_2 = -Q \operatorname{sn}[p(t - t_0)],$$

$$\omega_3 = R \operatorname{cn}[p(t - t_0)],$$

$$P^2 = \frac{H^2 - I_3 T}{I_1(I_1 - I_3)},$$

$$Q^2 = \frac{I_1 T - H^2}{I_2(I_1 - I_3)},$$

$$R^2 = \frac{I_1 T - H^2}{I_3(I_1 - I_3)},$$
$$p^2 = \frac{1}{I_1 I_2 I_3} (I_1 - I_2)(H^2 - I_3 T), \quad k^2 = \frac{(I_2 - I_3)(I_1 T - H^2)}{(I_1 - I_2)(H^2 - I_3 T)}$$

(iii) The cases $H^2 = I_1 T$, $H^2 = I_2 T$ are treated as before.

### On the Herpolhode motion

One can use the coordinates from [20, p. 436], arising from the geometric integration of Euler-Poinsot dynamics using the Hamiltonian formalism. The last one also adapted to conjugate points computations in Euler-Poinsot uncontrolled dynamics vs action-angle variables (Serret-Andoyer variables) which mix state and adjoint variable.

The Euler Poinsot dynamics is integrable using

$$H = \frac{1}{2} \left( \frac{H_1^2}{I_1} + \frac{H_2^2}{I_2} + \frac{H_3^2}{I_3} \right)$$

and the Hamiltonian lifts of the right-invariant vectors fields $A_i R$. Such a construction is equivalent to the momentum $m$ and kinetic energy first integrals used in [25]. Note that for geometric reason, for the momentum in the moving frame we switch the notation according to

$$M_i \rightarrow H_i.$$

The Euler equation is written using

$$H_i = I_i \Omega_i,$$

where $\Omega = (\Omega_1, \Omega_2, \Omega_3)^T$ is the angular velocity in the moving frame. We use the Euler angles as a chart to represent a rotation matrix with the convention

$$R = \exp(\Phi_1 A_3) \circ \exp(\Phi_2 A_2) \circ \exp(\Phi_3 A_3),$$

where the angles $\Phi_2, \Phi_3$ are given by

$$H_1 = -|H| \sin \Phi_2 \cos \Phi_3,$$

$$H_2 = -|H| \sin \Phi_2 \sin \Phi_3,$$

$$H_3 = |H| \cos \Phi_2$$

while the remaining angles are obtained integrating the equation

$$\frac{d\Phi_1}{dt} = \frac{|H|}{H_2 \sin \Phi_3 - H_1 \cos \Phi_3} \left( \sin \Phi_3 \frac{\partial H}{\partial H_2} - \cos \Phi_3 \frac{\partial H}{\partial H_1} \right).$$

Using this chart, the Hamiltonian is

$$H_n = \frac{1}{2I_1} \left( p_2 \sin \Phi_3 - \frac{\cos \Phi_3}{\sin \Phi_2} (p_1 - p_3 \cos \Phi_2) \right)^2$$
$$ + \frac{1}{2I_2} \left( p_2 \cos \Phi_3 + \frac{\sin \Phi_3}{\sin \Phi_2} (p_1 - p_3 \cos \Phi_2) \right)^2 + \frac{1}{2I_3} p_3^2, $$

where $p_i$, for $i = 1, 2, 3$ are the dual variables of the angles $\Phi_i$, for $i = 1, 2, 3$. The left-invariant Riemannian metric is denoted $g$ and is defined by the Hamiltonian.

The dynamics controlled with gas jets is obtained making the substitution

$$ b \rightarrow (I_i b_1, I_2 b_2, I_3 b_3)^T $$

to get

$$ \frac{dH_1}{dt} = H_2 H_3 \left( \frac{1}{I_3} - \frac{1}{I_2} \right) + u I_1 b_1 $$

$$ \frac{dH_2}{dt} = H_2 H_3 \left( \frac{1}{I_1} - \frac{1}{I_3} \right) + u I_2 b_2 $$

$$ \frac{dH_3}{dt} = H_1 H_2 \left( \frac{1}{I_3} - \frac{1}{I_2} \right) + u I_3 b_3. $$

**Theorem 7 ([32]) Closed herpolhode curves are given by an explicit formulae of the form**

$$ \psi(K(k)) = 2\pi\lambda, \quad \lambda \in \mathbb{N} $$

where the complete elliptic integral $K(k)$ is associated to the evolution of the precession Euler angle of the rigid body, in Euler angles representation.

### 2.4.2 The Euler controlled dynamics

To avoid the indices, we consider

$$ \begin{cases} \frac{dx}{dt} = a_1 y z, \\ \frac{dy}{dt} = a_2 x y, \\ \frac{dz}{dt} = a_3 x y, \end{cases} $$

where we have:

$$ a_1 = \frac{I_2 - I_3}{I_1}, \quad a_2 = \frac{I_3 - I_1}{I_2}, \quad a_3 = \frac{I_1 - I_2}{I_3}. $$

They are some resonant relations between the coefficients since the system admits two first integrals (the kinetic energy and the norm of the momentum), i.e.,

$$ E = I_1 x^2 + I_2 y^2 + I_3 z^2 $$

$$ M = E = I_1^2 x^2 + I_2^2 y^2 + I_3^2 z^2. $$

Hence, we have the resonant relations

$$ I \cdot a = 0, \quad I^2 \cdot a = 0. \qquad (12) $$
The system can be written as

$$ \frac{dv}{dt} = Q(v) + ub $$

and we add the attitude dynamics

$$ \frac{dR}{dt} = \sum_{i=1}^{3} \omega_i A_i R. $$

This led to a quadratic system embedded into $\mathbb{R}^{12}$ which has to be analyzed near $R = I_{12}$, $I_{12}$ denotes the identity in $M_{12}(\mathbb{R})$, since it is left-invariant.

**Theorem 8** The system is controllable provides $a_1b_3^2 - a_3b_1^2 \neq 0$ and at least two of the coefficients $b_i$ are nonzero.

*Proof* The proof and many generalizations follow from [8] since the dynamics is (for almost every initial condition) quasi-periodic and hence Poisson stable. □

# 3 Classification of Ternary cubic and the chord and tangent method. Application to Euler equation

## 3.1 Classification ternary cubic

We use the following notation for the equivalence class of a vector in the projective space:

$$ [v] = (x : y : z) $$

In the frame of geometric invariant theory, our presentation is based on [36] as a reference supplemented by [19]. See also [23] in relation with E. Cartan's method and the useful connection with real singularity theory [3]. Furthermore, a Maple package is associated to [23] to classify the cubic over $\mathbb{C}$.

Let $X_3 = \mathbf{A}^{10}$ the space of all ternary cubic

$$ U(x, y, z) = ax^3 + by^3 + cz^3 + a_2x^2y + a_3x^2z + b_1xy^2 + b_3y^2z + c_1xz^2 + c_2y^2 + mxyz. $$

So, the coordinate ring of $X_3$ is the polynomial ring

$$ K[a, b, c, a_2, a_3, b_1, b_3, c_1, c_2, m], $$

where $K = \mathbb{R}$ or $\mathbb{C}$ is the field of reference.

The Lie groups $Gl_3(K)$ and $Sl_3(K)$ act linearly by change of coordinates. Choosing the point $O = (0, 1, 0)$ as the base point on the line at infinity, after scaling it can be reduced to
$$y^2 z + a_1 xyz + a_3 yz^2 = x^3 + a_2 x^2 z + a_4 x z^2 + a_6 z^3.$$

Using non-homogeneous coordinates

$$x = \frac{X}{Z}, \quad y = \frac{Y}{Z}$$

it takes the form

$$E : y^2 + a_1 xy + a_3 y = x^3 + a_2 x^2 + a_4 x + a_6.$$

Further reduction gives an equation of the form

$$E : y^2 = 4x_3 + b_2 x^2 + 2b_4 x + b_6$$

where

$$b_2 = a_1^2 + 4a_2, \quad b_4 = 2a_4 + a_1 a_3, \quad b_6 = a_3^2 + 4a_6.$$

Define the quantities

$$b_8 = a_1^2 a_6 + 4 a_2 a_6 - a_1 a_3 a_4 + a_2 a_3^2 - a_4^2,$$

$$c_4 = b_2^2 - 24b_4,$$

$$c_6 = -b_2^3 + 36b_2 b_4 - 216b_6,$$

$$\Delta = -b_2^2 b_8 - 8b_4^2 - 27b_6^2 + 9b_2 b_4 b_6,$$

$$j = \frac{c_4^3}{\Delta},$$

$$\omega = dx/(2y + a_1 x + a_3) = dy/(3x^2 + 2a_2 x + a_4 - a_1 y).$$

One can verify that they satisfies the relations

$$4b_8 = b_2b_6 - b_4^2, \quad 1728\Delta = c_4^3 - c_6^2.$$

and, moreover, a further reduction gives the simpler equation

$$E : y^2 = x^3 - 27c_4x - 54c_6.$$

**Definition 4** Two cubics $U_1, U_2$ are properly projectively equivalent if there exist $g \in \text{Gl}_3(K)$ such that $U_2 = (\det g)^{-1}(gU_1)$.

**Definition 5** The quantity $\Delta$ is called the discriminant of the Weierstrass equation, $j$ is called the $j$-invariant of the elliptic curve $E$ and $\omega$ is the invariant differential associated to the Weierstrass equation.

**Theorem 9**
(i) The ring of invariants is a polynomial ring in two variables generated by $c_4$ and $c_6$.

(ii) A ternary cubic $U$ is non singular if $\Delta(U) \neq 0$.

(iii) The Hessian normalized as

$$H(x, y, z) = (-1/2) \det \begin{bmatrix} \frac{\partial^2 U}{\partial x^2} & \frac{\partial^2 U}{\partial x \partial y} & \frac{\partial^2 U}{\partial x \partial z} \\ \frac{\partial^2 U}{\partial x \partial y} & \frac{\partial^2 U}{\partial y^2} & \frac{\partial^2 U}{\partial y \partial z} \\ \frac{\partial^2 U}{\partial x \partial z} & \frac{\partial^2 U}{\partial y \partial z} & \frac{\partial^2 U}{\partial z^2} \end{bmatrix},$$

is a cubic covariant.

(iv) The weight of $c_4, c_6, \Delta, H$ are respectively 4, 6, 12, 2.

**Theorem 10** In the irreducible case two non singular cubic are equivalent on $\mathbb{C}$ if and only if they have the same invariant $c_4$ and $c_6$.

It is easy to complete the classification of the real locus when the coefficients are real. In particular in our study it is convenient to use the standard unfolding of a real cusp.

**Proposition 11** Real cubic curve with a singular point at zero is either one of the following.

(i) A real node where the origin is an ordinary point represented by $x^3 - x^2 + y^2 = 0$ with two tangents $x \pm y = 0$.

(ii) A complex node represented by $x^3 + x^2 + y^2 = 0$ with two complex tangents at 0: $x \pm iy = 0$ which decomposes on $\mathbb{R}$ onto the origin and a smooth branch.

(iii) The curve called a cusp represented by $x^3 - y^2 = 0$.

(iv) A node is given by $\Delta = c_4 \neq 0$, a cusp is given by $\Delta = c_4 = 0$.

(v) The three cases can be represented by the unfolding $x^3 + \lambda x^2 + y^2 = 0$, where $\lambda$ is a real parameter.

**Theorem 12** A reducible cubic is equivalent over $\mathbb{C}$ to one of the following form.

1. If it is the product of linear and quadratic factors then it takes the form

$$z(x^2 + yz), \quad z(x^2 + y^2 + z^2).$$

2. If it is the product of three linear factors then

(i) If the three factors are linearly independent we have: $xyz$

(ii) If three are linearly dependent but any pairs of them are linearly independent we have: $xy(x + y)$

(iii) Two factors are the same, we have: $x^2y$

(iv) Three factors are the same, we have: $x^3$.

**Remark 1** The real classification exists in the literature [5] but led to more cases.
It is useful to generate feedback invariant e.g. in relation of the bi-ratio of $L_1, L_2, L_3, L_4$. Let $(a : b : c : d)$ be 4 points in the complex projective space. The bi-ratio is

$$ \frac{a - c b - d}{a - d b - c} $$

and depends of the ordering of the points. The $j$-invariant does not depend on the order. For instance, taking $(m, 1, 0, \infty)$ one has

$$ j(m) = \frac{4 (1 - m + m^2)^3}{27 m^2 (1 - m)}. $$

Finally, we have the following result concerning the classification over $\mathbb{Q}$.

**Theorem 13** ([6]) In the non singular case every cubic can be written in the normal form

$$ E : y^2 = x^3 + Ax + B, $$

where $A$ and $B$ belong to $\mathbb{Z}$.

**Proposition 14** Consider a cubic equation normalized as

$$ y^2 = 4x^3 - g_2x - g_3. $$

The $\wp$-function satisfies the differential equation

$$ (\wp')^2 = 4\wp^3 - g_2\wp - g_3. $$

Furthermore, the solutions in homogeneous coordinates are $(\wp(z) : \wp'(z) : 1)$ and the real branches correspond to the intersection with the real plane.

## 3.2 Transforming a cubic to Weierstrass equation

We use [18] where a Sage implementation of the algorithm is provided.

Consider a general cubic equation $E$ with rational coefficients and take a rational point $P$ on this cubic.

**Step 1.** If $P$ is a flex the algorithm terminates as in [18].

**Step 2** If $P$ is not a flex (generic case). let $Q$ the intersection of the tangent line to the cubic $E$ with the cubic and observe that $Q$ is rational. We compute the tangent to the cubic at $Q$ and we denote by $R$ the intersection of this tangent with the cubic $E$. The triangle $PQR$ is the golden triangle. In such coordinates the cubic can be reduced to

$$ UV^2 + aU^2W + bUVW + cUW^2 + dVW^2. $$

Introduce the quadratic Cremona transformation :

$$ \sigma : (U : V : W) \rightarrow (K^2 : LN : KN). $$
The inverse transformation is defined as:

$$ (K : L : M) \to (UW : UV : W^2), $$

where we have divided by $UW^2$. Finally after normalization, the cubic can be written

$$ L^2 M + a' K L M + b' L M^2 = K^3 + c' K^2 M, $$

and, setting $M = 1$ one gets

$$ L^2 + a'KL + b'L = K^3 + c'K^2. $$

From [36], replacing

$$ L \to \frac{1}{2}(L - a'x - b') $$

the equation can be written in Weierstrass form

$$ y^2 = P(x), $$

where $P$ is a polynomial equation of degree 3. Finally using one additional Ferrari's normalization, we end up with the equation

$$ y^2 = x^3 - 27c_4x - 54c_6, $$

where $c_4$ and $c_6$ are the invariant.

We have the following result, which extend the previous Proposition 14 to the real case.

**Proposition 15** Consider a cubic equation normalized as

$$ y^2 = 4x^3 - g_2x - g_3. $$

The $\wp$-function satisfies the differential equation

$$ (\wp')^2 = 4\wp^3 - g_2\wp - g_3, $$

and, the solutions in homogeneous coordinates are $(\wp(z) : \wp'(z) : 1)$. Moreover, the real branches correspond to the intersection with the real plane and one can restrict $z$ to the period parallelogram.

Altogether one get generically (excluding the reductive cases) one or three real branches on the two sphere $S^2$.

## 3.3 Application to Euler equations

We consider a general control system in $\mathbb{R}^n$:

$$ \frac{dx}{dt} = Q(x) + ub, $$
and we introduce for a given vector field X the Hamiltonian lift $H = p \cdot X(x)$ where $p$ denotes the adjoint vector, and $\mathbf{H}$ is the associated Hamiltonian vector field.

Assuming $u \in \mathbb{R}$, from the Maximum Principle [34] the time extremal geodesics are given by the constrained Hamiltonian dynamics with $z = (x, p)$ being the coordinates on the cotangent bundle

$$ \begin{cases} \frac{dz}{dt} = \mathbf{H}_Q(z) + u\mathbf{H}_b \\ H_b(z) = 0. \end{cases} $$

Deriving the constraints along the trajectories leads to the following relations

$$ \begin{cases} H_b = H_{[Q,b]} = 0 \\ H_{[[b,Q],Q]} + u_s H_{[[b,Q],b]} = 0 \end{cases} $$

where the Lie bracket is computed with the convention

$$ [X, Y] = D_x X \cdot Y - D_x Y \cdot X. $$

## Corollary 16

1. In the 3d-case the adjoint vector can be eliminated and this led generically to compute the geodesic control as the feedback

$$ u_s = -\frac{D'(x)}{D(x)} $$

with

$$ D = \det(b, [b, Q], [[b, Q], b]), \quad D' = \det(b, [b, Q], [[b, Q], Q]). $$

The singular dynamics can be reparameterized as

$$ \frac{dv}{d\sigma} = DQ - D'Q $$

and is cubic.

2. Denoting $D'' = \det(b, [b, Q], Q)$, the abnormal geodesics defined by the zero level set of the Hamiltonian $H_Q = 0$ are contained in the cubic set $D'' = 0$.

Computing for the Euler dynamics one gets the following result.

**Proposition 17** In the Euler dynamics one has:

$$ D'' = a_2 a_3 b_1 x^2 (b_2 z - b_3 y) + a_1 a_3 b_2 y^2 (b_3 x - b_1 z) + a_1 a_2 b_3 z^2 (b_1 y - b_2 x), $$

$$ D = 2a_1 b_2 b_3 (a_3 b_2^2 - a_2 b_3^2)x + 2a_2 b_1 b_3 (a_1 b_3^2 - a_3 b_1^2)y + 2a_3 b_1 b_2 (a_2 b_1^2 - a_1 b_2^2)z, $$

$$ D' = a_3^2 b_1 b_2 (b_1 y + b_2 x) (a_2 x^2 - a_1 y^2) + a_1^2 b_2 b_3 (b_2 z + b_3 y) (a_3 y^2 - a_2 z^2) \\ + a_2^2 b_1 b_3 (b_1 z + b_3 x) (a_1 z^2 - a_3 x^2). $$
On the projective space we set

1. $L_1 = [b]$,

2. $L_2 = [w]$, with $w := [Q, b](w)$ collinear to $b$.

3. $L_3, L_4$ the real or complex lines defined by : $L_3, L_4 \neq L_2$ and solutions of $D' = D = 0$.

Then simple computations lead to.

**Proposition 18** Assume $L_1$ and $L_2$ be distinct. Then, we have:

(i) $D = 0$ is a plane generated by $L_1$ and $L_2$;

(ii) $L_1$ and $L_2$ belong to $D'' = D = 0$ and moreover the line $e_2 = (0, 1, 0)$ belongs to the cubic.

Moreover, from a further normalization, in such a case one can choose coordinates such that $b = e_1, w = e_3$ so that $D = 0$ is the plane $y = 0$.

### 3.3.1 Weierstrass representation of the cubic $D'' = 0$.

We choose $I_1 > I_2 > I_3 > 0$ as integers, so the coefficients $a_1, a_2, a_3$ are rational and the cubic $D''$ has rational coefficients.

Assume that the cubic is regular. We apply the chord and tangent algorithm starting with the base point $P = [b]$. In the general case, this point is not a flex. The tangent space to $D'' = 0$ is then given by $D = 0$ and intersects again the cubic at $L_2$. Hence $Q = [L_2]$ in the algorithm and we can compute the third point $R$ to form the golden triangle.

*Remark 2* Such an algorithm works in general with a system $(Q, b)$ with rational coefficients.

Altogether in Euler case (or in general quadratic case) this led, in the non singular case, to the following:

1. one real oval $C_1$ containing $L_1$ and $L_2$ and intersecting the equator $D = 0$;

2. two real symmetric ovals $C_2, C_3$ located respectively in the northern and southern hemisphere.

Additionally, in the Euler case each of the lines $e_i$ given by the canonical basis are in $D'' = 0$. Further analysis is provided in the appendix B and it will lead to easily distinguishing the reducible cases from the regular cases, in terms of $b$.

# 4 More about Classification of Quadratic Systems and Euler Dynamics

## 4.1 Generalities

First, we introduce the following notation to classify the Euler dynamics $Q = (a_1yz, a_2xz, a_3xy)^T$:
* GEuler case: $a_1, a_3 > 0, a_2 < 0$

* JEuler case: $a_1 = 1, a_2 = -1, a_3 = -k^2$, with $0 < k < 1$

* NEuler case: $a_1 = \frac{I_2 - I_3}{I_1}, a_2 = \frac{I_3 - I_1}{I_2}, a_3 = \frac{I_1 - I_2}{I_3}$ with $I_1 > I_2 > I_3$

Lemma 2 (L. Markus [27])

1. Let $u = \sin t$, $v = \cos t$, $w = \operatorname{dn} t$ be the Jacobi elliptic functions with modulus $k$. Then they are solutions of:

$$ \frac{du}{dt} = vw, \quad \frac{dv}{dt} = -uw, \quad \frac{dw}{dt} = -k^2 uv $$

2. For every $0 < k < 1$, the JEEuler cases are equivalent and moreover equivalent to the NEuler case ($k = 1$).

*Proof* Let us prove the first assertion. From [25, p.36], one has:

$$ (\sin t)' = (\cos t)(\sin t), \quad (\cos t)' = -(\sin t)(\cos t), \quad (\cos t)' = -k^2(\sin t)(\cos t). $$

Hence, Assertion 1 follows.

Let us prove the second assertion. The product is defined by the relations:

$$ e_1 \cdot e_1 = e_2 \cdot e_2 = e_3 \cdot e_3 = 0 $$

$$ e_1 \cdot e_2 = -\frac{1}{2}k^2e_3, \quad e_1 \cdot e_3 = -\frac{1}{2}e_2, \quad e_2 \cdot e_3 = \frac{1}{2}e_1 $$

Let

$$ E_1 = \frac{1}{k}e_1, \quad E_2 = \frac{1}{k}e_2, \quad E_3 = e_3 $$

Then

$$ E_1 \cdot E_1 = E_2 \cdot E_2 = E_3 \cdot E_3 = 0 $$

and

$$ E_1 \cdot E_2 = -\frac{1}{2}E_3, \quad E_1 \cdot E_3 = -\frac{1}{2}E_2, \quad E_2 \cdot E_3 = \frac{1}{2}E_1 $$

□

## 4.2 Euler Dynamics

To avoid indices, we write:

$$ \frac{dx}{dt} = a_1yz, \quad \frac{dy}{dt} = a_2xz, \quad \frac{dz}{dt} = a_3xy $$

where we have:

$$ a_1 = \frac{I_2 - I_3}{I_1}, \quad a_2 = \frac{I_3 - I_1}{I_2}, \quad a_3 = \frac{I_1 - I_2}{I_3} $$

There are resonant relations between the coefficients since the system admits two first integrals (the kinetic energy and the norm of the momentum):

$$ E = I_1x^2 + I_2y^2 + I_3z^2 $$
$$M = I_1^2 x^2 + I_2^2 y^2 + I_3^2 z^2$$

Hence we have the relations:

$$I \cdot a = 0, \quad I^2 \cdot a = 0 \qquad (13)$$

One has $a_1, a_3 > 0$ and $a_2 < 0$. The system evolves on the momentum sphere which is invariant and the dynamics is integrable.

Recall that the system is written in the principal axes of coordinates denoted $E_1, E_2, E_3$. We denote by $b = (b_1, b_2, b_3)$ the vector representing the applied torque associated to the position of the gas jets in the satellite moving frame. It can be taken as a unit vector and restricted to the upper hemisphere.

First, the following holds from computations.

**Proposition 19 1.** The associated algebra $\mathcal{A}$ is simple since $a_1a_2a_3 \neq 0$ with the following nontrivial subalgebras: algebras of dimension one corresponding to the axis $E_i$ and two subalgebras of dimension 2 given by the planes defined by the solutions of $a_1z^2 - a_3x^2 = 0$.

2. The controlled Euler equation is controllable if and only if $L^0 = \mathbb{R}^3$, that is when two of the coefficients $b_i$ are nonzero and $a_3b_1^2 - a_1b_3^2 \neq 0$.

Recall that geometrically it corresponds to $b$ not belonging to a subalgebra of dimension 1 or 2.

*Sketch of the proof.* One computes the three vectors $v_1 = b, v_2 = [[Q, v_1], v_1], v_3 = [[Q, v_2], b]$ and forms the determinant $\det(v_1, v_2, v_3)$. It is of full rank if and only if $a_i b_j^2 - a_j b_i^2 \neq 0$ for $i \neq j$.

### 4.2.1 Classification

Using projective coordinates: $Y = y/x, Z = z/x$, the dynamics write:

$$\frac{dY}{dt} = x[a_2Z - a_1Y^2Z],$$

$$\frac{dZ}{dt} = x[a_3Y - a_1YZ^2],$$

where we can use the time reparameterization $d\sigma = x dt$.

Limit cycles can be evaluated using polar coordinates. The case of Euler dynamics with the resonant relation corresponds to nonlinear centers around $E_1$ and $E_3$. Hence all solutions are periodic, except the separatrices of the unstable equilibria. This completes the rough linear classification of the Euler algebra using the parameter signs $a_1, a_2, a_3$ only.
Denoting:

$$
Q = \begin{pmatrix} a_1 y z \\ a_2 x z \\ a_3 x y \end{pmatrix}, \quad b = \begin{pmatrix} b_1 \\ b_2 \\ b_3 \end{pmatrix}
$$
$$
[Q, b] = \begin{pmatrix} a_1 b_2 z + a_1 b_3 y \\ a_2 b_1 z + a_2 b_3 x \\ a_3 b_1 y + a_3 b_2 x \end{pmatrix} = \begin{pmatrix} 0 & a_1 b_3 & a_1 b_2 \\ a_2 b_3 & 0 & a_2 b_1 \\ a_3 b_2 & a_3 b_1 & 0 \end{pmatrix} \begin{pmatrix} x \\ y \\ z \end{pmatrix} = A \begin{pmatrix} x \\ y \\ z \end{pmatrix}
$$

We have det A = 2a₁a₂a₃b₁b₂b₃ and we denote by adj A = com Aᵀ = A⁻¹ det A.
Computing additional Lie brackets leads to:

$$
D'' = \det(b, [b, Q], Q) = a_2 a_3 b_1 x^2 (b_2 z - b_3 y) + a_1 a_3 b_2 y^2 (b_3 x - b_1 z) + a_1 a_2 b_3 z^2 (b_1 y - b_2 x)
$$

$$
D = 2a_1 b_2 b_3 (a_3 b_2^2 - a_2 b_3^2) x + 2a_2 b_1 b_3 (a_1 b_3^2 - a_3 b_1^2) y + 2a_3 b_1 b_2 (a_2 b_1^2 - a_1 b_2^2) z
$$

$$
D' = a_3^2 b_1 b_2 (b_1 y + b_2 x) (a_2 x^2 - a_1 y^2) + a_1^2 b_2 b_3 (b_2 z + b_3 y) (a_3 y^2 - a_2 z^2) \\
\qquad + a_2^2 b_1 b_3 (b_1 z + b_3 x) (a_1 z^2 - a_3 x^2)
$$

We can represent the classification with the bifurcation diagram. It depends upon
the parameters $a_1, a_3, a_2$. The vector $b$ is of unit length and by symmetry one can
assume that $b$ belongs to the equator or to the upper hemisphere. This gives a 3D
representation if we normalize: $a_1 = a_3 = 1, a_2 = -1$.

The cubic homogeneous geodesic equations are given by the vector field $H(x) = D(x)Q(x) - D'(x)b$ which can be projected on the 2-sphere. The set $D = 0$ is identified to the equator and we represent the dynamics in the Northern hemisphere. In the non-singular case, branches of the cubic $D'' = 0$ on this hemisphere will lead to periodic abnormal trajectories for the projected dynamics while the branch intersecting the equator $D = 0$ will form polycycles. Limit cases of cubic $D'' = 0$ correspond to the limit cases $b = (b_1, b_2 = 0, b_3)$ and $(b_1, b_2, b_3 = 0)$ where the cubic is reducible.

Next we use [12].

4.2.2 Bifurcation Analysis

Case (b₁, b₂, b₃ = 0)

From the general formulae we have:

$$
D = 2a_3 b_1 b_2 (a_2 b_1^2 - a_1 b_2^2) z
$$

$$
D'' = a_3 b_1 b_2 (a_2 x^2 - a_1 y^2) z
$$

$$
D' = a_3^2 b_1 b_2 (b_1 y + b_2 x) (a_2 x^2 - a_1 y^2)
$$

We must analyze the cubic differential equation:

$$
\frac{dx}{dt} = DQ_1 - D'b_1
$$
$$ \begin{aligned} \frac{dy}{dt} &= DQ_2 - D'b_2 \\ \frac{dz}{dt} &= DQ_3 \end{aligned} $$

**Lemma 3** The cubic $D''$ is reducible and $D$ is a factor of $D''$. The set $D'' = 0$ is the union of the plane $z = 0$ and the straight line $x = y = 0$.

Rays and lines of singular points of the dynamics are analyzed in the general framework of [12]. Since the plane $D = 0$ is given by $z = 0$, one must compute for the cubic systems the rays or lines of singular points solving:

$$ \begin{aligned} DQ_1 - D'b_1 &= \lambda x \\ DQ_2 - D'b_2 &= \lambda y \\ DQ_3 &= \lambda z \end{aligned} $$

They lead to the singularities of the dynamics outside $z = 0$ using the projective coordinates: $u = \frac{x}{z}$, $v = \frac{y}{z}$. From the last equation with $D \sim z$, one has $\lambda \sim Q_3$ where $\sim$ means up to a scalar, with $\lambda = DQ_3/z$, quadratic. The basic property is to use that $D'$ is the product of a linear form and a quadratic form. Computations lead to the following.

**Lemma 4** Outside the plane $z = 0$, the singular points of the dynamics are:

1. The point $0$ on the cubic $D'' = 0$ is a degenerated focus taking in normal coordinates the form: $\frac{d\rho}{dt} = c\rho^2 + o(\rho^2)$ with $c$ nonzero.

2. The remaining singular rays define two saddle points denoted $S_1, S_2$.

**Case $b = (b_1, b_2, b_3 = \epsilon)$, $\epsilon$ nonzero and small**

By perturbation, the cubic becomes regular with three symmetric ovals when projecting on the spheres: one forming an abnormal polycycle passing through $L_1, L_2$ and the two others being symmetric abnormal cycles contained respectively in the Northern and Southern hemisphere where $D = 0$ is identified to the equator.

The main property is the following:

**Theorem 20** For $\epsilon \neq 0$, the ovals in each hemisphere can be seen in appropriate coordinates as circles whose center corresponds to a ray solution, approximating the cubic $D'' = 0$ by a quadric. Hence we have a semi-Hopf bifurcation with the birth of a small abnormal cycle which can be stable or unstable in relation with instability properties of the small oval cycle.

**Case $b = (b_1, b_2 = 0, b_3)$**
In this case the algebraic forms are similar but the analysis differs since $a_1, a_3 > 0$ and $a_2 < 0$. The set $D = 0$ becomes the plane $y = 0$ while $D'' = 0$ is the union of $y = 0$ with the two planes defined by $a_1z^2 - a_3x^2 = 0$ which intersect at $O: x = z = 0$.

**Lemma 5** The dynamics admit in each hemisphere four polycycles decomposing in each hemisphere into four triangles: $L_3 \to 0 \to L_4$. Two of such triangles pass through the lines $L_1, L_2$. Two others contain in their interior singular rays solutions which are foci.

Case $b = (b_1, b_2 = \epsilon, b_3)$, $\epsilon$ nonzero and small

**Theorem 21** By perturbation, the two foci and their stability or instability remain and one of the abnormal triangles becomes a large cycle surrounding a focus.

The remaining phenomenon in the bifurcation diagram is to analyze the case when $L_3 = L_4$ which is described by the algebraic invariants and a saddle-connection which has no simple algebraic formulation in terms of invariants. See Appendix B for numerical simulations illustrating the bifurcation analysis.

## 4.3 Abnormal gates in Attitude Control

**Step 1.** We use the cycle associated to the oval $C_2$ in the previous construction. From our analysis it correspond to a real branch which can be parameterized thanks to section 2.4 by Jacobi elliptic functions in the t parameter but where the k modulus and the amplitude of the cycle are defined by the abnormal feedback. This corresponds to a reshaping of the Riemannian metric of the free motion. The t-parameterization is regular provided this branch contains no singular point. We denote by $T$ the period.

**Step 2.** After one period we observe a change of the attitude from $I \to R_0$ defined by integrating the dynamics. From Euler theorem this corresponds to a rotation of an angle $\alpha_0$ with respect an axis denoted $\Delta$. After n cycles by invariance this leads to $I \to R_0^n$. This corresponds to a periodic or quasi-periodic trajectory for the full dynamics.

**Conclusion.** Different abnormal gates can be constructed with an internal reshaping of the satellite. One invariant is sufficient to control the reshaping parameter which can be taken as the j-parameter associated to the bi-ratio of the lines $L_i, i = 1, 2, 3, 4$.

# 5 $\alpha$ and $\omega$ Limit Sets of Time Minimal Geodesics of Single-Input Homogeneous Systems in Dimension 3

We consider on the 3-dimensional space systems of the form:

$$ \frac{dv}{dt} = P(x) + ub, \quad u \in \mathbb{R}, $$
where $P$ is a homogeneous polynomial of degree $k \ge 1$. We distinguish between the cases $k = 2q$ (even) and $k = 2q + 1$ (odd).

Pontryagin's Maximum Principle for time minimal geodesics leads to introduce:

$$ \frac{dz}{dt} = \mathbf{H}_P(z) + u\mathbf{H}_b, $$

where for a given vector field $X$, $H_X = p \cdot X(v)$ is the Hamiltonian lift and $\mathbf{H}_X$ is the associated Hamiltonian vector field. The Pontryagin's maximum principle tells us that candidates for time minimizers satisfy the Hamiltonian dynamics with the constraint:

$$ H_b = p \cdot b = 0. $$

Differentiating twice with respect to time yields:

$$ \begin{gathered} p \cdot [b, P] = 0, \\ p \cdot [[b, P], P] + u_s p \cdot [[b, P], b] = 0. \end{gathered} $$

Then, introducing the determinants

$$ \begin{align*} D &= \det[b, [b, P], [[b, P], b]], \\ D' &= \det[b, [b, P], [[b, P], P]], \\ D'' &= \det[b, [b, P], P], \end{align*} $$

and, by homogeneity, we have the degrees:

$$ (i) D: 0 + k - 1 + k - 2 = 2k - 3 $$

$$ (ii) D': 0 + k - 1 + 2k - 2 = 3k - 3 $$

$$ (iii) D'': 0 + k - 1 + k = 2k - 1 $$

$$ (iv) D'/D = 3k - 3 - 2k + 3 = k $$

**Lemma 6** The vector field

$$ X_s := P - \frac{D'}{D} b $$

is a homogeneous meromorphic vector field of degree $k$.

Using the reparameterization $dt = d\sigma \cdot D(v)$, we obtain the analytic homogeneous vector field:

$$ X'_s = D \cdot P - D' \cdot b \quad (14) $$

of degree $3k - 3$.

## 5.1 General Framework for Time Minimal Geodesics of Single-Input Control Systems

Consider the system:

$$ \frac{dx}{dt} = X(x) + uY(x), \quad u \in \mathbb{R}. \quad (15) $$
Remark 3 In general, the minimal time to steer $x_0$ to $x_1$ may fail to exist since $u$ is not uniformly bounded (e.g., linear case). In the homogeneous case, one must distinguish between odd and even degrees. See [21] for the case $k = 2m + 1$, and [11] for the quadratic case.

The geodesics dynamics is given by:

$$ \frac{dz}{dt} = \mathbf{H}_X + u\mathbf{H}_Y \qquad (16) $$

with the two constraints:

$$ H_Y = \{H_Y, H_X\} = 0. \qquad (17) $$

Then, the singular control is given as a dynamic feedback by differentiating the previous relation

$$ \{\{H_Y, H_X\}, H_X\} + u_s\{\{H_Y, H_X\}, H_Y\} = 0. \qquad (18) $$

The cotangent bundle is endowed with the standard symplectic form $\omega = dp \wedge dx$, and let $\Omega$ be the associated volume form. We denote by $\varphi$ the set of symplectomorphisms preserving the fiber associated to state diffeomorphisms $\varphi$ on $T^*M$, called Mathieu symplectomorphisms defined in canonical coordinates by $x = \varphi(y)$, $p^\top = \left(\frac{\partial \varphi^{-1}}{\partial y}(y)\right)^\top q$.

The extended time minimal Hamiltonian is given by:

$$ H' = H_X + uH_Y + p^0 \cdot t, \qquad (19) $$

where $p^0$ is constant. We consider both geodesics associated with the time minimal or maximal problem.

Lemma 7 Assuming the Legendre-Clebsch condition: $\{\{H_Y, H_X\}, H_Y\} \neq 0$ along the geodesics, then the geodesics dynamics is a Hamiltonian vector field on the reduced symplectic space using the two constraints. However, in such reduction the state space is not preserved, hence we shall always consider the geodesic systems in the original space with the two constraints.

## 5.2 Back to the Quadratic 3D Case

Using time reparameterization, one must analyze the cubic dynamics:

$$ \frac{dv}{dt} = DQ - D'b = H(v) \qquad (20) $$

which carries all the dynamical information, provided we take into account the change of time parameterization. In particular, in the domain $D < 0$, we must reverse the time orientation of the geodesics.
Lemma 8 ([9]) 1. The plane $D = 0$ is generated by the two lines $L_1, L_2$ (provided they are not distinct) and is an invariant plane.

2. The dynamics can be projected onto the 2-sphere using again an explicit time reparameterization:

$$ \frac{dr}{dt} = \langle x', H(x') \rangle r, \qquad (21) $$

$$ \frac{dx'}{d\sigma} = H(x') - \langle x', H(x') \rangle x', \qquad (22) $$

where the formulas come from Section 2.1.

3. The line $L_1 = [b]$ is a ray corresponding to a node, and the line $L_2$ consists of non-isolated equilibria; the associated singularity on the sphere is a stellar node.

4. The cubic $D'' = 0$ intersects the equator $D = 0$ at four points on the sphere, corresponding to two projective directions $L_1, L_2$.

Consider equation (20), one can apply the following result.

**Theorem 22** ([26]) Consider a smooth vector field on the 2-sphere with isolated singularities. Then limit sets of trajectories are either:

1. Isolated equilibria,

2. Limit cycles,

3. Monodromic polycycles.

## 5.3 Applications

This result can be applied to the regularized geodesics dynamics. To return to the control-theoretic interpretation, one must perform a double transfer:

(i) To the time minimal stratification, since the geodesics dynamics are only candidates as small time minimizers,

(ii) To the control system in $\mathbb{R}^3$ using the true time parameterization.

The first point is solved by the following high-order maximum principle [24].

**Proposition 23** ([24]) Outside the plane $D = 0$, the geodesics are either:

1. Small time hyperbolic minimizers if contained in the domain: $DD'' > 0$,

2. Small time elliptic maximizers contained in the domain: $DD'' < 0$,

3. Abnormal minimizers which are both small time minimizers and maximizers, contained in the set $D'' = 0$ stratified in the regular case by one or three branches.

The second point comes from our general analysis, proving the following.
**Proposition 24** ([9]) 1. At a point of $L_2$, the geodesics can cross smoothly the plane $D = 0$ (at finite distance).

2. At a point of $L_3, L_4$, there exists a single trajectory crossing smoothly the plane $D = 0$.

3. On the contrary, the line $L_1$ forms a set of asymptotic directions in the plane $D = 0$ where geodesics are attracted or repulsed, and the geodesic flow is not complete in those directions.

# 6 Conclusion

In this article we have completed and gathered in a single article a final analysis about the Euler-Poisot with a single torque, which were disseminated in different sources. Finally, it can be applied to the attitude control with a single reaction jet motivated by recent CubeSat projects. The 3-dimensional model can be generalized to the $n$-dimensional case but thanks to the dimension we were able to present a complete analysis in the frame of feedback classification of single input control systems where the invariants are described in relations with the classification of the limit uncontrolled case. The analysis leads to the introduction of commutative but non-associative algebra initiated in the seminal reference [27] and completed later by numerous articles see for instance [22]. Even in dimension 3, the problem is complex due the existences of cycles for such dynamics in relation with Hilbert's sixteen problem. Thanks to the case of geodesics dynamics algebraic abnormal cycles will occur in our situation due to the existence of a cubic abnormal foliation, which is related to the stationary rotations of the rigid body. Furthermore, we propose a bifurcation diagram for the geodesic dynamics which can be used to design new control techniques for modern under-actuated low cost satellites. It is an ideal setting to experiment with such complex feedback control laws on a new generation of satellites compared to previous attitude control projects, where only simple feedback control laws were implemented.

# References

[1] V. I. Arnold, *Sur la géométrie différentielle des groupes de Lie de dimension infinie et ses applications à l'hydrodynamique des fluides parfaits*, Ann. Inst. Fourier 16(1) (1966), pp. 319-361.

[2] V. I. Arnold, *Méthodes Mathématiques de la Mécanique Classique*, Editions Mir, Chap. 6 and Appendix 2, 1974.

[3] V. I. Arnold, *Singularity theory*, Birkhäuser, Basel-Boston-New York, 2000.

[4] J. Baillieul, *Controllability and observability of polynomial dynamical systems*. Nonlinear Anal., Theory Methods Appl. 5 (1981), pp. 543-552.

[5] M. Banchi, *Typical Ranks of Ternary Cubic Forms over $\mathbb{R}$*, PhD Thesis, Università degli Studi di Firenze, 2016.
[6] M. Bhargava and A. Shankar, *Binary quartic forms having bounded invariants, and the boundedness of the average rank of elliptic curves*, Ann. of Math. (2) **181** (2015), no. 1, 191–242.

[7] B. Bonnard, *Contrôlabilité des systèmes non linéaires*, C. R. Acad. Sci. Paris Sér. I Math. 292 (1981), pp. 535–537.

[8] B. Bonnard, *Contrôlabilité des systèmes mécaniques sur les groupes de Lie*, SIAM J. Control Optim. 22(5) (1984), pp. 711–722.

[9] B. Bonnard, *On Singular Extremals in the Time Minimal Control Problem in $\mathbb{R}^3$*, SIAM J. Control Optim. 23(5) (1985), pp. 794–802.

[10] B. Bonnard, *On Singular Extremals in the Time Minimal Control Problem*, In theory and applications of nonlinear control systems Ed. C.I. Byrnes, A. Lindquist, Northb-Holland (1986), pp.439-452.

[11] B. Bonnard, *Quadratic Control Systems*, in: Math. Control Signals Syst. 4(2) (1991), pp. 139–160, Springer-Verlag.

[12] B. Bonnard, *Feedback Equivalence for Nonlinear Systems and the Time Minimal Optimal Control Problem*, SIAM J. Control Optim. 29(6) (1991), pp. 1300–1321.

[13] B. Bonnard, O. Cots, N. Shcherbakova, *The Serret-Andoyer Riemannian Metric and Euler-Poinsot Rigid Body Motion*, Math. Control Relat. Fields 3(3) (2013), pp. 287–302.

[14] B. Bonnard, J. Rouot, *Feedback Classification and Optimal Control with Applications to the Controlled Lotka-Volterra Model*, Optimization, 74(5) (2024), 1173–1196.

[15] N. A. Chaturvedi, A. K. Sanyal, N. H. McClamroch, *Rigid Body Attitude Control Using Rotation Matrices for Continuous, Singularity-Free Control Laws*, IEEE Control Syst. 31(3) (2011), pp. 30–51.

[16] P. E. Crouch, *Spacecraft Attitude Control and Stabilization: Applications of Geometric Control Theory to Rigid Body Models*, IEEE Trans. Autom. Control 29(4) (1984), pp. 321-331.

[17] P.E. Crouch and B. Bonnard, *An Appraisal of Linear Analytic Systems Theory with Applications to Attitude Control*, ESA ESTEC Contract report, 1980.

[18] N. Duif, *Transforming a General Cubic Elliptic Curve Equation to Weierstrass Form: A Sage Implementation*, Technische Universiteit Eindhoven, 2011.

[19] T. Fisher, *Testing Equivalence of Ternary Cubics*, Algorithmic number theory, Lecture Notes in Comput. Sci. 4076, Springer, Berlin, 2006, pp. 333–345.
[20] V. Jurdjevic, *Geometric Control Theory*, Cambridge University Press, 1997.

[21] V. Jurdjevic and I. Kupka, *Polynomial control systems*, Math. Ann. **272** (1985), 361–368.

[22] M. K. Kinyon, A. A. Sagle, *Quadratic Dynamical Systems and Algebras*, J. Differ. Equations 117 (1995), pp. 67–126.

[23] I.A. Kogan and M.M. Maza, *Computations of Canonical Forms of Ternary Cubics*. Symposium Symbolic Algebr. Comput. (ISAAC), ACM Press, New York, NY, 2002, pp. 151–160.

[24] A. J. Krener, *The High Order Maximal Principle and Its Application to Singular Extremals*, SIAM J. Control Optim. 15(2) (1977), pp. 256–293.

[25] D. F. Lawden, *Elliptic Functions and Applications*, Applied Mathematical Sciences, vol. 80, Springer-Verlag, 1989.

[26] S. Lefschetz, *Differential Equations: Geometric Theory*, John Wiley and Sons, 1957.

[27] L. Markus, *Quadratic Differential Equations and Non-Associative Algebras: Contributions to the Theory of Nonlinear Oscillations*, Ann. Math. Stud., vol. 45, Princeton Univ. Press, 1960, pp. 185–213.

[28] J. Martinet, * Singularities of Smooth Functions and Maps*, Cambridge Univ. Press (1982).

[29] J. J. Moreau, *Mécanique Classique, Tome 1 et 2*, Masson, 1971.

[30] A. Mourtada, *Cyclicité Finie des Polycycles Hyperboliques de Champs de Vecteurs du Plan. Algorithme de Finitude*, Ann. Inst. Fourier 41(3) (1991), pp. 719–753.

[31] A. Mourtada, *Bifurcation de Cycles Limites au Voisinage de Polycycles Hyperboliques et Génériques à Trois Sommets*, Ann. Fac. Sci. Toulouse 6e série, 3(2) (1994), pp. 259–292.

[32] C. Murakami, *Analytical Solution of the Euler-Poisson Problem*, J. Geom. Symmetry Phys. 60 (2021), pp. 25-46.

[33] V. V. Nemitskii, V. V. Stepanov, *Qualitative Theory of Differential Equations*, Princeton, NJ, USA: Princeton Univ. Press, 1960.

[34] L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze and E. F. Mishchenko, *The Mathematical Theory of Optimal Processes*, Pergamon Press, Oxford, 1964, 362 pages.
[35] C. L. Siegel, J.K. Moser, *Lectures on Celestial Mechanics*, Springer-Verlag, Berlin, 1971, 302 pages.

[36] J. H. Silverman, *The Arithmetic of Elliptic Curves*, Graduate Texts in Mathematics, vol. 106, Springer-Verlag, 1986.

[37] R. J. Walker, *Algebraic Curves*, Springer-Verlag, New York, 1978.

[38] P. M. H. Wilson, *Some Geometry and Combinatory for the S-Invariant of Ternary Cubics*, Experiment. Math. 15, 479–490 (2006).

# A The 3D Case: Feedback Classification and Optimality

We restrict to the 3D-case; see [12] for the general case. Note that small time optimality can be discussed using [24]. In particular, normal-exceptional trajectories are candidates as time minimal or time maximal trajectories.

We consider a single-input affine control system of the form:

$$ \frac{dx}{dt} = X(x) + uY(x), $$

where $X$ and $Y$ are smooth vector fields.

**Definition 6** Define by $G_f$ the feedback group formed by the set of triples $g = \{(\varphi, \alpha, \beta)\}$ where $\varphi$ is a diffeomorphism and $u = \alpha(x) + \beta(x)u'$, which acts on the system $(X, Y)$ by changes of coordinates and feedback. The product rule is induced by the action.

Introduce the following determinantal mappings:

$$ D = \det\{Y, [Y, X], [[Y, X], Y]\}, $$

$$ D' = \det(Y, [Y, X], [[Y, X], X]), $$

$$ D'' = \det(Y, [Y, X], X). $$

**Definition 7** Let $V, W$ be two $\mathbb{R}$-vector spaces and let $G$ be a group acting linearly on $V$ and $W$.

* A *character* is a homomorphism $\chi$ from $G$ into $\mathbb{R} \setminus \{0\}$.

* A *semi-invariant* of weight $\chi$ of $G$ acting on $V$ is a map $\lambda: V \to \mathbb{R}$ such that for all $g \in G$ and for all $x \in V$, $\lambda(g \cdot x) = \chi(g)\lambda(x)$.

* If $\chi = 1$, a semi-invariant is called an *invariant*.

* A *semi-covariant* of weight $\chi$ for the actions of $G$ on $V, W$ is a map $\lambda: V \to W$ such that for all $x \in V$, $\lambda(g \cdot x) = \chi(g)g \cdot \lambda(x)$.

* It is a *covariant* if $\chi = 1$.
Lemma 9 If Z is a vector field, we denote $\varphi^*Z$ the image of Z by the diffeomorphism $\varphi$. If $(\varphi, \alpha, \beta) \in G_f$ then the following hold:

1. $D^{\varphi^*X, \varphi^*Y}(x) = (\det(\partial\varphi^{-1}/\partial x))D^{X,Y}(\varphi(x))$,

2. $D^{\prime\varphi^*X, \varphi^*Y}(x) = (\det(\partial\varphi^{-1}/\partial x))D^{\prime X,Y}(\varphi(x))$,

3. $D^{\prime\prime\varphi^*X, \varphi^*Y}(x) = (\det(\partial\varphi^{-1}/\partial x))D^{\prime\prime X,Y}(\varphi(x))$,

4. $D^{X+Y,\alpha,\beta} = \beta^4 D^{X,Y}$,

5. $D^{\prime X+Y,\alpha,\beta Y} = \beta^3(D^{\prime X,Y} + \alpha D^{X,Y})$,

6. $D^{\prime\prime X+Y,\alpha,\beta Y} = \beta^2 D^{\prime\prime X,Y}$.

**Proposition 25** Singular trajectories of minimal order are solutions of the dynamics:

$$ \frac{dx}{dt} = X_s(x) := X(x) + u_s(x)Y(x) \qquad (23) $$

where the singular control is given by the feedback:

$$ u_s(x) = -\frac{D'(x)}{D(x)}. \qquad (24) $$

The sets $D'' = 0$, $DD'' > 0$ and $DD'' < 0$ are invariant for the singular dynamics and correspond respectively to abnormal, hyperbolic and elliptic geodesics.

According to the high order maximum principle [24], we have:

**Proposition 26** Hyperbolic geodesics are small time minimizing and elliptic geodesics are small time maximizing trajectories, while abnormal geodesics are both.

**Theorem 27** For the action of the feedback group in the 3D-case where the action of mappings is reduced to the $\varphi$-action, we have the following:

1. $\lambda_s : (X, Y) \to X_s$ is a covariant,

2. $\lambda_s : (X, Y) \to D(x)X(x) - D'(x)Y(x)$ is a semi-covariant (coding the one-dimensional foliation),

3. $\lambda_2 : (X, Y) \to D^{X,Y}$ is a semi-covariant,

4. $\lambda_3 : (X, Y) \to D^{\prime\prime X,Y}$ is a semi-covariant,

5. $\lambda_5 : (X, Y) \to D^{X,Y} \cap D^{\prime X,Y}$ is a semi-covariant.

In particular, the sets $D = 0$, $D'' = 0$ and $D = D' = 0$ have an invariant meaning.

# B Symbolic and numerical computations

## B.1 Classification of the real ternary cubic $D''(x) = 0$

**Algorithm for computing the Weierstrass form**

Introduce $\delta_1 = a_1 b_2^2 - a_2 b_1^2$, $\delta_2 = a_2 b_3^2 - a_3 b_2^2$, $\delta_3 = a_3 b_1^2 - a_1 b_3^2$. The following lemmas specify the inputs of the algorithm.

**Lemma 10** The point $L_1$ is singular if and only if
* two elements among $a_1, a_2, a_3$ or among $b_1, b_2, b_3$ are zero,

* or $(a_i, b_i) = (0, 0)$ for some $i \in \{1, 2, 3\}$,

* or $b_2 = \delta_3 = 0$, or $b_1 = \delta_2 = 0$, or $b_3 = \delta_1 = 0$, or $\delta_1 = \delta_3 = 0$.

**Lemma 11** If $L_1$ is not a singular point of $C$, then the tangent to $C$ at $L_1$ contains $L_2$, and $L_2$ is not a singular point of $C$.

The following proposition is based on Appendix B.1 and [36].

**Proposition 28** (Weierstrass canonical form) Assume that for $i = 1, 2, 3$ we have $a_i, b_i \in \mathbb{Q}^*$ and that $\delta_1\delta_2\delta_3 \neq 0$. There exist projective coordinates in which the cubic $C$ has an equation of the form:

$$x_2^2 x_3 = x_1^3 + A x_1 x_3^2 + B x_3^3$$

where

$$A = -\frac{\sqrt{2}(a_1 a_2)^{2/3} a_3^2}{3(b_3^2 \delta_2^2 \delta_3^2)^{2/3}} (b_1^4 \delta_2^2 + b_2^4 \delta_3^2 + b_3^4 \delta_1^2) < 0, \qquad (25)$$

$$B = \frac{4a_1 a_2 a_3 (b_2^2 \delta_3 - b_3^2 \delta_1)(b_3^2 \delta_1 - b_1^2 \delta_2)(b_2^2 \delta_3 - b_1^2 \delta_2)}{27 b_3^2 \delta_2^2 \delta_3^2}.$$

The discriminant $\Delta := -(4A^3 + 27B^2)$ is

$$\Delta = \frac{16a_1^2 a_2^2 a_3^6 b_1^4 b_2^4 \delta_1^2}{\delta_2^2 \delta_3^2} > 0 \qquad (26)$$

and the invariant $j = -1728 \frac{4A^3}{\Delta}$ of the cubic $C$ is

$$j = \frac{256 (b_1^4 \delta_2^2 + b_2^4 \delta_3^2 + b_3^4 \delta_1^2)^3}{b_1^4 b_2^4 b_3^4 \delta_1^2 \delta_2^2 \delta_3^2}. \qquad (27)$$

*Proof* We cannot apply an existing algorithm such as those implemented in the library *Elliptic curves over the Rational Field* of SageMath<fnref n="1" /> because we want closed-form expressions for the invariants $A$ and $B$, whereas this library requires fixed parameters $a_i, b_i$, $i = 1, 2, 3$.

Therefore we implemented a Mathematica script<fnref n="2" /> to execute the following steps:

1. Under the assumption $a_1 a_2 a_3 b_1 b_2 b_3 \delta_1 \delta_2 \delta_3 \neq 0$, the cubic $C$ has a tangent at $L_1$ and $L_2$. Choose axes $(U, V, W)$ in the projective plane such that $W = 0$ is the tangent line to $C$ at $L_1$. This tangent line intersects $C$ at $L_2$, and we take the $U = 0$ axis to be tangent to $C$ at $L_2$. Finally, we choose $V = 0$ to be any line (other than $W = 0$) passing through $L_1$. Computing with $V = b_2 x_1 - b_1 x_2$, the determinant of the linear transformation $(x_1 : x_2 : x_3) \mapsto (U : V : W)$ is $-8 a_1 a_2 a_3^2 b_1^2 b_2^2 b_3 \delta_1^2 \delta_2 \delta_3$.

In this new frame the point $L_1$ has coordinates $(U : V : W) = (1 : 0 : 0)$, the point $L_2$ is $(0 : 1 : 0)$ and the list of monomials of $C$ (after division by $UV^2$) is $(U^2W, UV^2, UVW, UW^2, VW^2, W^3)$.

<footnote n="1">¹https://doc.sagemath.org/html/en/reference/arithmetric_curves/sage/schemes/ell_rational_field.html</footnote>

<footnote n="2">²The Mathematica script is available at https://github.com/jeremyrouot/cubsat/blob/main/compute-weierstrass_nform.m</footnote>
2. Apply the Cremona transformation $(U : V : W) \mapsto (K^2 : LN : KN)$ followed by the substitution $M \leftarrow 1$; the list of monomials of $C$ becomes $(K^3, K^2, KL, K, L^2, L)$.

3. Complete the square in the $L$ variable by substituting $L \leftarrow 1/2(L - (a'K + b'))$, where $a'$ and $b'$ are respectively the coefficients of the monomials $KL$ and $L$.

4. Normalize the coefficient of the $K^3$ monomial to $-1$ and the $L^2$ monomial to $1$ by scaling the variables $K$ and $L$.

5. Remove the $K^2$ monomial with a substitution of the form $K \leftarrow K - c'$ and a suitable choice of $c'$.

The algorithm returns the values of the invariants $A$ and $B$ such that the equation of $C$ is $y^2 = x^3 + Ax + B$, and yields the expressions of $A, B, \Delta$ and $j$ given in the statement of this proposition. The fact that $A < 0$ comes from the condition $\Delta > 0$, which implies $4A^3 < -27B^2$ and therefore $A < 0$. $\square$

**Remark 4** • When $a_1a_2a_3b_1b_2b_3\delta_1\delta_2\delta_3 = 0$, the cubic $D'' = 0$ is reducible and can be analyzed directly (see Section B.2). In this case, $C$ is singular and is either a parabola nodata or a parabola punctata in the Newton classification. We do not obtain the parabola pura case – later refined by Möbius – nor the parabola cuspidata case.

• We can initialize the algorithm with the base point $P = e_i$, $i = 1, 2, 3$, where $(e_i)_{i=1,2,3}$ is the canonical basis of $\mathbb{R}^3$, and the tangent of $C$ at $e_i$, $i = 1, 2, 3$ intersects $C$ at $L_1$. The assumption of Proposition B.1 can then be relaxed to $a_1a_2a_3b_1b_2b_3\delta_2 \neq 0$ with the choice $P = e_1$ and $Y = x_3$. The expressions of $A, B, \Delta$ and $j$ become:

$$ 
\begin{align} 
A &= -\frac{a_1^{2/3} (b_1^4 \delta_2^2 + b_2^4 \delta_3^2 + b_3^4 \delta_1^2)}{6a_2^2 a_3^2 b_1^4 (b_2 \delta_2)^{4/3}} < 0, \nonumber \\ 
B &= \frac{a_1 (b_2^2 \delta_3 - b_3^2 \delta_1) (b_3^2 \delta_1 - b_1^2 \delta_2) (b_2^2 \delta_3 - b_1^2 \delta_2)}{27 a_2^3 a_3^3 b_2^2 b_1^6 \delta_2^2}, \tag{28} \\ 
\Delta &= \frac{a_1^2 b_3^4 \delta_1^2 \delta_3^2}{a_2^6 a_3^6 b_1^8 \delta_2^2}, \quad j = \frac{256 (b_1^4 \delta_2^2 + b_2^4 \delta_3^2 + b_3^4 \delta_1^2)^3}{b_1^4 b_2^4 b_3^4 \delta_1^2 \delta_2^2 \delta_3^2}. \nonumber 
\end{align} 
$$

**Proposition 29** *Under the notations and assumptions of Proposition B.1, the cubic curve $C$ is a nonsingular parabola campaniformis cum ovali, and its three flex points all lie on the unbounded component of $C$.*

*Proof* From Proposition B.1 we know that $\Delta > 0$, so the depressed cubic $x^3 + Ax + B$ has three distinct real roots $r_1 < r_2 < r_3$, with $r_1 < 0$, $\text{sign}(r_2) = \text{sign}(B)$, and $r_3 > 0$. The curve $C$ then consists of an unbounded component for $x \ge r_3$ and a bounded oval for $r_1 \le x \le r_2$. The flex points are given by the roots of the 3-division polynomial

$$ \psi_3(x) = 3x^4 + 6Ax^2 + 12Bx - A^2. $$

Since its discriminant is $\text{Disc}_x(\psi_3) = -6912\Delta^2 < 0$, $\psi_3$ has exactly two distinct real roots $\gamma_- < \gamma_+$. For $i = 1, 2, 3$, we compute

$$ \psi_3(r_i) = 3r_i(-Ar_i - B) + 6Ar_i^2 + 12Br_i - A^2 = 3Ar_i^2 + 9Br_i - A^2, $$
which is negative when $\Delta > 0$. Moreover, since $y^2 = \psi'_3(x)$ is the (affine) equation of $C$, we deduce that $\gamma_- < r_1$ and $\gamma_+ > r_3$.

Therefore $\gamma_+$ yields two affine flex points on the unbounded branches

$$y = \pm\sqrt{(x - r_1)(x - r_2)(x - r_3)}, \quad x \ge r_3,$$

while the third flex is the point at infinity $(0 : 1 : 0)$, which also lies on the unbounded component. Hence $C$ has exactly three flex points, all belonging to its unbounded component. $\square$

<figure><img src="image_2.png" /></figure>

## B.2 Numerical results

We complement the previous bifurcation analysis with numerical computations. In each case, we display the orthogonal projection onto the plane $D = 0$. Specifically, each figure includes the projection of:

* the lines $L_i$ (for $i = 1, \dots, 4$) (if both of the latter exist),

* the images of the determinantal varieties $D, D',$ and $D''$.

* the flex points of the cubic curve $C: D''(x) = 0$,

* the tangential singular dynamics $X_s^T: \mathbb{R}^3 \ni v \mapsto X_s(v) - v \cdot X_s(v)v$, and

* the equilibria of $T_{X_s}$.

### Reducible cases

* $b = (2, 0, -1)$

<figure><img src="image_3.png" /><figcaption>Fig. 1: Projection on the plane $D = 0$ of the tangential singular flow. The computations were performed for $(a_1, a_2, a_3) = (1, -1, 1)$ and $(b_1, b_2, b_3) = (2, 0, -1)$. We are in the case where $D$ divides $D''$ and any point of $D = 0$ is a flex point of $D'' = 0$. The equilibria $\pm(0.3, -0.38)$ are focus, $(0, 0)$, $L_3$ and $L_4$ are saddle points, while the other equilibria $L_1, L_2$ are nodes.</figcaption></figure>

* $b_3 = 0$
The eigenvalues of the Jacobian of tangential singular dynamics at 0 are $\pm 2ib_1b_2(b_1^2+b_2^2)$ and direct computations of the first Lyapunov coefficient $\ell_1$ yield $\ell_1 = b_1^2b_2^2/2$. Therefore, $(0,0)$ is an unstable degenerate focus.

<figure><img src="image_4.png" /><figcaption>Fig. 2: Projection on the plane $D = 0$ of the tangential singular flow. The computations were performed for $(a_1, a_2, a_3) = (1, -1, 1)$ and $(b_1, b_2, b_3) = (2, -1, 0)$. We are in the reducible case where $D$ divides $D''$. In this case, any point of $D = 0$ is a flex point of $D'' = 0$. The equilibrium $(0,0)$ is an unstable degenerate focus, $\pm(-0.18, 0.52)$ are saddle points while the other equilibria $L_1$ and $L_2$ are nodes.</figcaption></figure>

• Uncontrollable case, $b = (2, 1, 2)$.

<figure><img src="image_5.png" /><figcaption>Fig. 3: Projection on the plane $D = 0$ of the tangential singular flow. We are in the case where $D$ divides $D'$ so any point of $D = 0$ is an equilibrium of $T_{X_s}$ and $D$ divides $D''$ so any point of $D = 0$ is a flex point of the curve $D'' = 0$. The computations were performed for $(a_1, a_2, a_3) = (1, -1, 1)$ and $(b_1, b_2, b_3) = (2, 1, 2)$. This is the uncontrollable case.</figcaption></figure>
Irreducible cases

$$
\bullet b = (2, 1/5, -1)
$$

<figure><img src="image_6.png" /><figcaption>Fig. 4: Projection on the plane $D = 0$ of the tangential singular flow. The computations were performed for $(a_1, a_2, a_3) = (1, -1, 1)$ and $(b_1, b_2, b_3) = (2, 1/5, -1)$. We are in the case where $D''$ is irreducible. The equilibria $(0.33, -0.46)$ and $(-0.28, 0.29)$ are focus, the point $(0.0002, -0.07)$ is a saddle point, while the other equilibria $L_1, L_2$ are nodes.</figcaption></figure>

$$
\bullet b = (2, -1, 0.15)
$$

<figure><img src="image_7.png" /><figcaption>Fig. 5: Projection on the plane $D = 0$ of the tangential singular flow. The computations were performed for $(a_1, a_2, a_3) = (1, -1, 1)$ and $(b_1, b_2, b_3) = (2, -1, 0.15)$. We are in the case where $D''$ is irreducible. The equilibrium $(0, 0)$ is a center, the points $(-0.22, 0.57)$ and $(0.14, -0.48)$ are saddle points, $(-0.03, 0.01)$ is a center while the other equilibria on $D = 0$ are nodes.</figcaption></figure>
C Further numerical simulations

C.1 Variation of b<sub>3</sub> for L<sub>2</sub> = L<sub>3</sub> transition.

Projection on the plane $D = 0$ of the tangential singular flow, computed for $(a_1, a_2, a_3) = (-2, -1, -0.207107)$ with $(b_1, b_2)$ fixed at $(-1, -1)$ and $b_3$ varying.

<figure><img src="image_8.png" /><figcaption>Fig. 6: Case $b_3 = -0.68$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7: Case $b_3 = -0.78$ (transition $L_2 = L_3$).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 8: Case $b_3 = -0.88$.</figcaption></figure>

## C.2 Variation of $a_3$ for $L_2 = L_3$ transition.

<figure><img src="image_11.png" /><figcaption>Fig. 9: Case $a_3 = -0.107107$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 10: Case $a_3 = -0.207107$ (transition $L_2 = L_3$).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11: Case $a_3 = -0.307107$.</figcaption></figure>
### C.3 Variation of $b_3$ for $L_3 = L_4$ transition.

<figure><img src="image_14.png" /><figcaption>Fig. 12: Case $b_3 = 197/30$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 13: Case $b_3 = 20/3$ (transition $L_3 = L_4$).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 14: Case $b_3 = 203/30$.</figcaption></figure>