TBCC组合动力发动机基于区域极点配置的推力鲁棒控制

Robust Thrust Control of TBCC Combined Cycle Engine Based on Regional Pole Assignment

摘要: 研究了某涡轮基组合动力发动机XTER在给定航迹下的全速域推力跟踪问题。基于模型扩展思想，在原有的基于系统辨识的小偏差模型上考虑了外扰和喷油惯性环节，建立了新的面向控制的LPV模型，将原推力跟踪问题转化为推力误差镇定问题。为有效抑制外界干扰并使发动机可以快速无超调跟踪目标推力，将区域极点配置引入到状态反馈鲁棒H_∞控制中，设计了一种兼顾稳态特性和暂态性能的推力跟踪控制器。最后将得到的控制算法加入到非线性模型上验证其控制效果。与已有的LADRC推力跟踪控制算法相比，本文所设计的控制器跟踪速度更快，鲁棒性能更好。

Abstract: The full-speed range thrust tracking problem of a certain TBCC combined cycle engine XTER under a given flight path was studied. Based on the model extension idea, the external disturbance and fuel injection inertia link were considered on the original small deviation model based on system identification, and a new LPV model oriented to control was established. The original thrust tracking problem was transformed into the thrust error stabilization problem. To effectively suppress external disturbances and enable the engine to quickly and without overshoot track the target thrust, the regional pole assignment was introduced into the state feedback robust control, and a thrust tracking controller that takes into account both steady-state characteristics and transient performance was designed. Finally, the obtained control algorithm was added to the nonlinear model to verify its control effect. Compared with the existing LADRC thrust tracking control algorithm, the controller designed in this paper has a faster tracking speed and better robust performance.

1. 引言

高超声速飞行被誉为继螺旋桨和喷气推进飞行后航空史上的第三次“革命”[1]，是21世纪航空航天技术制高点，其核心价值在于突破传统飞行器的速度与空域限制，实现快速全球打击、侦察与太空接轨等战略目标[2]。这一技术实现的关键是动力系统的创新，而涡轮基组合循环发动机(Turbine-Based Combined Cycle, TBCC)有着宽速域、适应性强、高比冲、可重复使用等优势，是当前最接近工程化应用的高超声速飞行动力解决方案。

TBCC组合动力发动机按布局分为串联式和并联式。并联式TBCC相较于串联式TBCC而言研制难度较低，可在现有的涡轮发动机技术以及成熟的冲压发动机技术的基础上进行实现[3]。目前已有多个国家针对并联式TBCC的建模、控制进行了研究。截至目前，并联式TBCC还未发展至成熟的工程应用阶段。飞发一体化、模态转换推力鸿沟、跨声速推力裕度不足等关键技术，一直是其工程化需要突破的难点[4]。为解决并联式TBCC发动机所存在的“推力鸿沟”的问题，厦门大学针对0~6 Ma超声速飞行条件，设计了一种以涡喷发动机为基本动力，集成引射–亚燃冲压发动机、超燃冲压发动机的四通道三动力组合发动机(Xiamen Turbine Ejector-Ramjet Combined Cycle, XTER) [5]。该方案以技术成熟度较高、模态转换性能较好且研制周期较短等特点，为我国涡轮基组合动力领域的发展提供一种新思路。目前该方案已完成整机建模、最优航迹规划和最优模态转换点求解等科学问题研究以及整机风洞实验。

在控制方面，TBCC发动机与传统的单一动力发动机相比，要兼顾水平起降、爬升和模态转换以及高速巡航航程的推力跟踪需求，因此更加复杂。如何改善发动机的安全性能，实现发动机在整个工作速域内稳定工作是TBCC发动机控制研究的重点问题。近年来，许多学者针对这些TBCC发动机的控制问题，进行了相应的学术研究。聂聆聪研究了一种串联式TBCC发动机模态转换的多变量控制方法，实现了模态转换中推力及流量的连续变化[6]。伍谦研究了TBCC涡轮模式的起动问题，对TBCC的工作特性进行分析，确定模态转换点，将工作模式划分五类[7]。对起动模式设计PI控制器，并对涡轮级间温度进行控制，避免了“超温启动”；其他模式采用增广LQR控制算法，合理地控制发动机内推力和涡轮前总温。许海琦针对TBCC组合进气道的喘振不起动问题进行了数值模拟和仿真分析，总结了其通用的耦合触发及影响规律，和发动机模型进行联合仿真得到进气道的安全边界。设计了TBCC模态转换过程中可以保证推力输出和进气道稳定裕度的闭环控制器[8]。仇小杰针对小型串联式TBCC发动机提出一种基于序列二次规划方法的模态转换控制器，实现了平滑推力下的模态转换[9]。对于XTER组合动力发动机的控制问题，文献[5]利用LQ最优控制理论完成了XTER发动机模态转换点的求解与各通道分流板开度规律的设计，得到了最佳模态切换点以及涡喷、引射–亚燃分流板开度的变化规律。同时对各通道设计了LADRC推力跟踪控制器，但在控制器设计时存在着未考虑模型不确定性的问题且控制器跟踪速度较慢。

