CN115114731A - 一种基于伪雅可比矩阵的航空发动机动态建模方法 - Google Patents

Abstract

本发明公开了一种基于伪雅可比矩阵的航空发动机动态建模方法。本发明针对双向差分计算雅可比矩阵的部件级模型实时性较差的问题，基于非线性模型动态线性化策略，将涡轴发动机共同方程进行线性化，构建伪雅可比矩阵评价函数，通过最优估计近似计算雅可比矩阵，获得了伪雅可比矩阵的递推计算方法，避免了对部件模型的重复调用，大大提高了动态模型的实时性。

Description

一种基于伪雅可比矩阵的航空发动机动态建模方法

技术领域

本发明属于航空宇航推进理论与工程中的***控制与仿真领域，具体涉及一种航空发动机动态实时模型建模方法。

背景技术

航空发动机数值仿真技术在航空发动机研制及控制***设计、验证过程中具有极其重要的作用。以数学模型代替真实发动机开展数值仿真或半物理仿真，可以对控制算法、故障诊断技术、健康管理技术等进行初步的验证，节约实验成本、降低实验风险。自上个世纪80年代末以来，西方国家相继开展对航空发动机仿真技术的研究，并开发出多种航空发动机数值仿真***。经过多年的研究发展与应用，目前相关理论技术己经成熟，仿真置信度与精度己经达到相当高的水平。随着模型基控制技术的发展，对模型精度和实时性的要求进一步提高，也开展了大量的研究工作。

航空发动机数学模型在引入真实发动机各部件基本参数的基础上，沿发动机气路流程，按照部件特性建立发动机各部件气动热力学模型，并通过求解部件间共同工作方程，获得各部件匹配当前工作状态的参数，而共同工作方程的收敛速度及精度将直接关系到发动机数学模型的准确性和实时性。

通常求解航空发动机共同工作方程的方法有牛顿-拉夫逊(N-R)法和Broyden拟牛顿法。其中N-R法由于需要重复调用发动机部件模型计算雅可比矩阵，而严重影响了模型的实时性；而拟Newton法是将雅可比矩阵作逼近处理，并通过递推获得逼近雅可比矩阵，降低了计算复杂性，但是算法迭代收敛要求初始解的偏离程度更小，并受计算步长影响，收敛能力弱于N-R法。

基于以上两种方法，西北工业大学陈玉春等人^[1]通过改变N-R法计算步长及限制初值的方法改善了求解的收敛性问题，但仍需差分计算雅可比矩阵。南京航空航天大学廖光煌等人^[2]提出了基于神经网络的发动机共同工作方程求解方法，采集基于N-R法计算的离线训练数据，以残差为输入，共同工作方程猜值修正量作为输出，利用神经网络对其进行训练，避免了雅可比矩阵计算和方程迭代求解过程，有效提高了模型的实时性，但是神经网络采取了离线训练方式，模型精度受到网络泛化能力的影响，对不确定性的适应能力弱。南京航空航天大学王元、伍谦等人^[3-4]将N-R法平方收敛与拟Newton法超线性收敛特性结合，提出自校正Broyden拟牛顿法，通过判断模型收敛趋势调整计算步长，收敛精度及速度均有所提高，但是在收敛趋势较差的工作点，仍需要通过差分进行校正矩阵的初始化，导致个别工作点的实时性低。南京航空航天大学庞淑伟^[5]提出基于精确偏导数的发动机部件级模型建立方法，以精确偏导数代替差分计算，构造雅可比矩阵，提高了模型的实时性，但是算法计算过程复杂，实现的难度较大。

发明内容

本发明所要解决的技术问题在于克服现有技术不足，提供一种基于伪雅可比矩阵的航空发动机动态建模方法，可避免差分算法对模型的反复调用，有效提高模型的实时性。

本发明具体采用以下技术方案解决上述技术问题：

一种基于伪雅可比矩阵的航空发动机动态建模方法，通过在航空发动机工作过程中构建并求解航空发动机的动态共同工作平衡方程，获得航空发动机各部件匹配当前工作状态的参数；在求解航空发动机的动态共同工作平衡方程时，通过以下方法对动态共同工作平衡方程的猜值进行迭代修正，以使得当前的动态共同工作平衡方程满足收敛条件：

步骤1、线性化k时刻的动态共同工作平衡方程；

步骤2、基于以下的优化目标函数，计算线性化动态共同工作平衡方程的雅可比矩阵J(k)的最优估计

(伪雅可比矩阵)：

其中，ε(k)、ε(k-1)分别表示k时刻、k-1时刻的动态共同工作平衡方程残差，Δx(k)＝x(k)-x(k-1)表示k时刻与k-1时刻动态共同工作平衡方程的猜值x(k)、x(k-1)的差值，μ为惩罚因子，

