Aircraft online trajectory optimization method based on parallel real-time interior point method of nested partitioning

Abstract

The application provides an aircraft online trajectory optimization method based on a parallel real-time interior point method of nested partitioning, which firstly converts the aircraft online trajectory optimization problem into a standard second-order cone programming problem; then solves the SOCP problem by using a primal-dual interior point method based on homogeneous self-dual embedding and scaling, rearranges the coefficient matrix of the KKT sparse linear system by using a nested partitioning algorithm, and obtains a balanced elimination tree; based on the elimination tree, the LDL matrix decomposition and triangular system solving tasks are divided into sub-tree / trunk sub-tasks, and each sub-tree calculation and trunk task calculation are distributed to each CPU calculation core, so that the calculation tasks of each CPU calculation core are very balanced, thereby improving the calculation efficiency; the application has high parallelism and good parallel performance, can realize a calculation speedup ratio of 3.2-3.5 times on a DSP platform, is suitable for medium / large-scale problems, has low memory occupation, and is suitable for rocket-borne / airborne embedded environments.

Description

Technical Field

[0001] This invention relates to the field of aircraft technology, and in particular to an online trajectory optimization method for aircraft based on a nested partitioning parallel real-time interior point method. Background Technology

[0002] Onboard trajectory optimization (ATO) refers to the technology of dynamically calculating the flight trajectory of an aircraft based on real-time environmental and status information during flight. ATO offers significant advantages over traditional offline trajectory optimization, enabling real-time responses to uncertainties such as changes in weather conditions, system failures, and mission adjustments.

[0003] In recent years, online trajectory optimization technology has been widely applied in various aircraft fields. For example, it has been used in the defense industry to enable atmospheric reentry, terminal guidance, and short-term trajectory correction for various missiles to accurately hit targets. It has also been applied to UAV arm patch and formation trajectory optimization to maintain formation and avoid collisions during patrol, mapping, and search and rescue missions.

[0004] Trajectory optimization is essentially an optimal control problem. Methods for solving optimal control problems are mainly divided into indirect and direct methods. Indirect methods, based on Pontryagin's maximum principle, first derive the necessary analytical conditions and then solve the two-point boundary value problem. Direct methods discretize the continuous-time optimal control problem, representing the trajectory (state and control variables) as finite-dimensional optimization variables, thus transforming the problem into a finite-dimensional nonlinear programming (NLP) problem. To meet real-time solution requirements, lossless convexity and sequential convexity techniques are often used to further transform the problem into solving a series of convex subproblems. The application of convex optimization theory is a major breakthrough in achieving online trajectory optimization. Numerous projects and tasks both domestically and internationally have adopted online real-time trajectory planning technology based on convex optimization to achieve millisecond-second-level solutions for online trajectories.

[0005] Interior point methods have become the mainstream approach for solving convex optimization problems due to their robustness and fast convergence. However, because interior point methods require iterative solutions to the KKT system using Newton's method, they suffer from high computational cost and difficulty in utilizing multi-core architectures.

[0006] Currently widely used commercial solvers such as Mosek and Gurobi, while advantageous for solving large-scale problems, suffer from limitations due to their complex implementations, numerous dependencies, and closed-source nature. These limitations prevent them from being customized for specific online trajectory optimization problems and make them unsuitable for solving convex problems on shipborne platforms. Based on the diverse needs of aerospace missions, solvers suitable for solving various convex problems on shipborne platforms have emerged, such as ECOS, OSQP, and QOCO. However, current solvers for rocket-borne platforms generally employ the primal-dual interior-point method, a second-order method. A challenge of this method lies in obtaining the search direction by solving a sparse linear system.

[0007] When solving sparse linear systems using the primal-dual interior-point method to optimize the search direction, most embedded solvers currently employ sparsity mining and utilization strategies, exemplified by the Approximate Minimum Degree (AMD) algorithm, along with customized explicit code generation strategies. However, the AMD algorithm leads to imbalanced elimination tree structures, making parallelization difficult; customized code is prone to instruction cache misses in medium to large-scale problems, resulting in significant performance degradation. Furthermore, existing methods largely focus on parallelizing matrix factorization, neglecting the parallel potential of solving triangular systems and failing to fully utilize the computational power of multi-core CPUs. Summary of the Invention

[0008] To address the shortcomings of existing technologies, this invention provides a method for online trajectory optimization of aircraft based on a parallel real-time interior point method using nested partitioning. This method aims to solve the technical problems of insufficient parallelism, imbalanced elimination trees, poor adaptation for medium to large scales, and high synchronization overhead in existing onboard embedded convex optimization solvers. It achieves efficient parallelization of the KKT system solution and improves the real-time performance of online trajectory optimization for aircraft.

[0009] The technical solution of this invention is: an online trajectory optimization method for aircraft based on a parallel real-time interior point method with nested partitioning, comprising the following steps:

[0010] S1) Establish the dynamic model and constraints of the aircraft;

[0011] S2) Reconstruct the online trajectory optimization problem of the aircraft thrust reduction fault into a convex optimization problem, and transform the linear trajectory optimization problem into a standard SOCP problem that can be directly solved by a convex solver;

[0012] S3) The final search direction is obtained by iteratively solving a series of perturbation KKT condition equations for the primal and dual problems using the prediction-correction primal-dual path tracking interior point method based on homogeneous self-dual embedding HSD and nested partitioning NT scaling techniques.

[0013] S4) The permutation matrix is ​​calculated using the nested partitioning ND algorithm. And based on the direct method of LDL matrix decomposition, the sparse KKT linear system is solved, and the system coefficient matrix of the sparse KKT linear system is... Rearrange;

[0014] S5), based on rearrangement matrix Elimination tree structure The LDL matrix decomposition task is divided into subtree decomposition tasks and trunk decomposition tasks, and the subtree decomposition tasks and trunk decomposition tasks are statically allocated to different computing cores of a multi-core CPU.

[0015] S6) Each core executes the subtree decomposition task in parallel, independently updating the trunk nodes through the submatrix replica mechanism, and then the designated core executes the trunk decomposition task to complete the process. break down;

