PCE and lu-nscp based robust trajectory fast optimization method for rbcc vehicle ascending phase

Abstract

The application particularly relates to a RBCC aircraft ascending segment robust trajectory fast optimization method based on PCE and LU-NSCP. First, considering initial state and dynamic parameter uncertainty, a RBCC aircraft ascending segment robust trajectory optimization problem model is constructed; second, a chaos polynomial expansion method is adopted to convert the robust trajectory optimization problem into a high-dimensional deterministic trajectory optimization problem; then, an improved Newton dimension reduction sequence convex optimization method is adopted for solution, and through Newton iteration, the extended state variable and dynamic constraint are eliminated, so that the problem scale is greatly reduced; in the convex subproblem construction process, an improved strategy based on LU decomposition is introduced to replace the direct inversion of a large Jacobian matrix, so that the calculation efficiency is further improved; finally, the robust optimal trajectory is obtained by iteratively solving the convex programming subproblem. The application can effectively deal with the strong nonlinearity and multi-dimensional uncertainty of the RBCC ascending segment, and realizes fast solution while guaranteeing the robustness and reliability of the trajectory.

Description

Technical Field

[0001] This invention relates to the field of aircraft trajectory optimization technology, specifically to a robust and rapid trajectory optimization method for the ascent phase of an RBCC aircraft based on PCE and LU-NSCP. Background Technology

[0002] Ascent phase trajectory optimization is one of the core technologies for ensuring the successful execution of flight missions by RBCC-powered aircraft. However, the strong coupling between the RBCC propulsion system and its state, coupled with the characteristics of high-speed cross-domain flight, makes the uncertainties of initial conditions and models particularly significant in actual flight. Traditional methods typically only focus on improving trajectory performance under nominal conditions, without considering the robustness and reliability of the trajectory under the influence of uncertainties. Furthermore, the increasing frequency of missions places higher demands on rapid response, necessitating the development of robust trajectory optimization methods capable of quickly generating trajectories. Therefore, it is necessary to consider the impact of uncertainties in trajectory optimization, taking into account the characteristics of RBCC ascent phase flight, and to study efficient robust trajectory optimization methods, thereby balancing trajectory optimality, reliability, and robustness while also ensuring rapid solution.

[0003] Typical robust trajectory optimization problem-solving strategies involve two key techniques: uncertainty quantification and deterministic trajectory optimization. Among these, the Polynomial Chaos Expansion (PCE) method has become one of the mainstream methods for handling uncertainty quantification problems due to its low computational cost and considerable accuracy. Direct methods, including pseudospectral methods, heuristic algorithms, and sequential convex optimization, are widely used in solving deterministic trajectory optimization problems because they easily handle the complex constraints of the ascent phase. However, although robust trajectory optimization methods based on PCE and direct methods have been gradually applied to robust trajectory optimization solutions in scenarios such as Mars entry and gliding reentry, research on strongly nonlinear and coupled scenarios such as the ascent phase of RBCC spacecraft remains limited. Some researchers have only conducted studies considering the uncertainty of dynamic parameters, failing to incorporate the uncertainty of the initial state and neglecting the robustness of performance indicators. Furthermore, as the dimension of the system input random variables increases, the dimension of the PCE state expansion grows exponentially, significantly increasing the number of discrete state nodes. This results in high computational costs for robust trajectory optimization methods combining PCE and direct methods, a problem that remains unsolved. It is worth noting that the Newton dimensionality reduction convex optimization method can eliminate the state variables in the problem, which is expected to become an effective means to deal with the PCE dimension explosion problem. However, at the same time, in the context of the explosion of problem size dimensions, how to efficiently handle the matrix inversion operation in the Newton dimensionality reduction convex optimization still needs to be studied.

[0004] It should be noted that the information disclosed in the background section above is only used to enhance the understanding of the background of the present invention, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0005] To address the challenges posed by the numerous uncertainties during the ascent phase of RBCC (Reverse Radiation Computational Carrier) vehicles and the highly nonlinear and tightly coupled nature of trajectory optimization problems, existing robust trajectory optimization methods suffer from dimensionality explosion due to chaotic polynomial expansion, resulting in low computational efficiency. Furthermore, the inversion of large matrices in Newtonian dimensionality reduction convex optimization is excessively time-consuming, making it difficult to balance trajectory robustness with fast solution processing. This invention provides a fast robust trajectory optimization method for the ascent phase of RBCC vehicles based on PCE (Potentially Parallel Equation of Circular Curves) and LU-NSCP (Low-Number of Circular Curves).

[0006] Other features and advantages of the invention will become apparent from the following detailed description, or may be learned in part by practice of the invention.

[0007] According to a first aspect of the present invention, a method for rapid optimization of robust trajectory during the ascent phase of an RBCC vehicle based on PCE and LU-NSCP is provided, the method comprising: Step S1: Construct a robust trajectory optimization model for the ascent phase of a rocket-based combined cyclic RBCC powered aircraft, considering initial state uncertainties and dynamic parameter uncertainties; Step S2: Using the chaotic polynomial expansion (PCE) method based on the stochastic response surface methodology, the uncertainty in the robust trajectory optimization problem model is quantified, transforming it into a high-dimensional deterministic trajectory optimization problem involving extended state variables and control variables. ; Step S3: Apply the improved Newton dimensionality reduction sequential convex optimization method to the high-dimensional deterministic trajectory optimization problem. The solution is performed to obtain the robust optimal trajectory; specifically, it includes: Step S31: Solve the high-dimensional deterministic trajectory optimization problem Discretize the problem and use Newton's iteration method to eliminate the constraints of the extended state variables and dynamic equations, resulting in a dimension-reduced nonconvex optimization problem containing only control variables and normalized time. ; Step S32: For the dimension-reduced nonconvex optimization problem By performing convexification, it is transformed into a convex programming subproblem that can be solved efficiently. In the process of constructing the convex programming subproblem, an improved strategy based on LU decomposition is adopted to replace the direct inversion of the Jacobian matrix; Step S33: Iteratively solve the convex programming subproblem The trajectory is updated until the convergence condition is met, and the robust optimal trajectory is output.

