航空发动机非线性控制系统建模方法与流程

本发明属于航空发动机控制技术领域，特别是一种航空发动机非线性控制系统建模方法。

背景技术：

航空发动机是一种高度复杂和精密的装备，是深度融合机械、热力学、电气、控制、材料等学科的产物。控制系统是航空发动机的“神经中枢”，控制系统的可靠性和安全性对航空发动机的安全运行至关重要。模型是控制的基础，航空发动机由于工作环境恶劣且工作条件多变，导致其被控对象是时变的，因此更加需要一种快速的被控对象建模方法。

传统的建模方法主要集中于对部件机理模型的研究。在分析被控系统构成的基础上，通过已知的物理和数学规律，推导出其输出变量与输入变量之间的数学关系，并以此作为其机理模型。在此基础上，考虑到实际发动机设计制造的差异性，以及某些实际存在但尚未通过数学公式进行描述的未建模特性，利用系统辨识方法对航空发动机控制系统进行建模的研究也层出不穷。传统系统辨识方法主要基于线性系统理论，无法有效描述实际系统中存在的一些非线性特性如饱和、死区、滞环等，这降低了系统辨识的精确度。利用非线性系统辨识方法进行航空发动机控制系统建模，充分考虑实际系统中存在的非线性特性，对于建立航空发动机控制系统的精确模型至关重要。另一方面，目前的非线性系统辨识方法主要集中于神经网络、支持向量机、hammerstein-wiener模型等方面的研究，在系统辨识的快速性和鲁棒性上表现较差，且难以在线辨识。航空发动机控制系统工况多变、可靠性要求极高，因此急需一种简单、快速和可靠的非线性系统在线辨识方法。

在背景技术部分中公开的上述信息仅仅用于增强对本发明背景的理解，因此可能包含不构成在本国中本领域普通技术人员公知的现有技术的信息。

技术实现要素：

针对现有技术中存在的问题，本发明提出一种航空发动机非线性控制系统建模方法，根据被控系统的实测输入和输出数据，利用一种非线性自适应滤波器进行在线系统辨识。本滤波方法由线性滤波部分和非线性滤波部分串联而成，其中线性部分是一个横向自适应fir滤波器，用于描述实际系统的线性特性；非线性部分是一个自适应三次样条插值器，用于描述实际系统的非线性特性。利用滤波器输出与实际系统输出之间的误差对滤波器参数进行自适应，采用动量随机梯度下降法进行参数迭代，从而加速算法收敛过程，提高非线性系统建模的快速性。

本发明的目的是通过以下技术方案予以实现，一种航空发动机非线性控制系统建模方法包括以下步骤：

第一步骤中，设定用于辨识控制系统输入和输出的非线性自适应滤波器初值，其中，设定用于描述控制系统线性特性的线性fir滤波器初始权值、设定用于描述非线性特性的样条插值器的样条插值点初始值及其两者的学习步长和动量因子，

第二步骤中，记录控制系统当前及历史时刻的输入，基于fir滤波器输入向量与fir滤波器权值进行内积运算，得到fir滤波输出，，

第三步骤中，基于所述fir滤波输出得到其在样条插值器中的局部区间索引以及局部归一化坐标，

第四步骤中，基于所述局部区间索引以及局部归一化坐标生成归一化向量坐标，并提取局部区间插值点的纵坐标向量，

第五步骤中，基于所述归一化向量坐标和纵坐标向量生成样条插值器输出，

第六步骤中，利用传感器采集控制系统当前时刻的输出，并与所述样条插值器输出进行对比，得到滤波误差，

第七步骤中，根据滤波误差分别计算当前时刻fir滤波器权值及样条插值点的梯度，

第八步骤中，根据梯度分别计算fir滤波器权值的梯度动量以及样条插值点的梯度动量，

第九步骤中，分别更新fir滤波器权向量和样条插值点坐标，

第十步骤中，继续从第二步骤开始进行下一时刻的自适应过程。

所述的方法中，被辨识的航空发动机控制系统包括三阶线性模型和非线性函数，其中三阶线性模型的传递函数如下式所示：

非线性函数为y＝5arctan(0.5s)，其中s为非线性函数输入，y为非线性函数输出，建模信号服从[-0.1，0.1]的均匀分布。

