# The Earth-Mars Transfer Trajectory Optimization of Solar Sail Based onhp-Adaptive Pseudospectral Method

Authors: Yufei Guo; Dongzhu Feng; Xin Wang; Cong Li; Yunzhao Liu
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## 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. Thehp-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 thehp-adaptive pseudospectral method is introduced, and the steps of solving the transfer trajectory optimization problem byhp-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 ofhp-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.


---

## Body

## 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 N, 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 tf of the N 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]. Thehp-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, thehp-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 ofhp-adaptive pseudospectral method and Gauss pseudospectral method are compared and analyzed in this paper. The results show that thehp-adaptive pseudospectral method is an effective trajectory optimization method.This paper is organized as follows. In Section2 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 thehp-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 by(1)Fi=PsAεcos⁡αn-sin⁡αtFt=PsAεcos⁡αn+sin⁡αtwhere Ps is the SRP, and Aε is the effective area of solar sail, Aε=Acos⁡α, and A is the sail area. The pitch angle α is the angle between the unit vector of incident rays r and the unit normal vector of the solar sail n, and n is perpendicular to the sail and opposite to the Sun. The unit tangent vector of the solar sail t is perpendicular to the normal n and its positive direction is counterclockwise. Then, the SRP resultant force acting on a solar sail can be described as [16](2)Fs=Fi+Ft=2PsAcos2αn=2PsAr·n2nFigure 1
Solar-radiation-pressure force model of an ideal solar sail.The SRP force is always along the direction of sail normal.f=n, where f 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 as(3)ap=βμr2r·n2nwhere μ=GM=1.3275×1020m3/s2 is solar gravitational coefficient, G is the gravitational constant, M is the mass of the Sun, and r 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 as(4)β=σ∗σwhere σ=m/A is the inverse of the area-to-mass ratio of the solar sail (area density), for the ideal sail, σ∗=1.53g/m2. The design parameter β plays a significant role in the structure and performance of the solar sail, and the SRP force can be represented as(5)Fs=map=βμmr2r·n2n


### 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(6)r¨+μr2r=βμr2r·n2nAs shown in Figure2, in the coplanar orbit polar coordinate system, the equations of motion are as follows:(7)r¨=rθ˙2-μr2+arrθ¨=aθ-2r˙θ˙where ar=βμ/r2cos3α and aθ=βμ/r2cos2αsin⁡α denote the magnitude of the acceleration of the SRP in the radial direction and tangential direction, respectively. ar provides radially outward acceleration, and aθ can change orbital angular momentum.Figure 2
Coplanar polar coordinate system centered on the Sun.It can be seen from Figure2 that, in the process of orbital transfer, r and θ are the radial position and polar angle of the solar sail relative to the Sun, respectively. The angle between the velocity vector v and the tangential velocity vector vθ 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 Fs.(  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 byx=r,θ,vr,vθT, where vr and vθ are the radial velocity and tangential velocity of the solar sail transfer orbit, respectively. Then, (7) can be expressed as [20, 21] (8)r˙=vrθ˙=vθrv˙r=-μr2+vθ2r+g1x,αv˙θ=-vrvθr+g2x,αwhere g1x,α=βμ/r2cos3α and g2x,α=βμ/r2cos2αsin⁡α 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]:(9)J=Φxt0,t0,xtf,tf+∫t0tfgxt,ut,tdtwhere Φ is the performance index function of end value type (Mayer type), and g is the performance index function of integral type (Lagrange type); the state xt∈Rn, the initial time t0, and the terminal time tf are subject to the dynamic constraints(10)x˙t=fxt,ut,tt∈t0,tf,the boundary conditions(11)ϕxt0,t0,xtf,tf=0,and the inequality path constraints(12)Cxt,ut,t⩽0Cxt,ut,t∈RlThe control constraints are denoted byu∈U⊂Rm,U=u∈Rm:uL⩽u⩽uR, where uL and uR are the lower and upper bounds of the control variables, respectively.The multiple-interval nonlinear optimal control problem divides time intervalt0,tf in the above problem into K mesh intervals and uses t0<t1<⋯<tK=tf to represent K+1 mesh points (the mesh point is the connecting point of the adjacent mesh interval). Nkk=1,⋯,K-1 represents the number of collocation points in the kth mesh interval tk-1,tk. 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 t∈tk-1,tk is transformed to the variable τ∈-1,1 via the time domain transformation(13)τ=2t-tk-tk-1tk-tk-1,tk-1<tkMoreover, the following relationship remains:(14)dτdt=2tk-tk-1,k=1,⋯,KAfter transformation, the cost functional in (9) can be rewritten as(15)J=Φx1-1,t0,xK+1,tK+∑k=1Ktk-tk-12∫-1+1gxkτ,ukτ,τdτThe system differential equations, boundary constraints, and path constraints are given, respectively, as(16)dxkτdτ=tk-tk-12fxkτ,ukτ,τϕx1-1,t0,xK+1,tK=0tk-tk-12Cxkτ,ukτ,τ⩽0In addition, there are linkage constraints that need to be satisfied(17)xk+1-xk+1-1=0


