如何实现F16的纵向增稳控制

今天，我们简要学习一下如下飞机的飞行控制问题。

1903莱特兄弟第一次有动力的飞行表明，飞机除了自身的稳定性外，驾驶员对于飞机操纵的可控性，同样具有重要意义。可以说，飞行控制技术是与飞机设计技术同步发展的，下面我们就来看一看飞行控制技术是如何演变的。

早期的有人驾驶飞机的飞行控制系统，往往是由拉杆、摇臂或钢索、滑轮等组成的机械传动结构。利用该机构，驾驶员的操纵指令（驾驶杆和脚蹬机构）直接与飞机的气动舵面（升降舵、副翼、方向舵等）相连接，通过体力实现对飞机的控制。

随着飞机包线的扩展、飞行性能的提升，各气动舵面载荷的增长，纵使飞行员体质再好，靠“体力”驾驶飞机已经显得力不从心。和汽车驾驶一样，慢慢演变出了助力装置来帮助驾驶员实现对气动舵面的控制，一个典型的机械操纵系统的示意图如下：

简单来说，就是驾驶员通过助力器来操纵舵面。

什么是增稳系统（Stability Augmentation）呢？——举个简单的例子，小区的门。以前，我们每次拉开了关猛地一关，门都会哐哐唧唧好多次才比上，噪音很大。现在小区的门都不会，因为大都安装了阻尼器。对于我们常见的二阶系统来说，合适的阻尼，是保证良好操纵性的基本条件。

飞机纵向运动，通过合适简化，也可以看成是某种二阶系统，随着飞机飞行高度、飞行速度的扩大，飞机的动态阻尼特性使下降的，飞机会变得更难操控。

所谓增稳系统，是在上述机械操纵系统的基础上，应用反馈原理而设计，旨在提高飞机动态稳定性的一种飞行控制系统，目的是为飞机提供符合要求的阻尼和自然频率（或固有频率）。典型的增稳系统示意图如下：

基本原理是在保持原机械操纵系统的同时，使用传感器测量飞机角速率、迎角和过载等信号，馈送至计算机，通过预定的控制律解算舵面的运动指令，驱动飞机的气动控制面，因此产生的气动力矩为飞机提供额外的运动阻尼和稳定性。采用增稳系统的飞机包括：麦道F-4，诺普T38等。

所谓控制增稳系统，是在增稳系统仅仅引入了飞机运动参量反馈的同时，又将驾驶员的操纵指令输送至增稳回路的一种增加了驾驶员控制指令的增稳系统。

对于大多数老式飞机来说，在人工操纵下，增稳系统就足以完成飞行控制的目的和满意的飞行品质。不过，对于高性能的飞机而言，随着飞行包线的扩大，飞行员应既能操纵飞机机动飞行至其性能极限，有要求其完成诸如精确跟踪目标等任务。众所周知，飞机的稳定性和机动性相互制约，控制增稳的概念，正是解决这一矛盾，由于在增加稳定性的同时，接受驾驶员的控制指令，如何合理的设置控制和增稳的信息量的时间匹配和量值协调，可将机动性和稳定性均发挥的恰到好处。控制增稳的示意图如下：

采用控制增稳的飞机有：

电传控制系统FBW（fly by wire）是英文的意译。一般而言，是指利用电气信号的形式，通过电缆实现飞行员对飞机进行操纵的系统。电传控制系统不仅取消了笨重的机械传动杆系，减小了机械传动的间隙、摩擦等非线性影响，最重要的是，电气信号易于综合、校正和转换，给飞行控制的设计带来极大的灵活性，使飞行控制系统从“机械操纵”跨入了“电气控制”。

电传控制系统，一般都设计成控制增稳系统，采用电传控制的飞机如下所示，其中F16是最早采用电传控制飞机。

我们今天就从简化的模型来看看F16的纵向增稳和控制增稳是如何实现的。