所述的方法中，第一步骤中，设定fir滤波器长度为m＝100，初始权值为w(0)＝[1，0，...，0]^t，样条插值点的横坐标为q^x(0)＝[-10，-9，...，0，...，9，10]^t，插值点均匀分布，插值点横坐标间隔δx＝1，插值点横坐标长度为21，样条插值点纵坐标的初始值与横坐标一致，即q^y(0)＝q^x(0)，线性fir滤波器学习步长μw为0.2，样条插值器的学习步长μq为0.2，动量因子β＝0.9。

所述的方法中，第二步骤中，所述输入向量为x(n)＝[x(n)，x(n-1)，…，x(n-m+1)]t，其中x(n)为n时刻的系统输入量，m是fir滤波器的长度，fir滤波器的输出为滤所述输入向量与fir滤波器权值向量的内积，即s(n)＝w^t(n)·x(n)，其中，n表示n时刻，上标t表示向量的转置，粗体变量为向量，w(n)为fir滤波器权值向量。

所述的方法中，第三步骤中，在样条插值点中查表寻找fir滤波器的输出s(n)所在的局部区间索引i(n)及归一化坐标u(n)，其中，局部区间索引通过求得，其中δx是均匀分布的样条插值点横坐标之间的间隔，是向下取整运算符，n+1是样条插值点的个数，归一化坐标通过求得，最终得到的区间索引和归一化坐标的取值范围分别为2≤i(n)≤n-2和0≤u(n)＜1。

所述的方法中，第五步骤中，样条插值器输出其中c表示样条插值基矩阵，归一化坐标向量u(n)＝[u³(n)，u²(n)，u(n)，1]^t以及局部区间插值点的纵坐标向量其中，i(n)表示n时刻输入数据所属的样条插值区间的索引，u(n)表示n时刻输入数据所属的样条插值区间的归一化横坐标，上标y表示样条插值点的纵坐标。

所述的方法中，第七步骤中，滤波误差e(n)＝d(n)-y(n)，其中，d(n)为传感器采集控制系统当前时刻的输出，样条插值器输出y(n)，当前时刻fir滤波器权值w(n)的梯度为及样条插值点的梯度其中，u(n)为归一化坐标向量，为其导数，c表示样条插值基矩阵，x(n)为滤波器输入向量，δx为样条插值点横坐标间隔。

所述的方法中，第八步骤中，fir滤波器权值w(n)的梯度动量mw(n)，以及样条插值点的梯度动量其中，β为动量因子，和分别为n时刻fir滤波器权值w(n)以及样条插值点坐标的梯度。

所述的方法中，第九步骤中，fir滤波器向量按照w(n+1)＝w(n)-μwmw(n)进行更新，样条插值点向量按照进行更新，其中μw和μq分别为fir滤波器向量与样条插值点向量的自适应更新学习步长。

所述的方法中，记录控制系统当前及历史时刻的输入中，控制系统测量噪声为零均值的高斯白噪声序列，信噪比为30db，控制系统采样时间为0.5ms。

和现有技术相比，本发明具有以下优点：

本发明充分考虑航空发动机控制系统中常见的非线性特性，利用本发明中的非线性自适应滤波器可以简单、快速、可靠地对航空发动机控制系统进行在线建模；利用动量随机梯度下降法进行滤波器参数更新，加快了滤波器收敛的速度，提升了在线建模的速度。

附图说明

通过阅读下文优选的具体实施方式中的详细描述，本发明各种其他的优点和益处对于本领域普通技术人员将变得清楚明了。说明书附图仅用于示出优选实施方式的目的，而并不认为是对本发明的限制。显而易见地，下面描述的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。而且在整个附图中，用相同的附图标记表示相同的部件。

在附图中：

图1为本发明的一种数据驱动的航空发动机非线性控制系统建模方法的结构图；

图2为本发明的一个实施例的控制系统辨识的输出结果图；

图3(a)为本发明的一个实施例的辨识结果中线性滤波器的权值图；

图3(b)为本发明的一个实施例的辨识结果中非线性曲线的拟合图；

图4为本发明的一个实施例中，本发明所述方法与两种传统的建模方法的性能对比图。

以下结合附图和实施例对本发明作进一步的解释。

