基于精确偏导数的航空发动机状态变量模型在线建立方法与流程

本发明属于航空宇航推进理论与工程中的系统仿真与控制领域，具体涉及一种基于精确偏导数的航空发动机状态变量模型在线建立方法。

背景技术

航空发动机数学模型在航空发动机控制、仿真及工程应用中有着巨大的价值。状态变量模型由于其计算简单，又有一定理论支撑，在航空发动机数学模型中具有代表性。

鉴于状态变量模型的重要性，目前已经发展出了两种主要的建立航空发动机状态变量模型方法，即拟合法和偏导数法。拟合法包括利用遗传算法、最小二乘法等优化算法直接拟合相应的实验数据。偏导数法主要是利用差分方法和小扰动方法，通过部件级模型来实现相应系数矩阵的求解。

但是，这两种建模方法存在比较明显的缺陷。在实时性方面，拟合法往往需要离线实现，偏导数法则由于差分计算需要反复调用原函数，且因变量越多，重复调用次数越多，直接制约了实时性。这就使得二者难以在线更新状态变量模型。在建模精度方面，由于航空发动机的强非线性特征和宽包线工作范围，这两种建模方法难以对每一个工作点实现建模，而只能在有限个工作点建立模型，使得在未建模点只能通过增益调度、线性变参数技术等插值方式获得相应模型，对模型的精度造成了比较明显的影响。

然而，随着航空发动机控制与状态监测技术的发展，先进的性能寻优控制(psc)、模型预测控制(mpc)以及卡尔曼滤波器设计，对状态变量模型的建模提出了实时性和精度上的要求。具体而言，在这些应用中，都需要采用能表现当前工作点特性的模型以保证估计的准确性。这就要求状态变量模型具有在线更新的能力以满足实时更新的需求。例如卡尔曼滤波器估计相应的状态量时，在算法中需要使用当前工作点的状态变量模型系统矩阵以实现递推。当状态变量模型不能及时更新而不能表示当前工作点时就会使估计过程产生误差。在模型的精度方面，由于在这些应用中状态变量模型的作用是提供某些量的估计，所以建模精度越高，则估计必然越准确。例如mpc的优势就在于可以通过相应的模型提供对未来时刻输出的估计，进而优化输出序列，使得系统获得最优的评价指标。若采用的模型足够准确，就能够有效地估计当前时刻之后数个时刻的输出变化曲线，若采用有限个状态变量模型线性插值方式或者参数调度方式获得预测模型，在未建模点的精度难以保证，且难以适应发动机的诸多不确定性，也就破坏了mpc控制算法的优势。

因此，为了满足新的控制技术应用要求，需要一种能够实现快速在线更新状态变量模型参数的方法。当建模方法能够满足实时性要求时，能够实现在发动机任意工作点的建模，也就提高了建模精度。

技术实现要素：

本发明所要解决的技术问题在于克服现有技术不足，提供一种基于精确偏导数的航空发动机状态变量模型在线建立方法，利用该方法，能够通过部件级模型和偏导数模型的计算获得状态变量模型所需的各个参数，从而在任意工作点在线计算出准确的状态变量模型。

本发明基于精确偏导数的航空发动机状态变量模型在线建立方法，包括以下步骤：

步骤a、确定航空发动机状态变量模型的表达形式；

步骤b、根据状态变量模型表达形式，基于泰勒展开，获得能够表示稳态点和动态点的统一离散状态变量形式；

步骤c、根据所需要建立的离散状态模型中的变量，确定偏导数计算过程中的中间变量，基于航空发动机部件级模型和链式求导法则，采用解析法建立相应的偏导数模型；

步骤d、联合部件级模型和偏导数模型利用多次通过算法进行共同计算，获得当前工作点的状态变量模型系数矩阵和初始值，获得相应工作点的状态变量模型。

优选地，所述步骤b具体包括：

步骤b1、将航空发动机部件级模型表示成非线性状态空间形式：

式中，x表示状态量，u为控制量；

步骤b2、记k时刻的发动机模型输入量和状态量为xk,0和uk,0，式(1)在此工作点进行泰勒展开，并省略高阶项，得到如式(2)所示的状态变量模型：

