机电控制技术试验指导书(基于matlab的控制系统设计与仿真)--大学毕业设计论文.doc

机电控制技术实验指导书
基于MATLAB的控制系统设计与仿真

王正杰 杨喆

北京理工大学机电学院
2012年4月

目录
1. MATLAB简介 1
1.1 MATLAB界面 1
1.2 文件和变量的命名规则 1
1.3 标点符号 2
1.4 数值计算 2
1.4.1 复数 2
1.4.2 矩阵和数组 3
1.4.3 基本运算函数 4
1.5 控制语句 5
1.5.1 循环语句 5
1.5.2 条件转移语句 5
1.5.3 break语句 6
1.6 函数绘图 6
1.6.1 曲线的颜色、线型和数据点形 6
1.6.2 绘制二维曲线 6
1.6.3 绘制三维图形 6
1.6.4 多个图形的绘制 7
1.6.5 交互式图形命令 8
1.6.6 图形的标注 8
2. 控制系统的数学模型 9
2.1 Laplace变换和反变换 9
2.1.1 直接调用MATLAB所提供的函数 9
2.1.2部分分式展开 9
2.2 传递函数描述法 11
2.3 零极点描述法 11
2.4 模型之间的转换 12
2.5 结构图的化简 13
2.5.1 串联结构 13
2.5.2 并联结构 13
2.5.3 反馈结构 14
3. 控制系统的时域分析 16
3.1系统稳定性分析 16
3.2 系统动态特性分析 17
4. 控制系统的频域分析 23
4.1 频域特性 23
4.2 绘制Bode图 24
4.3 绘制Nyquist曲线 29
4.4 幅值裕度和相位裕度 31
5. 控制系统的根轨迹分析 34
5. 1 绘制根轨迹 34
5.2 开环零点和极点对系统性能的影响 39
5.2.1 开环零点对系统性能的影响 39
5.2.2 开环极点系统性能的影响 42
5.2.3 开环偶极子对系统性能的影响 45

II

1. MATLAB简介
1.1 MATLAB界面
（1）命令窗口（Command Window）
进行各种MATLAB指令各种操作，可输入各种MATLAB指令、函数、表达式、变量等，并显示除图形外的所有运算结果。
（2）工作空间窗口（Workspace）
显示目前保存在内存中所有变量的变量名、类型，最大值和最小值等。
（3）当前目录窗口（Current Directory或Current Folder）
显示当前路径。
（4）历史命令窗口（Command History）
记录在命令窗口中已经运行过的指令、函数。
1.2 文件和变量的命名规则
（1）文件名和变量名必须是数字或英文字母，第一个字符不能是数字，最多可包含31个字符；
（2）文件名和变量名中不得包含空格、标点，可以包含下划线；
（3）文件名和变量名需区分大小写；
（4）文件名和函数名不得与MATLAB固有函数、特殊变量和常量名重复；
（5）文件存储路径必须为英文。
表1.1MATLAB的特殊变量与常量
变量名
取值
ans
默认变量名，以
pi
圆周率π
i或j
虚数单位，i=j=
inf
无穷大，如
NaN或nan
不定值，如，，
eps
计算机的最小数
realmax
最大的正实数
realmin
最小的正实数
nargin
函数的输入变量数目
narout
函数的输出变量数目
1.3 标点符号
表1.2 MATLAB中的常用标点符号
名称
符号
功能
空格

