GCKontrol状态空间模块及典型应用

摘要:状态空间是控制工程中的一个名词，是一种多输入多输出线性系统的描述方式。状态空间可以视为一个以状态变数为坐标轴的空间，因此系统的状态可表示为此空间中的一个向量。本文将介绍如何选择合适的状态变量，以及求出相应的四个状态矩阵，使用GCKontrol中的状态空间模块，以搭建飞行动力学模型为例，介绍状态空间模块的典型应用。

一、状态空间

1.1 理论知识

在一个实际问题当中，我们应该抽象出实际问题对应的数学模型，这个数学模型通常是用相应的微分方程来进行描述的。当我们所研究的控制系统是一个多输入多输出的系统时，用传统的传递函数的方法来进行分析，就有点显得力不从心了。如果对应的微分方程是单一的，那既可以使用传递函数的方法，也可以使用状态空间的方法；如果是微分方程组的形式，一般情况下会使用状态空间方法。这时我们可以使用状态空间法来建立系统的数学模型，该系统模型的方程如下：

其中是u输入，x是相应的状态变量，y是系统的输出。这里的u和y既可以是标量，也可以是向量。若是标量，就是单输入单输出系统；若是向量，就是多输入多输出系统。

首先，建立状态空间模型最重要的是选择合适的状态变量。其次是要求解出相应的A、B、C和D四个状态矩阵。本文以一个四阶常系数微分方程为例来展示如何求解相应的A、B、C和D四个状态矩阵，以及如何使用状态空间模块来求解相应的常系数微分方程。

线性常系数微分方程如下：

对于上面的线性常系数微分方程来说，一般情况下选择的状态变量，都令它等于要求解的原函数本身，也就是如下的公式，设：

此时，将原来的线性常系数微分方程转换成了状态空间的形式，接下来使用状态空间这个模块，来求解这个线性常系数微分方程。

1.2 求解过程

使用GCKontrol模型库中时间函数库里的状态空间模块，来对此线性常系数微分方程进行求解。

状态空间模块内部求解器使用四阶龙格库塔方法(Runge-Kutta)求解各状态值，四阶龙格库塔法(Runge-Kutta)的每一步需要四次计算函数值f，其截断误差更低，计算的精度更高。支持用户自定义配置状态空间模块内部求解器步数，当状态空间频率比GCKontrol调度周期高时，通过设置求解器步数，有效保持了计算精度。

以上就是状态空间模块的具体使用方法，使用难点就在于如何去求解A、B、C和D四个状态矩阵，以及如何选定相应的状态变量。

二、飞行动力学模型

2.1 稳定状态飞行

在此提供的数学模型所使用的是F-16飞机缩比模型在NASA兰利风洞试验中获得的数据。该数据是用于马赫数Ma=0.6以下的速度范围，曾用于NASA驾驶模拟，研究放宽静稳定性飞机的机动和失速/过失速特性。为了便于数据自动采集和使用，该数据覆盖的迎角α范围（-20°~90°）和侧滑角β范围（-30°~30°）很宽。但在目前技术条件下，还无法对过失速区进行精确动力学模拟，此外，迎角超过25°后，飞机机动的俯仰操纵力矩不足。因此，为了在此使用该模型，我们将数据范围缩小到-10°≤α≤45°，β角因事而定。

F-16飞机具有前缘襟翼，它作为迎角α和马赫数的函数自动操纵，在飞行机动过程中随α变化而迅速响应。在数据有效的速度范围内，襟翼随马赫数的变化很小，因此，我们没有考虑这种依赖关系。忽略襟翼作动器的动力学特性，假定襟翼仅依赖α变化，将所有与襟翼运动无关的数据都融合到启动数据表中。在简化数据中保留了襟翼偏角限制（而不是速率限制）的影响。

对于稳态飞行来说，我们希望在发动机功率限定范围内能够指定飞行高度和速度矢量（即速度和航迹角），对于传统飞机，我们期待存在唯一组控制输入和其余状态变量的组合。所有控制变量（即油门位置、升降舵、副翼、方向舵偏角）仅通过气动力数据表的形式进入模型，我们一般无法确定控制输入的解析约束。因此，必须采用数值算法对4个控制输入进行调整。但状态变量不会出现这种情况。

由于只与正切坐标系位置矢量的NED高度分量有关，且它可以预先指定，因此可以暂时不考虑三个位置状态。首先考虑稳定平动飞行。状态变量ϕ、P、Q、R都等于0，而方位ψ可自由指定，剩下来需要考虑的只有V_T、α、β（或U、V、W）和θ。侧滑角不能随便指定，必须采用配平算法计算，以使侧向力为零。还有三个变量V_T、α、θ，前两个变量通过飞机重量与总升力平衡相互关联，因此只有两个变量（θ和V_T，或θ和α）可能独立指定。我们通常希望在稳定状态条件上增加航迹角（γ）的限制，所以我们终将选择指定V_T和γ。