式中：

下标k表明这些变量及状态变量模型矩阵随工作点k变化；

步骤b3、将式(2)所示的模型进行离散化处理：

式中：

k时刻代表的是进行状态变量建模时发动机的工作点所处的时刻，m是连续模型进行离散化以及应用离散化模型时的采样时刻；

步骤b4、将建模工作点的初始导数用离散方式表达：

代入式(4)中化简之后，得到离散状态变量模型：

式中，下标d表示该矩阵为离散状态变量模型的系数矩阵，各变量满足如下定义：

优选地，所述步骤c具体包括：

步骤c1、分析状态变量模型涉及的系数矩阵对应的因变量和自变量及其所属部件；

步骤c2、记所需的自变量为v，且v∈r^p,p≥1，因变量为f，且f(u)∈r^q,q≥1，u∈r^s,s≥1为所有中间变量；

将自变量v的微分向量转化为对角矩阵：

式中，dv为v的微分向量，为dv对角化后的矩阵；

步骤c3、根据链式法则，通过解析推导，建立中间变量的热力学参数u对自变量v的偏微分关系和因变量f对中间变量u的偏微分关系，即：

步骤c4、获得因变量f对自变量v的部件级偏导数模型：

优选地，所述步骤d具体包括：

步骤d1、初始化模型；

步骤d2、给定k时刻的飞行条件和发动机输入，并初始化各变量微分初始化值；

步骤d3、联合各部件气动热力学模型和偏导数模型进行共同计算，更新部件级模型猜值；

步骤d4、判断是否达到相应迭代次数，若是，则执行步骤d5，否则返回步骤d3；

步骤d5、进行转子动力学模型和偏导数模型的共同计算，并更新高低压转子转速作为k+1时刻的转速；

步骤d6、从部件级模型中获得相应变量的值作为状态变量模型各个变量的初始值，从偏导数模型获得状态变量模型对应的系数矩阵，构成状态变量模型并输出；

步骤d7、判断动态过程是否结束，若是，则执行步骤d8，否则返回步骤d2；

步骤d8、程序结束。

进一步优选地，步骤d2中初始化各变量微分值时，各自变量的微分初始化值均为1，其余变量全部取0。

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

(1)具有在线建立状态变量模型能力：本发明结合部件级模型，基于解析推导进行偏导数计算，能够在线获得当前工作点的状态变量模型，避免了经典拟合法和偏导数法的不足；

(2)具有任意工作点建模能力：本发明只需要联合部件级模型和偏导数模型进行在线计算即可，对稳态工作点和动态工作点均能够适用，能够实现对发动机任意工作点的状态变量模型建模；

(3)建模精度高：偏导数计算过程中大量采用解析推导，利用了变量之间的解析表达式，能够反映物理量之间的运算关系提高了建模精度，且有效避免了已有技术只能在有限个工作点建模，需通过插值获得其余工作点状态变量模型导致的模型精度下降问题。

(4)通用性和可移植性：本发明基于航空发动机部件级模型实施，对能够建立部件级模型的各类航空发动机均适用。

附图说明

图1为典型的双轴混排涡扇发动机结构图；

图2为部件级模型的动态计算流程图；

图3为本发明在线建立状态变量模型流程图；

图4为稳态点主燃油流量wfb阶跃1％时状态变量模型和部件级模型的低压转子转速响应；

图5为稳态点主燃油流量wfb阶跃1％时状态变量模型和部件级模型的压气机出口总压响应；

图6为稳态点尾喷口面积a8阶跃1％时状态变量模型和部件级模型的低压转子转速响应；

图7为稳态点尾喷口面积a8阶跃1％时状态变量模型和部件级模型的压气机出口总压响应；

图8为动态点主燃油流量wfb阶跃1％时状态变量模型和部件级模型的低压转子转速响应；

图9为动态点主燃油流量wfb阶跃1％时状态变量模型和部件级模型的压气机出口总压响应；

图10为动态点尾喷口面积a8阶跃1％时状态变量模型和部件级模型的低压转子转速响应；

图11为动态点尾喷口面积a8阶跃1％时状态变量模型和部件级模型的压气机出口总压响应。

