matlab在自控原理中的应用

控制系统的模拟试验与MATLAB仿真

1 MATLAB简介

（略）

2 MATLAB基本操作命令

本节简单介绍与本书内容相关的一些基本知识和操作命令。

（1）简单矩阵的输入

MATLAB是一种专门为矩阵运算设计的语言，所以在MATLAB中处理的所有变量都是矩阵。这就是说，MATLAB只有一种数据形式，那就是矩阵，或者数的矩形阵列。标量可看作为1×1的矩阵，向量可看作为n×1或1×n的矩阵。这就是说，MATLAB语言对矩阵的维数及类型没有限制，即用户无需定义变量的类型和维数，MATLAB会自动获取所需的存储空间。

输入矩阵最便捷的方式为直接输入矩阵的元素，其定义如下： （1） 元素之间用空格或逗号间隔；

（2） 用中括号（[]）把所有元素括起来； （3） 用分号（；）指定行结束。

例如，在MATLAB的工作空间中，输入：

5　　6　　9] >> a?[2 3 4；　则输出结果为：

a?234

569矩阵a被一直保存在工作空间中，以供后面使用，直至修改它。

MATLAB的矩阵输入方式很灵活，大矩阵可以分成n行输入，用回车符代替分号或用续行符号（?）将元素续写到下一行。例如：

a?[1，　2，　3；　4，　5，　6；　7，　8，　9]

a?[ 1 2 3 4 5 6 7 8 9]a?[1，　2，　3；　4，　5，　???

　6；　7，　8，　9]以上三种输入方式结果是相同的。一般若长语句超出一行，则换行前使用续行符号（?）。

在MATLAB中，矩阵元素不限于常量，可以采用任意形式的表达式。同时，除了直接输入方式之外，还可以采用其它方式输入矩阵，如：

（1） 利用内部语句或函数产生矩阵； （2） 利用M文件产生矩阵；

- 1 -

（3） 利用外部数据文件装入到指定矩阵。

（2）复数矩阵输入

MATLAB允许在计算或函数中使用复数。输入复数矩阵有两种方法： （1） a=[12;34]+i*[56;78]

（2） a=[1+5i 2+6i;3+7i 4+8i]

注意，当矩阵的元素为复数时，在复数实部与虚部之间不允许使用空格符。如1 ＋5i将被认为是1和5i两个数。另外，MATLAB表示复数时，复数单位也可以用j。

（3） MATLAB语句和变量

MATLAB是一种描述性语言。它对输入的表达式边解释边执行，就象BASIC语言中直接执行语句一样。

MATLAB语句的常用格式为：

变量＝表达式[;]

或简化为：

表达式[;]

表达式可以由操作符、特殊符号、函数、变量名等组成。表达式的结果为一矩阵，它赋给左边的变量，同时显示在屏幕上。如果省略变量名和“＝”号，则MATLAB自动产生一个名为ans的变量来表示结果，如：

1900∕81 结果为：

ans? 23.4568

ans 是MATLAB提供的固定变量，具有特定的功能，是不能由用户清除的。常用的固定

变量还有eps、pi、Inf、NaN等。其特殊含义可以用7.2.10节介绍的方法查阅帮助。

MATAB允许在函数调用时同时返回多个变量,而一个函数又可以由多种格式进行调用，语句的典型格式可表示为：

[返回变量列表]＝fun-name(输入变量列表)

例如用bode()函数来求取或绘制系统的Bode图，可由下面的格式调用：

[mag,phase]?bode(num,den,W)

其中变量num、den表示系统传递函数分子和分母，W表示指定频段，mag为计算幅值，phase为计算相角。

（4）语句以“％”开始和以分号“；”结束的特殊效用

在MATLAB中以“％”开始的程序行，表示注解和说明。符号“％”类似于C++中的“//”。这些注解和说明是不执行的。这就是说，在MATLAB程序行中，出现“％”以后的一切内容都是可以忽略的。

分号用来取消打印，如果语句最后一个符号是分号，则打印被取消，但是命令仍在执行，而结果不再在命令窗口或其它窗口中显示。这一点在M文件中大量采用，以抑制不必要的信息显示。

（5）获取工作空间信息

MATLAB开辟有一个工作空间，用于存储已经产生的变量。变量一旦被定义，MATLAB系统会自动将其保存在工作空间里。在退出程序之前，这些变量将被保留在存储器中。

为了得到工作空间中的变量清单，可以在命令提示符>>后输入who 或 whos 命令，当

- 2 -

前存放在工作空间的所有变量便会显示在屏幕上。

命令clear能从工作空间中清除所有非永久性变量。如果只需要从工作空间中清除某个特定变量，比如“x”,则应输入命令clear x。

（6）常数与算术运算符

MATLAB采用人们习惯使用的十进制数。如： 3 –99 0.0001 9.6397238 1.60210e?20 6.62252e23 2i -3.14159i 3e5i 其中 i??1。

数值的相对精度为eps，它是一个符合IEEE标准的16位长的十进制数，其范围为：

10?308~10308。

MATLAB提供了常用的算术运算符：+，－，?，∕（﹨），＾（幂指数）。 应该注意：（∕）右除法和（﹨）左除法这两种符号对数值操作时，其结果相同，其斜线下为分母，如1∕4与4﹨1,其结果均为0.25，但对矩阵操作时，左、右除法是有区别的。

（7）选择输出格式

输出格式是指数据显示的格式，MATLAB提供format命令可以控制结果矩阵的显示，而不影响结果矩阵的计算和存储。所有计算都是以双精度方式完成的。

（1） 如果矩阵的所有元素都是整数，则矩阵以不带小数点的格式显示。 如输入：

0 1] x?[?1 则显示：

x? -1 0 1

（2） 如果矩阵中至少有一个元素不是整数，则有多种输出格式。常见格式有以下四

种：

① format short（短格式,也是系统默认格式） ② format short e（短格式科学表示） ③ format long（长格式）

④ format long e（长格式科学表示） 如：

1.2345e?6] x?[4/3 对于以上四种格式，其显示结果分别为：

x?

1.3333 0.0000x? 短格式5位表示

1.3333e?00 1.2345e?06

短格式科学表示

x? 1.33333333333333 0.00000123450000- 3 -

长格式16位表示

x? 1.33333333333333 e?00 1.234500000000000e?06（8）MATLAB图形窗口

长格式科学表示

一旦调用了某种格式，则这种被选用的格式将保持，直到对格式进行了改变为止。

当调用了一个产生图形的函数时，MATLAB会自动建立一个图形窗口。这个窗口还可分裂成多个窗口，并可在它们之间选择，这样在一个屏上可显示多个图形。

图形窗口中的图形可通过打印机打印出来。若想将图形导出并保存，可用鼠标点击菜单File|Export,导出格式可选emp、bmp、jpg等。命令窗口的内容也可由打印机打印出来：如果事先选择了一些内容，则可打印出所选择的内容；如果没有选择内容，则可打印出整个工作空间的内容。

（9）剪切板的使用

利用Windows的剪切板可在MATLAB与其它应用程序之间交换信息。

（1） 要将MATLAB的图形移到其它应用程序，首先按Alt-Print Screen键，将图形复制到剪切板中，然后激活其它应用程序，选择edit（编辑）中的paste（粘贴），就可以在应用程序中得到MATLAB中的图形。当然还可以借助于copy to Bitmap 或copy to Metafile 选项来传递图形信息。

（2） 要将其它应用程序中的数据传递到MATLAB，应先将数据放入剪切板，然后在MATLAB中定义一个变量来接收。

如键入： q=[

然后选择Edit中的paste，最后加上“]”,这样可将应用程序中的数据送入MATLAB的q变量中。

（10）MATLAB编程指南

MATLAB的编程效率比BASIC、C、FORTRAN和PASCAL等语言要高，且易于维护。在编写小规模的程序时，可直接在命令提示符>>后面逐行输入，逐行执行。对于较复杂且经常重复使用的程序，可按7.1.3介绍的方法进入程序编辑器编写M文件。

M文件是用MATLAB语言编写的可在MATLAB环境中运行的磁盘文件。它为脚本文件（Script File）和函数文件（Function File）,这两种文件的扩展名都是.m。

（1） 脚本文件是将一组相关命令编辑在一个文件中，也称命令文件。脚本文件的语句可以访问MATLAB工作空间中的所有数据，运行过程中产生的所有变量都是全局变量。例如下述语句如果以.m为扩展名存盘，就构成了M脚本文件,我们不妨将其文件名取为“Step_Response”。

% 用于求取一阶跃响应。

num=[1 4]; den=[1 2 8]; step(num,den)

当你键入help Step_Response时，屏幕上将显示文件开头部分的注释： 用于求取一阶跃响应。

很显然，在每一个M文件的开头，建立详细的注释是非常有用的。由于MATLAB提供了大量的命令和函数，想记住所有函数及调用方法一般不太可能，通过联机帮助命令help可容易地对想查询的各个函数的有关信息进行查询。该命令使用格式为：

help 命令或函数名

注意：若用户把文件存放在自己的工作目录上，在运行之前应该使该目录处在MATLAB

- 4 -

的搜索路径上。当调用时，只需输入文件名，MATLAB就会自动按顺序执行文件中的命令。

（2） 函数文件是用于定义专用函数的，文件的第一行是以function作为关键字引导的，后面为注释和函数体语句。

函数就像一个黑箱，把一些数据送进去，经加工处理，再把结果送出来。在函数体内使用的除返回变量和输入变量这些在第一行functon语句中直接引用的变量外，其它所有变量都是局部变量，执行完后，这些内部变量就被清除了。

函数文件的文件名与函数名相同（文件名后缀为.m）,它的执行与命令文件不同，不能键入其文件名来运行函数，M函数必须由其它语句来调用，这类似于C语言的可被其它函数调用的子程序。M函数文件一旦建立，就可以同MATLAB基本函数库一样加以使用。

例1 求一系列数的平均数，该函数的文件名为“mean.m”

function y=mean(x)

% 这是一个用于求平均数的函数

w=length(x); % length函数表示取向量x的长度 y=sum(x)/w; % sun函数表示求各元素的和

该文件第一行为定义行，指明是mean函数文件，y 是输出变量，x是输入变量，其后的％开头的文字段是说明部分。真正执行的函数体部分仅为最后二行。其中变量w是局部变量，程序执行完后，便不存在了。

在MATLAB命令窗口中键入

>> r=1:10; % 表示r变量取1到10共10个数 mean(r)

运行结果显示 ans =

5.5000

该例就是直接使用了所建立的M函数文件，求取数列r的平均数。

3 MATLAB在控制系统中的应用

MATLAB是国际控制界目前使用最广的工具软件，几乎所有的控制理论与应用分支中都有MATLAB工具箱。本节结合前面所学自控理论的基本内容，采用控制系统工具箱（Control Systems Toolbox）和仿真环境（Simulink），学习MATLAB的应用。

（1） 用MATLAB建立传递函数模型

1．有理函数模型

线性系统的传递函数模型可一般地表示为：

b1sm?b2sm?1?????bms?bm?1 G(s)? n?m (1)

sn?a1sn?1?????an?1s?an将系统的分子和分母多项式的系数按降幂的方式以向量的形式输入给两个变量num和den，就可以轻易地将传递函数模型输入到MATLAB环境中。命令格式为：

num?[b1,b2,???,bm,bm?1]; （2） den?[1,a1,a2,???,an?1,an]; （3）

在MATLAB控制系统工具箱中，定义了tf() 函数，它可由传递函数分子分母给出的变量构造出单个的传递函数对象。从而使得系统模型的输入和处理更加方便。

该函数的调用格式为：

- 5 -

G＝tf(num，den); （4）

例2 一个简单的传递函数模型:

G(s)?s?5 432s?2s?3s?4s?5可以由下面的命令输入到MATLAB工作空间中去。 >> num=[1，5];

den=[1，2，3，4，5]; G=tf(num，den) 运行结果：

Transfer function:

s + 5

----------------------------- s^4 + 2s^3 + 3s^2 + 4s + 5

这时对象G可以用来描述给定的传递函数模型,作为其它函数调用的变量。 例3 一个稍微复杂一些的传递函数模型：

G(s)?6(s?5)

(s2?3s?1)2(s?6)该传递函数模型可以通过下面的语句输入到MATLAB工作空间。 >> num=6*[1，5];

den=conv(conv([1，3，1]，[1，3，1])，[1，6]); tf(num,den) 运行结果

Transfer function: 6 s + 30

----------------------------------------- s^5 + 12 s^4 + 47 s^3 + 72 s^2 + 37 s + 6

其中conv()函数（标准的MATLAB函数）用来计算两个向量的卷积，多项式乘法当然也可以用这个函数来计算。该函数允许任意地多层嵌套，从而表示复杂的计算。 2.零极点模型

线性系统的传递函数还可以写成极点的形式：

G(s)?K(s?z1)(s?z2)???(s?zm) （5）

(s?p1)(s?p2)???(s?pn)将系统增益、零点和极点以向量的形式输入给三个变量KGain、Z和P，就可以将系统的零极点模型输入到MATLAB工作空间中，命令格式为：

KGain?K; (6)

(7) Z?[?z1；?z2；???；?zm]； (8) P?[?p1；?p2；???；?pn]；- 6 -

在MATLAB控制工具箱中，定义了zpk()函数，由它可通过以上三个MATLAB变量构造出零极点对象，用于简单地表述零极点模型。该函数的调用格式为：

G=zpk(Z,P,KGain) (9)

例4 某系统的零极点模型为： G(s)?6(s?1.9294)(s?0.0353?0.9287j)

(s?0.9567?1.2272j)(s?0.0433?0.6412j)该模型可以由下面的语句输入到MATLAB工作空间中。

>> KGain=6;

z=[-1.9294；-0.0353+0.9287j；-0.0353-0.9287j];

p=[-0.9567+1.2272j；-0.9567-1.2272j；0.0433+0.6412j；0.0433-0.6412j]; G=zpk(Z，P，KGain)

运行结果:

Zero/pole/gain:

6 (s+1.929) (s^2 + 0.0706s + 0.8637)

---------------------------------------------- (s^2 - 0.0866s + 0.413) (s^2 + 1.913s + 2.421)

注意：对于单变量系统，其零极点均是用列向量来表示的，故Z、P向量中各项均用分号（;）隔开。

3. 反馈系统结构图模型

设反馈系统结构图如图5所示。

控制系统工具箱中提供

了feedback()函数，用来求

图5 反馈系统结构图

取反馈连接下总的系统模

型，该函数调用格式如下：

G=feedback(G1,G2,sign); （10）

其中变量sign用来表示正反馈或负反馈结构，若sign=-1表示负反馈系统的模型，若省略sign变量，则仍将表示负反馈结构。G1和G2分别表示前向模型和反馈模型的LTI(线性时不变)对象。

例5 若反馈系统图5中的两个传递函数分别为： G1(s)?11G(s)?， 2s?1(s?1)2则反馈系统的传递函数可由下列的MATLAB命令得出

>> G1=tf(1,[1,2,1]); G2=tf(1,[1,1]); G=feedback(G1,G2) 运行结果:

- 7 -

Transfer function:

s + 1

--------------------- s^3 + 3 s^2 + 3 s + 2

若采用正反馈连接结构输入命令

>> G=feedback(G1,G2,1) 则得出如下结果:

Transfer function: s + 1

----------------- s^3 + 3 s^2 + 3 s

例6 若反馈系统为更复杂的结构如图6所示。其中

10s?51s3?7s2?24s?24G2(s)?H(s)?G1(s)?4，， 32s0.01s?1s?10s?35s?50s?24则闭环系统的传递函数可以由下面的MATLAB命令得出：

>> G1=tf([1,7,24,24],[1,10,35,50,24]);

G2=tf([10,5],[1,0]); H=tf([1],[0.01,1]); G_a=feedback(G1*G2,H)

得到结果：

Transfer function:

0.1 s^5 + 10.75 s^4 + 77.75 s^3 + 278.6 s^2 + 361.2 s + 120

-------------------------------------------------------------------- 0.01 s^6 + 1.1 s^5 + 20.35 s^4 + 110.5 s^3 + 325.2 s^2 + 384 s + 120

4. 有理分式模型与零 极点模型的转换

图6 复杂反馈系统

有了传递函数的

有理分式模型之后，求取零极点模型就不是一件困难的事情了。在控制系统工具箱中，可以由zpk()函数立即将给定的LTI对象G转换成等效的零极点对象G1。该函数的调用格式为： G1=zpk(G) (11)

例7 给定系统传递函数为：

6.8s2?61.2s?95.2 G(s)?4 32s?7.5s?22s?19.5s对应的零极点格式可由下面的命令得出 >> num=[6.8, 61.2, 95.2];

- 8 -

den=[1, 7.5, 22, 19.5, 0]; G=tf(num,den); G1=zpk(G) 显示结果：

Zero/pole/gain:

6.8 (s+7) (s+2)

------------------------- s (s+1.5) (s^2 + 6s + 13) 可见，在系统的零极点模型中若出现复数值，则在显示时将以二阶因子的形式表示相应的共轭复数对。

同样，对于给定的零极点模型，也可以直接由MATLAB语句立即得出等效传递函数模型。调用格式为：

G1=tf(G) (12) 例8 给定零极点模型：

G(s)?6.8(s?2)(s?7))

