Aircraft path planning method for space non-cooperative target

Abstract

The invention relates to the technical field of aircraft path planning, and particularly provides an aircraft path planning method for a space non-cooperative target. The method comprises the following steps: constructing an aircraft dynamics model and constraint conditions, and obtaining a non-convex optimal control problem; performing discretization processing on a non-convex optimal control problem through a non-equal-interval pseudo-spectrum discrete optimization model, and converting the non-convex optimal control problem into a fixed-time nonlinear programming problem about discrete nodes; according to the method, a fixed time nonlinear programming problem is convexed to form a convex sub problem, a local optimal solution of a non-convex optimal control problem is obtained through sequence iteration, under the multi-state constraint, path planning of an aircraft approaching a space non-cooperative target is achieved, the problem of no feasible solution caused by non-accurate linearization is solved, and the method is suitable for being applied to a non-cooperative target. And the convergence rate of the algorithm is accelerated.

Description

Technical Field

[0001] The present invention relates to the technical field of aircraft path planning, and in particular to an aircraft path planning method for a non-cooperative target in space. Background Art

[0002] Non-cooperative space targets are those in erratic motion that lack the ability to actively exchange information, sensitive markers for identification and positioning, or even attitude and orbit control capabilities. These typically include launch vehicle upper stages, defunct satellites, and enemy aircraft. Whether performing on-orbit space servicing tasks such as module assembly, refueling, and troubleshooting, or conducting surveillance, reconnaissance, and close-in capture in space situational awareness and offensive and defensive confrontation, all require pre-planned maneuvering trajectories to guide the vehicle toward its target.

[0003] The process of approaching a non-cooperative target in space is subject to numerous nonlinear path constraints, such as collision avoidance constraints, field-of-view pointing constraints, and thruster plume avoidance constraints. The main methods for designing vehicle maneuvering trajectories under multiple path constraints include the Gaussian pseudospectral method (GPOPS), artificial potential field methods, and convex programming. GPOPS is a nonlinear programming algorithm developed based on pseudospectral discretization and sequential quadratic programming (SQP). When both path constraints and dynamic constraints exist, the first-order optimality condition (KKT) of the nonlinear programming problem obtained through Gaussian pseudospectral discretization is equivalent to the discretized form of the first-order optimality condition of the original optimal control problem. Considering the task of approaching a tumbling non-cooperative target, GPOPS was used to solve the collision-free translational trajectory planning problem under time and energy optimization. GPOPS was also used to achieve simultaneous planning of the multi-stage flight trajectory of a hypersonic vehicle during ascent, cruise, and dive. Existing techniques have proposed the adaptive HP-pseudospectral method, also known as the local collocation method. This method divides the entire problem time domain into several segments (grids) based on an error criterion and locally adjusts the order of the interpolation polynomial within each segment to improve algorithm performance with less computational effort. Although its performance has been verified in practical applications, the Gaussian pseudospectral method still has shortcomings. For problems with non-smooth control inputs, the state and control trajectory represented by the polynomial are prone to chattering at discrete nodes.

[0004] The computational complexity of numerical optimization depends on the convex or non-convex nature of the original optimal control problem. With the deepening of research in interior-point optimization theory, convex programming algorithms such as second-order cone programming and semidefinite programming have been widely used in trajectory planning for rocket descent, asteroid landing, and approaching non-cooperative targets in space. Second-order cone programming was used to solve the optimal guidance problem for Mars landing fuel, and a lossless relaxation method was proposed to transform non-convex control amplitude constraints into convex constraints. Similarly, constraint relaxation or variable substitution methods have been used to achieve accurate lossless convexification of various non-convex state constraints. However, many optimal control problems involving nonlinear dynamics and path constraints are difficult to transform losslessly. To continue using convex programming to solve non-convex optimal control problems, numerous researchers have proposed and developed sequential convex programming (SCP). Similar to the SQP algorithm, which recursively approximates the original problem using a series of quadratic programming subproblems, the SCP algorithm approaches the local optimal solution of the original non-convex problem by sequentially solving approximate convex programming subproblems. Considering the trajectory planning problem for autonomous rendezvous and docking of aircraft under collision avoidance, docking corridor, and thrust pointing constraints, the non-convex constraints are first convexified using a linearization method. The fuel-optimal maneuvering trajectory is then designed using the SCP method, and the convergence of the sequential optimization algorithm is analyzed. Prior art proposes rotating hyperplane and dual hyperplane methods to convexify obstacle avoidance constraints. Based on second-order cone programming, a solution is provided for optimal control problems containing concave inequality constraints and nonlinear terminal equality constraints. Concave inequality constraints are approximated through sequential linearization, while nonlinear equality constraints are not only linearized but also corrected using second-order expansion terms during the sequential solution process to enhance the robustness of the algorithm's convergence. Currently, linearization and convexification methods can cause artificial infeasibility problems. That is, when the linearized model built around the reference trajectory differs significantly from the original nonlinear model, the convex subproblem may have no feasible solution. Summary of the Invention

[0005] In view of this, the present invention provides an aircraft path planning method for non-cooperative targets in space, which is used to realize path planning for aircraft approaching non-cooperative targets in space under multi-state constraints, alleviate the problem of no feasible solution caused by inexact linearization, and accelerate the convergence rate of the algorithm.

[0006] In a first aspect, the present invention provides a method for aircraft path planning for a non-cooperative space target, the method comprising:

[0007] Step 1: Construct the aircraft dynamics model and constraints to obtain the non-convex optimal control problem;