## 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 by(1)Fi=PsAεcos⁡αn-sin⁡αtFt=PsAεcos⁡αn+sin⁡αtwhere Ps is the SRP, and Aε is the effective area of solar sail, Aε=Acos⁡α, and A is the sail area. The pitch angle α is the angle between the unit vector of incident rays r and the unit normal vector of the solar sail n, and n is perpendicular to the sail and opposite to the Sun. The unit tangent vector of the solar sail t is perpendicular to the normal n and its positive direction is counterclockwise. Then, the SRP resultant force acting on a solar sail can be described as [16](2)Fs=Fi+Ft=2PsAcos2αn=2PsAr·n2nFigure 1
Solar-radiation-pressure force model of an ideal solar sail.The SRP force is always along the direction of sail normal.f=n, where f 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 as(3)ap=βμr2r·n2nwhere μ=GM=1.3275×1020m3/s2 is solar gravitational coefficient, G is the gravitational constant, M is the mass of the Sun, and r 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 as(4)β=σ∗σwhere σ=m/A is the inverse of the area-to-mass ratio of the solar sail (area density), for the ideal sail, σ∗=1.53g/m2. The design parameter β plays a significant role in the structure and performance of the solar sail, and the SRP force can be represented as(5)Fs=map=βμmr2r·n2n


## 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(6)r¨+μr2r=βμr2r·n2nAs shown in Figure2, in the coplanar orbit polar coordinate system, the equations of motion are as follows:(7)r¨=rθ˙2-μr2+arrθ¨=aθ-2r˙θ˙where ar=βμ/r2cos3α and aθ=βμ/r2cos2αsin⁡α denote the magnitude of the acceleration of the SRP in the radial direction and tangential direction, respectively. ar provides radially outward acceleration, and aθ can change orbital angular momentum.Figure 2
Coplanar polar coordinate system centered on the Sun.It can be seen from Figure2 that, in the process of orbital transfer, r and θ are the radial position and polar angle of the solar sail relative to the Sun, respectively. The angle between the velocity vector v and the tangential velocity vector vθ 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 Fs.(  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 byx=r,θ,vr,vθT, where vr and vθ are the radial velocity and tangential velocity of the solar sail transfer orbit, respectively. Then, (7) can be expressed as [20, 21] (8)r˙=vrθ˙=vθrv˙r=-μr2+vθ2r+g1x,αv˙θ=-vrvθr+g2x,αwhere g1x,α=βμ/r2cos3α and g2x,α=βμ/r2cos2αsin⁡α 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]:(9)J=Φxt0,t0,xtf,tf+∫t0tfgxt,ut,tdtwhere Φ is the performance index function of end value type (Mayer type), and g is the performance index function of integral type (Lagrange type); the state xt∈Rn, the initial time t0, and the terminal time tf are subject to the dynamic constraints(10)x˙t=fxt,ut,tt∈t0,tf,the boundary conditions(11)ϕxt0,t0,xtf,tf=0,and the inequality path constraints(12)Cxt,ut,t⩽0Cxt,ut,t∈RlThe control constraints are denoted byu∈U⊂Rm,U=u∈Rm:uL⩽u⩽uR, where uL and uR are the lower and upper bounds of the control variables, respectively.The multiple-interval nonlinear optimal control problem divides time intervalt0,tf in the above problem into K mesh intervals and uses t0<t1<⋯<tK=tf to represent K+1 mesh points (the mesh point is the connecting point of the adjacent mesh interval). Nkk=1,⋯,K-1 represents the number of collocation points in the kth mesh interval tk-1,tk. 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 t∈tk-1,tk is transformed to the variable τ∈-1,1 via the time domain transformation(13)τ=2t-tk-tk-1tk-tk-1,tk-1<tkMoreover, the following relationship remains:(14)dτdt=2tk-tk-1,k=1,⋯,KAfter transformation, the cost functional in (9) can be rewritten as(15)J=Φx1-1,t0,xK+1,tK+∑k=1Ktk-tk-12∫-1+1gxkτ,ukτ,τdτThe system differential equations, boundary constraints, and path constraints are given, respectively, as(16)dxkτdτ=tk-tk-12fxkτ,ukτ,τϕx1-1,t0,xK+1,tK=0tk-tk-12Cxkτ,ukτ,τ⩽0In addition, there are linkage constraints that need to be satisfied(17)xk+1-xk+1-1=0


