An analytical optimization method for low-thrust circular orbit out-of-plane rendezvous

Abstract

The present invention discloses an analytical low-thrust circular orbit out-of-plane rendezvous optimization method, which first obtains the input spacecraft initial orbit and the target orbit that needs to be rendezvoused, as well as the specified rendezvous time, and calculates the orbital root difference between the spacecraft and the target orbit that needs to be rendezvoused at the given rendezvous time without thrust; obtains a simplified model that describes the low-thrust circular orbit out-of-plane rendezvous trajectory; then establishes a fixed thrust startup coefficient, optimizes to obtain the optimal thrust startup coefficient and the semi-major axis, inclination, and right ascension change of the ascending node; obtains the phase difference of the rendezvous time relative to the target orbit that needs to be rendezvoused; and then obtains the new inclination and right ascension change of the ascending node; and finally outputs the corresponding low-thrust control rate and rendezvous transfer orbit. The present invention solves the problems of time-consuming integration and too many variables to be solved in the parameter optimization method in the prior art, and the problem that it is not easy to obtain the optimal rendezvous trajectory.

Description

Technical Field

[0001] The invention belongs to the technical field of aerospace navigation control, and in particular relates to an analytical small-thrust circular orbit out-of-plane rendezvous optimization method. Background Art

[0002] In space operations or debris removal missions, spacecraft need to rendezvous with targets and carry out various operations in the most fuel-efficient way. Electric propulsion has the advantage of high specific impulse compared to traditional chemical propulsion, but its thrust is small, the orbit transfer time is long, and there are certain difficulties in optimization calculations. Since the indirect optimization method can only meet the necessary conditions for optimality, it is easy to fall into local optimality, and it takes a long time to repeatedly integrate the orbit. The existing parameter optimization methods also have the problems of time-consuming integration and too many variables to solve, and it is not easy to obtain the optimal rendezvous trajectory. Summary of the invention

[0003] The purpose of the present invention is to provide an analytical low-thrust circular orbit out-of-plane rendezvous optimization method, which solves the problems of time-consuming integration and too many variables to be solved in the prior art parameter optimization method and the difficulty in obtaining the optimal rendezvous trajectory.

[0004] The technical solution adopted by the present invention is an analytical low-thrust circular orbit out-of-plane intersection optimization method, which is specifically implemented according to the following steps:

[0005] Step 1: First, obtain the input initial orbit of the spacecraft, the target orbit to be rendezvoused, and the specified rendezvous time, and calculate the orbital element difference between the spacecraft and the target orbit to be rendezvoused at the given rendezvous time without thrust;

[0006] Step 2, design an optimal thrust on / off strategy, which is divided into three stages during the entire rendezvous flight time: scheduled on / off at the beginning of the mission, no on / off in the middle of the mission, and scheduled on / off before the rendezvous moment; define a thrust on / off coefficient to represent the ratio of the thrust on / off time in the first and third stages to the duration of the first and third stages themselves; obtain a simplified model to describe the low-thrust circular orbit out-of-plane rendezvous trajectory;

[0007] Step 3: According to the dynamic law and the on / off strategy and model defined in step 2, a fixed thrust start-up coefficient is established, and the semi-major axis, inclination angle, and right ascension of the ascending node generated by the first-stage thrust defined in step 2 are used as optimization variables, and an equation constraint optimization model with the minimum small thrust equivalent speed increment as the goal is to discretize the thrust start-up coefficient between 0 and 1 and solve the corresponding equation constraint optimization model respectively, and find the optimal thrust start-up coefficient and the change of semi-major axis, inclination angle, and right ascension of the ascending node;

[0008] Step 4: Convert the optimal thrust start-up coefficient and the changes in the semi-major axis, inclination angle, and right ascension of the ascending node into the actual control rate of small thrust and substitute them into the orbit prediction model to predict and calculate the orbital elements of the spacecraft at the rendezvous time. Obtain the phase difference of the rendezvous time relative to the target orbit that needs to be rendezvoused;

[0009] Step 5, calculate the correction amount of the semi-major axis change according to the phase difference, then fix the semi-major axis change to the corrected value, re-solve the equation constraint optimization model in step 3, and obtain new changes in inclination and ascending node right ascension;

[0010] Step 6, return to step 4 to re-predict the orbital elements at the rendezvous time and the phase difference relative to the target orbit, and then determine whether the phase difference meets the given threshold range. If not, go to step 5 to continue to solve the new semi-major axis, inclination, and right ascension change of the ascending node and predict the phase difference; until the phase difference is less than the threshold, the corresponding low thrust control rate and rendezvous transfer orbit are output.

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

[0012] Step 1 is implemented as follows:

[0013] Assume that the input initial spacecraft orbit elements and the rendezvous target orbit elements are [a0, e0, i0, Ω0, ω0, M0] and [a f ,e f ,i f ,Ω f ,ω f ,M f ], where a is the semi-major axis, e is the eccentricity, i is the inclination, ω is the argument of perigee, Ω is the right ascension of the ascending node, M is the mean anomaly, subscript 0 indicates the initial orbit, and f indicates the target orbit. The difference between the two is expressed as [△a0,△i0,△Ω0,△u0] = [a f -a0,i f -i0,Ω f -Ω0,(ω f +M f -ω0-M0)].

[0014] Step 2 is implemented as follows:

[0015] Step 2.1, timed on / off means that the ratio of each power-on time to the single-circuit period is fixed as a variable to be optimized k thrust (values ​​between 0 and 1), the thrust direction only contains tangential and normal components, and the tangential component of the thrust direction remains unchanged, and the normal component switches its sign every half a circle.

[0016] u,u on ,u offThey are the equivalent action midpoint phase of half-turn thrust, the start-up point phase and the shutdown point phase, and have the following relationship: Assume α is the angle between the thrust direction and the tangential direction, then the tangential component of the thrust F is Fcosα, the normal component is F sinα in the first half circle, and -F sinα in the second half circle, and assume that the spacecraft period and mass remain approximately unchanged, recorded as constants T and m respectively, then the cumulative change of the thrust per circle to the orbital element number is written as:

[0017]

[0018] where △a T Indicates the change in the semi-major axis caused by the thrust in time T, △i T Indicates the change in inclination angle caused by thrust within T time, △Ω T It indicates the change in the right ascension of the ascending node caused by the thrust in time T. t ,△v n are the tangential and normal velocity increments generated by the thrust per revolution, a0 and i0 are the initial semi-major axis and inclination of the spacecraft, and V is the initial orbital average velocity.

[0019] remember represents the arc loss coefficient of the normal velocity increment accumulation, and assuming that the total time length of the timed on / off machine in the initial stage of the intersection is △t1, the equivalent midpoint phase is u1, the total time length of the timed on / off machine in the stage approaching the intersection is △t2, and the equivalent midpoint phase is u2, then according to formula (2), the change in the number of orbital elements corresponding to △t1 and △t2 is:

[0020]

[0021] Among them, △a1, △i1, △Ω1 represent the changes of thrust on semi-major axis, inclination and right ascension of ascending node in the initial stage of rendezvous respectively; △a2, △i2, △Ω2 represent the changes of thrust on semi-major axis, inclination and right ascension of ascending node in the stage approaching rendezvous respectively; thus, the relationship between the designed thrust on / off strategy and the change of orbital elements is obtained. This model will be used in the optimization model in the next section, among which △a1, △i1, △Ω1, △a2, △i2, △Ω2 and k thrust , α1, △t1, u1, α2, △t2, u2 are all values ​​to be optimized.

[0022] Step 3 is implemented as follows:

[0023] Step 3.1: Given k thrust In order to meet the intersection requirements, the orbital element changes △a1, △i1, △Ω1, △a2, △i2, △Ω2 calculated in step 3 satisfy the following constraints:

[0024]

[0025] in, J2,R is the drift rate of the spacecraft ascending node right ascension relative to the target orbit ascending node right ascension after the first stage of power on and off, e , μ are the J2 term of the non-spherical perturbation of the Earth, the radius of the Earth and the gravitational constant respectively.

[0026] Step 3.2: For a fixed k thrust , an equality constraint optimization model is established with △a1,△i1,△Ω1,△Ω2 as variables and the minimum total thrust start-up time as the goal:

[0027]

[0028] The parameters have been defined above. Combined with equation (4), the equality constraint function g is defined as:

[0029]

[0030] in is the right ascension drift rate of the ascending node corresponding to the initial orbit.

[0031] Step 3.3, order y=△i1,z1=△Ω1,z2=△Ω2, then according to the minimum principle, the optimal x,y,z1,z2 should satisfy the following relationship:

[0032]

[0033] where λ represents the constraint multiplier, Respectively represent the derivatives of △t1,△t2 with respect to x, y, z1, z2:

[0034]

[0035] Solve this equation, k in equation (7) thrust is a fixed value. In order to obtain the optimal k thrust , respectively take k thrust Let f be 50 numbers with equal intervals between 0 and 1 and solve equation (7) to find the k corresponding to the minimum f. thrust And x, y, z1, z2 are the optimal solutions. The thrust switch control variables △t1, △t2, α1, α2, u1, u2 are calculated by formula (3).

[0036] Step 4 is implemented according to the following steps:

[0037] Step 4.1, Substitute △t1,△t2,α1,α2,u1,u2 into the small thrust precise dynamic equation corresponding to equation (9) to predict the orbital elements at the rendezvous time. The dynamic equation is:

[0038]

[0039] Among them, r and v are the spacecraft position velocity vectors, which can be converted to orbital roots; a p is the sum of the perturbation accelerations, e t ,e n are the tangent and normal unit vectors in the orbital coordinate system, respectively, and the calculation method is common sense; c is the thrust on and off flag that changes with time, 0 is off, and 1 is on; α is the thrust yaw angle that changes with time. The specific values ​​of c and α are determined by the parameters △t1, △t2, α1, α2, u1, and u2. The specific relationship is divided into two cases as follows.

[0040] When time t≤t0+△t1:

[0041]

[0042] When t ≥ t f -△t2:

[0043]

[0044] The variables in the formula have been defined;

[0045] Step 4.2: Use the numerical integration method to calculate the position velocity and orbital element [a] of the spacecraft at the rendezvous moment according to the initial orbital element of the spacecraft and equation (9). p ,e p ,i p ,Ω p ,ω p ,M p ], calculate the phase difference △u with the target orbital element p =ω f +M f -ω p -M p .

[0046] Step 5 is implemented according to the following steps:

[0047] According to △u p Calculate the correction amount for △a1

[0048]

[0049] where μ is the Earth's gravitational constant. Then fix Resolve the remaining three variables y, z1, z2 in the optimization model (6); then, according to formula (3), the values ​​of the parameters △t1, △t2, α1, α2, u1, u2 can be recalculated.

[0050] Step 6 is implemented according to the following steps:

[0051] Substitute △t1,△t2,α1,α2,u1,u2 obtained in step 6 back into step 5, re-predict the orbital elements of the spacecraft and calculate △u p .

[0052] Determine △u p Is it less than a given threshold? If so, output △t1,△t2,α1,α2,u1,u2 and the corresponding predicted orbital elements, which are the position, velocity and orbital elements of the spacecraft during flight, which are the optimal rendezvous trajectory with low thrust. Otherwise, return to step 6 to continue to correct △a1 and continue to calculate the new △u p , until △u p is less than a given threshold.

[0053] The beneficial effect of the present invention is that an analytical low-thrust circular orbit out-of-plane rendezvous optimization method first calculates the orbital root difference between the spacecraft and the target at the rendezvous time when the spacecraft is not maneuvering; then designs an analytical low-thrust on-off mode, establishes an equation-constrained optimization model for the low-thrust startup parameters; and obtains the optimal startup parameters and transfer trajectory of the low-thrust rendezvous by solving a nonlinear equation group and performing iteration and local optimization. The method disclosed by the present invention solves the problem of long time consumption and large amount of calculation in solving long-term low-thrust orbit rendezvous, is applicable to the trajectory optimization of circular orbit out-of-plane rendezvous, and significantly improves the solution efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 It is a schematic diagram of the thrust on / off strategy;

[0055] Figure 2 is a flow chart of the method of the present invention;

[0056] Figure 3 It is a graph of thrust variation with time during rendezvous flight;

[0057] Figure 4 is the graph of the variation of the semi-major axis difference over time during the intersection process;

[0058] Figure 5 is the graph of the change of the inclination difference over time during the rendezvous process;

[0059] Figure 6 This is the graph of the right ascension difference of the ascending node during the conjunction process changing with time;

[0060] Figure 7 It is a graph showing the phase difference changing with time during the intersection process. DETAILED DESCRIPTION

[0061] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0062] In view of the characteristics of multi-turn low-thrust orbit transfer, a fast optimization method with fewer parameters and solveable using an equality constraint optimization model is established to solve the problem of fast solution of long-term orbit rendezvous. In order to achieve the above technical objectives, the technical solution adopted by the present invention is: an analytical low-thrust circular orbit out-of-plane rendezvous optimization method, that is, given the initial orbit elements of the spacecraft [a0, e0, i0, Ω0, ω0, M0] and the target orbit [a f ,e f ,i f ,Ω f ,ω f ,M f ], given the intersection time t f The process of rendezvous with a circular orbit with a fixed thrust of F and a duration of △t to obtain the optimal thrust on / off control rate and transfer trajectory (position speed change during the rendezvous process).

[0063] An analytical optimization method for low-thrust circular orbit out-of-plane rendezvous Figures 1-2 , follow the steps below:

[0064] Step 1: First, obtain the input initial orbit of the spacecraft, the target orbit to be rendezvoused, and the specified rendezvous time, and calculate the orbital element difference between the spacecraft and the target orbit to be rendezvoused at the given rendezvous time without thrust;

[0065] Step 1 is implemented as follows:

[0066] Assume that the input initial spacecraft orbit elements and the rendezvous target orbit elements are [a0, e0, i0, Ω0, ω0, M0] and [a f ,e f ,i f ,Ω f ,ω f ,M f ], where a is the semi-major axis, e is the eccentricity, i is the inclination, ω is the argument of perigee, Ω is the right ascension of the ascending node, M is the mean anomaly, subscript 0 indicates the initial orbit, and f indicates the target orbit. The difference between the two is expressed as [△a0,△i0,△Ω0,△u0] = [a f -a0,i f -i0,Ω f -Ω0,(ω f +M f -ω0-M0)].

[0067] Step 2: Design an optimal on-off strategy for the thrust. During the entire rendezvous flight time, it is divided into three stages: timed on-off at the beginning of the mission, no power-on in the middle stage of the mission, and timed on-off before approaching the rendezvous moment; define a thrust power-on coefficient to represent the ratio of the thrust power-on time to the duration of the first and third stages themselves in the first and third stages; obtain a simplified model describing the non-coplanar rendezvous trajectory of a small-thrust circular orbit.

[0068] Step 2 is specifically implemented according to the following steps:

[0069] Step 2.1: Timed on-off means that the ratio of the on-time duration per orbit to the single-orbit period is fixed as an optimization variable k thrust (with a value between 0 and 1). The thrust direction only includes tangential and normal components, and the tangential component of the thrust direction remains unchanged, while the normal component changes its sign every half orbit, as shown in the appendix Figure 1 . In the appendix Figure 1 u, u on , u off are the equivalent action midpoint phase, power-on point phase, and power-off point phase of the half-orbit thrust respectively, and there are the following relationships: Let α be the angle between the thrust direction and the tangent. Then the tangential component of the thrust F is F cosα, the normal component is F sinα in the upper half orbit and -F sinα in the lower half orbit. Assuming that the spacecraft period and mass are approximately constant, denoted as constant values T and m respectively, the cumulative change in orbital elements per orbit caused by the thrust is written as:

[0070]

[0071] where △a T represents the change in semi-major axis caused by the thrust within time T, △i T represents the change in inclination caused by the thrust within time T, △Ω T represents the change in right ascension of the ascending node caused by the thrust within time T. △v t , △v n are the tangential and normal velocity increments generated by the thrust per orbit respectively, a0 and i0 are the initial semi-major axis and inclination of the spacecraft, and V is the initial orbital average velocity.

[0072] Denote as the arc loss coefficient of the cumulative normal velocity increment. Let the total duration of timed on-off at the beginning of the defined rendezvous stage be △t1, the equivalent midpoint phase be u1, the total duration of timed on-off before approaching the rendezvous stage be △t2, and the equivalent midpoint phase be u2. Then, according to Equation (2), the changes in orbital elements corresponding to △t1 and △t2 are:

[0073]

[0074] Among them, △a1, △i1, △Ω1 represent the changes of thrust on semi-major axis, inclination and right ascension of ascending node in the initial stage of rendezvous respectively; △a2, △i2, △Ω2 represent the changes of thrust on semi-major axis, inclination and right ascension of ascending node in the stage approaching rendezvous respectively; thus, the relationship between the designed thrust on / off strategy and the change of orbital elements is obtained. This model will be used in the optimization model in the next section, among which △a1, △i1, △Ω1, △a2, △i2, △Ω2 and k thrust , α1, △t1, u1, α2, △t2, u2 are all values ​​to be optimized.

[0075] Step 3: According to the dynamic law and the on / off strategy and model defined in step 2, a fixed thrust start-up coefficient is established, and the semi-major axis, inclination angle, and right ascension of the ascending node generated by the first-stage thrust defined in step 2 are used as optimization variables, and an equation constraint optimization model with the minimum small thrust equivalent speed increment as the goal is to discretize the thrust start-up coefficient between 0 and 1 and solve the corresponding equation constraint optimization model respectively, and find the optimal thrust start-up coefficient and the change of semi-major axis, inclination angle, and right ascension of the ascending node;

[0076] Step 3 is implemented as follows:

[0077] Step 3.1: Given k thrust In order to meet the intersection requirements, the orbital element changes △a1, △i1, △Ω1, △a2, △i2, △Ω2 calculated in step 3 satisfy the following constraints:

[0078]

[0079] in, J2,R is the drift rate of the spacecraft ascending node right ascension relative to the target orbit ascending node right ascension after the first stage of power on and off, e , μ are the J2 term of the non-spherical perturbation of the Earth, the radius of the Earth and the gravitational constant respectively.

[0080] Step 3.2: For a fixed k thrust , an equality constraint optimization model is established with △a1,△i1,△Ω1,△Ω2 as variables and the minimum total thrust start-up time as the goal:

[0081]

[0082] The parameters have been defined above. Combined with equation (4), the equality constraint function g is defined as:

[0083]

[0084] in is the right ascension drift rate of the ascending node corresponding to the initial orbit.

[0085] Step 3.3, order y=△i1,z1=△Ω1,z2=△Ω2, then according to the minimum principle, the optimal x,y,z1,z2 should satisfy the following relationship:

[0086]

[0087] where λ represents the constraint multiplier, Respectively represent the derivatives of △t1,△t2 with respect to x, y, z1, z2:

[0088]

[0089] Solve this equation, k in equation (7) thrust is a fixed value. In order to obtain the optimal k thrust , respectively take k thrust Let f be 50 numbers with equal intervals between 0 and 1 and solve equation (7) to find the k corresponding to the minimum f. thrust And x, y, z1, z2 are the optimal solutions. The thrust switch control variables △t1, △t2, α1, α2, u1, u2 are calculated by formula (3).

[0090] Step 4: Convert the optimal thrust start-up coefficient and the changes in the semi-major axis, inclination angle, and right ascension of the ascending node into the actual control rate of small thrust and substitute them into the orbit prediction model to predict and calculate the orbital elements of the spacecraft at the rendezvous time. Obtain the phase difference of the rendezvous time relative to the target orbit that needs to be rendezvoused;

[0091] Step 4 is implemented according to the following steps:

[0092] Step 4.1, Substitute △t1,△t2,α1,α2,u1,u2 into the small thrust precise dynamic equation corresponding to equation (9) to predict the orbital elements at the rendezvous time. The dynamic equation is:

[0093]

[0094] Among them, r and v are the spacecraft position velocity vectors, which can be converted to orbital roots; a p is the sum of the perturbation accelerations, e t ,e n are the tangent and normal unit vectors in the orbital coordinate system, respectively, and the calculation method is common sense; c is the thrust on and off flag that changes with time, 0 is off, and 1 is on; α is the thrust yaw angle that changes with time. The specific values ​​of c and α are determined by the parameters △t1, △t2, α1, α2, u1, and u2. The specific relationship is divided into two cases as follows.

[0095] When time t≤t0+△t1:

[0096]

[0097] When t ≥ t f -△t2:

[0098]

[0099] The variables in the formula have been defined;

[0100] Step 4.2: Use the numerical integration method to calculate the position velocity and orbital element [a] of the spacecraft at the rendezvous moment according to the initial orbital element of the spacecraft and equation (9). p ,e p ,i p ,Ω p ,ω p ,M p ], calculate the phase difference △u with the target orbital element p =ω f +M f -ω p -M p .

[0101] Step 5, calculate the correction amount of the semi-major axis change according to the phase difference, then fix the semi-major axis change to the corrected value, re-solve the equation constraint optimization model in step 3, and obtain new changes in inclination and ascending node right ascension;

[0102] Step 5 is implemented according to the following steps:

[0103] According to △u p Calculate the correction amount for △a1

[0104]

[0105] where μ is the Earth's gravitational constant. Then fix Resolve the remaining three variables y, z1, z2 in the optimization model (6); then, according to formula (3), the values ​​of the parameters △t1, △t2, α1, α2, u1, u2 can be recalculated.

[0106] Step 6, return to step 4 to re-predict the orbital elements at the rendezvous time and the phase difference relative to the target orbit, and then determine whether the phase difference meets the given threshold range. If not, go to step 5 to continue to solve the new semi-major axis, inclination, and right ascension change of the ascending node and predict the phase difference; until the phase difference is less than the threshold, the corresponding low thrust control rate and rendezvous transfer orbit are output.

[0107] Step 6 is implemented according to the following steps:

[0108] Substitute △t1,△t2,α1,α2,u1,u2 obtained in step 6 back into step 5, re-predict the orbital elements of the spacecraft and calculate △u p .

[0109] Determine △u p Is it less than a given threshold? If so, output △t1,△t2,α1,α2,u1,u2 and the corresponding predicted orbital elements, which are the position, velocity and orbital elements of the spacecraft during flight, which are the optimal rendezvous trajectory with low thrust. Otherwise, return to step 6 to continue to correct △a1 and continue to calculate the new △u p , until △u p is less than a given threshold.

[0110] Example

[0111] Assume that the initial mass of the spacecraft is 500kg, the specific impulse is 1000s, and the thrust is 0.5N. It is required to rendezvous with the target orbit after 14 days. The elements of the local orbit and the target orbit are as follows

[0112] Table 1 Epochs for local and target orbits at initial moments

[0113] A(km) e i(degree) Ω(degrees) ω(degrees) M(degrees) This star 7100 0 98 0 0 0 Target 7140 0 98.1 1 0 0

[0114] The solution steps are as follows: (1) First, calculate the orbital difference between the local satellite and the target at the time of rendezvous: [△a0,△e x0 ,△e y0 ,△i0,△Ω0,△u0]=[40km,0,0,0.1 degree,1 degree,0,0]. (2) k thrust Discretize into 51 numbers with an interval of 0.02 between 1 and 0, solve the equations of formula (7) in sequence, and get k thrust = 0.38, the equivalent velocity increment is 93.78 m / s, which is the optimal value. (3) The calculated values ​​of u1, α1, u2, α2, △t1, and △t2 are [0.7931, 2.0756, 2.0050, 0.9559, 1.3019 days, and 1.5547 days], respectively. (4) The predicted calculated orbital elements at the time of rendezvous and the orbital elements of the target have a phase difference △u p =-0.8076. (5) According to formula (12), the corrected x=x+△a1 is calculated c , re-solve equation (7) to obtain p =0 for [u1,α1,u2,α2,△t1,△t2] = [0.7961, 2.0255, 2.0141, 0.9181, 1.3072 days, 1.5300 days].

[0115] (6) Output the orbital root number change and thrust on / off detailed information corresponding to the optimal solution according to formula (9). The change trend of thrust magnitude and direction is shown in the attached figure. Figure 3 The changes in the orbital element difference between the spacecraft and the target are shown in the attached Figures 4 to 7 Compared with other methods, the optimization time of the present invention is greatly reduced (the method of the present invention takes only less than 1 second, while the traditional indirect method takes more than 20 minutes and cannot stably converge to the global optimal solution); the fuel consumption of the two methods differs by only 5%; and the thrust direction angle α corresponding to the present invention is a constant value within the interval, which is easier to implement in engineering than the thrust direction angle calculated by the indirect method, which changes continuously.

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

Claims

1. An analytical low-thrust circular orbit out-of-plane intersection optimization method, characterized in that: Follow the steps below to implement it: Step 1: First, obtain the input initial orbit of the spacecraft, the target orbit to be rendezvoused, and the specified rendezvous time, and calculate the orbital element difference between the spacecraft and the target orbit to be rendezvoused at the given rendezvous time without thrust; Step 2, design an optimal thrust on / off strategy, which is divided into three stages during the entire rendezvous flight time: scheduled on / off at the beginning of the mission, no on / off in the middle of the mission, and scheduled on / off before the rendezvous moment; define a thrust on / off coefficient to represent the ratio of the thrust on / off time in the first and third stages to the duration of the first and third stages themselves; obtain a simplified model to describe the low-thrust circular orbit out-of-plane rendezvous trajectory; Step 3: According to the dynamics law and the on / off strategy and model defined in step 2, a fixed thrust start-up coefficient is established. The semi-major axis, inclination angle, and right ascension of the ascending node generated by the first-stage thrust defined in step 2 are used as optimization variables, and an equation-constrained optimization model is established with the goal of minimizing the small thrust equivalent speed increment. The thrust start-up coefficient is discretized between 0 and 1 and the corresponding equation-constrained optimization model is solved respectively to obtain the optimal thrust start-up coefficient and the semi-major axis, inclination angle, and right ascension of the ascending node. Step 4: Convert the optimal thrust start-up coefficient and the changes in the semi-major axis, inclination angle, and right ascension of the ascending node into the actual control rate of the small thrust and substitute them into the orbit prediction model to predict and calculate the orbital elements of the spacecraft at the rendezvous time, and obtain the phase difference of the rendezvous time relative to the target orbit that needs to be rendezvoused; Step 5, calculate the correction amount of the semi-major axis change according to the phase difference, then fix the semi-major axis change to the corrected value, re-solve the equation constraint optimization model in step 3, and obtain new changes in inclination and ascending node right ascension; Step 6, return to step 4 to re-predict the orbital elements at the rendezvous time and the phase difference relative to the target orbit, and then determine whether the phase difference meets the given threshold range. If not, go to step 5 to continue to solve the new semi-major axis, inclination, and right ascension change of the ascending node and predict the phase difference; until the phase difference is less than the threshold, the corresponding low thrust control rate and rendezvous transfer orbit are output.

2. The analytical low-thrust circular orbit out-of-plane intersection optimization method according to claim 1 is characterized in that: The step 1 is specifically implemented according to the following steps: Assume that the input initial spacecraft orbit elements and the rendezvous target orbit elements are [a0, e0, i0, Ω0, ω0, M0] and [a f ,e f ,i f ,Ω f ,ω f ,M f ], where a is the semi-major axis, e is the eccentricity, i is the inclination, ω is the argument of perigee, Ω is the right ascension of the ascending node, M is the mean anomaly, subscript 0 represents the initial orbit, and f represents the target orbit. The difference between the two is expressed as [△a0,△i0,△Ω0,△u0] = [a f -a0,i f -i0,Ω f -Ω0,(ω f +M f -ω0-M0)].

3. The analytical low-thrust circular orbit out-of-plane intersection optimization method according to claim 2 is characterized in that: The step 2 is specifically implemented according to the following steps: Step 2.1, timed on / off means that the ratio of each power-on time to the single-circuit period is fixed as a variable to be optimized k thrust , the thrust direction only contains the tangential component and the normal component, and the tangential component of the thrust direction remains unchanged, and the normal component switches its sign every half a circle, u,u on ,u off They are the equivalent action midpoint phase of half-turn thrust, the start-up point phase and the shutdown point phase, and have the following relationship: (1) Let α be the angle between the thrust direction and the tangential direction. Then the tangential component of the thrust F is Fcosα, and the normal component is Fsinα in the first half of the circle and -Fsinα in the second half of the circle. Assuming that the spacecraft period and mass remain approximately unchanged, denoted as constants T and m respectively, the cumulative change of the thrust per circle to the orbital element number can be written as: where △a T Indicates the change in the semi-major axis caused by the thrust in time T, △i T Indicates the change in inclination angle caused by thrust within T time, △Ω T It represents the change in the right ascension of the ascending node caused by the thrust within the time T, △v t ,△v n are the tangential and normal velocity increments generated by each rotation thrust, a0 and i0 are the initial semi-major axis and inclination of the spacecraft, V is the initial orbital average velocity, remember represents the arc loss coefficient of the normal velocity increment accumulation, and assuming that the total time length of the timed on / off machine in the initial stage of the intersection is △t1, the equivalent midpoint phase is u1, the total time length of the timed on / off machine in the stage approaching the intersection is △t2, and the equivalent midpoint phase is u2, then according to formula (2), the change in the number of orbital elements corresponding to △t1 and △t2 is: Among them, △a1, △i1, △Ω1 represent the changes of thrust on semi-major axis, inclination and right ascension of ascending node in the initial stage of rendezvous respectively; △a2, △i2, △Ω2 represent the changes of thrust on semi-major axis, inclination and right ascension of ascending node in the stage approaching rendezvous respectively; thus, the relationship between the designed thrust on / off strategy and the change of orbital elements is obtained. This model will be used in the optimization model in the next section, among which △a1, △i1, △Ω1, △a2, △i2, △Ω2 and k thrust , α1, △t1, u1, α2, △t2, u2 are all values ​​to be optimized.

4. The analytical low-thrust circular orbit out-of-plane intersection optimization method according to claim 3 is characterized in that: The step 3 is specifically implemented according to the following steps: Step 3.1: Given k thrust In order to meet the intersection requirements, the orbital element changes △a1, △i1, △Ω1, △a2, △i2, △Ω2 calculated in step 3 satisfy the following constraints: in, J2,R is the drift rate of the spacecraft ascending node right ascension relative to the target orbit ascending node right ascension after the first stage of power on and off, e ,μ are the earth's non-spherical perturbation J2 term, the earth's radius and the gravitational constant respectively; Step 3.2: For a fixed k thrust , an equality constraint optimization model is established with △a1,△i1,△Ω1,△Ω2 as variables and the minimum total thrust start-up time as the goal: The parameters have been defined above. Combined with equation (4), the equality constraint function g is defined as: in is the right ascension drift rate of the ascending node corresponding to the initial orbit; Step 3.3, order y=△i1,z1=△Ω1,z2=△Ω2, then according to the minimum principle, the optimal x,y,z1,z2 should satisfy the following relationship: where λ represents the constraint multiplier, Respectively represent the derivatives of △t1,△t2 with respect to x, y, z1, z2: Solve this equation, k in equation (7) thrust is a fixed value. In order to obtain the optimal k thrust , respectively take k thrust Let f be 50 numbers with equal intervals between 0 and 1 and solve equation (7) to find the k corresponding to the minimum f. thrust And x, y, z1, z2 are the optimal solutions; the thrust switch control variables △t1, △t2, α1, α2, u1, u2 are calculated by formula (3).

5. The analytical low-thrust circular orbit out-of-plane intersection optimization method according to claim 4 is characterized in that: The step 4 is specifically implemented according to the following steps: Step 4.1, Substitute △t1,△t2,α1,α2,u1,u2 into the small thrust precise dynamic equation corresponding to equation (9) to predict the orbital elements at the rendezvous time. The dynamic equation is: Among them, r and v are the spacecraft position velocity vectors, which can be converted to orbital roots; a p is the sum of the perturbation accelerations, e t ,e n are the tangential and normal unit vectors in the orbital coordinate system, respectively, and the calculation method is common sense; c is the thrust on and off flag that changes with time, 0 is off, and 1 is on; α is the thrust yaw angle that changes with time; the specific values ​​of c and α are determined by the parameters △t1, △t2, α1, α2, u1, and u2. The specific relationship is divided into two cases as follows; When time t≤t0+△t1: When t ≥ t f -△t2: The variables in the formula have been defined; Step 4.2: Use the numerical integration method to calculate the position velocity and orbital element [a] of the spacecraft at the rendezvous moment according to the initial orbital element of the spacecraft and equation (9). p ,e p ,i p ,Ω p ,ω p ,M p ], calculate the phase difference △u with the target orbital element p =ω f +M f -ω p -M p .

6. The analytical low-thrust circular orbit out-of-plane intersection optimization method according to claim 5 is characterized in that: The step 5 is specifically implemented according to the following steps: According to △u p Calculate the correction amount for △a1 where μ is the Earth's gravitational constant; then fix Resolve the remaining three variables y, z1, z2 in the optimization model (6); then, according to formula (3), the values ​​of the parameters △t1, △t2, α1, α2, u1, u2 can be recalculated.

7. The analytical low-thrust circular orbit out-of-plane intersection optimization method according to claim 6 is characterized in that: The step 6 is specifically implemented according to the following steps: Substitute △t1,△t2,α1,α2,u1,u2 obtained in step 6 back into step 5, re-predict the orbital elements of the spacecraft and calculate △u p ; Determine △u p Is it less than a given threshold? If so, output △t1,△t2,α1,α2,u1,u2 and the corresponding predicted orbital elements, which are the position, velocity and orbital elements of the spacecraft during flight, that is, the optimal rendezvous trajectory with low thrust; otherwise, return to step 6 to continue to correct △a1 and continue to calculate the new △u p , until △u p is less than a given threshold.

Images