Sliding-mode predictive control method for hypersonic flight vehicle

Abstract

The invention discloses a sliding mode prediction control method for a hypersonic aircraft, and relates to the technical field of aerospace. Compared with a previous sliding-mode predictive control method, the sliding-mode predictive control method solves the problem that sliding-mode predictive control design is carried out without considering system constraints; a sliding mode predictive control design method of a continuous time domain is adopted to carry out control design through a rigid hypersonic aircraft model; on the basis of uncertainty estimation, multi-step prediction is carried out on sliding mode dynamics, and an interior point is used for solving a constraint optimization problem. Theoretical derivation and result verification show that the designed control strategy makes full use of a predictive control rolling optimization mechanism, the chattering phenomenon of the control effect is weakened, the influence of uncertainty on the control system can be effectively restrained, and the robustness of the control system is greatly enhanced. The method is oriented to engineering practical problems, the achievement can provide reference for high-performance control of the hypersonic flight vehicle, and the method has important application value.

Description

Technical Field

[0001] The present invention relates to the field of aerospace technology, and in particular to a hypersonic aircraft sliding mode predictive control method. Background Art

[0002] Hypersonic vehicles are a new type of aircraft that fly at Mach 5 or above in near space. Their excellent concealment and defense against surprise attacks may make existing defense systems powerless. They have extremely important military value and have received widespread attention. However, hypersonic aircraft systems have strong uncertainty, and the flight environment is complex and the missions are diverse. In order to meet the guidance requirements, designing high-precision control systems faces huge challenges. Among them, improving anti-interference capabilities is the key to high-precision control.

[0003] Sliding mode control is completely robust to matching disturbances and is one of the most powerful tools for dealing with uncertain systems. It has been widely used in hypersonic vehicle control systems. However, the sliding mode control method also has some disadvantages. First, it is powerless against non-matching disturbances. In order to suppress the impact of non-matching disturbances on the control system, it needs to be further processed in combination with other methods, and the process is often complicated. Second, under the influence of time-varying disturbances, in order to maintain the robustness of the control system to disturbances, the control input often needs to be switched frequently, which is usually unfavorable to the actuator, thus limiting its scope of application.

[0004] Combining sliding mode control with predictive control methods is called sliding mode predictive control method. By predicting the sliding mode dynamics and then designing the control law, the control system's ability to suppress non-matching interference and reduce the jitter phenomenon of the control can be effectively enhanced. And sliding mode predictive control can also effectively deal with system constraint problems, which is crucial to ensuring the actual application performance of the control system. The existing Chinese patent with publication number CN115598973A discloses a nonlinear sliding mode predictive control method, device, equipment and storage medium, and the method includes: S11, constructing an object model group of a nonlinear system through model identification; the object model group includes a PID sub-object model, an actuator sub-object model and a controlled object sub-object model; the object model group is expressed in the form of a state space model; S12, generating a dynamic equation of the object model group, and calculating the disturbance characteristics of the object model group according to the dynamic equation; S13, respectively obtaining the actual output measurement value and state estimation value of the nonlinear system at the current moment; the acquisition of the state estimation value The method includes: taking the model parameters of the object model group as variables at the moment before the current moment, calculating and generating a state estimate value for estimating the state of the nonlinear system at the current moment according to the dynamic equation and its disturbance characteristics; S14, through an extended Kalman filter, recursively calculating the optimal state estimate value of the nonlinear system at the next moment according to the actual output measurement value and the state estimate value; S15, obtaining a corresponding switching function according to the optimal state estimate value; S16, bringing the switching function into the dynamic equation of the object model group and obtaining the optimal solution according to the sliding mode control algorithm; the optimal solution is used as the input of the object model group. This invention combines sliding mode control with predictive control, and uses the switching function as a new variable to substitute into the sliding mode predictive control for optimization and solution; wherein, the switching function and the control input are used as optimization variables in the objective function, and switching function constraints are added, and the control action as the optimal solution can be obtained by minimizing the designed objective function; the present invention uses the switching function in the sliding mode control as a new variable, and penalizes the deviation of the switching function and the control input relative to the equivalent control in the objective function, so that the objective function is zero, which means that the system state reaches the sliding surface, and the control law is equal to the equivalent control.

[0005] In view of the superiority of the sliding mode predictive control method, it has also been studied and applied in the control of hypersonic aircraft. Although the above method and the existing technology further improve the control performance of hypersonic aircraft, they also have the following shortcomings:

[0006] 1. The sliding mode predictive control design is performed without considering the system constraints.

[0007] 2. Use the sliding mode predictive control design method in the continuous time domain. However, most control systems are now implemented based on computers, and computer control systems are typical discrete time systems.

[0008] 3. Control design is performed based on a rigid hypersonic vehicle model; however, due to the slender shape of a hypersonic vehicle, elastic vibrations are inevitable during flight. Therefore, elastic modes cannot be ignored in the control design stage.

[0009] In order to solve the above problems, the present invention proposes a sliding mode predictive control method for a hypersonic aircraft. Summary of the invention

[0010] The purpose of the present invention is to propose a hypersonic aircraft sliding mode predictive control method to solve the problems raised in the background technology:

[0011] The sliding mode predictive control design is performed without considering the system constraints; the sliding mode predictive control design method in the continuous time domain is used to perform control design with a rigid hypersonic vehicle model.

[0012] In order to achieve the above object, the present invention adopts the following technical solutions:

[0013] A hypersonic vehicle sliding mode predictive control method comprises the following steps:

[0014] S1: Build a mathematical model of a hypersonic aircraft, including the longitudinal dynamics model of the hypersonic aircraft, the aerodynamic parameter expression of the hypersonic aircraft and related constraints;

[0015] S2: Design a sliding mode predictive controller combining constraint handling and uncertainty compensation control schemes; specifically including:

[0016] S2.1: Design sliding surface;

[0017] S2.2: Design the sliding mode predictive control law and perform corresponding constraint processing and control solution;

