Method for optimizing low-thrust trajectory near asteroid based on pseudo-spectral convex optimization

Abstract

The invention discloses a pseudo-spectral convex optimization-based method for optimizing a low-thrust trajectory near an asteroid, which comprises the following steps of: establishing a kinetic equation of an orbital motion of a spacecraft near the asteroid, discretizing the kinetic equation through a flipped Radau pseudo-spectral method, establishing a problem of no gravitational field and no anti-collision constraint, and solving by using a convex optimization method to obtain an initial nominal solution; constructing a gravitational field model in the kinetic equation based on a four-mass-point model, and establishing an asteroid nearby anti-collision model; establishing a convex optimization problem of an anti-collision small-thrust trajectory near the asteroid through convex methods such as nominal value replacement and Taylor expansion; establishing a successive solution model, and obtaining an optimal solution through sequence convex optimization; according to the method, a flipped Radau pseudo-spectrum method and convex optimization combined pseudo-spectrum convex optimization method is adopted, the motion low-thrust trajectory optimization of the spacecraft near the asteroid is realized, and the calculation efficiency is improved by using a four-mass-point gravitational field model.

Description

Technical Field

[0001] The present invention belongs to the field of aerospace technology, and specifically refers to a method for optimizing a small-thrust trajectory near an asteroid based on pseudo-spectral convex optimization. Background Art

[0002] Due to the irregular shape of asteroids, spacecraft are prone to collision with them during close-range detection. Therefore, mission designers need to rationally design and optimize the spacecraft orbit to meet mission requirements. In addition, the strength of the guidance capability of spacecraft for detecting asteroids is a key factor in achieving low-cost and lightweight missions in the future. When detecting some asteroids that are far away from the earth, there will be a large communication delay between the spacecraft and the earth, and it is difficult for the spacecraft to respond to the instructions issued by the earth in a timely manner. In recent years, many effective advanced guidance and control algorithms have been successfully constructed to achieve autonomous on-board control of spacecraft. Among them, the pseudo-spectral convex optimization method has been applied to the field of spacecraft trajectory control due to its high efficiency and simplicity.

[0003] The application of pseudo-spectral convex optimization methods is mostly in regular planetary exploration scenarios such as Mars landing, and the thrust magnitude is often large. There is a lack of small-thrust trajectory optimization scenarios near asteroids based on pseudo-spectral convex optimization methods, and most exploration trajectories near asteroids use convex optimization methods, which makes it difficult to achieve small-thrust trajectory optimization. Summary of the invention

[0004] In view of the deficiencies in the above-mentioned technology, the present invention provides a method for optimizing a small-thrust trajectory near an asteroid based on pseudo-spectral convex optimization.

[0005] To achieve the above-mentioned purpose, the present invention provides a method for optimizing a small-thrust trajectory near an asteroid based on pseudo-spectral convex optimization, comprising: step 1, setting parameters related to a spacecraft and its motion environment; step 2, establishing a dynamic equation of the spacecraft moving near an asteroid, discretizing the dynamic equation by the flipped Radau pseudo-spectral method, establishing a problem without a gravitational field and without anti-collision constraints, solving it using a convex optimization method, and obtaining an initial nominal solution; step 3, constructing a gravitational field model in the dynamic equation based on a four-particle model, and establishing an anti-collision model near an asteroid; step 4, establishing a convex optimization problem of a small-thrust trajectory for anti-collision near an asteroid by convexification methods such as nominal value replacement and Taylor expansion; step 5, establishing a successive solution model, and obtaining the optimal solution by sequential convex optimization.

[0006] Optionally, in step 1, the relevant parameters include one or more of the spacecraft flight time, the number of discrete points of the pseudo-spectral method, the four-particle gravitational field model parameters, the asteroid spin acceleration, the anti-collision model parameters, the spacecraft initial and final position velocity states, the maximum thrust, and the convergence tolerance.

[0007] Optionally, in step 2, the orbital dynamics equation is established as:

[0008] (twenty four);

[0009] Among them, x is the position and velocity state of the spacecraft, u is the control quantity provided by the spacecraft thruster, and t is the flight time.

[0010] Optionally, step 2 includes:

[0011] Step 21: Construct the Lagrange interpolation polynomial and use the Legendre Gauss-Radau integration point as the interpolation node to achieve the interpolation approximation of the orbital dynamics equation; let the interpolation node be τ, then the expression of the Lagrange basis function is:

[0012] (25);

[0013] Where n represents the degree of the polynomial, k represents the kth interpolation node, and i represents the ith basis function;

[0014] Step 22: Assume that the interpolation approximation function is x(τ), and the n-th Lagrange interpolation polynomial expression is:

[0015] (26);

[0016] Consider the n-order Legendre orthogonal polynomial sequence L n (τ), whose expression in the interval [-1,1] is:

[0017] (27);

[0018] in, represents differential;

[0019] The Legendre Gauss-Radau integration point is obtained by taking the root of equation (28);

[0020] (28);

[0021] By taking the derivative of formula (26), we can obtain:

[0022] (29);

[0023] Define the differential matrix D:

[0024] (30);

[0025] Where D represents an n×(n+1)-dimensional matrix, D j,i represents the j-th row and i-th column of matrix D;

[0026] Step 23: According to equation (24), the state space expression of the dynamic equation is:

[0027] (31);

[0028] (32);

[0029] Among them, A and B represent Jacobian matrices, I represents the identity matrix, and its subscript is the matrix dimension. lb and A rb The expression is as follows:

[0030] (33);

[0031] in, represents the asteroid's spin angular velocity;

[0032] Since the interpolation node τ is only generated within the interval [-1,1], the flight time needs to be Normalized to the interval [-1,1], we get the corresponding relationship between the two:

[0033] (34);

[0034] in, represents the initial flight time; represents the terminal flight time;

[0035] Define the coefficient s=(t f - t0) / 2, substituting equation (31) into equation (30) yields:

[0036] (35);

[0037] The discrete state vector X along n interpolation nodes is expressed as follows:

[0038] (36);

[0039] Among them, the superscript T represents transposition;

[0040] According to formula (35), the state correlation coefficient matrix A is defined as dyn and vector b dyn :

[0041] (37);

[0042] Among them, diag represents a diagonal matrix, and block represents a matrix block;

[0043] Matrix A dynThe diagonal block matrix and other block matrices of are defined as follows:

[0044] (38);

[0045] This results in the pseudo-spectral discretized dynamic equation:

[0046] (39);

[0047] Step 24: Take the integral of the control quantity u during the flight time as the optimization objective function to express the optimal fuel. Given the boundary conditions of the initial and final states, establish a problem without gravity field and anti-collision constraints; the terminal state constraint is expressed as:

[0048] (40);

[0049] Matrix A f and vector b f The expressions are:

[0050] (41);

[0051] The Gaussian integral formula is used to discretize the integral term in the objective function, and we get:

[0052] (42);

[0053] Among them, ω i represents the weight function component;

[0054] Step 25: After the order of weight function elements is reversed, the numerical stability and computational efficiency of Radau pseudospectral method are improved. The expression of flipped weight function is:

[0055] (43);

[0056] in, represents the weight function of the unflipped element;

[0057] After establishing the problem without gravitational field and anti-collision constraints, the MATLAB convex optimization toolbox CVX is used to solve and obtain the initial nominal value.

[0058] Optionally, step 3 includes:

[0059] Step 31: Establish a four-particle gravitational field model and add the gravitational acceleration term g(r) to the orbital dynamics equation; generate particle group data through the asteroid polyhedron model, and use the K-means clustering algorithm based on the data set to obtain four cluster centers to fit the gravitational acceleration at any position near the asteroid. The expression is:

[0060] (44);

[0061] in, is the gravitational constant corresponding to the ith particle, is the distance from the spacecraft to the i-th particle;

[0062] Step 32: Establish a collision avoidance model near the asteroid, construct an ellipsoid to completely cover the asteroid inside the sphere, and restrict the motion trajectory of the spacecraft outside the ellipsoid to achieve collision avoidance; assume that the position vector from the center of the ellipsoid to the spacecraft is r=[r x , r y , r z ] T , the lengths of the semi-major axes of the ellipsoid along the x, y, and z axes are a, b, and c respectively, then the anti-collision ellipsoid constraint can be expressed as:

[0063] (45);

[0064] Define the quadratic matrix R = diag(1 / a 2 ,1 / b 2 ,1 / c 2 ), then the vector expression of the anti-collision ellipsoid constraint is:

[0065] (46);

[0066] in, represents the transpose of the position vector from the center of the ellipsoid to the spacecraft;

[0067] Optionally, step 4 includes:

[0068] Step 41: Through the first-order Taylor expansion, the anti-collision ellipsoid constraint is converted into a convex constraint using the solved nominal value;

[0069] Step 42: Replace the gravitational acceleration term with a nominal value to achieve convexification of the dynamic equation;

[0070] Step 43: Apply thrust amplitude constraints and establish a convex optimization problem for collision avoidance and low-thrust trajectories near the asteroid.

[0071] Optionally, step 5 includes: using the nominal value obtained by solving the problem without gravitational field and anti-collision constraints as the initial value, using the CVX toolbox to solve the anti-collision small thrust trajectory convex optimization problem, calculating the difference between the solved trajectory and the nominal trajectory, and outputting the optimal solution if the given tolerance is met; if the given tolerance is not met, updating the output solution to the nominal value, and repeatedly iterating until the convergence condition is met to obtain the optimal solution.

[0072] The beneficial effects of the present invention compared with the prior art are as follows: the present small thrust trajectory optimization method is applicable to the small thrust trajectory optimization problem near an asteroid that can be convexified; the dynamic equations of the motion of a spacecraft near an asteroid are established, and the dynamic equations are discretized by the flipped Radau pseudo-spectral method; a problem without a gravitational field and without anti-collision constraints is constructed, and a nominal solution is obtained by using a convex optimization method; a four-particle gravitational field model is added to orbital dynamics to improve the computational efficiency based on the classical polyhedron gravitational field model; an ellipsoid constraint is imposed to achieve spacecraft trajectory anti-collision; a pseudo-spectral convex optimization problem of a small thrust trajectory is established by using a variety of convexification techniques; a successive solution model is applied, and the optimal solution is obtained through multiple solution iterations. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] Figure 1 It is a schematic diagram of the process of the present invention. DETAILED DESCRIPTION

[0074] In order to facilitate the understanding of those skilled in the art, the present invention is further described below in conjunction with embodiments and drawings. The contents mentioned in the implementation modes are not intended to limit the present invention.

[0075] Reference Figure 1 As shown, this embodiment provides a method for optimizing a small thrust trajectory near an asteroid based on pseudo-spectral convex optimization, and the steps are as follows:

[0076] 1) Set the parameters related to the spacecraft and its motion environment.

[0077] Specifically, the state parameters include spacecraft flight time, number of discrete points of pseudo-spectral method, four-particle gravitational field model parameters, asteroid spin acceleration, anti-collision model parameters, spacecraft initial and final position velocity state, maximum thrust, and convergence tolerance.

[0078] 2) The flipped Radau pseudospectral method is used to establish a problem without gravitational field and without anti-collision constraints, and the nominal solution is obtained.

[0079] Specifically, the dynamic equations of the spacecraft's orbital motion near the asteroid are established, and the flipped Radau pseudo-spectral method is used to globally approximate the dynamics using the Lagrange interpolation polynomial based on the Legendre Gauss-Radau integral point, and the weight function in the objective function is element-flipped. Compared with other pseudo-spectral methods, this method has the dual advantages of convergence speed and calculation accuracy, and the flipped processing simultaneously improves the numerical stability and calculation efficiency. The thrust amplitude constraint is imposed, thereby establishing the fuel-optimal low-thrust trajectory optimization problem without gravity field and anti-collision constraint. The nominal solution is obtained through single optimization using the MATLAB convex optimization toolbox CVX.

[0080] 3) Establish a four-particle gravitational field model and anti-collision model.

[0081] Specifically, a four-particle gravitational field model is used. Based on the polyhedron gravitational field model, four particles are used to fit the gravitational field, and the gravitational field model is added to the orbital dynamics. The asteroid is included in the ellipsoid model, and constraints are imposed so that the spacecraft can only move outside the ellipsoid, thereby establishing an anti-collision model.

