# Robust Controller Designing for an Air-Breathing Hypersonic Vehicle with an HOSVD-Based LPV Model

Authors: Xing He; Wei Jiang; Caisheng Jiang
Journal: International Journal of Aerospace Engineering
 (2021)

---

## Abstract

This paper focuses on the linear parameter varying (LPV) modeling and controller design for a flexible air-breathing hypersonic vehicle (AHV). Firstly, by selecting the measurable altitude and velocity as gain-scheduled variables, the original longitudinal nonlinear model for AHV is transformed into the LPV model via average gridding division, vertex trimming, Jacobian linearization, and multiple linear regression within the entire flight envelope. Secondly, using the tensor product model transformation method, the obtained LPV model is converted into the polytopic LPV model via high-order singular value decomposition (HOSVD). Third, the validity and applicability of the HOSVD-based LPV model are further demonstrated by designing a robust controller for command tracking control during maneuvering flight over a large envelope.


---

## Body

## 1. Introduction

Since the 50s and 60s of the 20th century, with the series of events such as the advent of long-range ballistic missiles, the successful return of manned spacecraft, and the X-15 test aircraft velocity faster than Mach 6, it marked that mankind formally entered the era of hypersonic velocity. The hypersonic vehicle has taken advantages of both spacecraft and aircraft and has become the new technological commanding point in the 21st century aerospace field. It has been hailed as the third epoch-making milestone in the history of aviation in the world after the invention of aircraft and the breakthrough of the sound barrier [1–3].Compared with a conventional aircraft, the complex aerodynamic characteristics of hypersonic vehicles and the application of technologies such as scramjet engines and integrated airframe/engine configuration make the coupling more powerful between elastic bodies, propulsion systems, aerodynamics, and structural dynamics of hypersonic vehicles; besides, the physical model is more nonlinear. At present, the main means of dynamic modeling for hypersonic vehicles are computational fluid dynamics (CFD), wind tunnel tests, aircraft tests, and other test methods. Because of the imperfection of hypersonic aerodynamic database and CFD software, the low degree of wind tunnel test approximation, the high cost and risk of flight tests, and the limitation of hypersonic vehicles modeling methods, the system modeling has great structural and parameter uncertainties. Second, the large-span flight envelope of hypersonic vehicles not only requires the flight control system to adapt to different flight environments but also requires to satisfy the strict constraints such as more complex dynamic pressure, overload, flight envelopes, and actual physical characteristics (e.g., actuator saturation). In summary, a hypersonic vehicle is a highly dynamic, strongly coupled, fast time-varying, and multivariable nonlinear system. Its control system design not only faces the influence of model uncertainty and external environment but also faces the constraints of actual existence and the time-delay dynamic limitation, which makes the control, stability, safety, rapidity, and accuracy of hypersonic vehicle systems face great challenges. It has become one of the key issues to be solved and is also a hot issue in the field of flight control, such as adaptive control strategies based on the multiple Lyapunov function method or using an auxiliary reaction control system for air-breathing hypersonic vehicles [4–8].Recently, the LPV controller has emerged as a reliable alternative to classical gain scheduling approaches and provides guaranteed stability, robustness, or performance [9, 10]. Gain-scheduled control strategies for LPV systems have been developed intensively and bring promising solutions in many fields, such as missiles, aircrafts, spacecrafts, winding systems, automated vehicles, and robotic systems [11–13, 14]. Zhang handled the reference tracking for a class of flexible hypersonic vehicles using a novel switched LPV framework [15]. Yang proposed a fault-tolerant controller based on robust model-predictive control and A polytopic LPV model for A hypersonic vehicle with external disturbances and actuator loss of effectiveness faults [16]. Hu proposes a novel passive fault-tolerant control method using weighted tube-based model predictive control via polytopic LPV for air-breathing hypersonic vehicles [17]. However, the formulation of the polytopic LPV model for nonlinear systems and controller design for LPV systems with uncertainty and external disturbance remain difficult, which motivates our study.This paper takes the flight control system design of the hypersonic vehicle as the background. In the process of modeling, it is necessary to deal with various factors that affect the performance of the control system. Aircraft modeling errors, flight parameter changes of large package flight, elastic vibration deformation of aircraft, and other factors are regarded as internal factors leading to system uncertainty. Factors such as the complex flight environment, dynamic pressure effect, aerodynamic thermal effect, and the frictional resistance effect are regarded as external disturbances. Then, the attitude stability and command tracking issues during the hypersonic vehicle maneuvering flight in large envelope are transformed into the robust controller synthesis and guaranteed performance control of the uncertain disturbance system. Firstly, for a nonlinear longitudinal elasticity model of air-breathing hypersonic vehicles, the altitude and velocity are selected as scheduling parameters, uniform meshing is performed within the range of parameter variation, and then, the LPV model of the aircraft is obtained by Jacobi linearization and multivariate linear fitting. Secondly, the LPV model of the aircraft is transformed into a polytopic LPV model by the tensor product model change method, which lays a foundation for the subsequent aircraft controller design. Finally, a robust controller is designed based on the obtained aircraft polytopic LPV model, which verifies the effectiveness of the aircraft polytopic LPV model in the command tracking control of aircraft maneuvering flight in large envelopes.


## 2. Problem Formulation

### 2.1. Control-Oriented Model

