Hybrid method of neural network and parameter optimization for low-thrust rendezvous trajectory optimization

Abstract

The present invention discloses a hybrid method of neural network and parameter optimization for low-thrust rendezvous trajectory optimization. First, a fast optimization method for low-thrust rendezvous control law of near-circular orbit with a flight time of 1 circle is obtained; then two neural network proxy models are established, corresponding to the given [Δa, Δi, ΔΩ, Δe x ,Δe y ] Output whether there is a feasible solution and the ratio of the arc length of the output thrust start to the entire cycle, where Δa is the change in semi-major axis, Δi is the change in inclination, ΔΩ is the change in right ascension of the ascending node, and Δe x is the change in eccentricity on the x-axis of the inertial system, Δe y The present invention solves the problem of excessive variables and low efficiency in existing technologies.

Description

Technical Field

[0001] The present invention belongs to the field of aerospace navigation control technology, and in particular relates to a hybrid method of neural network and parameter optimization for low-thrust rendezvous trajectory optimization. Background Art

[0002] Electric propulsion spacecraft has the characteristic of high specific impulse, and its fuel consumption is less than that of chemical propulsion under the same conditions. It is suitable for orbital rendezvous missions with long-term flights, such as on-orbit services. However, in low-Earth orbit, the acceleration of electric propulsion spacecraft is very small compared to the Earth's gravity, and many circles of flight are required to achieve a large range of orbital maneuvers, which also brings difficulties in optimizing the rendezvous trajectory. On the one hand, the orbit calculation time required to solve multi-circle flight trajectories increases, reducing efficiency; on the other hand, when the flight time is very long, the amount of maneuvering required to change the phase is very small, and there are local optimal solutions where the number of flight circles differs by an integer, which interferes with global optimization. The present invention proposes a hybrid optimization method, which approximates the complete rendezvous trajectory as a combination of multiple single-circle orbital maneuvers, and applies neural networks to improve the calculation of the single-circle optimal trajectory. Summary of the Invention

[0003] The purpose of the present invention is to provide a hybrid method of neural network and parameter optimization for low-thrust rendezvous trajectory optimization, which solves the problems of too many solution variables and low efficiency in the prior art.

[0004] The technical solution adopted by the present invention is a hybrid method of neural network and parameter optimization for low-thrust rendezvous trajectory optimization, which is specifically implemented according to the following steps:

[0005] Step 1: Obtain a fast optimization method for the low-thrust rendezvous control law for a near-circular orbit with a flight time of one circle;

[0006] Step 2: Establish two neural network agent models, corresponding to the given [Δa, Δi, ΔΩ, Δe x ,Δe y ] Output whether there is a feasible solution and the ratio of the arc length of the output thrust start to the entire cycle, where Δa is the change in semi-major axis, Δi is the change in inclination, ΔΩ is the change in right ascension of the ascending node, and Δe x is the change in eccentricity on the x-axis of the inertial system, Δe y is the change of eccentricity on the y-axis of the inertial system;

[0007] Step 3: Establish a parameterized optimization model for multi-circuit orbit rendezvous;

[0008] Step 4: Based on the given input, use the differential evolution algorithm to obtain the optimal solution and then obtain the optimal trajectory.

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

[0010] Step 1 is implemented as follows:

[0011] Given the spacecraft, the orbital element changes Δa, Δi, ΔΩ, and Δe that need to be achieved by a small thrust are: x ,Δe y , Δa is the change of semi-major axis, Δi is the change of inclination, ΔΩ is the change of right ascension of ascending node, Δe x is the change in eccentricity on the x-axis of the inertial system, Δe y is the change in eccentricity along the y-axis of the inertial system. The transfer time is fixed as one revolution, i.e., one orbital period. The time required for the spacecraft's latitude argument to change from 0 to 2π is used to solve the optimal small thrust control law.

[0012] According to the dynamic equation, for small eccentricity orbits, ignoring the coupling effect of eccentricity, we have:

[0013]

[0014] in are the three components of the unit vector of the thrust acceleration direction in the orbital coordinate system, ft, fn, fr are the tangential, normal and radial components respectively, c is a coefficient between 0 and 1, which represents the ratio of the thrust amplitude to the maximum thrust acceleration of the spacecraft itself, which changes with time and is recorded as c(t), α is the maximum amplitude of the thrust acceleration, [a, i, Ω, e x , e y ] are the spacecraft's instantaneous semi-major axis, inclination, right ascension of the ascending node, projection of the eccentricity vector on the x-axis of the inertial system, and projection of the eccentricity vector on the y-axis of the inertial system, respectively; t is time, V0 is the spacecraft's initial velocity, a0, i0 are the spacecraft's initial semi-major axis and inclination, respectively;

