CN103306822A - 一种基于喘振裕度估计模型的航空涡扇发动机控制方法 - Google Patents

Abstract

本发明公开了一种基于喘振裕度估计模型的航空涡扇发动机控制方法。本发明提出的喘振裕度估计模型的建模分为常规飞行时的无畸变模型与超机动飞行时的损失量模型两部分：无畸变模型是基于喘振裕度特征选择算法筛选最优模型输入，以非线性拟合方法建模实现；损失量模型则基于在线攻角预测模型实时评估航空发动机进口畸变度，进而计算获得畸变时喘振裕度损失量。而后利用上述估计模型对发动机的稳定性进行实时预测，在不改变发动机常规控制回路的基础上，对涡轮落压比控制指令进行喘振损失补偿，实现高稳定性控制。本发明可保证航空发动机稳定、高效地工作，实现航空发动机喘振裕度控制在11%~13%之间。

Description

一种基于喘振裕度估计模型的航空涡扇发动机控制方法

技术领域

本发明涉及一种航空涡扇发动机控制方法，尤其涉及一种基于喘振裕度估计模型的航空涡扇发动机控制方法，属于航空航天控制技术领域。

背景技术

超机动能力是现代战机的必备技术之一，但在过失速状态，航空发动机需工作在大攻角和侧滑角的条件下，此时进口畸变严重，致使航空发动机性能降低、稳定工作裕度减少，严重时甚至会产生喘振等不稳定工作现象。针对这一难题，美国NASA于1993年开展了航空发动机高稳定性控制(HISTEC,High Stability Engine Control)研究，主要利用压缩部件的压比损失来近似表征喘振裕度损失，而后通过设计先进控制规律及时修正风扇及压气机的瞬态工作线，以保证航空发动机足够的稳定裕度。近年来，主动失速/喘振控制作为一种非常规的控制方法也得到了一定的研究，该技术主要基于失速先兆与相关性测量进行航空发动机主动稳定性控制，如建立压力相关度测量值与喘振裕度之间的特定数学模型，在发现进喘趋势时通过快速的切油脉冲减少燃油，使得燃烧室温度及压力下降，清除压缩部件后的流动阻碍。这类控制方法在本质上属于限制保护控制的范畴，在与主控制回路协调工作以及充分发挥航空发动机性能潜力发面还存在一定的不足。

国内对HISTEC的研究较少，王立峰等在2000年“超机动飞行推进系统稳定性控制研究”中以航空发动机进口畸变时可通过风扇喘振压比损失进行补偿，使喘振裕度回到未畸变状态为思路，设计了航空发动机稳定性控制器，该方案虽能保证航空发动机的稳定性，但在机动过程中会有喘振裕度过大现象，航空发动机潜力不能得到充分的发挥。在超机动飞行条件下若能实时预测航空发动机的喘振裕度，通过HISTEC来限制喘振裕度的范围，控制航空发动机工作线与喘振边界线之间始终保持一定的距离，航空发动机即可始终处于高效、稳定的工作状态。袁春飞与陈霆昊等分别在2004年“飞行/推进系统综合优化控制模式及其关键技术”及2010年“基于攻角预测模型的航空发动机高稳定性控制”中提出了一种直接喘振裕度控制的方法，即通过实时估算喘振裕度并将其加入控制回路中作为被控量，来实现对喘振裕度的直接控制，很明显这种控制形态更为高效，但其先决条件是必须对喘振裕度进行快速准确的估计。建立航空发动机喘振裕度估计模型一直是防喘/消喘控制中的难题，陈霆昊等将喘振裕度分为畸变时的喘振裕度损失量和常规飞行时的未畸变值两部分来分别计算，但其采用等换算空气流量的思路求取喘振裕度(实际应取等换算转速线)，通过压缩部件压比来直接计算，即将喘振裕度看作可测量，这种方法在低速时存在一定的局限性。再者，传统的非线性拟合等方法均凭经验采用飞行高度、马赫数、主燃油量、喷口面积等作为模型输入，但这些参数与喘振裕度的相关性并非紧密，故很难实现喘振裕度的准确辨识。

发明内容

本发明所要解决的技术问题在于克服现有技术不足，提供一种基于喘振裕度估计模型的航空涡扇发动机控制方法，能够保证航空发动机稳定、高效地工作。

本发明的基于喘振裕度估计模型的航空涡扇发动机控制方法，发动机控制系统对下一时刻的发动机喘振裕度进行估计，具体按照以下方法：

步骤1、利用攻角预测模型对下一时刻的航空器攻角进行在线预测；

步骤2、利用预测的下一时刻攻角确定下一时刻发动机的压缩部件的喘振压比损失ΔP_rs；

步骤3、根据下式确定下一时刻发动机喘振裕度损失量：

ΔS_m(k+1)=ΔP_rs(1+S_m(k+1))，

其中，k、k+1分别表示当前时刻和当前时刻的下一时刻；ΔS_m(k+1)表示k+1时刻的发动机喘振裕度损失量；ΔP_rs表示k+1时刻的发动机的压缩部件的喘振压比损失；S_m(k+1)表示以当前时刻的发动机可测参数作为输入，利用无畸变喘振裕度模型估计出的k+1时刻的无畸变喘振裕度，所述无畸变喘振裕度模型为预先通过离线训练得到的神经网络；

步骤4、按照下式得到下一时刻的发动机喘振裕度估计值