用于输入变量之间的分隔符以及数组行元素之间的分隔符
逗号
,
用于要显示计算结果的命令之间的分隔符；用于输入变量之间的分隔符；用于数组行元素之间的分隔符
分号
;
用于不显示计算结果命令行的结尾；用于不显示计算结果命令之间的分隔符；用于数组元素行之间的分隔符
冒号
:
用于生成一维数值数组，表示一维数组的全部元素或多维数组的某一维的全部元素
百分号
%
用于注释的前面，在它后面的命令不需要执行
单引号
‘ ’
用于括住字符串
圆括号
( )
用于引用数组元素；用于函数输入变量列表；用于确定算术运算的先后次序
方括号
[ ]
用于构成向量和矩阵；用于函数输出列表
下划线
-
用于一个变量、函数或文件名中的连字符
续行号
…
用于把后面的行与该行连接以构成一个较长的命令
1.4数值计算
1.4.1 复数
复数由实部和虚部组成，MATLAB用特殊变量“i”和“j”表示虚数的单位。复数运算不需要特殊处理，可以直接进行。
（1）复数的表示
z=a+b*i或z=a+b*j
z=a+bi或z=a+bj(当b为标量时)
z=r*exp(i*theta)
（2）求一个复数的实部、虚部、幅值和相角
a=real(z) %计算实部
b=imag(z) %计算虚部
r=abs(z) %计算幅值
theta=angle(z) %计算相角
说明：
复数z的实部a=r*cos(θ)；
复数z的虚部b=r*sin(θ)；
复数z的幅值；
复数z的相角，以弧度为单位。
1.4.2 矩阵和数组
矩阵元素应用方括号括住，每行内的元素间用逗号或空格隔开，行与行之间用分号或回车键隔开，矩阵元素可以是数值或表达式。
（1）通过显式元素列表输入矩阵
表1.3通过显式元素列表输入矩阵
说明
例子
输入
结果
[ ]表示构成矩阵，分号分隔行，空格分隔元素
c=[1 2;3 4;5 3*2]
或c=[1,2;3,4;5,3*2]
c =
1 2
3 4
5 6
（2）通过语句生成矩阵
表1.4通过语句生成矩阵
函数名
说明
例子
输入
结果
from:step:to
from、step和to分别表示开始值、步长和结束值，step省略则默认为step=1
x1=1:2:6
x1 =
1 3 5
linspace(a,b,n)
a、b、n三个参数分别表示开始值、结束值和元素个数，生成从a到b之间线性分布的n个元素的行向量。n如果省略则默认值为100
x2=linspace(0,2*pi,3)
x2 =
0 3.1416 6.2832
logspace(a,b,n)
a、b、n三个参数分别表示开始值、结束值和数据个数，n如果省略则默认值为50。生成从10a到10b之间按对数等分的n个元素的行向量
x3=logspace(0,2,3)
x3 =
1 10 100
（3）由矩阵生成函数产生特殊矩阵
表1.5由矩阵生成函数产生特殊矩阵
函数名
功能
例子
输入
结果
zeros(m,n)
产生m×n的全0矩阵
zeros(2,3)
ans =
0 0 0
0 0 0
ones(m,n)
产生m×n的全1矩阵
ones(2,3)
ans =
1 1 1
1 1 1
rand(m,n)
产生均匀分布的随机矩阵，元素取值范围0.0～1.0
rand(2,3)
ans =
0.9501 0.6068 0.8913
0.2311 0.4860 0.7621
randn(m,n)
产生正态分布的随机矩阵
randn(2,3)
ans =
-0.4326 0.1253 -1.1465
-1.6656 0.2877 1.1909
magic(N)
产生N阶魔方矩阵（矩阵的行、列和对角线上元素的和相等）
magic(3)
ans =
8 1 6
3 5 7
4 9 2
eye(m,n)
产生m×n的单位矩阵
eye(3)
ans =
1 0 0
0 1 0
0 0 1
1.4.3 基本运算函数
表1.6MATLAB的基本运算函数
函数名
函数功能
函数名
函数功能
函数名
函数功能
sin
正弦运算
asin
反正弦
log
以e为底的对数
cos
余弦运算
acos
反余弦
log2
以2为底的对数
tan
正切运算
atan
反正切
log10
以10为底的对数
cot
余切运算
acot
反余切
sqrt
开平方
exp
指数函数
abs
求模
round
就近取整
floor
向下取整
ceil
向上取整
sum
求和
1.5控制语句
1.5.1 循环语句
（1）for
for 循环变量=表达式1:表达式2:表达式3
循环语句组
end
其中，表达式1,2,3分别为循环变量的初始值、增量、最终值。当增量为1时，表达式2可省略。
（2）while
while (条件式)
循环语句组
end
对于条件式来说，逻辑运算是必要的，MATLAB中使用“&”，“|”，“~”分别表示“与”，“或”，“非”。
1.5.2 条件转移语句
（1）if (条件式)
语句组
end
（2）if(条件式)
语句组1
else
语句组2
end
（3）if(条件式1)
语句组1
else if (条件式2)
语句组2
end
1.5.3 break语句
break用于终止包含它的最里面一层循环过程。
1.6函数绘图
1.6.1 曲线的颜色、线型和数据点形
表1.7曲线的颜色、线型和数据点形
选项
意义
选项
意义
选项
意义
y
黄色
-
实线（默认）
.
实点标记
m
紫色
：
点线
o
圆圈标记
c
青色
-.
点划线
x
叉号形
r
红色
--
虚线
+
十字形
g
绿色

