# Optimal control of a solar sail

Alesia Herasimenka, Lamberto Dell'Elce, Jean-Baptiste Caillau, Jean-Baptiste Pomet

► To cite this version:

Alesia Herasimenka, Lamberto Dell'Elce, Jean-Baptiste Caillau, Jean-Baptiste Pomet. Optimal control of a solar sail. 2023. hal-04158202v3

HAL Id: hal-04158202

https://inria.hal.science/hal-04158202v3

Preprint submitted on 3 Apr 2024

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
# Optimal control of a solar sail

Alesia Herasimenka, Lamberto Dell'Elce, Jean-Baptiste Caillau, Jean-Baptiste Pomet

► To cite this version:

Alesia Herasimenka, Lamberto Dell'Elce, Jean-Baptiste Caillau, Jean-Baptiste Pomet. Optimal control of a solar sail. Optimal Control Applications and Methods, 2024, 45, pp.2837-2855. 10.1002/oca.3189 . hal-04158202v4

HAL Id: hal-04158202

https://inria.hal.science/hal-04158202v4

Submitted on 3 Sep 2024

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
# ARTICLE TYPE

# Optimal control of a solar sail

Alesia Herasimenka*<fnref n="1" /> | Lamberto Dell'Elce<fnref n="2" /> | Jean-Baptiste Caillau<fnref n="1" /> | Jean-Baptiste Pomet<fnref n="2" />

## Correspondence

*Email: alesia.herasimenka@univ-cotedazur.fr

## Present Address

This is sample for present address text this is sample for present address text

## Summary

Optimal control of a solar sail orbiting around a fixed center of mass is considered. The sail is modelled as a plane surface with two sides having similar optical properties. The control is assumed to be the attitude of the sail and is represented as an element of the projective plane. Mapping this plane into the three dimensional ambient space to compute the actual force generated by a given attitude results in a generally non-convex set of admissible control values. A suitable convex relaxation is introduced to study the optimality system associated with maximising the change of the sail orbital parameters in a given direction along one orbit. Existence for both the relaxed and the original problems are deduced. Building on previous works<fnref n="1" />,<fnref n="2" />,³, a refined analysis of the control structure is given, proving rigorously that it contains only switchings, and a global bound of the number of switchings is obtained. In order to compute effective solutions, a multiple shooting approach is retained that is well suited for systems with a sequence of switchings between various control subarcs. Three issues are addressed: initialisation of shooting by means a tailored semi-definite relaxation that takes advantage of good convergence properties of modern convex optimisation algorithms; changes in the control structure that are accommodated by coupling shooting with differential continuation; implicit character of the Hamiltonian maximisation when using Pontrjagin maximum principle that is taken care of by incorporating the associated equation for the dynamical feedback into the shooting procedure. The method is illustrated on an example from NASA for which the sail inclination is optimally changed over one orbital period.

## KEYWORDS:

solar sailing, non-convex constraints, SDP relaxation, multiple shooting, differential continuation

# 1 | INTRODUCTION

## 1.1 | Solar radiation pressure

Solar sails are satellites that use solar radiation pressure (SRP) as propulsion for orbital manœuvres. SRP is caused by the interaction between photons and the surface of the sail, and its magnitude depends on the distance between the Sun and the sail, $r$. Specifically, denote by $\Phi_{SR} = 1367 \text{ W m}^{-2}$ the solar flux at $r_{\oplus} = 1 \text{ AU}$ distance (which is the average Sun-Earth distance), and by $c$ the speed of light; then, a simple model for the SRP is given by:⁴, Chap. 3

$$P_{SR} = \frac{\Phi_{SR}}{c} \left( \frac{r_{\oplus}}{r} \right)^2. \qquad (1)$$

<footnote n="1">¹Université Côte d'Azur, CNRS, Inria, LJAD, Nice, France</footnote>

<footnote n="2">²Inria, Sophia Antipolis, France</footnote>
# ARTICLE TYPE

# Optimal control of a solar sail

Alesia Herasimenka*<fnref n="1" /> | Lamberto Dell'Elce<fnref n="2" /> | Jean-Baptiste Caillau<fnref n="1" /> | Jean-Baptiste Pomet<fnref n="2" />

## Correspondence

*Email: alesia.herasimenka@univ-cotedazur.fr

## Present Address

This is sample for present address text this is sample for present address text

## Summary

Optimal control of a solar sail orbiting around a fixed center of mass is considered. The sail is modelled as a plane surface with two sides having similar optical properties. The control is assumed to be the attitude of the sail and is represented as an element of the projective plane. Mapping this plane into the three dimensional ambient space to compute the actual force generated by a given attitude results in a generally non-convex set of admissible control values. A suitable convex relaxation is introduced to study the optimality system associated with maximising the change of the sail orbital parameters in a given direction along one orbit. Existence for both the relaxed and the original problems are deduced. Building on previous works<fnref n="1" />,<fnref n="2" />,³, a refined analysis of the control structure is given, proving rigorously that it contains only switchings, and a global bound of the number of switchings is obtained. In order to compute effective solutions, a multiple shooting approach is retained that is well suited for systems with a sequence of switchings between various control subarcs. Three issues are addressed: initialisation of shooting by means a tailored semi-definite relaxation that takes advantage of good convergence properties of modern convex optimisation algorithms; changes in the control structure that are accommodated by coupling shooting with differential continuation; implicit character of the Hamiltonian maximisation when using Pontrjagin maximum principle that is taken care of by incorporating the associated equation for the dynamical feedback into the shooting procedure. The method is illustrated on an example from NASA for which the sail inclination is optimally changed over one orbital period.

## KEYWORDS:

solar sailing, non-convex constraints, SDP relaxation, multiple shooting, differential continuation

# 1 | INTRODUCTION

## 1.1 | Solar radiation pressure

Solar sails are satellites that use solar radiation pressure (SRP) as propulsion for orbital manœuvres. SRP is caused by the interaction between photons and the surface of the sail, and its magnitude depends on the distance between the Sun and the sail, $r$. Specifically, denote by $\Phi_{SR} = 1367 \text{ W m}^{-2}$ the solar flux at $r_{\oplus} = 1 \text{ AU}$ distance (which is the average Sun-Earth distance), and by $c$ the speed of light; then, a simple model for the SRP is given by:⁴, Chap. 3

$$P_{SR} = \frac{\Phi_{SR}}{c} \left( \frac{r_{\oplus}}{r} \right)^2. \qquad (1)$$

<footnote n="1">¹Université Côte d'Azur, CNRS, Inria, LJAD, Nice, France</footnote>

<footnote n="2">²Inria, Sophia Antipolis, France</footnote>
In this paper, similarly to$^5$, we consider a flat sail with surface A and mass m. Different optical and geometrical properties have impact on the resulting SRP force, which has components of the incoming, reflected, and thermal radiations, namely $f_a$, $f_r$, and $f_e$. Moreover, the reflected radiation has specular and diffuse contributions, $f_{rs}$ and $f_{ru}$, respectively. Each force component has different magnitude and direction, that can be identified through the Sun-sail direction, denoted as $\hat{s}$ and the unit vector normal to the sail having a positive component along $\hat{s}$, $\hat{n}$. In this paper we assume that both sides of the sail have the same optical properties, so only the (non-oriented) direction of the normal to the plane representing the plane will describe its attitude. We also assume that it is possible to control the attitude, and that the actual control will be the force generated by this attitude (see (6)). In this model, $\hat{n}$ belongs the projective plane $RP^2$ that one can describe as the union of one open hemisphere (whose axis is $\hat{s}$) with a circle whose antipodal points are identified (subset of directions perpendicular to $\hat{s}$). Fixing some basis $\hat{e}_1, \hat{e}_2$ of $\{\hat{s}\}^\perp$ in order that $(\hat{s}, \hat{e}_1, \hat{e}_2)$ is a direct orthogonal frame, one defines coordinates $(\beta, \delta) \in (-\pi/2, \pi/2) \times \mathbf{R}$ for $\hat{n}$ in the open hemisphere part setting as usual (see Figure 1)

$$ \hat{n} = \sin \beta(\cos \delta \hat{e}_1 + \sin \delta \hat{e}_2) + \cos \beta \hat{s}. $$

This is not a chart<fnref n="1" /> as no $\delta$, even restricted to $\mathbf{R}/\pi\mathbf{Z}$, can be uniquely associated with the direction $\hat{s}$. (See also Remark 1.) The angle $\beta$ is the so-called solar-sail cone angle. As shown in Fig. 1, let us introduce the direction of specular reflection given by $\hat{\xi}$, and the unit tangent vector $\hat{t}$ lying in the plane generated by $\hat{s}$ and $\hat{n}$. The force due to the incoming radiation, $f_a$, is directed along $\hat{s}$. The force provided by the specularly reflected radiation, $f_{rs}$, points along $\hat{\xi}$ and is caused by photons that are reflected symmetrically with respect to the normal of the sail, thus yielding an exchange of momentum. Diffuse reflection stems from the sail surface roughness, which causes photons to be uniformly reflected in all directions, yielding a component of the force toward the direction normal to the sail, $\hat{n}$. Finally, as the absorbed photons are re-radiated in all directions, the force $f_e$ is generated, which is orthogonal to the sail surface and points again along $\hat{n}$. We follow <fnref n="1" />, Chap. ² and express the unit vectors $\hat{s}$ and $\hat{\xi}$ in terms of $\hat{n}$ and $\hat{t}$,

$$ \hat{s} = \cos \beta \hat{n} + \sin \beta \hat{t}, \quad \hat{\xi} = \cos \beta \hat{n} - \sin \beta \hat{t}, \qquad (2) $$

so the forces presented above can be expressed as: ⁶

$$ \begin{aligned} f_a &= \varepsilon \cos \beta \hat{s} = \varepsilon \cos \beta(\cos \beta \hat{n} + \sin \beta \hat{t}), \\ f_{rs} &= \varepsilon \rho s \cos \beta \hat{\xi} = \varepsilon \rho s \cos \beta(\cos \beta \hat{n} - \sin \beta \hat{t}), \\ f_{ru} &= \varepsilon B_f \rho (1-s) \cos \beta \hat{n}, \\ f_e &= \varepsilon (1-\rho) \frac{\varepsilon_f B_f - \varepsilon_b B_b}{\varepsilon_b + \varepsilon_f} \cos \beta \hat{n}. \end{aligned} \qquad (3) $$

In (3), $\varepsilon$ is a constant parameter equal to $AP_{SR} m^{-1}$, which combines optical and physical parameters of the sail and has small magnitude, stressing the fact that the overall force is small; $\rho \in [0, 1]$ is the fraction of reflected radiation to total amount of radiation illuminating the sail; $s \in [0, 1]$ is the fraction of specularly reflected radiation to total reflected radiation; $\varepsilon_b$ and $\varepsilon_f$ are the back and front surface emissivity coefficients, respectively, and $B_b$ and $B_f$ are back and front non-Lambertian coefficients, respectively. The SRP force is the sum of these four contributions:

$$ f_{SRP} = f_a + f_{rs} + f_{ru} + f_e. \qquad (4) $$

Let us use the notations used in: ²

$$ b_1 = 1 - \rho s, \quad b_2 = 2\rho s, \quad b_3 = B_f \rho (1 - s) + (1 - \rho) \frac{\varepsilon_f B_f - \varepsilon_b B_b}{\varepsilon_f + \varepsilon_b}. $$

With these, the SRP force reads (note that because of the radial symmetry, its norm is independent of $\delta$)

$$ f_{SRP} = \varepsilon \cos \beta \begin{bmatrix} (b_1 + b_2) \cos^2 \beta + b_3 \cos \beta + b_1 \sin^2 \beta \\ b_2 \sin \beta \cos \beta \sin \delta + b_3 \sin \beta \sin \delta \\ b_2 \sin \beta \cos \beta \cos \delta + b_3 \sin \beta \cos \delta \end{bmatrix}. \qquad (5) $$

**Remark 1.** In our model, the magnitude of the SRP is continuous with respect to $\hat{n}$, going to zero when the Sun direction is contained in the sail plane (orthogonality of $\hat{s}$ and $\hat{n}$), but its direction is not: when going through $\hat{s} \perp \hat{n}$, the illuminated side of the sail (a thickless 2D object embedded into 3D space) is changed and the orientation of $\hat{n}$ is changed to opposite $(\beta = \pm \pi/2$ being changed to $-\beta$, still defining the same direction—a perpendicular to $\hat{s}$—in the projective plane). The resulting force, going

<footnote n="1">¹One actually retrieves the universal cover of the pointed open hemisphere by restricting to $(\beta, \delta)$ in $(0, \pi/2) \times \mathbf{R}$.</footnote>
In this paper, similarly to$^5$, we consider a flat sail with surface A and mass m. Different optical and geometrical properties have impact on the resulting SRP force, which has components of the incoming, reflected, and thermal radiations, namely $f_a$, $f_r$, and $f_e$. Moreover, the reflected radiation has specular and diffuse contributions, $f_{rs}$ and $f_{ru}$, respectively. Each force component has different magnitude and direction, that can be identified through the Sun-sail direction, denoted as $\hat{s}$ and the unit vector normal to the sail having a positive component along $\hat{s}$, $\hat{n}$. In this paper we assume that both sides of the sail have the same optical properties, so only the (non-oriented) direction of the normal to the plane representing the plane will describe its attitude. We also assume that it is possible to control the attitude, and that the actual control will be the force generated by this attitude (see (6)). In this model, $\hat{n}$ belongs the projective plane $RP^2$ that one can describe as the union of one open hemisphere (whose axis is $\hat{s}$) with a circle whose antipodal points are identified (subset of directions perpendicular to $\hat{s}$). Fixing some basis $\hat{e}_1, \hat{e}_2$ of $\{\hat{s}\}^\perp$ in order that $(\hat{s}, \hat{e}_1, \hat{e}_2)$ is a direct orthogonal frame, one defines coordinates $(\beta, \delta) \in (-\pi/2, \pi/2) \times \mathbf{R}$ for $\hat{n}$ in the open hemisphere part setting as usual (see Figure 1)

$$ \hat{n} = \sin \beta(\cos \delta \hat{e}_1 + \sin \delta \hat{e}_2) + \cos \beta \hat{s}. $$

This is not a chart<fnref n="1" /> as no $\delta$, even restricted to $\mathbf{R}/\pi\mathbf{Z}$, can be uniquely associated with the direction $\hat{s}$. (See also Remark 1.) The angle $\beta$ is the so-called solar-sail cone angle. As shown in Fig. 1, let us introduce the direction of specular reflection given by $\hat{\xi}$, and the unit tangent vector $\hat{t}$ lying in the plane generated by $\hat{s}$ and $\hat{n}$. The force due to the incoming radiation, $f_a$, is directed along $\hat{s}$. The force provided by the specularly reflected radiation, $f_{rs}$, points along $\hat{\xi}$ and is caused by photons that are reflected symmetrically with respect to the normal of the sail, thus yielding an exchange of momentum. Diffuse reflection stems from the sail surface roughness, which causes photons to be uniformly reflected in all directions, yielding a component of the force toward the direction normal to the sail, $\hat{n}$. Finally, as the absorbed photons are re-radiated in all directions, the force $f_e$ is generated, which is orthogonal to the sail surface and points again along $\hat{n}$. We follow <fnref n="1" />, Chap. ² and express the unit vectors $\hat{s}$ and $\hat{\xi}$ in terms of $\hat{n}$ and $\hat{t}$,

$$ \hat{s} = \cos \beta \hat{n} + \sin \beta \hat{t}, \quad \hat{\xi} = \cos \beta \hat{n} - \sin \beta \hat{t}, \qquad (2) $$

so the forces presented above can be expressed as: ⁶

$$ \begin{aligned} f_a &= \varepsilon \cos \beta \hat{s} = \varepsilon \cos \beta(\cos \beta \hat{n} + \sin \beta \hat{t}), \\ f_{rs} &= \varepsilon \rho s \cos \beta \hat{\xi} = \varepsilon \rho s \cos \beta(\cos \beta \hat{n} - \sin \beta \hat{t}), \\ f_{ru} &= \varepsilon B_f \rho (1-s) \cos \beta \hat{n}, \\ f_e &= \varepsilon (1-\rho) \frac{\varepsilon_f B_f - \varepsilon_b B_b}{\varepsilon_b + \varepsilon_f} \cos \beta \hat{n}. \end{aligned} \qquad (3) $$