具体实施方式

下面将参照附图1至附图4更详细地描述本发明的具体实施例。虽然附图中显示了本发明的具体实施例，然而应当理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更透彻地理解本发明，并且能够将本发明的范围完整的传达给本领域的技术人员。

需要说明的是，在说明书及权利要求当中使用了某些词汇来指称特定组件。本领域技术人员应可以理解，技术人员可能会用不同名词来称呼同一个组件。本说明书及权利要求并不以名词的差异来作为区分组件的方式，而是以组件在功能上的差异来作为区分的准则。如在通篇说明书及权利要求当中所提及的“包含”或“包括”为一开放式用语，故应解释成“包含但不限定于”。说明书后续描述为实施本发明的较佳实施方式，然所述描述乃以说明书的一般原则为目的，并非用以限定本发明的范围。本发明的保护范围当视所附权利要求所界定者为准。

为便于对本发明实施例的理解，下面将结合附图以具体实施例为例做进一步的解释说明，且各个附图并不构成对本发明实施例的限定。

为了更好地理解，如图1所示，一种航空发动机非线性控制系统建模方法，所述方法包括以下步骤：

第一步骤s1中，设定用于辨识控制系统输入和输出的非线性自适应滤波器初值，其中，设定用于描述控制系统线性特性的线性fir滤波器初始权值、设定用于描述非线性特性的样条插值器的样条插值点初始值及其两者的学习步长和动量因子，

第二步骤s2中，记录控制系统当前及历史时刻的输入，基于fir滤波器输入向量与fir滤波器权值进行内积运算，得到fir滤波输出，

第三步骤s3中，基于所述fir滤波输出得到其在样条插值器中的局部区间索引以及局部归一化坐标，

第四步骤s4中，基于所述局部区间索引以及局部归一化坐标生成归一化向量坐标，并提取局部区间插值点的纵坐标向量，

第五步骤s5中，基于所述归一化向量坐标和纵坐标向量生成样条插值器输出，

第六步骤s6中，利用传感器采集控制系统当前时刻的输出，并与所述样条插值器输出进行对比，得到滤波误差，

第七步骤s7中，根据滤波误差分别计算当前时刻fir滤波器权值及样条插值点的梯度，

第八步骤s8中，根据梯度分别计算fir滤波器权值的梯度动量以及样条插值点的梯度动量，

第九步骤s9中，分别更新fir滤波器权向量和样条插值点坐标，

第十步骤s10中，继续从第二步骤开始进行下一时刻的自适应过程。

本发明采用自适应三次样条插值器逼近实际系统的非线性特性，通过动量随机梯度下降法提高参数自适应的收敛速度。如图2所示为本实施例中的控制系统辨识的输出结果图。可以看出，在系统运行1s以后，滤波器输出与系统实际输出基本一致，滤波误差急剧减小并趋近于零。因此，本发明所述的非线性系统建模方法可以快速、精确地对航空发动机非线性控制系统进行辨识。

所述的方法的优选实施方式中，被辨识的航空发动机控制系统包括三阶线性模型和非线性函数，其中三阶线性模型的传递函数如下式所示：

非线性函数为y＝5arctan(0.5s)，其中s为非线性函数输入，y为非线性函数输出，建模信号服从[-0.1，0.1]的均匀分布。

所述的方法的优选实施方式中，第一步骤s1中，设定fir滤波器长度为m＝100，初始权值为w(0)＝[1，0，...，0]^t，样条插值点的横坐标为q^x(0)＝[-10，-9，...，0，...，9，10]^t，插值点均匀分布，插值点横坐标间隔δx＝1，插值点横坐标长度为21，样条插值点纵坐标的初始值与横坐标一致，即q^y(0)＝q^x(0)，线性fir滤波器学习步长μw为0.2，样条插值器的学习步长μq为0.2，动量因子β＝0.9。

所述的方法的优选实施方式中，第二步骤s2中，所述输入向量为x(n)＝[x(n)，x(n-1)，…，x(n-m+1)]^t，其中x(n)为n时刻的系统输入量，m是fir滤波器的长度，fir滤波器的输出为滤所述输入向量与fir滤波器权值向量的内积，即s(n)＝w^t(n)·x(n)，其中，n表示n时刻，上标t表示向量的转置，粗体变量为向量，w(n)为fir滤波器权值向量。