[0082] 4) The convex optimization problem of small-thrust collision avoidance trajectory near the asteroid is established through a variety of convexification methods.

[0083] Specifically, the nominal value of gravitational acceleration is used to replace the non-convex gravitational field model in the dynamic model, and the orbital dynamics equation is transformed into a convex constraint. The anti-collision ellipsoid constraint is convexified and transformed into a convex constraint through the first-order Taylor expansion. In this way, the pseudo-spectral convex optimization problem of small thrust trajectory is established.

[0084] 5) Based on the nominal solution obtained in step 2), a successive solution model is established and iterated until convergence to obtain the optimal solution.

[0085] Specifically, the solution obtained in step 2) is used as a reference value, and step 4) is repeated to solve the problem until the difference between the two iterative solutions meets the convergence condition, that is, the difference between the two consecutive solution trajectories is less than the convergence tolerance, thereby obtaining the optimal solution to the pseudo-spectral convex optimization problem of the small thrust trajectory.

[0086] The following is an example of the equilibrium point transfer trajectory of asteroid 1996 HW1:

[0087] Step 1. First, in this problem, it is assumed that the spacecraft is initially hovering at a certain equilibrium point near the asteroid 1996 HW1, and the spacecraft can maintain an initial velocity of 0 at the equilibrium point. Assume that the spacecraft is equipped with electric thrusters that can provide a maximum thrust of 1N. After 1.5 hours of flight, it reaches another equilibrium point near the asteroid and hovers for observation. At present, my country has developed Hall thrusters that can provide 4.6N continuous thrust, so the assumption of small thrust is reasonable. Based on the above background, set the flight time t f =1.5 h, number of integration points n=120, asteroid spin angular velocity ω=1.9929×10 -4 1 / s, the initial position of the spacecraft r0=[-3.26866,0.08414,-0.00103] T km, spacecraft initial speed v0=[0,0,0] T km / s, spacecraft terminal position r f =[3.21197,0.13383,-0.00233] T km, spacecraft terminal position v f =[0,0,0] T km / s, specific impulse I sp=3800 m / s, spacecraft mass m=1000 kg, Earth's gravitational acceleration g e =9.80665 m / s 2 , maximum thrust T max =1 N, convergence tolerance ε=0.1 m.

[0088] Step 2: Establish the orbital dynamics equation without gravitational field as follows:

[0089] (47);

[0090] Among them, the spacecraft position velocity state x = [r x ,r y ,r z ,v x ,v y ,v z ] T , the control quantity provided by the spacecraft thruster u = [u x ,u y ,u z ] T ,ω ast is the asteroid's spin angular velocity, and t is the flight time. The dynamic equation is linearized and discretized along n integral points using the flipped Radau pseudo-spectral method to obtain the pseudo-spectral discretized dynamic equation:

[0091] (48);

[0092] (49);

[0093] Apply thrust magnitude constraint:

[0094] (50);

[0095] Among them, u is the control quantity that the thruster can provide, which represents the acceleration in this problem. Since the electric thruster has a large specific impulse, the expression for the fuel consumption rate is obtained according to the Tsiolkovsky equation:

[0096] (51);

[0097] Where T is the thrust vector. According to the maximum thrust T max The fuel consumption of the spacecraft is only 0.14491 kg when the propulsion is carried out for 1.5 hours, which is almost negligible compared with the total mass of the spacecraft. Therefore, this paper does not consider the effect of fuel consumption on dynamics. The mass m in is taken as a fixed value.

[0098] The fuel-optimal low-thrust trajectory optimization problem without gravity field and anti-collision constraints is established. The objective function is discretized by Gaussian integral and the weight function ω with reversed element order is applied to obtain:

[0099] (52);

[0100] The nominal solution was then obtained through a single optimization using the MATLAB convex optimization toolbox CVX.

[0101] Step 3: Establish a four-particle gravitational field model and add the gravitational acceleration term g(r) to the orbital dynamics equation. Apply an ellipsoid anti-collision constraint to restrict the optimized trajectory outside the ellipsoid.