In (3), $\varepsilon$ is a constant parameter equal to $AP_{SR}m^{-1}$, which combines optical and physical parameters of the sail and has small magnitude, stressing the fact that the overall force is small; $\rho \in [0, 1]$ is the fraction of reflected radiation to total amount of radiation illuminating the sail; $s \in [0, 1]$ is the fraction of specularly reflected radiation to total reflected radiation; $\varepsilon_b$ and $\varepsilon_f$ are the back and front surface emissivity coefficients, respectively, and $B_b$ and $B_f$ are back and front non-Lambertian coefficients, respectively. The SRP force is the sum of these four contributions:

$$ f_{SRP} = f_a + f_{rs} + f_{ru} + f_e. \qquad (4) $$

Let us use the notations used in: ²

$$ b_1 = 1 - \rho s, \quad b_2 = 2\rho s, \quad b_3 = B_f \rho (1 - s) + (1 - \rho) \frac{\varepsilon_f B_f - \varepsilon_b B_b}{\varepsilon_f + \varepsilon_b}. $$

With these, the SRP force reads (note that because of the radial symmetry, its norm is independent of $\delta$)

$$ f_{SRP} = \varepsilon \cos \beta \begin{bmatrix} (b_1 + b_2) \cos^2 \beta + b_3 \cos \beta + b_1 \sin^2 \beta \\ b_2 \sin \beta \cos \beta \sin \delta + b_3 \sin \beta \sin \delta \\ b_2 \sin \beta \cos \beta \cos \delta + b_3 \sin \beta \cos \delta \end{bmatrix}. \qquad (5) $$

**Remark 1.** In our model, the magnitude of the SRP is continuous with respect to $\hat{n}$, going to zero when the Sun direction is contained in the sail plane (orthogonality of $\hat{s}$ and $\hat{n}$), but its direction is not: when going through $\hat{s} \perp \hat{n}$, the illuminated side of the sail (a thickless 2D object embedded into 3D space) is changed and the orientation of $\hat{n}$ is changed to opposite $(\beta = \pm \pi/2$ being changed to $-\beta$, still defining the same direction—a perpendicular to $\hat{s}$—in the projective plane). The resulting force, going

<footnote n="1">¹One actually retrieves the universal cover of the pointed open hemisphere by restricting to $(\beta, \delta)$ in $(0, \pi/2) \times \mathbf{R}$.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1 Components of the SRP force and orientation angles $\beta$ and $\delta$ of the solar sail.</figcaption></figure>

to zero in such cases, is continuous but not smooth. This singularity is inherent to the modeling and would be removed in a more realistic approach describing the sail as a genuine 3D object. This lack of smoothness is nonetheless not crucial here since, as will be clear from the optimality analysis in Section 2, an optimal force will have discontinuities, being either zero or with $\beta \in (-\beta^*, \beta^*)$ and $0 < \beta^* < \pi/2$ (if we exclude the ideal case for which $\beta^* = \pi/2$). So flips of illuminated side will not be encountered.

## 1.2 | Parametrisation of the control set

Controlling the sail attitude, i.e. the normal vector $\hat{n}$, allows to change the direction and magnitude of the resulting SRP. A reliable inference of optical coefficients is indeed mandatory to accurately estimate the mapping $f_{SRP}$. To carry out our analysis, solar sail dynamics is conveniently modeled as a nonlinear control-affine system (see Section 1.3), where the control variable is homogeneous to the renormalized force. The control set $U \subset \mathbb{R}^3$ is so given by:

$$ U = \left\{ u = \frac{f_{SRP}(\hat{n})}{\varepsilon} \in \mathbb{R}^3, \hat{n} \in \mathbb{R}^3 \right\}. \qquad (6) $$

It is a closed surface of revolution with axis $\hat{s}$ that bounds a domain of $\mathbb{R}^3$. Up to the scaling factor $\varepsilon$, this surface is the image of the projective plane by the mapping $f_{SRP}$. As a result, the surface minus the origin is parametrised by $(\beta, \delta)$ in $(-\pi/2, \pi/2) \times \mathbb{R}$, while the origin corresponds to the image of the circle of directions orthogonal to $\hat{s}$. In the rest of the paper, the control is defined by these two angles, $\beta$ and $\delta$, and we shall see that while $\delta$ is easily retrieved from Pontrjagin maximum principle, the computation of $\beta$ for an optimal control is more intricate. For a given pair $(\beta, \delta)$, the actual force is obtained through (5). Figure 2 shows the intersection of $U$ with the plane generated by $\hat{n}$ and $\hat{s}$ for various optical properties. The resulting curve now bounds a two-dimensional domain, non-convex unless the sail is ideal ($\rho = s = 1$), so the sail is perfectly reflective, a first extreme case). In addition to ideal sails, another extreme case can be identified: contrary to perfectly reflective sails, perfectly absorptive surfaces are the worst-case scenario ($\rho = 0$, $f_e$ neglected) because SRP is systematically parallel to $\hat{s}$, as shown in red in Figure 2 top left. Although sails are designed to be as close to ideal as possible, partial absorption of the energy is unavoidable in real-life applications and optical properties exhibit degradation in time. Hence, the fraction of reflected radiation decreases with the lifetime of the satellite as discussed in$^{7,8}$. Another significant optical coefficient that plays a crucial role in determining the control set is the specular reflectivity, denoted by $s$, as illustrated in Figure 2 top right. The red set represents a sail with a rough surface (corresponding to $s = 0$), where light is diffusely reflected in all directions. This results in a smaller control set, even for a high reflectivity value of $\rho = 0.8$. In Figure 2, the impact of other optical coefficients, namely emissivity and non-Lambertian coefficients, is shown. Although these coefficients are small, they have a noticeable effect on the shape and convexity of the control set. (The ideal solar sail is the only exception, exhibiting a vertical tangency at its vertex.) While
<figure><img src="image_2.png" /><figcaption>Figure 1 Components of the SRP force and orientation angles β and δ of the solar sail.</figcaption></figure>

to zero in such cases, is continuous but not smooth. This singularity is inherent to the modeling and would be removed in a more realistic approach describing the sail as a genuine 3D object. This lack of smoothness is nonetheless not crucial here since, as will be clear from the optimality analysis in Section 2, an optimal force will have discontinuities, being either zero or with $\beta \in (-\beta^*, \beta^*)$ and $0 < \beta^* < \pi/2$ (if we exclude the ideal case for which $\beta^* = \pi/2$). So flips of illuminated side will not be encountered.

## 1.2 | Parametrisation of the control set

Controlling the sail attitude, i.e. the normal vector $\hat{n}$, allows to change the direction and magnitude of the resulting SRP. A reliable inference of optical coefficients is indeed mandatory to accurately estimate the mapping $f_{SRP}$. To carry out our analysis, solar sail dynamics is conveniently modeled as a nonlinear control-affine system (see Section 1.3), where the control variable is homogeneous to the renormalized force. The control set $U \subset \mathbb{R}^3$ is so given by:

$$U = \left\{ u = \frac{f_{SRP}(\hat{n})}{\varepsilon} \in \mathbb{R}^3, \hat{n} \in \mathbb{R}^3 \right\}. \qquad (6)$$

It is a closed surface of revolution with axis $\hat{s}$ that bounds a domain of $\mathbb{R}^3$. Up to the scaling factor $\varepsilon$, this surface is the image of the projective plane by the mapping $f_{SRP}$. As a result, the surface minus the origin is parametrised by $(\beta, \delta)$ in $(-\pi/2, \pi/2) \times \mathbb{R}$, while the origin corresponds to the image of the circle of directions orthogonal to $\hat{s}$. In the rest of the paper, the control is defined by these two angles, $\beta$ and $\delta$, and we shall see that while $\delta$ is easily retrieved from Pontrjagin maximum principle, the computation of $\beta$ for an optimal control is more intricate. For a given pair $(\beta, \delta)$, the actual force is obtained through (5). Figure 2 shows the intersection of $U$ with the plane generated by $\hat{n}$ and $\hat{s}$ for various optical properties. The resulting curve now bounds a two-dimensional domain, non-convex unless the sail is ideal ($\rho = s = 1$), so the sail is perfectly reflective, a first extreme case). In addition to ideal sails, another extreme case can be identified: contrary to perfectly reflective sails, perfectly absorptive surfaces are the worst-case scenario ($\rho = 0$, $f_e$ neglected) because SRP is systematically parallel to $\hat{s}$, as shown in red in Figure 2 top left. Although sails are designed to be as close to ideal as possible, partial absorption of the energy is unavoidable in real-life applications and optical properties exhibit degradation in time. Hence, the fraction of reflected radiation decreases with the lifetime of the satellite as discussed in$^{7,8}$. Another significant optical coefficient that plays a crucial role in determining the control set is the specular reflectivity, denoted by $s$, as illustrated in Figure 2 top right. The red set represents a sail with a rough surface (corresponding to $s = 0$), where light is diffusely reflected in all directions. This results in a smaller control set, even for a high reflectivity value of $\rho = 0.8$. In Figure 2, the impact of other optical coefficients, namely emissivity and non-Lambertian coefficients, is shown. Although these coefficients are small, they have a noticeable effect on the shape and convexity of the control set. (The ideal solar sail is the only exception, exhibiting a vertical tangency at its vertex.) While
<figure><img src="image_3.png" /></figure>

**Figure 2** Deformation of the control set with respect to the various optical properties of the sail. Top left: control sets for different reflectivity coefficients and $s = 1$. The blue sets correspond to the perfectly reflective sail, referring to the ideal sail model, when the red set illustrates forces generated by a perfectly absorptive sail, often called a cannonball model. Top right: control sets for different specular reflectivity coefficients and $\rho = 0.8$. As shown, the plots of specular reflectivity plays a major role in the resulting shape of the control set. Bottom left: control sets for different front emissivity coefficients. Bottom right: control sets for different front non-Lambertian coefficients.

most control sets have a horizontal tangency at their vertices, a more comprehensive analysis is required to fully understand the dependence of the control set geometry on different optical coefficients of the sail surface.

## 1.3 | Equations of motion

The following assumptions are introduced:

1. Orbital period of the sail is much smaller than the one of the heliocentric orbit of the attractor, so that variations of the Sun direction $\hat{s}$ over a single orbit of the sail are neglected.

2. Solar eclipses are neglected. Introducing solar eclipses can be a way of improving the proposed algorithm.

3. Re-emitted radiation is neglected. In fact, this component of SRP can be reasonably regarded as a disturbance for control purposes.

Equations of motion are written in a set of Keplerian-like orbital elements, which leverages on the axial symmetry of the problem with respect to the Sun's direction. Namely, consider a reference frame $S$ with origin at the center of the planet, $\hat{X} = \hat{s}$, $\hat{Y}$ lies in the plane of the planet's orbit around the Sun and is orthogonal to $\hat{X}$, and $\hat{Z}$ is chosen to form a right-hand frame. Because this study focuses on short-time controllability (characteristic time is of the order of one orbital period), motion of this frame is neglected by virtue of the first assumption above. Figure 3 represents the vectors $h$, $e$ and $\hat{N}$, which denote the angular momentum, eccentricity and ascending node vectors, respectively. Let $\gamma_1, \gamma_2, \gamma_3$ be Euler angles orienting the eccentricity vector
<figure><img src="image_3.png" /></figure>

**Figure 2** Deformation of the control set with respect to the various optical properties of the sail. Top left: control sets for different reflectivity coefficients and $s = 1$. The blue sets correspond to the perfectly reflective sail, referring to the ideal sail model, when the red set illustrates forces generated by a perfectly absorptive sail, often called a cannonball model. Top right: control sets for different specular reflectivity coefficients and $\rho = 0.8$. As shown, the plots of specular reflectivity plays a major role in the resulting shape of the control set. Bottom left: control sets for different front emissivity coefficients. Bottom right: control sets for different front non-Lambertian coefficients.

most control sets have a horizontal tangency at their vertices, a more comprehensive analysis is required to fully understand the dependence of the control set geometry on different optical coefficients of the sail surface.

## 1.3 | Equations of motion

The following assumptions are introduced:

1. Orbital period of the sail is much smaller than the one of the heliocentric orbit of the attractor, so that variations of the Sun direction $\hat{s}$ over a single orbit of the sail are neglected.

2. Solar eclipses are neglected. Introducing solar eclipses can be a way of improving the proposed algorithm.

3. Re-emitted radiation is neglected. In fact, this component of SRP can be reasonably regarded as a disturbance for control purposes.

Equations of motion are written in a set of Keplerian-like orbital elements, which leverages on the axial symmetry of the problem with respect to the Sun's direction. Namely, consider a reference frame $S$ with origin at the center of the planet, $\hat{X} = \hat{s}$, $\hat{Y}$ lies in the plane of the planet's orbit around the Sun and is orthogonal to $\hat{X}$, and $\hat{Z}$ is chosen to form a right-hand frame. Because this study focuses on short-time controllability (characteristic time is of the order of one orbital period), motion of this frame is neglected by virtue of the first assumption above. Figure 3 represents the vectors $h$, $e$ and $\hat{N}$, which denote the angular momentum, eccentricity and ascending node vectors, respectively. Let $\gamma_1, \gamma_2, \gamma_3$ be Euler angles orienting the eccentricity vector
<figure><img src="image_4.png" /><figcaption>**Figure 3** Euler angles $γ_i$ orienting the orbit according to a $\hat{X}-\hat{Y}-\hat{X}$ rotation with respect to the reference frame $S$. Here, $h$ and $e$ denote the angular momentum and eccentricity vectors.</figcaption></figure>

according to a $X-Y-X$ rotation as depicted in Figure 3, so that $γ_2$ is the angle between the angular momentum of the orbit and the Sun direction, and $a, e$, and $f$ be semi-major axis, eccentricity and true anomaly, respectively. The motion of slow elements, $I = (γ_1, γ_2, γ_3, a, e) \in \mathcal{M}$, where $\mathcal{M}$ is the configuration manifold, and where the fast angle $f$ is governed by

$$ \frac{dI}{dt} = \varepsilon \sqrt{\frac{a(1-e^2)}{\mu}} G_0(I, f) R(I, f) u, \quad (7) $$

$$ \frac{df}{dt} = \omega(I, f) + \varepsilon F(I, f) R(I, f) u, $$

where $\mu$ is the gravitational constant of the Sun, and where components of $u$ are in the reference frame $S$, $R(I, f) = R_X(γ_3 + f)R_Y(γ_2)R_X(γ_1)$ is the rotation matrix from reference to local-vertical local-horizontal frames<fnref n="2" />,

$$ \omega(I, f) = \sqrt{\frac{\mu}{a(1-e^2)^3}} (1 + e \cos f)^2. \quad (8) $$

Both $F(I, f)$ and $G_0(I, f)$ can be deduced from Gauss variational equations (GVE) of classical elements, where

$$ F(I, f) = \sqrt{\frac{a(1-e^2)}{\mu}} \left[ \frac{\cos f}{e} - \frac{2+e \cos f \sin f}{1+e \cos f} \right], \quad (9) $$

$$ G_0(I, f) = \begin{bmatrix} 0 & 0 & \frac{\sin(\gamma_3 + f)}{\sin\gamma_2(1 + e \cos f)} \\ 0 & 0 & \frac{\cos(\gamma_3 + f)}{1 + e \cos f} \\ \frac{-\cos f}{e} & \frac{2 + e \cos f \sin f}{1 + e \cos f} & \frac{\cos(\gamma_3 + f)}{1 + e \cos f} \\ \frac{2 ae}{1 - e^2} \sin f & \frac{2 ae}{1 - e^2 (1 + e \cos f)} & 0 \\ \sin f & \frac{e \cos^2 f + 2 \cos f + e}{1 + e \cos f} & 0 \end{bmatrix}. \quad (10) $$

<footnote n="2">²Here, $R_A(f)$ denotes the rotation matrix of angle $f$ about the axis $\hat{A}$.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3** Euler angles $γ_i$ orienting the orbit according to a $\hat{X}-\hat{Y}-\hat{X}$ rotation with respect to the reference frame $S$. Here, $h$ and $e$ denote the angular momentum and eccentricity vectors.</figcaption></figure>

