3.4 建模实例

在本节中，我们将介绍一些系统建模的例子，以说明对于不同类型的系统，我们既可以建立微分方程模型，也可以建立差分方程模型。这些例子是特意从不同领域里挑选出来的，目的是想突出这样一个事实，即反馈和控制的概念可以广泛地应用于各种不同的系统。第4章将给出更为详尽的应用，以作为贯穿全书的运行实例。

3.4.1 运动控制系统

运动控制系统涉及用计算和反馈去控制一个机械系统的运动。运动控制系统的范围很广，从纳米定位系统（如原子力显微镜、自适应光学）到DVD播放机驱动的读/写头控制系统，从制造系统（如运输机械和工业机器人）到汽车控制系统（如反锁制动、悬浮控制、牵引控制），再到大气与空间飞行控制系统（如飞机、卫星、火箭以及行星登陆器）等。

例3.11 车辆转向

运动控制中的一个共同问题是如何用驱动器改变交通工具行驶方向，从而控制交通工具的运动轨迹。例如，汽车通过使用方向盘实现转向，自行车通过掌握前轮车把实现转向，船舶的转向以及飞机的俯仰动态控制也是同样的动态问题。在许多情况下，可以利用系统的基本动力学特征构建简单模型来理解和分析这类问题。下面以车辆转向为例。

考虑图3.17所示的具有固定后轴的前轮转向的车辆模型及其简化后的自行车模型。从转向角度来看，人们感兴趣的是描述车辆速度对转向角 δ 依赖关系的模型。具体来说，令 b 为轮轴距，设质心速度为 v ，质心到后轮的距离为 a 。令 x 和 y 为质心的坐标， θ 为方位角， α 为速度 v 与车辆中心线的夹角，点 O 为前、后轮法向的交点。

假定车轮不打滑，那么车辆的转向运动是绕着图3.17中的 O 点旋转产生的。假定从旋转中心 O 到后轮触地点的距离为 r _r ，那么从图3.17可得 b = r _r tan δ ， a = r _r tan α ，这意味着tan α =（ a / b ）tan δ ，所以可以得到 α 和转向角 δ 之间的关系为：

如果车辆质心处的速度是 v ，那么质心的运动由下式确定：

为了弄明白方位角 θ 如何受转向角 δ 的影响，观察图3.17，可以发现质心到旋转中心 O 的距离是 r _c = a /sin α 。因此车辆绕 O 点旋转的角速度是 v / r _c =（ v / a ）sin α 。因此有

式中的近似只对小的 δ 和 α 成立。

式（3.25）～式（3.27）可以用来对车轮不打滑且两个前轮可以用汽车中心上的单轮来近似的汽车进行建模。这个模型常称为 自行车模型 （bicycle model）。多添加一个状态变量以满足不打滑的要求，可得到一个用于描述船舶的转向动态以及飞机、导弹的俯仰动态的更实用的模型。也可以将坐标系的参考点选在后轮（对应于令 α =0），这样得到的模型常称为 杜宾斯汽车 （Dubins car）模型 ^[79] 。

图3.17表明，将速度的符号取反，那么这个模型也适用于后轮转向的场合。

例3.12 矢量推力飞机

现在考虑矢量推力飞机的运动，以图3.18a所示的鹞式战斗机为例。鹞式战斗机是一种垂直起降的喷气机，它通过向下调整推力，并利用机翼上的小型加速推进器实现垂直起降。鹞式战斗机的简化模型如图3.18b所示，这里重点考虑通过机翼的垂直平面内的机身运动。我们将朝下推进器和加速推进器产生的力分解为一对作用于机身下方 r 处的力 F ₁ 和 F ₂ （由推进器的几何布置决定）。

令（ x ， y ， θ ）表示飞机质心的位置和方向。令 m 为飞机的质量， J 为惯量， g 为重力加速度常数， c 为阻尼系数。那么机身的运动方程为：

为方便分析，将系统零输入时的某个平衡点作为原点，来重新定义输入。令 u ₁ = F ₁ 、 u ₂ = F ₂ -mg ，则方程变为：

这是由3个耦合的二阶微分方程构成的方程组，它描述了机身的运动。

