1 引言

航空发动机参数变化受到很多因素响，使得对其进行精确的数学模型变得非常难。由于模型简化、飞行条件的不确定性、结构参数的不确定性、模型结构的不确定性等原因，使得所建立的航空发动机模型在某种程度上是不确定的，这就使得发动机的控制器设计问题变得复杂。自适应控制方法能够有效改善因工作状况变化、外部干扰及建模误差等不确定性引起的动态跟踪问题，而航空发动机正是一个具有较大不确定性的复杂非线性系统，因此开展自适应控制方法设计及其在航空发动机的应用研究具有重要的意义^{[1, 2]}。航空发动机的大多数控制方法均基于某一基本控制器，进而实现常值指令的渐近跟踪，往往忽视了系统的不确定因素，容易导致控制品质不满足技术要求，甚至使控制系统出现不稳定等现象^{[3, 4]}。因此，开展基于基本控制器结构下的增广模型参考自适应控制方法及其在航空发动机控制中的应用研究，对提升发动机控制系统性能等方面具有重要的应用意义。

国外基于自适应控制方法对不确定性系统的跟踪控制问题进行了大量的研究且取得了一定的成果。Asha等^[5]对开环不稳定时变系统设计了一种基于李雅普诺夫的模型参考PD/PID控制器，提高了系统的跟踪性能。Stepanyan等^[6]针对状态和输出反馈相对阶数不匹配的非线性系统提出了一种间接自适应控制设计方法，降低了系统的输入输出跟踪误差。国内基于自适应控制方法对不确定性跟踪问题的相关研究亦开展了大量工作。陈志勇等^[7]针对慢变系统设计了基于增广法的自适应策略，实现了系统在参数未知下的运动轨迹渐进跟踪。周涛等^[8]提出了基于跟踪微分器的模型参考自适应控制方法来处理二阶系统模型参数的大范围不确定性问题，保证了系统的全局渐进稳定。针对航空发动机控制中的不确定性问题，丁凯峰等^[9]提出了一种基于Adaline网的发动机自适应控制方，该控制方案的稳定性无法证明。查旭等^[10]提出了一种基于Backstepping技术的鲁棒控制器。该方法通过将广义不确定项参数化，提出了状态反馈控制算法，但该方法只对SISO系统进行了讨论。Fakhari V等^[11]提出了一种基于δ修正的鲁棒自适应控制策略，并成功运用到发动机的控制中。该方法对于系统存在不确定性时具有鲁棒性，但未针对不确定性进行深入分析研究。张龙等^[12]针对航空发动机模型的结构和参数的不确定性进行了新的e修正鲁棒自适应控制方法设计及应用研究，仿真结果符合控制要求。

综上所述，国内外通过自适应控制方法对系统的不确定跟踪控制问题研究取得了大量的成果^[13~17]，但以航空发动机多变量系统为对象，针对不确定性补偿跟踪问题的深入研究很少。因此，本文采用增广模型参考自适应控制结构实现基于LQR控制器的不确定跟踪控制补偿，从而实现多变量系统的渐进一致稳定控制。

2 问题的描述

考虑如下形式的一类MIMO非线性系统

$ \begin{array}{l} {{\mathit{\boldsymbol{\dot x}}}_p} = {\mathit{\boldsymbol{A}}_p}{\mathit{\boldsymbol{x}}_p} + {\mathit{\boldsymbol{B}}_p}\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}\left[ {\mathit{\boldsymbol{u}} + \mathit{\boldsymbol{w}}\left( {{\mathit{\boldsymbol{x}}_p}} \right)} \right] + \mathit{\boldsymbol{\eta }}\left( t \right)\\ \mathit{\boldsymbol{y}} = {\mathit{\boldsymbol{C}}_p}{\mathit{\boldsymbol{x}}_p} + {\mathit{\boldsymbol{D}}_p}\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}\left[ {\mathit{\boldsymbol{u}} + \mathit{\boldsymbol{w}}\left( {{\mathit{\boldsymbol{x}}_p}} \right)} \right] \end{array} $

(1)

式中x_p ∈ R^n_p 是系统状态矢量，u ∈ R^m是控制输入，y ∈ R^m为系统输出。假设A_p ∈ R^{n_p × n_p}，B_p ∈ R^{n_p × m}，C_p ∈ R_{m × np}和D_p ∈ R^{m × m}是已知的常数矩阵。而Λ ∈ R^{m × m}为常量不确定性矩阵，它是一个具有严格正对角元素λ_i的未知对角矩阵。设(A_p, B_pΛ)可控。η(t) ≤ η_max(t)为有界未知干扰。w(x_p)代表系统匹配不确定性，假设可以表示为N个已知基本函数和未知常系数的线性组合，即

$ \mathit{\boldsymbol{w}}\left( {{\mathit{\boldsymbol{x}}_p}} \right) = {\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {{\mathit{\boldsymbol{x}}_p}} \right) $

(2)

式中Θ ∈ R^{m × n_p}是理想参数的未知矩阵，Φ(x_p) ∈ R^N代表已知的局部Lipschitz连续的回归矢量。

设r_cmd(t) ∈ R^m表示系统输出y ∈ R^m要跟随的有界指令。定义输出跟踪误差

$ {\mathit{\boldsymbol{e}}_y}\left( t \right) = \mathit{\boldsymbol{y}}\left( t \right) - {\mathit{\boldsymbol{r}}_{{\rm{cmd}}}}\left( t \right) $

(3)

定义e_yI为e_y的积分，则有

