A semi-analytical optimization method for out-of-plane rendezvous of multi-pulse elliptical orbits

Abstract

The present invention provides a semi-analytical multi-pulse elliptical orbit out-of-plane rendezvous optimization method. The method first uses a four-pulse rendezvous velocity increment rapid estimation method to obtain estimated four-pulse initial parameters. A precise dynamics model is then used to predict the difference between the orbit after the four-pulse maneuver and the target orbit to be rendezvoused. An analytical expression is then given to calculate a correction algorithm for the four-pulse size and timing. By iteratively repeating the correction and prediction steps, the optimal four-pulse vector and application timing are ultimately obtained, resulting in a zero difference between the orbit after the pulse maneuver and the target orbit to be rendezvoused. The method can be applied to low-orbit space debris removal and on-orbit servicing mission orbit optimization design, resolving the issues of multiple local optimal solutions and low computational efficiency during long-term out-of-plane rendezvous optimization of perturbed orbits.

Description

Technical Field

[0001] The invention belongs to the technical field of aerospace navigation control, and in particular relates to a semi-analytical multi-pulse elliptical orbit out-of-plane rendezvous optimization method. Background Art

[0002] Spacecraft in near-Earth space are significantly affected by perturbations, which must be accounted for during orbital rendezvous optimization. Furthermore, rendezvous on out-of-plane orbits is limited by the spacecraft's inherent thrust, requiring extended maneuvers to complete. These factors make traditional orbital rendezvous optimization methods computationally intensive and prone to falling into local optimal solutions. Summary of the Invention

[0003] The purpose of the present invention is to provide a semi-analytical multi-pulse elliptical orbit out-of-plane intersection optimization method, which solves the problem of multiple local optimal solutions and low computational efficiency in the optimization process in the prior art. The technical solution adopted by the present invention is a semi-analytical multi-pulse elliptical orbit out-of-plane intersection optimization method, characterized by being implemented in the following steps:

[0004] Step 1: Calculate the parameters that can represent the approximately optimal rendezvous four-pulse according to the input initial orbit element number and the target orbit element number to be rendezvoused, as well as the four-pulse estimated values ​​and corresponding four phase values ​​directly calculated from these parameters;

[0005] Step 2: Calculate the time to apply the four pulses based on the four phase values, and then predict the spacecraft at the specified rendezvous time t f The number of orbital elements;

[0006] Step 3: Calculate the difference between the orbital elements predicted in step 2 and the orbital elements of the given target orbit, and record it in the form of first-kind singularity-free elements;

[0007] Step 4: Calculate the correction amount for the semi-major axis control amount based on the first and sixth components of the orbital root difference;

[0008] Step 5: Calculate the corrections to the inclination and right ascension of the ascending node based on the fourth and fifth components of the orbital element difference.

[0009] Step 6: Calculate the correction amount for the eccentricity control amount based on the second and third components of the orbital element difference;

[0010] Step 7: Substitute the corrected control quantity back into step 1, recalculate the size and phase of the four pulses, and then substitute them back into steps 2 and 3. Predict the updated orbital element difference after parameter correction and determine whether it is less than the given threshold. If so, output the size and timing of the four pulses as the optimal control quantity. Otherwise, continue to use steps 4, 5, and 6 to correct the corresponding control quantity until the predicted orbital element difference is less than the given threshold.

[0011] The present invention is also characterized in that:

[0012] Step 1 is implemented as follows:

[0013] According to the input of the initial orbital elements of the spacecraft [a0, e0, i0, Ω0, ω0, M0], a represents the semi-major axis, e represents the eccentricity, i represents the inclination, Ω represents the right ascension of the ascending node, ω represents the argument of perigee, M represents the mean anomaly, subscript 0 represents the orbital elements of the spacecraft at the initial moment, and the orbital elements of the target orbit that needs to be rendezvoused, f is the subscript representing the target orbit, the orbital elements corresponding to the spacecraft and the target orbit at the rendezvous moment are predicted, and the difference is taken to obtain the orbital element difference described in the form of the first kind of non-singular point elements, which is recorded as [△a0,△e x0 ,△e y0 ,△i0,△Ω0,△u0], defined as:

[0014] △a0=a f -a0

[0015] △e x0 =e f cosω f -e0cosω0

[0016] △e y0 =e f sinω f -e0sinω0

[0017] Δi0=i f -i0

[0018] ΔΩ0=Ω f -Ω0

[0019] Δu0=ω f +M f -ω0-M0

[0020] According to [△a0,△e x0 ,Δe y0 ,Δi0,ΔΩ0,Δu0], the optimal four-pulse intersection size calculated is described as formula (1), where the first two pulses act on the first circle, the last two pulses act on the last circle, and the sizes of the four pulse velocity increments Δv1, Δv2, Δv3, and Δv4 are respectively:

[0021]

[0022] Among them, △v1, △v2, △v3, △v4 represent four pulse velocity increments, a=a0 is the initial semi-major axis, i=i0 is the initial inclination angle, represents the initial orbital average velocity, μ is the Earth's gravitational constant, Δa represents the change in the semi-major axis caused by the first two pulses, △i represents the change in the inclination caused by the first two pulses, △Ω / 2 represents the change in the right ascension of the ascending node caused by the first two pulses, and k r1 k represents the change in eccentricity of the radial component in the first two pulses. r2 represents the change of the eccentricity due to the radial component in the last two pulses, k0 represents the change of the eccentricity due to the tangential component in the first two pulses, and k f It represents the change of eccentricity due to the tangential component in the last two pulses;

