Real-time optimal control for attitude-constrained solar sailcrafts via
  neural networks
This work is devoted to generating optimal guidance commands in real time for
attitude-constrained solar sailcrafts in coplanar circular-to-circular
interplanetary transfers. Firstly, a nonlinear optimal control problem is
established, and necessary conditions for optimality are derived by the
Pontryagin's Minimum Principle. Under some assumptions, the attitude
constraints are rewritten as control constraints, which are replaced by a
saturation function so that a parameterized system is formulated to generate an
optimal trajectory via solving an initial value problem. This approach allows
for the efficient generation of a dataset containing optimal samples, which are
essential for training Neural Networks (NNs) to achieve real-time
implementation. However, the optimal guidance command may suddenly change from
one extreme to another, resulting in discontinuous jumps that generally impair
the NN's approximation performance. To address this issue, we use two co-states
that the optimal guidance command depends on, to detect discontinuous jumps. A
procedure for preprocessing these jumps is then established, thereby ensuring
that the preprocessed guidance command remains smooth. Meanwhile, the sign of
one co-state is found to be sufficient to revert the preprocessed guidance
command back into the original optimal guidance command. Furthermore, three NNs
are built and trained offline, and they cooperate together to precisely
generate the optimal guidance command in real time. Finally, numerical
simulations are presented to demonstrate the developments of the paper.
Kun Wang, Fangmin Lu, Zheng Chen, Jun Li
# Real-time optimal control for attitude-constrained solar sailcrafts via neural networks

###### Abstract

This work is devoted to generating optimal guidance commands in real time for attitude-constrained solar sailcrafts in coplanar circular-to-circular interplanetary transfers. Firstly, a nonlinear optimal control problem is established, and necessary conditions for optimality are derived by the Pontryagin's Minimum Principle. Under some mild assumptions, the attitude constraints are rewritten as control constraints, which are replaced by a saturation function so that a parameterized system is formulated. This allows one to generate an optimal trajectory via solving an initial value problem, making it efficient to collect a dataset containing optimal samples, which are essential for training Neural Networks (NNs) to achieve real-time implementation. However, the optimal guidance command may suddenly change from one extreme to another, resulting in discontinuous jumps that generally impair the NN's approximation performance. To address this issue, we use two co-states that the optimal guidance command depends on, to detect discontinuous jumps. A procedure for preprocessing these jumps is then established, thereby ensuring that the preprocessed guidance command remains smooth everywhere. Meanwhile, the sign of one co-state is found to be sufficient to revert the preprocessed guidance command back into the original optimal guidance command. Furthermore, three NNs are built and trained offline, and they cooperate together to precisely generate the optimal guidance command in real time. Finally, numerical simulations are presented to demonstrate the developments of the paper.

keywords: Solar sailing, Attitude constraints, Real-time optimal control, Neural networks

## 1 Introduction

Unlike conventional spacecrafts that use chemical or electric propellant to produce thrust during a flight mission, solar sailing exploits the Solar Radiation Pressure (SRP) to propel the solar sailcraft. Although the resulting propulsive force is smaller than that of chemical- or electric-based propulsion systems, the SRP offers an almost "infinite" energy source that
can be used to orient the solar sailcraft [1]. This makes it a very promising technology, especially for long-duration interplanetary transfer missions [2; 3]. While the study on solar sailing dates back to the 1920s, substantial progress in the development of solar sailing has not been achieved until recently. Following the first demonstrations of solar sailcraft in orbit by JAXA's IKAROS [4] and NASA's NanoSail-D2 [5] in 2010, the successes of LightSail 1 and 2 [6] have sparked renewed interests in this technology. These achievements have led to the planning of numerous solar sailcraft-based missions, including Solar Cruiser [6] and OKEANOS [7].

Designing the transfer trajectory for solar sailcrafts is one of the most fundamental problems. As the acceleration provided by the SRP is extremely small, the transfer time of the solar sailcraft in space is usually too long. Thus, it is fundamentally important to find the time-optimal control so that the solar sailcraft can be steered to its targeted orbit within minimum time. This is essentially equivalent to addressing a minimum-time control problem, which is a conventional Optimal Control Problem (OCP). Numerical methods for OCPs include indirect and direct methods [8]. The indirect method, based on the calculus of variations or Pontryagin's Minimum Principle (PMP), transforms the OCP into a Two-Point Boundary-Value Problem (TPBVP) according to the necessary conditions for optimality. The resulting TPBVP is then typically solved by Newton-like iterative methods. The indirect method has been applied to solve the time-optimal orbital transfer problem for solar sailcrafts in Refs. [9; 10; 11]. Mengali and Quarta [12] employed the indirect method to solve the three-dimensional time-optimal transfer problem for a non-ideal solar sailcraft with optical and parametric force models. Sullo _et al._[13] embedded the continuation method into the indirect method, so that the time-optimal transfer of the solar sailcraft was found from a previously designed low-thrust transfer trajectory. Recently, the solar sail primer vector theory, combined with the indirect method, was employed to design the flight trajectory that minimizes the solar angle over a transfer with a fixed time of flight [14]. On the other hand, direct methods reformulate the OCP (usually via direct shooting or pseudospectral methods) as a nonlinear programming problem, which can be solved using interior point or sequential quadratic programming method [8]. Solar sailcraft transfers to hybrid pole and quasi pole sitters on Mars and Venus were studied using a direct pseudospectral algorithm in Ref. [15]. In Ref. [16], the Gaussian quadrature pseudospectral optimization was combined with a slerped control parameterization method, and the transfer time was drastically reduced. Furthermore, a comparison study on application of indirect methods and direct methods to the time-optimal transfer for solar sailcrafts was conducted in Ref. [17]. In addition to these conventional methods, heuristic methods have also been utilized in the literature to solve the time-optimal transfer problem; see, e.g., Ref. [18].

Although the techniques mentioned in the previous paragraph are widely used, they are usually time-consuming and need appropriate initial guesses of co-state or state vector. Therefore, these methods are typically not suitable for onboard systems with limited computational resources. To overcome this issue, shape-based methods have been developed. The basic idea of the shape-based methods is to describe the geometric shape of trajectory using a set of mathematical expressions with some tunable parameters. The parameters are usually adjusted so that the mathematical expressions match the required boundary conditions.
In this regard, Peloni _et al._[19] expressed the trajectory as a function of four parameters for a multiple near-Earth-asteroid mission. A Radau pseudospectral method was leveraged to solve the resulting problem in terms of these four parameters. Since then, different shaping functions, such as Bezier curve [20] and Fourier series [21; 22], have been proposed for transfer problems. In order to further cope with constraints on the propulsive acceleration, some additional shape-based functions have been developed [23; 24]. It is worth mentioning that shape-based methods, albeit computationally efficient, often provide suboptimal solutions. Consequently, they are usually employed to provide initial guesses for direct methods.

According to the preceding review, existing approaches for solar sailcraft transfer suffer from computational burden, convergence issues, or solution suboptimality. In fact, solution optimality plays an important role in space exploration missions [25], and it is demanding for onboard systems to generate real-time solutions. Therefore, it is worth exploiting more effective methods capable of providing optimal solutions in real time. Thanks to the emerging achievements of artificial intelligence in recent decades, machine learning techniques have become a viable approach. Generally, there are two quite distinct machine learning based approaches. The first class, known as Reinforcement Learning (RL) involves training an agent to behave in a potentially changing environment via repeatedly interactions, and the goal is to maximize a carefully designed reward function. Owing to excellent generalization abilities of deep Neural Networks (NNs), deep RL algorithms have demonstrated great success in learning policies to guide spacecrafts in transfer missions; see, e.g., Refs [26; 27; 28]. Although the trained agent is generalizable to non-nominal conditions and robust to uncertainties, significant effort may be required to design an appropriate reward function, in order to derive a solution that is very close to optimal [29].