$ {{\mathit{\boldsymbol{\dot e}}}_{y{\rm{I}}}}\left( t \right) = {\mathit{\boldsymbol{e}}_y} = \mathit{\boldsymbol{y}} - {\mathit{\boldsymbol{r}}_{{\rm{cmd}}}} $

(4)

其增广开环动态可写为

$ \underbrace {\left[ \begin{array}{l} {{\mathit{\boldsymbol{\dot e}}}_{y{\rm{I}}}}\\ {{\mathit{\boldsymbol{\dot x}}}_p} \end{array} \right]}_{\mathit{\boldsymbol{\dot x}}} = \underbrace {\left[ {\begin{array}{*{20}{c}} 0&{{\mathit{\boldsymbol{C}}_p}}\\ 0&{{\mathit{\boldsymbol{A}}_p}} \end{array}} \right]}_\mathit{\boldsymbol{A}}\underbrace {\left[ \begin{array}{l} {\mathit{\boldsymbol{e}}_{y{\rm{I}}}}\\ {\mathit{\boldsymbol{x}}_p} \end{array} \right]}_\mathit{\boldsymbol{x}} + \underbrace {\left[ \begin{array}{l} {\mathit{\boldsymbol{D}}_p}\\ {\mathit{\boldsymbol{B}}_p} \end{array} \right]}_\mathit{\boldsymbol{B}}\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}\left[ {\mathit{\boldsymbol{u}} + \mathit{\boldsymbol{w}}\left( {{\mathit{\boldsymbol{x}}_p}} \right)} \right] + \\ \underbrace {\left[ \begin{array}{l} - {\mathit{\boldsymbol{I}}_{m \times m}}\\ {\mathit{\boldsymbol{0}}_{{n_p} \times m}} \end{array} \right]}_{{\mathit{\boldsymbol{B}}_{{\rm{ref}}}}}{\mathit{\boldsymbol{r}}_{{\rm{cmd}}}} + \underbrace {\left[ \begin{array}{l} \mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{\eta }}\left( t \right) \end{array} \right]}_{\mathit{\boldsymbol{\xi }}\left( t \right)} $

(5)

即

$ \mathit{\boldsymbol{\dot x}} = \mathit{\boldsymbol{A}}x + \mathit{\boldsymbol{B \boldsymbol{\varLambda} }}\left[ {\mathit{\boldsymbol{u}} + {\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {{\mathit{\boldsymbol{x}}_p}} \right)} \right] + {\mathit{\boldsymbol{B}}_{{\rm{ref}}}}{\mathit{\boldsymbol{r}}_{{\rm{cmd}}}} + \mathit{\boldsymbol{\xi }}\left( t \right) $

(6)

根据式（6）、式（2）、式（1）的调节输出y可写为

$ \mathit{\boldsymbol{y}} = \underbrace {\left[ {\begin{array}{*{20}{c}} 0&{{\mathit{\boldsymbol{C}}_p}} \end{array}} \right]}_\mathit{\boldsymbol{C}}\underbrace {\left[ \begin{array}{l} {\mathit{\boldsymbol{e}}_{y{\rm{I}}}}\\ {\mathit{\boldsymbol{x}}_p} \end{array} \right]}_\mathit{\boldsymbol{x}} + \underbrace {{\mathit{\boldsymbol{D}}_p}\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}}_\mathit{\boldsymbol{D}}\left[ {\mathit{\boldsymbol{u}} + {\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( \mathit{\boldsymbol{x}} \right)} \right] = \mathit{\boldsymbol{Cx}} + \mathit{\boldsymbol{D}}\left[ {\mathit{\boldsymbol{u}} + {\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( \mathit{\boldsymbol{x}} \right)} \right] $

(7)

因此，控制问题的目标是在存在系统常量参数不确定性Λ和Θ以及有界未知干扰ξ(t)的条件下实现有界跟踪。即设计控制输入u，使系统输出y以有界误差跟踪任意有界时变指令r_cmd，同时，系统中其他信号有界。

3 LQR基本控制器设计

令式（6）中Λ = I_{m × m}，Θ = 0_{N × m}，ξ(t) = 0得到线性标称开环动态，即

$ \begin{array}{l} \mathit{\boldsymbol{\dot x}} = \mathit{\boldsymbol{Ax}} + \mathit{\boldsymbol{Bu}} + {\mathit{\boldsymbol{B}}_{{\rm{ref}}}}{\mathit{\boldsymbol{r}}_{{\rm{cmd}}}}\\ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{Cx}} + \mathit{\boldsymbol{Du}} \end{array} $

(8)

设r_cmd为常量指令，利用具有比例+积分（PI）反馈连接LQR方法^{[14, 18]}，设计LQR最优控制律。问题转化为式（9）~式（10）的LQR问题。

$ \begin{array}{l} \mathit{\boldsymbol{\dot z}} = \mathit{\boldsymbol{Az}} + \mathit{\boldsymbol{Bv}}\\ \mathit{\boldsymbol{J}}\left( \mathit{\boldsymbol{v}} \right) = \int_0^\infty {\left( {{\mathit{\boldsymbol{z}}^{\rm{T}}}\mathit{\boldsymbol{Qz}} + {\mathit{\boldsymbol{v}}^{\rm{T}}}\mathit{\boldsymbol{Rv}}} \right){\rm{d}}t} \end{array} $

(9)

式中 $\mathit{\boldsymbol{z}} = \mathit{\boldsymbol{\dot x}} = {\left[ {{{\mathit{\boldsymbol{\dot e}}}_{y{\rm{I}}}}\; \; {{\mathit{\boldsymbol{\dot x}}}_{\rm{p}}}} \right]^{\rm{T}}}$ ， $\mathit{\boldsymbol{v}} = \mathit{\boldsymbol{\dot u}}$ ，Q = Q^T ≥ 0，R = R^T ≥ 0，控制输入v使得式（9）中J(v)最小化。

反馈形式的最优LQR解为

$ \mathit{\boldsymbol{v}} = \mathit{\boldsymbol{\dot u}} = - {\mathit{\boldsymbol{R}}^{ - 1}}{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{Pz}} = - \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{K}}_{\rm{I}}}}&{{\mathit{\boldsymbol{K}}_p}} \end{array}} \right]\left[ \begin{array}{l} {{\mathit{\boldsymbol{\dot e}}}_{y{\rm{I}}}}\\ {{\mathit{\boldsymbol{\dot x}}}_p} \end{array} \right] $

