An aircraft stability bifurcation analysis method considering flight state constraints

Abstract

The application relates to the technical field of aircraft control, and discloses an aircraft stability bifurcation analysis method considering flight state constraints, which comprises the following steps: establishing an aircraft six-degree-of-freedom nonlinear flight dynamics model, defining a system state parameter vector and a system control parameter vector, introducing a flight state constraint equation, selecting an extension parameter, constituting an extended state variable, constructing an extended bifurcation analysis model, adopting a pseudo-arc length continuation algorithm to obtain a trimming solution set and a trimming numerical mapping relationship, obtaining a physical dynamics system equation, calculating a Jacobian matrix and extracting eigenroots, judging stability and bifurcation types according to the eigenroots, and obtaining a bifurcation curve and using the bifurcation curve for flight control. The flight state constraints are converted into equivalent differential constraints and are solved together with the flight dynamics model, so that the equilibrium points obtained through the continuation calculation satisfy the flight state constraints, and the bifurcation curve is obtained in combination with the trimming numerical mapping relationship and the eigenroots, thereby providing a basis for the determination of the dynamics boundary of a dangerous flight mode.

Description

Technical Field

[0001] This invention relates to the field of aircraft control technology, specifically to an aircraft stability bifurcation analysis method that considers flight state constraints. Background Technology

[0002] Aircraft stability analysis is used to obtain the trim characteristics and stability changes of an aircraft under different flight states and control parameters. For aircraft with nonlinear aerodynamic characteristics and multiple control parameters, the equilibrium point and equilibrium point stability will change when the angle of attack, velocity, or control parameters change. Therefore, it is necessary to analyze the trim branches and stability bifurcation of the aircraft.

[0003] Existing aircraft stability analysis typically includes linearization analysis methods based on trim points and numerical simulation analysis methods based on nonlinear flight dynamics models. Linearization analysis based on trim points linearizes the nonlinear flight dynamics model near a given equilibrium point and determines the stability near that point by calculating the eigenvalues ​​of the linearized model. This method can obtain stability information near a specified equilibrium point, but for the trim branches formed by continuous changes in control parameters at the equilibrium point and the stability changes on these branches, further analysis using methods that can continuously track the equilibrium solution is needed.

[0004] Numerical simulation analysis methods based on nonlinear flight dynamics models can calculate the motion response of an aircraft under given initial conditions and control parameters. However, when analyzing the trim points and stability changes under different control parameters, it is usually necessary to calculate each parameter separately, making it difficult to directly generate trim branches and bifurcation curves that continuously change with the control parameters.

[0005] Numerical extension and bifurcation analysis can be used to track the trim branches formed by the equilibrium point as the extension parameters change, and can determine the stability and bifurcation type of the equilibrium point based on the changes in the eigenvalues ​​of the linearized model at the equilibrium point. However, when performing trim analysis on an aircraft under specific flight conditions such as stationary level flight, the constraints of the corresponding flight conditions must also be met. When the aircraft model contains multiple control parameters, if the flight condition constraints are not combined with the extension calculation process, the calculated trim solution and stability analysis results are difficult to correspond to the preset flight condition conditions.

[0006] Therefore, it is necessary to provide a bifurcation analysis method for aircraft stability that considers flight state constraints. Under flight state constraints, the method performs balancing solution extension calculations to obtain the balancing solution set that varies with the extension parameters. Based on the characteristic roots at each equilibrium point, the stability and bifurcation type are determined, thereby obtaining the bifurcation curves and dynamic boundary conditions of the dangerous flight modes under the corresponding flight states. Summary of the Invention

[0007] To address the shortcomings of existing technologies, this invention provides a bifurcation analysis method for aircraft stability that considers flight state constraints, thus solving the problems mentioned in the background section.

[0008] To achieve the above objectives, this invention provides a method for bifurcation analysis of aircraft stability considering flight state constraints, comprising the following steps:

[0009] Establish a six-degree-of-freedom nonlinear flight dynamics model for the aircraft, and define the system state parameter vector and the system control parameter vector;

[0010] By introducing flight state constraint equations, selecting one control parameter from the system control parameter vector as an extension parameter, and incorporating the control parameters other than the extension parameter as non-extension control parameters into the system state parameter vector to form extended-dimensional state variables, the flight state constraint equations are transformed into equivalent differential constraint equations, and the equivalent differential constraint equations are combined with the six-degree-of-freedom nonlinear flight dynamics model of the aircraft to construct an extended bifurcation analysis model.

[0011] Based on the extended bifurcation analysis model, at the known initial equilibrium point, the pseudo-arc length extension algorithm is used to alternately execute the prediction step and the correction step to perform extension calculation, obtain the balancing solution set, and extract the balancing numerical mapping relationship between the non-extension control parameters and the extension parameters.

[0012] Substitute the trim numerical mapping relationship into the system control parameter vector and substitute it back into the six-degree-of-freedom nonlinear flight dynamics model of the aircraft to obtain the physical dynamic system equations that satisfy the flight state constraint equations;

[0013] Substitute each equilibrium point in the balanced solution set into the physical dynamic system equations for linearization, calculate the Jacobian matrix, and extract the eigenvalues.

[0014] Based on the characteristic roots, the stability and bifurcation type of each equilibrium point are determined, and the stability determination results are marked on the balancing curve to obtain the bifurcation curve.

[0015] Determine the dynamic boundary conditions for dangerous flight modes and apply them to flight control.

[0016] In the above technical solution, the flight state constraint equations are used to define the geometric and kinematic boundaries of the aircraft in the trim state. Transforming these flight state constraint equations into equivalent differential constraint equations, and then combining these equivalent differential constraint equations with the six-degree-of-freedom nonlinear flight dynamics model of the aircraft, allows the equilibrium point obtained through extension calculations to satisfy the corresponding flight state constraints. The equivalent differential constraint equations are used for trim point search and extension calculations, and do not represent the actual time response of the non-extension control parameters.

[0017] Furthermore, the system state parameter vector includes flight speed, angle of attack, pitch rate, pitch angle, sideslip angle, roll rate, yaw rate, and roll angle;

[0018] The system control parameter vector includes elevator deflection angle, engine throttle, aileron deflection angle, and rudder deflection angle. Based on the net external forces and net external moments acting on the aircraft, translational dynamics equations characterizing the translational motion of the aircraft's center of mass, rotational dynamics equations characterizing the rotational motion of the aircraft around its center of mass, and kinematic equations characterizing the evolution of the aircraft's spatial attitude are constructed. The translational dynamics equations, rotational dynamics equations, and kinematic equations together constitute the six-degree-of-freedom nonlinear flight dynamics model of the aircraft.

[0019] Furthermore, the construction of the extended bifurcation analysis model includes:

[0020] The flight state constraint equations are set as generalized algebraic constraint equations characterizing the flight mission, which are used to limit the geometric and motion boundaries of the aircraft in the trim state.

[0021] Select one control parameter from the system control parameter vector as the extended parameter, and use the control parameters other than the extended parameter as the non-extended control parameters. Then, concatenate the non-extended control parameters with the system state parameter vector to form an eleven-dimensional extended state variable.

[0022] The vector composed of the time derivatives of the non-extended control parameters is set as the output vector of the generalized algebraic constraint equation to construct a three-dimensional equivalent differential constraint equation.

[0023] By combining the eight-dimensional basic flight dynamics differential equations with the three-dimensional equivalent differential constraint equations, an eleven-dimensional extended bifurcation analysis model is constructed.

[0024] Furthermore, when the generalized algebraic constraint equation is used to define the level flight state, the generalized algebraic constraint equation includes three constraint conditions: the sideslip angle is zero, the roll angle is zero, and the difference between the angle of attack and the pitch angle is zero.

[0025] The extended parameter is elevator deflection angle, and the non-extended control parameters include engine throttle, aileron deflection angle, and rudder deflection angle.

[0026] The equivalent differential constraint equations include:

[0027] The time derivative of the engine throttle is set to the sideslip angle, the time derivative of the aileron deflection angle is set to the roll angle, and the time derivative of the rudder deflection angle is set to the difference between the angle of attack and the pitch angle.

[0028] By adopting the above technical solution, the equivalent differential constraint equation can embed the flight state constraints corresponding to the stationary level flight state into the extended bifurcation analysis model during the trim point search and extension calculation process. It should be noted that setting the time derivative of the engine throttle as the sideslip angle, the time derivative of the aileron deflection angle as the roll angle, and the time derivative of the rudder deflection angle as the difference between the angle of attack and the pitch angle is one specific implementation of the equivalent differential constraint equation.

[0029] In other embodiments, the time derivatives of engine throttle, aileron deflection, and rudder deflection can also be correlated with sideslip angle, roll angle, and the difference between angle of attack and pitch angle according to a preset one-to-one correspondence. This preset one-to-one correspondence ensures that the time derivative of each non-extended control parameter corresponds to a different constraint condition; under the equilibrium point solution conditions of the extended bifurcation analysis model, the time derivatives of all non-extended control parameters are zero, thereby making the sideslip angle, roll angle, and the difference between angle of attack and pitch angle zero. The above equivalent differential constraint equations are only used for trim point search and extension calculations and do not represent the actual time response of engine throttle, aileron deflection, or rudder deflection.

[0030] Furthermore, the pseudo-arc length continuation algorithm is used to alternately execute the prediction step and the correction step for continuation calculation. The prediction step includes:

[0031] The initial equilibrium point is used as the starting point for the extension calculation;

[0032] At the initial equilibrium point, a first-order multivariate Taylor series expansion is performed on the extended bifurcation analysis model;

[0033] Based on the condition that the derivative of the extended state variable at the prediction point is zero, a system of linear equations between the state increment and the extension parameter increment is constructed using the state Jacobian matrix and control derivative vector calculated at the initial equilibrium point.

[0034] Solve the system of linear equations to obtain the tangent direction corresponding to the initial equilibrium point, and obtain the state space coordinates of the prediction point along the tangent direction according to the preset extension step size, which are used as the initial values ​​of the correction step.

[0035] Furthermore, the correction step includes:

[0036] At the predicted point, a first-order multivariate Taylor series expansion is performed on the extended bifurcation analysis model to establish the residual correction equation.

[0037] A pseudo-arc length constraint equation is introduced to ensure that the correction increment vector and the tangent direction obtained in the prediction step satisfy an orthogonal constraint.

[0038] By combining the pseudo-arc length constraint equation with the residual correction equation, an augmented Jacobian matrix equation is constructed.

[0039] The augmented Jacobian matrix equation is solved iteratively using Newton's method, and the convergence of the iteration is determined by the norm of the residual vector.

[0040] When the norm of the residual vector is less than the preset convergence tolerance threshold, the current iteration point is determined as the equilibrium point corresponding to the current extension step.

[0041] Furthermore, the physical dynamic system equations satisfying the flight state constraint equations are obtained, including:

[0042] Extract the balancing numerical mapping relationship between the non-extended control parameters and the extended parameters from the balancing solution set;