3.4.2 热流体系统

热流体系统通常用于过程控制、发电，以及建筑物和汽车的加热通风与空调。这些过程牵涉到流体的运动和能量的传递，典型的过程包括热交换器、蒸发器、冷却器和压缩机。由于两相流的存在，流体的动态通常是很复杂的，精确的建模往往需要用到偏微分方程和计算流体动力学。图3.19给出了两个例子。

例3.13 热水器

考虑图3.19a中的热水器，它是一个截面积为 A 的圆柱形水箱。设水的质量为 m ，温度为 T ，流入和流出的水流量分别为 q _in 和 q _out ，流入水的温度为 T _in ，流出水的温度为 T 。总质量 m = ρAh ，其中 ρ 为水的密度， h 为水位。 C 为水的比热容， mCT 是总能量。系统可以用质量平衡和能量平衡来建模，因此有：

式中， P 是来自电热器的功率。这里没有考虑能量的损失，并假定水箱中所有的水都具有相同的温度。

假定 C 恒定，将能量平衡的导数展开，可得：

从此方程中解出d T /d t ，并利用质量平衡消去d m /d t ，式（3.30）表示的质量平衡和能量平衡方程可以重写为：

式中，总质量 m 和温度 T 是状态变量，输入功率 P 和流入流量 q _in 为控制（输入）变量，流入水温度 T _in 和输出水流量 q _out 为干扰。

例3.14 汽包锅炉

汽包锅炉是一种用于产生蒸汽的设备，它可作为发电系统的一部分，产生蒸汽，以驱动与发电机相连的汽轮机。汽包锅炉中的汽包与电热水器有许多共同的特性，但是它有两点更复杂的地方：一是材料常数 ρ 和 C 与状态有关，二是在上升管和汽包中都存在着水和蒸汽的混合物。尽管仍然可以采用质量和能量平衡来建模，但是两相流的存在导致了显著的复杂性，以下对此进行简单的讨论。汽包锅炉如图3.19b所示。

汽包的液位控制是一个关键问题：液位过低，管子会烧穿；液位过高，水就可能进入汽轮机，从而对汽轮机的叶片造成损伤。这里重点对汽包的液位进行建模。系统的进水由进水阀控制，并以蒸汽的形式经蒸汽阀离开汽包。水在由汽包、下降管、上升管构成的回路中循环，并在上升管中得到加热。下降管和上升管中水的密度差异建立起水的自循环。图3.19b中仅绘出了一根上升管和一根下降管，但此处讨论的锅炉实际上有22根下降管和788根上升管，汽包体积为40 m ³ 。在下降管中以及在上升管的底部，装有纯水。蒸汽是通过加热上升管产生的，蒸汽量沿着上升管的高度增加。在汽包中是蒸汽和水的混合物。

考虑系统处于平衡状态并突然打开蒸汽阀的情况。此时，离开系统的蒸汽量增加，因而一般人都会以为汽包的液位将会下降。但这种情况并不会发生，因为当更多蒸汽离开系统时，汽包内的压力会降低，上升管和汽包中的气泡随之增大，液位就会开始升高。如果继续保持蒸汽阀打开，那么液位最终会开始下降。汽包液位与进水量之间的动态关系也具有类似的特性。如果进水流量增加，那么汽包内的水温将下降，气泡将塌缩，因而汽包液位将开始下降。这种效应，即所谓的 假水位 （shrink and swell）或 逆响应 （inverse response），使得汽包液位的控制变得很困难。

图3.20用一台中型锅炉的模拟数据和实验数据展示了逆响应效应。在图中可以清楚地看到逆响应特性。仿真所用的模型是一个基于质量、能量及动量之平衡的五阶模型，详情可参见文献[18]。

由于从进水到汽包液位的逆响应动态特性，使得汽包液位的控制相当困难。因此，在系统中设置了蒸汽阀和进水阀两个传感器，如图3.19b所示。多增加的两个传感器使得人们可以预测系统中水量和蒸汽量的增减。在14.4节中将讨论系统具有逆响动态时的后果。

3.4.3 信息系统

