A sequential convex optimization method for lunar landing braking phase 6-DoF trajectory

Abstract

The application relates to a sequence convex optimization method for a lunar landing braking phase 6-DoF trajectory, and belongs to the technical field of spaceflight. The method comprises the following steps: establishing a 6-DoF dynamics model of a lander; establishing an optimization problem model according to landing process constraints to obtain a continuous optimal control problem P1; performing convex processing on the 6-DoF dynamics model of the lander, and discretizing and converting the continuous optimal control problem P1 into a sequence convex programming problem P2; and solving the sequence convex programming problem P2 by using a sequence convex programming method. The method considers the moon rotation and a non-constant gravity field, establishes a lander attitude-orbit coupling dynamics, and can effectively solve a 6-DoF optimal trajectory planning problem of a lunar landing braking phase, so that the lunar landing precision is improved.

Description

Technical Field

[0001] This invention belongs to the field of aerospace technology and relates to a 6-DoF trajectory planning method for lunar landing braking phase, particularly a six-degree-of-freedom attitude-orbit coupled landing dynamics trajectory planning method that considers lunar rotation and non-constant gravitational field. Background Technology

[0002] Lunar soft landing typically consists of four phases: braking, rapid attitude adjustment, approach, and final descent. The braking phase is the longest and has the greatest impact on range, requiring maximum accuracy at the end of the braking process while minimizing fuel consumption. The braking phase has long employed a 3-DoF guidance system, where the lander's attitude is uniquely determined by the thrust vector, and the calculated trajectory only satisfies the attitude constraints represented by the thrust vector. Therefore, extensive simulations are required before mission execution to verify the feasibility of the calculated trajectory for a 6-DoF physical system. Furthermore, for landing missions in complex scenarios, more attitude constraints need to be applied to ensure the sensor field of view. Therefore, a 6-DoF powered descent trajectory that considers both translational and attitude dynamics needs to be designed. The paper "Lee U, Mesbahiy M. Optimal powered descent guidance with 6-DoF line of sight constraints via unit dual quaternions[C]. AIAA Guidance, Navigation, and Control Conference, Kissimmee, USA, January 5-9, 2015." discloses a model predictive control method. The accuracy of the dynamic equations described in this paper largely depends on the discrete time step. To ensure system stability, the prediction time domain length needs to be limited, leading to a complex optimization problem. Furthermore, additional conditions need to be added to ensure system stability and feasibility. The paper "Szmuk M, Eren U, B. Successive convexification for Mars 6-DoF powered descent landing guidance[C]. AIAA Guidance, Navigation, and Control Conference, Grapevine, USA, January 9-13, 2017. This paper proposes a sequential convex programming guidance algorithm for the six-DOF powered descent problem on Mars, considering nonlinear aerodynamic lift and drag. The method described in the literature studies the terminal phase of Mars landing but does not consider Mars' rotation and non-constant gravitational field. Summary of the Invention

[0003] The technical problem to be solved by this invention is:

[0004] To address the attitude-orbit coupling problem during the lunar landing braking phase and improve landing accuracy, this invention provides a sequential convex optimization method for the 6-DoF trajectory during the lunar landing braking phase.

[0005] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is as follows:

[0006] A sequential convex optimization method for the 6-DoF trajectory of a lunar landing braking phase, characterized by comprising:

[0007] Establish a 6-DoF dynamic model of the lander;

[0008] Based on the constraints of the landing process, an optimization problem model is established, resulting in the continuous optimal control problem P1.

[0009] The 6-DoF dynamic model of the lander is made convex, and the continuous optimal control problem P1 is discretized into a sequential convex programming problem P2. The sequential convex programming method is used to solve the sequential convex programming problem P2.

[0010] A further technical solution of the present invention: the 6-DoF dynamic model of the lander is as follows:

[0011]