[0016] S7) Based on the decomposition factor of the LDL matrix, each computational core is based on its local... The submatrix copy performs forward substitution, diagonal scaling, and backward substitution operations in parallel to obtain the solution of the sparse KKT linear system, and then completes the iteration of the original dual interior point method, finally outputting the online optimized trajectory of the aircraft.

[0017] Preferably, in step S1), the dynamic model is:

[0018]

[0019] In the formula, , , The velocity components of the spacecraft in the orbit-inertial coordinate system; , () is the first derivative of the position component of the spacecraft in the orbit-inertial coordinate system; , , These are the velocity components of the spacecraft in the orbit-inertial coordinate system. , , The first derivative; , This is the component of gravitational acceleration; For the mass of the aircraft; The magnitude of the thrust is constant. , These are the commanded pitch angle and yaw angle, respectively. This refers to the fuel consumption of an aircraft engine per unit time. This is a fuel consumption constant.

[0020] Preferably, in step S1), the constraints include:

[0021] Terminal constraints:

[0022]

[0023] In the formula, , These represent the terminal position components of the aircraft; , These represent the terminal velocity components of the aircraft; The gravitational constant of Earth;

[0024] Initial state and control constraints:

[0025]

[0026] In the formula, Indicates the initial time. State variables; Indicates the initial time. The state value; Indicates the initial time. Control variables; Indicates the initial time. The value of the control variable;

[0027] Path constraints for the rate of change of commanded attitude angles:

[0028]

[0029] In the formula, and These represent the rates of change of the commanded pitch and yaw angles, respectively. Upper limit.

[0030] Preferably, in step S1), the terminal constraint is relaxed and the relaxation term is penalized in the performance index. The relaxed terminal constraint is expressed as follows:

[0031]

[0032] In the formula, , These are slack variables; , These represent the spacecraft in orbit. Terminal position component in the direction; These represent the spacecraft in the orbital plane. Terminal velocity component in the direction; The gravitational constant of Earth;

[0033] The optimization objective is to maximize the altitude of the insertion point, and a weighted penalty is applied to all slack variables in the performance metrics. The augmented performance metric objective function is as follows. Represented as:

[0034]

[0035] In the formula, These are the weighting coefficients; , These represent the terminal position components of the aircraft. This indicates the speed of the aircraft in the vertical orbital plane; , These are slack variables.

[0036] Preferably, in step S3), in each iteration of the prediction-correction primary-dual path tracking interior point method, the final search direction is obtained by solving the sparse KKT linear system, wherein the coefficient matrix of the sparse KKT linear system is a matrix with... The sparse KKT linear system, in the form of a block-based sparse symmetric matrix, is represented as follows:

[0037]

[0038] In the formula, This represents the coefficient matrix of the original linear system; The right-hand term; It is a regularization parameter with a small value. , It is the identity matrix; , For regularization terms; Represent the solution of a linear system; These represent the various parts of the solution to the linear system; This represents the components of the right-hand side of a linear system.

[0039] Preferably, in step S4), the nested partitioning ND algorithm is used to calculate the permutation matrix. Specifically, it includes the following steps:

[0040] S41) The system coefficient matrix Using the corresponding adjacency graph as input, a vertex separator set is selected to divide the adjacency graph into two approximately equal-sized subgraphs. With the partition set ;

[0041] S42), recursively execute step S41), respectively on the subgraph The subgraphs are segmented until their size is less than a preset threshold. Subgraphs smaller than the threshold are then rearranged using the approximate minimum degree AMD algorithm to reduce filler elements.

[0042] S43), according to the subgraph Separated Sets Renumber the vertices of the adjacency graph in the following order to generate a permutation matrix. .

[0043] Preferably, in step S4), the direct method for solving sparse KKT linear systems based on LDL matrix decomposition includes three stages: offline preprocessing, online numerical decomposition, and online solution.

[0044] In the offline preprocessing stage, a permutation matrix is ​​sought. To the system coefficient matrix The rows and columns are rearranged; additionally, the sparsity pattern of the factorization factors is determined by performing symbolic decomposition and the necessary memory space is allocated in advance, and the rearranged matrix is ​​calculated in the online numerical decomposition stage. The factorization factors L and D are such that:

[0045]

[0046] In the formula, the decomposition factors L and D are the triangular matrix and the diagonal matrix, respectively; This represents the coefficient matrix of the original linear system; This represents the rearranged coefficient matrix; This represents the permutation matrix used for rearrangement;

[0047] The online solution phase obtains the solution by solving a series of equations whose computational cost is less than a cost threshold. ,Right now:

[0048]

[0049] In the formula, It is an intermediate variable in the calculation process; Represents the final solution of a linear system; Represents the right-hand side term of the original linear system; This represents the triangular and diagonal matrices after decomposition.

[0050] Preferably, in step S5), based on the rearrangement matrix... Elimination tree structure The LDL matrix decomposition task is divided into subtree decomposition tasks and trunk decomposition tasks, and the subtree decomposition tasks and trunk decomposition tasks are statically allocated to different computing cores of a multi-core CPU; specifically, the following steps are included:

[0051] S51), Rearrangement Matrix The elimination tree structure is then arranged in postorder tracing.

[0052] S52) Divide the nodes of the elimination tree after the post-order arrangement into subtree nodes and trunk nodes. The subtree nodes constitute independent subtree tasks, and the trunk nodes constitute trunk tasks that depend on the results of the subtree tasks.

[0053] S53) Perform two rounds of ND decomposition, in which the elimination tree is divided into 4 independent subtree tasks and 3 trunk tasks. The difference in the number of nodes in each subtree task is less than a preset ratio, thereby achieving load balancing of multi-core CPUs.

[0054] Preferably, in step S6), the LDL matrix decomposition is performed in parallel using a column-block-based left-look decomposition algorithm, specifically including the following steps:

[0055] S61) For each computing core, for its assigned subtree task, perform column block decomposition: traverse the corresponding column blocks of the subtree column by column, and initialize the column data as follows: For the corresponding column, update the current column with the calculated column based on the elimination tree dependency relationship, and then perform diagonal division on the non-zero elements of the current column;

[0056] S62) After the subtree decomposition is completed, each computing core performs column block updates in parallel: the subtree decomposition results are used to update the corresponding columns of the trunk nodes that it depends on;