In contrast to RL, Behavioral Cloning (BC) aims to clone the observed expert behavior through training an NN on optimal state-action pairs. Such pairs are usually obtained by solving deterministic OCPs via indirect or direct methods. Recently, BC has been widely used to generate optimal guidance commands in real time, or reduce the computational time in aerospace applications, such as spacecraft landing [30; 31], hypersonic vehicle reentry [32; 33], missile interception [34], low-thrust transfer [35], end-to-end human-Mars entry, powered-descent, and landing problem [36], minimum-fuel geostationary station-keeping [37], time-optimal interplanetary rendezvous [38], and solar sailcraft transfer [39]. Specifically, in Ref. [39], some minimum-time orbital transfer problems with random initial conditions were solved by the indirect method to formulate the dataset for the state-guidance command pairs. Then, the dataset was used to train NNs within a supervised learning framework. Finally, the trained NNs were able to generate the optimal guidance command in real time. However, constraints on the solar sailcraft attitude were not considered in that paper. In practice, factors such as power generation and thermal control generally reduce the admissible variation range of the cone angle, which is equal to the sail attitude for a perfectly flat ideal sail [40; 41; 42]. In addition, the optimal guidance command for the solar sailcraft transfer problem may be discontinuous, which generally degrades the approximation performance of NNs. In fact, approximating discontinuous controls via NNs, even with many hidden layers and/or neurons, can be quite challenging, as shown by the numerical simulations in Refs. [43; 44; 45].
As a continuation effort to the work in Ref. [39], this paper considers a scenario where the solar sailcraft's attitude is constrained, and an NN-based method will be developed for real-time generation of the optimal guidance command. Firstly, the time-optimal problem for the solar sailcraft with constraints on the attitude is formulated. Then, necessary conditions for optimal trajectories are established by employing the PMP. These conditions allow using a saturation function to approximate the optimal guidance law. Unlike conventional optimization-based approaches, such as indirect methods and direct methods which often suffer from convergence issues [46], we formulate a parameterized system to facilitate the generation of datasets for training NNs. In this method, an optimal trajectory can be generated by solving an initial value problem instead of a TPBVP. As a consequence, a large number of optimal trajectories can be readily obtained.

However, discontinuous jumps in the optimal guidance command may defer one from obtaining perfect approximation. To resolve this issue, one viable method is to use NNs to approximate smooth co-states instead, as was done in Ref. [47]. Unfortunately, using co-states to train NNs is not reliable as the relationship between co-states and the optimal guidance command is highly nonlinear. This may result in magnified propagation errors even for small errors from the predicted co-states [48]. To avoid this difficulty, we propose to employ two co-states to preprocess the guidance command, smoothing out the discontinuous jumps. After preprocessing, the guidance command can be reverted to its original optimal form by examining the sign of one specific co-state. To achieve real-time generation of the optimal guidance command, we employ a system comprising three NNs that predict the optimal time of flight, the preprocessed guidance command, and one co-state. These three NNs cooperate together to generate the optimal guidance command in real time.

The remainder of the paper is structured as follows. The OCP is formulated in Section 2. Section 3 presents the optimality conditions derived from the PMP, and a parameterized system for optimal solutions is established. In Section 4, procedures for generating the dataset and preprocessing discontinuous jumps are introduced, and the scheme for generating the optimal guidance command in real time is detailed. Section 5 provides some numerical examples to demonstrate the developments of the paper. This paper finally concludes by Section 6.

## 2 Problem Formulation

### Solar Sailcraft Dynamics

We consider a two-dimensional interplanetary transfer of an ideal solar sailcraft. Before proceeding, we present some units widely used in astrodynamics for better numerical conditioning. The distance is in units of Astronomical Unit (AU, the mean distance between the Earth and the Sun), and time is in units of Time Unit (TU, a time unit such that an object in a circular orbit at 1 AU would have a speed of 1 AU/TU). For clearer presentation, a heliocentric-ecliptic inertial frame \((O,X,Y)\) and a polar reference frame \((O,r,\theta)\) are used, as shown in Fig. 1. The origin \(O\) is located at the Sun's center; the \(Y\) axis points towards the opposite direction of the vernal equinox, and the \(X\) axis is specified by rotating the \(Y\) axis 90 degrees clockwise in the ecliptic plane. \(r\in\mathbb{R}_{0}^{+}\) denotes the distance from the
origin to the solar sailcraft, and \(\theta\in[0,2\pi]\) is the rotation angle of the Sun-solar sailcraft line measured counterclockwise in the polar reference frame from the positive \(X\) axis.

In the context of preliminary mission design, it is reasonable to make some assumptions. The Sun and solar sailcraft are treated as point masses. The motion of the solar sailcraft is only influenced by the Sun's gravitational attraction and the propulsive acceleration due to the SRP acting on the solar sailcraft. Other disturbing factors, such as the solar wind and light aberration are neglected. Then, the dimensionless dynamics of the solar sailcraft is given by [11]

\[\begin{cases}\dot{r}(t)=u(t),\\ \dot{\theta}(t)=\frac{v(t)}{r(t)};\\ \dot{u}(t)=\frac{\beta\cos^{3}\alpha(t)}{r^{2}(t)}+\frac{v^{2}(t)}{r(t)}- \frac{1}{r^{2}(t)},\\ \dot{v}(t)=\frac{\beta\sin\alpha(t)\cos^{2}\alpha(t)}{r^{2}(t)}-\frac{u(t)v(t) }{r(t)},\end{cases} \tag{1}\]

where \(t\geq 0\) is the time; \(u\in\mathbb{R}\) and \(v\in\mathbb{R}\) is the radial and transverse speed in units of AU/TU, respectively; \(\beta\) is the normalized characteristic acceleration; \(\alpha\) represents the cone angle, which is the angle between the normal direction of sail plane \(\hat{n}\) and the direction of incoming rays \(\hat{i}_{u}\), measured positive in the clockwise direction from the normal direction of sail plane \(\hat{n}\).

### Attitude Constraints

Note that the radial component \(a_{u}\) and transverse component \(a_{v}\) of the ideal propulsive acceleration vector acting on the solar sailcraft are given by [41]

\[\begin{cases}a_{u}(t):=\frac{\beta\cos^{3}\alpha(t)}{r^{2}(t)},\\ a_{v}(t):=\frac{\beta\sin\alpha(t)\cos^{2}\alpha(t)}{r^{2}(t)}.\end{cases} \tag{2}\]

Figure 1: Dynamics of an ideal solar sailcraft.
Based on the assumption that the ideal sail is approximated with a flat and rigid reflecting surface, the attitude constraints can be rewritten in terms of the cone angle as

\[\alpha\in[-\phi_{max},\phi_{max}], \tag{3}\]

in which \(\phi_{max}\in(0,\frac{\pi}{2})\) is a given parameter depending on the minimum acceptable level of the electric power generated by the solar sailcraft [41]. It essentially sets the limits for the solar sailcraft's orientation with respect to the incoming solar rays.

The "force bubble" of the ideal solar mentioned in Ref. [49] due to attitude constraints is displayed in Fig. 2. The colormap represents the value of \(\beta\) in \([0.01892,0.3784]\) (corresponding to a characteristic acceleration domain of \([0.1,2]\) mm/s\({}^{2}\)). Then, for a given \(r\), it is clear that \(\beta\) defines the size of the force bubble, whereas \(\phi_{max}\) determines its actual shape, as shown by the pink-dashed lines. Specifically, the propulsive acceleration is constrained to the radial direction if \(\phi_{max}=0\).

### Formulation of the OCP

Without loss of generality, consider an initial condition for the solar sailcraft at the initial time \(t_{0}=0\) given by

\[r(0)=r_{0},\theta(0)=\theta_{0},u(0)=u_{0},v(0)=v_{0}. \tag{4}\]

The mission purpose is, through orientating the cone angle \(\alpha\) subject to constraints in Eq. (3), to drive the solar sailcraft governed by Eq. (1), into a coplanar circular orbit of radius \(r_{f}\) with the terminal condition given by

\[r(t_{f})=r_{f},u(t_{f})=0,v(t_{f})=\frac{1}{\sqrt{r_{f}}}. \tag{5}\]

Figure 2: Shape of the ideal sail force bubble.
The functional cost \(J\) is to minimize the arrival time (final time) \(t_{f}\), i.e.,

\[J=\int_{0}^{t_{f}}1\ dt. \tag{6}\]

## 3 Parameterization of Optimal Trajectories

In this section, we first derive the necessary conditions for optimal trajectories. Then, we formulate a parameterized system that enables the generation of an optimal trajectory via solving an initial value problem. For simplicity of presentation, important results and claims are written in lemmas with their proofs postponed to the appendix.

### Optimality Conditions

Denote by \(\mathbf{\lambda}=[\lambda_{r},\lambda_{\theta},\lambda_{u},\lambda_{v}]\) the co-state vector of the state vector \(\mathbf{x}=[r,\theta,u,v]\). Then, the Hamiltonian for the OCP is expressed as