[0043] Substitute the balance numerical mapping relationship into the system control parameter vector to form a control parameter vector function that depends only on the extension parameter;

[0044] Substituting the control parameter vector function back into the original eight-dimensional state differential equations of the six-degree-of-freedom nonlinear flight dynamics model of the aircraft yields the physical dynamic system equations that satisfy the flight state constraint equations. These physical dynamic system equations do not include the control parameter derivative terms from the equivalent differential constraint equations; they are used for linearization and eigenvalue calculation at each equilibrium point corresponding to the balancing solution set.

[0045] Further, the physical dynamics system equations are linearized and the Jacobian matrix is ​​calculated, including:

[0046] For each equilibrium point in the balanced solution set, the physical dynamic system equations are linearized to the first order.

[0047] Using the central difference method, positive and negative numerical perturbations are applied to each state variable in the system state parameter vector. The partial derivatives of the physical dynamic system equations with respect to each state variable are calculated and assembled into a Jacobian matrix.

[0048] Using a preset minimum flight speed threshold greater than zero, the denominator of the division terms containing flight speed in the physical dynamics system equations is processed to prevent singularity.

[0049] Solve the characteristic equation corresponding to the Jacobian matrix to obtain the characteristic roots of each equilibrium point.

[0050] Furthermore, the stability and bifurcation type of each equilibrium point are determined based on the eigenvalues, including:

[0051] When the real parts of all characteristic roots corresponding to the equilibrium point are less than zero, the flight state corresponding to the equilibrium point is determined to be a stable state.

[0052] When there is at least one eigenvalue with a real part greater than zero among the eigenvalues ​​corresponding to the equilibrium point, the flight state corresponding to the equilibrium point is determined to be an unstable state.

[0053] During the continuous extension along the extension parameters, when the sign of the real part of the characteristic root corresponding to the adjacent extension step changes, it is determined that there is a stable bifurcation point at the corresponding position.

[0054] When the eigenvalue crossing the imaginary axis is a pure real root, the stability bifurcation point is determined to be a saddle-node bifurcation point.

[0055] When the eigenvalues ​​crossing the imaginary axis are a pair of conjugate complex roots, the stability bifurcation point is determined to be a Hopf bifurcation point. The stability determination result is plotted on the balancing curve to obtain the bifurcation curve.

[0056] This invention provides a bifurcation analysis method for aircraft stability considering flight state constraints. It has the following advantages:

[0057] 1. This invention transforms flight state constraint equations into equivalent differential constraint equations and combines these equations with a six-degree-of-freedom nonlinear flight dynamics model of the aircraft to construct an extended bifurcation analysis model. During the trim point search and extension calculation process, the equivalent differential constraint equations ensure that the obtained equilibrium points satisfy the corresponding flight state constraints, thereby enabling the analysis of trim branches and stability changes under a given flight state.

[0058] 2. This invention extracts the numerical mapping relationship between non-extended control parameters and extended parameters from the trim solution set, and substitutes the numerical mapping relationship back into the six-degree-of-freedom nonlinear flight dynamics model of the aircraft to obtain the physical dynamic system equations that satisfy the flight state constraint equations. At the same time, by alternately executing the prediction step and the correction step through the pseudo-arc length extension algorithm, the trim solution set that varies with the extension parameters can be obtained, and the computational basis can be provided for the linearization processing, characteristic root calculation and bifurcation curve drawing at each equilibrium point.

[0059] 3. This invention calculates a comprehensive hazard index based on the real part of the current largest eigenvalue root, the physical limit boundary of the state variable, and the rated safety buffer margin. The comprehensive hazard index data sequence that changes with the control parameters is compared with preset multi-level safety thresholds. When the comprehensive hazard index crosses the set highest level safety threshold, the critical control parameter vector that triggers the extreme dangerous flight mode is fed back to the upstream flight control computer to trigger the control law reconstruction strategy, thereby providing parameter basis for dangerous flight mode determination and flight control. Attached Figure Description

[0060] Figure 1 This is a flowchart of the process of the present invention;

[0061] Figure 2 The diagram shows the longitudinal aerodynamic characteristics of the UAV of this invention, where (a) is the lift coefficient curve; (b) is the drag coefficient curve; (c) is the polar curve; and (d) is the pitch moment coefficient curve.

[0062] Figure 3 The diagram shows the propeller dynamics model of the present invention, wherein (a) is the propeller thrust model diagram; and (b) is the propeller torque model diagram.

[0063] Figure 4 This is a trim curve diagram of the present invention under constant vertical flight state;

[0064] Figure 5 This is a bifurcation curve diagram showing the relationship between elevator deflection angle and angle of attack in the fixed-level flight state of the present invention.

[0065] Figure 6 This is a bifurcation curve diagram showing the relationship between angle of attack and flight speed in the constant vertical flight state of this invention. Detailed Implementation

[0066] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0067] See attached document Figure 1 As the core system architecture foundation of this invention, this invention provides an aircraft stability bifurcation analysis method considering flight state constraints, including the following main technical execution lines:

[0068] A six-degree-of-freedom nonlinear flight dynamics model for the aircraft is established. Starting from the general physical laws of flight dynamics, existing flight stability analysis methods mainly rely on small perturbation assumptions for local linearization near the trim point, or conduct numerical simulations for single operating conditions, failing to reflect the global evolution of aircraft stability as relevant flight parameters change. To overcome these limitations, this embodiment establishes a six-degree-of-freedom nonlinear flight dynamics equation that includes the aircraft's state and control parameters, serving as the reference physical space for subsequent global dynamic evolution mapping.

[0069] After establishing a six-degree-of-freedom nonlinear flight dynamics model of the aircraft, if the trim solution is directly performed, the eight state differential equations in the model form eight independent equilibrium equations under steady-state conditions. Simultaneously, the system variables to be solved include eight state variables and three control variables (excluding extension parameters), totaling eleven variables. Therefore, relying solely on these eight equilibrium equations is insufficient to uniquely determine the trim solution that satisfies the requirements of a specific flight mission. To match the number of equations with the number of variables to be solved and to limit the specific flight state corresponding to the trim solution, this embodiment introduces flight state constraint equations. These flight state constraint equations are used to express flight mission requirements such as level flight and coordinated turns as algebraic constraints.

[0070] After introducing the flight state constraint equations, a certain control parameter of the aircraft (as a preferred approach, elevator deflection or aileron deflection, which has a direct physical causal relationship with the pitch or roll path, is often chosen) is selected as the extended parameter, and the remaining non-extended control parameters are incorporated into the aircraft's state parameters to form new state variables. Simultaneously, the flight state constraint equations are incorporated into the aircraft's six-degree-of-freedom nonlinear flight dynamics equations to establish an extended bifurcation analysis model. The newly incorporated equations do not possess physical dynamic time evolution characteristics; they only satisfy the system's static equilibrium conditions. The technical purpose is to complete the system's analytical degrees of freedom, enabling the extended bifurcation analysis model to have the foundation for pure algebraic trim solutions.

[0071] Based on the aforementioned algebraic reconstruction model, this embodiment tracks the trim solution branches during the continuous variation of the extension parameters. Since the trim solution branches may have turning points or limit points, near these locations, the parametric extension method, with the extension parameters as the sole independent variable, is easily affected by local slope changes or changes in the Jacobian matrix conditions, leading to difficulties in the convergence of the prediction and correction process for adjacent equilibrium points. To improve the continuity of trim solution branch tracking, this embodiment uses a pseudo-arc-length extension algorithm based on an extended bifurcation analysis model to solve the trim curve of the aircraft under specific flight conditions. Near the known initial equilibrium point, a prediction step is performed to initially estimate the new equilibrium point position by performing a Taylor expansion of the system equations and ignoring higher-order terms. The pseudo-arc-length extension algorithm effectively avoids the computational divergence failure caused by traditional parametric extension when the system Jacobian matrix is ​​singular (i.e., the determinant of the matrix approaches zero at the extreme turning point) by introducing the arc-length parameter as an additional constraint space.

[0072] Due to the truncation error of higher-order terms in Taylor expansion, the coordinates obtained in the prediction step are only approximate values. Based on this, a correction step is performed, which establishes a system of nonlinear algebraic equations and iteratively solves for the correction increments of the state and control parameters to obtain an accurate equilibrium solution. In specific engineering implementations, a convergence tolerance threshold is set to control the computational accuracy of the numerical iteration. This convergence tolerance threshold is used to determine whether the residual vector in the correction step has reached the convergence condition, and its value ranges from 10. -6 Up to 10 -8 The order of magnitude, with specific values ​​determined based on the aircraft's attitude angular displacement sensitivity and the accuracy requirements of the trim calculation. To ensure numerical stability during subsequent central difference method calculations of the Jacobian matrix, the perturbation step size in the central difference method is matched with this convergence tolerance threshold, and the perturbation step size is not less than the square root of the convergence tolerance threshold. By iteratively executing continuous extension calculations of the prediction and correction steps, a global trim solution set is obtained, showing the state variables changing with the extension parameters.

[0073] After obtaining the global trim solution set from solving the steady-state algebraic equations, the trim numerical mapping relationship between the non-extended control parameters and the extended parameters obtained from the extended calculations is extracted. This relationship is then substituted into the initially established six-degree-of-freedom nonlinear flight dynamics model of the aircraft to replace the control parameters, reconstructing the physical dynamic system equations that strictly satisfy the flight state constraints. Based on the reconstructed equations, bifurcation analysis is performed to obtain the global stability information of the aircraft as parameters change under specific flight states.

[0074] At each precise equilibrium point obtained in the extended space, the nonlinear flight dynamics system of the aircraft under flight state constraints is linearized. The eigenvalues ​​of the Jacobian matrix are extracted by calculating the system's Jacobian matrix with respect to the physical state variables. During this process, the trajectory of the eigenvalues ​​of the Jacobian matrix as a function of the extension parameters is fully recorded, especially the evolution characteristics when the matrix experiences rank reduction or eigenvalues ​​cross the imaginary axis near the critical bifurcation point. This serves as a direct mathematical basis for subsequently judging the system's stability and bifurcation type.

[0075] The stability and bifurcation type of the equilibrium point are determined based on the calculated eigenvalues. The stability of the equilibrium point is determined by the real parts of all eigenvalues ​​of the Jacobian matrix at that point: when all eigenvalues ​​have negative real parts, the equilibrium point is considered stable; when any eigenvalue has a positive real part, the equilibrium point is considered unstable; when there is an eigenvalue with a zero real part and the real parts of the remaining eigenvalues ​​are not positive, the equilibrium point is considered a bifurcation point. For bifurcation points, based on whether the eigenvalue with a zero real part is a real root or a complex root, they are further distinguished as saddle-node bifurcation points and Hopf bifurcation points.

