Practical Optimization of Low-Thrust Minimum-Time Orbital Rendezvous in Sun-Synchronous Orbits

Jian Ma,Changxuan Wen and Chen Zhang

1University of Chinese Academy Sciences,Beijing,100049,China

2Key Laboratory of Space Utilization,Technology and Engineering Center for Space Utilization,Chinese Academy of Sciences,Beijing,100094,China

3School of Aerospace Engineering,Beijing Institute of Technology,Beijing,100081,China

ABSTRACT High-specific-impulse electric propulsion technology is promising for future space robotic debris removal in sun-synchronous orbits.Such a prospect involves solving a class of challenging problems of low-thrust orbital rendezvous between an active spacecraftand a free-flying debris.This study focuses on computing optimal lowthrust minimum-time many-revolution trajectories,considering the effects of the Earth oblateness perturbations and null thrust in Earth shadow.Firstly,a set of mean-element orbital dynamic equations of a chaser(spacecraft)and a target(debris)are derived by using the orbital averaging technique,and specifically a slow-changing state of the mean longitude difference is proposed to accommodate to the rendezvous problem.Subsequently,the corresponding optimal control problem is formulated based on the mean elements and their associated costate variables in terms of Pontryagin’s maximum principle,and a practical optimization procedure is adopted to find the specific initial costate variables,wherein the necessary conditions of the optimal solutions are all satisfied.Afterwards,the optimal control profile obtained in mean elements is then mapped into the counterpart that is employed by the osculating orbital dynamics.A simple correction strategy about the initialization of the mean elements,specifically the differential mean true longitude,is suggested,which is capable of minimizing the terminal orbital rendezvous errors for propagating orbital dynamics expressed by both mean and osculating elements.Finally,numerical examples are presented,and specifically,the terminal orbital rendezvous accuracy is verified by solving hundreds of rendezvous problems,demonstrating the effectiveness of the optimization method proposed in this article.

KEYWORDS Trajectory optimization;low thrust;many revolutions;orbital rendezvous;sun-synchronous orbits

1 Introduction

It is evidenced that the number of trackable in-orbit space debris has reached 14,790 among the total 20,475 space objects [1].Approximately 90% of the space debris are distributed in low Earth orbits (LEOs),and majority of them reside in sun-synchronous orbits (SSOs).The increase in the number of the debris may lead to mutual collisions known as the Kessler’s syndrome [2],which would result in the continuous increase of the debris population.Active debris removal seems to be the only effective means for cleaning orbital environment.Much attention on debris removal has been paid not only to debris manipulating approaches and proximity operations [3-6],but also to spacecraft far range rendezvous problems,such as those proposed by the 9th Global Trajectory Optimization Competition and the 8th Chinese Trajectory Optimization Competition.For the debris removal mission,the spacecraft is anticipated to be capable of rendezvousing with as many debris as possible using its limited propellant.High specific impulse of electric propulsion appears to be the most promising engine option,which can greatly decrease the propellent mass fraction of the spacecraft compared with chemical propulsion.As we know,electric propulsion only generates low thrust.Therefore,research on orbital rendezvous using low thrust between debris in SSOs is of great necessity.

Although extensive studies [7-15] have been devoted to low-thrust many-revolution orbital transfer problems (a spacecraft is steered to a given target orbit but the terminal phase angle on the target orbit is not constrained),few studies have addressed the low-thrust orbital rendezvous problems (LORPs) specifically in SSOs considering the practical use of solar electric propulsion(SEP).In LORPs,a chaser using solar electric propulsion would rendezvous with a free-flying target with least time of flight given a specific time epoch.The difficulties of solving such LORPs in LEOs or LLOs mainly lie in the following three aspects:(1) Both of the chaser and the target significantly experience Earth or lunar oblateness perturbations,and the final orbital positions at rendezvous may not be known in advance for minimum-time missions.(2) The practical SEP spacecraft only has low thrust-to-weight ratio,usually on the order of 10?5,resulting in the case that low-thrust orbital maneuvers take several days with a large number of orbital revolutions(e.g.,more than a hundred).(3) Low thrust must be turned off when the spacecraft encounters Earth or moon’s shadow (or solar eclipse) or actively saves fuel consumption,thus the discontinuity in thrust,compared to the purely continuous thrust,leads to different and complex optimal control profile that is much harder to model and optimize.This study aims to propose practical optimization techniques that resolve all above-mentioned difficulties.