具体实施方式

针对现有拟合法和偏导数法建立状态变量模型所存在的缺陷，本发明利用的航空发动机原理和气动热力计算方法，在部件级模型的基础上，通过将求导转化为求微分，结合链式求导法则建立相应的偏导数模型，通过部件级模型和偏导数模型共同运算，快速获得各工作点状态变量模型的初始值和系数矩阵，进而在线获得任意工作点的状态变量模型。其适用对象为能够建立复杂解析模型的各种动力机械，包括但不限于涡喷发动机、涡扇发动机、涡轴发动机、涡桨发动机、变循环发动机、涡轮基冲压组合发动机等。

为了便于公众理解，下面以图1所示的双轴混排涡扇发动机的状态变量模型建立过程为例来对本发明技术方案进行详细说明。

图1中的1截面为进气道进口；2截面为进气道出口和风扇进口；22截面为风扇出口；25和15截面为内外涵进口；3截面为压气机出口和燃烧室进口；4截面为燃烧室出口；41截面为高压涡轮进口；42截面为高压涡轮出口；45截面为低压涡轮进口；46截面为低压涡轮出口；16和6截面分别为外涵出口和内涵出口；8截面为尾喷管喉道；9截面为尾喷管出口。

本发明基于部件级模型实施，部件级模型的动态计算过程如图2所示，模型以风扇压比系数zf、压气机压比系数zc、高压涡轮落压比xpitg和低压涡轮落压比xpilg作为猜值，通过迭代算法修正4个猜值使高压涡轮进口流量连续方程e1、低压涡轮进口流量连续方程e2、内外涵出口静压平衡方程e3和尾喷管喉道处总压平衡e4方程四项的残差收敛。通过转子动力学方程更新每一时刻的高低压转速nf、nc。

根据步骤a，首先确定状态量、输入量以及输出量按x＝[nfnc]^t，u＝[wfba8]^t，选取，其中nf、nc为高低压转子转速，wfb为主燃油流量，a8为尾喷口喉道面积，为压气机出口总压。在输入量作用下的状态量和输出量均可由部件级模型进行计算。

根据步骤b1，把如图2所示的部件级模型，表示成非线性状态空间形式：

式中，x＝[nfnc]^t，u＝[wfba8]^t，

步骤b2.记发动机工作过程中k时刻所对应发动机模型输入量和状态量为xk,0和uk,0，式(1)在此工作点进行泰勒展开，并省略高阶项，可以得到如式(2)所示的状态变量模型：

式中：

式中，下标k表明这些变量及状态变量模型矩阵随工作点k变化，在发动机动态过程中，获得的是时变的状态变量模型。

步骤b3.根据后向差分计算方法，对当前工作点的状态变量模型进行离散化：

式中，m和m+1表示离散时刻m和离散时刻m+1，t为部件级模型的仿真步长。

将式(4)带入式(2)所示的模型，得到：

式中：

上式中，

k时刻代表的是进行状态变量建模时发动机的工作点所处的时刻，m是连续模型进行离散化以及应用离散化模型时的采样时刻；

步骤b4.同样采用后向差分，计算当前工作点初始状态的一阶导数。

化简之后，得到离散状态变量模型：

式中，下标d表示该矩阵为离散状态变量模型的系数矩阵，

由此，经过步骤b，得到了离散状态变量模型的表达形式，式(8)。所以接下来只要按照步骤c和步骤d继续实施本发明方案，获得相应的系数矩阵和初始值，就可以在任意工作点构建出上述模型。

根据步骤c1，按照上述假设选取的状态量、输入量和输出量，根据式(8)，有

根据式(9)和式(3)可知，所涉及到的自变量为高低压转子转速nf和nc、输入量wfb和a8，因变量为和其中，nf属于风扇部件，nc、属于压气机部件，wfb属于燃烧室部件，a8属于尾喷管部件。接下来，以求解对a8的偏导数(d22)为例，说明步骤c的具体实施过程。