s(s?3?j2)(s?1.5)可以用下面的MATLAB命令立即得出其等效的传递函数模型。输入程序的过程中要注意大小写。

>> Z=[-2,-7];

P=[0,-3-2j,-3+2j,-1.5]; K=6.8;

G=zpk(Z,P,K); G1=tf(G)

结果显示：

Transfer function:

6.8 s^2 + 61.2 s + 95.2

------------------------------- s^4 + 7.5 s^3 + 22 s^2 + 19.5 s

5. Simulink建模方法

在一些实际应用中，如果系统的结构过于复杂，不适合用前面介绍的方法建模。在这种情况下，功能完善的Simulink程序可以用来建立新的数学模型。Simulink是由Math Works 软件公司1990年为MATLAB提供的新的控制系统模型图形输入仿真工具。它具有两个显著的功能：Simul(仿真)与Link(连接)，亦即可以利用鼠标在模型窗口上“画”出所需的控制系统模型。然后利用SIMULINK提供的功能来对系统进行仿真或线性化分析。与MATLAB中逐行输入命令相比，这样输入更容易，分析更直观。下面简单介绍SIMULINK建立系统模型的基本步骤：

（1） SIMULINK的启动：在MATLAB命令窗口的工具栏中单击按钮

或者在命令提示符>>

下键入simulink命令，回车后即可启动Simulink程序。启动后软件自动打开Simullink模型库窗口，如图 7所示。这一模型库中含有许多子模型库，如Sources(输入源模块库)、Sinks(输出显示模块库)、Nonlinear(非线性环节)等。若想建立一个控制系统结构框图，则应该选择File| New菜单中的Model选项，或选择工具栏上new Model