*
星号形
b
蓝色

s
方框形
k
黑色

d
钻石形

p
五角星形
在系统默认的情况下，第一条曲线是蓝色，第二条是绿色，第三条是红色，第四条是青色，第五条是紫色，第六条是黄色，第七条是黑色，如此循环。
1.6.2绘制二维曲线
（1）绘制向量x和y的曲线
plot(x,’s’)%绘制以x为纵坐标的二维曲线，s是指定线型、色
%彩、数据点形的字符串。
plot(x,y,’s’) %绘制以x为横坐标y为纵坐标的二维曲线
（2）绘制矩阵x的曲线
plot(x,’s’) %当x为矩阵时，plot(x)为每个列向量画一条曲线
（3）绘制复向量的曲线
plot(z,’s’) %以复数z的实部为横坐标，虚部为纵坐标
（4）绘制多条曲线
plot(x1,y1,’s1’,x2,y2,’s2’,…) % MATLAB自动以不同的颜色绘制不同曲线
1.6.3绘制三维图形
（1）绘制三维曲线——plot3命令
plot3(x,y,z, 's') %绘制三维曲线
plot3(x1,y1,z1, 's1',x2,y2,z2, 's2',…) %绘制多条三维曲线
（2）绘制三维网线——mesh命令
mesh(z) %画三维网线图
mesh(x,y,z)
当只有参数z时，以z矩阵的行下标作为x坐标轴，把z的列下标当作y坐标轴；x、y分别为x、y坐标轴的自变量；当有x、y、z参数时，z是指定各点的数据。
1.6.4多个图形的绘制
（1）指定图形窗口
figure(n) %产生新图形窗口
如果该窗口不存在，则产生新图形窗口并设置为当前图形窗口，该窗口名为“Figure No.n”，不关闭其它窗口。
（2）同一窗口多个子图形
将plot命令改成subplot命令，将一个图形窗口划分为多个区域，每个区域一幅子图，即
subplot(m,n,k) %使(m×n)幅子图中的第k幅成为当前图
将图形窗口划分为m×n幅子图，k是当前子图的编号。子图的序号编排原则是：左上方为第1幅，先向右后向下依次排列，子图彼此之间独立。
（3）同一窗口多次叠绘
为了在一个坐标系中增加新的图形对象，可以用“hold”命令来保留原图形对象。
hold on %使当前坐标系和图形保留
hold off %使当前坐标系和图形不保留（默认）
hold %在以上两个命令中切换
在设置了“hold on”后，如果画多个图形对象，则在生成新的图形时保留当前坐标系中已存在的图形对象，MATLAB会根据新图形的大小，重新改变坐标系的比例。
1.6.5交互式图形命令
（1）从图上获取数据——ginput命令
[x,y]=ginput(n) %用鼠标从图形上获取n个点的坐标(x,y)
参数n应为正整数，是通过鼠标从图上获得数据点的个数；x、y用来存放所取点的坐标。由于鼠标所取点的位置可能有偏差，因此ginput命令获取的坐标并不是很精确。
（2）把字符串放置到图形中鼠标所指定的位置上——gtext命令
gtext(‘s’) %用鼠标把字符串放置到图形上
如果参数s是单个字符串或单行字符串矩阵，那么一次鼠标操作就可把全部字符以单行形式放置在图上；如果参数s是多行字符串矩阵，那么每操作一次鼠标，只能放置一行字符串，需要通过多次鼠标操作，把一行一行字符串放在图形的不同位置。
1.6.6 图形的标注
（1）使用grid命令显示分格线
grid on %显示分格线
grid off %不显示分格线（默认）
grid　 %在以上两个命令间切换
（2）文字标注
①添加图题
title(s) %书写图题
s为图题，为字符串，可以是英文或中文。
②添加坐标轴名
xlabel(s) %横坐标轴名
ylabel(s) %纵坐标轴名

