资源描述
《控制工程基础(经典控制部分)》的MATLAB分析
第一章 MATLAB的基本使用
§1-1 MATLAB语言简介
MATLAB是一种高级矩阵语言,它由Math Works公司于1984年正式推出,它的基本处理对象是矩阵,即使是一个标量纯数,MATLAB也认为它是只有一个元素的矩阵。随着MATLAB的发展,特别是它所包含的大量工具箱(应用程序集)的集结,使MATLAB已经成为带有独特数据结构、输入输出、流程控制语句和函数、并且面向对象的高级语言。
MATLAB语言被称为一种“演算纸式的科学计算语言”,它在数值计算、符号运算、数据处理、自动控制、信号处理、神经网络、优化计算、模糊逻辑、系统辨识、小波分析、图象处理、统计分析、甚至于金融财会等广大领域有着十分广泛的用途。
MATLAB语言在工程计算与分析方面具有无可比拟的优异性能。它集计算、数据可视化和程序设计于一体,并能将问题和解决方案以使用者所熟悉的数学符号或图形表示出来。
MATLAB语言和C语言的关系与C语言和汇编语言的关系类似。例如当我们需要求一个矩阵的特征值时,在MATLAB下只需由几个字符组成的一条指令即可得出结果,而不必去考虑用什么算法以及如何实现这些算法等低级问题,也不必深入了解相应算法的具体内容。就象在C语言下不必象汇编语言中去探究乘法是怎样实现的,而只需要采用乘积的结果就可以了。
MATLAB语言还有一个巨大的优点是其高度的可靠性。例如对于一个病态矩阵的处理,MATLAB不会得出错误的结果,而用C或其它高级语言编写出来的程序可能会得出错误的结果。这是因为MATLAB函数集及其工具箱都是由一些在该领域卓有研究成果,造诣很深的权威学者经过反复比较所得出来的最优方法,而且经过多年的实践检验被证明是正确可靠的。
§1-2 MATLAB的工作窗口
下面以MATLAB6.1为例介绍。
从实用的角度MATLAB的工作窗口包括命令窗口、M文件编辑器窗口、图形编辑窗口、数学函数庫、应用程序接口及在线窗口。下面首先介绍MATLAB的命令窗口及M文件编辑器。
一、命令窗口
启动MATLAB之后,屏幕上自动出现命令窗口MATLAB,它是MATLAB提供给用户的操作界面,用户可以在命令窗口内提示符“>>”之后(有的MATLAB版本命令窗口没有提示符)键入MATLAB命令,回车即获得该命令的答案。
命令窗口内有File、Edit、View、Web、Window、Help等菜单条。
二、M文件编辑窗口
M文件是MATLAB语言所特有的文件。用户可以在M文件编辑窗口内,编写一段程序,调试,运行并存盘,所保存的用户程序即是用户自己的M文件。MATLAB工具箱中大量的应用程序也是以M文件的形式出现的,这些M文件可以打开来阅读,甚至修改,但应注意,不可改动工具箱中的M文件!
1.进入M文件窗口有两种方法
1) 命令窗口 — File — New — M-File;
2) 命令窗口 — 点击“File”字样下面的图标。
M文件编辑窗口的标记是“Untitled”(无标题的)。当用户编写的程序要存盘时,Untitled作为默认文件名提供给用户,自然,用户可以,也应当自己命名。若用户不自己命名,则MATLAB会对Untitled进行编号。
2.M文件的执行:
返回命令窗口,在当前目录(Current Directory)内选择所要运行的M文件的目录,在命令窗口提示符“>>”后,直接键入文件名(不加后缀)即可运行。
注意:(1)机器默认路径为一级子目录MATLAB6 p1\work;
(2)MATLAB 6.1以前的版本,运行M文件的方法稍有不同,它必须在File菜单下,打开“Run Script…”子菜单,键入需要运行的文件路径及名称再回车,在这种情况下,work作为根目录对待,不出现在M文件的路径之中。本讲义的参考程序都是在M文件窗口下编制的。
三、在线帮助窗口
在命令窗口中键入Help(空格) 函数名,可以立即获得该函数的使用方法。
§1-3 MATLAB最基本的矩阵操作
作为命令窗口及M文件编辑器的应用实例,介绍几个最基本的矩阵运算命令。
一、矩阵的输入
在方括号内依次按行键入矩阵元素,在一行内的各元素之间用空格或逗号分开,每行之间用分号分开。
例如,在命令窗内输入
A=[2 2 3;4 5 4;7 8 9] (注意:方括号,分号为矩阵行标记)
B=[1,3,5;6,-4,2;3,5,1] (逗号与空格功能相同)
A= 2 2 3 B= 1 3 5
4 5 4 6 -4 2
7 8 9 3 5 1
同理:输入A1=[2 4 6]得到行矢量,
输入A2=[2;4;6]得到列矢量,
于是,当输入
C=[A;A1] 有
C= 1 2 3
4 5 6
7 8 9
2 4 6
A1作为矩阵C的最后一行,C和A相比,增加了一行。
二、矩阵的转置
矩阵A的转置用A′表示,显然,A1与A2互为转置,即A1'会得到以2,4,6为元素的列矢量。
思考一下输入 C1=[ A A2 ]
C2=[ A A1']
有什么结果?而输入[A;A1']有无意义?
三、矩阵的四则运算
1.矩阵的加减法:
当两个矩阵维数相同时可以直接进行“+”或“-”运算。
如 D1=A+B,D2=A-B
2. 矩阵的乘法:
当矩阵A,B维数相容时
C3=A﹡B:普通意义下的矩阵相乘
C4=A .﹡B:矩阵A与B的对应元素相乘
显然,A﹡B≠B﹡A(一般情况),而A .﹡B=B .﹡A。A .﹡B称为数列型乘法,它要求参加运算的矩阵或数列具有相同的行列数,这是MATLAB语言中的一种特殊运算,它在今后求取函数值等运算时是很重要的。实际上,前面所述的矩阵加、减法就是一种数列型运算。
3. 矩阵的除法
D4=A\B:表示A-1﹡B或inv(A)*B,即A的逆矩阵左乘矩阵B。
D5=B/A:表示B﹡A-1或B﹡inv(A),即A的逆矩阵右乘B。
D6=A .\B:表示B的每一个元素被A的对应元素除。
D7=A ./B:表示A的每一个元素被B的对应元素除。显然,A .\B与A ./B的各对应元素互为倒数。
读者可以思考一下,D6 .﹡D7等于什么?
D8=inv(A):A的逆矩阵。
打开M文件编辑窗口,将上述命令依次键入,得到fanli001如下:
参考程序fanli001:矩阵的四则运算
A=[2 2 3;4 5 4;7 8 9] % 三阶矩阵输入
B=[1,3,5;6,-4,2;3,5,1] % 三阶矩阵输入
A1=[2 4 6] %行向量
A2=[2;4;6] %列向量
C=[A;A1] %矩阵A增加一行
C1=[A A2] %矩阵A增加一列
C2=[A A1'] %矩阵A增加一列
D1=[A+B] %矩阵相加
D2=A-B %矩阵相减
C3=A*B %矩阵与矩阵相乘
C4=A.*B %矩阵的对应元素相乘
D3=A\B %A的逆左乘B
D4=B/A %A的逆右乘B
D6=A.\B %B的各元素被A的对应元素除
D7=A./B %A的各元素被B的对应元素除
D8=inv(A) %A的逆矩阵
语句后面的%为语句说明符。
MATLAB中矩阵运算的其它主要命令可通过在线帮助获得。
§1-4 MATLAB的符号运算操作
一、进入符号运算功能
在命令窗口键入
syms x y z t
此后,即可以使用x,y,z,t等作自变量定义函数。
syms x y z t real 规定所定义的变量为实型。
二、代数方程求解
使用命令solve 可以求解代数方程,如求下例方程组
的解,命令为
[x,y]=solve(‘x^2+x*y+y-3=0,x^2-4*x-2*y+3=0’)
程序见范例程序fanli002。
参考程序fanli002:代数方程求解
syms x y %进入符号运算功能
f1=x^2+x*y+y-3 %函数f1
f2=x^2-4*x-2*y+3 %函数f2
[x,y]=solve('x^2+x*y+y-3=0','x^2-4*x-2*y+3=0')
%求解方程组
f1a=simplify(subs(f1)) %用求解出的x,y检验方程1
f2a=simplify(subs(f2)) %用求解出的x,y检验方程2
x=double(x) %将符号变量转换成浮点数
y=double(y) %将符号变量转换成浮点数
f1=subs(f1) %用浮点数x,y检验方程1
f2=subs(f2) %用浮点数x,y检验方程2
由solve求出的根是根的符号表达形式,是准确解。
命令simplify(f)表示化简,subs(f1)表示将求解出x,y代回f1中;double是将符号变量转换成浮点数,是准确解x,y的近似值。这从程序运行f1a=0,f2a=0,f1≠0,f2≠0可以确认。
三、符号矩阵运算
符号矩阵可以和数值矩阵一样进行运算,例如:求矩阵特征值eig,求矩阵的逆inv等命令都支持符号运算。
设
计算其特征值eig(A),逆矩阵inv(A),程序见fanli003。
参考程序fanli003:符号矩阵的特征值
syms t real %定义为实型变量
A=[sin(t) -cos(t);cos(t) sin(t)] %定义矩阵A
B1=eig(A) %求矩阵A的特征值
B1=simple(B1) %化简A的特征值表达式
B2=inv(A) %求矩阵A的逆矩阵
B2=simple(B2) %化简A的逆矩阵表达式
C1=A*B2 %检验A的逆矩阵
C1=simple(C1) %C1为单位矩阵
注意函数的输入方法,自变量用圆括号括起来。
四、微积分运算
设函数,则MATLAB中微积分运算命令为
fiff(f):求函数f对自变量x的一阶导数;
diff(f,2):求函数f对自变量x的二阶导数;
int(f):求函数f的不定积分
例1.1,设
试计算其一阶,二阶导数,积分运算,并作出函数图象,见范例fanli004。
参考程序fanli004:函数的微分与积分
syms x
f=1/(5+4*cos(x))
ezplot(f) %函数f的曲线
f1=diff(f) %函数f的一阶导数
figure,ezplot(f1) %函数f一阶导数的曲线
f2=diff(f,2) %函数f的二阶导数
figure,ezplot(f2) %函数f二阶导数的曲线
g=int(int(f2)) %函数f的二阶导数f2的二重积分
figure,ezplot(g) %函数f2二重积分的曲线
e=f-g %二阶导数的二重积分与原函数的差
e=simple(e)
figure,ezplot(e)
程序中
ezplot(f):作函数的图形,x的取值范围默认值为-2π< x <2π
fanli004的函数曲线
fanli004的一阶导函数曲线
fanli004的二阶导函数曲线
fanli004的二阶导函数的二重积分曲线
ezplot(f)是一个很有用的作图命令,它的其它应用形式,请查在线帮助。
细心的读者会发现,一个函数求二阶导数后再对二阶导数进行二重积分,其结果与原函数相差一个常数。相当于纵坐标发生平移。
第二章 系统的时域特性
§2-1 传递函数
一、传递函数的两种形式
传递函数通常表达成s的有理分式形式及零极点增益形式。
设传递函数
1.有理分式形式
分别将分子、分母中s多项式的系数按降幂排列成行矢量,缺项的系数用0补齐。上述函数可表示为
num1=[2 1] %(注意:方括号,同一行的各元素间留空格或逗号)。
den1=[1 2 2 1]
syss1=tf(num1,den1)
运行后,返回传递函数的形式。这种形式不能直接进行符号运算!
2.零极点增益形式
[Z,P,K] = tf2zp(num1,den1)
sys2 = zpk(Z,P,K)
返回零、极点、增益表达式,其Z,P分别将零点和极点表示成列向量,若无零点或极点用[ ](空矩阵)代替。
运行得到的
点 Z = -0.5
极点 P= -1,-0.5±j0.866
增益 K = 2
指令zp2tf(Z,P,K)将零极点增益变换成有理分式形式,见程序fanli005。
参考程序fanli005:传递函数的有理分式及零极点增益模型
num1=[2 1] % 传递函数的分子系数向量
den1=[1 2 2 1] % 传递函数的分母系数向量
sys1=tf(num1,den1) % 传递函数的有理分式模型
[Z,P,K]=tf2zp(num1,den1)
% 有理分式模型转换成零极点增益模型
[num2,den2]=zp2tf(Z,P,K)
% 零极点增益模型转换成有理分式模型
sys2=zpk(Z,P,K) % 传递函数的零极点增益模型
[A1,B1,C1,D1]=tf2ss(num1,den1)
% 有理分式模型转换成状态空间模型
[A2,B2,C2,D2]=zp2ss(Z,P,K)
% 零极点及增益模型转换成状态空间模型
[num1,den1]=ss2tf(A1,B1,C1,D1)
% 状态空间模型转换成有理分式模型
[Z,P,K]=ss2zp(A2,B2,C2,D2)
% 状态空间模型转换成零极点增益模型
程序中,命令tf2ss,zp2ss及ss2tf,ss2zp是状态空间模型与有理分式及零、极点、增益模型之间的相互转换。
二、传递函数框图的处理
用框图可以方便地表示传递函数的并联,串联及反馈。为简洁,仅以有理分式模型为例。
G1
G1
G1+G2
1. 并联
sysp = parallel(sys1,sys2)
[num,den]=parallel(num1,den1,num2,den2)
2. 串联
G1(s)
G2(s)
G1(s)·G2(s)
syss = series(sys1,sys2)
[nums,dens] = series(num1,den1,num2,den2)
3. 反馈
G1(s)
G2(s)
G1(s)G2(s)
1+G1(s)G2(s) G3(s)
G3(s)
sysc=feedback(syss,sys3,±1) %默认值(-1)
[numc, denc] = feedback(nums, dens, num3, den3)
4. 单位反馈
G1(s)
G2(s)
G1(s)G2(s)
1+G1(s)G2(s)
sysd = feedback(syss, 1)
[numd, dend] = feedback(nums, dens, 1, 1) %(单位反馈)
上面给出了同一指令的两种形式,相当于两套平行指令。
对于零极点增益形式,书写稍复杂一些,可先用zpk转换成系统形式,或用zp2tf转折换成有理分式形式后再进行框图化简操作。
三、简单函数的拉普拉斯变换
在MATLAB的符号功能中,可以对简单函数进行拉普拉斯正、逆变换。
拉氏正变换:laplace(f(t))
拉氏逆变换:ilaplace(L(s))
其中为原函数,为象函数。命令格式参见fanli007。
参考程序fanli007:拉普拉斯变换
syms s t w a b c
f1=sqrt((b-a)^2+w^2)/w*exp(-a*t)*sin(w*t+atan(w/(b-a)))%原函数f1
L1=laplace(f1) %f1的拉氏变换(象函数)
L1=simple(L1) %化简
f2=ilaplace(L1) %L1的拉氏逆变换
f2=simple(f2) %化简
在MATLAB中使用laplace及ilaplace命令时,要注意象、原函数的符号,特别是对初相不等于零的振荡系统,运行结果常常同手册上的结果相差一个符号,这要注意函数表达式成立的条件。保险的办法是再使用拉氏变换的初值定理确定象、原函数的符号。
§2-2 系统时域特性曲线
在MATLAB中,当传递函数已知时,可以方便地求出系统的单位脉冲响应、单位阶跃响应等曲线。
一、系统的单位阶跃响应step
step有以下几种格式
step(sys):
直接作出sys的单位阶跃响应曲线。其中sys = tf(num, den) 或sys = zpk(z, p, k),MATLAB自动决定响应时间。
step(sys, t)
设定响应时间的单位阶跃响应。t可以设定为最大响应时间 t =
t终值(秒),也可以设置为一个向量
t = 0 : △t : t终值
注意冒号的使用。它产生一个从0到t终值的行矢量,元素之间的间隔为△t。
step(sys1, sys2, …, sysn)
在同一幅图上画出几个系统的单位阶跃响应。
[y, t] = step(sys);
命令输出对应时刻t的各个单位阶跃响应值,不画图。语句后的分号控制数据的屏幕显示。
如果要查看机器计算了多少个数据,可以使用命令
size(y)
得出的结果也表明数据作为列矢量的行数。要将计算出的[y, t]作成曲线,使用一般的作图命令
plot(t, y)
plot后面跟的两个参数横坐标在前,纵坐标在后。
参考程序见fanli008:
参考程序fanli008:系统的单位阶跃响应
num1=[4 2]
den1=[2 8 14 11 4]
sys1=tf(num1,den1) %系统G1(s)
num2=[2 1]
den2=[1 4 6 7 3]
sys2=tf(num2,den2) %系统G2(s)
[y1,t1]=step(sys1); %系统G1(s)的单位阶跃响应数据
[y2,t2]=step(sys2); %系统G2(s)的单位阶跃响应数据
step(sys1,sys2) %系统G1(s)、G2(s)的单位阶跃响曲线
figure,step(sys1,sys2,20)
%系统G1(s)、G2(s)在自选时间(20秒)内的单位阶跃响曲线
figure,plot(t1,y1) %系统G1(s)的单位阶跃响应曲线
figure,plot(t2,y2) %系统G2(s)的单位阶跃响应曲线
fanli008:step(sys1,sys2,t)单位阶跃曲线
二、系统的单位脉冲响应
impulse命令格式与单位阶跃响应step的命令格式完全相同,只需将语句中的step用impulse代替即可。针对同样的系统,其单位脉冲响应的参考程序见fanli009。
参考程序fanli009:系统的单位脉冲响应
num1=[4 2]
den1=[2 8 14 11 4]
sys1=tf(num1,den1)
num2=[2 1]
den2=[1 4 6 7 3]
sys2=tf(num2,den2)
[y1,t1]=impulse(sys1); %系统G1(s)的单位脉冲响应数据
[y2,t2]=impulse(sys2); %系统G2(s)的单位脉冲响应数据
impulse(sys1,sys2) %系统G1(s)、G2(s)的单位脉冲响应曲线
figure,impulse(sys1,sys2,20)
%系统G1(s)、G2(s)在自选时间(20秒)内的单位脉冲响应曲线
figure,plot(t1,y1) %系统G1(s)的单位脉冲响应曲线
figure,plot(t2,y2) %系统G2(s)的单位脉冲响应曲线
hold on,step(sys2) %系统G2(s)的单位阶跃和单位脉冲响应曲线
fanli009:impulse(sys1,sys2)单位脉冲响应曲线
程序的最后一句hold on是当前图形保护模式。当要将新图形作在当前图形上时,必须使用hold on。而figure的含意是另开一个新的图形窗口,如果不用figure或hold on,则新的图形会占用原图形窗口,始终只保留一个最新的图形窗口。
三、一阶系统及二阶系统的时域特性
一阶系统及二阶系统是最基本也是最重要的系统,高阶系统总可以视为由若干个一阶和(或)二阶系统组合构成。
1. 一阶系统(设增益为1)
影响系统特性的参数是其时间常数T,T越大,系统惯性越大,响应越慢。
参考程序fanli010给出了T=0.4, 1.2, 2.0, 2.8, 3.6, 4,4六条单位阶跃响应曲线。
参考程序fanli010:一阶系统的单位阶跃响应曲线
num=1;i=1;
for del=0.1:0.2:1.1 %一阶系统时间常数递增间隔
den=[4*del 1]; %一阶系统分母向量
step(tf(num,den)) %一阶系统单位阶跃响应曲线
hold on, %不同时间常数的一阶系统单位阶跃响应曲线簇
i=i+1;
end
同理,可以作出对应的单位脉冲响应曲线,参考程序fanli011。
参考程序fanli011:一阶系统的单位脉冲响应曲线
num=1;i=1;
fanli010一阶系统时间常数对单位阶跃响应的影响
fanli011一阶系统时间常数对单位脉冲响应的影响
for del=0.1:0.2:1.1
den=[4*del 1];
impulse(tf(num,den),10) %一阶系统单位阶脉冲应曲线
hold on, %不同时间常数的一阶系统单位脉冲响应曲线簇
i=i+1;
end
注意MATLAB中for语句的结构。读者可以改变不同的增益,看看图形有何变化。
2. 二阶系统(设0<ξ<1)
设二阶系统为
二阶系统的特征参数为固有频率及阻尼比ξ。当增大,系统振动频率加快,振荡加剧;而随着ξ减小,系统振荡加剧,振荡峰尖锐。
参考程序fanli012示出了当ξ=0.5, =1, 2, 3, 4, 5 rad/s时的间接阶跃曲线簇。
参考程序fanli012:
不同固有频率的二阶系统的单位阶跃响应曲线(ξ=0.5)
i=1;
for del=1:1:5; % 二阶系统固有频率递增间隔
num=del^2; % 二阶系统传递函数分子系数向量
den=[1 del del^2]; % 不同固有频率的二阶系统分母系数向量
step(tf(num,den),6) %二阶系统单位阶跃响应曲线
hold on, %不同固有频率的二阶系统单位阶跃响应曲线簇
i=i+1;
end
fanli012二阶系统固有频率对单位阶跃响应的影响
参考程序fanli013示出了同一二阶系统的单位脉冲响应曲线簇。
参考程序fanli013:
不同固有频率的二阶系统的单位脉冲响应曲线(ξ=0.5)
i=1;
for del=1:1:5;
num=del^2;
den=[1 del del^2]; % 不同固有频率的二阶系统分母系数向量
impulse(tf(num,den),6) %二阶系统单位脉冲响应曲线
hold on, %不同固有频率的二阶系统单位脉冲响应曲线簇
i=i+1;
end
fanli013二阶系统固有频率对单位脉冲响应的影响
参考程序fanlio14示出了当=1,ξ=0.1,0.3,0.5,0.7,0.9的二阶系统的单位阶跃响应曲线簇。
参考程序fanli014:
不同阻尼比的二阶系统的单位阶跃响应曲线(=1)
i=1;
for del=0.1:0.2:0.9; % 二阶系统阻尼比ξ递增间隔
num=1;
den=[1 2*del 1]; % 不同阻尼比的二阶系统分母系数向量
step(tf(num,den),30)
hold on, %不同阻尼比的二阶系统单位阶跃响应曲线簇
i=i+1;
end
famli014二阶系统阻尼比对单位阶跃响应的影响
参考程序fanli015示出了同一二阶系统当ξ取0.1,0.3,0.5,0.7,0.9的单位脉冲响应曲线簇。
参考程序fanli015:
不同阻尼比的二阶系统的单位脉冲响应曲线(=1)
i=1;
for del=0.1:0.2:0.9;
num=1;
den=[1 2*del 1];
impulse(tf(num,den),30)
hold on, %不同阻尼比的二阶系统单位脉冲响应曲线簇
i=i+1;
end
fanli015二阶系统阻尼比对单位脉冲响应的影响
四、带延时环节系统的典型响应
设具有纯延时环节的传递函数为
计算这种系统的单位阶跃响应不能使用一般方法。首先应使用Pade法对延时环节进行近似展开,
[numT, dent] = pade(T, 5) (表示使用5阶Pade级数)
得出的分子、分母系数的向量,然后将系统视为惯性环节与延迟环节的串联,再求其阶跃脉冲响应。
设K = 2.5,T=1,程序见fanli016。
参考程序fanli016: 带延时环节系统的单位阶跃响应
K=2.5;T=1; T1=0.5;
[numt,dent]=pade(T,5);%5阶Pade法近似延时环节
syst=tf(numt,dent); %延时环节的近似传递函数
num1=K;
den1=[T1 1];
sys1=tf(num1,den1);%惯性环节
sys=series(sys1,syst);% 带延时环节的惯性系统
step(sys,6),grid % 带延时环节惯性系统的单位阶跃响应
figure,impulse(sys,6),grid
% 带延时环节惯性系统的单位脉冲响应
fanli016带延迟的一阶系统的单位阶跃响应
fanli016带延迟的一阶系统的单位脉冲响应
从响应曲线上可见,具有延时环节的惯性系统当加上单位阶跃或单位脉冲输入前后,其响应初始有振荡特性,脉冲响应更为明显。这种振荡现象使得对这种系统进行校正需要使用特殊的方法。
在作图命令后,加上格线命令grid,机器自动给曲线图加上格线。
§2-3 响应曲线的动态分析
一、MATLAB图形的编辑
在MATLAB中,用于编辑图形的命令很多,因为可视化正是MATLAB的一种强大而优越的性能。从工程实用角度,仅在图形窗口进行编辑已经夠用了。
在图形窗口顶部的第二排依次排有下列命令按钮:
Edit Plot():图形编辑。点击进入图形编辑功能。点击选择编辑对象。例如,曲线,标题,纵、横坐标的说明等,可对选中项进行编辑。
Insert Text():插入文本。点击可在图区内插入文字,并可选择字体及大小等。
Insert Arrow():插入箭头。点击可在图区内画箭头,并可对该箭头进行编辑。
Insert Line():插入直线。点击可在图区内插入直线,并可对该直线进行编辑,改变线宽及颜色。
Zoom In():曲线放大。激活该按钮,可以逐次放大曲线,以观察曲线某些部分的细节,这对观察变化剧烈的曲线部分很有帮助。
Zoom Out():曲线缩小。激活该按钮,可以逐次从放大状态返回到原状态。
Rotate 3D():三维旋转。激活该按钮,在图区内按住鼠标左键,拖曳鼠标,原二维曲线可三维空间内任意旋转。
此外,还有平面旋转、三维动画、前后放大、平移、飞逸等。
例如,将曲线加粗,以利作图(曲线默认宽度为0.5)。按Edit Plot— 选择曲线双击左键 — 出现线宽选择框(line width)— 选择线宽— OK。
二、对曲线进行数据分析
不激活Edit Plot按钮(图形编辑不使能),
▲ 读取曲线上任一点的坐标值:鼠标指向选择点,单击左键,显示曲线编号,横坐标及纵坐标值;
▲ 研究曲线拟合情况:按Tools — Basic Fitting,出现Plot fitt(图形拟合)窗口,选择需要拟合的曲线,点击拟合方式,即可对所选图形进行曲线拟合。
MATLAB提供了样条插值、保形插值(shape-preserving interpolant)及直到10阶的多项式拟合。当选择多项式拟合时,机器会给出曲线及拟合方程。
在所给出的拟合方式中,三次样条插值(Cubic spline interpolant)及保形插值(shape-preserving interpolant)效果最好。
第三章 系统的频率特性
系统的频率特性就是,这是一个复变函数。
§3-1 MATLAB中的复数及复变量
在MATLAB中复数x的输入方法:
x = 2 + 3*i 或 2 + 3*j
也可以s = 2 + 3i 或 s = 3 + 2j
在符号运算功能下,复变函数表示为(设a,b都是实数):
syms a b real
x = a–b*I 或a–b*j(*不能省略)
在MATLAB中与复数有关的命令有
abs(x):求实数的绝对值及复数的模,支持符号运算;
angle(x):求复数的相角,单位为弧度,不支持符号运算。
conj(x):求复数的共轭复数,支持符号运算;
imag(x):求复数的虚部,支持符号运算;
real(x):求复数的实部,支持符号运算。
复数的运算见参考程序fanli017。
参考程序fanli017:复数运算(含符号功能)
syms a b c d real %符号均为实变量
X1=2-3i
X2=a+c/b*j
X3=c/a-2*d*j
S1=X2*X3/X1 %复数的乘除
S2=simple(abs(S1)) %复数的模
S3=simple(real(S1)) %复数的实部
S4=simple(imag(S1)) %复数的虚部
S5=simple(conj(S1)) %复数的共轭
S6=simple(expand(S1*S5)) %复数与其共轭的积
S7=simple(S3^2+S4^2) %复数实部、虚部的平方和
S8=simple(S6-S7) %检验复数运算的正确性
fanli018给出了系统的频率特性。
参考程序fanli018:系统频率特性的解析表达
syms w real %定义频率w 为实型变量
g=50*(0.6*w*i+1)/(i*w)^2/(4*w*j+1) %系统频率特性
g1=simple(conj(g)) %系统频率特性的复共轭
gr=simple(real(g)) %系统实频特性
gi=simple(imag(g)) %系统虚频特性
ga=simple(abs(g)) %系统幅频特性(模)
gb=simple(sqrt(gr^2+gi^2)) %系统幅频特性的另一计算方法
gc=simple(ga-gb) %检验系统幅频特性的两种计算方法
§3-2 频率特性的Nyquist图
将系统传递函数中的复变数s用纯虚数jω(ω为角频率,rad/s)代替,即得到系统的频率特性,又称为谐波传递函数,有三种表示方法:
式中
: 实频特性
: 虚频特性
: 幅频特性
: 相频特性
在研究控制系统的特性时,采用图示方法比采用解析方法直观、简单,便于对动态特性进行更深入的研究。以下讨论设系统传递函数为sys。
一、Nyquist图
1. Nyquist 图的有关命令
nyquist(sys)
直接返回系统sys的Nyquist图,机器自己确定频率范围(通常为 -∞ → +∞)
nyquist(sys, w)
返回在规定频率(rad/s)范围内的Nyquist图,频率范围的格式为
w = {wmin, wmax} (0<wmin<wmin)
nyquist(sys1,sys2,…,sysn)
nyquist(sys1,sys2,…,sysn,w)
在同一幅度图上画出几个系统的Nyquist图,可以由机器确定频率范围,也可以自选频率范围。
[re,im,w]=nyquist(sys)
分别返回系统sys对应频率w的实部(re)及虚部(im)值。
读取实部,虚部及对应频率点的方法:
re(1,1,k)
im(1,1,k)
w(k) (rad/s)
其中,k为机器计算Nyquist响应的序列编号,必须是正整数而且应当在计算范围内。可以使用命令size(w)查看。k值随问题的简、繁而变化。
这是一个非常有用的命令,当然也可以采用直接在Nyquist图上用鼠标读取的方法。
[re,im] = nyquist(sys,w)
返回在某一频率w(rad/s)或某一频率范围w ={ wmin, wmax }系统的实部及虚部值。当选定为频率范围时,调出实部,虚部及对应频率值的方法同前一指令。
例3.1 分析系统
的频率特性。
该传递函数是一个零极点及有理分式的混合形式。首先将子系统
变换成零极点增益形式[z1,p1,k1],再将系统变换成零极点增益模型(z,p,k)再作Nyquist图。
使用[re,im,w]=nyquist(sys)语句,可计算出若干频率点所对应的实部、虚部值,而且在
w(25)=0.8802rad/s
w(26)=0.8951rad/s
之间,曲线将从第三象限穿过负实轴进入第二象限,在w(0.8802,0.8951)之间使用语句
[re,im]=nyquist(sys,w)
[0.8802,Δw,0.8951]
进行更细致的计算,在需要的精度内求出穿越负实轴的参数。
参考程序见fanli019。
参考程序fanli019:用Nyquist图分析系统频率特性
num=5;den=[1 0.2 1];%子系统的传递函数
[Z1,P1,K1]=tf2zp(num,den) %子系统的极点
Z=[-2.5;-2] %系统的零点
P=[P1;1/1.2;-1/1.4] %系统的极点
K=num %系统的增益
sys=zpk(Z,P,K) %系统的零、极点、增益模型
nyquist(sys) %系统的Nyquist图
[re,im,w]=nyquist(sys) %计算系统若干频率点所对应的实部与虚部值
w(25) %w(25)=0.8802,曲线处于第三象限靠近负实轴的频率点
w(26) %w(26)=0.8951,曲线处于第二象限靠近负实轴的频率点
[re1,im1]=nyquist(sys,[0.8802:0.00005:0.8951])
%在(0.8802, 0.8951)(rad/s)范围内,以0.00005rad/s为间隔计算实部及虚部值
fanli019的Nyquist图
2. 带延时环节系统的Nyquist图
对带有延时环节的系统,由于不能表示成传递函数的有理分式,所以MATLAB中没有直接命令可用。但由于Nyquist图实际上是以频率ω为参变量,横坐标为实频特性,纵坐标为虚频特性的直角坐标系图。所以,可以通过计算实频特性与虚频特性,再使用plot(x,y)作出Nyquist图。
例3.2 求作具有延时环节的系统
的Ny
展开阅读全文