## 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 -1. 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 xk+1-xk+1-1=0 across mesh points by using the Radau scheme. For this reason, the RPM is adopted as the pseudospectral method ofhp-adaptive strategy in this paper.(  1) Collocation Points and Discretization. In mesh interval k, the RPM takes Nk order LGR points, which are the roots of the polynomial PNkτ+PNk-1τ, as the collocation points, where PNkτ is Nk order Legendre polynomial defined by(18)PNkτ=12NkNk!dNkdτNkτ2-1NkNk=0,1,2,⋯N k LGR points and τNk+1k=+1 are taken as the discrete nodes in the RPM, and Nk+1 Lagrange interpolation polynomials Ljkτj=1,⋯,Nk+1 are used as the basis functions to approximate the state variables xkτ within each mesh interval k∈1,⋯,K. So xkτ can be expressed as(19)xkτ≈Xkτ=∑j=1Nk+1LjkτXjkwhere the Lagrange interpolation basis function Ljkτ is(20)Ljkτ=∏i=1,i≠jNk+1τ-τikτjk-τikFurthermore, the interpolation basis functionsL~jkτj=1,⋯,Nk are used to approximate the control variables ukτ at the Nk LGR points in the kth mesh interval.(21)ukτ≈Ukτ=∑j=1NkL~jkτUjk(  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 as(22)x˙kτ≈X˙kτ=∑j=1Nk+1DijkXjkwhere Dijk is the differential matrix(23)Dijk=L˙jkτiki=1,⋯,Nkj=1,⋯,Nk+1Thus, the state constraints can be transformed into the equality constraints at theNk LGR collocation points(24)∑j=1Nk+1DijkXjk-tk-tk-12fXik,Uik,τik=0i=1,2,⋯,Nk(  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 as(25)min⁡J=ΦX11,t0,XNk+1K,tK+∑k=1K∑j=1Nktk-tk-12ωjkgXjk,Ujk,τjkwhere ωjk are the LGR weights at the kth mesh interval.(  4) Discretization of Constraint Functions. The boundary constraints of (11) are approximated as(26)ϕX11,t0,XNK+1K,tK=0The path constraints of (12) in mesh interval k∈1,⋯,K are enforced at the Nk LGR points as(27)CXik,Uik,τik⩽0i=1,⋯,NkFinally, continuity across the mesh points is maintained by the condition(28)XNk+1k=X1k+1k=1,⋯,K-1Here, 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 ofL points τ-1k,⋯,τ-Lk∈-1,1 is defined in the kth interval, the first-order differential residual of state variables at the sth sampling points ask will be(29)X˙kτsk-tk-tk-12fXsk,Usk,τsk=askThe residual of path constraintsbsk will be(30)CXsk,Usk,τsk=bskThen the maximum of approximation error in the current intervalemaxk can be described as(31)emaxk=max⁡asik,bsjkwhere s=1,…,L, i=1,…,n, where n is the dimension of state variables and j=1,…,l, where l is the number of the path constraints. It is limited by the error tolerance ε. If(32)emaxk≤ε,then 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 (p method) or refining mesh (h method) is determined by the curvature. According to the plane geometry, the curvature of the plane curve Xτ at the sth sampling point τsk is given by [22](33)κkτsk=X¨kτsk1+X˙kτsk23/2The relative curvature of thekth interval rkk is defined as(34)rkk=κmaxkκ-kwhere κmaxk and κ-k are the maximum curvature and the mean curvature for all L sampling points, respectively.Set the relative curvature threshold asrmax. If rkk≤rmax, the trajectory will be smoother, and a larger degree polynomial can be used to improve the approximation accuracy in interval k. If rkk>rmax, the subdivision time interval will be needed.(  1) Determination of New Polynomial Degree within a Mesh Interval. Suppose the polynomial degree within mesh interval k should be increased, the degree Nk of the polynomial in mesh interval k is given by(35)Nk=M+ceillog10emaxkε+Awhere M denotes the constant initial degree of the polynomial in this interval; ceil• is the rounding operator, and the adjustable factor A is an integer greater than 0. It can be seen from (35) that the larger the error emaxk, the greater the polynomial degree.(  2) Determination of Number and Placement of New Mesh Points. Suppose the kth mesh interval should be refined, the new number of mesh intervals, denoted nk, is given by(36)nk=ceilBlog10emaxkεwhere adjustable factor B is an integer greater than 0. It is seen from (36) that the larger the error emaxk, 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 asck of the kth mesh interval, where(37)ck∑s=1Nk+1κkτs=1Then theith subdivision grid point in the kth interval should be selected at the sampling point τP to satisfy(38)ck∑j=1Pκkτj≈ink1≤P≤L,1≤j≤P