The longitudinal elastic nonlinear model of the hypersonic vehicle used in this paper is mainly based on the winged-cone model provided by NASA’s Langley Research Center [18–20], which takes into account the influence of the elastic modality of the aircraft. The simplified longitudinal elastic nonlinear model of the aircraft is as follows:
(1)ḣ=Vsinθ−α,V̇=Tcosα−Dm−gsinθ−α,α̇=−Tsinα−LmV+Q+gcosθ−αV,θ̇=Q,Q̇=MIyy,η¨i=−2ςiωiη̇i−ωi2ηi+Ni,i=1,2.In equation (1), h represents the aircraft altitude, V represents the aircraft velocity, α indicates the angle of attack, θ denotes the pitch angle, Q represents the pitch rate, g indicates the gravitational constant, M is the pitching moment, Iyyrepresents the moment of inertia around the body axis, T is the engine thrust, D represents resistance, L represents lift, ηi(i=1,2) represents the elastic mode, ωi(i=1,2) represents the natural frequency of the elastic mode, ςi(i=1,2) represents the damping coefficient of the elastic mode, and Nidenotes the generalized force.The approximate expressions of the above aerodynamics, moments, thrusts, and generalized forces fitted near the equilibrium point are shown in equation (2). In equation (2), ρindicates air density, S is the reference area, CLα,δe represents the lift coefficient related to the angle of attack and elevator deflection, CDα,δe indicates the drag coefficient related to the angle of attack and elevator deflection, CM,αα indicates the torque coefficient related to the angle of attack, CM,δeδe indicates the torque coefficient related to the elevator deflection, zT represents the coupling coefficient between thrust and pitching moment, CTαj represents the thrust coefficient related to αj(i=1,2,3), Niαj represents the generalized force coefficients related to αj, and N2δedenotes the influence coefficient of elevator declination to generalized force N2, and CT0 and Ni0 represent constant items in thrust and generalized forces, respectively.
(2)L≈12ρV2SCLα,δe,D≈12ρV2SCDα,δe,M≈zTT+12ρV2Sc¯CM,αα+CM,δeδe,T≈CTα3α3+CTα2α2+CTαα+CT0,N1≈N1α2α2+N1αα+N10,N2≈N2α2α2+N2αα+N2δeδe+N20,(3)ρ=ρ0exp−h−h0hs,CL=CLαα+CLδeδe+CL0,CD=CDα2α2+CDαα+CDδe2δe2+CLδeδe+CD0,CM,α=CM,αα2α2+CM,ααα+CM,α0,CM,δe=ceδe,q¯=0.5ρV2,CTα3=β1h,q¯Φ+β2h,q¯,CTα2=β3h,q¯Φ+β4h,q¯,CTα=β5h,q¯Φ+β6h,q¯,CT0=β7h,q¯Φ+β8h,q¯.The relevant force and moment coefficients are shown in equation (3), where 1/hsis the air density decay rate, CL0 and CD0 are the constant term in CLα,δeand CDα,δe, respectively, CLαi and CDαi are the ith-order coefficient of α contribution to CLα,δe and CDα,δe, respectively, CLδe is the coefficient of δe contribution to CLα,δe, CDδei is the ith-order coefficient of δe contribution to CDα,δe, CM,ααi is the ith order coefficient of α contribution to CM,αα, CM,α0 is the constant term in CM,αα, q¯ is the dynamic pressure, and Ф is the stoichiometrically normalized fuel-to-air ratio, βih,q¯ is the ith thrust fit parameter. The parameter reference values of the elastic nonlinear model are detailed in Hu et al. [21] and listed as Table 1.Table 1
Nomenclature and nominal values.
SymbolNominal valuem300 slug/ftIyy5×105lb·ftc¯17 ftρ06.7429×10−5slugs·ft−3h085000 fths2.1358×104ftg31.92 ft·s−2S17 ft2·ft−1zT8.36 ftCLα4.6773 rad−1CLδe0.76224 rad−1CL0−1.8714×10−2CDα25.8224 rad−2CDα−4.5315×10−2rad−1CDδe28.1993×10−1rad−2CDδe2.7699×10−4rad−1CD01.0131×10−2CM,αα26.2926 rad−2CM,αα2.1335 rad−1CM,α00.18979ce−1.2897 rad−1β1−3.7693×105lb·ft−1·rad−3β2−3.7225×104lb·ft−1·rad−3β32.6814×104lb·ft−1·rad−3β4−1.7277×104lb·ft−1·rad−3β53.5542×104lb·ft−1·rad−3β6−2.4216×103lb·ft−1·rad−3β76.3785×103lb·ft−1β8−1.0090×102lb·ft−1N1α21.4013×103lb·ft−1·slug−0.5·rad−2N1α4.5737×103lb·ft−1·slug−0.5·rad−1N101.1752×102lb·ft−1·slug−0.5N2α2−5.0227×103lb·ft−1·slug−0.5·rad−2N2α2.8633×103lb·ft−1·slug−0.5·rad−1N2δe1.2465×103lb·ft−1·slug−0.5·rad−1N20−44.201×102lb·ft−1·slug−0.5Selectxt=h,V,α,θ,Q,η1,η̇1,η2,η̇2T as state variables, select ut=Φ,δeT as system inputs, and select zt=V,hT as system outputs. Substitute equations (2) and (3) into (1) to obtain the affine nonlinear model of the system as follows:
(4)ẋt=Fx,ut,t=f1x,t,f2x,t,f3x,t,⋯,f9x,tT.In equation (4), fix,t (i=1,2,⋯,9) is shown on the right side of equation (1). The equilibrium state of a hypersonic vehicle is a specific flight condition under which the line velocity, and angular velocity are constant or zero and the acceleration is zero under the given constraints. Therefore, the system equilibrium point sequence can be solved by equation (5) as follows:
(5)fix,t=0,i=1,2,⋯,9.Within the permissible range of parameters, the velocity and altitude space of the selected scheduling variables are gridded and each point on the grid is calculated and balanced according to equation (5) to obtain the equilibrium working point of the aircraft. A set of LTI systems is obtained by using the Jacobian linearization method. Finally, these linearized local systems are multivariate linearly fitted to obtain the LPV model. The equidistant gridding method is chosen in this paper. Five points are selected at equal intervals in the velocity interval 72009600ft/s, and four points are selected at equal intervals in the altitude interval 80000110000ft, totaling 20 points, as shown in Figure 1.Figure 1
Scheduling parameter mesh division.By means of parameter meshing, Jacobian linearization, and multivariate linear fitting, the LPV model of the hypersonic vehicle with respect to the scheduling parameters is shown in equation (6). The fitting results of partial LPV model parameters are shown in Figure 2, and its corresponding fitting results are listed as follows.
(6)ḣV̇α̇θ̇Q̇η̇1η¨1η̇2η¨2=00−VV00000A21A22A23−g00000A31A32A33010000000010000A51A52A5300000000000010000A7300−ω12−2ς1ω10000000000100A930000−ω22−2ς2ω2hVαθQη1η̇1η2η̇2+00B21B22B31B3200B51B520000000N2δeΦδe,(7)A21=0.005869328270053+0.000000350162259V−0.000000174371733h−0.000000000008442V·h+0.000000000022808V2+0.000000000001476h2,A22=−0.023088159404555−0.000002004091776V+0.000000734297646h+0.000000000065894V·h−0.000000000208165V2−0.000000000007750h2,A23=−4.281455848712256×102+0.000203656246563×102V+0.000109797017027×102h−0.000000002665318×102V·h+0.000000008045619×102V2−0.000000000711388×102h2,A31=0.108230120558901×10−5+0.000004011899188×10−5V−0.000003114516949×10−5h−0.000000000147221×10−5V·h+0.000000000376981×10−5V2+0.000000000030249×10−5h2,A32=−0.457306127416146×10−5+0.000041821380304×10−5V+0.000004391002007×10−5h+0.000000000788219×10−5V·h−0.000000005149628×10−5V2−0.000000000076799×10−5h2,A33=−0.003787681755457−0.000001119856725V+0.000000058391695h+0.000000000008971V·h−0.000000000041242V2−0.000000000000762h2,A51=0.296461521297026×10−4+0.000017734787744×10−4V−0.000008809870838×10−4h−0.000000000429950×10−4V·h+0.000000001166841×10−4V2+0.000000000074731×10−4h2,A52=−0.116805126192769×10−3−0.000010123132197×10−3V+0.000003712711998×10−3h+0.000000000334612×10−3V·h−0.000000001060704×10−3V2−0.000000000039242×10−3h2,A53=5.802500675695122−0.000054761923165V−0.000140973420646h+0.000000002348128V·h+0.000000030407243V2+0.000000000863671h2,A73=4.841368923917495×103−0.000005151942352×103V−0.000006725356453×103h−0.000000001382188×103V·h+0.000000005761907×103V2+0.000000000136397×103h2,A93=1.903890377392848×103+0.000018466182009×103V+0.000024105793090×103h+0.000000004954198×103V·h−0.000000020652486×103V2−0.000000000488889×103h2,B21=29.298643129461226−0.000509535839514V−0.000182979829453h−0.000000042211072V·h+0.000000177655838V2+0.000000004424883h2,B22=−2.492498015279544×102−0.000058628507583×102V+0.000067800689280×102h+0.000000001251828×102V·h−0.000000006600229×102V2−0.000000000504206×102h2,B31=−0.416859626570442×10−3+0.000021770424271×10−3V+0.000008073449568×10−3h+0.000000004341692×10−3V·h−0.000000020925782×10−3V2−0.000000000294067×10−3h2,B32≈−0.000000000019683×10−5−0.145625062719996×10−5V,B51=0.149606059591807−0.000002276698805V−0.000001007292996h−0.000000000224446V·h+0.000000000937140V2+0.000000000023361h2,B52=−0.000000167960940×10−7+0.000000000036265×10−7V+0.000000000000177×10−7h−0.251323593957002×10−7V2.Figure 2
The fitting effect of LPV model parameter.


