Deterministic Trajectory Design and Attitude Maneuvers of Gradient-Index Solar Sail in Interplanetary Transfers

Abstract

A refractive sail is a special type of solar sail concept, whose membrane exposed to the Sun’s rays is covered with an advanced engineered film made of micro-prisms. Unlike the well-known reflective solar sail, an ideally flat refractive sail is able to generate a nonzero thrust component along the sail’s nominal plane even when the Sun’s rays strike that plane perpendicularly, that is, when the solar sail attitude is Sun-facing. This particular property of the refractive sail allows heliocentric orbital transfers between orbits with different values of the semilatus rectum while maintaining a Sun-facing attitude throughout the duration of the flight. In this case, the sail control is achieved by rotating the structure around the Sun–spacecraft line, thus reducing the size of the control vector to a single (scalar) parameter. A gradient-index solar sail (GIS) is a special type of refractive sail, in which the membrane film design is optimized though a transformation optics-based method. In this case, the membrane film is designed to achieve a desired refractive index distribution with the aid of a waveguide array to increase the sail efficiency. This paper analyzes the optimal transfer performance of a GIS with a Sun-facing attitude (SFGIS) in a series of typical heliocentric mission scenarios. In addition, this paper studies the attitude control of the Sun-facing GIS using a simplified mathematical model, in order to investigate the effective ability of the solar sail to follow the (optimal) variation law of the rotation angle around the radial direction.

1. Introduction

The gradient-index solar sail (GIS) [1] is an innovative solar sail concept that can be considered as a kind of evolution of the refractive solar sail model [2,3], whose thin exposed membrane uses an advanced metamaterial, composed of metals and dielectrics, to convert sunlight into propulsive force. This interesting type of propellantless propulsion system was recently proposed by Firuzi et al. [1], who discussed a simplified thrust model for preliminary mission design and also analyzed the problem of trajectory optimization by considering spacecraft dynamics in a classical, Cartesian, heliocentric-ecliptic reference frame. From the thrust vector point of view, as in the case of the Swartzlander’s diffractive sail [4,5], a spacecraft whose primary propulsion system is GIS-based is able to create a transverse thrust component (which belongs to the nominal plane of the sail) even when the solar sail is in a Sun-facing condition, that is, in a configuration in which the direction normal to the plane of the sail coincides with the Sun–spacecraft line. In this case, the only control parameter is the sail roll angle around the radial direction, the value of which can be controlled using a robust sail attitude control system based, for example, on a roll stabilizer bar [6].

This paper analyzes the deterministic trajectories of an interplanetary spacecraft propelled by a GIS with a Sun-facing attitude (SFGIS) in a typical heliocentric transfer, considering a single control parameter, namely, the roll angle of the solar sail. In particular, two different guidance laws are considered: (1) a case where the sail roll angle can be freely chosen by the control system; and (2) a case where during flight, the roll angle can be chosen within a discrete set of values. In the latter scenario, the paper analyzes the characteristics of the attitude maneuver required to appropriately change the sail roll angle, using the simplified model proposed by Wie et al. [6].

2. Trajectory Design and Optimization

This section describes the mathematical model used to analyze the motion of an interplanetary spacecraft propelled by an SFGIS. Specifically, the spacecraft motion is studied in a three-dimensional heliocentric mission scenario, in which the deep space vehicle moves between two assigned Keplerian orbits. The characteristics of the parking and target orbits are described in Section 3 along with the optimal spacecraft transfer trajectory, while the next subsection describes the SFGIS thrust model used to schematize the propulsive acceleration vector. The SFGIS thrust vector used in this paper is a simplified version of the mathematical model recently proposed by Firuzi et al. [1]. The expression of the propulsive acceleration vector is then used in Section 2.2 to describe the heliocentric dynamics of the SFGIS-propelled spacecraft using a set of non-singular orbital elements [7] proposed by Walker et al. [8], i.e., the modified equinoctial orbital elements (MEOEs). In this regard, in order to simplify and speed up the discussion, details regarding the transition between classical and non-singular orbital elements are here omitted, while the interested reader can refer to the recent literature by the authors [9] for a more in-depth discussion of the use of MEOEs in describing the heliocentric dynamics of a spacecraft equipped with a propellantless propulsion system as the (photonic) solar sail or the Electric Solar Wind Sail [10,11,12].

2.1. Thrust Vector Model of the Gradient-Index Sail with a Sun-Facing Attitude

Consider an SFGIS and neglect the billowing of the sail membrane due to the interaction of the sail exposed surface with the solar radiation pressure [13,14,15]. In this specific configuration, we define the sail nominal plane $\mathcal{P}$ as the plane that (ideally) contains the sail membrane and observe that the $\mathcal{P}$ plane is perpendicular to the Sun–spacecraft line, because the sail attitude is Sun-facing; see the conceptual sketch in Figure 1.

Keeping in mind the scheme of Figure 1, we introduce the unit vector $\widehat{\mathit{n}}$ normal to the plane $\mathcal{P}$ in the direction opposite to the Sun, and the unit vector $\widehat{\mathit{m}}$ indicating a fixed direction along the sail nominal plane. In particular, the unit vector $\widehat{\mathit{n}}$ belongs to the Sun–spacecraft line, while the direction of $\widehat{\mathit{m}}$ (which is fixed in $\mathcal{P}$ by assumption) depends on the arrangement of the waveguide arrays covering the membrane of the solar sail [1]. According to Firuzi et al. [1], the SFGIS-induced propulsive acceleration vector ${\mathit{a}}_p$ can be written as

