CN102424119A - Interplanetary low-thrust transfer orbit design method based on polynomial approximation - Google Patents

Abstract

The invention relates to an interplanetary low-thrust transfer orbit design method based on polynomial approximation. Firstly, variable initial value guess is designed, the initial value guess of the transfer orbit design is given, then, the beginning and tail end boundary conditions of a detector are calculated, next, a second type Chebyshev polynomial is adopted for fitting the transfer orbit of the detector, the performance indexes and the constrained conditions are calculated, finally, whether the feasibility conditions are met or not is judged according to the calculated thrust restriction in accordance with whether the calculated performance indexes meet the optimum conditions or not, if all requirements are met, the optimization is successful, the optimum transfer orbit is obtained, and if one requirement is not met, the initial value guess of the design variable in the first step is regulated until the optimization is successful. The method utilizes the Chebyshev polynomial for approximating the low-thrust transfer orbit shape, the time is used as the independent variable, and the flight time restriction is avoided. The polynomial factor is determined through the beginning and tail end orbit state restriction of the detector, and the method has the advantage that the fast design on the low-thrust transfer orbits in different task types can be realized according to the given beginning and tail end boundary conditions.

Description

The interplanetary low thrust transfer orbit method of designing that approaches based on multinomial

Technical field

The present invention relates to a kind ofly, be applicable to the initial designs of interplanetary low thrust transfer orbit, belong to survey of deep space transfer orbit technical field based on polynomial interplanetary low thrust transfer orbit Fast design method.

Background technology

Be expert in the interstellar probe mission, how much fuel detector need consume could arrive the problem that the target celestial body is the task design first concern.Compare the traditional chemical propulsion system, thrustor has than leaps high, the characteristics of light weight, effectively fuel saving consumption when interplanetary probe utilizes it to realize the track transfer.Because need working long hours, thrustor just can reach the purpose that changes track; Make that the detector track is the non-Keplerian orbit of representative type strong nonlinearity; Many theories and method in the tradition pulse track design are no longer suitable, seek the focus that a kind of fast effective transfer orbit method of designing becomes present research.Based on the star-like curve approach method is to realize that at present transfer orbit designs valid approach the most fast; It is to use for reference the thought that conic section is described the pulse track; Utilize suitable function curve to describe the low-thrust trajectory shape; Obtain the parameters relationship between star-like curve and low thrust transfer orbit, and then effectively reduce the transfer orbit design difficulty.Choose which kind of star-like curve low-thrust trajectory is approached, how to resolve the low-thrust trajectory parameter on this basis, determined detector transfer orbit design correctness and design efficiency.Therefore the low thrust transfer orbit method of designing based on star-like curve is one of emphasis problem of current scientific and technical personnel's concern.

In the low thrust transfer orbit method of designing that has developed based on star-like curve; Formerly technological [1] (Petropoulos A E; Longuski J M.Shape-based algorithm for automated design of low-thrust, gravity-assist trajectories [J] .Journal of Spacecraft and Rockets, 2004; 41 (5): 87-796.); Adopt sinusoidal exponential curve that low-thrust trajectory is approached, this method is at first through supposing that thrust direction all the time along detector speed direction or reversing sense, has obtained the isoparametric analytical expression of thrustor acceleration/accel in the transfer process; Through going through, obtain to satisfy the optimal transfer orbit trajectory of task constraint then all over the sinusoidal exponential function parameter of search.Because sinusoidal exponential function free parameter quantity is few, the transfer orbit of trying to achieve can only satisfy the position constraint of detector end at the whole story, therefore can only be used to leap the design of type track; In addition, what this method adopted is to go through all over way of search design track, has shortcomings such as calculated amount is big, design efficiency is low.

Formerly technology [2] is (referring to Wall B.J.and Conway B.A.Shape-Based Approach to Low-ThrustRendezvous Trajectory Design [J] .Journal of Guidance; Control; And Dynamics, 2009,32 (1): 95-101.); Adopt contrary six order polynomials that low-thrust trajectory is approached; This method has equally also used thrust direction all the time along detector speed direction or reciprocal hypothesis, has tried to achieve the analytic solution of low-thrust trajectory parameter, adopts genetic algorithm to obtain best polynomial parameters then.Because contrary six order polynomials have seven free parameters, can satisfy the constraint of detector end position at whole story speed and transfer time simultaneously, so this method can be effective to the design of intersection type transfer orbit.But this method still can't break through the thrust direction hypothesis; And because it is to be the star-like curve independent variable to shift phase angle; Select under the uncomfortable situation at free parameter, the flight time constraint is difficult to satisfy, and this has directly influenced transfer orbit design-calculated robustness and efficient.