$S_m^{\′}(k+1)=S_m(k+1)-\mathrm{\Δ}{S}_m(k+1),$

其中，S_m(k+1)表示利用无畸变喘振裕度模型估计出的k+1时刻的无畸变喘振裕度，ΔS_m(k+1)表示步骤3所得到的k+1时刻的发动机喘振裕度损失量；

当所得到的

超出预设的工作裕度范围时，发动机控制系统对发动机涡轮落压比的控制指令π_t,r进行喘振压比损失补偿：当低于预设的工作裕度范围时，补偿后的发动机涡轮落压比的控制指令为π_t,r+Δπ′；当

高于预设的工作裕度范围时，补偿后的发动机涡轮落压比的控制指令为π_t,r+(-Δπ′)；Δπ′为所述喘振压比损失补偿的补偿量，根据下式得到：

Δπ'=(π_t/π_f)·ΔP_rs，

式中，ΔP_rs表示k+1时刻的发动机的压缩部件的喘振压比损失；π_t、π_f分别表示发动机的涡轮落压比、风扇落压比。

由于可作为无畸变喘振裕度模型的输入的发动机可测参数数量众多，为了提升模型的精度并降低计算量，有必要对可测参数进行筛选，作为本发明的进一步改进方案：作为所述无畸变喘振裕度模型的输入的发动机可测参数，是通过最小二乘支持向量回归机算法从所述航空涡扇发动机的所有N个可测参数中筛选出的与喘振裕度相关性最强的前M个可测参数，M＜N。

相比现有技术，本发明具有以下有益效果：

(1)本发明所采用的喘振裕度估计模型将喘振裕度分为无畸变喘振裕度与畸变时喘振裕度损失量，实现喘振裕度实时预测，并实时控制，可以准确估计出航空发动机进气畸变状况下的喘振裕度，提升航空发动机工作效率与安全性水平。

(2)本发明根据实时预测出的喘振裕度，在喘振裕度超出预设正常范围时，根据压缩部件喘振压比损失比例，计算出航空发动机涡轮压比损失量，进而在航空发动机原控制回路中直接对其控制指令进行补偿，如此便不必改变航空发动机原控制器，使该控制方案更具优势与工程应用性。

附图说明

图1是采用本发明方法的航空发动机控制系统结构图；

图2（a）、图2（b）分别为风扇和压气机喘振裕度输入特征选择过程示意图；

图3（a）-图3（d）是风扇喘振裕度预测模型精度测试效果对比，其中图3（a）、图3（b）分别为经特征选择后风扇喘振裕度预测模型训练的测试效果图和测试误差图；图3（c）、图3（d）分别为采用传统方法训练模型的测试效果图和测试误差图；

图4（a）-图4（d）是压气机喘振裕度预测模型精度测试效果对比，其中图4（a）、图4（b）分别为经特征选择后压气机喘振裕度预测模型训练的测试效果图和测试误差图；图4（c）、图4（d）分别为采用传统方法训练模型的测试效果图和测试误差图；

图5（a）、图5(b)是低空（高度H=5km，马赫数Ma=0.5）状态下的攻角预测模型精度测试结果，其中图5（a）为攻角实际值与预测值的对比曲线，图5(b)为攻角测试误差图；

图6（a）、图6(b)是高空（高度H=12km，马赫数Ma=0.8）状态下的攻角预测模型精度测试结果，其中图6（a）为攻角实际值与预测值的对比曲线，图6(b)为攻角测试误差图；

图7（a）、图7(b)是本发明的控制系统在高度H=6km，马赫数Ma=0.63状态下的仿真结果，图7(a)为风扇喘振裕度S_m,f的响应曲线，图中，“-no”为常规控制响应，“-histec”为高稳定性控制响应；图7(b)中从上到下依次为飞行攻角α和升降舵偏转角δ的变化曲线、进气畸变指标DC₆₀和喘振裕度损失的变化曲线、压气机喘振裕度S_m,c、涡轮落压比π_t、风扇相对换算转速n_f、压气机相对换算转速n_c、喷管面积A₈、主燃油量W_fb、加力燃油量W_fa的响应曲线；

图8（a）、图8(b)是本发明的控制系统在高度H=4km,，马赫数Ma=0.54状态下的仿真结果，图8(a)为风扇喘振裕度S_m,f的响应曲线，图中，“-no”为常规控制响应，“-histec”为高稳定性控制响应；图8(b)中从上到下依次为飞行攻角α和升降舵偏转角δ的变化曲线、进气畸变指标DC₆₀和喘振裕度损失的变化曲线、压气机喘振裕度S_m,c、涡轮落压比π_t、风扇相对换算转速n_f、压气机相对换算转速n_c、喷管面积A₈、主燃油量W_fb、加力燃油量W_fa的响应曲线。

具体实施方式

下面结合附图对本发明的技术方案进行详细说明：

本发明的思路是通过建立可准确估计航空发动机进气畸变工况下的喘振裕度的估计模型来对航空发动机的喘振裕度进行在线估计，并根据实时预测出的喘振裕度，在喘振裕度超出预设正常范围时，根据压缩部件喘振压比损失比例，计算出航空发动机涡轮压比损失量，进而在航空发动机原控制回路中直接对其控制指令进行补偿。从而在不改变常规控制回路基础上，即可实现中低空、大攻角低速飞行状态下航空发动机高稳定性控制。