专栏前面的文章中，我们介绍了飞机平飞过程中纵向小扰动的状态方程为： \begin {bmatrix} 1&0&0&0\\ 0&1-Z_{\dot w}&0&0\\ 0&M_{\dot w}&1&0\\ 0&0&0&1 \end {bmatrix} \begin {bmatrix} \dot u\\ \dot{w} \\ \dot q \\ \dot \theta \end {bmatrix} = \begin {bmatrix} X_u&X_w&0&-g\\ Z_u&Z_w&U+Z_q&0\\ M_u&M_w&M_q&0\\ 0&0&1&0 \end {bmatrix} \begin {bmatrix} u\\ w \\ q \\ \theta \end {bmatrix}

+ \begin {bmatrix} X_{\delta e}&X_{\delta T_h}\\ Z_{\delta e}&Z_{\delta T_h}\\ M_{\delta e}&M_{\delta {T_h}}\\ 0&0 \end {bmatrix} \begin{bmatrix} \delta e \\ \delta {T_h} \end{bmatrix}

为了应用现代控制理论的一些成果，我们需要写成标准的状态方程的形式： \begin {bmatrix} \dot u\\ \dot{w} \\ \dot q \\ \dot \theta \end {bmatrix} = \begin {bmatrix} X_u&X_w&0&-g\\ Z_u&Z_w&U+Z_q&0\\ M_u+M_{\dot w}Z_u&M_w+M_{\dot w}Z_w&M_q+M_{\dot w}U&0\\ 0&0&1&0 \end {bmatrix} \begin {bmatrix} u\\ w \\ q \\ \theta \end {bmatrix}

+ \begin {bmatrix} X_{\delta e}&X_{\delta T_h}\\ Z_{\delta e}&Z_{\delta T_h}\\ M_{\delta e}+M_{\dot w}Z_{\delta e}&M_{\delta {T_h}}+M_{\dot w}Z_{\delta T}\\ 0&0 \end {bmatrix} \begin{bmatrix} \delta e \\ \delta {T_h} \end{bmatrix}

当然，我们更习惯用迎角数据，可以用哪个下列形式进行转换：

\alpha\simeq\frac{w}{U}

然后代入纵向运动状态方程，可以得到：

\begin {bmatrix} \dot u\\ \dot{\alpha} \\ \dot q \\ \dot \theta \end {bmatrix} = \begin {bmatrix} X_u&X_{\alpha}&0&-g\\ Z_u/U&Z_{\alpha}/U&1+Z_q/U&0\\ M_u+M_{\dot \alpha}Z_u/U&M_{\alpha}+M_{\dot \alpha}Z_{\alpha}/U&M_q+M_{\dot \alpha}&0\\ 0&0&1&0 \end {bmatrix} \begin {bmatrix} u\\ \alpha \\ q \\ \theta \end {bmatrix}

+ \begin {bmatrix} X_{\delta e}&X_{\delta T_h}\\ Z_{\delta e}/U&Z_{\delta T_h}/U\\ M_{\delta e}+M_{\dot \alpha}Z_{\delta e}/U&M_{\delta {T_h}}+M_{\dot \alpha}Z_{\delta T}/U\\ 0&0 \end {bmatrix} \begin{bmatrix} \delta e \\ \delta {T_h} \end{bmatrix}

飞机的纵向在小扰动下可以分成两个运动的重新平衡过程：力的再平衡和力矩的再平衡。其中力矩的再平衡过程运动频率更高，结束更快，是飞行员最关注的的模态，一般称之为短周期模态。而力的平衡过程运动频率较低，周期更长，对飞行员来说不那么敏感，一般称之为长周期模态。短周期模态主要由迎角扰动 \alpha 和俯仰角速率扰动 q 两个状态变量决定，因此此时速度和俯仰角还没来得及变化。长周期模态主要由速度扰动 u 和俯仰角扰动 \theta ，因此此时迎角扰动和角速率扰动可认为已经达到稳态。

对于F16飞机，在飞行高度为3000ft，飞行速度为800ft/sec时，

飞机的纵向小扰动方程为：

\begin {bmatrix} \dot u\\ \dot{\alpha} \\ \dot q \\ \dot \theta \end {bmatrix} = \begin {bmatrix} -0.0067&42.1567&-0.2598&-31.6936\\ -0.0001&-0.7186&0.9645&0\\ -0.0000&13.9842&-0.0943&0\\ 0&0&1&0 \end {bmatrix} \begin {bmatrix} u\\ w \\ q \\ \theta \end {bmatrix}