[0057] S63) For trunk tasks, each computing core performs column block decomposition based on the updated trunk column data. If multiple cores are updating the same trunk node, the process is handled locally. The submatrix replica stores the update results, which are ultimately merged by the core responsible for trunk decomposition.

[0058] Preferably, in step S7), based on the decomposition factors of the LDL matrix, each computational core is based on its local... The submatrix replicas perform forward substitution, diagonal scaling, and backward substitution operations in parallel to obtain the solution of the sparse KKT linear system, thereby completing the iteration of the original-dual interior-point method and finally outputting the online optimized trajectory of the aircraft; specifically, the following steps are included:

[0059] S71) Divide the right-hand term of the sparse KKT linear system into segments according to the elimination tree nodes, allocate corresponding segmented data to each computational core, and initialize the local right-hand term.

[0060] S72) All computing cores execute Lsolve operations in parallel to solve the problem. Segment the right-hand side terms to obtain intermediate variables. ;

[0061] S73) Perform a multi-core synchronization, where each computing core updates the shared trunk node with its own Lsolve result. value;

[0062] S74) All computing cores execute the Dsolve operation in parallel to solve the problem. To obtain intermediate variables ;

[0063] S75) All computing cores execute the LTsolve operation in parallel to solve the problem. ; Obtain intermediate variables ;

[0064] S76), combining all computing cores Segmented data, processed by permutation matrix The inverse transform yields the final solution of the sparse KKT linear system.

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

[0066] 1. This invention has high parallelism and good parallel performance. It uses a nested partitioning algorithm to calculate the permutation matrix of a sparse linear system, which can obtain a highly balanced elimination tree structure. It distributes the calculation of each subtree and the trunk task to each CPU core, making the calculation tasks of each CPU core very balanced, thereby improving the calculation efficiency.

[0067] 2. The present invention has low synchronization overhead. By allocating a local submatrix of L to each computing CPU core, the synchronization overhead is further reduced.

[0068] 3. This invention achieves a speedup of 3.2-3.5 times on a 4-core DSP platform, reducing the KKT system solution time from 12ms to 3.6ms, meeting the real-time requirements of online optimization; the generated balanced elimination tree ensures that the difference in workload between cores is ≤10%, and the multi-core resource utilization rate reaches over 90%;

[0069] 4. This invention is not only applicable to aircraft fault reconstruction, but can also be extended to online trajectory optimization of other aircraft such as missiles and drones, and supports various convex problems such as QP and SOCP. Attached Figure Description

[0070] Figure 1 This is a schematic flowchart of the method of the present invention;

[0071] Figure 2 In this embodiment of the invention, the nested partitioning ND algorithm is used to calculate the permutation matrix. A flowchart; where (a) is... (b) is a schematic diagram of the block ND rearrangement matrix; (c) is a schematic diagram of the block decomposition factor; (d) is a schematic diagram of the elimination tree and its simplified block structure.

[0072] Figure 3This is a flowchart illustrating the decomposition of the tree structure elimination task in an embodiment of the present invention; where (a) is the rearrangement matrix. (a) is a schematic diagram of the decomposition factor L; (b) is a schematic diagram of the rearrangement matrix. (d) is a schematic diagram of the elimination tree structure; (d) is a schematic diagram of the elimination tree structure in the subsequent order.

[0073] Figure 4 This is a schematic diagram of the elimination tree after two rounds of ND partitioning and the sub-elimination tree structure obtained by copying the trunk node in an embodiment of the present invention;

[0074] Figure 5 This is a schematic diagram of the decomposition factor L after two rounds of ND partitioning and the submatrix obtained by copying the trunk column data in an embodiment of the present invention;

[0075] Figure 6 This is a schematic diagram of the computation process of the parallel LDL matrix decomposition algorithm in an embodiment of the present invention;

[0076] Figure 7 This is a schematic diagram of task scheduling based on the parallel LDL matrix factorization algorithm in an embodiment of the present invention. Detailed Implementation

[0077] The specific embodiments of the present invention will be further described below with reference to the accompanying drawings:

[0078] like Figure 1 As shown, this embodiment provides an online trajectory optimization method for aircraft based on a parallel real-time interior point method with nested partitioning, including the following steps:

[0079] S1) Establish the dynamic model and constraints of the aircraft; in this embodiment, the aircraft is a launch vehicle; specifically, it includes the following steps:

[0080] S11) Establish a three-degree-of-freedom dynamic model of the aircraft during its vacuum flight phase, namely:

[0081]

[0082] In the formula, , , The velocity components in the orbit-inertia coordinate system; , () is the first derivative of the position component in the orbit-inertia coordinate system; , , These are the velocity components in the orbit-inertia coordinate system. , , The first derivative; , This is the component of gravitational acceleration; For the mass of the aircraft; The magnitude of the thrust is constant. , These are the commanded pitch angle and yaw angle, respectively. This refers to the fuel consumption of an aircraft engine per unit time. This is a fuel consumption constant.

[0083] In this embodiment, the state variable is defined as: Control variables are defined as follows: The gravitational acceleration of the aircraft is based on the formula. The calculation of its components is expressed as follows:

[0084]

[0085] In the formula, is the Earth's gravitational constant.

[0086] The thrust reduction fault model is as follows:

[0087]

[0088] In the formula, This represents the percentage of thrust remaining after a thrust failure, with a value range of [value missing]. ; Indicates the nominal specific impulse of a rocket engine; This is a fuel consumption constant.

[0089] S12) Establish endpoint and path constraints;

[0090] In order to form a circular rescue trajectory coplanar with the target trajectory, the terminal constraints of the spacecraft at the entry point need to satisfy the following constraints:

[0091]

[0092] In the formula, , These represent the terminal position components of the aircraft; , These represent the terminal velocity components of the aircraft; The gravitational constant of Earth;

[0093] After experiencing a thrust failure, the aircraft begins executing a trajectory planning algorithm, therefore the following initial state and control constraints must be met:

[0094]

[0095] In the formula, Indicates the initial time. State variables; Indicates the initial time. The state value; Indicates the initial time. Control variables; Indicates the initial time. The value of the control variable;