2. 控制系统的数学模型
2.1 Laplace变换和反变换
2.1.1 直接调用MATLAB所提供的函数
F=laplace(f,t,s) %求时域函数f的Laplace变换F
f=ilaplace(F,s,t) %求F的Laplace反变换f
其中，返回结果F为s的函数，当参数s省略时，返回结果F默认为‘s’的函数；f为t的函数，当参数t省略时，默认自由变量为‘t’。
例2-1求的Laplace变换。
程序：
clear % 清除工作空间中的变量
clc % 清除命令窗口
syms t % 令t为符号变量
f=-exp(-4*t)+2*exp(-3*t)-exp(-2*t)-exp(-t)+1; % 这里的1相对于单位阶跃函数
% heaviside（1）
F=laplace(f) % 求f的Laplace变换

输出：
F =
2/(s + 3) - 1/(s + 2) - 1/(s + 1) - 1/(s + 4) + 1/s
为了便于阅读，可采用MATLAB提供的pretty()函数。
例2-1续求上例中F的Laplace反变换。
程序：
f2=ilaplace(F) % 求F的Laplace逆变换

结果：
f2 =
2/exp(3*t) - 1/exp(2*t) - 1/exp(t) - 1/exp(4*t) + 1
2.1.2部分分式展开
[r,p,k]=residue(num,den)
其中，num和den分别是分子多项式和分母多项式系数行向量：r是留数行向量；p为极点行向量；k为直项行向量。
例2-2将表达式进行部分分式展开。
程序：
clear
clc
num=[1 7 24 24]; % 分子多项式系数，降幂排列
den=[1 10 35 50 24 0]; % 分母多项式系数，注意：常数项为零（即s零
% 次方的系数为零）
[r,p,k]=residue(num,den)

结果：
r =
-1.0000
2.0000
-1.0000
-1.0000
1.0000
p =
-4.0000
-3.0000
-2.0000
-1.0000
0
k =
[]
由此可知，经部分分式展开后为

其Laplace反变换为

若留数r不是整数，可利用rat函数将其转换为分数形式：
[rnum,rden]=rat(r)
例如，[rnum,rden]=rat(0.8)的运行结果为rnum=4，rden=5，即将0.8化为。
2.2传递函数描述法
MATLAB中使用tf命令建立传递函数。
G=tf(num,den)
其中，num和den分别为传递函数分子多项式和分母多项式的系数行向量，按降幂排列。
在运算过程中，可能会用到多项式相乘和多项式求根。
p=conv(p1,p2) % p是多项式p1和p2的乘积多项式
r=roots(p) % r是多项式p的根，r以列向量形式保存
与roots相反，根据多项式的根来计算多项式的系数可以用poly函数来实现，即
p=poly (r)
例2-3已知传递函数的分子多项式系数行向量为num=[3 2 8]，分母多项式系数行向量den=[1 3 8 4 2]，求传递函数。
程序：
num=[3 2 8];
den=[1 3 8 4 2];
G=tf(num,den)

结果：
Transfer function:
3 s^2 + 2 s + 8
------------------------------------
s^4 + 3 s^3 + 8 s^2 + 4 s + 2
2.3零极点描述法
G=zpk(z,p,k) %由零点、极点和增益获得传递函数
其中，z为传递函数零点列向量，p为极点列向量，k为增益。

例2-4。
程序：
clear
clc
z=[-1;-2]; % 传递函数零点列向量
p=[0;-5;-10]; % 传递函数极点列向量
k=10; % 增益
G=zpk(z,p,k)

结果：
Zero/pole/gain:
10 (s+1) (s+2)
--------------------
s (s+5) (s+10)
2.4模型之间的转换
（1）零极点描述转换为传递函数
[num,den]=zp2tf(z,p,k)
例2-5将例2-4中零极点描述转换为传递函数
程序：
clear
clc
z=[-1;-2];
p=[0;-5;-10];
k=10;
[num,den]=zp2tf(z,p,k)