[0012] Where, r l ,v l ,P l ∈R 3 These represent the lander's position, velocity, and thrust vector in the lunar solid system, respectively, μ m Let ω be the gravitational parameter of the Moon, ω = [0, 0, ω m ] T Let ω be the lunar rotation angular velocity vector in the lunar-fixed system. m Let J be the lunar rotational angular velocity constant, m be the current mass of the lander, and J be the total mass of the lander. a ∈R 3×3 Let M be the lander inertial tensor matrix in this system. a ∈R 3 q is the torque acting on the lander. a←l =[q1,q2,q3,q4] T Let C be the quaternion parameter from the lunar solid system to this system. a←l For q a←l The corresponding rotation matrix, P a =[P,0,0] T I is the thrust vector in this system. sp For fuel specific impulse, g e This refers to the acceleration due to Earth's gravity.

[0013] A further technical solution of the present invention: The sequential convex programming problem P2 is:

[0014]

[0015] Where μ is the slack variable of the virtual control, and λ ν υ is the penalty factor for μ. f Let λ be the slack variable of the terminal constraint. f For υ f The penalty factor is M, F is the discretized dynamic equation matrix, and X is the variable to be optimized. Let m be the terminal state constraint function. dry h is the dry mass parameter of the lander. l R represents the lowest altitude of the lander above the lunar surface. m Let δ be the radius of the moon. x Let δ be the radius of the trust region of the state variable. u Let P be the radius of the trust region of the control variable, η be the radius of the combined trust region of the state variable and the control variable, υ be the virtual control variable, and P be the radius of the trust region of the control variable. l As the lower limit of thrust, P u M is the upper limit of thrust. u The upper limit of the torque modulus is given by the superscript k-1, which indicates the quantity related to the convexity expansion point, and the superscript k indicates the quantity to be solved in the k-th iteration.

[0016] A further technical solution of the present invention: the solution process of the sequential convex programming method includes:

[0017] Given the trust region update parameter β sh ,β gr ,ρ1,ρ2,0≤ρ0<ρ1,β sh ,β gr >1, convergence criterion parameter ε>0, μ tol >0, the number of discrete point intervals N, and the initial guesses of state and control variables. Let k = 1;

[0018] Step 1: With (x (k-1) ,u (k-1) Let be the point of convex expansion. Solving problem P2 yields (x) (k) ,u (k) );

[0019] Step 2: Check if the following termination conditions are met:

[0020] ||x k -x (k-1) ||2<εandμ<μ tol

[0021] If the condition is met, stop the iteration; otherwise, proceed to step 3.