[0008] In some exemplary embodiments, the initial state uncertainty and the dynamic parameter uncertainty are specifically as follows: (6) In the formula, Indicates the initial time; The initial state mean vector; This represents the initial state deviation. These are dynamic parameters; Represents a column vector where all elements are 1; This is the ratio of the deviation of the dynamic parameter to its nominal value; " indicates that two vectors are multiplied element by element; Let be the nominal value of the dynamic parameter; and Let them be independent random variables that all follow a normal distribution. For ease of representation, let them be denoted as... .

[0009] In some exemplary embodiments, the objective function of the robust trajectory optimization problem model is to minimize the negative value of the terminal residual mass, and a penalty term for the standard deviation of the terminal state is added to the objective function; the constraints include reliability process constraints based on the mean and standard deviation, terminal state mean constraints, and control variable boundary constraints; specifically expressed as: (12) in, For state vectors, For control vectors, It is a time variable; and To control the upper and lower boundaries of the quantity; , and These are the weighting coefficients to be designed; For stochastic dynamic equations containing random variables, For process constraints that include random variables, For terminal constraints; The negative value represents the remaining mass at the terminal. The mean, Standard deviation; These are the initial state constraints.

[0010] In some exemplary embodiments, the high-dimensional deterministic trajectory optimization problem Specifically, it is expressed as: (18) in, To expand the state variables, To control the quantity; Represents a sequence of state differential equations. For PCE projection matrix, To extend the state dynamics equations; This represents the average performance metric. The standard deviation of the performance index The standard deviation of the terminal state; For the first The mean of the process constraints For the first The standard deviation of the process constraints.

[0011] In some exemplary embodiments, step S31, which utilizes Newton's iteration method to eliminate the constraints of the extended state variables and dynamic equations, specifically includes: (1) Discretize the continuous extended state differential equation and express the discretized dynamic equation in the form of a nonlinear equation system; Take the number of discrete points as Let the discrete time point be... ( ).make , Then the discrete extended state sequence can be represented as The control sequence can be represented as The integral-type discrete form of the extended state differential equation shown in equation (18) is expressed as: (19) In the formula: ; For time-normalized variables; ,in For the constant integral matrix that is related to the discretization method, the Chebyshev pseudospectral discretization scheme is adopted; The constant matrix related to the sample is shown in Equation (20); Then it is a vector formed by the values ​​of the state equation at discrete points, as shown in equation (21); (20) (twenty one) Similarly, by substituting the discrete extended state and control into the process, terminal, and control constraints of equation (18), we can obtain its discrete form as follows: (twenty two) (twenty three) (twenty four) question The performance metrics are only related to the extended state at the terminal moment. For ease of representation, the following does not require additional discretization: (25) (2) Solve the nonlinear equation system by Newton iteration to obtain an expression for the extended state variables in explicit terms of the control variables; Building upon discrete problems, further utilizing Newton's iterative method for nonlinear equation systems, using... Indicating the number of iterations, the iterative solution of equation (19) can be expressed as: (26) In the formula, The inverse of the Jacobian matrix shown in equation (27); (27) (28) The right side of equation (26) can be simplified as follows: This formula explicitly gives the expression for the state quantity in terms of the control quantity; (3) Substitute equation (26) into the discretized constraints and performance indicators to eliminate state variables and achieve dimensionality reduction.

[0012] In some exemplary embodiments, the dimensionality reduction nonconvex optimization problem Specifically, it is expressed as: definition And satisfy ,Right now It can be considered as a discrete extended state sequence Extract Time-based expansion of state variables The operator; will Its introduction problem In the discrete form, the problem can be obtained by including only the control variable and the time-normalized variable. as follows: (29)

[0013] In some exemplary embodiments, the convexity process in step S32 includes: (1) The non-convex dynamic equations and process constraints are linearized using Taylor expansion. To solve the problem using convex optimization methods Furthermore, it is necessary to convexize the non-convex terms in the problem. Since both the performance index and constraints include extended states, the relationship between the extended states and the control sequence shown in equation (26) is first linearized; in the reference trajectory Taking the Taylor expansion of equation (26) above, we have: (30) In the formula, , and Represent equation (19) and its counterpart respectively. and The derivative of The value at that location, and The specific expressions are shown in equations (31) and (32); (31) (32) in, (33) For process constraints, since they contain nonlinear terms related to the state, they need to be linearized using Taylor expansion, as follows: (36) in, (37) (38) (2) Introduce second-order cone constraints to make the terminal state standard deviation term in the performance index convex. In the question Since the mean of the performance index and the terminal constraints only contain linear terms of the extended state, convexification can be achieved by directly substituting equation (30). However, for the performance index and the standard deviation of the terminal state, which are composed of the L2 norm of the extended state, a second-order cone constraint can be introduced to convexify them, i.e., let: (35) In the formula, and If the variable to be optimized is non-negative, adding it to the performance metric will complete the convexity of the performance metric. (3) Add trust region constraints to ensure the rationality of the linearization approximation near the reference trajectory. To ensure the boundedness of the problem, the following trust region constraint is added: (34) In the formula, and Indicates the size of the trust region.

[0014] In some exemplary embodiments, the solvable convex programming subproblem Specifically, it is expressed as: (39)

[0015] In some exemplary embodiments, the improved strategy based on LU decomposition in step S32 specifically includes: For the form in Newton's iteration process, such as A system of linear equations, in which For Jacobian matrices, For multiple right-hand terms, first, the Jacobian matrix is... Perform LU decomposition This yields the lower triangular matrix. and upper triangular matrix ; Solving the problem using forward substitution and backward substitution respectively. and This is to obtain solutions corresponding to multiple right-hand terms, thereby avoiding the need to modify the Jacobian matrix. The direct inverse of .