In literature,LORPs are modeled as orbital phasing problems or station-change problems,and the representative research documents are given as follows.Edelbaum [16] studied the timeoptimal phasing maneuver in the same circular orbit without considering any perturbation.The results suggest that the thrust directs tangentially in the first half of the transfer and then reverses to the opposite direction.Titus [17] transformed the problem into a maximum-phasing-change problem under fixed time of flight.Hall et al.[18] solved the time-optimal coplanar LORP using indirect optimization method.Shang et al.[19] extended the phasing problem on circular orbits to that on general elliptical orbits using the direct optimization method.Zhao et al.[20]proposed an analytical costate approximate method for the minimum-time station-change problem in GEOs.Wang et al.[21] studied a minimum-time coplanar LORP based on sequential convex programming,which was verified through an Earth-to-Mars rendezvous example with up to 50 orbital revolutions,with the terminal transfer angle set eitherπor 0.

More realistic work of LORPs considers theJ2perturbations that significantly affect orbital evolutions in Earth orbits.Zhao et al.[22] solved the fuel-optimal inclined LORP between space debris in SSOs by using indirect optimization method,where theJ2perturbation was considered.A similar scenario was also considered by Olympio et al.[23] and Li et al.[24],wherein multiple shooting methods and homotopic approach were employed to improve the convergence of numerical iteration.To avoid the difficulty of solving the TPBVP resulting from the optimal control theory,the direct method [25] is introduced to convert the primary problem into a nonlinear programming problem through parameterization.However,the computation workload rapidly increases when dealing with the many-revolution orbital transfers because a large number of nodes are required for numerical iteration.Zhang et al.[26] studied the time-optimal inclined LORP between geosynchronous satellites with fixed initial and terminal states based on the analytical orbit propagation.However,it is merely valid by imposing stringent constraints,i.e.,a constant yaw steering profile,J2perturbations only,and nearly zero eccentricity.

This work aims to develop a practical optimization method for solving the minimum-time many-revolution LORPs in SSOs,considering the practical use of SEP.To guarantee rendezvous,a phasing difference is introduced as an independent state into the dynamic equations,which contributes to the optimization process.The original problem is formulated as an optimal control problem by using the orbital averaging technique,which makes trajectory propagation easier to include the effects of perturbations and Earth shadow in which the thrust is off.The corresponding optimal control problem is then solved such that the necessary conditions derived from Pontryagin’s maximum principle are all satisfied.By employing an osculating-mean elements transformation,the resulting optimal control derived from the mean orbital dynamics and the corresponding costate variables is mapped to the control law that is used in osculating orbital dynamics.A simple correction strategy about the initialization of the mean elements is introduced to minimize the terminal rendezvous orbital errors.All the numerical iterations involved in trajectory optimization are based on the averaged orbital dynamics,greatly simplifying the orbital propagation complexity.The proposed method in this paper is capable to solving the minimumtime many-revolution LORPs in which the number of passive coast arcs in Earth shadow may be over one or two hundred.

The remainder of the paper is organized as follows:Section 2 states the minimum-time LORPs.Section 3 presents the optimal control problem formulation based on the mean-element dynamic equations.Section 4 and Section 5 describe the methods of solving the trajectory optimization problem and transforming the optimal control profile into the control law for propagating osculating elements.Section 6 shows representative numerical examples.Section 7 summarizes the conclusions.

2 Problem Statement

2.1 Assumptions

In this study,the following three assumptions are made:

Assumption 1:The chaser is assumed to be a SEP spacecraft,whereas the target is a freeflying debris.Supposing the chaser or the debris is moving in a circular SSO with an altitude of approximately 800 km.The zonal harmonics perturbations up to 20 orders is considered only,because they are much significant than the other perturbations as shown in Tab.1.

Assumption 2:The low thrust produced by electric propulsion is supposed to be directed to any direction without the change of thrust magnitude,which means that the maximum thrust available may keep constant during the entire transfer phase when the spacecraft is not in Earth shadow.This assumption is reasonable if the spacecraft’s solar panel may rotate at least in onedegree-of-freedom.