### 2.2. Polytopic LPV Model Transformation

Based on the LPV model of the hypersonic vehicle longitudinal elastic nonlinear model, the stability analysis and controller synthesis problem of the system can be transformed into a convex optimization problem within a set of LMI frameworks which depend on the scheduling parameters. If it can be ensured that the LMI is feasible on the entire continuous variable parameter trajectory, the resulting controller can guarantee the global stability of the closed-loop system and meet the desired performance index. Since the spatial continuity of the scheduling parameters means that an infinite number of LMIs are to be solved over the entire parameter trajectory, which is obviously difficult to achieve, the scheduling parameters can be spatially meshed to solve the problem of LMI at a finite number of mesh vertices but the performance of this method is more dependent on the mesh density; besides, the determination of the mesh density not only lacks theoretical proof but also is time-consuming. The polytopic set is a convex set, which could describe the entire system. Through transformation of the polytopic model, the problem of solving infinitely many LMI problems can be attributed to solving the LMI problem at a limited number of vertices, which can greatly reduce the amount of computation. The controller obtained by this method has continuous gain throughout the workspace and has global action features ([11] [12, 22]).The system equation (6) can be expressed as follows.
(8)Sρt≔AρtBρt.where Sρt is a parameter-varying object and ρt∈Ω is a time-varying N dimensional parameter vector within the closed hypercube Ω=a1b1×a2b2×⋯×aNbN⊂ℝN. Parameter ρt can also include some elements of xt. The N+2-dimensional tensor S is constructed from LTI vertex systems Si1i2,⋯,iN.
(9)ẋt=S^⊗n=1NAnρntxtut,where the row vector Anρnt∈ℝIn (n=1,⋯,N) contains one bounded variable and continuous weighting functions αn,inρn (in=1⋯In). The weighting function αn,inρnt∈0,1 is the in weighting function defined on the nth dimension of Ω, and ρnt is the nth element of vector ρt and ∑in=1Inαn,inρnt=1. In (n=1,⋯,N) denotes the number of the weighting functions used in the nth dimension of Ω. The symbol ⊗ denotes the Kronecker product. The tensor product model transformation method based on high-order singular value decomposition is an effective LPV system polytopic model transformation method. The brief process is as follows. Firstly, mesh the variable parameters in the LPV model, and then, discretize the system to form tensors. Finally, we perform high-order singular value decomposition. The related theory of high-order singular value decomposition can refer to the literature [23–25].On the basis of the LPV model of the hypersonic vehicle in equation (6), the tensor product model transformation method is used to obtain the polyhedral description of the LPV model of the hypersonic vehicle. By uniformly dividing the altitude variable h and velocity variable V in the LPV model described in equation (6), the original system is approached in the specified area (velocity interval 72009600 ft/s, altitude interval 80000110000 ft), where the number of sampling points is set to 200×200. Using the high-order singular value method decomposition, the singular value decomposition of the expansion matrix is arranged in a descending order, where N=2, I1=I2=3. Table 2 shows the top 5 largest singular values. All 6 (the remaining singular values are relatively small or equal to 0) of the singular values are preserved. The resulting numerical model is a convincing approximation of the original model, which has 9 polytopic vertex systems. The resulting approximation model is as follows:
(10)ẋt=∑i=19wiV,hSixtut=∑i=13∑j=13wiVwjhAi,jxt+Bi,jut=∑k=19wkV,hAkxt+Bkut.Table 2
HOSVD decomposition results of the LPV model.
Expand matrixThe 1st-order singular valueThe 2nd-order singular valueThe 3rd-order singular valueThe 4th-order singular valueThe 5th-order singular valueA12.62023e+0681445.31850.496.95697e−093.80140e−10A22.62101e+0650610.4916.2435.56458e−103.59899e−10In equation (10), wkV,h=wiVwjh, k∈1,2,⋯,9, i,j∈1,2,3. Based on the high-order singular value decomposition results, the maximum modeling error of the approximation model described in equation (8) can be calculated as 0.003679 and the average modeling error is 0.002499.In order to verify the approximate equivalence between the hypersonic elastic nonlinear model described in equation (1) and the polytopic LPV model obtained by the transformation of the tensor product model described in equation (8), an equilibrium state in the parameter variation interval in Table 3 is used as an example. The system approximation matrix at the equilibrium point is obtained by the Jacobian linearization method which can locally approximate the elastic nonlinear model and the tensor product model transformation method in equation (8). The zero-pole distribution map of the system is shown in Figure 3.Table 3
Trim condition of the hypersonic vehicle.
StatusValueStatusValueStatusValueh85000 ftη11.5122Φ0.2514V7202.0808 ft/sη̇10δe11.463 degα1.5153 degη21.2144θ1.5153 degη̇20Q0 deg/sFigure 3
Zero-pole comparison diagram.
(a)
Jacobi linearization(b)
Polytopic LPV modelComparing the zero-pole positions of the two methods at the same equilibrium point, we can find that the zero-pole positions of the two methods are very close, which indicates that the polytopic LPV model based on the transformation of the tensor product model accurately reflects the dynamic characteristics of longitudinal elasticity nonlinear models of the hypersonic vehicle.The shape of the time-varying weight functionwkV,h (k∈1,2,⋯,9) is shown in Figure 4.Figure 4
Time-varying weight function.The general distribution of the nine time-varying weight functions shows that when any time-varying weight is 1, the rest of the weights are 0, which indicates that the LPV model weight function obtained by the tensor product model transformation satisfies the requirement of convex coordinates.


