第5 章控制系统的Simulink建模Simulink建模与仿真
Simulink 最早解决的问题是控制系统仿真，控制问题仿真研究也是Simulink
发展与普及的一个巨大推动力。早在Simulink 出现之前，绝大部分控制系统仿真语
言都是编程式的。Simulink 的图形化、模块化建模方法使控制系统仿真研究变得更
加容易，所以Simulink 逐渐成为控制系统仿真研究的首选计算机工具。Simulink 提
供了大量的控制系统仿真模块，尽管如此，很多问题尚不能完美地解决，需要对其
中很多模块进一步扩展，得到适用性更广的控制系统建模与仿真模块。本章将系统
地介绍基于Simulink 的控制系统建模与仿真的内容。
5.1 节介绍一般连续线性系统的建模问题，首先介绍基于Simulink 底层模块的
建模方法，然后介绍奇异系统建模以及具有非零初值的传递函数建模，最后介绍基
于控制系统工具箱LTI 模块的复杂线性系统建模的通用方法。5.2 节介绍各类离散
线性系统建模与分析的方法，并介绍复杂离散模型的处理方法。5.3 节介绍基于查
表模块的单值或双值非线性环节的一般搭建方法，目标是描述任意的静态非线性
特性。该节还通过例子演示二维查表模块的使用。5.4 节介绍Simulink 仿真框图模
块的自动排序问题与准则，以便读者能够更好地理解Simulink 模型的仿真过程与
仿真机制。5.5 节介绍控制系统的近似仿真方法，例如非线性系统的线性化问题、延
迟模型的Padé 近似，演示仿真过程代数环的概念与消除问题，并给出一种消除代
数环的通用方法。当然，Simulink 理论上可以对任意复杂度的系统作精确的仿真分
析。如果不是特别必要，无须对系统进行近似的分析。
5.1 连续线性模型的Simulink建模
在传统的控制课程与控制工程中，描述线性连续系统通常有两大类方法系
统的传递函数描述方法和系统的状态方程描述方法。对多变量系统而言，传递函数
方法还可以进一步扩展为传递函数矩阵的描述方法。标准的MATLAB 与Simulink
提供了三组连续线性模型的Simulink 模块，包括标准组中的Continuous 模块组、
第5 章控制系统的Simulink 建模·127·
Simulink Extras 模块集中的Additional Linear（附加线性）模块组和控制系统工具
箱模块集。本节将主要探讨各类连续线性系统模块的Simulink 建模方法。特别地，
将构建并探讨具有非零初值的传递函数建模问题，这是用现有Simulink 模块难以
解决的问题。
5.1.1 传递函数模型
线性控制系统理论中，传递函数是描述控制系统最重要的数学模型。传递函数
模型是建立在常系数线性微分方程基础上的，式（2-3-1）给出了其一般形式。其数
学形式如下
G(s) =
b1sm + b2sm?1 + · · · + bms + bm+1
a1sn + a2sn?1 + · · · + ans + an+1
（5-1-1）
其中，ai，i = 1, 2, · · · , n+1，bi，i = 1, 2, · · · ,m+1 为常数，n称为系统的阶次。如果
m < n，则系统称为严格正则（strictly proper）系统；若m = n，则系统称为正则系
统；若n < m，则系统是物理不可实现的。
传递函数是系统输出信号的Laplace 变换表达式与输入信号Laplace 变换表
达式的比值，通常可以认为是系统在s 域上的放大倍数。值得指出的是，传递函
数的概念是在系统输入、输出信号及其各阶导数的初值均为零的前提下给出的。
MATLAB 控制系统工具箱提供了两种传递函数模型的输入方法。第一种方法
是将分子与分母多项式系数按照s 的降幂次序提取出来，构成分子与分母多项式系
数向量num 与den，这样，调用G=tf(num,den) 命令，则可以建立起传递函数对象
G。将变量名G填写到LTI System 模块的参数对话框中，则可以在Simulink 下直接
表示传递函数模型；另一种方法是由s=tf('s') 命令先定义算子s，然后用简单的
命令就可以建立传递函数对象。如果传递函数带有延迟时间常数T，则有三种方法
可以直接描述该模型。
G=tf(num,den,'ioDelay',T);
G=tf(num,den); G.ioDelay=T;
s=tf('s'); G=tf(num,den)*exp(?T*s);
　　例5-1 如果已知系统的传递函数模型如下，试将其用传递函数对象的形式输入给
MATLAB 环境。
G1(s) =
s3 + 7s2 + 24s + 24
s4 + 10s3 + 35s2 + 50s + 24
, G2(s) =
s2 + 4s + 3
s2(s + 4)??(s + 1)2 + 3
 　　解从给出的模型可见，G1(s) 的分子与分母多项式系数已知，所以可以考虑用