$${\mathit{a}}_p=a_c{\left({\displaystyle \frac{r_{\oplus }}{r}}\right)}^2\left[{\eta }_n\widehat{\mathit{n}}+{\eta }_m\widehat{\mathit{m}}\right]$$

(1)

where r is the spacecraft’s distance from the Sun, $r_{\oplus }=1\mathrm{AU}$ is a reference distance that coincides with the average Earth–Sun distance, and $a_c$ is the characteristic acceleration (which is the typical performance parameter in solar sail trajectory design) defined as the maximum value of $\parallel {\mathit{a}}_p\parallel$ when the spacecraft’s distance from the Sun is equal to $r_{\oplus }$. In Equation (1), the two terms $\{{\eta }_n,{\eta }_m\}$ are dimensionless coefficients that depend on the specific design of the waveguide array covering the sail membrane [1]. According to Ref. [1], the value of the pair $\{{\eta }_n,{\eta }_m\}$ can be obtained by using a semi-analytical procedure employing the ray tracing technique. The result is [1]

$${\eta }_n=0.6299,{\eta }_m=0.7767$$

(2)

with $\sqrt{{\eta }_n^2+{\eta }_m^2}=1$, according to the definition of the characteristic acceleration $a_c$.

The expression of the propulsive acceleration vector in Equation (1) allows us to make some important considerations. First, note that the thrust vector is fixed in a body (spacecraft) reference frame. Consequently, it is not possible to obtain the condition $\parallel {\mathit{a}}_p\parallel = 0$ during the flight. In other terms, the heliocentric trajectory of a spacecraft propelled by an SFGIS contains no coasting arcs. This interesting aspect is similar to the one found in the mathematical model describing the thrust vector of a diffractive sail [4,5,16,17] with a Sun-facing attitude, as discussed in Ref. [18]. More precisely, the main difference between the thrust model described by Equation (1) and the one proposed by Dubill and Swartzlander [18] is the value of the so-called “cone angle” $\alpha$, i.e., the angle between the direction of the propulsive acceleration vector and the Sun–spacecraft line. In fact, in a Sun-facing diffractive sail, the value of the cone angle is ${45}^{\circ }$, while for the case of an SFGIS, according to Equations (1) and (2), we have

$$\alpha =arctan\left({\displaystyle \frac{{\eta }_m}{{\eta }_n}}\right)\simeq {51}^{\circ }$$

(3)

which is therefore a constant of motion.

Second, since the sail nominal plane $\mathcal{P}$ is assumed perpendicular to the Sun–spacecraft line at each instant, the normal unit vector $\widehat{\mathit{n}}$ coincides with the (radial) Sun-spacecraft unit vector $\widehat{\mathit{r}}$. Therefore, Equation (1) indicates that the component of ${\mathit{a}}_p$ along (perpendicular to) the Sun–spacecraft direction, that is, the radial (or horizontal) component of the propulsive acceleration $a_{p_{\mathrm{R}}}$ (or $a_{p_{\mathrm{H}}}$) depends only on the distance of the spacecraft from the Sun, viz.,

$$\begin{array}{cc}\hfill a_{p_{\mathrm{R}}}& ={\mathit{a}}_p·\widehat{\mathit{n}}=a_c{\left({\displaystyle \frac{r_{\oplus }}{r}}\right)}^2{\eta }_n\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill a_{p_{\mathrm{H}}}& ={\mathit{a}}_p·\widehat{\mathit{m}}=a_c{\left({\displaystyle \frac{r_{\oplus }}{r}}\right)}^2{\eta }_m\hfill \end{array}$$

(5)

with the ratio $a_{p_{\mathrm{H}}}/a_{p_{\mathrm{R}}}={\eta }_m/{\eta }_n\simeq 1.233$ remaining fixed during flight.

2.2. Heliocentric Dynamics of the Spacecraft’s Center of Mass

Based on the mathematical model of the propulsive acceleration vector presented in the previous subsection, and taking into account the results obtained in Ref. [1], the thrust vectoring of an SFGIS is obtained by rotating the sail nominal plane around the Sun–spacecraft line. This is equivalent to rotating the sail nominal plane $\mathcal{P}$ around the direction of the unit vector $\widehat{\mathit{n}}$, that is, the direction of $\widehat{\mathit{r}}$.

The problem of spacecraft thrust vectoring can be more easily studied by introducing an orbital reference frame. For this purpose, consider a classical Radial–Transverse–Normal (RTN) reference frame ${\mathcal{T}}_{\mathrm{RTN}}$, whose origin coincides with the spacecraft’s center of mass C, in which ${\widehat{\mathit{i}}}_{\mathrm{R}}\equiv \widehat{\mathit{r}}$ is the radial unit vector, ${\widehat{\mathit{i}}}_{\mathrm{T}}$ is the transverse unit vector, and ${\widehat{\mathit{i}}}_{\mathrm{N}}$ is the normal unit vector. In particular, the direction of ${\widehat{\mathit{i}}}_{\mathrm{N}}$ coincides with that of the spacecraft angular momentum vector, while $({\widehat{\mathit{i}}}_{\mathrm{R}},{\widehat{\mathit{i}}}_{\mathrm{T}})$ coincides with the plane of the osculating orbit; see Figure 2.

