CN111268176A - Perturbation track four-pulse intersection rapid optimization method - Google Patents

Abstract

The invention discloses a perturbation track four-pulse intersection rapid optimization method, which comprises the following steps: according to known conditions of the spacecraft and the target, obtaining control quantity and minimum speed increment through a control model taking elimination of orbital root difference at the intersection time as an optimization result; establishing a four-pulse rendezvous optimization model according to the control quantity and the speed increment, taking values of the time and the size of the first two pulses in the neighborhood according to the estimated value of the step, and recalculating the last two pulses according to the actual orbit recursion result; and replacing the analytic perturbation kinetic model with a full perturbation accurate kinetic model, taking the optimal solution of the analytic kinetic model as an initial value, and solving again to obtain the optimal solution of the four-pulse intersection under the accurate kinetic model. The problem of low computing speed and the like in the prior art is solved, and computing efficiency is improved.

Description

Perturbation track four-pulse intersection rapid optimization method

Technical Field

The invention relates to the technical field of aerospace navigation control, in particular to a perturbation orbit four-pulse intersection rapid optimization method.

Background

On the near-earth orbit, the earth non-spherical gravitational perturbation, the atmospheric resistance and the like are the most main perturbation items. When pulse-cross orbit optimization is carried out, the influence of a perturbation term must be considered, but the corresponding orbit integral calculation amount is obviously increased compared with a two-body analysis dynamic model. In addition, the existing method for solving the optimal speed increment of the rail intersection lacks reasonable initial value setting, the solution space range needing to be searched when an evolutionary algorithm is used is large, the calculation time is long, and the convergence is slow particularly when the intersection transfer time is long. Most of the existing algorithms use analytic perturbation models, but the optimal calculation speed is still not ideal and has a certain difference with an accurate dynamic model.

Disclosure of Invention

The invention provides a perturbation orbit four-pulse intersection rapid optimization method which is used for overcoming the defects of long calculation time, slow convergence and the like in the prior art, realizing the accuracy of an iterative model and improving the calculation efficiency.

In order to achieve the purpose, the invention provides a perturbation track four-pulse intersection rapid optimization method, which comprises the following steps:

step 1, according to known conditions of a spacecraft and a target, obtaining control quantity and minimum speed increment through a control model taking elimination of orbital root difference at a meeting time as an optimization result;

changing the number of orbits when the spacecraft is controlled to run in the first circle of orbits according to the pulse which converts a part of control quantity and speed increment into two fixed periods at intervals; converting the residual control quantity and the speed increment into two other pulse control spacecrafts with fixed intervals to finish intersection when the last circle of the orbit runs; estimating the size of four pulses and components in three directions of radial direction, tangential direction and normal direction;

step 2, establishing a four-pulse rendezvous optimization model, taking values of the time and the size of the first two pulses in the neighborhood according to the estimated value in the step 1, and recalculating the last two pulses according to the actual orbit recursion result;

and 3, replacing the analytic perturbation kinetic model with a full perturbation accurate kinetic model, and solving the step 2 again by taking the optimal solution of the analytic kinetic model as an initial value to obtain the optimal solution of the four-pulse intersection under the accurate kinetic model.

The perturbation orbit four-pulse intersection rapid optimization method provided by the invention is suitable for the optimization problem of the orbit intersection with large ascension difference at the ascending intersection point and more similar numbers of other orbits. The method comprises the steps of firstly, providing a reasonable initial value for the four-pulse rendezvous optimization by adopting an estimation method according to the difference value of the number of orbits of a spacecraft and a rendezvous target and fixed transfer time, then designing an iterative model for quickly obtaining an accurate solution, obtaining the optimal solution of the four-pulse rendezvous through local search, and obtaining the accurate optimal solution of the four-pulse rendezvous by using less optimization variables; the method solves the problems of large kinetic recursion calculated amount and slow convergence when the transfer time is long in the traditional multi-pulse rendezvous optimization algorithm.

Drawings

In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly described below, it is obvious that the drawings in the following description are only some embodiments of the present invention, and for those skilled in the art, other drawings can be obtained according to the structures shown in the drawings without creative efforts.

FIG. 1 is a flow chart of a perturbation orbit four-pulse intersection fast optimization method proposed by the invention.

The implementation, functional features and advantages of the objects of the present invention will be further explained with reference to the accompanying drawings.

Detailed Description

The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.

It should be noted that all the directional indicators (such as up, down, left, right, front, and rear … …) in the embodiment of the present invention are only used to explain the relative position relationship between the components, the movement situation, etc. in a specific posture (as shown in the drawing), and if the specific posture is changed, the directional indicator is changed accordingly.

In addition, the descriptions related to "first", "second", etc. in the present invention are only for descriptive purposes and are not to be construed as indicating or implying relative importance or implicitly indicating the number of technical features indicated. Thus, a feature defined as "first" or "second" may explicitly or implicitly include at least one such feature. In the description of the present invention, "a plurality" means at least two, e.g., two, three, etc., unless specifically limited otherwise.

In the present invention, unless otherwise expressly stated or limited, the terms "connected," "secured," and the like are to be construed broadly, and for example, "secured" may be a fixed connection, a removable connection, or an integral part; the connection can be mechanical connection, electrical connection, physical connection or wireless communication connection; they may be directly connected or indirectly connected through intervening media, or they may be connected internally or in any other suitable relationship, unless expressly stated otherwise. The specific meanings of the above terms in the present invention can be understood by those skilled in the art according to specific situations.