(10)

式中P是对应矩阵黎卡提代数方程唯一的对称正定解，即

$ {\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}} + \mathit{\boldsymbol{PA}} + \mathit{\boldsymbol{Q}} - \mathit{\boldsymbol{PB}}{\mathit{\boldsymbol{R}}^{ - 1}}{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{P}} = 0 $

(11)

对式（10）积分得基本LQR-PI控制器

$ {\mathit{\boldsymbol{u}}_{bI}} = - \mathit{\boldsymbol{K}}_x^{\rm{T}}\mathit{\boldsymbol{x}} = - {\mathit{\boldsymbol{K}}_I}{\mathit{\boldsymbol{e}}_{yI}} - {\mathit{\boldsymbol{K}}_p}{\mathit{\boldsymbol{x}}_p} = {\mathit{\boldsymbol{K}}_I}\frac{{\left( {{\mathit{\boldsymbol{y}}_{{\rm{cmd}}}} - \mathit{\boldsymbol{y}}} \right)}}{\mathit{\boldsymbol{s}}} - {\mathit{\boldsymbol{K}}_p}{\mathit{\boldsymbol{x}}_p} $

(12)

其中最优控制增益矩阵

$ \mathit{\boldsymbol{K}}_x^{\rm{T}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{K}}_{\rm{I}}}}&{{\mathit{\boldsymbol{K}}_{\rm{p}}}} \end{array}} \right] $

(13)

分块为积分增益K_I和比例增益K_p。相应的基本LQR-PI控制框图如图 1所示。

4 控制器MRAC增广

在存在系统不确定性Λ和Θ以及有界未知干扰ξ(t)时，跟踪性能变差。为恢复系统预期跟踪性能，对上述控制系统进行自适应增广^[19~21]。这个过程包括参考模型的选取和自适应律的设计。

4.1 参考模型选取

选取参考模型以表示期望闭环系统动态，它可以由将基本控制器式（12）带入线性系统式（8）得到。参考模型动态变为

$ \begin{array}{l} {{\mathit{\boldsymbol{\dot x}}}_{{\rm{ref}}}} = {\mathit{\boldsymbol{A}}_{{\rm{ref}}}}{\mathit{\boldsymbol{x}}_{{\rm{ref}}}} + {\mathit{\boldsymbol{B}}_{{\rm{ref}}}}{\mathit{\boldsymbol{r}}_{{\rm{cmd}}}}\\ {\mathit{\boldsymbol{y}}_{{\rm{ref}}}} = {\mathit{\boldsymbol{C}}_{{\rm{ref}}}}{\mathit{\boldsymbol{x}}_{{\rm{ref}}}} \end{array} $

(14)

式中A_ref = A-ABK_x^T，AC_ref = C-DK_x^T，A_ref为Hurwitz矩阵。

4.2 自适应律设计

增广后可将总控制输入分解为基本控制器式（12）和自适应增广u_ad（待构建），即

$ \mathit{\boldsymbol{u}} = - \mathit{\boldsymbol{K}}_x^{\rm{T}}\mathit{\boldsymbol{x}} + {\mathit{\boldsymbol{u}}_{{\rm{ad}}}} = {\mathit{\boldsymbol{u}}_{{\rm{bl}}}} + {\mathit{\boldsymbol{u}}_{{\rm{ad}}}} $

(15)

将式（15）带入原系统动态式（6）得

$ \begin{array}{l} \mathit{\boldsymbol{\dot x}} = {\mathit{\boldsymbol{A}}_{{\rm{ref}}}}\mathit{\boldsymbol{x}} + \mathit{\boldsymbol{B \boldsymbol{\varLambda} }}\left[ {{\mathit{\boldsymbol{u}}_{{\rm{ad}}}} + \underbrace {\overbrace {\left( {{\mathit{\boldsymbol{I}}_{m \times m}} - {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}^{ - 1}}} \right)}^{\mathit{\boldsymbol{K}}_u^{\rm{T}}}{\mathit{\boldsymbol{u}}_{{\rm{bl}}}} + {\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {{\mathit{\boldsymbol{x}}_p}} \right)}_{{{\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\bar \varPhi} }}\left( {{\mathit{\boldsymbol{u}}_{{\rm{bl}}}},{\mathit{\boldsymbol{x}}_p}} \right)}} \right] + \\ \;\;\;\;\;{\mathit{\boldsymbol{B}}_{{\rm{ref}}}}{\mathit{\boldsymbol{y}}_{{\rm{cmd}}}} + \mathit{\boldsymbol{\xi }}\left( t \right) \end{array} $