在XTER发动机已有的控制研究基础上，为有效抑制发动机工作过程中受到的外部扰动，本文采用状态反馈鲁棒$H_{\infty }$ 控制方法进行推力跟踪控制。同时，为了兼顾XTER发动机在推力跟踪控制中的暂态特性，使发动机可以快速跟踪上目标推力，采用闭环系统区域极点配置技术。通过将这两种方法结合在一起，设计出一种既考虑鲁棒性又兼顾暂态特性的多目标控制器。这个多目标控制问题的求解比较困难，为便于求解，将其转化为线性矩阵不等式(Linear Matrix Inequality, LMI)的求解问题[10]。通过MATLAB中的LMI工具箱编写程序进行求解，得到控制器参数。利用线性系统理论设计了多目标控制器后，在XTER发动机气动热力学非线性模型上进行仿真验证。

2. 预备知识

本节先简单介绍XTER发动机的非线性模型和面向控制的LPV模型，以及鲁棒$H_{\infty }$ 状态反馈镇定和区域极点配置相关的定义和引理。XTER发动机的五种工作状态在文献[5]中均有介绍，本文选取较为复杂的涡喷到引射–亚燃模态转换阶段的LPV模型为例来说明控制器的设计过程。

2.1. XTER发动机气动热力学非线性仿真模型

XTER发动机是以涡喷发动机为基本动力，集成超燃冲压发动机、引射–亚燃冲压发动机的四通道三动力组合发动机，这四个通道共用一个公共的进气道和尾喷管，其结构如图1(a)所示。XTER发动机结构设计的关键在于引射火箭，该设计解决了推力鸿沟的问题且通过引射气流改善了亚燃发动机的入口气流参数。XTER发动机的公共进气道为三维内转式进气道，结构如图1(b)所示。入口处铰接有两个分流板，其中涡喷通道入口为涡喷分流板，设其开度为$\alpha$ 。同样，引射–亚燃通道入口铰接引射–亚燃分流板，设其开度为$\beta$ 。各个发动机的状态量均可采用传感器测量的方式来直接获取或通过模型进行合理估计。

Figure 1. (a) Structure diagram of XTER combined power engine; (b) Top view and side view of XTER combined air inlet

图1. (a) XTER组合动力发动机结构示意图；(b) XTER组合进气道俯视图及侧视图

文献[5]采用部件级建模方式，将气动热力学方程建立在各部件气动热力学关系的基础上，用这些方程代替发动机真正的工作部件。再由气动热力学和发动机工作必须遵守的转子动力学平衡关系串联起各个部件，构成一个共同的工作方程组。XTER发动机的所有参数都可以通过求解这些非线性方程组来得到，从而模拟出真实发动机的工作过程。其中，涡喷通道通过转子动力学效应与喷油惯性环节体现其动态特性；引射–亚燃通道的动态特性通过火箭燃烧室与亚燃燃烧室的容积动力学效应与喷油惯性环节描述；超燃通道的动态特性由超燃燃烧室的喷油惯性环节体现，最终得到的这个模型称它为非线性仿真模型。

完成非线性仿真模型的建立以后，在此基础上进一步研究XTER发动机在给定航迹下的推力跟踪控制问题。得到符合控制需求的控制算法后，非线性仿真模型作为验证平台来验证控制算法的正确性和合理性[11]。

2.2. XTER发动机面向控制的LPV模型

XTER发动机各个通道的工作区间如图2所示。其中$M{a}_{\text{c1}}、M{a}_{\text{c2}}、M{a}_{\text{c3}}、M{a}_{\text{c4}}$ 为各通道的模态转换点，且在文献[5]中已求得。

Figure 2. Working Mach range of each mode of XTER

图2. XTER各通道工作马赫数范围

LPV模型是指在被控对象每一个工作段上都选取若干个工作点，在每个工作点上用线性模型来代替非线性模型，最后将所有工作点处的线性模型对调度参数进行插值或拟合。XTER发动机基于控制的LPV模型的建模难点在于模态转换阶段，此阶段需对两个共同工作的子系统通过状态加维的方式建立。本文以涡喷到引射–亚燃模态转换阶段即工作马赫数区间为$\left[M{a}_{\text{c1}},M{a}_{\text{c2}}\right]$ 为例，在文献[5]的基础上考虑了XTER发动机的喷油惯性环节，给出其基于控制的LPV模型。

文献[5]中，该工作马赫数区间内涡喷和引射–亚燃子动力同时工作，LPV模型描述为

$\left[\begin{array}{c}\Delta {\dot{n}}_1\\ \Delta {\dot{T}}_1^{*}\\ \Delta {\dot{T}}_2^{*}\\ \Delta {\dot{F}}_2\end{array}\right]=\left[\begin{array}{cccc}{a}_1\left(Ma\right)& a_2\left(Ma\right)& 0& 0\\ a_3\left(Ma\right)& a_4\left(Ma\right)& 0& 0\\ 0& 0& a_5\left(Ma\right)& a_6\left(Ma\right)\\ 0& 0& a_7\left(Ma\right)& a_8\left(Ma\right)\end{array}\right]\left[\begin{array}{c}\Delta n_1\\ \Delta T_1^{*}\\ \Delta T_2^{*}\\ \Delta F_2\end{array}\right]+\left[\begin{array}{cc}{b}_1\left(Ma\right)& 0\\ b_2\left(Ma\right)& 0\\ 0& b_3\left(Ma\right)\\ 0& b_4\left(Ma\right)\end{array}\right]\left[\begin{array}{c}\Delta W{f}_1\\ \Delta W{f}_2\end{array}\right]$ (1)