Assumption 3:Initial orbital states of the chaser and the target are supposed to be not widely spaced,wherein the semi-major axis,eccentricity,and inclination are constrained by the concept of the near-circular SSOs.The right ascension of ascending nodes (RAAN) are supposed to be similar whereas the phase angle difference between the chaser and the target can be arbitrary.

Table 1:Orders of the perturbing accelerations in a typical SSO

Assumption 4:The navigation error and the control error are not considered currently.

The third assumption implies that the chaser may rendezvous with the target in its neighboring orbit only.For a widely RAAN spaced rendezvous case,the reasonable strategy is to divide the entire flight into two phases:Minimizing orbital plane difference and orbital rendezvous,as shown in Fig.1.The RAAN difference between the chaser and the target can be zeroed by adjusting the chaser’s orbital altitude,which leads to differential RAAN procession rate.Once the chaser reaches a target’s neighboring orbit,the phase for orbital rendezvous starts.

Figure 1:Illustration of two phases (left:widely spaced RAAN;middle:minimizing orbital plane difference;right:orbital rendezvous in this study)

2.2 Dynamic Equations

A set of non-singular modified equinoctial elements [27] (MEEs)is utilized to describe the states of the chaser and the target:

whereais the semi-major axis,eis the eccentricity,iis the inclination,Ω is the RAAN,ωis the argument of periapsis,andθis the true anomaly.

Letudenotes the acceleration resulting from the natural perturbations and the active control,the perturbed Keplerian motion governed by Gauss’ variational equations can be written as

whereμis the gravitational parameter of the Earth,uR,uT,uNare the components ofuin the radial-transverse-normal (RTN) reference frame centered in the spacecraft,andw=1+fcosL+gsinL,s2=1+h2+k2.

The accelerationuin Eq.(2) consists of two parts:The perturbing accelerationupfrom the natural perturbations and the active control accelerationuc

Based onAssumption 1(Section 2.1),the perturbations considered in this study include zonal harmonics perturbations only.Because they are conservative force,Lagrange planetary equation [27] is adopted to compute the zonal harmonics perturbationsfpinstead ofM·up,that is

In addition,the thrust amplitude of the chaser,T(xc),satisfies

whereTmaxis the maximum thrust amplitude available,L1andL2denote the exit and entry angles across the Earth shadow in each orbital revolution,respectively.BothL1andL2can be solved from the cylindrical Earth shadow model presented by [28].Meanwhile,the unit direction vector of the thrust forceαsatisfies

In this manner,the active control acceleration can be expressed as

wheremis the spacecraft mass.The mass flow rate is expressed as

whereIspis the specific impulse andge(=9.80665m/s2) is the gravitational acceleration at the sea level.

In Eq.(2),Drepresents the two-body Keplerian motion,which leads the sixth elementLto be a fast-changing variable.By contrast,the first five MEEs,p,f,g,h,k,are slow-changing variables defined as the following slow-changing vector

Then,the slow-changing vector is separated from the fast-changing variable to obtainx=Considering that no active control is imposed on the target,the osculating orbital dynamics of the chaser and the target can be readily obtained as

where the subscripts (c,d) represent the chaser and the target respectively,and

2.3 The LORP in SSOs

Assuming at the initial timet0,the state vectors of the chaser and the target are given by

The LORP discussed in this study aims to provide an optimal control profile for the chaser,so that at a future timethe chaser reaches the target,i.e.,satisfying

whereas the time of flight Δt=tf?t0is minimized simultaneously.The motion of the chaser and the target are governed by the osculating orbital dynamics expressed in Eq.(13).For convenience,this LORP described in osculating elements is stated asProblem 1.

The description inProblem 1is quite straightforward,but the problem is difficult to solve in reality because of the following two reasons.First,the practical thrust-to-weight ratio of the SEP spacecraft is generally small,which leads to a large number of orbital revolutions.As a result,direct numerical integration of Eq.(13) consumes an unacceptable amount of time during the numerical iteration for optimization.Second,the discontinuous thrust caused by Earth shadow is hard to formulate in an optimization problem using the osculating orbital dynamics Eq.(13).To avoid these difficulties,the orbital averaging technique is introduced in the following sections.

3 Optimal Control Problem Formulation Based on Orbital Averaging

3.1 Orbital Averaging

Referring to [9],by using the orbital averaging technique on Eq.(13),the averaged orbital dynamic equations can be obtained as:

where a symboldenotes the averaged value of variablexandis the mean of the first five MEEs,

and

Similarly,as to the sixth elementL,the mean longitude rate can be computed by

which represents the mean longitude difference between the chaser and the target.Thus,the full averaged orbital dynamics of the LORP can be formulated as

where

In this article,the averaged dynamic equations Eq.(23) are adopted to solve the optimal control problem for convenience.By using the orbital averaging technique,the discontinuous thrust caused by Earth shadow is easily handled by the definite integrals,which can be approximately computed by the numerical integration methods such as Gauss-Legendre quadrature.Moreover,all the mean elements on the left side of Eq.(23) change slowly and smoothly and the classical fourth-order Runge-Kutta method can be used for trajectory propagation.

3.2 Time Optimal Control Problem Formulation

Considering a time-optimal orbital rendezvous problem,the performance function is expressed as

The corresponding governing differential equations of costate variables can be written as

where the derivatives ofGc0,Gc1,Gc2,Gc3,Gc4are expressed as

According to Pontryagin’s maximum principle [29],the optimal controlu?c(=T?α?) is obtained if the Hamiltonian function takes the minimum,that is

where the admissible control domain is

Thus,the optimal thrust directionα?satisfies

while the optimal thrust amplitude should be maximum all the time,i.e.,The specific solution ofu?is attached in Appendix A.Consequently,if at the initial timet0,the initial conditionsare known,the time-optimal orbital rendezvous problem can be transformed into a two-point boundary-value problem,defined asProblem 2as follows:

Finding the specific value of the initial costate variablesas well as the final timetf,such that the final states obtained by propagating the Eqs.(23) and (27)-(29) fromt0totfcan meet the necessary conditions,namely the terminal conditions

and the transversality conditions

3.3 Initialization of the Mean Elements

The initial conditions of the chaser and the target are generally given in osculating MEEs by the navigation system,however,the averaged dynamic equation requires the mean elementswhich should be transformed from the osculating elements,Lc0,m0,,Ld0.The transformation is based on the singular value decomposition (SVD) algorithm used by Li [30],where the algorithm is applied to find the mean elements for station-keeping of GEO satellites.

First of all,according to the osculating dynamic equation Eq.(13),˙mis a piece-wise constant,which means that no short-term oscillation occurs form.Thus,

Considering that the thrust-to-weight ratio of the chaser is much smaller than the zonal harmonics perturbationsfpin SSOs,the effect of the active control can be reasonably neglected in a short time duration,for example only one orbital period.Thus,the elements’ transformations of the chaser are the same as that of the target.LetLbe the transformation,the details ofLare reported in Algorithm 1.

Algorithm 1:Initialization of the mean elements L Input:x0=■~x■0 L0■■images/BZ_187_556_532_577_578.png Output:x0=~x0■L0 images/BZ_187_747_532_768_578.png■.Step 1.Calculate the orbital period from the osculating initial element x0;P=2πimages/BZ_187_457_701_486_747.pngimages/BZ_187_727_701_755_747.png3 2(44)Step 2.Propagate the osculating orbital dynamic Eq.(7) on the time interval [t0,t0+P],where the initial conditions are the osculating x0 and the active control acceleration uc is set to be 0.Thus,a trajectory of the osculating elements without active control,S1,is obtained,which is expressed as S1:x(t;x0),t∈[t0,t0+P] (45)Step 3.Solve the averaged trajectory of S1,denoted by√p0 μimages/BZ_187_491_779_509_825.png1?f20?g20■S1,by using the SVD algorithm used by Li [30].As the time duration of S1 is only one orbital period,we only considered the secular terms in the SVD algorithm,which means that the S1 is a straight line,as shown in Fig.2.Step 4.Calculate the middle-value of x1(P/2).Step 5.Propagate the averaged orbital dynamic equations Eq.(17) on the same time interval[t0,t0+P],where the initial conditions are given by S1 at the time t=t0+P/2,which is denoted by x0=x0 and no thrust is included.Thus,a trajectory of the mean elements without active control,S0,is obtained,which is expressed as S0:x(t;x0),t∈[t0,t0+P] (46)Step 6.Calculate the middle-value of x0(P/2).Step 7.Calculate the initial mean elements as S0 at time t=t0+P/2,which is denoted by x0=x0+Δx=x0+x1(P/2)?x0(P/2) (47)