(16)

$ \mathit{\boldsymbol{y}} = {\mathit{\boldsymbol{C}}_{{\rm{ref}}}}\mathit{\boldsymbol{x}} + \mathit{\boldsymbol{D \boldsymbol{\varLambda} }}\left[ {{\mathit{\boldsymbol{u}}_{{\rm{ad}}}} + {{\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\bar \varPhi} }}\left( {{\mathit{\boldsymbol{u}}_{{\rm{bl}}}},{\mathit{\boldsymbol{x}}_p}} \right)} \right] $

定义回归矢量为

$ \mathit{\boldsymbol{ \boldsymbol{\bar \varPhi} }}\left( {{\mathit{\boldsymbol{u}}_{{\rm{bl}}}},{\mathit{\boldsymbol{x}}_p}} \right) = {\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{u}}_{{\rm{bl}}}^{\rm{T}}}&{{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}}\left( {{\mathit{\boldsymbol{x}}_p}} \right)} \end{array}} \right]^{\rm{T}}} $

(17)

未知/理想参数的增广矩阵为

$ \mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }} = {\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{K}}_u^{\rm{T}}}&{{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}^{\rm{T}}}} \end{array}} \right]^{\rm{T}}} $

(18)

选择自适应分量u_ad用于主导系统匹配不确定性 ${{\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\bar \varPhi} }}\left({{\mathit{\boldsymbol{u}}_{{\rm{b1}}}}, {\mathit{\boldsymbol{x}}_{\rm{p}}}} \right)$ ，即

$ \begin{array}{l} \mathit{\boldsymbol{\dot x}} = {\mathit{\boldsymbol{A}}_{{\rm{ref}}}}\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{B \boldsymbol{\varLambda} }}\Delta {{\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\bar \varPhi} }} + {\mathit{\boldsymbol{B}}_{{\rm{ref}}}}{\mathit{\boldsymbol{r}}_{{\rm{cmd}}}} + \mathit{\boldsymbol{\xi }}\left( t \right)\\ \mathit{\boldsymbol{y}} = {\mathit{\boldsymbol{C}}_{{\rm{ref}}}}\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{D \boldsymbol{\varLambda} }}\Delta {{\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\bar \varPhi} }} \end{array} $

(19)

式中 $\Delta \mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }} = \mathop {\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}}\limits^ \wedge -\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}$ 是参数估计误差矩阵。

引入状态跟踪误差

$ \mathit{\boldsymbol{e}} = \mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_{{\rm{ref}}}} $

(20)

由增广开环系统动态式（19）减去参考系统动态式（14）计算跟踪误差动态，即

$ \mathit{\boldsymbol{\dot e}} = {\mathit{\boldsymbol{A}}_{{\rm{ref}}}}\mathit{\boldsymbol{e}} - \mathit{\boldsymbol{B \boldsymbol{\varLambda} }}\Delta {{\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\bar \varPhi} }} + \mathit{\boldsymbol{\xi }}\left( t \right) $

(21)

为设计MRAC律并同时实现误差动态的闭环稳定性，考虑径向无界的二次李雅普诺夫候选函数

$ V\left( {\mathit{\boldsymbol{e}},\Delta \mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}} \right) = {\mathit{\boldsymbol{e}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_{{\rm{ref}}}}\mathit{\boldsymbol{e}} + tr\left( {\Delta {{\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_{\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}}^{ - 1}\Delta \mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} \boldsymbol{\varLambda} }}} \right) $

(22)

式中 ${\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_{\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}}} = \mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_{\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}}^{\rm{T}} > 0$ 的元素代表自适应率，P_ref = P_ref^T > 0是下面代数李雅普诺夫方程唯一的对称正定解

$ \mathit{\boldsymbol{A}}_{{\rm{ref}}}^{\rm{T}}{\mathit{\boldsymbol{P}}_{{\rm{ref}}}} + {\mathit{\boldsymbol{P}}_{{\rm{ref}}}}{\mathit{\boldsymbol{A}}_{{\rm{ref}}}} = - {\mathit{\boldsymbol{Q}}_{{\rm{ref}}}} $

(23)

其中Q_ref = Q_ref^T > 0。沿式（21）轨迹对V求时间微分得

$ \begin{array}{l} \dot V\left( {\mathit{\boldsymbol{e}},\Delta \mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}} \right) = - {\mathit{\boldsymbol{e}}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_{{\rm{ref}}}}\mathit{\boldsymbol{e}} + 2{\mathit{\boldsymbol{e}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_{{\rm{ref}}}}\mathit{\boldsymbol{\xi }}\left( t \right) + \\ 2{\rm{tr}}\left\{ {\Delta {{\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}}^{\rm{T}}}\left[ {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_{\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}}^{ - 1}\mathit{\boldsymbol{ \boldsymbol{\dot {\hat {\bar \varTheta}} } }} - \mathit{\boldsymbol{ \boldsymbol{\bar \varPhi} }}{\mathit{\boldsymbol{e}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_{{\rm{ref}}}}\mathit{\boldsymbol{B}}} \right]\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}} \right\} \end{array} $

(24)

定义基于射影理论的自适应律

