CN105676640A - 基于贝塞尔曲线的涡扇发动机加速过程控制律设计方法 - Google Patents

Abstract

本发明公开了一种基于贝塞尔曲线的涡扇发动机加速过程控制律设计方法，包括以下步骤：1、确定涡扇发动机在加速过程的起点时刻以及涡扇发动机各控制量在起点时刻和终点时刻的状态，并初始化涡扇发动机各控制量在时域空间的贝塞尔曲线阶数；2、以缩短控制量从初始值到目标值的耗时为目标，并结合相应约束条件，优化求解涡扇发动机各控制量在时域空间的贝塞尔曲线；3、改变涡扇发动机各控制量的贝塞尔曲线阶数，转至步骤2；4、反复迭代多次后，选择使得所述优化目标值最小的迭代步中所得到的涡扇发动机各控制量的最优贝塞尔曲线。相比现有技术，本发明从全局优化的角度出发，可以有效提高发动机在响应过程中的推力响应速率。

Description

基于贝塞尔曲线的涡扇发动机加速过程控制律设计方法

技术领域

本发明属于航空宇航推进理论与工程中的系统控制与仿真技术领域，具体涉及一种基于贝塞尔曲线的涡扇发动机加速过程控制律设计方法。

背景技术

在航空发动机控制系统设计中，对加速性能的要求越来越高，航空发动机是非常复杂的热力机械系统，其加速时的工作状态是强非线性，因而加速控制是航空发动机过渡态控制中最复杂的控制问题。在应急控制中，如紧急爬升、飞机复飞、追击、摆脱敌机等情况，发动机的加速性能无疑是一项极为重要指标。同时，它是研究其它过渡态控制，如进出加力、减速过程、起动过程的基础。因此，对于航空发动机加速过程的研究是航空发动机控制技术的一项重要研究内容。

近年来，非线性规划理论开始应用于航空发动机优化控制，并显现出其在求解复杂非线性对象优化问题方面的极大优势。这类算法可以直接对非线性对象进行优化控制计算，并且可以在优化控制过程中充分考虑各种约束条件，因此，很多学者在该方面做了大量的工作，成功地将罚函数法、约束变尺度法、遗传算法、二次序列规划法、可行二次序列规划法、功率提取法和动稳态法等方法应用于加速优化控制中，并取得良好效果。其研究的主要思路是，先把连续的优化问离散化，按时间顺序，在每一个离散点，建立能过提高加速性能的目标函数和能够保证发动机运行时的约束方程，运用优化算法对其进行求解，因此该方法存在一个致命缺陷，它只能寻求发动机在每个离散点取得最优值，而忽略全局离散点相互之间的联系，必然导致优化结果只能使发动机在某个阶段性能得到最优，然而，由于发动机的强非线性，加速过程中在某个阶段加速最快并不代表整个过程的加速时间最短。

发明内容

本发明所要解决的技术问题在于克服现有技术不足，提供一种基于贝塞尔曲线的涡扇发动机加速过程控制律设计方法，能够避免现有技术在优化过程单纯只寻求某个阶段性能最优的缺点，该方法从全局优化的角度出发，可以有效提高发动机在响应过程中的推力响应速率。

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

基于贝塞尔曲线的涡扇发动机加速过程控制律设计方法，包括以下步骤：

步骤1、确定涡扇发动机在加速过程的起点时刻以及涡扇发动机各控制量在起点时刻和终点时刻的状态，并初始化涡扇发动机各控制量在时域空间的贝塞尔曲线阶数；

步骤2、根据以下优化目标并结合包括涡扇发动机的风扇和压气机的喘振裕度和最高转子转速的限制、涡轮部件高温限制、富油熄火限制、控制量控制范围限制、最大控制量变化率限制的约束条件，优化求解涡扇发动机各控制量在时域空间的贝塞尔曲线，得到涡扇发动机各控制量在当前贝塞尔曲线阶数下的最优贝塞尔曲线：

minJ＝max(t_u1(n_u1),t_u2(n_u2),…,t_uM(n_uM))

其中，t_uj(n_uj)表示涡扇发动机第j个控制量uj的n_uj阶贝塞尔曲线起点到终点的时间，j＝1,2,…,M，M为涡扇发动机控制量的个数；

步骤3、改变涡扇发动机各控制量的贝塞尔曲线阶数，转至步骤2；