[0076] The local stability information extracted from the eigenvalues ​​is directly mapped and added to the trim curve, and a bifurcation curve containing trim information and stability evolution information is drawn. As a comprehensive physical judgment for flight safety assessment, based on the distribution of solid and dashed lines on the bifurcation curve and the location of singular points such as Hopf bifurcation points, stability mutation nodes are identified, and the dynamic boundary conditions for dangerous flight modes such as tailspin and wing rocking after the aircraft is disturbed are determined.

[0077] See attached document Figure 1 Based on the aforementioned system architecture, in order to accurately describe the motion evolution of the aircraft in space, this embodiment elaborates in detail the specific implementation path for establishing the generalized state space expression in the six-degree-of-freedom nonlinear flight dynamics model of the aircraft.

[0078] S101 establishes the fundamental physical assumptions and coordinate system reference for flight dynamics. To construct a mathematical model that balances computational efficiency and engineering accuracy, in this embodiment, the aircraft is equivalent to a rigid body with constant mass distribution, and the ground coordinate system is used as the inertial reference frame for spatial motion calculation. As a preferred engineering implementation method, the planar geodetic assumption is adopted, neglecting the influence of the Earth's rotation and curvature on the flight trajectory. Under this physical reference, the spatial trajectory and attitude changes of the aircraft are described by the Newton-Euler equations, thereby establishing the kinematic and dynamic analysis boundaries of the system.

[0079] S102 defines the system state parameter vector characterizing the instantaneous motion of the aircraft. Based on the general principles of multi-rigid-body dynamics, let the system state parameter vector be... To fully cover the longitudinal and lateral coupled motion modes of the aircraft in three-dimensional space and avoid the omission of stability information due to order reduction, the state parameter vector... It needs to possess complete physical dimensions, which can be represented in vector form as follows:

[0080] ;

[0081] In its implementation, the aforementioned state parameter vector is composed of multiple sub-dimensions of physical features, specifically covering parameters characterizing the relative relationship between the flight trajectory and the incoming flow, parameters characterizing the aircraft's rotational angular velocity, and parameters characterizing the aircraft's spatial attitude. Specifically, For flight speed, representing the magnitude of the velocity vector of the aircraft's center of mass relative to the air; The angle of attack reflects the angle of attack of the airflow between the longitudinal symmetry plane of the aircraft and the velocity vector; The sideslip angle reflects the airflow deflection angle of the velocity vector from the longitudinal symmetry plane of the aircraft. It is the pitch angular velocity; It is the roll angular velocity; Yaw angular velocity; The pitch angle; This refers to the roll angle. These eight physical quantities together constitute the basic set of variables describing the arbitrary spatial motion state of a rigid body aircraft.

[0082] It should be noted that in the subsequent actual aerodynamic calculations, the flight speed It is often used as a denominator in the numerical calculation of the derivatives of angle of attack and sideslip angle. To prevent calculation divergence and matrix singularities caused by the denominator approaching 0 when the aircraft is in a low-speed, high-angle-of-attack or deep stall state, this embodiment uses flight speed... Set a lower threshold. This threshold is pre-calibrated based on the minimum steady-state flight speed boundary of the specific aircraft, and is usually a constant greater than 0, in order to ensure the logical integrity of the extension algorithm under extreme conditions.

[0083] S103 defines the system control parameter vector that determines the direction of the aircraft's motion evolution. Let the system control parameter vector be... To support the bifurcation analysis of various control channels of the aircraft, control parameter vectors... It must include the main control input variables of the aircraft. Since the changes in the attitude and trajectory of the aircraft are essentially controlled by the deflection of aerodynamic control surfaces and the change in engine thrust, specifically, the above-mentioned system control parameter vector covers the longitudinal control surface deflection angle, the lateral control surface deflection angle, and the thrust parameters of the power system.

[0084] In this embodiment, based on the actuation mechanism configuration of a conventional aircraft, the control parameter vector is specifically expanded as follows: The physical meanings of each symbol in the formula are as follows: The elevator deflection angle is used as the main control surface for longitudinal control to adjust the pitch moment. For throttle, it represents the thrust output control amount of the engine power system; The aileron deflection angle is used to disrupt the lateral lift symmetry in order to adjust the rolling moment. This is the rudder deflection angle, used to adjust the yaw moment. This parameter matrix covers the core input variables for changing the aircraft's flight attitude.

[0085] S104, Constructing a generalized state-space differential expression for the nonlinear flight dynamics equations of the aircraft. After establishing the basic space of state variables and control inputs, and combining the aforementioned state parameter vectors and control parameter vectors, based on the principle of equilibrium of internal and external force systems, the following set of nonlinear ordinary differential equations characterizing the continuous evolution of the system state over time is constructed:

[0086] ;

[0087] In this mathematical equation, State parameter vector The time derivative vector represents the instantaneous rate of change of the state variable; This is a nonlinear function vector characterizing the physical-dynamic mapping relationship. This function vector implicitly includes the cross-linking effects of nonlinear aerodynamic forces of the fuselage, thrust of the propulsion system, gravitational components, and gyroscopic torque generated by the aircraft's rotation. This generalized state-space representation establishes the fundamental mathematical structure for the transformation of a nonlinear system from static control input to dynamic motion response. In practical engineering calculations, the aforementioned nonlinear function vector... The numerical calculation of internal forces and moments depends on the specific aerodynamic modeling logic. For the calculation of various aerodynamic and aerodynamic moment components in the body coordinate system, those skilled in the art can use an aerodynamic derivative database combined with current state parameters to perform interpolation fitting and polynomial calculations. The basic process of constructing a nonlinear aerodynamic model using aerodynamic derivatives is a well-known technique in this field and will not be elaborated upon here.

[0088] S201, based on the momentum theorem and relative derivative formulas in a moving coordinate system from theoretical mechanics, maps the spatial forces acting on the aircraft to the continuous evolution of velocity vectors and aerodynamic angles. As the core description of translational kinetics, in this embodiment, the translational kinetic equations, including the time derivatives of flight speed, angle of attack, and sideslip angle, are constructed as follows:

[0089] ;

[0090] ;

[0091] ;

[0092] In the above mechanical equations, The constant mass of the aircraft; It is the acceleration due to gravity; These represent the components of the net external force acting on the aircraft along the longitudinal, transverse, and vertical axes in the body coordinate system. This net external force essentially constitutes a vector space encompassing aerodynamic forces and engine thrust. The engineering purpose of this set of equations is to accurately describe how external excitations alter the aircraft's translational kinetic energy and the relative geometric configuration of the incoming flow and the aircraft body. To ensure the completeness of the algorithmic logic under extreme high-maneuver conditions, the derivative of the angle of attack needs to be adjusted. Singularity avoidance is performed on the division terms in the equation. When the sideslip angle... Approaching When, the denominator term Approaching zero will trigger numerical overflow in the computer. As a preferred engineering constraint, in this embodiment, an envelope constraint threshold is applied to the above differential equation, strictly limiting the absolute value of the sideslip angle within the computational domain to less than 85°. This range is determined based on the maximum sideslip angle flight envelope constraint boundary of conventional fixed-wing aircraft. Large sideslip flows exceeding this boundary are usually accompanied by strong asymmetric vortex breaking, leading to unpredictable nonlinear step jumps in aerodynamic parameters. This constraint condition provides dual protection, both mathematically and physically, for the convergence of the system evolution trajectory and its engineering effectiveness.

[0093] S202, establish the rotational dynamics equations characterizing the aircraft's rotational motion around its center of mass. Based on the angular momentum theorem, this step reveals the dynamic response characteristics of angular acceleration under the coupled effects of external aerodynamic torques and the aircraft's own gyroscopic effects. The rotational dynamics equations, including the time derivatives of roll, pitch, and yaw angular velocities, are constructed as follows:

[0094] ;

[0095] ;

[0096] ;

[0097] In the above equation, These are the resultant external moments of roll, pitch, and yaw in the body coordinate system, respectively. These represent the moments of inertia of the aircraft about its three coordinate axes; This represents the inertial product resulting from the asymmetric mass distribution of the organism. The technical purpose of this set of equations is to express the degree of nonlinear coupling between multiple axes, especially those involving... or The product term directly reflects the inertial cross-linking reaction between roll and dive. Furthermore, the denominator term of the equation... Based on the property that the positive definite matrix of the moment of inertia of a rigid body is always greater than zero, the continuity and non-singularity of the numerical solution of the rotational channel of the nonlinear dynamic system within the entire envelope are guaranteed.

[0098] S203 establishes the kinematic equations characterizing the evolution of the aircraft's spatial attitude. By establishing an instantaneous geometric mapping between the aircraft's following coordinate system and the ground inertial coordinate system, the aircraft's angular velocity is transformed into the Euler angle change rate characterizing spatial orientation. The kinematic differential equations containing the time derivatives of pitch and roll angles are constructed as follows:

[0099] ;

[0100] ;

[0101] By introducing Euler angles, the above kinematic equations relate the aircraft's rotational dynamics to a horizontal reference datum. Similarly, in the equations... The solution contains Item. When the aircraft enters a vertical climb or dive, resulting in a pitch angle Approaching ±90°, this term exhibits the gimbal lock singularity inherent in the Euler angle system. Therefore, before performing the extension calculation, it is necessary to determine whether the extension trajectory crosses this attitude lock boundary. For conventional civilian or non-post-stall maneuvering aircraft, the pitch angle... The analysis boundary is set in the range of [-80°, +80°]. This setting of the analysis boundary not only avoids the risk of data overflow caused by the denominator approaching unboundedness, but also fully covers the core analysis object of this invention, namely, conventional typical flight conditions such as constant straight flight and constant hovering.

[0102] S204, Establish the physical mapping relationship between the control input parameters and the dynamic differential equations. In this embodiment, the mechanism of the control parameters is essentially to drive the above-mentioned six-degree-of-freedom equations by changing the net external force and net external torque. In this embodiment, the aforementioned net external force and net external torque constitute the following mapping relationship:

[0103] ;

[0104] ;

[0105] in, Represents the vector of the combined aerodynamic and thrust functions of the three axes. This represents the composite function vector of the three-axis aerodynamic torques. The above physical mapping illustrates that the aerodynamic and thrust components are the functional results jointly driven by the state variables and control parameter vectors. Regarding the transmission of physical causality, changing the elevator deflection angle... Mainly causes pitching moment and the vertical axis force of the machine body Changes; Change the throttle It primarily determines the engine's thrust increment along the longitudinal axis; it also changes the aileron deflection angle. With rudder deflection This directly interferes with the lateral stress state. The physical causal relationship of this multi-source data determines that the control parameters can be extracted into independent extension parameters in bifurcation analysis, ensuring that the extension operation has real engineering control significance.