$ \mathit{\boldsymbol{ \boldsymbol{\dot {\hat {\bar \varTheta}} } }} = {\rm{Proj}}\left[ {\mathit{\boldsymbol{ \boldsymbol{\hat {\bar \varTheta}} }},{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_{\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}}}\mathit{\boldsymbol{ \boldsymbol{\bar \varPhi} }}\left( {{\mathit{\boldsymbol{u}}_{{\rm{bl}}}},{\mathit{\boldsymbol{x}}_p}} \right){\mathit{\boldsymbol{e}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_{{\rm{ref}}}}\mathit{\boldsymbol{B}}} \right] $

(25)

以保证自适应增益的一致有界性。其中，

$ {\rm{Proj}}\left( {\mathit{\boldsymbol{\theta }},\mathit{\boldsymbol{y}}} \right) = \left\{ \begin{array}{l} \mathit{\boldsymbol{y}} - \frac{{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}\nabla \mathit{\boldsymbol{f}}\left( \mathit{\boldsymbol{\theta }} \right){{\left( {\nabla \mathit{\boldsymbol{f}}\left( \mathit{\boldsymbol{\theta }} \right)} \right)}^{\rm{T}}}}}{{\left\| {\nabla \mathit{\boldsymbol{f}}\left( \mathit{\boldsymbol{\theta }} \right)} \right\|_\mathit{\Gamma }^2}}\mathit{\boldsymbol{y}}f\left( \theta \right),\left[ {\mathit{\boldsymbol{f}} > 0 \cap {\mathit{\boldsymbol{y}}^{\rm{T}}}\nabla \mathit{\boldsymbol{f}}\left( \mathit{\boldsymbol{\theta }} \right) > 0} \right]\\ \mathit{\boldsymbol{y}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{other}} \end{array} \right. $

式中Proj(θ, y)为射影算子，f为可微凸函数，|∇f||_Γ² =(∇f)^TΓ∇f为∇f的加权欧氏平方范数。因此

$ \begin{array}{l} \dot V\left( {\mathit{\boldsymbol{e}},\Delta \mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}} \right) \le - {\mathit{\boldsymbol{e}}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_{{\rm{ref}}}}\mathit{\boldsymbol{e}} + 2{\mathit{\boldsymbol{e}}^{\rm{T}}}\mathit{\boldsymbol{P\xi }}\left( t \right) \le {\lambda _{\min }}\left( {{\mathit{\boldsymbol{Q}}_{{\rm{ref}}}}} \right){\left\| \mathit{\boldsymbol{e}} \right\|^2} + \\ 2\left\| \mathit{\boldsymbol{e}} \right\|{\lambda _{\max }}\left( {{\mathit{\boldsymbol{P}}_{{\rm{ref}}}}} \right){\mathit{\boldsymbol{\xi }}_{\max }} = - {\lambda _{\min }}\left( {{\mathit{\boldsymbol{Q}}_{{\rm{ref}}}}} \right)\left\| \mathit{\boldsymbol{e}} \right\|\left( {\left\| \mathit{\boldsymbol{e}} \right\| - 2\frac{{{\lambda _{\max }}\left( {{\mathit{\boldsymbol{P}}_{{\rm{ref}}}}} \right){\mathit{\boldsymbol{\xi }}_{\max }}}}{{{\lambda _{\min }}\left( {{\mathit{\boldsymbol{Q}}_{{\rm{ref}}}}} \right)}}} \right) \end{array} $

(26)

定义紧集为

$ \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} = \left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}} \right) \in {R^n} \times {R^{\left( {m + N} \right) \times m}} $

式中 $\left\| \mathit{\boldsymbol{e}} \right\| \le 2\frac{{{\mathit{\boldsymbol{\lambda }}_{\max }}\left({{\mathit{\boldsymbol{P}}_{{\rm{ref}}}}} \right)}}{{{\mathit{\boldsymbol{\lambda }}_{\min }}\left({{\mathit{\boldsymbol{Q}}_{{\rm{ref}}}}} \right)}}{\mathit{\boldsymbol{\xi }}_{\max }}$ ， ${\left\| {\Delta \mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}} \right\|_F} \le \Delta {{\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}}_{\max }}$

式（26）在该紧集外部，有 $\mathit{\boldsymbol{\dot V}}\left({\mathit{\boldsymbol{e}}, \Delta \mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}} \right) < 0$ ，其中

$ \Delta {{\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}}_{\max }} = 2\underbrace {\left[ {\begin{array}{*{20}{c}} {\mathit{\bar \Theta }_1^{\max }}& \cdots &{\mathit{\bar \Theta }_{m + N}^{\max }} \end{array}} \right]}_{{{\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}}_{\max }}} = 2{{\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}}_{\max }} $

(27)

式中 $\mathit{\bar \Theta }_j^{\max }$ 为第j列 ${\widehat {\mathit{\bar \Theta }}_j}\left(t \right)$ 最大允许范围。这证明了闭环系统中所有信号的一致最终有界性。同时也证明了系统的调节输出y(t)能够以有界误差跟踪任意外部有界指令r_cmd(t)。

对于自适应参数矩阵 $\widehat {\mathit{\boldsymbol{ \boldsymbol{\bar \varTheta} }}} \in {R^{\left({N + m} \right) \times m}}$ 的第j列 ${\widehat {\mathit{\bar \Theta }}_j}$ ，引入射影公差 $\varepsilon _j^{\mathit{\bar \Theta }} > 0$ 并选择凸函数

$ {f_j} = f\left( {{{\mathit{\hat {\bar \Theta }}}_j}} \right) = \frac{{\left( {1 + \varepsilon _j^{\mathit{\bar \Theta }}} \right){{\left\| {{{\mathit{\hat {\bar \Theta }}}_j}} \right\|}^2} - {{\left( {\mathit{\bar \Theta }_j^{\max }} \right)}^2}}}{{\varepsilon _j^{\mathit{\bar \Theta }}{{\left( {\mathit{\bar \Theta }_j^{\max }} \right)}^2}}} $