## 2.1. Control-Oriented Model

The longitudinal elastic nonlinear model of the hypersonic vehicle used in this paper is mainly based on the winged-cone model provided by NASA’s Langley Research Center [18–20], which takes into account the influence of the elastic modality of the aircraft. The simplified longitudinal elastic nonlinear model of the aircraft is as follows:
(1)ḣ=Vsinθ−α,V̇=Tcosα−Dm−gsinθ−α,α̇=−Tsinα−LmV+Q+gcosθ−αV,θ̇=Q,Q̇=MIyy,η¨i=−2ςiωiη̇i−ωi2ηi+Ni,i=1,2.In equation (1), h represents the aircraft altitude, V represents the aircraft velocity, α indicates the angle of attack, θ denotes the pitch angle, Q represents the pitch rate, g indicates the gravitational constant, M is the pitching moment, Iyyrepresents the moment of inertia around the body axis, T is the engine thrust, D represents resistance, L represents lift, ηi(i=1,2) represents the elastic mode, ωi(i=1,2) represents the natural frequency of the elastic mode, ςi(i=1,2) represents the damping coefficient of the elastic mode, and Nidenotes the generalized force.The approximate expressions of the above aerodynamics, moments, thrusts, and generalized forces fitted near the equilibrium point are shown in equation (2). In equation (2), ρindicates air density, S is the reference area, CLα,δe represents the lift coefficient related to the angle of attack and elevator deflection, CDα,δe indicates the drag coefficient related to the angle of attack and elevator deflection, CM,αα indicates the torque coefficient related to the angle of attack, CM,δeδe indicates the torque coefficient related to the elevator deflection, zT represents the coupling coefficient between thrust and pitching moment, CTαj represents the thrust coefficient related to αj(i=1,2,3), Niαj represents the generalized force coefficients related to αj, and N2δedenotes the influence coefficient of elevator declination to generalized force N2, and CT0 and Ni0 represent constant items in thrust and generalized forces, respectively.
(2)L≈12ρV2SCLα,δe,D≈12ρV2SCDα,δe,M≈zTT+12ρV2Sc¯CM,αα+CM,δeδe,T≈CTα3α3+CTα2α2+CTαα+CT0,N1≈N1α2α2+N1αα+N10,N2≈N2α2α2+N2αα+N2δeδe+N20,(3)ρ=ρ0exp−h−h0hs,CL=CLαα+CLδeδe+CL0,CD=CDα2α2+CDαα+CDδe2δe2+CLδeδe+CD0,CM,α=CM,αα2α2+CM,ααα+CM,α0,CM,δe=ceδe,q¯=0.5ρV2,CTα3=β1h,q¯Φ+β2h,q¯,CTα2=β3h,q¯Φ+β4h,q¯,CTα=β5h,q¯Φ+β6h,q¯,CT0=β7h,q¯Φ+β8h,q¯.The relevant force and moment coefficients are shown in equation (3), where 1/hsis the air density decay rate, CL0 and CD0 are the constant term in CLα,δeand CDα,δe, respectively, CLαi and CDαi are the ith-order coefficient of α contribution to CLα,δe and CDα,δe, respectively, CLδe is the coefficient of δe contribution to CLα,δe, CDδei is the ith-order coefficient of δe contribution to CDα,δe, CM,ααi is the ith order coefficient of α contribution to CM,αα, CM,α0 is the constant term in CM,αα, q¯ is the dynamic pressure, and Ф is the stoichiometrically normalized fuel-to-air ratio, βih,q¯ is the ith thrust fit parameter. The parameter reference values of the elastic nonlinear model are detailed in Hu et al. [21] and listed as Table 1.Table 1
Nomenclature and nominal values.
SymbolNominal valuem300 slug/ftIyy5×105lb·ftc¯17 ftρ06.7429×10−5slugs·ft−3h085000 fths2.1358×104ftg31.92 ft·s−2S17 ft2·ft−1zT8.36 ftCLα4.6773 rad−1CLδe0.76224 rad−1CL0−1.8714×10−2CDα25.8224 rad−2CDα−4.5315×10−2rad−1CDδe28.1993×10−1rad−2CDδe2.7699×10−4rad−1CD01.0131×10−2CM,αα26.2926 rad−2CM,αα2.1335 rad−1CM,α00.18979ce−1.2897 rad−1β1−3.7693×105lb·ft−1·rad−3β2−3.7225×104lb·ft−1·rad−3β32.6814×104lb·ft−1·rad−3β4−1.7277×104lb·ft−1·rad−3β53.5542×104lb·ft−1·rad−3β6−2.4216×103lb·ft−1·rad−3β76.3785×103lb·ft−1β8−1.0090×102lb·ft−1N1α21.4013×103lb·ft−1·slug−0.5·rad−2N1α4.5737×103lb·ft−1·slug−0.5·rad−1N101.1752×102lb·ft−1·slug−0.5N2α2−5.0227×103lb·ft−1·slug−0.5·rad−2N2α2.8633×103lb·ft−1·slug−0.5·rad−1N2δe1.2465×103lb·ft−1·slug−0.5·rad−1N20−44.201×102lb·ft−1·slug−0.5Selectxt=h,V,α,θ,Q,η1,η̇1,η2,η̇2T as state variables, select ut=Φ,δeT as system inputs, and select zt=V,hT as system outputs. Substitute equations (2) and (3) into (1) to obtain the affine nonlinear model of the system as follows:
(4)ẋt=Fx,ut,t=f1x,t,f2x,t,f3x,t,⋯,f9x,tT.In equation (4), fix,t (i=1,2,⋯,9) is shown on the right side of equation (1). The equilibrium state of a hypersonic vehicle is a specific flight condition under which the line velocity, and angular velocity are constant or zero and the acceleration is zero under the given constraints. Therefore, the system equilibrium point sequence can be solved by equation (5) as follows:
(5)fix,t=0,i=1,2,⋯,9.Within the permissible range of parameters, the velocity and altitude space of the selected scheduling variables are gridded and each point on the grid is calculated and balanced according to equation (5) to obtain the equilibrium working point of the aircraft. A set of LTI systems is obtained by using the Jacobian linearization method. Finally, these linearized local systems are multivariate linearly fitted to obtain the LPV model. The equidistant gridding method is chosen in this paper. Five points are selected at equal intervals in the velocity interval 72009600ft/s, and four points are selected at equal intervals in the altitude interval 80000110000ft, totaling 20 points, as shown in Figure 1.Figure 1
Scheduling parameter mesh division.By means of parameter meshing, Jacobian linearization, and multivariate linear fitting, the LPV model of the hypersonic vehicle with respect to the scheduling parameters is shown in equation (6). The fitting results of partial LPV model parameters are shown in Figure 2, and its corresponding fitting results are listed as follows.
(6)ḣV̇α̇θ̇Q̇η̇1η¨1η̇2η¨2=00−VV00000A21A22A23−g00000A31A32A33010000000010000A51A52A5300000000000010000A7300−ω12−2ς1ω10000000000100A930000−ω22−2ς2ω2hVαθQη1η̇1η2η̇2+00B21B22B31B3200B51B520000000N2δeΦδe,(7)A21=0.005869328270053+0.000000350162259V−0.000000174371733h−0.000000000008442V·h+0.000000000022808V2+0.000000000001476h2,A22=−0.023088159404555−0.000002004091776V+0.000000734297646h+0.000000000065894V·h−0.000000000208165V2−0.000000000007750h2,A23=−4.281455848712256×102+0.000203656246563×102V+0.000109797017027×102h−0.000000002665318×102V·h+0.000000008045619×102V2−0.000000000711388×102h2,A31=0.108230120558901×10−5+0.000004011899188×10−5V−0.000003114516949×10−5h−0.000000000147221×10−5V·h+0.000000000376981×10−5V2+0.000000000030249×10−5h2,A32=−0.457306127416146×10−5+0.000041821380304×10−5V+0.000004391002007×10−5h+0.000000000788219×10−5V·h−0.000000005149628×10−5V2−0.000000000076799×10−5h2,A33=−0.003787681755457−0.000001119856725V+0.000000058391695h+0.000000000008971V·h−0.000000000041242V2−0.000000000000762h2,A51=0.296461521297026×10−4+0.000017734787744×10−4V−0.000008809870838×10−4h−0.000000000429950×10−4V·h+0.000000001166841×10−4V2+0.000000000074731×10−4h2,A52=−0.116805126192769×10−3−0.000010123132197×10−3V+0.000003712711998×10−3h+0.000000000334612×10−3V·h−0.000000001060704×10−3V2−0.000000000039242×10−3h2,A53=5.802500675695122−0.000054761923165V−0.000140973420646h+0.000000002348128V·h+0.000000030407243V2+0.000000000863671h2,A73=4.841368923917495×103−0.000005151942352×103V−0.000006725356453×103h−0.000000001382188×103V·h+0.000000005761907×103V2+0.000000000136397×103h2,A93=1.903890377392848×103+0.000018466182009×103V+0.000024105793090×103h+0.000000004954198×103V·h−0.000000020652486×103V2−0.000000000488889×103h2,B21=29.298643129461226−0.000509535839514V−0.000182979829453h−0.000000042211072V·h+0.000000177655838V2+0.000000004424883h2,B22=−2.492498015279544×102−0.000058628507583×102V+0.000067800689280×102h+0.000000001251828×102V·h−0.000000006600229×102V2−0.000000000504206×102h2,B31=−0.416859626570442×10−3+0.000021770424271×10−3V+0.000008073449568×10−3h+0.000000004341692×10−3V·h−0.000000020925782×10−3V2−0.000000000294067×10−3h2,B32≈−0.000000000019683×10−5−0.145625062719996×10−5V,B51=0.149606059591807−0.000002276698805V−0.000001007292996h−0.000000000224446V·h+0.000000000937140V2+0.000000000023361h2,B52=−0.000000167960940×10−7+0.000000000036265×10−7V+0.000000000000177×10−7h−0.251323593957002×10−7V2.Figure 2
The fitting effect of LPV model parameter.