按钮，打开一个

- 9 -

空白的模型编辑窗口如图 8所示。

(2) 画出系统的各个模块：打开相应的子模块库，选择所需要的元素，用鼠标左键点中后拖

图 7 simulink 模型库

图8 模型编辑窗口

到模型编辑窗口的合适位置。 (3) 给出各个模块参数：由于选中的各个模块只包含默认的模型参数，如默认的传递函数模

型为1/(s+1)的简单格式，必须通过修改得到实际的模块参数。要修改模块的参数，可以用鼠标双击该模块图标，则会出现一个相应对话框，提示用户修改模块参数。

(4) 画出连接线：当所有的模块都画出来之后，可以再画出模块间所需要的连线，构成完整

的系统。模块间连线的画法很简单，只需要用鼠标点按起始模块的输出端（三角符号），再拖动鼠标，到终止模块的输入端释放鼠标键，系统会自动地在两个模块间画出带箭头的连线。若需要从连线中引出节点，可在鼠标点击起始节点时按住Ctrl键，再将鼠标拖动到目的模块。

(5) 指定输入和输出端子：在Simulink下允许有两类输入输出信号，第一类是仿真信号，

可从source(输入源模块库)图标中取出相应的输入信号端子，从Sink(输出显示模块库)图标中取出相应输出端子即可。第二类是要提取系统线性模型，则需打开Connection(连接模块库)图标，从中选取相应的输入输出端子。

例9 典型二阶系统的结构图如图9所示。用SIMULINK对系统进行仿真分析。

图9 典型二阶系统结构图

按前面步骤，启动simulink并打开一个空白的模型编辑窗口。

（1） 画出所需模块，并给出正确的参数：

? 在sources子模块库中选中阶跃输入（step）图标，将其拖入编辑窗口，并用鼠标

左键双击该图标，打开参数设定的对话框，将参数step time(阶跃时刻)设为0。 ? 在Math(数学)子模块库中选中加法器（sum）图标，拖到编辑窗口中，并双击该图

标将参数List of signs(符号列表)设为|+-（表示输入为正，反馈为负）。

? 在continuous(连续)子模块库中、选积分器（Integrator）和传递函数（Transfer

Fcn）图标拖到编辑窗口中，并将传递函数分子（Numerator）改为〔900〕，分母

- 10 -

（Denominator）改为〔1，9〕。

? 在sinks(输出)子模块库中选择scope(示波器)和Out1(输出端口模块)图标并将之

拖到编辑窗口中。

（3）将画出的所有模块按图9用鼠标连接起来，构成一个原系统的框图描述如图10所示。 （4）选择仿真算法和仿真控制参数，启动仿真过程。

● 在编辑窗口中点击Simulation|Simulation parameters菜单，会出现一个参数对话框，

在solver模板中设置响应的仿真范围StartTime(开始时间)和StopTime(终止时间)，仿真步长范围Maxinum step size(最大步长)和Mininum step size(最小步长)。对于本例，StopTime可设置为2。最后点击Simulation|Start菜单或点击相应的热键启动仿真。双击示波器，在弹出的图形上会“实时地”显示出仿真结果。输出结果如图11所示。

在命令窗口中键入whos命令，会发现工作空间中增加了两个变量――tout和yout，这是因为Simulink中的Out1 模块自动将结果写到了MATLAB的工作空间中。利用MATLAB命令plot(tout,yout),可将结果绘制出来，如图12所示。比较11和12，可以发现这两种输

图10 二阶系统的simulink实现

出结果是完全一致的。

图11仿真结果示波器显示

图12 MATLAB命令得出的系统响应曲线

(2) 利用MATLAB进行时域分析

- 11 -

1. 线性系统稳定性分析

线性系统稳定的充要条件是系统的特征根均位于S平面的左半部分。系统的零极点模型可以直接被用来判断系统的稳定性。另外，MATLAB语言中提供了有关多项式的操作函数，也可以用于系统的分析和计算。 （1）直接求特征多项式的根

设p为特征多项式的系数向量，则MATLAB函数roots()可以直接求出方程p=0在复数范围内的解v,该函数的调用格式为：

v=roots(p) (13)

例10 已知系统的特征多项式为： x?3x?2x?x?1

特征方程的解可由下面的MATLAB命令得出。 >> p=[1,0,3,2,1,1]; v=roots(p) 结果显示： v =

0.3202 + 1.7042i 0.3202 - 1.7042i -0.7209 0.0402 + 0.6780i 0.0402 - 0.6780i

利用多项式求根函数roots(),可以很方便的求出系统的零点和极点，然后根据零极点分析系统稳定性和其它性能。 （2）由根创建多项式

如果已知多项式的因式分解式或特征根，可由MATLAB函数poly()直接得出特征多项式系数向量，其调用格式为：

p=poly(v) (14)

如上例中：

v=[0.3202+1.7042i;0.3202-1.7042i;

-0.7209;0.0402+0.6780i; 0.0402-0.6780i];

>> p=poly(v) 结果显示

p =

1.0000 -0.0000 3.0000 2.0000 1.0000 1.0000 由此可见，函数roots()与函数poly()是互为逆运算的。 （3）多项式求值

在MATLAB 中通过函数polyval()可以求得多项式在给定点的值，该函数的调用格式为： polyval(p,v) (15) 对于上例中的p值，求取多项式在x点的值，可输入如下命令：

>> p=[1,0,3,2,1,1]; x=1

polyval(p,x) 结果显示

532- 12 -

ans =

8

(4)部分分式展开 考虑下列传递函数：

M(s)numb0sn?b1sn?1?????bn ??N(s)dena0sn?a1sn?1?????an式中a0?0，但是ai和bj中某些量可能为零。 MATLAB函数可将函数的调用格式为：

?r,p,k??residue(num,den) （16）

M(s)展开成部分分式，直接求出展开式中的留数、极点和余项。该N(s)则

M(s)的部分分式展开由下式给出： N(s)

M(s)r(1)r(2)r(n)????????k(s) N(s)s?p(1)s?p(2)s?p(n)式中p(1)??p1,p(2)??p2,?, p(n)??pn,为极点，

r(1)??r1,r(2)??r2,?, r(n)??rn为各极点的留数，k(s)为余项。

例11 设传递函数为：

2s3?5s2?3s?6 G(s)?3 2s?6s?11s?6该传递函数的部分分式展开由以下命令获得： >> num=[2,5,3,6]; den=[1,6,11,6];

[r,p,k]=residue(num,den) 命令窗口中显示如下结果

r= p= k= -6.0000 -3.0000 2 -4.0000 -2.0000 3.0000 -1.0000

中留数为列向量r，极点为列向量p，余项为行向量k。 由此可得出部分分式展开式：

G(s)??6?43???2 s?3s?2s?1该函数也可以逆向调用，把部分分式展开转变回多项式

M(s)之比的形式，命令格式为： N(s)- 13 -

[num,den]=residue(r,p,k) (17) 对上例有：

>> [num,den]=residue(r,p,k) 结果显示 num=

2.0000 5.0000 3.0000 6.0000

den=

1.0000 6.0000 11.0000 6.0000

应当指出，如果p(j)=p(j+1)=?=p(j+m-1),则极点p(j)是一个m重极点。在这种情况下，部分分式展开式将包括下列诸项：

r(j)r(j?1)r(j?m?1) ??????2ms?p(j)?s?p(j)??s?p(j)?例12 设传递函数为：

s2?2s?3s2?2s?3 G(s)? ?332(s?1)s?3s?3s?1 则部分分式展开由以下命令获得：

>> v=[-1,-1,-1] num=[0,1,2,3];

den=poly(v);

[r,p,k]=residue(num,den)

结果显示

r=

1.0000 0.0000 2.0000 p=

-1.0000 -1.0000 -1.0000 k=

[ ]

其中由poly()命令将分母化为标准降幂排列多项式系数向量den, k=[]为空矩阵。 由上可得展开式为： G(s)?102???0 s?1(s?1)2(s?1)3（5）由传递函数求零点和极点。

在MATLAB控制系统工具箱中，给出了由传递函数对象G求出系统零点和极点的函数，其调用格式分别为：

Z=tzero(G) (18) P=G.P{1} (19) 注意：式19中要求的G必须是零极点模型对象，且出现了矩阵的点运算“.”和大括号

- 14 -

｛｝表示的矩阵元素，详细内容参阅后面章节。

例13 已知传递函数为：

6.8s2?61.2s?95.2 G(s)?4

s?7.5s3?22s2?19.5s输入如下命令：

num=[6.8,61.2,95.2]; den=[1,7.5,22,19.5,0]; G=tf(num,den); G1=zpk(G);

Z=tzero(G) P=G1.P{1}

结果显示 Z =

-7 -2 P =

0 -3.0000 + 2.0000i

-3.0000 - 2.0000i -1.5000

其结果与例8完全一致。 （6）零极点分布图。

在MATLAB中，可利用pzmap()函数绘制连续系统的零、极点图，从而分析系统的稳定性，该函数调用格式为：

pzmap(num,den) (20) 例 14 给定传递函数:

3s4?2s3?5s2?4s?6 G(s)?5

s?3s4?4s3?2s2?7s?2 利用下列命令可自动打开一个图形窗口，显示该系统的零、极点分布图，如图13所示。 >> num=[3,2,5,4,6]; den=[1,3,4,2,7,2];