[0016] In some exemplary embodiments, the convergence condition for the iterative solution in step S33 is that the norm of the difference between the control variable sequences obtained from two adjacent iterations is less than a preset threshold.

[0017] The robust trajectory optimization method for the ascent phase of an RBCC vehicle based on PCE and LU-NSCP provided in the embodiments of the present invention has the following advantages compared with the prior art: (1) High robustness and reliability: This invention fully considers various uncertainties faced by the RBCC aircraft during the ascent phase, such as the initial state and dynamic parameters. Uncertainty is quantified through chaotic polynomial expansion, and reliability constraints based on the mean and standard deviation, as well as terminal state penalty terms, are introduced into the trajectory optimization problem, which can effectively reduce the sensitivity of the trajectory to uncertainty. Monte Carlo simulation results show that, compared with traditional deterministic methods, the probability of process constraint violation of the trajectory obtained by this invention is reduced from 73.53% to 0.06%, which significantly improves the robustness and reliability of the flight trajectory.

[0018] (2) High computational efficiency, solving the "curse of dimensionality": To address the problem of the explosion of state variable dimensionality caused by chaotic polynomial expansion, this invention innovatively introduces the Newton dimensionality reduction sequential convex optimization method. By explicitly eliminating the extended state variables and dynamic equation constraints through Newton iteration, the number of variables to be optimized is significantly reduced. Compared with the traditional sequential convex optimization method, this invention significantly reduces the solution time while ensuring the solution accuracy.

[0019] (3) Further improves the construction efficiency of high-dimensional problems: This invention designs an improved strategy based on LU decomposition. In the process of constructing convex programming subproblems, LU decomposition is used to replace the direct inversion operation of large Jacobian matrices. Since LU decomposition has a smaller constant factor in actual calculations, it can significantly reduce the time and storage overhead of high-dimensional matrix operations, further improving the adaptability and computational efficiency of the method to high-dimensional robust trajectory optimization problems.

[0020] (4) Strong applicability: The method proposed in this invention is not only applicable to the trajectory optimization of the ascent phase of RBCC powered aircraft, but can also be extended to the trajectory design of other hypersonic aircraft with strong nonlinearity and various uncertainties, and has good versatility and engineering application value.

[0021] It should be understood that the above general description and the following detailed description are exemplary and explanatory only, and are not intended to limit the invention. Attached Figure Description

[0022] The accompanying drawings, which are incorporated in and constitute a part of this specification, illustrate embodiments consistent with the invention and, together with the description, serve to explain the principles of the invention. It is obvious that the drawings described below are merely some embodiments of the invention, and those skilled in the art can obtain other drawings based on these drawings without any inventive effort.

[0023] Figure 1 This is a schematic diagram of the method steps of the present invention; Figure 2 This is a schematic diagram of a specific embodiment of the present invention; Figure 3 The figures show the simulation results of the DO, RBO, and RO trajectories in an embodiment of the present invention. Figure 4 The figure shows the simulation results of the time consumed in a single iteration of solving the RO4 problem according to an embodiment of the present invention. Detailed Implementation

[0024] Exemplary embodiments will now be described more fully with reference to the accompanying drawings. However, these exemplary embodiments can be implemented in many forms and should not be construed as limited to the examples set forth herein; rather, they are provided so that the invention will be more comprehensive and complete, and will fully convey the concept of the exemplary embodiments to those skilled in the art. The described features, structures, or characteristics may be combined in any suitable manner in one or more embodiments.

[0025] Furthermore, the accompanying drawings are merely illustrative of the invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and therefore repeated descriptions of them will be omitted. Some block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. These functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor devices and / or microcontroller devices.

[0026] Based on previous research findings, "A Robust Trajectory Optimization Method for Aircraft Based on SRSM and NIPC", this invention conducts research on corresponding rapid optimization methods to meet the robust trajectory design requirements of RBCC-powered aircraft during the ascent phase.

[0027] To address the trajectory design problem of the ascending phase of an RBCC (Remotely Benchmarked Carrier) aircraft with uncertain parameters and initial states, this invention proposes a robust trajectory optimization method based on chaotic polynomial expansion and improved Newtonian dimensionality reduction convex optimization. First, considering various complex uncertainties, a robust trajectory optimization problem model for the ascending phase of an RBCC-powered aircraft is constructed. Second, an uncertainty quantification algorithm based on stochastic response surfaces and chaotic polynomials is designed to transform the robust trajectory optimization model into a deterministic trajectory optimization problem with increased dimension. Furthermore, based on pseudospectral discretization and Newtonian dimensionality reduction convex optimization, the problem is transformed into a convex programming subproblem by eliminating the extended state variables, avoiding the inefficiency of directly optimizing extended dimension problems. By introducing an improved strategy based on LU decomposition, the speed of constructing large-scale problems using Newtonian dimensionality reduction convex optimization is enhanced, thus achieving a rapid solution to the robust trajectory optimization problem of the ascending phase.

[0028] refer to Figure 1 As shown, the specific steps may include: Step S1: Construct a robust trajectory optimization model for the ascent phase of a rocket-based combined cyclic RBCC powered aircraft, considering initial state uncertainties and dynamic parameter uncertainties; Step S2: The uncertainty in the robust trajectory optimization problem model is quantified by using the chaotic polynomial expansion PCE method based on the stochastic response surface method, and it is transformed into a high-dimensional deterministic trajectory optimization problem with respect to the extended state variables and control variables. Step S3: Solve the high-dimensional deterministic trajectory optimization problem using an improved Newtonian dimensionality reduction sequence convex optimization method, specifically including: Step S31: Discretize the high-dimensional deterministic trajectory optimization problem and use Newton's iteration method to eliminate the extended state variables and dynamic equation constraints to obtain a dimension-reduced non-convex optimization problem containing only control variables and normalized time. Step S32: The dimension-reduced non-convex optimization problem is transformed into a convex programming subproblem that can be solved efficiently. In the process of constructing the convex programming subproblem, an improved strategy based on LU decomposition is used to replace the direct inversion of the Jacobian matrix. Step S33: Iteratively solve the convex programming subproblem and update the trajectory until the convergence condition is met, and output the robust optimal trajectory.