[0102] Step 4: Through the first-order Taylor expansion, use the nominal value obtained by the solution to transform the anti-collision ellipsoid constraint into a convex constraint; replace the gravitational acceleration term with the nominal value to realize the convexification of the dynamic equation; impose the thrust amplitude constraint and establish the convex optimization problem of the small thrust trajectory for anti-collision near the asteroid.

[0103] Step 5. A nominal solution can be obtained from step 2, including nominal values ​​of state, control, etc. Then, the CVX toolbox is used to solve the convex optimization problem of small-thrust trajectory for collision avoidance near the asteroid established in step 4, and the difference between the solved trajectory and the nominal trajectory is calculated. If the given tolerance is met, the optimal solution is output. If the given tolerance is not met, the output solution is updated to the nominal value, and the iteration is repeated until the convergence condition is met to obtain the optimal solution.

[0104] It can be seen from the above embodiments that the present invention establishes the dynamic equations of the orbital motion of the spacecraft near the asteroid, and discretizes the dynamic equations through the flipped Radau pseudo-spectral method; constructs a problem without gravitational field and anti-collision constraints, and uses a convex optimization method to obtain a nominal solution; adds a four-particle gravitational field model to orbital dynamics, and improves the calculation efficiency on the basis of the classical polyhedron gravitational field model; imposes an ellipsoid constraint to achieve spacecraft trajectory anti-collision; establishes a small thrust trajectory pseudo-spectral convex optimization problem through a variety of convexification techniques; applies a successive solution model, and obtains the optimal solution through multiple solution iterations.

[0105] The present invention has many specific application paths. The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements can be made without departing from the principle of the present invention. These improvements should also be regarded as the protection scope of the present invention.

Claims

1. A method for optimizing a small thrust trajectory near an asteroid based on pseudo-spectral convex optimization, characterized in that: include: Step 1: Set the spacecraft and its motion environment related parameters; Step 2: Establish the dynamic equations of the spacecraft moving near the asteroid, discretize the dynamic equations by the flipped Radau pseudo-spectral method, establish a problem without gravitational field and anti-collision constraints, solve it using the convex optimization method, and obtain the initial nominal solution; Step 3: Construct a gravitational field model in the dynamic equation based on the four-particle model and establish a collision avoidance model near the asteroid; Step 4: Establish the convex optimization problem of the small-thrust trajectory for collision avoidance near the asteroid through convexification methods such as nominal value replacement and Taylor expansion; Step 5: Establish a successive solution model and obtain the optimal solution through sequential convex optimization.

2. The method for optimizing a small-thrust trajectory near an asteroid according to claim 1, characterized in that: In step 1, the relevant parameters include one or more of the spacecraft flight time, the number of discrete points of the pseudo-spectral method, the four-particle gravitational field model parameters, the asteroid spin acceleration, the anti-collision model parameters, the spacecraft initial and final position velocity states, the maximum thrust, and the convergence tolerance.

3. The method for optimizing a small thrust trajectory near an asteroid according to claim 1, characterized in that: In step 2, the orbital dynamics equation is established as: （1）； Among them, x is the position and velocity state of the spacecraft, u is the control quantity provided by the spacecraft thruster, and t is the flight time.