[0018] S3: Apply the optimal control vector obtained by the control solution in S2.2 to the control input of the hypersonic aircraft.

[0019] Preferably, the longitudinal dynamics model of the hypersonic aircraft constructed in S1 is as follows:

[0020] The longitudinal dynamics model includes five rigid states, namely, speed, altitude, track angle, angle of attack, and pitch rate, and six elastic states, as follows:

[0021]

[0022] h=Vsinγ+d2

[0023]

[0024] Among them, V, h, γ, α, Q, are the speed, altitude, track angle, angle of attack, and pitch velocity of the aircraft respectively; are the rate of change of track angle and the rate of change of angle of attack respectively; m, g, I yy are mass, gravitational acceleration and moment of inertia respectively; ζ i (i=1,2,3) is the damping of the i-th elastic mode; ω i (i=1,2,3) is the natural vibration frequency of the i-th elastic mode; is the acceleration of the i-th elastic mode; is the velocity of the i-th elastic mode; η i is the displacement of the ith elastic mode; d1 is the elastic state of velocity; d2 is the elastic state of altitude; d3 is the elastic state of track angle; d4 is the elastic state of attack angle; d5 is the elastic state of angular velocity; L is lift; D is drag; T is thrust; M is pitch moment; N i is the generalized elastic force;

[0025]

[0026] in; is the dynamic pressure; ρ is the air density; S is the reference area; is the average aerodynamic chord length; is the combination vector of elastic modes; is the thrust coefficient related to air density and angle of attack; C T (α) is the thrust coefficient related to the angle of attack; is the thrust coefficient related to the elastic mode; C L (α, δ, η) are lift coefficients related to the angle of attack, control surface deflection angle and elastic mode; C D (α, δ, η) are the drag coefficients related to the angle of attack, the deflection angle of the control surface and the elastic mode; C M (α, δ, η) are the pitching moment coefficients related to the angle of attack, control surface deflection angle, and elastic mode; is the generalized force coefficient related to the square of the angle of attack; is the generalized force coefficient related to the angle of attack; is the generalized force coefficient related to the elevator deflection angle; δ e is the elevator deflection angle; is the generalized force coefficient related to the aileron deflection angle; δ c is the aileron deflection angle; is the generalized force coefficient associated with zero angle of attack and zero control surface deflection angle; are the generalized force coefficients associated with the elastic modes.

[0027] Preferably, the aerodynamic parameter expressions of the hypersonic aircraft constructed in S1 are specifically as follows:

[0028]

[0029] in, The fuel equivalent ratio The thrust coefficient related to the cube of the angle of attack α; The fuel equivalent ratio The thrust coefficient is related to the square of the angle of attack α; The fuel equivalent ratio The thrust coefficient is linearly related to the angle of attack α; The fuel equivalent ratio The associated constant term; is the thrust coefficient related to the cube of the angle of attack α; is the thrust coefficient related to the square of the angle of attack α; The thrust coefficient is linearly dependent on the angle of attack α; Constant term; is the pitch moment coefficient related to the square of the angle of attack; is the pitching moment coefficient related to the angle of attack; is the pitching moment coefficient related to the elevator deflection angle; is the pitching moment coefficient related to the aileron deflection angle; is the pitching moment coefficient associated with zero angle of attack and zero control surface deflection angle; is the pitching moment coefficient associated with the elastic mode; is the lift coefficient related to the square of the angle of attack; is the lift coefficient related to the angle of attack; is the lift coefficient related to the elevator deflection angle; is the lift coefficient related to the aileron deflection angle; is the lift coefficient associated with zero angle of attack and zero control surface deflection angle; is the pitch force coefficient associated with the elastic mode; is the drag coefficient related to the square of the angle of attack; is the drag coefficient related to the angle of attack; is the drag coefficient related to the square of the elevator deflection angle; is the drag coefficient related to the elevator deflection angle; is the drag coefficient related to the square of the aileron deflection angle; is the drag coefficient related to the aileron deflection angle; is the drag coefficient associated with zero angle of attack and zero control surface deflection angle; is the drag coefficient associated with the elastic mode; is the jth coefficient associated with the first elastic mode; is the jth coefficient associated with the second elastic mode; is the jth coefficient associated with the third elastic mode; is the generalized force coefficient associated with the first elastic mode; is the generalized force coefficient associated with the second elastic mode; is the generalized force coefficient associated with the third elastic mode;

[0030] δ e and δ c Satisfy the following formula:

[0031]

[0032] The hypersonic aircraft meets the following constraints:

[0033]

[0034] Among them, Δδ e is the change in the elevator deflection angle; is the change in fuel equivalence ratio.

[0035] Preferably, the sliding surface in S2.1 is designed as follows:

[0036]

[0037] in, is the sliding surface function; is the input vector, T is the transpose of the matrix; u=[Φ,δ e ] T is the control input; d=[d1,d2,d3,d4,d5] T is the disturbance input; f(x) is the evolution function of the system when there is no control input and disturbance input; g(x) is the influence function of the control input on the system state; g1(x) is the influence function of the disturbance input on the system state;

[0038] Discretize the sliding surface:

[0039] x(k+1)=Ax(k)+Bu(k)+Nd(k)

[0040]

[0041] Where, e is a mathematical constant; x(k) is the discretization of the input vector; u(k) is the discretization of the control input; d(k) is the discretization of the disturbance input; k is the discrete time step; A is the system matrix; B is the input matrix; N is the disturbance matrix; T s is the sampling period;

[0042] Define the linear sliding surface s(k) as:

[0043]

[0044] Where C∈R m×n is the switching matrix; c mn is the switching matrix factor, m and n are the vector dimensions;

[0045] In the nominal case, the ideal sliding mode dynamics is expressed as:

[0046]

[0047] Where O is a zero matrix;