信息系统的范围很广，从因特网这样的通信系统，到处理数据或管理企业级资源的软件系统，都是信息系统。反馈存在于所有这些系统中，为路由规划、流量控制及缓冲管理等功能设计策略是这种系统里很典型的问题。排队理论中的很多成果都源自电信系统的设计以及最新的因特网与计算机通信系统的发展 ^[43,149,218] 。管理排队以防止拥堵是这种系统的一个中心问题，因此我们将从排队系统的建模开始讨论。

例3.15 排队系统

图3.21是一个简单排队系统的原理图，请求到达后，先排队，再处理。到达速率和服务速率变动很大，当到达速率大于服务速率时，队列就变长。当队列变得太长时，就利用准入控制策略来拒绝服务。

该系统可以用很多不同的方法来建模。有个方法是对每个进入的请求进行建模，这样就得到了一个基于事件的离散状态模型，其状态是代表队列长度的一个整数。当一个请求到达时或一个请求获得服务时，队列就会改变。到达与服务的统计通常被建模为随机过程。在许多情况下，有可能确定出队列长度和服务时间等物理量的统计规律，但相应的计算会十分复杂。

用一个连续变量来近似表示离散队列的长度可以得到一个显著简化的模型。我们不再跟踪每个请求，而是将服务和请求看成连续流。这样得到的模型称为 流模型 （flow model），因为这种做法类似于流体力学中用连续流体来代替分子运动。因此，如果队列长度 x 是一个连续变量，到达和服务分别是流量为 λ 和 μ 的流体，那么该系统就可以用一阶微分方程来建模：

这个模型是Agnew提出来的（参见文献[5]）。服务速率 μ 取决于队列长度，如果没有容量的限制，将有 μ = x / T ，这里 T 是系统服务一个客户所花的时间。因此服务速率随队列长度线性增长。实际上，由于较长的队列需要更多的资源，因此增长率将随着队列的增长而变慢，服务的速率存在上限，即 μ _max 。将服务速率建模为 μ _max f （ x ）可以捕获这些影响，这里 f （ x ）为单调函数，当 x 较小时，它近似为线性，但 f （∞）=1。

对于具体的排队，函数 f （ x ）可以通过测量不同的到达速率和服务速率时的队列长度，用经验来确定。一个简单的选择是取 f （ x ）= x /（1+ x ），这可以得到下面的模型：

Tipper证明（参见文献[240]），如果到达过程和服务过程都是泊松过程，那么平均队列长度将由式（3.33）给出。

为了弄清楚式（3.33）这个模型的特点，下面先研究到达速率 λ 为常数时队列长度平衡点的数值。令式（3.33）中的导数d x /d t 为零，求解 x ，可以发现队列长度 x 趋近于以下的稳态值：

图3.22a所示为队列长度随 λ / μ _max （即服务强度）变化的函数关系。请注意，当 λ 趋向于 μ _max 时，队列长度迅速增大。要保持队列长度小于20，则要求 λ / μ _max <0.95。可以证明，响应一个请求的平均服务时间是 T _s =（ x +1）/ μ _max ，并且当 λ 趋向于 μ _max 时，这个时间急剧增长。

图3.22b所示为服务器在典型过负荷情况下的行为。该图表明队列增长得很快，但消减得很慢。由于响应时间正比于队列长度，因此这意味着长时间过负荷之后，服务质量将很差。这种行为称为 高峰时间效应 （rush-hour effect），已经在Web服务器及许多其他的排队系统（如汽车交通系统）中被观测到。

图3.22b中虚线所示为流体模型给出的平均队列长度曲线。可见，简化模型定性地抓取到了系统的行为特征，但在队列长度短的情况下，样本之间的一些变动在简化模型中难以反映出来。

许多复杂的系统采用离散控制作用。对于这种系统，可以对其中每个控制作用的具体情况进行分别描述，来实现整个系统的建模。下面举例说明。

例3.16 虚拟内存页面交换控制