## 2.2. Polytopic LPV Model Transformation

Based on the LPV model of the hypersonic vehicle longitudinal elastic nonlinear model, the stability analysis and controller synthesis problem of the system can be transformed into a convex optimization problem within a set of LMI frameworks which depend on the scheduling parameters. If it can be ensured that the LMI is feasible on the entire continuous variable parameter trajectory, the resulting controller can guarantee the global stability of the closed-loop system and meet the desired performance index. Since the spatial continuity of the scheduling parameters means that an infinite number of LMIs are to be solved over the entire parameter trajectory, which is obviously difficult to achieve, the scheduling parameters can be spatially meshed to solve the problem of LMI at a finite number of mesh vertices but the performance of this method is more dependent on the mesh density; besides, the determination of the mesh density not only lacks theoretical proof but also is time-consuming. The polytopic set is a convex set, which could describe the entire system. Through transformation of the polytopic model, the problem of solving infinitely many LMI problems can be attributed to solving the LMI problem at a limited number of vertices, which can greatly reduce the amount of computation. The controller obtained by this method has continuous gain throughout the workspace and has global action features ([11] [12, 22]).The system equation (6) can be expressed as follows.
(8)Sρt≔AρtBρt.where Sρt is a parameter-varying object and ρt∈Ω is a time-varying N dimensional parameter vector within the closed hypercube Ω=a1b1×a2b2×⋯×aNbN⊂ℝN. Parameter ρt can also include some elements of xt. The N+2-dimensional tensor S is constructed from LTI vertex systems Si1i2,⋯,iN.
(9)ẋt=S^⊗n=1NAnρntxtut,where the row vector Anρnt∈ℝIn (n=1,⋯,N) contains one bounded variable and continuous weighting functions αn,inρn (in=1⋯In). The weighting function αn,inρnt∈0,1 is the in weighting function defined on the nth dimension of Ω, and ρnt is the nth element of vector ρt and ∑in=1Inαn,inρnt=1. In (n=1,⋯,N) denotes the number of the weighting functions used in the nth dimension of Ω. The symbol ⊗ denotes the Kronecker product. The tensor product model transformation method based on high-order singular value decomposition is an effective LPV system polytopic model transformation method. The brief process is as follows. Firstly, mesh the variable parameters in the LPV model, and then, discretize the system to form tensors. Finally, we perform high-order singular value decomposition. The related theory of high-order singular value decomposition can refer to the literature [23–25].On the basis of the LPV model of the hypersonic vehicle in equation (6), the tensor product model transformation method is used to obtain the polyhedral description of the LPV model of the hypersonic vehicle. By uniformly dividing the altitude variable h and velocity variable V in the LPV model described in equation (6), the original system is approached in the specified area (velocity interval 72009600 ft/s, altitude interval 80000110000 ft), where the number of sampling points is set to 200×200. Using the high-order singular value method decomposition, the singular value decomposition of the expansion matrix is arranged in a descending order, where N=2, I1=I2=3. Table 2 shows the top 5 largest singular values. All 6 (the remaining singular values are relatively small or equal to 0) of the singular values are preserved. The resulting numerical model is a convincing approximation of the original model, which has 9 polytopic vertex systems. The resulting approximation model is as follows:
(10)ẋt=∑i=19wiV,hSixtut=∑i=13∑j=13wiVwjhAi,jxt+Bi,jut=∑k=19wkV,hAkxt+Bkut.Table 2
HOSVD decomposition results of the LPV model.
Expand matrixThe 1st-order singular valueThe 2nd-order singular valueThe 3rd-order singular valueThe 4th-order singular valueThe 5th-order singular valueA12.62023e+0681445.31850.496.95697e−093.80140e−10A22.62101e+0650610.4916.2435.56458e−103.59899e−10In equation (10), wkV,h=wiVwjh, k∈1,2,⋯,9, i,j∈1,2,3. Based on the high-order singular value decomposition results, the maximum modeling error of the approximation model described in equation (8) can be calculated as 0.003679 and the average modeling error is 0.002499.In order to verify the approximate equivalence between the hypersonic elastic nonlinear model described in equation (1) and the polytopic LPV model obtained by the transformation of the tensor product model described in equation (8), an equilibrium state in the parameter variation interval in Table 3 is used as an example. The system approximation matrix at the equilibrium point is obtained by the Jacobian linearization method which can locally approximate the elastic nonlinear model and the tensor product model transformation method in equation (8). The zero-pole distribution map of the system is shown in Figure 3.Table 3
Trim condition of the hypersonic vehicle.
StatusValueStatusValueStatusValueh85000 ftη11.5122Φ0.2514V7202.0808 ft/sη̇10δe11.463 degα1.5153 degη21.2144θ1.5153 degη̇20Q0 deg/sFigure 3
Zero-pole comparison diagram.
(a)
Jacobi linearization(b)
Polytopic LPV modelComparing the zero-pole positions of the two methods at the same equilibrium point, we can find that the zero-pole positions of the two methods are very close, which indicates that the polytopic LPV model based on the transformation of the tensor product model accurately reflects the dynamic characteristics of longitudinal elasticity nonlinear models of the hypersonic vehicle.The shape of the time-varying weight functionwkV,h (k∈1,2,⋯,9) is shown in Figure 4.Figure 4
Time-varying weight function.The general distribution of the nine time-varying weight functions shows that when any time-varying weight is 1, the rest of the weights are 0, which indicates that the LPV model weight function obtained by the tensor product model transformation satisfies the requirement of convex coordinates.