[0029] The steps in this exemplary embodiment will now be described in more detail with reference to the accompanying drawings and examples.

[0030] Step 1: Building a Robust Trajectory Optimization Problem Model Step 1.1: Center of Mass Motion Model Assuming the Earth is a non-rotating sphere and the spacecraft is considered a controllable point mass, the equations of motion for the center of mass in the longitudinal plane during the ascent phase of the RBCC spacecraft are established as follows: (1) In the formula, For state variables; For height; For speed; The trajectory inclination angle; For quality; For thrust; and These are lift and drag, respectively. For the angle of attack; It is the acceleration due to gravity; It is the distance from the Earth's center; This refers to the engine fuel flow rate per second.

[0031] Aerodynamics and It can be represented as: (2) In the formula, For dynamic pressure, Atmospheric density; For aircraft reference area; , The aerodynamic lift-drag coefficient can be expressed as the Mach number. and attack angle The function.

[0032] The RBCC-powered aircraft is propelled by both a rocket engine and a ramjet engine. The actual combination of these two engines is quite complex. To simplify the model, this invention expresses the thrust and fuel flow rate per second of the RBCC combined propulsion system as follows: (3) In the formula, and These represent the thrust of the rocket and the ramjet engine, respectively. and These represent the fuel flow rate per second for rockets and ramjet engines, respectively.

[0033] Rocket engine thrust can be expressed as: (4) In the formula, The specific impulse of a rocket engine can be characterized as a function of altitude. .

[0034] The thrust and fuel flow rate per second of a ramjet engine can be expressed as: (5) In the formula, The specific impulse of a ramjet engine can be characterized as a function of the Mach number. ; This represents the area of ​​the stamping flow path; The flow coefficient can be expressed as a function of Mach number and angle of attack. ; The fuel equivalence ratio is 1 in this invention; The mass ratio of fuel to air under complete combustion conditions is taken as 0.068 in this invention. It is evident that the introduction of the ramjet engine creates a strong coupling relationship between thrust magnitude and states and control variables such as altitude, speed, and angle of attack, further enhancing the nonlinearity of the dynamic model.

[0035] Step 1.2: Uncertainty Model During actual flight, aircraft face various uncertainties, primarily including input uncertainties and inherent uncertainties. Input uncertainties refer to the error between the aircraft's initial state and nominal state, while inherent uncertainties refer to errors caused by inaccurate modeling of atmospheric density, aerodynamic forces, etc. These uncertainties are generally modeled as random deviations relative to their nominal states and described using a certain probability distribution model. The uncertainties considered in this invention mainly include initial state, atmospheric density, aerodynamic coefficients, and thrust, and the model can be established as follows: (6) In the formula, Indicates the initial time; The initial state mean vector; This represents the initial state deviation. These are dynamic parameters; Represents a column vector where all elements are 1; This is the ratio of the deviation of the dynamic parameter to its nominal value; " indicates that two vectors are multiplied element by element; Let be the nominal value of the dynamic parameter; and Let them be independent random variables that all follow a normal distribution. For ease of representation, let them be denoted as... .

[0036] Step 1.3: Constraint Model For the design of the ascent trajectory of an RBCC aircraft, the main considerations are process, terminal, and control constraints.

[0037] (1) Process constraints To ensure flight safety, dynamic pressure must be met. Stagnation heat flux density and overload The process constraints are as follows: (7) In the formula, , and These are the maximum allowable values ​​for dynamic pressure, stagnation point heat flux density, and overload, respectively. Let be the heat flux density constant, the value of which is related to the aerodynamic shape of the aircraft. Equation (31) can be simplified as follows: .

[0038] (2) Control constraints Control parameters during the ascent phase of the RBCC spacecraft include the rocket fuel flow rate per second. and attack angle The corresponding constraints are: (8) In the formula, , and , These are the upper and lower limits of rocket fuel flow rate per second and angle of attack, respectively. If using... If we denote the control variable, then equation (32) can be simplified as follows: .

[0039] (3) Terminal constraints Terminal constraints include altitude, velocity, and trajectory inclination constraints, expressed as: (9) In the formula, For terminal time; , and These are the terminal height, velocity, and trajectory inclination constraint values, respectively. Equation (9) can be simplified as follows: .

[0040] Step 1.4: Trajectory Optimization Problem (1) Deterministic trajectory optimization model The dynamic equation shown in equation (1) can be simplified as follows: (10) In the formula, The function on the right side of the equation (25) is represented by the equal sign.

[0041] During the ascent phase of an RBCC aircraft, optimal fuel consumption is of paramount concern. Therefore, this invention selects maximizing the remaining mass at the terminal stage as the performance indicator, namely: (11) Based on the aforementioned dynamic equations, process constraints, and performance indicators, a deterministic trajectory optimization (DO) model can be established, as shown in equation (12).

[0042] (12) In the formula, These are the initial state constraints.

[0043] (2) Robust trajectory optimization model Consider the uncertainties in step 1.2 Then, the dynamic equations will be transformed into stochastic dynamic equations containing random variables, which can be expressed as: (13) Process constraints are also affected by input random variables and can be expressed as functions containing random variables. In uncertain environments, reliability process constraints (as shown in (14)) can be established based on the mean and standard deviation of the process constraints. If these constraints are met, it can be considered that there is a high probability that the process constraints will be met in actual flight. The probability of actually meeting the constraints is related to the constant. The value of is related to .