Summary of the invention

The present invention is directed to and have problems such as thrust direction hypothesis, poor robustness, efficient are low in the interplanetary low thrust transfer orbit design at present based on star-like curve; Provided a kind of interplanetary low thrust transfer orbit method of designing that approaches based on multinomial; Directly avoid the flight time constraint, improved robustness.

Be somebody's turn to do the interplanetary low thrust transfer orbit method of designing that approaches based on multinomial, may further comprise the steps:

The first step: the conjecture of design variable initial value, the initial value conjecture of given transfer orbit design variable;

Second step: hold boundary condition the whole story of calculating detector: through reading planet ephemeris file, according to preset t ₀Constantly obtain setting out celestial body day heart cartesian coordinate system position vector r _LWith velocity vector v _L, according to preset t _fConstantly obtain the position vector r of target celestial body _AWith velocity vector v _A, hold boundary condition to do the whole story that obtains detector $\begin{array}{c}{r}_0=r_L\\ v_0=v_L+V_{L\∞}\end{array},\begin{array}{c}{r}_f=r_A\\ v_f=v_A+V_{A\∞}\end{array};$

R wherein ₀Heliocentric place vector when launching for detector, v ₀Heliocentric velocity vector when launching for detector, V _{L ∞}The heliocentric velocity vector of the relative earth when launching for detector; r _fHeliocentric place vector when arriving the target celestial body for detector, v _fHeliocentric velocity vector during for arrival target celestial body, V _{A ∞}The heliocentric velocity vector of relative target celestial body when arriving the target celestial body for detector; Then the boundary condition in the cartesian coordinate system is transformed in the spherical coordinate system, and the phase angle is revised:

$Y_0^s=C_{\mathrm{sc}}^{-1}\left(t_0\right){[r_0^T,v_0^T]}^T$

$Y_f^s=C_{\mathrm{sc}}^{-1}\left(t_f\right){[r_f^T,v_f^T]}^T$

track condition in day bulbus cordis system of axes when launching wherein for detector, the track condition when

arrives the target celestial body for detector in day bulbus cordis system of axes;

The 3rd step: adopt the transfer orbit of Chebyshev polynomial of the second kind match detector, polynomial matrix can be expressed as

Wherein τ is normalized time variable, t ₀The moment corresponding τ=-1, t ₀+ t _fThe moment corresponding τ=1; Utilize the detector under the spherical coordinate system to hold boundary condition to calculate the Chebyshev polynomials coefficient whole story then;

The 4th step: calculation of performance indicators constraints:, calculate the performance figure and the constraint of thrustor thrust of transfer orbit according to the Chebyshev polynomials coefficient that obtains;

The 5th step: the performance figure according to calculating judge whether to satisfy optimality condition; Thrust constraint according to calculating judges whether to satisfy feasibility condition; If all satisfy, then optimize successfully, obtain optimal transfer orbit trajectory; If there is one not satisfy, the initial value conjecture of then adjusting design variable in the first step is until optimizing successfully.

The inventive method utilizes Chebyshev polynomials to approach low thrust transfer orbit shape, is that independent variable has been avoided the flight time constraint with time, has strengthened the algorithm robustness; Hold the track condition constraint to confirm multinomial coefficient the whole story through detector then, and directly obtained the acceleration/accel analytical expression in the transfer process, avoided the thrust direction hypothesis, improved computational efficiency; At last the transfer orbit design problem is converted into the mixed integer nonlinear programming problem that only has simple parameter, effectively reduces the track design difficulty.This method can design according to the low thrust transfer orbit of holding the given whole story boundary condition to the different task type fast.

Description of drawings

Fig. 1 is an interplanetary probe transfer orbit scheme drawing.

The specific embodiment

This interplanetary low thrust transfer orbit Fast design method is divided and is adopted iteration optimizing account form to find the solution, and solution procedure is divided into the conjecture of design variable initial value, the whole story and holds that boundary condition calculates, multinomial coefficient calculates, the performance figure constraints is calculated, five parts of design variable adjustment.

1) design variable initial value conjecture

The initial value conjecture Z of given transfer orbit design variable ₀

2) hold the whole story boundary condition to calculate