## 3. Robust Controller Design

The LPV model of a hypersonic vehicle is as follows.(11)ẋt=Aθtxt+Bθtut+Dwt,yt=Cxt.In equation (11), θt=VhT denotes the scheduling variables, Aθt and Bθt both are the corresponding parameter-dependent system matrix, and the numerical expression is described in equation (6). C is the observation matrix, D is the interference matrix, and wt denotes time-varying exogenous disturbance. The reference velocity and altitude of the aircraft are recorded as rt=VcomhcomT and the control output is recorded as yt=VthtT. The controller design goal of the hypersonic vehicle not only needs to ensure the stability of the aircraft but also requires that the output of the aircraft can track a given command signal in real time. Therefore, the design goal of the hybrid target robust controller is to design an effective controller within the whole flight envelope ensuring that the output tracking error of the aircraft is 0.
(12)limt⟶∞yt−rt=0.In order to evaluate the system’s steady-state error more effectively, the following tracking error integral term is introduced to eliminate the steady-state error in the feedback control strategy.(13)xet=∫0tyτ−rτdτ.According to the tracking error integral term and combining polytopic description and the LPV model of the hypersonic vehicle, the extended system can be obtained as follows.(14)Ẋt=A^θtXt+B1rt+B^θtut+D^wt,Yt=C^Xt,where
(15)Xt=xtxetT,A^θt=Aθt0C0=∑k=19wkV,hA^k,B^θt=Bθt0=∑k=19wkV,hB^k,A^k=Ak0C0,B^k=Bk0,C^=C0B1=0−IT,D^=DT0T.According to the parallel distributed compensation (PDC), the designed state feedback controller is as follows.(16)ut=∑k=19wkV,hKkXt=KθtXt.In equation (16), Kθt=∑k=19wkV,hKk. Substituting equation (16) into (14), the expression for the closed-loop system of a hypersonic vehicle can be obtained as follows.
(17)Ẋt=A^θt+B^θtKθtXt+B1rt+D^wt=∑i=19∑j=19wiV,hwjV,hA^i+B^iKjXt+B1rt+D^wt.The output tracking problem of the hypersonic vehicle is transformed into the robust stability problem of the above extended system. When the closed-loop system described in equation (17) is robustly stable, the command tracking target of the original system described in equation (9) can be achieved.


