放宽静稳定度设计是通过配置重心与焦点的相对位置以降低配平阻力和改善飞机机动性能的一种主动控制技术^[1]。现代战机采用放宽静稳定度设计方案能够提高其最大可用过载、大迎角机动性及过失速机动能力等关键性能^[2]。

飞机设计初步阶段的操稳评估是影响整个布局方案与飞行控制安全的关键环节^[3]，其传统作法主要是依据基于飞机本体动力学方程的静态评估准则(包括静稳定导数C_mα、倍幅时间T₂、最大低头力矩等)进行逐项考察^[4]。但在工程实践中，一方面这些准则的边界往往来自于经验，易致使偏于保守的评估结果；另一方面，对于即便满足静态准则的布局设计，在依据飞行品质的要求设计控制增稳系统后，仍可能出现因电传飞控系统整体时延带来的复杂动力学特性^[5](如分支与极限环等)所导致的控制震荡乃至发散不可控的情况。因此，一方面，在给定战机布局条件下，对其闭环控制最大允许时间延迟边界的确定，将有助于指导电传飞控系统的后续详细设计；而另一方面，随着电传飞控技术的持续发展，其固有时间延迟也在逐步减小，这对设计更为先进、本体更为静不稳定的飞机布局是有利的，在给定飞控系统最小可达时延能力时，确定其所能支撑的最大静不稳定布局边界，亦将有助于初步设计阶段布局方案的高效迭代。

针对战机纵向静稳定度所对应飞控系统时间延迟允许边界的确定问题，以纵向短周期方程为例，首先分析了静稳定度与短周期方程中各参数的关系，随后基于飞行品质要求设计纵向增稳控制律，同时引入等效输入时间延迟建立纵向闭环特征方程，利用根轨迹趋势理论并结合数值计算方法求解了不同静稳定度下系统的时间延迟稳定边界，最后分析了短周期方程中的动导不确定性及舵效不确定性对该边界的影响。该时间延迟边界求解方法，同样适用于对战机横航向的最大允许时间延迟的确定，这对于在初步设计阶段的飞机布局设计和飞控系统时间延迟指标确定，具有一定的工程指导意义。

1 战斗机电传飞控系统的时间延迟

战斗机电传飞控系统的主要构成如图 1所示，本节将描述在工程实践中各子环节具有的时延情况。

1) 传感器。用于飞行控制系统的信号主要包括大气数据(气流角、动静压等)、惯性运动数据(三轴速率、三轴过载等)，当前所采用的嵌入式大气数据传感系统^[6]和光学陀螺惯导方案，相比于原机械风标和机械陀螺，信号采集的稳定性和准确性都大幅提升，带来的信号测量时延极小，一般不超过2个测量周期。

2) 作动器。目前作动器的驱动方式较多，包括传统的液压、电动以及日渐成熟的电动静液作动器(Electro-Hydrostatic Actuator, EHA)^[7]等，从目前的工程实践来看，这一环节的时间延迟是显著的，一般来说在50~100 ms量级。

3) 总线传输。飞控系统总线传输的时间延迟与飞行控制周期相关。参照MIL-STD-1553标准^[8]，以15 ms周期为例，其对应的时延最大约22.5 ms，随着未来飞控系统总线传输技术向光传方案方向的发展^[9]，这部分时延将进一步降低。

4) 飞行员操纵时延。飞行员从接收到状态信息反馈到做出操纵的时间延迟主要来自飞行员的反应时间，一般在200 ms左右^[10]，加之飞行员指令输入单元(Pilot Input Unit, PIU)的40 ms左右时延，总时延约在240 ms量级。尽管飞行员操纵时延较大，但其输入是控制增稳系统(Control Augmentation System, CAS)的外部执行指令，因此在考虑CAS稳定性的时候，并不将该时延纳入系统时延。这一部分操纵时延主要涉及的问题是飞行员诱发震荡(Pilot Induced Oscillations, PIO)^[11]现象。