[0096] Considering the limitations of the aircraft's attitude control system, the following path constraint is applied to the rate of change of the commanded attitude angle:

[0097]

[0098] In the formula, and These represent the rates of change of the commanded pitch and yaw angles, respectively. Upper limit.

[0099] S13) Establish the optimization objective function;

[0100] Due to the uncertainty of the time of the aircraft failure and the percentage of thrust reduction, it is difficult to determine whether the nominal orbit is reachable. Therefore, it is necessary to jointly optimize the rescue orbit and the flight trajectory. For the circular target orbit, this embodiment adopts relaxed terminal constraints and penalizes the relaxation term in the performance index. The relaxed terminal constraints are expressed as follows:

[0101]

[0102] In the formula, , These are slack variables; , These represent the spacecraft in orbit. Terminal position component in the direction; These represent the spacecraft in the orbital plane. Terminal velocity component in the direction; The gravitational constant of Earth;

[0103] Since the trajectory optimization problem is described in the orbital coordinate system at the entry point, the terminal state variables in the z-axis direction can be directly constrained to make the rescue trajectory plane approximate the nominal trajectory plane. and The value is a relatively small value.

[0104] To maximize the likelihood of payload remaining on track after a fault, the optimization objective of this embodiment is to maximize the insertion point altitude. Weighted penalties are applied to all relaxed variables in the performance metrics, resulting in the objective function of the augmented performance metrics. The form of expression:

[0105]

[0106] In the formula, Here, the weighting coefficients are: and The value is used to adjust the shape of the reconstructed orbit, while the weighting coefficient... and The value of affects the orientation of the reconstructed orbital plane; by adjusting the values ​​of each weight in the performance indicators, the distribution direction of the remaining energy of the spacecraft can be changed, and the weighted optimization of the orbital plane, orbital altitude and orbital shape can be achieved, thereby reconstructing different rescue orbits.

[0107] S2) Reconstruct the online trajectory optimization problem for aircraft thrust reduction faults into a convex optimization problem and convert it into a standard SOCP problem. The convex optimization problem is expressed as follows:

[0108]

[0109] In the formula, This represents the state variables in the optimal control problem. The control input for the optimal control problem is represented by the first-order hold-and-discrete method. The state variables and control variables are discretized to various discrete points. Indicates the first The state of a discrete point, including position and velocity information. Indicates the first Control quantity at a discrete point; Represents the number of discrete points, for the terminal. , , , , These represent the orbital / inertial coordinate systems in the terminal state, respectively. Position and velocity in the direction;

[0110] The dynamic model of the original optimal control problem can be transformed into the following form of equality constraints:

[0111]

[0112] in, , , These represent the parameters of the dynamic equation.

[0113] The objective function of the original optimal control problem can be transformed as follows: The objective of maximizing the altitude of the entry point is equivalently transformed into minimizing the velocity of the entry point, that is... Transform into: Introducing slack variables A second-order cone constraint can be constructed. and will Add it to the objective function. Similarly, minimize it. , , , , and then introduce slack variables respectively , , More second-order cone constraints can be constructed to replace the quadratic objective function with a combination of second-order cone constraints and linear objective functions.

[0114] Initial value constraints can be established by specifying the initial state values ​​for the state and control. , , in The state variable representing the first discrete point. This represents the initial state value. The control variable representing the first discrete point, This represents initial control.

[0115] For relaxed circular orbit constraints, the original constraints are nonlinear constraints, which can be linearized into linearized constraints through sequence linearization, and their form is as follows:

[0116]

[0117] in, , , , Represents terminal location , Position and velocity variables in the direction; , It is a slack variable; , , , Representing the terminal location respectively Reference position and reference velocity parameters in the direction; , , , , , , , , , These represent other linear constraint parameters, respectively.

[0118] Establish constraints for control amplitude It can be converted into a second-order cone constraint: , in , , Representing discrete points Control quantities in all directions.

[0119] Establish constraints to control the rate of change:

[0120]

[0121] in , , Representing discrete points Control quantities in all directions; , , Representing discrete points The reference control parameters are the control parameters from the previous iteration. , , These represent the control constraint parameters.

[0122] Establish trust region constraints for the terminal's location and speed:

[0123]

[0124] in, , , , Represents terminal location , Position and velocity variables in direction, , , , Representing the terminal location respectively Reference position and reference velocity parameters in the direction, , , , These represent the trust region parameters.

[0125] For the optimization problem established in this embodiment, it is easy to see that this problem is a second-order cone programming problem, which involves converting the state variables of each discrete point... Control variables at each discrete point Various slack variables , , , , , , As optimization variables. The objective function is a linear function that satisfies In this form, the initial value constraints, dynamic constraints, and linearized circular trajectory constraints are linear equality constraints, satisfying... The control transform rate constraint and trust region constraint are linear inequality constraints, satisfying the form; The relaxed objective function and control amplitude constraint are second-order cone constraints, which can be easily written as... In the form of. For Introducing slack variables It can be transformed into: Therefore, we can obtain the standard SOCP problem that can be directly solved by a convex solver; that is:

[0126]

[0127] The dual problem of the second-order cone programming problem is expressed as:

[0128]

[0129] in, , , , , It is about optimizing the parameters of the problem; , These are the optimization variables of the original problem; , These are the optimization variables for the dual problem; This represents a cone composed of linear inequality constraints and various second-order cone constraints, where... , .

[0130] S3) The final search direction is obtained by iteratively solving a series of perturbation KKT condition equations for the primal and dual problems using the prediction-correction primal-dual path tracking interior point method based on homogeneous self-dual embedding HSD and nested partitioning NT scaling techniques.

[0131] S31) Given any initial point that satisfies the following conditions :

[0132]

[0133] In the formula, Let represent the initial values ​​chosen for the second-order cone programming problem; Indicates a convex cone internal; , These represent the optimization variables. dual variables The dimension; express , It is non-negative.

[0134] Define parameters Nonlinear KKT condition equations for perturbations, all The resulting solution set defines the central path of the HSD problem. The solution corresponding to this time is the optimal solution to the original problem; the nonlinear KKT condition equation is expressed as:

[0135]

[0136] Indicate the optimization variables of the primal problem and the dual problem. The value of the next iteration; , The self-dual variable represents the variable at the th... The next iteration value. Represents the NT scaling matrix; NT scaling matrix It is a linear transformation, which will and Mapped to the same point Make ; Represents a unit vector; operators It is defined within the Euclidean-Jordanian algebraic space; Indicates the duality gap at the th The measure of the next iteration , express The total number of cones;

[0137] Gradually decrease in each iteration Applying Newton's method to approximate equation (14), the solution variables are updated as follows:

[0138]

[0139] In the formula, This represents the step size obtained through line search; It follows the direction of Newton's search;

[0140] S32) In each iteration of the predict-correction original-dual path tracking interior point method, the final search direction is obtained by solving the sparse KKT linear system, wherein the coefficient matrix of the sparse KKT linear system is... The sparse KKT linear system, in the form of a block-based sparse symmetric matrix, is represented as follows:

[0141]

[0142] In the formula, This represents the coefficient matrix of the original linear system; The right-hand term; It is a regularization parameter with a small value. , It is the identity matrix; , For regularization terms; express; This represents the solution to a linear system.

[0143] S4) The permutation matrix is ​​calculated using the nested partitioning ND algorithm. And based on the direct method of LDL matrix decomposition, the sparse KKT linear system is solved, and the system coefficient matrix of the sparse KKT linear system is... Rearrange;

[0144] In this embodiment, the nested partitioning ND algorithm is used to calculate the permutation matrix. Specifically, it includes the following steps:

[0145] S41) The system coefficient matrix Using the corresponding adjacency graph as input, a vertex separator set is selected to divide the adjacency graph into two approximately equal-sized subgraphs. With the partition set ;

[0146] S42), recursively execute step S41), respectively on the subgraph The subgraphs are segmented until their size is less than a preset threshold. Subgraphs smaller than the threshold are then rearranged using the approximate minimum degree AMD algorithm to reduce filler elements.

[0147] S43), according to the subgraph Separated Sets Renumber the vertices of the adjacency graph in the following order to generate a permutation matrix. .

[0148] like Figure 2 As shown, Figure 2 (a) shows the sparse matrix after rearrangement using the ND algorithm. , Figure 2 (b) shows the decomposition The obtained matrix The block format, Figure 2 (c) shows the structure of the elimination tree after decomposition. Subtree 1 in the elimination tree corresponds to... and Subtree 2 corresponds to and tree trunk corresponds The block elimination tree reflects The computational dependencies of each block are described, and since triangular forward substitution and numerical decomposition have the same dependencies, the elimination tree also describes the dependencies in solving the trigonometric system. Figure 2The nodes in (c) are divided into two categories: nodes 1, 2, 4, and 5 in the second level are called subtree nodes; nodes 3 and 6 in the first level and node 7 in the 0th level are called trunk nodes.

[0149] In this embodiment, the direct method for solving sparse KKT linear systems based on LDL matrix factorization includes three stages: offline preprocessing, online numerical decomposition, and online solution.

[0150] In the offline preprocessing stage, a permutation matrix is ​​sought. To the system coefficient matrix The rows and columns are rearranged to reduce memory utilization and computational complexity in the online numerical decomposition and online solution stages; additionally, the sparsity pattern of the decomposition factors is determined by performing symbolic decomposition and the necessary memory space is allocated in advance, and the rearranged matrix is ​​calculated in the online numerical decomposition stage. The factorization factors L and D are such that:

[0151]

[0152] In the formula, the decomposition factor , These are triangular matrices and diagonal matrices, respectively. This represents the rearranged coefficient matrix; This represents the coefficient matrix of the original linear system; This represents the rearranged coefficient matrix; This represents the permutation matrix used for rearrangement;

[0153] The online solution phase obtains the solution by solving a series of equations whose computational cost is less than a cost threshold. ,Right now:

[0154]

[0155] In the formula, It is an intermediate variable in the calculation process. Represents the final solution of a linear system; Represents the right-hand side term of the original linear system; This represents the triangular and diagonal matrices after decomposition.

[0156] S5), based on rearrangement matrix Elimination tree structure The LDL matrix decomposition task is divided into subtree decomposition tasks and trunk decomposition tasks, and the subtree decomposition tasks and trunk decomposition tasks are statically allocated to different computing cores of a multi-core CPU; specifically, the following steps are included:

[0157] S51), Rearrangement Matrix The elimination tree structure is then arranged in post-order tracing.

[0158] S52) Divide the nodes of the elimination tree after the post-order arrangement into subtree nodes and trunk nodes. The subtree nodes constitute independent subtree tasks, and the trunk nodes constitute trunk tasks that depend on the results of the subtree tasks.

[0159] S53) Perform two rounds of ND decomposition, in which the elimination tree is divided into 4 independent subtree tasks and 3 trunk tasks. The difference in the number of nodes in each subtree task is less than a preset ratio, thereby achieving load balancing of multi-core CPUs.

[0160] In this embodiment, Figure 3 (a) gives a size of KKT example matrix Its decomposition factor L and eliminate trees Figure 3 (b) and Figure 3 As shown in (c); Figure 3 The elimination tree shown in (c) was rearranged to obtain... Figure 3 (d) shows an elimination tree structure where the post-order tracing places columns with similar sparse structures in the matrix close to each other. This allows for better data locality in column-based LDL decomposition computation, and such a data structure is cache-friendly; compared to the matrix In contrast, due to the system coefficient matrix decomposition... L The new non-zero elements introduced in the diagram are called "fill" elements, which are represented by red dots in the diagram.

[0161] In this embodiment, each column is treated as a computational unit during matrix operations. To facilitate efficient processor access and computation, sparse matrices... and its factors L All data is stored in a compressed sparse columnar format within the shared memory of the multi-core processor. This storage allocation method allows for easy access from any core of the multi-core processor. and L Any column.

[0162] Eliminate the matrix corresponding to each node in the tree L A certain column, for example, node 1 corresponds to L The first column. Eliminating the connections between nodes in the tree reflects the dependencies between columns during the calculation process; specifically, that is, the nodes... The calculation of the corresponding column must be completed before the node calculation is completed. The columns corresponding to all child nodes, for example in Figure 3 In the elimination tree shown in (d), the calculations must be completed first. L The calculation of column 8 can only begin after columns 6 and 7 have been completed.