[0023] The phase angles u1, u2, u3, and u4 of the four pulse maneuvering moments corresponding to △v1, △v2, △v3, and △v4 are:

[0024]

[0025] Among them, u1 and u2 correspond to a certain moment in the first circle after the start of orbit transfer, and u3 and u4 correspond to a certain moment in the last circle before rendezvous; according to the definition of orbital elements, u1, u2, u3, u4 are one-to-one corresponding to t1, t2, t3, t4, so the optimal four-pulse results are expressed by Δv1, Δv2, Δv3, Δv4 and u1, u2, u3, u4.

[0026] Step 2 is implemented as follows:

[0027] Step 2.1, calculate the time t1 when the first pulse is applied:

[0028] Assuming that the initial phase of the spacecraft at time t0 is u0=ω0+M0(3), the phase difference between u0 and u1 includes four cases:

[0029] Formula (4) ensures that Δu1 is always between 0 and π by treating different values ​​of Δu1 differently. According to the near-circular orbit assumption, the time of the first pulse relative to the initial moment is

[0030]

[0031] in is the average angular velocity of the orbit corresponding to the number of orbital elements at the initial moment;

[0032] Step 2.2, calculate the time t2 when the second pulse is applied:

[0033] The orbital elements of the spacecraft at the initial moment are predicted to the moment t1=t0+Δt1. At this moment, the velocity vector of the spacecraft is added with the velocity increment vector Δv1. The three-dimensional components of Δv1 in the radial, tangential and normal directions are:

[0034]

[0035] in: e t 、e r 、e n Represent the tangential unit vector, radial unit vector and normal unit vector respectively, r=[x,y,z], v=[v x ,v y ,v z ] are the position and velocity vectors of the spacecraft at time t1, r = |r|, v = |v| are the moduli of the position and velocity vectors;

[0036] At time t1, adding Δv1 to the spacecraft velocity vector gives the position velocity of the spacecraft after the first pulse is applied, which can be converted into six orbital numbers. According to formula (2), the second pulse phase u2 is: u2=u1+π(8),

[0037] Therefore, the time difference Δt2 of the second pulse relative to the first pulse time t1 is:

[0038]

[0039] Where n1 is the orbital angular velocity after the first pulse. The time when the second pulse is applied is then obtained as t2 = Δt2 + t1. The spacecraft orbit is predicted to time t2, and the second pulse Δv2 is added to the velocity vector:

[0040]

[0041] The definition of variables in formula (10) is the same as formula (6). The definitions of case 1, case 2, case 3, and case 4 are consistent with formula (4). At this time, the position velocity of the spacecraft after the first and second pulses are applied is obtained;

[0042] Step 2.3, calculate the third pulse time t3:

[0043] Since the third and fourth pulses are in the last circle, that is, the intersection time t f Nearby, first predict the orbit after applying the second pulse to t f At this moment, we get the orbital element number [a(t f ),e(t f ),i(t f ),Ω(t f ),ω(t f ),M(t f )], then t f The phase of the spacecraft at time u f for:

[0044] u f =ω f+M f (11)

[0045] According to formula (2), the phases of the third and fourth pulses are u3, u4 = u3 + π, and u4 and u f The difference between u3 and u f The phase difference Δu3 must be between π and 2π. Calculate Δu3 according to the following formula:

[0046]

[0047] Where mod(x, y) represents x modulo y; then the third pulse time t3 is:

[0048]

[0049] Predict the spacecraft orbit to time t3 and add the third pulse △v3 to the velocity vector:

[0050]

[0051] The definitions of Case 1, Case 2, Case 3, and Case 4 are consistent with those of Formula (12); e t ,e n ,e r The calculation of is the same as formula (6);

[0052] Step 2.4: Calculate the time t4 and t f The number of spacecraft orbital elements after four pulses have been applied at time t;

[0053] The time difference Δt4 of the fourth pulse relative to the third pulse is:

[0054] Where n3 is the orbital angular velocity after the third pulse. Predict the spacecraft orbit to time t4 and add the fourth pulse △v4 to the velocity vector:

[0055]

[0056] The definitions of Case 1, Case 2, Case 3, and Case 4 are consistent with those of Formula (12); e t ,e n ,e r The calculation is the same as formula (6); obtain the new position velocity and predict t f The number of spacecraft orbital elements at the time is denoted as [a p ,e p ,i p ,Ω p ,ω p ,M p Step 3 is implemented as follows:

[0057] Define the orbital element difference [Δa p ,Δe xp ,Δe yp ,Δi p ,ΔΩ p ,Δu p ] is the difference between the orbit of the spacecraft to which the four pulses are applied and the target orbit to be rendezvoused, expressed in the form of first-kind singularity-free elements:

[0058]

[0059] The smaller the difference is, the smaller the parameters [Δa, Δi, △Ω, k0, k f ,k r1 ,k r2 ] and the four-pulse estimated values ​​△v1, △v2, △v3, △v4 calculated directly from these seven parameters are closer to the exact value, the more accurate the value is, the better the result is. p ,Δe xp ,Δe yp ,Δi p ,ΔΩ p ,Δu p ] for [Δa,Δi,ΔΩ,k0,k f ,k r1 ,k r2 ] to make corrections.

[0060] Step 4 is implemented as follows:

[0061] According to Δa p ,Δu p ,Δi p ,ΔΩ p , calculate the correction amount for Δa as follows:

[0062] Assume that the correction amount of Δa is Δa1, then Δa1 and Δu p The following relationship must be satisfied:

[0063]

[0064] Where a0 and n0 have been defined in equation (5), Δt = t f -t0 represents the time difference from the start to the intersection;

[0065] At the same time, when recalculating △v1, △v2, △v3, and △v4, the term (△a0-△a) in equations (16) and (14) needs to be corrected to:

[0066] (△a0-△a p )-(△a+△a1) (19)

[0067] The results of this step are the correction amounts for Δa and (Δa-Δa0).

[0068] Step 5 is implemented as follows:

[0069] Keeping u1, u2, u3, and u4 unchanged, let △i in equations (6) and (10) be corrected to △i+△i1, and △i0-△i in equations (16) and (14) be corrected to △i0-△i+△i2, then △Ω in equations (6) and (10) will be proportionally changed to The △Ω in formula (16) and formula (14) will be proportionally changed to △i+△i1,△i0-△i+△i2, and △Ω should satisfy the following relationship:

[0070]

[0071] in It represents the right ascension drift rate of the ascending node after applying one or two pulses to the initial orbit. It represents the change in the right ascension drift rate of the ascending node caused by △a1 and △i1; according to the law of the J2 perturbation term in orbital mechanics, it is:

[0072]

[0073] Combining formula (2), formula (20) is transformed into:

[0074]

[0075] Formula (22) is a system of linear equations in two variables, which can be solved directly:

[0076]

[0077] Then, the correction values ​​of △i, △i0-△i and △Ω used in calculating △v1, △v2, △v3 and △v4 in equations (6), (10), (16) and (14) are obtained.

[0078] Step 6 is implemented as follows:

[0079] Let k0,k f ,k r1 ,k r2 The correction amount is Δk0, Δk f ,Δk r1 ,Δk r2 The eccentricity change caused by the correction should offset the Δe calculated by formula (17) xp ,Δe yp ,Right now:

[0080]

[0081] Where Δu1 represents the orbital perigee argument drift from the time the first and second pulses are applied until the last revolution:

[0082]

[0083] Among them, J2 is the second-order term of the non-spherical perturbation of the earth, R e is the Earth's equatorial radius,

[0084] Formula (24) can be written as:

[0085]

[0086] Then Δk0, Δk f ,Δk r1 ,Δk r2 The solution is considered as Δk0, Δk f ,Δk r1 ,Δk r2 To optimize the variables, use Equation (27) as the equality constraint to minimize the sum of △v1, △v2, △v3, and △v4:

[0087] min f=|△v1|+|△v2|+|△v3|+|△v4|

[0088]

[0089]

[0090] Obtain Δk0, Δk that minimizes Equation (28) f ,Δk r1 ,Δk r2 .

[0091] Step 7 is as follows:

[0092] The corrected Δa, Δi, ΔΩ, k0, k f ,k r1 ,k r2 Substitute back into step 1, recalculate Δv1, Δv2, Δv3, △v4 and u1, u2, u3, u4, and then substitute back into steps 2 and 3 to predict the updated orbital root difference [△a p ,△e xp ,△e yp ,△u p ,△i p ,△Ω p ], judge whether it is less than a given threshold, if so, output △v1,△v2,△v3,△v4 and t1,t2,t3,t4; otherwise, continue with steps 4, 5, and 6 to compare △a,△i,△Ω,k0,k f ,kr1 ,k r2 Correction is made until the forecast [△a p ,△e xp ,△e yp ,△u p ,△i p ,△Ω p ] is less than a given threshold and stops.

[0093] The beneficial effect of the present invention is a semi-analytical multi-pulse elliptical orbit out-of-plane rendezvous optimization method. It first uses a four-pulse rendezvous velocity increment rapid estimation method to obtain initial pulse values ​​with a certain error. Then, a precise dynamic model is used to predict the difference between the orbit after the pulse maneuver and the rendezvous target orbit. Finally, an analytical expression is given to calculate the pulse correction value, and finally the optimal four pulses required for precise rendezvous are obtained. This method is used to solve the problem of multiple local optimal solutions and low computational efficiency in pulse optimization for long-term out-of-plane rendezvous on perturbed orbits. BRIEF DESCRIPTION OF THE DRAWINGS

[0094] Figure 1 This is a flow chart of a semi-analytical multi-pulse elliptical orbit out-of-plane rendezvous optimization method of the present invention. DETAILED DESCRIPTION

[0095] The present invention will be described in detail below with reference to the accompanying drawings and specific embodiments.

[0096] The problem of the present invention is set as follows:

[0097] Given the six orbital elements of the spacecraft at the initial moment of the mission [a0, e0, i0, Ω0, ω0, M0] and the six orbital elements of the target [a f ,e f ,i f ,Ω f ,ω f ,M f ], given the intersection time t f (The precise rendezvous moment can be fine-tuned around this) and the flight duration △t, solve the four optimal pulse maneuver vectors △v1, △v2, △v3, △v4 and the corresponding four times t1, t2, t3, t4 required to make the spacecraft rendezvous with the target orbit.

[0098] A semi-analytical multi-pulse elliptical orbit out-of-plane rendezvous optimization method, the flow chart is as follows Figure 1 As shown, please follow the steps below:

[0099] Step 1: Calculate the parameters that can represent the approximately optimal rendezvous four-pulse according to the input initial orbit element number and the target orbit element number to be rendezvoused, as well as the four-pulse estimated values ​​and corresponding four phase values ​​directly calculated from these parameters;

[0100] Step 1 is implemented as follows:

[0101] Using the conventional orbit prediction method, according to the input of the initial orbit elements of the spacecraft [a0, e0, i0, Ω0, ω0, M0], a represents the semi-major axis, e represents the eccentricity, i represents the inclination, Ω represents the right ascension of the ascending node, ω represents the argument of perigee, M represents the mean anomaly, subscript 0 represents the orbit elements of the spacecraft at the initial moment, and the target orbit elements that need to be rendezvoused, f is the subscript indicating that it is the target orbit, in order to distinguish it from subscript 0, the orbit elements corresponding to the spacecraft and target orbits at the rendezvous moment are predicted, and the difference is taken to obtain the orbit element difference described in the form of the first kind of non-singular point elements, which is recorded as [△a0,△e x0 ,△e y0 ,△i0,△Ω0,△u0], defined as:

[0102] △a0=a f -a0

[0103] △e x0 =e f cosω f -e0cosω0

[0104] △e y0 =e f sinω f -e0sinω0

[0105] △i0=i f -i0

[0106] △Ω0=Ω f -Ω0

[0107] △u0=ω f +M f -ω0-M0

[0108] The orbit prediction method and the conversion between orbit elements and position velocity are common knowledge in this field and will not be described here or in the following sections. Taking the method in reference [1] as an example, according to [△a0,△e x0 ,△e y0 ,△i0,△Ω0,△u0], the optimal four-pulse intersection size is described as formula (1), where the first two pulses act on the first circle, the last two pulses act on the last circle, and the sizes of the four pulse velocity increments △v1, △v2, △v3, and △v4 are respectively:

[0109]

[0110] Among them, △v1, △v2, △v3, △v4 represent four pulse velocity increments, a=a0 is the initial semi-major axis, i=i0 is the initial inclination angle, represents the initial orbital average velocity, μ is the Earth's gravitational constant, △a represents the change in the semi-major axis caused by the first two pulses, △i represents the change in the inclination caused by the first two pulses, △Ω / 2 represents the change in the right ascension of the ascending node caused by the first two pulses, and k r1 k represents the change in eccentricity of the radial component in the first two pulses. r2 represents the change of the eccentricity due to the radial component in the last two pulses, k0 represents the change of the eccentricity due to the tangential component in the first two pulses, and k f It represents the change of eccentricity due to the tangential component in the last two pulses;

[0111] The phase angles u1, u2, u3, and u4 of the four pulse maneuvering moments corresponding to △v1, △v2, △v3, and △v4 are:

[0112]

[0113] Among them, u1 and u2 correspond to a certain moment in the first circle after the start of the orbit transfer, and u3 and u4 correspond to a certain moment in the last circle before the rendezvous; according to the definition of the orbital elements, u1, u2, u3, u4 and t1, t2, t3, t4 are one-to-one corresponding, so the optimal four-pulse result is represented by △v1, △v2, △v3, △v4 and u1, u2, u3, u4. [△a, △i, △Ω, k0, k f ,k r1 ,k r2 ] These seven parameters are calculated by the approximate method in the literature [1], so △v1, △v2, △v3, △v4 and u1, u2, u3, u4 calculated by formula (1) and formula (2) are all approximate values.

[0114] Step 2: Calculate the time to apply the four pulses based on the four phase values, and then predict the spacecraft at the specified rendezvous time t f The number of orbital elements;

[0115] Step 2 is implemented as follows:

[0116] Step 2.1, calculate the time t1 when the first pulse is applied:

[0117] Assuming that the initial phase of the spacecraft at time t0 is u0=ω0+M0(3), the phase difference between u0 and u1 includes four cases:

[0118] Formula (4) ensures that △u1 is always between 0 and π by treating different values ​​of △u1 differently. According to the near-circular orbit assumption, the time of the first pulse relative to the initial moment is

[0119]

[0120] in is the average angular velocity of the orbit corresponding to the number of orbital elements at the initial moment;

[0121] Step 2.2, calculate the time t2 when the second pulse is applied:

[0122] The orbital elements of the spacecraft at the initial moment are predicted to the moment t1=t0+△t1. At this moment, the velocity vector of the spacecraft is added with the velocity increment vector △v1. The three-dimensional components of △v1 in the radial, tangential and normal directions are:

[0123]

[0124] in: e t 、e r 、e n Represent the tangential unit vector, radial unit vector and normal unit vector respectively, r=[x,y,z], v=[v x ,v y ,v z ] are the position and velocity vectors of the spacecraft at time t1, r = |r|, v = |v| are the moduli (scalars) of the position and velocity vectors;

[0125] At time t1, adding △v1 to the spacecraft velocity vector gives the position velocity of the spacecraft after the first pulse is applied, which can be converted into six orbital numbers. According to formula (2), the second pulse phase u2 is: u2=u1+π(8),

[0126] Therefore, the time difference Δt2 of the second pulse relative to the first pulse time t1 is:

[0127]

[0128] Where n1 is the orbital angular velocity after the first pulse. The time when the second pulse is applied is then obtained as t2 = Δt2 + t1. The spacecraft orbit is predicted to time t2, and the second pulse Δv2 is added to the velocity vector:

[0129]

[0130] The definition of variables in formula (10) is the same as formula (6). The definitions of cases 1, 2, 3, and 4 are consistent with formula (4). At this time, the position velocity of the spacecraft after the first and second pulses are applied (which can be converted into six orbital numbers) is obtained;

[0131] Step 2.3, calculate the third pulse time t3:

[0132] Since the third and fourth pulses are in the last circle, that is, the intersection time t f Nearby, first predict the orbit after applying the second pulse to t f At this moment, we get the orbital element number [a(t f ),e(t f ),i(t f ),Ω(t f ),ω(t f ),M(t f )], then t f The phase of the spacecraft at time u f (Definition is the same as formula (3)) is:

[0133] u f =ω f +M f (11)

[0134] According to formula (2), the phases of the third and fourth pulses are u3, u4 = u3 + π, and u4 and u f The difference between u3 and u f The phase difference △u3 must be between π and 2π. Calculate △u3 according to the following formula:

[0135]

[0136] Where mod(x, y) represents x modulo y; then the third pulse time t3 is:

[0137]

[0138] Predict the spacecraft orbit to time t3 and add the third pulse △v3 to the velocity vector:

[0139]

[0140] The definitions of Case 1, Case 2, Case 3, and Case 4 are consistent with those of Formula (12); e t ,e n ,e r The calculation of is the same as formula (6);

[0141] Step 2.4: Calculate the time t4 and t f The number of spacecraft orbital elements after four pulses have been applied at time t;

[0142] The time difference Δt4 of the fourth pulse relative to the third pulse is:

[0143] Where n3 is the orbital angular velocity after the third pulse. Predict the spacecraft orbit to time t4 and add the fourth pulse △v4 to the velocity vector:

[0144]

[0145] The definitions of Case 1, Case 2, Case 3, and Case 4 are consistent with those of Formula (12); e t ,e n ,e r The calculation is the same as formula (6); obtain the new position velocity and predict t f The number of spacecraft orbital elements at the time is denoted as [a p ,e p ,i p ,Ω p ,ω p ,M p ](p is a meaningless subscript, only used to distinguish the previous orbital elements, indicating that it is the orbit obtained after adding 4 pulse control variables). Step 3. Calculate the difference between the orbital elements predicted in step 2 and the orbital elements of the given target orbit, and record it in the form of first-kind non-singular point elements;

[0146] Step 3 is implemented as follows:

[0147] Define the orbital root difference [△a p ,△e xp ,△e yp ,△i p ,△Ω p ,△u p ] is the difference between the orbit of the spacecraft to which the four pulses are applied and the target orbit to be rendezvoused, expressed in the form of first-kind singularity-free elements:

[0148]

[0149] The smaller the difference is, the smaller the parameters [△a,△i,△Ω,k0,k f ,k r1 ,k r2 ] and the four-pulse estimated values ​​△v1, △v2, △v3, △v4 calculated directly from these seven parameters are closer to the exact value, but generally not 0, and need to be calculated based on [△a p ,△e xp ,△e yp ,△i p ,△Ω p ,△u p ] for [△a,△i,△Ω,k0,k f ,k r1 ,k r2 ] to make corrections.

[0150] Step 4: Calculate the correction amount for the semi-major axis control amount based on the first and sixth components of the orbital root difference;

[0151] Step 4 is implemented as follows:

[0152] According to △a p ,△u p ,△i p ,△Ω p , calculate the correction amount for △a, as follows:

[0153] Assume that the correction amount of △a is △a1, then △a1 and △u p The following relationship must be satisfied:

[0154]

[0155] Where a0 and n0 have been defined in equation (5), △t=t f -t0 represents the time difference from the start to the intersection;

[0156] At the same time, when recalculating △v1, △v2, △v3, and △v4, the term (△a0-△a) in equations (16) and (14) needs to be corrected to:

[0157] (△a0-△a p )-(△a+△a1)(19)

[0158] The results of this step are the correction amounts for Δa and (Δa-Δa0).

[0159] Step 5: Calculate the corrections to the inclination and right ascension of the ascending node based on the fourth and fifth components of the orbital element difference.

[0160] Step 5 is implemented as follows:

[0161] Keeping u1, u2, u3, and u4 unchanged, let △i in equations (6) and (10) be corrected to △i+△i1, and △i0-△i in equations (16) and (14) be corrected to △i0-△i+△i2, then △Ω in equations (6) and (10) will be proportionally changed to The △Ω in formula (16) and formula (14) will be proportionally changed to △i+△i1,△i0-△i+△i2, and △Ω should satisfy the following relationship:

[0162]

[0163] in It represents the right ascension drift rate of the ascending node after applying one or two pulses to the initial orbit. It represents the change in the right ascension drift rate of the ascending node caused by △a1 and △i1; according to the law of the J2 perturbation term in orbital mechanics, it is:

[0164]

[0165] Combining formula (2), formula (20) is transformed into:

[0166]

[0167] Formula (22) is a system of linear equations in two variables, which can be solved directly:

[0168]

[0169] Then, the correction values ​​of △i, △i0-△i and △Ω used in calculating △v1, △v2, △v3 and △v4 in equations (6), (10), (16) and (14) are obtained.

[0170] Step 6: Calculate the correction amount for the eccentricity control amount based on the second and third components of the orbital element difference;

[0171] Step 6 is implemented as follows:

[0172] Let k0,k f ,k r1 ,k r2 The correction amount is △k0,△k f ,△k r1 ,△k r2 The eccentricity change caused by the correction should offset the △e calculated by formula (17) xp ,△e yp ,Right now:

[0173]

[0174] Where △u1 represents the orbital perigee argument drift from the time the first and second pulses are applied until the last revolution:

[0175]

[0176] Among them, J2 is the second-order term of the non-spherical perturbation of the earth, R e is the equatorial radius of the Earth. Since the formula (26) was used in the calculation in the literature [1]:

[0177]

[0178] Formula (24) can be written as:

[0179]

[0180] Then △k0,△kf ,△k r1 ,△k r2 The solution is considered as △k0,△k f ,△k r1 ,△k r2 To optimize the variables, use Equation (27) as the equality constraint to minimize the sum of △v1, △v2, △v3, and △v4:

[0181] min f=|△v1|+|△v2|+|△v3|+|△v4|

[0182]

[0183]

[0184] The present invention uses the differential evolution algorithm to solve the four-dimensional optimization problem, and can quickly obtain △k0, △k that minimizes formula (28). f ,△k r1 ,△k r2 The differential evolution algorithm is a known public method, and its detailed process will not be repeated here.

[0185] Step 7: Substitute the corrected control quantity back into step 1, recalculate the size and phase of the four pulses, and then substitute them back into steps 2 and 3. Predict the updated orbital element difference after parameter correction and determine whether it is less than the given threshold. If so, output the size and timing of the four pulses as the optimal control quantity. Otherwise, continue to use steps 4, 5, and 6 to correct the corresponding control quantity until the predicted orbital element difference is less than the given threshold.

[0186] Step 7 is as follows:

[0187] The corrected △a,△i,△Ω,k0,k f ,k r1 ,k r2 Substitute back into step 1, recalculate △v1,△v2,△v3,△v4 and u1,u2,u3,u4, then substitute back into steps 2 and 3, and predict the updated orbital root difference [△a p ,△e xp ,△e yp ,△u p ,△i p ,△Ω p ], judge whether it is less than a given threshold, if so, output △v1,△v2,△v3,△v4 and t1,t2,t3,t4; otherwise, continue with steps 4, 5, and 6 to compare △a,△i,△Ω,k0,k f ,k r1 ,k r2 Correction is made until the predicted [△ap ,△e xp ,△e yp ,△u p ,△i p ,△Ω p ] is less than a given threshold and stops.

[0188] Example:

[0189] The known initial spacecraft orbit and the number of orbital elements required for rendezvous are shown in Table 1. Rendezvous is required to be achieved after 25 days.

[0190] Table 1 Orbital elements of the spacecraft and the rendezvous target at the initial moment

[0191] a(m) e i (radians) Ω (radians) ω (radians) M (radians) spacecraft 7131917.466 0.01462 1.7304 3.0759 -0.2124 3.0366 Target 7058041.710 0.00872 1.7206 3.1247 -4.5384 6.0577

[0192] The calculation process according to the present invention is as follows:

[0193] In the first step, the approximate four pulse sizes △v1, △v2, △v3, and △v4 are calculated to be 74.55m / s, 0m / s, 64.58m / s, and 44.64m / s respectively. The four pulse phases u1, u2, u3, and u4 are 3.8274, 0.6858, 5.8232, and 2.6816 respectively. f ,k r1 ,k r2 They are 1, 0.5913, -0.0001851, and 0.0003683 respectively.

[0194] The second step is to use the approximate four-pulse result to calculate the orbital element difference between the spacecraft and the target at the rendezvous moment [△a p ,△e xp ,△e yp ,△i p ,△Ω p ,△u p ]=[-886.746m,9.78e-5,3.56e-4,1.19e-4,1.05e-3,-2.10].

[0195] The third step is to calculate △a1=4409.525m and △a2=-5296.271m according to formula (18).

[0196] The fourth step is to calculate △i1=4.99e-4,△i2=-3.78e-4,△Ω1=4.01e-4,△Ω2=1.87e-4 according to formula (23).

[0197] The fifth step is to solve the optimization model of formula (28) and obtain △k0, △k f ,△k r1 ,△k r2They are 0, 0.009334, 0.0001887, and 0 respectively. The minimum value of f is 0.052 m / s, which means that the total speed increment for the optimal rendezvous is 183.832 m / s.

[0198] Applying the corrected four-pulse results to the prediction process in step 2, the orbital element differences between the spacecraft and target at the rendezvous moment are recalculated as [-102.72m, -4.32e-4, 5.36e-4, 1.0e-4, -6.69e-5, 0.010]. The differences in each component are very small, demonstrating that high-precision orbit rendezvous can be achieved with a single correction. To further improve accuracy, the correction calculation process from steps 2 to 5 can be repeated with the corrected pulse size until the predicted orbital element difference is less than the desired threshold.

[0199] References

[0200] [1] Luo Yazhong, Huang Anyi, Li Hengnian, Wu Shenggang, Zhang Jin, Yang Zhen. A fast estimation method for optimal velocity increment of long-time orbit rendezvous under J2 perturbation[P]. Hunan Province: CN110789739A, 2020-02-14.

Claims

1. A semi-analytical multi-pulse elliptical orbit out-of-plane rendezvous optimization method, characterized in that: Please follow the steps below to implement: Step 1: Calculate the parameters that can represent the approximately optimal rendezvous four-pulse according to the input initial orbit element number and the target orbit element number to be rendezvoused, as well as the four-pulse estimated values ​​and corresponding four phase values ​​directly calculated from these parameters; Step 2: Calculate the time to apply the four pulses based on the four phase values, and then predict the spacecraft at the specified rendezvous time t f The number of orbital elements; Step 3: Calculate the difference between the orbital elements predicted in step 2 and the orbital elements of the given target orbit, and record it in the form of first-kind singularity-free elements; Step 4: Calculate the correction amount for the semi-major axis control amount based on the first and sixth components of the orbital root difference; Step 5: Calculate the corrections to the inclination and right ascension of the ascending node based on the fourth and fifth components of the orbital element difference. Step 6: Calculate the correction amount for the eccentricity control amount based on the second and third components of the orbital element difference; Step 7: Substitute the corrected control quantity back into step 1, recalculate the size and phase of the four pulses, and then substitute them back into steps 2 and 3. Predict the updated orbital element difference after parameter correction and determine whether it is less than the given threshold. If so, output the size and timing of the four pulses as the optimal control quantity. Otherwise, continue to use steps 4, 5, and 6 to correct the corresponding control quantity until the predicted orbital element difference is less than the given threshold.

2. A semi-analytical multi-pulse elliptical orbit out-of-plane intersection optimization method according to claim 1, characterized in that: The step 1 is specifically implemented according to the following steps: According to the input spacecraft initial orbital elements [a0, e0, i0, Ω0, ω0, M0], a represents the semi-major axis, e represents the eccentricity, i represents the inclination, Ω represents the right ascension of the ascending node, ω represents the argument of perigee, M represents the mean anomaly, subscript 0 represents the orbital elements of the spacecraft at the initial moment, and the target orbital elements that need to be rendezvoused, f is the subscript representing the target orbit, the orbital elements corresponding to the spacecraft and target orbits at the rendezvous moment are predicted, and the difference is taken to obtain the orbital element difference described in the form of the first kind of non-singular point elements, which is recorded as [Δa0, Δe x0 ,Δe y0 ,Δi0,ΔΩ0,Δu0], defined as: Δa0=a f -a0 No x0 =e f I'm sorry. f -e0 cosω0 No y0 =e f sin f -e0 sinω0 Δi0=i f -i0 ΔΩ0=Ω f -Ω0 Δu0=ω f +M f -ω0-M0 According to [Δa0,Δe x0 ,Δe y0 ,Δi0,ΔΩ0,Δu0], the optimal four-pulse intersection size calculated is described as formula (1), where the first two pulses act on the first circle, the last two pulses act on the last circle, and the sizes of the four pulse velocity increments Δv1, Δv2, Δv3, and Δv4 are respectively: Among them, Δv1, Δv2, Δv3, Δv4 represent four pulse velocity increments, a=a0 is the initial semi-major axis, i=i0 is the initial inclination angle, represents the initial orbital average velocity, μ is the Earth's gravitational constant, Δa represents the change in the semi-major axis caused by the first two pulses, Δi represents the change in the inclination caused by the first two pulses, ΔΩ / 2 represents the change in the right ascension of the ascending node caused by the first two pulses, and k r1 k represents the change in eccentricity of the radial component in the first two pulses. r2 represents the change of the eccentricity due to the radial component in the last two pulses, k0 represents the change of the eccentricity due to the tangential component in the first two pulses, and k f It represents the change of eccentricity due to the tangential component in the last two pulses; The phase angles u1, u2, u3, and u4 of the four pulse maneuvering moments corresponding to Δv1, Δv2, Δv3, and Δv4 are: Among them, u1 and u2 correspond to a certain moment in the first circle after the start of orbit transfer, and u3 and u4 correspond to a certain moment in the last circle before rendezvous; according to the definition of orbital elements, u1, u2, u3, u4 are one-to-one corresponding to t1, t2, t3, t4, so the optimal four-pulse results are expressed by Δv1, Δv2, Δv3, Δv4 and u1, u2, u3, u4.

3. A semi-analytical multi-pulse elliptical orbit out-of-plane intersection optimization method according to claim 1, characterized in that: The step 2 is specifically implemented according to the following steps: Step 2.1, calculate the time t1 when the first pulse is applied: Assuming that the initial phase of the spacecraft at time t0 is u0=ω0+M0(3), the phase difference between u0 and u1 includes four cases: Formula (4) ensures that Δu1 is always between 0 and π by treating different values ​​of Δu1 differently. According to the near-circular orbit assumption, the time of the first pulse relative to the initial moment is in is the average angular velocity of the orbit corresponding to the number of orbital elements at the initial moment; Step 2.2, calculate the time t2 when the second pulse is applied: The orbital elements of the spacecraft at the initial moment are predicted to the moment t1=t0+Δt1. At this moment, the velocity vector of the spacecraft is added with the velocity increment vector Δv1. The three-dimensional components of Δv1 in the radial, tangential and normal directions are: in: e t 、e r 、e n Represent the tangential unit vector, radial unit vector and normal unit vector respectively, r=[x,y,z], v=[v x ,v y ,v z ] are the position and velocity vectors of the spacecraft at time t1, r = |r|, v = |v| are the moduli of the position and velocity vectors; At time t1, adding Δv1 to the spacecraft velocity vector gives the position velocity of the spacecraft after the first pulse is applied, which can be converted into six orbital numbers. According to formula (2), the second pulse phase u2 is: u2=u1+π(8), Therefore, the time difference Δt2 of the second pulse relative to the first pulse time t1 is: Where n1 is the orbital angular velocity after the first pulse. The time when the second pulse is applied is then obtained as t2 = Δt2 + t1. The spacecraft orbit is predicted to time t2, and the second pulse Δv2 is added to the velocity vector: The definition of variables in formula (10) is the same as formula (6). The definitions of case 1, case 2, case 3, and case 4 are consistent with formula (4). At this time, the position velocity of the spacecraft after the first and second pulses are applied is obtained; Step 2.3, calculate the third pulse time t3: Since the third and fourth pulses are in the last circle, that is, the intersection time t f Nearby, first predict the orbit after applying the second pulse to t f At this moment, we get the orbital element number [a(t f ),e(t f ),i(t f ),Ω(t f ),ω(t f ),M(t f )], then t f The phase of the spacecraft at time u f for: u f =ω f +M f (11) According to formula (2), the phases of the third and fourth pulses are u3, u4 = u3 + π, and u4 and u f The difference between u3 and u f The phase difference Δu3 must be between π and 2π. Calculate Δu3 according to the following formula: Where mod(x, y) represents x modulo y; then the third pulse time t3 is: Predict the spacecraft orbit to time t3 and add the third pulse Δv3 to the velocity vector: The definitions of Case 1, Case 2, Case 3, and Case 4 are consistent with those of Formula (12); e t ,e n ,e r The calculation of is the same as formula (7); Step 2.4: Calculate the time t4 and t f The number of spacecraft orbital elements after four pulses have been applied at time t; The time difference Δt4 of the fourth pulse relative to the third pulse is: Where n3 is the orbital angular velocity after the third pulse. Predict the spacecraft orbit to time t4 and add the fourth pulse Δv4 to the velocity vector: The definitions of Case 1, Case 2, Case 3, and Case 4 are consistent with those of Formula (12); e t ,e n ,e r The calculation is the same as formula (7); get the new position velocity and predict t f The number of spacecraft orbital elements at the time is denoted as [a p ,e p ,i p ,Ω p ,ω p ,M p ].

4. A semi-analytical multi-pulse elliptical orbit out-of-plane intersection optimization method according to claim 3, characterized in that: The step 3 is specifically implemented according to the following steps: Define the orbital element difference [Δa p ,Δe xp ,Δe yp ,Δi p ,ΔΩ p ,Δu p ] is the difference between the orbit of the spacecraft to which the four pulses are applied and the target orbit to be rendezvoused, expressed in the form of first-kind singularity-free elements: The smaller the difference is, the smaller the parameters [Δa, Δi, ΔΩ, k0, k f ,k r1 ,k r2 ] and the four-pulse estimated values ​​Δv1, Δv2, Δv3, and Δv4 calculated directly from these seven parameters are closer to the exact values. p ,Δe xp ,Δe yp ,Δi p ,ΔΩ p ,Δu p ] for [Δa,Δi,ΔΩ,k0,k f ,k r1 ,k r2 ] to make corrections.

5. A semi-analytical multi-pulse elliptical orbit out-of-plane intersection optimization method according to claim 4, characterized in that: The step 4 is specifically implemented according to the following steps: According to Δa p ,Δu p ,Δi p ,ΔΩ p , calculate the correction amount for Δa as follows: Assume that the correction amount of Δa is Δa1, then Δa1 and Δu p The following relationship must be satisfied: Where a0 and n0 have been defined in equation (5), Δt = t f -t0 represents the time difference from the start to the intersection; At the same time, when recalculating Δv1, Δv2, Δv3, and Δv4, the term (Δa0-Δa) in equations (16) and (14) needs to be corrected to: (Δa0-Δa p )-(Δa+Δa1)(19) This step obtains the correction amount for Δa and (Δa-Δa0).

6. A semi-analytical multi-pulse elliptical orbit out-of-plane intersection optimization method according to claim 5, characterized in that: The step 5 is specifically implemented according to the following steps: Keeping u1, u2, u3, and u4 unchanged, let Δi in equations (6) and (10) be corrected to Δi+Δi1, and Δi0-Δi in equations (16) and (14) be corrected to Δi0-Δi+Δi2, then ΔΩ in equations (6) and (10) will be proportionally changed to The ΔΩ in formula (16) and formula (14) will be proportionally changed to Δi+Δi1, Δi0-Δi+Δi2, and ΔΩ should satisfy the following relationship: in It represents the right ascension drift rate of the ascending node after applying one or two pulses to the initial orbit. represents the change in the right ascension drift rate of the ascending node caused by Δa1 and Δi1; according to the law of the J2 perturbation term in orbital mechanics, it is: Combining formula (2), formula (20) is transformed into: Formula (22) is a system of linear equations in two variables, which can be solved directly: Then, the correction values ​​of Δi, Δi0-Δi, and ΔΩ used to calculate Δv1, Δv2, Δv3, and Δv4 in equations (6), (10), (14), and (16) are obtained.

7. A semi-analytical multi-pulse elliptical orbit out-of-plane intersection optimization method according to claim 6, characterized in that: The step 6 is specifically implemented according to the following steps: Let k0,k f ,k r1 ,k r2 The correction amount is Δk0, Δk f ,Δk r1 ,Δk r2 The eccentricity change caused by the correction should offset the Δe calculated by formula (17) xp ,Δe yp ,Right now: Where Δu1 represents the orbital perigee argument drift from the time the first and second pulses are applied until the last revolution: Where J2 is the second-order term of the non-spherical perturbation of the earth, R e is the Earth's equatorial radius, Formula (24) can be written as: Then Δk0, Δk f ,Δk r1 ,Δk r2 The solution is considered as Δk0, Δk f ,Δk r1 ,Δk r2 To optimize the variables, use Equation (26) as the equality constraint to minimize the sum of Δv1, Δv2, Δv3, and Δv4: minf=|Δv1|+|Δv2|+|Δv3|+|Δv4| Obtain Δk0, Δk that minimizes Equation (27) f ,Δk r1 ,Δk r2 .

8. A semi-analytical multi-pulse elliptical orbit out-of-plane intersection optimization method according to claim 7, characterized in that: The step 7 is specifically as follows: The corrected Δa, Δi, ΔΩ, k0, k f ,k r1 ,k r2 Substitute back into step 1, recalculate Δv1, Δv2, Δv3, Δv4 and u1, u2, u3, u4, then substitute back into steps 2 and 3, and predict the updated orbital root difference [Δa p ,Δe xp ,Δe yp ,Δu p ,Δi p ,ΔΩ p ], judge whether it is less than a given threshold, if so, output Δv1, Δv2, Δv3, Δv4 and t1, t2, t3, t4; otherwise, continue with steps 4, 5, and 6 to compare Δa, Δi, ΔΩ, k0, k f ,k r1 ,k r2 Correction is made until the forecast [Δa p ,Δe xp ,Δe yp ,Δu p ,Δi p ,ΔΩ p ] is less than a given threshold and stops.

Images