Figure 2:Initialization of mean elements using osculating elements

4 Trajectory Optimization with the Averaged Orbital Dynamics

Problem 2stated as a two-point boundary-value problem is still difficult to be solved directly using the shooting method because of the sensitivity of the initial costate variables.Thus,Problem 2is transformed into the parameter optimization problem that is solved by taking advantages of mathematical programming solvers already developed and a series of practical techniques.

It is observed that the state equations Eq.(17),the costate equations Eq.(23),as well as the optimal control law in Eqs.(39) and (40) are all homogeneous with respect to the costate variables (or Lagrange multipliers)except for the Hamiltonian function Eq.(26).Multiplying or dividing these multipliersby an arbitrary positive real number does not change the terminal states of the chaser and the target.Thus,if there are a set of Lagrange multipliers,whose initial value is,that satisfies Eq.(39) and

The necessary conditions left inProblem 2=0 can be easily achieved by multiplying a specific positive scaling factor,k(k＞0),that is

whereCis the corresponding terminal value of the Hamiltonian function with respect to.The new Lagrange multipliers

are then the specific solutions ofProblem 2,where all the necessary conditions involved are satisfied.

As a result,the approach for solvingProblem 2can be divided into two steps.First,the TPBVP is converted into a nonlinear programming problem (NLP),defined asProblem 3,where one of the transversality conditions,the equality constraint=0,is changed to an inequality constraint<1.Second,the specific scaling factor of the Lagrange multipliers is selected to satisfy the constraint=0.Problem 3is totally equivalent toProblem 2,but the convergence domain is enlarged.Details ofProblem 3are summarized in Tab.2.

Table 2:Nonlinear optimization model of Problem 3

The NLP is solved by sequential quadratic programming (SQP),which is a gradient based algorithm,and can be readily employed in many software tools.To improve the NLP convergence properties,the intelligent optimization algorithm,for example,the differential evolution (DE) [31]is employed to search good initial guesses.By using the normalization of the initial costate vector [32],the upper and lower bounds of the costate variables can be normalized to be [?1,1].As to the terminal time variabletf,the upper and lower bound can be estimated from the analytical phasing maneuver time Δt1[19] and the plane change time Δt2[24],and is obtained as follows:

whereσis the weighting function considering the Earth shadow,which represents the percentage of time the spacecraft is thrusting during one revolution:

SolvingProblem 3by both the DE and the SQP requires evolutions or iterations of optimization variables,and meanwhile propagating the state and costate equations.Owing to the use of orbital averaging,the computation complexity of such propagation is greatly simplified.

It is noted that DE usually costs much computational time to locate good initial guesses,while SQP obtains the converged solution much quickly if good initial guesses are available.

5 Trajectory Correction with the Osculating Orbital Dynamics

5.1 Control Law in Osculating Orbital Dynamic Equations

The optimal solution of orbital rendezvous is obtained in the solution space of mean elements.In the real mission,the osculating elements are obtained by the navigation system and considered to be the input to the control system.As a result,the optimal control profile expressed by mean elements needs to be transformed into the control law that steers the chaser using osculating elements.

First,the thrust amplitude,T(xc),satisfies

That is to say,the thrust keeps a maximum amplitude when the chaser is out of Earth shadow.The judgement of whether the chaser is in the shadow (the comparison betweenLcandL1,L2)is tested continuously during the propagation.

As to the thrust direction,the most straightforward way is to take the optimal control profile in mean elements as the control law in osculating elements,with the phase angleLis replaced,i.e.,:

whereLis the osculating phase angle and the mean costate variables,and the mean statesinvolved are computed by the interpolation of,andthat are optimization results fromProblem 3.The thrust direction is only related to the osculating phase angleL,while,andare pre-stored quantities.

An alternative approach is that the control law of the chaser is governed by all the osculating elements.As indicated in Reference [33],the influence of thrust on the mean elements is approximated as that acting on the corresponding osculating elements.It is reasonable that the osculating elementsand the mean elementsare deemed the same to form the control law.Thus,the thrust direction profile in osculating elements is simply obtained by substitutinginto Eq.(39)for:

where the mean costate variables,are computed by interpolation of,.In this case,the thrust direction is depending on the osculating elementsandL,whileandare pre-stored quantities that can be deemed the gains for the control law.