[0163] like Figure 3 As shown in (d), nodes in the elimination tree can be further divided into "subtree" nodes and "trunk" nodes. In a post-order elimination tree, all nodes belonging to a subtree or the trunk are numbered consecutively, for example, in... Figure 3 In (d), the elimination tree is divided into a trunk structure (nodes 11 and 12) and two subtree structures (subtree 1 contains nodes 1, 2, 3, 4, and 5, and subtree 2 contains nodes 6, 7, 8, 9, and 10).

[0164] The computation of subtree or trunk nodes is treated as an independent computational subtask in this paper. For subtree tasks, it is clear that they are computationally independent because no node in any subtree is a child or parent node of any other node in a different subtree. This allows subtree tasks to be executed in parallel on different threads. For example, in... Figure 3 In (d), the two threads running on two separate cores can simultaneously compute columns 1-5 in subtree 1 and columns 6-10 in subtree 2. For the trunk task, its computation appears to depend on the results of all its subtrees and can only be executed serially. However, there are also parallelizable parts within the trunk task. In fact, the computation of the trunk columns depends on the previously computed columns in different subtrees and the trunk itself. Therefore, the trunk columns can be updated in parallel first using the results from different subtrees, and then the remaining computations can be executed serially.

[0165] S6) Based on the Parallel Real-Time Interior Point Method (PIPM), each core executes the subtree decomposition task in parallel, independently updates the trunk nodes through a submatrix replication mechanism, and then a designated core executes the trunk decomposition task to complete the process. break down;

[0166] This embodiment employs a column-block-based left-look decomposition algorithm to perform LDL matrix decomposition in parallel, specifically including the following steps:

[0167] S61) For each computing core, for its assigned subtree task, perform column block decomposition: traverse the corresponding column blocks of the subtree column by column, and initialize the column data as follows: For the corresponding column, update the current column with the calculated column based on the elimination tree dependency relationship, and then perform diagonal division on the non-zero elements of the current column;

[0168] S62) After the subtree decomposition is completed, each computing core performs column block updates in parallel: the subtree decomposition results are used to update the corresponding columns of the trunk nodes that the subtree depends on.

[0169] S63) For trunk tasks, each computing core performs column block decomposition based on the updated trunk column data. If multiple cores are updating the same trunk node, the process is handled locally. The submatrix replica stores the update results, which are ultimately merged by the core responsible for trunk decomposition.

[0170] S7) Based on the decomposition factor of the LDL matrix, each computational core is based on its local... The submatrix replicas perform forward substitution, diagonal scaling, and backward substitution operations in parallel to obtain the solution of the sparse KKT linear system, thereby completing the iteration of the original-dual interior-point method and finally outputting the online optimized trajectory of the aircraft; specifically, the following steps are included:

[0171] S71) Divide the right-hand term of the sparse KKT linear system into segments according to the elimination tree nodes, allocate corresponding segmented data to each computational core, and initialize the local right-hand term.

[0172] S72) All computing cores execute Lsolve operations in parallel to solve the problem. Segment the right-hand side terms to obtain intermediate variables. ;

[0173] S73) Perform a multi-core synchronization, where each computing core updates the shared trunk node with its own Lsolve result. value;

[0174] S74) All computing cores execute the Dsolve operation in parallel to solve the problem. To obtain intermediate variables ;

[0175] S75) All computing cores execute the LTsolve operation in parallel to solve the problem. ; Obtain intermediate variables ;

[0176] S76), combining all computing cores Segmented data, processed by permutation matrix The inverse transform yields the final solution of the sparse KKT linear system. .

[0177] Figure 4 The diagram illustrates the elimination tree after two rounds of ND partitioning and the sub-elimination tree obtained by copying the trunk node, containing four subtrees and three trunks. The existence of four independent subtrees in the elimination tree indicates that four threads can be used to accelerate the decomposition process. Figure 4As shown, through the trunk node replication operation, the original elimination tree is divided into four independent sub-elimination trees, each assigned to one of four threads. Concurrent updates of trunk 3 by subtrees 1 and 2 can lead to write contention on trunk 3. To avoid the performance degradation caused by using synchronization mechanisms such as mutexes, two copies of trunk 3 are created, allowing subtrees 1 and 2 to update their respective copies independently, and then the results are aggregated. Based on this strategy, the elimination tree is divided into four sub-elimination trees, each corresponding to one computation thread. Figure 4 As shown, the data of root node 7 is copied in all child elimination trees. In the second level of the elimination tree (counting from top to bottom), the data of node 3 is copied to threads 1 and 2, while the data of node 6 is copied to threads 3 and 4. In the third level, the data of all nodes is assigned to a specific thread.

[0178] Figure 5 Decomposition factor after ND partitioning of the wheel L and the submatrix obtained by copying the data from the trunk columns, Figure 5 The left half gives the decomposition factor obtained based on two rounds of ND partitioning. L It is a Block-based lower triangular matrix. Figure 4 The elimination tree partitioning method shown here copies the column data corresponding to the trunk node, and the original The matrix is ​​divided into 4 groups with Block-based submatrices (denoted as respectively) Each thread owns one of the submatrices, such as Figure 5 As shown in the right half.

[0179] (1) By L Dividing it into four smaller submatrices results in a larger computational load. L The computational task can be decomposed into four threads for execution, thereby achieving parallel acceleration of matrix decomposition; (2) Each thread only needs to perform decomposition and update operations on its corresponding submatrix, thus effectively avoiding memory write conflicts that may be caused by parallel update operations; (3) Compared with the original Matrix operations, where each submatrix has a smaller dimension and fewer non-zero elements, are more effective. In embedded multi-core processors, each thread corresponds to a dedicated computational core. Operations on smaller submatrices improve data locality and are more conducive to data being cached by the processor cores, as embedded processors typically have limited cache capacity.