本具体实施方式以模型双轴混合排气带加力的航空发动机高稳定性控制方案为例，该控制系统如图1所示，包括航空发动机喘振裕度估计模型与飞机/发动机综合仿真模型，分别用于喘振裕度实时预测与航空发动机高稳定性控制、仿真与验证。该控制系统的建立包括以下步骤：

步骤A、通过LSSVR（最小二乘支持向量回归机）算法选择航空发动机的可测变量作为未畸变喘振裕度模型的输入量；

选取与喘振裕度相关性较强的可测参数作为喘振裕度模型输入，不仅能够提升模型的精度，而且较易实现喘振裕度的准确辨识。本发明采用LSSVR算法来进行可测参数的筛选，具体如下：

给定一组训练数据集

其中x_i∈Rⁿ,y_i∈R，N为训练样本的规模，则可得最小二乘支持向量回归机模型

式中：e_i为模型输出误差；w为高维特征空间的权矩阵；y∈R⁺为惩罚因子；

为输入空间Rⁿ到高维特征空间的非线性映射；b为偏置量。构建此模型的Lagrange函数

式中α_i为Lagrange乘子，其KKT（Karush-Kuhn-Tucker）条件为

在消去e_i和w后，可得线性系统

$\left[\begin{array}{cc}0& 1^T\\ 1& K\end{array}\right]\left[\begin{array}{c}b\\ \mathrm{\α}\end{array}\right]=\left[\begin{array}{c}0\\ Y\end{array}\right]---\left(4\right)$