Through reading planet ephemeris file, according to t ₀Obtain setting out celestial body day heart cartesian coordinate system position and speed vector r _LAnd v _L, according to t _fObtain the position and speed vector r of target celestial body _AAnd v _A, hold boundary condition to do the whole story that then can obtain detector

r ₀＝r _L

v ₀＝v _L+V _L∞

r _f＝r _A

v _f＝v _A+V _A∞

Boundary condition in the cartesian coordinate system is transformed in the spherical coordinate system does

$Y_0^s=C_{\mathrm{sc}}^{-1}\left(t_f\right){[r_0^T,v_0^T]}^T$

$Y_f^s=C_{\mathrm{sc}}^{-1}\left(t_f\right){[r_f^T,v_f^T]}^T$

Recomputate detector azimuth endways then and carry out the phase place correction

θ _f＝θ _f+2πP

3) multinomial coefficient calculates

Make τ=-1 and τ=1 respectively, utilize computes Chebyshev polynomial of the second kind matrix M (1) and M (1)

Utilize the detector under the spherical coordinate system to hold boundary condition to calculate the Chebyshev polynomials coefficient whole story then

Wherein s is

4) the performance figure constraints is calculated

According to the Chebyshev polynomials coefficient that obtains, calculate the performance figure of transfer orbit

$J={\∫}_{t_0}^{t_f}\left|\right|A\left(t\right)\left|\right|\mathrm{dt}$

Calculate the constraint of thrustor thrust

Φ＝max(||A(t)||)-A _max

5) design variable adjustment

Optimizing algorithm judges whether to satisfy optimality condition based on the performance indications of calculating, and judges whether to satisfy feasibility condition based on the thrust constraint of calculating.If all satisfy, then optimize successfully, obtain optimal transfer orbit trajectory, if there is one not satisfy, then adjust design variable, and repeating step 2)～5), until optimizing successfully.

The transfer orbit location status X of interplanetary probe can describe in day bulbus cordis system of axes that with the ecliptic plane is reference, and it is the function of flight time t, can be expressed as

R wherein; θ,

is respectively central meridian distance (CMD), azimuth and the elevation angle of detector in spherical coordinate system.

Thrustor thrust is very little, needs long-time acceleration just can reach the purpose that changes track, and this makes interplanetary transfer orbit be generally the multi-circle spiral flow wire shaped.Chebyshev polynomials have good characteristic aspect curve fitting, the Chebyshev polynomials of choosing the flight time and be independent variable are approached three location statuss respectively and obtained

$r\left(t\right)=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}{a}_n^r{U}_n\left(\mathrm{\τ}\right)$

$\mathrm{\θ}\left(t\right)=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}{a}_n^{\mathrm{\θ}}{U}_n\left(\mathrm{\τ}\right)$

Wherein N is the Chebyshev polynomials number of times, U _n(τ) be n rank Chebyshev polynomial of the second kind,

Be respectively the cooresponding Chebyshev polynomials coefficient in central meridian distance (CMD), azimuth and the elevation angle, τ ∈ [1,1] is normalized transfer time, and the relation of itself and actual flying time does

$\mathrm{\τ}=\frac{2(t-t_0)}{t_f-t_0}-1$

T wherein ₀Be initial time, t _fBe the end moment.Ask first derivative can obtain being expressed as of detector speed all directions component to above-mentioned three Chebyshev's approximating polynomials

Wherein the recurrence formula of n rank Chebyshev polynomial of the second kind first derivative

does

And have

By position of detector with speed Chebyshev is approximant can find out; If multinomial coefficient

confirms that then the track condition of detector is confirmed.When detector shifts between planet; Track must satisfy constraints such as the ephemeris constraint, flight time of planet; Since Chebyshev polynomials are independent variable with time, the flight time constraint is satisfied, so multinomial coefficient can only just can obtain through the ephemeris constraint.Assumed that the track detector state at the initial time constraint

orbital moment in the terminal state constraint is Order

Then M (τ) is 2 * (N+1) matrix, will hold the whole story track condition constraint to bring Chebyshev respectively into and can obtain in approximant

Can find out that by above-mentioned three formulas there is unique solution in set of equations when N=3, this is illustrated in and satisfies detector and hold the whole story under the condition of track condition constraint, and three times Chebyshev polynomials just can realize approaching of low thrust transfer orbit.Equation found the solution to obtain Chebyshev coefficient and do

Wherein s is r, θ or