\[\mathscr{H}=1+\lambda_{r}u+\lambda_{\theta}\frac{v}{r}+\lambda_{u}(\frac{ \beta\cos^{3}\alpha}{r^{2}}+\frac{v^{2}}{r}-\frac{1}{r^{2}})+\lambda_{v}(\frac {\beta\sin\alpha\cos^{2}\alpha}{r^{2}}-\frac{uv}{r}). \tag{7}\]

According to the PMP [50], we have

\[\begin{cases}\dot{\lambda}_{r}(t)=\lambda_{\theta}(t)\frac{v(t)}{r^{2}(t)}+ \lambda_{u}(t)[2\beta\frac{\cos^{3}\alpha(t)}{r^{3}(t)}+\frac{v^{2}(t)}{r^{2}( t)}-\frac{2}{r^{3}(t)}]+\lambda_{v}(t)[2\beta\cos^{2}\alpha(t)\frac{\sin \alpha(t)}{r^{3}(t)}-\frac{u(t)v(t)}{r^{2}(t)}],\\ \dot{\lambda}_{\theta}(t)=0,\\ \dot{\lambda}_{u}(t)=-\lambda_{r}(t)+\frac{\lambda_{v}(t)v(t)}{r(t)},\\ \dot{\lambda}_{v}(t)=-\frac{\lambda_{\theta}(t)}{r(t)}-2\frac{\lambda_{u}(t)v (t)}{r(t)}+\frac{\lambda_{v}(t)u(t)}{r(t)}.\end{cases} \tag{8}\]

Because \(t_{f}\) is free, it holds

\[\mathscr{H}(t)\equiv 0,\ \forall\ t\in[0,t_{f}]. \tag{9}\]

In addition, the terminal rotation angle \(\theta(t_{f})\) is not constrained, leading to

\[\lambda_{\theta}(t_{f})=0. \tag{10}\]

By the following lemma, we present the optimal guidance law.

**Lemma 1**: _The optimal guidance law that minimizes \(\mathscr{H}\) in Eq. (7) is_

\[\alpha^{*}=\text{Median}[-\phi_{max},\bar{\alpha},\phi_{max}],\text{where }\bar{\alpha}=\arctan\ \frac{-3\lambda_{u}-\sqrt{9\lambda_{u}^{2}+8\lambda_{v}^{2}}}{4\lambda_{v}}. \tag{11}\]

The proof of this lemma is postponed to A.
Notice that the optimal guidance law in Eq. (11) may not be smooth at some isolated points. In fact, the optimal guidance law in Eq. (11) can be approximated in a compact form as [51]

\[\alpha^{*}\approx\alpha^{*}(\delta)=\frac{1}{2}\ [\sqrt{\left(\bar{\alpha}+\phi_{ max}\right)^{2}+\delta}-\sqrt{\left(\bar{\alpha}-\phi_{max}\right)^{2}+\delta}]. \tag{12}\]

where \(\delta\) is a small non-negative number. It is clear that Eq. (12) is equivalent to Eq. (11) if \(\delta=0\). If \(\delta>0\) is sufficiently small, Eq. (12) acts like a smoothing filter function.

Note that the initial rotation angle \(\theta_{0}\) has no effect on the solutions because the transfer trajectory is rotatable [39]. For brevity, a triple \((r(t),u(t),v(t))\) for \(t\in[0,t_{f}]\) is said to be an optimal trajectory if all the necessary conditions in Eqs. (8), (9), (10), and (11) are met. In order to generate the optimal guidance command in real time via NNs, a training dataset containing a large number of optimal trajectories is required. In this regard, one viable approach is to employ some root-finding algorithms to solve the TPBVP

\[[r(t_{f})-r_{f},u(t_{f}),v(t_{f})-\frac{1}{\sqrt{r_{f}}},\mathscr{H}(t_{f})]= \mathbf{0}. \tag{13}\]

Nevertheless, this procedure is usually time-consuming and suffers from convergence issues. In the next subsection, we present a parameterized system so that an optimal trajectory can be readily obtained by solving an initial value problem instead.

### Parameterized System

Define a new independent variable \(\tau\) as below

\[\tau=t_{f}-t,t\in[0,t_{f}]. \tag{14}\]

Let us establish a first-order ordinary differential system

\[\begin{cases}\dot{R}(\tau)=-U(\tau),\\ \dot{U}(\tau)=-\frac{\beta\cos^{3}\hat{\alpha}(\tau)}{R^{2}(\tau)}-\frac{V^{2} (\tau)}{R(\tau)}+\frac{1}{R^{2}(\tau)},\\ \dot{V}(\tau)=-\frac{\beta\sin\hat{\alpha}(\tau)\cos^{2}\hat{\alpha}(\tau)}{R ^{2}(\tau)}+\frac{U(\tau)V(\tau)}{R(\tau)},\\ \dot{\lambda}_{R}(\tau)=-\lambda_{U}(\tau)[2\beta\frac{\cos^{3}\hat{\alpha}( \tau)}{R^{3}(\tau)}\frac{V^{2}(\tau)}{R^{2}(\tau)}-\frac{2}{R^{3}(\tau)}]- \lambda_{V}(\tau)[2\beta\cos^{2}\hat{\alpha}(\tau)\frac{\sin\hat{\alpha}(\tau )}{R^{3}(\tau)}-\frac{U(\tau)V(\tau)}{R^{2}(\tau)}],\\ \dot{\lambda}_{U}(\tau)=\lambda_{R}(\tau)-\frac{\lambda_{V}(\tau)V(\tau)}{R( \tau)},\\ \dot{\lambda}_{V}(\tau)=2\frac{\lambda_{U}(\tau)V(\tau)}{R(\tau)}-\frac{ \lambda_{V}(\tau)U(\tau)}{R(\tau)},\end{cases} \tag{15}\]

where \((R,U,V,\lambda_{R},\lambda_{U},\lambda_{V})\in\mathbb{R}_{0}^{+}\times\mathbb{ R}^{5}\), and \(\hat{\alpha}\) is defined as

\[\hat{\alpha}=\frac{1}{2}\ [\sqrt{\left(\bar{\alpha}+\phi_{max}\right)^{2}+ \delta}-\sqrt{\left(\bar{\alpha}-\phi_{max}\right)^{2}+\delta}]\ \text{with}\ \bar{\alpha}=\arctan\ \frac{-3\lambda_{\text{U}}-\sqrt{9\lambda_{\text{U}}^{2}+8\lambda_{\text{V} }^{2}}}{4\lambda_{\text{V}}}. \tag{16}\]
The initial value at \(\tau=0\) for the system in Eq. (15) is set as

\[R(0)=R_{0},U(0)=0,V(0)=\frac{1}{\sqrt{R_{0}}},\lambda_{R}(0)=\lambda_{R_{0}}, \lambda_{U}(0)=\lambda_{U_{0}},\lambda_{V}(0)=\lambda_{V_{0}}. \tag{17}\]

The value of \(\lambda_{V_{0}}\) satisfies the following equation

\[\begin{split} 1+\lambda_{R}(0)U(0)+\lambda_{U}(0)(\frac{\beta\cos^{3} \hat{\alpha}(0)}{R^{2}(0)}+\frac{V^{2}(0)}{R(0)}-\frac{1}{R^{2}(0)})\\ +\lambda_{V}(0)(\frac{\beta\sin\hat{\alpha}(0)\cos^{2}\hat{ \alpha}(0)}{R^{2}(0)}-\frac{U(0)V(0)}{R(0)})=0.\end{split} \tag{18}\]

Substituting Eq. (17) into Eq. (18) leads to

\[1+\lambda_{U_{0}}\frac{\beta\cos^{3}\hat{\alpha}(0)}{R_{0}^{2}}+\lambda_{V_{0} }\frac{\beta\sin\hat{\alpha}(0)\cos^{2}\hat{\alpha}(0)}{R_{0}^{2}}=0. \tag{19}\]

In view of Eqs. (16) and (19), it is clear that for a given \(\lambda_{U_{0}}\), the value for \(\lambda_{V_{0}}\) can be numerically determined by solving Eq. (19). By the following lemma, we shall show that an optimal trajectory can be generated by solving an initial value problem based on the parameterized system in Eq. (15).