[0106] S205, Determine the equilibrium dimension of the basic nonlinear flight dynamics model. According to the general principles of bifurcation theory for nonlinear dynamic systems, when performing aircraft trim and bifurcation analysis, the system's physical state must be in steady-state equilibrium, meaning the instantaneous rate of change of the state variables must be strictly zero. Let the aforementioned dynamic differential equations be... The system is transformed into a system of equilibrium equations consisting of nonlinear algebraic equations, i.e. The engineering purpose of this system of equations is to describe the static equilibrium of aerodynamic forces and moments in all directions of an aircraft. This algebraic system contains a total of eight independent constraint equations. Meanwhile, the unknowns involved in solving the system encompass eight state variables. and 4 control variables The total number of variables reaches 12 dimensions. This dimensionality calculation establishes the variable capacity of the basic model in the full parameter space.

[0107] S206 establishes the solution degree-of-freedom gap based on single-parameter extension. The core technical objective of the numerical extension algorithm in bifurcation analysis is to solve for and trace a continuous one-dimensional balancing manifold curve. A prerequisite for achieving this algorithm's objective is that a specific variable must be selected as an independent extension parameter in the control parameter vector U. As one implementation method, in this embodiment, the elevator deflection angle is selected. As an independent extension parameter, the elevator deflection angle, being one of the longitudinal control inputs of the aircraft, is chosen as an extension parameter to facilitate obtaining the trim solution and its stability as a function of the elevator deflection angle. After specifying this extension parameter, the number of unknowns to be solved in the system is reduced from 12 to 11. Comparing the number of basic equilibrium equations, it can be seen that the number of unknowns is strictly greater than the number of independent equations, thus placing the system in an underdetermined state. Specifically, there is a clear gap of 3 degrees of freedom between the 11 unknowns and the 8 equations. From the perspective of physical space evolution, this gap results in an infinite number of flight attitude solution spaces at a given control parameter point, making it impossible to lock onto a unique physical equilibrium point using conventional algebraic methods. For example, under the same elevator deflection angle and thrust input, due to the lack of lateral constraints, the aircraft may simultaneously correspond to multiple uncertain physical states such as straight level flight, sideslip flight, or spiral maneuvering.

[0108] S207 establishes the technical necessity of introducing flight state constraints. To eliminate the aforementioned dimensional contradictions and ensure the non-singularity and solvability of the equations, additional constraints must be added to the system to achieve mathematical closure. In real flight scenarios, combined with the actual control logic of the aircraft, specific flight missions (such as stationary level flight or coordinated turns) inevitably correspond to strict attitude or trajectory constraints on a macroscopic level. By transforming specific flight missions into algebraic geometric constraints, flight state constraint equations are formed.

[0109] S301 establishes generalized algebraic constraint equations characterizing a specific flight mission. To compensate for the three-dimensional degree-of-freedom gap in the basic model, equivalent constraints must be introduced based on the specific flight mission objectives. In aircraft systems engineering applications, whether it's stationary flight, coordinated turns, or stable hovering, specific flight trajectory or attitude requirements can be abstracted as algebraic functions of maneuvering state parameters. In this embodiment, let the generalized flight state constraint equation vector be... In this equation, The output vector characterizing the target constraint state; This is a vector of nonlinear spatial mapping functions determined by the specific flight mission; Let be the 8-dimensional system state parameter vector defined above. This constraint function vector contains three independent algebraic equations used to physically define the geometric and kinematic boundaries of the aircraft system under specific trim conditions, ensuring that the solution space strictly corresponds to the real engineering task.

[0110] S302, perform a dimension expansion operation on the system variable space based on specific extension parameters. After establishing the above algebraic constraint equations, the system parameter space needs to be repartitioned to meet the dimension alignment requirements of the numerical solver. In the control parameter vector... One control parameter is selected as the extension parameter. In this embodiment, the elevator deflection angle is selected. As extension parameters, the remaining control parameters are incorporated into the state variables as trim variables to be determined. This setup allows analysis of the trim point and its stability as a function of the elevator deflection angle. Subsequently, the control parameter vector... The remaining control parameters, excluding the elevator deflection angle, are extracted and compared with the original state parameter vector. Perform data concatenation. The new vector formed by this concatenation process is defined as the extended-dimensional state variable. .

[0111] Specifically, the remaining control parameters include engine throttle. Aileron deflection and rudder deflection After incorporating it into the original state space, the expanded state variables of the system are expanded into column vectors as follows: Through this dimensional expansion and reconstruction of the variable space, the control parameters, originally known external inputs, are physically transformed into unknown balancing state parameters to be solved within the system, thus enabling the expansion of the state variables. The total dimension is increased to 11. This redefinition of variable dimensions achieves strict equality between the number of unknowns to be solved (11) and the total number of independent equations available in the system (8 basic equations and 3 constraint equations) at the engineering algebra level, thereby eliminating the uncertainty of the solution space and laying the variable foundation for the subsequent construction of a full-rank Jacobian matrix.

[0112] S303, perform the equivalent differential form transformation of the constraint equations. This is to make the above algebraic constraint equations derived from the physical task... To solve the existing flight dynamics differential equations simultaneously within the same extended algorithm framework, an equivalent formal transformation of the algebraic equations is required. Conventional continuous numerical bifurcation analysis algorithms typically only accept standard first-order ordinary differential equation systems as input structures and cannot directly solve algebraic-differential hybrid systems. For this type of algorithm input format, this embodiment sets the extracted residual control parameter vector as... And construct equivalent differential constraint equations The aforementioned equivalent differential constraint equations are used to incorporate flight state constraints into a continuous extension solution framework; they do not represent the actual time response of the remaining control parameters. During the trim point solution process, the system's global derivative is set to zero, and the equivalent differential constraint equations enable the generalized algebraic constraint equations to... The synchronous satisfaction allows flight state constraints to be embedded into the extended bifurcation analysis model, providing an equational basis for subsequent trim branch extension calculations.

[0113] S304, Construct the global differential equation system of the extended bifurcation analysis model. Based on the principle of extended simultaneous equations, the 8-dimensional fundamental flight dynamics differential equations established in step S1 are mathematically expanded using vector dimension expansion with the previously constructed 3-dimensional equivalent differential constraint equations. In this embodiment, a nonlinear differential equation system is established with extended state variables as independent variables and elevator deflection angle as extension parameters. Its generalized expression is as follows:

[0114] ;

[0115] In this formula, This is a column vector of time derivatives of 11-dimensional extended-dimensional state variables; For the extended-dimensional state variables that include the original system state parameters and redundant control parameters, that is:

[0116] ;

[0117] Elevator deflection angle as an independent input item; This is the reconstructed extended nonlinear vector mapping function. Specifically, to achieve closed-loop computation of the system of equations, this nonlinear vector function exhibits a block matrix structure internally, i.e. Among them, the first 8 components Derived from the rigid body flight dynamics logic representing the balance of forces and moments, the latter three components This stems from the requirement that the system must satisfy the geometric or kinematic algebraic constraints of a specific flight mission. Through this system dimensional expansion and fusion, the originally non-closed underdetermined system of equations is reshaped into a dimension-aligned ordinary differential topology.

[0118] S305 explains the applicability of the extended bifurcation analysis model. This relates to the three additional equivalent differential equations added to the above system of equations, namely... Its function is to transform flight state constraints into differential forms that can be solved simultaneously with the flight dynamics equations. This equivalent differential equation is only used for trim point search and extension calculations, and does not represent the actual time response of engine throttle, aileron deflection, or rudder deflection.

[0119] In actual flight, changes in control parameters are determined by the flight control system or control commands. In this embodiment, the aforementioned equivalent differential equation is only used as a constraint embedding method in numerical solutions. When the derivative of the extended state variable is required to be zero during the extended solution process, this equivalent differential equation can ensure that the corresponding flight state constraints are satisfied synchronously. Therefore, the extended bifurcation analysis model is suitable for calculating continuous trim curves and tracking equilibrium solutions, but is not directly used for time integral simulation of the transient response of the aircraft.

[0120] From an engineering implementation perspective, when using numerical extension algorithms to search for the steady-state equilibrium point of a system, it is necessary to force the global derivative to be used. When this happens, the above system of equations must be rigorously derived. The algebraic constraints ensure that the calculated equilibrium solution fully meets the preset requirements for flight missions such as stationary and level flight. Based on the above logical analysis, this extended bifurcation analysis model, due to the inclusion of non-physical pseudo-dynamic variables, can only be strictly limited to the static calculation and manifold tracking process of continuous trim curves, and absolutely cannot be directly used for integral simulation of the system's dynamic time history or evaluation and analysis of transient stability response.

[0121] S306 establishes the numerical interface matching between the extended model and the bifurcation analysis solver. After establishing the above applicable boundaries, the 11th-order nonlinear differential equation system fully satisfies the input specifications of the continuous numerical extension algorithm in terms of algebraic characteristics. At a given elevator deflection angle... At the parameter point, the dimension of the unknown variable to be solved in the system (11-dimensional). The number of equations and the number of independent nonlinear equations (11) are in one-to-one correspondence. This strict alignment between the number of equations and the number of variables to be solved avoids the Newton iteration divergence problem caused by the non-full rank (singularity) of the Jacobian matrix at the algebraic level, thus satisfying the mathematical requirement of continuously pursuing the balanced solution in the extension algorithm.

[0122] S307 establishes the kinematic and geometric constraint equations for stationary level flight. Based on the fundamental kinematic principles of flight dynamics, in real flight scenarios, stationary level flight requires the aircraft to maintain a constant heading and no lateral motion at a given altitude and speed. This macroscopic physical task is transformed into strict algebraic characteristics, requiring that the aircraft's sideslip angle, roll angle, and track angle all remain zero. Specifically, under the assumption of no wind, the longitudinal track angle... Numerically, it strictly follows the kinematic geometric relationship. If the track angle is required To achieve level flight, it is necessary to deduce that the difference between the pitch angle and the angle of attack must be zero. Based on the above physical cause and effect, in this embodiment, three independent algebraic constraint equations are constructed for the fixed-line level flight state, specifically expressed as follows:

[0123] ;

[0124] ;

[0125] ;

[0126] In the above formula, It represents the sideslip angle of an aircraft and is used to constrain the lateral aerodynamic force symmetry of the aircraft and ensure sideslip-free flight. It indicates the roll angle, used to constrain the horizontal attitude of the wing to avoid lateral slippage; For the angle of attack, The pitch angle; Let be the aircraft's trajectory angle. The difference between these two values ​​is designed to ensure that the aircraft's longitudinal flight trajectory does not deviate vertically. These three equations, from both geometric and kinematic perspectives, close the physical boundaries of the stationary, level flight state.

[0127] S308 involves implementing an equivalent differential transformation and variable dimension expansion for specific constraints. To integrate the aforementioned purely algebraic constraints into a nonlinear differential equation analysis solver, an equivalent transformation is necessary. As one implementation method, the elevator deflection angle is selected. As extension parameters, they are used to obtain the trim solution and its stability extension curve as a function of elevator deflection angle. Subsequently, the remaining control parameters, namely engine throttle... Aileron deflection rudder deflection Transformed into undetermined balancing state parameters within the system, an 11-dimensional extended-dimensional state variable is constructed:

[0128] .

[0129] Based on the extended-dimensional state variables, this embodiment establishes three equivalent differential equations to replace the original algebraic constraints:

[0130] ;

[0131] ;

[0132] ;

[0133] In this equivalent construction, , , These represent the time derivatives of the three extended control variables introduced. The above equivalent differential equations are used to incorporate the level flight constraints into the extended bifurcation analysis model. Their effect is limited to trim point search and extension calculations, and does not indicate that the constraints between sideslip angle, roll angle, angle of attack, and pitch angle will actually drive the engine throttle, aileron deflection, or rudder deflection to change over time.

[0134] In the steady-state balancing solution, the derivatives of the extended-dimensional state variables are set to zero. Therefore, the equivalent differential equations described above enable the calculated equilibrium point to satisfy the corresponding fixed-axis and level-plane constraints. Furthermore, employing a first-order linear mapping ensures that the supplementary equations have definite partial derivatives at the corresponding constraint variables, thereby reducing the computational complexity of partial derivatives and minimizing matrix singularities or Newton iteration divergences caused by local derivative anomalies in the constraint equations.

[0135] S309 establishes a global extended analysis model for a fixed-axis level flight scenario. The eight basic rigid-body flight dynamics equations are mathematically combined with the aforementioned three equivalent differential constraint equations. In this embodiment, a model is established based on the elevator deflection angle... Extended analysis model for level flight with independent extension parameters:

[0136] ;

[0137] in, It is the derivative vector of the 11-dimensional state vector; The reconstructed composite mapping vector comprises 11 scalar nonlinear functions, encompassing both the actual flight dynamics derivatives and the aforementioned supplementary virtual trim derivatives. Through the specific model unfolding described above, the original underdetermined system with 8 unknowns and 12 variables is reconstructed into a defined system with 11 simultaneous equations and 11 unsolved variables aligned. This engineering transformation process, while eliminating redundancy in the system's degrees of freedom, implicitly embeds the physical boundaries of level flight, establishing a mathematical foundation with a clear physical task orientation for subsequent use of numerical pseudo-arc length algorithms to trace continuous trim manifold curves.

[0138] S401, Establish the initial equilibrium point of the system and its physical state constraints. As a prerequisite for executing the numerical extension algorithm, the system must obtain an accurate static balancing solution under a specific control parameter. In this embodiment, the known initial equilibrium point is defined as... At this coordinate point, the system must strictly satisfy the steady-state condition that the global derivative is zero, i.e. In this expression, The initial input reference value represents the elevator deflection angle as an extension parameter; This represents the elevator deflection angle at a given value. The system obtains the exact solution vector of the 11-dimensional extended-dimensional state variables through balancing calculations.

[0139] S402, based on the local Taylor expansion approximation of the Jacobian matrix. After obtaining the initial accurate equilibrium point, the local topological space of the system in the neighborhood of that point needs to be linearized and its dimensionality reduced. Since the nonlinear mapping function is difficult to directly provide an analytical solution over the entire envelope, performing a low-order truncation approximation near the known solution is a common engineering method for achieving continuous parameter tracking. In this embodiment, the reconstructed extended nonlinear mapping function... At the initial equilibrium point Perform a first-order multivariate Taylor series expansion at the initial point. After ignoring second-order and higher-order nonlinear terms, the local dynamic characteristics of the system approximately evolve into the following form:

[0140] ;

[0141] In this approximate differential formula, The column vector representing the infinitesimal state increments of the expanded-dimensional state variables of the system relative to the initial solution, with dimensions 11×1, is defined as follows: ; The scalar representing the small control increment of the elevator deflection angle relative to the initial parameters is defined as follows: .

[0142] Furthermore, the partial derivative matrix terms in the formula constitute the directional guide for the algorithm's evolution. Among them, This is the 11×11 order state Jacobian matrix evaluated at the initial equilibrium point of the system, whose internal elements characterize the local coupling gradients between the state parameters of the system. The 11×1 dimensional control derivative vector represents the gradient of the direct driving force of the extension parameter variation on the global mapping equation of the system.

[0143] For the aforementioned algebraic structure, the local singularity of the state Jacobian matrix in continuous computation can affect the stable operation of the algorithm. Although the aforementioned variable dimension expansion operation can make the 11×11 matrix consistent with the dimension of the unknowns, the determinant of the state Jacobian matrix approaches zero when the extension curve passes through critical extreme points such as saddle-knot bifurcation. At this time, if we directly let Inverting the Jacobian matrix to calculate the state increment can easily lead to numerical overflow or iterative divergence at critical bifurcation points.

[0144] See attached document Figure 1 Based on the aforementioned established local linear approximation model, this embodiment further discloses the algebraic derivation and calculation process of the prediction step based on the extension algorithm. In the numerical solution of nonlinear equations, if the initial value of the iteration deviates significantly from the true solution, it can easily lead to iteration divergence. Based on the above-mentioned objective observation, the core technical objective of the prediction step is to calculate the approximate state space coordinates of the next equilibrium point along the tangent direction of the known equilibrium solution, using a preset or adaptively adjusted extension step size based on the convergence situation, as a good initial value for subsequent correction steps.

[0145] S403, Construct the stable zero-derivative condition for the prediction step. In the context of bifurcation analysis of nonlinear systems, searching for a new equilibrium point means that after a small change in parameters, the global state derivative of the system needs to be reset to zero. In this embodiment, based on the aforementioned initial precise equilibrium point... The first-order Taylor approximation differential equation constructed nearby forces the system's dynamic rate of change at the prediction point to be... Therefore, it can be deduced that the initially estimated new equilibrium point must satisfy the following local algebraic constraints:

[0146] ;

[0147] In this equation, the partial derivative matrices of each term are all at the initial equilibrium point. A static evaluation is performed at this point. This formula represents the linear balance relationship that should be satisfied between the state increment and the extension parameter increment under the first-order approximation condition. By solving this relationship, the local tangent direction from the current equilibrium point to the adjacent prediction point can be obtained.

[0148] S404, perform the inverse algebraic matrix operation to solve the predicted state vector. Let the initially estimated new equilibrium point be... .in, The predicted elevator deflection angle for the next extension step; Let be the 11-dimensional predicted state vector of the system under a given deflection angle. Based on the zero derivative condition mentioned above, the state increment of the predicted point can be obtained by rearranging the matrix and solving the system of linear equations.

[0149] In practical numerical computation, although the above derivation can be expressed in matrix inversion form, the inverse of the Jacobian matrix is ​​usually not calculated directly. Instead, it is transformed into a system of linear equations to improve computational stability. When the state Jacobian matrix approaches singularity, direct solution can lead to numerical overflow or iterative divergence. In this case, an augmented matrix equation can be constructed by introducing pseudo-arc length constraints to reduce computational divergence caused by matrix rank reduction near critical bifurcation points.

[0150] S405 explains the approximate properties of the predicted solution. Obtain the predicted points. Subsequently, because the prediction step is based on a first-order Taylor expansion and ignores second-order and higher-order terms, the predicted point is usually not an exact equilibrium solution. Substituting it back into the original nonlinear mapping function will result in a non-zero residual, i.e. The predicted point is used to provide initial values ​​for subsequent Newton correction iterations, and the final balancing solution needs to be obtained through further solving in the correction step.

[0151] S406, Construct the Taylor expansion and residual equation based on the exact equilibrium solution. To establish the mathematical evolution direction for error elimination, the system assumes that an exact equilibrium solution exists within the currently given extended neighborhood. In this embodiment, let the coordinates of this exact equilibrium solution be... It must satisfy the steady-state constraints of flight dynamics, that is... The equilibrium point obtained by the original nonlinear mapping function in the prediction step. We perform a first-order multivariate Taylor series expansion to establish a local approximation equation that approximates the exact solution:

[0152] ;

[0153] In this approximate equation, This is the residual column vector calculated by the system at the prediction point, with dimensions 11x1. This term physically represents the virtual aerodynamic forces or moments that are not fully balanced due to model nonlinearity when the aircraft is in the prediction state. ; The tiny scalar representing the control parameter is defined as follows: The partial derivative matrices are all at the current predicted equilibrium point. Static assignment is performed at the point. The above formula transforms the nonlinear root-finding problem into a linear incremental solution problem. Its technical intention is to use the local state gradient information of the system to solve for the parameter increment that can produce an equal amount of reverse dynamic response, thereby neutralizing and zeroing the current prediction residual.

[0154] S407, construct the augmented Jacobian matrix equation to avoid dimensionality loss. By rearranging the terms in the above expansion equation, we obtain the residual correction equation for iterative calculation:

[0155] ;

[0156] For the residual correction equations, the system must perform algorithm-level data dimension calculations. Because the extension algorithm introduces additional control parameter increments, the above equation set contains only 11 independent equations (corresponding to 11-dimensional residual vectors) in algebraic structure, but simultaneously contains 12 unknowns. This typical underdetermined equation set has infinitely many solutions; directly calling a conventional inverse solver will trigger a singularity error in the computation program. To address the solution dilemma of the above underdetermined system and provide it with a unique analytical basis, as a preferred implementation, this embodiment introduces a scalar orthogonal constraint equation based on pseudo-arc length extension. Specifically, the system mandates that the correction increment vector of the correction step must remain orthogonal to the tangent advance vector of the prediction step. Let the normalized tangent direction vector obtained from the previous step be... ,in It is an 11×1 dimensional column vector of state tangents. To control the tangent components with scalars, we supplement with the following one-dimensional algebraic equation:

[0157] ;

[0158] In the formula, superscript This represents the matrix transpose operation. The geometric constraint of this formula ensures that the iterative correction path is always perpendicular to the current tangent prediction direction, thus guaranteeing that the iteration can quickly converge to the balanced manifold closest to the prediction point. By combining this orthogonal constraint equation with the original 11-dimensional residual correction equation, the system is reconstructed into a 12×12-order established algebraic equation system. This augmented matrix construction not only physically fills the dimensional gap in the underdetermined equations but also, due to the introduction of the tangent geometric vector, effectively avoids the risk of division overflow at critical bifurcation points caused by the rank reduction of the original 11×11-order state Jacobian matrix, where the determinant approaches zero.

[0159] S408 performs multi-dimensional iterative updates and comprehensive tolerance determination. After successfully constructing a full-rank augmented matrix system, the underlying solver calculates the current increment using conventional numerical algebra operations such as Gaussian elimination. and Subsequently, the equilibrium solution coordinates are iteratively updated according to the following rules:

[0160] ;

[0161] ;