In addition, the technical solutions in the embodiments of the present invention may be combined with each other, but it must be based on the realization of those skilled in the art, and when the technical solutions are contradictory or cannot be realized, such a combination of technical solutions should not be considered to exist, and is not within the protection scope of the present invention.

Example one

As shown in fig. 1, an embodiment of the present invention provides a method for rapidly optimizing perturbation orbit four-pulse intersection, which is characterized by comprising the following steps:

step 1, according to known conditions of a spacecraft and a target, obtaining control quantity and minimum speed increment through a control model taking elimination of orbital root difference at a meeting time as an optimization result;

changing the number of orbits when the spacecraft is controlled to run in the first circle of orbits according to the pulse which converts a part of control quantity and speed increment into two fixed periods at intervals; converting the residual control quantity and the speed increment into two other pulse control spacecrafts with fixed intervals to finish intersection when the last circle of the orbit runs; estimating the size of four pulses and components in three directions of radial direction, tangential direction and normal direction;

given that the spacecraft and the target are at the transfer start time t₀In the case of the number of tracks of (1), the calculation is at t_fTime-of-day intersection, transferring an optimal four-pulse (including the time and magnitude of application of each pulse) of duration Δ t, said step 1 comprising:

step 11, respectively calculating the orbit root difference of the spacecraft and the target which need to be eliminated by the pulse at the intersection moment according to the known conditions and the J2 analytic dynamic model;

step 12, designing an optimal control model for eliminating the track root number difference, and obtaining control quantity and minimum speed increment required by intersection after pulse elimination;

step 13, taking the solution obtained in the step 12 as an initial value, and calculating the minimum speed increment required by intersection under a high-order approximation model;

step 14, correcting the minimum speed increment obtained in the step 13 according to the phase difference to obtain the minimum speed increment required by intersection after the phase difference is corrected;

and step 15, correcting the minimum speed increment obtained in the step 14 according to the eccentricity vector difference, and obtaining the estimated size of each pulse and components in three directions, namely radial direction, tangential direction and normal direction.

Step 2, establishing a four-pulse rendezvous optimization model, taking values of the time and the size of the first two pulses in the neighborhood according to the estimated value in the step 1, and recalculating the last two pulses according to the actual orbit recursion result;

and 3, replacing the analytic perturbation kinetic model with a full perturbation accurate kinetic model, and solving the step 2 again by taking the optimal solution of the analytic kinetic model as an initial value to obtain the optimal solution of the four-pulse intersection under the accurate kinetic model.

The known spacecraft and the target are at a starting time t₀The number of tracks. Let t_fFor the meeting time, the transition duration from the transition start time to the meeting time is Δ t. The number of the tracks is respectively a semi-major axis a of the track, an eccentricity e of the track, an inclination angle i of the track, a right ascension channel omega of a rising intersection point, an angle omega of a near center point and an angle M of a horizontal near point.

Respectively calculating the intersection time t according to the J2 analytic dynamic model_fThe number of spacecraft orbits and the number of target orbits. Wherein, the J2 analytic kinetic model is as follows:

wherein:

the derivative of the number of tracks with respect to time, J₂Being a second order term of Earth's gravity, R_eIs the radius of the earth, p ═ a (1-e)²) The semi-diameter of the track is the track semi-diameter,

mu is the earth's gravitational constant, orbital angular velocity.

The orbit root difference [ Delta a ] of the spacecraft and the target at the moment of intersection, which needs to be eliminated by the pulse₀,Δe₀,Δi₀,ΔΩ₀,Δω₀,ΔM₀]Is t_fSubtracting the orbit number of the spacecraft from the target orbit number at the moment.

Elimination of Delta omega in step 12 of the invention₀The optimal control model is as follows: applying control to the spacecraft at the transfer starting moment to enable the number of the orbits of the spacecraft to generate variation delta a, delta i and delta omega ', wherein the delta a is the variation of the semi-major axis of the orbits of the spacecraft at the transfer starting moment, the delta i is the variation of the inclination angle of the orbits of the spacecraft at the transfer starting moment, and the delta omega' is the variation of the right ascension channel of the ascending nodes of the spacecraft at the transfer starting moment; applying control quantities with the same reverse direction to the spacecraft at the meeting moment; Δ a, Δ i, Δ Ω' are required to satisfy the condition that the right ascension of the spacecraft at the intersection relative to the target at the intersection is zeroThe beam, under this constraint, solves for the optimal Δ a, Δ i, Δ Ω' distribution in the control quantities, making the speed increment take a minimum value. The specific implementation method comprises the following steps:

for the semi-major axis a and the ascent point right ascension Ω of the spacecraft orbit, according to the J2 analytic dynamic model, that is, equation (1), when the orbit inclination angle of the spacecraft and the rendezvous target is close and the orbit eccentricity is close to 0, it is considered that:

the control quantity delta i and delta a of the orbit inclination angle and the orbit semi-major axis of the spacecraft can cause the right ascension at the lifting point of the spacecraft

Change of (2):

the pulse application mode is as follows: at the transfer start time t_fApplication of a single pulse Δ v at₀Changes the semi-major axis and the orbit inclination angle of the orbit of the spacecraft, and further changes the ascension point right ascension channel of the spacecraft

Then at the meeting time t_fApplying a single pulse Deltav_fRestoring the semi-major axis and the orbit inclination angle of the spacecraft orbit to the initial values when the control quantity is not applied, and then t_fThe right ascension difference of the intersection target relative to the spacecraft at the moment is as follows:

the speed increment required for eliminating delta omega' is divided into two and synthesized into double pulses, and then in order to eliminate delta omega₀The required speed increments for the spacecraft are:

Δ v in formula (5)_a，Δv_i，Δv_ΩRespectively the speed increment corresponding to the variation delta a of the semi-major axis of the orbit of the spacecraft, the variation delta i of the inclination angle of the orbit and the variation delta omega' of the right ascension of the ascending intersection point,

the orbit average velocity is the initial moment of the spacecraft and sini ≈ sin (i + Δ i) needs to be assumed.

The estimation problem of the optimal speed increment is converted into the solution of Δ v_sumI.e. Δ v₀The extreme value problem of (2). The estimation problem for the optimal speed increment can be described as:

will be provided with

Is calculated as a first order approximation, then

Wherein

The method for calculating the red channel drift rate of the lifting intersection point at the initial time of the spacecraft is obtained by only retaining a first order term according to a J2 analytic dynamic model, namely formula (1), in linear expansion of the number of orbits of the spacecraft at the initial time.

Considering the normalized variables x Δ a/a, y Δ i, z Δ Ω', f and g are equivalent to:

defining L ═ f + λ g, where λ is a constraint multiplier, then the extremum condition is:

equation (9) is a fourth order linear equation that can be directly solved:

only Δ Ω is eliminated₀The minimum speed increments required are:

further, the implementation method of step 13 of the present invention is as follows:

to at t_fTime of day elimination of Δ a₀And Δ i₀Δ v in formula (5)_fThe rewrite is:

when Δ a and Δ i are large, the calculation is performed by equation (7)

Is large, the calculation of formula (3) is adopted

And the semi-major axis of the spacecraft orbit can not be less than a in the process of meeting maneuver_minConstraint (a)_minTaken as 6600km), a more accurate solution model is obtained as:

wherein a is_minIs the minimum spacecraft orbit semi-major axis of constraint.

Likewise, let x be Δ a/a_cY ═ Δ i, z ═ Δ Ω', and L ═ f + λ g₁+χg₂Wherein λ and χ are constraint multipliers, respectively, and the corresponding extreme conditions are:

wherein

Calculated when a is replaced by a + ax and i is replaced by i + y in equation (2)

Namely, it is

At the moment, the equation set (15) is nonlinear, x, y, z and lambda solved in the formula (11) in the step 12 are used as initial values of x, y, z and lambda in the equation set (15), the initial value of chi in the equation set (15) is set to be 0, multi-step iterative correction is carried out, new x, y, z, lambda and chi are obtained, and then new x, y, z, lambda and chi are obtained

Further, the implementation method of step 14 of the present invention is as follows:

let Δ λ₀＝Δω₀+ΔM₀The spacecraft and the target are at t without control_fA phase difference of time, applying Δ v₀After, t_fThe time phase difference is changed to:

normalizing Δ λ to [ - π, π]To satisfy the convergence constraint, Δ a and Δ i need to be modified to produce a relative drift rate of the phase

Considering sin²(i + Δ i) is insensitive to variations in i, changes in Δ a are considered dominant; the correction amount is:

superimpose it on Δ v₀I.e. fixing x in formula (5)_fix＝(Δa+Δa_mod) A, the system of equations translates to:

using the new y, z, λ obtained by iterative correction in step 13 as the initial values of y, z, λ in equation (20), re-solving Δ i and Δ Ω' to obtain the minimum speed increment required for intersection after correcting the phase difference

Further, the implementation method of step 15 of the present invention is as follows:

applying Δ v to a spacecraft at an initial time₀Later, the drift rate of the spacecraft ω changes, corresponding to the crossing time t_fThe new ω is:

t_fthe difference of the orbit eccentricity of the spacecraft and the target at the moment is

According to Δ v₀、Δv_fThe normal component of the two pulses to the change quantity of the orbit inclination angle i and the ascension point right ascension omega of the spacecraft, and the control phase of the two pulses is known to be

Wherein at u₀To apply Δ v₀The variation of the eccentricity of the spacecraft orbit rotates along with the drift of omega, and t is considered to be_fThe phase for this change in time is:

if Δ v is to be reduced₀，Δv_fSplit into two pulses (defined as av) separated by half a period₁、Δv₂Equivalent to Δ v₀，Δv₃、Δv₄Equivalent to Δ v_f) The distribution proportion is different according to the pulse size, and the change quantity delta e of the eccentricity ratio of the spacecraft orbit is_v0And Δ e_vfComprises the following steps:

-Δa/a＜|Δe_v0|＜Δa/a

-(Δa-Δa₀)/a＜|Δe_vf|＜(Δa-Δa₀)/a (25)

let k₀、k_f∈[0,1]Are respectively Δ v₀，Δv_fThe delta v is obtained after the superposition of the variation coefficient of the eccentricity of the spacecraft orbit₀，Δv_fChange quantity delta e to orbit eccentricity of spacecraft_x,Δe_y]Comprises the following steps:

the existing control quantity is used for eliminating the eccentricity difference relative to the target, and the method is equivalent to solving the extreme value problem:

after the solution is obtained in min (f),if Δ e^remA value of 0, indicating that no additional control is required; if Δ e^remValues greater than zero require an additional amount of radial control: Δ v_r＝Δe^remV, Δ V in step (4) according to equation (28)_fIn the middle of superposition, the estimated optimal meeting speed increment is obtained