+ \begin {bmatrix}-0.2469&0\\-0.0013&0\\-0.1476&0\\ 0&0 \end {bmatrix} \begin{bmatrix} \delta e \\ \delta {T_h} \end{bmatrix}

对于短周期模态，可以认为速度还没来得及变化，即 u\simeq 0 ，可以将我们更关心的短周期运动进一步抽离：

\begin {bmatrix} \dot{\alpha} \\ \dot q \\ \end {bmatrix} = \begin {bmatrix} Z_{\alpha}/U&1+Z_q/U\\ M_{\alpha}+M_{\dot \alpha}Z_{\alpha}/U&M_q+M_{\dot \alpha} \end {bmatrix} \begin {bmatrix} \alpha \\ q \\ \end {bmatrix}

+ \begin {bmatrix} Z_{\delta e}/U&Z_{\delta T_h}/U\\ M_{\delta e}+M_{\dot \alpha}Z_{\delta e}/U&M_{\delta {T_h}}+M_{\dot \alpha}Z_{\delta T}/U \end {bmatrix} \begin{bmatrix} \delta e \\ \delta {T_h} \end{bmatrix}

从F16实际的数值来看，发动机对于短周期运动的影响也可以忽略，同时为简单期间，忽略 Z_q 的影响，方程可进一步简化：

\begin {bmatrix} \dot{\alpha} \\ \dot q \\ \end {bmatrix} = \begin {bmatrix} Z_{\alpha}/U&1\\ M_{\alpha}+M_{\dot \alpha}Z_{\alpha}/U&M_q+M_{\dot \alpha} \end {bmatrix} \begin {bmatrix} \alpha \\ q \\ \end {bmatrix}

+ \begin {bmatrix} Z_{\delta e}/U\\ M_{\delta e}+M_{\dot \alpha}Z_{\delta e}/U \end {bmatrix} \delta e

此时可以获得短周期的固有频率为：

\omega_{sp}=\sqrt{\frac{Z_{\alpha}M_q}{U}-M_{\alpha}}

阻尼比为：

\xi=-\frac{M_q+M_{\dot \alpha}+Z_{\alpha}/U}{2\omega_{sp}}

一般情况下， M_{\dot \alpha} 项也可省略：

\begin {bmatrix} \dot{\alpha} \\ \dot q \\ \end {bmatrix} = \begin {bmatrix} Z_{\alpha}/U&1\\ M_{\alpha}&M_q\end {bmatrix} \begin {bmatrix} \alpha \\ q \\ \end {bmatrix}+ \begin {bmatrix} Z_{\delta e}/U\\ M_{\delta e}\end {bmatrix} \delta e

此时，

\omega_{sp}\simeq\sqrt{-M_{\alpha}}

阻尼比为：

\xi=-\frac{M_q+Z_{\alpha}/U}{2\omega_{sp}}

对于F16飞机而言，其短周期运动方程为：

\begin {bmatrix} \dot{\alpha} \\ \dot q \\ \end {bmatrix} = \underbrace{\begin {bmatrix} -0.7186&0.9645\\ 13.9842&-0.0943\end {bmatrix}}_{A} \begin {bmatrix} \alpha \\ q \\ \end {bmatrix}+ \underbrace{\begin {bmatrix} -0.0013\\ -0.1476\end {bmatrix}}_{B} \delta e

搭个模型来看一下状态变量的响应情况：

迎角和角速率的响应为：

可见，迎角和角速率都趋向于无穷大，系统是不稳定的，这是因为F16本身就是一个静不稳定的飞机。当然，我们也可以从数学上来验证一下：

计算结果为：

可见极点有一个是正值，导致可系统的不稳定。因此，对于F16飞机，在增稳之前是无法正常飞行的，要想让飞机具有良好的飞行品质，就必须有合适的固有频率和阻尼比，为此，需要对飞机增加增稳系统（Stability Augmentation System）和/或控制增稳系统（Control Augmentation System）。

对于SAS系统，可以采用的反馈状态变量包括迎角 \alpha 和角速率 q ，可以使用任意一种状态变量或者他们之间的组合进行状态反馈。

当反馈为迎角 \alpha 时，有

