一种航空发动机宽范围模型预测控制方法

1.本发明涉及航空发动机控制方法领域，具体为一种航空发动机宽范围模型预测控制方法。


背景技术：

2.现有的控制方法为协调推力管理和约束管理，广泛采用的航空发动机控制架构是线性控制器加上min-max选择器的结构。其结构示意图如图1所示。
3.通过最小最大这两层选择器的操作，能够保证发动机控制器产生的控制信号不超过最小和最大约束的限制。由于min-max控制器具有控制逻辑严密，容易实现的特点，其在民用航空发动机中得到了广泛的运用。但这种传统的结构在本质上是保守的，这是由于在传统的约束管理中，采用的限制调节器都是简单的线性控制器，并且是相对独立设计的。而min-max选择器的采用，使得最终控制量不仅取决于各线性控制器输出值与设定值的偏差大小，还与偏差变化的速率有关。其会导致某一控制器在其输出值还未达到设定边界时就被提前激活，代替主控制器发挥作用，这种限制控制器在没有达到限制值时就开始工作的现象就是传统约束管理保守性的具体体现。这制约了发动机潜能的进一步发挥。在紧急情况下，比如需要在极短的时间内提供更大推力的时候，发动机过渡态控制中响应速度将被限制，进而导致事故的发生。


技术实现要素：