[0008] Step 2: Through the non-uniformly spaced pseudo-spectral discrete optimization model, the non-convex optimal control problem is discretized and converted into a fixed-time nonlinear programming problem with discrete nodes;

[0009] Step 3: Convexify the fixed-time nonlinear programming problem to form a convex subproblem, and obtain the local optimal solution of the non-convex optimal control problem through sequential iteration.

[0010] Optionally, building the aircraft dynamics model in step 1 includes:

[0011] To describe the relative motion of aircraft, we use represents the J2000 Earth Inertial Coordinate System, Represent the aircraft and target's local systems respectively, The orbital coordinate system of the target is represented by the following equations: The translational and rotational dynamic equations of the spacecraft approaching the non-cooperative target in space are expressed as follows:

[0012]

[0013]

[0014]

[0015] in, m C Indicates the mass of the aircraft; n ref for Relative to The angular velocity vector of the system is The weight in the system; express The translational motion state of the aircraft relative to the target, To express in The control quantity under the system, is the attitude control torque of the aircraft, is the moment of inertia of the aircraft, For aircraft relative The attitude of the system is corrected by the Rodrigues parameters MRPs, For the aircraft relative to The angular velocity of the system, For the aircraft relative to The relationship between the angular velocity and the speed is expressed as:

[0016]

[0017] in, Indicates The orbital angular velocity in the system, For Transfer to The coordinate transformation matrix of the system; using Indicates the rotational motion state of the aircraft.

[0018] Optionally, the constraints constructed in step 1 include:

[0019] The target spacecraft body is replaced by a sphere, and the solar panels installed on the target body are represented by an ellipsoid. The collision avoidance constraint is established, and its expression is:

[0020]

[0021]

[0022] in, r ball represents the radius of the collision sphere, is a constant matrix, C p =[I3,0 3×9 ], Represents the collision ellipsoid feature matrix, r x , r y , r z is the semi-major axis of the ellipsoid, Indicates from Tie to The coordinate transformation matrix of the system; Indicates from Tie to The coordinate transformation matrix of the system is determined according to the terminal attitude of the target aircraft;

[0023] Create a path constraint, the expression is:

[0024]

[0025] Among them, β fov is the half cone angle of the navigation equipment field of view, are the direction of the device optical axis and its position vector in this system;

[0026] The control output of the aircraft has an amplitude constraint, which is expressed as:

[0027]

[0028] Among them, F max and T max denote the magnitudes of the control resultant force and resultant moment, respectively;

[0029] Taking the optimal energy consumption as the indicator, the non-convex optimal control problem of autonomous approach is established, and its expression is:

[0030]

[0031] Among them, t0, t f Denote the initial and terminal moments of the maneuver, x0, x frepresents the corresponding nominal state at the initial and final moments; st(1)-(3), (5)-(8) represent the constraints imposed by formulas (1)-(3) and (5)-(8).

[0032] Optionally, step 2 includes:

[0033] The non-uniformly spaced pseudo-spectral discretization method is used, that is, the Gauss-Legendre points of the state variables and control variables are discretized. The Gauss-Legendre points are distributed on [-1, 1]. The time interval [t0, t f ] is converted to the interval [-1, 1], and a variable is introduced for this purpose, and its expression is:

[0034]

[0035] Among them, t0, t f denotes the initial and terminal moments of the maneuver;

[0036] The continuous optimal control problem after time domain transformation is expressed as:

[0037]

[0038]

[0039] Among them, L(u τ ) represents the index function to be optimized; f(x τ ,u τ ,τ;t0,t f ) represents the right side of differential equations (1) to (3); function C(x τ ,u τ ,τ;t0,t f ) represents various path constraints; is the terminal constraint, Represent the control force and torque respectively; x τ represents the state variable; τ represents the time factor; Indicates the initial state; Indicates terminal status;

[0040] Gauss-Legendre points are N-th degree Legendre polynomials P N The zero point of (τ) is expressed as:

[0041]

[0042] The Lagrange interpolation method is used to discretize the state variables in the continuous Bolza problem at the first N+1 LG nodes. The approximate expression of the state variables is:

[0043]

[0044] Among them, K represents the order of interpolation fitting, that is, how many points are used to fit the state curve; X i represents the control point, i.e. the state variable of the interpolation point; L i (τ) represents the Lagrange interpolation basis function;

[0045] The Lagrange interpolation of the control variables is performed at the discrete nodes, and the approximate expression of the discretized control variables is:

[0046]

[0047] Among them, U i represents the control points of the control variable fitting;

[0048] For the differential equation of the state variable, take the derivative of equation (12) and get:

[0049]

[0050] Among them, D ki equal X(τ i ) represents the column vector composed of control points;

[0051] Based on the differential approximation of the state variables mentioned above, the dynamic differential equation is configured at the discrete nodes, and its expression is:

[0052]

[0053] Among them, f(X(τ k ), U(τ k ), τ; t0, t f ) represents the system differential equation, and its corresponding state variables are the control points of the state variables and the control points of the control variables; X(τ k ) represents τ k The control point of the state variable at the moment; U(τ k ) represents τ k The control point of the moment control variable; τ k Represents a time value;

[0054] Using Gauss quadrature formula, we can quadrature the dynamic equation and get:

[0055]

[0056] Among them, ω i represents the fitting coefficient;

[0057] Then the terminal state of the system at each discrete point is approximately:

[0058]