Wherein c is₁For the radial pulse component Δ v_rAt Δ v₀、Δv_fAccording to the distribution ratio of Δ v obtained in step (4)₀、Δv_fValue calculation

Δv₁,Δv₂,Δv₃,Δv₄The sizes are respectively as follows:

wherein k is₀，k₀See equation (27), where the semimajor axis component corresponds to the tangential component of the pulse, the dip and lift-off right ascension components correspond to the normal component, and the remainder corresponds to the radial component. By the formula (23), Δ v₁,Δv₂,Δv₃,Δv₄The phase angles of (a) are respectively:

wherein u is₁And u₂First turn after task start, u₃And u₄At the last turn.

Preferably, the step 2 optimizes the four-pulse intersection optimization model by using a differential evolution algorithm.

Preferably, the step 2 includes:

step 21, directly applying the four-pulse estimation value obtained in the step 1 to a spacecraft, recurrently delivering the orbit to a known rendezvous moment, and defining a local search range coefficient of a speed increment and a local search range coefficient of a control quantity according to the deviation between the recurrently delivering orbit and a target orbit; the step 21 includes:

in step 211, the four-pulse estimation values obtained in step 1 are respectively: Δ v₁,Δv₂,Δv₃,Δv₄The phase angles are respectively: u. of₁,u₂,u₃,u₄Wherein u is₁,u₂For the first turn after the start of the rendezvous task, u₃,u₄The last circle after the rendezvous task starts;

at step 212, four pulses are applied directly to the spacecraft and the orbit is recurred until time t_fDeviation from the target trajectory, defining 4 coefficients [ x ]₁x₂x₃x₄]In a local search, where x₁,x₂,x₃Local search range coefficients, x, representing the tangential, normal and radial components of the velocity increment, respectively₄A local search range coefficient representing the corresponding eccentricity control quantity in the step 1; equations (29) and (30), i.e. the estimated values of the four pulses, are applied directly to the spacecraft and the orbit is recurred to time t_fThere will be some small deviation from the target trajectory. The invention considers that the difference between the estimated value and the actual optimal solution is small, and only local search is needed.

Step 213, convert Δ v₁、Δv₂The changes are as follows:

wherein, delta a is the variation of the semi-major axis of the orbit of the spacecraft at the moment of the transfer initiation, delta i is the variation of the inclination angle of the orbit of the spacecraft at the moment of the transfer initiation, delta omega is the variation of the right ascension channel of the ascending intersection point of the spacecraft at the moment of the transfer initiation, and delta v_rTo require an additional amount of radial control, a₁For a predetermined upper limit of the local change of the semi-major axis, k₀For controlling space flightThe coefficient of the variation of the eccentricity of the orbit, i is the inclination angle of the orbit,

the average orbit velocity of the spacecraft.

Step 22, determining the moments of four pulses through an optimization model established by the time of pulse intersection, the phase of the spacecraft and the orbital angular velocity when the spacecraft intersects; preferably, said step 22 comprises:

step 221, defining an optimization model, first, the number of orbits of the spacecraft is deduced to the starting moment t of the handover control_f- Δ t, set to [ a_s,e_s,i_s,Ω_s,ω_s,M_s]If the phase at that time is u ═ ω_s+M_sDue to the pulse Δ v₁Is u₁Then the pulse Deltav is known₁Time deltat relative to the start of a session₁Comprises the following steps:

Δt₁＝mod(u₁-u,2π)/n(32)

wherein mod represents the sum of₁-u is normalized to between 0 and 2 π,

for the orbital angular velocity at the beginning of a spacecraft rendezvous, namely, the spacecraft rendezvous needs to wait for delta t₁After that, the first pulse is applied, and Δ v is known from the formula (31)₁The three-dimensional components in the spacecraft orbit coordinate system are:

the number of orbits of the spacecraft is pushed to t_f-Δt+Δt₁At the moment, the number of tracks is converted into a position velocity vector, and delta v is superposed on the velocity₁Obtaining a spacecraft position and speed vector after the first pulse, and then converting to obtain the orbit number;

step 222, due to Δ v₁And Δ v₂Phase difference of u₂-u₁Pi, so that the second pulse is separated from the first by a time at₂Comprises the following steps:

Δt₂＝π/n₁(34)

wherein

To add Δ v₁New orbital angular velocity of the spacecraft after, i.e. passing at₂Applying a second pulse of formula (31), Δ v₂The three-dimensional components in the spacecraft orbit coordinate system are:

then the orbit number of the spacecraft after the second pulse is pushed to t_fTime is set as [ a ]_f,e_f,i_f,Ω_f,ω_f,M_f]The number of tracks [ a ] at the same time as the target_t(t_f),e_t(t_f),i_t(t_f),Ω_t(t_f),ω_t(t_f),M_t(t_f)]The difference shows that the right ascension angle Δ Ω at this time_fAnd the difference of inclination angle Δ i_f；

Step 223, according to the principle of equation (30), the phases u of the third and the fourth pulses are determined₃,u₄The change is:

u₃＝atan((Δi_f)/(ΔΩ_fsini/2))

u₄＝atan((Δi_f)/(ΔΩ_fsini/2))+π (36)

the time at which the fourth pulse is applied and t_fThe moments not necessarily coinciding, with a time difference Δ t₄Comprises the following steps:

Δt₄＝-(mod(u₄-(ω_f+M_f),2π)-2π)/n_f(37)

wherein

Track angular velocity for the rendezvous target, (mod (u)₄-(ω_f+M_f),2π)-2π) Indicating handle u₄-(ω_f+M_f) Normalized to between-2 pi and 0, i.e. the fourth pulse has a time t_f-Δt₄The time advance delta t of the third pulse relative to the fourth pulse₃Comprises the following steps:

Δt₃＝π/n₃(38)

wherein n is₃For the angular velocity of the spacecraft orbit after the third pulse is applied, since n₃Is difficult to be directly calculated, and introduces a correction coefficient x to be solved for compensating eccentricity₅Order:

Δt₃＝x₅π/n_f(39)

the timings of the third and fourth pulses are determined by equations (37) and (39). But its size has not yet been determined. Calculating Δ ν using Lambert's algorithm taking perturbation into account₃And Δ v₄The size of (2). The Lambert problem can be described as a given start-stop time instant (t)_inAnd t_out) And a corresponding start-stop position vector (r)_inAnd r_out) Calculating to obtain a start-stop velocity vector (v)_inAnd v_out) The position and speed state of the spacecraft is represented as r_inv_in]Time passes t_out-t_inThe long time orbit recursion can reach the state r_outv_out]. The method has an analytic calculation method for a two-body dynamics model, has more descriptions in open documents, and is regarded as a function in the invention, which is shown in a formula (40).

The magnitude of the last two pulses is calculated by the Lambert algorithm taking perturbation into account, step 23. Preferably, the step 23 includes:

step 231, for a given start-stop time t_in、t_outAnd corresponding start-stop position vector r_in、r_outCalculating a start-stop velocity vector v for a two-body dynamics model_in、v_out：

[v_inv_out]＝Lambert(r_in,r_out,t_in,t_out) (40)

Applying the Lambert algorithm, i.e. setting the starting time t_inIs set to t_f-Δt₄-Δt₃Terminate atTime t_outIs set to t_f-Δt₄At a starting position r_inFor the number of orbits of the spacecraft after the second pulse to be pushed to t_f-Δt₄-Δt₃Position vector obtained at time, end position r_outStepping to t for target trajectory_f-Δt₄A position vector obtained at a time;

step 232, firstly, obtaining the v corresponding to the two-body dynamic model by using an analytic calculation method_inAnd v_outAnd obtaining v corresponding to the J2 analytic perturbation dynamics model by differential correction_inAnd v_out；

Firstly, the v corresponding to the two-body dynamic model is obtained by using a public analytical calculation method_inAnd v_outThen, the differential correction is used to obtain v corresponding to the perturbation dynamics model formula (1)_inAnd v_out；

Step 233, set v₃ ^-For the number of orbits of the spacecraft after the second pulse to be pushed to t_f-Δt₄-Δt₃The velocity vector taken at the moment, the third pulse needs to change the spacecraft velocity to v_inThe meeting location constraint can be satisfied; v. the₄ ⁺Stepping to t for target trajectory_f-Δt₄The velocity vector obtained at the moment, the fourth pulse needs to drive the spacecraft velocity from v_outChange to v₄ ⁺The meeting speed constraint can be met; thus Δ v₃And Δ v₄The sizes are respectively as follows:

the total velocity increment for the four pulses is then:

Δv_sum＝Δv₁+Δv₂+Δv₃+Δv₄(42)

to sum up, the local optimization model includes 5 variables: x is the number of₁,x₂,x₃,x₄,x₅The optimization index is Δ v calculated by the following equations (31) to (42)_sumMinimum; x is the number of₁,x₂,x₃,x₄,x₅The value ranges are as follows:

solving by using a differential evolution algorithm to obtain optimal x₁,x₂,x₃,x₄,x₅And accurate values of the four pulse moments and the four pulse magnitudes under the analytic dynamics model can be obtained.

The differential evolution algorithm can be easily solved, and is a general algorithm for solving the optimization problem, and is not described herein again. Obtaining the optimal x₁,x₂,x₃,x₄,x₅And then, accurate values of the time and the size of the four pulses under the analytic kinetic model can be obtained.

Preferably, the step 3 comprises:

step 31, replacing the J2 analytic perturbation kinetic model of formula (1) with the following accurate kinetic equation, and applying to all processes involving orbit recursion calculation in step 2:

wherein r is the three-dimensional position vector of the spacecraft, v is the three-dimensional velocity vector, a_PerturbationThe acceleration caused by the gravity forces such as the non-spherical gravity of the earth, the atmospheric resistance, the sun-moon gravity, the sunlight pressure and the like is considered; the calculation formula can be obtained in public teaching materials, and is not described in detail herein. It is obvious that the optimal solution of step (6) cannot be directly used as the corresponding optimal solution of the exact kinetic equation (44).

Step 32, using x solved in step 2₂,x₃,x₄Re-optimizing x as an optimal solution for accurate kinetic models₁,x₅；

And step 33, recalculating the optimization models corresponding to the expressions (31) to (42) by using a differential evolution algorithm, so as to obtain the four pulse moments and the four pulse magnitudes with the optimal precise dynamics.