pzmap(num,den)

title(1Pole-Zero Map1) % 图形标题。

- 15 -

图13 MATLAB函数零、极点分布图

2. 系统动态特性分析。

（1）时域响应解析算法――部分分式展开法。

用拉氏变换法求系统的单位阶跃响应，可直接得出输出c(t)随时间t变化的规律，对于高阶系统，输出的拉氏变换象函数为：

C(s)?G(s)?1num1?? （21） sdens对函数c(s)进行部分分式展开，我们可以用num,[den,0]来表示c(s)的分子和分母。 例 15 给定系统的传递函数：

s3?7s2?24s?24 G(s)?4

s?10s3?35s2?50s?24用以下命令对

G(s)进行部分分式展开。 s>> num=[1,7,24,24] den=[1,10,35,50,24]

[r,p,k]=residue(num,[den,0])

输出结果为

r= p= k=

-1.0000 -4.0000 [ ] 2.0000 -3.0000 -1.0000 -2.0000

-1.0000 -1.0000

1.0000 0

输出函数c(s)为：

- 16 -

C(s)??12111?????0 s?4s?3s?2s?1s拉氏变换得：

c(t)??e?4t?2e?3t?e?2t?e?t?1

（2）单位阶跃响应的求法：

控制系统工具箱中给出了一个函数step()来直接求取线性系统的阶跃响应，如果已知传递函数为：

G(s)?num den则该函数可有以下几种调用格式：

step(num,den) (22) step(num,den,t) (23)

或

step(G) (24)

step(G,t) (25)

该函数将绘制出系统在单位阶跃输入条件下的动态响应图，同时给出稳态值。对于式23和25，t为图像显示的时间长度，是用户指定的时间向量。式22和24的显示时间由系统根据输出曲线的形状自行设定。

如果需要将输出结果返回到MATLAB工作空间中，则采用以下调用格式：

c=step(G) (26) 此时，屏上不会显示响应曲线，必须利用plot()命令去查看响应曲线。plot 可以根据两 个或多个给定的向量绘制二维图形，详细介绍可以查阅后面的章节。

例16 已知传递函数为： G(s)?25 2s?4s?25

利用以下MATLAB命令可得阶跃响应曲线如图14所示。

图14 MATLAB绘制的响应曲线

>> num=[0,0,25]; den=[1,4,25];

- 17 -

step(num,den)

grid % 绘制网格线。

title(1Unit-Step Response of G(s)=25/(s^2+4s+25) 1) % 图像标题 我们还可以用下面的语句来得出阶跃响应曲线 >> G=tf([0,0,25],[1,4,25]);

t=0:0.1:5; % 从0到5每隔0.1取一个值。 c=step(G,t); % 动态响应的幅值赋给变量c

plot(t,c) % 绘二维图形，横坐标取t,纵坐标取c。 Css=dcgain(G) % 求取稳态值。

系统显示的图形类似于上一个例子，在命令窗口中显示了如下结果

Css= 1

（3）求阶跃响应的性能指标

MATLAB提供了强大的绘图计算功能，可以用多种方法求取系统的动态响应指标。我们首先介绍一种最简单的方法――游动鼠标法。对于例16，在程序运行完毕后，用鼠标左键点击时域响应图线任意一点，系统会自动跳出一个小方框，小方框显示了这一点的横坐标（时间）和纵坐标（幅值）。按住鼠标左键在曲线上移动，可以找到曲线幅值最大的一点――即曲线最大峰值，此时小方框中显示的时间就是此二阶系统的峰值时间，根据观察到的稳态值和峰值可以计算出系统的超调量。系统的上升时间和稳态响应时间可以依此类推。这种方法简单易用，但同时应注意它不适用于用plot()命令画出的图形。

另一种比较常用的方法就是用编程方式求取时域响应的各项性能指标。与上一段介绍的游动鼠标法相比，编程方法稍微复杂，但通过下面的学习，读者可以掌握一定的编程技巧，能够将控制原理知识和编程方法相结合，自己编写一些程序，获取一些较为复杂的性能指标。

通过前面的学习，我们已经可以用阶跃响应函数step( )获得系统输出量，若将输出量返回到变量y中，可以调用如下格式

[y,t]=step(G) (27) 该函数还同时返回了自动生成的时间变量t，对返回的这一对变量y和t的值进行计算，可以得到时域性能指标。

① 峰值时间(timetopeak)可由以下命令获得：

[Y,k]=max(y); (28) timetopeak=t(k) (29) 应用取最大值函数max()求出y的峰值及相应的时间，并存于变量Y和k中。然后在变量t中取出峰值时间，并将它赋给变量timetopeak。

② 最大(百分比)超调量(percentovershoot)可由以下命令得到： C=dcgain(G);

[Y,k]=max(y); (30) percentovershoot=100*(Y-C)/C (31) dcgain( )函数用于求取系统的终值，将终值赋给变量C，然后依据超调量的定义，由Y和C计算出百分比超调量。

③ 上升时间(risetime)可利用MATLAB中控制语句编制M文件来获得。首先简单介绍一下循环语句while的使用。

while循环语句的一般格式为： while<循环判断语句> 循环体

- 18 -

end

其中,循环判断语句为某种形式的逻辑判断表达式。

当表达式的逻辑值为真时，就执行循环体内的语句；当表达式的逻辑值为假时，就退出当前的循环体。如果循环判断语句为矩阵时，当且仅当所有的矩阵元素非零时，逻辑表达式的值为真。为避免循环语句陷入死循环，在语句内必须有可以自动修改循环控制变量的命令。

要求出上升时间，可以用while语句编写以下程序得到： C=dcgain(G); n=1;

while y(n)

risetime=t(n)

在阶跃输入条件下，y 的值由零逐渐增大，当以上循环满足y=C时，退出循环，此时对应的时刻，即为上升时间。

对于输出无超调的系统响应，上升时间定义为输出从稳态值的10%上升到90%所需时间，则计算程序如下：

C=dcgain(G); n=1;

while y(n)<0.1*C n=n+1; end m=1;

while y(n)<0.9*C m=m+1; end

risetime=t(m)-t(n)

④ 调节时间(setllingtime)可由while语句编程得到： C=dcgain(G);

i=length(t);

while(y(i)>0.98*C)&(y(i)<1.02*C) i=i-1; end

setllingtime=t(i)

用向量长度函数length( )可求得t序列的长度，将其设定为变量i的上限值。 例 17 已知二阶系统传递函数为:

G(s)?3

(s?1?3i)(s?1?3i)利用下面的stepanalysis.m程序可得到阶跃响应如图 15及性能指标数据。 >> G=zpk([ ],[-1+3*i,-1-3*i ],3);

% 计算最大峰值时间和它对应的超调量。 C=dcgain(G) [y,t]=step(G);

plot(t,y)

- 19 -

grid

[Y,k]=max(y); timetopeak=t(k)

percentovershoot=100*(Y-C)/C % 计算上升时间。 n=1;

while y(n)

risetime=t(n)

% 计算稳态响应时间。 i=length(t);

while(y(i)>0.98*C)&(y(i)<1.02*C) i=i-1; end

setllingtime=t(i)

运行后的响应图如图 15,命令窗口中显示的结果为

C = timetopeak = 0.3000 1.0491 percentovershoot = risetime = 35.0914 0.6626 setllingtime = 3.5337

图 15 二阶系统阶跃响应

有兴趣的读者可以用本节介绍的游动鼠标法求取此二阶系统的各项性能指标。将它们与例18得出的结果相比较，会发现它们是一致的。 3. 利用MATLAB绘制系统根轨迹

假设闭环系统中的开环传递函数可以表示为：

- 20 -

sm?b1sm?1?????bm?1?bm(s?z1)(s?z2)???(s?zm)numGk(s)?Kn?K?K?KG0(s)n?1den(s?p1)(s?p2)???(s?pn)s?a1s?????an?1s?an 则闭环特征方程为： 1?Knum?0 den特征方程的根随参数K的变化而变化，即为闭环根轨迹。控制系统工具箱中提供了rlocus()函数，可以用来绘制给定系统的根轨迹，它的调用格式有以下几种：

rlocus(num，den) （32) rlocus(num，den，K) （33） 或者 rlocus(G) (34) rlocus(G，K) (35)

以上给定命令可以在屏幕上画出根轨迹图，其中G为开环系统G0(s)的对象模型,K为用户自己选择的增益向量。如果用户不给出K向量，则该命令函数会自动选择K向量。如果在函数调用中需要返回参数，则调用格式将引入左端变量。如

[R，K]=rlocus(G) (36)

此时屏幕上不显示图形，而生成变量R和K。

R为根轨迹各分支线上的点构成的复数矩阵，K向量的每一个元素对应于R矩阵中的一行。若需要画出根轨迹，则需要采用以下命令：

plot(R，11) （37）

plot()函数里引号内的部分用于选择所绘制曲线的类型，详细内容见表1。控制系统工具箱中还有一个rlocfind()函数，该函数允许用户求取根轨迹上指定点处的开环增益值，并将该增益下所有的闭环极点显示出来。这个函数的调用格式为：

[K，P]=rlocfind(G) (38)

这个函数运行后，图形窗口中会出现要求用户使用鼠标定位的提示，用户可以用鼠标左键点击所关心的根轨迹上的点。这样将返回一个K变量，该变量为所选择点对应的开环增益，同时返回的P变量则为该增益下所有的闭环极点位置。此外，该函数还将自动地将该增益下所有的闭环极点直接在根轨迹曲线上显示出来。