With reference to the sketch of Figure 2, we introduce the sail clock angle $\delta \in {[0,360]}^{\circ }$ defined as the angle between the direction of the transverse unit vector ${\widehat{\mathit{i}}}_{\mathrm{T}}$ and the direction of the (body-fixed) unit vector $\widehat{\mathit{m}}$. The clock angle is measured in the $\mathcal{P}$-plane counterclockwise from the direction of ${\widehat{\mathit{i}}}_{\mathrm{T}}$. According to Equation (1) and recalling the expressions of $\{a_{p_{\mathrm{R}}},a_{p_{\mathrm{H}}}\}$ in Equations (4) and (5), the angle $\delta$ allows the transverse ($a_{p_{\mathrm{T}}}$) and normal ($a_{p_{\mathrm{N}}}$) component of the propulsive acceleration vector to be written as

$$\begin{array}{cc}\hfill a_{p_{\mathrm{T}}}& ={\mathit{a}}_p·{\widehat{\mathit{i}}}_{\mathrm{T}}=a_c{\left({\displaystyle \frac{r_{\oplus }}{r}}\right)}^2{\eta }_{m}cos\delta \hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill a_{p_{\mathrm{N}}}& ={\mathit{a}}_p·{\widehat{\mathit{i}}}_{\mathrm{N}}=a_c{\left({\displaystyle \frac{r_{\oplus }}{r}}\right)}^2{\eta }_{m}sin\delta \hfill \end{array}$$

(7)

with $\sqrt{a_{p_{\mathrm{N}}}^2+a_{p_{\mathrm{T}}}^2}=a_{p_{\mathrm{H}}}$, where $a_{p_{\mathrm{H}}}$ is given by Equation (5), so that

$${\left[{\mathit{a}}_p\right]}_{\mathrm{RTN}}=a_c{\left({\displaystyle \frac{r_{\oplus }}{r}}\right)}^2\left[\begin{array}{c}{\eta }_n\\ {\eta }_{m}cos\delta \\ {\eta }_{m}sin\delta \end{array}\right]$$

(8)

The orientation of the thrust vector is therefore defined by a single (scalar) control variable, namely, the sail clock angle $\delta$, whose value indicates the direction of the horizontal component of the propulsive acceleration $a_{p_{\mathrm{H}}}$ with respect to the direction of ${\widehat{\mathit{i}}}_{\mathrm{T}}$. Note that, when $\delta ={\{0,180\}}^{\circ }$ (or $\delta ={\{90,270\}}^{\circ }$) the thrust vector belongs to a plane coincident with (or perpendicular to) the plane of the osculating orbit. The ratio of the generic component $\{a_{p_{\mathrm{R}}},a_{p_{\mathrm{T}}},a_{p_{\mathrm{N}}}\}$ to the characteristic acceleration $a_c$, when the solar distance is equal to $r_{\oplus }$, is drawn in Figure 3 as a function of the clock angle.

The expressions of the three components of ${\mathit{a}}_p$ in the RTN reference frame, see Equation (8), can be used to write the spacecraft equations of motion in a heliocentric scenario in which ${\mu }_{\odot }$ is the Sun’s gravitational parameter. To this end, the spacecraft state vector is described in terms of the six modified equinoctial orbital elements $\{p,f,g,h,k,L\}$, the geometric interpretation of which is discussed in Refs. [8]. According to the mathematical model proposed by Betts [19] and using Equations (9)–(12) of Ref. [9], the heliocentric spacecraft equations of motion are

$${\displaystyle \frac{\mathrm{d}p}{\mathrm{d}t}}={\displaystyle \frac{2p}{1+fcosL+gsinL}}\sqrt{{\displaystyle \frac{p}{{\mu }_{\odot }}}}{a}_{p_{\mathrm{T}}}$$

(9)

$$\begin{array}{c}{\displaystyle \frac{\mathrm{d}f}{\mathrm{d}t}}=sinL\sqrt{{\displaystyle \frac{p}{{\mu }_{\odot }}}}{a}_{p_{\mathrm{R}}}+{\displaystyle \frac{\left(2+fcosL+gsinL\right)cosL+f}{1+fcosL+gsinL}}\sqrt{{\displaystyle \frac{p}{{\mu }_{\odot }}}}{a}_{p_{\mathrm{T}}}+\hfill \\ -{\displaystyle \frac{g\left(hsinL-kcosL\right)}{1+fcosL+gsinL}}\sqrt{{\displaystyle \frac{p}{{\mu }_{\odot }}}}{a}_{p_{\mathrm{N}}}\hfill \end{array}$$

(10)