\delta e=u-k_{\alpha}\alpha

带入到短周期运动状态方程：

\begin {bmatrix} \dot{\alpha} \\ \dot q \\ \end {bmatrix} = \begin {bmatrix} Z_{\alpha}/U&1\\ M_{\alpha}&M_q\end {bmatrix} \begin {bmatrix} \alpha \\ q \\ \end {bmatrix}+ \begin {bmatrix} Z_{\delta e}/U\\ M_{\delta e}\end {bmatrix} (u-k_{\alpha}\alpha)

变换成标准形式：

\begin {bmatrix} \dot{\alpha} \\ \dot q \\ \end {bmatrix} = \begin {bmatrix} Z_{\alpha}/U-k_{\alpha}Z_{\delta e}/U&1\\ M_{\alpha}-k_{\alpha}M_{\delta e}&M_q\end {bmatrix} \begin {bmatrix} \alpha \\ q \\ \end {bmatrix}+ \begin {bmatrix} Z_{\delta e}/U\\ M_{\delta e}\end {bmatrix} u

此时，有

\omega'_{sp}=\sqrt{(M_{\delta e}-\frac{Z_{\delta e}}{U}M_q)k_{\alpha}+\frac{Z_{\alpha}M_q}{U}-M_{\alpha}}\simeq \sqrt{M_{\delta e}k_{\alpha}-M_{\alpha}}

可见，此时比无迎角反馈时相比，固有频率增加到了 \omega'_{sp} ，如果我们知道了期望的固有频率\omega'_{sp}，就可以反算出期望的增益 k_{\alpha} 。

k_{\alpha}=\frac{\omega'^2_{sp}+M_{\alpha}-Z_{\alpha}M_q/U}{M_{\delta e}-Z_{\delta e}M_q/U}

但需要注意的是，此时

\xi'=-\frac{M_q+Z_{\alpha}/U-k_{\alpha}Z_{\delta e}/U}{2\omega'_{sp}}

也就是说，阻尼比会减小。

我们以F16为例，我们用根轨迹法来看一下增益的变化对极点的影响。

可见，整体上来看，随着增益的增加，阻尼比是减小的，与我们的分析是一致的。对于在原点附近的根轨迹，我们可以通过放大图来仔细观察：

可见，极点的出发点（x所在的位置）是两个独立的点，其中一个还在正半轴，导致系统的不稳定，随着增益 k_{\alpha} 的增加，两个独立的极点逐渐变为一对共轭极点，并且这对共轭极点的实部处于坐标系的负半轴，系统开始变得稳定。

代码如下：

当增益 k_{\alpha}=-100 时，系统的极点和阻尼比为：

可见，此时两个独立的极点变成了一对共轭极点-4.98e-01+7.91e-01i，极点的实部是负值，系统的是稳定的，阻尼比达到了0.533。MATLAB代码如下：

搭个Simulink模型验证一下：

迎角扰动和角速率状态变量的响应为：

当反馈为迎角 q 时，有

\delta e=u-k_q q

带入到短周期运动状态方程：

\begin {bmatrix} \dot{\alpha} \\ \dot q \\ \end {bmatrix} = \begin {bmatrix} Z_{\alpha}/U&1\\ M_{\alpha}&M_q\end {bmatrix} \begin {bmatrix} \alpha \\ q \\ \end {bmatrix}+ \begin {bmatrix} Z_{\delta e}/U\\ M_{\delta e}\end {bmatrix} (u-k_qq)

变换成标准形式：

\begin {bmatrix} \dot{\alpha} \\ \dot q \\ \end {bmatrix} = \begin {bmatrix} Z_{\alpha}/U&1-k_qZ_{\delta e}/U\\ M_{\alpha}&M_q-k_qM_{\delta e}\end {bmatrix} \begin {bmatrix} \alpha \\ q \\ \end {bmatrix}+ \begin {bmatrix} Z_{\delta e}/U\\ M_{\delta e}\end {bmatrix}u

此时，有

\omega'_{sp}=\sqrt{(\frac{M_{\alpha}Z_{\delta e}}{U}-\frac{M_{\delta e}Z_{\alpha}}{U})k_{\alpha}+\frac{Z_{\alpha}M_q}{U}-M_{\alpha}}\simeq \sqrt{-M_{\alpha}}