例18 已知系统的开环传递函数模型为：

Gk(s)?K?KG0(s)

s(s?1)(s?2)利用下面的MATLAB命令可容易地验证出系统的根轨迹如图16所示。 >> G=tf(1,[conv([1,1],[1,2]),0]); rlocus(G); grid

title(1Root_Locus Plot of G(s)=K/[s(s+1)(s+2)]1) xlabel(1Real Axis1) % 给图形中的横坐标命名。 ylabel(1Imag Axis1) % 给图形中的纵坐标命名。 [K,P]=rlocfind(G)

用鼠标点击根轨迹上与虚轴相交的点，在命令窗口中可发现如下结果 select_point=0.0000+1.3921i K=

5.8142

- 21 -

p=

-2.29830

-0.0085+1.3961i -0.0085-1.3961i

所以，要想使此闭环系统稳定，其增益范围应为0

参数根轨迹反映了闭环根与开环增益K的关系。我们可以编写下面的程序，通过K的变化，观察对应根处阶跃响应的变化。考虑K=0.1,0.2,?,1,2,?,5,这些增益下闭环系统的阶跃响应曲线。可由以下MATLAB命令得到。

>> hold off; % 擦掉图形窗口中原有的曲线。

图16 系统的根轨迹

t=0:0.2:15; Y=[ ];

for K=[0.1:0.1:1,2:5]

GK=feedback(K*G,1); y=step(GK,t); Y=[Y,y];

end

plot(t,Y)

对于for循环语句，循环次数由K给出。系统画出的图形如图17所示。可以看出，当K的值增加时，一对主导极点起作用，且响应速度变快。一旦K接近临界K值，振荡加剧，性能变坏。

- 22 -

图17 不同K值下的阶跃响应曲线

4. MATLAB绘图的基本知识

通过以上实例的应用，我们已初步尝试了MATLAB的绘图功能。MATLAB具有丰富的获取图形输出的程序集。我们已用命令plot()产生线性x-y图形（用命令loglog、semilogx、semilogy或polar取代命令plot,可以产生对数坐标图和极坐标图）。所有这些命令的应用方式都是相似的，它们只是在如何给坐标轴进行分度和如何显示数据上有所差别。 （１）二维图形绘制

如果用户将X和Y轴的两组数据分别在向量x和y中存储，且它们的长度相同，则命令

plot(x,y) (39)

将画出y值相对于x值的关系图。

例19 如果想绘制出一个周期内的正弦曲线，则首先应该用t=0:0.01:2*pi（pi是系统自定义的常数,可用help命令显示其定义）命令来产生自变量t；然后由命令y=sin(t)对t向量求出正弦向量y,这样就可以调用plot(t,y)来绘制出所需的正弦曲线，如图18所示。

图18一个周期内的正弦曲线

（２）一幅图上画多条曲线。

- 23 -

利用具有多个输入变量的plot( )命令，可以在一个绘图窗口上同时绘制多条曲线，命令格式为：

plot(x1，Y1，x2，Y2，?，xn，Yn) (40)

x1、Y1、x2、Y2等一系列变量是一些向量对，每一个x-y对都可以用图解表示出来，因而可以在一幅图上画出多条曲线。多重变量的优点是它允许不同长度的向量在同一幅图上显示出来。每一对向量采用不同的线型以示区别。

另外，在一幅图上叠画一条以上的曲线时，也可以利用hold命令。hold命令可以保持当前的图形，并且防止删除和修改比例尺。因此，后来画出的那条曲线将会重叠在原曲线图上。当再次输入命令hold，会使当前的图形复原。也可以用带参数的hold命令――hold on 和hold off来启动或关闭图形保持。 （３）图形的线型和颜色

为了区分多幅图形的重叠表示，MATLAB提供了一些绘图选项，可以用不同的线型或颜色来区分多条曲线，常用选项见下表１

表１MATLAB绘图命令的多种选项

选项 ′-′ ′: ′ ′r′ ′b′ ′g′ ′y′

表１中绘出的各个选项有一些可以并列使用，能够对一条曲线的线型和颜色同时作出规定。例如′--g′表示绿色的短划线。带有选项的曲线绘制命令的调用格式为：

plot(X1，Y1,S1,X2,Y2,S2,?) (41) (4)加进网络线、图形标题、x轴和y轴标记

一旦在屏幕上显示出图形，就可以依次输入以下相应的命令将网络格线、图形标题、x、y轴标记叠加在图形上。命令格式如下：

grid(网络线) (42)

title(′图形标题′) (43) xlabel（′x轴标记′） (44) ylabel（′y轴标记′） (45)

函数引号内的字符串将被写到图形的坐标轴上或标题位置。 （５）在图形屏幕上书写文本。

如果想在图形窗口中书写文字，可以单击按钮

，选择屏幕上一点，点击鼠标，在光

意义 实线 虚线 红色 蓝色 绿色 黄色 选项 ′--′ ′-.′ ′*′ ′o′ ′.′ ′×′ 意义 短划线 点划线 用星号绘制各个数据点 用圆圈绘制各个数据点 用圆点绘制各个数据点 用叉号绘制各个数据点 标处输入文字。另一种输入文字的方法是用text()命令。它可以在屏幕上以（x,y）为坐标的某处书写文字，命令格式如下：

text（x，y，′text′） （46）

- 24 -

例如，利用语句

text(3，0.45，′sint′)

将从点（３，0.45）开始，水平的写出“sint”。 （６）自动绘图算法及手工坐标轴定标

在MATLAB图形窗口中，图形的横、纵坐标是自动标定的，在另一幅图形画出之前，这幅图形作为现行图将保持不变，但是在另一幅图形画出后，原图形将被删除，坐标轴自动地重新标定。关于瞬态响应曲线、根轨迹、伯德图、奈魁斯特图等的自动绘图算法已经设计出来，它们对于各类系统具有广泛的适用性，但是并非总是理想的。因此，在某些情况下，可能需要放弃绘图命令中的坐标轴自动标定特性，由用户自己设定坐标范围，可以在程序中加入下列语句：

v=[x-min x-max y-min y-max] （47）

axis(v) (48)

式中v是一个四元向量。axis(v)把坐标轴定标建立在规定的范围内。对于对数坐标图，v的元素应为最小值和最大值的常用对数。

执行axis(v)会把当前的坐标轴标定范围保持到后面的图中，再次键入axis可恢复系统的自动标定特性。

Axis(′sguare′)能够把图形的范围设定在方形范围内。对于方形长宽比，其斜率为１

o

的直线恰位于45上，它不会因屏幕的不规则形状而变形。Axis(′normal′)将使长宽比恢复到正常状态。

５、线性系统的频域分析。 （１）频率特性函数G(j?)。

设线性系统传递函数为：

b0sm?b1sm?1?????bm?1s?bm G(s)?nna0s?a1s?????an?1s?an则频率特性函数为：

b0(j?)m?b1(j?)m?1?????bm?1(j?)?bm G(jw)?nna0(j?)?a1(j?)?????an?1(j?)?an由下面的MATLAB语句可直接求出G(jw)。

i=sqrt(-1) % 求取-1的平方根 （49） GW=polyval(num，i*w)./polyval(den，i*w) （50）

其中（num，den）为系统的传递函数模型。而w为频率点构成的向量，点右除（./）运算符表示操作元素点对点的运算。从数值运算的角度来看，上述算法在系统的极点附近精度不会很理想，甚至出现无穷大值，运算结果是一系列复数返回到变量GW中。

（２）用MATLAB作奈魁斯特图。

控制系统工具箱中提供了一个MATLAB函数nyquist( )，该函数可以用来直接求解Nyquist 阵列或绘制奈氏图。当命令中不包含左端返回变量时，nyquist（）函数仅在屏幕上产生奈氏图，命令调用格式为：

nyquist(num,den) (51) nyquist(num,den,w) (52) 或者

nyquist(G) (53)

- 25 -

nyquist(G,w) (54) 该命令将画出下列开环系统传递函数的奈氏曲线： G(s)?num(s)

den(s)如果用户给出频率向量w,则w包含了要分析的以弧度/秒表示的诸频率点。在这些频率点上，将对系统的频率响应进行计算，若没有指定的w向量，则该函数自动选择频率向量进行计算。

对于式43和45用户不必给定频率向量，系统会自动选择频率向量进行计算。式44和46需要用户给出率向量w。w包含了用户要分析的以弧度/秒表示的诸频率点,MATLAB会自动计算这些点的频率响应。

当命令中包含了左端的返回变量时，即:

[re,im,w]=nyquist(G) (55) 或

[re,im,w]=nyquist(G,w) (56) 函数运行后不在屏幕上产生图形，而是将计算结果返回到矩阵re、im和w中。矩阵re和im分别表示频率响应的实部和虚部，它们都是由向量w中指定的频率点计算得到的。

在运行结果中，w数列的每一个值分别对应re、im数列的每一个值。 例20 考虑二阶典型环节: G(s)?1

s2?0.8s?1试利用MATLAB画出奈氏图。

利用下面的命令，可以得出系统的奈氏图，如图19所示。 >> num=[0,0,1];

den=[1,0.8,1]; nyquist(num,den) % 设置坐标显示范围 v=[-2,2,-2,2]; axis(v) grid

- 26 -

title(′Nyquist Plot of G(s)=1/(s^2+0.8s+1) ′)

图19 二阶环节奈氏图

（3）用MATLAB作伯德图