总体来说，主要考虑飞控系统中传感器、作动器及总线传输中存在的时间延迟，当前飞控系统总体时延大约在120 ms，将统一考虑为等效输入时延来探讨对闭环系统的影响。

2 静稳定度与短周期方程关系

以纵向短周期模态为例，分析静稳定度与短周期方程参数间的关系。静稳定度定义为

$ {K_n} = {h_n} - h $	（1）

式中：h_n为中性点位置；h为重心位置。在低中速飞行的常规迎角区域，h_n基本不变^[12]。式(1)表明，静稳定度与重心位置基本线性相关，在飞机构型不变的情况下，通过配置重心位置即可实现对其静稳定度的设计。纵向静稳定系数为

$ {C_{{m_\alpha }}} = {C_{{L_\alpha }}}(h - {h_n}) = - {C_{{L_\alpha }}}{K_n} $	（2）

式中：C_Lα为升力迎角系数。纵向短周期方程可以表达为^[13]

$ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} {\Delta \dot \alpha }\\ {\Delta \dot q} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - {Z_\alpha }}&1\\ {{{\bar M}_\alpha } - {{\bar M}_{\dot \alpha }}{Z_\alpha }}&{{{\bar M}_q} + {{\bar M}_{\dot \alpha }}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\Delta \alpha }\\ {\Delta q} \end{array}} \right] + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left[ {\begin{array}{*{20}{c}} { - {Z_{{\delta _e}}}}\\ {{{\bar M}_{{\delta _e}}} - {{\bar M}_{\dot \alpha }}{Z_{{\delta _e}}}} \end{array}} \right]\Delta {\delta _e} = \mathit{\boldsymbol{Ax}} + \mathit{\boldsymbol{Bu}} \end{array} $

（3）

式中:各参数具体表达式为

$ \left\{ \begin{array}{l} \begin{array}{*{20}{l}} {{Z_\alpha } = \frac{{({C_{{D^*}}} + {C_{{L_\alpha }}}){Q_*}{\kern 1pt} S}}{{m{V_*}}}}\\ {{{\bar M}_\alpha } = {C_{{m_\alpha }}}\frac{{{Q_*}Sc}}{{{I_y}}}} \end{array}\\ \begin{array}{*{20}{l}} {{{\bar M}_q} = {C_{{m_q}}}\left( {\frac{c}{{2{V_*}}}} \right)\frac{{{Q_*}{\kern 1pt} Sc}}{{{I_y}}}}\\ {{{\bar M}_{\dot \alpha }} = {C_{{m_{\dot \alpha }}}}\left( {\frac{c}{{2{V_*}}}} \right)\frac{{{Q_*}{\kern 1pt} Sc}}{{{I_y}}}} \end{array}\\ \begin{array}{*{20}{l}} {{Z_{{\delta _e}}} = {C_{{L_{{\delta _e}}}}}\frac{{{Q_*}{\kern 1pt} Sc}}{{m{V_*}}}}\\ {{{\bar M}_{{\delta _e}}} = {C_{{m_{{\delta _e}}}}}\frac{{{Q_*}{\kern 1pt} Sc}}{{{I_y}}}} \end{array} \end{array} \right. $

（4）

其中，C_D^*为配平阻力系数; Q_*为动压；S为机翼参考面积；m为飞机质量；V_*为空速；c为参考弦长；I_y为惯矩；C_mq为动导阻尼项; C_{m$\dot{\alpha }$}为洗流时差力矩导数; C_Lδe为舵面升力系数; C_mδe为舵面俯仰力矩系数。

以上各项系数中，配平阻力系数随着静稳定度的放宽而减小；舵面俯仰力矩系数因静稳定度放宽时导致舵面力臂缩短而减小；其余参数随静稳定度放宽时无明显变化。至此，分析了短周期方程中各参数随静稳定度K_n的变化关系。

以F-16飞机为对象^[14]，通过迎角α、油门位置c_τ及升降舵偏度δ_e对飞机进行配平，在状态点H=100 m, Ma=0.46处，得到飞机配平迎角及油门位置随静稳定度变化曲线如图 2所示。可以看出，随着静稳定度的放宽，飞机所需的配平迎角将减小，配平阻力将减小，相应的油门位置c_τ将减小。

短周期方程A、B矩阵中各参数随静稳定度变化曲线如图 3所示，A₁₁代表A矩阵中第一行第一列元素，其他类似。可以看到，随着静稳定度的放宽，A₂₁项将产生较大的变化，该变化主要来自于C_mα随静稳定度的变化；控制矩阵B中的B₂₁项也有明显变化，该变化主要来自于重心后移导致的舵面力臂减小从而引起舵面俯仰力矩系数C_mδe的降低。A、B阵中其他参数随静稳定度的变化不明显。

短周期方程开环根轨迹随静稳定度的变化曲线如图 4所示。可以看出，随着静稳定度的放宽，飞机在该处的开环根轨迹朝着右半平面的方向移动，其自然模态的频率和阻尼比逐步降低，即飞机的稳定性逐步降低。

以上研究分析了纵向短周期方程的开环特性随静稳定度的变化关系。为研究闭环系统的延迟稳定边界问题，依据飞行品质的要求，采用极点配置方法设计增稳反馈控制律^[15]，控制律设计为u=Kx，则得到其闭环系统方程为

$ \mathit{\boldsymbol{\dot x}} = \mathit{\boldsymbol{Ax}} + \mathit{\boldsymbol{BKu}} $	（5）

第1节中指出，现有飞控系统架构下，闭环系统的延迟主要考虑为等效输入延迟，引入输入时延τ后对式(5)作Laplace变换，闭环控制系统变为

$ \mathit{\boldsymbol{sx}} = \mathit{\boldsymbol{Ax}} + \mathit{\boldsymbol{BKx}}{\rm{exp}}( - \tau s) $	（6）

式中：s代表Laplace算子，其相应的特征方程为

$ {\rm{ det }}(s\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{A}} - \mathit{\boldsymbol{BK}}{\kern 1pt} {\rm{exp}}( - \tau s)) = 0 $	（7）