4. The method for optimizing a small thrust trajectory near an asteroid according to claim 3, characterized in that: The step 2 comprises: Step 21: Construct the Lagrange interpolation polynomial and use the Legendre Gauss-Radau integration point as the interpolation node to achieve the interpolation approximation of the orbital dynamics equation; let the interpolation node be τ, then the expression of the Lagrange basis function is: （2）； Where n represents the degree of the polynomial, k represents the kth interpolation node, and i represents the ith basis function; Step 22: Assume that the interpolation approximation function is x(τ), and the n-th Lagrange interpolation polynomial expression is: （3）； Consider the n-order Legendre orthogonal polynomial sequence L n (τ), whose expression in the interval [-1,1] is: （4）； in, represents differential; The Legendre Gauss-Radau integration points are given by the equation Find the root; （5）； By taking the derivative of formula (3), we can get: （6）； Define the differential matrix D: （7）； Where D represents an n×(n+1)-dimensional matrix, D j,i represents the j-th row and i-th column of matrix D; Step 23: According to equation (1), the state space expression of the dynamic equation is: （8）； （9）； Among them, A and B represent Jacobian matrices, I represents the identity matrix, and its subscript is the matrix dimension. lb and A rb The expression is as follows: （10） in, represents the asteroid's spin angular velocity; Since the interpolation node τ is only generated within the interval [-1,1], the flight time needs to be Normalized to the interval [-1,1], we get the corresponding relationship between the two: （11）； in, represents the initial flight time; represents the terminal flight time; Define the coefficient s=(t f - t0) / 2, substituting equation (8) into equation (7) yields: （12）； The discrete state vector X along n interpolation nodes is expressed as follows: （13）； Among them, the superscript T represents transposition; According to formula (12), the state correlation coefficient matrix A is defined as dyn and vector b dyn : （14）； Among them, diag represents a diagonal matrix, and block represents a matrix block; Matrix A dyn The diagonal block matrix and other block matrices of are defined as follows: （15）； This results in the pseudo-spectral discretized dynamic equation: （16）； Step 24: Take the integral of the control quantity u during the flight time as the optimization objective function to express the optimal fuel. Given the boundary conditions of the initial and final states, establish a problem without gravity field and anti-collision constraints; the terminal state constraint is expressed as: （17）； Matrix A f and vector b f The expressions are: （18）； The Gaussian integral formula is used to discretize the integral term in the objective function, and we get: （19）； Among them, ω i represents the weight function component; Step 25: After the order of weight function elements is reversed, the numerical stability and computational efficiency of Radau pseudospectral method are improved. The expression of flipped weight function is: （20）； in, represents the weight function of the unflipped element; After establishing the problem without gravitational field and anti-collision constraints, the MATLAB convex optimization toolbox CVX is used to solve and obtain the initial nominal value.

5. The method for optimizing a small-thrust trajectory near an asteroid according to claim 1, characterized in that: The step 3 comprises: Step 31: Establish a four-particle gravitational field model and add the gravitational acceleration term g(r) to the orbital dynamics equation; generate particle group data through the asteroid polyhedron model, and use the K-means clustering algorithm based on the data set to obtain four cluster centers to fit the gravitational acceleration at any position near the asteroid. The expression is: （21）； in, is the gravitational constant corresponding to the ith particle, is the distance from the spacecraft to the i-th particle; Step 32: Establish a collision avoidance model near the asteroid, construct an ellipsoid to completely cover the asteroid inside the sphere, and restrict the motion trajectory of the spacecraft outside the ellipsoid to achieve collision avoidance; assume that the position vector from the center of the ellipsoid to the spacecraft is r=[r x , r y , r z ] T , the lengths of the semi-major axes of the ellipsoid along the x, y, and z axes are a, b, and c respectively, then the anti-collision ellipsoid constraint can be expressed as: （22）； Define the quadratic matrix R = diag(1 / a 2 ,1 / b 2 ,1 / c 2 ), then the vector expression of the anti-collision ellipsoid constraint is: （23）； in, Represents the transpose of the position vector from the center of the ellipsoid to the spacecraft.

6. The method for optimizing a small thrust trajectory near an asteroid according to claim 1, characterized in that: The step 4 comprises: Step 41: Through the first-order Taylor expansion, the anti-collision ellipsoid constraint is converted into a convex constraint using the solved nominal value; Step 42: Replace the gravitational acceleration term with a nominal value to achieve convexification of the dynamic equation; Step 43: Apply thrust amplitude constraints and establish a convex optimization problem for collision avoidance and low-thrust trajectories near the asteroid.

7. The method for optimizing a small-thrust trajectory near an asteroid according to claim 1, characterized in that: Step 5 includes: using the nominal value obtained by solving the problem without gravitational field and anti-collision constraint as the initial value, using the CVX toolbox to solve the anti-collision small thrust trajectory convex optimization problem, calculating the difference between the solved trajectory and the nominal trajectory, and outputting the optimal solution if the given tolerance is met; if the given tolerance is not met, updating the output solution to the nominal value, and iterating repeatedly until the convergence condition is met to obtain the optimal solution.

Images