As a result,once theProblem 3is solved,the osculating trajectory of the chaser can be obtained by propagating the osculating orbital dynamic Eq.(13) from the initial timet0to the final timetf,where the initial conditions arexc0,m0and the active control satisfies Eq.(54) and Eq.(55) or Eq.(56).Based on our experiences of solving a large number of numerical examples,the two strategies of the control lawα?(t)perform almost the same in precisely propagating the osculating orbital dynamic equations with the same set of initial states.This fact implies that the replacement ofbyis empirically acceptable,and the second strategy in Eq.(56) is more closely linked to the osculating elements that are obtained by the navigation system.It is noted that Zhong et al.[34] proposed an onboard mean-elements estimator,but it is not used in this study and our experiences show that Eq.(56) makes the control law effective and simpler.

5.2 Accuracy Verification and Error Correction Strategy

When we numerically propagated the orbital dynamic equations expressed by osculating elements with the control laws in Eqs.(54) and (56) obtained by solvingProblem 3,we found that the terminal orbital rendezvous errors exist,namely the terminal states of the chaser may not precisely meet those of the target.These errors stem from a series of approximations,as listed in Tab.3.

Table 3:List of approximations

It is found that the most significant rendezvous error is the difference in the true longitude between the chaser and the target at the terminal time.Considering the true longitudeLis a fast-changing variable,a correction about the initial mean longitude differencecan effectively compensate for the rendezvous errors.Letηbe the correction factor,which is a single scalar;andbe the terminal position and velocity errors between the chaser and the target in osculating elements,the correction of the initial mean longitude difference is expressed as

As a result,the new problem is proposed,that is:finding the optimalη?that contributes to the minimum of objective function

whereαandβare the weight coefficients and are configured according to different emphasis on the terminal position and velocity errors.The optimization interval ofηis chosen as [?0.15,0.15]empirically,which means a slight correction ofplays an important role in eliminating the terminal rendezvous errors.This univariate optimization problem is readily solved by many mature numerical methods,such as the Golden section search and the Fibonacci search [35].For convenience,this strategy is termedη-correction in this work.

During the process ofη-correction,there are two layers of optimization iterations.The outer iteration is the univariate optimization with respect toη,and the inner iteration is the NLP optimization (solvingProblem 3) with the initial conditions changed (related toη),which can be readily solved by SQP from the solution with the initial conditions be=L(x0),whereη=0.WhenJ+reaches its minimum,Problem 3is solved to obtain the optimal control in the solution space of mean elements and orbital accuracy of the control laws is also met with osculating elements.In this study,the orbital rendezvous is deemed to be achieved ifJ+is small enough.

Theη-correction perturbs the initialand leads to adjustment of the costate variables or the gains of control law.This strategy can be applicable to other mean elements if necessary.Although empirical,it is effective for solving concerning problems,which is validated by several hundreds of examples shown in Section 6.

6 Numerical Example

In this section,two typical LORP examples in SSOs are provided in details,namely:(1) The quasi phasing maneuver problem and (2) The general LORP.In example (1),the first five MEEs of the chaser and the target are nearly the same initially,the chaser mainly implements the phasing maneuvers;in Example (2),the chaser and the target differs in the semi-major axis(100 km difference) and the phase location at the initial time.More examples are tested via continuation of the chaser’s initial states,which generated several hundreds of examples to validate the proposed method.

For simplicity,the initial time epoch is set at January 1,2025,00:00:00.000 UTCG for all simulation examples.The target is chosen as a real free-flying debris moving in an SSO,whose DISCOS ID is 23753.It is a rocket mission related object and the size is larger than 5 m2[36].The orbital states of the target at the initial time epoch can be obtained via the SGP4 propagation model.The chaser,a servicing SEP spacecraft,is supposed to be injected into a nearby orbit with the initial mass and the thruster’s configuration presented in Tab.4.The low thrust produced by electric propulsion is supposed to be directed to any direction and the maximum thrust available keeps constant during the entire transfer phase.The initial states,in terms of the classic orbital elements under the J2000 inertial frame,of the chaser and the target,are given in Tab.5.The perturbations include the zonal harmonics of the non-spherical Earth gravitational perturbations up to 20 orders of the EGM96.

Table 4:Configuration of the chaser

Table 5:Initial classical orbital elements of the chaser and the target α