所述的方法的优选实施方式中，第三步骤s3中，在样条插值点中查表寻找fir滤波器的输出s(n)所在的局部区间索引i(n)及归一化坐标u(n)，其中，局部区间索引通过求得，其中δx是均匀分布的样条插值点横坐标之间的间隔，是向下取整运算符，n+1是样条插值点的个数，归一化坐标通过求得，最终得到的区间索引和归一化坐标的取值范围分别为2≤i(n)≤n-2和0≤u(n)＜1。

所述的方法的优选实施方式中，样条插值器输出其中c表示样条插值基矩阵，归一化坐标向量u(n)＝[u³(n)，u²(n)，u(n)，1]^t以及局部区间插值点的纵坐标向量其中，i(n)表示n时刻输入数据所属的样条插值区间的索引，u(n)表示n时刻输入数据所属的样条插值区间的归一化横坐标，上标y表示样条插值点的纵坐标。

所述的方法的优选实施方式中，第七步骤s7中，滤波误差e(n)＝d(n)-y(n)，其中，d(n)为传感器采集控制系统当前时刻的输出，样条插值器输出y(n)，当前时刻fir滤波器权值w(n)的梯度为及样条插值点的梯度其中，u(n)为归一化坐标向量，为其导数，c表示样条插值基矩阵，x(n)为滤波器输入向量，δx为样条插值点横坐标间隔。

所述的方法的优选实施方式中，第八步骤s8中，fir滤波器权值w(n)的梯度动量mw(n)，以及样条插值点的梯度动量其中，β为动量因子，和分别为n时刻fir滤波器权值w(n)以及样条插值点坐标的梯度。

所述的方法的优选实施方式中，第九步骤s9中，fir滤波器向量按照w(n+1)＝w(n)-μwmw(n)进行更新，样条插值点向量按照进行更新，其中μw和μq分别为fir滤波器向量与样条插值点向量的自适应更新学习步长。所述的方法的优选实施方式中，记录控制系统当前及历史时刻的输入中，控制系统测量噪声为零均值的高斯白噪声序列，信噪比为30db，控制系统采样时间为0.5ms。

为了进一步理解本发明，在一个实施例中，建模方法包括如下步骤：

s1，设定非线性自适应滤波器初值。设置fir滤波器权值w(0)＝[w0，w1，...，wm-1]^t，样条插值点初值和给定滤波器学习步长μw、μq以及动量因子β。其中m是fir滤波器的长度，n+1是样条插值点的个数。

s2，记录控制系统当前时刻的输入x(n)及其历史时刻的输入，求fir滤波器输出s(n)。记滤波器输入向量为x(n)＝[x(n)，x(n-1)，…，x(n-m+1)]^t，则fir滤波器的输出可表示为滤波器输入向量与fir滤波器权值向量的内积，即s(n)＝w^t(n)·x(n)。其中括号索引中的n表示n时刻，上标t表示向量的转置，粗体变量为向量。

s3，在样条插值点中查表寻找s(n)所在的局部区间索引i(n)及归一化坐标u(n)。局部区间索引可通过求得，其中δx是均匀分布的样条插值点横坐标之间的间隔，是向下取整运算符。归一化坐标可通过求得。最终得到的区间索引和归一化坐标的取值范围分别为2≤i(n)≤n-2和0≤u(n)＜1。

s4，组成局部区间内的归一化坐标向量u(n)＝[u³(n)，u²(n)，u(n)，1]^t以及局部区间插值点的纵坐标向量

s5，计算样条插值器输出其中c表示样条插值基矩阵，具体为

s6，利用传感器采集控制系统当前时刻的输出d(n)，与滤波器输出y(n)进行比较并计算滤波误差e(n)＝d(n)-y(n)。

s7，分别计算当前时刻w(n)的梯度以及的梯度二者分别表示为和其中，

s8，分别计算当前时刻的梯度动量mw(n)和fir滤波器权向量w(n)的梯度动量为局部样条插值点向量的梯度动量为

s9，更新滤波器权向量。fir滤波器向量按照w(n+1)＝w(n)-μwmw(n)进行更新。样条插值点向量按照进行更新。