tf() 函数直接构造传递函数模型，而G2(s) 模型中，分母多项式的系数不是直接已知
的，如果用常规的tf() 函数调用方式，需要首先展开各个因式，再合并同类项。所以，
·128· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
应该考虑先定义s 算子，再用简单数学表达式构造传递函数模型。对这个具体的例子而
言，可以直接给出下面的语句试图描述传递函数模型。
>> G1=tf([1 7 24 24],[1 10 35 50 24])
s=tf('s'); G2=(s^2+4*s+3)/s^2/(s+4)/((s+1)^2+3)
如果系统的传递函数已知，则可以用Continuous 模块组的Transfer Fcn 模块直
接描述。不过，若传递函数是物理不可实现的、传递函数带有延迟环节或传递函
数对应非零初值，则不可采用这个模块描述系统模型。如果已经按照例5-1 描述的
方式建立了传递函数对象，则可以由[n,d]=tfdata(G,'v') 命令提取分子和分
母多项式系数向量，由变量n和d返回。如果模型含有延迟时间常数，则可以由语
句T=G.ioDelay 提取T 的值。提取之后，就可以将其填写到Transfer Fcn 模块中，
直接在Simulink 下描述传递函数模型。
5.1.2 状态方程模型
线性系统的状态方程模型是另一类常用的数学模型，对物理可实现的系统而
言，可以直接将其写成系统的状态方程模型。本节先给出状态方程模型的一般形式，
然后指出该模型的局限性，并引入线性奇异系统的状态方程模型。
　　定义5-1 连续线性定常系统的状态方程模型为
(x′(t) = Ax(t) + Bu(t)
y(t) = Cx(t) +Du(t)
（5-1-2）
其中，x(t) ∈ Rn×1 称为状态变量向量，且已知常数矩阵A ∈ Rn×n，B ∈ Rn×q，
C ∈ Rp×n，D ∈ Rp×q。状态方程模型的状态变量在t = t0 时刻的初值x(t0) 为已
知常数向量，默认为零向量。为避免与其他状态方程形式混淆，本书必要时将这种
状态方程模型称为常规状态方程模型。
状态方程的输入可以利用MATLAB 控制系统工具箱提供的ss() 函数直接实
现，该语句的调用格式为G=ss(A,B,C,D)。状态方程模型是可以直接描述多变
量系统的。
　　例5-2 多变量系统的状态方程模型可以用前面介绍的方法直接输入，无须再进行
特殊的处理。试将下面的双输入双输出系统的状态方程模型输入MATLAB 环境。
8>
>>><>>>>:
x′(t) =264
?12 ?17.2 ?16.8 ?11.9
6 8.6 8.4 6
6 8.7 8.4 6
?5.9 ?8.6 ?8.3 ?6
375
x(t) +264
1.5 0.2
1 0.3
2 1
0 0.5
375
u(t)
y(t) =  2 0.5 0 0.8
0.3 0.3 0.2 1 x(t)
　　解由给出的状态方程模型可以直接识别出A、B、C 和D矩阵，这样，由ss() 可
第5 章控制系统的Simulink 建模·129·
以直接将状态方程模型输入MATLAB 工作空间中。
>> A=[-12,-17.2,-16.8,-11.9; 6,8.6,8.4,6; 6,8.7,8.4,6;
-5.9,-8.6,-8.3,-6];
B=[1.5,0.2; 1,0.3; 2,1; 0,0.5];
C=[2,0.5,0,0.8; 0.3,0.3,0.2,1]; D=zeros(2,2); % 输入4 个矩阵
G=ss(A,B,C,D) % 输入并显示系统状态方程模型
如果已知系统的传递函数模型G，则可以由G1=ss(G) 函数直接变换出等效
的状态方程模型，转换完成后，可以用[A,B,C,D]=ssdata(G1) 命令提取各个
矩阵。若已知状态方程模型G1，则可以由G=tf(G1) 命令变换出等效的传递函数
模型G。
状态方程模型可以由Continuous 模块组中的State Space 模块直接表示，用户
需要在该模块的参数对话框中直接填写A、B、C 和D矩阵。如果状态方程含有非
零初值，则在Initial States（初始状态）编辑框中填写x(t0) 向量即可。
更一般地，MATLAB 控制系统工具箱支持带有内部延迟的状态方程模型[47]。
由于内部延迟的数学形式比较复杂，这里不给出其数学形式与框图，有兴趣的读者
可以参考相关的文献。
5.1.3 线性奇异系统的状态方程模型
从给出的数学表达式看，常规状态方程模型只适合描述正则的传递函数模型，
即分子多项式的阶次不高于系统的阶次。在某些特定场合下，常规的状态方程模型
是不够的。例如，如果已知某个严格正则的传递函数模型，则其逆系统是物理不可
实现的。这样，用常规的状态方程模型就不能描述逆系统模型了。因此，需要扩展现
有的状态方程模型。本节将介绍线性奇异系统的状态方程模型。
　　定义5-2 连续线性定常系统的奇异状态方程模型为
(Ex′(t) = Ax(t) + Bu(t)
y(t) = Cx(t) +Du(t)
（5-1-3）
其中，E ∈ Rn×n 称为描述符矩阵（descriptor matrix）。如果E为单位矩阵，则系统
模型为定义5-1 描述的形式；如果E为奇异矩阵，则系统称为奇异系统；如果E?1
为有限元素矩阵，则引入z(t) = Ex(t)，则可以建立起关于z(t) 状态变量的常规状
态方程模型。
MATLAB 控制系统工具箱函数支持奇异系统的输入与分析。如果已知系统的
状态方程模型，则可以使用G=dss(A,B,C,D,E) 命令输入模型。如果将这5 个
矩阵输入Descriptor State-Space 模块，则可以在Simulink 下描述奇异系统。
·130· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
　　例5-3 考虑例5-1 中给出的严格正则传递函数模型G1(s)，试写出其逆系统的状
态方程模型。
　　解由于G1(s) 模型是正则的，所以其逆系统G2(s) = 1/G1(s) 是物理不可实现的，
该模型对应的常规状态方程模型不存在。必须使用奇异系统模型描述逆系统的状态方
程模型。由下面的语句可以直接求取其状态方程模型。
>> G1=tf([1 7 24 24],[1 10 35 50 24]); G2=1/G1;
G3=ss(G2), step(G3) % 变换出状态方程模型并尝试获得阶跃响应
可以得出奇异状态方程模型为
26664
0 0.25 0 0 0
0 0 0.25 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
37775
x′(t) =26664
1 0 0 0 0
0 1 0 0 0
0 0 ?7 ?6 ?3
0 0 4 0 0
0 0 0 2 0
37775
x(t) +26664
0
0
8
0
0
37775
u(t)
且输出方程为y(t) = [2, 5, 4.375, 1.5625, 0.375]x(t)。
当然，由于这样的奇异系统是物理不可实现的，所以调用step(G3) 命令绘制系统
阶跃响应曲线时，将给出错误信息“Cannot simulate the time response of improper
(non-causal) models（非正则系统不能得出时域响应）”。
　　例5-4 在某些特定物理系统的建模中，如果状态变量选择不当或选择了冗余的状
态变量，很容易将其状态方程模型写成奇异的形式。考虑如图5-1 所示的简单电阻、电
容和电感串联电路。如果状态变量选作i(t)、ul(t)、uc(t) 和ur(t)[48]，试写出该电路的状
态方程模型。
?
ur(t)
+
?
u(t) R
+
L
+ ul(t) ?
i(t)
C
+ ? uc(t)
图5-1 电阻、电容与电感的串联电路
　　解正常情况下[25]，在建模时选择状态变量x1(t) = uc(t)，x2(t) = i(t)，则可以得
出常规的二阶微分方程。不过，若按照本例要求取状态变量，则可以写出下面4 个微分
方程（其中有代数方程）：
8>
>><>>>:
i′(t) = ul(t)
u′c(t) = i(t)/R
ur(t) = Ru(t)
u(t) = ul(t) + uc(t) + ur(t)
第5 章控制系统的Simulink 建模·131·
若将其写成矩阵形式，不难整理出如下的奇异微分方程。
264
L 0 0 0
0 0 1 0
0 0 0 0
0 0 0 0
375
264
i′(t)
u′l(t)
u′c(t)
u′r(t)
375
=264
0 1 0 0
1/C 0 0 0
?R 0 0 1
0 1 1 1
375
264
i(t)
ul(t)
uc(t)
ur(t)
375
+264
0
0
0
?1
375
u(t)
且输出方程为y(t) = uc(t)。
　　例5-5 考虑例5-4 给出的奇异系统，若R = L = C = 1，试绘制其阶跃响应曲线。
　　解将R = L = C = 1 代入模型，则可以得出如下的奇异状态方程模型。
264
1 0 0 0
0 0 1 0
0 0 0 0
0 0 0 0
375
x′(t) =264
0 1 0 0
1 0 0 0
?1 0 0 1
0 1 1 1
375
x(t) +264
0
0
0
?1
375
u(t), y(t) = x3(t)
这样，就可以由下面的语句输入奇异系统的数学模型，然后调用step() 函数直接
绘制系统的阶跃响应曲线，如图5-2 所示。其实，若将其变换成零极点模型，则可以得出
G1(s) = 1/(s2 + s + 1)，说明该模型事实上就是一个二阶传递函数模型，冗余的状态会
被转换函数自动剔除。
>> E=[1,0,0,0; 0,0,1,0; 0,0,0,0; 0,0,0,0];
A=[0,1,0,0; 1,0,0,0; -1,0,0,1; 0,1,1,1];
B=[0; 0; 0; -1]; C=[0,0,1,0]; % 输入各个矩阵
G=dss(A,B,C,0,E); step(G); % 绘制每个状态
G1=zpk(G) % 求出等效零极点模型
图5-2 奇异系统的阶跃响应曲线
如果想用Simulink 对该奇异系统仿真，则应该用Continuous 模块组的Descriptor
State-Space 模块搭建仿真框图，如图5-3 所示。注意，在仿真过程中应该使用刚性微分
方程求解算法，如ode15s，且相对误差限不宜选择太小。这里将相对误差限选择为10?8，
得出的结果与图5-2 完全一致。
·132· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
图5-3 奇异系统的Simulink 模型（文件名：c5mdss.slx）
5.1.4 带有非零初值的传递函数模型
前面介绍了传递函数模型，不过这样的传递函数模型有一个致命的局限性，就
是传递函数是定义在零初值基础上的。如果含有非零的初值，则不易处理。
Simulink Extras 模块集的Additional Linear（附加线性）模块组提供了一些可以
处理初值的模块，如图5-4 所示。不过从实际应用角度看，但只能输入y(t0)，不能处
理带有y′(t0)，y′′(t0) 等初值的传递函数，所以应该考虑建立通用的变换方法，将一
般的含有非零初值的传递函数模型转换成Simulink 可以描述的形式，构造一个能
处理这些初值的模块。本节将通过例子演示这类模型的绘制方法。第7 章还将介绍
构造封装模块的方法，得出通用的模块。
图5-4 Additional Linear 模块组
　　定义5-3 一般线性常系数高阶微分方程模型为
y(n)(t) +a1y(n?1)(t) + a2y(n?2)(t) + · · · + an?1y′(t) + any(t)
= b1u(n?1)(t) + b2u(n?2)(t) · · · + bn?1u′(t) + bnu(t)
（5-1-4）
其中，已知输入信号u(t)，并已知初值y(t0)，y′(t0)，· · · ，y(n?1)(t0)。
注意，在前面的模型中有意假设u(t) 的最高阶次比y(t) 的阶次低1；如果u(t)
的最高阶次与y(t) 一致，则可以事先处理阶次，将其提取出来单独处理；如果u(t)
的最高阶次高于y(t)，则微分方程为物理不可实现的，本书不予进一步考虑。
这样的微分方程用现有的Simulink 模块是没有办法直接描述的，因此应该建
立一个模块描述该微分方程。其实，前面已经试图探讨过该微分方程的Simulink 建
模，例如，图4-26 给出的解决方案使用微分模块，严重影响仿真精度，而图4-28 给出
的解决方案需要事先求出u′(t) 等信号，这在现实应用中要求过于苛刻，因为输入信
号的数学表达式有时不是事先预知的。所以需要一个更好的Simulink 解决方案。
　　定理5-1 如果选择一组状态变量x1(t)，x2(t)，· · · ，xn(t)，则可以得出如下的状
第5 章控制系统的Simulink 建模·133·
态方程模型[49]
8>
>>>><>>>>>:
x′1(t) = x2(t)
x′2(t) = x3(t)
...
x′n?1(t) = xn(t)
x′n(t) = ?anx1(t) ? an?1x2(t) ? · · · ? a1xn(t) + u(t)
y(t) = bnx1(t) + bn?1x2(t) + · · · + b2xn?1(t) + b1xn(t)
（5-1-5）
根据状态方程模型，可以搭建起如图5-5 所示的框图实现[49]。该微分方程可以
通过MATLAB 命令的方式绘制仿真模型。
- Σ - 1
s
- 1
s
-
· · ·
- 1
s
- 1
s
?an
?an?1
?a2
?a1

6

6

6 6
- b1 -

Σ -
b2
bn?1
bn
-
?
-
?
-
? u(t) y(t)
x1(t)
x2(t)
xn?1(t)
xn(t)
x′n(t)
图5-5 线性微分方程的框图表示
由于已知y(t) 及其各阶导数的初值，未知状态变量的初值，所以使用普通的状
态方程模块是不能表示这种传递函数模型的，必须进行适当的底层变换，求出状态
变量的等效初值，然后利用状态方程或框图的结构，获得微分方程的Simulink 仿真
模型，分析这样的微分方程。不过，这样的变换与绘制方法是很烦琐的。另一种方法
是直接获得状态方程模型，然后用State-Space 模块在Simulink 下表示该微分方程
模型。
　　定理5-2 定理5-1 中状态方程的矩阵形式为
A=26664
0 1 0 · · · 0
0 0 1 · · · 0
......
...
... ...
?an ?an?1 ?an?2 · · · ?a1
37775
, B=26664
0
0...
1
37775
, CT=26664
bn
bn?1 ...
b1
37775
（5-1-6）
从给出的状态方程模型看，由已知的传递函数可以直接写出A、B、C 和D矩
·134· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
阵，不过状态变量的初值x(t0) 是不能直接写出的。下面尝试由已知的y(t0)，y′(t0)
等条件求出等效初始状态变量值的方法。
记β1
i = bi，i = 1, 2, · · · , n，则由输出方程可以直接写出
y(t0) = β1n
x1(t0) + β1 n?1x2(t0) + · · · + β1
2xn?1(t0) + β1
1xn(t0) （5-1-7）
对输出方程两边求导，并将状态方程代入，则可以直接写出
y′(t0) = β1n
x′1(t0) + β1 n?1x′2(t0) + · · · + β1
2x′n?1(t0) + β1
1x′n(t0)
= β1n
x2(t0) + β1 n?1x3(t0) + · · · + β1
2xn(t0)+
β1
1??? anx1(t0) ? an?1x2(t0) ? · · · ? a1xn(t0) + u(t0)
 = β2n
x1(t0) + β2 n?1x2(t0) + · · · + β2
1xn(t0) + β2
1u(t0)
（5-1-8）
其中，β2n
= ?β1
1an，β2
i = β1
i+1 ? β1
1ai，i = 1, 2, · · · , n ? 1。
继续求导，则可以推导出
y′′(t0) = ?β2n
x′1(t0) + β2 n?1x′2(t0) + · · · + β2
1x′n(t0) + β1
1u′(t0)
= β3,nx1(t0) + β3,n?1x2(t0) + · · ·+
β3,1xn(t0) + β2,1u(t0) + β1,1u′(t0)
（5-1-9）
其中，β3n
= ?β2
1an，β3
i = β2
i+1 ? β2
1ai，i = 1, 2, · · · , n ? 1。
根据上面的推导可见
y(k?1)(t0) = βkn
x1(t0) + βk n?1x2(t0) + · · · + βk
1 xn(t0)+
βk?2
1 u(k?2)(t0) + · · · + β2
1u′(t0) + β1
1u(t0)
（5-1-10）
其中，βkn
= ?βk?1
1 an，βk
i = βk?1
i+1 ? βk?1
1 ai，i = 1, 2, · · · , n ? 1，k = 2, 3, · · · , n。
　　定理5-3 如果已知输入、输出信号与各阶导数的初值，则可以如下构造出等效
的状态变量初值。
266664
β1n
β1 n?1 · · · β1
1
β2n
β2 n?1 · · · β2
1 ...
...
... ...
βn
n βn
n?1 · · · βn
1
377775
266664
x1(t0)
x2(t0)
...
xn(t0)
377775
=
26666664 y(t0)
y′(t0) ? β1
1u(t0)
...
y(n?1)(t0) ?
n?1 Pi=1
βi1
u(i?1)(t0)
37777775
（5-1-11）
若方程的系数矩阵满秩，则方程（5-1-11）有唯一解，可以计算出等效的状态变
量初值x(t0)。
　　定理5-4 可以证明，系数矩阵是系统的可观测性判定矩阵
266664
β1n
β1 n?1 · · · β1
1
β2n
β2 n?1 · · · β2
1 ...
...
... ...
βn
n βn
n?1 · · · βn
1
377775
=
266664
C
CA
...
CAn?1
377775
（5-1-12）
第5 章控制系统的Simulink 建模·135·
所以，若系统可观测，则可以唯一地确定状态变量的等效初值。
由前面给出的算法，不难编写出如下计算状态变量等效初值的MATLAB 函
数。注意，如果将βk
i 用矩阵表示，比较直接的表示方法即b(k, i) = βk
i 。
function [x0,G1,T,Y]=find_x0(G,y0,u0)
[num,den]=tfdata(G,'v'); n=length(den)-1;
d0=den(1); den=den/d0; num=num/d0; Y=y0(1); D=num(1);
y0=[y0(:); zeros(n-length(y0),1)];
u0=[u0(:); zeros(n-length(u0)-1,1)];
num=num-D*den; num=num(2:end); b=num; a=den(2:end);
A=diag(ones(1,n-1),1); A(end,:)=-den(end:-1:2);
B=[zeros(n-1,1); 1]; C=num(n:-1:1); G1=ss(A,B,C,0);
for k=2:n
b(k,n)=-b(k-1,1)*a(n);
for i=1:n-1, b(k,i)=b(k-1,i+1)-b(k-1,1)*a(i); end
Y=[Y; y0(k)-b(1:k-1,1).'*u0(1:k-1)];
end
T=fliplr(b); x0=inv(T)*Y;
函数的调用格式为[x0,G1,T , ?y0]=find_x0(G,y0,u0)，其中，G为微分方
程等效的传递函数模型；y0 = [y(t0), y′(t0), · · · , y(n?1)(t0)] 与u0 = [u(t0), u′(t0),
· · · , u(n?2)(t0)] 分别为输出信号与输入信号及其各阶导数的初始值构成的向量。返
回的x0 为状态方程模型的等效初值向量；G1 为等效的状态方程模型；T 为方程的
系数矩阵，也是系统的可观测性判定矩阵；?y0 为方程右侧的向量。
　　例5-6 试用简洁的方式为例4-11 建立微分方程的Simulink 仿真模型，并评价其
仿真方法与精度。
　　解例4-11 使用了微分器方法，使得仿真速度慢、精度低。虽然例4-12 引入了简洁
的建模方法，不过方法过于苛刻，因为需要提供输入信号导数的解析表达式。这里，将尝
试使用等效状态方程的仿真方法，用Simulink 求解微分方程。
由给定的微分方程模型，不难写出等效的传递函数为
G(s) =
2s + 1
s4 + 5s3 + 9s2 + 7s + 2
将传递函数模型输入MATLAB 工作空间中，并将输入、输出的向量也输入MATLAB
工作空间中，这时，调用find_x0() 函数，就可以得出等效的状态方程模型与状态
变量的等效初值。
>> num=[2,1]; den=[1,5,9,7,2]; G=tf(num,den);
u0=[0;1;0]; y0=[1;0;0;0];
[x0,G1,T]=find_x0(G,y0,u0) % 计算等效的状态变量初值
·136· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
得出的等效状态方程模型如下，且状态变量的等效初值为x0=[?29, 16,?8, 4]T/3。
x′(t) =264
0 1 0 0
0 0 1 0
0 0 0 1
?2 ?7 ?9 ?5
375
x(t) +264
0
0
0
1
375
u(t), y(t) = [1, 2, 0, 0]x(t)
将这些信息输入给State Space 模块，则可以直接建立如图5-6 所示的Simulink 仿
真模型。对模型进行仿真，并与微分方程解的理论值进行比较
>> [A,B,C,D]=ssdata(G1);
tic, [t1,~,y1]=sim('c5mlin1'); toc, length(t1)
y0=7*exp(-t1)/4-2*exp(-2*t1)/5+t1.*exp(-t1)+...
3*t1.^2.*exp(-t1)/4-sqrt(2)*cos(t1-atan(1/7))/4;
max(abs(y0-y1)), plot(t1,y1,'--',t1,y0)
则可以发现，这时得出的最大误差为1.4211×10?14，计算点数为2142，耗时0.2219 s。可
以看出，这里介绍的仿真方法比前面介绍的方法更易于使用，精度更高。
图5-6 非零初值传递函数的Simulink 模型（文件名：c5mlin1.slx）
第6 章将对这里建立的非零初值传递函数模块进行封装，建立一个可重用的带
有非零初值的传递函数模块。
5.1.5 传递函数矩阵的Simulink 建模
传递函数矩阵是多变量线性控制系统中常用的数学模型。本节将给出传递函
数矩阵的定义，再探讨传递函数矩阵的Simulink 建模方法，并给出一种简洁、通用
的建模方法。
　　定义5-4 连续多变量系统传递函数矩阵的一般形式为
G(s) =26664
g11(s) g12(s) · · · g1p(s)
g21(s) g22(s) · · · g2p(s)
...
...
... ...
gq1(s) gq2(s) · · · gqp(s)
37775
（5-1-13）
其中，gij(s) 定义为第i 路输出信号对第j 路输入信号的放大倍数，称为(i, j) 子传递
函数。该矩阵为q × p 矩阵，表明系统有q 路输出信号，p 路输入信号。
代入系统的输入向量与输出向量，则可以写出
26664 y1(s)
y2(s)
...
yq(s)
37775
=26664
g11(s) g12(s) · · · g1p(s)
g21(s) g22(s) · · · g2p(s)
...
...
... ...
gq1(s) gq2(s) · · · gqp(s)
37775
26664
u1(s)
u2(s)
...
up(s)
37775
（5-1-14）
第5 章控制系统的Simulink 建模·137·
可以看出，如果想从底层刻画系统的传递函数矩阵模型，则是一件很烦琐的
事。下面通过例子演示一个2 × 2 传递函数矩阵的底层建模方法。
　　例5-7 试用Simulink 描述下面的带有时间延迟的多变量传递函数矩阵[50]。
G(s) =2664
0.1134e?0.72s
1.78s2 + 4.48s + 1
0.924
2.07s + 1
0.3378e?0.3s
0.361s2 + 1.09s + 1
?0.318e?1.29s
2.93s + 1
3775
　　解前面介绍了传递函数与延迟环节的Simulink 模块，不过这些模块只能用于单
变量系统的描述，不能直接描述多变量系统。若想在底层描述这里给出的2 × 2 传递函
数矩阵，则需要由式（5-1-14）推导，可以得出下面的显式公式。
y1(s) = g11(s)u1(s) + g12(s)u2(s)
y2(s) = g21(s)u1(s) + g22(s)u2(s)
如果每个子传递函数gij(s) 都用Transfer Fcn 和Transport Delay 模块从底层搭建，
可以建立起如图5-7 所示的Simulink 仿真框图。可以看出，这里介绍的底层建模方法是
比较烦琐的，也是容易出错的，不宜采用这样的方法表示多变量系统。
图5-7 多变量系统的底层建模（文件名：c5mmimo1.slx）
对传递函数矩阵或其他线性系统模型，可以考虑采用MATLAB 控制系统工具
箱支持的LTI（linear time invariant，线性时不变）对象的形式直接输入。MATLAB
支持使用tf() 函数建立传递函数对象，也支持使用ss() 函数建立状态方程甚至
描述符状态方程等对象。
Simulink 提供了控制系统工具箱模块集，如图5-8 所示，其中LTI System 模块
直接表示LTI 模型。可以先用上述函数构造LTI 对象，将建立起来的LTI 对象的变
量名填写到模块的参数对话框中，这样就可以直接使用这些模块了。下面将通过例
子演示基于LTI System 模块的Simulink 建模方法。
　　例5-8 试利用Simulink 的LTI System 模块直接表示例5-7 中的多变量系统模型。
·138· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
图5-8 控制系统工具箱模块集
　　解若想输入例5-7 中的带有延迟的传递函数矩阵模型，一种更直接的方法是先输
入4 个子传递函数模型，然后用普通矩阵的输入方式，直接构造传递函数矩阵模型G，
然后将变量G填写到LTI System 的参数对话框，这样，看起来比较复杂的带有延迟的
传递函数矩阵用这样一个模块就可以表示了，无须建立图5-7 中的底层模型。
>> g11=tf(0.1134,[1.78 4.48 1],'ioDelay',0.72);
g12=tf(0.924,[2.07 1]);
g21=tf(0.3378,[0.361 1.09 1],'ioDelay',0.3);
g22=tf(-0.318,[2.93 1],'ioDelay',1.29);
G=[g11, g12; g21, g22]; % 和矩阵定义一样，这样可以输入传递函数矩阵
还可以用这种方法表示更复杂的模型，包括不能用传递函数模块描述的复合
模型。MATLAB 控制系统工具箱支持带有内部延迟的状态空间模型[51]，还支持带
有可变参数的模型等。这里将通过例子演示复杂线性时不变模型的表示方法。
　　例5-9 考虑例2-7 中给出的复合传递函数模型。该例曾使用Simulink 搭建了复合
传递函数的底层模型，不过方法比较麻烦，也比较容易出错。试用更简单的方法重新建
立其Simulink 模型。
　　解如果用手工转换的方法变换原有的传递函数模型，可以看出，该模型是不能转
换成式（5-1-1）的常规传递函数模型，不过，利用MATLAB 控制系统工具箱提供的模
型处理方法，可以将其直接转换为带有内部延迟的状态方程模型。这样，就可以使用LTI
System 模块直接表示该传递函数模型了。该复合传递函数模型可以通过下面的命令直
接输入MATLAB 工作空间中。然后，将变量G直接填写到LTI System 的编辑框中，这
样，就可以用一个模块直接表示复杂的系统模型了。
>> s=tf('s'); % 定义s 算子
G=(1+3*exp(-s)/(s+1))/(s+2); % 直接输入复合传递函数模型
5.1.6 可变参数线性系统的Simulink 建模
　　定义5-5 线性可变参数状态方程模型的数学表达式为
(x′(t) = A(p)x(t) + B(p)u(t)
y(t) = C(p)x(t) +D(p)u(t)
（5-1-15）
其中，p = [p1, p2, · · · , pm] 为可变参数向量。
在MATLAB 控制系统工具箱中，可以使用realp() 函数定义可变参数。
这种模型在Simulink 下，可以通过如图5-8 所示的控制系统工具箱模块集的
第5 章控制系统的Simulink 建模·139·
Linear Parameter Varying（线性变参数）模块组中的模块直接表示。双击该模块组的
图标，则得出如图5-9 所示的变参数模块与控制器的下级模块组。
图5-9 变参数模块与控制器模块组
5.2 离散线性模型的建模与仿真
前面介绍了连续系统的理论与模型表示方法。所谓连续系统是指系统中信号
在连续变化的系统。如果系统中信号不连续，则系统可能带有离散环节，需要引入
专门的处理方法来建模与仿真。本节先给出离散线性系统的各种数学模型，介绍离
散系统实现模型输入MATLAB 环境的方法，并介绍连续与离散线性模型之间的转
换方法，然后介绍离散线性模型的Simulink 建模方法。最后，还通过例子介绍仿真
控制参数或其他因素对仿真结果的影响。
5.2.1 离散线性系统的数学模型
和连续线性系统类似，离散线性系统同样有自己的差分方程、传递函数与状态
方程表现形式。这里首先给出离散线性系统的数学模型，然后介绍基于MATLAB
控制系统工具箱的模型输入方法。
　　定义5-6 离散线性差分方程的数学形式为
a1y((k + n)T)+ a2y((k + n ? 1)T) + · · · + +an+1y(kT)
= b1u((k+n)T) + b2u((k+n?1)T) + · · · + bn+1u(kT)
（5-2-1）
式中，T 为离散系统的采样周期。更简洁地，简记t + n = (k + n)T，在描述差分方
·140· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
程时可以略去T。
a1y(t + n)+ a2y(t + n ? 1) + · · · + an+1y(t)
= b1u(t + n) + b2u(t + n ? 1) + · · · + bn+1u(t)
（5-2-2）
或更简洁地
a1yt+n + a2yt+n?1 + · · · + an+1yt = b1ut+n + b2ut+n?1 + · · · + bn+1ut （5-2-3）
　　定义5-7 离散线性系统的传递函数模型为
G(z) =
b1zn + b2zn?1 + · · · + bnz + bn+1
a1zn + a2zn?1 + · · · + anz + an+1
（5-2-4）
　　定义5-8 离散线性状态方程模型的数学形式为
(xt+1 = Fxt + Gut
yt = Cxt +Dut
（5-2-5）
MATLAB 的控制系统工具箱提供了tf() 函数和ss() 函数，允许用户直接输
入传递函数与状态方程对象。这些函数的调用格式与连续模型差不多，不同的是，
需要指定采样周期的值。
G=tf(num,den,'Ts',T);
z=tf('z',T);
G=ss(A,B,C,D,'Ts',T);
　　例5-10 如果系统的采样周期为T = 0.1 s，试输入下面两个传递函数模型。
G1(z) =
2z + 1
z4 + 2.6z3 + 2.52z2 + 1.08z + 0.1728
, G2(z) =
2z + 1
(z + 0.6)3(z + 0.8)
　　解从传递函数给出的数学模型看，前者已知分子和分母多项式系数，所以适合于
第一种输入方法。而后者是关于z 的运算，因此，应该先定义z 变换算子，然后用数学表
达式的输入方法输入后者的模型。
>> G1=tf([2 1],[1 2.6 2.52 1.08 0.1728],'Ts',0.1)
z=tf('z',0.1); G2=(2*z+1)/(z+0.6)^3/(z+0.8)
5.2.2 连续离散模型的相互转换
连续线性模型与离散线性模型是可以相互转换的。MATLAB 控制系统工具
箱提供的c2d() 与d2c() 函数可以直接进行这样的转换。如果已知连续LTI 模型
G，且给定采样周期T，则由G1=c2d(G,T) 命令就可以直接获得等效的离散模型
G1。如果G为传递函数模型，则转换后的G1 为离散传递函数；如果G为状态方程模
型，则离散化的G1 为离散状态方程模型。
如果已知离散模型G1，则由G=d2c(G1) 就可以直接获得等效的连续模型G。
在语句调用中无须给出采样周期T，因为该信息已经包含在模型G1 中了。
第5 章控制系统的Simulink 建模·141·
　　例5-11 考虑例5-6 给出的连续传递函数模型，试选择采样周期T = 0.1 s 与T =
0.01 s，对系统进行仿真，并比较离散化系统的阶跃响应。
　　解可以先输入连续系统的传递函数模型，然后在不同采样周期下对系统进行离
散化，这时，可以绘制连续模型与离散模型的阶跃响应曲线，如图5-10 所示。可以看出，
采样周期越小，离散化模型的阶跃响应越接近连续系统模型。此外，可以看出，离散系
统的阶跃响应曲线是以阶梯线的形式给出的。在每个采样周期内，输出信号的值保持不
变。从效果上看，相当于输出信号后面接了一个零阶保持器模块。
>> num=[2,1]; den=[1,5,9,7,2]; G=tf(num,den);
G1=c2d(G,0.1), G2=c2d(G,0.01), step(G,G1,'--',G2,':')
图5-10 连续传递函数离散化的阶跃响应比较
上述语句得出的离散传递函数模型为
G1(z) =
0.000298z3 + 0.0007821z2 ? 0.0007757z ? 0.0002263
z4 ? 3.533z3 + 4.679z2 ? 2.752z + 0.6065
G2(z) =
3.296×10?7z3 + 9.757×10?7z2 ? 9.749×10?7z ? 3.207×10?7
z4 ? 3.95z3 + 5.852z2 ? 3.853z + 0.9512
可以用下面的命令对离散化模型进行连续化
>> G3=d2c(G1)
得出的模型如下所示。可以看出，这里的转换基本上可以还原原始的连续传递函数模
型，只在某些系数上引入微小的误差。
G3(s) =
1.291×10?17s3 ? 3.248×10?16s2 + 2s + 1
s4 + 5s3 + 9s2 + 7s + 2
5.2.3 离散模型的Simulink 建模
前面介绍的Discrete 模块组包含了离散传递函数、离散状态方程等离散系统模
块，也包含了零阶保持器等一般的离散系统模块，可以用这样的模块直接搭建离散
系统的Simulink 仿真模型。
·142· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
由于传递函数模型的天然缺陷，相应的模块只能处理零初值的传递函数模块，
如果需要处理非零初值的传递函数模块，则可以考虑使用Simulink Extras 模块集提
供的Additional Discrete（附加离散）模块组，如图5-11 所示。当然，这样的模块有时
也不利于建模处理，可以仿照前面介绍的连续非零初始条件的传递函数转换方法，
推导出适用于离散传递函数的数学公式与MATLAB 实现方法。
图5-11 Additional Discrete 模块组
除了这两个模块组之外，可以直接使用控制系统工具箱自带模块集中给出的
LTI System 模块，该模块可以直接处理离散的LTI 对象，无须用户将其转换成连续
模型。直接将离散对象赋给LTI System 模块即可直接仿真。
　　例5-12 考虑例5-11 中的离散化问题。试用Simulink 建模并得出连续与离散系统
的阶跃响应。
　　解有两种方法可以描述离散模型，对应的Simulink 仿真模型如图5-12 所示。其
中，分别使用Discrete Transfer Function 和LTI System 模块描述离散环节。为方便建模
与仿真，可用下面的命令写入模型的Model Properties 属性。
num=[2,1]; den=[1,5,9,7,2]; G=tf(num,den);
G1=c2d(G,0.1); [numd,dend]=tfdata(G1,'v');
（a）方法一（文件名：c5mc2d1.slx ） （b）方法二（c5mc2d2.slx）
图5-12 离散化建模与仿真
在Discrete Transfer Fcn 模块中，将变量名numd 和dend 分别填写到分子与分母参
数编辑框中，采样周期填写为0.1。由于这里给定的离散传递函数的参数存于变量中，所
以在模型中不显示传递函数的表达式。在LTI System 模块中，在模型对象编辑框填写
G1，无须填写采样周期的值，因为该值已经包含在G1 对象中了。
对这两个系统进行仿真，得出的仿真结果与图5-10 的完全一致。注意，在图形绘制
时，连续系统的输出信号由曲线表示，而离散模块的输出信号应该由阶梯信号表示。
在一个Simulink 模型中，如果需要使用离散模块，则应该在模块的参数对话框
第5 章控制系统的Simulink 建模·143·
中填写采样周期的值。采样周期可以填写实际的数值，也可以填写?1，表示继承模
块输入信号的采样周期。如果系统中多个离散模块的采样周期是一致的，则没有必
要逐个模块填写采样周期的值，可以统一填写为?1。这样，在描述输入信号或使用
第一个离散模块时填写正确的采样周期即可，其余模块会自动继承这个采样周期，
这使得Simulink 建模更简单。
5.2.4 仿真控制参数对离散系统的影响
一般情况下，如果待仿真的系统只包含离散环节，不包含连续环节，且系统各
个部件的采样周期的值是完全一致的，则应该采用定步长的仿真方法实现仿真过
程。如果系统中含有连续环节，则不能采用离散仿真算法，必须采用变步长仿真方
法，否则，连续模块的仿真精度难以保证。这里将通过例子演示变步长方法中，控制
参数对仿真结果的影响。
　　例5-13 假设连续系统模型由例5-11 给出。选择采样周期T = 0.1 s，则可以得出离
散化模型。试分析仿真控制参数对仿真结果的影响。
　　解考虑例5-11 给出的仿真模型。由于系统中既有离散化模块，也有连续模块，所
以不能采用定步长离散系统的仿真方法，必须采用变步长仿真算法实现系统的仿真。如
果采用默认的控制参数，则可以得出如图5-13 所示的仿真结果。
图5-13 默认参数下的仿真结果
为使得结果更容易观察，这里有意将仿真终止时间设置为2 s。从得出的结果看，尽
管在采样点上仿真结果是正确的，但这显然不是一般期望的离散系统仿真结果，因为正
常情况下，离散系统的输出曲线是阶梯形的。为解决这个问题，应该将Simulink 的默认
误差限设置为很小的值，如前面介绍的3×10?14。如果将绝对与相对误差限设置为这
样的值，则可以直接得出如图5-10 所示的阶梯曲线。从这个例子可以看出，仿真控制参
数对离散系统的仿真结果是有显著影响的。如果想得到正确的结果，应该在允许的条件
下，尽可能设置小的误差限。
·144· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
5.2.5 离散参数精度对仿真结果的影响
如果离散传递函数的极点接近单位圆，则传递函数分母多项式系数发生微小
的变化，则可能影响整个系统的仿真结果。这里将通过例子演示离散系统分母多项
式系数对仿真结果的影响。
　　例5-14 假设连续系统模型由例5-11 给出。选择采样周期T = 0.1 s，则可以得出离
散化模型。试分析离散化模型分母多项式系数小数点保留位数对仿真结果的影响。
　　解可以用下面的语句直接输入连续传递函数模型，并对其离散化。由tfdata()
函数提取离散化模型的分子与分母多项式系数。如果想保留分母多项式系数的小数点
后两位数字，则可以将参数先乘以100 后取整，再乘以0.01，这样就可以舍去多余的位，
并将变换后的分母多项式赋给G1a。通过这样的方法可以保留分母多项式不同的位数，
然后比较模型的阶跃响应曲线，如图5-14 所示，这里，只给出保留4 位小数与保留6 位
小数的离散化系统，因为G1a 与G1b 是不稳定的。可以看出，保留分母多项式小数点后
不同的位数会对仿真结果有巨大的影响。位数保留过少则可能导致系统不稳定，保留4
位（系数具有万分之一误差）则得出的响应明显是错误的，稳态值也有显著差异；保留6
位也有比较明显的偏差。这在使用离散系统传递函数仿真的时候应该格外注意。
>> num=[2,1]; den=[1,5,9,7,2]; G=tf(num,den);
G1=c2d(G,0.1); [n,d]=tfdata(G1,'v');
G1a=G1; G1b=G1; G1c=G1; G1d=G1;
d0=0.01*round(d*100); G1a.den{1}=d0;
d0=0.001*round(d*1000); G1b.den{1}=d0;
d0=0.0001*round(d*10000); G1c.den{1}=d0;
d0=1e-6*round(d*1e6); G1d.den{1}=d0; step(G,G1c,G1d)
保留小数点后4 位的错误结果
连续系统响应
保留小数点后6 位的偏差结果
图5-14 离散化的仿真结果比较
分析离散化模型G1，可以发现，在0.9048 处有三重极点，比较接近单位圆，另
一个极点位于0.8781。如果保留离散化模型小数点后两位，则模型G1a 的极点位置为
第5 章控制系统的Simulink 建模·145·
1.0931 ± 0.2490j，0.6719 ± 0.1841j，前一对极点的幅值为1.1211，位于单位圆外，所以离
散化的模型是不稳定的。如果保留小数点后4 位有效数字，则G1c 模型的极点为0.9759，
0.8858 ± 0.0796j，0.7857，第一个极点离单位圆比较近。可以看出，该极点与离散化模型
G1 有巨大差距，所以仿真曲线截然不同。
>> eig(G1), d1=max(abs(ans))
eig(G1a), d2=max(abs(ans)), eig(G1c), d3=max(abs(ans))
5.3 非线性环节的Simulink构造
从Discontinuities 模块组现有模块可见，该模块组提供的非线性模块是很有限
的，而例2-8 给出的非线性模块等效变换方法比较麻烦，可表示的非线性模块也是
很有限的，因此需要更一般的非线性模块描述方法。本节将介绍利用查表模块构造
任意单值非线性静态环节的方法，并介绍基于查表与开关模块构造双值非线性静
态模块搭建方法。
5.3.1 查表模块
双击图2-36 中给出的Lookup Tables 模块组中的1-D Lookup Table 模块图标，
则打开的参数对话框如图5-15 所示。其中，Number of table dimensions（表格的维
数）是查表维数，这里给出的是一维，如果给出其他数值，则该模块可以同样用于二
维甚至高维查表功能。
图5-15 一维查表模块的参数对话框
该模块有两个重要的参数向量，Breakpoints 1（一维转折点）和Table data（表
格数据），其值分别为转折点的横坐标向量u = [u1, u2, · · · , um] 与纵坐标向量y =
[y1, y2, · · · , ym]，它们的值需要用户自行填写。有了这两个向量，就可以由模块的输
·146· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
入信号唯一地确定模块的输出信号。在实际应用中，应该事先给定这两个向量。
默认情况下，若模块输入信号满足uk ? u(t) < uk+1，则输出信号y(t) 为
y(t) =
yk+1 ? yk
uk+1 ? uk
(u(t) ? uk) （5-3-1）
这种计算方法实际上是线性插值。如果用户想用其他算法计算输出信号，则可
以选择图5-15 中的Algorithm（算法）标签，在得到的对话框中再设置Interpolation
（内插插值）和Extrapolation（外插）算法。
如果想编辑这些转折点的坐标，还可以单击Edit table and breakpoints（编辑表
格与转折点）按钮，弹出如图5-16 所示的表格编辑界面。用户可以选中横坐标表格
项修改横坐标的值，选中纵坐标表格项则可以修改纵坐标的值，用这种可视的方法
更方便地修改转折点参数。
图5-16 表格参数编辑对话框
模块参数设定完成后，如果横坐标、纵坐标向量的长度一致，则Simulink 机制
会根据输入的数值将函数的形状自动绘制出来，作为该模块的图标。
5.3.2 单值非线性模块的搭建
为叙述方便起见，这里将静态非线性函数分为单值非线性函数与双值非线性
函数分别进行探讨，并介绍Simulink 的一般建模方法。
　　定义5-9 单值非线性函数是指输入信号值u(t) 在每个t 时刻对应唯一输出信
号值y(t) 的函数，其数学表达式为y(t) = f(u(t))。
现在考虑分段线性的静态非线性函数模块的建模方法。单值非线性静态模块
可以由一维查表模块构造出来。考虑如图5-17 所示的分段线性非线性静态特性，已
知非线性特性的转折点为(x1, y1)，(x2, y2)，· · · , (xn?1, yn?1), (xn, yn)。如果想用
第5 章控制系统的Simulink 建模·147·
Simulink 的查表模块表示此非线性模块，则需要在x1 点之前任意选择一个x0 点，
即x0 < x1，这样可以根据非线性函数本身求出该点对应的y0 值。同样还应该任意
选择一个xn+1 点，使得xn+1 > xn，并根据延长线求出yn+1 的值，这样就可以构造
两个向量xx 和yy，使得
xx=[x0,x1,x2,· · · ,xn,xn+1]; yy=[y0,y1,y2,· · · ,yn,yn+1];
-
6
x
y
x1 x2 x3
y1
y2
y3
xn?1 xn
xn?1
yn
图5-17 分段线性单值非线性模块的一般形式
将xx 填写到Lookup Tables 模块的Breakpoints 1 编辑框中，并将yy 变量名填
写到Table data 编辑框中，则可以使用查表模块在Simulink 中直接描述非线性特
性。下面将通过例子描述单值非线性模块的Simulink 建模方法。
　　例5-15 试用Simulink 的查表模块重新描述例2-8 的非线性特性。
　　解由图2-30（a）给出的非线性环节可见，现有的转折点为(?1.2, 1)，(?0.2, 0)，
(0.2, 0) 和(1.2, 1)。如果想用查表模块描述该非线性特性，则应该在左侧和右侧的延长
线上各加入一个点，例如，添加(?1.5, 1) 和(2, 1)，这样，就可以分别构造出两个向量，
xx=[-1.5,-1.2,-0.2,1.2,2]，yy=[?1,?1,0,0,1,1]。将这两个向量输入MATLAB
工作空间中，再将这两个变量填写到Lookup Tables 模块的参数对话框，则该模块可以
反映所需的非线性特性。可以看出，这里给出的方法比例2-8 介绍的方法更直观，其适
用性也更广，可以用这样的方法描述任意复杂的单值静态非线性环节。
5.3.3 双值非线性模块的搭建
双值非线性函数是指同一个u(t) 在不同的外部条件（例如，u(t) 的走行方向）
下，对应的输出信号y(t) 可能取不同值的函数。双值非线性模块的构造就没有这么
简单了，这里用简单例子来演示如何对双值非线性静态环节进行Simulink 建模，并
总结一般的建模方法。
由前面的叙述可以看出，任何单值非线性函数均可以采取该方式来建立或近
似，但如果非线性中存在回环或双值属性，则简单地采用这样的方法是不能构造的，
解决这类问题则需要使用开关模块。
·148· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
　　例5-16 试用Simulink 搭建如图5-18 所示的双值非线性特性。
-
6
-
-
?
 6

(2, 0)
(2, 1)
(1, 0)
图5-18 带有回环的继电非线性环节
　　解可以看出，该特性不是单值的，该模块中输入在增加时走一条折线，减小时走
另一条折线。将这个非线性函数分解成如图5-19 所示的两个单值函数，当然这两个单
值函数是有条件的，它们是区分输入信号上升还是下降的。
-
6
-
6
(1, 0)
(1, 1)
(?2, 0)
(?2,?1)
(2, 1)
(2, 0)
(?1, 0)
(?1,?1)
（a）当输入量增加时（b）当输入量减小时
图5-19 回环函数分解为单值函数
建立双值非线性环节的一个关键问题是如何判断输入信号是上升信号还是下降
信号。如果能识别出输入信号的走向，则可以通过开关模块控制上分支和下分支，使
得它们能按要求导通。Simulink 的连续模块组中提供了一个Memory 模块，该模块记忆
前一个计算步长上的信号值，所以可以按照图5-20 中所示的格式构造一个Simulink 模
型。在该框图中使用了一个比较符号来比较当前的输入信号与上一步输入信号的大小，
u(t) 上升分支
u(t) 下降分支
图5-20 非线性模块的Simulink 表示（文件名：c5mloop1.slx）
其输出是逻辑变量，在上升时输出的值为1，下降时输出的值为0。由该信号可以控制后
面的开关模块，设开关模块参数对话框的Threshold（阈值）为0.5，则当输入信号为上升
时由上面的通路计算整个系统的输出，而下降时由下面的通路计算输出。
第5 章控制系统的Simulink 建模·149·
两个查表模块的输入输出转折点向量分别为
x1 = [?3,?1,?1 + ?, 2, 2 + ?, 3], y1 = [?1,?1, 0, 0, 1, 1]
x2 = [?3,?2,?2 + ?, 1, 1 + ?, 3], y2 = [?1,?1, 0, 0, 1, 1]
其中，? 可以取一个很小的数值。例如，可以取10?10。注意，该值不宜选择为过小的值，
如MATLAB 保留的常数eps，否则不能正确绘制图标。
　　例5-17 试建立如图5-21 所示的带有饱和的继电非线性环节。
-
6
(2, 0)
(2, 1) (3, 1)
(1, 0)  -
-

图5-21 带有饱和的继电非线性环节
　　解从本质上看，这里的非线性环节与前面介绍的是完全一致的，都需要将回环非
线性分解成两个单值的非线性环节。不同的是，两个分解的单值函数与前面的单值函
数转折点不一致，需要针对新的非线性特性重新提取。仍可以利用前面建立的Simulink
模型，只需修改两个查表函数参数对话框的向量即可。
x1 = [?3,?2,?1, 2, 3, 4]，y1 = [?1,?1, 0, 0, 1, 1]
x2 = [?4,?3,?2, 1, 2, 3]，y2 = [?1,?1, 0, 0, 1, 1]
从而就能得出整个系统的Simulink 仿真框图，如图5-22 所示。
u(t) 上升分支
u(t) 下降分支
图5-22 非线性模块的Simulink 表示（文件名：c5mloop2.slx）
从前述的分析结果可以看出，任意的非线性静态环节，无论是单值非线性还是
双值非线性，均可以用类似的方法用Simulink 搭建起模块，直接用于仿真。
5.3.4 多维查表模块
Simulink 不但提供了一维查表的模块，还提供了二维和多维查表模块。提供的
多维查表可以使用线性插值的算法求出输出。下面以二维查表模块为例演示多维
查表的应用。
·150· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
　　例5-18 考虑二维函数z = f(x, y) = (x2 ? 2x)e?x2?y2?xy。假设能构造较稀疏的
x-y 平面上的网格点，并可以计算出在这些点上的高度值：
>> xx=linspace(-3,3,15); yy=linspace(-2,2,15);
[x,y]=meshgrid(xx,yy); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
surf(xx,yy,z);
其样本点对应的三维曲面如图5-23 所示。试根据这些样本点构造Simulink 仿真模型。
图5-23 样本点数据的三维曲面图
　　解已知这些点的数据xx、yy 和z，其中，xx 与yy 分别为x 轴与y 轴的刻度向量，z
为二元函数的取值矩阵。用两个Random Number（随机数发生器）模块作为输入信号，
后面接查表模块。现在试图用均匀分布随机数的方式生成一组x和一组y 变量作为输
入信号，输入二维查表模块，就可以构造出如图5-24（a）所示的Simulink 框图。注意，在
两个随机数发生模块中，伪随机数的种子不能选择相同的值，否则两路输入信号将是一
致的，不能很好演示这里的例子。可以将第二路随机数种子选择为不同的值，如123456，
另外，假设两个随机数发生模块的范围分别为(?3, 3) 和(?2, 2)。
双击二维查表模块图标打开如图5-24（b）所示的对话框，读者可以将这3 个样本点
变量分别填写到相应的位置，从而正确地建立二维查表模块。
在仿真中如果选择定步长为0.001 s，并设定终止仿真时间为10，则可以用仿真的
（a）Simulink 框图（文件名：c5m2dtab.slx） （b）查表模块参数对话框
图5-24 二维查表模块设置与仿真
第5 章控制系统的Simulink 建模·151·
方法计算出10000 个点，利用这些数据，用描点的方式就可以绘制输出信号的三维图，
如图5-25 所示。可以看出其形状类似于已知数据给出的三维曲面。
>> plot3(yout(:,2),yout(:,1),yout(:,3),'.')
图5-25 仿真结果的散点图
默认设置下，这个模块采用的是线性插值的方法。如果读者需要不同的插值算
法，则可以选择Algorithm 标签，在得出的对话框中选择其他插值算法。除了二维表
格模块之外，Simulink 还允许进行多维查表运算。
5.4 模块运行的自动排序
Simulink 模型是由模块构成的，所以整个仿真过程是由这些模块每步依次运
行而推演出来的。这里提及“依次”，是指Simulink 模型运行时，模块的运行应该有
一个次序。这个次序是如何决定的呢？
Simulink 模型建立完成后，其运行机制就会自动地为模块的运行进行排序预
处理。一般情况下，不依赖自身输入端信号计算其输出信号的模块会自动排在前面
执行，而其他的模块按照连接的关系接着执行。不依赖自身输入端信号的模块主要
包括信号源模块和积分器模块。信号源是没有输入端子的模块，依赖其本身的参数
和内在的时间，就可以直接计算出模块的输出信号；对积分器模块而言，在一个仿
真步长内，由积分器的初值就可以直接计算出模块的输出信号，而整个模型运行起
来以后，积分器的输入信号就可以依次计算出来了。积分器的输入信号会更新积分
器的初值，以备下一步仿真之用。
从模块的自动排序而言，如果信号源与积分器的输出信号已经计算出来了，则
可以利用这些模块直接计算其连接的其他模块，从而构造出所需的仿真次序。
由传递函数或状态方程模型看，如果传递函数为正则的，或状态方程的D = 0，
则可以依赖其底层积分器的初值先计算输出信号，然后用随后计算出来的输入信
·152· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
号更新底层积分器的值，以备下一步仿真，所以这类模块也应该自动排在前面。如
果系统传递函数的分子阶次与分母阶次相等，或状态方程中D ?= 0，则计算模块的
输出信号中，有一部分是需要已知输入信号，然后叠加到模块的输出信号上的，所
以这时的模块就不能排在前面，只有等其输入信号计算出来后才能计算模块的输
出信号。这种输入信号在模块输出信号中直接反映出来的模块又称为直馈（direct
feedthrough）模块，其DirectFeedthrough 属性值为1。模块的属性名和属性值的
问题后面还将详细介绍。
Simulink 模型窗口的快捷菜单Other Displays-Blocks-Sorted Execution Order（
模块执行排序）如图5-26 所示，允许用户显示模型中的模块执行顺序。下面将
通过例子演示该菜单项，并解释排序的情况。
图5-26 显示模块执行顺序的快捷菜单
　　例5-19 试显示例4-23 中给出的延迟微分方程组模块执行顺序。
　　解打开c4mdde1.slx 文件，将在模型窗口中打开该仿真模型。选择图5-26 中给出
的快捷菜单项，则在模型中会自动显示模块的执行顺序，如图5-27 所示。注意，排序后
图5-27 模型中模块的执行顺序
模块的背景颜色被自动变换成浅色，但个别模块（如这里的Fcn 模块）的前景颜色不会
发生变化。
从给出的顺序显示看，因为Transfer Fcn 模块是正则的传递函数，模型首先计算出
其输出信号，该信号直接馈入Out1 模块，所以，输出端子排在第2 位。另外，Step 模块是
信号源，因此被自动排在第3 位。由于信号1 是已知的，因此Transport Delay1 模块可
第5 章控制系统的Simulink 建模·153·
以直接计算，被排在第4 位。这时，暂时不能计算Add 模块的输出，因为其输入信号不
全是已知的。现在考虑积分器模块，由于其初值是已知的，所以其信号可以直接流入
Transport Delay、Fcn 模块，也可以同时流入Gain 模块。当然，Integrator 自身也可以计
算。因此，这些模块排在第5 ～ 8 位。现在，Add 模块的输入各路信号都是已知的了，由
此，可以计算其输出信号。这样，Add 模块就排在最后了。
在实际计算中，第4、5、7 与8 位的实际计算次序不会影响仿真的进程，第1 和3 位
交换次序也不会影响仿真的进程。但很多模块的相对顺序还是很重要的，例如，模块1
已知才能计算模块2；模块5 已知才能计算模块6；模块3、4、6 和8 已知才能计算模块9，
这些模块的相对次序必须正确排列，否则仿真的结果将值得怀疑。Simulink 的机制可以
生成可靠的自动排序结果，得出正确的仿真结果。
　　例5-20 试对例4-17 中给出的隐式Simulink 仿真模型进行排序分析。
　　解打开该方程的仿真模型c4mimp.slx，并选择前面提及的快捷菜单项，则将显示
出如图5-28 所示的自动排序结果。
图5-28 隐式微分方程的模块执行顺序
从给出的模型排序结果看，积分器、时钟信号应该排在最前面。如果右侧的积分器
已经计算出来，则可以直接计算出输出端子的信号，因此该模块也排在前面。这些模块
的排序没有关系，只要右侧积分器排到输出端子之前就行了。
现在看Fcn 模块与Algebraic Constraint 模块，这两个模块的顺序关系比较奇怪：
（1）如果想先计算Fcn 模块的输出，但Fcn 模块的输入端包含了u(3)，而u(3) 是
Algebraic Constraint 模块的输出，所以，这个计算顺序不行。
（2）如果想先计算Algebraic Constraint 模块的输出，但其输入端为Fcn 的输出，因
此这个次序也是不行的。
这就出现了一个悖论。先计算哪个都不行，那就同时计算吧。这些相关的信号满足
代数方程，因此这个现象又称为代数环。后面还将专门探讨代数环的问题。如果想正确
地求出这两个信号，则在每步仿真中需要求解代数方程。从排序上看，Simulink 将这两
个模块排为5.1 和5.2 是合理的。
·154· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
5.5 非线性控制系统的近似方法
在传统的非线性控制理论中，由于没有合适的仿真工具支持，所以分析非线性
系统的手段受到了极大的限制。人们通常采用线性化方法，将非线性系统在某个
小的范围内近似成线性系统。此外，还经常采用描述函数（describing function）方
法对某些非线性环节直接近似，得出近似的仿真结果。如果一个系统带有延迟环
节，或含有无理环节，也可以使用Padé 近似技术将其近似为不含有延迟的线性系
统，对其进行近似分析。本节介绍非线性系统的近似仿真方法，如延迟模型的Padé
近似方法与一般非线性系统的线性化方法等，还介绍代数环的概念与代数环消除
的近似方法等。
值得指出的是，这些近似方法通常可能带来很大的误差。这样近似的目的是期
望套用成熟的线性系统分析工具对其进行理论分析。如果不追求这样的理论分析，
当然，有了Simulink 这样的强有力仿真工具，则可以对复杂的非线性系统进行精确
的仿真分析，并评价传统近似方法的求解精度。
5.5.1 延迟环节的Padé 近似
对连续系统而言，纯延迟环节的数学形式为e?Ts，其中，T 为延迟时间常数。延
迟环节是不能用线性环节直接表示的，所以在线性化过程中，经常采用Padé 近似
对其进行逼近。Padé 近似是以法国数学家Henri Padé（1863?1953）命名的，其最原
始的目标是用有理式去逼近任意函数，而s 的有理式正巧是线性系统需要的传递函
数模型。
对一般函数f(s) 进行Taylor 幂级数展开，则可以得出
f(s) = c1 + c2s + c3s2 + · · · + cmsm?1 + · · · （5-5-1）
而利用Padé 有理式近似，则可以得出
Gr/k(s) =
β1 + β2s + · · · + βrsr?1 + βr+1sr
α1 + α2s + · · · + αksk?1 + αk+1sk
（5-5-2）
其中，r 和k 分别为用户选定的分子和分母阶次。如果想让有理式保留原函数幂
级数展开式的前r + k + 1 个系数，根据适当运算，可以推导出计算βi、αi 的计
算式，从而得出函数的Padé 近似。MATLAB 控制系统工具箱提供了pade() 函数，
计算延迟环节e?Ts 的Padé 近似表达式Gr/r(s)，G=pade(T,r)。本书作者编写了
任意函数、任意阶次组合的Padé 近似表达式求取函数pade_app()，其调用格式
为G=pade_app(c,r,k)，其中，c 是函数Taylor 级数系数向量[52]。
一般地，如果G模型为LTI 对象，即使模型含有复杂的内部延迟，仍可以使用
第5 章控制系统的Simulink 建模·155·
控制系统工具箱通过的G1=pade(G,r) 函数变换出近似模型，其中，G1 模型不含
有延迟环节。
　　例5-21 考虑例5-7 中给出的多变量传递函数矩阵模型。试通过仿真比较对延迟环
节作2 阶Padé 近似后的误差。
　　解MATLAB 提供的pade() 不但能对纯延迟环节求取Padé 近表达式，还可以对
传递函数矩阵做相应的近似。因此由下面的语句可以直接绘制出原始延迟传递函数近
似与Padé 近似模型的阶跃响应曲线，如图5-29 所示。可以看出，除了g22(s) 模型之外，
逼近效果是很理想的，几乎看不出二者的区别。g22(s) 自身的延迟时间常数相对较大，
在初始时刻稍有误差，不过一般是可以接受的。
>> g11=tf(0.1134,[1.78 4.48 1],'ioDelay',0.72);
g12=tf(0.924,[2.07 1]);
g21=tf(0.3378,[0.361 1.09 1],'ioDelay',0.3);
g22=tf(-0.318,[2.93 1],'ioDelay',1.29);
G=[g11, g12; g21, g22];
step(G,pade(G,2),'--') % 比较原系统与近似模型的阶跃响应
图5-29 Padé 近似的比较
5.5.2 Simulink 模型的线性化
在控制理论中，线性系统分析与设计工具是很完备的，可以通过不同的方法，
例如，时域方法与频域方法等分析线性系统的各种性质。相比之下，除了仿真方法
之外，只有很少的分析方法可以用于一般的非线性系统。因此，可以考虑在某些固
定的工作点（或称平衡点）附近对非线性系统进行线性化近似。在介绍非线性系统
线性化之前，首先介绍系统平衡点的概念与求解方法。
　　定义5-10 如果非线性系统状态方程为
x′(t) = f??x(t)
 （5-5-3）
·156· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
与t 无显式关系，则系统称为自治系统。
　　定义5-11 如果定义5-10 中自治系统的状态变量的导数x′(t) 理解成该状态变
量的速度，或理解成质点在各个坐标轴上的速度分量。如果速度的分量都等于0，则
状态方程描述的质点就会停止下来。满足状态变量导数等于零的点称为系统的平
衡点。
　　定理5-5 自治状态方程x′(t) = f??x(t)
的平衡点可以由下面的代数方程直
接求出。
f??x(t)
= 0 （5-5-4）
这里给出的平衡点搜索方法只局限于自治微分方程，对非自治微分方程而言，
还需要考虑时间变量t，由于方程的个数小于未知数的个数，因此平衡点是不能求
取的。
如果有Simulink 模型需要得出其平衡点，读者应该首先在某些输入端连接输
入端子，在输出端连接输出端子。做了这些准备工作之后，则可以调用trim() 函数
直接计算平衡点：
[x0,u0,y,Δx]=trim(模型名)
[x0,u0,y,Δx]=trim(模型名,x0,u0) % 指定平衡点的搜索初值
其中，返回的x0 与u0 为平衡点处的状态变量和输入信号的值；y 为平衡点处的输
出信号；Δx为平衡点处状态变量的增量，应该是接近零的值。在后一种调用格式
中，允许用户指定状态变量和输入信号的搜索初值。
有了微分方程的平衡点x0 和u0，可以对微分方程进行线性化近似。本节介绍
微分方程的线性化并给出MATLAB 求解方法。
首先，考虑一个简单的二维自治非线性系统。
( x′(t) = F??x(t), y(t)
 y′(t) = G??x(t), y(t)
 （5-5-5）
已知平衡点x0 和y0。在平衡点附近选择一个点(x0+Δx, y0+Δy)，假设|Δx|?1，
|Δy| ? 1，则可以写出
8>
><>>:
F??x0 + Δx, y0 + Δy
=
?
?x
F(x0, y0)Δx +
?
?y
F(x0, y0)Δy + · · ·
G??x0 + Δx, y0 + Δy
=
?
?x
G(x0, y0)Δx +
?
?y
G(x0, y0)Δy + · · ·
（5-5-6）
其中，略去了Δx 和Δy 的高次项，并略去了F(x0, y0) 与G(x0, y0)，因为平衡点处这
第5 章控制系统的Simulink 建模·157·
两个值均为0。由式（5-5-6）不难写出
8>
><>>:
Δx′(t) ≈
?
?x
F(x0, y0)Δx(t) +
?
?y
F(x0, y0)Δy(t)
Δy′(t) ≈
?
?x
G(x0, y0)Δx(t) +
?
?y
G(x0, y0)Δy(t)
（5-5-7）
令z1(t) = Δx(t)，z2(t) = Δy(t)，则可以写出关于z(t) 的线性状态方程。
z′(t) = Jz(t) （5-5-8）
其中，J 为Jacobi 矩阵，其数学形式为
J =2664
?
?x
F(x0, y0)
?
?y
F(x0, y0)
?
?x
G(x0, y0)
?
?y
G(x0, y0)
3775
（5-5-9）
　　定理5-6 更一般地，x′(t) = f??x(t)
给出的非线性一阶显式微分方程组可以
线性化为式（5-5-8）的线性形式，其中
J =266664
?f1/?x1 ?f1/?x2 · · · ?f1/?xn
?f2/?x1 ?f2/?x2 · · · ?f2/?xn ...
...
... ...
?fn/?x1 ?fn/?x2 · · · ?fn/?xn
377775
(
x1(0),x2(0),·
·
·
,xn(0))
（5-5-10）
微分方程的线性化可以通过符号运算实现，可以直接使用符号运算工具箱中
的jacobian() 函数计算Jacobi 矩阵。
如果非线性系统由Simulink 模型描述，且指定了输入与输出端子，则可以利
用G=linearize(模型名) 函数直接计算系统的线性化模型G，其中，G为LTI 对
象。如果系统中不含有非线性环节，则可以直接对Simulink 描述的线性模型进行化
简，得出系统的LTI 对象。
还可以用linmod2() 等底层函数直接提取Simulink 模型的线性化状态方程矩
阵A、B、C 和D。
[A,B,C,D]=linmod2(模型名,x0,u0)
如果系统模型只由线性模块构成，线性化实际上就是线性模型的等效化简，没
有模型近似。在线性化语句调用时，可以忽略掉输入变元x0 和u0，因为线性系统的
等效模型提取与平衡点无关。
对一般连续非线性Simulink 模型而言，推荐使用linmod2() 函数。如果系统
中含有延迟环节，则建议使用linmod() 函数进行线性化计算。该函数使用Padé 近
似去逼近延迟环节；如果系统模型中含有离散环节，则建议使用dlinmod() 函数，
获得离散系统的状态方程模型。
·158· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
　　例5-22 考虑例3-1 中给出的非线性模型，试对其进行线性化近似，并比较线性化
模型的近似效果。
　　解假设饱和非线性的宽度为2，死区非线性的宽度为0.1，则可以建立如图5-30 所
示的仿真模型。与图3-8 给出的模型不同的是，因为线性化过程需要模型含有输入端子，
所以这里为该模型的输入端添加了一个输入端子。这样，由下面的命令可以直接对原非
线性系统线性化近似。
>> G=linearize('c5mlinr1'), zpk(G)
图5-30 非线性系统的仿真模型（文件名：c5mlinr1.slx）
得出的线性化模型是以状态方程形式给出的，其数学形式从略。可以将其转换为零
极点增益模型为
G(s) =
(s + 0.3)
(s + 0.2185)(s2 + 0.4498s + 0.4221)(s2 + 3.332s + 3.253)
如果想评价线性化模型对原系统的逼近效果，则可以先对原系统进行仿真分析，然
后给出下面的命令。得出的系统响应对比如图5-31 所示。可以看出，这样的线性化模型
与原来的非线性系统偏差较大。
非线性系统
线性化模型
图5-31 非线性系统的阶跃响应逼近
>> [y1 t1]=step(G,25); plot(tout,yout,t1,y1,'--')
如果非线性系统死区的宽度变成±0.3，则得出的线性化模型没有变化，但两个系
统的阶跃响应曲线如图5-32 所示。可以看出，随着死区变宽（更确切地说，是非线性程
度增高），线性化模型对原系统的逼近效果变得更差，不能很好地描述原来的性质。
第5 章控制系统的Simulink 建模·159·
非线性系统
线性化模型
图5-32 死区变宽后非线性系统的阶跃响应逼近
从上面的例子可以看出，线性化逼近会带来较大的仿真误差，这种误差和原始
模型的非线性程度有关。非线性程度越高，误差越大，甚至可能由误差得出误导性
的结果。所以单纯从数值仿真角度看，实际应用中应该慎用或不用非线性系统的线
性化技术。
　　例5-23 考虑例2-7 中给出的复合传递函数模型，试利用线性化方法得出其线性近
似模型，并比较线性化模型的阶跃响应近似效果。
　　解因为需要线性化，所以应该使用输入、输出端子描述系统的输入与输出信号；
因为需要求取系统的阶跃响应，所以需要用阶跃输入信号激励系统。考虑到这两方面因
素，可以建立起如图5-33 所示的仿真模型。
图5-33 复合系统的仿真模型（文件名：c5mlinr2.slx）
如果想对含有延迟环节的模型做线性化处理，应该双击延迟模块，在打开的对话框
下部可以填写Pade order（Padé 近似的阶次），如分别选择为2 和3，这样，由下面语句
可以得出线性化模型，并绘制原模型与线性化模型的阶跃响应比较，如图5-34 所示，其
中，实线为延迟模型，虚线为三阶Padé 近似模型，点线为二阶Padé 近似模型。可以看
出，Simulink 仿真模型的线性化近似还是比较精确的。此外，可以看出，三阶近似比二阶
近似模型的逼近效果更精确一些。
>> G=linearize('c5mlinr2'); G1=zpk(G)
[y1 t1]=step(G); plot(tout,yout,t1,y1,'--')
得出二阶和三阶Padé 近似的线性化传递函数模型如下，可见，这两个模型都不含
有延迟。
G2(s) =
(s + 10.44)(s2 ? 0.4404s + 4.598)
(s + 2)(s + 1)(s2 + 6s + 12)
·160· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
图5-34 原模型与线性化模型比较
G3(s) =
(s2 ? 0.4289s + 4.443)(s2 + 10.43s + 108)
(s + 4.644)(s + 2)(s + 1)(s2 + 7.356s + 25.84)
事实上，对这样的线性系统而言，即使不使用线性化方法，也可以得出同样的线性
化模型，得出的结果是完全一致的。
>> s=tf('s'); G=(1+3*exp(-s)/(s+1))/(s+2);
G1=zpk(pade(G,2)), G2=zpk(pade(G,3))
5.5.3 代数环现象与代数环消除
例5-20 已经演示了代数环的现象。代数环是指在一个仿真模型中，某个或某些
模块的输出信号取决于其输入信号，而输入信号同时又取决于其输出信号的现象。
这需要在每步仿真过程中，求解一次代数方程。由于代数方程的求解是比较耗时的，
带有代数环的仿真模型势必增加仿真过程的计算量。
在如图5-35 所示的仿真参数设置对话框中，如果选择了Diagnostics（诊断）栏
目，则第一项为Algebraic loop（代数环）。该项有3 个选项：none（不处理）、warning
（警告）与error（错误），默认的是warning 选项。如果仿真模型中有代数环，则给出警
告信息，但仿真过程继续执行，建议保持这个选项；如果设置为error，一旦发现代数
环则仿真过程终止，给出错误信息，这样的方式是不可取的；none 选项会忽略代数
环的存在，直接求解代数方程，完成仿真过程，不给出任何提示。
下面首先通过一个例子来演示Simulink 仿真的代数环现象及问题，然后探讨
消除代数环的几种方法，并给出一种有效的方法。
　　例5-24 考虑图5-36（a）中给出的一个简单的线性反馈系统模型。试建立Simulink
仿真模型，并观察其中的代数环。
　　解对该模型而言，因为开环模块的分子阶次与分母阶次是相同的，所以计算输出
信号y(t)，首先需要已知模块的输入信号u(t)，而想计算u(t) = r(t) ? y(t)，又需要事先
已知y(t)，这就出现了一个怪圈。这个怪圈就是所谓的代数环现象。在实际仿真中，如果
第5 章控制系统的Simulink 建模·161·
图5-35 诊断对话框
- - s + 3
s + 1
-
6
r(t) u(t) y(t)
?
（a）反馈系统框图（b）仿真模型（文件名：c5malg1.slx）
图5-36 简单反馈系统模型
遇到代数环，则需要在每步仿真过程中求解一次代数方程，解出满足代数方程的u(t) 和
y(t) 信号。
由相应的模块可以绘制出如图5-36（b）所示的Simulink 仿真框图。对这样的系统
进行仿真，将在命令窗口中给出如下警告信息：
Warning: Block diagram 'c5malg1' contains 1 algebraic loop(s).
To see more details about the loops use the command
Simulink.BlockDiagram.getAlgebraicLoops('c5malg1') or the
command line Simulink debugger by typing sldebug('c5malg1')
in the MATLAB command window. To eliminate this message,
set Algebraic loop to "none".
说明该模型有一个代数环，可以尝试不同的命令分析代数环。如果不想显示代数环提示
信息，也可以将图5-35 中的Algebraic loop 警告设置为none。
　　例5-25 试用控制理论推导的方法消除例5-24 中系统的代数环。
　　解其实，如果利用负反馈系统等效模型计算公式，就可以把这个回路的等效传递
函数直接计算出来
G1(s) =
G(s)
1 + G(s)
=
(s + 3)/(s + 1)
1 + (s + 3)/(s + 1)
=
s + 3
2s + 4
·162· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
这样，用Step 模块直接驱动G1(s) 模块就可以彻底消除代数环了。
虽然这里给出的方法完美地避免了代数环，不过这种方法的适用性不强，因为
很多模型是不能这样化简的。例如，如果回路中含有非线性环节，则这样的等效变
换方法是不可行的，必须探讨其他的通用方法。
代数环之所以成为代数环，是因为输入和输出信号是同时发生的，这才形成了
代数环。如果能让这两个信号不同时发生，则有望避免极限环。一种显然的方法是
在输出信号后面接一个延迟环节，将其作为输入信号馈入环中的另一个模块。如果
这个延迟的时间常数足够小，则有望用这种方法避开代数环；另一个方法，是在输
出模块后接一个低通滤波器，1/(Ts + 1)，避免直馈现象。如果T 足够小，则用这样
的方法有可能写出代数环。
　　例5-26 试引入小延迟e?τs 消除例5-24 系统中的代数环，并评价处理的效果。
　　解在图5-36（b）的传递函数模型后面加一个常数延迟模块，则可以构造出如图
5-37 所示的Simulink 仿真模型。
图5-37 人为引入延迟环节的仿真模型（文件名：c5malg2a.slx）
假设延迟时间常数为0.01，将仿真算法设置为ode15s，并将相对误差限设置为
10?8，则得出的仿真结果如图5-38 所示，且仿真过程极其耗时，结果存在巨大的误差。
进一步减小延迟时间常数，则可能得出更杂乱无章的结果。显然，这样的方法虽然避开
了代数环，但引入了过多的问题，不适合采用。
图5-38 系统的数值解（结果是错误的）
　　例5-27 试引入低通滤波器消除例5-24 系统的代数环，并评价处理的效果。
第5 章控制系统的Simulink 建模·163·
　　解由给出的等效传递函数可以反推回系统的微分方程
2y′(t) + 4y(t) = u′(t) + 3u(t), y(0) = 0.5
其中，u(t) 为阶跃函数或Heaviside 函数。这样，可以使用下面的方法直接求取微分方程
的解析解。
>> syms t y(t); u=heaviside(t);
y0=dsolve(2*diff(y)+4*y==diff(u)+3*u, y(0)==0.5)
latex(simplify(y0))
可以得出输出信号的解析解为
y0(t) =
1
8
e?2t +
3
8
sign (t) ?
3
8
e?2tsign (t) +
3
8
在图5-36（b）的传递函数模型后面加一个传递函数模块，并将分子多项式设置为1，
分母多项式设置为[T,1]，则可以构造出如图5-39 所示的Simulink 仿真模型。
图5-39 人为引入低通滤波器的仿真模型（文件名：c5malg2b.slx）
选择T = 0.001，对其进行仿真，并计算理论值，则可以得出仿真结果与理论值曲
线，如图5-40 所示。从显示的仿真结果看，二者似乎没有什么区别。
>> T=0.001; tic, [t,~,y1]=sim('c5malg2b'); toc
y0=exp(-2*t)/8+3*sign(t)/8-3*exp(-2*t).*sign(t)/8+3/8;
plot(t,[y1 y0])
图5-40 时域响应曲线的比较
事实上，在初始时刻，二者是有显著区别的。由于原始模型直馈现象的存在，输入信
号在t = 0 时刻就会在输出信号中体现出来，这时输出信号的初值为0.5。如果采用低通
滤波器，则这样的直馈现象被回避了，输出信号在短时间内由0 平缓地变化到0.5，后续
的响应与理论值很接近。这个过程是很短的，如果T 足够小，这个过程是可以忽略不计
·164· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
的。除去这个短暂的时间段，可以得出最大误差为4.6089×10?5，计算点数为6690。
>> ii=t>1e-5; max(abs(y(ii)-y1(ii))), length(t)
减小T 值，则可以得出不同T 值的仿真结果对比（见表5-1）。可以看出，不加低通滤
波器（即表中的“直馈”）虽然导致代数环现象，耗时稍长，其仿真结果是最好的。一般
情况下不必理会代数环，让Simulink 去自行求解代数方程即可。在某些特定的场合下，
Simulink 求解出现困难时，可以考虑引入低通滤波器消除代数环。可以看出，如果T 足
够小，逼近的效果还是比较好的。
表5-1 不同T 值的仿真结果对比
滤波常数T 直馈0.001 0.0001 10?6 10?10
最大误差1.6875×10?14 4.6089×10?5 4.5995×10?6 4.5985×10?8 4.4973×10?12
计算点数823 6690 60675 1353 1128
耗时/s 0.4133 0.1722 0.5211 0.1640 0.1961
对复杂的代数环系统而言，仍然可以使用前面介绍的低通滤波器近似方法，很
好地避开代数环的问题。除了这里介绍的添加环节的近似方法之外，还可以尝试
Simulink 求解代数环的类似于直馈的求解方法。默认情况下，Simulink 采用基于置
信域（trust region）的算法。如果该算法不能很好地解决代数环问题，读者可以尝试
直线搜索（line search）的算法，设置的命令为
set_param(模型名,'AlgebraicLoopSolver','LineSearch')
本章习题
5.1 例5-22 给出了一个控制系统的Simulink 模型c5mlinr1.slx。试由该模型的组成
部分判定模块的执行顺序，并用Simulink 自动排序的方法验证排序结果。
5.2 含有磁滞回环非线性环节的控制系统框图如图5-41 所示，试建立Simulink 仿真
模型，并绘制不同输入幅值下系统的阶跃响应曲线。
1 + 0.8s
s
10
1 + 0.1s
1
s
- - - - - -
?6
u(t) 6- y(t)
图5-41 非线性系统方框图
5.3 试绘制出如图5-42（a）、图5-42（b）所示的非线性反馈系统[23] 的Simulink 仿真框
图，得出系统的阶跃响应曲线，并探讨非线性环节参数对仿真结果的影响。
5.4 建立起如图5-43 所示的非线性系统[53] 的Simulink 框图，并观察在单位阶跃信号
输入下系统的输出曲线和误差曲线。另外，本系统中涉及两个非线性环节的串联，
试问这两个非线性环节可以互换吗？试从仿真结果上加以解释。
第5 章控制系统的Simulink 建模·165·
1
s2
s
6-
1
1
x3
2
4 ?x2
- - - - -
6
 
? 6?
x2
r(t) y(t)
（a）
y(t)
s
6-
1
1
K
s(s + 1)
- -
-
?
6
?1
- -
?
?
r(t)
 
5
s + 5
6
-
（b）
图5-42 习题5.3 图
10
s(s + 1)2
- - - - -
6

-
6
-
6
6 -
0
1
2
1
2
1
1
r(t) y(t)
图5-43 习题5.4 图
5.5 试建立Simulink 模型描述图5-44 中给出的非线性系统，找出其平衡点并求出线
性化模型。试比较原始模型与线性化模型的阶跃响应曲线。
3s2
s3+6.5s2+6.8s+2.5
6 -
1
1 1 + s
s2
- -
6
y(t) - -
?
r(t)-
图5-44 习题5.5 的非线性系统方框图
5.6 建立起如图5-45 所示的非线性系统[54] 的Simulink 框图，并观察在单位阶跃信号
输入下系统的输出曲线和误差曲线。
5.7 已知某系统的Simulink 仿真框图如图5-46 所示，试由该框图写出系统的数学模
型公式。
·166· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
6 -
?0.5
?2.5
0.5
2.5
e?0.4s
30
0.8s + 1
1
s
30
s2 + 6.5s + 1
r(t)- - - - - y-(t)
6

图5-45 习题5.6 的系统方框图
图5-46 习题5.7 的Simulink 仿真框图
5.8 假设线性系统由下面的常微分方程给出：
8>
><>>:
x′1(t) = ?x1(t) + x2(t)
x′2(t) = ?x2(t) ? 3x3(t) + u1(t)
x′3(t) = ?x1(t) ? 5x2(t) ? 3x3(t) + u2(t)
且y(t) = ?x2(t) + u1(t) ? 5u2(t)，式中有两个输入信号u1(t) 与u2(t)，请用两种
方法在Simulink 下将该模型表示出来。
5.9 试将下面给出的多变量系统传递函数矩阵用Simulink 表示出来：
G(s) = 2664
0.806s + 0.264
s2 + 1.15s + 0.202
?15s ? 1.42
s3 + 12.8s2 + 13.6s + 2.36
1.95s2 + 2.12s + 0.49
s3 + 9.15s2 + 9.39s + 1.62
7.15s2 + 25.8s + 9.35
s4 + 20.8s3 + 116.4s2 + 111.6s + 18.8
3775
5.10 对线性传递函数模型G(s) =
6(s + 1)
(s + 2)(s + 3)(s + 4)
在白噪声信号γ(t) 激励下的
行为进行仿真，求出其输出函数的概率密度曲线，并与由仿真结果求出的结果进
行比较；另外，对此系统进行自相关函数和概率谱密度分析，并比较仿真结果和理
论结果是否一致。
5.11 考虑图5-47 中给出的非线性系统框图，如果随机输入信号γ(t) 是Gauss 白噪
声信号，其均值为0、方差为4，试通过仿真方法求出误差信号e(t) 的概率密度
分布。
第5 章控制系统的Simulink 建模·167·
s + 1
s3+1.9s2+1.1s+0.18
-
6
- - - ?- -
6
r(t)=0 y(t)
Δ = 2
γ(t)
e(t)
?
图5-47 习题5.11 非线性系统框图
5.12 假设系统的受控对象与控制器模型分别为
G(s) =
e?s
s(s + 1)4 , Gc(s) = Kp = 0.173 + 0.348s
试得出闭环系统模型，并得出其Padé 近似模型。选择一个合适的Padé 近似阶次，
使得闭环模型的阶跃响应可以由Padé 近似模型很好地逼近。
5.13 已知一级倒立摆的数学模型为：
8>
><>>:
y′′ =
f/m + lθ2 sin θ ? g sin θ cos θ
M/m + sin2 θ
θ′′ = ?f cos θ/m + (M + m)g sin θ/m ? lθ2 sin θ cos θ
l(M/m + sin2 θ)
其中，θ 为摆体与垂直方向的夹角（单位为rad），y 为小车的位移（单位为m），f
为电机对小车的作用力（单位为N），M 和m分别为小车和摆体的质量（单位为
kg），l 为摆长的一半（单位为m），g 为重力加速度（9.81 m/s2）。试建立起倒立摆的
Simulink 模型。若取m = 0.21 kg，M = 0.455 kg，l = 0.61/2 m，并取f 为系统的输
入信号，试在平衡点y = θ = 0 处对该系统进行线性化，并比较原系统和线性化系
统的阶跃响应曲线。
5.14 文献[55] 将F-14 战斗机模型作为基准测试问题，用来测试控制系统仿真软件的
功能。MATLAB 提供了该模型的Simulink 实现，如图5-48 所示。试求出输出端1
对输入端1 的等效线性模型，并得出该模型全部的零极点。
5.15 试将图5-49 中给出的非线性系统模型输入Simulink 环境中，并得出该系统在阶
跃输入下的工作点及线性化模型。
5.16 某模型参考自适应系统的框图如图5-50 所示[56]。其中，模型与控制器参数为
b0 = 0.5，a1 = 0.447，a2 = 0.1，a3 = 2，d0 = 1，d1 = 0.5，k1 = 0.03，k2 = 1，
?b
0(0) = 0.2。该框图使用了两个乘法器及相关运算来改变自适应控制的增益
?b
0(t)，使得系统的输出可以跟踪参考模型的输出。试根据给出的系统框图建立系
统的Simulink 模型，并得出方波信号激励下系统的响应曲线，并绘制自适应增益
?b
0(t) 的曲线。
·168· 薛定宇教授大讲堂（卷VI）：Simulink 建模与仿真
图5-48 F-14 战斗机的Simulink 模型（文件名：f14.slx）
3 1
s + 1 6
s
s2 + 2
4s + 2
(s + 1)2
1
s2
50
s2 + 2
s3 + 14
- - - - - - -
 

6

6
r(t) y(t)
?
图5-49 习题5.15 图
1
a2s + a1
1
s
- - -
6
1
a2s + a3
1
s
- - -
6
?
6
-
b0
d1
- d0
?
6 - ?
k2s + k1
s
  
6
-
-
6
? ?b
0(0)
ym(t)
ys(t)
u(t)
?(t)
?′(t)
?
?
?
?
?b
0(t)
图5-50 模型参考自适应控制系统