[0180] Figure 6 This is a schematic diagram illustrating the parallel LDL matrix factorization algorithm used in this embodiment. The algorithm is based on... Figure 4The sub-elimination tree structure is shown, and calculations are performed in ascending order. First, data initialization is performed. Based on the node data copying method, for the trunk node, only the column corresponding to the submatrix with the smallest thread number is used. Initialize the columns in the first thread, while columns in other threads are initialized to 0. For subtree nodes, all columns of the submatrix are initialized to 0. initialization.

[0181] After initialization, the third-level calculations are performed first. Each of the four threads performs column block decomposition operations on its respective subtree nodes (1, 2, 4, 5). After the decomposition calculations are complete, column block update operations are immediately performed on all trunk nodes. Next, the second-level calculations are performed, with threads 1 and 2 performing the calculations. The sum is used to update the original value, i.e. Thread 3 and Thread 4 compute and replace with and and The value is then determined. After this, threads 1 and 2 perform column block decomposition on node 3, while threads 3 and 4 perform column block decomposition on node 6. Subsequently, thread 1 updates node 7 by performing a column block update operation on node 3, and thread 3 updates node 7 by performing a similar operation on node 6. Finally, the first-level computation is performed, with all four threads calculating simultaneously. The calculated sum is then used to replace the original value. Then, all threads perform column block decomposition operations on their respective nodes (nodes 7), and finally, all four threads complete their calculations to obtain their respective decomposition factors. and .

[0182] Each thread has a lower triangular factorization factor. and diagonal matrix Their data are respectively composed of matrices Decomposition factor L and D It consists of different parts. 4 threads and The data are as follows:

[0183] Thread 1:

[0184] Thread 2:

[0185] Thread 3:

[0186] Thread 4:

[0187] In step S7, the original vector is used. The right-hand side of each thread is initialized with segmented data. In the formula The subscript numbers represent the segmented data corresponding to the node numbers in the elimination tree:

[0188]

[0189] 4 threads execute in parallel Lsovle Calculation and solution simultaneously get :

[0190] ;

[0191] exist Lsolve After the calculation is complete, the four threads perform a synchronization operation. Then, thread 1 and thread 2 use the sum... renew and Thread 3 and Thread 4 use and renew and All 4 threads use and renew That is:

[0192]

[0193] All threads execute in parallel. Dsolve Calculation and solution simultaneously get :

[0194]

[0195] By combination Different data segments yield vectors Finally, the permutation matrix is ​​used. arrangement The elements are used to obtain the final solution vector. .

[0196] This embodiment employs parallel processing. The matrix factorization algorithm achieves high parallel computing performance. The computational complexity of different parallel subtasks in each thread is roughly estimated using the number of multiplication and division operations involved in the calculation process; the statistical results are given in Table 1. The results show that all four threads achieve good load balancing. Based on Table 1, the following can be concluded: Figure 7The task scheduling diagram shown illustrates that most of the computational burden can be accelerated through parallelism, confirming the effectiveness of the above design.

[0197] Table 1. Theoretical estimates of computational cost for different subtasks in parallel LDL decomposition.

[0198]

[0199] The convex subproblems obtained by the sequential convexification method all have the same problem structure. Therefore, the task can be statically allocated to 4 CPUs before the actual flight, and a task scheduling graph as shown in the figure below can be established.

[0200] Applying parallel computing methods to subsequent... Lsolve , Dsolve and LTsolve The operation can also achieve a good theoretical speedup. Using four computing CPU cores, the load on each core is evenly distributed, and the computational workload of each core is shown in Table 2.

[0201] Table 2. Rough statistics on the computational cost of different subtasks in solving the parallel triangular system.

[0202]

[0203] When using four CPU cores, Lsolve , Dsolve and LTsolve The operation can achieve a theoretical speedup of approximately 3.63.

[0204] In addition, when the number of discrete points of FOH is selected to be between 35 and 95 in this embodiment, the acceleration data of the KKT system solution measured on the onboard computer platform is shown in Table 3:

[0205] Table 3 Performance verification of the parallel solution method for the KKT system

[0206]

[0207] The ND-PIPM algorithm in this embodiment is used to solve the KKT system. Matrix decomposition and triangular system solution (Triangular, , , The parallel computing operation is compared with the serial sparse algebra SparseLDL. As can be seen from the table above, the ND-PIPM algorithm can achieve a computational speedup of 3.2-3.5 on 4 DSP computing cores, which proves the effectiveness of the algorithm proposed in this embodiment.

[0208] The embodiments and descriptions above are merely illustrative of the principles and preferred embodiments of the present invention. Various changes and modifications may be made to the present invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed.

Claims

1. A method for online trajectory optimization of aircraft based on a parallel real-time interior point method with nested partitioning, characterized in that, Includes the following steps: S1) Establish the dynamic model and constraints of the aircraft; S2) Reconstruct the online trajectory optimization problem of the aircraft thrust reduction fault into a convex optimization problem, and transform the linear trajectory optimization problem into a standard SOCP problem that can be directly solved by a convex solver; S3) The final search direction is obtained by iteratively solving a series of perturbation KKT condition equations for the primal and dual problems using the prediction-correction primal-dual path tracking interior point method based on homogeneous self-dual embedding HSD and nested partitioning NT scaling techniques. S4) The permutation matrix is ​​calculated using the nested partitioning ND algorithm. Specifically, it includes the following steps: S41) The system coefficient matrix Using the corresponding adjacency graph as input, a vertex separator set is selected to divide the adjacency graph into two approximately equal-sized subgraphs. With the partition set ; S42), recursively execute step S41), respectively on the subgraph The subgraphs are segmented until their size is less than a preset threshold. Subgraphs smaller than the threshold are then rearranged using the approximate minimum degree AMD algorithm to reduce filler elements. S43), according to the subgraph Separated Sets Renumber the vertices of the adjacency graph in the following order to generate a permutation matrix. ; Furthermore, the sparse KKT linear system is solved using a direct method based on LDL matrix decomposition, and the system coefficient matrix of the sparse KKT linear system is analyzed. Rearrange; S5), based on rearrangement matrix Elimination tree structure The LDL matrix decomposition task is divided into subtree decomposition tasks and trunk decomposition tasks, and the subtree decomposition tasks and trunk decomposition tasks are statically allocated to different computing cores of a multi-core CPU. The specific steps include the following: S51), Rearrangement Matrix The elimination tree structure is then arranged in post-order tracing. S52) Divide the nodes of the elimination tree after the post-order arrangement into subtree nodes and trunk nodes. The subtree nodes constitute independent subtree tasks, and the trunk nodes constitute trunk tasks that depend on the results of the subtree tasks. S53) Perform two rounds of ND decomposition, in which the elimination tree is divided into 4 independent subtree tasks and 3 trunk tasks. The difference in the number of nodes in each subtree task is less than a preset ratio, thereby achieving load balancing of multi-core CPUs. S6) Each core executes the subtree decomposition task in parallel, independently updating the trunk nodes through the submatrix replica mechanism, and then the designated core executes the trunk decomposition task to complete the process. break down; S7) Based on the decomposition factor of the LDL matrix, each computational core is based on its local... The submatrix copy performs forward substitution, diagonal scaling, and backward substitution operations in parallel to obtain the solution of the sparse KKT linear system, and then completes the iteration of the original dual interior point method, finally outputting the online optimized trajectory of the aircraft.