[0048] Through linear transformation Z(k) = Λx(k), Z(k) is the state vector after linear transformation; and B2 represents a non-singular matrix, and Λ represents a non-singular transformation matrix;

[0049] Convert the above formula to:

[0050]

[0051] Among them, Z(k+1) is the state vector after linear transformation at time k+1; Z1(k+1) and Z2(k+1) are the two state components of Z(k+1); A 11 and A 12 A is the system submatrix transferred from the state component Z1(k) at time k to Z1(k+1); 21 and A 22 is the system submatrix transferred from the state component Z2(k) at time k to Z2(k+1);

[0052] Based on the above transformation results, the ideal sliding surface is expressed as:

[0053] s(k)=C z1 Z1(k)+C z2 Z2(k)=O m×1

[0054] Among them, C z2 is a reversible matrix, C z1 is the adjustment parameter matrix;

[0055] The converted ideal sliding mode dynamics and ideal sliding mode surface are integrated to obtain:

[0056]

[0057] Take Z2(k) as the input of Z1(k+1), and define

[0058] The expression of Z1(k+1) can be further simplified:

[0059] Select C z2 is the unit matrix, then K=C z1 , the switching matrix is ​​expressed as:

[0060] C=[C z1 C z2 ]Λ=[C z2 KC z2 ]Λ=[KI]Λ

[0061] Where I represents the unit matrix.

[0062] Preferably, the S2.2 is as follows:

[0063] The discretization of the sliding surface is expressed as a one-step prediction form:

[0064] x(k+1|k)=Ax(k|k)+Bu(k|k)+Nd(k|k)

[0065] Among them, x(k+1k) is the state vector at time k+1 predicted at time k; x(kk) is the state vector at time k; u(kk) is the control input vector at time k; d(kk) is the interference vector at time k;

[0066] Through iterative calculation, the full time domain state prediction is obtained as follows:

[0067] X=Νx(k)+ΩU+ΓΕ

[0068] in, It is the full time domain state prediction; is the iterative system matrix; Input matrix for iteration; is the full time domain control input vector; is the iterative interference matrix; is the full time domain interference vector; H p is the prediction time domain; H u To control the time domain;

[0069] The sliding mode dynamic full-time prediction is expressed as:

[0070]

[0071] Among them, s(k+1k) is the sliding surface at time k+1 predicted at time k; s(k+H p k) is the k+H predicted at time k p The sliding surface at the moment; Ψ is a block diagonal matrix;

[0072] Simplify the full time domain prediction of sliding mode dynamics:

[0073]

[0074] in, is the simplified sliding surface prediction vector; Θ = ΨΩ and P = ΨΓ are simplified matrices transformed with respect to the block diagonal matrix Ψ;

[0075] The performance function J(k) is defined as:

[0076]

[0077] Where Q(i) and R(i) are weighted matrices; s(k+i|k) is the sliding surface at time k+i predicted at time k; u(k+i|k) is the control input vector at time k+i predicted at time k; definition The performance function J(k) is further expressed as:

[0078] J(k)=ξ T Qξ+2U T Θ T Qξ+U T [Θ T QΘ+R]U

[0079] in, is a diagonal matrix about Q(i), R = diag([R(1), R(2), ..., R(H u -1)]) is a diagonal matrix about R(i), and diag(·) is a diagonal matrix;

[0080] Define L = 2θ T Qξ and H=[Θ T QΘ+R], the performance function J(k) further expressed in the above formula is further simplified to:

[0081] J(k)=U T HU+U T L+const

[0082] Where const = ξ T Qξ is a term independent of the control vector U;

[0083] The system constraints are processed and converted into inequality constraints of the control vector U; the increment limit in the constraint condition is expressed as:

[0084] E[ΔU(k|k) … ΔU(k+H u -1|k) 1] T ≤vec(0)

[0085] Among them, ΔU(k|k) is the increment of the control vector U at time k; ΔU(k+H u -1|k) is the k+H predicted by the control vector U at time k u -1; vec(0) represents a column vector, each element of which is 0; E is represented as e0 is the last column, expand the above formula:

[0086]

[0087] Rewrite the above formula into vector product form:

[0088]

[0089] The vector product form of the incremental limit in the constraint condition is further expressed as:

[0090]

[0091] in,

[0092] Similarly, the amplitude limit in the constraint condition is expressed as:

[0093] F[U(k|k) … U(k+H u -1|k 1] T ≤vec(0)

[0094] Let F be f is the last column, and the amplitude limit is further expressed as:

[0095]

[0096] Rewrite the above formula into vector product form:

[0097]

[0098] The vector product of the magnitude restrictions in the constraints is further expressed as:

[0099]

[0100] in,

[0101] For the further simplified performance function J(k), remove the term const that is irrelevant to the control vector U, and transform the optimization of the performance function J(k) into minimizing the following performance function J'(k):

[0102] J'(k)=1 / 2(U T 2HU)+UT L

[0103] The performance function J'(k) satisfies the following inequality constraints:

[0104]

[0105] The above performance function J'(k) and the corresponding inequality constraints are solved by the interior point method; the following estimator is used to estimate the unknown interference d in the solution process before optimization:

[0106]

[0107] Where z1(k) represents the estimation of the system state, z2(k) represents the estimation of the disturbance d; Υ, β0, β1 are the estimator design parameters;

[0108] When solving the control problem, let d(k|k)=d(k+H u -1|k) for further solution operation.

[0109] Compared with the prior art, the present invention provides a hypersonic aircraft sliding mode predictive control method, which has the following beneficial effects:

[0110] Based on uncertainty estimation, the present invention uses interior point generation to solve the constrained optimization problem by performing multi-step prediction of sliding mode dynamics. Theoretical derivation and result verification show that the designed control strategy fully utilizes the mechanism of predictive control rolling optimization, which not only weakens the chattering phenomenon of the control action, but also effectively suppresses the influence of uncertainty on the control system, greatly enhancing the robustness of the control system. The present invention is aimed at practical engineering problems, and the results will provide a reference for high-performance control of hypersonic aircraft, which has important application value. BRIEF DESCRIPTION OF THE DRAWINGS

[0111] Figure 1 This is a flow chart of the method mentioned in Example 1 of the present invention;

[0112] Figure 2 The system input and input response diagram under the external interference condition mentioned in Example 1 of the present invention;

[0113] Figure 3 This is a response diagram of other internal variables under the external interference condition mentioned in Example 1 of the present invention;

[0114] Figure 4 This is a sliding surface response diagram under the external interference condition mentioned in Example 1 of the present invention. DETAILED DESCRIPTION

[0115] The technical solutions in the embodiments of the present invention will be described clearly and completely below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments.

[0116] Based on uncertainty estimation, the present invention uses interior point generation to solve the constrained optimization problem by performing multi-step prediction of sliding mode dynamics. Theoretical derivation and result verification show that the designed control strategy makes full use of the mechanism of predictive control rolling optimization, which not only weakens the chattering phenomenon of the control action, but also effectively suppresses the influence of uncertainty on the control system, greatly enhancing the robustness of the control system. The present invention is aimed at practical engineering problems, and the results will provide a reference for high-performance control of hypersonic aircraft, and have important application value. Specifically, it includes the following contents.

[0117] Embodiment 1:

[0118] See also Figure 1-4 The present invention provides a hypersonic aircraft sliding mode predictive control method, comprising:

[0119] S1: Build a mathematical model of a hypersonic aircraft, including the longitudinal dynamics model of the hypersonic aircraft, the aerodynamic parameter expression of the hypersonic aircraft, and related constraints; the details are as follows:

[0120] The longitudinal dynamics model is accurately constructed, covering the motion equations of five rigid states, namely speed, altitude, track angle, angle of attack, and pitch angle rate, and six elastic states. For example, the speed equation needs to comprehensively consider the influence of various factors such as thrust, drag, and gravity on the change of aircraft speed; the altitude equation should accurately describe the motion law of the aircraft in the vertical direction, involving factors such as lift, gravity, and track angle change; the track angle equation should reflect the component effects of thrust, lift, gravity and other forces in the track angle direction; the angle of attack equation should consider the influence of aerodynamic force, pitch moment, and the change of the aircraft's own attitude on the angle of attack; the pitch angle rate equation needs to be related to factors such as pitch moment and moment of inertia; the elastic state equation should be able to accurately describe the vibration characteristics of the elastic mode, including the influence of parameters such as elastic mode damping and natural vibration frequency on elastic deformation. The coefficients in these equations, such as the lift coefficient, drag coefficient, thrust coefficient, pitch moment coefficient, etc., should be accurately determined through theoretical calculations, wind tunnel tests or flight data identification based on factors such as the aircraft's aerodynamic shape and flight conditions, to ensure that the dynamic model can truly reflect the changes in the aircraft's flight state.

[0121] The longitudinal dynamics model is specifically constructed as follows:

[0122]

[0123] h=Vsinγ+d2

[0124]

[0125] Among them, V, h, γ, α, Q, are the speed, altitude, track angle, angle of attack, and pitch velocity of the aircraft respectively; are the rate of change of track angle and the rate of change of angle of attack respectively; m, g, I yy are mass, gravitational acceleration and moment of inertia respectively; ζ i (i=1,2,3) is the damping of the i-th elastic mode; ω i (i=1,2,3) is the natural vibration frequency of the i-th elastic mode; is the acceleration of the i-th elastic mode; is the velocity of the i-th elastic mode; η i is the displacement of the ith elastic mode; d1 is the elastic state of velocity; d2 is the elastic state of altitude; d3 is the elastic state of track angle; d4 is the elastic state of attack angle; d5 is the elastic state of angular velocity; L is lift; D is drag; T is thrust; M is pitch moment; N i is the generalized elastic force.

[0126]

[0127] in; is the dynamic pressure; ρ is the air density; S is the reference area; is the average aerodynamic chord length; is the combination vector of elastic modes; is the thrust coefficient related to air density and angle of attack; C T (α) is the thrust coefficient related to the angle of attack; is the thrust coefficient associated with the elastic mode; C L (α, δ, η) are lift coefficients related to the angle of attack, control surface deflection angle and elastic mode; C D (α, δ, η) are the drag coefficients related to the angle of attack, the deflection angle of the control surface and the elastic mode; C M (α, δ, η) are the pitching moment coefficients related to the angle of attack, control surface deflection angle, and elastic mode; is the generalized force coefficient related to the square of the angle of attack; is the generalized force coefficient related to the angle of attack; is the generalized force coefficient related to the elevator deflection angle; δ e is the elevator deflection angle; is the generalized force coefficient related to the aileron deflection angle; δ c is the aileron deflection angle; is the generalized force coefficient associated with zero angle of attack and zero control surface deflection angle; are the generalized force coefficients associated with the elastic modes.

[0128] Completely establishing the aerodynamic parameter expressions can clarify the expressions of various aerodynamic coefficients related to fuel equivalence ratio, angle of attack, control surface deflection angle (such as elevator deflection angle, aileron deflection angle) and elastic mode. For example, the thrust coefficient should be accurately expressed as a polynomial function relationship with the fuel equivalence ratio and angle of attack to reflect the changes in engine performance with fuel supply and flight attitude; the lift coefficient, drag coefficient, pitch moment coefficient, etc. should reasonably consider the influence of angle of attack, control surface deflection angle and elastic mode, and determine their function form and coefficient value through a large amount of test data and theoretical analysis to ensure that the aerodynamic parameters can accurately describe the aerodynamic force and torque characteristics of the aircraft under different flight conditions. At the same time, clarify the constraint relationship between the aircraft control input (such as fuel equivalence ratio and elevator deflection angle), such as the linear transformation relationship between the two, and the constraints that the aircraft must meet during flight, such as the limit range of the change in elevator deflection angle and the change in fuel equivalence ratio. These constraints are important bases for ensuring the safe and stable flight of the aircraft and must be accurately set during the model building process.

[0129] The aerodynamic parameter expressions for building a hypersonic aircraft are as follows:

[0130]

[0131] in, The fuel equivalent ratio The thrust coefficient related to the cube of the angle of attack α; The fuel equivalent ratio The thrust coefficient is related to the square of the angle of attack α; The fuel equivalent ratio The thrust coefficient is linearly related to the angle of attack α; The fuel equivalent ratio The associated constant term; is the thrust coefficient related to the cube of the angle of attack α; is the thrust coefficient related to the square of the angle of attack α; The thrust coefficient is linearly dependent on the angle of attack α; Constant term; is the pitch moment coefficient related to the square of the angle of attack; is the pitching moment coefficient related to the angle of attack; is the pitching moment coefficient related to the elevator deflection angle; is the pitching moment coefficient related to the aileron deflection angle; is the pitching moment coefficient associated with zero angle of attack and zero control surface deflection angle; is the pitching moment coefficient associated with the elastic mode; is the lift coefficient related to the square of the angle of attack; is the lift coefficient related to the angle of attack; is the lift coefficient related to the elevator deflection angle; is the lift coefficient related to the aileron deflection angle; is the lift coefficient associated with zero angle of attack and zero control surface deflection angle; is the pitch force coefficient associated with the elastic mode; is the drag coefficient related to the square of the angle of attack; is the drag coefficient related to the angle of attack; is the drag coefficient related to the square of the elevator deflection angle; is the drag coefficient related to the elevator deflection angle; is the drag coefficient related to the square of the aileron deflection angle; is the drag coefficient related to the aileron deflection angle; is the drag coefficient associated with zero angle of attack and zero control surface deflection angle; is the drag coefficient associated with the elastic mode; is the jth coefficient associated with the first elastic mode; is the jth coefficient associated with the second elastic mode; is the jth coefficient associated with the third elastic mode; is the generalized force coefficient associated with the first elastic mode; is the generalized force coefficient associated with the second elastic mode; is the generalized force coefficient associated with the third elastic mode.

[0132] δ e and δ c Satisfy the following formula:

[0133]

[0134] Hypersonic vehicles meet the following constraints:

[0135]

[0136] Among them, Δδ e is the change in the elevator deflection angle; is the change in fuel equivalence ratio.

[0137] S2: Design a sliding mode predictive controller; the details are as follows:

[0138] Design the sliding surface as follows:

[0139]

[0140] in, Sliding surface function; is the input vector, T is the transpose of the matrix; u=[Φ,δ e ] T is the control input; d=[d1,d2,d3,d4,d5] T is the disturbance input; f(x) is the evolution function of the system when there is no control input and disturbance input; g(x) is the influence function of the control input on the system state; g1(x) is the influence function of the disturbance input on the system state; through this form, the relationship between the system state change and the control input and disturbance input can be clearly sorted out, providing a basic framework for sliding surface design.

[0141] Discretize the sliding surface:

[0142] x(k+1)=Ax(k)+Bu(k)+Nd(k)

[0143]

[0144] Where, e is a mathematical constant; x(k) is the discretization of the input vector; u(k) is the discretization of the control input; d(k) is the discretization of the disturbance input; k is the discrete time step; A is the system matrix; B is the input matrix; N is the disturbance matrix; T s is the sampling period;

[0145] Define the linear sliding surface s(k) as:

[0146]

[0147] Where C∈R m×n is the switching matrix; c mn is the switching matrix factor, m and n are the vector dimensions; the design of the switching matrix should comprehensively consider the control objectives and performance requirements of the aircraft to ensure that the sliding surface can guide the system state to the desired stable state.

[0148] In the nominal case, the ideal sliding mode dynamics is expressed as:

[0149]

[0150] Where O is a zero matrix;

[0151] Through linear transformation Z(k) = Λx(k), and making B2 represents a non-singular matrix, and Λ represents a non-singular transformation matrix;

[0152] Convert the above formula to:

[0153]

[0154] Among them, Z(k+1) is the state vector after linear transformation at time k+1; Z1(k+1) and Z2(k+1) are the two state components of Z(k+1); A 11 and A 12 A is the system submatrix transferred from the state component Z1(k) at time k to Z1(k+1); 21 and A 22 is the system submatrix transferred from the state component Z2(k) at time k to Z2(k+1);

[0155] Based on the above transformation results, the ideal sliding surface is expressed as:

[0156] s(k)=C z1 Z1(k)+C z2 Z2(k)=O m×1

[0157] Among them, C z2 is a reversible matrix, C z1 is the adjustment parameter matrix;

[0158] The converted ideal sliding mode dynamics and ideal sliding mode surface are integrated to obtain:

[0159]

[0160] Take Z2(k) as the input of Z1(k+1), and define

[0161] The expression of Z1(k+1) can be further simplified:

[0162] Select C z2 is the unit matrix, then K=C z1 , the switching matrix is ​​expressed as:

[0163] C=[C z1 C z2 ]Λ=[C z2 KC z2 ]Λ=[KI]Λ

[0164] Where I represents the unit matrix. The specific form of the sliding surface is finally determined so that it can adapt to the complex flight state changes of the hypersonic vehicle and realize the effective switching and control of the system state.

[0165] Then design the sliding mode predictive control law and perform corresponding constraint processing and control solution:

[0166] The discretization of the sliding surface is expressed as a one-step prediction form:

[0167] x(k+1|k)=Ax(k|k)+Bu(k|k)+Nd(k|k)

[0168] Among them, x(k+1k) is the state vector at time k+1 predicted at time k; x(kk) is the state vector at time k; u(kk) is the control input vector at time k; d(kk) is the interference vector at time k. Through this one-step prediction form, based on the system state, control input and interference estimation at the current moment, the system state at the next moment is predicted, providing basic data for the subsequent control law design.

[0169] Through iterative calculation, the full time domain state prediction is obtained as follows:

[0170] X=Νx(k)+ΩU+ΓΕ

[0171] in, It is the full time domain state prediction; is the iterative system matrix; Input matrix for iteration; is the full time domain control input vector; is the iterative interference matrix; is the full time domain interference vector; H p is the prediction time domain; H u is the control time domain. p and H u The choice of should be determined according to the flight mission and control performance requirements of the hypersonic vehicle to ensure effective prediction and optimization of the system state and control input within a sufficiently long time frame.

[0172] The sliding mode dynamic full-time prediction is expressed as:

[0173]

[0174] Among them, s(k+1k) is the sliding surface at time k+1 predicted at time k; s(k+H p k) is the k+H predicted at time k p The sliding surface at time ; Ψ is a block diagonal matrix.

[0175] Simplify the full time domain prediction of sliding mode dynamics:

[0176]

[0177] in, is the simplified sliding surface prediction vector; Θ=ΨΩ and Ρ=ΨΓ are simplified matrices transformed with respect to the block diagonal matrix Ψ. Through these steps, the sliding mode dynamic prediction is transformed into a form that is convenient for analysis and control, laying the foundation for performance function definition and constraint processing.

[0178] The performance function J(k) is defined as:

[0179]

[0180] Where Q(i) and R(i) are weighted matrices; s(k+i|k) is the sliding surface at time k+i predicted at time k; u(k+i|k) is the control input vector at time k+i predicted at time k; definition The performance function J(k) is further expressed as:

[0181] J(k)=ξ T Qξ+2U T Θ T Qξ+U T [Θ T QΘ+R]U

[0182] in, is a diagonal matrix about Q(i), R = diag([R(1), R(2), ..., R(H u -1)]) is a diagonal matrix with respect to R(i), and diag(·) is a diagonal matrix.

[0183] Define L = 2θ T Qξ and H=[Θ T QΘ+R], the performance function J(k) further expressed in the above formula is further simplified to:

[0184] J(k)=U T HU+U T L+const

[0185] Where const = ξ T Qξ is a term independent of the control vector U. By rationally designing the performance function and comprehensively considering the sliding surface state and the changes in the control input, a clear objective function is provided for the optimization of the control law.

[0186] The system constraints are processed and converted into inequality constraints of the control vector U; the incremental limit in the constraint condition is expressed as:

[0187] E[ΔU(k|k) … ΔU(k+H u -1|k) 1] T ≤vec(0)

[0188] Among them, ΔU(k|k) is the increment of the control vector U at time k; ΔU(k+H u -1|k) is the k+H predicted by the control vector U at time k u -1; vec(0) represents a column vector, each element of which is 0; E is represented as e0 is the last column, expand the above formula:

[0189]

[0190] Rewrite the above formula into vector product form:

[0191]

[0192] The vector product form of the incremental limit in the constraint condition is further expressed as:

[0193]

[0194] in,

[0195] Similarly, the amplitude limit in the constraint condition is expressed as:

[0196] F[U(k|k) LU(k+H u -1|k) 1] T ≤vec(0)

[0197] Let F be f is the last column, and the amplitude limit is further expressed as:

[0198]

[0199] Rewrite the above formula into vector product form:

[0200]

[0201] The vector product of the magnitude restrictions in the constraints is further expressed as:

[0202]

[0203] in,

[0204] Similarly, for the amplitude limit constraint condition, similar processing is performed to convert it into the corresponding inequality constraint form. Through these constraint processing steps, it is ensured that the solution of the control law is carried out under the condition of satisfying the actual operation constraints of the aircraft, ensuring the feasibility and effectiveness of the control strategy.

[0205] For the further simplified performance function J(k), remove the term const that is irrelevant to the control vector U, and transform the optimization of the performance function J(k) into minimizing the following performance function J'(k):

[0206] J'(k)=1 / 2(U T 2HU)+U T L

[0207] The performance function J'(k) satisfies the following inequality constraints:

[0208]

[0209] The above performance function J'(k) and the corresponding inequality constraints are solved by the interior point method. The application of the interior point method requires the reasonable setting of parameters such as the initial point and iteration step size according to the characteristics of the performance function and the inequality constraints, and the optimal solution is gradually approached through iterative calculations. In the solution process, it is necessary to estimate the unknown interference before optimization. The following estimator is used to estimate the unknown interference d in the solution process before optimization:

[0210]

[0211] Where z1(k) represents the estimation of the system state, and z2(k) represents the estimation of the disturbance d; Υ, β0, β1 are the design parameters of the estimator. According to the system state, control input and disturbance estimation at the current moment, the system state and disturbance at the next moment are predicted and estimated to provide accurate disturbance information for the optimization of the control law. At the same time, when solving the control problem, let d(k|k)=d(k+H u -1|k); this is also reasonable in practice, because the control vector U is obtained according to the current disturbance value d(k|k), but only the natural control component u(k|k) is implemented.

[0212] S3: Apply the optimal control vector obtained by the control solution in S2.2 to the control input of the hypersonic aircraft. The details are as follows:

[0213] Through the above steps, the optimal control vector is finally solved and applied to the control input of the hypersonic aircraft to achieve precise control of the aircraft and ensure that the aircraft operates stably and efficiently in complex flight environments.

[0214] Assume that the aircraft performs maneuvers during cruising, the speed increases by 500ft / s and the altitude rises by 800ft. The design parameters are shown in Table 1:

[0215] Table 1 Control parameters

[0216]

[0217]

[0218] The method of this embodiment is represented as SMPC, and the general sliding mode control method is represented as SMC. Considering the external disturbance suffered by the system, specifically the influence of d1=0.7sin(0.9πt), d2=0.2sin(0.9πt), d3=0.4sin(0.2πt), d4=0.6sin(0.3πt), d5=0.4sin(0.8πt), the results are as follows: Figure 2-4 As shown. Figure 2 As shown in the figure, both the control systems based on SMPC and SMC can make the output tracking trajectory good. However, the SMPC method proposed in this embodiment has a faster response speed and a shorter system adjustment time. In addition, in terms of system input, the high-frequency oscillation amplitude of the control system based on SMPC is smaller, which is helpful for actual execution. Figure 3 It can be seen that the control system based on SMPC requires a smaller angle of attack, which can increase the maneuverability of the aircraft. Finally, from the response of the sliding surface, the control method SMPC proposed in this embodiment has a great advantage in buffeting suppression. For details, please refer to Figure 4 In summary, the control method proposed in this embodiment has strong robustness.

[0219] The above description is only a preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any technician familiar with the technical field can make equivalent replacements or changes according to the technical scheme and inventive concept of the present invention within the technical scope disclosed by the present invention, which should be covered by the protection scope of the present invention.

Claims

1. A hypersonic vehicle sliding mode predictive control method, characterized in that: The steps include: S1: Build a mathematical model of a hypersonic aircraft, including the longitudinal dynamics model of the hypersonic aircraft, the aerodynamic parameter expression of the hypersonic aircraft and related constraints; S2: Design a sliding mode predictive controller combining constraint handling and uncertainty compensation control schemes; specifically including: S2.1: Design sliding surface; S2.2: Design the sliding mode predictive control law and perform corresponding constraint processing and control solution; S3: Apply the optimal control vector obtained by the control solution in S2.2 to the control input of the hypersonic aircraft.

2. A hypersonic vehicle sliding mode predictive control method according to claim 1, characterized in that: The longitudinal dynamics model of the hypersonic aircraft constructed in S1 is as follows: The longitudinal dynamics model includes five rigid states, namely, speed, altitude, track angle, angle of attack, and pitch rate, and six elastic states, as follows: h=Vsinγ+d2 Among them, V, h, γ, α, Q, are the speed, altitude, track angle, angle of attack, and pitch velocity of the aircraft respectively; are the rate of change of track angle and the rate of change of angle of attack respectively; m, g, I yy are mass, gravitational acceleration and moment of inertia respectively; ζ i (i=1,2,3) is the damping of the i-th elastic mode; ω i (i=1,2,3) is the natural vibration frequency of the i-th elastic mode; is the acceleration of the i-th elastic mode; is the velocity of the i-th elastic mode; η i is the displacement of the ith elastic mode; d1 is the elastic state of velocity; d2 is the elastic state of altitude; d3 is the elastic state of track angle; d4 is the elastic state of attack angle; d5 is the elastic state of angular velocity; L is lift; D is drag; T is thrust; M is pitch moment; N i is the generalized elastic force; in; is the dynamic pressure; ρ is the air density; S is the reference area; c is the average aerodynamic chord length; is the combination vector of elastic modes; is the thrust coefficient related to air density and angle of attack; C T (α) is the thrust coefficient related to the angle of attack; is the thrust coefficient associated with the elastic mode; C L (α, δ, η) are lift coefficients related to the angle of attack, control surface deflection angle and elastic mode; C D (α, δ, η) are the drag coefficients related to the angle of attack, the control surface deflection angle and the elastic mode; C M (α, δ, η) are the pitching moment coefficients related to the angle of attack, control surface deflection angle, and elastic mode; is the generalized force coefficient related to the square of the angle of attack; is the generalized force coefficient related to the angle of attack; is the generalized force coefficient related to the elevator deflection angle; δ e is the elevator deflection angle; is the generalized force coefficient related to the aileron deflection angle; δ c is the aileron deflection angle; is the generalized force coefficient associated with zero angle of attack and zero control surface deflection angle; are the generalized force coefficients associated with the elastic modes.

3. A hypersonic vehicle sliding mode predictive control method according to claim 2, characterized in that: The aerodynamic parameter expressions of the hypersonic aircraft constructed in S1 are as follows: in, The fuel equivalent ratio The thrust coefficient related to the cube of the angle of attack α; The fuel equivalent ratio The thrust coefficient is related to the square of the angle of attack α; The fuel equivalent ratio The thrust coefficient is linearly related to the angle of attack α; The fuel equivalent ratio The associated constant term; is the thrust coefficient related to the cube of the angle of attack α; is the thrust coefficient related to the square of the angle of attack α; The thrust coefficient is linearly dependent on the angle of attack α; Constant term; is the pitching moment coefficient related to the square of the angle of attack; is the pitching moment coefficient related to the angle of attack; is the pitching moment coefficient related to the elevator deflection angle; is the pitching moment coefficient related to the aileron deflection angle; is the pitching moment coefficient associated with zero angle of attack and zero control surface deflection angle; is the pitching moment coefficient associated with the elastic mode; is the lift coefficient related to the square of the angle of attack; is the lift coefficient related to the angle of attack; is the lift coefficient related to the elevator deflection angle; is the lift coefficient related to the aileron deflection angle; is the lift coefficient associated with zero angle of attack and zero control surface deflection angle; is the pitch force coefficient associated with the elastic mode; is the drag coefficient related to the square of the angle of attack; is the drag coefficient related to the angle of attack; is the drag coefficient related to the square of the elevator deflection angle; is the drag coefficient related to the elevator deflection angle; is the drag coefficient related to the square of the aileron deflection angle; is the drag coefficient related to the aileron deflection angle; is the drag coefficient associated with zero angle of attack and zero control surface deflection angle; is the drag coefficient associated with the elastic mode; is the jth coefficient associated with the first elastic mode; is the jth coefficient associated with the second elastic mode; is the jth coefficient associated with the third elastic mode; is the generalized force coefficient associated with the first elastic mode; is the generalized force coefficient associated with the second elastic mode; is the generalized force coefficient associated with the third elastic mode; δ e and δ c Satisfy the following formula: The hypersonic aircraft meets the following constraints: Among them, Δδ e is the change in the elevator deflection angle; is the change in fuel equivalence ratio.

4. A hypersonic vehicle sliding mode predictive control method according to claim 3, characterized in that: The sliding surface design in S2.1 is as follows: in, is the sliding surface function; is the input vector, T is the transpose of the matrix; u=[Φ,δ e ] T is the control input; d=[d1,d2,d3,d4,d5] T is the disturbance input; f(x) is the evolution function of the system when there is no control input and disturbance input; g(x) is the influence function of the control input on the system state; g1(x) is the influence function of the disturbance input on the system state; Discretize the sliding surface: x(k+1)=Ax(k)+Bu(k)+Nd(k) Where, e is a mathematical constant; x(k) is the discretization of the input vector; u(k) is the discretization of the control input; d(k) is the discretization of the disturbance input; k is the discrete time step; A is the system matrix; B is the input matrix; N is the disturbance matrix; T s is the sampling period; Define the linear sliding surface s(k) as: Where C∈R m×n is the switching matrix; c mn is the switching matrix factor, m and n are the vector dimensions; In the nominal case, the ideal sliding mode dynamics is expressed as: Where O is a zero matrix; Through linear transformation Z(k) = Λx(k), Z(k) is the state vector after linear transformation; and B2 represents a non-singular matrix, and Λ represents a non-singular transformation matrix; Convert the above formula to: Among them, Z(k+1) is the state vector after linear transformation at time k+1; Z1(k+1) and Z2(k+1) are the two state components of Z(k+1); A 11 and A 12 A is the system submatrix transferred from the state component Z1(k) at time k to Z1(k+1); 21 and A 22 is the system submatrix transferred from the state component Z2(k) at time k to Z2(k+1); Based on the above transformation results, the ideal sliding surface is expressed as: s(k)=C z1 Z1(k)+C z2 Z2(k)=O m×1 Among them, C z2 is a reversible matrix, C z1 is the adjustment parameter matrix; The converted ideal sliding mode dynamics and ideal sliding mode surface are integrated to obtain: Take Z2(k) as the input of Z1(k+1), and define The expression of Z1(k+1) can be further simplified: Select C z2 is the unit matrix, then K=C z1 , the switching matrix is ​​expressed as: C=[C z1 C z2 ]Λ=[C z2 KC z2 ]Λ=[KI]Λ Here, I represents the unit matrix.

5. A hypersonic vehicle sliding mode predictive control method according to claim 4, characterized in that: The S2.2 is as follows: The discretization of the sliding surface is expressed as a one-step prediction form: x(k+1|k)=Ax(k|k)+Bu(k|k)+Nd(k|k) Among them, x(k+1|k) is the state vector at time k+1 predicted at time k; x(k|k) is the state vector at time k; u(k|k) is the control input vector at time k; d(k|k) is the interference vector at time k; Through iterative calculation, the full time domain state prediction is obtained as follows: X=Νx(k)+ΩU+ΓΕ in, It is the full time domain state prediction; is the iterative system matrix; Input matrix for iteration; is the full time domain control input vector; is the iterative interference matrix; is the full time domain interference vector; H p is the prediction time domain; H u To control the time domain; The sliding mode dynamic full-time prediction is expressed as: Among them, s(k+1|k) is the sliding surface at time k+1 predicted at time k; s(k+H p |k) is the k+H predicted at time k p The sliding surface at the moment; Ψ is a block diagonal matrix; Simplify the full time domain prediction of sliding mode dynamics: in, is the simplified sliding surface prediction vector; Θ = ΨΩ and P = ΨΓ are simplified matrices transformed with respect to the block diagonal matrix Ψ; The performance function J(k) is defined as: Where Q(i) and R(i) are weighted matrices; s(k+i|k) is the sliding surface at time k+i predicted at time k; u(k+i|k) is the control input vector at time k+i predicted at time k; definition The performance function J(k) is further expressed as: J(k)=ξ T Qξ+2U T I T Qξ+U T [I T QΘ+R]U Where, Q = diag([Q(1), Q(2), ..., Q(H p )]) is a diagonal matrix about Q(i), R = diag([R(1), R(2), …, R(H u -1)]) is a diagonal matrix about R(i), and diag(·) is a diagonal matrix; Define L = 2θ T Qξ and H=[Θ T QΘ+R], the performance function J(k) further expressed in the above formula is further simplified to: J(k)=U T HU+U T L+const Where const = ξ T Qξ is a term independent of the control vector U; The system constraints are processed and converted into inequality constraints of the control vector U; the increment limit in the constraint condition is expressed as: E[ΔU(k|k)…ΔU(k+H u -1|k) 1] T ≤vec(0) Among them, ΔU(k|k) is the increment of the control vector U at time k; ΔU(k+H u -1|k) is the k+H predicted by the control vector U at time k u -1; vec(0) represents a column vector, each element of which is 0; E is represented as e0 is the last column, expand the above formula: Rewrite the above formula into vector product form: The vector product form of the incremental limit in the constraint condition is further expressed as: in, Similarly, the amplitude limit in the constraint condition is expressed as: F[U(k|k)…U(k+H u -1|k) 1] T ≤vec(0) Let F be f is the last column, and the amplitude limit is further expressed as: Rewrite the above formula into vector product form: The vector product of the magnitude restrictions in the constraints is further expressed as: in, For the further simplified performance function J(k), remove the term const that is irrelevant to the control vector U, and transform the optimization of the performance function J(k) into minimizing the following performance function J'(k): J'(k)=1 / 2(U T 2HU)+U T L The performance function J'(k) satisfies the following inequality constraints: The above performance function J'(k) and the corresponding inequality constraints are solved by the interior point method; the following estimator is used to estimate the unknown interference d in the solution process before optimization: Where z1(k) represents the estimation of the system state, z2(k) represents the estimation of the disturbance d; γ, β0, β1 are the estimator design parameters; When solving the control problem, let d(k|k)=d(k+H u -1|k) for further solution operation.

Images