因为大气密度随着高度变化，因此严格地讲，稳定状态飞行条件并不包括非零航迹角条件。然而，在任意给定高度，通常可以确定非零航迹角的配平条件，由于随后会确定爬升率（ROC）并得到非零航迹角的线化动力学模型。因此，我们将导出一般爬升率约束条件，该限制条件允许出现非零滚转角，这样就可以用与稳定盘旋飞行状态。

现在考虑稳定盘旋飞行，变量ϕ、P、Q和R不再为0。用欧拉角度率ψ描述转弯，即机头朝向的变化速率（机头初始朝向仍然可随意指定）。然后，给定姿态角ϕ（滚转角）和θ（俯仰角）利用运动方程确定状态变量P、Q和R之值。如果已知ϕ，则可根据ROC约束条件确定所需的θ值。下面我们考虑如何确定θ值。

稳定盘旋飞行的滚转角ϕ可以自由指定，但是，一般情况下，会出现较大的侧滑角，“打滑”转弯。在“协调”转弯中飞机会滚转一定的角度，使得不出现气动力沿机体y轴方向的分量。该条件是后面推导协调转弯约束条件的基础。可以发现，在协调转弯约束条件中包含有θ和ϕ因此，必须与ROC约束条件同时求解。

飞机的动力学特性由相对风速 （-v_rel）和大气高度所决定，因此取决于变量V_T、α、β和h。这种特性基本不受风速v_(W/e)影响。因此，如果为了进行动力学研究、希望确定稳定状态飞行条件时，我们可以将风速设为0。

无风时方程可写为：

在平面地球方程中，爬升率为sin，这是速度在NED坐标系负z方向的分量。因此，上述方程可化为：

其中，星号表示“不关注”的分量。如果将该方程展开，并整理成关于θ的方程，则结果为：

sinγ= asinθ-bcosθ

其中

a = cosαcosβ，b = sinϕsinβ   cosϕsinαcosβ

求θ解可得：

作为检验，在无侧滑水平飞行时，该方程退化为θ=α γ。

通过求解非线性状态方程，使所有状态导数值U ̇、V ̇、W ̇（或V ̇_T 、α ̇、β ̇）和P ̇、Q ̇、R ̇都为0，求得状态矢量和控制矢量，确定稳态飞行条件。比较方便的方法是采用合适的数值算法，以上述导数的平方和作为标量代价函数。然后采用函数最小化计算程序，调整控制变量和适当的状态变量，使该标量代价函数之值最小。下图说明了整个配平算法的逻辑流程针对具体飞机或状态方程，只需要对代价函数进行裁剪。

图1 稳态配平流程图

将F-16飞机平衡到最小阻力状态，在模型的有效速度范围内水平直线飞行，升降舵的配平偏角变化非常小，且会出现反常变化。在极低的速度下，动压很小，需要较高的升力系数来支撑飞机重量。这会使诱导阻力变大，由于迎角较大，支撑飞机重量的很大一部分力来自于发动机推力。因此，低速飞行时必须增大油门位置。在飞机速度接近声速时，由于阻力增加，油门位置也要增大，因此，油门位置随速度的变化曲线必然要经过一个最小值。

2.2 增稳

大多数高性能军用和民用飞机都需要以某种方式增强稳定性。有些军用飞机实际上是不稳定的，如果没有自动控制系统，飞机就不能飞。SAS通常利用传感器测量飞行器的体轴角速率，并对这些信号进行处理，反馈到伺服机械系统，驱动气动操纵面。在这种方式中，角速度与气动力矩成正比，得到其导数，产生运动阻尼效应。如果基本模态是不稳定的，或者，如果希望分别改变阻尼和固有频率，则需要增加反馈信号。

传统增稳系统设计将纵向和横航向分开，在大多数飞行条件下，都可以对飞机动力特性进行解耦。

俯仰增稳系统旨在提供短周期模态令人满意的固有频率和阻尼。该模态涉及的变量包括迎角和俯仰速率，将这些变量反馈到升降舵，修改频率和阻尼。如下图所示流程，如果短周期模态是小阻尼的、但其他方面都是合适的，则只需要俯仰速率反馈。如果频率和阻尼都不满意，或该模态是不稳定的，则必要进行迎角反馈。这种反馈不应影响长周期模态。例如，外反馈控制回路往往环绕俯仰增稳系统，提供自动驾驶功能。可以在外反馈回路连接后再安排增稳（内）回路反馈增益自动调节，以使整体性能最优。

图3 俯仰轴增稳

迎角反馈效应的物理原理如下：一个静态不稳定飞机，俯仰力矩曲线在某些迎角范围的斜率为正值；如果将感受到的迎角扰动反馈到升降舵舵机，产生俯仰恢复力矩，则在工作迎角附近范围可使俯仰力矩曲线的斜率变负更大。而且，总体俯仰力矩曲线和升降舵配平偏角不会受到影响，从而保持配平阻力和设计师赋予基本飞机的设计机动性特点。