[0059] Among them, ω k Represents the Gauss integral coefficient, and its expression is:

[0060]

[0061] in, Indicates P N In τ k The derivative at ;

[0062] The inequality path constraints of state variables and control variables are approximated at each discrete point as follows:

[0063] C(X k , U k , τ k ;t0,t f )≤0(k=1,....N)(19)

[0064] The terminal constraints at each discrete point are approximated as:

[0065] φ(X(τ0),t0,X(τ f ), t f )=0(20)

[0066] The objective function is approximately:

[0067]

[0068] Among them, g(X k , U k ) represents the indicator function;

[0069] In summary, the non-convex optimal control problem is discretized and converted into a fixed-time nonlinear programming problem about discrete nodes, which is expressed as follows:

[0070]

[0071] Among them, st(15), (19), (20) means being constrained by formulas (15), (19), (20).

[0072] Optionally, step 3 includes:

[0073] Given a reference trajectory Will The nonlinear term f(X k , U k , τ k ;t0,t f ) is expanded by Taylor around it, and we get:

[0074]

[0075] in,

[0076] matrix is the Jacobian matrix of the nonlinear term about the discrete nodes, and its expression is:

[0077]

[0078] Bring it in In, we get:

[0079]

[0080] The nonlinear path constraint function C(X k , U k , τ k ;t0,t f ) in the reference trajectory Expand near the linearized inequality constraint, and its expression is:

[0081]

[0082] in, is the gradient of the constraint function at the reference trajectory;

[0083] A soft trust region constraint is additionally introduced into the optimization problem to solve the problem of no feasible solution, that is, an inequality constraint is imposed, which is expressed as:

[0084]

[0085] Among them, η k Represents the trust region radius, which is the variable to be optimized. In order to promote the continuous shrinkage of the trust region radius during the iteration process, a related penalty term is added to the indicator function, and its expression is:

[0086]

[0087] In summary, the expression of the convex subproblem after convexification is:

[0088] min J+J tr

[0089] st(20),(24),(25),(26)(27)

[0090] Use z to uniformly represent the optimization variable z=[X T , U T ] T , define p m=z m+1 -z m Indicates the search direction, z m The equality and inequality constraint violations and the index function are expressed by the exact penalty function, which is expressed as follows:

[0091]

[0092] Among them, μ = [μ1, μ2] is the penalty factor, f E，k (z m ) represents a nonlinear equality constraint, which is expressed as:

[0093]

[0094] Define the step factor α∈[0,1], the line search method starts from α=1, and shrinks in sequence α′=cα, c∈(0,1), along the direction p m Search for new iteration points Until it satisfies the following non-exact one-dimensional search Armijo condition, which is expressed as:

[0095] φ1(z m +αp m ;μ)≤φ1(z m ;μ)+c1αD(φ1(z m ;μ);p m );

[0096] Among them, 0<c1<1, D(φ1(z m ;μ);p m ) represents the function φ1(z m ; μ) in direction p m The directional derivative on is expressed as:

[0097]

[0098] In summary, when solving the non-convex optimal control problem, the initial reference trajectory is first obtained by linear interpolation. Then, in each iteration, the iteration step size and penalty factor are updated based on line search, and the coefficient matrix related to the linearized dynamics and convex path constraints is calculated. Then, the newly established convex subproblem, formula (27), is solved based on convex optimization. When the difference between adjacent iterative trajectories meets the stopping condition, the sequential optimization process is stopped. The expression of the stopping condition is:

[0099] ΔX=||X m -X m-1 || ∞≤ε. In a second aspect, an embodiment of the present invention provides a computer-readable storage medium, wherein the computer-readable storage medium includes a stored program, wherein when the program is run, the device where the computer-readable storage medium is located is controlled to execute the aircraft path planning method for a non-cooperative space target in the first aspect or any possible implementation of the first aspect.

[0100] In a third aspect, an embodiment of the present invention provides an electronic device comprising: one or more processors; a memory; and one or more computer programs, wherein the one or more computer programs are stored in the memory, and the one or more computer programs include instructions that, when executed by the device, enable the device to execute the aircraft path planning method for a non-cooperative target in space in the first aspect or any possible implementation of the first aspect.

[0101] In the technical solution provided by the present invention, the method includes constructing an aircraft dynamics model and constraints to obtain a non-convex optimal control problem; through a non-uniformly spaced pseudo-spectral discrete optimization model, the non-convex optimal control problem is discretized and converted into a fixed-time nonlinear programming problem about discrete nodes; the fixed-time nonlinear programming problem is convexified to form a convex subproblem, and the local optimal solution of the non-convex optimal control problem is obtained through sequence iteration. Under multi-state constraints, this method realizes the path planning of the aircraft approaching a non-cooperative target in space, alleviates the problem of no feasible solution caused by imprecise linearization, and accelerates the convergence rate of the algorithm. BRIEF DESCRIPTION OF THE DRAWINGS

[0102] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0103] Figure 1 A flow chart of an aircraft path planning method provided by an embodiment of the present invention;

[0104] Figure 2 A schematic diagram of a three-dimensional maneuvering trajectory provided by an embodiment of the present invention;

[0105] Figure 3 A distance curve diagram for collision avoidance provided by an embodiment of the present invention;

[0106] Figure 4 A graph showing the angle change of the field of view provided by an embodiment of the present invention;

[0107] Figure 5 A schematic diagram of a trust region convergence process provided by an embodiment of the present invention;