[0022] Step 3: Using (x (k) ,u (k) Calculate the value of ρ using the following formula:

[0023]

[0024] Adjust the trust domain radius according to the trust domain update mechanism;

[0025] If ρ < ρ0, return directly to step 1, where ρ0 is the initial trust region radius; otherwise, let k = k + 1 and return to step 1.

[0026] A further technical solution of the present invention: the trust region update mechanism is specifically as follows:

[0027] If ρ < ρ1, it indicates that the convex model is inaccurate, and the trust region radius should be reduced by letting η = η / β. sh If ρ1≤ρ<ρ2, it indicates that the convex model is more accurate, and the trust region radius remains unchanged; if ρ≥ρ2, it indicates that the convex model is conservative, and the trust region radius should be increased, letting η=β. gr η; ρ1 and ρ2 are given values.

[0028] A computer system is characterized by comprising: one or more processors, and a computer-readable storage medium for storing one or more programs, wherein when the one or more programs are executed by the one or more processors, the one or more processors cause the one or more processors to implement the method described above.

[0029] A computer-readable storage medium is characterized by storing computer-executable instructions, which, when executed, are used to implement the above-described method.

[0030] The beneficial effects of this invention are as follows:

[0031] This invention provides a sequential convex optimization method for the 6-DoF trajectory of the lunar landing braking phase. It considers the lunar rotation and non-constant gravitational field, establishes the attitude-orbit coupling dynamics of the lander, and can effectively solve the optimal trajectory planning problem of the 6-DoF phase of the lunar landing braking phase, thereby improving the accuracy of lunar landing. Attached Figure Description

[0032] The accompanying drawings are for illustrative purposes only and are not intended to limit the invention. Throughout the drawings, the same reference numerals denote the same parts.

[0033] Figure 1 This is a flowchart of the trajectory optimization problem solving method of the present invention.

[0034] Figure 2 This is a three-dimensional schematic diagram of the landing process obtained by the method of this invention.

[0035] Figure 3 x is obtained by the method of this invention. l -Oy l Landing trajectory in the plane.

[0036] Figure 4 This is the distance-to-lunar-surface-height curve obtained by the method of this invention.

[0037] Figure 5 This is the thrust magnitude curve obtained by the method of this invention.

[0038] Figure 6 This is the torque magnitude curve obtained by the method of this invention.

[0039] Figure 7 This is the curve showing the variation of the body angular velocity components obtained by the method of this invention.

[0040] Figure 8 The ontology x obtained by the method of this invention is a Curve showing the change in the angle between the axis and the lunar surface.

[0041] Figure 9 This is the unit quaternion and its modulus variation curve obtained by the method of this invention. Detailed Implementation

[0042] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention. Furthermore, the technical features involved in the various embodiments of this invention described below can be combined with each other as long as they do not conflict with each other.

[0043] This invention provides a sequential convex optimization method for the 6-DoF trajectory of the lunar landing braking phase. To establish a 6-DoF dynamic model of the lander, the method first assumes that the lander relies on the main engine to adjust thrust and is equipped with an independent attitude engine to adjust torque. Unit quaternions are used to describe the attitude motion relationship of the lander relative to the lunar center, which is fixed. Based on this, the lunar rotation and non-constant gravitational field model are considered. In order to quickly solve the braking phase trajectory optimization problem, linearization, relaxation, virtual control and other techniques are combined to convexize the dynamic model and related flight constraints, thereby modeling the braking phase trajectory optimization problem as a sequential convex programming problem. A trust region radius update mechanism is introduced to improve the iterative convergence, so as to solve this complex nonlinear model to obtain the 6-DoF trajectory.

[0044] like Figure 1 As shown, it includes the following steps:

[0045] Step 1: Establish a 6-DoF dynamic model of the lander.

[0046] A 6-DoF dynamic model is established by selecting the lunar center-fixed linkage and the lander's own system, as defined below:

[0047] (1) Monthly fixed system Ox l y l z l The origin is located at the center of the moon, x l The axis points to the prime meridian in the lunar equatorial plane, z l The axis is perpendicular to the lunar equatorial plane and aligns with the direction of the rotational angular velocity vector. l The axis and the two axes form a right-handed rectangular coordinate system.

[0048] (2)This systemox a y a z a The origin is located at the lander's center of mass, x a The axis points in the opposite direction to the main engine nozzle, and the axis is located in the longitudinal plane of the lander and points upwards. a The axis and the two axes form a right-handed right-angled system.

[0049] To simplify the calculations, we assume that the lander's center of mass position and inertial tensor remain constant throughout the flight. First, the translational dynamics model in the lunar-fixed frame is as follows:

[0050]

[0051] In the formula: r l ,v l ,P l ∈R 3 These represent the lander's position, velocity, and thrust vector in the lunar solid system, respectively, μ m Let ω be the gravitational parameter of the Moon, ω = [0, 0, ω m ] T Let ω be the lunar rotation angular velocity vector in the lunar-fixed system. m Let m be the lunar rotational angular velocity constant, and m be the current mass of the lander.

[0052] This invention uses quaternions to describe the attitude relationship between the lunar-fixed system and the local system. The lander attitude dynamics are as follows:

[0053]

[0054] In the formula: ω a ∈R 3 The antisymmetric matrix describes the angular velocity vector between two coordinate systems within this system.

[0055]

[0056] Among them, Ja ∈R 3×3 Let M be the lander inertial tensor matrix in this system. a ∈R 3 q is the torque acting on the lander. a←l =[q1,q2,q3,q4] T Let be the quaternion parameters from the lunar solid system to the local system, then the corresponding rotation matrix is ​​described as follows:

[0057]

[0058] And there are

[0059] In summary, the 6-DoF dynamic model of the lander is as follows:

[0060]

[0061] In the formula, P a =[P,0,0] T I is the thrust vector in this system. sp For fuel specific impulse, g e Let gravitational acceleration be the acceleration due to Earth's gravity. Define state variables. Control quantity

[0062] Step 2: Establish an optimization problem model.

[0063] The landing process needs to satisfy a series of state and control constraints, defined as follows:

[0064] (1) Terminal constraints

[0065] The initial state of the braking phase is given as a known condition. To ensure a smooth transition to the subsequent approach phase, the terminal position, velocity, and quaternions are constrained as follows:

[0066] r l (t f ) = r lf ,v l (t f ) = v lf ,q a←l (t f )=q f ,ω a (t f )=ω f (6)

[0067] In short:

[0068] Ψ(x(t f ),t f )=0 (7)

[0069] (2) State constraints

[0070] The inherent constraints on state quantities include dry mass constraints:

[0071] m(t)≥m dry (8)

[0072] Lander body rotation angular velocity magnitude constraint:

[0073] ||ω a ||≤ω max (9)

[0074] And quaternion modulus constraints:

[0075] ||q a←l ||2=1 (10)

[0076] To avoid hitting terrain obstacles, the braking section height is maintained at h. l above:

[0077] ||r l ||2-R m ≥h l (11)

[0078] In the formula, R m The radius is the radius of the moon.

[0079] (3) Control constraints

[0080] Limit the magnitude of thrust and torque to within permissible limits:

[0081] P l ≤P≤P u (12)

[0082] ||M a ||≤M u (13)

[0083] Taking into account the above dynamic model and related constraints, and using fuel optimization as the performance index, the following optimal control problem can be formed:

[0084]

[0085] Step 3: Sequential convex planning method.

[0086] (1) Model convexity processing

[0087] The dynamic model (5) has strong nonlinear characteristics. For ease of description, it is abbreviated as:

[0088]

[0089] Then perform successive linearization approximation on it:

[0090]

[0091] in,

[0092]

[0093] Linearized expansion point x (k-1) and u (k-1) The known solution obtained in the previous iteration is called the convex expansion point. In the early stages of the sequential solution (when the iteration number k is small), the possibility of solution failure due to the infeasibility of linearization dynamics (15) is relatively high. To avoid this problem and enhance the robustness of the solution process, a virtual control υ∈R is introduced into the linearization dynamics equation (15). 14 :

[0094]

[0095] As the sequence solution progresses, the size of the virtual control should gradually decrease. To achieve this, a slack variable μ is introduced and constraints are applied:

[0096] ||υ||≤μ (18)

[0097] A penalty term for μ is added to the performance metrics. To ensure the accuracy of the linearized dynamics in approximating the original dynamics, and to keep the sequential convex programming problem bounded throughout the solution process, trust region constraints are defined:

[0098] δ x =||xx (k-1) ||2,δ u =||uu (k-1) ||2,δ x +δ u ≤η (19)

[0099] In the formula: η is the trust region radius.

[0100] Furthermore, considering that the terminal constraints of the equations may be too stringent and lead to solution failure, a relaxation process is performed on them:

[0101] Ψ(x(t f ),t f )=υ f (20)

[0102] Add information about ||υ to the performance metrics. f The penalty term for ||. State constraints (10) and (11) are non-convex constraints. Constraint (10) is relaxed as follows:

[0103] ||q a← l||2≤1 (21)

[0104] The safety height constraint (11) is approximated using a linearization method:

[0105]

[0106] (2) Discretization

[0107] To transform the continuous optimal control problem P1 into a finite-dimensional optimization problem, in the finite time domain [0, t f The area is divided into N equal intervals, forming a total of N+1 discrete points, with a discrete point interval Δt = t. f / N, the corresponding discrete state variables, control variables, and virtual control variables are denoted as: x i =x(t) i ), u i =u(t) i ), υ i =υ(t) i ), i = 0, 1, ..., N. According to Euler's method, the dynamic equation (17) can be discretized as:

[0108] x i =x i-1 +Δt(A i-1 (k-1) x i-1 +B i-1 (k-1) u i-1 +d i-1 (k-1) +υ i-1 ) (twenty three)

[0109] The optimization variables include state variables {x0, x1, ..., x}. N}、Control variables {u0,u1,...,u N} and virtual control quantities {υ0,υ1,...,υ N} can be uniformly listed as:

[0110]

[0111] After simplification, the constraint (23) at each discrete point can be expressed as:

[0112] MX = F (25)

[0113] In the formula:

[0114]

[0115] Finally, the sequential convex programming (sub-problem) for the k-th (k=1,2,...) iteration requires a solution as follows:

[0116]

[0117] In the formula, the superscript k-1 corresponds to the convex expansion point related quantity, the superscript k represents the quantity to be solved in the k-th iteration, and λ ν ∈R + and λ f ∈R + These are the penalty factors corresponding to equations (18) and (20), respectively.

[0118] (3) Trust region radius update

[0119] Define the linearization accuracy of the dynamics as δ = [δ0,...,δ N ] T :

[0120] δ i =x i+1 -(x i +ΔtF(t i ,x i ,u i (26)

[0121] That is, δ i The state variable x obtained at a certain discrete point i The result of positive integration over a discrete interval and the solved state variable x i+1 The difference. If this difference is greater than a certain threshold, it indicates that the obtained solution is infeasible for the original nonconvex dynamics. The nonlinear augmented cost function is further defined as:

[0122] J λ (x,u)=-m(t f )+Ψ(x(t f ),t f )+||δ|| ∞ (27)

[0123] And the linear augmentation cost function is:

[0124] L λ (x,u)=-m(t f )+||M|| ∞ (28)

[0125] In the formula, M = [υ0,...,υ N ] T The convexity accuracy of the dynamics can then be measured by the following formula:

[0126]

[0127] In the formula, the superscript k-1 represents the solution of the previous iteration, k represents the solution of the current iteration, the denominator represents the change in the predicted performance index, and the numerator represents the change in the actual performance index. Based on the value of ρ, the following trust region radius update mechanism can be formed: Setting the parameter β...sh ,β gr >1 and 0 < ρ1 < ρ2 ≤ 1. If ρ < ρ1, it indicates that the convex model is inaccurate, and the trust region radius should be reduced. Let η = η / β sh If ρ1≤ρ<ρ2, it indicates that the convex model is more accurate, and the trust region radius remains unchanged; if ρ≥ρ2, it indicates that the convex model is conservative, and the trust region radius should be increased, letting η=β. g The values ​​of rη, ρ1, and ρ2 can be selected based on the solution requirements; typically, ρ2 is taken to be around 1.

[0128] (4) Solution process

[0129] Given the trust region update parameter β sh ,β g r, ρ1, ρ2, 0 ≤ ρ0 < ρ1, and convergence criterion parameter ε > 0, μ tol >0, the number of discrete point intervals N, and the initial guesses of state and control variables. Let k = 1. The iterative solution process for the sequential convex optimization of problem P2 is as follows:

[0130] Step 1: With (x (k-1) ,u (k-1) Let be the point of convex expansion. Solving problem P2 yields (x) (k) ,u (k) ).

[0131] Step 2: Check if the following termination conditions are met:

[0132] ||x k -x (k-1) ||2<εandμ<μ ol

[0133] If the condition is met, stop the iteration; otherwise, proceed to step 3.

[0134] Step 3: Using (x (k) ,u (k) Calculate the value of ρ according to (29) and adjust the radius of the trust region according to the trust region update mechanism. If ρ < ρ0, return directly to step 1; otherwise, let k = k + 1 and return to step 1.

[0135] The following simulation scenario is used, assuming the flight mission starts deceleration from the lunar perigee of a 200×15km elliptical orbit with an inclination of 0°, the target landing point is located at the lunar latitude and longitude, and the desired terminal altitude is 2.5km. Specific simulation parameters are shown in Table 1. The specified terminal time t... f =500s, the number of discrete points is 301, and the initial guess at each discrete point is obtained by linear interpolation using the specified initial and terminal states.

[0136] Table 1

[0137]

[0138]

[0139] Simulation results are as follows Figures 2-9 As shown, the results indicate a clear bang-bang characteristic in the thrust magnitude profile, with both thrust and torque magnitudes within the constrained amplitude range. During landing, the lander primarily orbits the x-axis. a The shaft rotates to change the direction of the thrust. Figure 8 Display of the braking phase terminal lander body x a The angle between the axis and the lunar surface reaches the constrained 45° angle. Figure 9 The changes in the unit quaternion components and their moduli are shown, indicating that the convex constraint (21) plays an active role throughout the landing process.

[0140] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any person skilled in the art can easily conceive of various equivalent modifications or substitutions within the scope of the technology disclosed in the present invention, and such modifications or substitutions should all be covered within the scope of protection of the present invention.

Claims

1. A sequential convex optimization method for the 6-DoF trajectory of the lunar landing braking phase, characterized in that, include: Establish a 6-DoF dynamic model of the lander; Based on the constraints of the landing process, an optimization problem model is established, resulting in the continuous optimal control problem P1. The 6-DoF dynamics model of the lander is made convex, and the continuous optimal control problem P1 is discretized into a sequential convex programming problem P2; the sequential convex programming method is used to solve the sequential convex programming problem P2. The sequential convex programming problem P2 is: in, For slack variables in virtual control, for The penalty factor For the slack variables of the terminal constraints, for The penalty factor and The discretized dynamic equation matrix, For the variable to be optimized, The terminal state constraint function, For the dry mass parameters of the lander, This represents the lowest altitude of the lander above the lunar surface. The radius of the moon, Let the radius of the trust region of the state variable be denoted as . To control the radius of the trust region, The radius of the combined trust region for state and control variables. For virtual control variables, As the lower limit of thrust, This is the upper limit of thrust. The upper limit of the torque modulus, indicated by the superscript. k -1 indicates quantities related to the convexity expansion point, superscript k Indicates the first k The quantity to be solved in the next iteration.

2. The sequential convex optimization method for the 6-DoF trajectory of the lunar landing braking segment according to claim 1, characterized in that, The 6-DoF dynamics model of the lander is as follows: in, These represent the lander's position, velocity, and thrust vector in the lunar-fixed system. For lunar gravitational parameters, The vector of the Moon's rotational angular velocity in the lunar-fixed system. Let ω be the constant of the Moon's rotational angular velocity. m For the current mass of the lander, Let be the lander inertial tensor matrix in this system. The torque acting on the lander, The quaternion parameters from the lunar solid system to this system are: for The corresponding rotation matrix, This is the thrust vector in this system. For fuel specific impulse, This refers to the acceleration due to Earth's gravity.

3. The sequential convex optimization method for the 6-DoF trajectory of the lunar landing braking segment according to claim 1, characterized in that, The solution process of the sequential convex programming method includes: Given trust region update parameters , Convergence criterion parameters Number of discrete point intervals N Initial guesses of state variables and control variables , ,remember k = 1; Step 1: With To find the convex expansion point, solve problem P2 to obtain... ; Step 2: Check if the following termination conditions are met: If the condition is met, stop the iteration; otherwise, proceed to step 3. Step 3: Utilize Calculate according to the following formula value: Adjust the trust domain radius according to the trust domain update mechanism; like Return directly to step 1. Let be the initial trust region radius; otherwise, let k = k +1, return to step 1.

4. The sequential convex optimization method for the 6-DoF trajectory of the lunar landing braking segment according to claim 1, characterized in that, The trust domain update mechanism is specifically as follows: like This indicates that the convexized model is inaccurate, and the trust region radius should be reduced. ;like This indicates that the convexized model is more accurate, and the trust region radius remains unchanged; if This indicates that the convexized model is conservative, and the trust region radius should be increased, letting... ; and For a given value.

5. A computer system, characterized in that... include: One or more processors, a computer-readable storage medium for storing one or more programs, wherein, when the one or more programs are executed by the one or more processors, the one or more processors cause the one or more processors to perform the method of any one of claims 1-4.

6. A computer-readable storage medium, characterized in that... The device stores computer-executable instructions, which, when executed, are used to implement the method described in any one of claims 1-4.

Images