[0044] (14) In the formula, The mean, The standard deviation is denoted as .

[0045] For state endpoint constraints, after considering the influence of random variables, they are transformed into constraints on the mean values ​​of the initial and terminal states, which can be expressed as: (15) In the formula, Represents the terminal state constraints, where As a unit array, It is a zero vector.

[0046] Therefore, by combining the aforementioned stochastic dynamic equations and process constraints, a robust trajectory optimization (RO) model as shown in equation (16) can be established. Under the influence of uncertainty, the system state... Performance indicators and process constraints All are transformed into random quantities related to the input random variables. The RO problem, on the one hand, reduces the probability of violation of process constraints and improves the reliability of the trajectory by transforming process constraints into reliability constraints composed of their mean and standard deviation; on the other hand, it adds the performance index and the standard deviation of the terminal state as a penalty term to the objective function and imposes constraints on the mean of the terminal state to reduce the dispersion under random disturbances and enhance the robustness of the trajectory.

[0047] (16) In the formula, , and All are normal numbers.

[0048] In the RO problem, parameters and The robustness level of the obtained trajectory is determined. and The smaller the value, the worse the robustness of the obtained trajectory. When both parameters are 0, the RO problem degenerates into a reliability trajectory optimization (RBO) problem, as shown in equation (17). This problem only considers the reliability of process constraints and does not consider the robustness of terminal states and performance indicators.

[0049] (17) Step 2: Quantifying Uncertainty and Transforming the Problem This invention builds upon previous research findings, "A Robust Trajectory Optimization Method for Aircraft Based on SRSM and NIPC," to quantify the uncertainty of the robust trajectory optimization problem, transforming it into a deterministic trajectory optimization problem with extended state dimensions. Specifically, it introduces a Latin hypercube sampling strategy and uses the Stochastic Response Surface Method (SRSM) to quantify the uncertainty of the state vector containing uncertain parameters. The problem is expanded into a high-dimensional ordinary differential equation with respect to the coefficients of the PC polynomials. The mean and standard deviation of the states, constraints, and objective function in the Rotametry problem are then calculated using the PC polynomial coefficients, transforming the robust trajectory optimization problem for an aircraft into an equivalent high-dimensional deterministic optimization problem. This step mainly includes state variables. The invention comprises three parts: PC expansion, process constraints, and PC expansion and statistical calculation of the objective function, as detailed in "A Robust Trajectory Optimization Method for Aircraft Based on SRSM and NIPC". Here, the invention only uses the RO problem as an example, presenting its transformation result. After transformation, the RO problem can be expressed as a problem concerning the extended state. and control quantity High-dimensional deterministic optimization problem ,question The solution can be obtained using optimization algorithms.

[0050] (18) It is worth noting that, regarding the question This can be solved using a direct method. However, compared to deterministic trajectory optimization problems, The number of state variables is the same as the original. This means that with the same number of discrete points, [times]. The number of variables to be solved in the discrete parameter optimization problem is close to the original number. The time difference is several times greater, regardless of whether the pseudospectral method or the sequential convex optimization method is used to solve the problem. The process takes a long time, and this phenomenon is particularly noticeable when the dimension d of the input random variable is large.

[0051] Step 3: Solving the robust trajectory optimization problem To solve the problem efficiently This step presents an improved Newton dimensionality reduction convex optimization method. First, based on Newton iteration to handle nonlinear dynamics, it explicitly eliminates state variables and state equations from the problem to improve the solution efficiency of convex programming problems. Second, it designs a matrix inversion strategy based on LU decomposition to improve the shortcomings of Newton dimensionality reduction convex optimization in terms of long time consumption in the construction of large problems.

[0052] Step 3.1: Problem Discretization and Transformation Take the number of discrete points as Let the discrete time point be... ( ).make , Then the discrete extended state sequence can be represented as The control sequence can be represented as The integral-type discrete scheme of the extended state differential equation shown in equation (18) can be expressed as: (19) In the formula: ; For time-normalized variables; ,in For the constant integral matrix related to the discretization method, this invention adopts the Chebyshev pseudospectral discretization scheme; The constant matrix related to the sample is shown in Equation (20); Then it is a vector formed by the values ​​of the state equation at discrete points, as shown in equation (21).

[0053] (20) (twenty one) Similarly, by substituting the discrete extended state and control into the process, terminal, and control constraints of equation (18), we can obtain its discrete form as follows: (twenty two) (twenty three) (twenty four) question The performance metrics are only related to the extended state at the terminal moment. For the sake of convenience, this invention is referred to as follows: (The following text appears to be a separate, unrelated section about a different invention: "Regarding the absence of additional discretization, for the sake of convenience, this invention is referred to as:") (25) Building upon discrete problems, further utilizing Newton's iterative method for nonlinear equation systems, using... Indicating the number of iterations, the iterative solution of equation (19) can be expressed as: (26) In the formula, The inverse of the Jacobian matrix is ​​shown in expression (27).

[0054] (27) (28) The right side of equation (26) can be simplified as follows: This formula explicitly gives the expression for the state variables in terms of control variables, which can be substituted into the discretized constraints and indices to eliminate the state variables.

[0055] For ease of explanation, the present invention defines And satisfy ,Right now It can be considered as a discrete extended state sequence Extract Time-based expansion of state variables The operator. Its introduction problem In the discrete form, the problem can be obtained by including only the control variable and the time-normalized variable. as follows: (29) It is worth noting that, assuming the same number of discrete points, let the dimension of the control quantity be... Given the same number of discrete points, the problem There are a total of variables to be optimized in the data. One, and the problem The number of variables to be optimized is only One. Generally speaking. Therefore, the introduction of equation (26) significantly reduces the number of variables to be solved in the extended problem, solves the "curse of dimensionality" caused by the dimensionality expansion of the PCE state, and lays the foundation for fast solution of the RO problem. Furthermore, compared to the problem... ,question It also eliminates dynamic constraints, resulting in fewer equality constraints and easier convergence.