根据步骤c2，自变量为a8，因变量为此外，当采用多次通过算法时，涉及到多次迭代过程，这就涉及到所有的部件，因此除这两个量以外的有关变量均应当视为中间变量。具体地，按部件级模型动态计算流程，a8的影响需要通过迭代算法修正猜值的方式，才能传递给因此，要计算对a8的偏导数，就需要初始化a8的微分da8。值得指出，在式(10)所示的自变量包括nf、nc、a8和wfb，因此在实际计算过程中有dnf、dnc、da8和dwfb，然后将四个微分构成一个四维微分向量并对角化后用于完整的偏导数计算。但下文仅以a8到的计算为例，所以仅给出da8有关的计算过程。

根据步骤c3，此处以da8为例，需要建立a8到平衡方程e4的微分模型cm1、平衡方程e4到猜值的微分模型cm2、当前猜值到下一次迭代猜值的微分模型cm3、猜值到风扇出口流量m22、出口总压和出口总温的微分模型cm4和风扇出口参数到的微分模型cm5。值得指出的是，在微分模型cm3中，包括了cm1、cm4和cm5，因为cm3涉及到了从风扇部件到尾喷管部件所有热力参数的微分计算。因此只需建立cm3即可，但是下文给出cm3的同时还给出cm1、cm4和cm5，以保证整个计算逻辑的完整性，使计算过程便于公众理解。

为了便于后文对偏导数模型的表述，此处给出风扇部件模型cm6、压气机部件模型cm7和尾喷管部件模型cm8。

模型cm6.风扇气动热力学部件模型：

步骤cm6.1风扇部件的气动热力计算，计算风扇折合转速

式中，为风扇进口设计总温，为风扇进口总温。

步骤cm6.2根据和风扇压比系数zf，由相应的特性图插值得到进口折合流量m2c、风扇压比πf和风扇效率ηf：

式中，g1、g2和g3为相应的风扇特性图插值函数。

步骤cm6.3.计算风扇实际进口流量m2：

式中，和为风扇进口总压和总温，为风扇进口设计总压。

步骤cm6.4.根据风扇进口总温和油气比f2由热力学关系式得到进口熵s2和进口焓h2：

f2＝0(16)

式中，g4和g5为相应的热力学函数。

步骤cm6.5.计算出口截面理想熵s22i：

s22i＝s2+log(πf)(19)

步骤cm6.6.根据油气比f22和理想熵s22i由热力学关系式g6得到出口截面理想焓h22i：

f22＝f2(20)

h22i＝g6(f22,s22i)(21)

步骤cm6.7.计算出口截面实际焓h22和出口截面总温

h22＝h2+(h22i-h2)/ηf(22)

式中，g7为相应的热力学函数。

步骤cm6.8.计算风扇出口总压和出口流量m22：

m22＝m2(25)

由于压气机部件的热力计算流程与风扇基本相似，简要地给出如下计算过程以方便后文描述。

模型cm7.压气机部件模型：

步骤cm7.1.计算压气机进口总温和总压

式中，σ25为风扇出口到压气机进口总压恢复系数。

步骤cm7.2.计算压气机出口总压

式中，g8为与风扇气动热力计算过程类似的压气机气动热力计算过程。

模型cm8.尾喷管部件模型：已知尾喷口总压恢复系数σ8和大气静压p0。

步骤cm8.1计算尾喷管进口流量m8、进口总压进口总温和进口油气比f8；

m8＝m7(29)

式中，m7、为掺混室出口燃气流量、总压和总温，wfb为主燃油流量。

步骤cm8.2.根据f8计算得到进口燃气多变系数k8和气体常数r8；

式中，g9和g10表示相应的热力学函数关系。

步骤cm8.3.计算尾喷管临界压比π8；

步骤cm8.4.当时，继续步骤cm8.5，否则转步骤cm8.6；

步骤cm8.5.计算8截面的流量系数，并转step7：

步骤cm8.6.计算8截面的流量系数

q(λ8)＝1(38)

步骤cm8.7.根据流量公式计算8截面总压

步骤cm8.8.根据进口参数求得的总压和由流量公式求解得到的应该相等，满足压力平衡条件，则有

此外，为了后面表述方便，在此假设步骤d中的迭代算法采用的是newton-raphson算法，根据此算法，采用式(41)对猜值进行修正：