according to a $X-Y-X$ rotation as depicted in Figure 3, so that $γ_2$ is the angle between the angular momentum of the orbit and the Sun direction, and $a, e$, and $f$ be semi-major axis, eccentricity and true anomaly, respectively. The motion of slow elements, $I = (γ_1, γ_2, γ_3, a, e) \in M$, where $M$ is the configuration manifold, and where the fast angle $f$ is governed by

$$ \frac{dI}{dt} = \varepsilon \sqrt{\frac{a(1-e^2)}{\mu}} G_0(I, f) R(I, f) u, \quad (7) $$

$$ \frac{df}{dt} = \omega(I, f) + \varepsilon F(I, f) R(I, f) u, $$

where $\mu$ is the gravitational constant of the Sun, and where components of $u$ are in the reference frame $S$, $R(I, f) = R_X(γ_3 + f)R_Y(γ_2)R_X(γ_1)$ is the rotation matrix from reference to local-vertical local-horizontal frames<fnref n="2" />,

$$ \omega(I, f) = \sqrt{\frac{\mu}{a(1-e^2)^3}} (1 + e \cos f)^2. \quad (8) $$

Both $F(I, f)$ and $G_0(I, f)$ can be deduced from Gauss variational equations (GVE) of classical elements, where

$$ F(I, f) = \sqrt{\frac{a(1-e^2)}{\mu}} \left[ \frac{\cos f}{e} - \frac{2+e \cos f \sin f}{1+e \cos f} \right], \quad (9) $$

$$ G_0(I, f) = \begin{bmatrix} 0 & 0 & \frac{\sin(\gamma_3 + f)}{\sin\gamma_2(1 + e \cos f)} \\ 0 & 0 & \frac{\cos(\gamma_3 + f)}{1 + e \cos f} \\ \frac{-\cos f}{e} & \frac{2 + e \cos f \sin f}{1 + e \cos f} & \frac{\cos(\gamma_3 + f)}{1 + e \cos f} \\ \frac{2 ae}{1 - e^2} \sin f & \frac{2 ae}{1 - e^2} (1 + e \cos f) & 0 \\ \sin f & \frac{e \cos^2 f + 2 \cos f + e}{1 + e \cos f} & 0 \end{bmatrix}. \quad (10) $$

<footnote n="2">²Here, $R_A(f)$ denotes the rotation matrix of angle $f$ about the axis $\hat{A}$.</footnote>
## 1.4 | Summary of contributions and outline of the paper

In Section 2 we introduce the optimal control problem of interest in this paper: it does not consist in minimising time or consumption, as is customary in space mechanics, but consists instead in maximising the displacement of the sail in a given direction over one orbit. A characteristic feature of this problem is the non-convexity of its control set; we define a suitable convex relaxation and review existence and necessary conditions of the original problem, giving precise bounds of the number of switchings of the control. A tailored method of resolution that combines convex optimisation and shooting is presented in Section 3. While a bounded polyhedral cone allows to approximate the set of admissible controls, a semi-infinite formulation of the control problem is proposed. Using a sum-of-squares approach à la Nesterov, this problem is recast as an SDP (semidefinite program) on the cone of positive-semidefinite matrices. For such a program, there are efficient convex optimisation algorithms with guaranteed convergence properties. The resulting solution is expected to provide an initial guess for shooting. In particular, we rely on the SDP step to capture the switching structure of the optimal control on a bounded cone. Then the solution on this cone is connected to the solution on the original non convex control set by means of differential continuation, coupled with multiple shooting. As the structure of the control may change during continuation (appearance or disappearance of control subarcs), a callback procedure is used to monitor the process and restart it with an updated shooting function. An additional difficulty comes from the fact that there is no explicit expression of the control that maximises the Hamiltonian coming from Pontrjagin maximum principle applied to the problem. This issue is accommodated by incorporating the equation implicitly defining the control into the shooting procedure. The last section is devoted to the numerical treatment by the described approach of an example coming from the NASA Jet Propulsion Lab: a non-ideal solar sail is considered, and the target is to increase the inclination of the initial orbit. Convex optimisation turns to provide a very precise solution to initialise the continuation. A change of structure is nonetheless observed during the early steps of homotopy as the number of arcs drops from three 5 to 3, and the callback is used to eventually obtain the zero-bang-zero optimal control of the sail.

While the literature on solar sailing is now well established, including contributions in optimal control, the paper fills some gaps in the analysis of the problem; in particular, existence is proved by using an auxiliary convexification, and a detailed analysis of the extremals is provided: the control only has isolated switchings (no singular arc), and a global bound on the number of switchings is given. On the algorithmic side, a tailored convex approximation based on trigonometric polynomials is derived that is crucial to initialise robustly the shooting procedure. To the best of our knowledge, the numerical results obtained by coupling shooting and continuation are new. In contrast with other papers in the literature, optimal control strategies are obtained that maximise the considered cost (change of orbital elements in a prescribed direction) globally on an orbit and not only instantaneously. Moreover, two genuinely original techniques are introduced: implicit computation of the control (crucial for this application since the equation verified by one of the two control angles cannot be solved explicitly) on one had, callback technique monitor changes in the structure of the control (disappearance of a subarc) during differential continuation. These two approaches are discussed in the two last sections of the paper.

# 2 | CONTROL OVER ONE ORBITAL PERIOD

## 2.1 | Optimal control formulation

We are interested in moving solar sail in the desired direction after one orbital period. Therefore, it is interesting to rewrite System (7) in terms of displacement of the slow state elements, denoted $\delta I$. Therefore, its dynamics is given by

$$ \delta I' = \varepsilon G(I, f)u \quad (11) $$

where $'$ := $d/df$ and with

$$ G(I, f) := \frac{a(1 - e^2)^2}{\mu(1 + e \cos f)^2} G_0(I, f) R(I, f). $$

We emphasise that the independent variable is now the true anomaly, and angle taken as the new time of the system over one period: $f \in [0, 2\pi]$. As mentioned earlier, SRP has a very small magnitude, this is why it is usually considered as a perturbation. Thus, changes on slow variables are very small over one orbital period so that $I$ will be assumed constant in the rest of the paper. For the same reason, we also neglect the order one term in $\varepsilon$ in the dynamics of the fast variable in (7). The goal is to maximise the size of the displacement in a given direction fixed by a unit vector, $d_I$, so that the final value of $\delta I$ is parallel to $d_I$. This problem can be written in Mayer form as follows (note the simple form of the dynamics, given by an explicit integral,
## 1.4 | Summary of contributions and outline of the paper

In Section 2 we introduce the optimal control problem of interest in this paper: it does not consist in minimising time or consumption, as is customary in space mechanics, but consists instead in maximising the displacement of the sail in a given direction over one orbit. A characteristic feature of this problem is the non-convexity of its control set; we define a suitable convex relaxation and review existence and necessary conditions of the original problem, giving precise bounds of the number of switchings of the control. A tailored method of resolution that combines convex optimisation and shooting is presented in Section 3. While a bounded polyhedral cone allows to approximate the set of admissible controls, a semi-infinite formulation of the control problem is proposed. Using a sum-of-squares approach à la Nesterov, this problem is recast as an SDP (semidefinite program) on the cone of positive-semidefinite matrices. For such a program, there are efficient convex optimisation algorithms with guaranteed convergence properties. The resulting solution is expected to provide an initial guess for shooting. In particular, we rely on the SDP step to capture the switching structure of the optimal control on a bounded cone. Then the solution on this cone is connected to the solution on the original non convex control set by means of differential continuation, coupled with multiple shooting. As the structure of the control may change during continuation (appearance or disappearance of control subarcs), a callback procedure is used to monitor the process and restart it with an updated shooting function. An additional difficulty comes from the fact that there is no explicit expression of the control that maximises the Hamiltonian coming from Pontrjagin maximum principle applied to the problem. This issue is accommodated by incorporating the equation implicitly defining the control into the shooting procedure. The last section is devoted to the numerical treatment by the described approach of an example coming from the NASA Jet Propulsion Lab: a non-ideal solar sail is considered, and the target is to increase the inclination of the initial orbit. Convex optimisation turns to provide a very precise solution to initialise the continuation. A change of structure is nonetheless observed during the early steps of homotopy as the number of arcs drops from three 5 to 3, and the callback is used to eventually obtain the zero-bang-zero optimal control of the sail.

While the literature on solar sailing is now well established, including contributions in optimal control, the paper fills some gaps in the analysis of the problem; in particular, existence is proved by using an auxiliary convexification, and a detailed analysis of the extremals is provided: the control only has isolated switchings (no singular arc), and a global bound on the number of switchings is given. On the algorithmic side, a tailored convex approximation based on trigonometric polynomials is derived that is crucial to initialise robustly the shooting procedure. To the best of our knowledge, the numerical results obtained by coupling shooting and continuation are new. In contrast with other papers in the literature, optimal control strategies are obtained that maximise the considered cost (change of orbital elements in a prescribed direction) globally on an orbit and not only instantaneously. Moreover, two genuinely original techniques are introduced: implicit computation of the control (crucial for this application since the equation verified by one of the two control angles cannot be solved explicitly) on one had, callback technique monitor changes in the structure of the control (disappearance of a subarc) during differential continuation. These two approaches are discussed in the two last sections of the paper.

# 2 | CONTROL OVER ONE ORBITAL PERIOD

## 2.1 | Optimal control formulation

We are interested in moving solar sail in the desired direction after one orbital period. Therefore, it is interesting to rewrite System (7) in terms of displacement of the slow state elements, denoted $\delta I$. Therefore, its dynamics is given by

$$ \delta I' = \varepsilon G(I, f)u \quad (11) $$

where $'$ := d/df and with

$$ G(I, f) := \frac{a(1 - e^2)^2}{\mu(1 + e \cos f)^2} G_0(I, f) R(I, f). $$

We emphasise that the independent variable is now the true anomaly, and angle taken as the new time of the system over one period: $f \in [0, 2\pi]$. As mentioned earlier, SRP has a very small magnitude, this is why it is usually considered as a perturbation. Thus, changes on slow variables are very small over one orbital period so that $I$ will be assumed constant in the rest of the paper. For the same reason, we also neglect the order one term in $\varepsilon$ in the dynamics of the fast variable in (7). The goal is to maximise the size of the displacement in a given direction fixed by a unit vector, $d_I$, so that the final value of $\delta I$ is parallel to $d_I$. This problem can be written in Mayer form as follows (note the simple form of the dynamics, given by an explicit integral,
<figure><img src="image_5.png" /><figcaption>**Figure 4** Convexification of the control set (in blue) by a convex cone (in red). Left: definition of the convex cone $K_\alpha$, the bounded cone $\hat{K}_\alpha$ and the critical angle $\beta^*$. Right: view in 3D.</figcaption></figure>

as the right-hand side does not depend on $\delta I$ in our approximation):

$$ \max_{u(f) \in U} (\delta I(2\pi)|d_I) \quad \text{subject to} \quad \delta I' = \varepsilon \sum_{i=1}^{3} u_i G_i(I, f), \quad \delta I(0) = 0, \quad \delta I(2\pi) \text{ parallel to } d_I. \quad (12) $$

Hereabove, and in the rest of the paper, the notation $(\cdot|\cdot)$ stands for the standard Euclidean product, $(x|y) = \sum_i x_i y_i$. Building upon results in$^{9,10}$, we have access to an effective test (related to the convex SDP approximation discussed in Section 3) to check that it is indeed possible to move in the direction $d_I$ after one revolution. So we assume in the sequel that the problem is controllable.

## 2.2 | Existence and necessary conditions for optimality

We first consider the relaxation of (12) obtained by replacing the control set $U$ by its convex hull: $u(f) \in \text{conv}(U)$ (see Figure 4). As the control set is now compact and convex, and since we have assumed controllability using controls valued in $U \subset \text{conv}(U)$, Filippov theorem entails that

**Proposition 1.** The relaxed problem has a solution.

Let $K_\alpha$ be the convex cone generated both by $U$ and by its convex hull, $\alpha$ denoting the half-angle at the cone vertex. To formulate the necessary optimality conditions for the problem on $\text{conv}(U)$ we introduce the costate $p_{\delta I}$ of $\delta I$, a covector of dimension 5. The Hamiltonian associated with the dynamics is

$$ H(I, f, p_{\delta I}, u) = \varepsilon p_{\delta I} G(I, f) u. \quad (13) $$

Remember that $I$ is a constant, and note that the Hamiltonian does not depend on the state $\delta I$ because of the very simple form of the dynamics. (The ODE defines a mere quadrature, here.) Clearly, $p_{\delta I}$ is constant and transversality conditions write

$$ (p_{\delta I}|d_I) = -p^0 \|d_I\|^2 = -p^0 \quad (14) $$

where $p^0$ is the nonpositive multiplier associated with the cost. In particular, $p_{\delta I}$ is not zero, since otherwise both $p^0$ and $p_{\delta I}$ would vanish. By homogeneity in $(p^0, p_{\delta I})$ there are two cases: (i) the abnormal case ($p^0 = 0$) when $(p_{\delta I}|d_I) = 0$ and where one can normalise setting $\|p_{\delta I}\| = 1$; (ii) the normal case ($p^0 < 0$) when $(p_{\delta I}|d_I) > 0$ and where one can normalise setting $(p_{\delta I}|d_I) = 1$. Let us set $\psi := p_{\delta I} G(I, f)$.

**Lemma 1.** For any $I$, the matrix formed by $G(I, f)$ and $\partial G(I, f)/\partial f$ has maximum rank for all $f \in [0, 2\pi]$.
<figure><img src="image_5.png" /><figcaption>**Figure 4** Convexification of the control set (in blue) by a convex cone (in red). Left: definition of the convex cone $K_\alpha$, the bounded cone $\hat{K}_\alpha$ and the critical angle $\beta^*$. Right: view in 3D.</figcaption></figure>

as the right-hand side does not depend on $\delta I$ in our approximation):

$$ \max_{u(f) \in U} (\delta I(2\pi)|d_I) \quad \text{subject to} \quad \delta I' = \varepsilon \sum_{i=1}^{3} u_i G_i(I, f), \quad \delta I(0) = 0, \quad \delta I(2\pi) \text{ parallel to } d_I. \quad (12) $$

Hereabove, and in the rest of the paper, the notation $(\cdot|\cdot)$ stands for the standard Euclidean product, $(x|y) = \sum_i x_i y_i$. Building upon results in⁹,¹⁰, we have access to an effective test (related to the convex SDP approximation discussed in Section 3) to check that it is indeed possible to move in the direction $d_I$ after one revolution. So we assume in the sequel that the problem is controllable.

## 2.2 | Existence and necessary conditions for optimality

We first consider the relaxation of (12) obtained by replacing the control set $U$ by its convex hull: $u(f) \in \text{conv}(U)$ (see Figure 4). As the control set is now compact and convex, and since we have assumed controllability using controls valued in $U \subset \text{conv}(U)$, Filippov theorem entails that

**Proposition 1.** The relaxed problem has a solution.

Let $K_\alpha$ be the convex cone generated both by $U$ and by its convex hull, $\alpha$ denoting the half-angle at the cone vertex. To formulate the necessary optimality conditions for the problem on $\text{conv}(U)$ we introduce the costate $p_{\delta I}$ of $\delta I$, a covector of dimension 5. The Hamiltonian associated with the dynamics is

$$ H(I, f, p_{\delta I}, u) = \varepsilon p_{\delta I} G(I, f) u. \quad (13) $$

Remember that $I$ is a constant, and note that the Hamiltonian does not depend on the state $\delta I$ because of the very simple form of the dynamics. (The ODE defines a mere quadrature, here.) Clearly, $p_{\delta I}$ is constant and transversality conditions write

$$ (p_{\delta I} | d_I) = -p^0 \|d_I\|^2 = -p^0 \quad (14) $$