表示k-1时刻的线性化动态共同工作平衡方程的雅可比矩阵的最优估计；

步骤3、通过迭代计算更新k+1时刻的动态共同工作平衡方程猜值x(k+1)：

其中，λ∈(0,1]是迭代步长。

优选地，最优估计

通过下式计算得到：

其中，η∈(0,1]是步长因子，Δε(k)＝ε(k)-ε(k-1)表示k时刻、k-1时刻的动态共同工作平衡方程残差ε(k)、ε(k-1)的差值，上标“T”表示矩阵转置，“|| ||²”表示向量2范数的平方。

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

(1)实时性高：本发明以递推伪雅可比矩阵计算替代差分形式的雅可比矩阵计算，大大减小了部件级模型调用次数，大幅度提高了动态模型计算的实时性。

(2)通用性强：伪雅可比矩阵的计算基于动态线性化方法，本身是一种数据驱动的方法，对发动机特性变化的适应能力强，可以用于多种型号的发动机动态实时模型计算。

附图说明

图1为双转子涡轴发动机的结构示意图；

图2为基于伪雅可比矩阵的涡轴发动机动态建模流程图；

图3为双向差分算法与本发明伪雅可比矩阵算法的关键发动机输出参数对比(1km高度)；

图4为双向差分算法与本发明伪雅可比矩阵算法的关键发动机输出参数对比(3km高度)。

具体实施方式

针对现有基于差分方法计算雅可比矩阵实时性差的问题，本发明结合非线性模型动态线性化方法，利用最优估计进行伪雅可比矩阵计算，提高了动态模型的实时性。

本发明具体采用以下技术方案解决上述技术问题：

一种基于伪雅可比矩阵的航空发动机动态建模方法，通过在航空发动机工作过程中构建并求解航空发动机的动态共同工作平衡方程，获得航空发动机各部件匹配当前工作状态的参数；在求解航空发动机的动态共同工作平衡方程时，通过以下方法对动态共同工作平衡方程的猜值进行迭代修正，以使得当前的动态共同工作平衡方程满足收敛条件：

步骤1、线性化k时刻的动态共同工作平衡方程；

步骤2、基于以下的优化目标函数，计算线性化动态共同工作平衡方程的雅可比矩阵J(k)的最优估计

(伪雅可比矩阵)：

其中，ε(k)、ε(k-1)分别表示k时刻、k-1时刻的动态共同工作平衡方程残差，Δx(k)＝x(k)-x(k-1)表示k时刻与k-1时刻动态共同工作平衡方程的猜值x(k)、x(k-1)的差值，μ为惩罚因子，

表示k-1时刻的线性化动态共同工作平衡方程的雅可比矩阵的最优估计；

步骤3、通过迭代计算更新k+1时刻的动态共同工作平衡方程猜值x(k+1)：

其中，λ∈(0,1]是迭代步长。

本发明动态建模方法的适用对象包括但不限于涡轴发动机、涡桨发动机、涡扇发动机、变循环发动机、组合发动机等。

为了便于公众理解，下面通过一个具体实施例并结合附图来对本发明的技术方案进行详细说明：

本实施例以图1所示双转子涡轴发动机为研究对象，该双转子涡轴发动机的燃气涡轮和压气机共轴，动力涡轮和主旋翼共轴。在工作过程中，空气从进气道进入，经压气机流向燃烧室，与燃油混合燃烧，产生高温高压燃气，燃气流经燃气涡轮时提取一部分能量来驱动压气机，离开燃气涡轮的气体流经动力涡轮，动力涡轮提取剩余能量来驱动主旋翼，满足直升机的工作需求。

对于图1所示双转子涡轴发动机，本发明所提出的动态建模方法，如图2所示，其中t为动态模型仿真时刻，k为平衡方程求解的迭代次数，T为仿真步长。该方法具体包括以下步骤：

步骤A、获取t时刻的飞行条件(直升机飞行高度H、前飞速度V_x)和燃油流量W_f；

步骤B、沿着气流流程从进气道到尾喷管依次计算各个部件热力参数及旋翼负载功率需求；

步骤C、构建涡轴发动机动态共同工作平衡方程；

在已知飞行条件、燃油流量的条件下，根据双转子涡轴发动机动态工作过程中需满足的流量连续、压力平衡条件进行分析，通常选择3个共同方程，记ε_i,i＝1,2,3为动态共同工作方程残差，包括

(1)燃气涡轮进口流量连续方程φ₁(x)：

φ₁(x)＝(W_41xs-Q_41xs)/Q_41xs＝ε₁ (1)

其中，x为动态共同工作方程中的猜值，W_41xs为由燃气涡轮部件特性曲线计算出的当前工作条件下燃气涡轮进口相似流量，Q_41xs为按气路流程计算获得的从燃气涡轮导向器进入燃气涡轮的相似流量，ε₁为方程φ₁(x)残差；

