A small-thrust trajectory optimization method near asteroids based on pseudospectral convex optimization

Abstract

The present invention discloses a small-thrust trajectory optimization method near an asteroid based on pseudospectral convex optimization, including: establishing the dynamic equation of the spacecraft's orbital motion near the asteroid, discretizing the dynamic equation by the flipped Radau pseudospectral method, establishing a problem of no-gravitational-field and no-collision constraint, and using a convex optimization method to solve and obtain an initial nominal solution; constructing a gravitational field model in the dynamic equation based on a four-particle model, and establishing a collision avoidance model near the asteroid; establishing a convex optimization problem of small-thrust trajectory for collision avoidance near the asteroid through convexification methods such as nominal value substitution and Taylor expansion; establishing a successive solution model, and obtaining the optimal solution through sequential convex optimization; the present invention adopts a pseudospectral convex optimization method combining the flipped Radau pseudospectral method and convex optimization, realizes the optimization of the small-thrust trajectory of the spacecraft's motion near the asteroid, and improves the calculation efficiency by using the four-particle gravitational field model.

Description

Technical Field

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

[0002] Due to the irregular shape of asteroids, spacecraft are prone to collide with them during close-range detection. Therefore, mission designers need to reasonably design and optimize the spacecraft orbit to meet mission requirements. In addition, the strength of the guidance ability of spacecraft for detecting asteroids is a key factor for future low-cost and lightweight missions. When detecting some asteroids far from the Earth, there will be a large communication delay between the spacecraft and the Earth, and it is difficult for the spacecraft to give timely feedback to the instructions sent from the Earth end. In recent years, people have successfully constructed many effective advanced guidance and control algorithms to achieve on-board autonomous control of spacecraft. Among them, the pseudospectral convex optimization method has been applied to the field of spacecraft trajectory control due to its high efficiency and simplicity.

[0003] The application of the pseudospectral convex optimization method is mostly in regular planet detection scenarios such as Mars landing, and the thrust level is often large. There is a lack of small-thrust trajectory optimization scenarios near asteroids based on the pseudospectral convex optimization method, and convex optimization methods are mostly used for detection trajectories near asteroids, making it difficult to achieve small-thrust trajectory optimization. Summary of the Invention

[0004] Aiming at the deficiencies in the above technologies, the present invention provides a small-thrust trajectory optimization method near an asteroid based on pseudospectral convex optimization.

[0005] To achieve the above object, the present invention provides a small-thrust trajectory optimization method near an asteroid based on pseudospectral convex optimization, including: Step 1, setting relevant parameters of the spacecraft and its motion environment; Step 2, establishing the dynamic equation of the spacecraft's motion near the asteroid, discretizing the dynamic equation by the flipped Radau pseudospectral method, establishing a problem of no-gravitational-field and no-collision-prevention constraint, and using the convex optimization method to solve it to obtain the initial nominal solution; Step 3, constructing the gravitational field model in the dynamic equation based on the four-particle model, and establishing a collision-prevention model near the asteroid; Step 4, establishing a convex optimization problem of the small-thrust trajectory for collision prevention near the asteroid through convexification methods such as nominal value substitution and Taylor expansion; Step 5, establishing a successive solution model, and obtaining the optimal solution through sequential convex optimization.

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

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

[0008] (24);

[0009] where 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 a Lagrange interpolation polynomial, using the Legendre Gauss-Radau integration points as interpolation nodes to achieve 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 polynomial degree, k represents the k-th interpolation node, and i represents the i-th basis function;

[0014] Step 22: Let the interpolated approximate function be x(τ), and the expression of the n-th order Lagrange interpolation polynomial is:

[0015] (26);

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

[0017] (27);

[0018] where represents differentiation;

[0019] The Legendre Gauss-Radau integration points are obtained by solving the roots of equation (28);

[0020] (28);

[0021] Differentiating equation (26) gives:

[0022] (29);

[0023] Define the differential matrix D:

[0024] (30);

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

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

[0027] (31);

[0028] (32);

[0029] where, A and B represent the Jacobian matrix, I represents the identity matrix, and their subscripts are the matrix dimensions. A lb and A rb are expressed as follows:

[0030] (33);

[0031] where, represents the spin angular velocity of the asteroid;

[0032] Since the interpolation nodes τ are only generated inside the interval [-1, 1], the flight time needs to be normalized to the interval [-1, 1] to obtain their corresponding relationship:

[0033] (34);

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

[0035] Define the coefficient s=(t f - t0) / 2, and substitute Equation (31) into Equation (30) to get:

[0036] (35);

[0037] The expression of the state vector X discretized along n interpolation nodes is as follows:

[0038] (36);