反馈用于计算机系统中的一个早期例子，是用在IBM370计算机的OS/VS操作系统中 ^[54,67] 。该系统使用了虚拟内存技术，以使程序可以定址访问的内存要多于实际可得的快速内存。当数据位于当前快速内存（即随机存取内存）时，则进行直接访问，而当数据位于较慢的存储体（磁盘）时，则先被自动加载到快速内存。系统采用了特殊的实现方法，使得在编程者眼中系统仅拥有单一的大块内存。该系统在许多场合都运行得相当出色，但在过载情况下，就会遇到执行时间变长的问题，如图3.23a中“ ”标记的点所示。这个问题通过采用简单的离散反馈系统得到了解决。系统同时对中央处理器（CPU）的负荷以及快速内存与慢速存储体之间的页面交换数进行监测。系统运行的区域被划分为三种状态：正常状态、低负荷状态和过负荷状态。正常状态的特点是CPU活动性高；低负荷状态的特点是CPU活动性低、页面交换少；过负荷状态的特点是具有中等到偏低的CPU负荷，但有大量的页面交换（见图3.23b）。三个区域间的边界以及负荷测量所用时间是基于对典型负荷的仿真确定的。相应的控制策略是这样的：在正常负荷下不作任何控制；在过负荷情况下则从内存中移出进程；在低负荷情况下则容许载入新的进程或载入以前移出的进程。图3.23a中“ ”标记的点表明，该简单反馈系统对所仿真的负荷是有效的。在许多其他的场合，例如通信系统、Web服务器控制等场合，也使用了基于状态粗糙量化、简单启发式算法等类似的原理。

例3.17 传感器网络中的共识协议

传感器网络在许多场合得到了应用，人们希望利用经通信网络连接在一起的多个传感器来对一个空间区域中的信息进行采集和聚合。这种例子包括，监测地理区域（或一个建筑物）内的环境状况，监测动物或车辆的运动，监测一组计算机上所加载的资源，等等。在很多传感器网络中，计算资源是随着传感器分布的，因此这些分布式代理能否就某个属性（例如，一个区域的平均温度，一组计算机的平均计算负荷等）达成共识是极其重要的。

为了说明这样的共识是如何达成的，考虑计算一组数据平均值的问题，这些数据对于各个代理都是本地可用的。我们要设计一个“协议”（或算法），使所有代理都同意该平均值。此处假设通信网络是连接的（没有任何两组代理是彼此完全隔离的），但不是所有代理都能直接相互通信。图3.24a显示了这种情况的一个简单例子。

用图来对传感器网络的连接情况进行建模，用节点表示传感器，用边表示两个节点之间是否存在直接的通信链。对于任何这样的图，都可以建立一个 邻接矩阵 （adjacency matrix），该矩阵的每一行和每一列都对应一个节点，元素值为1表示相应行和列所对应的两个节点是相连的。对于图3.24a所示的网络，对应的邻接矩阵为：

用记号 N _i 来表示节点 i 的邻点的集合。例如，在图3.24a所示的网络中，有 N ₂ ={1，3，4，5}及 N ₃ ={2，4}。

为了求解共识问题，令 x _i 为第 i 个传感器的状态，它对应于传感器 i 对想要计算的平均值给出的估算值。将状态初始化为各个传感器测得的数值。这样一来，共识协议（算法）就可以作为一个局部更新法则来实现：

这个协议试图基于各个代理的邻居的值，来更新各个代理的局部状态，从而算出平均值。所有代理的组合动态可以写成以下形式：

式中， A 是邻接矩阵； D 是对角矩阵，其对角元素的值等于相应节点的邻居个数。常数 γ 是更新速率，节点将根据来自邻近节点的信息，以这一速率来更新平均值的估算值。定义矩阵 L = D-A ，称作该图的 拉普拉斯矩阵 （Laplacian）。

式（3.36）的 平衡点 （equilibrium point）是满足 x _e [ k +1]= x _e [ k ]的状态的集合。可以证明，如果网络是连通的，那么 x _e =（ α ， α ，…， α ）是系统的一个平衡状态，它对应于每个传感器对平均值的估算值都为同一个数值 α 。更进一步，可以证明 α 实际上是初始状态的平均值。由于图中可能存在 圈 （circle），因此系统的状态有可能进入无穷循环而永远无法收敛到希望的共识状态。正式的分析需要用到本书后面介绍的一些工具，不过现在就可以证明，对于任何连通图，总可以找到一个数 γ ，使各个代理的状态收敛到平均值。图3.24b给出的仿真结果说明了这个性质。尽管这里集中讨论的只是对一组测量值的平均值的共识，但通过选择适当的反馈机制，也可以实现其他共识状态。例如，求网络中的最大值或最小值，计算网络中的节点数，或计算一个分布量的高阶统计矩 ^[64,197] 。