[0162] After completing a coordinate update, the system determines whether the current iteration has converged based on the updated residual vector. Specifically, the residual vector is calculated... norm When the norm is less than the preset convergence tolerance threshold At that time, the current iteration is determined to have converged. Based on the aircraft trim accuracy and the double-precision floating-point truncation characteristics of the computer, this convergence tolerance threshold is determined. In this embodiment, it is set to 10. -8 Magnitude. After the convergence condition is met, the updated... That is, it is confirmed as the exact equilibrium solution at that parameter step. If the residual norm is still higher than the threshold, the system recalculates the partial derivative matrix based on the updated coordinates and returns to perform a new round of incremental solution. By establishing this continuous prediction and correction cycle alternation mechanism, the system can reliably track and characterize the global static balancing manifold curve of the state variables as the elevator deflection angle evolves in the multidimensional parameter space.

[0163] S409 defines the multidimensional physical boundaries and iterative derivation mechanism for parameter extension. The underlying solver successfully obtains the exact equilibrium solution at the current parameter step. Then, the system assigns its entire data coordinates to the initial reference point of the next iteration cycle, that is, lets To avoid invalid iterations of the algorithm in non-physical regions that have no engineering significance, the system must perform envelope limit checks on multi-dimensional synthesis logic. As a preferred implementation, the system not only provides extension parameters (i.e., elevator deflection angle)... Pre-loaded mandatory physical action boundary thresholds At the same time, it is the core aerodynamic state variable (such as angle of attack). Flight speed Aerodynamic limit thresholds were set. and The aforementioned action boundary thresholds are typically determined based on the mechanical stop limit deflection angle of the aircraft's actual elevator and service actuators, generally ranging from -25° to +25°. The aerodynamic limit thresholds correspond to the aircraft's stall angle of attack and maximum Mach number limits. After each correction step achieves tolerance convergence, the algorithm is forced to extract the aforementioned multiple-dimensional parameters from the current solution for joint verification. Only when all parameters are within their respective safety boundaries does the system restart the alternating prediction and correction solution based on a new round of tangent directions. If any dimension parameter is detected to reach or exceed its corresponding limit threshold—for example, the angle of attack exceeding the stall redline or the deflection angle reaching the mechanical limit—the system immediately triggers a buffer termination interruption, stopping the algebraic operations along the current evolution branch. Through this closed-loop iterative evolution, the algorithm accumulates a series of densely arranged and continuous discrete-precise equilibrium solution sets in the extended-dimensional state space.

[0164] S410, After obtaining the balanced solution set covering the entire envelope, the system needs to expand the dimension of the state variables. Perform a reverse dimensionality decomposition. Based on the dimensionality reduction and replacement operation performed when establishing the extended bifurcation analysis model in the previous steps, the auxiliary control quantity (throttle) of the original physical system is... Aileron deflection rudder deflection The 11-dimensional extended state variables have been pre-incorporated. The system participates in the root-finding process. At this point, the system independently extracts these three parameter components from the point-state data and establishes their principal extension parameters. Mapping relationship of changing balance values:

[0165] ;

[0166] ;

[0167] ;

[0168] The underlying physical significance lies in the fact that, due to the dynamic symmetry of the aircraft's longitudinal and lateral directions in a fixed-level flight state, the trim angles of the ailerons and rudder are algebraically constrained to a theoretical solution of zero. Furthermore, to counteract aerodynamic drag fluctuations caused by changes in angle of attack and velocity, the throttle command exhibits a specific nonlinear variation law regarding the elevator deflection angle. For the continuous mapping of the aforementioned discrete data points, those skilled in the art can employ conventional reconstruction algorithms such as cubic spline interpolation or least-squares polynomial fitting. The underlying multidimensional grid addressing and fitting smoothing principles are well-known techniques in this field and will not be elaborated upon here.

[0169] S411 integrates the physical attitude vector and outputs a globally continuous balancing curve. In addition to separating auxiliary control parameters, the system simultaneously extracts extended-dimensional state variables. Using conventional 8-dimensional physical flight state parameters, establish a global mapping relationship between the aircraft's state space and control commands:

[0170] ;

[0171] In the formula, The 8×1 dimensional pure flight state column vector recovered after removing the three auxiliary control parameters mentioned above has the following internal elements: ; The global state evolution continuous function matrix is ​​established. Representing this multidimensional functional relationship in the parametric coordinate system generates a continuous trim curve. This curve characterizes the relationship between the control inputs and physical state outputs required for the aircraft to maintain steady-state flight under specific kinematic constraints. The resulting trim numerical mapping relationship is used to characterize the trim values ​​of the non-extended control parameters as the extended parameters change. Substituting this trim numerical mapping relationship back into the original eight-dimensional flight dynamics equations yields a physical dynamics model that satisfies the flight state constraints. This physical dynamics model no longer contains the virtual control derivatives in the equivalent differential constraints and can be used for subsequent linearization at each trim point, as well as for calculating the Jacobian matrix and its eigenvalues.

[0172] S501, Extract the trim function relationship and perform algebraic substitution of the control parameters. The system extracts the algebraic mapping relationship between the auxiliary control parameters and the main extension parameters from the continuous trim manifold generated in the previous steps. In the constant-throttle level flight condition corresponding to this embodiment, the established trim numerical mapping group is as follows:

[0173] ;

[0174] ;

[0175] ;

[0176] In the formula, For throttle command parameters; For aileron deflection angle parameters; This refers to the rudder deflection parameter; Elevator deflection parameter as an independent extension variable; This is a continuous numerical mapping function characterizing the longitudinal trim law. The physical basis of this algebraic substitution is that, in a stationary, level flight state, the aircraft maintains dynamic symmetry in the lateral direction, therefore the trim angle of the ailerons and rudder is always zero; while in longitudinal motion, to maintain the set altitude and speed, the engine thrust must be matched in real time with the aerodynamic drag changes caused by the elevator deflection. Through the above mapping relationship, the system directly substitutes it into the original aircraft control parameter vector constructed in the initial steps. In the middle, after algebraic substitution, the original multidimensional independent control vector is reduced in dimension and converged to a vector that depends on only a single parameter. Vector functions:

[0177] ;

[0178] At the engineering implementation level, this operation essentially materializes the flight state constraints of level flight into the underlying linkage commands of the control system, providing a prerequisite for the subsequent elimination of redundant degrees of freedom in the model.

[0179] S502 reconstructs nonlinear flight dynamics equations with real physical meaning. As a preferred implementation, the system reduces the dimensionality of the control parameter vector... Substituting the entire equations back into the original eight-dimensional state differential equations, we can then determine the net forces acting on the aircraft in the body coordinate system (including aerodynamic components). , , (and engine thrust) and net external torque They are no longer functions of multiple independent control variables, but have been transformed into functions completely controlled by the current state variable. Compared to a single elevator deflection angle The associated variables. To ensure the accuracy of aerodynamic calculations during dynamic evolution, the system needs to utilize the flight velocity in the current state vector when reconstructing the aerodynamic forces and moments model. Real-time Mach number and dynamic pressure are calculated and used as a reference for real-time alignment and three-dimensional interpolation lookup of multi-source aerodynamic data tables. Through the above function substitution and parameter solidification, the original flight dynamics equations are reconstructed into a constrained physical dynamics model, the generalized mathematical representation of which is:

[0180] ;

[0181] In the formula, To restore to the pure physical dimension The flight state parameter vector contains elements that represent flight speed, angle of attack, pitch rate, pitch angle, sideslip angle, roll rate, yaw rate, and roll angle, respectively. This is the original nonlinear implicit function operator; This is the nonlinear explicit dynamic mapping operator formed after substituting static constraints. The reconstructed equations eliminate the artificially added balancing terms in the preceding extended model, and the function structure on its right side can realistically and accurately reflect the system transient recovery or divergence gradient when the aircraft is subjected to external atmospheric disturbances.

[0182] S503 performs algebraic dimension verification and deterministic validation of the engineering system. After establishing the reconstructed model, the system must undergo rigorous alignment checks of mathematical dimensions and degrees of freedom to avoid non-full-rank singularity errors during subsequent Jacobian matrix solving. In this embodiment, the newly reconstructed model consists of a total of 11 fundamental equations from both physical and algebraic perspectives, specifically including 8 purely physical nonlinear differential equations (i.e., The corresponding 8 expansions) and 3 algebraic function relationships (i.e., explicit trim equations for throttle, aileron, and rudder).

[0183] Further analysis of the system's unknown variable pool revealed that the system's overall design includes 12 fundamental unknowns: 8 state variables and 4 control variables. In the engineering algebraic analytical logic, the system forcibly specifies the elevator deflection angle. As the only independent external driving parameter of the system, the remaining 11 unknowns and the 11 established independent equations form a strict one-to-one correspondence in dimensionality. This dimensional balance result proves that the dynamic system after incorporating flight state constraints is in a deterministic state that is exactly solvable on the engineering mathematical level. This state effectively eliminates redundant degrees of freedom interference, fundamentally ensuring that the system state transition matrix will not cause irreversible computational divergence due to dimensional mismatch when performing local linearization and eigenvalue extraction on the reconstructed model.

[0184] S504 is a numerical linearization calculation of the Jacobian matrix for discrete exact equilibrium points. After obtaining the continuous balancing solution set, the system extracts each exact equilibrium solution one by one. As a preferred implementation, the algorithm reconstructs the nonlinear mapping operator at a given equilibrium point. Perform a Taylor expansion and truncate the first-order terms to construct the physical state vector. 8×8 Jacobian matrix :

[0185] ;

[0186] In the formula, the subscript The reference point for calculating partial derivatives is anchored to the current static balancing operating point. Considering that the reconstruction function embeds a multi-dimensional discrete aerodynamic data table, directly deriving analytical partial derivatives is not feasible in the underlying computer architecture. Therefore, this embodiment uses the central difference method to perform the numerical calculation of the partial derivative matrix. In the specific computer execution logic, the system processes the state vector... Each independent element in Apply small numerical perturbations in both positive and negative directions respectively. To balance nonlinear truncation error, residual convergence error, and floating-point rounding error, the perturbation step size is... The dimensionless range is set in

[10] . -4 10 -3 Within the range, and not less than the square root of the aforementioned convergence tolerance threshold. This range is related to 10. -6 Up to 10 -8 The convergence tolerance thresholds are matched to the order of magnitude, which avoids the central difference results being affected by residual convergence errors or floating-point rounding errors due to excessively small perturbation step sizes. Simultaneously, it ensures that the perturbed state variables remain within a local linear neighborhood of the current trim point. Furthermore, it addresses the issues present in the original flight dynamics equations... For the division operation term, the system incorporates denominator anti-singularity processing logic into the computational architecture. This is based on the minimum flight speed threshold applied during boundary verification in the preceding physical theory. The constraint (this threshold corresponds to the stall boundary, and the value is strictly greater than zero) cuts off the path that the denominator approaches 0 and causes numerical divergence from the parameter input end.

[0187] S505, construct the characteristic equation of the state transition matrix and extract local eigenvalues. Complete the Jacobian matrix. After numerical assembly, the system is then subjected to algebraic analysis of its underlying dynamic evolution modes. By solving the characteristic equation of the linearized system, the multidimensional coupled physical states are decoupled into mutually independent natural motion modes. The mathematical expression of the characteristic equation is:

[0188] ;

[0189] In the formula, This refers to the operation of determining the determinant of a matrix. It is an 8×8 identity matrix with the same dimensions as the Jacobian matrix; Let be the characteristic roots of the complex field system to be solved. Based on the fundamental theorem of linear algebra, the above 8th-order equation will analytically yield 8 corresponding characteristic roots. For numerically solving for the eigenvalues ​​of a matrix, those skilled in the art can employ mature matrix iterative algorithms such as QR decomposition. The underlying orthogonalization and eigenvector extraction are well-known techniques in the field and will not be elaborated upon here. After extracting each complex eigenvalue, the system uniformly defines them in standard algebraic form:

[0190] ;

[0191] In the formula, For the first The real part of each characteristic root is the damping coefficient corresponding to the mode, including the attenuation or divergence coefficient, which reflects the transient convergence trend of the system after being disturbed by external factors. The imaginary part of the characteristic root represents the natural oscillation frequency of the dynamic mode; The imaginary unit. Explicitly extracting implicit dynamic information into... The data pairs not only completed the qualitative assessment of the asymptotic stability of a specific equilibrium point, but also established the bifurcation topology quantification criterion for distinguishing between pure exponential divergence (such as saddle knot instability) and periodic oscillations (such as Hopf instability).

[0192] S506 relies on the analytical local asymptotic stability of the real parts of the eigenvalues. The system extracts the eight complex eigenvalues ​​obtained from the preceding solution. (in ), and for the real parts of all eigenvalues Perform a traversal comparison. In this embodiment, the system sets a strict absolute stability criterion. If the real parts of all eigenvalues ​​calculated for the current precise equilibrium point satisfy... If the system determines that the flight state corresponding to the equilibrium point is stable, then from the perspective of physical evolution, after the aircraft is subjected to any small aerodynamic disturbance in this state, the inherent aerodynamic damping of the system can effectively dissipate the disturbance energy, causing the various flight attitude parameters to decay exponentially over time and eventually automatically recover to the initial trim reference.

[0193] Conversely, if at least one real part is detected among the eight eigenvalues ​​satisfying... The system then marks this equilibrium point as an unstable state. This mathematical characteristic physically represents the absence or negative damping of local aerodynamic damping in the system. Under this condition, small disturbances will be continuously amplified by the dynamic coupling mechanism, causing the aircraft attitude to diverge exponentially and deviate from the original balance envelope.

[0194] S507, Numerical capture and topology type quantization classification of critical bifurcation points. (Along the elevator deflection angle) In the continuous parameter extension process, the critical algebraic point where the system state transitions from stable to unstable is called the bifurcation point. In this embodiment, the system locates the critical point by monitoring the real part of the eigenvalue crossing the imaginary axis of the complex plane, i.e., the sign of the real part changes. Considering the inherent truncation and rounding errors in low-level floating-point arithmetic, absolute zero is difficult to solve precisely in the numerical solution domain. To avoid the risk of missed detections, the system introduces a zero-value tolerance threshold in its low-level logic. When the sign of the real part of the eigenvalues ​​changes in adjacent extension steps, and the absolute value is satisfied... When the sign of the real part changes, the system determines that a stability bifurcation has been triggered. This necessitates the introduction of a zero-tolerance threshold. The specific range of values ​​is usually adaptively set based on the condition number of the current system's Jacobian matrix. As a preferred method, it is generally fixed in the range

[10] in conventional double-precision solving environments. -6 10 -4 Within the range.

[0195] After identifying the critical bifurcation point, the imaginary part of the corresponding eigenvalue is further extracted. Determine the type of bifurcation: (1) Saddle-node bifurcation determination: If the characteristic root that satisfies the critical crossing condition is a pure real root, that is, its imaginary part satisfies Then the critical state is determined to be the saddle point bifurcation point. The occurrence of saddle point bifurcation indicates that the aircraft has lost its static stability, and the original equilibrium branch of the system folds away and disappears at this point. In engineering, this usually manifests as an irreversible change in the aircraft's attitude or a transition to another undesirable equilibrium state that is far away, such as a sudden nose pitch or deep stall. (2) Hopf bifurcation determination: If the characteristic roots that satisfy the critical crossing condition are a pair of conjugate complex roots, that is, their imaginary parts are obviously not zero ( If the critical state is reached (e.g., a Hopf bifurcation), then the Hopf bifurcation is determined to be a Hopf bifurcation point. The occurrence of a Hopf bifurcation signifies that the aircraft has lost its dynamic stability, and the system will produce limiting cycle oscillations near the critical point. In engineering, this typically manifests as periodic dynamic oscillations of equal amplitude or divergent amplitude, such as wing rocking or Dutch roll instability.

[0196] S508 constructs a continuous bifurcation map of the global state evolution under control. As a preferred implementation, the system extracts all exact equilibrium solutions along the extension path. The system uses its corresponding stability labels to perform curve fitting and state rendering in a two-dimensional or three-dimensional Cartesian coordinate system. In the conventional longitudinal motion evaluation, the system selects the elevator deflection angle. As the x-axis, select the angle of attack. or pitch angle The ordinate is used as the vertical axis. The mathematical definition of this set of bifurcation curves is:

[0197] ;