(2)动力涡轮进口流量连续方程φ₂(x)：

φ₂(x)＝(W_44xs-Q_44xs)/Q_44xs＝ε₂ (2)

其中，W_44xs为由动力涡轮部件特性曲线计算出的当前工作条件下动力涡轮进口相似流量，Q_44xs为按气路流程计算获得的从动力涡轮导向器进入动力涡轮的相似流量，ε₂为方程φ₂(x)残差；

(3)尾喷管出口压力平衡方程φ₃(x)：

φ₃(x)＝(p_c7-p₇)/p₇＝ε₃ (3)

其中，p_c7为进入尾喷管出口气流总压，p₇为喷口背压，ε₃为方程φ₃(x)残差；

所述部件级模型共同工作平衡方程中，选取猜值x＝[Z_CZ_GZ_P]^T；其中，Z_C为压气机压比系数，Z_G为燃气涡轮压比系数，Z_P为动力涡轮压比系数；以Z_c为例，压比系数定义为：

其中，π_c为当前工作点处的压气机压比，下标min和max分别代表当前转速线上压比的最小值和最大值。

步骤D、判断动态共同工作平衡方程是否满足收敛条件，如果是flag＝1，则跳到步骤F，否则flag＝0执行步骤E；

其中，k_max为最大允许的迭代次数。

本实施例中的收敛条件如下：

ε_i＜10^-5,i＝1,2,3，or k＞k_max。

步骤E、计算伪雅可比矩阵，修正猜值，返回步骤B；具体包括以下子步骤：

步骤E1、线性化k时刻的动态共同工作平衡方程：

ε(k)＝φ(x(k)) (6)

假设φ(x(k))对x的偏导数连续，且满足||ε(k)-ε(k-1)||≤b||x(k)-x(k-1)||，既如果猜值变化有限时，方程残差的变化也是有限的，当模型收敛时满足此条件；

对式(6)进行线性化，得到如下形式的线性化方程：

其中，Δε(k)＝ε(k)-ε(k-1)，Δx(k)＝x(k)-x(k-1)，

为伪雅可比矩阵；式(7)为线性化的平衡方程；

步骤E2、采用最优估计方法来计算伪雅可比矩阵

具体采用如下的评价指标：

其中，μ为惩罚因子，J(k)为雅可比矩阵；

式(8)表明当L(J(k))＝0时，ε(k)-ε(k-1)＝J(k)Δx(k)，

即线性化方程(6)成立，且对雅可比矩阵的变化进行约束；

雅可比矩阵如式(9)所示：

由于航空发动机工作过程是强非线性的复杂气动热力过程，其数学建模依赖于部件特性图和气动热力参数的插值和迭代运算，因此φ_i(x),i＝1,2,3是包含发动机各个部件计算过程的隐式方程，无法直接进行偏导数计算，只能通过数值差分的方法获得。为了提高模型精度，通常采用中间差分法来计算雅可比矩阵中的偏导数，即：

在式(10)的偏导数计算过程中，需对猜值x_i分别进行向上和向下的小扰动，调用部件级数学模型进行差分计算，共有三个猜值，故共需进行6次计算，严重影响了动态模型的实时性。为此本发明采取递推的方式进行伪雅可比矩阵的计算。

对式(8)求关于J(k)的极小值，以获得J(k)的最优估计

有：

移项整理有

在等式两边同时减去

所以有

其中，η∈(0,1]是步长因子；

当t＝0时，

通过差分进行初始化，否则，

此处下标t代表动态模型仿真时刻，即当前仿真时刻的伪雅可比矩阵初始值为上一仿真时刻迭代结束时候的伪雅可比矩阵值；

式(14)以递推的方式进行伪雅可比矩阵的计算，计算过程中仅需当前方程残差的增量和猜值的增量，不需要额外调用部件模型进行差分计算，计算过程大大简化；

步骤E3、通过迭代计算更新猜值：

其中，λ∈(0,1]是迭代步长。

步骤F、计算转子加速度，更新转速；

其中，n_G、n_P是燃气涡轮和动力涡轮转速，J_G、J_P是燃气涡轮转子和动力涡轮转子的转动惯量，η_G、η_P是燃气涡轮转子和动力涡轮转子的机械效率，W_GT、W_PT是燃气涡轮和动力涡轮功率，W_C、W_R是压气机和直升机旋翼功率。

为了验证本发明采用的伪雅可比矩阵算法在涡轴发动机求解动态共同工作方程过程中的优势，在T700涡轴发动机部件级模型上进行仿真，设置涡轴发动机初始飞行高度为H＝1km，初始前飞速度V_x＝20m/s，仿真过程中模拟发动机动态工作过程，对两种算法的模型各调用1万次。模型基于VC++6.0开发，计算机***为Windows 10家庭版，CPU i5-8250U1.60GHz,8GB RAM。伪雅可比矩阵算法模型1万次调用的平均计算时间为3.19s，而双向差分计算雅可比矩阵模型平均计算时间为6.97s。可见，伪雅可比矩阵算法所用时间约为双向差分算法的3/7，验证了算法在实时性方面的优越性。