$\Delta y_{t1}=\Delta y_1+\Delta y_2=\left[\begin{array}{cccc}{c}_1\left(Ma\right)& c_2\left(Ma\right)& 0& 1\end{array}\right]\left[\begin{array}{c}\Delta n_1\\ \Delta T_1^{*}\\ \Delta T_2^{*}\\ \Delta F_2\end{array}\right]$ (2)

式中，$\Delta n_1$ 为涡轮转速增量、$\Delta T_1^{*}$ 为涡轮前总温增量，$\Delta T_2^{*}$ 为亚燃燃烧室总温增量，$\Delta F_2$ 引射–亚燃通道推力增量为输出量，这4个量为状态量。$\Delta W{f}_1$ 为涡喷通道燃油流量增量，$\Delta W{f}_2$ 为引射–亚燃通道燃油流量增量，这2个量为模型的输入量。$\Delta y_1$ 为涡喷通道推力输出增量，$\Delta y_2$ 为引射–亚燃通道推力输出增量，模型的输出量$\Delta y_{t1}=\Delta y_1+\Delta y_2$ 是涡喷到引射–亚燃过渡模态XTER发动机的总推力增量。

将涡喷和引射–亚燃发动机喷油装置的时间常数考虑为0.1 s。令涡喷发动机和引射–亚燃发动机的喷油变化量$\Delta W{f}_1$ 和$\Delta W{f}_2$ 也作为状态变量[12]，得到XTER发动机的增广LPV模型为