[0198] In the formula, This is the set of geometric manifolds representing the global equilibrium state of the system. For continuously changing control input parameters; The static equilibrium state vector that is strictly matched to the control input; This is the reconstructed bounded physics dynamics mapping operator. To intuitively distinguish the dynamic evolution attributes of different branches, the system introduces a linear mapping mechanism in the rendering logic. Specifically, in this operation, the real parts of the corresponding eigenvalues ​​are all less than zero (i.e., the maximum real part). The stable equilibrium sequence of ) is rendered as a solid line; while the real part of ) is greater than zero (i.e. The sequence of unstable equilibrium points is rendered as dashed lines. Saddle-node bifurcation points and Hopf bifurcation points where the system state polarity reverses are anchored to discrete primitives with different geometries or colors. This graphical topological representation directly shows engineers the aircraft's handling safety margin and instability critical boundaries in the current configuration.

[0199] Specific application examples:

[0200] The technical solution of the present invention will be described below using the stability bifurcation analysis of a vector flying wing UAV in a fixed straight-line level flight state as an example.

[0201] Step 1: Establish a nonlinear flight dynamics model. Based on aerodynamic data obtained from wind tunnel tests and calculation results using the vortex lattice method, a nonlinear aerodynamic model of the UAV is established. This nonlinear aerodynamic model is used to determine the aerodynamic coefficients and aerodynamic moment coefficients in the six-degree-of-freedom nonlinear flight dynamics equations of the UAV, and their expressions are as follows:

[0202] ;

[0203] ;

[0204] ;

[0205] ;

[0206] ;

[0207] ;

[0208] in, These represent the lift coefficient, drag coefficient, side force coefficient, roll moment coefficient, pitch moment coefficient, and yaw moment coefficient, respectively. This represents the static aerodynamic coefficient obtained by interpolation from static wind tunnel test data. The dynamic derivative is represented by the vortex lattice method, whose data are obtained and used in the model through interpolation. The longitudinal aerodynamic characteristic curve of the UAV is shown below. Figure 2 As shown.

[0209] Furthermore, a propeller dynamic model was established based on wind tunnel test data of propeller dynamic thrust. This propeller dynamic model is used to determine the thrust and torque generated by the power system, and its expression is as follows:

[0210] ;

[0211] ;

[0212] ;

[0213] In the formula, Indicates the propeller speed. Indicates the throttle position. This indicates the inflow velocity of the propeller. Indicates propeller thrust. This represents the propeller torque. The relationship between propeller thrust and torque as a function of inflow velocity and throttle position is as follows: Figure 3 As shown.

[0214] Based on the characteristics of the aircraft, the forces and torques acting on the UAV under the airframe are first determined through aerodynamic and dynamic models.

[0215] Specifically, the forces and torques acting on the drone can be expressed as:

[0216] ;

[0217] ;

[0218] ;

[0219] ;

[0220] ;

[0221] ;

[0222] in, This indicates the net external force acting on the drone within the drone system. Components on the axis, This indicates the net external force acting on the drone within the drone system. Components on the axis, This indicates the net external force acting on the drone within the drone system. Components on the axis; This indicates the rolling torque acting on the drone. This indicates the pitching moment experienced by the drone. This indicates the yaw moment experienced by the drone; Indicates dynamic pressure. Indicates the wing reference area. Indicates the length of the exhibition. Indicates the mean aerodynamic chord length. Indicates pitch angle, Indicates the roll angle. Indicates the quality of the drone. Represents gravitational acceleration; Indicates the determination of the dynamic model axial thrust component Indicates the determination of the dynamic model axial thrust component; This is the axial force coefficient. This is the lateral force coefficient. The normal force coefficient, , , These are the roll moment coefficient, pitch moment coefficient, and yaw moment coefficient of the UAV, respectively. , , , , , These represent the corresponding torque components in the dynamic model. The force coefficients and torque coefficients mentioned above are determined based on the UAV aerodynamic model, and their expressions are as follows:

[0223] ;

[0224] ;

[0225] in, The lift coefficient, The drag coefficient, This is the angle of attack for flight. After obtaining the above relationships between forces, moments, and aerodynamic coefficients, [the following will be implemented / determined]... , , , , , Substituting the rigid body's six-degree-of-freedom kinematic and dynamic equations, and combining them with the UAV's inertia parameters, the following equations for the UAV's six-degree-of-freedom nonlinear flight dynamics model are obtained for subsequent trim calculations and bifurcation analysis under stationary level flight conditions:

[0226] ;

[0227] ;

[0228] ;

[0229] ;

[0230] ;

[0231] ;

[0232] ;

[0233] ;

[0234] In the formula, , , , , , , , These represent the speed, angle of attack, pitch rate, pitch angle, sideslip angle, roll rate, yaw rate, and roll angle of the UAV, respectively. , , These represent the drone orbiting system. axis, shaft and Moment of inertia of the shaft, Indicates the drone orbiting system shaft and The product of inertia of the axis. The above equations of rotational dynamics retain the product of inertia. Item, not with Zero is taken as a modeling premise; under a specific symmetric configuration, when the inertial product... When negligible, the above equations can degenerate into corresponding simplified forms. The aforementioned six-degree-of-freedom nonlinear flight dynamics model of the UAV is then used for trim calculations and stability bifurcation analysis under constant-level flight conditions.

[0235] Step 2: Construct a constrained bifurcation analysis model for the stationary level flight state, and obtain the trim curve and bifurcation curve for the stationary level flight state based on this model. Building upon the aforementioned six-degree-of-freedom nonlinear flight dynamics model of the UAV, flight state constraint equations corresponding to the stationary level flight state are introduced to ensure that the trim point obtained from the extension calculation satisfies the conditions for the stationary level flight state.

[0236] In this embodiment, the flight state constraint equation under constant vertical flight state is expressed as:

[0237] ;

[0238] ;

[0239] ;

[0240] In the formula, Indicates the sideslip angle. Indicates the roll angle. Indicates the angle of attack. This represents the pitch angle. The above constraints are used to define the relationship between the sideslip angle, roll angle, and angle of attack and pitch angle in stationary level flight.

[0241] Based on the equivalent differential transformation and extended bifurcation analysis model construction method of the aforementioned flight state constraints, the above-mentioned level flight state constraints are incorporated into the extended bifurcation analysis model, and the elevator deflection angle is selected. Continuous extension calculations are performed using these extension parameters to obtain the trim curve under constant-height level flight conditions, such as... Figure 4 As shown.

[0242] Figure 4 The trim curves shown include the trim relationships between speed and elevator deflection, angle of attack and elevator deflection, throttle and elevator deflection, and speed and angle of attack. Based on these trim curves, the flight speed can be obtained under the constraint of constant-speed level flight. Angle of attack and accelerator With elevator deflection The result of the change in balancing.

[0243] After obtaining the trim branches under constant vertical flight conditions, the bifurcation curves under constant vertical flight conditions are plotted based on the stability assessment results at each trim point, such as... Figure 5 and Figure 6 As shown. Among them, Figure 5 The horizontal axis represents the elevator deflection angle, and the vertical axis represents the angle of attack. Figure 6 The horizontal axis represents the angle of attack, and the vertical axis represents the flight speed.

[0244] Depend on Figure 5 and Figure 6 It can be seen that before reaching the stall angle of attack, the trim branch corresponding to the UAV's stationary level flight state is in a stable state; within a certain angle of attack range after stall, the stability of the trim branch changes. For the bifurcation points identified in the bifurcation curve, their stability and bifurcation type can be determined by combining the changes in the eigenvalues ​​of the linearized system at the corresponding trim points. Since some bifurcation points in the bifurcation curve may be affected by numerical calculation errors, the determination results of the bifurcation points should be confirmed in conjunction with the eigenvalue calculation results.

[0245] Through the above calculations, this embodiment obtains the trim curve and bifurcation curve of the UAV in the fixed-level flight state under the constraint of introducing a fixed-level flight state. Based on the stability changes of the trim branches in the bifurcation curve and the location of the bifurcation point, the stability variation range of the UAV in the fixed-level flight state can be analyzed.

[0246] Therefore, this embodiment obtains the constrained bifurcation curve and stability change results under specific flight conditions by establishing a six-degree-of-freedom nonlinear flight dynamics model for the UAV, introducing constraint equations for specific flight states, constructing an extended bifurcation analysis model, and performing continuation calculations. The above analysis process demonstrates that this invention can analyze the equilibrium point and stability boundary of the UAV under given flight state constraints, providing a basis for flight state assessment and control parameter boundary analysis.

Claims

1. A method for bifurcation analysis of aircraft stability considering flight state constraints, characterized in that, Includes the following steps: Establish a six-degree-of-freedom nonlinear flight dynamics model for the aircraft, and define the system state parameter vector and the system control parameter vector; By introducing flight state constraint equations, selecting one control parameter from the system control parameter vector as an extension parameter, and incorporating the control parameters other than the extension parameter as non-extension control parameters into the system state parameter vector to form extended-dimensional state variables, the flight state constraint equations are transformed into equivalent differential constraint equations, and the equivalent differential constraint equations are combined with the six-degree-of-freedom nonlinear flight dynamics model of the aircraft to construct an extended bifurcation analysis model. Based on the extended bifurcation analysis model, at the known initial equilibrium point, the pseudo-arc length extension algorithm is used to alternately execute the prediction step and the correction step to perform extension calculation, obtain the balancing solution set, and extract the balancing numerical mapping relationship between the non-extension control parameters and the extension parameters. Substitute the trim numerical mapping relationship into the system control parameter vector and substitute it back into the six-degree-of-freedom nonlinear flight dynamics model of the aircraft to obtain the physical dynamic system equations that satisfy the flight state constraint equations; Substitute each equilibrium point in the balanced solution set into the physical dynamic system equations for linearization, calculate the Jacobian matrix, and extract the eigenvalues. Based on the characteristic roots, the stability and bifurcation type of each equilibrium point are determined, and the stability determination results are marked on the balancing curve to obtain the bifurcation curve; the dynamic boundary conditions of the dangerous flight mode are determined and used for flight control.

2. The aircraft stability bifurcation analysis method considering flight state constraints according to claim 1, characterized in that, The system state parameter vector includes flight speed, angle of attack, pitch rate, pitch angle, sideslip angle, roll rate, yaw rate, and roll angle; The system control parameter vector includes elevator deflection angle, engine throttle, aileron deflection angle, and rudder deflection angle. The establishment of the six-degree-of-freedom nonlinear flight dynamics model of the aircraft includes: based on the net external forces and net external torques acting on the aircraft, constructing translational mechanics equations characterizing the translational motion of the aircraft's center of mass, rotational mechanics equations characterizing the rotational motion of the aircraft around its center of mass, and kinematic equations characterizing the evolution of the aircraft's spatial attitude. The translational mechanics equations, the rotational mechanics equations, and the kinematic equations together constitute the six-degree-of-freedom nonlinear flight dynamics model of the aircraft.

3. The aircraft stability bifurcation analysis method considering flight state constraints according to claim 2, characterized in that, The construction of the extended bifurcation analysis model includes: The flight state constraint equations are set as generalized algebraic constraint equations characterizing the flight mission, which are used to limit the geometric and motion boundaries of the aircraft in the trim state. Select one control parameter from the system control parameter vector as the extended parameter, and use the control parameters other than the extended parameter as the non-extended control parameters. Then, concatenate the non-extended control parameters with the system state parameter vector to form an eleven-dimensional extended state variable. The vector composed of the time derivatives of the non-extended control parameters is set as the output vector of the generalized algebraic constraint equation to construct a three-dimensional equivalent differential constraint equation. By combining the eight-dimensional basic flight dynamics differential equations with the three-dimensional equivalent differential constraint equations, an eleven-dimensional extended bifurcation analysis model is constructed. The equivalent differential constraint equation is used for balancing point search and extension calculation, but does not represent the actual time response law of the non-extension control parameters.

4. The aircraft stability bifurcation analysis method considering flight state constraints according to claim 3, characterized in that, When the generalized algebraic constraint equation is used to define the level flight state, the generalized algebraic constraint equation includes three constraint conditions: the sideslip angle is zero, the roll angle is zero, and the difference between the angle of attack and the pitch angle is zero. The extended parameter is elevator deflection angle, and the non-extended control parameters include engine throttle, aileron deflection angle, and rudder deflection angle. The equivalent differential constraint equations include: The time derivative of the engine throttle is set to the sideslip angle, the time derivative of the aileron deflection angle is set to the roll angle, and the time derivative of the rudder deflection angle is set to the difference between the angle of attack and the pitch angle.

5. The aircraft stability bifurcation analysis method considering flight state constraints according to claim 1, characterized in that, The pseudo-arc length continuation algorithm is used to alternately execute the prediction step and the correction step for continuation calculation. The prediction step includes: The initial equilibrium point is used as the starting point for the extension calculation; At the initial equilibrium point, a first-order multivariate Taylor series expansion is performed on the extended bifurcation analysis model; Based on the condition that the derivative of the extended state variable at the prediction point is zero, a system of linear equations between the state increment and the extension parameter increment is constructed using the state Jacobian matrix and control derivative vector calculated at the initial equilibrium point. Solve the system of linear equations to obtain the tangent direction corresponding to the initial equilibrium point, and obtain the state space coordinates of the prediction point along the tangent direction according to the preset extension step size, which are used as the initial values ​​of the correction step.

6. The aircraft stability bifurcation analysis method considering flight state constraints according to claim 5, characterized in that, The correction step includes: At the predicted point, a first-order multivariate Taylor series expansion is performed on the extended bifurcation analysis model to establish the residual correction equation. A pseudo-arc length constraint equation is introduced to ensure that the correction increment vector and the tangent direction obtained in the prediction step satisfy an orthogonal constraint. By combining the pseudo-arc length constraint equation with the residual correction equation, an augmented Jacobian matrix equation is constructed. The augmented Jacobian matrix equation is solved iteratively using Newton's method, and the convergence of the iteration is determined by the norm of the residual vector. When the norm of the residual vector is less than the preset convergence tolerance threshold, the current iteration point is determined as the equilibrium point corresponding to the current extension step.

7. The aircraft stability bifurcation analysis method considering flight state constraints according to claim 3, characterized in that, Obtaining the physical dynamic system equations that satisfy the flight state constraint equations includes: Extract the balancing numerical mapping relationship between the non-extended control parameters and the extended parameters from the balancing solution set; Substitute the balance numerical mapping relationship into the system control parameter vector to form a control parameter vector function that depends only on the extension parameters; Substituting the control parameter vector function back into the original eight-dimensional state differential equations in the six-degree-of-freedom nonlinear flight dynamics model of the aircraft, we obtain the physical dynamic system equations that satisfy the flight state constraint equations. The physical dynamics system equations do not include the control parameter derivative terms in the equivalent differential constraint equations, and are used for linearization and characteristic root calculation at each equilibrium point corresponding to the balanced solution set.

8. The aircraft stability bifurcation analysis method considering flight state constraints according to claim 1, characterized in that, The physical dynamics system equations are linearized and the Jacobian matrix is ​​calculated, including: For each equilibrium point in the balanced solution set, the physical dynamic system equations are linearized to the first order. Using the central difference method, positive and negative numerical perturbations are applied to each state variable in the system state parameter vector. The partial derivatives of the physical dynamic system equations with respect to each state variable are calculated and assembled into a Jacobian matrix. Using a preset minimum flight speed threshold greater than zero, the denominator of the division terms containing flight speed in the physical dynamics system equations is subjected to anti-singularity processing. Solve the characteristic equation corresponding to the Jacobian matrix to obtain the characteristic roots of each equilibrium point.

9. The aircraft stability bifurcation analysis method considering flight state constraints according to claim 2, characterized in that, The stability and bifurcation type of each equilibrium point are determined based on the eigenvalues, including: When the real parts of all characteristic roots corresponding to the equilibrium point are less than zero, the flight state corresponding to the equilibrium point is determined to be a stable state. When there is at least one eigenvalue with a real part greater than zero among the eigenvalues ​​corresponding to the equilibrium point, the flight state corresponding to the equilibrium point is determined to be an unstable state. During the continuous extension along the extension parameters, when the sign of the real part of the eigenvalue corresponding to the adjacent extension step changes, it is determined that there is a stable bifurcation point at the corresponding position. When the eigenvalue crossing the imaginary axis is a pure real root, the stability bifurcation point is determined to be a saddle-node bifurcation point. When the eigenvalues ​​crossing the imaginary axis are a pair of conjugate complex roots, the stability bifurcation point is determined to be a Hopf bifurcation point.

Images