4.本发明的目的在于提供一种航空发动机宽范围模型预测控制方法。该方法可以在飞行包线内较宽范围实现多输入多输出控制；同时在保证控制发动机安全性的前提下能够避免传统控制器由于存在min-max选择架构导致的控制保守性问题。为实现上述目的，本发明提供如下技术方案：
5.一种航空发动机宽范围模型预测控制方法，包括以下步骤：
6.第一步，构建单点多变量模型预测控制器。使得控制器能够在控制器设计点附近局有较好的控制效果。
7.第二步，并列多个单点多变量模型预测控制器，使用分段控制器调度模块调度多个单点多变量模型预测控制器扩展控制器应用范围，形成标称点控制器。
8.第三步，使用k-means算法在航空发动机飞行包线中划分控制区域，在每个控制区域中设计标称点控制器，最终实现全飞行包线内的宽范围模型预测控制器。
9.作为对本发明一种航空发动机宽范围模型预测控制方法的改进，整个控制器为三层嵌套组合形式，最底层的控制器为单点多变量控制器，其能够在折合转速点附近控制发动机状态；第二层控制器为标称点控制器，通过设计的折合转速调度方案，能够在选取的控制区域所确定的高度马赫数标称点附近控制发动机状态；顶层控制器即为设计的宽范围mpc控制器，通过标称点控制器调度方案能够根据进口高度，进口马赫数以及当前发动机低压轴折合转速来确定发动机控制量，整个飞行包线控制区域内的发动机状态控制。
10.在第一步中，作为对本发明一种航空发动机宽范围模型预测控制方法的改进，制定控制器性能指标及参考轨迹的步骤如下：
11.将发动机输出量分为跟踪输出量y
t
和限制输出量y
l
分别进行计算。性能指标可具体定义为：发动机跟踪输出量的预测输出与设计的发动机响应轨迹值之间的偏差。同时在约束条件中，添加对控制量以及限制输出量的限制。此外，将控制增量变化量也作为性能指标中的一项，最终制定的发动机性能指标为：
[0012][0013]
式中，δpyr(k j)是从当前时刻起未来j步的参考轨迹，对应发动机的跟踪输出量归一化增量的期望响应轨迹；δpy
t
(k j)和δpy
lcor
(k j)分别表示当前时刻起未来j步跟踪输出量归一化增量以及反馈校正后限制输出量归一化增量；δpu(k i)表示从当前时刻起未来i步的预测控制量的归一化增量，其权系数λ为正定矩阵，求和符号上的n
p
和nc分别表示预测时域和控制时域，通常来说，n
p
≥nc。约束条件中，δpu
min
和δpu
max
分别表示控制量归一化增量的最大和最小值限制，用于表征实际执行机构极限位置，和表示控制量归一化变化速率的最大和最小值限制，用于表征实际执行机构变化速率的限制；δpy
lmin
和δpy
lmax
代表限制输出量归一化增量的最大最小值限制。
[0014]
参考轨迹即模型输出的期望轨迹，通过优化得到的控制量使得系统能够沿着参考轨迹最终趋近设定值。本发明采取一种根据设定值r和输出测量值确定参考轨迹的方法，其可表示为如下的指数形式：δpyr(k j)＝αjδpy
t
(k) (1-αj)δpr，式中，ts为采样周期，t为参考轨迹的时间常数，δpr为设定值的归一化增量。在每个采样时刻，构成了δpyr(k)＝[δpyr(k 1) δpyr(k 2)
…
δpyr(k n
p
)]
t
形式的期望轨迹，其组成元素每个时刻随着δpy
t
(k)更新。
[0015]
在第二步中，作为对本发明一种航空发动机宽范围模型预测控制方法的改进，设计了折合转速调度方案。采用当前转速作为调度量，目标是在相邻的两个控制器之间切换时，控制器的输出变量不应出现不连续的变化。对于一个工作点而言，其转速若位n
i-n
i 1
折合转速之间，则控制器最终输出的值是ni控制器和n
i 1
控制器输出值插值之后的值。其公式如下所示：
[0016][0017]
式中，w
f_cmd
是分段控制器最终输出的控制量，其大小根据当前转速n
l
和相邻控制器的设计点转速的插值乘以两个控制器的当前输出控制量得到。
[0018]
在第三步中，作为对本发明一种航空发动机宽范围模型预测控制方法的改进，确定了控制区域的划分原则，定义了广义距离d
*
。
[0019]
控制区域的划分原则，可描述为：归为同一区域的飞行点之间，表征其状态的状态空间模型之间的差值应尽可能的小，以便基于标称点设计的控制器能够发挥理想的控制作用。针对发动机对象而言，对于某一稳态工作点，在一定的控制规律作用下，供油量确定时，
其状态可以由飞行高度和飞行马赫数决定。可以利用进口总温和总压来具体反映不同飞行点状态空间模型的差异。引入广义距离d
*
的定义，并以此为评价准则对飞行包线进行划分。
[0020][0021]
式中，p
t2_d
，t
t2_d
分别代表标称点的进口总压和总温；d
*
表示不同飞行点时，进口总温和总压与标称点的差异，其大小可综合反映由于飞行高度，飞行马赫数不同而引起的与标称点的发动机性能差异；ε则反映了包线划分时对区域内点所期望的一致性程度，其值越小，各点的模型一致性越好，反之，则一致性越差。式中不等号的方向表明，ε表征区域内部差异的最大程度，当包线内某一状态点的广义距离满足以上不等式时，可将其划分在以该标称点为代表的子区域内。
[0022]
在第三步中，作为对本发明一种航空发动机宽范围模型预测控制方法的改进，设计了一种使用kmeans算法的控制区域划分方法。
[0023]
首先利用随机数覆盖控制区域，并将之作为数据对象，其为飞行高度和飞行马赫数组成的向量。从数据对象中随机选定初始的聚类中心，之后，利用广义距离的概念，计算其余数据对象与各聚类中心的广义距离，并将它们按照广义距离的大小进行分类。在这之后，新聚类中心的值由每个聚类的平均值计算而来。以上过程不断重复，直至聚类中心变动满足平方和差小于特定值或者迭代次数满足要求为止。具体要求可由下式表示：
[0024][0025]
式中，x，y代表聚类中心的坐标，下标new，old代表新聚类中心和旧聚类中心。o为用户定义的循环结束条件。其流程图如图2所示。
[0026]
与现有技术相比，本发明的有益效果如下：
[0027]
1、采用mpc控制方法避免了传统min-max架构导致的控制保守性问题，在保证航空发动机参数不超限的前提下，能够充分发挥航空发动机的性能。
[0028]
2、采用控制器调度方案，大大减少了需要设计的单点控制器数量，通过在不同单点控制器之间进行调度，控制器工作范围可以覆盖整个飞行包线的。
[0029]
3、使用kmeans聚类分析方法对控制区域进行划分，能够在控制效果与需要的控制区域数量之间达到最优化划分方法。
附图说明
[0030]
图1为背景技术介绍中传统的基于线性控制器和min-max的控制架构
[0031]
图2为用于飞行包线划分的k-means算法流程图
[0032]
图3为单点多变量模型预测控制器的工作结构图
[0033]
图4为本发明一种实施例中提出的分段控制器的实现图
[0034]
图5为本发明一种实施例中提出的包线划分参数选取流程图
[0035]
图6为本发明一种实施例中提出的高度马赫数调度方案示意图
[0036]
图7为本发明一种实施例中提出的宽范围控制器的实现示意图
具体实施方式
[0037]
以一型双轴涡扇分开排气发动机为例.
[0038]
第一步构建单点多变量模型预测控制器。
[0039]
首先，利用小偏差线性化方法，可以得到发动机在某一稳态点的状态变量模型：
[0040][0041]
由于系数矩阵a，b，c，d中组成元素的数量级差距很大，通过设计点参数对各变量进行归一化处理。经过处理之后的发动机状态变量模型如下所示，px表示归一化后的参数：
[0042][0043]
为仿真模拟需要，再对其进行离散化，其离散形式可表示为：
[0044][0045]
为增强控制系统跟踪变化参考指令和抗干扰的能力，引入积分控制的概念。具体操作为，定义增广状态δpxa(k)＝[δpx(k) δpu(k-1)]，则离散情况下的矩阵可变换为：
[0046]ce
＝[c d]，de＝d
ꢀꢀ
(11)
[0047]
则方程组形式可以变为以下形式，此处δδpu指控制量增量的增量：
[0048][0049]
其紧凑形式为：
[0050][0051]
至此，得到了需要的预测模型。
[0052]
然后，确定控制器的性能指标和参考轨迹。将发动机输出量分为跟踪输出量y
t
和限制输出量y
l
分别进行计算。性能指标可具体定义为：发动机跟踪输出量的预测输出与设计的发动机响应轨迹值之间的偏差。同时在约束条件中，添加对控制量以及限制输出量的限制。此外，将控制增量变化量也作为性能指标中的一项，最终制定的发动机性能指标为：
[0053][0054]
式中，δpyr(k j)是从当前时刻起未来j步的参考轨迹，对应发动机的跟踪输出量归一化增量的期望响应轨迹；δpy
t
(k j)和δpy
lcor
(k j)分别表示当前时刻起未来j步跟踪输出量归一化增量以及反馈校正后限制输出量归一化增量；δpu(k i)表示从当前时刻起未来i步的预测控制量的归一化增量，其权系数λ为正定矩阵，选择的控制量为燃油量wf和风扇进口导叶角vsv，求和符号上的n
p
和nc分别表示预测时域和控制时域，通常来说，n
p
≥nc。约束条件中，δpu
min
和δpu
max
分别表示控制量归一化增量的最大和最小值限制，用于表征实际执行机构极限位置，和表示控制量归一化变化速率的最大和最小值
限制，用于表征实际执行机构变化速率的限制；δpy
lmin
和δpy
lmax
代表限制输出量归一化增量的最大最小值限制。
[0055]
参考轨迹即模型输出的期望轨迹，通过优化得到的控制量使得系统能够沿着参考轨迹最终趋近设定值。本文采取一种根据设定值r和输出测量值确定参考轨迹的方法，其可表示为如下的指数形式：
[0056]
δpyr(k j)＝αjδpy
t
(k) (1-αj)δprꢀꢀ
(15)
[0057]
式中，ts为采样周期，t为参考轨迹的时间常数，δpr为设定值的归一化增量。在每个采样时刻，构成了δpyr(k)＝[δpyr(k 1) δpyr(k 2)
…
δpyr(k n
p
)]
t
形式的期望轨迹，其组成元素每个时刻随着δpy
t
(k)更新。
[0058]
至此，确定了控制器的性能指标和参考轨迹。
[0059]
再然后，设计反馈矫正。由于预测模型选取的是某一稳态点附近的线性化模型，就导致在实际的控制过程中，会出现模型失匹的情况。针对这一问题，需要对预测模型的输出进行反馈校正。反馈校正的方法是多样的本文中采用根据实际输出与预测模型输出的误差来对预测输出进行校正的方法。
[0060]
每个时刻的预测校正输出可表示为：
[0061][0062]
可以简化表示为：
[0063][0064]
上式中，代表当前时刻发动机实际输出与预测模型输出的误差的归一化形式，代表当前时刻发动机实际输出值的增量归一化形式，δpy(k)是相应预测模型的输出。下标“cor”表示校正输出。h
cor_l
为校正系数矩阵，其元素选取在0-1之间选择，在本文中，选择取h1＝1。在目标函数表达式中，其代表了误差的权重系数。
[0065]
至此，反馈校正设计完毕。
[0066]
最后，设计控制算法。mpc算法过程可简单描述为：在当前k时刻，以目标函数j在满足约束条件的情况下达到最小为目标，对控制输入进行优化求解。具体解释，就是保证跟踪输出量δpy
t
(k j)，j＝1，2，...，n
p
的预测值与期望响应的误差及控制量变化值之和最小的前提下，计算当前时刻起控制时域内最优的控制量序列δpu(k i)，i＝0，1，...，n
u-1。在下一个时刻到来时，重复以上过程。
[0067]
将跟踪输出和限制输出都用输入量的归一化形式δpu(k)表示。
[0068][0069]
预测输出是根据预测模型的迭代得到的。这里因为控制时域nu的关系，当i≥nu时，δδpu(k i)＝0。
[0070]
则状态量归一化增量形式的预测输出可以写作以下形式：
[0071][0072]
同理，跟踪量归一化增量形式的预测输出为：
[0073][0074]
限制量归一化增量形式的预测输出与跟踪量归一化增量形式的预测输出类似。出于表达简洁的考虑，将复杂的表达式表示为紧凑格式。
[0075][0076]
式中，预测状态量的归一化形式为：
[0077]
δpxa(k)＝[δpxa(k 1) δpxa(k 2)
…
δpxa(k n
p
)]
t
；
[0078]
预测输出跟踪量的归一化增量形式为：
[0079]
δpy
t
(k)＝[δpy
t
(k 1) δpy
t
(k 2)
…
δpy
t
(k n
p
)]
t
；
[0080]
预测输出限制量的归一化增量形式为：
[0081]
δpy
l
(k)＝[δpy
l
(k 1) δpy
l
(k 2)
…
δpy
l
(k n
p
)]
t
；
[0082]
预测控制量增量的归一化增量形式为：
[0083]
δδpu(k)＝[δδpu(k 1) δδpu(k 2)
…
δδpu(k n
c-1)]
t
。
[0084]
矩阵有如下形式：
[0085][0086]
[0087][0088][0089][0090][0091]
根据上面的推导过程，性能指标可变为如下形式：
[0092][0093]
令则上式可被表示为：
[0094][0095]
约束条件可变化为以δδpu(k)表示的形式。
[0096][0097][0098]
上式可化简为：
[0099][0100]
对于输出量限制，与式δpy
lcor
(k)＝p
l
δpxa(k) h
l
δδpu(k) h
cor_l
δpe(k)联立可得：
[0101][0102]
化简得到如下形式：
[0103][0104]
式中：
[0105]
f＝[cmꢀ‑cm i
ꢀ‑
i h
l
ꢀ‑hl
]
t
ꢀꢀ
(30)
[0106]
d(k)＝[dumax(k) dumin(k) deltaumax deltaumin dymax(k) dymin(k)]
t
ꢀꢀ
(31)
[0107]
至此，将该优化问题转变成了一个典型的带约束条件的二次规划问题，存在有效的数值解。在matlab中，存在成熟的二次规划优化函数“quadprog”。在本文中利用此函数对δδpu(k)的值进行反复求解，并取控制序列的第一个量作用于被控对象，在下一时刻重新进行优化，并不断重复，实现滚动优化。控制算法设计完毕。
[0108]
第一步设计完成的单点多变量模型预测控制器的工作结构图如图3所示。
[0109]
由上面推导的算法可知dumax(k)，dumin(k)，dymax(k)，dymin(k)，m(k)与当前或上一时刻的变量值有关，因此，必须要进行在线更新。f，deltaumax，deltaumin这样不随时间变化的量可以提前进行离线计算，由此降低在线计算量。
[0110]
在第二步中，构建标称点控制器的步骤如下。
[0111]
首先，利用分段设计控制器的思想，对于80％-110％这一段大范围转速范围，选定10％折合转速作为阶梯点，针80％，90％，100％，110％折合转速这个四个点按照第一步所阐述的方法设计单点多变量控制器。
[0112]
其次，设计针对四个控制器的调度方案，采用当前转速作为调度量，目标是在相邻的两个控制器之间切换时，控制器的输出变量不应出现不连续的变化。下面举例说明本文采取的切换规律。对于一个工作点而言，其转速若位于80％-90％折合转速之间，则控制器最终输出的值是80％控制器和90％控制器输出值插值之后的值。其公式如(3.1)所示：
[0113][0114]
式中，w
f_cmd
是分段控制器最终输出的控制量，其大小根据当前转速n
l
和相邻控制器的设计点转速的插值乘以两个控制器的当前输出控制量得到。对公式的分析可以观察到，当转速从80％逐渐增加到90％时，控制量的变化在两个控制器的输出间是平滑切换的。
[0115]
具体架构图如图4所示。
[0116]
在第三步中，构建宽范围控制器的具体实施的方法为：
[0117]
首先，采取分段式控制器设计，先将飞行包线划分为多个区域，在其中各选择一个标称点。
[0118]
对于飞行包线的划分原则为，归为同一区域的飞行点之间，表征其状态的状态空间模型之间的差值应尽可能的小，以便基于标称点设计的控制器能够发挥理想的控制作用。对于某一稳态工作点，在一定的控制规律作用下，供油量确定时，其状态可以由飞行高度和飞行马赫数决定。由于飞行高度和飞行马赫数可转化为进口总温和总压的函数，由此，可以利用进口总温和总压来具体反映不同飞行点状态空间模型的差异。引入广义距离d*的定义，并以此为评价准则对飞行包线进行划分。
[0119][0120]
式中，分别代表标称点的进口总压和总温；d
*
表示不同飞行点时，进口总温和总压与标称点的差异，其大小可综合反映由于飞行高度，飞行马赫数不同而引起的与标称点的发动机性能差异；ε则反映了包线划分时对区域内点所期望的一致性程度，其值越小，各点的模型一致性越好，反之，则一致性越差。式中不等号的方向表明，ε表征区域内部差异的最大程度，当包线内某一状态点的广义距离满足以上不等式时，可将其划分在以该标称点为代表的子区域内。
[0121]
根据广义距离概念对实际飞行包线的划分主要取决于两个参数，一是ε的大小，二是标称点的个数。参数的选取可以由图5表示。
[0122]
确定好初始参数之后，本方法采用k-means聚类分析方法，对区域进行划分。先利用随机数覆盖控制区域，并将之作为数据对象，其为飞行高度和飞行马赫数组成的向量。从数据对象中随机选定初始的聚类中心，之后，利用广义距离的概念，计算其余数据对象与各聚类中心的广义距离，并将它们按照广义距离的大小进行分类。在这之后，新聚类中心的值由每个聚类的平均值计算而来。以上过程不断重复，直至聚类中心变动满足平方和差小于特定值或者迭代次数满足要求为止。具体要求可由式(34)表示：
[0123][0124]
式中x，y代表聚类中心的坐标，下标new，old代表新聚类中心和旧聚类中心。o为用户定义的循环结束条件。
[0125]
其流程图如图2所示。
[0126]
然后，设计一种调度方案来覆盖飞行包线内的宽飞行范围。
[0127]
为避免区域切换时所带来的控制量突变问题，在不同飞行高度和马赫数下，需要借助相邻控制器的输出进行插值。不同高度马赫数可用标称点来代替，即需要对相邻标称点控制器的输出进行插值选择。又因为飞行高度和马赫数是进口总温和总压的函数。由此，可通过进口总温和总压来代替转速对相邻控制器的输出进行插值选择。为进一步减少变量数，更接近转速调度方案情况，利用广义距离定义，将进口总温和总压进一步转化为“距离”定义。根据当前状态点与相邻标称点的“距离”来对控制器输出量进行插值选择。
[0128]
为明确说明本文采用的高度马赫数调度方案，以0.23马赫数，11000m高度处，低压转子折合转速为85％的双轴涡扇发动机状态点为例。具体方案的示意图如图6所示。
[0129]
0.23马赫数，11000m高度这个状态点，位于两个标称点(在此例中为标称点2和3)
之间，可以通过进口总温总压公式与广义距离公式计算得到其与相邻两个标称点的“距离”。定义如下：
[0130][0131]
式中，下标n表示其是与n号标称点的广义距离。
[0132]dn*
的大小表征了选取的状态点和n号标称点工作状态的接近程度，d
n*
的值越小，则状态点与n号标称点的状态越接近。至此，可将d
n*
类比式32中的n
l-n
l80％
，则选取的状态点对应的控制量可表示为：
[0133][0134]
对上述公式进行分析。若此时状态点离2号标称点较近，则d
2*
值较小，(w
f_3-w
f_2
)前系数较小，2号标称点的控制量w
f_2
在最后输出的控制量中占主导作用，反之，当状态点接近3号标称点时，(w
f_3-w
f_2
)前系数接近1，则最后输出的控制量也越接近3号标称点控制器的输出量w
f_3
。且在经过两个标称点的过程中，w
f_cmd
的值是在w
f_2
和w
f_3
之间平滑过渡的。需要注意的是，以上只是对该调度方案的一个举例说明，并不意味着其只适用于2号标称点和3号标称点的切换过程，其适用于任意标称点之间的切换。
[0135]
最终宽范围模型预测控制器的结构如图7所示。