结果：
num =
0 10 30 20
den =
1 15 50 0
（2）传递函数转换为零极点描述
[z,p,k]=tf2zp(num,den)

例2-6将例2-5中传递函数转换为零极点描述
程序：
clear
clc
num=[10 30 20];
den=[1 15 50 0];
[z,p,k]=tf2zp(num,den)

结果：
z =
-2
-1
p =
0
-10
-5
k =
10
2.5结构图的化简
2.5.1 串联结构
串联结构如图2.1所示。

图2.1串联结构
实现串联结构传递函数的命令为：
G=G1*G2，或G=series(G1,G2)

2.5.2 并联结构
并联结构如图2.2所示。

图2.2并联结构
实现并联结构传递函数的命令为：
G=G1+G2或G=parallel(G1,G2)。
2.5.3 反馈结构
反馈结构如图2.3所示。

图2.3反馈结构
实现反馈结构传递函数的命令为：
G=feedback(G1,G2,sign)
其中，sign用来表示正反馈或负反馈，Sign=-1或省略则表示为负反馈。
例2-7根据系统的结构框图求出整个系统的传递函数，结构框图如图2.4所示，其中，，，。

图2.4系统结构框图
程序：
clear
clc
G1=tf(1,[1 2 1]);
G2=tf(1,[1 1]);
G3=tf(1,[2 1]);
G4=tf(1,[1 0]);
G12=G1+G2; %并联结构
G34=G3-G4; %并联结构
G=feedback(G12,G34,-1) %反馈结构

结果：
Transfer function:
2 s^4 + 7 s^3 + 7 s^2 + 2 s
--------------------------------------------
2 s^5 + 7 s^4 + 8 s^3 + s^2 - 4 s – 2

3. 控制系统的时域分析
3.1系统稳定性分析
根据Routh-Hurwitz稳定性判据，我们可以用roots()函数、pzmap()函数来判断系统的稳定性。
p=roots(den) % 求多项式的跟，den为分母多项式系数向量
pzmap(G) % 绘制传递函数G的零极点分布图
[p,z]=pzmap(G) % 求G的极点和零点
例3-1分析如图3.1所示闭环系统的稳定性。

图3.1闭环系统
程序：
clear
clc
G=tf([1],[1 1 2 23]); % 求前向通道传递函数
Gc=feedback(G,1,-1) % 求系统闭环传递函数
pzmap(Gc) % 绘制传递函数G的零极点分布图
[p,z]=pzmap(Gc) % 求G的极点和零点

结果：
Transfer function:
1
-------------------------
s^3 + s^2 + 2 s + 24
p =
-3.0000
1.0000 + 2.6458i
1.0000 - 2.6458i
z =
Empty matrix: 0-by-1

图3.2系统零极点分布图
可以看出，该系统中存在不稳定的极点。
3.2 系统动态特性分析
阶跃响应可以用step命令来实现。
step(G, tf) %绘制从0到tf时间段内的系统阶跃响应曲线
[y,t]=step(G,tf) %求阶跃响应在（t,y）时间点上的值，不绘制曲线
y=step(G,t) %求阶跃响应在制定时间点t上的值，不绘制曲线
其中，G为系统模型，可以是传递函数、状态方程、零极点增益的形式；y为时间响应；t为时间向量；tf为时间区间终点，可省略。
稳态值可使用dcgain命令来得出。
k=dcgain(G) %获得稳值

例3-2系统的传递函数为，绘制系统的阶跃响应曲线。
程序：
clear;
clc;
num=1;
den=conv([1 0.1 5],[1 2 3 4]);
step(num,den,40); % 绘制系统从0到40s的阶跃响应曲线
hold on
t1=0:1:40;
y1=step(num,den,t1); %求对应时间点t1的响应值
plot(t1,y1,'r');
grid on;

结果：