Lemma 2: For any given constants \(\lambda_{R_{0}}\) and \(\lambda_{U_{0}}\), a fixed final time \(t_{f}\), and a parameter \(\lambda_{V_{0}}\) to be determined by solving Eq. (19), denote by

\[\begin{split}\mathcal{F}:=(R(\tau,\lambda_{R_{0}},\lambda_{U_{0} }),U(\tau,\lambda_{R_{0}},\lambda_{U_{0}}),V(\tau,\lambda_{R_{0}},\lambda_{U_ {0}}),\\ \lambda_{R}(\tau,\lambda_{R_{0}},\lambda_{U_{0}}),\lambda_{U}( \tau,\lambda_{R_{0}},\lambda_{U_{0}}),\lambda_{V}(\tau,\lambda_{R_{0}}, \lambda_{U_{0}}))\in\mathbb{R}^{+}\times\mathbb{R}^{5}.\end{split}\]

the solution of the \((\lambda_{R_{0}},\lambda_{U_{0}})\)-parameterized system in Eq. (15) with the initial value specified in Eq. (17). Define \(\mathcal{F}_{p}\) as

\[\begin{split}&\mathcal{F}_{p}:=\{R(\tau,\lambda_{R_{0}},\lambda_{U_ {0}}),U(\tau,\lambda_{R_{0}},\lambda_{U_{0}}),V(\tau,\lambda_{R_{0}},\lambda_{U _{0}}),\tau|\\ &(R(\tau,\lambda_{R_{0}},\lambda_{U_{0}}),U(\tau,\lambda_{R_{0}},\lambda_{U_{0}}),V(\tau,\lambda_{R_{0}},\lambda_{U_{0}}),\tau)\in\mathcal{F} \}.\end{split}\]

Then \(\mathcal{F}_{p}\) represents the solution space of an optimal trajectory starting from a circular orbit with a radius of \(R_{0}\).

The proof of this lemma is postponed to A.

With the parameterized system established, obtaining an optimal trajectory becomes a straightforward process that involves solving an initial value problem, rather than tackling a TPBVP. Once the optimal trajectory \(\mathcal{F}_{p}\) is obtained, the corresponding optimal guidance command \(\hat{\alpha}\) can be determined through the definition of \(\mathcal{F}\) and reference to Eq. (16). To maintain clarity and prevent any ambiguity in notation, we continue to employ \((r,u,v)\) to represent the flight state and we use \(\tau\) to represent the optimal time of flight along the optimal trajectory \(\mathcal{F}_{p}\). Since the optimal time of flight is of great importance for preliminary
mission design, evaluating the optimal time of flight without the need for the exact guidance law is quite attractive, as was done in Refs. [31; 35]. In this paper, we include \(\tau\) into the flight state, and it can be evaluated via a trained NN, as shown by our scheme in Subsection 4.3. Likewise, the optimal guidance command is again denoted by \(\alpha\). Then we shall define \(f\) as the nonlinear mapping from the flight state \((r,u,v,\tau)\) to the optimal guidance command \(\alpha\), i.e.,

\[f:=(r,u,v,\tau)\rightarrow\alpha.\]

According to the universal approximation theorem [52], employing a large number of sampled data that map the flight state \((r,u,v,\tau)\) to the optimal guidance command \(\alpha\) to train feedforward NNs enables the trained networks to accurately approximate the nonlinear mapping \(f\). Before generating the dataset, a _nominal trajectory_ is typically required. Then, by applying perturbations over \(\lambda_{R_{0}}\) and \(\lambda_{U_{0}}\) on the _nominal trajectory,_ lots of optimal trajectories containing optimal state-guidance command pairs can be aggregated. With this dataset, the optimal guidance command can be generated in real time within the supervised learning framework, as we will detail in the next section.

## 4 Generation of Optimal Guidance Commands in Real Time via NNs

In this section, we begin by outlining the procedure for generating the training dataset. As our analysis progresses, we observe a notable characteristic of the optimal guidance command, that is, it may abruptly change from one extreme to another, resulting in discontinuous jumps. This behavior poses a challenge, as the approximation performance of the NN typically deteriorates under such circumstances, a phenomenon seen in Refs. [43; 44; 45]. To address this issue, we introduce a preprocessing method that employs two co-states to detect the discontinuous jump. With the saturation function defined in Eq. (12), the pre-processed guidance command becomes smooth everywhere. Following this procedure, we present a scheme for generating optimal guidance commands in real time via NNs.

### Generation of the Training Dataset

#### 4.1.1 Nominal Trajectory

Consider a solar sailcraft subject to attitude constraints \(\alpha\in[-\frac{\pi}{3},\frac{\pi}{3}]\), starting from Earth's orbit around the Sun, which is a circular orbit with radius of 1 AU (1 AU = \(1.4959965\times 10^{11}\) m), i.e.,

\[r(0)=1,\theta(0)=0,u(0)=0,v(0)=1, \tag{20}\]

and the mission is to fly into Mars' orbit around the Sun, which is a circular orbit with radius of 1.524 AU, i.e.,

\[r(t_{f})=1.524,u(t_{f})=0,v(t_{f})=\frac{1}{\sqrt{1.524}}. \tag{21}\]
The normalized characteristic acceleration \(\beta\) is set as a constant of 0.16892, corresponding to a characteristic acceleration of 1 mm/s\({}^{2}\). To demonstrate the effects of the saturation function in Eq. (12), the optimal guidance laws in Eqs. (11) and (12) are embedded into the shooting function defined by Eq. (13). To enhance the convergence of the indirect method, we initially consider cases without attitude constraints and subsequently apply a homotopy technique to accommodate attitude constraints.

Fig. 3 shows the optimal guidance command profiles for the unconstrained and constrained cases with two smoothing constants. It is evident that the constrained cases contain rapid changes in the guidance command. Additionally, the guidance command with \(\delta=1\times 10^{-3}\) displays smoother behavior compared to the case with \(\delta=0\). Besides, the arrival time for the unconstrained case is 7.01204 TU ( 1 TU \(=5.0226757\times 10^{6}\) s), whereas this time span extends to 7.01818 TU for the constrained case with \(\delta=0\). Regarding \(\delta=1\times 10^{-3}\), the arrival time is 7.01838 TU, indicating that a smoothing constant of \(\delta=1\times 10^{-3}\) leads to a minor functional penalty of merely \(2\times 10^{-4}\) TU. Consequently, we designate the transfer trajectory with \(\delta=1\times 10^{-3}\) as the _nominal trajectory_ henceforth.

#### 4.1.2 Dataset

Denote by \(\lambda_{r_{f}}^{*}\), \(\lambda_{u_{f}}^{*}\) the final values along the _nominal trajectory_ for the co-state \(\lambda_{r}\) and \(\lambda_{u}\), respectively. They are computed as \(\lambda_{r_{f}}^{*}=-21.51\) and \(\lambda_{u_{f}}^{*}=8.54\). Consider a new set of values given by

\[\lambda_{R_{0}}^{\prime}=\lambda_{r_{f}}^{*}+\delta\lambda_{R}, \lambda_{U_{0}}^{\prime}=\lambda_{u_{f}}^{*}+\delta\lambda_{U}, \tag{22}\]

where the perturbations \(\delta\lambda_{R}\) and \(\delta\lambda_{U}\) are chosen such that the resulting \(\lambda_{R_{0}}^{\prime}\) and \(\lambda_{U_{0}}^{\prime}\) are uniformly distributed in the interval \([-23,-20]\) and \([5,11]\), respectively. Set \(\lambda_{R_{0}}=\lambda_{R_{0}}^{\prime}\),\(\lambda_{U_{0}}=\lambda_{U_{0}}^{\prime}\), and calculate the value of \(\lambda_{V_{0}}\) by solving Eq. (19). Then, \(\mathcal{F}\) can be obtained by propagating the \((\lambda_{R_{0}},\lambda_{U_{0}})\)-parameterized system with the initial value specified in Eq. (17). Define an empty set \(\mathcal{D}\), and insert \(\mathcal{F}_{p}\) depending on the pair \((\lambda_{R_{0}},\lambda_{U_{0}})\), into \(\mathcal{D}\) until the

Figure 3: Optimal guidance command profiles for the unconstrained and constrained cases.
perturbation process described in Eq. (22) ends. In this context, we can amass a dataset \(\mathcal{D}\) comprising the essential optimal state-guidance command pairs required for NN training.