$$\begin{array}{c}{\displaystyle \frac{\mathrm{d}g}{\mathrm{d}t}}=-cosL\sqrt{{\displaystyle \frac{p}{{\mu }_{\odot }}}}{a}_{p_{\mathrm{R}}}+{\displaystyle \frac{\left(2+fcosL+gsinL\right)sinL+g}{1+fcosL+gsinL}}\sqrt{{\displaystyle \frac{p}{{\mu }_{\odot }}}}{a}_{p_{\mathrm{T}}}+\hfill \\ +{\displaystyle \frac{f\left(hsinL-kcosL\right)}{1+fcosL+gsinL}}\sqrt{{\displaystyle \frac{p}{{\mu }_{\odot }}}}{a}_{p_{\mathrm{N}}}\hfill \end{array}$$

(11)

$${\displaystyle \frac{\mathrm{d}h}{\mathrm{d}t}}={\displaystyle \frac{\left(1+h^2+k^2\right)cosL}{2\left(1+fcosL+gsinL\right)}}\sqrt{{\displaystyle \frac{p}{{\mu }_{\odot }}}}{a}_{p_{\mathrm{N}}}$$

(12)

$${\displaystyle \frac{\mathrm{d}k}{\mathrm{d}t}}={\displaystyle \frac{\left(1+h^2+k^2\right)sinL}{2\left(1+fcosL+gsinL\right)}}\sqrt{{\displaystyle \frac{p}{{\mu }_{\odot }}}}{a}_{p_{\mathrm{N}}}$$

(13)

$${\displaystyle \frac{\mathrm{d}L}{\mathrm{d}t}}={\displaystyle \frac{hsinL-kcosL}{1+fcosL+gsinL}}\sqrt{{\displaystyle \frac{p}{{\mu }_{\odot }}}}{a}_{p_{\mathrm{N}}}+\sqrt{{\mu }_{\odot }p}{\left({\displaystyle \frac{1+fcosL+gsinL}{p}}\right)}^2$$

(14)

where $\{a_{p_{\mathrm{R}}},a_{p_{\mathrm{T}}},a_{p_{\mathrm{N}}}\}$ are given by Equation (8) as a function of the performance characteristics ($a_c$, ${\eta }_n$, and ${\eta }_m$), the dimensionless control term ($\delta$), and the solar distance (r). The latter can be written as a function of the modified equinoctial orbital elements as [20]

$$r={\displaystyle \frac{p}{1+fcosL+gsinL}}$$

(15)

The initial conditions of the differential Equations (9)–(14) depend on the characteristics of the spacecraft parking orbit at the initial time instant $t_0=0$. In this regard, as discussed in detail in Section 3, we assume that the shape of the spacecraft parking orbit coincides with that of the Earth’s heliocentric orbit. Specifically, the orbital elements of Earth’s heliocentric orbit were retrieved from JPL’s Horizon system [21] on 1 August 2024. The control law $\delta =\delta \left(t\right)$ is discussed in the next two subsections.

2.3. Trajectory Optimization: Unconstrained Case

The time variation in the sail clock angle during the interplanetary transfer is obtained by solving an optimization problem in which the performance index to maximize is $J=-\mathsf{\Delta }t$, where $\mathsf{\Delta }t=t_f-t_0\equiv t_f$ is the flight time ($t_f$ is the final time instant) required to complete the heliocentric transfer between two assigned Keplerian orbits. In this case, the angular position of the spacecraft both along the parking orbit (at time $t_0$) and along the target orbit (at time $t_f$) is left free, so that the problem considered here is a classical obit-to-orbit optimal transfer without ephemeris constraints.

The procedure used to solve the optimization problem parallels the approach discussed in Ref. [9]. More precisely, an indirect method [22,23,24] is used to determine the spacecraft rapid transfer trajectory, while the Pontryagin Maximum Principle (PMP) [25,26,27] is employed to obtain the optimal control law in terms of the variation in the sail clock angle during the flight. The approach used is standard and has been discussed in several articles by authors. For this reason, only the calculation of the optimal control law is discussed here, since the analytical expression obtained is new and extends the literature results related to the (optimal) performance of a gradient-index solar sail [1]. In this subsection, we assume that the sail clock angle is freely selectable within the range ${[0,360]}^{\circ }$, while the next subsection discusses the case where $\delta$ can take a finite set of values.

In order to use the PMP to obtain the optimal control law, we first determine the expression of the part of the Hamiltonian function ${\mathcal{H}}_c$ that depends explicitly on the control parameter $\delta$. To this end, firstly note that the clock angle appears in the expressions of $a_{p_{\mathrm{T}}}$ and $a_{p_{\mathrm{N}}}$, see Equations (6) and (7), and observe that the two components $\{a_{p_{\mathrm{T}}},a_{p_{\mathrm{N}}}\}$ appear in the expressions of the time-derivative of the generic modified equinoctial orbital element, see Equations (9)–(14). Introducing the six adjoint variables $\{{\lambda }_p,{\lambda }_f,{\lambda }_g,{\lambda }_h,{\lambda }_k,{\lambda }_L\}$ [22,23], we obtain the following expression of ${\mathcal{H}}_c$

$${\mathcal{H}}_c={\displaystyle \frac{a_c{\left({\displaystyle \frac{r_{\oplus }}{r}}\right)}^2\sqrt{{\displaystyle \frac{p}{{\mu }_{\odot }}}}}{1+fcosL+gsinL}}\left[Tcos\delta +Nsin\delta \right]$$

(16)