从图中可以看出，仿真时间的步长不宜取得太大，否则将影响曲线的光滑程度。
例3-3某典型二阶系统的传递函数为，求和时，系统阶跃响应。
程序：
close all
clear
clc
wn=1;
num=[wn^2];
zeta=[0 0.1 0.5 0.7 1 1.5 2];
fori=1:7
den=[1 2*zeta(i)*wn wn^2];
G=tf(num,den);
step(G,10)
hold on
end
grid on
title('step responses with different zeta')
figure
zeta=0.7;
wn=[1 2 4];
fori=1:3
num=wn(i)^2;
den=[1 2*zeta*wn(i) wn(i)^2];
G=tf(num,den);
step(G,5)
hold on
end
grid on
title('step responses with different wn')

结果：

例3-4已知系统传递函数为，计算系统瞬态响应指标（稳态误差、峰值、超调量、上升时间、调节时间）。
程序：
clc;
clear
% 系统模型建立
num=[3];
den=conv([1 1+3j],[1 1-3j]);
% 求系统的单位阶跃响应
step(num,den);
grid on;
[y,x,t]=step(num,den);
% 求响应的稳态值
finalvalue=dcgain(num,den)
% 求响应的峰值及对应的下标
[yss,n]=max(y);
% 计算超调量及峰值时间
overshoot=100*(yss-finalvalue)/finalvalue
time_to_peak=t(n)
% 计算上升时间
n=1;
while y(n)<0.1*finalvalue
n=n+1;
end
m=1;
while y(m)<0.9*finalvalue
m=m+1;
end
rise_time=t(m)-t(n)
% 计算调节时间
er=input('Pleasr input expect error range ='); % 输入误差带范围
k=length(t);
while (y(k)>(1-er)*finalvalue)&(y(k)<(1+er)*finalvalue)
k=k-1;
end
setting_time=t(k)

结果：
finalvalue =
0.3000
overshoot =
35.0913
time_to_peak =
1.0492
rise_time =
0.4197
Pleasr input expect error range =0.02
setting_time =
3.5147

4. 控制系统的频域分析
4.1 频域特性
频率特性函数与传递函数之间的关系为

频率特性可由下式求出
Gw=polyval(num,j*w)/polyval(den,j*w)
mag=abs(Gw) %幅频特性
pha=angle(Gw) %相频特性
其中，函数polyval(p,s)用来计算多项式在给定变量时的值，即将s带入多项式p。
例4-1求二阶系统传递函数在时的频域特性。
程序：
clear
clc
num=1;
den=[1 1.414 1];
w=1 ;
Gw=polyval(num,j*w)/polyval(den,j*w)
mag=abs(Gw)
pha=angle(Gw)

结果：
Gw =
0 - 0.7072i
mag =
0.7072
pha=
-1.5708
4.2 绘制Bode图
bode(G,w) %绘制bode图，w可省略
[mag,pha]=bode(G,w) %得出w对应的幅值和相角
[mag,pha,w]=bode(G) %得出幅值、相角和频率
例4-2绘制系统的Bode图。
程序：
clear
clc
num=1;
den=conv([1 1],[1 ,2]);
bode(num,[den 0])
grid on

结果：

例4-2续使用semilogx命令绘制对数幅频和相频特性。
w=logspace(a,b,n) % 用对数尺度（a，b）间频区间上取n个点
semilogx与plot命令类似，但在画图时对横轴所表示的量取以10为底的对数。

程序：
figure
w=logspace(-2,2,100);
[mag,pha]=bode(num,[den 0],w);
subplot(2,1,1) % 绘图窗口分成两行一列（即上下两张图），
% 选择第一张图
semilogx(w,20*log10(mag))
grid on
subplot(2,1,2) % 绘图窗口分成两行一列（即上下两张图），
% 选择第二张图
semilogx(w,pha)
grid on

结果：

例4-3某典型二阶系统的传递函数为，绘制和时，系统的Bode图。

程序：
close all
clear
clc
wn=1;
num=[wn^2];
zeta=[0 0.1 0.3 0.7 1 1.5 2];
fori=1:7
den=[1 2*zeta(i)*wn wn^2];
G=tf(num,den);
Bode(G)
hold on
end
grid on
title('Bode diagrams with different zeta')
figure
zeta=0.7;
wn=[1 2 4];
fori=1:3
num=wn(i)^2;
den=[1 2*zeta*wn(i) wn(i)^2];
G=tf(num,den);
Bode(G)
hold on
end
grid on
title('Bode diagrams with different wn')