For the low thrust acceleration/accel of calculating detector in transfer process, detector spherical coordinate system lower railway state exchange is arrived cartesian coordinate system

Y ^c(t)＝C _sc(t)Y ^s(t)

Y wherein ^c=[r _x, r _y, r _z, v _x, v _y, v _z] ^T, C _ScBe the transition matrix of spherical coordinate system to cartesian coordinate system, then the thrustor acceleration/accel can be expressed as

Wherein μ is the gravitational acceleration constant of the sun.

Can find out by following formula; Acceleration magnitude and the direction of thrustor in transfer process determined by the Chebyshev polynomials coefficient; The acceleration/accel direction need not to suppose in advance to limit; And acceleration magnitude will receive the restriction of maximum thrust usually, must consider the thrust size constraint when therefore transfer orbit being designed

max(||A(t)||)≤A _max

A wherein _MaxThe maximum thrust that the thrustor that adopts for task allows.

In sum, the variable that influences the Chebyshev polynomials coefficient is t launch time of detector ₀With hyperbola hypervelocity V _{L ∞}, the time t of arrival target planet _fWith hyperbola hypervelocity V _{A ∞}, the whole number of turns P of detector low thrust transfer orbit, then the track design problem can be summed up as

Design parameters

Z＝[t ₀，t _f，V _L∞，V _A∞，P] ^T

Performance figure

$J\left(Z\right)={\∫}_{t_0}^{t_f}\left|\right|A\left(t\right)\left|\right|\mathrm{dt}\→\min$

Constraint condition

Φ(Z)＝max(||A(t)||)-A _max≤0

The problem of more than summing up is simple mixed integer nonlinear programming problem, can adopt the differential evolution scheduling algorithm directly to find the solution, the promptly corresponding best low thrust transfer orbit that satisfies the task constraint of optimal solution.

Claims (1)

1. the interplanetary low thrust transfer orbit method of designing that approaches based on multinomial is characterized in that may further comprise the steps:

The first step: the conjecture of design variable initial value, the initial value conjecture of given transfer orbit design variable;

Second step: hold boundary condition the whole story of calculating detector: through reading planet ephemeris file, according to preset t ₀Constantly obtain setting out celestial body day heart cartesian coordinate system position vector r _LWith velocity vector v _L, according to preset t _fConstantly obtain the position vector r of target celestial body _AWith velocity vector v _A, hold boundary condition to do the whole story that obtains detector $\begin{array}{c}{r}_0=r_L\\ v_0=v_L+V_{L\∞}\end{array},\begin{array}{c}{r}_f=r_A\\ v_f=v_A+V_{A\∞}\end{array};$

R wherein ₀Heliocentric place vector when launching for detector, v ₀Heliocentric velocity vector when launching for detector, V _{L ∞}The heliocentric velocity vector of the relative earth when launching for detector; r _fHeliocentric place vector when arriving the target celestial body for detector, v _fHeliocentric velocity vector during for arrival target celestial body, V _{A ∞}The heliocentric velocity vector of relative target celestial body when arriving the target celestial body for detector; Then the boundary condition in the cartesian coordinate system is transformed in the spherical coordinate system, and the phase angle is revised:

$Y_0^s=C_{\mathrm{sc}}^{-1}\left(t_0\right){[r_0^T,v_0^T]}^T$

$Y_f^s=C_{\mathrm{sc}}^{-1}\left(t_f\right){[r_f^T,v_f^T]}^T$

track condition in day bulbus cordis system of axes when launching wherein for detector, the track condition when

arrives the target celestial body for detector in day bulbus cordis system of axes;

The 3rd step: adopt the transfer orbit of Chebyshev polynomial of the second kind match detector, polynomial matrix can be expressed as

Wherein τ is normalized time variable, t ₀The moment corresponding τ=-1, t ₀+ t _fThe moment corresponding τ=1; Utilize the detector under the spherical coordinate system to hold boundary condition to calculate the Chebyshev polynomials coefficient whole story then;

The 4th step: calculation of performance indicators constraints:, calculate the performance figure and the constraint of thrustor thrust of transfer orbit according to the Chebyshev polynomials coefficient that obtains;

The 5th step: the performance figure according to calculating judge whether to satisfy optimality condition; Thrust constraint according to calculating judges whether to satisfy feasibility condition; If all satisfy, then optimize successfully, obtain optimal transfer orbit trajectory; If there is one not satisfy, the initial value conjecture of then adjusting design variable in the first step is until optimizing successfully.

Images