3 时间延迟稳定边界

针对式(7)所代表的特征方程，研究的思路为：通过求解出其特征根穿越虚轴时对应的时间延迟和穿越虚轴的方向，来研究闭环特征方程随时间延迟变化所引起的稳定性变化趋势。代入短周期方程式(3)展开后可以得到特征方程为

$ \begin{array}{*{20}{c}} {{s^2} - {\rm{tr}} (\mathit{\boldsymbol{A}})s - {\rm{tr}} (\mathit{\boldsymbol{BK}}) s{\rm{exp}} ( - \tau s) + }\\ {\eta {\rm{exp}}( - \tau s) + {\rm{det}} (\mathit{\boldsymbol{A}}) = 0} \end{array} $	（8）

式中：tr(·)代表矩阵的迹；det(A)为A的行列式；η=A₁₁BK₂₂+A₂₂BK₁₁－A₁₂BK₂₁－A₂₁BK₁₂。将闭环系统的特征方程重新表达为

$ \mathit{\boldsymbol{CE}}(s, \tau ) = \sum\limits_{l = 0}^n {{P_l}} ({s^\gamma }){\rm{exp}}( - l\tau s) = 0 $	（9）

式中：P_l(s^γ)表示关于s的多项式，l表示时延τ的阶次，文献[16]证明了式(9)在τ∈$\mathbb{R}$⁺范围内仅存在有限个穿越虚轴的特征根。穿越虚轴时，其特征根为纯虚根，即s=±jω_ck, j=－1，其中ω_ck∈$\mathbb{R}$⁺为穿越处的频率，k=0, 1, …, m代表有限特征根的个数。依据Rekasius变换^[17]，可以得到