where $\{T,N\}$ are two auxiliary functions defined as

$$\begin{array}{c}T\triangleq 2p{\lambda }_p{\eta }_m+\left[\left(2+fcosL+gsinL\right)cosL+f\right]{\lambda }_f{\eta }_m+\hfill \\ +\left[\left(2+fcosL+gsinL\right)sinL+f\right]{\lambda }_g{\eta }_m\hfill \end{array}$$

(17)

$$\begin{array}{c}N\triangleq -g\left(hsinL-kcosL\right){\lambda }_f{\eta }_n+f\left(hsinL-kcosL\right){\lambda }_g{\eta }_n+cosL{\displaystyle \frac{1+h^2+k^2}{2}}{\lambda }_h{\eta }_n+\hfill \\ +sinL{\displaystyle \frac{1+h^2+k^2}{2}}{\lambda }_k{\eta }_n+\left(hsinL-kcosL\right){\lambda }_L{\eta }_n\hfill \end{array}$$

(18)

According to the PMP, the optimal value of the clock angle is the one that maximizes at each time instant the value of ${\mathcal{H}}_c$. Bearing in mind Equation (16), this corresponds to selecting the value of $\delta$ that maximizes $\left(Tcos\delta +Nsin\delta \right)$, viz.,

$$cos\delta ={\displaystyle \frac{T}{\sqrt{T^2+N^2}}},sin\delta ={\displaystyle \frac{N}{\sqrt{T^2+N^2}}}$$

(19)

where $\{T,N\}$ are given by Equations (17) and (18).

2.4. Trajectory Optimization: Case of Constrained Clock Angle

Consider now the case in which the clock angle can take only a finite number of values in the admissible interval ${[0,360]}^{\circ }$. Assume, for example, that the angle $\delta$ can take four values only, that is, $\delta \in \mathcal{I}$, where $\mathcal{I}$ is a discrete set. In this case, at a given point of the transfer trajectory, the function ${\mathcal{H}}_c$ can take only four values depending on the selected sail clock angle. Therefore, the optimal clock angle can be easily obtained using a (numerical) vector sorting procedure [28] that finds the value of $\delta$ that maximizes ${\mathcal{H}}_c$.

3. Mission Application and Trajectory Simulations

The optimal control law of Equation (19) was used to obtain the minimum-time orbit-to-orbit transfer in some typical heliocentric mission scenarios. In all the numerical simulations, the differential equation was integrated using a PECE solver [29] based on the Adams–Bashforth method with a tolerance of ${10}^{-12}$, while the associated boundary value problem was solved by a shooting procedure [28] with a tolerance of ${10}^{-6}$. The initial guess was obtained by adapting the procedure recently described in Ref. [30].

The optimization procedure was first validated using, as a comparison, the results reported in Ref. [1] which, however, refer to ephemeris-constrained transfer scenarios. Therefore, it was reasonable to expect that the flight times obtained with the technique proposed in this paper (which refer to a less constrained orbit-to-orbit transfer) would be lower than those reported in Ref. [1]. In the case of an unconstrained sail clock angle, consider an Earth–Venus mission scenario, and assume a characteristic acceleration $a_c=0.175\mathrm{mm}/{\mathrm{s}}^2$ [1]. The proposed procedure gave a minimum flight time of about $434.8\mathrm{days}$, while Ref. [1] indicated a (constrained) flight time of $456.3\mathrm{days}$ with a departure on 29 April 2021 and arrival on 27 July 2022. The time variation in the sail clock angle during the flight is shown in Figure 4, while Figure 5 shows the optimal transfer trajectory in a heliocentric-ecliptic Cartesian reference frame [31].

The second mission scenario illustrated in Ref. [1] is a heliocentric orbit transfer from Earth to asteroid 433 Eros. In that case, considering again a characteristic acceleration of $0.175\mathrm{mm}/{\mathrm{s}}^2$, Firuzi et al. [1] indicated an optimal (constrained) flight time of $1207.3\mathrm{days}$, departing on 11 November 2024 and arriving on 2 March 2028. The procedure proposed in this paper gave a minimum flight time (in an orbit-to-orbit transfer), of $1125\mathrm{days}$, with the optimal time variation in the clock angle shown in Figure 6 and the optimal transfer trajectory shown in Figure 7. Note that the function $\delta =\delta \left(t\right)$ drawn in Figure 6 is only apparently discontinuous, due to the range of variation chosen to describe the sail clock angle.

The simulations performed demonstrated that the optimization method and control law described above provided numerical results consistent with the literature and therefore could be used to analyze new mission scenarios, as discussed in the next subsection.

3.1. Case of Unconstrained Clock Angle

Consider the general case of an unconstrained sail clock angle, in which the optimal control law is given by Equation (19). The proposed optimization procedure was used to parametrically analyze the Earth–Venus transfer as a function of the characteristic acceleration value $a_c\in [0.1,0.2]\mathrm{mm}/{\mathrm{s}}^2$. The results are summarized in Figure 8, which shows that the flight time $\mathsf{\Delta }t$ depends strongly on $a_c$. Note, in fact, that changing from a value of $a_c=0.175\mathrm{mm}/{\mathrm{s}}^2$ to a value of $a_c=0.2\mathrm{mm}/{\mathrm{s}}^2$ results in a reduction in the flight time of more than $60\mathrm{days}$.