[0108] Figure 6 A schematic diagram of an electronic device provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0109] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts shall fall within the scope of protection of the present invention.

[0110] It should be understood that the embodiments described are only a portion of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by persons of ordinary skill in the art without creative work are within the scope of protection of the present invention.

[0111] The terms used in the embodiments of the present invention are only for the purpose of describing specific embodiments and are not intended to limit the present invention. The singular forms "a", "the" and "the" used in the embodiments of the present invention are also intended to include plural forms, unless the context clearly indicates other meanings.

[0112] It should be understood that the term "and / or" as used herein simply describes a relationship between associated objects, indicating that three possible relationships exist. For example, "A and / or B" can represent: A alone, A and B together, or B alone. Furthermore, the character " / " in this document generally indicates an "or" relationship between the associated objects.

[0113] The word "if," as used herein, may be interpreted as "at the time of" or "when" or "in response to determining" or "in response to detecting," depending on the context. Similarly, the phrases "if it is determined" or "if (stated condition or event) is detected" may be interpreted as "when it is determined" or "in response to the determination" or "when detecting (stated condition or event)" or "in response to detecting (stated condition or event)," depending on the context.

[0114] Figure 1 A flow chart of an aircraft path planning method provided by an embodiment of the present invention is shown in FIG. Figure 1 As shown, the method includes:

[0115] Step 1: Construct the aircraft dynamics model and constraints to obtain the non-convex optimal control problem;

[0116] In the embodiment of the present invention, constructing the aircraft dynamics model in step 1 includes:

[0117] To describe the relative motion of aircraft, we use represents the J2000 Earth Inertial Coordinate System, Represent the aircraft and target's local systems respectively, The orbital coordinate system of the target is represented by the following equations: The translational and rotational dynamic equations of the spacecraft approaching the non-cooperative target in space are expressed as follows:

[0118]

[0119]

[0120]

[0121] in, m C Indicates the mass of the aircraft; n ref for Relative to The angular velocity vector of the system is The weight in the system; express The translational motion state of the aircraft relative to the target, To express The control quantity under the system, is the attitude control torque of the aircraft, is the moment of inertia of the aircraft, For aircraft relative The attitude of the system is corrected by the Rodrigues parameters MRPs, For the aircraft relative to The angular velocity of the system, For the aircraft relative to The relationship between the angular velocity and the speed is expressed as:

[0122]

[0123] in, Indicates The orbital angular velocity in the system, For Transfer to The coordinate transformation matrix of the system; using Indicates the rotational motion state of the aircraft.

[0124] In the embodiment of the present invention, constructing the constraint conditions in step 1 includes:

[0125] During the approach process, the spacecraft should avoid collision with the target spacecraft and its payload. The target spacecraft body is replaced by a sphere, and the solar panels installed on the target body are represented by an ellipsoid. The collision avoidance constraint is established, and its expression is:

[0126]

[0127]

[0128] in, r ball represents the radius of the collision sphere, is a constant matrix, C p =[I3,0 3×9 ], Represents the collision ellipsoid feature matrix, r x , r y , r z is the semi-major axis of the ellipsoid, Indicates from Tie to The coordinate transformation matrix of the system. It can be seen from this matrix that the ellipsoid collision avoidance constraint couples the target's posture information; Indicates from Tie to The coordinate transformation matrix of the system is determined according to the terminal attitude of the target aircraft;

[0129] Since non-cooperative targets cannot actively provide their own status information, the spacecraft needs to rely on its own navigation equipment to measure its relative status online. Therefore, the spacecraft's measurement equipment should always point to the target spacecraft. Path constraints are established, and their expression is:

[0130]

[0131] Among them, β fov is the half cone angle of the navigation equipment field of view, are the direction of the device optical axis and its position vector in this system;

[0132] The control output of the aircraft has an amplitude constraint, which is expressed as:

[0133]

[0134] Among them, F max and T max denote the magnitudes of the control resultant force and resultant moment, respectively;

[0135] Taking the optimal energy consumption as the indicator, the non-convex optimal control problem of autonomous approach is established, and its expression is:

[0136]

[0137] Among them, t0, t f Denote the initial and terminal moments of the maneuver, x0, x f represents the corresponding nominal state at the initial and final moments; st(1)-(3), (5)-(8) represent the constraints imposed by formulas (1)-(3) and (5)-(8).

[0138] Step 2: Through the non-uniformly spaced pseudo-spectral discrete optimization model, the non-convex optimal control problem is discretized and converted into a fixed-time nonlinear programming problem with discrete nodes;

[0139] In this embodiment of the present invention, step 2 includes:

[0140] For the optimal control problem containing nonlinear equality and inequality constraints, that is, formula (9), a direct method is usually used to solve it, that is, by discretizing in the time domain, the continuous optimization problem is transformed into a parameter programming problem. In order to reduce the number of optimization variables and avoid the Runge phenomenon of equidistant dispersion, a non-equidistant pseudo-spectral discretization method is used, that is, the Gauss-Legendre points of the state variables and control variables are discretized. The Gauss-Legendre points are distributed on [-1, 1]. The time interval [t0, t f ] is converted to the interval [-1, 1], and a variable is introduced for this purpose, and its expression is:

[0141]

[0142] Among them, t0, t f denotes the initial and terminal moments of the maneuver;

[0143] The continuous optimal control problem after time domain transformation is expressed as:

[0144]