$ \exp ( - {\tau _{ck}}s){|_{s = {\rm{j}}{\omega _{ck}}}} = \frac{{1 - {\rm{j}}{T_{ck}}}}{{1 + {\rm{j}}{T_{ck}}}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {T_{ck}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \in {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathbb{R} $	（10）

其中，

$ \left\{ {\begin{array}{*{20}{l}} {{T_{ck}} = {\rm{tan}}\left( {\frac{{{\tau _{ck}}{\omega _{ck}}}}{2}} \right)}\\ {{\tau _{ck}} = \frac{2}{{{\omega _{ck}}}}({\rm{ta}}{{\rm{n}}^{ - 1}}(T) + p\pi ){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} p = 0, 1, \cdots , \infty } \end{array}} \right. $

（11）

相应地，可以得到

$ {\rm{exp}}( - l{\tau _{ck}}s){|_{s = {\rm{j}}{\omega _{ck}}}} = {\left( {\frac{{1 - {\rm{j}}{T_{ck}}}}{{1 + {\rm{j}}{T_{ck}}}}} \right)^l}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {T_{ck}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \in {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathbb{R} $

（12）

将式(9)乘以(1+jT_ck)ⁿ，并代入式(12)可得

$ \sum\limits_{l = 0}^n {{P_l}} ({s^\gamma }){(1 - {\rm{j}}{T_{ck}})^l}{(1 + {\rm{j}}{T_{ck}})^{(n - l)}} $	（13）

由于虚轴上，P_l(s^γ)仅与穿越频率ω_ck相关，则式(13)的解仅与ω_ck、T_ck相关，按照虚部和实部分别整理，得到

$ h(s, {T_{ck}}){|_{s = {\rm{j}}{\omega _{ck}}}} = {h_R}({\omega _{ck}}, {T_{ck}}) + {\rm{j}}{h_{\rm{I}}}({\omega _{ck}}, {T_{ck}}) $	（14）

当式(14)为零成立时，其虚部和实部均为零，即

$ \left\{ {\begin{array}{*{20}{l}} {{h_{\rm{R}}} = \sum\limits_{i = 0}^n {{a_i}} ({\omega _{ck}})T_{ck}^i = 0}\\ {{h_{\rm{I}}} = \sum\limits_{i = 0}^n {{b_i}} ({\omega _{ck}})T_{ck}^i = 0} \end{array}} \right. $	（15）

根据文献[18]，式(15)成立的必要条件为其对应的sylvester矩阵的行列式为零，即式(20)成立。当式(20)成立时，以下条件中至少有一项成立：①存在相应的(ω_ck, T_ck)，使得式(15)成立；② h_R、h_I的最高次系数a_n(ω)=b_n(ω)=0；③ h_R的系数a_i(ω_ck)=0；④ h_I的系数b_i(ω_ck)=0。

若式(20)无解，则说明在所有的时间延迟条件下，时延闭环系统均不存在稳定性的变化，即保持稳定或不稳定的状态。若式(20)存在解，通过排除与检查条件②③④则可以在条件①成立的前提下，求出所有的ω_ck及与其相对应的T_ck的值，且ω_ck与T_ck为一一对应关系，对应的时间延迟量通过式(11)给出。

至此，得到所有虚轴穿越处的频率ω_ck及其对应的T_ck，但此处对应的时间延迟τ_ck为周期解，即其个数为无限个。判断τ_ck对应的虚轴穿越引起的系统稳定性变化，将是求取时延稳定边界需要解决的问题。

方程特征根在虚轴穿越处对时延τ的导数定义为

$ S_\tau ^s{|_{s = {\rm{j}}{\omega _{ck}}, \tau = {\tau _{ck}}}} = {\left. {\frac{{{\rm{d}}s}}{{{\rm{d}}\tau }}} \right|_{s = {\rm{j}}{\omega _{ck}}, \tau = {\tau _{ck}}}} $	（16）

根据隐函数的求导原则，并结合特征方程式(9)可以得到

$ \begin{array}{l} S_\tau ^s{|_{s = {\rm{j}}{\omega _{ck}}, \tau = {\tau _{ck}}}} = {\left. { - \frac{{\frac{{\partial C}}{{\partial \tau }}}}{{\frac{{\partial C}}{{\partial s}}}}} \right|_{s = {\rm{j}}{\omega _{ck}}, \tau = {\tau _{ck}}}} = \\ {\left. {\frac{{\sum\limits_{l = 0}^n {{P_l}} ({s^\gamma })ls{\rm{exp}}( - l\tau s)}}{{\sum\limits_{l = 0}^n {\left| {\frac{{{\rm{d}}{P_l}({s^\gamma })}}{{{\rm{d}}s}} - l{P_l}({s^\alpha })\tau } \right|{\rm{exp}}( - {\rm{j}}\tau s)} }}} \right|_{s = {\rm{j}}{\omega _{ck}}, \tau = {\tau _{ck}}}} \end{array} $

（17）

对稳定性的分析，关心的是根轨迹穿越是从$\mathbb{C}$⁺→$\mathbb{C}$^－，即从右半平面到左半平面的稳定穿越；或从$\mathbb{C}$^－→$\mathbb{C}$⁺的不稳定穿越。因此研究τ_ck变化时，方程特征根s的实部的变化趋势，定义根轨迹趋势(Root Tendency, RT)为^[19]

$ {\rm{RT}}{|_{s = {\rm{j}}{\omega _{ck}}, \tau = {\tau _{ck}}}} = {\rm{sgn}} \left[ {{\rm{Re}} \left( {{{\left. {\frac{{{\rm{d}}s}}{{{\rm{d}}\tau }}} \right|}_{s = {\rm{j}}{\omega _{ck}}, \tau = {\tau _{ck}}}}} \right)} \right] $	（18）

结合式(17)得到

$ \begin{array}{l} {\rm{RT}}{|_{s = {\rm{j}}{\omega _{ck}}, \tau = {\tau _{ck}}}} = {\rm{sgn}} \left[ {{\rm{Re}} \left( {{{\left. {\frac{{{\rm{d}}s}}{{{\rm{d}}\tau }}} \right|}_{s = {\rm{j}}{\omega _{ck}}, \tau = {\tau _{ck}}}}} \right)} \right] = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{sgn}} \left[ {{\rm{Re}} \left( {{{\left. {\frac{{\sum\limits_{l = 0}^n {{P_l}} ({s^\gamma })s{\rm{exp}}( - l\tau s)}}{{\sum\limits_{l = 0}^n {\left| {\frac{{{\rm{d}}{P_l}({s^\gamma })}}{{{\rm{d}}s}} - l{P_l}({s^\gamma })\tau } \right|{\rm{exp}}( - l\tau s)} }}} \right|}_{s = {\rm{j}}{\omega _{ck}}, \tau = {\tau _{ck}}}}} \right)} \right] = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{sgn}} \left[ {{\rm{Im}} \left( {{{\left. {\frac{{\sum\limits_{l = 0}^n {\frac{{{\rm{d}}{P_l}({s^\gamma })}}{{{\rm{d}}s}}} {\rm{exp}}( - l\tau s)}}{{\sum\limits_{l = 0}^n {{P_l}({s^\gamma }){\rm{exp}}( - l\tau s)} }}} \right|}_{s = {\rm{j}}{\omega _{ck}}, \tau = {\tau _{ck}}}}} \right)} \right] \end{array} $

（19）

$ \begin{array}{l} {\rm{RT}} ({h_{\rm{R}}}, {h_{\rm{I}}}) = {\left| {\begin{array}{*{20}{c}} {{a_n}(\omega )}&{{a_{n - 1}}(\omega )}&{{a_{n - 2}}(\omega )}&{{a_{n - 3}}(\omega )}& \cdots & \cdots & \cdots &0\\ 0&{{a_n}(\omega )}&{{a_{n - 1}}(\omega )}&{{a_{n - 2}}(\omega )}&{{a_{n - 3}}(\omega )}& \cdots & \cdots &0\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0& \cdots & \cdots & \cdots &{{a_3}(\omega )}&{{a_2}(\omega )}&{{a_1}(\omega )}&{{a_0}(\omega )}\\ {{b_n}(\omega )}&{{b_{n - 1}}(\omega )}&{{b_{n - 2}}(\omega )}&{{b_{n - 3}}(\omega )}& \cdots & \cdots & \cdots &0\\ 0&{{b_n}(\omega )}&{{b_{n - 1}}(\omega )}&{{b_{n - 2}}(\omega )}&{{b_{n - 3}}(\omega )}& \cdots & \cdots &0\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0& \cdots & \cdots & \cdots &{{b_3}(\omega )}&{{b_2}(\omega )}&{{b_1}(\omega )}&{{b_{n0}}(\omega )} \end{array}} \right|_{2n \times 2n}} = 0\\ \end{array} $

（20）

结合式(12)及在穿越处特征根为纯虚数的特点，可看出根轨迹趋势的正负与时延τ_ck无关，即在穿越点处的RT仅与ω_ck相关^[19]。穿越处对应的时延τ_ck可通过式(11)给出，通过分析在每个时延区间对应的不稳定根个数即可给出闭环系统的稳定时延区间，具体细节将在第4节的算例中给出。值得一提的是，随着时延的增加，特征方程式(9)的根个数将增多，其动态特性将变得更加复杂，这点将在后面的根轨迹图中得到展现。

4 时间延迟稳定边界算例与应用

现以F-16飞机的短周期方程为例，在不同静稳定度情况下，求解其时间延迟稳定边界，求得的关系也可以用于在飞控系统时间延迟已知的条件下，确定可放宽静稳定度的边界。

此外，在初步设计阶段，通过计算流体力学(Computational Fluid Mechanics, CFD)等方法能给出较为可靠的静态气动导数，但动态导数C_{m$\dot{\alpha }$}、C_mq通常缺乏足够的精度，舵效在实际飞行中也存在较大的不确定性，因此舵效和动导不确定性与时延边界的关系也将在此进行分析。

4.1 静稳定度与时延边界的关系

以前文提到的F-16在马赫数Ma=0.46、高度H=100 m、配平迎角α₀=2.08°处的状态为例，其基本的重心位置为重心在焦点前方0.05c位置，即静稳定度为K_n=0.05，对应的短周期矩阵为

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} { - 1.0386}&1\\ { - 2.7206}&{ - 1.1132} \end{array}} \right]}\\ {\mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{l}} { - 0.1424}\\ { - 11.7839} \end{array}} \right]} \end{array}} \right. $	（21）