### 3.3. Processes of hp-Adaptive Pseudospectral Method

The iterative steps to solve the transfer trajectory optimization problem of solar sail by thehp-adaptive pseudospectral method are shown in Table 1, and the iterative flow chart is shown in Figure 3.Table 1
Iterative steps ofhp-adaptive pseudospectral method.
Step  1 Initialization. Given the initial states of the solar sail, set the number of intervalsK, the number of initial collocation points M for each interval, the error threshold ε, and the relative curvature threshold rmax.  Step  2 The optimal control problem is discretized into the NLP problem with the selected nodes and is solved.  Step  3 Determine the relation between the maximum residualemaxk and ε in the current interval k: if emaxk≤ε, then continue to check the next interval; otherwise, proceed to step 4 or step 5.  Step  4 Ifrkk<rmax, determine the new number collocation points Nk of this mesh interval by (35), and then check the k+1th mesh interval.  Step  5 Ifrkk≥rmax, then refine the kth mesh interval into nk subintervals by (36), and set the number of collocation points on each of the subintervals to be M and proceed to the next interval.  Step  6 After all mesh intervals have been updated, return to Step 2 and proceed to the next iteration.Figure 3
Iterative process ofhp-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 -1. 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 xk+1-xk+1-1=0 across mesh points by using the Radau scheme. For this reason, the RPM is adopted as the pseudospectral method ofhp-adaptive strategy in this paper.(  1) Collocation Points and Discretization. In mesh interval k, the RPM takes Nk order LGR points, which are the roots of the polynomial PNkτ+PNk-1τ, as the collocation points, where PNkτ is Nk order Legendre polynomial defined by(18)PNkτ=12NkNk!dNkdτNkτ2-1NkNk=0,1,2,⋯N k LGR points and τNk+1k=+1 are taken as the discrete nodes in the RPM, and Nk+1 Lagrange interpolation polynomials Ljkτj=1,⋯,Nk+1 are used as the basis functions to approximate the state variables xkτ within each mesh interval k∈1,⋯,K. So xkτ can be expressed as(19)xkτ≈Xkτ=∑j=1Nk+1LjkτXjkwhere the Lagrange interpolation basis function Ljkτ is(20)Ljkτ=∏i=1,i≠jNk+1τ-τikτjk-τikFurthermore, the interpolation basis functionsL~jkτj=1,⋯,Nk are used to approximate the control variables ukτ at the Nk LGR points in the kth mesh interval.(21)ukτ≈Ukτ=∑j=1NkL~jkτUjk(  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 as(22)x˙kτ≈X˙kτ=∑j=1Nk+1DijkXjkwhere Dijk is the differential matrix(23)Dijk=L˙jkτiki=1,⋯,Nkj=1,⋯,Nk+1Thus, the state constraints can be transformed into the equality constraints at theNk LGR collocation points(24)∑j=1Nk+1DijkXjk-tk-tk-12fXik,Uik,τik=0i=1,2,⋯,Nk(  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 as(25)min⁡J=ΦX11,t0,XNk+1K,tK+∑k=1K∑j=1Nktk-tk-12ωjkgXjk,Ujk,τjkwhere ωjk are the LGR weights at the kth mesh interval.(  4) Discretization of Constraint Functions. The boundary constraints of (11) are approximated as(26)ϕX11,t0,XNK+1K,tK=0The path constraints of (12) in mesh interval k∈1,⋯,K are enforced at the Nk LGR points as(27)CXik,Uik,τik⩽0i=1,⋯,NkFinally, continuity across the mesh points is maintained by the condition(28)XNk+1k=X1k+1k=1,⋯,K-1Here, 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 ofL points τ-1k,⋯,τ-Lk∈-1,1 is defined in the kth interval, the first-order differential residual of state variables at the sth sampling points ask will be(29)X˙kτsk-tk-tk-12fXsk,Usk,τsk=askThe residual of path constraintsbsk will be(30)CXsk,Usk,τsk=bskThen the maximum of approximation error in the current intervalemaxk can be described as(31)emaxk=max⁡asik,bsjkwhere s=1,…,L, i=1,…,n, where n is the dimension of state variables and j=1,…,l, where l is the number of the path constraints. It is limited by the error tolerance ε. If(32)emaxk≤ε,then 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 (p method) or refining mesh (h method) is determined by the curvature. According to the plane geometry, the curvature of the plane curve Xτ at the sth sampling point τsk is given by [22](33)κkτsk=X¨kτsk1+X˙kτsk23/2The relative curvature of thekth interval rkk is defined as(34)rkk=κmaxkκ-kwhere κmaxk and κ-k are the maximum curvature and the mean curvature for all L sampling points, respectively.Set the relative curvature threshold asrmax. If rkk≤rmax, the trajectory will be smoother, and a larger degree polynomial can be used to improve the approximation accuracy in interval k. If rkk>rmax, the subdivision time interval will be needed.(  1) Determination of New Polynomial Degree within a Mesh Interval. Suppose the polynomial degree within mesh interval k should be increased, the degree Nk of the polynomial in mesh interval k is given by(35)Nk=M+ceillog10emaxkε+Awhere M denotes the constant initial degree of the polynomial in this interval; ceil• is the rounding operator, and the adjustable factor A is an integer greater than 0. It can be seen from (35) that the larger the error emaxk, the greater the polynomial degree.(  2) Determination of Number and Placement of New Mesh Points. Suppose the kth mesh interval should be refined, the new number of mesh intervals, denoted nk, is given by(36)nk=ceilBlog10emaxkεwhere adjustable factor B is an integer greater than 0. It is seen from (36) that the larger the error emaxk, 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 asck of the kth mesh interval, where(37)ck∑s=1Nk+1κkτs=1Then theith subdivision grid point in the kth interval should be selected at the sampling point τP to satisfy(38)ck∑j=1Pκkτj≈ink1≤P≤L,1≤j≤P


## 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 ofL points τ-1k,⋯,τ-Lk∈-1,1 is defined in the kth interval, the first-order differential residual of state variables at the sth sampling points ask will be(29)X˙kτsk-tk-tk-12fXsk,Usk,τsk=askThe residual of path constraintsbsk will be(30)CXsk,Usk,τsk=bskThen the maximum of approximation error in the current intervalemaxk can be described as(31)emaxk=max⁡asik,bsjkwhere s=1,…,L, i=1,…,n, where n is the dimension of state variables and j=1,…,l, where l is the number of the path constraints. It is limited by the error tolerance ε. If(32)emaxk≤ε,then 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 (p method) or refining mesh (h method) is determined by the curvature. According to the plane geometry, the curvature of the plane curve Xτ at the sth sampling point τsk is given by [22](33)κkτsk=X¨kτsk1+X˙kτsk23/2The relative curvature of thekth interval rkk is defined as(34)rkk=κmaxkκ-kwhere κmaxk and κ-k are the maximum curvature and the mean curvature for all L sampling points, respectively.Set the relative curvature threshold asrmax. If rkk≤rmax, the trajectory will be smoother, and a larger degree polynomial can be used to improve the approximation accuracy in interval k. If rkk>rmax, the subdivision time interval will be needed.(  1) Determination of New Polynomial Degree within a Mesh Interval. Suppose the polynomial degree within mesh interval k should be increased, the degree Nk of the polynomial in mesh interval k is given by(35)Nk=M+ceillog10emaxkε+Awhere M denotes the constant initial degree of the polynomial in this interval; ceil• is the rounding operator, and the adjustable factor A is an integer greater than 0. It can be seen from (35) that the larger the error emaxk, the greater the polynomial degree.(  2) Determination of Number and Placement of New Mesh Points. Suppose the kth mesh interval should be refined, the new number of mesh intervals, denoted nk, is given by(36)nk=ceilBlog10emaxkεwhere adjustable factor B is an integer greater than 0. It is seen from (36) that the larger the error emaxk, 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 asck of the kth mesh interval, where(37)ck∑s=1Nk+1κkτs=1Then theith subdivision grid point in the kth interval should be selected at the sampling point τP to satisfy(38)ck∑j=1Pκkτj≈ink1≤P≤L,1≤j≤P


## 3.3. Processes of hp-Adaptive Pseudospectral Method

The iterative steps to solve the transfer trajectory optimization problem of solar sail by thehp-adaptive pseudospectral method are shown in Table 1, and the iterative flow chart is shown in Figure 3.Table 1
Iterative steps ofhp-adaptive pseudospectral method.
Step  1 Initialization. Given the initial states of the solar sail, set the number of intervalsK, the number of initial collocation points M for each interval, the error threshold ε, and the relative curvature threshold rmax.  Step  2 The optimal control problem is discretized into the NLP problem with the selected nodes and is solved.  Step  3 Determine the relation between the maximum residualemaxk and ε in the current interval k: if emaxk≤ε, then continue to check the next interval; otherwise, proceed to step 4 or step 5.  Step  4 Ifrkk<rmax, determine the new number collocation points Nk of this mesh interval by (35), and then check the k+1th mesh interval.  Step  5 Ifrkk≥rmax, then refine the kth mesh interval into nk subintervals by (36), and set the number of collocation points on each of the subintervals to be M and proceed to the next interval.  Step  6 After all mesh intervals have been updated, return to Step 2 and proceed to the next iteration.Figure 3
Iterative process ofhp-adaptive pseudospectral method.


## 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(39)J=∫t0tf1dt=tf-t0The orbital radiusr is in units of AU, where 1AU=1.496×1011m represents the average distance from the Earth to the Sun. The unit of time 1TU=AU3/μ≈58.1255day. The problem is solved for two solar sails with different lightness number of β1=0.17 and β2=0.1.Initial conditions:(40)r0=1,θ0=0.1,vr0=0,vθ0=1Terminal constraint:(41)rf=1.524,vrf=0,vθf=11.524Control variable constraints:(42)0≤α≤π2The 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 thehp-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 forhp-adaptive method is set to 10-6. The optimized results are shown in Figures 4–9.Figure 4
The optimized trajectory of the solar sail from Earth orbit to Mars orbit.Figure 5
The variation curve of orbit transfer radius.Figure 6
The radial velocity of transfer orbit.Figure 7
The angular position of transfer orbit.Figure 8
The tangential velocity of transfer orbit.Figure 9
The required pitch angle of the optimized trajectory.For the solar sail with the lightness number ofβ1=0.17, the orbit transfer times optimized byhp-adaptive method and GPM algorithm are 406.641 days and 406.638 days, respectively. For β2=0.1, the orbit transfer times optimized byhp-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 β1=0.17 and β2=0.1 have experienced a total spiral motion of about 4.426rad(253.59∘) and 5.530rad(317.01∘), 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 1.7∘/day. 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 byhp-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 byhp-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 byhp-adaptive method and GPM algorithm are shown in Table 2.Table 2
Average residuals between the state variable values and the real values.
β Algorithms State Variables 1 0 - 7 Δ r / A U Δ θ / r a d Δ v r / k m ⋅ s - 1 Δ v θ / k m ⋅ s - 1 β 1 = 0.17 hp-adaptive method 12.75 11.40 55.69 22.07 GPM 36.58 73.60 145.2 40.07 β 2 = 0.1 hp-adaptive method 3.59 2.29 14.38 1.69 GPM 4.52 46.93 44.52 14.63It is visible that the average residuals ofhp-adaptive method are small. Compared with GPM algorithm, the accuracies of the four state variables calculated byhp-adaptive algorithm for β1=0.17 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 thehp-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, thehp-adaptive pseudospectral method satisfies the first-order necessary conditions of optimal control theory.Figure 10
The curve of Hamiltonian-time.(  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, thehp-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 thehp-adaptive method has a good ability of adaptively adjusting the collocation points. The performance comparison of the two methods is shown in Table 3.Table 3
Performance comparison betweenhp-adaptive method and GPM.
β Algorithms The number of collocation points The number of variables The number of iterations Time consuming /s β 1 = 0.1 7 hp-adaptive method 52 265 1045 7.611 GPM 60 309 1147 18.402 β 2 = 0.1 hp-adaptive method 54 275 1050 8.131 GPM 78 399 1168 21.150The 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 Table3 show that thehp-adaptive method uses fewer collocation points and has a higher efficiency. This is because thehp-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, thehp-adaptive method has greater potential in online trajectory optimization.


## 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(39)J=∫t0tf1dt=tf-t0The orbital radiusr is in units of AU, where 1AU=1.496×1011m represents the average distance from the Earth to the Sun. The unit of time 1TU=AU3/μ≈58.1255day. The problem is solved for two solar sails with different lightness number of β1=0.17 and β2=0.1.Initial conditions:(40)r0=1,θ0=0.1,vr0=0,vθ0=1Terminal constraint:(41)rf=1.524,vrf=0,vθf=11.524Control variable constraints:(42)0≤α≤π2The 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 thehp-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 forhp-adaptive method is set to 10-6. The optimized results are shown in Figures 4–9.Figure 4
The optimized trajectory of the solar sail from Earth orbit to Mars orbit.Figure 5
The variation curve of orbit transfer radius.Figure 6
The radial velocity of transfer orbit.Figure 7
The angular position of transfer orbit.Figure 8
The tangential velocity of transfer orbit.Figure 9
The required pitch angle of the optimized trajectory.For the solar sail with the lightness number ofβ1=0.17, the orbit transfer times optimized byhp-adaptive method and GPM algorithm are 406.641 days and 406.638 days, respectively. For β2=0.1, the orbit transfer times optimized byhp-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 β1=0.17 and β2=0.1 have experienced a total spiral motion of about 4.426rad(253.59∘) and 5.530rad(317.01∘), 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 1.7∘/day. 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 byhp-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 byhp-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 byhp-adaptive method and GPM algorithm are shown in Table 2.Table 2
Average residuals between the state variable values and the real values.
β Algorithms State Variables 1 0 - 7 Δ r / A U Δ θ / r a d Δ v r / k m ⋅ s - 1 Δ v θ / k m ⋅ s - 1 β 1 = 0.17 hp-adaptive method 12.75 11.40 55.69 22.07 GPM 36.58 73.60 145.2 40.07 β 2 = 0.1 hp-adaptive method 3.59 2.29 14.38 1.69 GPM 4.52 46.93 44.52 14.63It is visible that the average residuals ofhp-adaptive method are small. Compared with GPM algorithm, the accuracies of the four state variables calculated byhp-adaptive algorithm for β1=0.17 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 thehp-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, thehp-adaptive pseudospectral method satisfies the first-order necessary conditions of optimal control theory.Figure 10
The curve of Hamiltonian-time.(  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, thehp-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 thehp-adaptive method has a good ability of adaptively adjusting the collocation points. The performance comparison of the two methods is shown in Table 3.Table 3
Performance comparison betweenhp-adaptive method and GPM.
β Algorithms The number of collocation points The number of variables The number of iterations Time consuming /s β 1 = 0.1 7 hp-adaptive method 52 265 1045 7.611 GPM 60 309 1147 18.402 β 2 = 0.1 hp-adaptive method 54 275 1050 8.131 GPM 78 399 1168 21.150The 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 Table3 show that thehp-adaptive method uses fewer collocation points and has a higher efficiency. This is because thehp-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, thehp-adaptive method has greater potential in online trajectory optimization.


## 5. Conclusions

In this paper, the application ofhp-adaptive pseudospectral method in the solar sail Earth-Mars transfer trajectory optimization is studied. The basic principle and the iterative process ofhp-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 thehp-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.



---

*Source: 6916848-2018-09-02.xml*