The second mission scenario analyzed an Earth–Mars interplanetary transfer. In that case, we again assumed the same characteristic acceleration reported in Ref. [1], that is, $a_c=0.175\mathrm{mm}/{\mathrm{s}}^2$. The minimum flight time for an orbit-to-orbit transfer was about $752\mathrm{days}$, and the SFGIS-propelled spacecraft completed one revolution around the Sun before reaching the target planet orbit, as shown in Figure 9. The optimal control law $\delta =\delta \left(t\right)$ is shown in Figure 10, where the clock angle remains close to the value $\delta =0$ (or $\delta ={360}^{\circ }$, which is equivalent in terms of the direction of $\widehat{\mathit{m}}$ in the RTN reference frame) during the entire flight.

The last scenario analyzed in this subsection was an Earth–Mercury transfer with $a_c=0.175\mathrm{mm}/{\mathrm{s}}^2$. The low-thrust transfer to Mercury is a challenging problem, typically requiring a long flight time and numerous revolutions around the Sun to complete [32,33,34,35]. In this case, an SFGIS-propelled spacecraft required about $780\mathrm{days}$ to reach Mercury’s heliocentric orbit, with four complete revolutions around the Sun. The spacecraft transfer trajectory and the time variation in the sail clock angle are shown in Figure 11 and Figure 12, respectively.

3.2. Case of Constrained Clock Angle

The performance of an SFGIS-powered spacecraft in a case where the sail clock angle is constrained to assume only the values included in a discrete set $\mathcal{I}$ was analyzed considering an Earth–Venus mission scenario with $a_c=0.175\mathrm{mm}/{\mathrm{s}}^2$. The results obtained in the same transfer, but with the unconstrained clock angle, summarized in Figure 4 and Figure 5, gave a minimum flight time of $434.8\mathrm{days}$.

Bearing in mind the shape of the graph $\delta =\delta \left(t\right)$ of Figure 4, let us consider two possible sets of admissible values of the sail clock angle, viz.,

$${\mathcal{I}}_{①}={\{180,210,240,270,300\}}^{\circ },{\mathcal{I}}_{②}={\{180,240,300\}}^{\circ }$$

(20)

In other terms, when $\mathcal{I}={\mathcal{I}}_{①}$ (or $\mathcal{I}={\mathcal{I}}_{②}$) the sail clock angle can only take five (or three) different values during the transfer.

In this $\delta$-constrained mission scenario, the optimization procedure gave a minimum flight time of $\mathsf{\Delta }{t}_{①}=437.7\mathrm{days}$ when $\mathcal{I}={\mathcal{I}}_{①}$, and $\mathsf{\Delta }{t}_{②}=441.9\mathrm{days}$ when $\mathcal{I}={\mathcal{I}}_{②}$. Note that in both cases, the optimal flight time was (slightly) higher than the value obtained in the unconstrained case, with $\mathsf{\Delta }{t}_{①}<\mathsf{\Delta }{t}_{②}$ as expected, while the difference was only a few days. According to Figure 13, the optimal control law in the $\delta$-constrained mission scenario presented a number of attitude maneuvers in which the direction of the unit vector $\widehat{\mathit{m}}$ in the RTN frame was suitably changed. The attitude maneuvers were performed at a solar distance r that could be retrieved from Figure 14. Specifically, there were 11 maneuvers when $\mathcal{I}={\mathcal{I}}_{①}$ and 6 maneuvers when $\mathcal{I}={\mathcal{I}}_{②}$. From the point of view of trajectory design, these attitude maneuvers were considered instantaneous, although obviously, the rotation of the SFGIS around the Sun–spacecraft direction required a finite time interval. This analysis is the subject of the next section.

4. Spacecraft Attitude Maneuvers

To address the problem of attitude maneuvers, it was assumed that the SFGIS was equipped with a pair of rectangular control vanes, mounted at the end of the sides of the SFGIS symmetrically to the diagonal of the sail; see Figure 15. These control vanes were reflective surfaces (essentially small solar sails) capable of generating a force orthogonal to the plane of their surface. They were also capable of generating a torque component along the axis normal to the sail nominal plane by rotating the control angle of the vane, i.e., ${\beta }_1$ for surface 1 and ${\beta }_2$ for surface 2, as shown in Figure 15. By construction, the rotations were positive (or negative) if the surfaces rotated toward the positive (or negative) semi-axis z.

Using the reference system described in Figure 15 and with the sign convention adopted for the vane control angles, the vector arms of the forces due to the solar radiation pressure acting on the two vanes are

$${\mathit{p}}_1={\displaystyle \frac{1}{2}}\left[\begin{array}{c}l+bcos{\beta }_1\\ l\left(1-\gamma \right)\\ bsin{\beta }_1\end{array}\right],{\mathit{p}}_2=-{\displaystyle \frac{1}{2}}\left[\begin{array}{c}l+bcos{\beta }_2\\ l\left(1-\gamma \right)\\ -bsin{\beta }_2\end{array}\right]$$

(21)

where l is the length of the SFGIS side, b is the width of the two control surfaces, and $\gamma l$ (with $\gamma \in (0,0.5]$) represents their length. The two force vectors due to the solar radiation pressure acting on the control surfaces are given by