where $p^0$ is the nonpositive multiplier associated with the cost. In particular, $p_{\delta I}$ is not zero, since otherwise both $p^0$ and $p_{\delta I}$ would vanish. By homogeneity in $(p^0, p_{\delta I})$ there are two cases: (i) the abnormal case ($p^0 = 0$) when $(p_{\delta I}|d_I) = 0$ and where one can normalise setting $\|p_{\delta I}\| = 1$; (ii) the normal case ($p^0 < 0$) when $(p_{\delta I}|d_I) > 0$ and where one can normalise setting $(p_{\delta I}|d_I) = 1$. Let us set $\psi := p_{\delta I} G(I, f)$.

**Lemma 1.** For any $I$, the matrix formed by $G(I, f)$ and $\partial G(I, f)/\partial f$ has maximum rank for all $f \in [0, 2\pi]$.
*Proof.* This computation is actually equivalent to the rank condition that can be verified in terms of Lie brackets (and, e.g., Cartesian coordinates) in¹¹ (check Lemma 1). □

As a result, the zeros of the dimension three covector $\psi$ (as a function of the true anomaly $f$) are isolated on $[0, 2\pi]$. Indeed, the previous lemma implies that $\psi$ and $d\psi/df$ cannot vanish simultaneously as then, $p_{\delta I}$ would be orthogonal to all columns of $G(I, f)$ and of its derivative, so $p_{\delta I}$ would be zero (a contradiction). So there are only finitely many such zeros on $[0, 2\pi]$, defining a locus of codimension greater than one in the $(I, f)$ space. For the sake of simplicity, we assume in the sequel that $\psi$ actually never vanishes. For a detailed discussion on the associated singularities of the dynamics, see¹².

The polar cone $K_{\alpha}^0$ is the set of directions having a nonpositive scalar product with those in $K_{\alpha}$. The drop-shaped curve obtained when intersecting the control set with a plane is parametrised by the angle $\beta$ alone, and we denote $\beta^* \in (0, \pi/2)$ the parameter associated with the tangency point of this curve with its conical hull (see Figure 4). Thus, angle $\alpha$ of the convex cone can be retrieved from the critical angle $\beta^*$. In the sequel, we recall and complete the analysis from², providing precise bounds on the number of switchings on the control.

**Lemma 2.**² The critical angle $\beta^*$ is solution of

$$ \cos \beta^* = \frac{-b_1 b_3 - 2 b_2 b_3 + \sqrt{b_1^2 b_3^2 - 4 b_1 b_2 b_3^2 + 8 b_1^2 b_2^2 + 4 b_1 b_2^3}}{4 b_1 b_2 + 2 b_2^2}, $$

and is related to the cone angle $\alpha$ by the following relation:

$$ \tan \alpha = \frac{(b_3 + b_2 \cos \beta^*) \sin \beta^*}{b_1 + b_2 \cos^2 \beta^* + b_3 \cos \beta^*}. $$

**Proposition 2.** An optimal control $u$ verifies the following: (i) when $\psi$ belongs to the interior of $K_{\alpha}^0$, $u$ is zero; (ii) when $\psi$ does not belong to $K_{\alpha}^0$, the coordinates $(\beta, \delta)$ of the control verify the following relations:

$$ \psi_1 \sin \beta(b_1 + 3b_2 \cos^2 \beta + 2b_3 \cos \beta) = \sqrt{\psi_2^2 + \psi_3^2} (\cos^2 \beta(b_2 \cos \beta + b_3) - \sin^2 \beta(2b_2 \cos \beta + b_3)), \quad \beta \in (-\beta^*, \beta^*), \quad (15) $$

and

$$ \delta = \pi/2 - \arg(\psi_2 + i\psi_3) \mod \pi. \quad (16) $$

Moreover, any optimal control is made of finitely many subarcs corresponding to case (i) or (ii), and has at most 8 switchings (transverse contacts with $\partial K_{\alpha}^0$) over one period.

*Proof.* According to Pontrjagin's maximum principle (PMP) and to the expression (13) of the Hamiltonian, for almost all true anomaly $f$ an optimal control must be a maximiser of the scalar product $(\psi|u)$ for $u$ in conv($U$). Clearly, when $\psi$ belongs to the interior of the polar cone of $K_{\alpha}$, this scalar product is negative for any nonzero $u$, so $u = 0$ is the only maximiser. Conversely, when $\psi$ belongs to the open complement of $K_{\alpha}^0$, maximisers must annihilate the gradient of the Hamiltonian with respect to the chosen coordinates of the control,

$$ \frac{\partial H}{\partial \beta} = 0, \quad \frac{\partial H}{\partial \delta} = 0, $$

which gives the expressions in alternative (ii) of the statement. See Figure 5 for the geometrical interpretation. Moreover, $\psi$ belongs to the boundary of $K_{\alpha}^0$ if and only if $\psi_1 \cos \alpha + \sqrt{\psi_2^2 + \psi_3^2} \sin \alpha = 0$, implying that

$$ \psi_1^2 \cos^2 \alpha - (\psi_2^2 + \psi_3^2) \sin^2 \alpha = 0. \quad (17) $$

Every component of $\psi$ is trigonometric in $f$, and this (nontrivial) equation results in a trigonometric polynomial of degree 4. As it has isolated zeros, there are finitely many zeros (at most eight, see Remark 3) defining isolated contacts with $\partial K_{\alpha}^0$. □

**Remark 2.** Equation (15) defines the optimal solution $\beta$ as an implicit function of $\psi$ and optical parameters of the sail, notably, $b_1, b_2, b_3$. This relation becomes explicit when considering an ideal solar sail model, what is widely used in the literature for different preliminary analysis design. In this case, $b_1 = b_3 = 0$, and $b_2 = 2$. Therefore, Equation (15) becomes:

$$ 6 \psi_1 \sin \beta \cos^2 \beta = \sqrt{\psi_2^2 + \psi_3^2} (2 \cos^2 \beta \cos \beta - 4 \sin^2 \beta \cos \beta), \quad \beta \in (-\beta^*, \beta^*), \quad (18) $$
*Proof.* This computation is actually equivalent to the rank condition that can be verified in terms of Lie brackets (and, e.g., Cartesian coordinates) in¹¹ (check Lemma 1). □

As a result, the zeros of the dimension three covector $\psi$ (as a function of the true anomaly $f$) are isolated on $[0, 2\pi]$. Indeed, the previous lemma implies that $\psi$ and $d\psi/df$ cannot vanish simultaneously as then, $p_{\delta I}$ would be orthogonal to all columns of $G(I, f)$ and of its derivative, so $p_{\delta I}$ would be zero (a contradiction). So there are only finitely many such zeros on $[0, 2\pi]$, defining a locus of codimension greater than one in the $(I, f)$ space. For the sake of simplicity, we assume in the sequel that $\psi$ actually never vanishes. For a detailed discussion on the associated singularities of the dynamics, see¹².

The polar cone $K_{\alpha}^{0}$ is the set of directions having a nonpositive scalar product with those in $K_{\alpha}$. The drop-shaped curve obtained when intersecting the control set with a plane is parametrised by the angle $\beta$ alone, and we denote $\beta^{*} \in (0, \pi/2)$ the parameter associated with the tangency point of this curve with its conical hull (see Figure 4). Thus, angle $\alpha$ of the convex cone can be retrieved from the critical angle $\beta^{*}$. In the sequel, we recall and complete the analysis from², providing precise bounds on the number of switchings on the control.

**Lemma 2.**² The critical angle $\beta^{*}$ is solution of

$$ \cos \beta^* = \frac{-b_1 b_3 - 2b_2 b_3 + \sqrt{b_1^2 b_3^2 - 4b_1 b_2 b_3^2 + 8b_1^2 b_2^2 + 4b_1 b_2^3}}{4b_1 b_2 + 2b_2^2}, $$

and is related to the cone angle $\alpha$ by the following relation:

$$ \tan \alpha = \frac{(b_3 + b_2 \cos \beta^*) \sin \beta^*}{b_1 + b_2 \cos^2 \beta^* + b_3 \cos \beta^*}. $$

**Proposition 2.** An optimal control $u$ verifies the following: (i) when $\psi$ belongs to the interior of $K_{\alpha}^{0}$, $u$ is zero; (ii) when $\psi$ does not belong to $K_{\alpha}^{0}$, the coordinates $(\beta, \delta)$ of the control verify the following relations:

$$ \psi_1 \sin \beta(b_1 + 3b_2 \cos^2 \beta + 2b_3 \cos \beta) = \sqrt{\psi_2^2 + \psi_3^2} (\cos^2 \beta(b_2 \cos \beta + b_3) - \sin^2 \beta(2b_2 \cos \beta + b_3)), \quad \beta \in (-\beta^*, \beta^*), \quad (15) $$

and

$$ \delta = \pi/2 - \arg(\psi_2 + i\psi_3) \mod \pi. \quad (16) $$

Moreover, any optimal control is made of finitely many subarcs corresponding to case (i) or (ii), and has at most 8 switchings (transverse contacts with $\partial K_{\alpha}^{0}$) over one period.

*Proof.* According to Pontrjagin's maximum principle (PMP) and to the expression (13) of the Hamiltonian, for almost all true anomaly $f$ an optimal control must be a maximiser of the scalar product $(\psi|u)$ for $u$ in conv($U$). Clearly, when $\psi$ belongs to the interior of the polar cone of $K_{\alpha}$, this scalar product is negative for any nonzero $u$, so $u = 0$ is the only maximiser. Conversely, when $\psi$ belongs to the open complement of $K_{\alpha}^{0}$, maximisers must annihilate the gradient of the Hamiltonian with respect to the chosen coordinates of the control,

$$ \frac{\partial H}{\partial \beta} = 0, \quad \frac{\partial H}{\partial \delta} = 0, $$

which gives the expressions in alternative (ii) of the statement. See Figure 5 for the geometrical interpretation. Moreover, $\psi$ belongs to the boundary of $K_{\alpha}^{0}$ if and only if $\psi_1 \cos \alpha + \sqrt{\psi_2^2 + \psi_3^2} \sin \alpha = 0$, implying that

$$ \psi_1^2 \cos^2 \alpha - (\psi_2^2 + \psi_3^2) \sin^2 \alpha = 0. \quad (17) $$

Every component of $\psi$ is trigonometric in $f$, and this (nontrivial) equation results in a trigonometric polynomial of degree 4. As it has isolated zeros, there are finitely many zeros (at most eight, see Remark 3) defining isolated contacts with $\partial K_{\alpha}^{0}$. □

**Remark 2.** Equation (15) defines the optimal solution $\beta$ as an implicit function of $\psi$ and optical parameters of the sail, notably, $b_1, b_2, b_3$. This relation becomes explicit when considering an ideal solar sail model, what is widely used in the literature for different preliminary analysis design. In this case, $b_1 = b_3 = 0$, and $b_2 = 2$. Therefore, Equation (15) becomes:

$$ 6 \psi_1 \sin \beta \cos^2 \beta = \sqrt{\psi_2^2 + \psi_3^2} (2 \cos^2 \beta \cos \beta - 4 \sin^2 \beta \cos \beta), \quad \beta \in (-\beta^*, \beta^*), \quad (18) $$
<figure><img src="image_6.png" /><figcaption>Figure 5 Geometrical interpretation of the PMP.</figcaption></figure>

leading to the classical well-known near-optimal steering law for an ideal sail¹:

$$ \beta = \tan^{-1} \left( \frac{-3|\psi_1| + \sqrt{9\psi_1^2 + \psi_2^2 + \psi_3^2}}{4(\psi_2^2 + \psi_3^2)} \right). \quad (19) $$

Remark 3. Roots of a trigonometric polynomial can be found using companion-matrix methods¹³. Consider the degree 4 polynomial

$$ \mathcal{T}(f) = \sum_{j=0}^{4} a_j \cos(jf) + \sum_{j=1}^{4} b_j \sin(jf). $$

Fourier-Frobenius companion matrix elements are

$$ B_{jk} = \begin{cases} \delta_{j,k-1}, & j = 1, \dots, 7, k = 1, \dots, 8, \\ (-1)^{k-1} \frac{h_{k-1}}{a_4 - ib_4}, & j = 8, k = 1, \dots, 8, \end{cases} \quad (20) $$

where $\delta_{jk}$ are the Kronecker functions such that $\delta_{jk} = 0$ if $j \neq k$ and $\delta_{jj} = 1$, and $h_k$ are

$$ h_k = \begin{cases} a_{4-k} + ib_{4-k}, & k = 0, \dots, 3, \\ 2a_0, & k = 4, \\ a_{k-4} - ib_{k-4}, & k = 5, \dots, 8. \end{cases} \quad (21) $$

The roots of $\mathcal{T}(f)$ are obtained from eigenvalues $z_k$ of the matrix defined in Eq. (20) as

$$ f_{k,m} = \arg(z_k) - i \log(|z_k|) \mod (2\pi), \quad k = 1, \dots, 8. $$

Real-valued roots of $\mathcal{T}(f)$ are such that $|z_k| = 1$. Therefore, this technique allows to find roots of the switch function and, thus, find out the structure of the solution for a given costate. It is important to stress that the trigonometric polynomial is of degree 4, which means that the switching function can have at most 8 roots.

**Corollary 1.** The original optimal control problem (12) has a solution.

*Proof.* The relaxed problem has at least one solution (Proposition 1), and any control solution actually belongs to $U$ by virtue of Proposition 2. Such controls must be optimal for the original problem, whence existence. □

# 3 | SOLUTION USING CONVEX OPTIMISATION AND CONTINUATION

## 3.1 | Convex approximation for a reliable initial guess

In order to use indirect shooting methods for solving optimal control problem, we need first a reliable initial guess for the costate $p_{\delta I}$. We propose an approximation by a convex mathematical program similar to the one used in⁹ for controllability check
<figure><img src="image_6.png" /><figcaption>Figure 5 Geometrical interpretation of the PMP.</figcaption></figure>

leading to the classical well-known near-optimal steering law for an ideal sail¹:

$$ \beta = \tan^{-1} \left( \frac{-3|\psi_1| + \sqrt{9\psi_1^2 + \psi_2^2 + \psi_3^2}}{4(\psi_2^2 + \psi_3^2)} \right). \quad (19) $$

Remark 3. Roots of a trigonometric polynomial can be found using companion-matrix methods¹³. Consider the degree 4 polynomial

$$ \mathcal{T}(f) = \sum_{j=0}^{4} a_j \cos(jf) + \sum_{j=1}^{4} b_j \sin(jf). $$

Fourier-Frobenius companion matrix elements are

$$ B_{jk} = \begin{cases} \delta_{j,k-1}, & j = 1, \dots, 7, k = 1, \dots, 8, \\ (-1)^{k-1} \frac{h_{k-1}}{a_4 - ib_4}, & j = 8, k = 1, \dots, 8, \end{cases} \quad (20) $$

where $\delta_{jk}$ are the Kronecker functions such that $\delta_{jk} = 0$ if $j \neq k$ and $\delta_{jj} = 1$, and $h_k$ are

$$ h_k = \begin{cases} a_{4-k} + ib_{4-k}, & k = 0, \dots, 3, \\ 2a_0, & k = 4, \\ a_{k-4} - ib_{k-4}, & k = 5, \dots, 8. \end{cases} \quad (21) $$

The roots of $\mathcal{T}(f)$ are obtained from eigenvalues $z_k$ of the matrix defined in Eq. (20) as

$$ f_{k,m} = \arg(z_k) - i \log(|z_k|) \mod (2\pi), \quad k = 1, \dots, 8. $$

Real-valued roots of $\mathcal{T}(f)$ are such that $|z_k| = 1$. Therefore, this technique allows to find roots of the switch function and, thus, find out the structure of the solution for a given costate. It is important to stress that the trigonometric polynomial is of degree 4, which means that the switching function can have at most 8 roots.

**Corollary 1.** The original optimal control problem (12) has a solution.

*Proof.* The relaxed problem has at least one solution (Proposition 1), and any control solution actually belongs to $U$ by virtue of Proposition 2. Such controls must be optimal for the original problem, whence existence. □

# 3 | SOLUTION USING CONVEX OPTIMISATION AND CONTINUATION

## 3.1 | Convex approximation for a reliable initial guess

In order to use indirect shooting methods for solving optimal control problem, we need first a reliable initial guess for the costate $p_{\delta I}$. We propose an approximation by a convex mathematical program similar to the one used in⁹ for controllability check
purposes. To this end, define the bounded cone $\hat{K}_{\alpha}$ obtained by truncating the $K_{\alpha}$ at its tangency points with $U$ (check Figure 6). This cone is bounded by a disk denoted $D_{\alpha}$. This new control set is a subset of the convex hull of $U$, in order that any solution of