控制系统工具箱里提供的bode()函数可以直接求取、绘制给定线性系统的伯德图。 当命令不包含左端返回变量时，函数运行后会在屏幕上直接画出伯德图。如果命令表达式的左端含有返回变量，bode()函数计算出的幅值和相角将返回到相应的矩阵中，这时屏幕上不显示频率响应图。命令的调用格式为：

[mag,phase,w]=bode(num,den) (57) [mag,phase,w]=bode(num,den,w) (58) 或

[mag,phase,w]=bode(G) （59） [mag,phase,w]=bode(G,w) (60) 矩阵mag、phase包含系统频率响应的幅值和相角，这些幅值和相角是在用户指定的频率点上计算得到的。用户如果不指定频率w,MATLAB会自动产生w向量，并根据w向量上各点计算幅值和相角。这时的相角是以度来表示的，幅值为增益值，在画伯德图时要转换成分贝值，因为分贝是作幅频图时常用单位。可以由以下命令把幅值转变成分贝：

magdb=20﹡log10(mag) (61)

绘图时的横坐标是以对数分度的。为了指定频率的范围，可采用以下命令格式： logspace(d1，d2) （62) 或

logspace(d1，d2，n) （63）

公式（62）是在指定频率范围内按对数距离分成50等分的，即在两个十进制数?1?10和

2d1?2?10d之间产生一个由50个点组成的分量，向量中的点数50是一个默认值。例如

要在?1?0.1弧度/秒与?2?100弧度/秒之间的频区画伯德图，则输入命令时，

- 27 -

返回到工作d1?log10(?1),d2?log10(?2)在此频区自动按对数距离等分成50个频率点，空间中，即

w=logspace(-1，2)

要对计算点数进行人工设定，则采用公式（63）。例如，要在?1?1与?2?1000之间产生100个对数等分点，可输入以下命令：

w=logspace(0，3，100)

在画伯德图时，利用以上各式产生的频率向量w，可以很方便地画出希望频率的伯德图。

由于伯德图是半对数坐标图且幅频图和相频图要同时在一个绘图窗口中绘制，因此，要用到半对数坐标绘图函数和子图命令。

（1） 对数坐标绘图函数

利用工作空间中的向量x，y绘图，要调用plot函数，若要绘制对数或半对数坐标图，只需要用相应函数名取代plot即可，其余参数应用与plot完全一致。命令公式有：

semilogx(x，y，s) (64)

上式表示只对x轴进行对数变换，y轴仍为线性坐标。

semilogy(x，y，s) (65)

上式是y轴取对数变换的半对数坐标图。

Loglog(x，y，s) (66)

上式是全对数坐标图，即x轴和y 轴均取对数变换。 （2） 子图命令

MATLAB允许将一个图形窗口分成多个子窗口，分别显示多个图形，这就要用到subplot()函数，其调用格式为：

subplot(m，n，k)

该函数将把一个图形窗口分割成m×n个子绘图区域，m为行数，n为列数，用户可以通过参数k调用各子绘图区域进行操作，子图区域编号为按行从左至右编号。对一个子图进行的图形设置不会影响到其它子图，而且允许各子图具有不同的坐标系。例如，subplot(4，3，6)则表示将窗口分割成4×3个部分。在第6部分上绘制图形。 MATLAB最多允许9×9的分割。

例21 给定单位负反馈系统的开环传递函数为：

G(s)?10(s?1)

s(s?7)试画出伯德图。

利用以下MATLAB程序，可以直接在屏幕上绘出伯德图如图20。 >> num=10*[1,1]; den=[1,7,0]; bode(num,den) grid

title(′Bode Diagram of G(s)=10*(s+1)/[s(s+7)] ′) 该程序绘图时的频率范围是自动确定的，从0.01弧度/秒到30弧度/秒，且幅值取分贝值，?轴取对数，图形分成2个子图，均是自动完成的。

- 28 -

图20 自动产生频率点画出的伯德图

如果希望显示的频率范围窄一点，则程序修改为： >> num=10*[1,1]; den=[1,7,0];

w=logspace(-1,2,50); % 从0.1至100，取50个点。 [mag, phase, w]=bode(num, den, w);

magdB=20*log10(mag) % 增益值转化为分贝值。

% 第一个图画伯德图幅频部分。 subplot(2,1,1);

semilogx(w,magdB, ′-r′) % 用红线画 grid

title(′Bode Diagram of G(s)= 10*(s+1)/[s(s+7)] ′) xlabel(1Frequency(rad/s)1) ylabel(1Gain(dB)1)

% 第二个图画伯德图相频部分。 subplot(2,1,2);

semilogx(w,phase, 1-r1); grid

xlabel(1Frequency(rad/s)1) ylabel(′Phase(deg) ′)

- 29 -

修改程序后画出的伯德图如21所示。

图21 用户指定的频率点画出的伯德图

4. 用MATLEB求取稳定裕量

同前面介绍的求时域响应性能指标类似，由MATLAB里bode()函数绘制的伯德图也可以采用游动鼠标法求取系统的幅值裕量和相位裕量。例如，我们可以在图20的幅频曲线上按住鼠标左键游动鼠标，找出纵坐标（Magnitude）趋近于零的点，从提示框图中读出其频率约为7.25dB。然后在相频曲线上用同样的方法找到横坐标（Frequence）最接近7.25dB的点，可读出其相角为-53.9度，由此可得，此系统的相角裕量为126.1度。幅值裕量的计算

方法与此类似。

此外，控制系统工具箱中提供了margin()函数来求取给定线性系统幅值裕量和相位裕量，该函数可以由下面格式来调用：

[Gm, Pm, Wcg, Wcp]=margin(G); (67)

可以看出，幅值裕量与相位裕量可以由LTI对象G求出，返回的变量对（Gm, Wcg）为幅值裕量的值与相应的相角穿越频率，而（Pm, Wcp）则为相位裕量的值与相应的幅值穿越频率。若得出的裕量为无穷大，则其值为Inf，这时相应的频率值为NaN(表示非数值)，Inf和NaN均为MATLAB软件保留的常数。

如果已知系统的频率响应数据，我们还可以由下面的格式调用此函数。

[Gm, Pm, Wcg, Wcp]=margin(mag, phase, w);

其中(mag, phase, w)分别为频率响应的幅值、相位与频率向量。

例22 已知三阶系统开环传递函数为：

G(s)?7 322(s?2s?3s?2)利用下面的MATLAB程序，画出系统的奈氏图，求出相应的幅值裕量和相位裕量，并求出闭环单位阶跃响应曲线。

- 30 -

>> G=tf(3.5,[1,2,3,2]);

subplot(1,2,1); % 第一个图为奈氏图 nyquist(G); grid

xlabel('Real Axis') ylabel('Imag Axis') % 第二个图为时域响应图 [Gm,Pm,Wcg,Wcp]=margin(G) G_c=feedback(G,1); subplot(1,2,2); step(G_c) grid

xlabel(′Time(secs) ′) ylabel(′Amplitude′) 显示结果为：

ans=1.1429 1.1578 1.7321 1.6542

图22 三阶系统的奈氏图和阶跃响应图

画出的图形如图22 所示。由奈氏曲线可以看出，奈氏曲线并不包围（-1，j0）点，故闭环系统是稳定的。由于幅值裕量虽然大于1，但很接近1，故奈氏曲线与实轴的交点离临

o

界点（-1，j0）很近，且相位裕量也只有7.1578，所以系统尽管稳定，但其性能不会太好。观察闭环阶跃响应图，可以看到波形有较强的振荡。

o

如果系统的相角裕量γ>45,我们一般称该系统有较好的相角裕量。 例23 考虑一个新的系统模型，开环传递函数为：

100(s?5)2 G(s)?2(s?1)(s?s?9)- 31 -

由下面MATLAB语句可直接求出系统的幅值裕量和相位裕量:

>> G=tf(100*conv([1,5],[1,5]), conv([1,1],[1,1,9])); [Gm, Pm, Wcg, Wcp]=margin(G)

结果显示 Gm = Pm =

Inf 85.4365

Wcg = Wcp =

NaN 100.3285 再输入命令

>> G_c=feedback(G,1);

step(G_c) grid

xlalel(′Time(sec) ′) ylalel(′Amplitude′)

o

可以看出，该系统有无穷大幅值裕量，且相角裕量高达85.4365。所以系统的闭环响应是较理想的，闭环响应图如图23.

图23 较理想的系统响应

5 时间延迟系统的频域响应

(1) 时间延迟系统的传递函数模型

-Ts

带有延迟环节e的系统不具有有理函数的标准形式，在MATLAB中，建立这类系统的模型。要由一个属性设置函数set()来实现。该函数的调用格式为：

set(H, ′属性名′, ′属性值′) （68）

其中H为图形元素的句柄（handle）。在MATLAB中，当对图形元素作进一步操作时，只需对该句柄进行操作即可。例如以下调用格式

h=plot(x,y) G=tf(num,den)

Plot()函数将返回一个句柄h，tf()函数返回一个句柄G，要想改变句柄h所对应曲线的颜色，则可以调用下面命令：

Set(h,color,[1,0,0]);

即对“color”参数进行赋值，将曲线变成红色（由[1,0,0]决定）

- 32 -

同样，要想对G句柄所对应模型的延迟时间’Td’进行修改，则可调用下面命令

Set(G, ′Td′,T)

其中T为延迟时间。由此修改后，模型G就已具有时间延迟特性。 (2) 时间延迟系统的频域响应

含有一个延迟环节的系统，其开环频域响应为 G(j?)e?jT??G(j?)ej?(?)?T?