[0015] Make the following assumptions:

[0016] (1) The change in the number of orbital elements caused by a small thrust within one cycle is very small, and the cycle is also approximately unchanged. Rewriting dt in equation (1) as du = n0dt yields:

[0017]

[0018] Where n0 is the initial orbital angular velocity;

[0019] (2) Since the orbital angular velocity remains unchanged, the optimal fuel index function for a single lap is written as:

[0020]

[0021] The constraint equation for the change in the number of orbital elements is written as:

[0022]

[0023] In summary, according to the minimum principle, let:

[0024]

[0025] Among them, the definition is the derivative of the orbital element, λ σ =[λ a ,λ i ,λ Ω ,λ ex ,λ ey ] is the undetermined Lagrange multiplier, also known as the co-state variable, then according to the minimum principle, the λ corresponding to the optimal fuel trajectory is σ And c should make L take the minimum value, and the terms related to f in formula (5) are combined into:

[0026]

[0027] Then, remember In order to minimize L, χ and f should be in opposite directions, that is:

[0028]

[0029] Substituting into formula (5), we get In order to minimize L, we have:

[0030]

[0031] At the same time, according to the minimum principle, the covariate variable λ σ Need to meet:

[0032]

[0033] Where σ=[a,i,Ω,e x ,e y ] is the number of orbital elements, indicating that λ σ is a constant vector. When the number of orbital elements is changed in a single cycle, the optimization model of the optimal small thrust is: the undetermined variables are the 5-dimensional co-state variables plus the coefficient λ0, which is 6-dimensional in total. The shooting equation is also 6-dimensional, that is, based on formula (4), a constraint formula is added, such as formula (10):

[0034]

[0035] The solution method of this model uses random guessing of the unknown variables as the initial value and uses nonlinear solution tools to solve it. In this way, the optimal λ is obtained. σ And the corresponding f and c, we get the optimal J; According to the definition of J (3), the ratio of the arc length of the thrust start to the entire cycle is

[0036] Step 2 is implemented as follows:

[0037] The first neural network structure is as follows: the input layer is Δa, Δi, ΔΩ, Δe x ,Δe y , the output layer is a 1-dimensional Boolean variable, 0 means no solution, 1 means there is a solution, the middle hidden layer is set to 2 layers, each layer has 60 nodes;

[0038] The second neural network structure is designed to be the same as the first one except for the output layer. The output layer of the second neural network structure is a 1-dimensional real variable, which represents the ratio of the arc length of the start-up to the entire cycle when there is a solution.

[0039] The training process of the neural network in step 2 is as follows,

[0040] First, according to formula (2), the maximum semi-major axis change Δa that can be achieved when all the thrust is tangential thrust is obtained: max for:

[0041]

[0042] Similarly, when all thrust is normal thrust, the maximum inclination angle and right ascension of the ascending node that can be achieved are Δi max and ΔΩ max They are:

[0043]

[0044] When all the thrust is radial thrust, the two components of the maximum eccentricity change Δe that can be achieved xmax and Δe ymax They are:

[0045]

[0046] Randomly generate 100,000 sets of 5-dimensional random variables between -1 and 1 {x i}, i = 1, 2, 3, 4, 5, corresponding to Δa, Δi, ΔΩ, Δe x ,Δe y The sampling value is x1Δa max ,x2Δi max ,x3ΔΩ max ,x4Δe xmax ,x5Δe ymax , substitute into model (10) to find out whether there is an optimal solution and the ratio of the corresponding startup arc length to the entire cycle;

[0047] The public neural network training framework Keras is used to train the neural network using the generated sample data to obtain the given Δa, Δi, ΔΩ, Δe expressed in the form of a neural network. x ,Δe y , quickly determine whether there is a solution and calculate the corresponding Proxy model.

[0048] Step 3 is implemented as follows:

[0049] The small thrust multi-turn orbit rendezvous problem is formulated as follows: given the initial orbit elements [a0, i0, Ω0, e x0 ,e y0 ,u0], where u0 is the latitude argument, given the target orbital element number [a f ,i f ,Ω f ,e xf ,e yf ,u f ], given the flight time Δt, solve c and f that satisfy the intersection constraint and minimize equation (15):

[0050]

