matlab在控制中的应用.docx
《matlab在控制中的应用.docx》由会员分享，可在线阅读，更多相关《matlab在控制中的应用.docx（62页珍藏版）》请在冰豆网上搜索。
matlab在控制中的应用
《控制工程基础（经典控制部分）》
的MATLAB分析
机电工程系
郭天石 建明
轻化工学院
二零零二年九月
前 言
在学习《工程控制基础》及其它相关的控制理论课程时，许多计算问题越来越离不开计算机辅助分析及计算机辅助设计。
在计算机性能迅速提高和各种软件大量涌现的今天，机辅分析与机辅设计也越来越方便。
我们希望，在进行机辅分析及机辅设计时，把主要精力放在对计算机计算处理后的数据、曲线的控制特性分析上，而不在编程及调试程序上花费太多的时间，使计算机真正成为学习相关课程和解决实际问题的工具。
因此，我们编写这本在MATLAB语言支持下进行控制方面的机辅分析及机辅设计的讲义，力图从实用的角度，介绍MATLAB语言在经典控制方面的一些典型应用。
第一章 MATLAB的基本使用
§1-1 MATLAB语言简介
MATLAB是一种高级矩阵语言，它由MathWorks公司于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文件的执行：
返回命令窗口，在当前目录（CurrentDirectory）选择所要运行的M文件的目录，在命令窗口提示符“>>”后，直接键入文件名（不加后缀）即可运行。
注意：
（1）机器默认路径为一级子目录MATLAB6p1\work；
（2）MATLAB6.1以前的版本，运行M文件的方法稍有不同，它必须在File菜单下，打开“RunScript…”子菜单，键入需要运行的文件路径及名称再回车，在这种情况下，work作为根目录对待，不出现在M文件的路径之中。
本讲义的参考程序都是在M文件窗口下编制的。
三、在线帮助窗口
在命令窗口中键入Help（空格）函数名，可以立即获得该函数的使用方法。
§1-3MATLAB最基本的矩阵操作
作为命令窗口及M文件编辑器的应用实例，介绍几个最基本的矩阵运算命令。
一、矩阵的输入
在方括号依次按行键入矩阵元素，在一行的各元素之间用空格或逗号分开，每行之间用分号分开。
例如，在命令窗输入
A＝[223;454;789]↙（注意：
方括号，分号为矩阵行标记）
B＝[1,3,5;6,-4,2;3,5,1]↙（逗号与空格功能相同）
A＝2 23B＝135
4546-42
789351
同理：
输入A1＝［2 4 6］得到行矢量，
输入A2＝［2；4；6］得到列矢量，
于是，当输入
C＝［A；A1］↙ 有
C＝123
456
789
246
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=[223;454;789]%三阶矩阵输入
B=[1,3,5;6,-4,2;3,5,1]%三阶矩阵输入
A1=[246]%行向量
A2=[2;4;6]%列向量
C=[A;A1]%矩阵A增加一行
C1=[AA2] %矩阵A增加一列
C2=[AA1'] %矩阵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的符号运算操作
一、进入符号运算功能
在命令窗口键入
symsxyzt↙
此后，即可以使用x，y，z，t等作自变量定义函数。
symsxyztreal↙规定所定义的变量为实型。
二、代数方程求解
使用命令solve可以求解代数方程，如求下例方程组
的解，命令为
[x,y]=solve（‘x^2+x*y+y-3=0,x^2-4*x-2*y+3=0’）
程序见例程序fanli002。
参考程序fanli002:
代数方程求解
symsxy %进入符号运算功能
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:
符号矩阵的特征值
symstreal %定义为实型变量
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中微积分运算命令为
diff（f）：
求函数f对自变量x的一阶导数；
diff（f，2）：
求函数f对自变量x的二阶导数；
int（f）：
求函数f的不定积分
例1.1，设
试计算其一阶，二阶导数，积分运算，并作出函数图象，见例fanli004。
参考程序fanli004:
函数的微分与积分
symsx
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πfanli004的函数曲线
fanli004的一阶导函数曲线
fanli004的二阶导函数曲线
fanli004的二阶导函数的二重积分曲线
ezplot（f）是一个很有用的作图命令，它的其它应用形式，请查在线帮助。
细心的读者会发现，一个函数求二阶导数后再对二阶导数进行二重积分，其结果与原函数相差一个常数。
相当于纵坐标发生平移。
第二章 系统的时域特性
§2-1 传递函数
一、传递函数的两种形式
传递函数通常表达成s的有理分式形式及零极点增益形式。
设传递函数
1．有理分式形式
分别将分子、分母中s多项式的系数按降幂排列成行矢量，缺项的系数用0补齐。
上述函数可表示为
num1=[21]（注意：
方括号，同一行的各元素间留空格或逗号）。
den1=[1221]
syss1=tf（num1,den1）
运行后，返回传递函数
的形式。
这种形式不能直接进行符号运算！
2．零极点增益形式
[Z，P，K]=tf2zp（num1，den1）
sys2=zpk（Z，P，K）
返回零、极点、增益表达式，其Z，P分别将零点和极点表示成列向量，若无零点或极点用[]（空矩阵）代替。
运行得到
的
1点Z=-0.5
极点 P=-1，-0.5±j0.866
增益 K=2
指令zp2tf（Z，P，K）将零极点增益变换成有理分式形式，见程序fanli005。
参考程序fanli005:
传递函数的有理分式及零极点增益模型
num1=[21]%传递函数的分子系数向量
den1=[1221]%传递函数的分母系数向量
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是状态空间模型与有理分式及零、极点、增益模型之间的相互转换。
二、传递函数框图的处理
用框图可以方便地表示传递函数的并联，串联及反馈。
为简洁，仅以有理分式模型为例。
1．并联
sysp=parallel（sys1，sys2）
[num，den]＝parallel（num1，den1，num2，den2）
2．串联
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．单位反馈
sysd=feedback（syss,1）
[numd,dend]=feedback（nums,dens,1,1）（单位反馈）
上面给出了同一指令的两种形式，相当于两套平行指令。
对于零极点增益形式，书写稍复杂一些，可先用zpk转换成系统形式，或用zp2tf转折换成有理分式形式后再进行框图化简操作。
三、简单函数的拉普拉斯变换
在MATLAB的符号功能中，可以对简单函数进行拉普拉斯正、逆变换。
拉氏正变换：
laplace（f（t））
拉氏逆变换：
ilaplace（L（s））
其中
为原函数，
为象函数。
命令格式参见fanli007。
参考程序fanli007:
拉普拉斯变换
symsstwabc
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=[42]
den1=[2814114]
sys1=tf（num1,den1）%系统G1（s）
num2=[21]
den2=[14673]
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=[42]
den1=[2814114]
sys1=tf（num1,den1）
num2=[21]
den2=[14673]
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）的单位脉冲响应曲线
holdon,step（sys2）%系统G2（s）的单位阶跃和单位脉冲响应曲线
fanli009:
impulse（sys1,sys2）单位脉冲响应曲线
程序的最后一句holdon是当前图形保护模式。
当要将新图形作在当前图形上时，必须使用holdon。
而figure的含意是另开一个新的图形窗口，如果不用figure或holdon，则新的图形会占用原图形窗口，始终只保留一个最新的图形窗口。
三、一阶系统及二阶系统的时域特性
一阶系统及二阶系统是最基本也是最重要的系统，高阶系统总可以视为由若干个一阶和（或）二阶系统组合构成。
1．一阶系统（设增益为1）
影响系统特性的参数是其时间常数T，T越大，系统惯性越大，响应越慢。
参考程序fanli010给出了T＝0.4,1.2,2.0,2.8,3.6,4,4六条单位阶跃响应曲线。
参考程序fanli010:
一阶系统的单位阶跃响应曲线
num=1;i=1;
fordel=0.1:
0.2:
1.1%一阶系统时间常数递增间隔
den=[4*del1];%一阶系统分母向量
step（tf（num,den））%一阶系统单位阶跃响应曲线
holdon,%不同时间常数的一阶系统单位阶跃响应曲线簇
i=i+1;
end
同理，可以作出对应的单位脉冲响应曲线，参考程序fanli011。
参考程序fanli011:
一阶系统的单位脉冲响应曲线
num=1;i=1;
fanli010一阶系统时间常数对单位阶跃响应的影响
fanli011一阶系统时间常数对单位脉冲响应的影响
fordel=0.1:
0.2:
1.1
den=[4*del1];
impulse（tf（num,den）,10）%一阶系统单位阶脉冲应曲线
holdon,%不同时间常数的一阶系统单位脉冲响应曲线簇
i=i+1;
end
注意MATLAB中for语句的结构。
读者可以改变不同的增益，看看图形有何变化。
2．二阶系统（设0<ξ<1）
设二阶系统为
二阶系统的特征参数为固有频率
及阻尼比ξ。
当
增大，系统振动频率加快，振荡加剧；而随着ξ减小，系统振荡加剧，振荡峰尖锐。
参考程序fanli012示出了当ξ＝0.5,
=1,2,3,4,5rad/s时的间接阶跃曲线簇。
参考程序fanli012:
不同固有频率的二阶系统的单位阶跃响应曲线（ξ=0.5）
i=1;
fordel=1:
1:
5;%二阶系统固有频率递增间隔
num=del^2;%二阶系统传递函数分子系数向量
den=[1deldel^2];%不同固有频率的二阶系统分母系数向量
step（tf（num,den）,6）%二阶系统单位阶跃响应曲线
holdon,%不同固有频率的二阶系统单位阶跃响应曲线簇
i=i+1;
end
fanli012二阶系统固有频率对单位阶跃响应的影响
参考程序fanli013示出了同一二阶系统的单位脉冲响应曲线簇。
参考程序fanli013:
不同固有频率的二阶系统的单位脉冲响应曲线（ξ=0.5）
i=1;
fordel=1:
1:
5;
num=del^2;
den=[1deldel^2];%不同固有频率的二阶系统分母系数向量
impul
升级会员 