As mentioned before in the initialization of the mean elements (Section 3.3),there are two cases for a rendezvous problem:leading and trailing,whereinis set positive and negative separately.It is uncertain in advance that the transfer time of which case will be shorter,as a result,in the following examples,both cases were solved.The integral interval for propagating the averaged dynamic equations is set as 2～3 revolutions,and the integral precision adopted to calculate the terminal errors in osculating elements by propagating the osculating dynamic equations is defined by a relative error tolerance of 1×10?10and an absolute error tolerance of 1×10?10.The weight coefficientsαandβare both set to be 1 to just weight the terminal errorsδrandδvequally.Meanwhile,the stopping criteria of the univariate optimization forη-correction is set to be 1×10?4.During the computation,all the unit of the range and velocity are normalized by Earth’s radiusReandseparately to get a better convergence.

6.1 Example 1:Quasi Phasing Maneuver Problem

In this example,the chaser and the target are supposed to be located in nearby SSOs with the orbital altitudes of around 810 km.WhenProblem 3is solved and theη-correction is made,the solutions of the leading case and trailing case are summarized in Tab.6.

Table 6:Solutions of Example 1

whereη?is the optimal correction factor,Δt?is the minimum flight time,Δrand Δvrepresent the corresponding terminal differences of position and velocity after correction,respectively.In addition,Δtonand Δmdenote the working time of the thruster and the fuel consumption during the maneuver.It shows that the trailing case results in a shorter flight time,thus,the trailing case is discussed further below.

The time history of the osculating and mean MEEs (p,f,g,h,k) and the longitude difference ΔLare shown in Fig.3.These results show that the chaser gradually approaches the target.In addition,the entire maneuver turns out to be a “catch-up” process,as the phase difference increases from negative to 0 monotonically.

Figure 3:Evolutions of p,f,g,h,k,Δ L in Example 1

The time histories of the thrust directionα?under the RTN reference frame are presented in Fig.4.It shows that the thrust direction changes periodically along the transfer phase.At the first half transfer,thecomponent of the thrust is negative,thus,decreasing the semi-major axis to catch up with the phase;while at the second half,it turns to positive to increase the semi-major axis to the desired target value.Thecomponent of the thrust keeps a relative high value to compensate for the orbital plane change throughout the flight as the change of semi-major axis leads to a continuous drift difference of RAAN between the chaser and the target.In addition,thecomponent of the thrust is imposed mainly at the middle part of the maneuver.

Figure 4:History of the control direction α?under RTN reference frame in Example 1

Under the J2000 inertial frame,trajectories of the chaser and the target are shown in Fig.5.At the initial time,the chaser and the target are located at the points represented by the black and green dots,respectively.At the final time,terminal positions of the chaser and the target are denoted by the black triangle and green pentagon.The overlapping of the terminal positions indicates that the chaser rendezvouses with the target.The chaser passes the Earth shadow in each revolution,as the blue lines represent the chaser enters the shadow in which the thruster is off.The total number of the complete revolutions for the chaser to rendezvous with the target is 127.

Figure 5:Trajectories of the chaser for the rendezvous in Example 1

6.2 Example 2:General LORP

In this example,the chaser and the target are supposed to be located in SSOs with the altitude of 900 km and 810 km separately.WhenProblem 3is solved and theη-correction is made,the solutions of the leading case and trailing case are summarized in Tab.7.Different from Example 1,the leading case in Example 2 takes a better result,which is discussed further below.

Table 7:Solutions of Example 2

The time history of the osculating and mean MEEs (p,f,g,h,k) and the longitude difference ΔLare shown in Fig.6.The entire maneuver turns out to be a process including both “waitingfor”and “catch-up,”as the phase difference changes from positive to negative and eventually to 0.

Figure 6:Evolutions of p,f,g,h,k,ΔL in Example 2

The time histories of the thrust directionα?under the RTN reference frame are presented in Fig.7.The thrust direction during the rendezvous still displays a periodic change along the transfer phase,but the whole evolution is quite different with that in Example 1.The sign of thecomponent is not exactly changed at the half and thecomponent is imposed mainly at the end part.Moreover,thecomponent is small at the beginning while in the later,the proportion is large where the main plane change is conducted.

Figure 7:History of the control direction α?under RTN reference frame in Example 2