$$ \max_{u(f) \in \hat{K}_{\alpha}} (\delta I(2\pi)|d_I|) \quad \text{subject to} \quad \delta I' = \varepsilon \sum_{i=1}^{3} u_i G_i(I, f), \quad \delta I(0) = 0, \quad \delta I(2\pi) \text{ parallel to } d_I, \qquad (22) $$

will define an admissible control for the convex relaxation of the original control problem. Note that existence holds for this new problem (Filippov again, as $\hat{K}_{\alpha}$ is convex and bounded) and that any solution will also have a bang-bang structure. A similar analysis to the one of Section 2.2 on conv($U$) indeed allows to prove the result below.

**Proposition 3.** An optimal control $u$ of problem (22) on $\hat{K}_{\alpha}$ verifies the following: (i) when $\psi$ belongs to the interior of $K_{\alpha}^0$, $u$ is zero; when $\psi$ does not belong to $K_{\alpha}^0$, (ii-a) the control is uniquely determined and belongs to the circle $\partial D_{\alpha}$, unless (ii-b) $\psi$ is colinear to the axis $\hat{s}$ of the cone $K_{\alpha}$ in which case the control still belongs to the $\partial D_{\alpha}$ but is not uniquely determined, as shown in Figure 6. Moreover, any optimal control is made of finitely many subarcs corresponding to case (i) or (ii-a) over one period.

*Proof.* As $\hat{K}_{\alpha}$ and $K_{\alpha}$ have the same polar cone, (i) is clear. Conversely, when $\psi$ belongs to the open complement of $K_{\alpha}^0$, the colinearity condition $\psi \wedge \hat{s} = 0$ reduces to checking a polynomial condition in $f$ and has only isolated zeros corresponding to case (ii-b). When $\psi$ is not colinear to $\hat{s}$, the unique maximiser of $(\psi|u)$ for $u$ in $\hat{K}_{\alpha}$ indeed belongs to the circle $\partial D_{\alpha}$, which is case (ii-a). □

This structure being analogous to that of solutions of the original problem, one hopes to retrieve a reasonable approximation to be used to initiate a differential continuation (see Section 3.2). In particular, we note that the original problem (12) on $U$ and problem (22) on $\hat{K}_{\alpha}$ share the same switching function associated with contacts with $\partial K_{\alpha}^0$ and given by (17).

<figure><img src="image_7.png" /><figcaption>**Figure 6** Left: approximation of the control set by a cone bounded at the points of tangency. Middle: optimal solution will be given by the controls on the boundary of the cone or a zero. Controls on the boundary of the cone are actually situated on a circle, illustrated here by two dots projected on the two-dimensional plane. Right: possible optimal solutions on the real control set.</figcaption></figure>

Consider the following discretisation of (22): the control set $\hat{K}_{\alpha}$ is approximated by a polyhedral cone $\hat{K}_{\alpha}^g \subset \hat{K}_{\alpha}$ generated as the convex hull of $g$ vertices $V_1, \dots, V_g$ chosen in $\partial \hat{K}_{\alpha}$, as shown on Figure 7. (Note that the 3D cone $\hat{K}_{\alpha}$ is not finitely generated.) Any control in $\hat{K}_{\alpha}^g$ is given by a bounded conical combination

$$ u(f) = \sum_{j=1}^{g} v_j(f)V_j, \quad v_j(f) \ge 0, \quad \sum_{j=1}^{g} v_j(f) \le 1, \quad f \in \mathbf{S}^1 \simeq \mathbf{R}/2\pi\mathbf{Z}, \quad j = 1, \dots, g. \qquad (23) $$

The functions $v_j$ are modeled using an $N$-dimensional basis of trigonometric polynomials, $\Phi(f) = (1, e^{if}, e^{2if}, \dots, e^{(N-1)if})$:

$$ v_j(f) = (\Phi(f) | c_j)_H \qquad (24) $$

with $c_j \in \mathbf{C}^N$ complex-valued coordinates of $v_j$ in $\Phi(f)$, and $(\cdot|\cdot)_H$ stands for the Hermitian product on $\mathbf{C}^N$. To enforce the positivity constraint, we leverage on the formalism of squared functional systems outlined in¹⁴. It allows to recast continuous positivity constraints as Linear Matrix Inequalities (LMI), that can be solved using convex optimisation. Actually, $\Phi(f)$ has a corresponding squared functional system given by $S^2(f) = \Phi(f)\Phi^H(f)$, with $\Phi^H(f)$ being the conjugate transpose of $\Phi(f)$. According to¹⁴, let us define a linear operator $\Lambda_H : \mathbf{C}^N \to \mathbf{C}^{N \times N}$ that maps coefficients of a polynomial in $\Phi(f)$ to its squared base, and its adjoint operator $\Lambda_H^* : \mathbf{C}^{N \times N} \to \mathbf{C}^N$ such that

$$ (Y|\Lambda_H(c))_H = (\Lambda_H^*(Y)|c)_H, \quad Y \in \mathbf{C}^{N \times N}. \qquad (25) $$
purposes. To this end, define the bounded cone $\hat{K}_\alpha$ obtained by truncating the $K_\alpha$ at its tangency points with $U$ (check Figure 6). This cone is bounded by a disk denoted $D_\alpha$. This new control set is a subset of the convex hull of $U$, in order that any solution of

$$ \max_{u(f) \in \hat{K}_\alpha} (\delta I(2\pi)|d_I|) \quad \text{subject to} \quad \delta I' = \epsilon \sum_{i=1}^{3} u_i G_i(I, f), \quad \delta I(0) = 0, \quad \delta I(2\pi) \text{ parallel to } d_I, \qquad (22) $$

will define an admissible control for the convex relaxation of the original control problem. Note that existence holds for this new problem (Filippov again, as $\hat{K}_\alpha$ is convex and bounded) and that any solution will also have a bang-bang structure. A similar analysis to the one of Section 2.2 on conv($U$) indeed allows to prove the result below.

**Proposition 3.** An optimal control $u$ of problem (22) on $\hat{K}_\alpha$ verifies the following: (i) when $\psi$ belongs to the interior of $K_\alpha^0$, $u$ is zero; when $\psi$ does not belong to $K_\alpha^0$, (ii-a) the control is uniquely determined and belongs to the circle $\partial D_\alpha$, unless (ii-b) $\psi$ is colinear to the axis $\hat{s}$ of the cone $K_\alpha$ in which case the control still belongs to the $\partial D_\alpha$ but is not uniquely determined, as shown in Figure 6. Moreover, any optimal control is made of finitely many subarcs corresponding to case (i) or (ii-a) over one period.

*Proof.* As $\hat{K}_\alpha$ and $K_\alpha$ have the same polar cone, (i) is clear. Conversely, when $\psi$ belongs to the open complement of $K_\alpha^0$, the colinearity condition $\psi \wedge \hat{s} = 0$ reduces to checking a polynomial condition in $f$ and has only isolated zeros corresponding to case (ii-b). When $\psi$ is not colinear to $\hat{s}$, the unique maximiser of $(\psi|u)$ for $u$ in $\hat{K}_\alpha$ indeed belongs to the circle $\partial D_\alpha$, which is case (ii-a). □

This structure being analogous to that of solutions of the original problem, one hopes to retrieve a reasonable approximation to be used to initiate a differential continuation (see Section 3.2). In particular, we note that the original problem (12) on $U$ and problem (22) on $\hat{K}_\alpha$ share the same switching function associated with contacts with $\partial K_\alpha^0$ and given by (17).

<figure><img src="image_7.png" /><figcaption>**Figure 6** Left: approximation of the control set by a cone bounded at the points of tangency. Middle: optimal solution will be given by the controls on the boundary of the cone or a zero. Controls on the boundary of the cone are actually situated on a circle, illustrated here by two dots projected on the two-dimensional plane. Right: possible optimal solutions on the real control set.</figcaption></figure>

Consider the following discretisation of (22): the control set $\hat{K}_\alpha$ is approximated by a polyhedral cone $\hat{K}_\alpha^g \subset \hat{K}_\alpha$ generated as the convex hull of $g$ vertices $V_1, \dots, V_g$ chosen in $\partial \hat{K}_\alpha$, as shown on Figure 7. (Note that the 3D cone $\hat{K}_\alpha$ is not finitely generated.) Any control in $\hat{K}_\alpha^g$ is given by a bounded conical combination

$$ u(f) = \sum_{j=1}^{g} v_j(f)V_j, \quad v_j(f) \ge 0, \quad \sum_{j=1}^{g} v_j(f) \le 1, \quad f \in \mathbf{S}^1 \simeq \mathbf{R}/2\pi\mathbf{Z}, \quad j = 1, \dots, g. \qquad (23) $$

The functions $v_j$ are modeled using an $N$-dimensional basis of trigonometric polynomials, $\Phi(f) = (1, e^{if}, e^{2if}, \dots, e^{(N-1)if})$:

$$ v_j(f) = (\Phi(f) | c_j)_H \qquad (24) $$

with $c_j \in \mathbf{C}^N$ complex-valued coordinates of $v_j$ in $\Phi(f)$, and $(\cdot|\cdot)_H$ stands for the Hermitian product on $\mathbf{C}^N$. To enforce the positivity constraint, we leverage on the formalism of squared functional systems outlined in¹⁴. It allows to recast continuous positivity constraints as Linear Matrix Inequalities (LMI), that can be solved using convex optimisation. Actually, $\Phi(f)$ has a corresponding squared functional system given by $S^2(f) = \Phi(f)\Phi^H(f)$, with $\Phi^H(f)$ being the conjugate transpose of $\Phi(f)$. According to¹⁴, let us define a linear operator $\Lambda_H : \mathbf{C}^N \to \mathbf{C}^{N \times N}$ that maps coefficients of a polynomial in $\Phi(f)$ to its squared base, and its adjoint operator $\Lambda_H^* : \mathbf{C}^{N \times N} \to \mathbf{C}^N$ such that

$$ (Y|\Lambda_H(c))_H = (\Lambda_H^*(Y)|c)_H, \quad Y \in \mathbf{C}^{N \times N}. \qquad (25) $$
<figure><img src="image_8.png" /><figcaption>Figure 7 Approximation of the bounded convex cone $\hat{K}_{\alpha}$ by the polyhedral cone $\hat{K}_{\alpha}^{g}$.</figcaption></figure>

The theory of squared functional systems states that, for the trigonometric polynomial $(\Phi(f)|c)$ to be non-negative, it is necessary and sufficient that there exists a Hermitian semidefinite positive matrix $Y \ge 0$, such that

$$ \Lambda_H^*(Y) = c. $$

The fact that there is no gap between nonnegativeness of the polynomials and reduction to a sum of squares (that can be expressed in SDP terms) is remarkable and peculiar to the univariate case, without restriction on the degree (actually, it remains true for bivariate trigonometric polynomials, see¹⁵; this could allow us to include the second angle, $\delta$, instead of using a polyhedral approximation). All in all, this approach permits to translate the semi-infinite constraint (25) into an SDP one that can be very efficiently treated by modern convex optimisation solvers. The operator $\Lambda_H^*$ is defined by means of Toeplitz matrices:

$$ \Lambda_H^*(Y) = \begin{bmatrix} \mathrm{tr}(Y|T_0) \\ \vdots \\ \mathrm{tr}(Y|T_{N-1}) \end{bmatrix} \qquad (26) $$

with $kl$-coefficients of $T_j$ such that

$$ T_0 = I, \\ T_j^{(k,l)} = \begin{cases} 2 & \text{if } k-l=j \\ 0 & \text{otherwise} \end{cases} \quad j=1, \dots, N-1. \qquad (27) $$

For an admissible control $u$ valued in $\hat{K}_{\alpha}^{g}$, one has

$$ \int_{0}^{2\pi} \sum_{i=1}^{3} u_i(f) G_i(I, f) \, df = \sum_{j=1}^{g} (L_j c_j + \bar{L}_j \bar{c}_j) \qquad (28) $$

with $L_j(I)$ in $C^{5 \times N}$ defined by

$$ L_j(I) = \frac{1}{2} \sum_{i=1}^{3} \int_{0}^{2\pi} V_{ij} G_i(I, f) \Phi^H(f) \, df, \qquad (29) $$

where $V_j = (V_{ij})_{i=1,...,3}$. We note that the components of $L_j(I)$ are Fourier coefficients of the function $\sum_{i=1}^{3} V_{ij} G_i(I, f)$. $L_j(I)$ are approximated using a Discrete Fourier Transform. Since vector fields $G_i$ are smooth, truncation of the series is justified by the fast decrease of the coefficients. Finally, for a control $u$ valued in $\hat{K}_{\alpha}^{g}$, coefficients $v_j$ are truncated Fourier series of order
<figure><img src="image_8.png" /><figcaption>Figure 7 Approximation of the bounded convex cone $\hat{K}_{\alpha}$ by the polyhedral cone $\hat{K}_{\alpha}^{g}$.</figcaption></figure>

The theory of squared functional systems states that, for the trigonometric polynomial $(\Phi(f)|c)$ to be non-negative, it is necessary and sufficient that there exists a Hermitian semidefinite positive matrix $Y \ge 0$, such that

$$ \Lambda_H^*(Y) = c. $$

The fact that there is no gap between nonnegativeness of the polynomials and reduction to a sum of squares (that can be expressed in SDP terms) is remarkable and peculiar to the univariate case, without restriction on the degree (actually, it remains true for bivariate trigonometric polynomials, see¹⁵; this could allow us to include the second angle, $\delta$, instead of using a polyhedral approximation). All in all, this approach permits to translate the semi-infinite constraint (25) into an SDP one that can be very efficiently treated by modern convex optimisation solvers. The operator $\Lambda_H^*$ is defined by means of Toeplitz matrices:

$$ \Lambda_H^*(Y) = \begin{bmatrix} \mathrm{tr}(Y|T_0) \\ \vdots \\ \mathrm{tr}(Y|T_{N-1}) \end{bmatrix} \qquad (26) $$

with $kl$-coefficients of $T_j$ such that

$$ T_0 = I, \\ T_j^{(k,l)} = \begin{cases} 2 & \text{if } k-l=j \\ 0 & \text{otherwise} \end{cases}, \quad j=1, \dots, N-1. \qquad (27) $$

For an admissible control $u$ valued in $\hat{K}_{\alpha}^{g}$, one has

$$ \int_{0}^{2\pi} \sum_{i=1}^{3} u_i(f) G_i(I, f) \, df = \sum_{j=1}^{g} (L_j c_j + \bar{L}_j \bar{c}_j) \qquad (28) $$

with $L_j(I)$ in $C^{5 \times N}$ defined by

$$ L_j(I) = \frac{1}{2} \sum_{i=1}^{3} \int_{0}^{2\pi} V_{ij} G_i(I, f) \Phi^H(f) \, df, \qquad (29) $$

where $V_j = (V_{ij})_{i=1,...,3}$. We note that the components of $L_j(I)$ are Fourier coefficients of the function $\sum_{i=1}^{3} V_{ij} G_i(I, f)$. $L_j(I)$ are approximated using a Discrete Fourier Transform. Since vector fields $G_i$ are smooth, truncation of the series is justified by the fast decrease of the coefficients. Finally, for a control $u$ valued in $\hat{K}_{\alpha}^{g}$, coefficients $v_j$ are truncated Fourier series of order
$N-1$. As a result, for a given vector $d_I$, the SDP approximation is

$$ 
\begin{array}{ll}
\max_{c_j \in \mathbb{C}^N, Y_b, Y_I \in \mathbb{C}^{N \times N}} & (\delta I | d_I) \quad \text{subject to} \quad \delta I = \varepsilon \sum_{j=1}^{g} (L_j c_j + \bar{L}_j \bar{c}_j) \quad \text{parallel to } d_I \\
& Y_b \ge 0, \quad \Lambda^*(Y_b) = c_j, \quad j = 1, \dots, g, \\
& Y_b \ge 0, \quad \Lambda^*(Y_b) = (1, 0, \dots, 0) - \sum_{j=1}^{g} c_j.
\end{array}
\quad (30) $$