(28)

对于每个j = 1, 2, …, m + N，定义两个凸集

$ \begin{array}{l} \mathit{\Omega }_0^j = \left\{ {{{\mathit{\hat {\bar \Theta }}}_j} \in {R^{\left( {m + N} \right) \times 1}}:f\left( {{{\mathit{\hat {\bar \Theta }}}_j}} \right) \le 0} \right\} = \\ \left\{ {{{\mathit{\hat {\bar \Theta }}}_j} \in {R^{\left( {m + N} \right) \times 1}}:\left\| {{{\mathit{\hat {\bar \Theta }}}_j}} \right\| \le \frac{{\mathit{\bar \Theta }_j^{\max }}}{{\sqrt {1 + \varepsilon _j^{\mathit{\bar \Theta }}} }}} \right\}\\ \mathit{\Omega }_1^j = \left\{ {{{\mathit{\hat {\bar \Theta }}}_j} \in {R^{\left( {m + N} \right) \times 1}}:f\left( {{{\mathit{\hat {\bar \Theta }}}_j}} \right) \le 1} \right\} = \\ \left\{ {{{\mathit{\hat {\bar \Theta }}}_j} \in {R^{\left( {m + N} \right) \times 1}}:\left\| {{{\mathit{\hat {\bar \Theta }}}_j}} \right\| \le \mathit{\bar \Theta }_j^{\max }} \right\} \end{array} $

(29)

第j个凸函数式（29）的梯度为

$ \nabla {f_j} = \frac{{\left( {1 + \varepsilon _j^{\mathit{\bar \Theta }}} \right)}}{{\varepsilon _j^{\mathit{\bar \Theta }}{{\left( {\mathit{\bar \Theta }_j^{\max }} \right)}^2}}}\nabla \left[ {{{\left\| {{{\mathit{\bar \Theta }}_j}} \right\|}^2}} \right] = \frac{{2\left( {1 + \varepsilon _j^{\mathit{\bar \Theta }}} \right)}}{{\varepsilon _j^{\mathit{\bar \Theta }}\mathit{\bar \Theta }_j^{\max }}}{{\mathit{\hat {\bar \Theta }}}_j} $

(30)

因此，式（25）的自适应律变为

$ {{\mathit{\dot {\hat {\bar \Theta} } }}_j} = {\mathit{\Gamma }_{\mathit{\bar \Theta }}}\left\{ \begin{array}{l} {\left( {\mathit{\bar \Phi }{\mathit{e}^{\rm{T}}}{P_{{\rm{ref}}}}B} \right)_j} - \frac{{\nabla {f_j}\nabla f_j^{\rm{T}}}}{{\left\| {\nabla {f_j}} \right\|_{{\mathit{\Gamma }_{\mathit{\bar \Theta }}}}^2}}{\mathit{\Gamma }_{\mathit{\bar \Theta }}}{\left( {\mathit{\bar \Phi }{\mathit{e}^{\rm{T}}}{P_{{\rm{ref}}}}B} \right)_j}{f_j}\\ \left[ {{f_i} > 0 \cap \left( {\mathit{\bar \Phi }{\mathit{e}^{\rm{T}}}{P_{{\rm{ref}}}}B} \right)_j^{\rm{T}}{\mathit{\Gamma }_{\mathit{\bar \Theta }}}\nabla {f_j}} \right]\\ {\left( {\mathit{\bar \Phi }{\mathit{e}^{\rm{T}}}{P_{{\rm{ref}}}}B} \right)_j},\;\;\;\;\;{\rm{other}} \end{array} \right. $

(31)

上述方法，确保了自适应过程时变矩阵 ${\widehat {\mathit{\bar \Theta }}_j}\left(t \right)$ 的一致有界性。即

$ \begin{array}{l} \left\{ {\left\| {{{\mathit{\hat {\bar \Theta }}}_j}\left( 0 \right)} \right\| \le \frac{{\mathit{\bar \Theta }_j^{\max }}}{{\sqrt {1 + \varepsilon _j^{\mathit{\bar \Theta }}} }}} \right\} \Rightarrow \\ \left\{ {\left\| {{{\mathit{\hat {\bar \Theta }}}_j}\left( t \right)} \right\| \le \mathit{\bar \Theta }_j^{\max },\forall t \ge 0,1 \le j \le m + N} \right\} \end{array} $

(32)

自适应增广分量可以写为

$ {\mathit{\boldsymbol{u}}_{{\rm{ad}}}} = - {{\mathit{\boldsymbol{ \boldsymbol{\hat {\bar \varTheta}} }}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\bar \varPhi} }}\left( {{\mathit{\boldsymbol{u}}_{{\rm{bl}}}},{\mathit{\boldsymbol{x}}_p}} \right) $

(33)

则总控制输入（LQR PI基本控制器+自适应）为