The main difference between the analytic kinetic equation and the accurate kinetic equation is that accumulated phase difference is generated after long-time recursive calculation and can be eliminated through semi-long axis correction. Therefore, in order to further improve the calculation efficiency, x obtained by the solution of the step (6) is used₂,x₃,x₄Re-optimizing only x as an optimal solution for accurate kinetic models₁,x₅And (4) finishing. The number of the optimized variables is reduced to 2, and the solving speed is higher. Recalculating the optimization models (fixed x) corresponding to equations (31) to (42) by applying a differential evolution algorithm₂,x₃,x₄) Then, four pulse time moments and sizes optimal for precise dynamics can be obtained.

Next, the method of the present invention is simulated to obtain the optimal speed increment through a specific embodiment.

In this embodiment: spacecraft and target initial time t₀The number of orbitals is 0, as shown in table 1:

TABLE 1 spacecraft and target initial time orbital radical

	A(km)	e	i (degree)	Omega (degree)	w (degree)	M (degree)
							Spacecraft	6638.14	0.009039	42.05	172.6	120	1
Target	6720.14	0.00001	42	169.2	100	145

To calculate the optimal four-pulse control required for a 14 th day rendezvous starting at the initial time compared to the track transfer. I.e. Δ t equals 14 days, t_fDay 14.

Firstly, calculating the orbit root difference of the spacecraft and the target at the meeting time.

First extrapolate the number of tracks to t_fThe track root difference [ Delta a ] is obtained at the moment₀,Δi₀,ΔΩ₀]Is [82000 m0.0008730.00585 ]]. Δ i here₀,ΔΩ₀All are radians.

And secondly, calculating the control quantity required by the intersection according to a first-order approximate model.

According to Δ Ω₀Solving the sum formula (11) to obtain [ x y z]＝[0.00103 6.61e-5 1.04e-4]. Then the optimal control quantity [ delta a, delta i, delta omega ] at the initial moment required for eliminating the ascension crossing point declination under the first-order approximation model']Are respectively [6814.99m 6.61 e-51.04 e-4 ]]。

And thirdly, calculating the control quantity required by intersection under the high-order approximation model by taking the solution of the second step as an initial value.

Using it as an initial value, substituting it into formula (15), wherein a_minThe radius of the earth plus 200km is set as the solution to [ x y z]＝[0.00104 -7.32e-5 0]. Substitution of formula (14) to eliminate [ Delta a ]₀,Δi₀,ΔΩ₀]Required optimum control amount Δ v₀And Δ v_fRespectively 4.06m/s and 44.27 m/s.

And step four, correcting the control quantity obtained in the step three according to the phase difference.

According to the number of orbits of the spacecraft and delta v₀Can be recalculated at t_fThe phase difference of the time target relative to the spacecraft is-3.31 rad, and the semi-long axis correction quantity required by corresponding delta t is 10381.17 m. Substitute it into x_fix＝(Δa+Δa_mod) (ii) when/a is 0.0026, the formula (20) is solved to obtain [ y z%]＝[-0.0014 -0.0067]. Re-substitution (14), i.e. correcting phase difference to obtain optimum control quantity delta v₀And Δ v_f22.75m/s and 41.79m/s, respectively.

And fifthly, correcting the control quantity obtained in the fourth step according to the eccentricity ratio vector difference.

According to the number of orbits of the spacecraft and delta v₀Calculated as t_fThe eccentricity difference of the target at the moment relative to the spacecraft is as follows:

u is represented by the formula (23)₀And u_fRespectively-0.26 and-1.36, and solving the formula (27) to obtain Delta e^rem0.0052, i.e. Δ v_r＝Δe^remV was 40.38 m/s. k is a radical of₀,k_fAre respectively 1 and-0.64, c₁Is 0.35. The four pulse sizes thus estimated are: 26.70m/s, 0.13m/s, 8.87m/s, 40.42 m/s.

And sixthly, establishing a four-pulse intersection optimization model and applying a differential evolution algorithm for optimization.

According to the estimated pulse and the local optimization model in the fifth step, solving X ═ X₁x₂x₃x₄x₅]Is [ -0.21,0.8,1.0, -0.28,0.56 [)]The optimal four pulse sizes when the analytic kinetic equation is used are respectively as follows: 23.33m/s, 9.02m/s, 15.72m/s, 26.64m/s, the sum being 65.68 m/s. The time is respectively as follows: 1760.71 seconds, 4461.18 seconds, 1206291.32 seconds, 1207800.56 seconds. Local optimization takes only 1 second or less.

And seventhly, replacing the analytic perturbation dynamic model with a full perturbation accurate dynamic model (10-order gravitational field, atmospheric resistance (coefficient 2.2, face value ratio 0.003) and daily and monthly gravitational perturbation), and solving the sixth step again by taking the optimal solution of the analytic dynamic model as an initial value.

The optimal four pulse sizes of the accurate dynamic model are respectively as follows: 20.89m/s, 8.07m/s, 9.98m/s, 38.70m/s, the sum being 69.57 m/s. The time is respectively as follows: 1760.71 seconds, 4461.18 seconds, 1205335.53 seconds, 1207465.19 seconds. Local optimization takes less than 2 minutes.

The method has the advantages of small error of the estimated speed increment and high solving speed. The method solves the problem of fast solving of the four-pulse optimal intersection.

The above description is only a preferred embodiment of the present invention, and is not intended to limit the scope of the present invention, and all modifications and equivalents of the present invention, which are made by the contents of the present specification and the accompanying drawings, or directly/indirectly applied to other related technical fields, are included in the scope of the present invention.

Claims (9)

1. A perturbation orbit four-pulse intersection rapid optimization method is characterized by comprising the following steps:

step 1, according to known conditions of a spacecraft and a target, obtaining control quantity and minimum speed increment through a control model taking elimination of orbital root difference at a meeting time as an optimization result;

changing the number of orbits when the spacecraft is controlled to run in the first circle of orbits according to the pulse control method of converting a part of control quantity and speed increment into two fixed intervals; converting the residual control quantity and the speed increment into two other pulse control spacecrafts with fixed intervals to finish intersection when the last circle of the orbit runs; estimating the size of four pulses and components in three directions of radial direction, tangential direction and normal direction;

step 2, establishing a four-pulse rendezvous optimization model, taking values of the time and the size of the first two pulses in the neighborhood according to the estimated value in the step 1, and recalculating the last two pulses according to the actual orbit recursion result;

and 3, replacing the analytic perturbation kinetic model with a full perturbation accurate kinetic model, and solving the step 2 again by taking the optimal solution of the analytic kinetic model as an initial value to obtain the optimal solution of the four-pulse intersection under the accurate kinetic model.

2. The perturbation orbit four-pulse intersection rapid optimization method according to claim 1, wherein the step 2 adopts a differential evolution algorithm to optimize a four-pulse intersection optimization model.

3. A perturbation orbit four-pulse intersection fast optimization method as claimed in claim 2, wherein said step 2 comprises:

step 21, directly applying the four-pulse estimation value obtained in the step 1 to a spacecraft, recurrently delivering the orbit to a known rendezvous moment, and defining a local search range coefficient of a speed increment and a local search range coefficient of a control quantity according to the deviation between the recurrently delivering orbit and a target orbit;

step 22, determining the moments of four pulses through an optimization model established by the time of pulse intersection, the phase of the spacecraft and the orbital angular velocity when the spacecraft intersects;

the magnitude of the last two pulses is calculated by the Lambert algorithm taking perturbation into account, step 23.

4. A perturbation orbit four-pulse intersection fast optimization method as claimed in claim 3, wherein said step 21 comprises:

in step 211, the four-pulse estimation values obtained in step 1 are respectively: Δ v₁,Δv₂,Δv₃,Δv₄The phase angles are respectively: u. of₁,u₂,u₃,u₄Wherein u is₁,u₂For the first turn after the start of the rendezvous task, u₃,u₄The last circle after the rendezvous task starts;

at step 212, four pulses are applied directly to the spacecraft and the orbit is recurred until time t_fDeviation from the target trajectory, defining four coefficients [ x ]₁x₂x₃x₄]Searching locallyWherein x is₁,x₂,x₃Local search range coefficients, x, representing the tangential, normal and radial components of the velocity increment, respectively₄A local search range coefficient representing the corresponding eccentricity control quantity in the step 1;

step 213, convert Δ v₁、Δv₂The changes are as follows:

wherein, delta a is the variation of the semi-major axis of the orbit of the spacecraft at the moment of the transfer initiation, delta i is the variation of the inclination angle of the orbit of the spacecraft at the moment of the transfer initiation, delta omega is the variation of the right ascension channel of the ascending intersection point of the spacecraft at the moment of the transfer initiation, and delta v_rTo require an additional amount of radial control, a₁For a predetermined upper limit of the local change of the semi-major axis, k₀I is the orbit inclination angle for the coefficient of the variation of the control quantity to the orbit eccentricity of the spacecraft,

the average orbit velocity of the spacecraft.

5. A perturbation orbit four-pulse intersection fast optimization method as claimed in claim 4, wherein said step 22 comprises:

step 221, the number of orbits of the spacecraft is pushed to the starting moment t of the handover control_f- Δ t, set to [ a_s,e_s,i_s,Ω_s,ω_s,M_s]If the phase at that time is u ═ ω_s+M_sDue to the pulse Δ v₁Is u₁Then the pulse Deltav is known₁Time deltat relative to the start of a session₁Comprises the following steps:

Δt₁＝mod(u₁-u,2π)/n (32)

wherein mod represents the sum of₁-u is normalized to between 0 and 2 π,

for the orbital angular velocity at the beginning of a spacecraft rendezvous, namely, the spacecraft rendezvous needs to wait for delta t₁After that, the first pulse is applied, and Δ v is known from the formula (31)₁The three-dimensional components in the spacecraft orbit coordinate system are:

the number of orbits of the spacecraft is pushed to t_f-Δt+Δt₁At the moment, the number of tracks is converted into a position velocity vector, and delta v is superposed on the velocity₁Obtaining a spacecraft position and speed vector after the first pulse, and then converting to obtain the orbit number;

step 222, due to Δ v₁And Δ v₂Phase difference of u₂-u₁Pi, so that the second pulse is separated from the first by a time at₂Comprises the following steps:

Δt₂＝π/n₁(34)

wherein

To add Δ v₁New orbital angular velocity of the spacecraft after, i.e. passing at₂Applying a second pulse of formula (31), Δ v₂The three-dimensional components in the spacecraft orbit coordinate system are:

then the orbit number of the spacecraft after the second pulse is pushed to t_fTime is set as [ a ]_f,e_f,i_f,Ω_f,ω_f,M_f]The number of tracks [ a ] at the same time as the target_t(t_f),e_t(t_f),i_t(t_f),Ω_t(t_f),ω_t(t_f),M_t(t_f)]The difference shows that the right ascension angle Δ Ω at this time_fAnd the difference of inclination angle Δ i_f；

Step 223, the phases u of the third and fourth pulses₃,u₄The change is:

u₃＝atan((Δi_f)/(ΔΩ_fsini/2))

u₄＝atan((Δi_f)/(ΔΩ_fsini/2))+π (36)

the time at which the fourth pulse is applied and t_fThe moments not necessarily coinciding, with a time difference Δ t₄Comprises the following steps:

Δt₄＝-(mod(u₄-(ω_f+M_f),2π)-2π)/n_f(37)

wherein

Track angular velocity for the rendezvous target, (mod (u)₄-(ω_f+M_f) 2 π) -2 π) represents u₄-(ω_f+M_f) Normalized to between-2 pi and 0, i.e. the fourth pulse has a time t_f-Δt₄The time advance delta t of the third pulse relative to the fourth pulse₃Comprises the following steps:

Δt₃＝π/n₃(38)

wherein n is₃For the angular velocity of the spacecraft orbit after the third pulse is applied, since n₃Is difficult to be directly calculated, and introduces a correction coefficient x to be solved for compensating eccentricity₅Let us order

Δt₃＝x₅π/n_f(39)

Thereby determining the time instants of the third and fourth pulses.

6. A perturbation orbit four-pulse intersection fast optimization method as claimed in claim 5, wherein said step 23 comprises:

step 231, for a given start-stop time t_in、t_outAnd corresponding start-stop position vector r_in、r_outCalculating a start-stop velocity vector v for a two-body dynamics model_in、v_out：

[v_inv_out]＝Lambert(r_in,r_out,t_in,t_out) (40)

Applying the Lambert algorithm, i.e. setting the starting time t_inIs set to t_f-Δt₄-Δt₃End time t_outIs set to t_f-Δt₄At a starting position r_inFor the number of orbits of the spacecraft after the second pulse to be pushed to t_f-Δt₄-Δt₃Position vector obtained at time, end position r_outStepping to t for target trajectory_f-Δt₄A position vector obtained at a time;

step 232, firstly, obtaining the v corresponding to the two-body dynamic model by using an analytic calculation method_inAnd v_outAnd obtaining v corresponding to the J2 analytic perturbation dynamics model by differential correction_inAnd v_out；

Step 233, set v₃ ^-For the number of orbits of the spacecraft after the second pulse to be pushed to t_f-Δt₄-Δt₃The velocity vector taken at the moment, the third pulse needs to change the spacecraft velocity to v_inThe meeting location constraint can be satisfied; v. the₄ ⁺Stepping to t for target trajectory_f-Δt₄The velocity vector obtained at the moment, the fourth pulse needs to drive the spacecraft velocity from v_outChange to v₄ ⁺The meeting speed constraint can be met; thus Δ v₃And Δ v₄The sizes are respectively as follows:

the total velocity increment for the four pulses is then:

Δv_sum＝Δv₁+Δv₂+Δv₃+Δv₄(42)

to sum up, the local optimization model includes 5 variables: x is the number of₁,x₂,x₃,x₄,x₅The optimization index is Δ v calculated by the following equations (31) to (42)_sumMinimum; x is the number of₁,x₂,x₃,x₄,x₅The value ranges are as follows:

solving by using a differential evolution algorithm to obtain optimal x₁,x₂,x₃,x₄,x₅And accurate values of the four pulse moments and the four pulse magnitudes under the analytic dynamics model can be obtained.

7. A perturbation orbit four-pulse intersection fast optimization method as claimed in claim 6, wherein said step 3 comprises:

step 31, replacing the J2 analytic perturbation dynamics model with an accurate dynamics equation of the following formula, and applying the equation to all the processes related to the orbit recursion calculation in step 2:

wherein r is the three-dimensional position vector of the spacecraft, v is the three-dimensional velocity vector, a_PerturbationThe acceleration caused by the gravity forces such as the non-spherical gravity of the earth, the atmospheric resistance, the sun-moon gravity, the sunlight pressure and the like is considered;

step 32, using x solved in step 2₂,x₃,x₄Re-optimizing x as an optimal solution for accurate kinetic models₁,x₅；

And step 33, recalculating the optimization models corresponding to the expressions (31) to (42) by using a differential evolution algorithm, so as to obtain the four pulse moments and the four pulse magnitudes with the optimal precise dynamics.

8. A perturbation orbit four-pulse intersection fast optimization method as claimed in any one of claims 1 to 7, wherein in step 1, the known conditions include: initial time of the spacecraft and the target, the number of orbits of the spacecraft and the target at the initial time and rendezvous time;

the predetermined period of the interval between the two pulses is half a period.

9. A perturbation orbit four-pulse intersection fast optimization method as claimed in claim 8, wherein said step 1 comprises:

step 11, respectively calculating the orbit root difference of the spacecraft and the target which need to be eliminated by the pulse at the intersection moment according to the known conditions and the J2 analytic dynamic model;

step 12, designing an optimal control model for eliminating the track root number difference, and obtaining control quantity and minimum speed increment required by intersection after pulse elimination;

step 13, taking the solution obtained in the step 12 as an initial value, and calculating the minimum speed increment required by intersection under a high-order approximation model;

step 14, correcting the minimum speed increment obtained in the step 13 according to the phase difference to obtain the minimum speed increment required by intersection after the phase difference is corrected;

and step 15, correcting the minimum speed increment obtained in the step 14 according to the eccentricity vector difference, and obtaining the estimated size of each pulse and components in three directions, namely radial direction, tangential direction and normal direction.

Images