式中：Y=[y₁,y₂,…,y_N]^T,α=[α₁,α₂,…,α_N]^T,1[1₁,1₂,…,1_N]^T,K中元素 $K_{i,j}=k(x_i,x_j)+{\mathrm{\δ}}_{i,j}/\mathrm{\γ},{\mathrm{\δ}}_{i,j}=\left\{\begin{array}{cc}1& i=j\\ 0& i\≠j\end{array}\right.,k(\·,\·)$ 为核函数。解式(4)即可求得α,b，进而给定一个新的输入变量x，可得：

$f\left(x\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_{i}k(x,x_i)+b---\left(5\right)$

对无畸变喘振裕度估计模型，式(5)中的f(x)即为其估计值，x为模型输入，因而x的维数必须加以确定。首先，结合式(2)和式(3)，得到如下不带约束的Wolfe对偶优化问题：

$\min \{\stackrel{\_}{L}(b,{\mathrm{\α}}_i)=[\frac{1}{2}\underset{i,j=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_i{\mathrm{\α}}_{j}k\left(x_i{,x}_j\right)+\frac{1}{2\mathrm{\γ}}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_i^2-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_i{y}_i+b\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_i\left]\right\}---\left(6\right)$

进而，对式(6)进行如下等价变形：

$\stackrel{\_}{L}(b,\mathrm{\α})=[\frac{1}{2}\left[\begin{array}{cc}b& {\mathrm{\α}}^T\end{array}\right]\left[\begin{array}{cc}0& 1^T\\ 1& K\end{array}\right]\left[\begin{array}{c}b\\ \mathrm{\α}\end{array}\right]-\left[\begin{array}{cc}b& {\mathrm{\α}}^T\end{array}\right]\left[\begin{array}{c}0\\ Y\end{array}\right]\left]\right\}---\left(7\right)$

再令 $\frac{d\stackrel{\‾}{L}}{d\stackrel{\‾}{\mathrm{\α}}}=0,$ 其中 $\stackrel{\‾}{\mathrm{\α}}={[b,{\mathrm{\α}}^T]}^T,$ 则有：

$\left[\begin{array}{c}b\\ \mathrm{\α}\end{array}\right]={\left[\begin{array}{cc}0& 1^T\\ 1& K\end{array}\right]}^{-1}\left[\begin{array}{c}0\\ Y\end{array}\right]---\left(8\right)$

将式(8)入式(7)，则可得到最优值

${\stackrel{\‾}{L}}^{*}=-\frac{1}{2}\left[\begin{array}{cc}0& Y^T\end{array}\right]{\left[\begin{array}{cc}0& 1^T\\ 1& K\end{array}\right]}^{-1}\left[\begin{array}{c}0\\ Y\end{array}\right]=-\frac{1}{2}\left[\begin{array}{cc}0& Y^T\end{array}\right]\left[\begin{array}{c}b\\ \mathrm{\α}\end{array}\right]=-\frac{1}{2}{Y}^T\mathrm{\α}---\left(9\right)$

在计算过程中，若

取得了最优值

则此时的喘振裕度估计精度最高，若将输入向量x中的第i个参数移除，则最优值变为

${\stackrel{\‾}{L}}^{*}\left(i\right)=-\frac{1}{2}{Y}^T\mathrm{\α}\left(i\right)$

式中α(i)为移除第i个变量后式(11)的解。由此，可设计如下准则，即根据对喘振裕度的影响程度来对x中各个参数进行排序：

Δ(i)=|Y^T(α-α(i))|       (10)

对于x中第i和第j两个不同参数，如果Δ(i)的值小于Δ(j)，则说明第i个参数对最优值

的贡献程度不及第j个参数。通过Δ(i)序列，即可每次剔除一个对喘振裕度影响最小的参数。最终，在不同飞行高度和马赫数下通过对油门杆角度充分激励，通过采集3816组动态数据用于喘振裕度特征选择，并将与风扇及压气机喘振裕度相关的21个航空发动机可测参数作为筛选范围。该筛选过程如图2（a）、图2（b）所示，图中“特征选择次序”对应图中一个黑点，该黑点表示其对应的纵轴参数在该步筛选中被舍弃，越早舍弃说明该参数的价值越小，其中H为飞行高度，M_a为马赫数，n_f和n_c分别为风扇和压气机相对换算转速，W_fb、W_fa、A₈分别为主燃油量、加力燃油量和尾喷口喉道面积，T、P、τ、π分别为总温、总压、温比和压比，下标2、22、23、25、3分别表示风扇进口截面、风扇出口截面、风扇引气截面、压气机进口截面、压气机出口截面。从而依次排出对喘振裕度贡献最大的参数队列，为尽量减少传感器数目以及简化训练模型，在上述两组参数中仅选用前5个相关性较强的参数分别作为模型的输入，即风扇喘振裕度模型的输入参数包括T₂₂,P₂,n_f,P₃,W_fb，压气机喘振裕度模型则包括n_c,T₃,π_c,P₃,T₂₂。

步骤B、采用对非线性映射具有任意逼近能力的BP神经网络设计无畸变喘振裕度模型；

本发明采用BP三层神经网络进行模型设计，其数学表达为：

Y=f(X)      （11）

对于风扇喘振裕度，有

$\left\{\begin{array}{c}X=[T_{22}(k-1),T_{22}(k-2),\·\·\·,T_{22}(k-N_1),\\ P_2(k-1),P_2(k-3),\·\·\·,P_2(k-N_2),\\ n_f(k-1),n_f(k-2),\·\·\·,n_f(k-N_3),\\ P_3(k-1),P_3(k-2),\·\·\·,P_3(k-N_4),\\ W_{\mathrm{fb}}(k-1),W_{\mathrm{fb}}(k-2),\·\·\·,W_{\mathrm{fb}}(k-N_5)]\\ Y=\left[S_{m,f}\left(k\right)\right]\end{array}\right.---\left(12\right)$

由于航空发动机模型一般可简化为一个二阶对象，为了保证S_m模型具有合适的动静态精度且输入参数尽可能少，将N₁,N₂,N₃,N₄,N₅均设置为3。根据该模型，即可实时预测下一时刻的喘振裕度值，即：

Y_k+1=f(X_k+1)

此外，由于神经网络类似于一个非线性插值器，在内插时精度较高，而外插时精度较低，因而所选的训练数据应尽量包含各变量可能出现的最大最小值。在模型训练过程中，首先在3～7km,0.3～0.7M_a的包线范围内通过对油门杆角度进行充分激励，采集18227组动态数据用作模型训练样本，另外随机采集3717组数据作为测试样本(仿真步长20ms)；而后以调试法来确定映射网络的神经元个数，最终选择S_m,f模型的神经元个数为14，该模型的测试结果如图3（a）-图3（d）所示。同理，也可对压气机喘振裕度模型进行设计，最终确定S_m,c模型的神经元个数为11，其测试结果如图4（a）-图4（d）所示。两组图中，(a)和(b)分别为经特征选择后模型训练的测试效果图和测试误差图，(c)和(d)则为采用传统方法训练模型的测试效果图和测试误差图(即模型输入参数为传统的H,M_a,W_fb,W_fa,A₈)，其中测试误差均为相对误差，表达如下：

$\mathrm{error}=\left|\frac{{\tilde{s}}_m-s_m}{s_m}\right|\×100\%$

式中：为喘振裕度估计值，S_m为航空发动机部件级模型实际输出，即目标值。

步骤C、基于OSP-LSSVR（在线滚动最小二乘支持向量回归机）算法设计在线攻角预测模型；

攻角预测模型，采用OSP-LSSVR算法设计，具体如下：

基于LVSSR算法，迭代策略来减轻训练代价，即假定在第n步：

${R^n=\left[\begin{array}{cc}0& 1_{\left|P\right|}^T\\ 1_{\left|P\right|}& {\stackrel{\_}{K}}_{\mathrm{PP}}\end{array}\right]}^{-1}---\left(13\right)$

式中：P是第n步时支持向量的索引集合，|·|代表集合的势。则可进一步求得：

$R^{n+1}=\left[\begin{array}{cc}{R}^n& 0\\ 0^T& 0\end{array}\right]+\mathrm{\λ}\left[\begin{array}{c}\mathrm{\β}\\ -1\end{array}\right]\left[\begin{array}{cc}{\mathrm{\β}}^T& -1\end{array}\right]---\left(14\right)$

式中： $\mathrm{\β}=R^n\left[\begin{array}{c}1\\ {\stackrel{\‾}{k}}_s\end{array}\right],\mathrm{\λ}={({\stackrel{\‾}{k}}_{\mathrm{ss}}-{\left[\begin{array}{cc}1& {\stackrel{\‾}{k}}_s^T\end{array}\right]}^T{\mathrm{\β}}^T)}^{-1},$ $s\∈P,{\stackrel{\_}{k}}_{\mathrm{ss}}=k(x_s,x_s)+1/\mathrm{\γ}.$ 当第n步的 $\left[\begin{array}{c}{b}^n\\ a_P^n\end{array}\right]=R^n\left[\begin{array}{c}0\\ Y_P\end{array}\right]$ 已知，则第n+1步：

$\left[\begin{array}{c}{b}^{n+1}\\ a_p^{n+1}\\ a_s^{n+1}\end{array}\right]=R^{n+1}\left[\begin{array}{c}0\\ Y_P\\ y_s\end{array}\right]=\left[\begin{array}{cc}{R}^n& \left[\begin{array}{c}0\\ Y_P\end{array}\right]\\ & 0\end{array}\right]+\mathrm{\λ}({\mathrm{\β}}^T\left[\begin{array}{c}0\\ Y_P\end{array}\right]-y_s)\left[\begin{array}{c}\mathrm{\β}\\ -1\end{array}\right]=\left[\begin{array}{c}{b}^n\\ a_P^n\\ 0\end{array}\right]+\mathrm{\λ}({\mathrm{\β}}^T\left[\begin{array}{c}0\\ Y_P\end{array}\right]-y_s)\left[\begin{array}{c}\mathrm{\β}\\ -1\end{array}\right]---\left(15\right)$

如此，R,α,b就可以有效地迭代计算出来，进而可实现攻角预测模型的构建。但此时支持向量集的容量没有上限，还不能满足攻角预测模型的实时性需求，下面结合滚动窗法进行改进。设滚动窗口样本数据的最大容量为个，若第m+1时刻加入输入向量x_m+1后，窗口数据不超过

个，则将x_m+1直接加入窗口中作为支持向量，否则根据下述推导先将最远时刻的输入向量

删掉后再加入。删除数据时，首先根据Sherman-Morrison定理将m时刻的

分解如下：

$K_m^{-1}={\left[\begin{array}{cc}{g}_m& G_m^T\\ G_m& {\stackrel{\‾}{K}}_m\end{array}\right]}^{-1}=\left[\begin{array}{cc}0& 0\\ 0& {\stackrel{\‾}{K}}_m^{-1}\end{array}\right]+{\stackrel{\‾}{z}}_m{\stackrel{\‾}{z}}_m^T{\stackrel{\‾}{r}}_m=\left[\begin{array}{cc}0& 0\\ 0& {\stackrel{\‾}{K}}_m^{-1}\end{array}\right]+\left[\begin{array}{cc}{U}^{-1}& V^T\\ V& U^{-1}{\mathrm{VV}}^T\end{array}\right]---\left(16\right)$

式中：

表示m时刻删除最远样本后的K阵，g_m=k(x₁,x₁)+1/γ, ${\stackrel{\‾}{r}}_m={(g_m-G_m^T{\stackrel{\‾}{K}}_m^{-1}{G}_m)}^{-1},$ ${G_m=[k(x_1,x_2),k(x_1,x_3),\·\·\·,k(x_1,x_{\stackrel{\_}{l}})]}^T,$ ${{\stackrel{\‾}{z}}_m=[-1,G_m^T{\stackrel{\‾}{K}}_m^{-1}]}^T,$ $U={{\stackrel{\_}{r}}_m}^{-1},V=-{{\stackrel{\_}{r}}_m}^{-1}{{\stackrel{\_}{K}}_M}^{-1}{G}_m.$

通过式(16)可提取出

而后将式(4)改写为：

$\left[\begin{array}{c}{b}_m\\ a_m\end{array}\right]={\left[\begin{array}{cc}0& {\stackrel{\→}{1}}^T\\ \stackrel{\→}{1}& K_m\end{array}\right]}^{-1}\left[\begin{array}{c}0\\ Y_m\end{array}\right]---\left(17\right)$

不难求得：

$b_m=\frac{1^T{K}_m^{-1}{Y}_m}{1^T{K}_m^{-1}1},{\mathrm{\α}}_m=K_m^{-1}[Y_m-\frac{{11}^T{K}_m^{-1}{Y}_m}{1^T{K}_m^{-1}1}]$

进而，即可利用式(14)、（15)计算出m+1时刻的R,α,b，以实现对攻角状态的实时在线滚动预测，而始终保持有效支持向量的数量恒定且最接近当前时刻。

上述在线滚动LSSVR算法仍将所有的输入数据

都作为支持向量，由于滚动窗口容量有限，一些冗余的、重复的输入信息必然使得攻角预测精度大打折扣。所以，加入阈值判断准则：事先设定一个较小正数ε，当预测精度小于ε时，则认为该输入信息为冗余信息并舍弃，反之则作为支持向量执行式(7-10)，进而更新预测模型，如此便实现了OSP-LSSVR算法。

攻角预测模型使用最近时刻且最有效的攻角信息建模，可以很好地保留系统动态特性，而由于相邻时刻的攻角特性相近，故退一步采用k时刻的攻角预测误差代替k+1时刻的误差，即以

作为预测模型更新的判断准则，设计3阶的攻角预测模型如下：

Y′=f(X′)   (18)

式中， $\{X^{\′}=\begin{array}{c}\left[\begin{array}{ccc}{\mathrm{\α}}_{k-3}& {\mathrm{\α}}_{k-2}& {\mathrm{\α}}_{k-1}\end{array}\right]\\ Y\′=\left[\begin{array}{c}{\mathrm{\α}}_k\end{array}\right]\end{array}$

由此可推导出在k+1时刻模型输入-输出关系：

$Y_{k+1}^{\′}=f\left(X_{k+1}^{\′}\right)---\left(19\right)$

式中， $\left\{\begin{array}{c}{X}_{k+1}^{\′}=\left[\begin{array}{ccc}{\mathrm{\α}}_{k-2},& {\mathrm{\α}}_{k-1},& {\mathrm{\α}}_k\end{array}\right]\\ Y_{k+1}^{\′}=\left[\begin{array}{c}{\mathrm{\α}}_{k+1}\end{array}\right]\end{array}\right.$

如此设计，由于当前时刻的α_k可测，即可通过式(18)求取α_k后，经

进行阈值判断，决定对现有模型的支持向量集是否进行更新，而后进行式(19)的迭代计算。为了提升模型精度，在设计过程中对输入参数均作了归一化处理，因而对阈值ε也进行了改进，即采用k时刻模型输出的相对偏差量e_k作为目标值，如下：

$e_k=\left|\frac{{\widetilde{\mathrm{\&alpha;}}}_k-{\mathrm{\&alpha;}}_k}{{\mathrm{\&alpha;}}_k}\right|\&times;100\%<\mathrm{\&epsiv;}\text{'}$

在阈值判断过程中，若e_k小于ε，则认为此时的攻角预测模型可以满足精度要求，不需要更新支持向量集便可准确预测出下一时刻的攻角α_k+1，反之则认为此时预测模型精度不足，需要加入k时刻数据对支持向量集进行实时更新后，方可进行k+1时刻的攻角预测。

最终，通过在OSP-LSSVR算法中使用高斯核函数k(x_i,x_j)=exp(-||x_i-x_j||²/2v²)设计在线攻角预测模型，经调试设定模型相关参数如下：

滚动窗口容量

个，判断阈值ε′=0.005，核参数v=1.6，惩罚因子γ=2^20。

其中

的选择原则为滚动窗口容量尽量小以优化在线模型的实时性，判断阈值ε′=0.005表示模型训练精度可达到5‰，进而可使得模型测试精度较高约为1%，核参数与惩罚因子则通过交叉验证的方法确定。该攻角预测模型的设计具有很好的实时性与鲁棒性，在全包线范围均可适用。

分别以H=5km,Ma=0.5和H=12km,Ma=0.8状态为例，说明上述模型的有效性。图5（a）、图5（b）是低空（高度H=5km，马赫数Ma=0.5）状态下的攻角预测模型精度测试结果，其中图5（a）为攻角实际值与预测值的对比曲线，图5(b)为攻角测试误差图；图6（a）、图6（b）高空（高度H=12km，马赫数Ma=0.8）状态下的攻角预测模型精度测试结果，其中图6（a）为攻角实际值与预测值的对比曲线，图6(b)为攻角测试误差图。在测试过程中，通过为对升降舵角δ进行充分激励，产生攻角的动态变化，其变化范围为-20°～+70°。如图5(a)所示，为k+1时刻的攻角实际值与预测值对比曲线，包含2016组动态数据，经OSP-LSSVR算法在线运算仅需219ms，且从图中可以看出两条曲线的吻合程度较高，说明本发明所采用的攻角预测模型具有很好的响应速度及预测精度；再从图5(b)的攻角测试误差图可以看出，其相对测试误差在1%以内，比起采用离线方法训练攻角预测模型达到3°的预测误差，具有更高的精度优势。其次如图6（a）、图6（b）所示，在高空对其预测效果进行测试，仿真结果与H=5km,Ma=0.5状态下类似，该模型同样能够快速并准确地预测出k+1时刻的攻角状态，其测试精度也达到1%，实现了很好的预测效果。

步骤D、通过由进气畸变度转化的喘振压比损失进而建立喘振裕度损失量模型；

喘振裕度损失量模型的建立主要分为三个步骤：

(a)利用预测的攻角状态确定下一时刻的进气畸变量。

基于步骤C建立的攻角预测模型，本具体实施方式中采用DC₆₀为畸变指标衡量参数，则可求得：

DC₆₀(k+1)=DC₆₀(M_a,α_k+1)

式中，DC₆₀(k+1)表示k+1时刻的DC₆₀进气畸变量，M_a表示飞行马赫数，α_k+1表示所述攻角预测模型在线预测出的k+1时刻的攻角。

(b)由进气畸变量转换为喘振压比损失。

在进气畸变条件下，航空发动机对畸变的敏感性（即压缩部件的喘振压比损失）与进气畸变指标DC₆₀成一定的线性关系，经验公式如下：

ΔP_rs=π_s(K_circ·DC₆₀+Const)

式中，π_s为压缩部件的喘振边界压比，K_circ为航空发动机对周向畸变的敏感系数，常数Const可认为是径向畸变和平面波的影响因素。

(c)求取喘振裕度损失量。

在航空发动机部件级模型中，通常采用如下公式来表示航空发动机压缩部件的喘振裕度：

$S_m=\frac{{\mathrm{\&pi;}}_s/W_{\mathrm{cor},s}-{\mathrm{\&pi;}}_o/W_{\mathrm{cor},o}}{{\mathrm{\&pi;}}_o/W_{\mathrm{cor},o}}=\frac{{\mathrm{\&pi;}}_s}{{\mathrm{\&pi;}}_o}\&CenterDot;\frac{W_{\mathrm{cor},o}}{W_{\mathrm{cor},s}}-1---\left(20\right)$

式中，S_m表示压缩部件喘振裕度，π表示压缩部件增压比，W_cor表示换算空气流量，下标s,o分别代表喘振边界值和工作线上的值。在航空发动机使用过程中遇到进气道出口流场不均匀或总温突升等情况，则会使得航空发动机喘振边界下移或工作线上移，引起喘振裕度损失，以喘振边界下移为例，则有

$\mathrm{\&Delta;}{S}_m=\frac{{\mathrm{\&pi;}}_s/W_{\mathrm{cor},s}-{\mathrm{\&pi;}}_1/W_{\mathrm{cor},1}}{{\mathrm{\&pi;}}_o/W_{\mathrm{cor},o}}---\left(21\right)$

式中，ΔS_m表示喘振裕度损失量，下标1表示下移的稳定边界值。而此时，由于稳定边界线下移，航空发动机压缩部件的压比也必然产生了一定的损失，该损失量ΔP_rs的表达如下：

$\mathrm{\&Delta;}{P}_{\mathrm{rs}}=\frac{{\mathrm{\&pi;}}_s/W_{\mathrm{cor},s}-{\mathrm{\&pi;}}_1/W_{\mathrm{cor},1}}{{\mathrm{\&pi;}}_s/W_{\mathrm{cor},s}}---\left(22\right)$

结合式(20)～(22)，喘振裕度损失量可表示为：

ΔS_m=ΔP_rs(1+S_m)

由上式即可计算k+1时刻的航空发动机喘振裕度损失量：

ΔS_m(k+1)=ΔP_rs(1+S_m(k+1))    （23）式中，S_m(k+1)通过步骤B建立的无畸变喘振裕度模型来获取。

步骤E、计算航空发动机的预测喘振裕度，若判定航空发动机安全裕度超出预设范围，则切入高稳定性控制回路，改变涡轮落压比使得航空发动机裕度工作在设定范围内；

通过公式(20)及公式(23)便可直接求得处于畸变状态下的k+1时刻航空发动机喘振裕度，如下：

$S_m^{\&prime;}(k+1)=S_m(k+1)-\mathrm{\&Delta;}{S}_m(k+1)$

基于上述航空发动机喘振裕度估计模型，根据图1所示控制结构，以风扇压缩部件为例，阐述本发明的高稳定性控制方案（HISTEC）。在超机动飞行的控制过程中，认为航空发动机在无畸变状态下保留有足够的喘振裕度，但其效率并不一定较高，故设定航空发动机高效稳定的工作裕度范围S′为11%～13%。

若航空发动机发生进气畸变使得喘振裕度预测值小于11%，则判定航空发动机安全裕度过小，进而切入高稳定性控制回路，为了不改变航空发动机原闭环控制回路特性，仅在原闭环回路的基础上对涡轮落压比π_t的控制指令π_t,r进行相应的喘振压比损失补偿，其补偿量由风扇部件的压比损失按比例进行计算，如下：

Δπ'=p·ΔP_rs=(π_t/π_f)·ΔP_rs

进而可根据控制率给出新的涡轮落压比控制指令π_t,r+Δπ′，使得风扇喘振裕度达到11%以上，以保证航空发动机安全稳定工作。相反，若喘振裕度预测值大于13%，则判断航空发动机效率较低，安全裕度过剩，则给出压比指令π_t,r+(-Δπ′)，使得风扇喘振裕度低于13%，保持较高的发动机效率。再者，若喘振裕度预测值在S′范围内，则认为航空发动机可高效稳定工作，从而保持常规控制不变。

在图1所示的飞机/发动机综合仿真模型中，飞机模型为某型飞行纵向平面内的非线性动力学仿真模型，可模拟平飞直行、爬升、加速和下滑等飞行任务；该飞机由两台某型航空发动机提供推力，航空发动机模型为采用部件法建立的非线性气动热力学模型，航空发动机加力控制器为采用增广LQR方法设计的双变量控制器，控制量为A₈和W_fb，被控量为π_t和n_c(或者π_t和n_f，根据航空发动机进口温度选择相应的控制模式)，W_fa为开环给定；该综合仿真模型的仿真步长为20ms，其中使用到的攻角数据来源于NASA Dryden飞行研究中心F/A-18A飞机的飞行试验数据。

为了验证本发明控制方法的有效性，在主频2.53GHz/内存2G的仿真环境下，基于上述飞机/发动机综合仿真模型，在中低空、大攻角低速飞行条件下，对超机动飞行任务进行仿真模拟，操作过程中将油门杆拉至最大加力状态（Pla=110°），其仿真结果如图7（a）-图7（b）和图8（a）-图8（b）所示。两组图中的图（a）为风扇喘振裕度S_m,f的响应曲线，“-no”为常规控制响应，“-histec”为高稳定性控制响应，图(b)依次为飞行攻角α和升降舵偏转角δ的变化曲线、进气畸变指标DC₆₀和喘振裕度损失的变化曲线、压气机喘振裕度S_m,c、涡轮落压比π_t、风扇相对换算转速n_f、压气机相对换算转速n_c、喷管面积A₈、主燃油量W_fb、加力燃油量W_fa的响应曲线。

图7（a）-图7（b）为H=6km,Ma=0.63状态下的大攻角机动飞行仿真，在t=0s时分别采用常规控制方法和本发明中HISTEC方法进行控制，从图7(b)可以看出，在t=2s,5s,8s,12s,15s时刻，分别对升降舵偏转角δ设置阶跃(负阶跃表示向上偏转升降舵)，当δ减小时，攻角增大，进气畸变量DC₆₀增大，而DC₆₀决定了喘振裕度损失量(ΔS_m,f,ΔS_m,c)的大小，从而致使航空发动机喘振裕度降低，当机动操作完成后，进气畸变量则回到原来的水平。此外，从攻角响应曲线可以看出当前攻角α_k与预测值α_k+1非常接近，一方面说明了相邻时刻的攻角特性相似，另一方面也体现出本发明中攻角预测模型具有很好的精度。

再结合图7(a)和图7(b)，通过高稳定性控制与常规控制的响应曲线对比，可知该控制过程的被控量为π_t和n_c。当S_m,f预测值小于11%时，高稳定性控制对涡轮落压比指令进行相应的喘振损失补偿，涡轮落压比π_t上升，喷管面积A₈扩张程度明显大于常规控制，进而使得图7(a)中的S_m,f保持在11%以上；而当S_m,f预测值大于13%时，涡轮落压比指令也会得到相应的修正，喷管面积缩小，最终S_m,f保持在13%以内，致使航空发动机具有较高的效率，相比之下常规控制方法则完全不能达到这种效果。另外，由于飞行条件的改变导致W_fb和W_fa产生变化，而两种控制方法的飞行条件变化一致，故W_fb和W_fa在两种控制方法中的响应曲线几乎重合。

图8（a）-图8（b）为H=4km,Ma=0.54状态下的大攻角机动飞行仿真，各图物理意义与图7相同。从图8(b)可以看出，在t=2s时通过操纵升降舵，使飞机由平飞进入大攻角机动，而后分别在t=5s,8s,12s,15s对升降舵偏转角δ设置阶跃，各仿真段均能基于喘振裕度估计模型很好地实现航空发动机高稳定性控制。从图8(a)结果可明显看出，高稳定性控制方案能够准确地把S_m,f控制在11%～13%之间，使航空发动机在机动过程中稳定性得到保证，且具有较高效率，这正体现了本发明方案的优越性及有效性。

Claims (6)

1.一种基于喘振裕度估计模型的航空涡扇发动机控制方法，其特征在于，发动机控制系统对下一时刻的发动机喘振裕度进行估计，具体按照以下方法：

步骤1、利用攻角预测模型对下一时刻的航空器攻角进行在线预测；

步骤2、利用预测的下一时刻攻角确定下一时刻发动机的压缩部件的喘振压比损失                                               

；

步骤3、根据下式确定下一时刻发动机喘振裕度损失量：

  ，

其中，k、k+1分别表示当前时刻和当前时刻的下一时刻；

表示k+1时刻的发动机喘振裕度损失量；

表示k+1时刻的发动机的压缩部件的喘振压比损失；表示以当前时刻的发动机可测参数作为输入，利用无畸变喘振裕度模型估计出的k+1时刻的无畸变喘振裕度，所述无畸变喘振裕度模型为预先通过离线训练得到的神经网络；

步骤4、按照下式得到下一时刻的发动机喘振裕度估计值

：

  ，

其中，表示利用无畸变喘振裕度模型估计出的k+1时刻的无畸变喘振裕度，

表示步骤3所得到的k+1时刻的发动机喘振裕度损失量；

当所得到的超出预设的工作裕度范围时，发动机控制系统对发动机涡轮落压比的控制指令

进行喘振压比损失补偿：当低于预设的工作裕度范围时，补偿后的发动机涡轮落压比的控制指令为

；当

高于预设的工作裕度范围时，补偿后的发动机涡轮落压比的控制指令为

；

为所述喘振压比损失补偿的补偿量，根据下式得到：

  ，

式中，

表示k+1时刻的发动机的压缩部件的喘振压比损失；

、

分别表示发动机的涡轮落压比、风扇落压比。

2.如权利要求1基于喘振裕度估计模型的航空涡扇发动机控制方法，其特征在于，作为所述无畸变喘振裕度模型的输入的发动机可测参数，是通过最小二乘支持向量回归机算法从所述航空涡扇发动机的所有N个可测参数中筛选出的与喘振裕度相关性最强的前M个可测参数，M<N。

3.如权利要求1所述基于喘振裕度估计模型的航空涡扇发动机控制方法，其特征在于，所述无畸变喘振裕度模型为三层BP神经网络。

4.如权利要求1所述基于喘振裕度估计模型的航空涡扇发动机控制方法，其特征在于，所述攻角预测模型，采用在线滚动最小二乘支持向量回归机算法在线训练得到。

5.如权利要求1所述基于喘振裕度估计模型的航空涡扇发动机控制方法，其特征在于，所述预设的工作裕度范围为11%~13%。

6.如权利要求1所述基于喘振裕度估计模型的航空涡扇发动机控制方法，其特征在于，所述利用预测的下一时刻攻角确定下一时刻发动机的压缩部件的喘振压比损失

，具体如下：首先根据下式利用预测的下一时刻攻角确定下一时刻的

进气畸变量：

 ，

式中，

表示k+1时刻的

进气畸变量，

表示飞行马赫数，

表示所述攻角预测模型在线预测出的k+1时刻的攻角；

然后根据以下公式将进气畸变量转换为喘振压比损失

：

式中，

为发动机压缩部件的喘振边界压比，

为发动机对周向畸变的敏感系数， Const为常数。

Images