Trajectory optimization nonlinear programming solving method supporting parallel computing

Abstract

The invention discloses a trajectory optimization nonlinear programming solving method supporting parallel computing. The method comprises the steps of 1, dispersing a trajectory optimization problem in a general form into a nonlinear programming problem; 2, constructing a first-order partial derivative matrix of a nonlinear programming problem and identifying a sparse type of the first-order partial derivative matrix; step 3, constructing a second-order partial derivative matrix of the nonlinear programming problem and identifying a sparse type of the second-order partial derivative matrix; 4, calculating discrete values of a first-order / second-order partial derivative matrix and a dynamic system by adopting a parallel mode; and step 5, evaluating and optimizing an acceleration effect, and determining an optimal parallel kernel number. According to the method, the first-order / second-order partial derivative matrix needing to be frequently evaluated in the nonlinear programming solving process and the discrete value of a dynamic system are efficiently calculated in a parallel mode, the trajectory optimization efficiency is effectively improved, the calculation time consumption is reduced, and the online application potential is enhanced.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of aircraft trajectory optimization and optimal control, and proposes a trajectory optimization nonlinear programming solving method supporting parallel computing, which has application potential in the fields of trajectory fast optimization and online optimization. BACKGROUND

[0002] Trajectory optimization is one of the key technologies in aircraft mission planning, and is widely used in the mission design and guidance control system research of various aircrafts. Such problems are essentially optimal control problems, which can usually be transformed into nonlinear programming problems (NLP) for numerical solution.

[0003] With the increasing demand for online optimization and real-time decision-making, trajectory optimization methods are required to provide higher precision solutions in a shorter time, which poses more stringent performance challenges to existing serial NLP solvers (such as IPOPT, SNOPT, etc.). In recent years, parallel computing technology has developed rapidly, and parallel algorithms based on multi-core CPUs, GPUs and distributed computing platforms have shown significant advantages in scientific and engineering calculations. By introducing parallel computing strategies into the NLP solving process of trajectory optimization, it is expected to break through the bottleneck of traditional serial framework in terms of computing speed, thereby significantly improving the solving efficiency of trajectory optimization, especially for application scenarios with high real-time requirements. However, the current mainstream NLP solvers (such as IPOPT, SNOPT) do not support direct use of parallel computing resources during problem solving, which limits their application potential in real-time trajectory optimization.

[0004] Therefore, the present application proposes a trajectory optimization nonlinear programming solving method that can utilize parallel computing resources. Through structural reconstruction of the solving process, parallel acceleration of the calculation of first and second derivative matrices and the evaluation of dynamic equations and constraints is achieved, thereby effectively improving the overall solving efficiency. SUMMARY

[0005] To achieve the above objectives, the present application proposes a nonlinear programming fast solving method supporting parallel computing. The core of this method is as follows: first, accurately identify the positions of non-zero elements in the first and second derivative matrices of the nonlinear programming problem; on this basis, evaluate the above non-zero elements and the values of the dynamic system at each discrete node simultaneously in parallel; finally, integrate this parallel computing module into the nonlinear programming solver to achieve parallel acceleration of trajectory optimization. Specifically, the present application includes the following steps:

[0006] Step 1: Discretize the general form of trajectory optimization problem into a nonlinear programming problem.

[0007] (1) Standard form of trajectory optimization problem

[0008] Trajectory optimization problem is essentially an optimal control problem. Depending on the form of performance index, optimal control problem can be classified into Mayer problem, Lagrange problem and Bolza problem. The first two can be regarded as special forms of the third. Therefore, without loss of generality, the present invention takes Bolza problem as the research object.

[0009] The general form of Bolza-type trajectory optimization / optimal control problem can be described as: solving the optimal continuous control variable u(t)∈R m , so that the objective function of the following form is minimized

[0010]

[0011] and satisfies the state equation

[0012]

[0013] End-point condition

[0014] Φ(x(t0),t0,x(t f ),t f )=0 (3)

[0015] Path constraint

[0016] C(x(t),u(t),t)≤0,t∈[t0,t f ] (4)

[0017] In formula (1)-(4), the definitions of variables M, L, f, Φ and C are as follows

[0018]

[0019] The continuous optimal control problem described by formula (1)-(4) is called Bolza-type optimal control problem.