[0056] Step 3.2: Problem Convexification To solve the problem using convex optimization methods Furthermore, it is necessary to convexize the non-convex terms in the problem. Since both the performance index and constraints include extended states, we first linearize the relationship between the extended states and the control sequence shown in equation (26). In the reference trajectory Taking the Taylor expansion of equation (26) above, we have: (30) In the formula, , and Represent equation (19) and its counterpart respectively. and The derivative of The value at that location, and The specific expressions are shown in equations (31) and (32).

[0057] (31) (32) in, (33) To ensure the boundedness of the problem, the following trust region constraint is added: (34) In the formula, and Indicates the size of the trust region.

[0058] In the question Since the mean of the performance index and the terminal constraints only contain linear terms of the extended state, convexification can be achieved by directly substituting equation (30). For the performance index and the standard deviation of the terminal state, which are composed of the L2 norm of the extended state, convexification can be achieved by introducing a second-order cone constraint, i.e., let: (35) In the formula, and If the variable to be optimized is non-negative, adding it to the performance metrics will complete the convexity of the performance metrics.

[0059] For process constraints, since they contain nonlinear terms related to the state, they need to be linearized using Taylor expansion, as follows: (36) in, (37) (38) Thus far, the present invention has addressed the extended problem. Transformed into a convex programming problem where the variables to be optimized consist only of the control sequence and time-normalized variables. This can be solved using an efficient convex optimization solver. The solution process for the Roma problem based on Newton's dimension reduction convex optimization can be summarized as follows: In each iteration, first start with a convex programming subproblem with reduced space. Solving new control and new time normalization variables Then, the new state is recovered using the Newton iteration method given in equation (26). If the convergence condition is met. ,but If the solution is not converged, proceed to the next iteration and adjust the trust region according to certain rules.

[0060] (39) It is worth noting that in traditional sequential convex optimization algorithms, in order to ensure the feasibility of the subproblems after convexification, it is often necessary to add virtual control variables as relaxation terms for the constraints of the linearized state equations. However, the Newton-Simplified Convex Optimization Algorithm (NSCP) used in this invention not only eliminates state variables in the problem construction but also avoids the introduction of virtual control variables, enabling fast solution.

[0061] Furthermore, as shown in equation (30), the Newton dimensionality reduction convex optimization method requires adjusting the Jacobian matrix during the construction of the convex programming problem. Performing inversion operations, for robust optimization problems, the matrix... It is a size of The square array, when taken , , hour, Since inverting large matrices consumes a significant amount of time and storage space, using Newton's dimension reduction convex optimization method to solve high-dimensional Roma problems will require a substantial amount of time during the problem construction process compared to solving deterministic optimization problems.

[0062] Step 3.3: An Improved Method Based on LU Decomposition Generally, the purpose of matrix inversion is to solve a system of linear equations. In the problem studied in this invention, equation (30) can be used... This term is considered a solution to the linear system of equations (40), while the other terms in equation (30) can be considered solutions to equations with the same coefficient matrix but different right-hand side terms as the system of equations (40). Note that LU decomposition requires only one decomposition when solving linear systems of equations with multiple right-hand side terms, exhibiting high efficiency. Therefore, to further reduce the time consumption of the Newton dimensionality reduction convex optimization method in the problem construction process, this step employs an LU decomposition-based method to improve the algorithm, making it more suitable for solving RO problems.

[0063] (40) Without loss of generality, assume the matrix All principal minors are non-zero. For the matrix... LU decomposition yields: (41) Substituting equation (41) into equation (40), we get: (42) Summarized as follows: (43) Similarly, the same operation can be performed on the other terms in equation (30). Substituting the LU decomposition result back into equation (30) yields: (44) It is worth noting that although the time complexity of LU decomposition and matrix inversion is the same, LU decomposition Compared to finding the inverse The constant is smaller. This means that LU decomposition is actually faster in computation, and the efficiency improvement is very significant when the problem size is large, which can also be verified by the subsequent simulations.

[0064] The improved method based on LU decomposition proposed in this step can effectively enhance the efficiency of Newton's dimensionality reduction convex optimization method in solving large-scale problems. Furthermore, since this improved method is relatively independent, its effectiveness can be verified separately. (In the matrix...) Size is For problems involving dimensional square matrices, explicit inversion, Krylov subspace-based methods, and LU decomposition-based methods are used to calculate the expression in equation (30) and... The three items took 18.83s, 26.55s and 6.35s respectively, which proves the effectiveness of the improved method proposed in this invention.

[0065] This step proposes a method for solving robust trajectory optimization problems based on polynomial chaos expansion (PCE) and improved Newton type-sequential convex programming (LU-Newton type-Sequential Convex Programming), which is abbreviated as the LUNSCP-PCE method. This method first transforms the RO problem into a deterministic optimization problem with state expansion using the PCE method. Based on this, the problem was eliminated using Newton's iterative approach. The extended state and state equation constraints in the problem form a convex programming subproblem with a reduced space. This reduces the number of variables and constraints to be optimized in the problem, while improving the solution efficiency and convergence. Subsequently, an improved strategy based on LU decomposition was designed to further reduce the time spent constructing the problem. Finally, efficient solution of the RO problem was achieved.