[0051] Assume that the optimal thrust form is to divide Δt into three parts: Δt1, Δt-Δt1-Δt2, Δt2, where Δt1 and Δt2 correspond to the time when there is thrust at both ends of the rendezvous, Δt-Δt1-Δt2 is the time when the thrust is fixed at 0, and within the Δt1 time, the change in the number of orbital elements caused by the small thrust is [Δa1, Δi1, ΔΩ1, Δe x1 ,Δe y1 ], Δa1 is the change of the semi-major axis within Δt1, Δi1 is the change of the inclination within Δt1, ΔΩ1 is the change of the right ascension of the ascending node within Δt1, Δe x1 is the change in eccentricity on the x-axis of the inertial system during Δt1, Δe y1 is the change of eccentricity on the y-axis of the inertial system during Δt1. During Δt2, the change of orbital elements caused by the small thrust is [Δa2, Δi2, ΔΩ2, Δe x2 ,Δe y2 ], Δa2 is the change of the semi-major axis within Δt2, Δi2 is the change of the inclination within Δt2, ΔΩ2 is the change of the right ascension of the ascending node within Δt2, Δe x2 is the change in eccentricity on the x-axis of the inertial system during Δt2, Δe y2 is the change of eccentricity on the y-axis of the inertial system during the time Δt2, then [Δa1, Δi1, ΔΩ1, Δe x1 ,Δe y1] and [Δa2,Δi2,ΔΩ2,Δe x2 ,Δe y2 ] satisfies the following constraints:

[0052]

[0053] in, ΔΩ d is the change in right ascension of the ascending node caused by the J2 perturbation, calculated according to the following formula:

[0054]

[0055] in, is the right ascension drift rate of the ascending node of the initial orbit, is the right ascension drift rate of the ascending node after adding Δa1 and Δi1 to the initial orbit, is the right ascension drift rate of the ascending node of the target orbit. The calculation formula of the right ascension drift rate of the ascending node is:

[0056]

[0057] Where J2=0.0010826, r E =6378137, a and i are the semi-major axis and inclination, calculate as well as When using the initial orbit, the initial orbit plus Δa1 and Δi1, and the target orbit, respectively;

[0058] Δω d is the change in the argument of perigee caused by the J2 perturbation, calculated according to the following formula:

[0059]

[0060] in, is the perigee argument drift rate of the initial orbit, is the perigee argument drift rate after adding Δa1 and Δi1 to the initial orbit, is the perigee argument drift rate of the target orbit. The calculation formula for the perigee argument drift rate is:

[0061]

[0062] The definition and calculation of the parameters substituted in are the same as ΔΩ in formula (17) d same;

[0063] Δu d is the phase angle change, calculated as follows:

[0064]

[0065] in, is the phase drift rate of the initial orbit, is the phase drift rate after adding Δa1 and Δi1 to the initial orbit, is the phase drift rate of the target orbit, and the calculation formula of the phase drift rate is:

[0066]

[0067] The definition and calculation of the parameters substituted in are the same as ΔΩ in formula (17) d same;

[0068] In summary, the optimal thrust corresponding to the time Δt1 is As well as Δt1 and Δt2 are optimization variables, given any values, the actual thrust start-up time is calculated;

[0069] The process is, first Substituting into equations (2) and (4), we can calculate the change in the number of orbital elements generated during the time Δt1 [Δa1, Δi1, ΔΩ1, Δe x1 ,Δe y1 ] and the corresponding thrust startup time

[0070] Then [Δa1,Δi1,ΔΩ1,Δe x1 ,Δe y1 ] into formula (15), calculate the change in the number of orbital elements corresponding to the time Δt2 [Δa2, Δi2, ΔΩ2, Δe x2 ,Δe y2 ];

[0071] [Δa2,Δi2,ΔΩ2,Δe x2 ,Δe y2 ] is converted into the change in the number of orbital elements of a single circle, that is,

[0072] Substitute the first neural network proxy model obtained in step 2 to determine whether there is a solution. If there is no solution, add a large constant penalty term to the actual thrust startup time. If there is a solution, continue to use the second neural network to calculate the corresponding The total thrust start time is That is, the optimal fuel index function;

[0073] Finally, we get the optimization variable [λ a λ i λ Ω λ ex λ ey λ0] and Δt1, Δt2 to calculate the objective function Optimization model.

[0074] Step 4 is implemented as follows:

[0075] Get the optimization variable [λ a λ i λ Ω λ ex λ ey After obtaining the optimal solution of λ0], Δt1, and Δt2, the thrust control law and trajectory during the entire rendezvous process are calculated according to Equations (2), (4), and (10) to complete the trajectory optimization.

[0076] The beneficial effect of the present invention is that the hybrid method of neural network and parameter optimization for low-thrust rendezvous trajectory optimization can quickly obtain the low-thrust optimal rendezvous trajectory considering J2 perturbation at a relatively low computational cost. BRIEF DESCRIPTION OF THE DRAWINGS