$\{\begin{array}{l}{\dot{x}}_p=A_{\text{p}}\left(Ma\right)x_{\text{p}}+B_{\text{p}}\left(Ma\right)u_{\text{p}}\\ y_{\text{p}}=C_{\text{p}}\left(Ma\right)x_{\text{p}}+D_{\text{p}}\left(Ma\right)u_{\text{p}}\end{array}$ (3)

其中

$x_{\text{p}}={\left[\begin{array}{cccccc}\Delta n_1& \Delta T_1^{*}& \Delta T_2^{*}& \Delta F_2& \Delta W{f}_1& \Delta W{f}_2\end{array}\right]}^{\text{T}},$ $y_{\text{p}}=\Delta y_{t1},$ $u_{\text{p}}=\left[\begin{array}{c}\Delta W{f}_{\text{r1}}\\ \Delta W{f}_{\text{r2}}\end{array}\right],$ $W{f}_{\text{r1}}，\text{}W{f}_{\text{r2}}$ 为喷油指令

$A_{\text{p}}=\left[\begin{array}{cccccc}{a}_1\left(Ma\right)& a_2\left(Ma\right)& 0& 0& b_1\left(Ma\right)& 0\\ a_3\left(Ma\right)& a_4\left(Ma\right)& 0& 0& b_2\left(Ma\right)& 0\\ 0& 0& a_5\left(Ma\right)& a_6\left(Ma\right)& 0& b_3\left(Ma\right)\\ 0& 0& a_7\left(Ma\right)& a_8\left(Ma\right)& 0& b_4\left(Ma\right)\\ 0& 0& 0& 0& -10& 0\\ 0& 0& 0& 0& 0& -10\end{array}\right],$ $B_{\text{p}}={\left[\begin{array}{cccccc}0& 0& 0& 0& 10& 0\\ 0& 0& 0& 0& 0& 10\end{array}\right]}^{\text{T}},$ $C_{\text{p}}=\left[\begin{array}{cccc}{c}_1\left(Ma\right)& c_2\left(Ma\right)& 0& 1\end{array}\begin{array}{cc}0& 0\end{array}\right],$ $D_{\text{p}}=[\begin{array}{cc}0& 0\end{array}]$ 。

2.3. 鲁棒H_∞状态反馈镇定

鲁棒$H_{\infty }$ 控制是一种基于频域优化的控制策略，通过最小化系统传递函数的$H_{\infty }$ 范数(即最坏扰动下的最大能量增益)，确保系统在模型不确定性、参数摄动和外部干扰下仍能保持稳定性和预设性能指标。其本质是在频域内设计控制器，使闭环系统对各类扰动的敏感性降至最低。

鲁棒$H_{\infty }$ 控制的广义系统如图3所示：

Figure 3. Generalized systems of robust $H_{\infty }$ control

图3. 鲁棒$H_{\infty }$ 控制的广义系统

被控系统$P\left(s\right)$ 为如式(4)所示的线性时不变连续时间系统

$P\left(s\right):\{\begin{array}{l}\dot{x}=Ax+B_{1}w+B_{2}u\\ z=C_{1}x+D_{11}w+D_{12}u\\ y=C_{2}x+D_{21}w+D_{22}u\end{array}$ (4)

其中，$x\in R^{\text{n}}$ 是状态向量，$u\in R^{\text{m}}$ 是控制输入，$y\in R^{\text{p}}$ 是测量输出，$z\in R^{\text{r}}$ 是感兴趣的系统被调输出，$w\in R^{\text{q}}$ 是外部扰动，这里考虑的外部扰动是不确定的，但具有有限能量。

在使用鲁棒$H_{\infty }$ 控制器镇定系统时，需要先将被控系统转化为上述标准鲁棒模型的形式，然后根据控制需求进行控制器的设计。

对系统(4)设计状态反馈控制器，需要求解出状态反馈矩阵$K$ 。此时控制量$u$ 为

$u=Kx$ (5)

外部干扰$w$ 到被控输出$z$ 的传递函数$T_{zw}\left(s\right)$ 的$H_{\infty }$ 范数可以表示为

${\Vert T_{zw}\left(s\right)\Vert }_{\infty }={\Vert \left(C_1+D_{12}K\right){\left(sI-A-B_{2}K\right)}^{-1}{B}_1+D_{11}\Vert }_{\infty }$ (6)

$\gamma$ 次优$H_{\infty }$ 控制标准问题(P1)：给定一个标量$\gamma >0$ ，对系统(4)设计如式(5)所示的控制器，使传递函数$T_{zw}\left(s\right)$ 的$H_{\infty }$ 范数满足如下条件

${\Vert T_{zw}\left(s\right)\Vert }_{\infty }<\gamma$ (7)

称满足式(7)的控制器$u=Kx$ 为一个$\gamma$ 次优状态反馈鲁棒$H_{\infty }$ 控制器。

引理1 [13]：对于被控系统(4)，$\gamma$ 次优$H_{\infty }$ 控制标准问题(P1)可解的充要条件是存在一个对称正定矩阵$W$ 和任意合适维数的矩阵$V$ ，使得如下LMI

$\left[\begin{array}{ccc}AW+B_{2}V+{\left(AW+B_{2}V\right)}^{\text{T}}& B_1& {\left(C_{1}W+D_{12}V\right)}^{\text{T}}\\ *& -\gamma I& D_{11}^{\text{T}}\\ *& *& -\gamma I\end{array}\right]<0$ (8)

成立，使得$\gamma$ 次优$H_{\infty }$ 控制标准问题(P1)可解的控制器为$u=V{W}^{-1}x$ 。

2.4. 区域极点配置

定义1 [13]：对复平面中的区域$D$ ，如果存在一个对称矩阵$L\in R^{\text{m}\times \text{m}}$ 和矩阵$M\in R^{\text{m}\times \text{m}}$ ，使得

$D=\left\{s\in C:L+sM+\overline{s}{M}^{\text{T}}<\text{0}\right\}$ (9)

则称$D$ 是一个线性矩阵不等式区域(简记为LMI区域)。矩阵值函数

$f_{\text{D}}\left(z\right)=L+sM+\overline{s}{M}^{\text{T}}$ (10)

称为LMI区域$D$ 的特征函数。

定义2 [13]：对复平面中给定的LMI区域$D$ 和矩阵$A\in R^{\text{n}\times \text{n}}$ ，如果矩阵$A$ 的所有特征值都位于区域$D$ 中，即$\sigma \left(A\right)\subset D$ ，则称矩阵$A$ 是D-稳定的。

引理2 [13]：给定由式(8)描述的LMI区域$D$ ，则矩阵$A\in R^{\text{n}\times \text{n}}$ 是D-稳定的充分必要条件是存在一个对称正定矩阵$X\in R^{\text{n}\times \text{n}}$ ，使得

$M_{\text{D}}\left(A,X\right)<0$ (11)

其中，$M_{\text{D}}\left(A,X\right)=L\otimes X+M\otimes \left(AX\right)+M^{\text{T}}\otimes {\left(AX\right)}^{\text{T}}$

引理3 [13]：给定两个LMI区域$D_1$ 和$D_2$ ，矩阵$A$ 同时是$D_1$ -稳定和$D_2$ -稳定的充分必要条件是存在一个对称正定矩阵$X$ ，使得

$M_{{\text{D}}_{\text{1}}}\left(A,X\right)<0,M_{{\text{D}}_{\text{2}}}\left(A,X\right)<0$ (12)

3. 控制器设计

本文所研究的问题是在考虑XTER发动机喷油环节、模型的不确定性和有限能量外扰的情况下，进行XTER发动机在已知航迹下的推力跟踪控制。同时，为了保证控制系统具有一定的暂态性能，将每个工作点处闭环系统的极点配置到如图4中阴影区域所示位置，称这个位置为极点配置目标区域[14]。本文每隔0.1 Ma作为一个工作点设计控制器，最后通过插值的方法实现全速域范围内控制器的调度。

Figure 4. Pole assignment target area

图4. 极点配置目标区域

为了完成这个控制目标，采用积分补偿的方法，将偏差的积分也作为增广系统的状态变量，建立新的LPV模型，将推力跟踪问题转化为推力误差镇定问题。

针对图5所示的控制系统结构进行控制器设计。其中$r$ 为参考输出，$d$ 为外部干扰信号，输出偏差为$e=r-\left(y_{\text{p}}+d\right)$ ，偏差的积分为$\xi ={\displaystyle {\int }_0^t\left(r\left(\tau \right)-\left(y_{\text{p}}\left(\tau \right)+d\left(\tau \right)\right)\right)\text{d}\tau }$ ，考虑外扰的LPV模型为

$\{\begin{array}{l}\dot{\overline{x}}=\overline{A}\left(Ma\right)\overline{x}+{\overline{B}}_1\left(Ma\right)\overline{w}+{\overline{B}}_2\left(Ma\right)\overline{u}\\ \overline{z}={\overline{C}}_1\left(Ma\right)\overline{x}+{\overline{D}}_{11}\left(Ma\right)\overline{w}+{\overline{D}}_{12}\left(Ma\right)\overline{u}\\ \overline{y}={\overline{C}}_2\left(Ma\right)\overline{x}+{\overline{D}}_{21}\left(Ma\right)\overline{w}+{\overline{D}}_{22}\left(Ma\right)\overline{u}\end{array}$ (13)

其中

$\overline{x}=\left[\begin{array}{c}{x}_{\text{p}}\\ \xi \end{array}\right],$$\overline{z}=y_{\text{p}}$为被控输出，$\overline{u}=u_{\text{p}}，$$\overline{w}=\left[\begin{array}{c}d\\ r\end{array}\right],$ $\overline{y}=\left[\begin{array}{c}\xi \\ y_{\text{d}}\end{array}\right]$，$y_{\text{d}}=y_{\text{p}}+d$ ，

$\overline{A}\left(Ma\right)=\left[\begin{array}{cc}{A}_{\text{p}}\left(Ma\right)& 0\\ -C_{\text{p}}\left(Ma\right)& 0\end{array}\right],$${\overline{B}}_1\left(Ma\right)=\left[\begin{array}{cc}0& 0\\ -I& I\end{array}\right],$ ${\overline{B}}_2\left(Ma\right)=\left[\begin{array}{c}{B}_{\text{p}}\left(Ma\right)\\ -D_{\text{p}}\left(Ma\right)\end{array}\right],$${\overline{C}}_1\left(Ma\right)\overline{=}\left[\begin{array}{cc}{C}_{\text{p}}\left(Ma\right)& 0\end{array}\right]，$${\overline{D}}_{11}\left(Ma\right)=\left[0\right],$ ${\overline{D}}_{12}\left(Ma\right)=\left[D_{\text{p}}\left(Ma\right)\right],$${\overline{C}}_2\left(Ma\right)=\left[\begin{array}{cc}0& I\\ C_{\text{p}}\left(Ma\right)& 0\end{array}\right],$${\overline{D}}_{21}\left(Ma\right)=\left[\begin{array}{cc}0& 0\\ I& 0\end{array}\right],$ ${\overline{D}}_{22}\left(Ma\right)=\left[\begin{array}{c}0\\ D_{\text{p}}\left(Ma\right)\end{array}\right]$。

Figure 5. Block diagram of state feedback control system

图5. 状态反馈控制系统框图

这样本文待求解的推力跟踪控制问题的数学描述如下。

推力跟踪控制问题(P2)：考虑LPV模型(13)，基于图5所示的控制系统框图，设计如下状态反馈控制器

$u\left(Ma\right)=K\left(Ma\right)\overline{x}\left(Ma\right)$ (14)

使得在有限能量范围的外扰下，(13)和(14)构成的闭环系统

$\{\begin{array}{l}{\dot{\overline{x}}}_{\text{cl}}={\overline{A}}_{\text{cl}}\left(Ma\right){\overline{x}}_{\text{cl}}+{\overline{B}}_{\text{cl}}\left(Ma\right)\overline{w}\\ {\overline{z}}_{\text{cl}}={\overline{C}}_{\text{cl}}\left(Ma\right){\overline{x}}_{\text{cl}}+{\overline{D}}_{\text{cl}}\left(Ma\right)\overline{w}\end{array}$ (15)

满足

(1) $H_{\infty }$ 性能指标：闭环系统(15)在每个工作点处渐近稳定，且从$\overline{w}$ 到${\overline{z}}_{\text{cl}}$的闭环传递函数$T_{wz\text{cl}}\left(s\right)$ 满足

${\Vert T_{wz\text{cl}}\left(s\right)\Vert }_{\infty }<\gamma$ (16)

即输出推力误差对于外界干扰具有指定的抑制度$\gamma$ 。

(2) 极点配置约束：闭环系统(15)在每个工作点处的极点都位于极点配置目标区域$S\left(\alpha ,r,\theta \right)$ 内，即保证系统在每个工作点处都具有最小的衰减度$\alpha 、$ 最小阻尼比$\zeta =\cos \theta$ 和一个最大自然频率${\omega }_n=r\sin \theta$ 。

其中，${\overline{x}}_{\text{cl}}=\overline{x}，$${\overline{z}}_{\text{cl}}=\overline{z}，$${\overline{A}}_{\text{cl}}\left(Ma\right)=\overline{A}\left(Ma\right)+{\overline{B}}_2\left(Ma\right)K\left(Ma\right)，$ ${\overline{B}}_{\text{cl}}\left(Ma\right)={\overline{B}}_1\left(Ma\right)，$${\overline{C}}_{\text{cl}}\left(Ma\right)={\overline{C}}_1\left(Ma\right)+{\overline{D}}_{12}\left(Ma\right)K\left(Ma\right)，$${\overline{D}}_{\text{cl}}\left(Ma\right)={\overline{D}}_{11}\left(Ma\right)$。

对XTER某个工作点处推力跟踪控制问题P2的求解，根据引理1，$\gamma$ 次优状态反馈鲁棒$H_{\infty }$ 控制器可以通过求解如下关于对称正定矩阵$X_1$ 和任意合适维数的矩阵$Y_1$ 的LMI得到

$\left[\begin{array}{ccc}\overline{A}{X}_1+{\overline{B}}_2{Y}_1+{\left(\overline{A}{X}_1+{\overline{B}}_2{Y}_1\right)}^{\text{T}}& {\overline{B}}_1& {\left({\overline{C}}_1{X}_1+{\overline{D}}_{12}{Y}_1\right)}^{\text{T}}\\ *& -\gamma I& {\overline{D}}_{11}^{\text{T}}\\ *& *& -\gamma I\end{array}\right]<0$ (17)

且控制器为$u=Y_1{X}_1^{-1}\overline{x}$ 。

考虑到要兼顾被控系统的暂态特性，将闭环系统配置到如图4所示的位置中。这个极点配置的目标区域可以看成是一个具有$\alpha$ 稳定度的半平面区域$D_{\alpha }$ 、一个半径为$r$ 圆盘和一个圆锥扇形$S\left(0,0,\theta \right)$ 这三个LMI区域的交。依据引理2和引理3，若存在适当维数的对称正定矩阵$X_2$ ，使得不等式(18)~(20)可解，那么闭环系统(15)的极点都位于该目标区域内。

即闭环系统极点配置问题转化为矩阵不等式(18)~(20)的求解。

$\left(\overline{A}+{\overline{B}}_{2}K\right)X_2+X_2{\left(\overline{A}+{\overline{B}}_{2}K\right)}^{\text{T}}+2\alpha X_2<0$ (18)

$\left[\begin{array}{cc}-r{X}_2& \left(\overline{A}+{\overline{B}}_{2}K\right)X_2\\ X_2{\left(\overline{A}+{\overline{B}}_{2}K\right)}^{\text{T}}& -r{X}_2\end{array}\right]<0$ (19)

$[\begin{array}{c}\sin \theta \left(\left(\overline{A}+{\overline{B}}_{2}K\right)X_2+X_2{\left(\overline{A}+{\overline{B}}_{2}K\right)}^{\text{T}}\right)\\ \cos \theta \left(X_2{\left(\overline{A}+{\overline{B}}_{2}K\right)}^{\text{T}}-\left(\overline{A}+{\overline{B}}_{2}K\right)X_2\right)\end{array}\begin{array}{c}\cos \theta \left(\left(\overline{A}+{\overline{B}}_{2}K\right)X_2-X_2{\left(\overline{A}+{\overline{B}}_{2}K\right)}^{\text{T}}\right)\\ \sin \theta \left(\left(\overline{A}+{\overline{B}}_{2}K\right)X_2+X_2{\left(\overline{A}+{\overline{B}}_{2}K\right)}^{\text{T}}\right)\end{array}]<0$ (20)

这一组矩阵不等式是非线性的。为便于求解，参考文献[13]采用变量替换法将它们转化为关于新变量的LMI。

令$X_2=X_1$ ，定义替换变量$Y_1=K{X}_2=K{X}_1$ ，则$Y_1^{\text{T}}=X_1^{\text{T}}{K}^{\text{T}}=X_1{K}^{\text{T}}$ ，将这个变量替换公式带入到非线性矩阵不等式(18)~(20)中，这三个不等式分别转化为

$\overline{A}{X}_1+{\overline{B}}_2{Y}_1+X_1{\overline{A}}^{\text{T}}+Y_1{}^{\text{T}}{\overline{B}}_2^{\text{T}}+2\alpha X_1<0$ (21)

$\left[\begin{array}{cc}-r{X}_1& \overline{A}{X}_1+B_2{Y}_1\\ X_1{\overline{A}}^{\text{T}}+Y_1{}^{\text{T}}{\overline{B}}_2^{\text{T}}& -r{X}_1\end{array}\right]<0$ (22)

$[\begin{array}{c}\sin \theta \left(\overline{A}{X}_1+{\overline{B}}_2{Y}_1+X_1{\overline{A}}^{\text{T}}+Y_1^{\text{T}}{\overline{B}}_2^{\text{T}}\right)\\ \cos \theta \left(X_1{\overline{A}}^{\text{T}}+Y_1{}^{\text{T}}{\overline{B}}_2^{\text{T}}-\overline{A}{X}_1-{\overline{B}}_2{Y}_1\right)\end{array}\begin{array}{c}\cos \theta \left(\overline{A}{X}_1+{\overline{B}}_2{Y}_1-X_1{\overline{A}}^{\text{T}}-Y_1^{\text{T}}{\overline{B}}_2^{\text{T}}\right)\\ \sin \theta \left(\overline{A}{X}_1+{\overline{B}}_2{Y}_1+X_1{\overline{A}}^{\text{T}}+Y_1^{\text{T}}{\overline{B}}_2^{\text{T}}\right)\end{array}]<0$ (23)

至此，非线性矩阵不等式(18)~(20)均已转化为关于矩阵变量$X_1$ 和$Y_1$ 的线性矩阵不等式(21)~(23)。

对于推力跟踪控制问题P2的求解，给出如下定理。

定理1：若关于对称正定矩阵$X$ 和适当维数矩阵$Y$ 的LMI (24)~(27)同时有解，那么XTER某工作点处推力跟踪控制问题P2可解且控制器为$u=Y{X}^{-1}\overline{x}$ 。

$\left[\begin{array}{ccc}\overline{A}X+{\overline{B}}_{2}Y+{\left(\overline{A}X+{\overline{B}}_{2}Y\right)}^{\text{T}}& {\overline{B}}_1& {\left({\overline{C}}_{1}X+{\overline{D}}_{12}Y\right)}^{\text{T}}\\ *& -\gamma I& {\overline{D}}_{11}^{\text{T}}\\ *& *& -\gamma I\end{array}\right]<0$ (24)

$\overline{A}X+{\overline{B}}_{2}Y+X{\overline{A}}^{\text{T}}+Y^{\text{T}}{\overline{B}}_2^{\text{T}}+2\alpha X<0$ (25)

$\left[\begin{array}{cc}-rX& \overline{A}X+B_{2}Y\\ X{\overline{A}}^{\text{T}}+Y^{\text{T}}{\overline{B}}_2^{\text{T}}& -rX\end{array}\right]<0$ (26)

$[\begin{array}{c}\sin \theta \left(\overline{A}X+{\overline{B}}_{2}Y+X{\overline{A}}^{\text{T}}+Y^{\text{T}}{\overline{B}}_2^{\text{T}}\right)\\ \cos \theta \left(X{\overline{A}}^{\text{T}}+Y^{\text{T}}{\overline{B}}_2^{\text{T}}-\overline{A}X-{\overline{B}}_{2}Y\right)\end{array}\begin{array}{c}\cos \theta \left(\overline{A}X+{\overline{B}}_{2}Y-X{\overline{A}}^{\text{T}}-Y^{\text{T}}{\overline{B}}_2^{\text{T}}\right)\\ \sin \theta \left(\overline{A}X+{\overline{B}}_{2}Y+X{\overline{A}}^{\text{T}}+Y^{\text{T}}{\overline{B}}_2^{\text{T}}\right)\end{array}]<0$ (27)

注1：该定理求解得到的控制器属于γ-次优控制器，为了得到最优控制器，还需要把γ极小化。

依据定理1求得原问题的γ-次优控制器后，为了得到最优控制器，把γ极小化。即求解如下优化问题

$\begin{array}{l}\min \gamma \\ s.t.\left(24\right)~\left(27\right)\end{array}$ (28)

如果该优化问题有解，则推力跟踪控制问题P2的最优控制器存在且控制器为$u^{*}=Y^{*}{X}^{*-1}\overline{x}$ 。

通过上述方法，可将XTER某个工作点推力跟踪控制问题P2转化为相应的LMI的求解。通过求解LMI，对该工作点设计兼顾鲁棒性和闭环系统的暂态特性的状态反馈控制器。对每个工作点都设计了控制器后，对控制器参数采用对马赫数插值的方法实现全速域控制器的调度。

4. 仿真验证

为验证基于LPV模型所设计的控制器的控制效果，将得到的控制器加入到非线性仿真模型上进行推力跟踪验证，并将控制效果和原LADRC控制算法进行对比。

4.1. 单个工作点处推力跟踪验证

对于单个工作点处的推力跟踪控制，以2.1 Ma为例进行说明。在2.1 Ma时，XTER发动机的飞行高度为11.592 km，推力需求为3132.1 N，涡喷发动机平衡点燃油流量为0.0693 kg/s，亚燃发动机的平衡点燃油流量为0.038 kg/s。2.1 Ma处XTER发动机考虑外扰的线性控制模型为

$\begin{array}{c}\left[\begin{array}{c}\Delta {\dot{n}}_1\\ \Delta {\dot{T}}_1^{*}\\ \Delta {\dot{T}}_2^{*}\\ \Delta {\dot{F}}_2\\ \Delta \dot{W}{f}_1\\ \Delta \dot{W}{f}_2\\ \dot{\xi }\end{array}\right]=\left[\begin{array}{ccccccc}\text{4}\text{.7312}& -\text{1}\text{.4586}& 0& 0& -1.4955& 0& 0\\ \text{19}\text{.5260}& -\text{5}\text{.9502}& 0& 0& -\text{1}\text{.4228}& 0& 0\\ 0& 0& -\text{5}\text{.8080}& -\text{5}\text{.9088}& 0& \text{7}\text{.1000}& 0\\ 0& 0& \text{9}\text{.7467}& -\text{5}\text{.4440}& 0& \text{7}\text{.0080}& 0\\ -\text{7}\text{.2058}& \text{1}\text{.8310}& 0& -\text{1}\text{.8080}& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\Delta n_1\\ \Delta T_1^{*}\\ \Delta T_2^{*}\\ \Delta F_2\\ \Delta W{f}_1\\ \Delta W{f}_2\\ \xi \end{array}\right]\\ +\left[\begin{array}{cc}0& 0\\ 0& 0\\ 0& 0\\ 0& 0\\ -1& 1\end{array}\right]\left[\begin{array}{c}d\\ r\end{array}\right]+\left[\begin{array}{cc}0& 0\\ 0& 0\\ 0& 0\\ 0& 0\\ 10& 0\\ 0& 10\\ 0& 0\end{array}\right]\left[\begin{array}{c}\Delta W{f}_{1r}\\ \Delta W{f}_{2r}\end{array}\right]\end{array}$ (29)

$\begin{array}{l}\Delta y_{t1}=\Delta y_1+\Delta y_2=\left[\begin{array}{ccccccc}\text{7}\text{.2058}& -\text{1}\text{.8310}& 0& 1& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\Delta n_1\\ \Delta T_1^{*}\\ \Delta T_2^{*}\\ \Delta F_2\\ \Delta W{f}_1\\ \Delta W{f}_2\\ \xi \end{array}\right]\\ \end{array}$ (30)

$\left[\begin{array}{c}\xi \\ y_d\end{array}\right]=\left[\begin{array}{ccccccc}0& 0& 0& 0& 0& 0& 1\\ \text{7}\text{.2058}& -\text{1}\text{.8310}& 0& 1& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\Delta n_1\\ \Delta T_1^{*}\\ \Delta T_2^{*}\\ \Delta F_2\\ \Delta W{f}_1\\ \Delta W{f}_2\\ \xi \end{array}\right]+\left[\begin{array}{cc}0& 0\\ 1& 0\end{array}\right]\left[\begin{array}{c}d\\ r\end{array}\right]$ (31)

令$\gamma =2$ ，$r=2$ ，$\alpha =0.6$ ，$\theta =\frac{\pi }{3}$ ，利用MATLAB中的LMI工具箱编写程序即可求得状态反馈矩

阵$K=\left[\begin{array}{ccccc}\text{1}\text{.6355}& -\text{0}\text{.4205}& \text{43}\text{.1429}& -\text{27}\text{.5760}& -\text{0}\text{.0003}\\ -\text{0}\text{.0005}& \text{0}\text{.0000}& \text{9}\text{.6016}& -\text{6}\text{.1336}& \text{0}\text{.0000}\end{array}\right]$ 。此时的闭环系统伯德图和闭环阶跃响应如

图6(a)和图6(b)所示。

将在线性模型上设计好的控制器加入到XTER发动机的非线性模型上进行验证。图7为XTER发动机非线性仿真模型的仿真示意图。模型中加入幅度为0.3的外部扰动，如图8所示。

相比LADRC控制算法而言，基于极点配置的LPV鲁棒$H_{\infty }$ 控制算法有更好的跟踪速度和跟踪精度。在2.1 Ma工作点处，基于极点配置的LPV鲁棒$H_{\infty }$ 控制器能使发动机推力输出在0.05 s以内跟踪上目标推力，而LADRC控制器需要大约1 s的时间使推力输出跟踪上目标推力(图9)；基于极点配置的LPV鲁棒$H_{\infty }$ 控制器的推力跟踪相对误差始终保持在0.1%以内，而LADRC控制器的相对误差相对较大(图10)。进一步提取对比各通道燃油量曲线(图11)，发动机的燃油输入也在正常工作的范围内。

因此采用基于极点配置的LPV鲁棒$H_{\infty }$ 控制算法相较于原LADRC控制算法，可以改善控制系统的动态特性并且减小稳态误差。

(a) (b)

Figure 6. (a) Bode diagram of 2.1 Ma closed-loop system; (b) Step response of 2.1 Ma closed-loop system

图6. (a) 2.1 Ma闭环系统伯德图；(b) 2.1 Ma闭环系统阶跃响应

Figure 7. Schematic diagram of nonlinear simulation of XTER engine

图7. XTER发动机非线性仿真示意图

Figure 8. External disturbance

图8. 外部扰动

Figure 9. 2.1 Ma thrust tracking relative error

图9. 2.1 Ma推力跟踪相对误差曲线

Figure 10. Relative error of thrust output of XTER engine

图10. XTER发动机推力输出相对误差

Figure 11. Fuel flow of each channel

图11. 各通道燃油流量曲线

4.2. 鲁棒性验证

为验证所设计控制器的鲁棒性，对XTER发动机的环境量马赫数和高度分别进行±5%和±10%的拉偏，图12和图13为拉偏后的推力输出和原环境条件下推力输出的比值。从图12和图13可以看出，系统具有一定的鲁棒性，对于合理范围内的环境量变化，能够维持原来的推力输出基本不变。相对于马赫数而言，飞行高度的变化对于XTER发动机的推力输出影响更大。

Figure 12. The relative error of thrust output of XTER engine after Mach number deviation

图12. 马赫数拉偏后XTER发动机推力输出相对误差

Figure 13. The relative error of thrust output of XTER engine after flight height deviation

图13. 飞行高度拉偏后XTER发动机推力输出相对误差

4.3. XTER发动机全速域推力跟踪验证

通过仿真，验证XTER发动机在全速域内对目标航迹推力的跟踪情况。XTER发动机的工作马赫数和飞行高度随时间的变化关系如图14(a)和图14(b)所示。发动机每个工作点处的燃油输入由设计点的平衡燃油通过马赫数插值得到，涡喷和亚燃分流板开度的输入参考文献[5]。将设计的控制器叠加到非线性仿真模型上，非线性仿真模型的推力跟踪曲线以及推力跟踪相对误差如图15和图16所示。可以看出，在飞行条件不断变化的情况下，控制器能使发动机在全速域内以较快速度跟踪上目标推力，且推力跟踪相对误差始终保持在0.5%以内，满足XTER发动机推力跟踪需求。

(a) (b)

Figure 14. (a) Variation curve of Mach number; (b) Variation curve of flying height

图14. (a) 马赫数变化曲线；(b) 飞行高度变化曲线

Figure 15. XTER engine thrust tracking curve

图15. XTER发动机推力跟踪曲线

Figure 16. Relative error of XTER engine thrust tracking

图16. XTER发动机推力跟踪相对误差

5. 结论

本文针对XTER发动机的推力跟踪控制中存在的不足，设计出一种既考虑鲁棒性又兼顾暂态特性的多目标控制器。通过求解LMI的方式对控制器进行求解。最后将得到的控制器加入到非线性仿真模型中进行验证。仿真结果表明：在存在扰动的情况下，基于状态反馈的考虑区域极点配置的LPV鲁棒$H_{\infty }$ 控制器能够使XTER发动机在0.05 s以内跟踪上目标推力，且对目标推力跟踪的相对误差不超过0.1%。在外部环境量变化的情况下能保持较小的推力变化。相较于原LADRC控制器而言，具有很强的鲁棒性，同时具有良好的稳态和动态性能。

致 谢

感谢何增明和张俊龙师兄对XTER发动机建模的指导，感谢孙洪飞教授对本文的修改完善。

基金项目

本文受航空科学基金资助(20220058068001)。