The Lagrange variable of the discretisation of the equality constraint that $\delta I$ is parallel to $d_I$ from the convex program is expected to be a fair approximation of the costate $p_{\delta I}$ of (22). More importantly, it is hoped that the bang-bang control structure associated with this $p_{\delta I}$ is indeed the same as for the solution of the problem defined on $\hat{K}_\alpha$.

## 3.2 | Multiple shooting, differential continuation and callback

Homotopy, also referred to as continuation, allows to solve a complex problem by connecting it continuously to a simpler problem. The idea is then to follow the path (assumed to be regular enough) of solutions from the simpler problem towards the targeted one. See, e.g., references$^{16,17}$ for applications in optimal control. In our case, a parameter $\lambda$ defined between 0 and 1 allows to connect the problem with control set the bounded convex cone $\hat{K}_\alpha$ at $\lambda = 0$, to the original problem with the non-convex drop-like control set $U$ at $\lambda = 1$. In order to be able to solve the problem for $\lambda = 0$, we rely on the solution of the convex program on $\hat{K}_\alpha^g$ to provide an admissible solution. This solution is used not only to compute an educated guess for the initial costate but also to devise the appropriate multiple shooting function. To do so, we use the control structure corresponding to the approximation of $p_{\delta I}$ provided by the convex optimisation and described at Proposition 3. This proposition tells us that, when $\psi$ (a function of $p_{\delta I}$ and $f$) belongs to the open complement of the polar cone $K_\alpha^0$, the control must be equal to the dynamical feedback described in case (ii-a) (apart for some isolated points that correspond to case (ii-b) that we can neglect); we denote $u_b^0(f, p_{\delta I})$ this control. Similarly, for such values of $\psi$, Proposition 2 for the problem on $\text{conv}(U)$—and actually $U$, check Corollary 1—, implies that the control must be a solution of (15)-(16). (While these equations provide an explicit solution for the coordinate $\delta$ of the control, $\beta$ is only implicitly defined and we discuss its actual computation in Section 3.3.) We assume that this solution is unique and denote it $u_b^1(f, p_{\delta I})$. Then, for $\lambda$ in $[0, 1]$ and $\psi$ outside the polar cone, we define

$$ u_b(f, p_{\delta I}, \lambda) := (1 - \lambda)u_b^0(f, p_{\delta I}) + \lambda u_b^1(f, p_{\delta I}) \quad (31) $$

as the convex combination of the dynamical feedbacks for $\lambda = 0$ and $\lambda = 1$. (Such convex combinations are a common way to craft continuous in the parameter connections between nonlinear functions.) Conversely, for any $\lambda$ in $[0, 1]$ and $\psi$ in the interior of the polar cone, the control is set to zero.

For a given $\lambda$, one has a finite sequence of arcs with either $u = u_b$ (bang arcs), or $u = 0$ (zero arcs). Contacts with $\partial K_\alpha^0$ are characterized by (17) whose left-hand side defines the switching function, denoted $\varphi(f, p_{\delta I})$ (not depending on $\lambda$ in our particular setting). To this finite sequence of arcs is associated a multiple shooting function in a standard fashion. Assume for instance that the structure is bang-zero-bang. Then the shooting function has three arguments: the (constant) value of the costate, $p_{\delta I}$, and the two switchings times (true anomalies) bounding the central zero arc, $f_1$ and $f_2$. (So that $(p_{\delta I}, f_1, f_2)$ belong to $\mathbb{R}^7$.) Plugging $u = u_b(f, p_{\delta I}, \lambda)$ into the dynamics of $\delta I$ and integrating on $[0, f_1]$ from $\delta I(0) = 0$ allows to compute $\delta I_1 := \delta I(f_1)$. As the control is zero on $[f_1, f_2]$, $\delta I$ remains constant on the corresponding arc and we set $\delta I_2 := \delta I_1$. The control $u = u_b(f, p_{\delta I}, \lambda)$ is eventually plugged again on $[f_2, 2\pi]$ to compute $\delta I_f := \delta I(2\pi)$, starting from $\delta I_2$. The associated value of the shooting function is obtained by concatenating the left-hand side of the four equations below, forming a vector of dimension $4+1+2=7$ (note that the first colinearity equation indeed has dimension $5-1=4$):

$$ \delta I_f \wedge d_I = 0, \quad (32) $$

$$ (p_{\delta I} | d_I) - 1 = 0, \quad (33) $$

$$ \varphi(f_1, p_{\delta I}) = 0, \quad (34) $$

$$ \varphi(f_2, p_{\delta I}) = 0. \quad (35) $$

This defines a shooting function $S(\xi, \lambda)$ with, for this bang-zero-bang structure, $\xi := (p_{\delta I}, f_1, f_2)$. Once the first solution for $\lambda = 0$ is obtained, the path of zeros is followed by differential continuation, typically using a parametrisation by its curvilinear
$N-1$. As a result, for a given vector $d_I$, the SDP approximation is

$$ \begin{array}{ll} \max_{c_j \in \mathbb{C}^N, Y_b, Y_I \in \mathbb{C}^{N \times N}} & (\delta I | d_I) \quad \text{subject to} \quad \delta I = \varepsilon \sum_{j=1}^g (L_j c_j + \bar{L}_j \bar{c}_j) \quad \text{parallel to } d_I \\ & Y_b \ge 0, \quad \Lambda^*(Y_b) = c_j, \quad j = 1, \dots, g, \\ & Y_b \ge 0, \quad \Lambda^*(Y_b) = (1, 0, \dots, 0) - \sum_{j=1}^g c_j. \end{array} \quad (30) $$

The Lagrange variable of the discretisation of the equality constraint that $\delta I$ is parallel to $d_I$ from the convex program is expected to be a fair approximation of the costate $p_{\delta I}$ of (22). More importantly, it is hoped that the bang-bang control structure associated with this $p_{\delta I}$ is indeed the same as for the solution of the problem defined on $\hat{K}_\alpha$.

## 3.2 | Multiple shooting, differential continuation and callback

Homotopy, also referred to as continuation, allows to solve a complex problem by connecting it continuously to a simpler problem. The idea is then to follow the path (assumed to be regular enough) of solutions from the simpler problem towards the targeted one. See, e.g., references$^{16,17}$ for applications in optimal control. In our case, a parameter $\lambda$ defined between 0 and 1 allows to connect the problem with control set the bounded convex cone $\hat{K}_\alpha$ at $\lambda = 0$, to the original problem with the non-convex drop-like control set $U$ at $\lambda = 1$. In order to be able to solve the problem for $\lambda = 0$, we rely on the solution of the convex program on $\hat{K}_\alpha^g$ to provide an admissible solution. This solution is used not only to compute an educated guess for the initial costate but also to devise the appropriate multiple shooting function. To do so, we use the control structure corresponding to the approximation of $p_{\delta I}$ provided by the convex optimisation and described at Proposition 3. This proposition tells us that, when $\psi$ (a function of $p_{\delta I}$ and $f$) belongs to the open complement of the polar cone $K_\alpha^0$, the control must be equal to the dynamical feedback described in case (ii-a) (apart for some isolated points that correspond to case (ii-b) that we can neglect); we denote $u_b^0(f, p_{\delta I})$ this control. Similarly, for such values of $\psi$, Proposition 2 for the problem on $\text{conv}(U)$—and actually $U$, check Corollary 1—, implies that the control must be a solution of (15)-(16). (While these equations provide an explicit solution for the coordinate $\delta$ of the control, $\beta$ is only implicitly defined and we discuss its actual computation in Section 3.3.) We assume that this solution is unique and denote it $u_b^1(f, p_{\delta I})$. Then, for $\lambda$ in $[0, 1]$ and $\psi$ outside the polar cone, we define

$$ u_b(f, p_{\delta I}, \lambda) := (1 - \lambda)u_b^0(f, p_{\delta I}) + \lambda u_b^1(f, p_{\delta I}) \quad (31) $$

as the convex combination of the dynamical feedbacks for $\lambda = 0$ and $\lambda = 1$. (Such convex combinations are a common way to craft continuous in the parameter connections between nonlinear functions.) Conversely, for any $\lambda$ in $[0, 1]$ and $\psi$ in the interior of the polar cone, the control is set to zero.

For a given $\lambda$, one has a finite sequence of arcs with either $u = u_b$ (bang arcs), or $u = 0$ (zero arcs). Contacts with $\partial K_\alpha^0$ are characterized by (17) whose left-hand side defines the switching function, denoted $\varphi(f, p_{\delta I})$ (not depending on $\lambda$ in our particular setting). To this finite sequence of arcs is associated a multiple shooting function in a standard fashion. Assume for instance that the structure is bang-zero-bang. Then the shooting function has three arguments: the (constant) value of the costate, $p_{\delta I}$, and the two switchings times (true anomalies) bounding the central zero arc, $f_1$ and $f_2$. (So that $(p_{\delta I}, f_1, f_2)$ belong to $\mathbb{R}^7$.) Plugging $u = u_b(f, p_{\delta I}, \lambda)$ into the dynamics of $\delta I$ and integrating on $[0, f_1]$ from $\delta I(0) = 0$ allows to compute $\delta I_1 := \delta I(f_1)$. As the control is zero on $[f_1, f_2]$, $\delta I$ remains constant on the corresponding arc and we set $\delta I_2 := \delta I_1$. The control $u = u_b(f, p_{\delta I}, \lambda)$ is eventually plugged again on $[f_2, 2\pi]$ to compute $\delta I_f := \delta I(2\pi)$, starting from $\delta I_2$. The associated value of the shooting function is obtained by concatenating the left-hand side of the four equations below, forming a vector of dimension $4+1+2=7$ (note that the first colinearity equation indeed has dimension $5-1=4$):

$$ \delta I_f \wedge d_I = 0, \quad (32) $$

$$ (p_{\delta I} | d_I) - 1 = 0, \quad (33) $$

$$ \varphi(f_1, p_{\delta I}) = 0, \quad (34) $$

$$ \varphi(f_2, p_{\delta I}) = 0. \quad (35) $$

This defines a shooting function $S(\xi, \lambda)$ with, for this bang-zero-bang structure, $\xi := (p_{\delta I}, f_1, f_2)$. Once the first solution for $\lambda = 0$ is obtained, the path of zeros is followed by differential continuation, typically using a parametrisation by its curvilinear
abscissa:

$$s \mapsto (\lambda(s), \xi(s)) \text{ with } S(\xi(s), \lambda(s)) = 0.$$

We refer, e.g., to¹⁸ for the assumptions needed to do so. Note that, according to (33), we look for normal extremals (compare with (14)). Although abnormal extremals with $p^0 = 0$ in (14) are not expected to fill a large (in the sense of measure) subset of the state space, it would be possible to encompass them in the search by extending our set of variables to $\xi := (p^0, p_{\delta I}, f_1, f_2)$ (scalar $p^0$ added), and replacing (33) by two new equations:

$$ (p_{\delta I} | d_I) + p^0 = 0, \\ \|p_{\delta I}\|^2 + (p^0)^2 = 1. $$

The second equation is used to take into account homogeneity in $(p^0, p)$, and one has to check a *posteriori* that $p^0$ is nonpositive.

One important issue in practice is that it might not be possible to reach $\lambda = 1$ because, at some $\lambda(\bar{s})$ in (0, 1), the structure of the solution changes; for instance because one subarcs disappears. It is crucial to be able to detect such a change during homotopy since then, the shooting function has to be redefined according to the new structure. This is achieved using a standard callback mechanism along with differential continuation. On the previous bang-zero-bang example, the continuation is monitored and, at each step of the path following procedure, a simple test is performed: if the exit time of the zero arc, $f_2$, becomes inferior to the entry time $f_1$ (this is detected by a sign change on $f_2 - f_1$, as going forward in time makes sense mathematically but is not allowed to obtain admissible trajectories), the continuation is stopped. And restarted at $\lambda(\bar{s})$ with a new shooting function (in this case, a single shooting one, as only one bang arc would be left), using $\xi(\bar{s})$ as initial guess. More elaborated tests can be constructed to detect a new arc appearing, etc. In our case, a callback is used to detect a structure change from 5 subarcs to 3 (see Section 4). This methodology can readily be extended to any other structure with a finite number of arcs, and even to a more complicated structure change (e.g. bang-bang towards bang-singular-bang) provided one is able to write the appropriate test in the callback.

## 3.3 | Implicit treatment of the Hamiltonian maximisation

Regarding the computation of $u_b^1(f, p_{\delta I})$, we know after Proposition 2 that the control is either zero, either solution of (15-16). The first equation for the coordinate $\beta$ of $u$ has no closed form solution. There is a preliminary numerical discussion of the number of solutions in² (we actually look for a global maximiser of the Hamiltonian over $U$, which may allow to eliminate some strictly local minimiser that also verify (15)) for a particular set of values of the sail parameters. More generally, while maximisation of the Hamiltonian often yields an explicit expression of the control as a dynamics feedback function of the state and the costate, it is not always the case. In such a situation, we advocate an implicit treatment of this maximisation, incorporating the stationarity equation of the Hamiltonian into the shooting procedure. We sketch below a simple way to do so in a general setting. Note that this approach also makes sense for other numerical methods to solve the boundary value problem resulting from the maximum principle (collocation, e.g.); it allows not to care about explicit Hamiltonian maximisation, for instance in convoluted examples coming from penalty methods.¹⁹

Assume that, after applying Pontrjagin maximum principle, one has to integrate the following system (x denoting the state, p the costate):

$$ \dot{x}(t) = \nabla_p H(x(t), p(t), u(t)), \quad \dot{p}(t) = -\nabla_x H(x(t), p(t), u(t)), \qquad (36) $$

where, at each time t, the m-dimensional control $u(t)$ verifies

$$ \nabla_u H(x(t), p(t), u(t)) = 0. \qquad (37) $$

The last stationarity equation corresponds to an unconstrained situation—whereas a Lagrangian, plus an additional finite dimensional multiplier, should be considered in the presence of constraints—, and defines a semi-explicit differential-algebraic-equation. Assume that the strong Legendre-Clebsch condition holds in an open neighbourhood of the reference extremal times the open control set, $\nabla_{uu}^2 H \le -cI_m$ for some positive constant $c$. Then the Hamiltonian has a unique maximiser, that satisfies $\nabla_u H = 0$, and the previous differential-algebraic equation (DAE) is of index 1 (differentiating once (37) allows to solve for $\dot{u}$). In particular, one can extend the Hamiltonian system (36) by adding the equation

$$ \dot{u} = -\nabla_{uu}^2 H^{-1}(\nabla_{ux} H \cdot \nabla_p H - \nabla_{up} H \cdot \nabla_x H)(x, p, u) := g(x, p, u), \qquad (38) $$
abscissa:

$$s \mapsto (\lambda(s), \xi(s)) \text{ with } S(\xi(s), \lambda(s)) = 0.$$

We refer, e.g., to¹⁸ for the assumptions needed to do so. Note that, according to (33), we look for normal extremals (compare with (14)). Although abnormal extremals with $p^0 = 0$ in (14) are not expected to fill a large (in the sense of measure) subset of the state space, it would be possible to encompass them in the search by extending our set of variables to $\xi := (p^0, p_{\delta I}, f_1, f_2)$ (scalar $p^0$ added), and replacing (33) by two new equations:

$$\begin{gathered} (p_{\delta I} | d_I) + p^0 = 0, \\ \|p_{\delta I}\|^2 + (p^0)^2 = 1. \end{gathered}$$

The second equation is used to take into account homogeneity in $(p^0, p)$, and one has to check a *posteriori* that $p^0$ is nonpositive.