## 4. Numerical Simulation

Under conditions of cruising flight of the hypersonic vehicle (altitude 85000 ft, velocity 7702 ft/s), a robust model predictive controller for state feedback of mixed performance indicators was designed to achieve tracking control of velocity commands and altitude commands [26]. In the actual control of the aircraft, the command must comprehensively consider the physical constraints of the aircraft such as structural stiffness and payload. In general, it is assumed that the given velocity and altitude reference command signals are generated by the second-order filters.
(18)Hs=wn2s2+2ζwns+wn2.In equation (18), the velocity instruction signal filtering system natural frequency wV=0.0195 and damping coefficient ζV=0.9. The altitude command signal filtering system natural frequency wh=0.0219 and damping coefficient ζh=0.9. It is assumed that the initial equilibrium state of the longitudinal dynamics of the aircraft is as follows. The rigid state V0=7702ft/s, h0=85000ft, α0=1.6262deg, θ0=1.6262deg, and Q=0. The elastic mode η1=0.9127, η̇1=0, η2=0.7623, and η̇2=0. The control input Φ=0.2666 and δe=12.4978deg. The control input satisfies the hard constraints of the actuator; the engine fuel equivalence ratio and the range of the elevator declination angle are 0.1≤Φ≤1.4 and ‐30≤δe≤30.The discretization sampling timeTs of the continuous system is 0.01 s. The aerodynamic parameters in the aircraft model are mainly obtained by computational fluid dynamics simulations or wind tunnel experiments. Besides, the unknown external environment and complex dynamic influences such as dynamic pressure effects, viscous effects, lean gas effects, and atmospheric turbulence may be encountered in the flight process. Therefore, the abovementioned factors are assumed as uncertainties of model inertial parameters and aerodynamic parameters. Similar to the literature [27], the value of parameters lies within ±10% of the nominal value such as m=m1+△m where △m≤0.07 and ρ=ρ1+△ρ where △ρ≤0.08. Compared with the ones in Qin et al. [27], the uncertainties are more aggressive and energy-bounded external disturbance is considered. When the desired velocity Ve=9200ft/s and desired altitude he=104000ft, the simulation results using the state feedback controller design method are shown in Figures 5–9. The simulation results are compared with the results obtained through Qin’s method [27] and Cai’s method [28]. In the sequence, results obtained through the proposed method, Qin’s method, and Cai’s method are depicted by a red line marked by circles, blue dashed line marked by circles, and green dashed line marked by asterisks. Since the tracking curves are quite close, thus, only the tracking errors are given.Figure 5
Velocity tracking performance and tracking error.Figure 6
Altitude tracking performance and tracking error.Figure 7
Angle of attack.Figure 8
Pitch angle.Figure 9
Control input.In Figures5 and 6, it can be seen that the robust controller can ensure that a hypersonic vehicle achieves the fast and accurate velocity and altitude command tracking during the maneuvering flight of the large flight envelope. The above controller has a smaller tracking error and a relatively smooth state, which can effectively suppress the influence of serious aerodynamic parameter uncertainties and elastic modes in the flight process. Figures 7 and 8 show the variation curves of the angle of attack and pitch angle in the whole process, respectively. In order to achieve a high degree of climb while increasing the velocity, the hypersonic vehicle in the initial trim state adds engine fuel input to the controller output signal and the elevator is pulled up until the stable command tracking is achieved again to reach the equilibrium state. The angle of attack and the pitch angle will be stabilized in the new trim state. The curve of the control input shown in Figure 9 indicates that the engine fuel equivalence ratio and the elevator deflection angle all satisfy the system control input hard constraints.In summary, in the case of external disturbances, parameter uncertainties, and elastic modal couplings during the maneuvering flight, the designed controller has strong robustness and ensures that the system can quickly track the command signal.


## 5. Conclusion

Aiming at a type of hypersonic vehicle longitudinal elasticity model disclosed in the existing literature, altitude and velocity are selected as the gain-scheduled variables, through a series of sequential steps such as meshing within the range of parameter variation, finding equilibrium points, Jacobian linearization, and multivariate linear fitting; the continuous-time LPV model of the aircraft is established. Then, using the tensor product model transformation method, the general LPV model is transformed into a polytopic LPV system model with 9 vertices and the rationality of the LPV model is verified. At the same time, regarding the maneuvering flight command tracking problem for hypersonic vehicles, a state feedback controller is used to translate the command tracking control problem into a robust stability problem for the corresponding extended closed-loop system combined with the obtained polytopic LPV system model. Simulation results demonstrate that the obtained polytopic LPV model based on high-order singular value decomposition in this paper is effective in solving the maneuvering flight command tracking and attitude stabilization control of general complex nonlinear systems, such as hypersonic vehicles. However, since the linearization process lose some beneficial nonlinearities, some gain perturbations due to finite word length or time delays of state or input might bring unpredictable complexities. LPV nonfragile control of air-breathing hypersonic vehicles with time delays needs further studies.



---

*Source: 7570059-2021-12-03.xml*