Abstract
Solar sails have many advantages over traditional chemical propulsion spacecraft, such as needlessness of fuel, high payload ratio, long service life, and great application potential in deep space exploration, interstellar voyage, and other aerospace application fields. However, the period of solar sail transfer from the initial orbit to the desired orbit is relatively long. Thus, it is necessary to optimize the transfer trajectory of solar sail according to specific tasks. The hp-adaptive pseudospectral method combines the global pseudospectral method with the finite element method, which adopts double layer optimization strategy to solve the optimal control problem and has higher computational efficiency and accuracy than the traditional global pseudospectral method. In this paper, the pressure of solar light acting on the solar sail is analyzed, and the kinematic equation of the solar sail is established in polar coordinate system first; then the basic principle of the hp-adaptive pseudospectral method is introduced, and the steps of solving the transfer trajectory optimization problem by hp-adaptive pseudospectral method are proposed; finally, the trajectory optimization of the solar sail from Earth orbit to Sun-centered Mars orbit is simulated as an example to demonstrate the effectiveness of hp-adaptive pseudospectral method in the orbit transfer optimization problem of solar sail. The simulation results show that the adopted method is not sensitive to the initial values and has more reasonable distribution and less computational cost than the Gauss pseudospectral method.
1. Introduction
With the increasing complexity and space distance of deep space exploration missions, the conventional propulsion methods hardly meet the demands. Solar sail is a new type of spacecraft that produces propulsion by a large area of light film that reflects solar photons, where the power source is natural solar energy, without consuming the chemical fuel. Therefore, the lifespan of the solar sail in space is not limited by fuel. The impact force of solar photons is very weak, but they continuously act on the solar sail, which will generate continuous small propulsion and persistent acceleration. The speed of the solar sail can reach up to 93km/s, which is 5 to 7 times the speed of the current fastest spacecraft [1]. In addition, the characteristics of light structure, low cost, and small launching risk make the solar sail have great advantages in long period space missions such as deep space exploration and interplanetary navigation.
The key technologies of solar sail include structure design, material selection, folding and deployment of solar sail, attitude control, and trajectory optimization. After nearly a century of development, solar sail technologies have advanced greatly. On the one hand, many scholars have made in-depth theoretical analysis of solar sail; on the other hand, a number of space agencies have carried out some ground or space tests on the key technologies of solar sail. Trajectory optimization for spacecraft can prolong the orbital life and increase the ability of performing tasks [2]. Thus, it is necessary to optimize the transfer trajectory of solar sail to shorten the orbital transfer time during deep space exploration missions.
Trajectory optimization of solar sail is essentially a nonlinear dynamic optimal control problem with state constraints and control constraints and is very difficult to solve. At present, the research on the trajectory optimization of solar sail is still in the early stage, which needs further theoretical research and experimental verification. Carl G. Sauer [3] (in 1976) used the Calculus of Variations of the classical indirect method to optimize the solar sail rendezvous trajectories for each of the terrestrial planets (Venus, Mercury, and Mars) and an asteroid roundtrip mission of a solar sail. The results proved the versatility of the solar sail spacecraft. However, the indirect method has some disadvantages including the difficulty of deducing the optimal condition, the small radius of convergence, and the difficulty of predicting the initial value of the covariate variable. Michiel Otten [4] (in 2001) obtained the nearly minimum transfer time for Earth-Mars coplanar trajectory of solar sail by using the direct discrete method. The author divided the whole orbital transfer process into equal periods of , and the control variables in each segment (i.e., the solar sail steering angle ) were constants, so that the optimized parameters were the values of and the terminal time of the time periods. Although the steering law obtained by this method is easier to implement for a realistic solar sail, the orbit transfer time of the optimal solution is longer than that of the indirect method. Nassiri [5] (in 2005) introduced the force model and orbital dynamics model for an ideal solar sail and used the direct collocation method to cast the solar sail trajectory optimization problem as a nonlinear programming problem. It is concluded that the direct collocation method is capable of finding accurate solutions to the optimal control problem of solar sail interplanetary trajectories. Gau [6] (in 2015) proposed a method using a solar sail to change the orbit of a near Earth asteroid. He derived the dynamic model for a tethered system formed by an asteroid and a solar sail and utilized the Legendre pseudospectral method to solve the optimal control problem of its trajectory. As an improvement of the direct collocation method, the pseudospectral method can obtain higher solution accuracy with less computational cost. Moreover, some scholars have applied intelligent optimization algorithms to solve the interplanetary optimal transfer trajectory problem for a solar sail, such as Evolutionary Neurocontrol [7] and Generalized Extremal Optimization [8, 9]. The intelligent optimization algorithms have the advantages of strong robustness and global convergence, but they are prone to produce premature phenomenon, and their local optimization ability is poor.
Since the collocation points of the pseudospectral method are fixed (dense near the boundaries and sparse in the middle) [10], the interpolation polynomial with higher dimensions is often needed to obtain an ideal approximate solution [11]. The hp-adaptive pseudospectral method [12–14] proposed by C. L. Darby and A. V. Rao et al. has the characteristics of flexibility of collocation points distribution, sparsity of calculation, and fast convergence and can adjust adaptively the number of the units or the degree of the polynomial in the corresponding segments according to the requirement of calculation accuracy. This method combines the advantages of the global pseudospectral method and the finite element method and can greatly reduce the computational cost. Here, the hp-adaptive pseudospectral method is adopted to solve the trajectory optimization problem for a solar sail during the transfer from Earth orbit to Mars orbit; furthermore, the collocation points distribution and time-consuming of hp-adaptive pseudospectral method and Gauss pseudospectral method are compared and analyzed in this paper. The results show that the hp-adaptive pseudospectral method is an effective trajectory optimization method.
This paper is organized as follows. In Section 2 we establish the kinematic equation of the solar sail. In Section 3 we present the steps for solving the transfer trajectory optimization problem based on the hp-adaptive pseudospectral method. In Section 4 we take the solar sail Earth-Mars transfer trajectory optimization as an example for the simulation and provide a discussion and analysis of the results. In Section 5 we provide concluding remarks.
2. Kinematics Model of the Ideal Solar Sail
We assume that the solar sail is an ideal plane and perfectly reflecting, and the orbital inclination of the desired planetary orbit is ignored. Therefore, the three-dimensional trajectory optimization problem can be simplified into a two-dimensional trajectory optimization problem.
2.1. Solar-Radiation-Pressure Model
The driving force of the solar sail is produced by momentum exchange between the solar photons and the highly reflective sail, so the solar-radiation-pressure (SRP) plays a dominant role on solar sail [15]. As shown in Figure 1, considering an ideal solar sail, the solar pressure energy is generated by the momentum transfer of incident light and reflected light with the sail. Since the solar photons are ideally reflected on the ideal plane sails, the magnitudes of the forces generated by the incident and reflected photons are equal, and they are given bywhere is the SRP, and is the effective area of solar sail, , and is the sail area. The pitch angle is the angle between the unit vector of incident rays and the unit normal vector of the solar sail , and is perpendicular to the sail and opposite to the Sun. The unit tangent vector of the solar sail is perpendicular to the normal and its positive direction is counterclockwise. Then, the SRP resultant force acting on a solar sail can be described as [16]
The SRP force is always along the direction of sail normal. , where is the thrust unit vector, and it is always in the direction of the SRP force and deviates from the Sun.
According to McInnes [17], The SRP acceleration can be expressed aswhere is solar gravitational coefficient, is the gravitational constant, is the mass of the Sun, and is the distance between the solar sail and the Sun. The dimensionless quantity corresponds to the sail lightness number (also known as the SRP factor) and is defined as the ratio between the SRP acceleration and the solar gravitational attraction [18, 19]. is defined aswhere is the inverse of the area-to-mass ratio of the solar sail (area density), for the ideal sail, . The design parameter plays a significant role in the structure and performance of the solar sail, and the SRP force can be represented as
2.2. Kinematics Equations
(1) Motion Equations in the Polar Coordinate System. Consider the two-body problems of an ideal solar sail in the Sun-centered ecliptic movement. The kinematic equation is
As shown in Figure 2, in the coplanar orbit polar coordinate system, the equations of motion are as follows:where and denote the magnitude of the acceleration of the SRP in the radial direction and tangential direction, respectively. provides radially outward acceleration, and can change orbital angular momentum.
It can be seen from Figure 2 that, in the process of orbital transfer, and are the radial position and polar angle of the solar sail relative to the Sun, respectively. The angle between the velocity vector and the tangential velocity vector is . Because the solar sail has the characteristics of continuous small propulsion and the SRP acceleration is inversely proportional to the square of the orbital radius, the interplanetary orbit transfer of solar sail can be achieved by using a spiral orbit. The orbit transfer can be realized by controlling the pitch angle of solar sail to adjust the SRP resultant force .
(2) State Differential Equation of Motion. Consider the optimal trajectory transfer problem of coplanar circular orbits centered on the Sun: by optimizing the control variables, the time to arrive to the target orbit is minimized while satisfying the orbital motion conditions and terminal constraints.
The state vector is given by , where and are the radial velocity and tangential velocity of the solar sail transfer orbit, respectively. Then, (7) can be expressed as [20, 21] where and are, respectively, the radial acceleration and the tangential acceleration produced by the SRP force, and the control variable is the pitch angle .
2.3. Description of the Transfer Trajectory Optimization Problem of Solar Sail
The transfer trajectory optimization problem for a solar sail can be described as a general form of optimal control problem: determining control variables with the consideration of dynamic constraints, path constraints, and boundary constraints, so as to minimize the cost functional in general Bolza form [22]:where is the performance index function of end value type (Mayer type), and is the performance index function of integral type (Lagrange type); the state , the initial time , and the terminal time are subject to the dynamic constraintsthe boundary conditionsand the inequality path constraints
The control constraints are denoted by , where and are the lower and upper bounds of the control variables, respectively.
The multiple-interval nonlinear optimal control problem divides time interval in the above problem into mesh intervals and uses to represent mesh points (the mesh point is the connecting point of the adjacent mesh interval). represents the number of collocation points in the th mesh interval . Therefore, the number of nodes in the whole interval is equal to the sum of the mesh points number and the collocation points number. The variable is transformed to the variable via the time domain transformationMoreover, the following relationship remains:
After transformation, the cost functional in (9) can be rewritten as
The system differential equations, boundary constraints, and path constraints are given, respectively, as
In addition, there are linkage constraints that need to be satisfied
3. hp-Adaptive Pseudospectral Method
3.1. Radau Pseudospectral Method
The basic principle of the global pseudospectral method for solving optimal control problems is that the state variables and the control variables are discretized on a series of collocation points, and the Lagrange interpolation polynomials are constructed by these discrete points as nodes to approximate the state variables and the control variables. Then, the time derivatives of the state variables are approximated by the derivatives of the global interpolation polynomials, which convert the constraints of the trajectory differential equations into a set of algebraic constraints. Next, the integral term in the cost function is calculated by the Gauss integral, and the terminal state is obtained from the initial state and the integral of the right function [2]. Through the above transformation, the optimal control problem can be transformed into a parameter optimization problem governed by a series of algebraic constraints, that is, a nonlinear programming (NLP) problem.
There are many kinds of pseudospectral methods. Different collocation points and node positions are adopted to make various pseudospectral methods different in numerical approximation methods. Among the pseudospectral methods applied in Aeronautics and Astronautics field, we find the Legendre pseudospectral method (LPM) [23], the Gauss pseudospectral method (GPM) [10], and the Radau pseudospectral method (RPM) [24]. The Radau pseudospectral method uses the Legendre–Gauss–Radau (LGR) collocation points that contain the point . Besides, the discrete form of the optimal condition of NLP problem is consistent with the discrete form of the optimal condition of the original time-continuous problem. Hence, it is particularly easy to implement the continuity conditions across mesh points by using the Radau scheme. For this reason, the RPM is adopted as the pseudospectral method of hp-adaptive strategy in this paper.
(1) Collocation Points and Discretization. In mesh interval , the RPM takes order LGR points, which are the roots of the polynomial , as the collocation points, where is order Legendre polynomial defined by
LGR points and are taken as the discrete nodes in the RPM, and Lagrange interpolation polynomials are used as the basis functions to approximate the state variables within each mesh interval . So can be expressed aswhere the Lagrange interpolation basis function is
Furthermore, the interpolation basis functions are used to approximate the control variables at the LGR points in the th mesh interval.
(2) Transformation of the State Equations. After parameterizing the state variables with the interpolation polynomials, the differential operations of states can be approximated by the differential operations of the interpolation basis functions, such aswhere is the differential matrix
Thus, the state constraints can be transformed into the equality constraints at the LGR collocation points
(3) Transformation of the Cost Functional. The integral term in the cost functional of (15) can be approximated by the Gauss integral. Hence, we can obtain the approximate performance index function of the RPM aswhere are the LGR weights at the th mesh interval.
(4) Discretization of Constraint Functions. The boundary constraints of (11) are approximated as
The path constraints of (12) in mesh interval are enforced at the LGR points as
Finally, continuity across the mesh points is maintained by the condition
Here, the parameter optimization problem that arises from the RPM to minimize the cost function of (25) is subject to the algebraic constraints of (24) and (26)~(28) and can be solved by appropriate NLP algorithms.
3.2. hp-Adaptive Strategy
3.2.1. Evaluation Criteria of Approximation Error in a Mesh Interval
Generally, in the process of numerical solution, the exact solution of the original continuous time optimal control problem cannot be obtained. However, if the numerical method can approximate the original problem infinitely, then the differential-algebraic constraint equation should be satisfied at any point. Therefore, the satisfaction degree of the differential-algebraic constraint equation among the collocation points can be used as the approximate error for the solution. Suppose a set of points is defined in the th interval, the first-order differential residual of state variables at the th sampling points will be
The residual of path constraints will be
Then the maximum of approximation error in the current interval can be described aswhere , , where is the dimension of state variables and , where is the number of the path constraints. It is limited by the error tolerance . Ifthen the state and control values meet the accuracy tolerance and the iterative computation in the interval will be stopped. Otherwise, the current interval needs to be modified by either increasing the degree of the approximating polynomial or dividing the current interval into more intervals to improve accuracy until the accuracy tolerance is satisfied in all intervals.
3.2.2. Determination of Increasing Collocation Points or Refining Mesh
Increasing polynomial degree ( method) or refining mesh ( method) is determined by the curvature. According to the plane geometry, the curvature of the plane curve at the th sampling point is given by [22]
The relative curvature of the th interval is defined aswhere and are the maximum curvature and the mean curvature for all sampling points, respectively.
Set the relative curvature threshold as . If , the trajectory will be smoother, and a larger degree polynomial can be used to improve the approximation accuracy in interval . If , the subdivision time interval will be needed.
(1) Determination of New Polynomial Degree within a Mesh Interval. Suppose the polynomial degree within mesh interval should be increased, the degree of the polynomial in mesh interval is given bywhere denotes the constant initial degree of the polynomial in this interval; is the rounding operator, and the adjustable factor is an integer greater than 0. It can be seen from (35) that the larger the error , the greater the polynomial degree.
(2) Determination of Number and Placement of New Mesh Points. Suppose the th mesh interval should be refined, the new number of mesh intervals, denoted , is given bywhere adjustable factor is an integer greater than 0. It is seen from (36) that the larger the error , the more the number of subdivision intervals.
In general, the larger mean curvature of a trajectory in a section indicates the trajectory of the segment more oscillatory. Under the same interpolation polynomial degree, the time span of the interval should be shortened to ensure the approximation accuracy. Based on this principle, the interval subdivision position can be determined.
Let us define the curvature density constant factor as of the th mesh interval, where
Then the th subdivision grid point in the th interval should be selected at the sampling point to satisfy
3.3. Processes of hp-Adaptive Pseudospectral Method
The iterative steps to solve the transfer trajectory optimization problem of solar sail by the hp-adaptive pseudospectral method are shown in Table 1, and the iterative flow chart is shown in Figure 3.
4. Numerical Simulations
4.1. Simulation Conditions
Deep space exploration focuses on the moon and Mars. The purpose and significance of Mars exploration is exploring the existence of life in Mars, detection of useful resources, and exploring the transformation of Mars and immigration prospects [25]. In this paper, the transfer of a solar sail from Earth orbit to Mars orbit is taken as an example, and the minimum transfer time is taken as the optimal index. Suppose both Earth orbit and Mars orbit are circular coplanar orbits centered on the Sun, the kinematic model of the solar sail is shown in (8). The objective function is given by
The orbital radius is in units of , where represents the average distance from the Earth to the Sun. The unit of time . The problem is solved for two solar sails with different lightness number of and .
Initial conditions:
Terminal constraint:
Control variable constraints:
The simulation has been performed with Matlab R2014a, running on a laptop with a CPU i5-3230m and a clock frequency of 2.60 GHz.
4.2. Analysis of Simulation Results
The Gauss pseudospectral method and the hp-adaptive pseudospectral method introduced in Section 2 have been used to solve the optimal control problem, respectively, and the SNOPT toolkit is used to optimize the calculation. Moreover, the accuracy requirement of collocation points for hp-adaptive method is set to . The optimized results are shown in Figures 4–9.
For the solar sail with the lightness number of , the orbit transfer times optimized by hp-adaptive method and GPM algorithm are 406.641 days and 406.638 days, respectively. For , the orbit transfer times optimized by hp-adaptive method and GPM algorithm are approximately 505.056 days. The smaller the lightness number of the solar sails is (i.e., the larger the area density is), the smaller the obtained characteristic acceleration is and the longer the orbital transfer time will be. The orbital radius variation curve of the solar sail transferring from Earth orbit to Mars orbit is shown in Figure 5, and it can be seen that the solar sail enters the Martian target orbit eventually and remains on the Mars orbit. We can see from Figure 6 that the radial velocity tends to zero when the solar sail reaches the target orbit. The angular position of transfer orbit and the tangential velocity variation curve are given in Figures 7 and 8, respectively, which shows that the solar sails with the lightness number of and have experienced a total spiral motion of about and , respectively, when they complete the orbital transfer. Moreover, it is easy to observe that the change of pitch angle is slow from Figure 9 and the maximum rate of change is about . Therefore, the attitude control of the solar sail is easy to achieve during the transfer process.
(1) Analysis of Simulation Accuracy. In order to analyze the accuracy of the optimal trajectory calculated by hp-adaptive method, the optimized optimal control quantities are substituted into the original differential equation to compute the forward integral of trajectory by ode45 command of Matlab. And then the real flight trajectory under the optimized control law is obtained. The results are shown as the red solid line in Figures 5–8. It can be seen that the optimal trajectory obtained by hp-adaptive method and the real trajectory are in good agreement, and the terminal conditions are basically satisfied. The average residuals between the real values and the state variable values obtained by hp-adaptive method and GPM algorithm are shown in Table 2.
It is visible that the average residuals of hp-adaptive method are small. Compared with GPM algorithm, the accuracies of the four state variables calculated by hp-adaptive algorithm for are increased by 65.14%, 84.51%, 61.65%, and 44.92%, respectively. This is because the fitting accuracy of the lower order polynomial for the smooth trajectory is higher than that of the high-order polynomial for the nonsmooth trajectory.
(2) Analysis of Optimality. The Hamilton function variation curve given by the hp-adaptive method is shown in Figure 10. It is known by the minimum principle that the Hamiltonian for this time optimal control problem should remain constant and is -1 (because the Hamiltonian is not an explicit function of time). The results in Figure 10 show that the Hamiltonian is basically -1 in the entire simulation process. Therefore, the hp-adaptive pseudospectral method satisfies the first-order necessary conditions of optimal control theory.
(3) Analysis of Collocation Points and Simulation Time. From the simulation results of Figures 5 and 7, it is visible that the distribution of collocation points of GPM algorithm shows the characteristics of being dense at both ends and being sparse in the middle. As can be seen from Figures 6, 8, and 9, the hp-adaptive method can automatically adjust the number of mesh intervals and the order of the basis function during the optimization process; i.e., there are fewer collocation points in places with small gradient change and more collocation points in places with large gradient change. It indicates that the hp-adaptive method has a good ability of adaptively adjusting the collocation points. The performance comparison of the two methods is shown in Table 3.
The fewer the number of collocation points, the fewer the number of variables to be optimized (i.e., the smaller the size of the NLP) and the fewer the number of iterations. The comparison results in Table 3 show that the hp-adaptive method uses fewer collocation points and has a higher efficiency. This is because the hp-adaptive method divides the oscillation area into several relatively smooth intervals, and the low-order interpolation polynomial is used to approximate the optimal trajectory in each interval. Compared with the high-order global polynomial approximation adopted by GPM algorithm, there is a decrease in the computational burden. Hence, the hp-adaptive method has greater potential in online trajectory optimization.
5. Conclusions
In this paper, the application of hp-adaptive pseudospectral method in the solar sail Earth-Mars transfer trajectory optimization is studied. The basic principle and the iterative process of hp-adaptive pseudospectral method are described in detail, and the mathematical simulation is carried out. By comparing with the results of Gauss pseudospectral method, it is shown that the hp-adaptive pseudospectral method has fewer collocation points, more reasonable distribution of points, and faster convergence rate and can effectively solve the real-time optimal control problem.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest.
Acknowledgments
This work was supported in part by the National Natural Science Foundation of China (no. 61473228).