[0039] where, the superscript T represents the transpose;

[0040] According to Equation (35), define the state - dependent coefficient matrix A dyn and the vector b dyn :

[0041] (37);

[0042] where, diag represents the diagonal matrix and block represents the matrix block;

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

[0044] (38);

[0045] From this, the dynamic equation with discrete pseudo-spectrum is obtained:

[0046] (39);

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

[0048] (40);

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

[0050] (41);

[0051] Adopt the Gaussian integration formula to discretize the integral term in the objective function, and obtain:

[0052] (42);

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

[0054] Step 25: After reversing the order of the weight function elements, the numerical stability and calculation efficiency of the Radau pseudo-spectrum method are improved. The flipped weight function expression is:

[0055] (43);

[0056] Among them, represents the weight function without element flipping;

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

[0058] Optionally, the said 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 clustering centers to fit the gravitational acceleration at any position near the asteroid. Its expression is:

[0060] (44);

[0061] Among them, is the gravitational constant corresponding to the \(i\)-th particle, is the distance from the spacecraft to the \(i\)-th particle;

[0062] Step 32: Establish an anti-collision model near the asteroid, construct an ellipsoid to completely cover the asteroid inside the sphere, and limit the motion trajectory of the spacecraft outside the ellipsoid to achieve anti-collision; assume that the position vector from the center of the ellipsoid to the spacecraft is \(\mathbf{r}=[r_{ x , r_{ y , r_{ z \) T , and the semi-major axis lengths 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 form matrix \(\mathbf{R}=\text{diag}(1 / a_{ 2 ,1 / b_{ 2 ,1 / c_{ 2 ), then the vector expression of the anti-collision ellipsoid constraint is:

[0065] (46);

[0066] Among them, 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 first-order Taylor expansion, use the obtained nominal value to transform the anti-collision ellipsoid constraint into a convex constraint;

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

[0070] Step 43: Apply the thrust amplitude constraint to establish a convex optimization problem for the anti-collision low-thrust trajectory near the asteroid.

[0071] Optionally, step 5 includes: Using the nominal value obtained by solving the problem of no gravitational field and no anti-collision constraint as the initial value, using the CVX toolbox to solve the convex optimization problem of the anti-collision low-thrust trajectory, calculating the difference between the obtained trajectory and the nominal trajectory. If the given tolerance is satisfied, output the optimal solution. If the given tolerance is not satisfied, update the current output solution to the nominal value and iterate repeatedly 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: This small-thrust trajectory optimization method is applicable to the small-thrust trajectory optimization problem near asteroids that can be convexified; establish the dynamic equations of the spacecraft's motion near the asteroid, and discretize the dynamic equations by the flipped Radau pseudospectral method; construct a problem of no gravitational field and no collision avoidance constraint, and obtain the nominal solution using the convex optimization method; add the four-particle gravitational field model to the orbital dynamics to improve the calculation efficiency based on the classical polyhedron gravitational field model; impose an ellipsoidal constraint to achieve collision avoidance of the spacecraft trajectory; establish a small-thrust trajectory pseudospectral convex optimization problem through various convexification techniques; apply the successive solution model, and obtain the optimal solution through multiple solution iterations. Brief Description of the Drawings

[0073] Figure 1 It is a schematic flow chart of the method of the present invention. Detailed Embodiments

[0074] For the convenience of understanding by those skilled in the art, the present invention will be further described below in conjunction with the embodiments and the drawings. The content mentioned in the embodiments does not limit the present invention.

[0075] Refer to Figure 1 As shown, this embodiment provides a small-thrust trajectory optimization method near asteroids based on pseudospectral convex optimization, and the steps are as follows:

[0076] 1) Set the relevant parameters of the spacecraft and its motion environment.

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

[0078] 2) Establish a problem of no gravitational field and no collision avoidance constraint through the flipped Radau pseudospectral method to obtain the nominal solution.

[0079] Specifically, establish the dynamic equations of the spacecraft's orbital motion near the asteroid, and use the flipped Radau pseudospectral method. Use the Lagrange interpolation polynomial based on the Legendre Gauss-Radau integration points for global approximate fitting of the dynamics, and perform element flipping on the weight function in the objective function. This method has the dual advantages of convergence speed and calculation accuracy compared with other pseudospectral methods, and synchronously improves the numerical stability and calculation efficiency through the flipped processing. Impose a thrust amplitude constraint, and thus establish a small-thrust trajectory optimization problem of fuel optimization without gravitational field and without collision avoidance constraint. Obtain the nominal solution through a single optimization using the MATLAB convex optimization toolbox CVX.

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

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

[0082] 4) Establish a convex optimization problem for the small-thrust trajectory for collision avoidance near the asteroid through various convexification methods.

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

[0084] 5) Based on the nominal solution obtained in step 2), establish a successive solution model and iterate 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 for solution 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 takes the equilibrium point transfer trajectory of asteroid 1996HW1 as an example for illustration:

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

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

[0089] (47);

[0090] Among them, the spacecraft position and velocity state x = [r x , r y , r z , v x , v y , v z T , the control quantity u provided by the spacecraft thruster = [u x , u y , u z T , ω ast is the spin angular velocity of the asteroid, and t is the flight time. Linearize the dynamics equation and discretize it along n integration points by the flipped Radau pseudospectral method to obtain the pseudospectral discretized dynamics equation:

[0091] (48);

[0092] (49);

[0093] Apply the thrust amplitude constraint:

[0094] (50);

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

[0096] (51);

[0097] Among them, T is the thrust vector. Calculated according to the maximum thrust T max for 1.5 hours of propulsion, the fuel consumption of the spacecraft is only 0.14491 kg, which is almost negligible compared to the total mass of the spacecraft. Therefore, this paper does not consider the impact of fuel consumption on the dynamics, and the mass m in Equation takes a fixed value.

[0098] ​​The problem of optimizing the small-thrust trajectory without gravitational field and collision avoidance constraints with optimal fuel is established. The objective function is discretized by Gaussian quadrature and a weight function ω with reversed element order is applied to obtain:

[0099] (52);

[0100] Subsequently, the nominal solution is obtained through a single optimization using the CVX toolbox in MATLAB.

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

[0102] Step 4: Through first-order Taylor expansion and using the obtained nominal values, transform the collision avoidance ellipsoid constraint into a convex constraint; replace the gravitational acceleration term with the nominal value to convexify the dynamics equation; apply the thrust magnitude constraint to establish the convex optimization problem of the small-thrust trajectory for collision avoidance near the asteroid.

[0103] Step 5: A nominal solution including nominal values of states, controls, etc. can be obtained from Step 2. Subsequently, use the CVX toolbox to solve the convex optimization problem of the small-thrust trajectory for collision avoidance near the asteroid established in Step 4, calculate the difference between the obtained trajectory and the nominal trajectory. If the given tolerance is satisfied, output the optimal solution. If the given tolerance is not satisfied, update the current output solution to the nominal value and iterate repeatedly 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 dynamics equation of the spacecraft's orbital motion near the asteroid, discretizes the dynamics equation by the flipped Radau pseudospectral method; constructs the problem of no gravitational field and no collision avoidance constraints, and uses the convex optimization method to obtain the nominal solution; adds the four-particle gravitational field model to the orbital dynamics to improve the calculation efficiency on the basis of the classical polyhedron gravitational field model; applies the ellipsoidal constraint to achieve collision avoidance of the spacecraft trajectory; establishes the pseudospectral convex optimization problem of the small-thrust trajectory through various convexification techniques; applies the successive solution model, and obtains the optimal solution through multiple solution iterations.

[0105] There are many specific application ways of the present invention. The above are only the preferred embodiments of the present invention. It should be noted that for those of ordinary skill in the art in this technical field, without departing from the principle of the present invention, several improvements can still be made, and these improvements should also be regarded as the protection scope of the present invention.

Claims

1. A small-thrust trajectory optimization method near an asteroid based on pseudospectral convex optimization, characterized in that Including: Step 1: Set parameters related to the spacecraft and its motion environment; the related parameters include one or more of the spacecraft flight time, the number of discretization points of the pseudospectral method, the four-particle gravitational field model parameters, the asteroid spin acceleration, the anti-collision model parameters, the initial and final position and velocity states of the spacecraft, the maximum thrust, and the convergence tolerance; Step 2: Establish the dynamic equation of the spacecraft's motion near the asteroid. Discretize the dynamic equation by the flipped Radau pseudospectral method, establish the problem of no gravitational field and no anti-collision constraint, and use the convex optimization method to solve it to obtain the initial nominal value; Step 3: Construct the four-particle gravitational field model in the dynamic equation based on the four-particle model, and establish the anti-collision model near the asteroid; Including: Step 31: Establish the four-particle gravitational field model, and add the gravitational acceleration term g(r) to the dynamic equation; generate the particle group data through the asteroid polyhedron model, and use the K-means clustering algorithm based on the data set to obtain four clustering centers to fit the gravitational acceleration at any position near the asteroid. Its expression is: （21）； Among them, is the gravitational constant corresponding to the i-th 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 limit 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 . If the semi-major axis lengths of the ellipsoid along the x, y, and z axes are a, b, and c respectively, the collision avoidance ellipsoid constraint can be expressed as:​ （22）； Define the quadratic form matrix \(R = \text{diag}(1 / a 2 , 1 / b 2 , 1 / c 2 ). Then the vector expression of the anti-collision ellipsoid constraint is: （23）； Among them, represents the transpose of the position vector from the center of the ellipsoid to the spacecraft; Step 4: Establish the convex optimization problem of the anti-collision small-thrust trajectory near the asteroid by the convex optimization method of nominal value substitution and Taylor expansion; Step 5: Establish a successive solution model and obtain the optimal solution through sequential convex optimization.

2. The small-thrust trajectory optimization method near an asteroid according to claim 1, wherein In step 2, the dynamic equation is established as: （1）； where 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.

3. The small-thrust trajectory optimization method near an asteroid according to claim 2, wherein Step 2 includes: Step 21: Construct the Lagrange interpolation polynomial, use the Legendre Gauss-Radau integration points as the interpolation nodes, and realize the interpolation approximation of the dynamic equation; let the interpolation node be τ, then the expression of the Lagrange basis function is: （2）； where n represents the polynomial degree, k represents the k-th interpolation node, and i represents the i-th basis function; Step 22: Let the interpolated approximate function be x(τ), and the expression of the n-th Lagrange interpolation polynomial is: （3）； Consider the sequence of n - th Legendre orthogonal polynomials \(L\) n (τ), whose expression in the interval \([-1,1]\) is given by: （4）； Among them, represents differentiation; The Legendre Gauss-Radau quadrature points are obtained by solving the equation for its roots; （5）； Taking the derivative of equation (3) gives: （6）； Define the differential matrix D: （7）； where D represents an n×(n + 1)-dimensional matrix, and D j,i represents the element in the i-th column and j-th row of matrix D; Step 23: According to equation (1), the state-space expression of the dynamic equation is: （8）； （9）； where A and B denote Jacobian matrices, I denotes the identity matrix, and their subscripts are the matrix dimensions, A lb and A rb The expressions are as follows: （10） Among them, represents the spin angular velocity of the asteroid; Since the interpolation nodes τ are only generated inside the interval [-1, 1], the time of flight needs to be normalized to the interval [-1, 1] to obtain the corresponding relationship between the two: （11）； Among them, represents the initial flight time; represents the end flight time; Define the coefficient s = (t f - t0) / 2, substitute Equation (8) into Equation (7) to obtain: （12）； The expression of the state vector X discretized along n interpolation nodes is as follows: （13）； where the superscript T represents the transpose; Define the state-dependent coefficient matrix A according to Equation (12) dyn and the vector b dyn : （14）； where diag represents the diagonal matrix and block represents the matrix block; Matrix A dyn The diagonal block matrix and other block matrices are defined as follows: （15）； Thus, the pseudospectral-discretized dynamic equation is obtained: （16）； Step 24: Take the integral of the control quantity u during the flight time as the optimization objective function to represent the fuel optimization, and given the boundary conditions of the initial and final states, establish the problem of no gravitational field and no anti-collision constraint; among them, the terminal state constraint is expressed as: （17）； Matrix A f and vector b f The expressions are respectively: （18）； Use the Gaussian integration formula to discretize the integral term in the objective function to obtain: （19）； where ω i represents the weight function component; After reversing the order of the weight function elements, the numerical stability and calculation efficiency of the Radau pseudospectral method are improved, and the flipped weight function expression is: （20）； Among them, represents the weight function for non-element flipping; After establishing the problem of no gravitational field and no anti-collision constraint, use the MATLAB convex optimization toolbox CVX to solve and obtain the initial nominal value.

4. The small-thrust trajectory optimization method near an asteroid according to claim 1, wherein Step 4 includes: Step 41: Through the first-order Taylor expansion, use the obtained nominal value to transform the anti-collision ellipsoid constraint into a convex constraint; Step 42: Replace the gravitational acceleration term with the nominal value to achieve the convex optimization of the dynamic equation; Step 43: Apply the thrust magnitude constraint to establish a convex optimization problem for the small-thrust trajectory for collision avoidance near the asteroid.

5. The small-thrust trajectory optimization method near an asteroid according to claim 1, wherein Step 5 includes: Using the nominal value obtained by solving the problem without gravitational field and without collision avoidance constraint as the initial value, solving the convex optimization problem for the small-thrust trajectory for collision avoidance using the CVX toolbox, calculating the difference between the obtained trajectory and the nominal trajectory. If the given tolerance is satisfied, output the optimal solution. If the given tolerance is not satisfied, update the current output solution to the nominal value, and iterate repeatedly until the convergence condition is met to obtain the optimal solution.

Images