Because we are specifically dealing with orbit transfer missions from Earth to Mars in this study, any optimal trajectory with \(R(\tau)>1.54\) AU for \(\tau\in[0,t_{f}]\) (\(t_{f}\) is fixed as \(10\)) is classified as beyond the region of interest and excluded from \(\mathcal{D}\). Moreover, to explore the NN's generalization ability, the _nominal trajectory_ is excluded from \(\mathcal{D}\), as was done in Ref. [35]. Ultimately, a total of \(18,324\) optimal trajectories are acquired, and they are evenly sampled, resulting in \(6,416,613\) pairs.

### Preprocessing Discontinuous Jumps

While generating the dataset, we encounter instances where the optimal guidance command displays discontinuous jumps. Remember that these jumps can be detected using \(\lambda_{u}\) and \(\lambda_{v}\), as elaborated in A. For a representative optimal trajectory defined in \([0,t_{f}]\), spotting a discontinuous jump point \(t_{j}\) is straightforward, as illustrated in Fig. 4.

The optimal trajectory is then split into two segments, i.e., \([0,t_{j}]\) and \([t_{j},t_{f}]\). Assume that there is only one jump in this scenario, we can derive a preprocessed smooth guidance command, denoted by \(\alpha_{pre}\), using the following approach

\[\alpha_{pre}=\begin{cases}-\alpha,\text{for }t\in[0,t_{j}]\\ \phantom{-}\alpha,\text{for }t\in[t_{j},t_{f}]\end{cases} \tag{23}\]

A representative example is provided in Fig. 5. It is evident that, thanks to the saturation function in Eq. (12) and the preprocessing procedure outlined in Eq. (23), \(\alpha_{pre}\) now demonstrates smooth behavior along the entire optimal trajectory. Conversely, if the discontinuous jump does not appear along the optimal trajectory, it is clear that \(\alpha_{pre}=\alpha\) always holds for \(t\in[0,t_{f}]\). Moreover, we observe that multiple discontinuous jumps are rare and primarily occur in orbital transfers with unusually long time of flight. Hence, in this study, we limit our consideration to orbital transfers with a maximum of one discontinuous jump.

Figure 4: Two co-state profiles along an optimal trajectory.
With the parameterized system, a sample of dataset \(\mathcal{D}_{\alpha}:=\{(r,u,v,\tau),\alpha\}\) can be extracted from the dataset \(\mathcal{D}\). Thanks to the preprocessing procedure aforementioned, we can transform the dataset \(\mathcal{D}_{\alpha}\), which could contain discontinuous jumps, into a refined dataset \(\mathcal{D}_{\alpha_{pre}}:=\{(r,u,v,\tau),\alpha_{pre}\}\) exclusively consisting of smooth guidance commands. Consequently, a well-trained NN \(\mathcal{N}_{\alpha_{pre}}\) based on \(\mathcal{D}_{\alpha_{pre}}\) is anticipated to yield reduced approximation errors compared to an NN \(\mathcal{N}_{\alpha}\) trained on \(\mathcal{D}_{\alpha}\)[25].

For the sake of completeness, we will elucidate the process of reverting the output \(\alpha_{pre}\) from \(\mathcal{N}_{\alpha_{pre}}\) back into the original optimal guidance command \(\alpha\) for practical implementation. It can be demonstrated that the sign of the original optimal guidance command \(\alpha\) is always opposite to that of \(\lambda_{v}\), i.e.,

\[4\lambda_{v}\cdot\arctan\alpha=-3\lambda_{u}-\sqrt{9\lambda_{u}^{2}+8\lambda_ {v}^{2}}<0. \tag{24}\]

Hence, the original optimal guidance command \(\alpha\) can be obtained by simply comparing the sign of the preprocessed guidance command \(\alpha_{pre}\) and \(\lambda_{v}\), i.e.,

\[\alpha=\begin{cases}-\alpha_{pre},\text{if }\alpha_{pre}\cdot\lambda_{v}>0\\ \alpha_{pre},\ \text{ otherwise}\end{cases} \tag{25}\]

### Scheme for Generating Optimal Guidance Commands in Real Time

To facilitate real-time implementation, a system composed of three NNs is established, as depicted in Fig. 6. To elaborate, \(\mathcal{D}\) is divided into three samples, i.e., \(\mathcal{D}_{\tau}:=\{(r,u,v),\tau\}\), \(\mathcal{D}_{\alpha_{pre}}:=\{(r,u,v,\tau),\alpha_{pre}\}\), and \(\mathcal{D}_{\lambda_{v}}:=\{(r,u,v,\tau),\lambda_{v}\}\); and the corresponding NN is denoted by \(\mathcal{N}_{\tau}\), \(\mathcal{N}_{\alpha_{pre}}\), and \(\mathcal{N}_{\lambda_{v}}\), respectively. \(\mathcal{N}_{\tau}\) is designed to forecast the optimal time of flight \(\tau\) given a flight state \((r,u,v)\). Additionally, this network facilitates preliminary mission design by providing time of flight evaluations without necessitating an exact solution. The output of \(\mathcal{N}_{\tau}\) is subsequently used as part of the input for \(\mathcal{N}_{\alpha_{pre}}\), which yields the preprocessed guidance command \(\alpha_{pre}\). Concurrently, \(\mathcal{N}_{\lambda_{v}}\) predicts the co-state \(\lambda_{v}\), in which its sign

Figure 5: Preprocessed and original optimal guidance command profiles along an optimal trajectory.
\(sgn(\lambda_{v})\) is used to revert the preprocessed guidance command \(\alpha_{pre}\) back into the original guidance command \(\alpha\). Once these three NNs are trained offline, they enable the generation of the optimal guidance command \(\alpha\) given an initial condition \((r_{0},u_{0},v_{0})\). Consequently, the trained NNs offer closed-form solutions to the OCP. This endows the proposed method with robustness and generalization abilities, as shown in Subsections 5.2 and 5.3.

### NN Training

Now, our focus shifts to the implementation of NN training algorithm. In addition to the three aforementioned NNs, i.e., \(\mathcal{N}_{\tau}\), \(\mathcal{N}_{\lambda_{v}}\), and \(\mathcal{N}_{\alpha_{pre}}\), we also include the training of \(\mathcal{N}_{\alpha}\) based on the sample \(\mathcal{D}_{\alpha}=\{(r,u,v,\tau),\alpha\}\). This is done to highlight the enhancement in approximation resulting from the preprocessing procedure. All the networks considered are feedforward NNs with multiple hidden layers. _Prior_ to training, the dataset samples are shuffled to establish a split of 70% for training, 15% for validation, and 15% for testing sets. All input and output data are normalized by the min-max scaling and scaled to the range [0, 1]. Selecting an appropriate NN structure, particularly the number of hidden layers and neurons per layer, is a non-trivial task. A structure that is too simplistic tends to lead to underfitting. Conversely, overfitting often arises when the structure is overly complex. Making a balance between time consumption and the overfitting issue, we adopt a structure with three hidden layers, each containing 20 neurons. Subsequently, the sigmoid function serves as the activation function. For the output layer, a linear function is utilized. The crux of the training lies in minimizing the loss function, quantified as the mean squared error (MSE) between the predicted values from the trained NNs and the actual values within the dataset samples. We employ the 'Levenberg-Marquardt' in Ref. [53] for training NNs, and the training is terminated after 1,000 epochs or when the loss function drops below \(1\times 10^{-8}\). To prevent overfitting, we employ an early stopping criteria based on the concept of patience. Other hyperparameters utilized in training are set to their default values.

Figure 6: Generation of optimal guidance commands in real time via NNs.
Fig. 7 illustrates the training progression of four NNs. As seen in Fig. 6(a), it is evident that the MSEs for the training, validation, and test sets decrease to \(1\times 10^{-8}\) in less than 300 epochs. Table 1 displays the MSEs of the four NNs upon completion of the training process. Notably, the validation and test errors of \(\mathcal{N}_{\alpha_{pre}}\) are smaller than those of \(\mathcal{N}_{\alpha}\), indicating that our preprocessing procedure enhances the approximation accuracy. This is further confirmed through simulations detailed in the next section.

## 5 Numerical Simulations

In this section, we present some numerical simulations to demonstrate the efficacy and performance of the proposed methodology. We commence by showing its optimality performance. Subsequently, robustness against perturbations and generalization ability of the proposed method are investigated. In addition, we evaluate the real-time execution of the