迎角测量可以通过皮托（pitot）静压大气数据系统或安装在前机身侧边的小型“风标”获得，测量迎角应能覆盖很宽的飞行条件范围，其安装位置必须经过反复试验和校准。可将两个传感器分别安装在飞机两侧，以提供冗余度，还有可能通过信号平均消除由侧滑引起的测量误差。此外，为了将机体相对于来流的迎角与传感器位置的流场方向联系起来，还必须根据“指示迎角”、空速和马赫数（实时）计算出“真”迎角。迎角传感器信号通常因湍流而包含噪声，需要利用噪声滤波器减小进入控制系统的噪声量。

尽可能避免迎角反馈，因为要得到准确、响应快速、无噪声的迎角测量是很困难的，而且传感器容易受到机械损伤。迎角传感器噪声会使精确定位（如寻靶）成为困难，因此，通常要限制迎角反馈量。

俯仰速率传感器通常是机械陀螺装置，用以测量绕俯仰轴的（惯性）角速率。选择陀螺安装位置必须非常谨慎，避免飞机结构振动的影响。在理想化结构振荡的节点上，有角运动但没有位移；在波腹上，有位移没有角运动。因此，速率陀螺安装的首选位置是在最重要的结构模态的波腹上。然后，根据飞行试验调整陀螺安装位置。陀螺位置选择不好，会对操纵品质产生不利影响，极端情况下，可能导致飞行控制系统振荡。上图中的陀螺滤波器通常是必要的，以消除噪声和/或抵消结构模态振动。

选取F-16的飞行配平条件及结果如下表所示：

表1 F-16的配平飞行条件

上表给出了 F-16模型另外的配平条件，F-16模型气动力数据的x方向参考点在0.35c ̅处，该点称为“名义”重心位置。表中的数据包括名义重心、前重心、后重心、前重心稳态盘旋、前重心拉起等条件的数据。重心位置对稳定性的影响都包括在前重心和后重心的情况里。针对两种机动情况，用前重心位置说明不附加不稳定动力特性的机动影响。

F-16模型在表1中的名义飞行条件下的纵向（4状态）雅可比矩阵为:

单输入是升降舵偏角δ_e ，单位为度；两个输出是适当的反馈信号：迎角和俯仰速率。为了与输入保持一致，C矩阵将单位转换为度。

由系数矩阵得到的两个单输入单输出传递函数，任何一个都能体现此飞行条件下的运动模态。升降舵到迎角的传递函数为：

与稳定重心位置（如x_cg=0.3c ̅）的传递函数有所不同，这个传递函数不会出现通常的长周期和短周期模态的极点。s=0. 098表明，不稳定指数模态的时间常数约为10s。对应于复共辘极点的振荡模态的周期为33s、阻尼比为0.79，这像是具有短周期阻尼比的长周期模态。这种模式称为静不稳定飞机的“第三振荡模态”。

2.3 仿真结果

使用GCKontrol展示了F-16飞机俯仰角的稳态控制，通过升降舵控制俯仰角速率使其达到稳态，如下图：

图2 GCKontrol主页面全局图

工程中配置的仿真时间为10秒，仿真步长和采样步长为0.001秒，曲线结果展示中配置了3条输出曲线，分别为一维插值表elevator.y和状态空间FlightDynamics.y，其中FlightDynamics.y输出F-16飞机的俯仰角和俯仰角速率。

图4 Yt曲线

当我们将模型的闭环取消时，会发现此时的俯仰角结果发散，如下：

三、其他应用场景

3.1 模型算法控制（MPC）

通过模型来预测系统在某一未来时间段内的表现来进行优化控制。多用于数位控制，系统模型通过离散状态空间表达x_(k 1)=Ax_k Bu_k。MPC一般分三步。

步骤一：估计或测量读取当前系统状态；

步骤二：基于最优控制策略得到的控制器输出u_k,u_(k 1)...u_(k N)来进行最优化，对应的被控对象模型预测输y_(k 1),y_(k 2)...y_(k N)；

步骤三：在执行下一时刻控制时，只取u_k，然后回到步骤一进行滚动优化控制。

根据上面的过程可以看出，MPC在执行每一步都需要求解最优化问题（步骤二）获取控制输入u，因此对控制器的计算能力要求很高。

3.2 故障诊断

故障诊断就是将根据状态预报得到的下一步应该处于的理想状态与测控系统测量获得的实际状态进行匹配和比较，得出关于模型运行是否正常的结论。如果运行不正常，还需进一步确定模型处于哪一种故障模式。这些判断的过程在状态空间描述法的基础上将易于实现。

3.3 化工加工系统

化工过程通常是多变量、非线性的，将过程多个变量反映到一个多维状态数据空间里，采用状态空间方法分析可以更全面、直观和有效的描述过程状态的变化。多变量分析方法中的主成分分析PCA（Principle Component Analysis）是目前常用的方法。独立成分分析ICA（Independent Component Analysis）是近几年在PCA基础上发展起来的一种新，分解出的各分量之间相互独立的特点，使ICA在信号处理领域，尤其是盲源信号分离方面受到广泛的关注。化工过程的状态空间是确定的，可以通过过程的性质作出确定的判断。通过ICA算法从状态空间中计算出对应的独立分量，实现了过程状态控件确定。通过仿真，证实了这种状态空间确定方法的合理性。