[0077] Figure 1 It is a flow chart of the hybrid method of neural network and parameter optimization for low-thrust rendezvous trajectory optimization of the present invention. DETAILED DESCRIPTION

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

[0079] The present invention is used for the hybrid method of neural network and parameter optimization for low-thrust rendezvous trajectory optimization, the flow chart is as follows Figure 1 As shown, please follow the steps below:

[0080] Step 1: Obtain a fast optimization method for the low-thrust rendezvous control law for a near-circular orbit with a flight time of one circle;

[0081] Step 1 is implemented as follows:

[0082] Given the spacecraft, the orbital element changes Δa, Δi, ΔΩ, and Δe that need to be achieved by a small thrust are: x ,Δe y , Δa is the change of semi-major axis, Δi is the change of inclination, ΔΩ is the change of right ascension of ascending node, Δe x is the change in eccentricity on the x-axis of the inertial system, Δe y is the change in eccentricity along the y-axis of the inertial system. The transfer time is fixed as one revolution, i.e., one orbital period. The time required for the spacecraft's latitude argument to change from 0 to 2π is used to solve the optimal small thrust control law.

[0083] According to the dynamic equation, for small eccentricity orbits, ignoring the coupling effect of eccentricity, we have:

[0084]

[0085] in are the three components of the unit vector of the thrust acceleration direction in the orbital coordinate system, ft, fn, fr are the tangential, normal and radial components respectively, c is a coefficient between 0 and 1, which represents the ratio of the thrust amplitude to the maximum thrust acceleration of the spacecraft itself, which changes with time and is recorded as c(t), α is the maximum amplitude of the thrust acceleration, [a, i, Ω, e x , e y ] are the spacecraft's instantaneous semi-major axis, inclination, right ascension of the ascending node, projection of the eccentricity vector on the x-axis of the inertial system, and projection of the eccentricity vector on the y-axis of the inertial system, respectively; t is time, V0 is the spacecraft's initial velocity, a0, i0 are the spacecraft's initial semi-major axis and inclination, respectively;

[0086] Make the following assumptions:

[0087] (1) The change in the number of orbital elements caused by a small thrust within one cycle is very small, and the cycle is also approximately unchanged. Rewriting dt in equation (1) as du = n0dt yields:

[0088]

[0089] Where n0 is the initial orbital angular velocity;

[0090] (2) Since the orbital angular velocity remains unchanged, the optimal fuel index function for a single lap is written as:

[0091]

[0092] The constraint equation for the change in the number of orbital elements is written as:

[0093]

[0094] In summary, according to the minimum principle, let:

[0095]

[0096] Among them, the definition is the derivative of the orbital element, is the undetermined Lagrange multiplier, also known as the co-state variable, then according to the minimum principle, the λ corresponding to the optimal fuel trajectory is σ And c should make L take the minimum value, and the terms related to f in formula (5) are combined into:

[0097] Then, remember In order to minimize L, χ and f should be in opposite directions, that is:

[0098]

[0099] Substituting into formula (5), we get In order to minimize L, we have:

[0100]

[0101] At the same time, according to the minimum principle, the covariate variable λ σ Need to meet:

[0102]

[0103] Where σ=[a,i,Ω,e x ,e y ] is the number of orbital elements, indicating that λ σ is a constant vector. When the number of orbital elements is changed in a single cycle, the optimization model of the optimal small thrust is: the undetermined variables are the 5-dimensional co-state variables plus the coefficient λ0, which is 6-dimensional in total. The shooting equation is also 6-dimensional, that is, based on formula (4), a constraint formula is added, such as formula (10):

[0104]

[0105] The solution method of this model uses random guessing of the unknown variables as the initial value and uses nonlinear solution tools to solve it. In this way, the optimal λ is obtained. σ And the corresponding f and c, we get the optimal J; According to the definition of J (3), the ratio of the arc length of the thrust start to the entire cycle is It should be noted that if the given orbital element number change is too large and the small thrust amplitude is too small, Equation (10) may have no solution, corresponding to the situation where the small thrust and full thrust cannot be achieved.

[0106] Step 2: Establish two neural network agent models, corresponding to the given [Δa, Δi, ΔΩ, Δe x ,Δe y ] Output whether there is a feasible solution and the ratio of the arc length of the output thrust start to the entire cycle, where Δa is the change in semi-major axis, Δi is the change in inclination, ΔΩ is the change in right ascension of the ascending node, and Δe x is the change in eccentricity on the x-axis of the inertial system, Δe y is the change of eccentricity on the y-axis of the inertial system;

[0107] Step 2 is implemented as follows:

[0108] The first neural network structure is as follows: the input layer is Δa, Δi, ΔΩ, Δe x ,Δe y , the output layer is a 1-dimensional Boolean variable, 0 means no solution, 1 means there is a solution, the middle hidden layer is set to 2 layers, each layer has 60 nodes;

[0109] The second neural network structure is designed to be the same as the first one except for the output layer. The output layer of the second neural network structure is a 1-dimensional real variable, which represents the ratio of the arc length of the start-up to the entire cycle when there is a solution.

[0110] The training process of the neural network in step 2 is as follows,

[0111] First, according to formula (2), the maximum semi-major axis change Δa that can be achieved when all the thrust is tangential thrust is obtained: max for:

[0112]

[0113] Similarly, when all thrust is normal thrust, the maximum inclination angle and right ascension of the ascending node that can be achieved are Δi max and ΔΩ max They are:

[0114]

[0115] When all the thrust is radial thrust, the two components of the maximum eccentricity change Δe that can be achieved xmax and Δe ymax They are:

[0116]

[0117] Randomly generate 100,000 sets of 5-dimensional random variables between -1 and 1 {x i}, i = 1, 2, 3, 4, 5, corresponding to Δa, Δi, ΔΩ, Δe x ,Δe y The sampling value is x1Δa max ,x2Δi max ,x3ΔΩ max ,x4Δe xmax ,x5Δe ymax , substitute into model (10) to find out whether there is an optimal solution and the ratio of the corresponding startup arc length to the entire cycle;

[0118] The public neural network training framework Keras is used to train the neural network using the generated sample data. Keras is common knowledge in the field and will not be described here. The given Δa, Δi, ΔΩ, Δe expressed in the form of a neural network are obtained. x ,Δe y , quickly determine whether there is a solution and calculate the corresponding Proxy model.

[0119] Step 3: Establish a parameterized optimization model for multi-circuit orbit rendezvous;

[0120] Step 3 is implemented as follows:

[0121] The small thrust multi-turn orbit rendezvous problem is formulated as follows: given the initial orbit elements [a0, i0, Ω0, e x0 ,e y0 ,u0], where u0 is the latitude argument, (a0,i0,Ω0,e x0 ,e y0 As defined above), given the target orbital elements [a f ,i f ,Ω f ,e xf ,e yf ,u f ], given the flight time Δt, solve c and f that satisfy the intersection constraint and minimize equation (15):

[0122]

[0123] Assume that the optimal thrust form is to divide Δt into three parts: Δt1, Δt-Δt1-Δt2, Δt2, where Δt1 and Δt2 correspond to the time when there is thrust at both ends of the rendezvous, Δt-Δt1-Δt2 is the time when the thrust is fixed at 0, and within the Δt1 time, the change in the number of orbital elements caused by the small thrust is [Δa1, Δi1, ΔΩ1, Δe x1 ,Δe y1 ], Δa1 is the change of the semi-major axis within Δt1, Δi1 is the change of the inclination within Δt1, ΔΩ1 is the change of the right ascension of the ascending node within Δt1, Δe x1 is the change in eccentricity on the x-axis of the inertial system during Δt1, Δe y1 is the change of eccentricity on the y-axis of the inertial system during Δt1. During Δt2, the change of orbital elements caused by the small thrust is [Δa2, Δi2, ΔΩ2, Δe x2 ,Δe y2 ], Δa2 is the change of the semi-major axis within Δt2, Δi2 is the change of the inclination within Δt2, ΔΩ2 is the change of the right ascension of the ascending node within Δt2, Δe x2 is the change in eccentricity on the x-axis of the inertial system during Δt2, Δe y2 is the change of eccentricity on the y-axis of the inertial system during the time Δt2, then [Δa1, Δi1, ΔΩ1, Δe x1 ,Δe y1 ] and [Δa2,Δi2,ΔΩ2,Δe x2 ,Δe y2 ]) satisfies the following constraints:

[0124]

[0125] in, ΔΩ d is the change in right ascension of the ascending node caused by the J2 perturbation, calculated according to the following formula:

[0126]

[0127] in, is the right ascension drift rate of the ascending node of the initial orbit, is the right ascension drift rate of the ascending node after adding Δa1 and Δi1 to the initial orbit, is the right ascension drift rate of the ascending node of the target orbit. The calculation formula of the right ascension drift rate of the ascending node is:

[0128]

[0129] Where J2=0.0010826, r E =6378137, a and i are the semi-major axis and inclination, calculate as well as When using the initial orbit, the initial orbit plus Δa1 and Δi1, and the target orbit, respectively;

[0130] Δω d is the change in the argument of perigee caused by the J2 perturbation, calculated according to the following formula:

[0131]

[0132] in, is the perigee argument drift rate of the initial orbit, is the perigee argument drift rate after adding Δa1 and Δi1 to the initial orbit, is the perigee argument drift rate of the target orbit. The calculation formula for the perigee argument drift rate is:

[0133]

[0134] The definition and calculation of the parameters substituted in are the same as ΔΩ in formula (17) d same;

[0135] Δu d is the phase angle change, calculated as follows:

[0136]

[0137] in, is the phase drift rate of the initial orbit, is the phase drift rate after adding Δa1 and Δi1 to the initial orbit, is the phase drift rate of the target orbit, and the calculation formula of the phase drift rate is:

[0138]

[0139] The definition and calculation of the parameters substituted in are the same as ΔΩ in formula (17) d same;

[0140] In summary, the optimal thrust corresponding to the time Δt1 is As well as Δt1 and Δt2 are optimization variables, given any values, the actual thrust start-up time is calculated;

[0141] The process is, first Substituting into equations (2) and (4), we can calculate the change in the number of orbital elements generated during the time Δt1 [Δa1, Δi1, ΔΩ1, Δe x1 ,Δe y1 ] and the corresponding thrust startup time Then [Δa1,Δi1,ΔΩ1,Δe x1 ,Δe y1 ] into formula (15), calculate the change in the number of orbital elements corresponding to the time Δt2 [Δa2, Δi2, ΔΩ2, Δe x2 ,Δe y2 ];

[0142] [Δa2,Δi2,ΔΩ2,Δe x2 ,Δe y2 ] is converted into the change in the number of orbital elements of a single circle, that is,

[0143] Substitute the first neural network proxy model obtained in step 2 to determine whether there is a solution. If there is no solution, add a large constant penalty term to the actual thrust startup time. If there is a solution, continue to use the second neural network to calculate the corresponding The total thrust start time is That is, the optimal fuel index function;

[0144] Finally, we get the optimization variable [λ a λ i λ Ω λ ex λ ey λ0] and Δt1, Δt2 to calculate the objective function Optimization model.

[0145] Step 4: Based on the given input, use the differential evolution algorithm to obtain the optimal solution and then obtain the optimal trajectory.

[0146] Step 4 is implemented as follows:

[0147] The differential evolution algorithm is common knowledge in the field and will not be described in detail in this invention. Its basic process is to initialize a set of random optimization variables as the initial population and calculate the objective function corresponding to each set of variables; then perform operator operations such as selection, crossover, and mutation according to the value of the objective function; then use the generated new individuals as the next generation population and recalculate the objective function value; then continuously repeat the selection, crossover, mutation and other processes until the individuals converge to the optimal solution.

[0148] Get the optimization variable [λ a λ i λ Ω λ ex λ ey After obtaining the optimal solution of λ0], Δt1, and Δt2, the thrust control law and trajectory during the entire rendezvous process are calculated according to Equations (2), (4), and (10) to complete the trajectory optimization.

[0149] Example

[0150] This invention provides a specific embodiment of a low-thrust rendezvous control law optimization method based on simplified parameter optimization, which is used to optimize the optimal control law for low-thrust rendezvous. The inputs are as follows: the initial and target orbital elements of the spacecraft are shown in Table 1, the rendezvous flight time is 10.0 days, the thrust is 0.35N, the initial mass is 100kg, and the acceleration is 3.5e-3m / s. 2 .

[0151] Table 1 Spacecraft orbital elements

[0152]

[0153]

[0154] The steps and results are as follows:

[0155] (1) First, a single-loop orbit optimization model is established, where a0 and i0 are taken as the root-mean-square semi-major axis and inclination of the initial orbit, respectively.

[0156] (2) Generate 100,000 sets of random vectors Δa, Δi, ΔΩ, Δe x ,Δe y , solve the equation corresponding to formula (10) to obtain whether there is an optimal solution, and the ratio of the arc length of the startup corresponding to the optimal solution to the entire cycle.

[0157] (3) Substitute the initial orbit, target orbit, acceleration, and rendezvous flight time into the objective function to obtain the optimization model.

[0158] (4) Set the dimension of the differential evolution algorithm to 7, the population size to 50, the evolutionary generations to 2000, and the crossover mutation probability to 0.8. The optimal solution is: Δt1 = 2.345 days, Δt2 = 0.836 days, That is, the spacecraft first flies in a specific thrust mode for 2.345 days, then switches to no thrust, and flies in another thrust mode for the last 0.836 days, finally achieving rendezvous with the target orbit. The actual thrust-on time in Δt1 and Δt2 totals 0.841 days, and the equivalent speed increment is 254.4 m / s.

[0159] The optimization calculation of the present invention takes only 30 seconds, and the obtained trajectory can be used as a reference trajectory for orbit rendezvous mission analysis.

Claims

1. A hybrid method of neural network and parameter optimization for low-thrust rendezvous trajectory optimization, characterized in that: Please follow the steps below to implement it: Step 1: Obtain a fast optimization model for the low-thrust rendezvous control law in a near-circular orbit with a flight time of one circle; Step 2: Establish two neural network agent models, corresponding to the given [ ] Output whether there is a feasible solution and the ratio of the arc length of the output thrust start to the entire cycle, where, is the change of the semi-major axis, is the tilt angle change, is the change in right ascension of the ascending node, is the change in eccentricity on the x-axis of the inertial system, is the change in eccentricity along the y-axis of the inertial system; Step 3: Establish a parameterized optimization model for multi-circuit orbit rendezvous; Step 4: Based on the given input, use the differential evolution algorithm to obtain the optimal solution and then obtain the optimal trajectory.

2. The hybrid method of neural network and parameter optimization for low-thrust rendezvous trajectory optimization according to claim 1, characterized in that: The step 1 is specifically implemented according to the following steps: The change in orbital elements that a given spacecraft needs to achieve with a small thrust , is the change of the semi-major axis, is the tilt angle change, is the change in right ascension of the ascending node, is the change in eccentricity on the x-axis of the inertial system, is the change in eccentricity along the y-axis of the inertial system. The transfer time is fixed as one revolution, i.e., one orbital period. The time required for the spacecraft's latitude argument to change from 0 to 2π is used to solve the optimal small thrust control law. According to the dynamic equation, for small eccentricity orbits, ignoring the coupling effect of eccentricity, we have: （1） in are the three components of the unit vector of the thrust acceleration direction in the orbital coordinate system, ft, fn, fr are the tangential, normal and radial components respectively, The coefficient between 0 and 1 represents the ratio of the thrust amplitude to the maximum thrust acceleration of the spacecraft itself, which changes with time and is recorded as , is the maximum magnitude of thrust acceleration, [ a,i,Ω,e x ，e y ] are the spacecraft's instantaneous semi-major axis, inclination, right ascension of the ascending node, and eccentricity vector in the inertial system. x Axis projection, eccentricity vector in inertial system y The projection of the axis, t For time, is the initial velocity of the spacecraft, , are the initial semi-major axis and inclination of the spacecraft, respectively; Make the following assumptions: (1) The change of the number of orbital elements caused by a small thrust in one cycle is very small, and the cycle is also approximately unchanged. t Rewrite as ,get: （2） in is the initial orbital angular velocity; (2) Since the orbital angular velocity remains unchanged, the optimal fuel index function for a single lap is written as: （3） The constraint equation for the change in the number of orbital elements is written as: （4） In summary, according to the minimum principle, let: （5） Among them, the definition is the derivative of the orbital element, is the undetermined Lagrange multiplier, also known as the co-state variable, then according to the minimum principle, the optimal fuel trajectory corresponds to as well as c should be Take the minimum value, and the formula (5) is Related items are combined into: （6） Then, remember , in order to make L Take the minimum value, and The direction should be opposite, that is: （7） Substituting into formula (5), we get , in order to L Taking the minimum value, we have: （8） At the same time, according to the minimum principle, the covariate variables Need to meet: （9） in, is the number of orbital elements, indicating is a constant vector, then when the number of orbital elements in a single cycle is changed, the optimal small thrust optimization model is: the undetermined variable is the 5-dimensional co-state variable plus the coefficient There are 6 dimensions in total, and the shooting equation is also 6-dimensional, that is, based on formula (4), a constraint formula is added, such as formula (10): （10） The solution method of this model uses random guessing of the unknown variables as the initial value and uses nonlinear solution tools to solve it. In this way, the optimal and the corresponding and c, we get the optimal ;according to The definition of formula (3) shows that the ratio of the arc length of thrust start to the entire cycle is .

3. The hybrid method of neural network and parameter optimization for low-thrust rendezvous trajectory optimization according to claim 2, characterized in that: The step 2 is specifically implemented according to the following steps: The first neural network structure is as follows: the input layer is , the output layer is a 1-dimensional Boolean variable, 0 means no solution, 1 means there is a solution, the middle hidden layer is set to 2 layers, each layer has 60 nodes; The second neural network structure is designed to be the same as the first one except for the output layer. The output layer of the second neural network structure is a 1-dimensional real variable, which represents the ratio of the arc length of the start-up to the entire cycle when there is a solution. .

4. The hybrid method of neural network and parameter optimization for low-thrust rendezvous trajectory optimization according to claim 3, characterized in that: The training process of the neural network in step 2 is as follows: First, according to formula (2), the maximum semi-major axis change that can be achieved when all the thrust is tangential thrust is obtained: for: （11） Similarly, when all thrust is normal thrust, the maximum inclination angle and right ascension of the ascending node that can be achieved are and They are: （12） When all thrust is radial thrust, the two components of the maximum eccentricity change that can be achieved and They are: （13） Randomly generate 100,000 sets of 5-dimensional random variables between -1 and 1 , corresponding to the generation The sampling value is , substitute into model (10) to find out whether there is an optimal solution and the ratio of the corresponding startup arc length to the entire cycle; The public neural network training framework Keras is used to train the neural network using the generated sampling data to obtain the given , quickly determine whether there is a solution and calculate the corresponding Proxy model.

5. The hybrid method of neural network and parameter optimization for low-thrust rendezvous trajectory optimization according to claim 4, characterized in that: The step 3 is specifically implemented according to the following steps: The low-thrust multi-turn orbit rendezvous problem is formulated as: given the number of initial orbit elements ,in, is the latitude argument, given the number of target orbit elements , given the flight time Δ t , solve the problem of satisfying the intersection constraint c and f, and minimize equation (15): （14） Assume that the optimal thrust form is to t Divided into three parts: , , ,in and The time when there is thrust at both ends of the intersection, The time when the thrust is fixed at 0, and During this time, the change in the number of orbital elements caused by the small thrust is , for The change of the semi-major axis in time, for The change in inclination angle over time, for The change in right ascension of the ascending node during this time, for The change in eccentricity on the x-axis of the inertial system during time, for The change in eccentricity on the y-axis of the inertial system during time is During this time, the change in the number of orbital elements caused by the small thrust is , for The change of the semi-major axis in time, for The change in inclination angle over time, for The change in right ascension of the ascending node during this time, for The change in eccentricity on the x-axis of the inertial system during time, for The change in eccentricity on the y-axis of the inertial system during time is and Satisfy the following constraints: （15） in, , , is the change in right ascension of the ascending node caused by the J2 perturbation, calculated according to the following formula: （16） in, is the right ascension drift rate of the ascending node of the initial orbit, Add the initial track and The right ascension drift rate of the ascending node after is the right ascension drift rate of the ascending node of the target orbit. The calculation formula of the right ascension drift rate of the ascending node is: （17） in =0.0010826, =6378137, a and i are the semi-major axis and inclination, calculate 、 as well as When using the initial orbit, the initial orbit plus and , and the value of the target orbit; is the change in the argument of perigee caused by the J2 perturbation, calculated according to the following formula: （18） in, is the perigee argument drift rate of the initial orbit, Add the initial track and The perigee argument drift rate after is the perigee argument drift rate of the target orbit. The calculation formula for the perigee argument drift rate is: （19） The definition and calculation of the parameters substituted in are the same as those in formula (17). same; is the phase angle change, calculated as follows: （20） in, is the phase drift rate of the initial orbit, Add the initial track and The phase drift rate after is the phase drift rate of the target orbit, and the calculation formula of the phase drift rate is: （21） The definition and calculation of the parameters substituted in are the same as those in formula (17). same; In summary, corresponding Optimal thrust within a certain time ,as well as and To optimize the variables, given any value, the actual thrust start-up time is calculated; The process is, first Substituting into equations (2) and (4), we can calculate The change in the number of orbital elements produced over time And the corresponding thrust startup time ; then Substitute into formula (15) and calculate Change in the number of orbital elements corresponding to time ;Will The change in the number of orbital elements converted into a single circle is , substitute the first neural network proxy model obtained in step 2 to determine whether there is a solution. If there is no solution, add a large constant penalty term to the actual thrust startup time. If there is a solution, continue to use the second neural network to calculate the corresponding ; The total thrust start time is , which is the optimal fuel index function; Finally, we get the optimized variables and 、 Calculate the objective function Optimization model.

6. The hybrid method of neural network and parameter optimization for low-thrust rendezvous trajectory optimization according to claim 5, characterized in that: The step 4 is specifically implemented according to the following steps: Get optimization variables and 、 After obtaining the optimal solution, the thrust control law and trajectory during the entire rendezvous process are calculated according to Equations (2), (4), and (10), completing the trajectory optimization.

Images