3.4.4 生物系统

生物系统也许为反馈和控制提供了最为丰富的实例来源。在分子机器、细胞、生物体以及生态系统等各个层次的生物系统中，存在着许多复杂的反馈相互作用，像温度、血糖水平等物理量围绕着某个固定值进行调节的动态平衡问题，就是这种复杂反馈相互作用的例子。

例3.18 转录调控

转录是从DNA片段产生信使RNA（mRNA）的过程。基因的启动子区域允许其他的蛋白质控制转录，这种蛋白质被称作 转录因子 （transcription factor），它们与启动子区域结合，抑制或激活RNA聚合酶——一种从DNA转录产生mRNA的酶。然后基于mRNA的核苷酸序列转换出蛋白质。这一过程如图3.25所示。

转录调控过程的一个简单模型是采用Hill函数 ^[70,186] 。考虑某蛋白质A的调控问题，设该蛋白质的浓度为 p _a ，相应mRNA的浓度为 m _a 。设B为另一个蛋白质，其浓度为 p _b ，它通过转录调控来抑制蛋白质A的产生。这样一来， p _a 和 m _a 的动态方程可以写成：

式中， α _ab + α _a0 是未受调控的转录速率； δ _a 表示mRNA的降解率； α _ab 、 k _ab 及 n _ab 是描述B对A的抑制作用的参数； κ _a 表示从相应的mRNA产生蛋白质的速率； γ _a 表示蛋白质A的降解率。参数 α _a0 描述启动子的“泄漏”， n _ab 称为Hill系数，它跟启动子的协调性有关。

当一种蛋白质在另一种蛋白质的产生中起激活作用（激活子）而非抑制作用（抑制子）时，也可以使用类似的模型。此时的方程具有以下形式：

其中的变量跟前面描述的相同。请注意，在激活子的情况下，若 p _b 为零，则生产率就是 α _a0 ≪ α _ab （抑制子的为 α _ab + α _a0 ）。随着 p _b 增大，表达式 中右边第一项趋向于 α _ab ，而转录速率则变成 α _ab + α _a0 （抑制子的为 α _a0 ）。由此可见，激活子和抑制子的作用方式正好相反。

为了说明如何使用这些模型，考虑由Elowitz和Leibler首先提出的 压缩振荡子 （repressilator）模型 ^[85] 。压缩振荡子是一个人工回路，其中有三个蛋白质以一个抑制另一个的方式形成一个环，其原理如图3.26a所示，其中的三个蛋白质分别是TetR、 λ cI和LacI。压缩振荡子的基本思想是这样的，如果TetR存在，那它就抑制 λ cI的产生。如果 λ cI不存在，那么就以不受调控的转录速率产生LacI，这进而抑制TetR。一旦TetR受到抑制，那么 λ cI就不再受抑制。如果该回路的动态特性设计得当，所获得的蛋白质浓度将发生振荡。

可以将式（3.37）复制三份，来对该系统建模，但是在各个方程中，要用TetR、 λ cI及LacI的恰当组合来代替A和B。这样一来，系统的状态就是 x =（ m _TetR ， p _TetR ， m _{λ cI} ， p _{λ cI} ， m _LacI ， p _LacI ）。图3.26b显示了参数取 n =2、 α =0.5、 k =6.25×10 ^-4 、 α ₀ =5×10 ^-4 、 δ =5.8×10 ^-3 、 κ =0.12、 γ =1.2×10 ^-3 以及初始条件取 x （0）=（1，200，0，0，0，0）时，三种蛋白质浓度的变化轨迹（基于文献[85]）。

例3.19 神经细胞