[0145] Among them, L(u τ ) represents the index function to be optimized; f(x τ ,u τ ,τ;t0,t f ) represents the right side of differential equations (1) to (3); function C(x τ ,u τ ,τ;t0,t f ) represents various path constraints; is the terminal constraint, Represent the control force and torque respectively; x τ represents the state variable; τ represents the time factor; Indicates the initial state; Indicates terminal status;

[0146] Gauss-Legendre points are N-th degree Legendre polynomials P N The zero point of (τ) is expressed as:

[0147]

[0148] The Lagrange interpolation method is used to discretize the state variables in the continuous Bolza problem at the first N+1 LG nodes. The approximate expression of the state variables is:

[0149]

[0150] Among them, K represents the order of interpolation fitting, that is, how many points are used to fit the state curve; X i represents the control point, i.e. the state variable of the interpolation point; L i (τ) represents the Lagrange interpolation basis function;

[0151] The Lagrange interpolation of the control variables is performed at the discrete nodes, and the approximate expression of the discretized control variables is:

[0152]

[0153] Among them, U i represents the control points of the control variable fitting;

[0154] For the differential equation of the state variable, take the derivative of equation (12) and get:

[0155]

[0156] Among them, D ki equal X(τ i ) represents the column vector composed of control points;

[0157] Based on the differential approximation of the state variables mentioned above, the dynamic differential equation is configured at the discrete nodes, and its expression is:

[0158]

[0159] Among them, f(X(τ k ), U(τ k ), τ k ;t0,t f ) represents the system differential equation, and its corresponding state variables are the control points of the state variables and the control points of the control variables; X(τ k ) represents τ k The control point of the state variable at the moment; U(τk ) represents τ k The control point of the moment control variable; τ k Represents a time value;

[0160] Using Gauss quadrature formula, we can quadrature the dynamic equation and get:

[0161]

[0162] Among them, ω i represents the fitting coefficient;

[0163] Then the terminal state of the system at each discrete point is approximately:

[0164]

[0165] Among them, ω k Represents the Gauss integral coefficient, and its expression is:

[0166]

[0167] in, Indicates P N In τ k The derivative at ;

[0168] The inequality path constraints of state variables and control variables are approximated at each discrete point as follows:

[0169] C(X k , U k , τ k ;t0,t f )≤0(k=1,....N)(19)

[0170] The terminal constraints at each discrete point are approximated as:

[0171] φ(X(τ0),t0,X(τ f ), t f )=0 (20)

[0172] The objective function is approximately:

[0173]

[0174] Among them, g(X k , U k ) represents the indicator function;

[0175] In summary, the non-convex optimal control problem is discretized and converted into a fixed-time nonlinear programming problem about discrete nodes, which is expressed as follows:

[0176]

[0177] Among them, st(15), (19), (20) means being constrained by formulas (15), (19), (20).

[0178] Step 3: Convexify the fixed-time nonlinear programming problem to form a convex subproblem, and obtain the local optimal solution of the non-convex optimal control problem through sequential iteration.

[0179] Before using convex optimization to establish a solution algorithm, it is necessary to convexify the motion equations and nonlinear path constraints of the fixed-time nonlinear programming problem, i.e., the equations of motion and the nonlinear path constraints in formula (25), to form a convex subproblem; then, the local optimal solution of the non-convex optimal control problem is approximated by sequential iteration.

[0180] In this embodiment of the present invention, step 3 includes:

[0181] Given a reference trajectory Will The nonlinear term f(X k , U k , τ k ;t0,t f ) is expanded by Taylor around it, and we get:

[0182]

[0183] in,

[0184] matrix is the Jacobian matrix of the nonlinear term about the discrete nodes, and its expression is:

[0185]

[0186] Bring it in In, we get:

[0187]

[0188] The nonlinear path constraint function C(X k , U k , τ k ;t0,t f ) in the reference trajectory Expand near the linearized inequality constraint, and its expression is:

[0189]

[0190] in, is the gradient of the constraint function at the reference trajectory;

[0191] Since the terminal moment equality constraint is a convex constraint, convexification is not required. In addition, the linearization operation in the convexification process may cause the solution set to be unbounded, resulting in no feasible solution to the optimization problem. In order to solve the problem of no feasible solution, a soft trust region constraint is additionally introduced into the optimization problem, that is, an inequality constraint is imposed, which is expressed as:

[0192]

[0193] Among them, η k Represents the trust region radius, which is the variable to be optimized. In order to promote the continuous shrinkage of the trust region radius during the iteration process, a related penalty term is added to the indicator function, and its expression is:

[0194]

[0195] In summary, the expression of the convex subproblem after convexification is:

[0196]

[0197] It should be noted that traditional sequence convexification algorithms usually use the linear interpolation trajectory between the initial and terminal states for initialization, and directly use the optimal trajectory X of the current iteration during the sequence optimization process. m , U m To update the reference trajectory for the next iteration Thus we can get the optimal solution X of the corresponding convex subproblem m+1 , U m+1 , m is the number of iterations. m , U m Compared with the subsequent optimal trajectory X m+1 , U m+1 Although this can further reduce the optimization index, it does not guarantee that the violation of the original non-convex path constraint will continue to decrease. There may still be another trajectory between two adjacent optimization trajectories that minimizes the sum of the optimization index and the path constraint violation. This application uses an inexact line search to obtain this trajectory.

[0198] Use z to uniformly represent the optimization variable z=[X T , U T ] T , define p m =z m+1 -z m Indicates the search direction, z m The equality and inequality constraint violations and the index function are expressed by the exact penalty function, which is expressed as follows:

[0199]

[0200] Among them, μ = [μ1, μ2] is the penalty factor, fE，k (z m ) represents a nonlinear equality constraint, which is expressed as:

[0201]

[0202] Define the step factor α∈[0,1], the line search method starts from α=1, and shrinks in sequence α′=cα, c∈(0,1), along the direction p m Search for new iteration points Until it satisfies the following non-exact one-dimensional search Armijo condition, which is expressed as:

[0203] φ1(z m +αp m ;μ)≤φ1(z m ;μ)+c1αD(φ1(z m ;μ);p m );

[0204] Among them, 0<c1<1, D(φ1(z m ;μ);p m ) represents the function φ1(z m ; μ) in direction p m The directional derivative on is expressed as:

[0205]

[0206] In summary, when solving the non-convex optimal control problem, the initial reference trajectory is first obtained by linear interpolation. Then, in each iteration, the iteration step size and penalty factor are updated based on line search, and the coefficient matrix related to the linearized dynamics and convex path constraints is calculated. Then, the newly established convex subproblem, formula (27), is solved based on convex optimization. When the difference between adjacent iterative trajectories meets the stopping condition, the sequential optimization process is stopped. The expression of the stopping condition is:

[0207] ΔX=||X m -X m-1 || ∞ ≤ε.

[0208] Simulation analysis of the present invention:

[0209] The effectiveness of the proposed method is illustrated by a simulation example. The number of optimized discrete nodes is 80, the maneuvering time is set to 150s, the algorithm convergence threshold is ε=1e-5, and the other simulation parameters are shown in Table 1.

[0210] Table 1 Simulation parameters

[0211] Parameter (unit) Value Initial position (m) <![CDATA[[10,-20,2] T ]]> Target spacecraft docking port (m) <![CDATA[[0.00;5.00;2.00] T ]]> Initial speed (m / s) <![CDATA[[0,0,0] T ]]> Terminal velocity (m / s) <![CDATA[[0,0,0] T ]]> Initial posture (-) <![CDATA[[0.34,0.41,0.37] T ]]> Field of view constraint (°) 25 Semi-major axis of ellipsoid obstacle zone (m) <![CDATA[[3.00,22.00,3.00] T ]]> Radius of obstacle sphere (m) 12.50 Initial angular velocity (rad / s) <![CDATA[[0,0,0] T ]]> Terminal angular velocity (rad / s) <![CDATA[[0,0,0] T ]]>

[0212] In the embodiment of the present invention, Figure 2 As shown in the figure, the aircraft starts to maneuver from the initial position, bypasses the circular and ellipsoidal obstacles and arrives at the terminal position. The distance curve between the aircraft and the obstacles is shown in the figure. Figure 3 As shown in , the two curves are always greater than zero, indicating that there is no collision between the two during the maneuver. In addition, the field of view angle change curve of the aircraft pointing to the target is shown in Figure 4 As shown in the figure, the field of view angle changes by less than 25° in the entire time domain, indicating that the planned translational and rotational maneuvering trajectories can ensure that the aircraft's field of view always points to the target. This also shows that the convexification method of non-convex inequality constraints given in Section 4 is effective.

[0213] The trust region convergence process after adding line search is as follows Figure 5 As shown in the figure, the algorithm finally converges after 12 iterations, obtaining a local optimal solution to the nonconvex optimal control problem. In contrast, the iterative algorithm without line search requires at least 20 iterations to converge, which is slower. This suggests that adjusting the search direction and step size between two consecutive iterative trajectories can help improve the algorithm's convergence.

[0214] Aiming at the scenario of an aircraft autonomously approaching a non-cooperative target with field-of-view pointing and collision avoidance constraints, the present invention proposes a solution to the non-convex optimal control problem based on convex programming and a time-domain non-equidistance dispersion method. In addition, to further improve the solution efficiency, a line search strategy is introduced to accelerate the convergence process of the algorithm by dynamically adjusting the step size.

[0215] The present invention uses a convex optimization method to study the path planning problem of an aircraft approaching a non-cooperative target under state constraints such as field of view pointing, control amplitude, and collision avoidance. In the process of converting the nonlinear optimization problem into a convex subproblem, in order to reduce the number of optimization variables and improve the discretization accuracy and avoid the Runge phenomenon of equidistant discretization, the present invention uses a non-uniformly spaced pseudo-spectral discretization method to discretize the non-convex optimal control problem, uses the Taylor expansion method to convexify the nonlinear motion equations and path constraints, and introduces soft trust region constraints to alleviate the infeasible solution problem that may be caused by inexact linearization. The local optimal solution of the non-convex optimal control problem is approximated by solving the convex subproblem through sequential iteration. Line search is also introduced in the iterative process to adjust the adjacent iteration step size and search direction to accelerate the convergence rate of the algorithm. The simulation example shows that this method can effectively improve the search efficiency, and the solution trajectory obtained can satisfy the non-convex inequality constraints.

[0216] In the technical solution provided by the present invention, the method includes constructing an aircraft dynamics model and constraints to obtain a non-convex optimal control problem; through a non-uniformly spaced pseudo-spectral discrete optimization model, the non-convex optimal control problem is discretized and converted into a fixed-time nonlinear programming problem about discrete nodes; the fixed-time nonlinear programming problem is convexified to form a convex subproblem, and the local optimal solution of the non-convex optimal control problem is obtained through sequence iteration. Under multi-state constraints, this method realizes the path planning of the aircraft approaching a non-cooperative target in space, alleviates the problem of no feasible solution caused by imprecise linearization, and accelerates the convergence rate of the algorithm.