结果：

从图中可以看出，当时，对数幅频特性在角频率略低于处有一个谐振峰，阻尼系数越小，谐振峰越高。相频特性在处有，在处，，在的频段有。相频特性关于，点中心对称。整条相频特性的形状因的不同而有很大的不同，越小，相频特性在附近变化就越剧烈，曲线越陡。
4.3 绘制Nyquist曲线
nyquist (G,w) %绘制nyquist曲线，w可省略
[Re,Im]= nyquist (G,w) %由w得出对应的实部和虚部
[Re,Im,w]= nyquist (G) %得出实部、虚部和频率
其中，G为系统模型；w为频率向量，也可以用{wmin,wmax}表示频率的范围；Re为频率特性的实部，Im为频率特性的虚部。
例4-4单位负反馈系统的开环传递函数为，绘制其Nyquist曲线。
程序：
close all
clear
clc
G=tf(1000,conv([1 1],conv([1 2],[1 5])));
nyquist(G)
figure
nyquist(G)
axis([-12,1,-8,8]) % 设置坐标轴范围

结果：

从图中可以看出，Nyquist顺时针包围（-1，j0）点两周，而开环传递函数在右半平面不存在极点，故系统不稳定。求系统闭环极点为
ans =
-12.8196
2.4098 + 8.5427i
2.4098 - 8.5427i
显然闭环系统右两个不稳定的极点。
4.4 幅值裕度和相位裕度
margin(G) %绘制bode图并标出幅值裕度和相角裕度
[mag,map,wcg,wcp]=margin(G) %得出幅值裕度和相角裕度
其中：mag为幅值裕度，wcg为幅值裕度对应的频率；map为相角裕度，wcp为相角裕度对应的频率（穿越频率）。如果wcg或wcp为nan为Inf，则对应的Gm或Pm为无穷大。
例4-5已知传递函数，求K分别为2和20时的幅值裕度和相位裕度。
程序：
close all
clear
clc
num1=2;
num2=20;
den=conv([1 0],conv([1 1],[0.2 1]));
figure(1)
[mag1,map1,wcg1,wcp1]=margin(num1,den)
margin(num1,den)
grid on
figure(2)
w=logspace(-1,2);
[mag2,pha2]=bode(num2,den,w);
margin(mag2,pha2,w);
grid on

结果：
mag1 =
3.0000
map1 =
25.3898
wcg1 =
2.2361
wcp1 =
1.2271

将用margin函数求出的增益裕量以dB表示，20*log10(mag1)= 9.5424，即与图中所示相一致。

5. 控制系统的根轨迹分析
5. 1绘制根轨迹
rlocus(G,k) %根据给定的k变化范围绘制根轨迹，k可省略
[r,k]=rlocus(G) %得出闭环极点和对应的K，不画出根轨迹
r= rlocus(G,k) %根据K得出对应的闭环极点，不画出根轨迹
[k,p]=rlocfind(G) % 找出给定的根（闭环极点）所对应的根轨迹增益
该函数要求在屏幕上预先画好根轨迹图，然后，此命令将产生一个光标，供用户选择希望的闭环极点。命令的执行结果为：给出选择点所对应的根轨迹开环增益，并给出在该根轨迹增益下的闭环特征根。
例5-1已知系统，绘制其根轨迹。
程序：
close all
clear
clc
num=[2 4];
den=[8 3 1 0];
G=tf(num,den);
rlocus(num,den);
grid on

结果：

例5-2已知开环系统的数学模型为，绘制其根轨迹。
程序：
close all
clear
clc
num=[1];
den=conv([0.01 1 0],[0.02 1]);
rlocus(num,den);
[k1,p1]=rlocfind(num,den)
[k2,p2]=rlocfind(num,den)

结果：
Select a point in the graphics window
selected_point =
-14.6327 +34.9379i
k1 =
34.7393
p1 =
1.0e+002 *
-1.2046
-0.1477 + 0.3498i
-0.1477 - 0.3498i

Select a point in the graphics window
selected_point =

-6.1019 +55.1242i
k2 =
85.1077
p2 =
1.0e+002