步骤4、反复迭代多次后，选择使得所述优化目标值最小的迭代步中所得到的涡扇发动机各控制量的最优贝塞尔曲线，即为涡扇发动机各控制量在加速过程的最优控制曲线。

优选地，所述约束条件被N点离散化为以下约束表达式：

$s.t.\left\{\begin{array}{c}1-\frac{N_f\[k+1\]\left(u\[k\]\right)}{N_{f,\max }}\≥0;1-\frac{N_c\[k+1\]\left(u\[k\]\right)}{N_{c,\max }}\≥0;\\ \frac{S_{mf}\[k+1\]\left(u\[k\]\right)}{S_{mf,\min }}-1\≥0;\frac{S_{mc}\[k+1\]\left(u\[k\]\right)}{S_{mc,\min }}-1\≥0;\\ 1-\frac{T_4\[k+1\]\left(u\[k\]\right)}{T_{4,\max }}\≥0;1-\frac{{\mathrm{FAR}}_4\[k+1\]\left(u\[k\]\right)}{{\mathrm{FAR}}_{4,\max }}\≥0;\\ u_{\min }\≤u\[k\]\≤u_{\max };\mathrm{\Δ}u\[k\]\≤{\mathrm{\Δu}}_{\max };k=1,2,\mathrm{...},N;\\ t_{uj}\[i\]\≤t_{uj}\[i+1\];i=0,1,\mathrm{...},n_{uj}-1;j=1,2,\mathrm{...},M\end{array}\right.$

其中，N_f[k+1](u[k])、N_c[k+1](u[k])、S_mf[k+1](u[k])、S_mc[k+1](u[k])、T₄[k+1](u[k])、FAR₄[k+1](u[k])分别表示第k+1个离散点处的风扇转速、压气机转子转速、风扇喘振裕度、压气机喘振裕度、涡轮前温度、燃烧室燃气比，其为第k个离散点处的控制量u[k]的函数；N_f,max、N_c,max、S_mf,min、S_mc,min、T_4,max、FAR_4,max分别为所允许的风扇最大转速、压气机转子最大转速、最小风扇喘振裕度、最小压气机喘振裕度、最大涡轮前温度、最大燃烧室燃气比；u_min、u_max分别为预设的最小、最大控制量；Δu[k]为第k个离散点处的控制量变化率；Δu_max为控制量变化率上限；M为涡扇发动机控制量的个数；t_uj[i]表示涡扇发动机第j个控制量uj的n_uj阶贝塞尔曲线的起点到第i个控制点的时间。

优选地，使用FSQP算法求解优化求解涡扇发动机各控制量在时域空间的贝塞尔曲线。

优选地，在3～7的范围内改变涡扇发动机各控制量的贝塞尔曲线阶数。

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

本发明技术方案是基于Bezier(贝塞尔)控制量曲线的优化，以缩短加速过程中发动机控制量最后稳定不变的时间为优化目标，可以有效地克服常规优化控制方法只考虑某个阶段最优的缺点，从而提高发动机推力的整体响应速度。

附图说明

图1是具体实施例中所针对的发动机界面图；

图2是推力F的响应曲线；

图3是燃油Wfb的优化结果曲线；

图4是尾喷管喉道面积A₈的优化结果曲线；

图5是涡轮前温度T₄的响应曲线；

图6是压气机转子转速N_c的响应曲线；

图7是风扇喘振裕度Smf的响应曲线；

图8是风扇转子转速N_f的响应曲线；

图9是压气机喘振裕度Smc的响应曲线；

图10是主燃烧室的油气比FAR₄的响应曲线；

图2～图10中的Wfb、A₈、F、和T₄作了归一化处理。

具体实施方式

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

针对现有技术在优化过程单纯只寻求某个阶段性能最优的缺点，本发明的思路是利用Bezier曲线设计涡扇发动机控制量在加速过程中的曲线，将发动机加速控制优化问题转化为优化发动机各控制量的Bezier曲线的数学问题，综合考虑了每个时刻控制量输入对发动机最终稳定工作的影响，从而得到全局优化的发动机控制曲线。

本发明具体包括以下步骤：

步骤1、确定涡扇发动机在加速过程的起点时刻以及涡扇发动机各控制量在起点时刻和终点时刻的状态，并初始化涡扇发动机各控制量在时域空间的贝塞尔曲线阶数。

发动机不同的加速起点和加速终点，其加速限制必然不同，因此，在进行加速优化前需先确定加速阶段的起点和终点。对于发动机的某一加速过程而言，其起点时刻以及该时刻各控制量的状态已知，终点时刻的各控制量状态也是已知的。由于不同阶数的Bezier曲线其对应的曲线方程不一样，因而需要初步设置各发动机控制变量的Bezier曲线阶数。

步骤2、根据以下优化目标并结合包括涡扇发动机的风扇和压气机的喘振裕度和最高转子转速的限制、涡轮部件高温限制、富油熄火限制、控制量控制范围限制、最大控制量变化率限制的约束条件，优化求解涡扇发动机各控制量在时域空间的贝塞尔曲线，得到涡扇发动机各控制量在当前贝塞尔曲线阶数下的最优贝塞尔曲线：

minJ＝max(t_u1(n_u1),t_u2(n_u2),…,t_uM(n_uM))

其中，t_uj(n_uj)表示涡扇发动机第j个控制量uj的n_uj阶贝塞尔曲线起点到终点的时间，j＝1,2,…,M，M为涡扇发动机控制量的个数。

由于Bezier曲线具有良好可控性，使得它在工程上得到广泛的应用。本发明即利用Bezier曲线设计航空发动机控制量的曲线，并以此作为优化对象。对于给定空间n+1个控制点P_i(i＝0,1,…n)，其Bezier曲线由以下计算公式可得：

$P\left(t\right)=\underset{i=0}{\overset{n}{\Σ}}{P}_i{B}_{i,n}\left(t\right),t\∈\[0,1\]$

其中B_i,n(t)是n次Bernstein基函数：

$B_{i,n}\left(t\right)=C_n^i{t}^i{(1-t)}^{n-1},C_n^i=\frac{n!}{(n-1)!i!},(i=0,1,...n)$

假设对象发动机共存在M个控制量，则其中第j个控制量uj的贝塞尔曲线的控制点可表示为(t_uj(0),uj(0)),(t_uj(1),uj(1)),…,(t_uj(n_uj),uj(n_uj))，n_uj为当前迭代步中控制量uj的贝塞尔曲线阶数。则各控制量的贝塞尔曲线如下：

$P_{uj}\left(t\right)=\underset{i=0}{\overset{n_{uj}}{\Σ}}{P}_i^{uj}{B}_{i,n_{uj}}\left(t\right),t\∈\[0,1\]$

其中是n_uj阶Bernstein基函数，表示涡扇发动机第j个控制量uj的控制点，j＝1,2,…,M，M为涡扇发动机控制量的个数。

发动机在运行过程中，如果控制量越早接近平稳状态，则发动机越早地平稳地运行，这意味着越短的加速时间，因此为了提高发动机的响应速度，本发明直接以缩短控制量从初始值到目标值的耗时为目标函数，可得：

$\underset{x}{min}J=max\left(t_{u1}\right(n_{u1}),t_{u2}(n_{u2}),...,t_{uM}(n_{uM}\left)\right)---\left(2\right)$

其中x为优化变量，由于加速起点固定，终点状态固定，则(t_uj(0),uj(0))和uj(n_uj)也确定，只需优化其它控制点和终点时间，所以：

x＝[(t_u1(1),u1(1)),…,(t_u1(n_u1-1),u1(n_u1-1)),

(t_u2(1),u2(1)),…,(t_u2(n_u2-1),u2(n_u2-1)),

…，

(t_uM(1),uM(1)),…,(t_uM(n_uM-1),uM(n_uM-1)),

t_u1(n_u1),t_u2(n_u2),…,t_uM(n_uM)]^H

为了保证发动机安全平稳地工作，在加速优化过程中，除了要提高发动机推力的响应速度，同时还得满足包括涡扇发动机的风扇和压气机的喘振裕度和最高转子转速的限制、涡轮部件高温限制、富油熄火限制、控制量控制范围限制、最大控制量变化率限制的约束，在优化过程中这些都应是曲线约束方程，为了减小计算量，本发明进一步对这些约束进行离散化，从而得到以下约束表达：

$s.t.\left\{\begin{array}{c}1-\frac{N_f\[k+1\](u\[k\])}{N_{f,max}}\≥0;1-\frac{N_c\[k+1\](u\[k\])}{N_{c,max}}\≥0;\\ \frac{S_{mf}\[k+1\](u\[k\])}{S_{mf,\min }}-1\≥0;\frac{S_{mc}\[k+1\](u\[k\])}{S_{mc,\min }}-1\≥0;\\ 1-\frac{T_4\[k+1\](u\[k\])}{T_{4,\max }}\≥0;1-\frac{{\mathrm{FAR}}_4\[k+1\](u\[k\])}{{\mathrm{FAR}}_{4,max}}\≥0;\\ u_{\min }\≤u\[k\]\≤u_{\max };\mathrm{\Δ}u\[k\]\≤\mathrm{\Δ}{u}_{\max };k=1,2,...,N\end{array}\right.$

其中，N_f[k+1](u[k])、N_c[k+1](u[k])、S_mf[k+1](u[k])、S_mc[k+1](u[k])、T₄[k+1](u[k])、FAR₄[k+1](u[k])分别表示第k+1个离散点处的风扇转速、压气机转子转速、风扇喘振裕度、压气机喘振裕度、涡轮前温度、燃烧室燃气比，其为第k个离散点处的控制量u[k]的函数；N_f,max、N_c,max、S_mf,min、S_mc,min、T_4,max、FAR_4,max分别为所允许的风扇最大转速、压气机转子最大转速、最小风扇喘振裕度、最小压气机喘振裕度、最大涡轮前温度、最大燃烧室燃气比；u_min、u_max分别为预设的最小、最大控制量；Δu[k]为第k个离散点处的控制量变化率；Δu_max为控制量变化率上限，N为离散点数量。

在以上约束条件基础上再结合贝塞尔曲线本身的约束，可得完整的约束条件

如下： $s.t.\left\{\begin{array}{c}1-\frac{N_f\[k+1\]\left(u\[k\]\right)}{N_{f,\max }}\≥0;1-\frac{N_c\[k+1\]\left(u\[k\]\right)}{N_{c,\max }}\≥0;\\ \frac{S_{mf}\[k+1\]\left(u\[k\]\right)}{S_{mf,\min }}-1\≥0;\frac{S_{mc}\[k+1\]\left(u\[k\]\right)}{S_{mc,\min }}-1\≥0;\\ 1-\frac{T_4\[k+1\]\left(u\[k\]\right)}{T_{4,\max }}\≥0;1-\frac{{\mathrm{FAR}}_4\[k+1\]\left(u\[k\]\right)}{{\mathrm{FAR}}_{4,\max }}\≥0;\\ u_{\min }\≤u\[k\]\≤u_{\max };\mathrm{\Δ}u\[k\]\≤{\mathrm{\Δu}}_{\max };k=1,2,\mathrm{...},N;\\ t_{uj}\[i\]\≤t_{uj}\[i+1\];i=0,1,\mathrm{...},n_{uj}-1;j=1,2,\mathrm{...},M\end{array}\right.---\left(3\right)$

其中最后一行为贝塞尔曲线本身的约束，t_uj[i]表示涡扇发动机第j个控制量uj的n_uj阶贝塞尔曲线的起点到第i个控制点的时间。

从式(2)和式(3)所构成的优化模型可以看出，该优化问题可以简化为如下具有不等式约束的非线性规划问题：

minf(x)

$s.t.\left\{\begin{array}{c}{g}_i\left(x\right)\≤0,i=1,2,...m_1\\ h_i\left(x\right)=0,i=m_1+1,...m\end{array}\right.---\left(4\right)$

式(4)中，x＝[x₁,x₂,…,x_n]^T，f(x)、g(x)和h(x)为连续可微的非线性函数。

以上优化问题可采用现有的多种算法进行求解，例如线性规划(LP)和序列二次规划(SQP)，但LP是采用分段线性优化的方法，对于强非线性的发动机加速过程来说存在精度低的缺点,SQP虽然在精度上得到改进，但其迭代点往往不可行，而可行性的序列二次规划算法(FSQP)，该算法不仅具有全局收敛和局部超线性收敛的特性，而且直接以目标函数为效益函数，可以避免罚因子的影响。为此，本发明优选采用FSQP算法进行优化求解。FSQP算法为现有技术，详细内容可参考文献[BüskensC,MaurerH.SQP-MethodsforSolvingOptimalControlProblemswithControlandStateConstraints:AdjointVariables,SensitivityAnalysisandReal-TimeControl[J].Journalofcomputationalandappliedmathematics,2000,120(1):85-108.]，其基本原理如下：

在每一步(即迭代点x_k)构造以下二次规划(QP)子问题：

$\begin{array}{cc}min& \frac{1}{2}{d}_k^T{H}_k{d}_k+{\mathrm{\γ}}_k\end{array}$

$s.t.\left\{\begin{array}{c}\▿f{\left(x_k\right)}^T{d}_k\≤{\mathrm{\γ}}_k\\ g_i\left(x_k\right)+\▿g_i{\left(x_k\right)}^T{d}_k\≤{\mathrm{\γ}}_k{\mathrm{\η}}_k,i=1,...,m_1\\ h_i\left(x_k\right)+\▿h_i{\left(x_k\right)}^T{d}_k\≤{\mathrm{\γ}}_k{\mathrm{\η}}_k,i=m_1+1,...,m\end{array}\right.---\left(5\right)$

式(5)中，H_k为矩阵，且 $\begin{array}{cc}Hessian& 0<H_k=H_k^T\&Element;R^{n\&times;n},0\&le;{\mathrm{\&eta;}}_k\&Element;R,\end{array}$ d_k为可行性搜索方向，η_k是通过求解如下QP问题来更新：

$\begin{array}{cc}min& \frac{1}{2}{\tilde{d}}_k^T{H}_k{\tilde{d}}_k+\&dtri;f{\left(x_k\right)}^T{\tilde{d}}_k^T\end{array}$

$s.t.\left\{\begin{array}{c}{g}_j\left(x_k\right)+\&dtri;g_j{\left(x_k\right)}^T{\tilde{d}}_k=0,j\&Element;I_k\\ h_i\left(x_k\right)+\&dtri;h_i{\left(x_k\right)}^T{\tilde{d}}_k=0,i\&Element;I_k\\ {\mathrm{\&eta;}}_k=C_k\&CenterDot;\left|\right|{\tilde{d}}_k|{|}^2\end{array}\right.---\left(6\right)$

式(6)中，I_k称为积极约束集，0＜C_k∈R。

步骤3、改变涡扇发动机各控制量的贝塞尔曲线阶数，转至步骤2。

步骤4、反复迭代多次后，选择使得所述优化目标值最小的迭代步中所得到的涡扇发动机各控制量的最优贝塞尔曲线，即为涡扇发动机各控制量在加速过程的最优控制曲线。

根据大量实验发现，在3～7的范围内改变涡扇发动机各控制量的贝塞尔曲线阶数，即利用穷举法在3～7的范围内调试各控制量Bezier曲线的不同阶数组合，可快速获得较优的加速控制方案。

为了验证本发明所设计的加速控制方法具有全局最优性，以普拉特·惠特尼公司生产的F100-PW-200涡扇发动机为仿真对象，进行了仿真对比实验。

图1给出了该涡扇发动机截面的示意图，图1中，Wfb表示主燃烧室燃油流量，A₈表示尾喷管喉道面积，1截面为发动机进口，2截面为风扇进口，22截面为风扇出口，13和23截面分别为外涵道、内涵道进口，25截面为压气机进口，3截面为燃烧室进口，4截面为高压涡轮进口，42截面为高压涡轮出口，45截面为低压涡轮转子进口，46截面为低压涡轮出口，16截面为外涵道出口，6截面为内涵道出口及掺混室进口，7截面为加力燃烧室进口，75截面为加力燃烧室出口，8截面为尾喷管喉道，9截面为尾喷管出口。

采用本发明方法设计该涡扇发动机的加速控制方案，具体包括以下步骤：

步骤1、确定涡扇发动机在加速过程的起点时刻以及涡扇发动机各控制量在起点时刻和终点时刻的状态，并初始化涡扇发动机各控制量在时域空间的贝塞尔曲线阶数。

步骤2、优化求解涡扇发动机各控制量在时域空间的贝塞尔曲线，得到涡扇发动机各控制量在当前贝塞尔曲线阶数下的最优贝塞尔曲线：

对于该发动机的两个控制量：主燃烧室燃油流量Wfb和尾喷管喉道面积A8，其贝塞尔曲线的控制点可以分别表示为：

(t_Wfb(0),Wfb(0)),(t_Wfb(1),Wfb(1)),…,(t_Wfb(n_Wfb),Wfb(n_Wfb))；

$(t_{A_8}\left(0\right),A_8\left(0\right)),(t_{A_8}\left(1\right),A_8\left(1\right)),...,(t_{A_8}\left(n_{A_8}\right),A_8\left(n_{A_8}\right));$

其中n_Wfb和分别表示Wfb和A₈的Bezier曲线的阶数。则相应的优化模型具体如下：

$\underset{x}{min}J=max\left(t_{Wfb}\right(n_{Wfb}),t_{A_8}(n_{A_8}\left)\right)$

$s.t.\left\{\begin{array}{c}1-\frac{N_f\&lsqb;k+1\&rsqb;\left(u\&lsqb;k\&rsqb;\right)}{N_{f,\max }}\&GreaterEqual;0;1-\frac{N_c\&lsqb;k+1\&rsqb;\left(u\&lsqb;k\&rsqb;\right)}{N_{c,\max }}\&GreaterEqual;0;\\ \frac{S_{mf}\&lsqb;k+1\&rsqb;\left(u\&lsqb;k\&rsqb;\right)}{S_{mf,\min }}-1\&GreaterEqual;0;\frac{S_{mc}\&lsqb;k+1\&rsqb;\left(u\&lsqb;k\&rsqb;\right)}{S_{mc,\min }}-1\&GreaterEqual;0;\\ 1-\frac{T_4\&lsqb;k+1\&rsqb;\left(u\&lsqb;k\&rsqb;\right)}{T_{4,\max }}\&GreaterEqual;0;1-\frac{{\mathrm{FAR}}_4\&lsqb;k+1\&rsqb;\left(u\&lsqb;k\&rsqb;\right)}{{\mathrm{FAR}}_{4,\max }}\&GreaterEqual;0;\\ t_{W_{fb}}\&lsqb;i\&rsqb;\&le;t_{W_{fb}}\&lsqb;i+1\&rsqb;;i=0,1,\mathrm{...},n_{W_{fb}}-1\\ t_{A_8}\&lsqb;j\&rsqb;\&le;t_{A_8}\&lsqb;j+1\&rsqb;;j=0,1,\mathrm{...},n_{A_8}-1\\ u_{\min }\&le;u\&lsqb;k\&rsqb;\&le;u_{\max };\mathrm{\&Delta;}u\&lsqb;k\&rsqb;\&le;{\mathrm{\&Delta;u}}_{\max };k=1,2,\mathrm{...},N\end{array}\right.$

其中x为优化变量，由于加速起点固定，终点状态固定，则(t_Wfb(0),Wfb(0)),Wfb(n_Wfb)，也确定，只需优化其它控制点和终点时间，所以 $x=\&lsqb;\left(t_{W_{fb}}\right(1),W_{fb}\left(1\right)),...,\left(t_{W_{fb}}\right(n_{W_{fb}}-1),W_{fb}(n_{W_{fb}}-1)),\left(t_{A_8}\right(1),A_8\left(1\right))...,$ $\left(t_{A_8}\right(n_{A_8}-1),A_8(n_{A_8}-1\left)\right),t_{W_{fb}}\left(n_{W_{fb}}\right),t_{A_8}\left(n_{A_8}\right){\&rsqb;}^T;$ N_f[k+1](u[k])、N_c[k+1](u[k])、S_mf[k+1](u[k])、S_mc[k+1](u[k])、T₄[k+1](u[k])、FAR₄[k+1](u[k])分别表示第k+1个离散点处的风扇转速、压气机转子转速、风扇喘振裕度、压气机喘振裕度、涡轮前温度、燃烧室燃气比，其为第k个离散点处的控制量u[k]的函数；N_f,max、N_c,max、S_mf,min、S_mc,min、T_4,max、FAR_4,max分别为所允许的风扇最大转速、压气机转子最大转速、最小风扇喘振裕度、最小压气机喘振裕度、最大涡轮前温度、最大燃烧室燃气比；u_min、u_max分别为预设的最小、最大控制量；Δu[k]为第k个离散点处的控制量变化率；Δu_max为控制量变化率上限，N为离散点数量。

利用FSQP算法对上述优化模型进行优化求解。

步骤3、在3～7的范围内改变涡扇发动机各控制量的贝塞尔曲线阶数，转至步骤2。

步骤4、反复迭代多次后，选择使得所述优化目标值最小的迭代步中所得到的涡扇发动机各控制量的最优贝塞尔曲线，即为涡扇发动机各控制量在加速过程的最优控制曲线。

本发明方法与传统方法的仿真环境都都在发动机工作在高度H＝0km、马赫数Ma＝0、标准大气压的条件下。发动机加速起点都为油门杆角度PLA＝20°时对应稳定工作状态。进行优化仿真前，还需要确定各限制因素的限制值。N_f和N_c的最大值都是设计点的102％，Smf和Smc最小值均为5％，T₄的最大值为设计点T₄的101.35％，FAR₄的最大值为0.035。燃烧室每秒最大供油量为120％Wfb_,ds，Wfb_,ds为在设计点时Wfb的值,其最大每秒增速限制值为60％Wfb_,ds，A₈的变化范围为95％～122％A_8,ds，A_8,ds为在设计点时A₈的值，A₈每秒最大增量与A_8,ds相等。

本发明算法所输出的Wfb和A₈的Bezier曲线的阶数均为5阶，其控制点的优化结果如表1所示，本发明方法与传统方法的对比结果如图2～图10所示，图中用“新方法”表示本发明方法，“传统方法”表示现有的分段优化的方法(具体可参见文献[黄伟.航空发动机高稳定性控制及其在加速控制中的应用[D],南京：南京航空航天大学，2013])。

表1Bezier控制点的优化结果

	P0	P1	P2	P3	P4	P5
							Wfb	(0,0.165)	(0.99,1.62)	(1.63,0.44)	(3.05,0.71)	(3.58,1.02)	(3.6,1)
A8	(0,1)	(0.36,1.3)	(0.48,1.35)	(1.47,0.97)	(3.47,1.01)	(3.55,1)

从表1、图3和图4可以看出，本发明的控制方案中，Wfb和A₈最后保持不变的时间分别为3.6sec和3.55sec，而常规优化最后保持不变的时间分别为5.375sec和5.425sec。控制量的越早稳定也意味着发动机推力越快的响应速度，最后推力增加到95％的最大状态推力，常规控制方法和本发明分别需要5.125sec和3.8sec。

从图2中可以看出，加速前半段常规优化方法的响应速度大于基于控制量曲线优化方法，而随着加速的运行，在加速后半段，基于控制量曲线优化方法的响应速度快于常规优化控制方法。这主要因为常规优化控制只寻求发动机在某个阶段性能最优，而忽视整体响应速度，在加速前半段，常规优化控制为了使推力快速增加，Wfb的增加速度大于本发明提出的方法，但随着加速的进行，如图9所示，越快的Wfb增加速度使发动机越快地接近喘振限制边界，进而使得Wfb在加速后半段的增加速度受阻，因而影响推力整体的响应速度。而基于控制量曲线优化方法由于充分考虑到加速前后控制量输入大小对加速时间的影响，虽然在加速前半段响应速度稍微慢于常规方法，但在加速后期起响应速度却远快于常规方法，因此可以使得在整个加速性能得以提高。

如图5～图10所示,加速优化过程中，两种控制方法都使发动机工作点紧贴着喘振裕度保护限制边界，并且不出现不超温、超转和燃烧室富油现象。

上述对比实验表明，在加速前半段，常规优化方法的推力响应速度大于本发明提出的方法，但在加速后半段，本发明提出的方法快于常规优化控制，算例中本发明提出的方法使发动机推力增加到95％的最大状态推力减少了1.325sec，因而验证了本发明比现有技术更有效。

Claims (5)

1.基于贝塞尔曲线的涡扇发动机加速过程控制律设计方法，其特征在于，包括以下步骤：

步骤1、确定涡扇发动机在加速过程的起点时刻以及涡扇发动机各控制量在起点时刻和终点时刻的状态，并初始化涡扇发动机各控制量在时域空间的贝塞尔曲线阶数；

步骤2、根据以下优化目标并结合包括涡扇发动机的风扇和压气机的喘振裕度和最高转子转速的限制、涡轮部件高温限制、富油熄火限制、控制量控制范围限制、最大控制量变化率限制的约束条件，优化求解涡扇发动机各控制量在时域空间的贝塞尔曲线，得到涡扇发动机各控制量在当前贝塞尔曲线阶数下的最优贝塞尔曲线：

minJ＝max(t_u1(n_u1),t_u2(n_u2),…,t_uM(n_uM))

其中，t_uj(n_uj)表示涡扇发动机第j个控制量uj的n_uj阶贝塞尔曲线起点到终点的时间，j＝1,2,…,M，M为涡扇发动机控制量的个数；

步骤3、改变涡扇发动机各控制量的贝塞尔曲线阶数，转至步骤2；

步骤4、反复迭代多次后，选择使得所述优化目标值最小的迭代步中所得到的涡扇发动机各控制量的最优贝塞尔曲线，即为涡扇发动机各控制量在加速过程的最优控制曲线。

2.如权利要求1所述涡扇发动机加速过程控制律设计方法，其特征在于，所述约束条件被N点离散化为以下约束表达式：

$s.t.\left\{\begin{array}{c}1-\frac{N_f\&lsqb;k+1\&rsqb;(u\&lsqb;k\&rsqb;)}{N_{f,\max }}\&GreaterEqual;0;1-\frac{N_c\&lsqb;k+1\&rsqb;(u\&lsqb;k\&rsqb;)}{N_{c,\max }}\&GreaterEqual;0;\\ \frac{S_{mf}\&lsqb;k+1\&rsqb;(u\&lsqb;k\&rsqb;)}{S_{mf,\min }}-1\&GreaterEqual;0;\frac{S_{mc}\&lsqb;k+1\&rsqb;(u\&lsqb;k\&rsqb;)}{S_{mc,\min }}-1\&GreaterEqual;0;\\ 1-\frac{T_4\&lsqb;k+1\&rsqb;(u\&lsqb;k\&rsqb;)}{N_{4,\max }}\&GreaterEqual;0;1-\frac{{\mathrm{FAR}}_4\&lsqb;k+1\&rsqb;(u\&lsqb;k\&rsqb;)}{{\mathrm{FAR}}_{4,\max }}\&GreaterEqual;0;\\ u_{\min }\&le;u\&lsqb;k\&rsqb;\&le;u_{\max };\mathrm{\&Delta;}u\&lsqb;k\&rsqb;\&le;{\mathrm{\&Delta;u}}_{\max };k=1,2,\mathrm{...},N;\\ t_{uj}\&lsqb;i\&rsqb;\&le;t_{uj}\&lsqb;i+1\&rsqb;;i=0,1,\mathrm{...},n_{uj}-1;j=1,2,\mathrm{...},M\end{array}\right.$

其中，N_f[k+1](u[k])、N_c[k+1](u[k])、S_mf[k+1](u[k])、S_mc[k+1](u[k])、T₄[k+1](u[k])、FAR₄[k+1](u[k])分别表示第k+1个离散点处的风扇转速、压气机转子转速、风扇喘振裕度、压气机喘振裕度、涡轮前温度、燃烧室燃气比，其为第k个离散点处的控制量u[k]的函数；N_f,max、N_c,max、S_mf,min、S_mc,min、T_4,max、FAR_4,max分别为所允许的风扇最大转速、压气机转子最大转速、最小风扇喘振裕度、最小压气机喘振裕度、最大涡轮前温度、最大燃烧室燃气比；u_min、u_max分别为预设的最小、最大控制量；Δu[k]为第k个离散点处的控制量变化率；Δu_max为控制量变化率上限；M为涡扇发动机控制量的个数；t_uj[i]表示涡扇发动机第j个控制量uj的n_uj阶贝塞尔曲线的起点到第i个控制点的时间。

3.如权利要求1所述涡扇发动机加速过程控制律设计方法，其特征在于，使用FSQP算法求解优化求解涡扇发动机各控制量在时域空间的贝塞尔曲线。

4.如权利要求1所述涡扇发动机加速过程控制律设计方法，其特征在于，所述涡扇发动机的控制量包括主燃烧室燃油流量和尾喷管喉道面积。

5.如权利要求1所述涡扇发动机加速过程控制律设计方法，其特征在于，在3～7的范围内改变涡扇发动机各控制量的贝塞尔曲线阶数。