[0217] Each step of the embodiment of the present invention may be performed by an electronic device, including but not limited to a tablet computer, a portable PC, a desktop computer, etc.

[0218] An embodiment of the present invention provides a computer-readable storage medium, which includes a stored program, wherein when the program is running, the electronic device where the computer-readable storage medium is located is controlled to execute the above-mentioned embodiment of the aircraft path planning method for non-cooperative targets in space.

[0219] Figure 6 A schematic diagram of an electronic device provided by an embodiment of the present invention is shown in FIG. Figure 6 As shown, the electronic device 21 includes: a processor 211, a memory 212, and a computer program 213 stored in the memory 212 and executable on the processor 211. When the computer program 213 is executed by the processor 211, the aircraft path planning method for non-cooperative space targets in the embodiment is implemented. To avoid repetition, they are not described here one by one.

[0220] The electronic device 21 includes, but is not limited to, a processor 211 and a memory 212. Those skilled in the art will understand that Figure 6 It is only an example of the electronic device 21 and does not constitute a limitation of the electronic device 21. It may include more or fewer components than shown in the figure, or a combination of certain components, or different components. For example, the electronic device may also include input and output devices, network access devices, buses, etc.

[0221] The processor 211 may be a central processing unit (CPU), or other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. A general-purpose processor may be a microprocessor or any conventional processor.

[0222] The memory 212 can be an internal storage unit of the electronic device 21, such as a hard drive or memory of the electronic device 21. The memory 212 can also be an external storage device of the electronic device 21, such as a plug-in hard drive, a Smart Media Card (SMC), a Secure Digital (SD) card, a flash memory card, etc. equipped on the electronic device 21. Furthermore, the memory 212 can include both an internal storage unit of the electronic device 21 and an external storage device. The memory 212 is used to store computer programs and other programs and data required by the network device. The memory 212 can also be used to temporarily store data that has been output or is about to be output.

[0223] Those skilled in the art will clearly understand that, for the convenience and brevity of description, the specific working processes of the systems, devices and units described above can refer to the corresponding processes in the aforementioned method embodiments and will not be repeated here.

[0224] The above description is only a preferred embodiment of the present invention and is not intended to limit 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. A method for aircraft path planning for non-cooperative targets in space, characterized in that: The method comprises: Step 1: Construct the aircraft dynamics model and constraints to obtain the non-convex optimal control problem; Step 2: Through the non-uniformly spaced pseudo-spectral discrete optimization model, the non-convex optimal control problem is discretized and converted into a fixed-time nonlinear programming problem with discrete nodes; Step 3: Convexify the fixed-time nonlinear programming problem to form a convex subproblem, and obtain the local optimal solution of the non-convex optimal control problem through sequential iteration.

2. The method according to claim 1, characterized in that The construction of the aircraft dynamics model in step 1 includes: To describe the relative motion of aircraft, we use represents the J2000 Earth Inertial Coordinate System, Represent the aircraft and target's local systems respectively, The orbital coordinate system of the target is represented by the following equations: The translational and rotational dynamic equations of the spacecraft approaching the non-cooperative target in space are expressed as follows: in, m C Indicates the mass of the aircraft; n ref for Relative to The angular velocity vector of the system is The weight in the system; express The translational motion state of the aircraft relative to the target, To express in The control quantity under the system, is the attitude control torque of the aircraft, is the moment of inertia of the aircraft, For aircraft relative The attitude of the system is corrected by the Rodrigues parameters MRPs, For the aircraft relative to The angular velocity of the system, For the aircraft relative to The relationship between the angular velocity and the speed is expressed as: in, Indicates The orbital angular velocity in the system, For Transfer to The coordinate transformation matrix of the system; using Indicates the rotational motion state of the aircraft.

3. The method according to claim 2, characterized in that The constraints constructed in step 1 include: The target spacecraft body is replaced by a sphere, and the solar panels installed on the target body are represented by an ellipsoid. The collision avoidance constraint is established, and its expression is: in, r ball represents the radius of the collision sphere, is a constant matrix, C p =[I3,0 3×9 ], Represents the collision ellipsoid feature matrix, r x ,r y ,r z is the semi-major axis of the ellipsoid, Indicates from Tie to The coordinate transformation matrix of the system; Indicates from Tie to The coordinate transformation matrix of the system is determined according to the terminal attitude of the target aircraft; Create a path constraint, the expression is: Among them, β fov is the half cone angle of the navigation equipment field of view, are the direction of the device optical axis and its position vector in this system; The control output of the aircraft has an amplitude constraint, which is expressed as: Among them, F max and T max denote the magnitudes of the control resultant force and resultant moment, respectively; Taking the optimal energy consumption as the indicator, the non-convex optimal control problem of autonomous approach is established, and its expression is: Among them, t0, t f Denote the initial and terminal moments of the maneuver, x0, x f represents the corresponding nominal state at the initial and final moments; st(1)-(3), (5)-(8) represent the constraints imposed by formulas (1)-(3) and (5)-(8).