按照极点配置的方法，将其闭环极点配置到－3±j 3处，增稳控制律得到的相应增益为K=[0.869 30.316 1]，得到式(8)中各参数为

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{det}} (\mathit{\boldsymbol{A}}) = 3.875{\kern 1pt} {\kern 1pt} {\kern 1pt} 8}\\ {\eta = 14.124{\kern 1pt} {\kern 1pt} 2}\\ {{\rm{tr}} (\mathit{\boldsymbol{A}}) = - 2.151{\kern 1pt} {\kern 1pt} {\kern 1pt} 8}\\ {{\rm{tr}} (\mathit{\boldsymbol{BK}}) = - 3.848{\kern 1pt} {\kern 1pt} {\kern 1pt} 2} \end{array}} \right. $	（22）

相应的式(14)展开为

$ \begin{array}{l} h(s, {T_{ck}}) = {\rm{CE}}(s, {\tau _{ck}})(1 + {\rm{j}}{T_{ck}}) = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} [ {\rm{tr}} (\mathit{\boldsymbol{A}}){\omega _{ck}} - {\rm{tr}} (\mathit{\boldsymbol{BK}}){\omega _{ck}}] \cdot T + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} [\eta + {\rm{det}} (\mathit{\boldsymbol{A}}) - \omega _{ck}^2] + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{j}}\{ ( - \omega _{ck}^2 + {\rm{det}} (\mathit{\boldsymbol{A}}) - \eta ) \cdot T + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} [ - {\rm{tr}} (\mathit{\boldsymbol{A}}){\omega _{ck}} - {\rm{tr}} (\mathit{\boldsymbol{BK}}){\omega _{ck}}]\} \end{array} $