$ \mathit{\boldsymbol{u}} = {\mathit{\boldsymbol{u}}_{{\rm{bl}}}} + {\mathit{\boldsymbol{u}}_{{\rm{ad}}}} = \underbrace { - \mathit{\boldsymbol{K}}_x^{\rm{T}}\mathit{\boldsymbol{x}}}_{{\mathit{\boldsymbol{u}}_{{\rm{bl}}}}} + \underbrace {\left[ {\mathit{\boldsymbol{\hat K}}_u^{\rm{T}}{\mathit{\boldsymbol{u}}_{{\rm{bl}}}} - {{\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {{\mathit{\boldsymbol{x}}_p}} \right)} \right]}_{{\mathit{\boldsymbol{u}}_{{\rm{ad}}}}} $

(34)

即

$ \begin{array}{l} \mathit{\boldsymbol{u}} = \left( {{\mathit{\boldsymbol{I}}_{m \times m}} - \mathit{\boldsymbol{\hat K}}_u^{\rm{T}}} \right){\mathit{\boldsymbol{u}}_{{\rm{bl}}}} - {{\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {{\mathit{\boldsymbol{x}}_p}} \right) = \\ - \left( {{\mathit{\boldsymbol{I}}_{m \times m}} - \mathit{\boldsymbol{\hat K}}_u^{\rm{T}}} \right)\mathit{\boldsymbol{K}}_x^{\rm{T}}\mathit{\boldsymbol{x}} - {{\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {{\mathit{\boldsymbol{x}}_p}} \right) = \\ \left( {{\mathit{\boldsymbol{I}}_{m \times m}} - \mathit{\boldsymbol{\hat K}}_u^{\rm{T}}} \right)\left( {{\mathit{\boldsymbol{K}}_{\rm{I}}}\frac{{{\mathit{\boldsymbol{r}}_{{\rm{cmd}}}} - \mathit{\boldsymbol{y}}}}{\mathit{\boldsymbol{s}}} - {\mathit{\boldsymbol{K}}_{\rm{p}}}{\mathit{\boldsymbol{x}}_p}} \right){{\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {{\mathit{\boldsymbol{x}}_p}} \right) \end{array} $

(35)

整体控制框图如图 2所示。图中虚线表示增益 ${{\mathit{\boldsymbol{\hat K}}}_u}$ 和 ${\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}}$ 的自适应性，其动态由跟踪误差e = x-x_ref驱动，其中x是增广系统式（5）的状态。

对于每个设计，尽管存在匹配不确定性和有界未知干扰，控制器依然迫使系统输出y渐近跟踪参考模型输出r_cmd。

5 仿真验证

航空发动机的动态可以用一个二阶微分方程逼近

$ \begin{array}{l} \underbrace {\left[ \begin{array}{l} {{\dot n}_{\rm{H}}}\\ {{\dot n}_{\rm{L}}} \end{array} \right]}_{{{\mathit{\boldsymbol{\dot x}}}_p}} = \underbrace {\left[ {\begin{array}{*{20}{c}} { - 2.148}&{ - 0.403}\\ {2.862}&{ - 3.017} \end{array}} \right]}_{{\mathit{\boldsymbol{A}}_p}}\underbrace {\left[ \begin{array}{l} {n_{\rm{H}}}\\ {n_{\rm{L}}} \end{array} \right]}_{{\mathit{\boldsymbol{x}}_p}} + \\ \underbrace {\left[ {\begin{array}{*{20}{c}} {0.298}&{0.594}\\ {0.285}&{1.708} \end{array}} \right]}_{{\mathit{\boldsymbol{B}}_p}}\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}\left( {\underbrace {\left[ \begin{array}{l} {q_{{\rm{m}},{\rm{f}}}}\left( t \right)\\ {A_8}\left( t \right) \end{array} \right]}_\mathit{\boldsymbol{u}} + \mathit{\boldsymbol{w}}\left( {{\mathit{\boldsymbol{x}}_p}} \right)} \right) + \mathit{\boldsymbol{\eta }}\left( t \right) \end{array} $

(36)

式中燃油供油量q_{m, f}(t)和尾喷口面积A₈(t)为控制量，是系统的控制输入，为了便于仿真验证，这里假设w(x_p)是状态高压转子转速n_H和低压转子转速n_L的线性组合。定义低压转子转速n_L和涡轮前温度T₅为输出变量

$ \mathit{\boldsymbol{y}}\left( t \right) = \left[ \begin{array}{l} {n_{\rm{L}}}\left( t \right)\\ {T_5}\left( t \right) \end{array} \right] = \underbrace {\left[ {\begin{array}{*{20}{c}} 0&1\\ { - 0.013}&{ - 0.364} \end{array}} \right]}_{{\mathit{\boldsymbol{C}}_p}}{\mathit{\boldsymbol{x}}_p} + \left[ {\begin{array}{*{20}{c}} 0&0\\ {0.365}&{ - 0.232} \end{array}} \right]\mathit{\boldsymbol{u}} $

(37)

假设Λ = diag(1, 1), Θ = 0_{2 × 2}, η(t)= 0利用式（5）对其增广得增广开环动态为

$ \begin{array}{l} \underbrace {\left[ {\begin{array}{*{20}{c}} {{{\dot e}_{{n_{\rm{H}}}I}}}\\ {{{\dot e}_{{n_{\rm{L}}}I}}}\\ {{{\mathit{\boldsymbol{\dot x}}}_p}} \end{array}} \right]}_x = \underbrace {\left[ {\begin{array}{*{20}{c}} {{0_{2 \times 2}}}&{{\mathit{\boldsymbol{C}}_p}}\\ {{0_{2 \times 2}}}&{{\mathit{\boldsymbol{A}}_p}} \end{array}} \right]}_\mathit{\boldsymbol{A}}\underbrace {\left[ {\begin{array}{*{20}{c}} {{e_{{n_{\rm{H}}}I}}}\\ {{e_{{n_{\rm{L}}}I}}}\\ {{\mathit{\boldsymbol{x}}_p}} \end{array}} \right]}_\mathit{\boldsymbol{x}} + \underbrace {\left[ \begin{array}{l} {0_{2 \times 2}}\\ {\mathit{\boldsymbol{B}}_p} \end{array} \right]}_\mathit{\boldsymbol{B}}\underbrace {\left[ \begin{array}{l} {q_{{\rm{m}},{\rm{f}}}}\left( t \right)\\ {A_8}\left( t \right) \end{array} \right]}_\mathit{\boldsymbol{u}} + \underbrace {\left[ \begin{array}{l} - {I_{2 \times 2}}\\ {0_{2 \times 2}} \end{array} \right]}_{{\mathit{\boldsymbol{B}}_{{\rm{ref}}}}}\underbrace {\left[ \begin{array}{l} {n_{{\rm{H}},{\rm{cmd}}}}\\ {n_{{\rm{L}},{\rm{cmd}}}} \end{array} \right]}_{{\mathit{\boldsymbol{r}}_{{\rm{cmd}}}}}\\ \;\;\;\;\;\mathit{\boldsymbol{y}} = \underbrace {\left[ {\begin{array}{*{20}{c}} {{0_{2 \times 2}}}&{{\mathit{\boldsymbol{C}}_p}} \end{array}} \right]}_\mathit{\boldsymbol{C}}x + {\mathit{\boldsymbol{D}}_p}\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}\left[ {\mathit{\boldsymbol{u}} + \mathit{\boldsymbol{w}}\left( \mathit{\boldsymbol{x}} \right)} \right] \end{array} $

(38)

式（9）中，Q，R取如下值

$ \mathit{\boldsymbol{Q}} = \left[ {\begin{array}{*{20}{c}} {100}&0&0&0\\ 0&{150}&0&0\\ 0&0&0&0\\ 0&0&0&0 \end{array}} \right],\mathit{\boldsymbol{R = }}\left[ {\begin{array}{*{20}{c}} 1&0\\ 0&1 \end{array}} \right] $

(39)

求得LQR基本控制器为

$ \begin{array}{l} {\mathit{\boldsymbol{u}}_{{\rm{b1}}}} = - \left[ {\begin{array}{*{20}{c}} {4.7137}&{10.8015}\\ {8.8194}&{ - 5.7731} \end{array}} \right]\int {\left[ {\begin{array}{*{20}{c}} {{n_{{\rm{L}},{\rm{cmd}}}} - {n_{\rm{L}}}}\\ {{T_{5,{\rm{cmd}}}} - {T_5}} \end{array}} \right]{\rm{d}}t} - \\ \;\;\;\;\;\;\;\left[ {\begin{array}{*{20}{c}} {0.4652}&{0.5718}\\ {0.9590}&{1.9212} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{W_{\rm{f}}}}\\ {{A_8}} \end{array}} \right] \end{array} $