进一步在包线范围内对模型的精度进行验证，这里给出了2组仿真结果，如图3、图4所示，其它飞行区域和发动机工作状态下的仿真结果类似。图中PJ(Pseudo-Jacobian)代表采用伪雅可比矩阵算法的部件模型仿真效果；BDD(Bi-directional differential)代表双向差分计算雅可比矩阵算法的模型仿真效果。图中给出了涡轴发动机关键参数的变化曲线，其中sfc为耗油率，T₄₁为燃气涡轮进口总温，e代表PPD算法与差分算法模型相比的输出偏差，其偏差计算公式为

(1)飞行高度1km下的仿真测试

设置直升机飞行高度H＝1km，前飞速度V_x由20m/s增至50m/s，输出参数变化曲线如图3所示。

如图3可见，相比双向差分算法，基于伪雅可比矩阵的部件级模型输出的转速偏差小于0.15％，温度偏差小于0.25％，燃油流量由于本身数量级较小，闭环控制中的最大偏差小于3％，由于燃油流量偏差引起的耗油率计算偏差小于1.2％，详细误差信息见表1，由表1和图3可见，基于伪雅可比矩阵的涡轴发动机数学模型取得了较高的精度，平均误差小于0.14％，验证了本文算法在动态模型精度方面的有效性。

表1伪雅可比矩阵模型相比双向差分模型的误差

(2)飞行高度3km下的仿真测试

设置直升机飞行高度H＝3km，前飞速度V_x由20m/s增至50m/s，输出参数变化曲线如图4所示。

由图4可见，相比双向差分算法，基于伪雅可比矩阵的部件级模型在高度3km处输出的转速偏差小于0.3％，温度偏差小于0.4％，燃油流量由于本身数量级较小，闭环控制中的最大偏差小于1.1％，由于燃油流量偏差引起的耗油率计算偏差也小于1.1％，最大建模误差优于高度1km的情况，验证了其在包线内的适用性。

参考文献：

[1]陈玉春，徐思远，杨云凯，等.改善航空发动机特性计算收敛性的方法[J].航空动力学报，2008,23(12):2242-2248.

[2]廖光煌，焦洋，李秋红，等.涡轴发动机高精度实时部件级模型研究[J].推进技术，2016，37(01):25-33.

[3]王元，李秋红，黄向华.基于自校正Broyden拟牛顿法的航空发动机模型数值计算[J].航空动力学报，2016,31(1):249-256.

[4]伍谦，李秋红，单睿斌.航空发动机气动热力***模型求解方法研究[J].计算机仿真，2019，36(01):76-81.

[5]Pang,S.；Li,Q.；Ni,B.Improved nonlinear MPC for aircraft gas turbineengine based on semi-alternative optimization strategy[J].Aerosp.Sci.Technol.2021,118,106983.

Claims (2)

1.一种基于伪雅可比矩阵的航空发动机动态建模方法，通过在航空发动机工作过程中构建并求解航空发动机的动态共同工作平衡方程，获得航空发动机各部件匹配当前工作状态的参数；其特征在于，在求解航空发动机的动态共同工作平衡方程时，通过以下方法对动态共同工作平衡方程的猜值进行迭代修正，以使得当前的动态共同工作平衡方程满足收敛条件：

步骤1、线性化k时刻的动态共同工作平衡方程；

步骤2、基于以下的优化目标函数，计算线性化动态共同工作平衡方程的雅可比矩阵J(k)的最优估计

其中，ε(k)、ε(k-1)分别表示k时刻、k-1时刻的动态共同工作平衡方程残差，Δx(k)＝x(k)-x(k-1)表示k时刻与k-1时刻动态共同工作平衡方程的猜值x(k)、x(k-1)的差值，μ为惩罚因子，

表示k-1时刻的线性化动态共同工作平衡方程的雅可比矩阵的最优估计；

步骤3、通过迭代计算更新k+1时刻的动态共同工作平衡方程猜值x(k+1)：

其中，λ∈(0,1]是迭代步长。

2.如权利要求1所述基于伪雅可比矩阵的航空发动机动态建模方法，其特征在于，最优估计

通过下式计算得到：

其中，η∈(0,1]是步长因子，Δε(k)＝ε(k)-ε(k-1)表示k时刻、k-1时刻的动态共同工作平衡方程残差ε(k)、ε(k-1)的差值，上标“T”表示矩阵转置，“||||²”表示向量2范数的平方。

Images