4. The method according to claim 2, characterized in that Step 2 includes: The non-uniformly spaced pseudo-spectral discretization method is used, that is, the Gauss-Legendre points of the state variables and control variables are discretized. The Gauss-Legendre points are distributed on [-1,1]. The time interval [t0,t f ] is converted to the interval [-1,1], and a variable is introduced for this purpose. The expression is: Among them, t0, t f denotes the initial and terminal moments of the maneuver; The continuous optimal control problem after time domain transformation is expressed as: Among them, L(u τ ) represents the index function to be optimized; f(x τ ,u τ ,τ;t0,t f ) represents the right side of differential equations (1) to (3); function C(x τ ,u τ ,τ;t0,t f ) represents various path constraints; is the terminal constraint, Represent the control force and torque respectively; x τ represents the state variable; τ represents the time factor; Indicates the initial state; Indicates terminal status; Gauss-Legendre points are N-th degree Legendre polynomials P N The zero point of (τ) is expressed as: The Lagrange interpolation method is used to discretize the state variables in the continuous Bolza problem at the first N+1 LG nodes. The approximate expression of the state variables is: Among them, K represents the order of interpolation fitting, that is, how many points are used to fit the state curve; X i represents the control point, i.e. the state variable of the interpolation point; L i (τ) represents the Lagrange interpolation basis function; The Lagrange interpolation of the control variables is performed at the discrete nodes, and the approximate expression of the discretized control variables is: Among them, U i represents the control points of the control variable fitting; For the differential equation of the state variable, take the derivative of equation (12) and get: Among them, D ki equal X(τ i ) represents the column vector composed of control points; Based on the differential approximation of the state variables mentioned above, the dynamic differential equation is configured at the discrete nodes, and its expression is: Among them, f(X(τ k ),U(τ k ),τ k ;t0,t f ) represents the system differential equation, and its corresponding state variables are the control points of the state variables and the control points of the control variables; X(τ k ) represents τ k The control point of the state variable at the moment; U(τ k ) represents τ k The control point of the moment control variable; τ k Represents a time value; Using Gauss quadrature formula, we can quadrature the dynamic equation and get: Among them, ω i represents the fitting coefficient; Then the terminal state of the system at each discrete point is approximately: Among them, ω k Represents the Gauss integral coefficient, and its expression is: in, Indicates P N In τ k The derivative at ; The inequality path constraints of state variables and control variables are approximated at each discrete point as follows: C(X k ,U k ,τ k ;t0,t f )≤0(k=1,....N) (19) The terminal constraints at each discrete point are approximated as: φ(X(τ0),t0,X(τ f ),t f )=0 (20) The objective function is approximately: Among them, g(X k ,U k ) represents the indicator function; In summary, the non-convex optimal control problem is discretized and converted into a fixed-time nonlinear programming problem about discrete nodes, which is expressed as follows: Among them, st(15), (19), (20) means being constrained by formulas (15), (19), (20).

5. The method according to claim 1 or 4, characterized in that Step 3 includes: Given a reference trajectory Will The nonlinear term f(X k ,U k ,τ k ;t0,t f ) is expanded by Taylor around it, and we get: in, matrix is the Jacobian matrix of the nonlinear term about the discrete nodes, and its expression is: Bring it in In, we get: The nonlinear path constraint function C(X k ,U k ,τ k ;t0,t f ) in the reference trajectory Expand near the linearized inequality constraint, and its expression is: Among them, X C is the gradient of the constraint function at the reference trajectory; A soft trust region constraint is additionally introduced into the optimization problem to solve the problem of no feasible solution, that is, an inequality constraint is imposed, which is expressed as: Among them, η k Represents the trust region radius, which is the variable to be optimized. In order to promote the continuous shrinkage of the trust region radius during the iteration process, a related penalty term is added to the indicator function, and its expression is: In summary, the expression of the convex subproblem after convexification is: Use z to uniformly represent the optimization variable z=[X T ,U T ] T , define p m =z m+1 -z m Indicates the search direction, z m The equality and inequality constraint violations and the index function are expressed by the exact penalty function, which is expressed as follows: Among them, μ = [μ1, μ2] is the penalty factor, f E,k (z m ) represents a nonlinear equality constraint, which is expressed as: Define the step factor α∈[0,1], the line search method starts from α=1, and shrinks in sequence α′=cα, c∈(0,1), along the direction p m Search for new iteration points Until it satisfies the following non-exact one-dimensional search Armijo condition, which is expressed as: φ1(z m +αp m ;μ)≤φ1(z m ;μ)+c1αD(φ1(z m ;μ);p m ); Among them, 0<c1<1, D(φ1(z m ;μ);p m ) represents the function φ1(z m ; μ) in direction p m The directional derivative on is expressed as: In summary, when solving the non-convex optimal control problem, the initial reference trajectory is first obtained by linear interpolation. Then, in each iteration, the iteration step size and penalty factor are updated based on line search, and the coefficient matrix related to the linearized dynamics and convex path constraints is calculated. Then, the newly established convex subproblem, formula (27), is solved based on convex optimization. When the difference between adjacent iterative trajectories meets the stopping condition, the sequential optimization process is stopped. The expression of the stopping condition is: ΔX=||X m -X m-1 || ∞ ≤ε。 6. A computer-readable storage medium, characterized in that The computer-readable storage medium includes a stored program, wherein when the program is executed, the device where the computer-readable storage medium is located is controlled to execute the aircraft path planning method for a non-cooperative space target according to any one of claims 1 to 5.

7. An electronic device, characterized in that: include: one or more processors; Memory; and one or more computer programs, wherein the one or more computer programs are stored in the memory, and the one or more computer programs include instructions that, when executed by the device, cause the device to perform the aircraft path planning method for a non-cooperative target in space as described in any one of claims 1 to 5.

Images