式中，x＝[zfzcxpitgxpilg]^t为4个猜值构成的向量，e＝[e1e2e3e4]^t为4个平衡方程构成的残差向量，λ为迭代步长，j为e对x的雅克比矩阵。

首先建立a8到平衡方程e4的微分模型cm1。

步骤cm1.1.初始化a8的微分量da8＝1.0；

步骤cm1.2.根据式(39)，求的偏微分

步骤cm1.3.根据式(40)，求e4的偏微分de4：

接着建立平衡方程e4到猜值的微分模型cm2：

根据式(41)，求猜值x的偏微分dx：

dx＝-λj^-1(:,4)de4(44)

式中，j^-1(:,4)是j^-1的第四列，dx＝[dzfdzcdxpitgdxpilg]^t。

然后建立当前猜值到下一次迭代猜值的微分模型cm3。由于cm3涉及到的是整个部件级模型，且具体中间某一变量到另一变量的微分模型建立方法与cm1、cm2和cm4相似。因此，为简便起见，将整个从当前猜值x到下一次迭代猜值x的热力计算过程表示成式(45)和式(41)两式：

e＝t(x(a8),a8)(45)

式中，t表示从风扇部件起进行部件热力计算，直至尾喷管部件计算完毕，求得四个平衡方程残差e的热力计算过程。

据此建立cm3：

dx＝dx-λj^-1de(47)

在多次通过算法下，在k时刻的气动热力计算过程中，第一次迭代执行偏导数模型cm2，若要求的迭代次数超过1次，则后续迭代反复执行偏导数模型cm3直至达到迭代次数。

建立猜值到风扇出口参数的微分模型cm4：

步骤cm4.1.根据式(12)-(14)，求m2c、πf和ηf的偏微分dm2c、dπf和dηf：

步骤cm4.2.根据式(15)，计算风扇进口流量m2的偏微分dm2：

步骤cm4.3.根据式(19)，计算风扇出口理想焓s22i的偏微分ds22i：

步骤cm4.4.根据式(21)，计算风扇出口理想焓h22i的偏微分dh22i：

步骤cm4.5.根据式(22)-(23)，计算风扇出口实际焓h22和出口总温的偏微分dh22和

步骤cm4.6.根据式(24)-(25)，计算风扇出口总压和出口流量m22的偏微分和dm22：

dm22＝dm2(57)

最后建立风扇出口参数到的微分方程cm5：

步骤cm5.1.根据式(26)-(27)，计算压气机进口总温和总压的偏微分和

步骤cm5.2.根据式(28)，计算压气机出口总压的偏导数

由于c3设置了da8为1，步骤c4被省略。由此，就建立了a8到偏导数模型。相同的方法可以建立其余变量的偏导数模型。根据上述的建立过程可知，偏导数模型在计算过程中大量复用了部件模型热力计算参数，因此可以与部件级模型进行同步计算。例如偏导数模型cm4可以同模型cm6进行同步计算，cm5可以与cm7进行同步计算。由此可以看出，在利用部件级模型计算出相应的气动热力参数值时，可以同步计算相应的导数值。

特别指出，根据式(9)，离散形式的系数矩阵是根据连续形式的系数矩阵计算得到的，因此，对于a阵和b阵而言，可以先求得转子加速度对相应状态量和输入量的导数，即得到ak和bk，再由式(9)计算获得ak,d和bk,d。

在完成相应的偏导数模型建立后，即可以根据步骤d在线计算出任意工作点的状态变量模型。具体而言：

执行步骤d1.初始化模型；

执行步骤d2.给定k时刻的飞行条件和发动机输入，并初始化各自变量微分；

执行步骤d3.联合各部件气动热力学模型和部件偏导数模型进行共同计算，更新部件级模型猜值；

执行步骤d4.判断是否达到相应迭代次数，若是，则执行步骤d5，否则返回步骤d3；

执行步骤d5.进行转子动力学模型和转子动力学偏导数模型的共同计算，并更新高低压转子转速作为k+1时刻的转速；