$${\mathit{F}}_1=2PA{cos}^2{\beta }_1\left[\begin{array}{c}sin{\beta }_1\\ 0\\ -cos{\beta }_1\end{array}\right],{\mathit{F}}_2=2PA{cos}^2{\beta }_2\left[\begin{array}{c}-sin{\beta }_2\\ 0\\ -cos{\beta }_2\end{array}\right]$$

(22)

where

$$P\triangleq P_{\oplus }{\left({\displaystyle \frac{r_{\oplus }}{r}}\right)}^2\mathrm{with}{P}_{\oplus }=4.5391\times {10}^{-6}\mathrm{Pa}$$

(23)

is the solar radiation pressure at a distance r from the Sun, $A\triangleq \gamma lb$ is the area of each control surface, while the factor two in Equation (22) means that the control vanes are assumed perfectly reflective. Note that, in principle, ${\beta }_1$ and ${\beta }_2$ can take different values. The torque vectors are therefore

$$\begin{array}{cc}\hfill {\mathit{M}}_1& ={\mathit{p}}_1\times {\mathit{F}}_1=PA{cos}^2{\beta }_1\left[\begin{array}{c}-cos{\beta }_{1}l\left(1-\gamma \right)\\ b{sin}^2{\beta }_1+cos{\beta }_1\left(l+bcos{\beta }_1\right)\\ -sin{\beta }_{1}l\left(1-\gamma \right)\end{array}\right]\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill {\mathit{M}}_2& ={\mathit{p}}_2\times {\mathit{F}}_2=PA{cos}^2{\beta }_2\left[\begin{array}{c}cos{\beta }_{2}l\left(1-\gamma \right)\\ -b{sin}^2{\beta }_2-cos{\beta }_2\left(l+bcos{\beta }_2\right)\\ -sin{\beta }_{2}l\left(1-\gamma \right)\end{array}\right]\hfill \end{array}$$

(25)

from which it is easily concluded that the first two components of $\mathit{M}\triangleq {\mathit{M}}_1+{\mathit{M}}_2$ are zero only if ${\beta }_1={\beta }_2=\beta$. Under this assumption, the total torque becomes

$$\mathit{M}=2PA{cos}^2\beta \left[\begin{array}{c}0\\ 0\\ -sin\beta l\left(1-\gamma \right)\end{array}\right]$$

(26)

Therefore, a rotation of the same angle of the two control vanes generates a pure windmill torque around the z-axis (orthogonal to the plane of the SFGIS) given by

$$M\triangleq \mathit{M}·\widehat{\mathit{k}}=-2P_{\oplus }{\left({\displaystyle \frac{r_{\oplus }}{r}}\right)}^{2}b{l}^2\gamma \left(1-\gamma \right){cos}^2\beta sin\beta$$

(27)

where $\beta \in [-{\beta }_{max},{\beta }_{max}]$ is the vane control angle. It is assumed that ${\beta }_{max}$ coincides with the value of $\beta$ that maximizes the torque magnitude, viz.,

$${\beta }_{max}\triangleq arcsin\left({\displaystyle \frac{\sqrt{3}}{3}}\right)\simeq {35.26}^{\circ }$$

(28)

Equivalently, Equation (27) may be rewritten by introducing the dimensionless control variable $u\in [-1,1]$, defined as

$$u\triangleq {\displaystyle \frac{{cos}^2\beta sin\beta }{max\left({\cos }^2\beta sin\beta \right)}}$$

(29)

from which

$$M={\displaystyle \frac{4}{3\sqrt{3}}}{P}_{\oplus }{\left({\displaystyle \frac{r_{\oplus }}{r}}\right)}^{2}b{l}^2\gamma \left(1-\gamma \right)u$$

(30)

With reference to the simplified model discussed by Wie et al. [6], to estimate the time required for an SFGIS to perform a reorientation maneuver around the z-axis, it is sufficient to describe the attitude dynamics in the form

$$I\ddot{\delta }=M$$

(31)

where $I\triangleq m_{\mathrm{sail}}{l}^2/6$ is the moment of inertia of the SFGIS around z. Substituting Equation (30) into Equation (31) gives

$$\ddot{\delta }=ku,\mathrm{with}k\triangleq {\displaystyle \frac{8P_{\oplus }{\left({\displaystyle \frac{r_{\oplus }}{r}}\right)}^{2}b\gamma \left(1-\gamma \right)}{\sqrt{3}{m}_{\mathrm{sail}}}}$$

(32)

where k is a design parameter that encompasses information about the geometry of the control surfaces, the mass of the sail, and the distance from the Sun. The problem is to determine a control law $u=u\left(t\right)$ (with $\left|u\right|\le 1$) capable of tracking a generic (constant) reference value of $\delta$, i.e., ${\delta }_{\mathrm{ref}}$. To this end, consider the block diagram shown in Figure 16, where $G_{\mathrm{c}}$ represents the controller, $G_{\mathrm{s}}$ represents the system (i.e., the SFGIS attitude dynamics), $e\triangleq {\delta }_{\mathrm{ref}}-\delta$ is the error signal, and $\overline{u}$ is the (unsaturated) input to the system.