One important issue in practice is that it might not be possible to reach $\lambda = 1$ because, at some $\lambda(\bar{s})$ in (0, 1), the structure of the solution changes; for instance because one subarcs disappears. It is crucial to be able to detect such a change during homotopy since then, the shooting function has to be redefined according to the new structure. This is achieved using a standard callback mechanism along with differential continuation. On the previous bang-zero-bang example, the continuation is monitored and, at each step of the path following procedure, a simple test is performed: if the exit time of the zero arc, $f_2$, becomes inferior to the entry time $f_1$ (this is detected by a sign change on $f_2 - f_1$, as going forward in time makes sense mathematically but is not allowed to obtain admissible trajectories), the continuation is stopped. And restarted at $\lambda(\bar{s})$ with a new shooting function (in this case, a single shooting one, as only one bang arc would be left), using $\xi(\bar{s})$ as initial guess. More elaborated tests can be constructed to detect a new arc appearing, etc. In our case, a callback is used to detect a structure change from 5 subarcs to 3 (see Section 4). This methodology can readily be extended to any other structure with a finite number of arcs, and even to a more complicated structure change (e.g. bang-bang towards bang-singular-bang) provided one is able to write the appropriate test in the callback.

## 3.3 | Implicit treatment of the Hamiltonian maximisation

Regarding the computation of $u_b^1(f, p_{\delta I})$, we know after Proposition 2 that the control is either zero, either solution of (15-16). The first equation for the coordinate $\beta$ of $u$ has no closed form solution. There is a preliminary numerical discussion of the number of solutions in² (we actually look for a global maximiser of the Hamiltonian over $U$, which may allow to eliminate some strictly local minimiser that also verify (15)) for a particular set of values of the sail parameters. More generally, while maximisation of the Hamiltonian often yields an explicit expression of the control as a dynamics feedback function of the state and the costate, it is not always the case. In such a situation, we advocate an implicit treatment of this maximisation, incorporating the stationarity equation of the Hamiltonian into the shooting procedure. We sketch below a simple way to do so in a general setting. Note that this approach also makes sense for other numerical methods to solve the boundary value problem resulting from the maximum principle (collocation, e.g.); it allows not to care about explicit Hamiltonian maximisation, for instance in convoluted examples coming from penalty methods.¹⁹

Assume that, after applying Pontrjagin maximum principle, one has to integrate the following system (x denoting the state, p the costate):

$$\dot{x}(t) = \nabla_p H(x(t), p(t), u(t)), \quad \dot{p}(t) = -\nabla_x H(x(t), p(t), u(t)), \qquad (36)$$

where, at each time t, the m-dimensional control $u(t)$ verifies

$$\nabla_u H(x(t), p(t), u(t)) = 0. \qquad (37)$$

The last stationarity equation corresponds to an unconstrained situation—whereas a Lagrangian, plus an additional finite dimensional multiplier, should be considered in the presence of constraints—, and defines a semi-explicit differential-algebraic-equation. Assume that the strong Legendre-Clebsch condition holds in an open neighbourhood of the reference extremal times the open control set, $\nabla_{uu}^2 H \le -c I_m$ for some positive constant $c$. Then the Hamiltonian has a unique maximiser, that satisfies $\nabla_u H = 0$, and the previous differential-algebraic equation (DAE) is of index 1 (differentiating once (37) allows to solve for $\dot{u}$). In particular, one can extend the Hamiltonian system (36) by adding the equation

$$\dot{u} = -\nabla_{uu}^2 H^{-1}(\nabla_{ux} H \cdot \nabla_p H - \nabla_{up} H \cdot \nabla_x H)(x, p, u) := g(x, p, u), \qquad (38)$$
<figure><img src="image_9.png" /><figcaption>Figure 8 Algorithm for solving optimal control problem (OCP).</figcaption></figure>

with initial condition $\nabla_u H(x(0), p(0), u(0)) = 0$. The new system remains Hamiltonian as is clear setting $\hat{x} := (x, u)$, $\hat{p} := (p, p_u)$ and

$$ \hat{H}(x, u, p, p_u) := H(x, p, u) + (p_u|g(x, p, u)) $$

with $p_u(0) = 0$. (One can obviously eliminate the trivial equation on $p_u$, which is an extra but identically zero costate.) In the case of a shooting approach, the value of $u(0)$ is an additional shooting variable. Keeping the system in Hamiltonian form is convenient in the algorithmic framework described in Section 4, but other approaches for DAE such as predictor-corrector ones can of course be considered. In our case, we use this approach with $x = \delta I$, $p = p_{\delta I}$ to deal with the implicit equation (15) on $\beta$ (while we use (16) to solve it explicitly for $\delta$). The combination of this implicit approach with multiple shooting, homotopy and callback is described in the last section.

## 4 | NUMERICAL EXAMPLE: THE NASA SQUARE SAIL

This section presents the results for a specific example of sail orbiting around a planet for which a change in the solution structure occurs during differential continuation. It is important to note that, in our simulations, such a change is more of an exception rather than the norm as the majority of the simulated trajectories exhibited the same structure both for the bounded cone at $\lambda = 0$ and the original drop-like control set at $\lambda = 1$. This particular example highlights the need to use a callback procedure to detect the change of structure and restart the continuation with an updated number of arcs that define the multiple shooting function, following the algorithm highlighted in Figure 8. We consider the Square Sail from NASA Jet Propulsion Lab defined in ¹, Table 2.¹ whose optical properties are as follows: $\rho = 0.88$, $s = 0.94$, $\epsilon_b = 0.55$, $\epsilon_f = 0.05$, $B_b = 0.55$, $B_f = 0.79$. The orbit is determined by $I = (10^\circ, 50^\circ, 30^\circ, 1, 0.1)$. The desired manoeuvre to be performed is an increase of the inclination $\gamma_2$, so that
<figure><img src="image_9.png" /><figcaption>Figure 8 Algorithm for solving optimal control problem (OCP).</figcaption></figure>

with initial condition $\nabla_u H(x(0), p(0), u(0)) = 0$. The new system remains Hamiltonian as is clear setting $\hat{x} := (x, u), \hat{p} := (p, p_u)$ and

$$ \hat{H}(x, u, p, p_u) := H(x, p, u) + (p_u|g(x, p, u)) $$

with $p_u(0) = 0$. (One can obviously eliminate the trivial equation on $p_u$, which is an extra but identically zero costate.) In the case of a shooting approach, the value of $u(0)$ is an additional shooting variable. Keeping the system in Hamiltonian form is convenient in the algorithmic framework described in Section 4, but other approaches for DAE such as predictor-corrector ones can of course be considered. In our case, we use this approach with $x = \delta I, p = p_{\delta I}$ to deal with the implicit equation (15) on $\beta$ (while we use (16) to solve it explicitly for $\delta$). The combination of this implicit approach with multiple shooting, homotopy and callback is described in the last section.

## 4 | NUMERICAL EXAMPLE: THE NASA SQUARE SAIL

This section presents the results for a specific example of sail orbiting around a planet for which a change in the solution structure occurs during differential continuation. It is important to note that, in our simulations, such a change is more of an exception rather than the norm as the majority of the simulated trajectories exhibited the same structure both for the bounded cone at $\lambda = 0$ and the original drop-like control set at $\lambda = 1$. This particular example highlights the need to use a callback procedure to detect the change of structure and restart the continuation with an updated number of arcs that define the multiple shooting function, following the algorithm highlighted in Figure 8. We consider the Square Sail from NASA Jet Propulsion Lab defined in ¹, Table 2.1 whose optical properties are as follows: $\rho = 0.88, s = 0.94, \epsilon_b = 0.55, \epsilon_f = 0.05, B_b = 0.55, B_f = 0.79$. The orbit is determined by $I = (10^\circ, 50^\circ, 30^\circ, 1, 0.1)$. The desired manoeuvre to be performed is an increase of the inclination $\gamma_2$, so that
<figure><img src="image_10.png" /><figcaption>**Figure 9** Results from convex optimisation (in red) on $\hat{K}_\alpha^g$ compared to multiple shooting on $\hat{K}_\alpha$ with $\lambda = 0$ (in blue). Only the norm of the controls are represented. The solution of the problem on the polyhedral approximation of the cone (here with $g = 18$) does capture the 5-arc structure of the control, and provides a very accurate initial guess for the multiple shooting on the bounded cone. Small Gibbs phenomenon is observed on the red curve for a truncation at order $N = 80$ of the Fourier series.</figcaption></figure>

gives $d_I = (0, 1, 0, 0, 0)$. We first solve a convex optimisation problem on the bounded polyhedral cone, $\hat{K}_\alpha^g$. This result serves as a reliable initial guess for the optimal control problem considered. This initial guess consists of five arcs, with the first, third and fifth arc begin zero or coast arcs. The costate approximation retrieved from the Lagrange multiplier associated with the target direction $d_I$ of the convex program (see Section 3.1) is

$$p_{\delta I} = (-0.0837, 1, -0.0052, 0.0398, 0.0852),$$

and switchings between zeros and bangs occurs at $f_1 = 49.4^\circ$, $f_2 = 237.9^\circ$, $f_3 = 265.6^\circ$ and $f_4 = 286.9^\circ$. These values allow to find a zero $\xi := (p_{\delta I}, f_1, f_2, f_3, f_4)$ of the multiple shooting function on the bounded cone $\hat{K}_\alpha$ for $\lambda = 0$. Both controls, on $\hat{K}_\alpha^g$ and $\hat{K}_\alpha$, are portrayed Figure 9. The remarkable precision of the solution provided by convex optimisation is observed: not only does it give the good control structure (here with five arcs), but also the control values on each subarc and the switching times are very accurately initialised.

During the continuation, initiated with a $\gamma_0\gamma_+\gamma_0\gamma_+$ 5-arc structure (a $\gamma_0$ denotes an arc with zero control, corresponding to the alternative (i) of Proposition 3—see also discussion after (31) for $\lambda$ in [0, 1]—, while $\gamma_+$ corresponds to alternative (ii-a) of the same proposition), the callback function detects that $f_4 < f_3$, which means that the fourth $\gamma_+$ arc vanishes. On this example, the change of structure takes place around $\lambda \approx 0.0256$. Continuation is stopped there, then restarted with a multiple shooting function devised for a $\gamma_0\gamma_+\gamma_0$ structure whose argument is now $\xi := (p_{\delta I}, f_1, f_2)$ (a pair of switching times have been removed). This structure is preserved until $\lambda = 1$, allowing to solve the targeted problem (check Figure 10). The solution adjoint state is

$$p_{\delta I} = (-0.1637, 1, -0.0972, 0.0712, 1.6037).$$

Figure 11 shows the solution on the real control set in terms of control angles, defining attitude of the solar sail. When $\beta = \frac{\pi}{2}$, the sail is aligned with the solar rays, so no part at all of surface is exposed to the solar light, which corresponds to $\|u\| = 0$. (In this case the angle $\delta$ is not defined.)

Convex optimisation is performed using CVX.20,21 Multiple shooting with differential continuation is done thanks to the control-toolbox package, leveraging the Hamiltonian setting discussed in the previous section to compute the shooting function as a concatenation of Hamiltonian flows. More precisely, we use two Hamiltonians to compute extremals of the system with extended state $(\delta I, \beta)$ (the first angle of the control has been added and is computed according to the implicit procedure described Section 3.3) and costate $(p_{\delta I}, p_\beta)$: first, the trivial Hamiltonian associated with the zero control (coast arc), second the Hamiltonian that corresponds to Proposition 2. The resulting flows are then automatically computed using automatic differentiation and numerical integration, and concatenated to evaluate the multiple shooting function. For instance, for a $\gamma_0\gamma_+\gamma_0$ structure, three Hamiltonian flows are evaluated. Note that in order to compute the right-hand side for $p_\beta$ (extra adjoint associated with the
<figure><img src="image_10.png" /><figcaption>**Figure 9** Results from convex optimisation (in red) on $\hat{K}_\alpha^g$ compared to multiple shooting on $\hat{K}_\alpha$ with $\lambda = 0$ (in blue). Only the norm of the controls are represented. The solution of the problem on the polyhedral approximation of the cone (here with $g = 18$) does capture the 5-arc structure of the control, and provides a very accurate initial guess for the multiple shooting on the bounded cone. Small Gibbs phenomenon is observed on the red curve for a truncation at order $N = 80$ of the Fourier series.</figcaption></figure>

gives $d_I = (0, 1, 0, 0, 0)$. We first solve a convex optimisation problem on the bounded polyhedral cone, $\hat{K}_\alpha^g$. This result serves as a reliable initial guess for the optimal control problem considered. This initial guess consists of five arcs, with the first, third and fifth arc begin zero or coast arcs. The costate approximation retrieved from the Lagrange multiplier associated with the target direction $d_I$ of the convex program (see Section 3.1) is

$$p_{\delta I} = (-0.0837, 1, -0.0052, 0.0398, 0.0852),$$

and switchings between zeros and bangs occurs at $f_1 = 49.4^\circ$, $f_2 = 237.9^\circ$, $f_3 = 265.6^\circ$ and $f_4 = 286.9^\circ$. These values allow to find a zero $\xi := (p_{\delta I}, f_1, f_2, f_3, f_4)$ of the multiple shooting function on the bounded cone $\hat{K}_\alpha$ for $\lambda = 0$. Both controls, on $\hat{K}_\alpha^g$ and $\hat{K}_\alpha$, are portrayed Figure 9. The remarkable precision of the solution provided by convex optimisation is observed: not only does it give the good control structure (here with five arcs), but also the control values on each subarc and the switching times are very accurately initialised.

During the continuation, initiated with a $\gamma_0\gamma_+\gamma_0\gamma_+\gamma_0$ 5-arc structure (a $\gamma_0$ denotes an arc with zero control, corresponding to the alternative (i) of Proposition 3—see also discussion after (31) for $\lambda$ in [0, 1]—, while $\gamma_+$ corresponds to alternative (ii-a) of the same proposition), the callback function detects that $f_4 < f_3$, which means that the fourth $\gamma_+$ arc vanishes. On this example, the change of structure takes place around $\lambda \approx 0.0256$. Continuation is stopped there, then restarted with a multiple shooting function devised for a $\gamma_0\gamma_+\gamma_0$ structure whose argument is now $\xi := (p_{\delta I}, f_1, f_2)$ (a pair of switching times have been removed). This structure is preserved until $\lambda = 1$, allowing to solve the targeted problem (check Figure 10). The solution adjoint state is

$$p_{\delta I} = (-0.1637, 1, -0.0972, 0.0712, 1.6037).$$

Figure 11 shows the solution on the real control set in terms of control angles, defining attitude of the solar sail. When $\beta = \frac{\pi}{2}$, the sail is aligned with the solar rays, so no part at all of surface is exposed to the solar light, which corresponds to $\|u\| = 0$. (In this case the angle $\delta$ is not defined.)

Convex optimisation is performed using CVX.20,21 Multiple shooting with differential continuation is done thanks to the control-toolbox package, leveraging the Hamiltonian setting discussed in the previous section to compute the shooting function as a concatenation of Hamiltonian flows. More precisely, we use two Hamiltonians to compute extremals of the system with extended state $(\delta I, \beta)$ (the first angle of the control has been added and is computed according to the implicit procedure described Section 3.3) and costate $(p_{\delta I}, p_\beta)$: first, the trivial Hamiltonian associated with the zero control (coast arc), second the Hamiltonian that corresponds to Proposition 2. The resulting flows are then automatically computed using automatic differentiation and numerical integration, and concatenated to evaluate the multiple shooting function. For instance, for a $\gamma_0\gamma_+\gamma_0$ structure, three Hamiltonian flows are evaluated. Note that in order to compute the right-hand side for $p_\beta$ (extra adjoint associated with the
<figure><img src="image_11.png" /><figcaption>**Figure 10** Differential continuation and callback, control norms. Left: homotopy from the 5-arc $\gamma_0\gamma_+\gamma_0\gamma_+\gamma_0$ solution from $\lambda = 0$ (black) to $\lambda \approx 0.0256$ (red), just before a change of structure (vanishing of the fourth $\gamma_+$ arc) is detected by the callback. Intermediate solutions are represented in grey. Right: homotopy from $\lambda \approx 0.0256$ (red) to $\lambda = 1$ (black). The continuation is restarted at $\lambda \approx 0.0256$ with a 3-arc $\gamma_0\gamma_+\gamma_0$ structure for the multiple shooting. As $\lambda$ increases towards 1, no further change of structure is detected. The convex combination (31) of controls becomes essentially valued in $U$ and gets larger. Intermediate solutions are again represented in grey.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11** Control for the target problem ($\lambda = 1$), $(\beta, \delta)$ coordinates. The 3-arc $\gamma_0\gamma_+\gamma_0$ is clearly observed with $\beta = 90^\circ$ on the initial and final zero arcs (and $\delta$ undefined there).</figcaption></figure>

implicit computation of the angle $\beta$ of the control), equation (38) must be differentiated once (actually twice since the linearisation of the ODE, or *variational system*, is also assembled so as to compute efficiently the derivative of the shooting function), which involves high order derivatives of the original Hamiltonian (up to order four, in practice). On nontrivial examples (here Gauss equations (10) are used), this is only possible by means of automatic differentiation. A predictor-corrector method is used to follow the path of zeros with, in practice, strong requirements on the prediction (ODE integration) and few correction steps. Regarding the implicit computation of the maximising control discussed in Section 3.3, we monitor the computation of the Hessian of the Hamiltonian wrt. $\beta$ on bang arcs along the continuation to ensure that at least of strict local solution is obtained. As we start from an initial guess provided by convex optimisation, we actually expect that a global maximiser is preserved along
<figure><img src="image_11.png" /><figcaption>**Figure 10** Differential continuation and callback, control norms. Left: homotopy from the 5-arc $\gamma_0\gamma_+\gamma_0\gamma_+\gamma_0$ solution from $\lambda = 0$ (black) to $\lambda \approx 0.0256$ (red), just before a change of structure (vanishing of the fourth $\gamma_+$ arc) is detected by the callback. Intermediate solutions are represented in grey. Right: homotopy from $\lambda \approx 0.0256$ (red) to $\lambda = 1$ (black). The continuation is restarted at $\lambda \approx 0.0256$ with a 3-arc $\gamma_0\gamma_+\gamma_0$ structure for the multiple shooting. As $\lambda$ increases towards 1, no further change of structure is detected. The convex combination (31) of controls becomes essentially valued in $U$ and gets larger. Intermediate solutions are again represented in grey.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11** Control for the target problem ($\lambda = 1$), $(\beta, \delta)$ coordinates. The 3-arc $\gamma_0\gamma_+\gamma_0$ is clearly observed with $\beta = 90^\circ$ on the initial and final zero arcs (and $\delta$ undefined there).</figcaption></figure>

implicit computation of the angle $\beta$ of the control), equation (38) must be differentiated once (actually twice since the linearisation of the ODE, or *variational system*, is also assembled so as to compute efficiently the derivative of the shooting function), which involves high order derivatives of the original Hamiltonian (up to order four, in practice). On nontrivial examples (here Gauss equations (10) are used), this is only possible by means of automatic differentiation. A predictor-corrector method is used to follow the path of zeros with, in practice, strong requirements on the prediction (ODE integration) and few correction steps. Regarding the implicit computation of the maximising control discussed in Section 3.3, we monitor the computation of the Hessian of the Hamiltonian wrt. $\beta$ on bang arcs along the continuation to ensure that at least of strict local solution is obtained. As we start from an initial guess provided by convex optimisation, we actually expect that a global maximiser is preserved along
<figure><img src="image_13.png" /></figure>