2. The online trajectory optimization method for aircraft based on the parallel real-time interior point method with nested partitioning as described in claim 1, characterized in that: In step S1), the dynamic model is: ； （1） In the formula, , , The velocity components of the spacecraft in the orbit-inertial coordinate system; , () is the first derivative of the position component of the spacecraft in the orbit-inertial coordinate system; , , These are the velocity components of the spacecraft in the orbit-inertial coordinate system. , , The first derivative; , This is the component of gravitational acceleration; For the mass of the aircraft; The magnitude of the thrust is constant. , These are the commanded pitch angle and yaw angle, respectively. This refers to the fuel consumption of an aircraft engine per unit time. This is a fuel consumption constant.

3. The online trajectory optimization method for aircraft based on the parallel real-time interior point method with nested partitioning as described in claim 2, characterized in that: In step S1), the constraints include: Terminal constraints: ； （4） In the formula, , These represent the terminal position components of the aircraft; , These represent the terminal velocity components of the aircraft; The gravitational constant of Earth; Initial state and control constraints: ； （5） In the formula, Indicates the initial time. State variables; Indicates the initial time. The state value; Indicates the initial time. Control variables; Indicates the initial time. The value of the control variable; Path constraints for the rate of change of commanded attitude angles: ； （6） In the formula, and These represent the rates of change of the commanded pitch and yaw angles, respectively. Upper limit.

4. The online trajectory optimization method for aircraft based on the parallel real-time interior point method with nested partitioning as described in claim 3, characterized in that: In step S1), the terminal constraints are relaxed and the relaxation term is penalized in the performance index. The relaxed terminal constraints are expressed as follows: ； （7） In the formula, , These are slack variables; The optimization objective is to maximize the altitude of the insertion point, and a weighted penalty is applied to all slack variables in the performance metrics. The augmented performance metric objective function is as follows. Represented as: ； （8） In the formula, These are the weighting coefficients.

5. The online trajectory optimization method for aircraft based on the parallel real-time interior point method with nested partitioning as described in claim 4, characterized in that: In step S4), the direct method for solving sparse KKT linear systems based on LDL matrix decomposition includes three stages: offline preprocessing, online numerical decomposition, and online solution. In the offline preprocessing stage, a permutation matrix is ​​sought. To the system coefficient matrix The rows and columns are rearranged; additionally, the sparsity pattern of the factorization factors is determined by performing symbolic decomposition and the necessary memory space is allocated in advance, and the rearranged matrix is ​​calculated in the online numerical decomposition stage. Decomposition factor and , so that: ； （17） In the formula, the decomposition factor and These are triangular matrices and diagonal matrices, respectively. The online solution phase obtains the solution by solving a series of equations whose computational cost is less than a cost threshold. ,Right now: ； （18） In the formula, It is an intermediate variable in the calculation process; Represents the final solution of a linear system; This represents the right-hand side term of the original linear system.

6. The online trajectory optimization method for aircraft based on the parallel real-time interior point method with nested partitioning as described in claim 5, characterized in that: In step S6), the LDL matrix decomposition is performed in parallel using a column-block-based left-look decomposition algorithm, specifically including the following steps: S61) For each computing core, for its assigned subtree task, perform column block decomposition: traverse the corresponding column blocks of the subtree column by column, and initialize the column data as follows: For the corresponding column, update the current column with the calculated column based on the elimination tree dependency relationship, and then perform diagonal division on the non-zero elements of the current column; S62) After the subtree decomposition is completed, each computing core performs column block updates in parallel: the subtree decomposition results are used to update the corresponding columns of the trunk nodes that it depends on; S63) For trunk tasks, each computing core performs column block decomposition based on the updated trunk column data. If multiple cores are updating the same trunk node, the process is handled locally. The submatrix replica stores the update results, which are ultimately merged by the core responsible for trunk decomposition.

7. The online trajectory optimization method for aircraft based on the parallel real-time interior point method with nested partitioning as described in claim 6, characterized in that: In step S7), based on the decomposition factors of the LDL matrix, each computational core is based on its local... The submatrix replicas perform forward substitution, diagonal scaling, and backward substitution operations in parallel to obtain the solution of the sparse KKT linear system, thereby completing the iteration of the original-dual interior-point method and finally outputting the online optimized trajectory of the aircraft; specifically, the following steps are included: S71) Divide the right-hand term of the sparse KKT linear system into segments according to the elimination tree nodes, allocate corresponding segmented data to each computational core, and initialize the local right-hand term. S72) All computing cores execute Lsolve operations in parallel to solve the problem. Segment the right-hand side to obtain intermediate variables. ; S73) Perform a multi-core synchronization, where each computing core updates the shared trunk node with its own Lsolve result. value; S74) All computing cores execute the Dsolve operation in parallel to solve the problem. To obtain intermediate variables ; S75) All computing cores execute the LTsolve operation in parallel to solve the problem. ; Obtain intermediate variables ; S76), combining all computing cores Segmented data, processed by permutation matrix The inverse transform yields the final solution of the sparse KKT linear system.

Images