Using Equation (32), the transfer function that describes the dynamics of the system is

$$G_{\mathrm{s}}={\displaystyle \frac{k}{s^2}}$$

(33)

where s is the complex variable, while the assumption is made that $G_{\mathrm{c}}$ has the form of a proportional–derivative controller with a first-order filter on the derivative term (PDF), viz.,

$$G_{\mathrm{c}}=K_p+{\displaystyle \frac{K_{d}s}{T_{f}s+1}}$$

(34)

The choice of control law parameters ($K_p$, $K_d$ and $T_f$) was achieved through a numerical optimization procedure, in which the objective function to be minimized was the settling time of the response of the closed-loop system in Figure 16. More precisely, the value of the gain crossover frequency (${\omega }_c$) of the open-loop system was varied iteratively, and for each value of ${\omega }_c$, the calculation of $K_p$, $K_d$, and $T_f$ was performed with the help of MATLAB R2024a’s built-in function pidtune, which automatically adjusts the controller parameters to balance performance (i.e., response time) against robustness (i.e., stability margins) [36]. The optimization problem was solved by using MATLAB’s built-in function fminsearch.

5. Numerical Results of the Spacecraft Attitude Maneuvers

This section discusses the results of the spacecraft attitude control problem, with the goal of estimating the optimal time required for the SFGIS to complete a reorientation maneuver. Recall that, with reference to the results obtained in the cases of Earth–Venus transfers with constrained clock angles, the reorientation maneuvers were on the order of 60 deg each. Based on that, the minimum settling time $t_s$, defined as the time instant such that $\left|e\right|\le 0.01\mathrm{rad}$ for $t\ge t_s$, was estimated with respect to a rotation maneuver of $1\mathrm{rad}$ around z. Numerical simulations were performed for values of k ranging between $5\times {10}^{-8}\mathrm{rad}/{\mathrm{s}}^2$ and $1.2\times {10}^{-6}\mathrm{rad}/{\mathrm{s}}^2$. To give an idea of the order of magnitude of the numbers used, consider a reorientation maneuver that occurs when the Sun–spacecraft distance is $r=r_{\oplus }$, assuming $b=0.2\mathrm{m}$, $\gamma =0.125$, and $m_{\mathrm{sail}}=5\mathrm{kg}$. In this case, we obtain $k\simeq 9.1724\times {10}^{-8}\mathrm{rad}/{\mathrm{s}}^2$. The simulation results for this particular case are shown in Figure 17, from which it can be seen that the maneuver takes about $120\min$ to complete.

Note that the angle $\delta$ reached the reference value through a time response having a maximum peak slightly greater than one, while the input to the system revealed that the control law was essentially of bang–off–bang type. These results are in agreement with what Wie found using a nonlinear PID control logic [37]. In all cases analyzed for different values of k, control laws characterized by a limited maximum peak and with bang–off–bang inputs were always found. Figure 18 shows the minimum settling time $t_s$ to perform a rotation of ${\delta }_{\mathrm{ref}}=1\mathrm{rad}$ around z as a function of k.

As might be expected, the time of the optimal reorientation maneuver decreased as k (and thus control angular acceleration) increased. For example, $t_s\simeq \{136,56,28\}\min$ were obtained when $k=\{5,30,120\}\times {10}^{-8}\mathrm{rad}/{\mathrm{s}}^2$, respectively. Finally, Figure 19 shows the variation in the optimal control parameters of the PDF controller with k.

6. Conclusions

A gradient-index solar sail is able to generate a nonzero thrust component along the nominal plane of the sail that is greater than that in the direction normal to the plane of the sail. This allows for heliocentric orbital transfers while keeping the plane of the sail constantly in a Sun-facing condition, that is, orthogonal to the direction of the Sun’s rays. This feature is very interesting from an application point of view because the Sun-facing attitude can be maintained passively, greatly simplifying sail maneuvering. In practice, the direction of thrust depends solely on a single control variable (the clock angle), i.e., the angle of rotation of the sail with respect to the radial direction. The performance of a gradient-index solar sail in a Sun-facing condition was studied from an optimal point of view by preliminarily determining the control law that minimized the total transfer time between two heliocentric orbits. Several optimal trajectories were studied, including Earth–Venus, Earth–Mars, Earth–Mercury transfers, and the optimal trajectory from Earth to asteroid 433 Eros. To further simplify the sail control, the optimal transfer problem was also studied in the case where the clock angle was constrained to take only a finite number of values (three or five depending on the cases studied). In an Earth–Venus transfer scenario, simulations revealed that the constrained clock angle case provided similar performance to the unconstrained problem, with an increase in flight time of about 1% only. These results were obtained under the assumption of instantaneous maneuvers. However, in the last part of the paper, the problem of estimating the time required for the sail to perform such maneuvers was studied. Assuming that the maneuvers were performed with the aid of two control vanes generating a pure windmill torque, the results were expressed as a function of a single parameter that took into account the surface of the control vanes, the mass of the sail, and the distance from the Sun. The simulations obtained showed that a sail rotation maneuver of 1 rad took about 2 h at a distance of $1\mathrm{AU}$ from the Sun. Since the characteristic transfer times are several hundred days, the assumption of instantaneous maneuvering is amply justified.