\begin{table}
\begin{tabular}{c c c c c} \hline  & \(\mathcal{N}_{\tau}\) & \(\mathcal{N}_{\lambda_{v}}\) & \(\mathcal{N}_{\alpha_{pre}}\) & \(\mathcal{N}_{\alpha}\) \\ \hline Training & \(1\times 10^{-8}\) & \(4.94\times 10^{-6}\) & \(7.53\times 10^{-7}\) & \(1.72\times 10^{-6}\) \\ Validation & \(1\times 10^{-8}\) & \(5.02\times 10^{-6}\) & \(6.82\times 10^{-7}\) & \(3.53\times 10^{-6}\) \\ Test & \(1\times 10^{-8}\) & \(5.07\times 10^{-6}\) & \(7.70\times 10^{-7}\) & \(2.55\times 10^{-6}\) \\ \hline \end{tabular}
\end{table}
Table 1: MSEs of four NNs

Figure 7: Training histories of four NNs.
guidance command generation. It is worth noting that the proposed method has been implemented on a desktop equipped with an Intel Core i9-10980XE CPU @3.00 GHz and 128 GB of RAM.

### Optimality Performance

This subsection is devoted to examining whether or not the proposed approach can produce optimal solutions compared with existing methods. Table 2 outlines the initial conditions for two orbital transfers. For comparison, two strategies are employed to steer the solar sailcraft. The first strategy employs three NNs, i.e., \(\mathcal{N}_{\tau}\), \(\mathcal{N}_{\alpha_{pre}}\), and \(\mathcal{N}_{\lambda_{v}}\), while the second one uses two NNs, \(\mathcal{N}_{\tau}\) and \(\mathcal{N}_{\alpha}\). The indirect method, resolving the shooting function in Eq. (13) will be used for comparison.

For Case 1, the solar sailcraft embarks on a journey from Earth's orbit around the Sun, and the mission is to steer the solar sailcraft into Mars' orbit around the Sun. Fig. 8 represents the outcomes yielded by the trained NNs and the indirect method. The optimal time of flight obtained via the indirect method is 7.0184 TU, while the prediction from \(\mathcal{N}_{\tau}\) is 7.0192 TU for the given initial condition. Remarkably, \(\mathcal{N}_{\tau}\) exhibits precise approximation of the optimal time of flight throughout the entire transfer, as shown by Fig. 8a. On the other hand, Fig. 8b shows less accurate approximation of the optimal co-state \(\lambda_{v}\), especially at the end of the transfer. Fortunately, the sign of the predicted \(\lambda_{v}\), rather than its exact value, is sufficient for reverting the preprocessed guidance command back to the original guidance command. This observation is depicted in Fig. 8c, where the guidance commands derived from the first strategy closely align with the indirect method, except at the end of the transfer. Since the chattering guidance command derived from the proposed method at the end of the transfer is not practically implementable, the numerical simulation is terminated once the predicted time of flight \(\tau\) from \(\mathcal{N}_{\tau}\) drops below 0.005 TU (equivalent to 0.2907 days) hereafter. In contrast, the guidance command derived from the second strategy exhibits lower accuracy, especially in the presence of rapid changes in the optimal guidance command. The degeneration of approximation performance, caused by the second strategy, is demonstrated in Fig. 8d. It is worth emphasizing that this problem of performance degeneration may not be well fixed by simply adjusting the hyperparameters for training NNs, as shown by the numerical simulations in Ref. [43]. In that work, a random search for the hyperparameters was employed, and the trained NNs, even with the best approximation performance, were not able to capture the rapid change in the guidance command. In contrast, thanks to the preprocessing procedure in this paper, the well-trained NNs are capable of precisely generating the optimal guidance command.

\begin{table}
\begin{tabular}{c c c c c} \hline \hline  & \(r_{0}\) & \(\theta_{0}\) & \(u_{0}\) & \(v_{0}\) \\ \hline Case 1 & 1 & 0 & 0 & 1 \\ Case 2 & 1.05 & 0 & 0.15 & 1.03 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Initial conditions for two orbital transfer cases
To assess the performance of fulfilling the terminal condition, we define the flight error \(\Delta\Phi\) as the difference between the NN-driven terminal flight state \((r^{\mathcal{N}}(t_{f}),u^{\mathcal{N}}(t_{f}),v^{\mathcal{N}}(t_{f}))\) and desired terminal flight state \((1.524,0,\frac{1}{\sqrt{1.524}})\). This error is computed as follows

\[\Delta\Phi=(|r^{\mathcal{N}}(t_{f})-1.524|,|u^{\mathcal{N}}(t_{f})|,|v^{ \mathcal{N}}(t_{f})-\frac{1}{\sqrt{1.524}}|).\]

Then, the flight error caused by the first strategy is \(\Delta\Phi=(5.3703\times 10^{-7},3.0121\times 10^{-5},3.6943\times 10^{-5})\), which is lower than that of the second strategy with \(\Delta\Phi=(2.2220\times 10^{-5},4.1852\times 10^{-4},2.5968\times 10^{-4})\). Although both strategies can guide the solar sailcraft into Mars' orbit with negligible flight errors, the guidance command derived from the first strategy is more feasible and accurate than the second strategy, as shown in Figs. (c)c and (d)d. In addition, the transfer trajectories derived from the first strategy and the indirect method are illustrated in Fig. 9.

Given the superiority of the first strategy over the second, our focus will exclusively be on the first strategy, which will be called NN-based approach hereafter. Shifting to Case 2, the output of \(\mathcal{N}_{\tau}\) for the initial condition stands at 6.1643 TU, while the optimal solution attained via the indirect method is 6.1552 TU, resulting in a marginal functional penalty of merely 0.5290 days. Furthermore, Fig. 10 illustrates that the NN-based approach precisely captures the discontinuous jump and effectively approximates the optimal guidance command, except at the end of the transfer. The transfer trajectories derived from the NN-based approach

Figure 8: Comparison of results from the trained NNs and the indirect method for Case 1.
and the indirect method are depicted in Fig. 11. As a consequence, the flight error caused by the NN-based approach is \(\Delta\Phi=(5.2011\times 10^{-6},3.3526\times 10^{-4},1.5393\times 10^{-4})\), implying that the NN-based approach capably steers the solar sailcraft into Mars' orbit with an acceptable flight error.

### Robustness Analysis

In reality, the normalized characteristic acceleration \(\beta\) is influenced by diverse factors, such as solar sailing efficiency and solar radiation pressure [54]. Consequently, the solar sailcraft must continually adjust its flight trajectory based on its current flight state to rectify any flight errors. Unfortunately, this process is quite challenging for onboard systems with limited computational resources. Additionally, solving OCPs featuring parameters under perturbations tends to be problematic for conventional direct and indirect methods

Figure 10: Guidance command profiles derived from the NN-based approach and the indirect method for Case 2.

Figure 9: Transfer trajectories derived from the first strategy and the indirect method for Case 1.
that rely on precise system models. For this reason, we evaluate the proposed method's robustness against perturbations concerning \(\beta\).

The initial condition for the solar sailcraft is set as

\[r_{0}=1.1,\theta_{0}=\frac{\pi}{2},u_{0}=0.18,v_{0}=0.93,\]

and the mission is to fly into Mars' orbit. The normalized characteristic acceleration \(\beta\) is under perturbations following a standard uniform distribution within the range of \((-15\%,15\%)\).

Figs.12 and 13 depicts the transfer trajectory and the corresponding guidance command derived from the NN-based approach, respectively.

It can be seen that the colormap-identified speed exhibits a gradual reduction for most of the transfer, followed by a gradual increase towards the end of the transfer, as depicted in Fig.12. Reg

Figure 11: Transfer trajectories derived from the NN-based approach and the indirect method for Case 2.

Figure 12: NN-based transfer trajectory with \(\beta\) under perturbations.
ance command is generally smooth during the initial phase of the transfer. Subsequently, to uphold precise flight, the NN-based approach acts like an error corrector by dynamically adapting the guidance command in the latter phase of the transfer. A remarkable observation is that even when the solar sailcraft's dynamics is under perturbations, the attitude constraints remain unviolated throughout the entire transfer. Furthermore, it takes 5.9152 TU for the solar sailcraft to accomplish its journey with a minor flight error of \(\Delta\Phi=(2.0160\times 10^{-5},9.8708\times 10^{-5},8.4738\times 10^{-5})\).

### Generalization Analysis