[0020] (2) Discretization of trajectory optimization problem

[0021] Suppose that N+1 discrete points on the unit interval [0,1] are

[0022]

[0023] In the formula: τ i is called node, τ i can be uniformly distributed or non-uniformly distributed on [0,1].

[0024] The present invention uses trapezoidal format (TR format) to discretize the trajectory optimization problem, and x i = x(τ i ), ui =u(τ i The NLP optimization variables obtained by discretization are (x0,···,x N ;u0,···,u N ;γ;t0,t f ).

[0025] The discretized objective function is

[0026] J=M(x0,t0,x f ,t f )+γ (7)

[0027] The constraints for discretization are

[0028]

[0029] C i ≤0,(i=0,1,…,N) (9)

[0030]

[0031] Φ(x0,t0,x f ,t f )=0 (11)

[0032] Where: f i =f(x i ,u i ,t i ), C i =C(x i ,u i ,t i ), L i =L(x i ,u i ,t i ).

[0033] In order to facilitate the application of pseudospectral method to discretize the above general form of optimal control problem, it is necessary to discretize the time interval t∈[t0,t f ] is transformed to τ∈[-1,+1], and the transformation method is as follows

[0034]

[0035] Step 2: Construct the first-order partial derivative matrix of the nonlinear programming problem and identify its sparsity.

[0036] (1) First-order partial derivative matrix of nonlinear programming

[0037] In order to describe the sparse property of the NLP Jacobian matrix, the present application records the values of the same component of the variables or constraints at different nodes as a vector. Taking the discrete residual of the state equation as an example,

[0038] ξ :,j =(ξ 0,j ,ξ 1,j ,…,ξ N-1,j ) T ,(j=1,2,…,n) (13)

[0039] In the formula, the first subscript of ξ i,j indicates the i-th node, and the second subscript indicates the j-th component. It is known that ξ :,j is an N×1 vector. Similarly, the optimization variables C :,j , x :,j , u :,j , t : , τ : , and L : may be defined by using this method.

[0040] According to this variable recording method, the nonlinear programming problem obtained by using the trapezoidal format for discretization can be described again as follows: solving the optimization variable z∈R (N+1)·(n+m)+2 , so as to minimize the following objective function

[0041] J(z) (14) and satisfying the constraint condition

[0042] F min ≤F(z)≤F max (15)

[0043] In the formula, the expression of J(z) is shown in formula (7), and the definitions of the optimization variable z and the constraint function F(z) are as follows

[0044]

[0045] The first-order derivative matrix of the nonlinear programming is composed of the first-order derivative of the objective function with respect to the independent variable (referred to as the gradient vector) and the first-order derivative of the constraint with respect to the independent variable (referred to as the Jacobian matrix).

[0046] The gradient vector J1 of the objective function is

[0047]

[0048] The Jacobian matrix J2 of the constraint is

[0049]

[0050] (2) Sparse type of the first-order derivative matrix of the nonlinear programming

[0051] Each element in the Jacobian matrix shown in equation (18) is still a matrix block. The matrix is a sparse matrix, in which there are a large number of zero elements. In order to reduce the amount of calculation, it is necessary to identify the position of the zero elements in the matrix, and only the non-zero elements are calculated. For the state equation, the discrete residual constraint is as follows

[0052]

[0053] In the formula: the matrices D1, D2 and h are defined as follows

[0054]

[0055] The partial derivative of both sides of equation (19) with respect to is obtained

[0056]

[0057] The partial derivative of both sides of equation (19) with respect to is obtained

[0058]

[0059] The partial derivative of both sides of equation (19) with respect to γ is obtained

[0060]

[0061] The partial derivative of both sides of equation (19) with respect to time t0 is obtained

[0062]

[0063] Where

[0064]

[0065] The partial derivative of both sides of equation (19) with respect to t f is obtained

[0066]

[0067] Where

[0068]

[0069] According to equations (20)-(24), the first-order partial derivative matrix of the NLP is a function of the first-order partial derivative of the original trajectory optimization problem. Taking the state equation of the original trajectory optimization problem as an example, its partial derivative is as follows

[0070]

[0071] In the formula: each item of G1 is still a matrix, and the form is as follows

[0072]

[0073] To describe the sparsity pattern of G1, define the following struct function

[0074]

[0075] Let

[0076] S1 = struct (G1) (27)

[0077] where struct(G1) means struct operation on each element of G1; S1 means the sparsity pattern of G1.

[0078] According to the sparsity pattern of G1, the sparsity pattern of the corresponding element block in the Jacobian matrix can be obtained by means of equations (20)-(24). Similar methods can be used to obtain the expressions of other element blocks in the Jacobian matrix.

[0079] Step 3: Construct the second-order derivative matrix of the nonlinear programming problem and identify its sparsity pattern.

[0080] (1) The second-order derivative matrix of the nonlinear programming

[0081] The Lagrangian function of NLP is

[0082]

[0083] Taking the second-order derivative of it with respect to the independent variable, the second-order derivative matrix H of NLP can be obtained as

[0084]

[0085] The second-order derivative matrix of NLP shown in equation (29) is called the Hessian matrix. Since this matrix is a symmetric matrix, only the lower left part is listed here. Each element in equation (29) is still a matrix block.

[0086] (2) The sparsity pattern of the second-order derivative matrix of the nonlinear programming

[0087] The Lagrangian function L can be divided into two parts: the endpoint function term L E and the integral function term L I

[0088]

[0089] where

[0090]

[0091] In the integral function term L I ​In the equation (31), the corresponding item of the state equation discrete residual error is

[0092]

[0093] The partial derivative of the equation (31) is obtained :,i

[0094]

[0095] The partial derivative of the equation (32) is obtained

[0096]

[0097] According to the equations (31)-(36), the second-order partial derivative matrix of the NLP is the function of the second-order partial derivative of the original trajectory optimization problem. i (x, u, t), and the second-order partial derivative is as follows

[0098]

[0099] In the formula, the elements of struct (H ) are still matrix blocks, and the form is as follows

[0100]

[0101] The traditional method identifies the second-order partial derivative sparsity which is too conservative, and the application adopts the super-dual number method to accurately identify the sparsity of the second-order partial derivative of the original trajectory optimization problem. i The sparsity of H

[0102] S i (2) = struct (H i ) (38)

[0103] In the formula, struct (H i ) represents that the struct operation is performed on each element of H i .

[0104] According to the sparsity of H i , the sparsity corresponding to the discrete residual error component of the state equation in the Hessian matrix can be obtained by means of the equations (32)-(36). The similar method can be used to obtain the sparsity corresponding to other components.

[0105] Step 4: the first-order / second-order partial derivative matrix and the discrete value of the dynamic system are calculated in a parallel mode.

[0106] ​​The NLP solver needs to frequently calculate the non-zero elements in the first-order and second-order partial derivative matrices of the NLP and frequently evaluate the values of the dynamic equation, the constraint condition and the objective function at different discrete nodes when solving the NLP discretized from the trajectory optimization problem. The application efficiently calculates these values by using a parallel method, thereby improving the solving efficiency of the NLP and further improving the trajectory optimization efficiency.

[0107] Step 5: Evaluate the optimization acceleration effect and determine the optimal number of parallel cores.

[0108] The trajectory optimization parallel computing effect is related to the calculation amount of the trajectory optimization problem and the calculation capacity of the platform. For a specific trajectory optimization problem, the calculation effect in the case of different numbers of parallel cores needs to be tested on the adopted calculation platform to determine the number of parallel cores with the optimal calculation efficiency.

[0109] The method of the application uses a parallel computing technology to perform parallel acceleration on the key calculation links in the nonlinear programming, thereby improving the trajectory optimization efficiency. The method is particularly suitable for application scenarios with high real-time calculation requirements, such as online optimization of aircraft trajectories, guidance control in high dynamic environments and the like. Meanwhile, the method of the application improves the solving efficiency of the trajectory optimization problem.

[0110] The application provides a trajectory optimization nonlinear programming solving method supporting parallel computing, which has the advantage that the first-order / second-order partial derivative matrix and the values of the dynamic system at each discrete point, which need to be frequently evaluated in the nonlinear programming solving process, are efficiently calculated in a parallel manner, thereby improving the solving speed of the nonlinear programming. The method can effectively improve the trajectory optimization efficiency, reduce the calculation time consumption and improve the online application potential of the trajectory optimization. BRIEF DESCRIPTION OF DRAWINGS

[0111] Figure 1 Fig. 1 is a schematic diagram of the trajectory optimization parallel computing scheme in the embodiment of the application;

[0112] Figure 2 Fig. 2 is a height-time curve of the optimal reentry trajectory in the embodiment of the application;

[0113] Figure 3 Fig. 3 is a speed-time curve of the optimal reentry trajectory in the embodiment of the application;

[0114] Figure 4 Fig. 4 is a track angle-time curve of the optimal reentry trajectory in the embodiment of the application;

[0115] Figure 5 Fig. 5 is an attack angle-time curve of the optimal reentry trajectory in the embodiment of the application;

[0116] Figure 6 Fig. 6 is a roll angle-time curve of the optimal reentry trajectory in the embodiment of the application;

[0117] Figure 7 Three-dimensional curve of the optimal reentry trajectory and its ground projection in the embodiment of the present application;

[0118] Figure 8 Heat flow curve of the optimal reentry trajectory in the embodiment of the present application. DETAILED DESCRIPTION

[0119] In order to make the objectives, technical solutions and effects of the present application clearer and more explicit, the present application is further described in detail below by citing examples. It should be noted that the specific implementation described herein is only used to explain the present application and does not limit the present application.

[0120] The present application is specifically described below in combination with the drawings and specific implementation cases.

[0121] Step 1: Discretize the trajectory optimization problem in general form into a nonlinear programming problem.

[0122] (1) Specific form of the trajectory optimization problem

[0123] The maximum lateral range problem of a space shuttle is a classic example in the field of optimal control. Different versions of this problem have been extensively studied. Let x = (r, θ, φ, v, ψ, γ) represent the radial, longitudinal, latitudinal, velocity, heading angle and track angle respectively, then the equation set describing the motion of the mass center of the space shuttle is T

[0124]

[0125] In the formula: g is the gravitational acceleration, g = μ / r 2 , μ = 3.98603195 × 10 14 m 3 / s 2 . The tangential, normal and lateral components of the aerodynamic force acceleration a s , a n , a w are respectively

[0126]

[0127] σ is the velocity roll angle, m is the mass of the aircraft, m = 102204.6 kg. The lift and drag are respectively

[0128] L = q d AC L (α), D = q d AC D (α) (41)

[0129] Reference area A = 250 m 2 , dynamic pressure q d ​=0.5ρ(r)v 2 The lift coefficient and drag coefficient are

[0130]

[0131] The angle of attack α is in radians. The atmospheric density is calculated using the exponential model.

[0132]

[0133] Where: ρ0 = 1.225 kg / m 3 , r0=6371200.4m,

[0134] The maximum lateral range trajectory optimization problem of the space shuttle is described as: finding the optimal control u(t) = [α(t)δ(t)] T , so that the latitude at the terminal moment is minimized under the premise of satisfying dynamic constraints and boundary conditions

[0135] J=φ(t f ) (44)

[0136] Satisfying the differential equation constraints (39), initial conditions: h0 = 79248 m, θ0 = 0°, φ0 = 0°, v0 = 7802.88 m / s, ψ0 = 0°, γ0 = -1.064°; terminal conditions (TAEM window): h f =24384m,v f =762m / s,γ f = -5°; and the aerodynamic heating rate constraint q≤79.441w / cm 2 . The expression of q is as follows

[0137]

[0138] Where: N = 0.5, M = 3.07, C = 1.7827 × 10 -4 J.s 2.07 / m 3.57 / kg 0.5 , h0=1.067, h1=-1.101, h2=0.6988, h3=-0.1903.

[0139] (2) Discretization of trajectory optimization problem

[0140] Assume that the N+1 discrete points on the unit interval [0,1] are

[0141]

[0142] Where: τ i is called a node, τi It can be uniformly distributed or non-uniformly distributed on [0, 1].

[0143] The present application adopts a trapezoidal format (TR format) to discretize the trajectory optimization problem, and records x i = x (τ i ), u i = u (τ i ). The optimization variables of the NLP obtained by discretization are (x0, ···, x N ; u0, ···, u N ; γ; t0, t f ).

[0144] The discretized objective function is

[0145] J = M (x0, t0, x f , t f ) + γ (47)

[0146] The discretized constraint condition is

[0147]

[0148] C i ≤ 0, (i = 0, 1, …, N) (49)

[0149]

[0150] Φ (x0, t0, x f , t f ) = 0 (51)

[0151] In the formula: f i = f (x i , u i , t i ), C i = C (x i , u i , t i ), L i = L (x i , u i , t i ).

[0152] In order to facilitate the application of the pseudospectral method to the discretization of the optimal control problem of the above general form, it is necessary to transform the time interval t ∈ [t0, t f ] of the optimal control problem to τ ∈ [-1, +1], and the transformation is as follows

[0153]

[0154] Step 2: Construct the first-order partial derivative matrix of the nonlinear programming problem and identify its sparsity.

[0155] (1) First-order partial derivative matrix of nonlinear programming

[0156] In order to facilitate the description of the sparse characteristics of the NLP partial derivative matrix, the present invention records the values ​​of the same component of the variable or constraint at different nodes as a vector. Taking the discrete residual of the state equation as an example,

[0157] ξ :,j =(ξ 0,j ,ξ 1,j ,…,ξ N-1,j ) T ,(j=1,2,…,n) (53)

[0158] Where: i,j The first subscript of represents the i-th node, and the second subscript represents the j-th component. It is easy to see that ξ :,j is an N×1 vector. Similarly, this method can be used to define the optimization variable C :,j ,x :,j ,u :,j ,t : ,τ : ,L : .

[0159] According to this variable notation, the nonlinear programming problem obtained by trapezoidal format discretization can be re-described as: solving the optimization variable z∈R (N+1)·(n+m)+2 , so that the following objective function is minimized

[0160] J(z) (54) and satisfies the constraints

[0161] F min ≤F(z)≤F max (55)

[0162] Where: The expression of J(z) is given by equation (47), and the optimization variable z and constraint function F(z) are defined as follows:

[0163]

[0164] The first-order partial derivative matrix of nonlinear programming consists of two parts: the first-order partial derivatives of the objective function with respect to the independent variable (called the gradient vector) and the first-order partial derivatives of the constraints with respect to the independent variable (called the Jacobian matrix).

[0165] The gradient vector J1 of the objective function is

[0166]

[0167] The Jacobian matrix J2 of the constraint is

[0168]

[0169] (2) Nonlinear programming, first order partial derivative matrix sparsity

[0170] Each element in the Jacobian matrix shown in equation (58) is still a matrix block. The matrix is a sparse matrix, in which there are a large number of zero elements. In order to reduce the amount of calculation, it is necessary to identify the position of the zero elements in the matrix, and only the non-zero elements are calculated. For the state equation, the discrete residual constraint is as follows

[0171]

[0172] In the formula: matrix D1, D2 and h are defined as follows

[0173]

[0174] Taking the partial derivative of both sides of equation (59) with respect to , we get

[0175]

[0176] Taking the partial derivative of both sides of equation (59) with respect to , we get

[0177]

[0178] Taking the partial derivative of both sides of equation (59) with respect to γ, we get

[0179]

[0180] Taking the partial derivative of both sides of equation (59) with respect to time t0, we get

[0181]

[0182] Where

[0183]

[0184] Taking the partial derivative of both sides of equation (59) with respect to t f , we get

[0185]

[0186] Where

[0187]

[0188] According to equations (60)-(64), the first-order partial derivative matrix of NLP is a function of the first-order partial derivative of the original trajectory optimization problem. Take the state equation of the original trajectory optimization problem as an example, its partial derivative is as follows

[0189]

[0190] In the formula, each term of G1 is still a matrix, and the form is as follows,

[0191]

[0192] In order to describe the sparse type of G1, the following struct function is defined

[0193]

[0194] Let

[0195] S1 = struct (G1) (67)

[0196] In the formula, struct (G1) represents that the struct operation is performed on each element of G1; S1 represents the sparse type of G1.

[0197] According to the sparse type of G1, the sparse type of the corresponding element block in the Jacobian matrix can be obtained by means of equations (60)-(64). By using a similar method, the expression of other element blocks in the Jacobian matrix can be obtained.

[0198] Step 3: Construct the second-order partial derivative matrix of the nonlinear programming problem and identify its sparse type.

[0199] (1) Second-order partial derivative matrix of nonlinear programming

[0200] The Lagrangian function of NLP is

[0201]

[0202] The second-order partial derivative of the Lagrangian function with respect to the independent variable can obtain the second-order partial derivative matrix H of NLP as

[0203]

[0204] The second-order partial derivative matrix of NLP shown in equation (69) is called the Hessian matrix. Since the matrix is a symmetric matrix, only the lower left part is listed here. Each element in equation (69) is still a matrix block.

[0205] (2) Sparse type of second-order partial derivative matrix of nonlinear programming

[0206] The Lagrangian function L can be divided into two parts, the endpoint function term L E and the integral function term L I ​

[0207]

[0208] where

[0209]

[0210]

[0211] In the integral function term L I , the state equation discrete residual error corresponding term is

[0212]

[0213] The partial derivative of equation (71) with respect to x :,i is obtained

[0214]

[0215] The partial derivative of equation (72) with respect to and t0is obtained

[0216]

[0217] According to equations (72)-(76), the second-order partial derivative matrix of the NLP is a function of the second-order partial derivative of the original trajectory optimization problem. For the i-th state equation f i (x, u, t), its second-order partial derivative is as follows

[0218]

[0219] In the formula: The elements of H are still matrix blocks, and the form is as follows

[0220]

[0221] The traditional method identifies the second-order partial derivative sparsity type too conservatively, and the present application accurately identifies the second-order partial derivative sparsity type of the original trajectory optimization problem by using the hyperdual number method. Denote the sparsity type of H i as

[0222] S i (2) = struct (H i ) (78)

[0223] In the formula: struct (H i ) represents performing a struct operation on each element of H i .

[0224] According to H iThe sparse form of the state equation residual component in the Hessian matrix can be obtained by means of equations (72)-(76). The sparse forms of other components can be obtained by using a similar method.

[0225] Step 4: Calculate the first-order / second-order partial derivative matrix and the discrete value of the dynamic system in a parallel manner.

[0226] When solving the NLP obtained by discretizing the trajectory optimization problem, the nonlinear programming (NLP) solver needs to frequently calculate the values of the dynamic equation, the constraint condition and the objective function at different discrete nodes, and calculate the non-zero elements in the first-order and second-order partial derivative matrices of the NLP. The present application efficiently calculates these values by using a parallel method, thereby improving the solving efficiency of the NLP, and further improving the trajectory optimization efficiency, as shown in Figure 1

[0227] Step 5: Evaluate the optimization acceleration effect and determine the optimal number of parallel cores.

[0228] The parallel computing effect of trajectory optimization is related to the calculation amount of the trajectory optimization problem and the computing capacity of the platform. The parallel computing method described in the present application is tested on a desktop computer. The computer is configured as: memory 16 GB, processor 12th Gen Intel(R) Core(TM) i7-12700 2.10 GHz. Table 1 gives the comparison of the time consumption of serial and parallel computing under different discrete node numbers, and Table 2 gives the comparison of the corresponding optimal objective function. Table 1 shows that the parallel computing method proposed in the present application can significantly reduce the total calculation time of trajectory optimization under different discrete node numbers, and the maximum calculation time can be compressed to 40% to 50% of the original serial time consumption. In the case of using 2 to 8 threads, the speedup gradually increases with the increase of the number of threads, showing good parallel scalability. However, when the number of threads is further increased, the speedup decreases, which is because the thread scheduling and synchronization overhead in parallel computing increases, offsetting part of the performance benefits brought by parallel computing. For this example, the optimal speedup effect can be obtained by using 8 threads. In addition, as the number of discrete nodes increases, the speedup also increases. This is because the more nodes there are, the larger the calculation scale of the trajectory optimization problem, thereby making the advantage of parallel computing in load distribution more significant, and the overall parallel efficiency improves. Table 2 further verifies that the method of the present application improves the calculation efficiency without affecting the objective function value of the trajectory optimization problem, indicating that the proposed parallel strategy effectively improves the solving efficiency without sacrificing the optimization accuracy.

[0229] Table 1 Comparison of time consumption of serial and parallel computing of space shuttle reentry problem

[0230] ​

[0231]

[0232] Table 2 Target function comparison between serial and parallel computation of space shuttle reentry problem

[0233]

[0234] Figures 2-8 The optimal trajectory curve is given, wherein Figure 2 is a height-time curve, Figure 3 is a velocity-time curve, Figure 4 is a flight-path angle-time curve, Figure 5 is an angle-of-attack-time curve, Figure 6 is a roll angle-time curve, Figure 7 is a three-dimensional trajectory curve, Figure 8 is a heat flow-time curve. The trajectory is smooth and continuous, and the heat flow strictly meets the constraint requirement, which indicates that the trajectory optimization method has good engineering applicability and precision.

[0235] The above case verifies the feasibility of the method of the present application in improving the trajectory optimization efficiency by parallel computation of first-order and second-order partial derivative matrices of nonlinear programming and dynamics equations at different discrete points.

[0236] The above only describes the preferred embodiments of the present application, and it should be noted that, for those skilled in the art, several improvements can be made without departing from the principles of the present application, and these improvements should also be considered as the protection scope of the present application.

Claims

1. A trajectory optimization nonlinear programming solution method supporting parallel computing, characterized in that: The method described is: Step 1: Discretize the general trajectory optimization problem into a nonlinear programming problem; Step 2: Construct the first-order partial derivative matrix of the nonlinear programming problem and identify its sparsity; Step 3: Construct the second-order partial derivative matrix of the nonlinear programming problem and identify its sparsity; Step 4: Calculate the first-order / second-order partial derivative matrix and the discrete values ​​of the dynamic system in parallel; Step 5: Evaluate the optimization acceleration effect and determine the optimal number of parallel cores.

2. The trajectory optimization nonlinear programming solution method supporting parallel computing according to claim 1, characterized in that: The step 1 is specifically as follows: 1.1: Standard form of trajectory optimization problem (Bolza problem); The general form of the Bolza-type trajectory optimization / optimal control problem can be described as: solving the optimal continuous control variable u(t)∈R m , so that the objective function of the following form is minimized And satisfy the state equation Endpoint conditions Φ(x(t0),t0,x(t f ),t f )=0 (3) Path Constraints C(x(t),u(t),t)≤0,t∈[t0,t f ] (4) In equations (1)-(4), the variables M, L, f, Φ, and C are defined as follows: The continuous optimal control problem described by equations (1)-(4) is called the Bolza type optimal control problem. 1.2: Discretization of trajectory optimization problem; Assume that the N+1 discrete points on the unit interval [0,1] are Where: τ i is called a node, τ i It can be uniformly distributed or non-uniformly distributed on [0,1].

3. The trajectory optimization nonlinear programming solution method supporting parallel computing according to claim 2, characterized in that: The trajectory optimization problem is discretized using the trapezoidal format (TR format), and x is i =x(τ i ),u i =u(τ i ); The NLP optimization variables obtained by discrete are (x0,···,x N ;u0,···,u N ;γ;t0,t f ); The discretized objective function is J=M(x0,t0,x f ,t f )+γ (7) The constraints for discretization are C i ≤0,(i=0,1, …, N) (9) Φ(x0,t0,x f ,t f )=0 (11) where: f i = f(x i , u i , t i ), C i = C(x i , u i , t i ), L i = L(x i , u i , t i ); In order to facilitate the application of pseudospectral method to discretize the above general form of optimal control problem, it is necessary to discretize the time interval t∈[t0,t f ] is transformed to τ∈[-1,+1], and the transformation method is as follows 4. The trajectory optimization nonlinear programming solution method supporting parallel computing according to claim 1, characterized in that: The step 2 is specifically as follows: 2.1: First-order partial derivative matrix for nonlinear programming; The values ​​of the same component of a variable or constraint at different nodes are recorded as a vector, and the discrete residual of the state equation is: x :,j =(ξ 0,j ,x 1,j ,…,x N-1,j ) T ,(j=1,2, …, n) (13) Where: i,j The first subscript of represents the i-th node, the second subscript represents the j-th component, :,j is the j component at different nodes; similarly, the variable C can be defined in this way :,j ,x :,j ,u :,j ,t : ,τ : ,L : ; According to this variable notation, the nonlinear programming problem obtained by trapezoidal format discretization can be re-described as: solving the optimization variable z∈R (N+1)·(n+m)+2 , so that the following objective function is minimized J(z)(14) and satisfy the constraints F min ≤F(z)≤F max (15) Where: The expression of J(z) is given by equation (7), and the optimization variable z and constraint function F(z) are defined as follows: The first-order partial derivative matrix of nonlinear programming consists of two parts: the first-order partial derivative of the objective function with respect to the independent variable (called the gradient vector) and the first-order partial derivative of the constraint with respect to the independent variable (called the Jacobian matrix); The gradient vector J1 of the objective function is The Jacobian matrix J2 of the constraint is 2.2: Sparse form of first-order partial derivative matrix for nonlinear programming; Each element in the Jacobian matrix shown in equation (18) is still a matrix block; The matrix is ​​a sparse matrix, in which a large number of elements are zero. In order to reduce the amount of calculation, it is necessary to identify the position of the zero elements in the matrix and only calculate the non-zero elements. For the state equation, the discrete residual constraint is as follows Where: matrices D1, D2 and h are defined as follows Compare both sides of equation (19) to Taking the partial derivative we get Compare both sides of equation (19) to Taking the partial derivative we get Taking partial derivatives of both sides of equation (19) with respect to γ, we can obtain Taking partial derivatives of both sides of equation (19) with respect to time t0, we get in Compare both sides of equation (19) to t f Taking the partial derivative we get in According to equations (20)-(24), the first-order partial derivative matrix of NLP is a function of the first-order partial derivative of the original trajectory optimization problem; Taking the state equation of the original trajectory optimization problem as an example, its partial derivative is as follows Where: Each item of G1 is still a matrix, in the following form To describe the sparse type of G1, define the following struct function remember S1=struct (G1) (27) Where: struct(G1) means performing struct operation on each element of G1; S1 represents the sparse type of G1; According to the sparse type of G1, the sparse type of the corresponding element block in the Jacobian matrix can be obtained with the help of equations (20)-(24); the expressions of other element blocks in the Jacobian matrix can be obtained according to the above method.

5. The trajectory optimization nonlinear programming solution method supporting parallel computing according to claim 1, wherein the step 3 is specifically: 3.1: Second-order partial derivative matrix for nonlinear programming; The Lagrangian function of NLP is Taking the second-order partial derivative of the independent variable, we can get the second-order partial derivative matrix H of NLP: The NLP second-order partial derivative matrix shown in formula (29) is called the Hessian matrix. Since this matrix is ​​a symmetric matrix, formula (29) only lists the lower left part; each element in formula (29) is still a matrix block; 3.2: Sparse form of the second-order partial derivative matrix for nonlinear programming; The Lagrangian function L can be divided into the endpoint function term L E and the integral function term L I Two parts in In the integral function term L I In the equation of state, the discrete residual corresponds to Reverse equation (31) to x :,i Taking the partial derivative we get Convert equation (32) to Taking the partial derivative of t0 we get According to equations (31)-(36), the second-order partial derivative matrix of NLP is a function of the second-order partial derivatives of the original trajectory optimization problem; for the i-th state equation f i (x,u,t), its second-order partial derivative is as follows Where: The elements are still matrix blocks, in the following form The superdual number method is used to accurately identify the sparse form of the second-order partial derivatives of the original trajectory optimization problem; denoted by H i The sparse type is S i (2) =struct (H i ) (38) Where: struct(H i ) indicates H i Perform struct operation on each element of According to H i The sparse type of the discrete residual component of the state equation in the Hessian matrix can be obtained by using equations (32)-(36); the sparse types corresponding to other components can be obtained according to the above method.

6. The trajectory optimization nonlinear programming solution method supporting parallel computing according to claim 1, characterized in that: Specifically, step 4 is as follows: when the NLP solver solves the NLP obtained by discretizing the trajectory optimization problem, it is necessary to frequently calculate the non-zero elements in the first-order and second-order partial derivative matrices of the NLP, and frequently evaluate the values ​​of the dynamic equations, constraints, and objective functions at different discrete nodes.

7. The trajectory optimization nonlinear programming solution method supporting parallel computing according to claim 1, characterized in that: The specific step 5 is as follows: the parallel computing effect of trajectory optimization is related to the computational complexity of the trajectory optimization problem and the computing power of the platform; for a specific trajectory optimization problem, the computing effect of using different numbers of parallel cores is actually tested on the computing platform used to determine the number of parallel cores with the best computing efficiency.

Images