**Figure 12** Integration of the dynamics of the initial system using the optimal control solution. As anticipated, all orbital parameters remain constant throughout the orbital period, with the exception of $γ_2$, which continues to increase.

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

**Figure 13** Resulting trajectory of the sail using the optimal control solution. Orbital inclination increases while other parameters stay constant. Red arrows show direction of the SRP force generated by the sail.

the entire continuation, but this would require further analysis. The code reproducing the presented results is available and executable online. <fnref n="3" /> To verify the accuracy of the solution, we integrate the trajectory of the initial system using the provided solution. Figure 12 depicts the dynamics of each of the five orbital elements over one orbit. The magnitude of the displacement, $\delta I$, depends on the parameter $\epsilon$ of the sail. However, since it does not impact the interpretation of the results, we normalize the displacement. As shown, all orbital parameters remain constant after the orbital period, except for $γ_2$, which continues to increase. It is important to note that our algorithm assumes a fixed state for one orbital period, but the final integration is performed using the actual dynamics. As expected, divergence becomes noticeable after integrating for ten or more orbits, which is consistent with the fact that the optimal control is provided for the initial orbit. Ideally, the proposed methodology should be utilized as a feedback control system to determine the optimal trajectory for each orbital period, requiring reinitialization after each orbit. Finally, Figure 13 shows the trajectory of the sail around a body after one orbital period. The inclination of the orbit changes indeed, according to the required direction of the displacement.

<footnote n="3">³Check ct.gitlabpages.inria.fr/gallery</footnote>
<figure><img src="image_13.png" /></figure>

**Figure 12** Integration of the dynamics of the initial system using the optimal control solution. As anticipated, all orbital parameters remain constant throughout the orbital period, with the exception of $γ_2$, which continues to increase.

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

**Figure 13** Resulting trajectory of the sail using the optimal control solution. Orbital inclination increases while other parameters stay constant. Red arrows show direction of the SRP force generated by the sail.

the entire continuation, but this would require further analysis. The code reproducing the presented results is available and executable online.<fnref n="3" /> To verify the accuracy of the solution, we integrate the trajectory of the initial system using the provided solution. Figure 12 depicts the dynamics of each of the five orbital elements over one orbit. The magnitude of the displacement, $\delta I$, depends on the parameter $\epsilon$ of the sail. However, since it does not impact the interpretation of the results, we normalize the displacement. As shown, all orbital parameters remain constant after the orbital period, except for $γ_2$, which continues to increase. It is important to note that our algorithm assumes a fixed state for one orbital period, but the final integration is performed using the actual dynamics. As expected, divergence becomes noticeable after integrating for ten or more orbits, which is consistent with the fact that the optimal control is provided for the initial orbit. Ideally, the proposed methodology should be utilized as a feedback control system to determine the optimal trajectory for each orbital period, requiring reinitialization after each orbit. Finally, Figure 13 shows the trajectory of the sail around a body after one orbital period. The inclination of the orbit changes indeed, according to the required direction of the displacement.

<footnote n="3">³Check ct.gitlabpages.inria.fr/gallery</footnote>
# ACKNOWLEDGEMENTS

This work was partially supported by ESA (contract no. 4000134950/21/NL/GLC/my). The authors wish to express their gratitude to S. Bennani from ESA (ESTEC, Noordwijk) for many fruitful exchanges during the preparation of this work. The authors also thank T. Calendini from Polytech Nice Sophia, O. Cots from Université de Toulouse, and S. Sed from Inria Centre of Université Côte d'Azur for their help during the preparation of the numerical simulations of Section 4.

# References

1. McInnes C. *Solar Sailing*. Springer. 1999. doi: 10.1007/978-1-4471-3992-8

2. Mengali G, Quarta AA. Optimal Three-Dimensional Interplanetary Rendezvous Using Non-Ideal Solar Sail. *Journal of Guidance, Control, and Dynamics* 2005; 28(1): 173-177. doi: 10.2514/1.8325

3. Macdonald M, McInnes C, Dachwald B. Heliocentric Solar Sail Orbit Transfers with Locally Optimal Control Laws. *J Spacecraft Rockets* 2007; 44: 273-276. doi: 10.2514/1.17297

4. Montenbruck O, Gill E. *Satellite Orbits*. Springer. 2000. doi: 10.1007/978-3-642-58351-3

5. Herasimenka A, Dell'Elce L, Caillau JB, Pomet JB. Controllability Properties of Solar Sails. *Journal of Guidance, Control, and Dynamics* 2023; 45(9): 900-909. doi: 10.2514/1.g007250

6. Rios-Reyes L, Scheeres DJ. Generalized Model for Solar Sails. *J Spacecraft Rockets* 2005; 42(1): 182-185. doi: 10.2514/1.9054

7. Dachwald B, Macdonald M, McInnes CR, Mengali G, Quarta AA. Impact of Optical Degradation on Solar Sail Mission Performance. *J Spacecraft Rockets* 2007; 44(4): 740-749. doi: 10.2514/1.21432

8. Niccolai L, Quarta AA, Mengali G. Trajectory Approximation of a Solar Sail With Constant Pitch Angle and Optical Degradation. *IEEE Transactions on Aerospace and Electronic Systems* 2022; 58(4): 3643-3649. doi: 10.1109/taes.2021.3124867

9. Herasimenka A, Caillau JB, Dell'Elce L, Pomet JB. Controllability Test for Systems with Constrained Control. Application to Solar Sailing. *Proc. 20th European Control Conf.* 2022: 2143-2148. doi: 10.23919/ECC55457.2022.9838146

10. Caillau JB, Dell'Elce L, Herasimenka A, Pomet JB. On the Controllability of Nonlinear Systems with a Periodic Drift. HAL preprint hal-03779482, https://inria.hal.science/hal-03779482; 2022.

11. Caillau JB, Daoud B. Minimum time control of the restricted three-body problem. *SIAM J. Control Optim.* 2012; 50(6): 3178-3202.

12. Caillau JB, Féjoz J, Orieux M, Roussarie R. On singularities of min time affine control systems. *SIAM J. Control Optim.* 2022; 60(2): 1143-1162.

13. Boyd JP. Computing the zeros, maxima and inflection points of Chebyshev, Legendre and Fourier series: solving transcendental equations by spectral interpolation and polynomial rootfinding. *J Eng Math* 2007; 56(3): 203-219. doi: 10.1007/s10665-006-9087-5

14. Nesterov Y. Squared Functional Systems and Optimization Problems. In: Pardalos PM, Hearn D, Frenk H, Roos K, Terlaky T, Zhang S., eds. *High Performance Optimization*. 33. Boston, MA: Springer. 2000 (pp. 405-440)

15. Dumitrescu B. *Positive Trigonometric Polynomials and Signal Processing Applications*. Springer. 2017. doi: 10.1007/978-3-319-53688-0.

16. Gergaud J, Haberkorn T. Homotopy method for minimum consumption orbit transfer problem. *ESAIM Control Optim. and Calc. Var.* 2006; 12(2): 294-310.
# ACKNOWLEDGEMENTS

This work was partially supported by ESA (contract no. 4000134950/21/NL/GLC/my). The authors wish to express their gratitude to S. Bennani from ESA (ESTEC, Noordwijk) for many fruitful exchanges during the preparation of this work. The authors also thank T. Calendini from Polytech Nice Sophia, O. Cots from Université de Toulouse, and S. Sed from Inria Centre of Université Côte d'Azur for their help during the preparation of the numerical simulations of Section 4.

# References

1. McInnes C. *Solar Sailing*. Springer. 1999. doi: 10.1007/978-1-4471-3992-8

2. Mengali G, Quarta AA. Optimal Three-Dimensional Interplanetary Rendezvous Using Non-Ideal Solar Sail. *Journal of Guidance, Control, and Dynamics* 2005; 28(1): 173-177. doi: 10.2514/1.8325

3. Macdonald M, McInnes C, Dachwald B. Heliocentric Solar Sail Orbit Transfers with Locally Optimal Control Laws. *J Spacecraft Rockets* 2007; 44: 273-276. doi: 10.2514/1.17297

4. Montenbruck O, Gill E. *Satellite Orbits*. Springer. 2000. doi: 10.1007/978-3-642-58351-3

5. Herasimenka A, Dell'Elce L, Caillau JB, Pomet JB. Controllability Properties of Solar Sails. *Journal of Guidance, Control, and Dynamics* 2023; 45(9): 900-909. doi: 10.2514/1.g007250

6. Rios-Reyes L, Scheeres DJ. Generalized Model for Solar Sails. *J Spacecraft Rockets* 2005; 42(1): 182-185. doi: 10.2514/1.9054

7. Dachwald B, Macdonald M, McInnes CR, Mengali G, Quarta AA. Impact of Optical Degradation on Solar Sail Mission Performance. *J Spacecraft Rockets* 2007; 44(4): 740-749. doi: 10.2514/1.21432

8. Niccolai L, Quarta AA, Mengali G. Trajectory Approximation of a Solar Sail With Constant Pitch Angle and Optical Degradation. *IEEE Transactions on Aerospace and Electronic Systems* 2022; 58(4): 3643-3649. doi: 10.1109/taes.2021.3124867

9. Herasimenka A, Caillau JB, Dell'Elce L, Pomet JB. Controllability Test for Systems with Constrained Control. Application to Solar Sailing. *Proc. 20th European Control Conf.* 2022: 2143-2148. doi: 10.23919/ECC55457.2022.9838146

10. Caillau JB, Dell'Elce L, Herasimenka A, Pomet JB. On the Controllability of Nonlinear Systems with a Periodic Drift. HAL preprint hal-03779482, https://inria.hal.science/hal-03779482; 2022.

11. Caillau JB, Daoud B. Minimum time control of the restricted three-body problem. *SIAM J. Control Optim.* 2012; 50(6): 3178-3202.

12. Caillau JB, Féjoz J, Orieux M, Roussarie R. On singularities of min time affine control systems. *SIAM J. Control Optim.* 2022; 60(2): 1143-1162.

13. Boyd JP. Computing the zeros, maxima and inflection points of Chebyshev, Legendre and Fourier series: solving transcendental equations by spectral interpolation and polynomial rootfinding. *J Eng Math* 2007; 56(3): 203-219. doi: 10.1007/s10665-006-9087-5

14. Nesterov Y. Squared Functional Systems and Optimization Problems. In: Pardalos PM, Hearn D, Frenk H, Roos K, Terlaky T, Zhang S., eds. *High Performance Optimization*. 33. Boston, MA: Springer. 2000 (pp. 405-440)

15. Dumitrescu B. *Positive Trigonometric Polynomials and Signal Processing Applications*. Springer. 2017. doi: 10.1007/978-3-319-53688-0.

16. Gergaud J, Haberkorn T. Homotopy method for minimum consumption orbit transfer problem. *ESAIM Control Optim. and Calc. Var.* 2006; 12(2): 294-310.
17. Zhu J, Trélat E, Cerf M. Geometric Optimal Control and Applications to Aerospace. *Pacific Journal of Mathematics for Industry* 2017; 9(8). doi: 10.1186/s40736-017-0033-4

18. Caillau JB, Cots O, Gergaud J. Differential pathfollowing for regular optimal control problems. *Optim. Methods Softw.* 2012; 27(2): 177-196.

19. Malisani P, Chaplais F, Petit N. An interior penalty method for optimal control problems with state and input constraints of nonlinear systems. *Optim. Control Appl. Meth.* 2016; 37: 3-33.

20. Grant M, Boyd S. CVX: Matlab Software for Disciplined Convex Programming, version 2.1. http://cvxr.com/cvx; 2014.

21. Grant M, Boyd S. Graph implementations for nonsmooth convex programs. In: Blondel V, Boyd S, Kimura H., eds. *Recent Advances in Learning and Control* Lecture Notes in Control and Information Sciences. Springer-Verlag Limited. 2008 (pp. 95–110). http://stanford.edu/~boyd/graph_dcp.html.

**How to cite this article:** A. Herasimenka, L. Dell'Elce, J.-B. Caillau, and J.-B. Pomet (xxxx), xxxx, JJJJ, xxxx;yy:zz–zz.
17. Zhu J, Trélat E, Cerf M. Geometric Optimal Control and Applications to Aerospace. *Pacific Journal of Mathematics for Industry* 2017; 9(8). doi: 10.1186/s40736-017-0033-4

18. Caillau JB, Cots O, Gergaud J. Differential pathfollowing for regular optimal control problems. *Optim. Methods Softw.* 2012; 27(2): 177-196.

19. Malisani P, Chaplais F, Petit N. An interior penalty method for optimal control problems with state and input constraints of nonlinear systems. *Optim. Control Appl. Meth.* 2016; 37: 3-33.

20. Grant M, Boyd S. CVX: Matlab Software for Disciplined Convex Programming, version 2.1. http://cvxr.com/cvx; 2014.

21. Grant M, Boyd S. Graph implementations for nonsmooth convex programs. In: Blondel V, Boyd S, Kimura H., eds. *Recent Advances in Learning and Control* Lecture Notes in Control and Information Sciences. Springer-Verlag Limited. 2008 (pp. 95–110). http://stanford.edu/~boyd/graph_dcp.html.

**How to cite this article:** A. Herasimenka, L. Dell'Elce, J.-B. Caillau, and J.-B. Pomet (xxxx), xxxx, JJJJ, xxxx;yy:zz–zz.