可见，该系统的幅频特性不变，只加大了相位滞后。 例24 考虑系统的开环模型为: G(s)?1?Tse s?1当T=1时，我们可以由下面的MATLAB命令绘出系统的奈氏图，如图24所示,此系统对应的时域响应图为25。

>> G=tf(1,[1, 1]); T=[1];

w=[0, logspace(-3, 1, 100), logspace(1,2,200)]; set(G,‘Td', T); % 延迟1秒。 nyquist(G,w) grid

figure % 建立一个新的绘图窗口 step(G)

图24 时间延迟系统奈氏图

图25 时间延迟系统的阶跃响应

4 频域法串联校正的MATLAB方法

利用MATLAB可以方便的画出Bode图并求出幅值裕量和相角裕量。将MATLAB应用到经典理论的校正方法中，可以方便的校验系统校正前后的性能指标。通过反复试探不同校正参数对应的不同性能指标，能够设计出最佳的校正装置。

例25 给定系统如图26 所示，试设计一个串联校正装置，使系统满足幅值裕量大于

o

10分贝，相位裕量≥45

- 33 -

解：为了满足上述要求，我们试探地采用超前校正装置Gc(s)，使系统变为图27的结构。

图26 校正前系统

图27 校正后系统

我们可以首先用下面地MATLAB语句得出原系统的幅值裕量与相位裕量。 >> G=tf(100, [0.04, 1, 0]);

[Gw, Pw, Wcg, Wcp]=margin(G); 在命令窗口中显示如下结果

w = Pw =

Inf 28.0243

Wcg = Wcp =

Inf 46.9701

o

可以看出，这个系统有无穷大的幅值裕量，并且其相位裕量?＝28，幅值穿越频率Wcp=47rad/sec。

引入一个串联超前校正装置：

Gc(s)?0.025s?1

0.01s?1我们可以通过下面的MATLAB语句得出校正前后系统的Bode图如图28,校正前后系统的阶跃响应图如图29。其中?1、?1、ts1分别为校正前系统的幅值穿越频率、相角裕量、调节时间，?2、?2、ts2分别为校正后系统的幅值穿越频率、相角裕量、调节时间。

>> G1=tf(100,[0.04,1,0]); % 校正前模型

G2=tf(100*[0.025,1],conv([0.04,1,0],[0.01,1])) % 校正后模型

- 34 -

% 画伯德图,校正前用实线，校正后用短划线。 bode(G1) hold

bode(G2, ′--′)

% 画时域响应图，校正前用实线，校正后用短划线。 figure

G1_c=feedback(G1,1) G2_c=feedback(G2,1) step(G1_c) hold

step(G2_c, ′--′)

图28 校正前后系统的Bode图

图29 校正前后系统的阶跃响应图

可以看出，在这样的控制器下，校正后系统的相位裕量由28增加到48，调节时间由

o

。0.28s减少到0.08s。系统的性能有了明显的提高，满足了设计要求。

5 自动控制理论模拟实验

《自动控制理论》是一门理论性和实践性很强的专业基础课，前面我们通过计算机仿真，可以方便地研究系统性能，验证理论的正确性，加深对理论知识的理解。本节我们再通过电子模拟实验，学习和掌握系统模拟电路的构成和测试技术，进一步培养学生的实际动手能力和分析、研究问题的能力。

在控制理论课程中，大部分院校目前拥有的实验设备是电子模拟学习机。这种专为教学实验制造的电子模拟学习机，体积较小，使用方便，实验箱中备有多个运算放大器构成的独立单元，再加上常用的电阻、电容等器件，通过手工连线、可以构成多种特性的被控对象和控制器。

在基础训练阶段，实验手段采用模拟方法，除了灵活方便之外，还具有以下两个优点： （1） 电子模拟装置可建立较准确的数学模型，从而可以避免实际系统中常碰到的各种复

杂因素，使初学者能够根据所学理论知识循序渐进地完成各项实验。

（2） 在工程实践中，电子模拟方法有一定的实用价值，也是实验室常用的一种实验方法。

当然，对于将来从事实际工作的学生来说，仅仅掌握模拟实验方法还是不够的，应

- 35 -

在此基础上进行一些以实际系统为主要设备的实验训练。 1. 常用的实验设备和仪器

以自控理论电子模拟学习机为核心的一组基本实验设备和仪器，共同完成对各种实验对象的模拟和测试任务，传统的测试手段下，构成基本实验必备仪器有以下几种：

（1） 电子模拟学习机。 （2） 超低频双踪示波器。 （3） 超低频信号发生器。 （4） 万用表。

按照被测系统的数学模型，在电子模拟学习机上用基本运放单元模拟出相应的电路模型，然后按图30所示的方法进行模拟实验测试。

图30 传统仪器组合

随着计算机软、硬件的快速发展，人们越来越多地利用计算机实现的虚拟仪器代替传统仪器。目前，大多数实验室都是用计算机来实现信号的产生、测量与显示、系统的控制及数据处理，使实验过程更加方便，功能更强大。现在的模拟实验组件是按图31来实现的。

图31 计算机仿真模拟实验

A/D、D/A卡起模拟信号与数字信号的转换作用，还可产生不同的输入信号（阶跃、三角、正弦等），供实验时选用。使用时用RS232串口电缆将A/D、D/A卡与计算机连接起来。如果配备打印机，则可在实验的同时将实验结果打印输出。由于计算机可以方便地输入数据、观察数据，初学者可以在屏幕的提示下进行实验过程，使学习变的更加轻松。

2. 自控理论的基本实验

实验一. 典型环节及阶跃响应测试

1. 实验的基本原理

控制系统的模拟实验是采用复合网络法来模拟各种典型环节，即利用运算放大器和RC

- 36 -

组成的不同输入网络和反馈网络组合，模拟出各种典型环节，然后按照给定系统的结构图将这些模拟环节连接起来，便得到了相应的模拟系统。然后将输入信号加到模拟系统的输入端，使系统产生动态响应。这时，可利用计算机或示波器等测试仪器，测量系统的输出，便可观测到系统的动态响应过程，并进行性能指标的测量。若改变系统的某一参数，还可进一步分析研究参数对系统性能的影响。

在以下的实验过程中，为了更好地检验实验结果，避免过多地出现错误操作，我们将每一环节的正确结果，通过Simulink仿真软件绘出正确的图形，以便于读者检验实验结果的正确性。

2. 时域性能指标测量方法 （1）最大超调量?%

利用示波器或计算机显示器上测到的输出波形，读出响应最大值和稳态值所具有的刻度值，代入下式算出超调量：

?%?Ymax?Y??100% （60） Y?（2）峰值时间tp

根据示波器或显示器上输出的波形最大值，找出这一点在水平方向上所具有的刻度值，即可换算出或读出峰值时间tp。

（3）调节时间ts

同样，读出水平方向上对应输出从零到进入5％或2％误差带时所占的刻度值，即可得到调节时间ts。

3. 实验内容

（1）比例环节 G(s)??K 模拟线路如图32

图32 比例放大电路

RRC(s)??2，K?2 R(s)R1R1一般可通过改变电阻R2来调整放大倍数K。

- 37 -

由于输入信号r(t)是从运算放大器的反相输入端输入，所以输出信号和输入信号在相位上正好相反，传递函数中出现负号。为了观测方便，可以从输入端输入负阶跃信号。也可以在输出端连接一个反相器，如图33。

图33

取R1?100K，R2?200K，将模拟学习机上手动阶跃信号（或信号发生器置于“手动阶跃”）引入环节输入端，观测输出波形，并作记录。（为便于比较，应将输入信号与输出信号同时送入双踪示波器或计算机，两路信号同时在一个坐标系下显示。绘制曲线时，也用这种形式）。

（a）仿真模块 (b) 仿真输出

图34 比例环节的Simulink仿真

图34为Simulink 的仿真模块。为便于观察，阶跃信号输入时间设置为1s（系统默认值），后面的各个例题也都适当调整输入时间。增益（Gain）模块的增益放大倍数设为2。另外，也可以用鼠标双击各模块，设置适合其它参数。

- 38 -

（2）积分环节 G(s)?1 Ts

模拟线路如图35 所示。

C(s)?1, T?R1C ?R(s)R1CS

图35 积分电路

积分时间常数可通过改变电阻R1或电容C来选择。

取R1?100K，C?1?f，按上述同样方法观测阶跃响应波形。用Simulink仿真的环节模块图如图36 (a)。由于积分环节附带的增益比较大（积分时间常数T＝0.1），Scope(示波器)绘出图形的辐值显示范围并不是很理想。我们可以在Scope的显示图中点击鼠标右键，选Axes properties菜单，在弹出的对话框中设置Y-max属性为100，则输出结果如图36 (b)所示。

（a）仿真模块 (b) 仿真输出

图36 积分环节的Simulink仿真

（3）微分环节 G(s)??Ts 模拟线路图如图37 。

R2C1sC(s)????R2C1s R(s)R2C2s?1其中： C2??C1，T?R2C1

图37 微分电路

微分时间常数T可通过改变R2和C1来选

- 39 -

取。令R2?100K，C1?1?f，C2?0.01?f，按上述同样步骤进行模拟和测试，观察微分环节的阶跃响应波形。用Simulink仿真的模块图为图38(a),在Scope绘出的图形中调整横纵坐标，得出的时域响应图如图38（b）所示。 （4）惯性环节

?K

G(s)?Ts?1

（a）仿真模块 （b）仿真输出

图38 微分环节的Simulink仿真

模拟线路如图39。

图39 惯性环节电路图

R2RR1C(s)，K?2，T?R2C ?R1R(s)R2Cs?1?取R1?100KR2?200K，C?1?f，观测其阶跃响应输出，测出ts，并与理论值

ts?3T（或4T）相比较，用Simulink仿真结果如图40。

- 40 -

（a）仿真模块 （b）仿真输出

图40 惯性环节的Simulink仿真

（5）振荡环节

2??n G(s)?2 2s?2??ns??n其模拟线路如图41 。

图41 振荡环节

对应的结构图如图42。

图42 振荡环节结构图

- 41 -

C(s)?R(s)?(s2?12)R1C1112s?()R2C1R1C1，?n?R1，??1 R1C12R2改变C1可改变?n的大小，改变R2可改变?的大小。按表2给出的参数测量阶跃响应，并记入表中。

表2 不同的?n和?所对应的时域性能指标 记 录 参 数 ?p% 100% tp(ms) ts(ms) 阶跃响应波形 R2?? ??0 R2?200K ?n?10rad/s （R1?100K ??0.25 R2?100K C1?1?f） ??0.5 R2?50K ??1 R2?200K ?n?100rad/s （R1?100K ??0.25 R2?100K C1?0.01?f）

??0.5 当?n?10，??0.25时用Simulink仿真结果如图43。

- 42 -

（a）仿真模块 （b）仿真输出

图43 振荡环节的Simulink 仿真

4. 实验报告

① 画出实验线路图，记录原始数据、测试数据及波形。 ② 讨论惯性环节（一阶系统）：

按实验给出的响应曲线，求出K，T，ts，与理论计算值比较，并得出由实验结果求惯性环节传递函数的方法。

③ 讨论振荡环节（二阶系统）：

按实验给出的欠阻尼下的响应曲线，求出tp和ts，与理论值相比较，并确定参数?，?n，最终可得出由实验结果求振荡环节传递函数的方法。

讨论振荡环节性能指标与?，?n的关系。

实验二 系统频率特性测量

利用简单仪器测量频率特性，测量精度是较差的，但物理意义明显，波形直观是其特点。本实验通过“李沙育图形”法进行频率特性测试，可以使学生通过实验观测到物理系统的频率响应，并根据测量值算出频率特性的幅值和相角，通过实验可以掌握测试频率特性的基本原理和方法。

1. 原理

一个稳定的线性系统，在正弦信号的作用下，它的稳态输出将是一个与输入信号同频率的正弦信号，但其振幅和相位却随输入信号的频率不同而变化，测取不同频率下系统的输出、输入信号的振幅比及相位差?，即可求得这个系统的幅频特性和相频特性。

设线性系统输入和稳态响应分别为以下两个正弦信号：

?x(?t)?Xmsin?t ?y(?t)?Ysin(?t??)m?- 43 -

幅频特性 |G(j?|?Ym(?) （69）

Xm(?)

（70）

相频特性 ?G(j?)??(?)若以x(?t)为横轴，以y(?t)为纵轴，而以?t作为参变量，则随?t的变化，x(?t)和

y(?t)所确定的点的轨迹，将在x－y平面上描绘出一条封闭的曲线（通常是一个椭圆）。这

就是所谓的“李沙育图形”，如图44。

图44 “李沙育图形”原理图

不断的改变x(?t)的频率，就可以获得一系列形状不同的李沙育图形。由此求出各个频率所对应的相位差和幅值比，就可获得系统的频率特性。

幅值比由测量数据按式（69）直接求出；而相位差的具体求法如下：

令?t＝0，则 ??x(0)?0

?y(0)?Ymsin?即得

??sin?102yy(0)?sin?10 （71） Ym2Ym0显然上式仅当0???90时成立，“李沙育图形”在四个象限的形状如图45 所示，注意箭头方向。

实际的控制系统一般为相位滞后系统，即频率特性的相频是负的角度，相频特性滞后角

- 44 -

按“李沙育图形”法，应按下式确定：

第四象限： （逆时针） 第三象限： ?4??sin?12y0 （72） 2Ym2y0) （73） ?3??(??sin?12Ym（逆时针） 第二象限： ?12y02??(??sin?2Y) m（顺时针） 第一象限： ?2y01??(2??sin?12Y) m（顺时针）

图45 “李沙育图形”形状

2. 实验内容

（1）给出三阶系统模拟电路如图46 所示

- 45 -

74） （75）

（

图46 三阶系统模拟图 对应的系统结构图如图47 所示。

图47 系统结构图

选元件R1?100K,C1?1?f，R2?R3?100K，C2?C3?0.1?f 则开环传递函数为：

Gk(s)?10 2(0.1s?1)(0.001s?0.1s?1)（2）断开闭环系统模拟电路图46中主反馈线路，按开环三阶系统在学习机上接好线路，并将有关测试仪器按图48连接。

图48 测试电路

- 46 -

将超低频正弦输入信号输入系统，调节输入信号幅度使被测对象在避免饱和的情况下，输出幅度尽可能大，以便于测量。然后调节示波器Y轴增益（量程范围），使在所取信号幅度下，图像达到满刻度。

（3）在示波器上测量此时输入信号幅值（用2Xm表示），并记录在表3中，此后在输出幅度能有效测出时，一般不再改变输入信号的幅度。

按表中给定的测点依次改变输入信号频率，测试并记录于表3中。 为了提高读数精度，对示波器的X，Y轴增益可随时调节，以获得较好的“李沙育图形”。注意在X轴与Y轴增益不一致时，“李沙育图形”的形状可能会有所变化。读数后按相应的增益正确折算出Xm，Ym值。另外，在转折频率附近以及穿越频率?c附近应多测几点。

表3 频率特性测试结果记录 ?(1s) f(Hz) 2Xm 2Ym 8 1.27 10 1.6 15 2.4 20 3.2 30 4.8 40 6.4 60 9.6 80 12.7 100 16 200 32 20lg2Ym(db)2Xm 2y0 sin?12y0 2Ym ? 李沙育图形 3. 实验报告

（1） 按被测对象的传递函数，画出模拟电路图

（2） 整理表3中的实验数据，在半对数坐标纸上作出被测系统的对数幅频特性和相频特

性。

（3） 采用MATLAB语言，画出被测对象的Bode图，与上图进行比较，特别验证在测试

点处的结果是否一致。

（4） 讨论“李沙育图形”法测试频率特性的优缺点，有效频率范围及测试精度。 附：（1）用以下MATLAB命令绘制的开环系统伯德图见图49。

>> G=tf(10,[conv([0,1,1],[0.001,0.1,1])]); Bode(G)

- 47 -

grid

49 开环系统伯德图

（2）采用以下MATLAB命令绘出的“李沙育图形”见图50 >> G=tf(10,[conv([0.1,1],[0.001,0.1,1])]); %画频率为8的图。

t=0:0.01:2*pi; u=sin(8*t); y=lsim(G,u,t); plot(u,y)

%画频率为60的图。

t=0:0.1:2*pi; % 设置新的时间间隔。 figure

u=sin(60*t); y=lsim(G,u,t); plot(u,y)

其中，MATLAB函数lsim()是求任意输入下的响应，调用格式与step()函数基本一致，只是输入变量中必须包含任意输入u，向量u表示系统输入在各个时刻的值。

如果想求取y0和ym，可输入如下程序 >> i=length(u);

while(u(i)>0&u(i)<0) i=i-1; end

y0=abs(y(i)) ym=max(y)

- 48 -

对频率为8时给定的u和y，可得如下结果: y0 = ym =

8.9149 9.2467

（a）?=8时的“李沙育图形” （b）?=60时的“李沙育图形”

图50 MATLAB仿真“李沙育图形”

（3）应用Simulink仿真工具输出的响应曲线及画出的“李沙育图形”如图51所示。

图51（a） Simulink仿真图

图51（b）响应波形

图51（c）Simulink绘出的“李沙育图形”

- 49 -

在图51(a)所示的编辑窗口中，点击Simulation|Simulation parameters菜单，将Simulation time设为：Start time设为0.0s，Stop time 设为3.0s。双击Sine Wave控件，将其Frequency设为8，Amplitude设为0.5，双击XY Graph控件，将x-min，x-max，y-min, y-max, sample time(采样时间间隔)分别设为-0.5，0.5，-4，4，-1。输出响应波形为图51（b），XY Graph绘出的“李沙育图形”为51(c).

实验三 连续系统的频率法串联校正

频率法串联校正，主要是根据工程上提出的频率指标?，?c ，Kg等，在频率特性曲线上进行计算测量，最终得出欲增加的校正装置的功能，然后验证校正后的实际结果，经过反复调整，直至满意为止。

本实验的校正原理已在前面章节介绍，理论计算过程可以由读者自己完成。本节只通过对校正前后系统性能的测量，来进一步体会串联校正的实际效果。

1. 实验目的

（1）加深理解串联校正装置对系统动态性能的校正作用。

（2）对给定系统进行串联校正设计，并通过模拟实验检验设计的正确性。 2. 实验内容

对于给定的单位反馈闭环系统如图52 ，其对应的开环传递函数为：

Gk?其中 K?K1K2K3

K

s(Ts?1)

图52 原闭环系统结构图

（1） 串联超前校正

系统的模拟电路图如图53 所示，图中开关S断开时对应校正前的系统，当S接通时串入超前校正环节。

- 50 -

本文来源：https://www.bwwdw.com/article/8nzo.html