可见，此时比无角速率反馈时相比，固有频率 \omega'_{sp} 基本不发生变化，此时阻尼比为：

\xi'=-\frac{M_q+Z_{\alpha}/U-k_qM_{\delta e}}{2\omega'_{sp}}=\xi+\frac{k_qM_{\delta e}}{2\omega'_{sp}}

可见，阻尼比增加了。如果我们知道了期望的阻尼比\xi'，就可以反算出期望的增益 k_q 。

k_q=\frac{2\xi'\sqrt{-M_{\alpha}}+M_q+Z_{\alpha}/U}{M_{\delta e}}

B747在某飞行状态下纵向短周期方程为：

\begin {bmatrix} \dot{\alpha} \\ \dot q \\ \end {bmatrix} = \begin {bmatrix} -0.5242&0.9735\\-0.5382&-0.3767\end {bmatrix} \begin {bmatrix} \alpha \\ q \\ \end {bmatrix}+ \begin {bmatrix} -0.0286\\-0.4240\end {bmatrix} u

其固有频率和阻尼比为：

可见，由于B747是静稳定的飞机，在未增稳的情况下，固有频率为0.849Hz，阻尼比为0.531。假如我们现在加入角速率反馈，来看一下B747飞机短周期模态随着增益 k_q 变化有什么变化，其根轨迹如下图所示。

可见，随着反馈增益的增加，固有频率略有增高，阻尼比明显增大。

前面B747是静稳定飞机，在增加的角速率的反馈后，阻尼比增大，那对于静不稳定的F16飞机，如果只有角速率反馈，会出现什么情况呢？我们来看一下F16飞机短周期模态随着增益 k_q 变化有什么响应。

下图是只有角速率反馈的根轨迹图，可以看出，随着增益的增大，两个独立的极点仍为独立的极点，当增益大大一定程度后，会将坐标系右半轴的那个极点校正到左半轴，此时系统将获得稳定性，这时系统的响应为两个惯性环节的叠加，飞机响应将会变得很慢，阻尼比过大。

比如，当 k_q=-150 时，迎角和角速率的响应为：

可见，对于静稳定的飞机，角速率反馈增稳效果明显，而对于静不稳定的飞机，单单使用角速率反馈，效果并不是很好。

前面的分析可以看出，采用迎角反馈可以增加固有频率，采用角速率反馈可以增加阻尼比，那如果同时采用迎角反馈和角速率反馈呢？——这其实就是现代控制理论的全状态反馈（State Feedback），理论上可以将极点配置到任意位置，也就是任意的固有频率和阻尼比。下面我们就来分析一下。

当迎角和角速率同时作为反馈变量时，状态方程变为：

\begin {bmatrix} \dot{\alpha} \\ \dot q \\ \end {bmatrix} = \begin {bmatrix} Z_{\alpha}/U&1\\ M_{\alpha}&M_q\end {bmatrix} \begin {bmatrix} \alpha \\ q \\ \end {bmatrix}+ \begin {bmatrix} Z_{\delta e}/U\\ M_{\delta e}\end {bmatrix} (u-k_{\alpha}\alpha -k_qq)

转换为标准形式：

\begin {bmatrix} \dot{\alpha} \\ \dot q \\ \end {bmatrix} = \begin {bmatrix} Z_{\alpha}/U-k_{\alpha}Z_{\delta e}/U&1-k_qZ_{\delta e}/U\\ M_{\alpha}-k_{\alpha}M_{\delta e}&M_q-k_qM_{\delta e}\end {bmatrix} \begin {bmatrix} \alpha \\ q \\ \end {bmatrix}+ \begin {bmatrix} Z_{\delta e}/U\\ M_{\delta e}\end {bmatrix}u

这样，当我们知道期望的固有频率 \omega'_{sp} 和阻尼比 \xi' 时，就可以得到两个状态变量的各自的增益 k_{\alpha} 和 k_q ，其实就是解线性方程组，也不复杂，但是我们还是更习惯用MATLAB来做这个事情。比如，我们想把F16在30000ft和800ft/sec状态下的固有频率和阻尼比校正到2Hz和0.8，即极点校正到 p=-1.6\pm1.2i ，使用MATLAB的place函数可以很容获得状态增益，代码如下：

此时计算的反馈增益为： [k_{\alpha} \quad k_q]=[ -108.9878 -14.8517] ，我们将该增益代入验证一下：

计算结果为：

可见，和我们的设定是完全一致的。迎角和角速率状态变量响应为：

前面设计的SAS，主要是以飞机的状态变量进行反馈，以期望获得的“等效飞机”有较为满意的固有频率和阻尼比，这主要的增加了飞机的稳定性。实际中，我们还需要控制增稳CAS，也就是在增稳的同时，还要增加一定操纵性。

也就是说，在CAS的加持下，驾驶杆/脚蹬控制的不在是舵面，而是飞机的姿态，我可以根据需要确定控制指令类型，比如对过载（加速度的表征）敏感，就可以过载作为指令；当然还有更流行的C*指令，后面说。

根据上面的框图可以看出，这其实是一个跟踪问题。我们在J Pan：如何入门现代控制理论一文中说了，跟踪问题需要增加一个新状态变量（state）来表征输出和输入之间误差。如果这个状态变量为零，那输入和输出就一致了。

状态变量定义如下：

\upsilon=\int {y}-{r}dt

两边都微分一下，就是：

\begin{equation} \begin{aligned}{\dot \upsilon}&= {y}-{r}\\&=C\bm{x}-{r}\end{aligned} \end{equation}

把这个新变量扩充到原来的状态方程

\left\{ \begin{array}{lcl} \bm{\dot x}=A\bm{x}+B{u}\\ {y}=C\bm{x}+D{u} \end{array} \right.

可以得到：

\underbrace{ \begin{bmatrix} \bm{\dot x}\\ {\dot \upsilon}\end{bmatrix} }_{\bm{\dot x}_a} = \underbrace{ \begin{bmatrix} A &0\\C & 0\end{bmatrix} }_{A_a} \begin{bmatrix} \bm{ x} \\ {\upsilon}\end{bmatrix} + \underbrace{\begin{bmatrix} B\\ 0\end{bmatrix}}_{{B_a}} {u} + \underbrace{\begin{bmatrix} 0\\ -1\end{bmatrix}}_{{B_r}} {r}

我们把扩展的矩阵起个新的名字，分别为 A_a , B_a , B_r 。在这里我们先选用法向过载作为指令，此时有

\underbrace{n_z}_{y}=\frac{U}{g}\dot \gamma=\frac{U}{g}(\dot \theta-\dot \alpha)=\frac{U}{g}(q-\dot \alpha)

其中 \gamma为轨迹角，将状态变量代入，即可以的到;

\underbrace{n_z}_{y}=\underbrace{\begin{bmatrix} {-Z_{\alpha}}/{g} & -{Z_q}/{g}\end{bmatrix}}_{C}\begin{bmatrix} \alpha\\ q\end{bmatrix}+ \underbrace{\begin{bmatrix} -{Z_{\delta e}}/{g}\end{bmatrix}}_{D}\underbrace{\delta e}_{u}

代入F16的短周期运动方程，可以得到增广后的状态方程为：

\begin{bmatrix} \dot\alpha \\ \dot q \\ \dot \upsilon\end{bmatrix} =\begin{bmatrix} Z_{\alpha}/U & 1+Z_q/U &0\\ M_{\alpha} & M_q&0\\ {-Z_{\alpha}}/{g} & -{Z_q}/{g} &0\end{bmatrix} \begin{bmatrix} \alpha \\ q \\ \upsilon\end{bmatrix} +\begin{bmatrix} Z_{\delta e}/U\\ M_{\delta e}\\-{Z_{\delta e}}/{g}\end{bmatrix}u +\begin{bmatrix} 0 \\ 0\\-1\end{bmatrix}r

画成框图就容易多了：

接下来的任务就是该如何去确定反馈增益 k_i 。

根据上面的框图，可以得到输入 {u} 的表达式为：

{u}=-K_c\bm{x}-k_i{\nu} =-\begin{bmatrix} K_c &k_i\end{bmatrix} \begin{bmatrix} \bm {x} \\{\nu}\end{bmatrix}

令 \bm{x_a}= \begin{bmatrix} \bm{x} \\{\nu}\end{bmatrix} 以及 K_a=\begin{bmatrix} K_c &k_i\end{bmatrix}

则可以得到有输出反馈的系统状态方程：

\bm{\dot x_a}=\underbrace{(A_a-B_aK_a)}_{\text{new A}} \bm{x_a}+B_r{r}

因此，闭环控制系统的极点分布为：

\left| \lambda I-(A_a-B_aK_a) \right|=0

也就是说，我们只要配置矩阵 A_a-B_aK_a 的极点就可以了，具体代码如下：

注意，由于增加了一个积分器，因此，在极点配置的时候，需要增加一个负的实值极点。计算结果为：

代入Simulink模型仿真一下看看：

过载的响应曲线为：

迎角和角速率的响应为：

代入F16的短周期运动方程，可以得到增广后的状态方程为：

\begin{bmatrix} \dot\alpha \\ \dot q \\ \dot \upsilon\end{bmatrix} =\begin{bmatrix} Z_{\alpha}/U & 1+Z_q/U &0\\ M_{\alpha} & M_q&0\\ {-Z_{\alpha}}/{g} & -{Z_q}/{g} &0\end{bmatrix} \begin{bmatrix} \alpha \\ q \\ \upsilon\end{bmatrix} +\begin{bmatrix} Z_{\delta e}/U\\ M_{\delta e}\\-{Z_{\delta e}}/{g}\end{bmatrix}u +\begin{bmatrix} 0 \\ 0\\-1\end{bmatrix}r

画成框图就容易多了：

所谓C*控制指令，是用法向过载和角速率融合在一起的一个变量。用C*作为反馈信号，在飞机高速飞行时，飞行员主要感受到法向过载的变化；在低速飞行时，飞行员主要感受到的是俯仰角速率的变化。C*借助俯仰角速率反馈提高短周期阻尼比，改善动态特性；通过法向果子啊，提高短周期频率，改善操纵性。

C*可定义成：

C*=n_z+12.4q

此时

\underbrace{C*}_{y}=\underbrace{\begin{bmatrix} {-Z_{\alpha}}/{g} & -{Z_q}/{g}+12.4\end{bmatrix}}_{C}\begin{bmatrix} \alpha\\ q\end{bmatrix}+ \underbrace{\begin{bmatrix} -{Z_{\delta e}}/{g}\end{bmatrix}}_{D}\delta e

通过极点配置，计算增益，MATLAB代码如下：

计算结果为：

用Simulink验证一下，模型如下：

C*指令的响应曲线如下，可见可以很好的实现很好的跟随。

法向过载的响应如下：

迎角和角速率的响应如下：

当然，实际的飞控系统，舵面是靠作动器驱动的，需要假如作动器模型；迎角信号由噪音，需要滤波；角速率信号也需要滤波（此处暂不考虑），我们可以得到完整模型如下:

由前面介绍， 可知： C*=\underbrace{\begin{bmatrix} -Z_{\alpha}/g&-Z_q/g+12.4&-Z_{\delta e}/g&0&0\end{bmatrix}}_{H} \begin{bmatrix} \alpha \\ q \\\delta e \\ \alpha_f\\ \upsilon\end{bmatrix}

对于F16飞机而言， a=10 , b=20 ,代入F16的短周期运动方程，可以得到增广后的状态方程为：

\begin{bmatrix} \dot\alpha \\ \dot q \\ \dot \delta e \\ \dot \alpha_f\\ \dot \upsilon\end{bmatrix} =\underbrace{\begin{bmatrix} Z_{\alpha}/U & 1+Z_q/U&Z_{\delta e}/U&0 &0\\ M_{\alpha} & M_q&M_{\delta e}&0&0\\ 0&0&-20&0&0\\10&0&0&-10&0\\ -{Z_{\alpha}}/{g} & -{Z_q}/{g} +12.4&-Z_{\delta e}/g&0&0\end{bmatrix} }_{A_a}\begin{bmatrix} \alpha \\ q \\\delta e \\ \alpha_f\\ \upsilon\end{bmatrix}

+\underbrace{\begin{bmatrix} 0\\ 0\\-20\\ 0\\ 0\end{bmatrix}}_{B_a}u +\underbrace{\begin{bmatrix} 0 \\ 0\\0\\0\\-1\end{bmatrix}}_{G_a}r

输出状态方程为：

\begin{bmatrix} \alpha_f\\q\\e\\\upsilon \end{bmatrix} =\underbrace{\begin{bmatrix} 0&0&0&1&0\\ 0&1&0&0&0\\ -{Z_{\alpha}}/{g} & -{Z_q}/{g} +12.4&-Z_{\delta e}/g&0&0\\ 0&0&0&0&1\\ \end{bmatrix}}_{C_a}\begin{bmatrix} \alpha \\ q \\\delta e \\ \alpha_f\\ \upsilon\end{bmatrix}

+\underbrace{\begin{bmatrix} 0 \\ 0\\-1\\0\end{bmatrix}}_{D_a}r

当输入方程为：

u=-\underbrace{\begin{bmatrix} k_{\alpha}&k_q&k_p&k_i \end{bmatrix}}_{K} \underbrace{\begin{bmatrix} \alpha_f\\q\\e\\\upsilon \end{bmatrix}}_{y}

则状态方程则变为：

\bm{\dot x}=A_a\bm{x}-B_aK(C_a \bm{x}+D_a r)

即

\bm{\dot x}=\underbrace{(A_a-B_aKC_a)}_{A_c}\bm{x}+\underbrace{(G_a-B_aKD_a)}_{B_c} r

显然，要同时整定这些增益难度是比较大的，我们需要进 行简化，否则将无从下手。首先作动器的简化，由于作动器的时间常数为 1/20=0.05s ，而前面我们整定完后的飞机的短周期时间常数为0.6s左右，因此作动器环节可以简化成一个纯增益环节，即：

\frac{-20}{s+20}\Rightarrow -1

角速率的滤波器，同样可以进行简化：

\frac{10}{s+10}\Rightarrow 1

这样，系统的阶数就可以从5阶降到3阶，具体如下：

\begin{bmatrix} \dot\alpha \\ \dot q \\ \dot \upsilon\end{bmatrix} =\begin{bmatrix} Z_{\alpha}/U & 1+Z_q/U &0\\ M_{\alpha} & M_q&0\\ {-Z_{\alpha}}/{g} & -{Z_q}/{g} &0\end{bmatrix} \begin{bmatrix} \alpha \\ q \\ \upsilon\end{bmatrix} +\begin{bmatrix} Z_{\delta e}/U\\ M_{\delta e}\\-{Z_{\delta e}}/{g}\end{bmatrix}(-u) +\begin{bmatrix} 0 \\ 0\\-1\end{bmatrix}r

u=-\begin{bmatrix} k_{\alpha}&k_q&k_i \end{bmatrix} \begin{bmatrix} \alpha\\q\\\upsilon \end{bmatrix}

这就回到了之前我们分析的情况。

我们仍选取之前的固有频率和阻尼比，这样就可以计算出需要的状态反馈增益：

现在我们用简化模型得到的增益带到完整的模型里面，看看结果如何。其Simulink模型如下：

其C*响应为：

注意此时 k_p=0 ，结果和之前完全一样。我们看看假如作动器模型和滤波器模型后极点和阻尼比的变化情况。MATLAB代码如下：

计算结果如下：

可见极点从三个变成了五个，主导极点由-1.6±1.2i变为了-2.13±1.02i，阻尼比和固有频率都增加了。惯性关节的极点由-10变到了-7.27，时间常数变长。同时还增加了一对快速衰减的共轭极点-9.64+17.4i。可见，当采用简化模型设计增益时，可以适当将期望极点的阻尼比选小一点，固有频率也小一点。

现在考察一下 k_p 的影响，当然，最好的方法就是根轨迹法了。代码如下：

可见对于非主导的共轭极点，阻尼比减小，固有频率增加。对于主导的共轭极点，刚好相反，阻尼比增加，固有频率减小。

这样，我们就可以根据简化的模型先初选 k_a,k_q,k_i,在通过根轨迹法去选择 k_p ，从而设计出整个纵向短周期控制律。

注：本人只说明设计方法，具体参数值不具有参考意义。