神经元是所有人类和动物控制系统的关键元件。神经元有多种类型：感觉神经元感受刺激，运动神经元控制肌肉和其他器官，中间神经元在其他神经元之间的信号传递中起中介作用。神经元往往连成网络，人脑有近1000亿个神经元。

神经元有三个部分： 细胞体 （soma）、轴突和树突，如图3.27a所示。细胞体的大小从4 μ m到100 μ m不等，轴突的长度从1 mm～1 m不等。细胞有一层细胞膜将其与外界环境（细胞外空间）隔开，膜上有分子尺度的离子通道让离子通过细胞膜，从而在细胞膜内外两侧产生电位差。当电位差达到临界值时将产生一个电脉冲（尖峰）。产生的脉冲频率范围为1 Hz～1 kHz，它沿着轴突往轴突终末传播。

神经元通过树突接收来自其他神经元的信号。在轴突和另一个细胞的树突之间的交界面上存在电化学反应，这使得神经元之间的信息传递成为可能。轴突终末有包含神经递质的囊泡，当轴突受到电脉冲刺激时，神经递质就会被释放到突触间隙中，如图3.27b所示。神经递质刺激细胞膜上的离子通道，使它们打开。离子通道的种类很多，两种常见的是钠离子（Na ⁺ ）通道和钾离子（K ⁺ ）通道。钾离子通道具有缓慢的兴奋作用，钠离子通道则具有快速的兴奋和缓慢的抑制作用。

神经细胞的动态是理解细胞中信号过程的一个基础机制。Hodgkin-Huxley方程是研究神经细胞动态的模型。它将细胞膜当成电容来建模：

式中， V 是细胞膜的电位， C 是电容， 和 是钠离子和钾离子通过细胞膜传输而引起的电流， I _leak 是漏电流， I _input 是细胞外部的激励。各个电流都满足欧姆定律，即

式中，电导 g _Na 、 g _K 、 g _leak 通过变量 m 、 n 、 h 与电压 V 相关， g _Na 正比于 m ³ h ， g _K 正比于 n ⁴ ， g _leak 为常数。变量 m 、 n 、 h 由以下微分方程给定：

式中，函数 m _a 、 h _a 、 n _a 、 τ _m 、 τ _h 和 τ _n 来自实验数据；函数 m _a 和 n _a 随 V 单调递增，引起兴奋性行为；函数 h _a 单调递减，产生抑制行为；时间常数 τ _m 比时间常数 τ _h 和 τ _n 几乎小一个数量级。

平衡电位 E _Na+ 和 E _K+ 由 能斯脱定律 （Nernst's law）给出：

式中， R 是玻耳兹曼（Boltzmann）常数， T 是绝对温度， F 是法拉第（Faraday）常数， n 是离子的电荷（或化合价）， c _i 和 c _e 分别是细胞内和细胞外液体的离子浓度。在20℃时，有 RT / F =20 mV， E _Na+ =55 mV， E _K+ =-92 mV。

Hodgkin-Huxley方程相当复杂，其中包含许多不同的时间尺度，为此人们提出了许多近似。其中一个近似是FitzHugh-Nagumo模型（习题3.11）。用该模型仿真了神经元对外部电流刺激的行为，结果如图3.28所示。系统最初处于当 I =0、 V =0的静止状态。在 t =5 ms时加入一个短电流脉冲，神经元被激励，发出一个脉冲作为响应。然后在 t =30 ms时神经元被持续激励，神经元开始发出尖峰。

Hodgkin-Huxley模型建立的初衷是为了定量地预测鱿鱼巨轴突的行为特性 ^[120] 。Hodgkin和Huxley因“对神经细胞放电的电化学事件的分析”而与J.C.Eccles共享了1963年的诺贝尔生理学奖。函数 m _a （ V ）、 n _a （ V ）和 h _a （ V ）的确定用到了1.4节介绍的电压钳。基于Hodgkin-Huxley模型的神经元动态模型还有许多其他的变种，例如文献[202]所介绍的模型。有些模型结合了常微分方程和离散转换，即所谓的 整合发放模型 （integrate-and-fire model）或混合系统。 zVnsFuuYQGAwg6ukrY2U8JNTL3zaHLAJfBP0au9dO6uy+Qgua56hOIy0GgouABVA