The assessment of generalization ability pertains to evaluating the predictive performance for inputs that lie outside the training dataset. Define a state space \(\mathcal{A}\), that is not inside \(\mathcal{D}\), as below

\[r_{0}\in[1,1.15],\theta_{0}\in[0,2\pi],u_{0}\in[0,0.1],v_{0}\in[0.8,1.2].\]

Within this space, the solar sailcraft's initial conditions are randomly chosen. Subsequently, the proposed approach is employed to guide the solar sailcraft from the selected initial condition to Mars' orbit. A total of 30 tests are conducted, resulting in 20 successful transfers, as depicted in Fig. 14, in which the small circles denote the solar sailcraft's initial positions. Among these successful cases, the average flight error is \(\Delta\Phi=(4.5091\times 10^{-5},1.4726\times 10^{-4},1.0147\times 10^{-4})\). This indicates that the proposed method is able to steer the spacecraft to Mars' obit accurately even with some initial states that are not in the training set.

### Real-Time Performance

Since the output of a feedforward NN is essentially a composition of linear mappings of the input vector, a trained NN can produce the output within a constant time. Recall that our method requires three trained NNs with simple structures. 10,000 trials of the proposed

Figure 13: NN-based guidance command profile with \(\beta\) under perturbations.
method across various flight states are run in a C-based computational environment, and the mean execution time for generating a guidance command is 0.0177 ms. This translates to approximately 0.5310 ms on a typical flight processor operating at 100 MHz [55]. Conversely, indirect methods typically necessitate time-consuming processes of integration and iteration. We utilize _fsolve_ function to solve the relevant TPBVP and observe a convergence time of approximately 1.22 seconds, even when an appropriate initial guess is provided. Consequently, the computational burden associated with optimization-based indirect methods can be substantially alleviated through adopting the proposed method.

## 6 Conclusions

Real-time optimal control for attitude-constrained solar sailcrafts in coplanar circular-to-circular interplanetary transfers was studied in the paper. The necessary conditions governing optimal trajectories were deduced through the PMP. Different from the conventional methods that rely on optimization-based solvers to collect a training dataset, we formulated a parameterized system. This enabled us to generate an optimal trajectory by simply solving an initial value problem, and the optimal state-guidance command pairs could be readily obtained. Regarding the challenge posed by potential discontinuities in the optimal guidance command, which could undermine the NN's approximation performance, we developed a technique to preprocess the guidance command. With the saturation function introduced in the optimal guidance law, this preprocessing technique smoothened the guidance command. In this way, the approximation performance of the preprocessed guidance command was enhanced, and another NN was used to predict one co-state, whose sign was adopted to revert the preprocessed guidance command into the potentially discontinuous original guidance command. As a result, the well-trained NNs were capable of not only generating the optimal guidance command precisely in real time but also predicting the optimal time of flight given a flight state.

Figure 14: NN-based transfer trajectories with initial conditions randomly chosen in \(\mathcal{A}\).
## Acknowledgement

This research was supported by the National Natural Science Foundation of China under Grant No. 62088101.

## Appendix A Proof of Lemmas in Section 2

Proof of Lemma 1. In view of \(\dot{\lambda}_{\theta}(t)=0\) in Eqs. (8) and (10), it is easy to see that \(\lambda_{\theta}(t)\) remains zero for \(t\in[0,t_{f}]\) along an optimal trajectory. Then, Eq. (7) reduces to

\[\mathscr{H}=1+\lambda_{r}u+\lambda_{u}(\frac{\beta\cos^{3}\alpha}{r^{2}}+\frac {v^{2}}{r}-\frac{1}{r^{2}})+\lambda_{v}(\frac{\beta\sin\alpha\cos^{2}\alpha}{r ^{2}}-\frac{uv}{r}). \tag{10}\]

To obtain stationary values of \(\mathscr{H}\), we have

\[\frac{\partial\mathscr{H}}{\partial\alpha}=\frac{\beta}{r^{2}}\cos^{3}\alpha [-3\lambda_{u}\tan\alpha+\lambda_{v}(1-2\tan^{2}\alpha)]=0. \tag{11}\]

Assume that \(\lambda_{u}\) and \(\lambda_{v}\) are not both zero [10]. Because \(\beta>0\), and \(\cos\alpha\neq 0\) for \(\alpha\in[-\phi_{max},\phi_{max}]\), two different roots, denoted by \(\alpha_{1}\) and \(\alpha_{2}\), can be obtained by solving Eq. (11) as

\[\begin{cases}\alpha_{1}=\arctan\ \frac{-3\lambda_{u}+\sqrt{9\lambda_{u}^{2}+8 \lambda_{v}^{2}}}{4\lambda_{v}},\\ \alpha_{2}=\arctan\ \frac{-3\lambda_{u}-\sqrt{9\lambda_{u}^{2}+8\lambda_{v}^{2}}}{4 \lambda_{v}}.\end{cases} \tag{12}\]

Now we proceed to check for the Hessian's positive definiteness. According to Eq. (11), taking the second partial derivative of \(\mathscr{H}\) w.r.t. \(\alpha\) yields

\[\frac{\partial^{2}\mathscr{H}}{\partial\alpha^{2}}=\frac{\beta}{r^{2}}\{-3 \cos^{2}\alpha\sin\alpha[-3\lambda_{u}\tan\alpha+\lambda_{v}(1-2\tan^{2} \alpha)]+\cos^{3}\alpha[-3\lambda_{u}\sec^{2}\alpha-4\lambda_{v}\tan\alpha \sec^{2}\alpha]\}. \tag{13}\]

Substituting Eq. (12) into Eq. (13) leads to

\[\begin{cases}\frac{\partial^{2}\mathscr{H}}{\partial\alpha^{2}}=\frac{\beta }{r^{2}}\cos\alpha(-\sqrt{9\lambda_{u}^{2}+8\lambda_{v}^{2}}),\\ \frac{\partial^{2}\mathscr{H}}{\partial\alpha_{2}^{2}}=\frac{\beta}{r^{2}} \cos\alpha(+\sqrt{9\lambda_{u}^{2}+8\lambda_{v}^{2}}).\end{cases} \tag{14}\]

Clearly, \(\frac{\partial^{2}\mathscr{H}}{\partial\alpha^{2}}>0\) always holds if and only if \(\alpha=\alpha_{2}\). Therefore, the local minimum solution \(\bar{\alpha}\) is given by

\[\bar{\alpha}=\alpha_{2}=\arctan\ \frac{-3\lambda_{u}-\sqrt{9\lambda_{u}^{2}+8 \lambda_{v}^{2}}}{4\lambda_{v}}. \tag{15}\]

Notice that Eq. (15) will become singular if \(\lambda_{v}(t)=0\). Then \(\lambda_{v}(t)=0\) can hold at some isolated points for \(t\in[0,t_{f}]\) or \(\lambda_{v}(t)\equiv 0\) holds in a time interval \([t_{a},t_{b}]\in[0,t_{f}]\). We now
analyze the first case.

If \(\lambda_{u}(t)<0\) for \(t\in[t_{a},t_{b}]\), then

\[\begin{split}\lim_{\lambda_{u}<0,\lambda_{v}\to 0}\frac{-3 \lambda_{u}-\sqrt{9\lambda_{u}^{2}+8\lambda_{v}^{2}}}{4\lambda_{v}}=\lim_{ \lambda_{u}<0,\lambda_{v}\to 0}\frac{-3\lambda_{u}+3\lambda_{u}\sqrt{1+\frac{8 \lambda_{v}^{2}}{9\lambda_{u}^{2}}}}{4\lambda_{v}}\\ \approx\lim_{\lambda_{u}<0,\lambda_{v}\to 0}\frac{-3 \lambda_{u}+3\lambda_{u}(1+\frac{1}{2}\frac{8\lambda_{v}^{2}}{9\lambda_{u}^{2} })}{4\lambda_{v}}=\lim_{\lambda_{u}<0,\lambda_{v}\to 0}\frac{\lambda_{v}}{3 \lambda_{u}}=0.\end{split}\] (A7)

Analogously, if \(\lambda_{u}(t)>0\) for \(t\in[t_{a},t_{b}]\), we have