[0066] Example like Figure 2 The diagram shown is a flowchart of a specific embodiment of the present invention, including the following steps: Step 1: Based on the requirements of the RBCC ascending phase task, and taking into account the uncertainties of dynamic parameters and initial state, establish a robust trajectory optimization (RO) problem model; Step 2: Using the chaotic polynomial expansion (PCE) method, the RO problem is transformed into a high-dimensional deterministic optimization problem with respect to the extended state. ; Step 3: For The problem is discretized, and Newton's dimensionality reduction strategy is applied to eliminate the extended state variables, resulting in a problem that is only concerned with the control variable and the time-normalized variable. ; Step 4: For the non-convex problem after dimensionality reduction By performing linearization and convexization, it is transformed into a convex programming subproblem that can be solved efficiently. ; Step 5: In the construction of the convex subproblem, an improved strategy based on LU decomposition is adopted to replace the direct inversion of large Jacobian matrices, so as to improve computational efficiency. Step 6: Iteratively solve the convex subproblem and update the trajectory, and determine whether it has converged. If it has converged, end the calculation and output the robust optimal trajectory; otherwise, update the reference trajectory and jump to step 3 until the convergence condition is met.

[0067] The invention will be further illustrated below with reference to simulation: Simulations were conducted using an air-launched RBCC concept aircraft. The initial mass of the aircraft was 4 tons, carrying 3.2 tons of fuel. The aircraft was horizontally launched from a carrier aircraft, and the initial and terminal constraints of the mission are shown in Table 1. The angle of attack constraint was set as follows. The rocket engine fuel flow rate per second constraint is Process constraints are set as follows: , , This invention considers 4-dimensional dynamic parameters and 4-dimensional initial state uncertainties, assuming they all follow a zero-mean normal distribution. In the RBO problem, the parameters are taken as follows: In the RO problem, the standard deviation of terminal height is added to the performance metric as a penalty term, with the parameter value set to... , , .

[0068] Table 1 Initial and Terminal State Constraints

[0069] All simulations were performed on a desktop computer equipped with an Intel Core i7-8700 3.20 GHz processor, using MATLAB 2020b as the simulation platform. All convex optimization subproblems were solved using the MOSEK solver.

[0070] To verify the effectiveness of the LUNSCP-PCE algorithm proposed in this invention, this step involves conducting a robust trajectory optimization simulation for the ascent phase of an RBCC spacecraft. The optimization results for the three problems are as follows: Figure 3 As shown.

[0071] Figure 3 The optimization results for the three problems were compared. The control variables from the optimization results were used for 3000 Monte Carlo (MC) simulations. The MC simulation results are shown in Table 2. This represents the probability of process constraint violation.

[0072] To further demonstrate the speed of the proposed LUNSCP-PCE method, this invention compares and simulates the proposed algorithm with the Newton-Non-Simplified Convex Optimization-Chaotic Polynomial Method (NSCP-PCE) and the Sequential Convex Optimization-Chaotic Polynomial Method (SCP-PCE). The time consumption of each of the three methods for solving the Roma problem in a single iteration is as follows: Figure 4 As shown, the blue, red, and yellow bars in the bar chart represent the time spent on problem construction, problem solving, and post-processing, respectively. Furthermore, the number of iterations for the three methods are 79, 68, and 68, respectively.

[0073] The simulation results above show that the NSCP-PCE method can reduce the number of variables to be designed in the robust trajectory optimization problem and shorten the solution time of the optimization problem; the LUNSCP-PCE method designed in this invention can further reduce the time required for problem construction by LU decomposition of large Jacobian matrices.

[0074] Furthermore, the above figures are merely illustrative of the processes included in the method according to exemplary embodiments of the present invention, and are not intended to be limiting. It is readily understood that the processes shown in the above figures do not indicate or limit the temporal order of these processes. Additionally, it is readily understood that these processes may be executed synchronously or asynchronously, for example, in multiple modules.

[0075] Other embodiments of the invention will readily occur to those skilled in the art upon consideration of the specification and practice of the invention herein. This application is intended to cover any variations, uses, or adaptations of the invention that follow the general principles of the invention and include common knowledge or customary techniques in the art not disclosed herein. The specification and embodiments are to be considered exemplary only, and the true scope and spirit of the invention are indicated by the claims.

[0076] It should be understood that the present invention is not limited to the precise structure described above and shown in the accompanying drawings, and various modifications and changes can be made without departing from its scope. The scope of the invention is defined only by the appended claims.

Claims

1. A robust trajectory optimization method for the ascent phase of an RBCC vehicle based on PCE and LU-NSCP, characterized in that, The method includes: Step S1: Construct a robust trajectory optimization model for the ascent phase of a rocket-based combined cyclic RBCC powered aircraft, considering initial state uncertainties and dynamic parameter uncertainties; Step S2: Using the chaotic polynomial expansion (PCE) method based on the stochastic response surface methodology, the uncertainty in the robust trajectory optimization problem model is quantified, transforming it into a high-dimensional deterministic trajectory optimization problem involving extended state variables and control variables. ; Step S3: Apply the improved Newton dimensionality reduction sequential convex optimization method to the high-dimensional deterministic trajectory optimization problem. The solution is performed to obtain the robust optimal trajectory; specifically, it includes: Step S31: Solve the high-dimensional deterministic trajectory optimization problem Discretize the problem and use Newton's iteration method to eliminate the constraints of the extended state variables and dynamic equations, resulting in a dimension-reduced nonconvex optimization problem containing only control variables and normalized time. ; Step S32: For the dimension-reduced nonconvex optimization problem By performing convexification, it is transformed into a convex programming subproblem that can be solved efficiently. In the process of constructing the convex programming subproblem, an improved strategy based on LU decomposition is adopted to replace the direct inversion of the Jacobian matrix; Step S33: Iteratively solve the convex programming subproblem The trajectory is updated until the convergence condition is met, and the robust optimal trajectory is output.