执行步骤d6.从部件级模型中获得相应变量的值作为状态变量模型各个变量的初始值，从偏导数模型中获得状态变量模型对应的系数矩阵，构成状态变量模型并输出；

执行步骤d7.判断动态过程是否结束，若是，则执行步骤d8，否则返回步骤d2；

执行步骤d8.程序结束。

相应地，如图2所示的部件级模型计算过程也发生了变化。图3给出了本发明的基于精确偏导数的航空发动机状态变量模型在线建立方法的流程图。从图3中可以看出，在计算部件热力关系的同时，通过偏导数模型获得了相应的偏导数值。最后结合部件模型的计算结果和偏导数模型计算结果直接得到状态变量模型。

为了验证本发明的有效性，分别在稳态点和动态点建立式(10)所示的状态变量模型。

本文采用的验证方法是在某一稳态点或者动态点建立状态变量模型。然后对状态变量模型输入量进行一定幅值的阶跃，获得此工作点状态变量模型的输出响应，即状态变量模型响应。在发动机部件模型的同一工作点，对输入量做相同幅值的阶跃，获得部件级模型响应。比较两个模型的响应是否一致来确定建立模型的精度。

稳态点的仿真以飞行条件0高度、0马赫数、主燃油wfb为0.99kg/s和a8为0.28m²为例。动态点仿真以发动机从0高度、0马赫数爬升至高度10km、马赫数1.2，主燃油流量wfb从1.46kg/s变化至0.6kg/s、尾喷口喉道面积a8从0.2888m²变化至0.27m²的过渡过程为例，仿真步长取20毫秒，建模工作点取动态过程中的第8秒状态。此时，对于稳态仿真点有即已经处于稳态。对于动态仿真点有表明尚处于动态过程。

仿真环境为dellt5810win7旗舰版，cpu为intelxeon(tm)1650v43.6ghz。程序运行平台为vs2010旗舰版release模式。部件级模型共同工作方程采用六次通过算法求解，迭代步长λ取1。

图4到图7给出了在稳态点得到的状态变量模型响应形式，分别作wfb和a8的1％幅值的阶跃仿真，得到相应的仿真曲线。图4和图5表示主燃油流量wfb阶跃1％时状态变量模型和部件级模型的输出响应，图6和图7表示尾喷口面积a8阶跃1％时状态变量模型和部件级模型的输出响应。图中svmresponse表示状态变量模型的响应，clmresponse表示部件级模型的响应。从图中可以看出，在线建立的状态变量模型在各输入量阶跃1％的响应能够与部件级模型的响应相吻合。

表1给出了在稳态点得到的状态变量模型分别作wfb和a8的1％和2％幅值的阶跃仿真时响应的稳态误差，表中稳态误差单位为百分比。从表中可以看出，在2％的阶跃幅值内，响应的稳态误差能够小于1％。综合图4到图7以及表1，表明了本发明在线建模方法在稳态点的有效性。

表1wfb和a8阶跃时稳态误差

图8到图11给出了在动态点得到的状态变量模型，分别作wfb和a8的1％幅值的阶跃仿真，得到相应的仿真曲线。图8和图9表示主燃油流量wfb阶跃1％时状态变量模型和部件级模型的输出响应，图10和图11表示尾喷口面积a8阶跃1％时状态变量模型和部件级模型的输出响应。从图中可以看出，在各输入量分别阶跃1％的条件下，状态变量模型的响应与部件级模型的响应能够很好地吻合。

表2给出了在动态点得到的状态变量模型分别作wfb和a8的1％和2％幅值的阶跃仿真时响应的稳态误差，表中稳态误差单位为百分比。从表中可以看出，在2％的阶跃幅值内，响应的稳态误差能够小于2％。综合图8到图11以及表2，表明了本发明能够适用于动态点的状态变量模型建立。

表2wfb和a8阶跃时稳态误差

为了验证本发明建立状态变量模型过程的实时性，对每个工作点重复100次建模计算，以避免计算机自身性能变动的影响，取100次建模耗时的平均值作为最终耗时。测试结果为稳态点平均耗时1.42毫秒，动态点平均耗时1.68毫秒，远小于20ms的采样步长，满足实时性需求。