(40)

当Λ = 0.95， $\mathit{\boldsymbol{ \boldsymbol{\varTheta} = }}1.8 \times \left[ {\begin{array}{*{20}{c}} { -{A_p}\left({1, 1} \right)} & {{A_p}\left({1, 2} \right)}\\ { -{A_p}\left({2, 1} \right)} & { -{A_p}\left({2, 2} \right)} \end{array}} \right]$ 时，选择如表 1所示自适应控制参数。

分别对该系统的原LQR基本控制器和增广自适应控制补偿方法进行仿真分析，仿真结果如图 3~8所示。

由图 3可以看出，系统在没有不确定性时，原LQR控制器保证了系统无超调且调节时间t_s ≤ 1.5s，符合控制要求。图 4所示给出了系统存在不确定性时，原LQR控制器无法满足控制要求，系统超调量很大且调节时间t_s ≥ 5s，同时，系统在10s后出现大幅度振荡。

图 5~8给出了引入增广模型参考自适应方法的仿真结果。可以看出，补偿方法很好地实现了系统存在不确定性时的跟踪控制，控制误差小于0.25%，超调量小于0.5%，调节时间t_s ≤1.5s，且输入量主燃油和喷口喉部面积等变化均符合执行机构要求。

图 8给出了控制器增益K_x的2范数和自适应参数Θ的范数的变化情况。由图 8可以看出，所选取的自适应参数在控制周期内能够保持有界。综上所述，采用的增广模型参考自适应控制方法能够实现原LQR基本控制器对不确定性控制问题的优化。

图 9、图 10给出了增广模型参考自适应控制方法在发动机非线性模型上的仿真结果。在不同高度、马赫数、低压转子转速条件下，调整控制器自适应参数进行仿真验证。由图 9、图 10可以看出，由该方法所设计的控制器在发动机不同状态下控制效果均满足要求。限于篇幅原因，其他状态点仿真结果不一一示出。

6 结论

采用增广自适应控制方法对多变量系统的不确定性跟踪补偿问题进行了研究，得到了以下结论：

（1）系统存在不确定性时，采用增广模型参考自适应控制方法能够实现LQR基本控制器的跟踪控制优化，解决了单一LQR控制器在系统不确定性时引起的大幅度振荡现象。

（2）分别基于单一LQR控制方法和增广模型参考自适应控制补偿方法对航空发动机多变量系统进行了控制器设计与应用。结果表明：补偿方法对发动机系统存在不确定时的控制误差小于0.25%，超调量小于0.5%，且调节时间小于1.5s，控制效果符合航空发动机控制要求。