Under the J2000 inertial frame,trajectories of the chaser and target are shown in Fig.8.The total number of complete revolutions for the chaser to rendezvous with the target is 234.

Figure 8:Trajectories of the chaser for the rendezvous in Example 2

6.3 Continuation of the Initial States of the Chaser

In this subsection,more examples are tested to verify the correction strategy numerically.Considering it is time-consuming to do the Monto-Carlo simulation because of the DE involved in solvingProblem 3,the continuation method is adopted instead.The parameter of the continuation is set to be the chaser’s initial orbital states,namelya,e,i,Ω andθ,whereas the target remains unchanged.The starting point of the continuation is the same with the chaser’s initial states stated in Example 1.Figs.9 and 10 show the minimum rendezvous time obtained,and the intervals between the discrete points (i.e.,the continuation step) are set as Δa=2 km,Δe=2×10?5,Δi=0.02°,ΔΩ=0.02°,Δθ=1°seperately.

Figure 9:Continuation on a

Figure 10:Continuation on e,i,Ω,θ

Obviously,when the initial states of the chaser changed,the resulting difference of both orbital plane and orbital size lead to a different minimum rendezvous time.The red dots represent the leading case,whereas the blue represent the trailing case.The cross point between the blue line and red line in Fig.9 indicates that at this specific point the minimum rendezvous time of both cases meet a same result,whereas the transfer orbits are totally different.In addition,compared to other elements,the change of eccentricity (nearly zero) has little effect on the minimum rendezvous time.

The continuation generates 750 examples.The flight time correction and correction factors for all these examples are summarized in Fig.11,showing that the resulting changes of the flight timeδtare less than 0.2 day,which is no more than 3 revolutions compared to the total revolutions up to hundreds.Tab.8 represents the statistics ofη?,demonstrating that the corrections made are relatively slight.In addition,as shown in Fig.12,the terminal orbital errors in position and velocity can be reduced to the order of 1 km and 1 m/s simultaneously,which are acceptable by the relative navigation system.Theη-correction strategy,which is empirical but simple,demonstrates its effectiveness to generate the control laws that precisely steer the chaser to the target.

Figure 11:Flight time correction and correction factors during the continuation

Table 8:Statistics of correction factors η?

Figure 12:Terminal errors before and after correction

7 Conclusions

Considering the prospect of the solar electric propulsion spacecraft for debris removal missions,a practical optimization method is proposed to compute the optimal low-thrust minimumtime many-revolution orbital rendezvous in sun-synchronous orbits,which usually consists of a large number of ballistic orbital arcs in Earth shadow.Different from the existing works,the prominent advantage of the proposed method are twofold:(1) The orbital averaging technique is used with a new slowing-changing variable (the mean longitude difference) proposed to accommodate to the rendezvous problem,which makes it easily to model the perturbations and thrust discontinuity in Earth shadow,and the computational complexity is greatly simplified in the aspect of numerical iterations for trajectory optimization.(2) By using a simple correction strategy,the resulting optimal control derived from the mean orbital dynamics and the corresponding costate variables is mapped to the control law that is used in osculating orbital dynamics to guarantee the rendezvous describing in osculating elements.The effectiveness of the proposed optimization method was successfully demonstrated by two typical low-thrust-orbit-rendezvous missions and hundreds of running with the continuation of the initial states.

Funding Statement:This work was supported by the National Key Research and Development Project (Grant No.2018YFB1900605) and the Key Research Program of Chinese Academy of Sciences (Grant No.ZDRW-KT-2019-1).

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.

Appendix A:Deduction of Optimal Control u?Substituting Eq.(23) into Eq.(26),the Hamiltonian function can be rewritten as

where Λ(uc) is the specific items related to the active controluc.The expressions of Λ(uc) andΥare

were

It is observed that

Thus,to minimize(uc),Γ should be minimized,resulting to that the optimal unit direction vectorα?is opposite to the vectorsatisfies

Substituting the optimal thrust directionα?into Eq.(28),it is readily to obtain≤0.Meanwhile,the terminal costate of the mass satisfies the transversality condition=0.Thus,is a non-negative real number during the entire transfer phase.In addition,substitutingα?into Γ,the coefficient beforeTis rewritten as

which is a non-positive real number.As a result,the optimal thrust amplitude should be maximum all the time,i.e.,: