Indirect optimization method for small celestial body landing anti-interference trajectory

Abstract

The invention discloses an indirect optimization method for a landing anti-interference trajectory of a small celestial body, and relates to the field of deep space exploration. Aiming at the problem that a small celestial body landing trajectory is sensitive to uncertain parameters, the following solving method is provided: firstly, establishing a stochastic dynamics model of small celestial body landing, and constructing an anti-interference trajectory optimization problem; secondly, establishing a relation between an anti-interference trajectory optimization problem and a gravity-energy-free optimal control problem through a performance functional, a kinetic equation and a homotopy function of process noise, and estimating a co-state initial value of the gravity-energy-free optimal control problem; and finally, gradually converting the zero-gravity energy optimal control problem into an anti-interference trajectory optimization problem by changing homotopy coefficients, and in the conversion process, converting the optimization problem corresponding to each homotopy function into a two-point boundary value problem by using a Pontryagin minimum principle, and solving the two-point boundary value problem by using a targeting method. The sensitivity of the trajectory to uncertainty is reduced, so that the anti-interference capability of the landing trajectory is improved.

Description

Technical Field

[0001] The invention belongs to the field of deep space exploration technology, and in particular relates to an indirect optimization method for the anti-disturbance trajectory of a small celestial body landing Background Art

[0002] Small body exploration is of great significance to human understanding of the origin and evolution of the solar system, defense against near-Earth small body impacts, and verification of related deep space exploration technologies. Landing on the surface of a small body is an important means of exploration, and landing trajectory optimization is an important technical basis for the analysis, guidance, and control of landing exploration missions. Trajectory optimization methods mainly include direct methods based on parameter optimization and indirect methods based on the Pontryagin minimum principle. Trajectory optimization refers to optimizing an important indicator, such as flight time or fuel consumption, under the premise of satisfying terminal state constraints, path constraints, and control constraints. However, in actual deep space missions, the state of the spacecraft is difficult to accurately estimate, and the estimation of the dynamic parameters of small bodies may also have errors. If the influence of these uncertainties is not considered in trajectory design, the actual trajectory will deviate greatly from the nominal trajectory. Therefore, the design of anti-disturbance trajectory is crucial. At the same time, the direct method of trajectory optimization based on parameter optimization has a large amount of calculation, especially the Gaussian pseudo-Pu method, and the accuracy of its solution depends on the initial value guesses of state variables and control variables.

[0003] Based on the above analysis, the present invention aims to provide an indirect optimization method for the anti-disturbance landing trajectory of small celestial bodies. Aiming at the problem that the trajectory is sensitive to uncertainty, a performance functional weighted by covariance and fuel consumption is constructed, and the solution is obtained through the indirect method of trajectory optimization to reduce the sensitivity of the trajectory to uncertainty and improve the anti-disturbance capability. Summary of the invention

[0004] In view of this, the present invention aims to provide an indirect optimization method for the anti-disturbance landing trajectory of a small celestial body. Aiming at the problem that the trajectory is sensitive to uncertainty, a performance functional weighted by covariance and fuel consumption is constructed, and the solution is obtained through an indirect trajectory optimization method to reduce the sensitivity of the trajectory to uncertainty and improve the anti-disturbance performance.

[0005] In order to achieve the above-mentioned object of the invention, the present invention adopts the following technical scheme:

[0006] An indirect optimization method for the anti-disturbance trajectory of a small celestial body landing, characterized in that the method comprises the following steps:

[0007] S1: Establish a stochastic dynamics model for small body landing, construct a performance functional weighted by covariance and fuel consumption, and construct an anti-disturbance trajectory optimization problem

[0008] S2: The anti-disturbance trajectory optimization problem is linked to the zero-gravity energy optimal control problem by defining the homotopy coefficient performance functional, the dynamic equations, and the homotopy function of the process noise, and the initial value of the co-state of the zero-gravity energy optimal control problem is estimated.

[0009] S3: The zero-gravity energy optimal control problem is converted into an anti-disturbance trajectory optimization problem by changing the homology coefficient. For each trajectory optimization problem corresponding to the homology coefficient, the Pontryagin minimum principle is used to convert it into a two-point boundary value problem and then solved by the shooting method. In the process, the optimal control that minimizes the Hamiltonian function is solved using a one-dimensional search algorithm.

[0010] Furthermore, in step S1, the stochastic dynamics equation of the landing of the small celestial body is:

[0011]

[0012] In the above formula, X a =[r T v T mω z ] T is the state variable, r is the position of the detector, v is the velocity of the detector, m is the mass of the detector, and W is the process noise. The state covariance matrix of the above formula can be calculated along the nominal trajectory by the following formula:

[0013]

[0014] In the above formula, A is the Jacobian matrix of the dynamics, and Q is the power spectral density of the process noise W.

[0015] Furthermore, in step S1, the following covariance and fuel consumption weighted performance functional is constructed:

[0016]

[0017] In the above formula is the equation of the ith state variable at the terminal moment, C i are the weight coefficients of the weights of each state variable, c0 is the weight coefficient of the distribution covariance and fuel consumption, λ0 is a performance index coefficient greater than zero, and the subscript N represents the nominal quantity. The different elements of the covariance matrix are reorganized into vector form and expanded into the state of the state equation, together with the performance functional, equality constraints and inequality constraints to form the following anti-disturbance trajectory optimization problem:

[0018]

[0019] In the above formula, Vector(P) is the vector form of the upper right triangular elements of the covariance matrix.

[0020] Furthermore, in step S2, the following performance functionals, dynamic equations, and homotopy functions of process noise are defined respectively:

[0021]

[0022]

[0023]

[0024] When ε1=1 and ε0=0 in J1, J1 is the performance functional of the anti-disturbance trajectory optimization problem. When ε1=0 and ε2=0 in J1, J1 is the performance functional of the fuel optimal control problem. When ε1=0 and ε2=1 in J1, J1 is the performance functional of the energy optimal control problem.

[0025] When f h1 (ε3) and f h2 When ε3=0 in (ε3), use f h1 (ε3), f h2 (ε3) replaces the dynamic equation of velocity and process noise, and when ε1=0 and ε2=1 in J1, J1 is the performance functional of the optimal control problem of zero-gravity energy. h1 (ε3) and f h2 When ε3=1 in (ε3), ε1=0 and ε2=1 in J1, J1 is the performance functional of the energy optimal control problem. The connection between the anti-disturbance trajectory optimization problem and the zero-gravity energy optimal control problem is established through the homotopy function.

[0026] Furthermore, in step S2, the initial co-state value of the zero-gravity energy optimal control problem is estimated, and the initial co-state value of the position and velocity of the zero-gravity energy optimal control problem is estimated by solving the analytical solution of the initial co-state value of the zero-gravity control unconstrained energy optimal control problem:

[0027]

[0028]

[0029] The optimal solution of the optimal control problem of zero-gravity energy is used to approximate the initial value of the mass co-state of the optimal control problem of zero-gravity control and unconstrained energy:

[0030]

[0031] The initial value of the covariance of the optimal control problem of unconstrained energy without gravity control is:

[0032]

[0033] Furthermore, in step S3, the zero-gravity energy optimal control problem is converted into an anti-disturbance trajectory optimization problem by changing the homotopy coefficient. The optimization problem corresponding to each homotopy coefficient is converted into a two-point boundary value problem and solved by the shooting method. First, let f h1 (ε3) and f h2 ε3 in (ε3) changes from 0 to 1, ε1=0, ε2=1 in J1, this process converts the weightless energy optimal control problem into the energy optimal control problem. Secondly, let ε1=0, ε2 in J1 change from 1 to 0. This process converts the energy optimal control problem into the fuel optimal control problem. Finally, let ε2=0, ε1 in J1 change from 0 to 1. This process converts the fuel optimal control problem into the anti-disturbance trajectory optimization problem. In the above process, the optimization problem corresponding to each homotopy coefficient is converted into a two-point boundary value problem using the Pontryagin minimum principle. The guessed value of the costate initial value used to solve the target shooting equation uses the costate initial value obtained from the previous two-point boundary value problem.

[0034] Furthermore, step S3 uses a one-dimensional search algorithm to search for the optimal control that minimizes the Hamiltonian function, and uses a fixed-step fourth-order Runge-Kutta (RK4) algorithm to discretize the differential equation. For each time step k, the Hamiltonian function H at the current moment is k (u) Minimize the optimal control thrust ratio A one-dimensional search algorithm is used for searching, and the initial search interval is defined as [a, b], and the two split points are:

[0035]

[0036] In the above formula is the split ratio, and then the search interval is continuously updated according to the search interval update algorithm until the search interval is small enough. The optimal control thrust ratio is approximately:

[0037]

[0038] Compared with the prior art, the advantages and positive effects of the present invention are:

[0039] The present invention establishes a random dynamic model of small celestial body landing and constructs an anti-disturbance trajectory optimization problem; then, the anti-disturbance trajectory optimization problem is linked to the zero-gravity energy optimal control problem through performance functionals, dynamic equations, and homotopy functions of process noise, and the co-state initial value of the zero-gravity energy optimal control problem is estimated; finally, the zero-gravity energy optimal control problem is gradually converted into the anti-disturbance trajectory optimization problem by changing the homotopy coefficient. During the conversion process, the optimization problem corresponding to each homotopy function is converted into a two-point boundary value problem using the Pontryagin minimum principle and solved using the shooting method, thereby reducing the sensitivity of the trajectory to uncertainty and improving the anti-disturbance capability of the landing trajectory. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 A flowchart of an indirect optimization method for anti-disturbance track of small celestial body landing proposed by the present invention;

[0041] Figure 2 This is a flow chart of a homotopy algorithm in an indirect optimization method for anti-disturbance trajectory of small celestial body landing proposed by the present invention;

[0042] Figure 3 A one-dimensional search for solving the optimal control in an indirect optimization method for anti-disturbance trajectory of small celestial body landing proposed by the present invention;

[0043] Figure 4 The small celestial body landing simulation result in the embodiment of the present invention includes position, speed, mass and thrust;

[0044] Figure 5 is a position covariance variation curve of the fuel optimal trajectory and the anti-disturbance trajectory in an embodiment of the present invention;

[0045] Figure 6 It is the speed covariance variation curve of the optimal fuel trajectory and the anti-disturbance trajectory in the embodiment of the present invention. DETAILED DESCRIPTION

[0046] The present invention will be described in further detail below in conjunction with the accompanying drawings.

[0047] The present invention aims at the uncertainty and sensitivity of trajectory in deep space environment and proposes an indirect optimization method for anti-disturbance trajectory of small celestial body landing. The present invention will be described in detail below with reference to the accompanying drawings.

[0048] See also Figure 1 The present invention provides an indirect optimization method for a small celestial body landing anti-disturbance trajectory, the method comprising the following steps:

[0049] Step S1: Establish a stochastic dynamics model for the landing of a small celestial body and construct an anti-disturbance trajectory optimization problem:

[0050] Taking the center of mass of the small celestial body as the origin, the three coordinate axes of space x, y, and z are along the three principal axes of the small celestial body, and the fixed coordinate system of the small celestial body is established. Then the landing dynamics equation of the probe is:

[0051]

[0052] For ease of expression, rewrite the above formula as follows:

[0053]

[0054] In the above formula, r = [xyz] T is the position of the detector in the coordinate system, v = [v x v y v z ] T is the velocity of the detector, m is the mass of the detector, ω=[00ω z ] T is the spin angular velocity of the small celestial body, T max is the maximum thrust of the probe engine, u∈[0,1] is the thrust ratio, α=[α x α y α z ] T is the direction of thrust, and a=1, · represents the Euclidean norm, I sp is the specific impulse of the engine, g0 is the gravity acceleration at sea level, V = [V x V y V z ] T is the first-order partial derivative of the gravitational potential function of the small celestial body with respect to position. The uncertainty of the gravitational field is expressed as:

[0055]

[0056] In the above formula, C Vx , C Vy , C Vz is a random process that follows a standard normal distribution, σ Vx , σ Vy , σ Vz are the standard deviations of the gravitational field modeling errors in the three-axis directions. The thrust errors in the three-axis directions are expressed as:

[0057]

[0058] In the above formula, T=[T x T y T z ] T is the thrust in the three-axis direction, and T = T max uα,δ T =C Tσ T is the error of thrust size, C T is a random process that follows a standard normal distribution, σ T is the standard deviation of the thrust error. The uncertainty of the gravitational field and the uncertainty of the thrust are modeled as process noise in the following form:

[0059]

[0060] The uncertainty of the angular velocity of a small celestial body is: ω z =(1+C ω σ ωz )ω z,N , C ω In order to obey the random process of standard normal distribution, the spin angular velocity of small celestial bodies is expanded into the state of the state equation, and process noise is added to the velocity state equation to construct the random dynamics equation:

[0061]

[0062] In the above formula, X a =[r T v T mω z ] T ,and W=[000W x W y W z 00] T The state covariance matrix of the above formula can be calculated along the nominal trajectory as follows:

[0063]

[0064] In the above formula is the dynamic Jacobian matrix, Q is the power spectral density of the process noise W, and has the following relationship with the expectation of the process noise:

[0065]

[0066] In the above formula, δ is the Dirac function. Since the covariance matrix is ​​a symmetric matrix, the n-order covariance matrix has (n+1) / 2 different elements. The elements of the covariance matrix can be expanded to the state of the state equation, and the following covariance expansion dynamics equation is formed:

[0067]

[0068] In the above formula Vector (P) is the vector form of the upper right triangular element of the covariance matrix. Then, the performance functional weighted by covariance and fuel consumption is constructed, and the desensitized trajectory optimization problem is formed with the state equation, equality constraints and inequality constraints. It can be ignored when solving the optimization problem, so the anti-disturbance trajectory optimization problem is as follows:

[0069]

[0070]

[0071] In the above formula is the variance of the terminal time of the ith state variable, c0 is the weight coefficient of the distribution covariance and fuel consumption, C i is the weight coefficient for assigning weights to each state variable, for:

[0072]

[0073] Step S2: The anti-disturbance trajectory optimization problem is linked to the zero-gravity energy optimal control problem through the performance functional, the dynamic equation and the homotopy function of the process noise. And the approximate analytical solution of the initial value of the co-state of the zero-gravity energy optimal control problem is obtained.

[0074] First, define the homotopy function of the following performance functional:

[0075]

[0076] In the above formula, ε1, ε2∈[0,1] are homotopy coefficients, λ0>0 is the performance functional coefficient, and when ε1=1, ε2=0, J1 is the performance functional of the anti-disturbance trajectory optimization problem. When ε1=0, ε2=0, J1 is the performance functional of the fuel optimal control problem. When ε1=0, ε2=1, J1 is the performance functional of the energy optimal control problem. Then, define the homotopy function of the following dynamic equation:

[0077]

[0078] In the above formula, ε3∈[0,1] is the homotopy coefficient. The velocity dynamics equation in is:

[0079]

[0080] In the above formula, G(r,v,ω) is generalized gravity, expressed as:

[0081] G(r,v,ω)=-2ω×v-ω×(ω×r)+V

[0082] At the same time, the homotopy function of the process noise is defined as follows:

[0083]

[0084] Define the new process noise as: W′=f h2 (ε3). The rate spectrum density of the new process noise is recorded as Q'. The new covariance state equation is recorded as dVector(P') / dt, and Replace the covariance state equation in , and let And ε3=0 in W′, ε1=0, ε2=1 in J1, the generalized gravity in the dynamic equation disappears, and J1 is a gravity-free energy optimal control problem.

[0085] Secondly, in order to estimate the initial value of the co-state of the optimal control problem of zero-gravity energy, the following simplified dynamic equation is defined:

[0086]

[0087] In the above formula, a is the control acceleration, and the process noise is defined as follows:

[0088]

[0089] In the above formula, δ a is the error of acceleration (only contains the error of acceleration magnitude). The process noise is applied to the velocity state equation of simplified dynamics to form the following random dynamics equation:

[0090]

[0091] Derive the state covariance matrix equation of the above formula, expand it to the state, and ignore the terms with zero derivatives in the simplified dynamic equation The following non-gravity control unconstrained energy optimal control problem is constructed:

[0092]

[0093] In the above formula, Vector(P1) is the vector form of the upper right triangular element of the random dynamics covariance matrix. Using the Pontryagin minimum principle, the two-point boundary value problem (TPBVP) corresponding to the above formula, the Hamiltonian function is:

[0094]

[0095] The optimal control acceleration minimized by equation H3 can be calculated as follows:

[0096]

[0097] The co-state equation is:

[0098]

[0099]

[0100]

[0101] The transversality condition is:

[0102]

[0103] Substituting the transversal condition into the covariance covariance equation, we can obtain:

[0104]

[0105] From the above formula, we can get that the covariance is always equal to zero, and substituting the above formula into We can get:

[0106]

[0107] Combined with the above formula Integrate from the initial time t0 to the terminal time t f The co-state of velocity and position can be obtained:

[0108]

[0109]

[0110] Let the homotopy function J1 of the performance functional be ε1=0, ε2=1, ε3=0 in the homotopy function of dynamics and process noise, and 1 / m0≈1 / m, so that:

[0111]

[0112] By assigning the performance functional coefficients, the zero-gravity energy optimal control problem has the same performance functional as the zero-gravity control unconstrained energy optimal control problem. Since the zero-gravity control unconstrained energy optimal control problem does not provide a mass co-state, its value can be approximately calculated from the optimal solution of the zero-gravity energy optimal control problem:

[0113]

[0114] The initial value of the mass co-state for the optimal control problem of non-gravity control and unconstrained energy is:

[0115]

[0116] Step S3: The zero-gravity energy optimal control problem is converted into an anti-disturbance trajectory optimization problem by adjusting the homotopy coefficients in the performance functional. The trajectory optimization problem corresponding to each homotopy coefficient is converted into a two-point boundary value problem and solved using the shooting method. The optimal control that minimizes the Hamiltonian function is solved using a one-dimensional search algorithm.

[0117] First, let ε1=0, ε2=1 in J1, and let ε3 in the homotopy function of the dynamic equation and process noise change from 0 to 1, and the performance functional J1 becomes:

[0118]

[0119] The Hamiltonian function at this time is:

[0120]

[0121] In the above formula is the covariate variable, and the covariate equation is:

[0122]

[0123] In the above formula α * are the optimal control thrust ratio and optimal control thrust direction that minimize the Hamiltonian function, By using a one-dimensional search algorithm to search for u that minimizes the Hamiltonian function at the current moment N,k To obtain, the specific process is as follows:

[0124] First, the fourth-order Runge-Kutta (RK4) algorithm is used to discretize the differential equations of the two-point boundary value problem. For each time step k, the Hamiltonian function H at the current moment is k (u) Minimize the optimal control thrust ratio A one-dimensional search algorithm is used for searching, and the initial search interval is defined as [a, b], and the two split points are:

[0125]

[0126] In the above formula is the segmentation ratio; then, the interval is reduced according to the search interval update algorithm until the search interval is small enough; finally, the optimal control thrust ratio can be approximated as:

[0127]

[0128] The optimal control thrust direction can be obtained by the following method:

[0129] The Hamiltonian function is reformulated as follows:

[0130]

[0131] In the above formula, h(X P ,λ,u N ) is a function including state variables, co-state variables and control variables, p = [p (1) p (2) p (3) ] T is the primer vector, and its expression is as follows:

[0132]

[0133]

[0134]

[0135] Then the optimal control thrust direction is parallel to p and opposite to it, that is, The state equation and the co-state equation constitute a differential equation of a two-point boundary value problem, and the boundary conditions are:

[0136]

[0137] Since the quality and state covariance at the terminal moment are unknown, the following transversal conditions are also required:

[0138]

[0139] The shooting equation used to solve the two-point boundary value problem is:

[0140] Φ(λ r0 ,λ v0 ,λ m0 ,λ Vector(P′)0 )=[r(t f )-r f ,v(t f ),λ m (t f ),λ Vector(P′) (t f )] T

[0141] In the above formula, λ r0 , v0 , m0 and λ Vector(P′)0 is the guessed value of the initial value of the co-state. The process of converting the zero-gravity energy optimal control problem into the energy optimal control problem is as follows:

[0142] First, the initial value of the costate of the optimal control problem of unconstrained energy under gravity-free control is taken as the guess value of the initial value of the costate when ε3=0. Then, ε3 is increased by a fixed step size dε3 in each step. For each two-point boundary value problem corresponding to ε3, the initial value of the costate obtained by the previous ε3 is used as the guess value of the initial value of the costate to solve the targeting equation. If the solution fails, dε3 is reduced to 0.8dε3 to continue solving until ε3=1. The above conversion process is recorded as Problem 1.

[0143] Secondly, ε2 in J1 = 0, ε1 changes from 1 to 0, and J1 is:

[0144]

[0145] The Hamiltonian function at this time is:

[0146]

[0147] The co-state equation is:

[0148]

[0149] The optimal control thrust ratio that minimizes the Hamiltonian function in the above formula is The same method as in Problem 1 is used to obtain the optimal control thrust direction α * The same as problem 1, boundary conditions, transversal conditions and target equations are also the same as process 1. The process of converting the energy optimal control problem into the fuel optimal control problem is as follows: First, the initial value of the co-state obtained in problem 1 (ε3 = 1) is used as the guess value of the initial value of the co-state of the two-point boundary value problem when ε2 = 1 to solve the target equation, and then ε2 is reduced by a fixed step size dε2 in each step. For each ε2, the initial value of the co-state obtained by the previous ε2 is used as the guess value of the co-state of the target equation to solve. If the solution fails, dε2 is reduced to 0.8dε2 to continue solving until ε2 = 0. The above conversion process is recorded as problem 2.

[0150] Finally, let ε2 in J1 = 0, and ε1 changes from 0 to 1. At this time, J1 is:

[0151]

[0152] The Hamiltonian function at this time is:

[0153]

[0154] The co-state equation is the same as that in Problem 2. The optimal control thrust ratio that minimizes the Hamiltonian function is The same method as in Problem 1 is used to obtain the optimal control thrust direction α *Same as Problem 1, the boundary conditions are the same as Process 1, and the transversal conditions are:

[0155]

[0156] The shooting equation is:

[0157]

[0158] The conversion process of the fuel optimal control problem into the disturbance rejection trajectory optimization problem is as follows:

[0159] First, the initial value of the costate obtained in process 2 (ε2 = 0) is used as the guess value of the initial value of the costate when ε1 = 0 to solve the target equation. Then, ε1 is increased by a fixed step size dε1 in each step. For each ε1, the initial value of the costate obtained by the previous ε1 is used as the guess value of the costate of the target equation to solve. If the solution fails, dε1 is reduced to 0.8dε1 to continue solving until ε1 = 1. The above conversion process is recorded as Problem 1.

[0160] In this simulation, the landing asteroid 433Eros was selected as the background and the model was built in the MATLAB environment. The nominal spin angular velocity of the asteroid is 3.314×10 -4 rad / s. The maximum thrust of the detector's engine is 20N, and the specific impulse is 20s. The initial position of the detector is r0 = [1011769568256] T , the initial velocity is v0 = [-25-12-17] T , the initial probe mass is 300kg, and the landing position is r f =[853501045] T , the landing speed is v f =

[000] T The thrust uncertainty parameter σ of the probe is T Satisfies Gaussian distribution N(0,(0.5%T max u N ) 2 ), the uncertainty parameter of the spin angular velocity σ ωz Obey Gaussian distribution N(0,1%ω z,N ), the uncertainty parameter σ of the gravitational field Vx , σ Vx , σ Vz Obey Gaussian distribution N(0,(15%) 2 ), the initial position error follows the uniform distribution U(-1,1), and the initial velocity error follows the uniform distribution U(-0.1,0.1).

[0161] The present invention is used for a small celestial body landing trajectory optimization method under uncertain conditions. Figure 4a, 4b, 4c and 4d are respectively the position, speed, fuel consumption of the optimal fuel trajectory and the anti-disturbance trajectory, and the nominal thrust of the optimal fuel trajectory and the anti-disturbance trajectory obtained by the method of the present invention. The fuel consumption trajectory can land at the target position at zero speed under the nominal state, and the fuel consumption is increased by 0.0082% compared with the optimal fuel consumption trajectory. Figure 5 a, 5b, and 5c are the landing position variances of the x, y, and z axes, respectively. Compared with the fuel optimal trajectory, the anti-disturbance trajectory reduces the landing position variances of the three axes by 22.17%, 2.08%, and 22.16%, respectively, and the landing speed method of the three axes reduces them by 9.37%, 0.9%, and 10.75%, respectively.

[0162] The specific descriptions, technical solutions and beneficial effects described above are further explained in detail. It should be understood that the embodiments are only used to explain the present invention and are not used to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. An indirect optimization method for the anti-disturbance trajectory of a small celestial body landing, characterized in that: The method comprises the following steps: S1: Establish a stochastic dynamics model for small body landing, construct a performance functional weighted by covariance and fuel consumption, and construct an anti-disturbance trajectory optimization problem S2: The anti-disturbance trajectory optimization problem is linked to the zero-gravity energy optimal control problem by defining the homotopy coefficient performance functional, the dynamic equations, and the homotopy function of the process noise, and the initial value of the co-state of the zero-gravity energy optimal control problem is estimated. S3: The zero-gravity energy optimal control problem is converted into an anti-disturbance trajectory optimization problem by changing the homology coefficient. For each trajectory optimization problem corresponding to the homology coefficient, the Pontryagin minimum principle is used to convert it into a two-point boundary value problem and then solved by the shooting method. In the process, the optimal control that minimizes the Hamiltonian function is solved using a one-dimensional search algorithm.

2. The indirect optimization method for the anti-disturbance landing trajectory of a small celestial body according to claim 1, characterized in that: In step S1, the stochastic dynamics equation of the landing of the small celestial body is: In the above formula, X a =[r T v T mω z ] T is the state variable, r is the position of the detector, v is the velocity of the detector, m is the mass of the detector, and W is the process noise. The state covariance matrix of the above formula can be calculated along the nominal trajectory by the following formula: In the above formula, A is the Jacobian matrix of the dynamics, and Q is the power spectral density of the process noise W.

3. The indirect optimization method for the anti-disturbance landing trajectory of a small celestial body according to claim 1, characterized in that: In step S1, the following covariance and fuel consumption weighted performance functional is constructed: In the above formula is the equation of the ith state variable at the terminal moment, C i are the weight coefficients of the weights of each state variable, c0 is the weight coefficient of the distribution covariance and fuel consumption, λ0 is a performance index coefficient greater than zero, and the subscript N represents the nominal quantity. The different elements of the covariance matrix are reorganized into vector form and expanded into the state of the state equation, together with the performance functional, equality constraints and inequality constraints to form the following anti-disturbance trajectory optimization problem: In the above formula, Vector(P) is the vector form of the upper right triangular elements of the covariance matrix.

4. The indirect optimization method for the anti-disturbance landing trajectory of a small celestial body according to claim 1, characterized in that: In step S2, the following performance functionals, kinetic equations, and homotopy functions of process noise are defined respectively: When ε1=1 and ε0=0 in J1, J1 is the performance functional of the anti-disturbance trajectory optimization problem. When ε1=0 and ε2=0 in J1, J1 is the performance functional of the fuel optimal control problem. When ε1=0 and ε2=1 in J1, J1 is the performance functional of the energy optimal control problem. When f h1 (ε3) and f h2 When ε3=0 in (ε3), use f h1 (ε3), f h2 (ε3) replaces the dynamic equation of velocity and process noise, and when ε1=0 and ε2=1 in J1, J1 is the performance functional of the optimal control problem of zero-gravity energy. h1 (ε3) and f h2 When ε3=1 in (ε3), ε1=0 and ε2=1 in J1, J1 is the performance functional of the energy optimal control problem. The connection between the anti-disturbance trajectory optimization problem and the zero-gravity energy optimal control problem is established through the homotopy function.

5. The indirect optimization method for the anti-disturbance landing trajectory of a small celestial body according to claim 1, characterized in that: In step S2, the initial value of the co-state of the optimal control problem of zero-gravity energy is estimated. The initial value of the co-state of the position and velocity of the optimal control problem of zero-gravity energy is estimated by solving the analytical solution of the initial value of the co-state of the optimal control problem of zero-gravity energy without constraints: The optimal solution of the optimal control problem of zero-gravity energy is used to approximate the initial value of the mass co-state of the optimal control problem of zero-gravity control and unconstrained energy: The initial value of the covariance of the optimal control problem of unconstrained energy without gravity control is:

6. The indirect optimization method for the anti-disturbance landing trajectory of a small celestial body according to claim 1, characterized in that: In step S3, the zero-gravity energy optimal control problem is converted into an anti-disturbance trajectory optimization problem by changing the homotopy coefficient. The optimization problem corresponding to each homotopy coefficient is converted into a two-point boundary value problem and solved by the shooting method. First, let f h1 (ε3) and f h2 ε3 in (ε3) changes from 0 to 1, ε1=0, ε2=1 in J1, this process converts the weightless energy optimal control problem into the energy optimal control problem. Secondly, let ε1=0, ε2 in J1 change from 1 to 0. This process converts the energy optimal control problem into the fuel optimal control problem. Finally, let ε2=0, ε1 in J1 change from 0 to 1. This process converts the fuel optimal control problem into the anti-disturbance trajectory optimization problem. In the above process, the optimization problem corresponding to each homotopy coefficient is converted into a two-point boundary value problem using the Pontryagin minimum principle. The guessed value of the costate initial value used to solve the target shooting equation uses the costate initial value obtained from the previous two-point boundary value problem.

7. The indirect optimization method for the anti-disturbance landing trajectory of a small celestial body according to claim 1, characterized in that: Step S3 uses a one-dimensional search algorithm to search for the optimal control that minimizes the Hamiltonian function, and uses a fixed-step fourth-order Runge-Kutta (RK4) algorithm to discretize the differential equation. For each time step k, the Hamiltonian function H at the current moment is k (u) Minimize the optimal control thrust ratio A one-dimensional search algorithm is used for searching, and the initial search interval is defined as [a, b], and the two split points are: In the above formula is the split ratio, and then the search interval is continuously updated according to the search interval update algorithm until the search interval is small enough. The optimal control thrust ratio is approximately:

Images