（23）

$ \begin{array}{l} {\rm{RT}} ({h_R}, {h_I}) = \\ \left| {\begin{array}{*{20}{c}} { {\rm{tr}} (\mathit{\boldsymbol{A}}){\omega _{ck}} - {\rm{tr}} (\mathit{\boldsymbol{BK}}){\omega _d}}&{\eta + {\rm{det}} (\mathit{\boldsymbol{A}}) - \omega _d^2}\\ { - \omega _{dd}^2 + {\rm{det}} (\mathit{\boldsymbol{A}}) - \eta }&{ - {\rm{tr}} (A){\omega _{ck}} - {\rm{tr}} (\mathit{\boldsymbol{BK}}){\omega _{ck}}} \end{array}} \right| = \\ 0 \end{array} $

（24）

由式(24)可以得到

$ \left\{ {\begin{array}{*{20}{l}} {{\omega _{c1}} = 5.023{\kern 1pt} {\kern 1pt} {\kern 1pt} 8}\\ {{T_{c1}} = 0.849{\kern 1pt} {\kern 1pt} 4} \end{array}} \right. $

说明该闭环系统仅存在一个穿越频率，即只存在一个根轨迹趋势值。

相应的延迟为

$ {\tau _{c1}} = 0.280{\kern 1pt} {\kern 1pt} {\kern 1pt} 3 + 0.398{\kern 1pt} {\kern 1pt} {\kern 1pt} 1p\pi {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} p = 0, 1, \cdots , \infty $

相应的根轨迹趋势用到的参数为

$ \left\{ {\begin{array}{*{20}{l}} {{P_0}({s^\gamma }) = {s^2} - {\rm{tr}} (\mathit{\boldsymbol{A}})s + {\rm{det}} (\mathit{\boldsymbol{A}})}\\ {{P_1}({s^\gamma }) = - s + \eta - {\rm{tr}} (\mathit{\boldsymbol{BK}})} \end{array}} \right. $	（25）

在该处根轨迹趋势为RT=+1，表明该穿越频率ω_c1对应的穿越为从$\mathbb{C}$^－→$\mathbb{C}$⁺的不稳定穿越，根据时延区间与根轨迹趋势关系得到的分段时间延迟区间的稳定性情况如表 1所示，从表中可以看出，仅在τ∈[0, 0.280 3]s区间内，该闭环系统稳定，即该静稳定度条件下，其闭环时间延迟稳定边界为τ=0.280 3s。

改变飞机的静稳定度，按照此方法可以求出K_n在一定范围内变化时所对应的时间延迟边界，如图 5所示。图中标注了两条边界，左边一条为K_n=0时对应的时间延迟边界值，针对算例中的飞机，其值τ=0.260 3s，该值为保证飞机可控时，放宽飞机到静不稳定条件下，飞控系统可接受的最大时间延迟。值得一提的是，在电传飞控增稳系统出现之前，考虑到飞行员的操纵负担与数百毫秒量级的操纵延迟，当时战斗机是不允许出现静不稳定布局的，该计算结果也对该现象提供了佐证。第二条边界对应为当前的飞控系统的120 ms延迟，通过该值，可以确定该飞机在此状态点处，目前可放宽的最大静不稳定度大约在25%。

为了从根轨迹上直观展示特征式(8)在该状态点时，其特征根随时间延迟的变化趋势，并验证上述方法计算结果的正确性。这里使用柯西多项式匹配寻根(Quasi-Polynomial mapping based Rootfinder, QPmR)函数求解特征多项式(9)在该状态点处，不同时间延迟分别对应的多项式特征根。该函数通过在指定区间搜索迭代以给出式(9)的数值结果，具体细节见文献[20-21]。需要指出的是，不同于以QPmR为代表的解空间遍历检索的数值计算方法，文中的求解方法是基于理论方法结合数值计算给出的精确结果，结果更可靠且需要的数值计算量很小。

该状态点上根轨迹随时延的变化曲线如图 6所示，图中将实轴和虚轴进行了压缩以便于观察。可以看到，随着延迟量的增加，特征方程的根将增多，从－3±3i出发的特征值首先穿越了虚轴，其穿越时对应的时间延迟与表 1中所求的结果一致。

4.2 时延边界对舵效的敏感性

考虑4.1节中的状态点和相应的增益，现在将舵效进行±30%的拉偏，分析舵效的变化对时延边界的影响，得到结果如图 7所示。

舵效偏大时，将导致时延稳定边界减小，且在小范围内，时延边界与拉偏系数成准线性关系，舵效偏小时，时延边界将增大，且其关系呈现出非线性的特征，即偏大的控制量将对时间延迟更为敏感，这对追求性能的战斗机控制律而言将更加明显。数值上，以原静稳定度为例，30%舵效拉偏将使边界从0.280 3 s降低到0.231 2 s，变化幅度为-17.5%，-30%舵效拉偏将使边界提高到0.366 1 s，变化幅度为30.6%。

4.3 时延边界对动导数的敏感性

动导数主要体现气动的非定常特性，对飞机过失速区的运动特性有显著影响，其数值虽然可以基于准定常的原则对其进行估计，但在初步设计阶段，仅通过经验公式或简单CFD计算等方式无法获得足够精度的动导数数据。因此需要研究时延边界对动导数的敏感性问题。

如图 8所示，动导数的影响主要体现在C_mq、C_m上，考虑短周期方程中，将C_mq、C_m进行±30%拉偏，其对短周期方程的影响近似于将_q+进行相应拉偏，得到时延边界与动导数敏感性的关系。

由于动导数本身为小量，在静稳定度足够的区域，其对时延边界的影响几乎可以忽略，只有在静稳定度本身很小的区域，其影响才得以体现，数值上，30%动导拉偏将使边界从0.280 3 s提高到0.298 2 s，变化幅度为6.4%，-30%动导数拉偏将使边界降低到0.263 8 s，变化幅度为-5.9%。

4.4 不确定性因素的综合影响

对静稳定度与时延稳定边界的影响，舵效不确定性相对动导数不确定性更为明显。在实际设计中应综合考虑两者的影响，以前文所述状态点为例，考虑两者综合影响的结果如图 9所示。此时，最小时延边界变为了0.181 7 s，变化幅度为-35.2%，因此，综合不确定性因素分析是有必要的，但可能导致过于保守的结果，这需要设计人员的权衡。

5 结论

1) 分析了目前电传飞控系统中时间延迟的来源，指出在初步设计阶段进行可控性评估时可以主要考虑输入时间延迟，分析了静稳定度与短周期方程参数间的关系。

2) 根据时延系统根轨迹穿越虚轴时为纯虚数的特点，并结合根轨迹穿越时其穿越方向仅与穿越频率相关而与具体时延无关的特点，给出了基于短周期方程的闭环时延系统的稳定边界的精确求解方法。

3) 针对短周期方程参数中的主要不确定性因素，分析了其对时延边界的影响，研究表明舵效不确定性是影响时间延迟稳定边界主要因素，但在静稳定度较低甚至静不稳定的情况下，也应考虑到动导不确定性对结果的影响，实际设计中，应综合两者的影响进行设计。

4) 解决了在飞机设计初步阶段，给定静稳定度布局情况下，基于可控性考虑的精确时间延迟稳定边界求解问题，该边界对飞控系统的设计具有指导意义；另一方面，该求解方法同时提供了一种在飞控系统时间延迟已知情况下，对静稳定度布局边界的约束指标。