s10，继续从s2开始进行下一个时刻的自适应过程。

在一个实施例中，被辨识的航空发动机控制系统由一个三阶线性模型和一个非线性函数组成。其中三阶线性模型的传递函数如下式所示：

非线性函数表示为y＝5arctan(0.5s)，其中s为非线性函数输入，y为非线性函数输出。建模信号x(n)服从[-0.1，0.1]的均匀分布，系统测量噪声为零均值的高斯白噪声序列，信噪比为30db。系统采样时间为0.5ms。

在该实施例中，第一步骤s1中，时刻n＝0。设定fir滤波器长度为m＝100，初始权值为w(0)＝[1，0，...，0]¹。样条插值点的横坐标为q^x(0)＝[-10，-9，...，0，...，9，10]¹，插值点均匀分布，间隔δx＝1，n＝20。样条插值点纵坐标的初始值与横坐标一致，即q^y(0)＝q^x(0)。滤波器学习步长为μw＝μq＝0.2，动量因子β＝0.9。

在该实施例中，第二步骤s2中，时刻n＞0。记滤波器输入向量为x(n)＝[x(n)，x(n-1)，…，x(n-m+1)]^t，则fir滤波器的输出可表示为滤波器输入向量与fir滤波器权值向量的内积，即s(n)＝w^t(n)·x(n)。其中括号索引中的n表示n时刻，上标t表示向量的转置，粗体变量为向量。

存该实施例中，第三步骤s3中，存样条插值点中查表寻找s(n)所在的局部区间索引i(n)及归一化坐标u(n)。局部区间索引可通过求得，其中δx是均匀分布的样条插值点横坐标之间的间隔，是向下取整运算符。归一化坐标可通过求得。最终得到的区间索引和归一化坐标的取值范围分别为2≤i(n)≤n-2和0≤u(n)＜1。

在该实施例中，第四步骤s4中，组成局部区间内的归一化坐标向量u(n)＝[u³(n)，u²(n)，u(n)，1]^t以及局部区间插值点的纵坐标向量

在该实施例中，第五步骤s5中，计算样条插值器输出其中c表示样条插值基矩阵。

在该实施例中，第六步骤s6中，利用传感器采集控制系统当前时刻的输出d(n)，与滤波器输出y(n)进行比较并计算滤波误差e(n)＝d(n)-y(n)。

在该实施例中，第七步骤s7中，分别计算当前时刻w(n)的梯度以及的梯度二者分别表示为和其中，

在该实施例中，第八步骤s8中，分别计算当前时刻的梯度动量mw(n)和fir滤波器权向量w(n)的梯度动量为局部样条插值点向量的梯度动量为

在该实施例中，第九步骤s9中，更新滤波器权向量。fir滤波器向量按照w(n+1)＝w(n)-μwmw(n)进行更新。样条插值点向量按照进行更新。

在该实施例中，第十步骤s10中，继续从第二步骤s2开始进行下一个时刻的自适应过程。

如图3(a)和图3(b)所示为本实施例中，本发明提出非线性自适应滤波器在收敛后，其fir滤波器权值与样条插值点与真实系统参数的对比。图3(a)是fir滤波器最终收敛得到的权值与真实系统权值的对比，可以看出，二者是一致的。图3(b)是滤波器中样条插值点所形成的曲线与实际系统非线性曲线的对比，可以看出，二者在数据范围内也是一致的。因此，本发明所述方法可以精确地描述实际系统参数，具有优良的泛化性能。

如图4所示为本实施例中，在同一系统、同一参数条件下，本发明的一种数据驱动的航空发动机非线性系统快速建模方法与两种传统建模方法的误差曲线对比情况。这一结果是经过100次蒙特卡罗试验之后得到的结果。可以看出，相比于其他两种方法，本发明所述方法在收敛速度和辨识精度上，均具有更加优良的特性。

尽管以上结合附图对本发明的实施方案进行了描述，但本发明并不局限于上述的具体实施方案和应用领域，上述的具体实施方案仅仅是示意性的、指导性的，而不是限制性的。本领域的普通技术人员在本说明书的启示下和在不脱离本发明权利要求所保护的范围的情况下，还可以做出很多种的形式，这些均属于本发明保护之列。