\[\begin{split}\lim_{\lambda_{u}>0,\lambda_{v}\to 0}\frac{-3 \lambda_{u}-\sqrt{9\lambda_{u}^{2}+8\lambda_{v}^{2}}}{4\lambda_{v}}=\lim_{ \lambda_{u}>0,\lambda_{v}\to 0}\frac{-3\lambda_{u}-3\lambda_{u}\sqrt{1+ \frac{8\lambda_{v}^{2}}{9\lambda_{u}^{2}}}}{4\lambda_{v}}\\ \approx\lim_{\lambda_{u}>0,\lambda_{v}\to 0}\frac{-3 \lambda_{u}-3\lambda_{u}(1+\frac{1}{2}\frac{8\lambda_{v}^{2}}{9\lambda_{u}^{2} })}{4\lambda_{v}}=\lim_{\lambda_{u}>0,\lambda_{v}\to 0}-(\frac{3 \lambda_{u}}{2\lambda_{v}}+\frac{\lambda_{v}}{3\lambda_{u}})=\lim_{\lambda_{u} >0,\lambda_{v}\to 0}-\frac{3\lambda_{u}}{2\lambda_{v}}.\end{split}\] (A8)

Therefore, from Eq. (A7), it is clear that the local minimum solution \(\bar{\alpha}\) in Eq. (A6) will automatically reduce to zero if \(\lambda_{u}(t)<0\) and \(\lambda_{v}(t)=0\) at some isolated points, indicating that Eq. (A6) still holds in such case. On the other hand, from Eq. (A8), the sign of the local minimum solution \(\bar{\alpha}\) will change if \(\lambda_{v}(t)\) crosses zero and \(\lambda_{u}(t)>0\) at the isolated points, which will result in the appearance of the discontinuous jump. As a result, a discontinuous jump can be detected using the signs of \(\lambda_{v}\) and \(\lambda_{u}\).

Now we prove that the second case, that is, \(\lambda_{v}(t)\equiv 0\) in a time interval \(t\in[t_{a},t_{b}]\), does not hold along an optimal trajectory. By contradiction, we assume that \(\lambda_{v}(t)\equiv 0\) for \(t\in[t_{a},t_{b}]\). Recall that \(\lambda_{\theta}(t)\equiv 0\) along an optimal trajectory. Then, in view of Eq. (8), we have

\[\dot{\lambda}_{v}(t)=-2\frac{\lambda_{u}v(t)}{r(t)}\equiv 0,\forall\ t \in[t_{a},t_{b}].\] (A9)

Note that if \(v(t)\equiv 0\) in any time interval, it will result in \(\theta\) being constant in such time interval, as shown by Eq. (1). This is obviously impossible during an orbital transfer. Thus, to make Eq. (A9) true, \(\lambda_{u}(t)\) must be kept zero for \(t\in[t_{a},t_{b}]\). In this case, it implies

\[\lambda_{r}(t)\equiv 0,\forall\ t\in[t_{a},t_{b}].\] (A10)

Clearly, if \(\lambda_{v}(t)\equiv 0\) for \(t\in[t_{a},t_{b}]\) does hold, the equation as follows will be valid

\[\boldsymbol{\lambda}(t)=[\lambda_{r}(t),\lambda_{\theta}(t),\lambda_{u}(t), \lambda_{v}(t)]\equiv\boldsymbol{0},\forall\ t\in[t_{a},t_{b}],\] (A11)

which contradicts the PMP. Thus, \(\lambda_{v}(t)\) can only be zero at some isolated points.

Remember that the optimal guidance law is still ambiguous for the case that \(\lambda_{v}(t)\) crosses zero and \(\lambda_{u}(t)>0\), as shown in Eq. (A8). Because a nonlinear function takes its global
minimum at one of its local minima or at one of the endpoints of its feasible domain [42], we have

\[\alpha^{*}=\operatorname*{argmin}_{\alpha\in\{-\phi_{max},\bar{\alpha},\phi_{max} \}}\mathscr{H}. \tag{101}\]

To further resolve the ambiguity in terms of \(\alpha\), rewrite Eq. (100) as

\[\mathscr{H}(\alpha,\mathbf{\lambda},\mathbf{x})=\mathscr{H}_{1}(\alpha,\mathbf{\lambda}, \mathbf{x})+\mathscr{H}_{2}(\mathbf{\lambda},\mathbf{x}), \tag{102}\]

where \(\mathscr{H}_{1}\) is part of \(\mathscr{H}\) related to \(\alpha\), and \(\mathscr{H}_{2}\) is the rest of \(\mathscr{H}\) independent from \(\alpha\). Then, \(\mathscr{H}_{1}\) satisfies

\[\mathscr{H}_{1}(\alpha,\mathbf{\lambda},\mathbf{x})=\frac{\beta}{r^{2}}(\lambda_{u} \cos^{3}\alpha+\lambda_{v}\sin\alpha\cos^{2}\alpha). \tag{103}\]

Recall Eq. (3) and \(\phi_{max}\in(0,\frac{\pi}{2})\), thus we obtain

\[\mathscr{H}_{1}(\frac{\pi}{2},\mathbf{\lambda},\mathbf{x})=\mathscr{H}_{1}(-\frac{ \pi}{2},\mathbf{\lambda},\mathbf{x})=0. \tag{104}\]

Substituting Eq. (102) into Eq. (103) yields

\[\begin{cases}\mathscr{H}_{1}(\alpha_{1})=\frac{\beta}{4r^{2}}\cos^{3}\alpha \left(\lambda_{u}+\sqrt{9\lambda_{u}^{2}+8\lambda_{v}^{2}}\right)\geq 0,\\ \mathscr{H}_{1}(\alpha_{2})=\frac{\beta}{4r^{2}}\cos^{3}\alpha\left(\lambda_{u }-\sqrt{9\lambda_{u}^{2}+8\lambda_{v}^{2}}\right)\leq 0,\end{cases} \tag{105}\]

which indicates that \(\alpha_{1}\) and \(\alpha_{2}\) is the global maximum and global minimum for \(\mathscr{H}\) in \([-\frac{\pi}{2},\frac{\pi}{2}]\), respectively. Define a variable \(\Delta\mathscr{H}\) as

\[\Delta\mathscr{H}=\mathscr{H}_{1}(\phi_{max},\mathbf{\lambda},\mathbf{x})-\mathscr{H} _{1}(-\phi_{max},\mathbf{\lambda},\mathbf{x})=\frac{2\beta}{r^{2}}\cos^{2}\phi_{max} \sin\phi_{max}\lambda_{v}. \tag{106}\]

Without loss of generality, we consider the case that \(\lambda_{v}<0\), which leads to a negative local maximum point \(\alpha_{1}\) and a positive local minimum point \(\alpha_{2}\). If \(\alpha_{2}<\phi_{max}\), then \(\alpha_{2}\) is the global minimum point in \([-\phi_{max},\phi_{max}]\). If \(\alpha_{2}>\phi_{max}\), then \(\phi_{max}\) will be the global minimum in \([-\phi_{max},\phi_{max}]\) since \(\Delta\mathscr{H}<0\), as shown in Fig. 15. Hence, in view of Eq. (103), we have that the optimal guidance law \(\alpha^{*}\) minimizing the Hamiltonian \(\mathscr{H}\) w.r.t. \(\alpha\) is

\[\alpha^{*}=\text{Median}[-\phi_{max},\bar{\alpha},\phi_{max}],\text{where }\bar{\alpha}=\arctan\ \frac{-3\lambda_{u}-\sqrt{9\lambda_{u}^{2}+8\lambda_{v}^{2}}}{4\lambda_{v}}. \tag{107}\]

Proof of Lemma 2. Taking Eq. (17) into consideration, it is easy to see that the solution to the parameterized system in Eq. (15) at \(\tau=0\) represents a circular orbit with a radius of \(R_{0}\). Moreover, regarding Eq. (10), we have that Eq. (9) holds at \(\tau=0\) for the pair \((\lambda_{R_{0}},\lambda_{U_{0}})\) and \(\lambda_{V_{0}}\) determined by Eq. (19). By propagating the parameterized system in Eq. (15) for \(\tau\in[0,t_{f}]\), \(\mathcal{F}\) meets all the necessary conditions for an optimal trajectory. By the definition
of \(\mathcal{F}_{p}\), it is obvious that \(\mathcal{F}_{p}\) represents the solution space of an optimal trajectory, and \(\tau\in[0,t_{f}]\) defines the time of flight. In other words, an optimal trajectory can be readily obtained simply by arbitrarily choosing two parameters \(\lambda_{R_{0}}\) and \(\lambda_{U_{0}}\) and solving an initial value problem governed by Eqs. (15) and (17).