2. The method according to claim 1, characterized in that, The initial state uncertainty and the dynamic parameter uncertainty are specifically as follows: (6) In the formula, Indicates the initial time; The initial state mean vector; This represents the initial state deviation. These are dynamic parameters; Represents a column vector where all elements are 1; This is the ratio of the deviation of the dynamic parameter to its nominal value. " indicates that two vectors are multiplied element by element; Let be the nominal value of the dynamic parameter; and Let them be independent random variables that all follow a normal distribution. For ease of representation, let them be denoted as... .

3. The method according to claim 2, characterized in that, The objective function of the robust trajectory optimization problem model is to minimize the negative value of the terminal residual mass, and a penalty term for the standard deviation of the terminal state is added to the objective function; the constraints include reliability process constraints based on the mean and standard deviation, terminal state mean constraints, and control variable boundary constraints; specifically expressed as follows: (12) in, For state vectors, For control vectors, It is a time variable; and To control the upper and lower boundaries of the quantity; , and These are the weighting coefficients to be designed; For stochastic dynamic equations containing random variables, For process constraints that include random variables, For terminal constraints; The negative value represents the remaining mass at the terminal. The mean, Standard deviation; These are the initial state constraints.

4. The method according to claim 3, characterized in that, The high-dimensional deterministic trajectory optimization problem Specifically, it is expressed as: (18) in, To expand the state variables, To control the quantity; Represents a sequence of state differential equations. For PCE projection matrix, To extend the state dynamics equations; This represents the average performance metric. The standard deviation of the performance index The standard deviation of the terminal state; For the first The mean of the process constraints For the first The standard deviation of the process constraints.

5. The method according to claim 4, characterized in that, Step S31, which uses Newton's iteration method to eliminate the constraints of the extended state variables and dynamic equations, specifically includes: (1) Discretize the continuous extended state differential equation and express the discretized dynamic equation in the form of a nonlinear equation system; Take the number of discrete points as Let the discrete time point be... ( ).make , Then the discrete extended state sequence can be represented as The control sequence can be represented as The integral-type discrete form of the extended state differential equation shown in equation (18) is expressed as: (19) In the formula: ; For time-normalized variables; ,in For the constant integral matrix that is related to the discretization method, the Chebyshev pseudospectral discretization scheme is adopted; The constant matrix related to the sample is shown in Equation (20); Then it is a vector formed by the values ​​of the state equation at discrete points, as shown in equation (21); (20) (21) Similarly, by substituting the discrete extended state and control into the process, terminal, and control constraints of equation (18), we can obtain its discrete form as follows: (22) (23) (24) question The performance metrics are only related to the extended state at the terminal moment. For ease of representation, the following does not require additional discretization: (25) (2) Solve the nonlinear equation system by Newton iteration to obtain an expression for the extended state variables in explicit terms of the control variables; Building upon discrete problems, further utilizing Newton's iterative method for nonlinear equation systems, using... Indicating the number of iterations, the iterative solution of equation (19) can be expressed as: (26) In the formula, The inverse of the Jacobian matrix shown in equation (27); (27) (28) The right side of equation (26) can be simplified as follows: This formula explicitly gives the expression for the state quantity in terms of the control quantity; (3) Substitute equation (26) into the discretized constraints and performance indicators to eliminate state variables and achieve dimensionality reduction.

6. The method according to claim 5, characterized in that, The dimension reduction nonconvex optimization problem Specifically, it is expressed as: definition And satisfy ,Right now It can be considered as a discrete extended state sequence Extract Time-based expansion of state variables The operator; will Its introduction problem In the discrete form, the problem can be obtained by including only the control variable and the time-normalized variable. as follows: (29) 7. The method according to claim 6, characterized in that, The convexity process in step S32 includes: (1) The non-convex dynamic equations and process constraints are linearized using Taylor expansion. To solve the problem using convex optimization methods Furthermore, it is necessary to convexize the non-convex terms in the problem. Since both the performance index and constraints include extended states, the relationship between the extended states and the control sequence shown in equation (26) is first linearized; in the reference trajectory Taking the Taylor expansion of equation (26) above, we have: (30) In the formula, , and Represent equation (19) and its counterpart respectively. and The derivative of The value at that location, and The specific expressions are shown in equations (31) and (32); (31) (32) in, (33) For process constraints, since they contain nonlinear terms related to the state, they need to be linearized using Taylor expansion, as follows: (36) in, (37) (38) (2) Introduce second-order cone constraints to make the terminal state standard deviation term in the performance index convex. In the question Since the mean of the performance index and the terminal constraints only contain linear terms of the extended state, convexification can be achieved by directly substituting equation (30). However, for the performance index and the standard deviation of the terminal state, which are composed of the L2 norm of the extended state, a second-order cone constraint can be introduced to convexify them, i.e., let: (35) In the formula, and If the variable to be optimized is non-negative, adding it to the performance metric will complete the convexity of the performance metric. (3) Add trust region constraints to ensure the rationality of the linearization approximation near the reference trajectory. To ensure the boundedness of the problem, the following trust region constraint is added: (34) In the formula, and Indicates the size of the trust region.

8. The method according to claim 7, characterized in that, The solvable convex programming subproblem Specifically, it is expressed as: (39) 9. The method according to claim 1, characterized in that, The improved strategy based on LU decomposition in step S32 specifically includes: For the form in Newton's iteration process, such as A system of linear equations, in which For Jacobian matrices, For multiple right-hand terms, first, the Jacobian matrix is... Perform LU decomposition This yields the lower triangular matrix. and upper triangular matrix ; Solving the problem using forward substitution and backward substitution respectively. and This is to obtain solutions corresponding to multiple right-hand terms, thereby avoiding the need to modify the Jacobian matrix. The direct inverse of .

10. The method according to claim 1, characterized in that, The convergence condition for the iterative solution in step S33 is that the norm of the difference between the control variable sequences obtained from two adjacent iterations is less than a preset threshold.

Images