控制系统的MATLAB仿真算法
1 MATLAB简介编程
MATLAB是Mathworks公司开发的一种集数值计算、符号计算和图形可视化三大基本功能于一体的功能强大、操做简单的优秀工程计算应用软件。MATLAB不只能够处理代数问题和数值分析问题,并且还具备强大的图形处理及仿真模拟等功能。从而可以很好的帮助工程师及科学家解决实际的技术问题。数组
MATLAB的含义是矩阵实验室(Matrix Laboratory),最初主要用于方便矩阵的存取,其基本元素是无需定义维数的矩阵。通过十几年的扩充和完善,现已发展成为包含大量实用工具箱(Toolbox)的综合应用软件,不只成为线性代数课程的标准工具,并且适合具备不一样专业研究方向及工程应用需求的用户使用。浏览器
MATLAB最重要的特色是易于扩展。它容许用户自行创建完成指定功能的扩展MATLAB函数(称为M文件),从而构成适合于其它领域的工具箱,大大扩展了MATLAB的应用范围。目前,MATLAB已成为国际控制界最流行的软件,控制界不少学者将本身擅长的CAD方法用MATLAB加以实现,出现了大量的MATLAB配套工具箱,如控制系统工具箱(control systems toolbox),系统识别工具箱(system identification toolbox),鲁棒控制工具箱(robust control toolbox),信号处理工具箱(signal processing toolbox)以及仿真环境SIMULINK等。网络
(1) MATLAB的安装jsp
本节将讨论操做系统为Microsoft Windows环境下安装MATLAB6的过程。编辑器
将MATLAB6的安装盘放入光驱,系统将自动运行auto-run.bat文件,进行安装;也能够执行安装盘内的setup.exe文件启动MATLAB的安装程序。启动安装程序后,屏幕将显示安装MATLAB的初始界面,根据Windows安装程序的常识,不断单击[Next],输入正确的安装信息,具体操做过程以下:ide
输入正确的用户注册信息码;函数
选择接收软件公司的协议;工具
输入用户名和公司名;
选择MATLAB组件(Toolbox);
选择软件安装路径和目录;
单击[Next]按钮进入正式的安装界面。安装过程界面如图1所示。
图1 MATLAAB安装过程界面
图2 MATLAAB启动过程界面
安装完毕后,选择[Restart my computer now]选项以从新启动计算机。
从新启动计算机后,用户就能够点击图标使用MATLAB6了。MATLAB启动过程界面如图2所示。
(2) MATLAB桌面系统
MATLAB的桌面系统由桌面平台以及桌面组件共同构成,如图3。桌面平台是各桌面组件的展现平台,它提供了一系列的菜单操做以及工具栏操做,而不一样功能的桌面组件构成了整个MATLAB操做平台。其组件主要包含以下8个组件部分:
①命令窗口(Command Window)②历史命令窗口(Command History)③组件平台(Launch Pad)④路径浏览器(Current Directory Browser)⑤帮助浏览器(Help Browser)⑥工做空间浏览器(Workspace Browser)⑦数组编辑器(Array Editor)⑧M文件编辑调试器(Editor-Debugger)。
用户能够在View菜单下选择打开或关闭某个窗口。
图3 MATLAB桌面平台
(3) MATLAB命令窗口
MATLAB能够认为是一种解释性语言。在MATLAB命令窗口中,标志>>为命令提示符,在命令提示符后面键入一个MATLAB命令时,MATLAB会当即对其进行处理,并显示处理结果。
这种方式简单易用,但在编程过程当中要修改整个程序比较困难,而且用户编写的程序不容易保存。若是想把全部的程序输入完再运行调试,能够用鼠标点击快捷或File|New|M-file菜单,在弹出的编程窗口中逐行输入命令,输入完毕后点击Debug|Run(或F5)运行整个程序。运行过程当中的错误信息和运行结果显示在命令窗口中。整个程序的源代码能够保存为扩展名为".m"的M文件。
在介绍MATLAB的强大计算和图象处理功能前,咱们能够先运行一个简单的程序。
设系统的闭环传递函数为:
求系统的时域响应图,可输入下面的命令:
>> num=[1,4];
den=[1,2,8];
step(num,den)
图4 动态响应时域图
程序运行后会在一个新的窗口中显示出系统的时域动态响应曲线,如图4。用鼠标左键点击动态响应曲线的某一点,系统会提示其响应时间和幅值。按住左键在曲线上移动鼠标的位置能够很容易的根据幅值观察出上升时间、调节时间、峰值及峰值时间,进而求出超调量。若是想求根轨迹,可将程序的第三行变为rlocus(num,den),求伯德图可改成bode(num,den)。所不一样的是,在根轨迹和伯德图中,G(s)为开环传递函数。
MATLAB的语法规则相似于C语言,变量名、函数名都与大小写有关,即变量A和a是两个彻底不一样的变量。应该注意全部的函数名均由小写字母构成。
MATLAB是一个功能强大的工程应用软件,它提供了至关丰富的帮助信息,同时也提供了多种得到帮助的方法。若是用户第一次使用MATLAB,则建议首先在>>提示符下键入DEMO命令,它将启动MATLAB的演示程序。用户能够在此演示程序中领略MATLAB所提供的强大的运算和绘图功能。
2 MATLAB基本操做命令
本节简单介绍与本书内容相关的一些基本知识和操做命令。
(1)简单矩阵的输入
MATLAB是一种专门为矩阵运算设计的语言,因此在MATLAB中处理的全部变量都是矩阵。这就是说,MATLAB只有一种数据形式,那就是矩阵,或者数的矩形阵列。标量可看做为1×1的矩阵,向量可看做为n×1或1×n的矩阵。这就是说,MATLAB语言对矩阵的维数及类型没有限制,即用户无需定义变量的类型和维数,MATLAB会自动获取所需的存储空间。
输入矩阵最便捷的方式为直接输入矩阵的元素,其定义以下:
例如,在MATLAB的工做空间中,输入:
>>
则输出结果为:
矩阵a被一直保存在工做空间中,以供后面使用,直至修改它。
MATLAB的矩阵输入方式很灵活,大矩阵能够分红n行输入,用回车符代替分号或用续行符号(…)将元素续写到下一行。例如:
以上三种输入方式结果是相同的。通常若长语句超出一行,则换行前使用续行符号(…)。
在MATLAB中,矩阵元素不限于常量,能够采用任意形式的表达式。同时,除了直接输入方式以外,还能够采用其它方式输入矩阵,如:
(2)复数矩阵输入
MATLAB容许在计算或函数中使用复数。输入复数矩阵有两种方法:
(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 是MATLAB提供的固定变量,具备特定的功能,是不能由用户清除的。经常使用的固定变量还有eps、pi、Inf、NaN等。其特殊含义能够用7.2.10节介绍的方法查阅帮助。
MATAB容许在函数调用时同时返回多个变量,而一个函数又能够由多种格式进行调用,语句的典型格式可表示为:
[返回变量列表]=fun-name(输入变量列表)
例如用bode()函数来求取或绘制系统的Bode图,可由下面的格式调用:
其中变量num、den表示系统传递函数分子和分母,W表示指定频段,mag为计算幅值,phase为计算相角。
(4)语句以"%"开始和以分号";"结束的特殊效用
在MATLAB中以"%"开始的程序行,表示注解和说明。符号"%"相似于C++中的"//"。这些注解和说明是不执行的。这就是说,在MATLAB程序行中,出现"%"之后的一切内容都是能够忽略的。
分号用来取消打印,若是语句最后一个符号是分号,则打印被取消,可是命令仍在执行,而结果再也不在命令窗口或其它窗口中显示。这一点在M文件中大量采用,以抑制没必要要的信息显示。
(5)获取工做空间信息
MATLAB开辟有一个工做空间,用于存储已经产生的变量。变量一旦被定义,MATLAB系统会自动将其保存在工做空间里。在退出程序以前,这些变量将被保留在存储器中。
为了获得工做空间中的变量清单,能够在命令提示符>>后输入who 或 whos 命令,当前存放在工做空间的全部变量便会显示在屏幕上。
命令clear能从工做空间中清除全部非永久性变量。若是只须要从工做空间中清除某个特定变量,好比"x",则应输入命令clear x。
(6)常数与算术运算符
MATLAB采用人们习惯使用的十进制数。如:
3 –99 0.0001 9.6397238
2i -3.14159i 3e5i
其中 。
数值的相对精度为eps,它是一个符合IEEE标准的16位长的十进制数,其范围为:。
MATLAB提供了经常使用的算术运算符:+,-,,∕(﹨),^(幂指数)。
应该注意:(∕)右除法和(﹨)左除法这两种符号对数值操做时,其结果相同,其斜线下为分母,如1∕4与4﹨1,其结果均为0.25,但对矩阵操做时,左、右除法是有区别的。
(7)选择输出格式
输出格式是指数据显示的格式,MATLAB提供format命令能够控制结果矩阵的显示,而不影响结果矩阵的计算和存储。全部计算都是以双精度方式完成的。
如输入:
则显示:
① format short(短格式,也是系统默认格式)
② format short e(短格式科学表示)
③ format long(长格式)
④ format long e(长格式科学表示)
如:
对于以上四种格式,其显示结果分别为:
短格式5位表示
短格式科学表示
长格式16位表示
长格式科学表示
一旦调用了某种格式,则这种被选用的格式将保持,直到对格式进行了改变为止。
(8)MATLAB图形窗口
当调用了一个产生图形的函数时,MATLAB会自动创建一个图形窗口。这个窗口还可分裂成多个窗口,并可在它们之间选择,这样在一个屏上可显示多个图形。
图形窗口中的图形可经过打印机打印出来。若想将图形导出并保存,可用鼠标点击菜单File|Export,导出格式可选emp、bmp、jpg等。命令窗口的内容也可由打印机打印出来:若是事先选择了一些内容,则可打印出所选择的内容;若是没有选择内容,则可打印出整个工做空间的内容。
(9)剪切板的使用
利用Windows的剪切板可在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。
% 用于求取一阶跃响应。
num=[1 4];
den=[1 2 8];
step(num,den)
当你键入help Step_Response时,屏幕上将显示文件开头部分的注释:
用于求取一阶跃响应。
很显然,在每个M文件的开头,创建详细的注释是很是有用的。因为MATLAB提供了大量的命令和函数,想记住全部函数及调用方法通常不太可能,经过联机帮助命令help可容易地对想查询的各个函数的有关信息进行查询。该命令使用格式为:
help 命令或函数名
注意:若用户把文件存放在本身的工做目录上,在运行以前应该使该目录处在MATLAB的搜索路径上。当调用时,只需输入文件名,MATLAB就会自动按顺序执行文件中的命令。
函数就像一个黑箱,把一些数据送进去,经加工处理,再把结果送出来。在函数体内使用的除返回变量和输入变量这些在第一行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.有理函数模型
线性系统的传递函数模型可通常地表示为:
(1)
将系统的分子和分母多项式的系数按降幂的方式以向量的形式输入给两个变量和,就能够轻易地将传递函数模型输入到MATLAB环境中。命令格式为:
; (2)
; (3)
在MATLAB控制系统工具箱中,定义了tf() 函数,它可由传递函数分子分母给出的变量构造出单个的传递函数对象。从而使得系统模型的输入和处理更加方便。
该函数的调用格式为:
G=tf(num,den); (4)
例2 一个简单的传递函数模型:
能够由下面的命令输入到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 一个稍微复杂一些的传递函数模型:
该传递函数模型能够经过下面的语句输入到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.零极点模型
线性系统的传递函数还能够写成极点的形式:
(5)将系统增益、零点和极点以向量的形式输入给三个变量、Z和P,就能够将系统的零极点模型输入到MATLAB工做空间中,命令格式为:
(6) (7) (8)
在MATLAB控制工具箱中,定义了zpk()函数,由它可经过以上三个MATLAB变量构造出零极点对象,用于简单地表述零极点模型。该函数的调用格式为:
G=zpk(Z,P,KGain) (9)
例4 某系统的零极点模型为:
该模型能够由下面的语句输入到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所示。
图5 反馈系统结构图
控制系统工具箱中提供了feedback()函数,用来求取反馈链接下总的系统模型,该函数调用格式以下:
G=feedback(G1,G2,sign); (10)
其中变量sign用来表示正反馈或负反馈结构,若sign=-1表示负反馈系统的模型,若省略sign变量,则仍将表示负反馈结构。G1和G2分别表示前向模型和反馈模型的LTI(线性时不变)对象。
例5 若反馈系统图5中的两个传递函数分别为:
,
则反馈系统的传递函数可由下列的MATLAB命令得出
>> G1=tf(1,[1,2,1]);
G2=tf(1,[1,1]);
G=feedback(G1,G2)
运行结果:
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所示。其中
,,
则闭环系统的传递函数能够由下面的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
图6 复杂反馈系统
4. 有理分式模型与零极点模型的转换
有了传递函数的有理分式模型以后,求取零极点模型就不是一件困难的事情了。在控制系统工具箱中,能够由zpk()函数当即将给定的LTI对象G转换成等效的零极点对象G1。该函数的调用格式为:
G1=zpk(G) (11)
例7 给定系统传递函数为:
对应的零极点格式可由下面的命令得出
>> num=[6.8, 61.2, 95.2];
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 给定零极点模型:
能够用下面的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按钮,打开一个空白的模型编辑窗口如图 8所示。
图 7 simulink 模型库
图8 模型编辑窗口
(2) 画出系统的各个模块:打开相应的子模块库,选择所须要的元素,用鼠标左键点中后拖到模型编辑窗口的合适位置。
(3) 给出各个模块参数:因为选中的各个模块只包含默认的模型参数,如默认的传递函数模型为1/(s+1)的简单格式,必须经过修改获得实际的模块参数。要修改模块的参数,能够用鼠标双击该模块图标,则会出现一个相应对话框,提示用户修改模块参数。
(4) 画出链接线:当全部的模块都画出来以后,能够再画出模块间所须要的连线,构成完整的系统。模块间连线的画法很简单,只须要用鼠标点按起始模块的输出端(三角符号),再拖动鼠标,到终止模块的输入端释放鼠标键,系统会自动地在两个模块间画出带箭头的连线。若须要从连线中引出节点,可在鼠标点击起始节点时按住Ctrl键,再将鼠标拖动到目的模块。
例9 典型二阶系统的结构图如图9所示。用SIMULINK对系统进行仿真分析。
图9 典型二阶系统结构图
按前面步骤,启动simulink并打开一个空白的模型编辑窗口。
(3)将画出的全部模块按图9用鼠标链接起来,构成一个原系统的框图描述如图10所示。
(4)选择仿真算法和仿真控制参数,启动仿真过程。
图10 二阶系统的simulink实现
在命令窗口中键入whos命令,会发现工做空间中增长了两个变量――tout和yout,这是由于Simulink中的Out1 模块自动将结果写到了MATLAB的工做空间中。利用MATLAB命令plot(tout,yout),可将结果绘制出来,如图12所示。比较11和12,能够发现这两种输出结果是彻底一致的。
图11仿真结果示波器显示
图12 MATLAB命令得出的系统响应曲线
(2) 利用MATLAB进行时域分析
1. 线性系统稳定性分析
线性系统稳定的充要条件是系统的特征根均位于S平面的左半部分。系统的零极点模型能够直接被用来判断系统的稳定性。另外,MATLAB语言中提供了有关多项式的操做函数,也能够用于系统的分析和计算。
(1)直接求特征多项式的根
设p为特征多项式的系数向量,则MATLAB函数roots()能够直接求出方程p=0在复数范围内的解v,该函数的调用格式为:
v=roots(p) (13)
例10 已知系统的特征多项式为:
特征方程的解可由下面的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)
结果显示
ans =
8
(4)部分分式展开
考虑下列传递函数:
式中,可是和中某些量可能为零。
MATLAB函数可将展开成部分分式,直接求出展开式中的留数、极点和余项。该函数的调用格式为:
(16)
则的部分分式展开由下式给出:
式中,,…, ,为极点,
,,…, 为各极点的留数,为余项。
例11 设传递函数为:
该传递函数的部分分式展开由如下命令得到:
>> 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。
由此可得出部分分式展开式:
该函数也能够逆向调用,把部分分式展开转变回多项式之比的形式,命令格式为:
[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重极点。在这种状况下,部分分式展开式将包括下列诸项:
例12 设传递函数为:
则部分分式展开由如下命令得到:
>> 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=[]为空矩阵。
由上可得展开式为:
(5)由传递函数求零点和极点。
在MATLAB控制系统工具箱中,给出了由传递函数对象G求出系统零点和极点的函数,其调用格式分别为:
Z=tzero(G) (18)
P=G.P{1} (19)
注意:式19中要求的G必须是零极点模型对象,且出现了矩阵的点运算"."和大括号{}表示的矩阵元素,详细内容参阅后面章节。
例13 已知传递函数为:
输入以下命令:
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 给定传递函数:
利用下列命令可自动打开一个图形窗口,显示该系统的零、极点分布图,如图13所示。
>> num=[3,2,5,4,6];
den=[1,3,4,2,7,2];
pzmap(num,den)
title(?Pole-Zero Map?) % 图形标题。
图13 MATLAB函数零、极点分布图
2. 系统动态特性分析。
(1)时域响应解析算法――部分分式展开法。
用拉氏变换法求系统的单位阶跃响应,可直接得出输出c(t)随时间t变化的规律,对于高阶系统,输出的拉氏变换象函数为:
(21)
对函数c(s)进行部分分式展开,咱们能够用num,[den,0]来表示c(s)的分子和分母。
例 15 给定系统的传递函数:
用如下命令对进行部分分式展开。
>> 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)为:
拉氏变换得:
(2)单位阶跃响应的求法:
控制系统工具箱中给出了一个函数step()来直接求取线性系统的阶跃响应,若是已知传递函数为:
则该函数可有如下几种调用格式:
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 已知传递函数为:
利用如下MATLAB命令可得阶跃响应曲线如图14所示。
图14 MATLAB绘制的响应曲线
>> num=[0,0,25];
den=[1,4,25];
step(num,den)
grid % 绘制网格线。
title(?Unit-Step Response of G(s)=25/(s^2+4s+25) ?) % 图像标题
咱们还能够用下面的语句来得出阶跃响应曲线
>> 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<循环判断语句>
循环体
end
其中,循环判断语句为某种形式的逻辑判断表达式。
当表达式的逻辑值为真时,就执行循环体内的语句;当表达式的逻辑值为假时,就退出当前的循环体。若是循环判断语句为矩阵时,当且仅当全部的矩阵元素非零时,逻辑表达式的值为真。为避免循环语句陷入死循环,在语句内必须有能够自动修改循环控制变量的命令。
要求出上升时间,能够用while语句编写如下程序获得:
C=dcgain(G);
n=1;
while y(n)<C
n=n+1;
end
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 已知二阶系统传递函数为:
利用下面的stepanalysis.m程序可获得阶跃响应如图 15及性能指标数据。
>> G=zpk([ ],[-1+3*i,-1-3*i ],3);
% 计算最大峰值时间和它对应的超调量。
C=dcgain(G)
[y,t]=step(G);
plot(t,y)
grid
[Y,k]=max(y);
timetopeak=t(k)
percentovershoot=100*(Y-C)/C
% 计算上升时间。
n=1;
while y(n)<C
n=n+1;
end
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绘制系统根轨迹
假设闭环系统中的开环传递函数能够表示为:
则闭环特征方程为:
特征方程的根随参数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,??) (37)
plot()函数里引号内的部分用于选择所绘制曲线的类型,详细内容见表1。控制系统工具箱中还有一个rlocfind()函数,该函数容许用户求取根轨迹上指定点处的开环增益值,并将该增益下全部的闭环极点显示出来。这个函数的调用格式为:
[K,P]=rlocfind(G) (38)
这个函数运行后,图形窗口中会出现要求用户使用鼠标定位的提示,用户能够用鼠标左键点击所关心的根轨迹上的点。这样将返回一个K变量,该变量为所选择点对应的开环增益,同时返回的P变量则为该增益下全部的闭环极点位置。此外,该函数还将自动地将该增益下全部的闭环极点直接在根轨迹曲线上显示出来。
例18 已知系统的开环传递函数模型为:
利用下面的MATLAB命令可容易地验证出系统的根轨迹如图16所示。
>> G=tf(1,[conv([1,1],[1,2]),0]);
rlocus(G);
grid
title(?Root_Locus Plot of G(s)=K/[s(s+1)(s+2)]?)
xlabel(?Real Axis?) % 给图形中的横坐标命名。
ylabel(?Imag Axis?) % 给图形中的纵坐标命名。
[K,P]=rlocfind(G)
用鼠标点击根轨迹上与虚轴相交的点,在命令窗口中可发现以下结果
select_point=0.0000+1.3921i
K=
5.8142
p=
-2.29830
-0.0085+1.3961i
-0.0085-1.3961i
因此,要想使此闭环系统稳定,其增益范围应为0<K<5.81。
参数根轨迹反映了闭环根与开环增益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值,振荡加重,性能变坏。
图17 不一样K值下的阶跃响应曲线
4. MATLAB绘图的基本知识
经过以上实例的应用,咱们已初步尝试了MATLAB的绘图功能。MATLAB具备丰富的获取图形输出的程序集。咱们已用命令plot()产生线性x-y图形(用命令loglog、semilogx、semilogy或polar取代命令plot,能够产生对数坐标图和极坐标图)。全部这些命令的应用方式都是类似的,它们只是在如何给坐标轴进行分度和如何显示数据上有所差异。
(1)二维图形绘制
若是用户将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一个周期内的正弦曲线
(2)一幅图上画多条曲线。
利用具备多个输入变量的plot( )命令,能够在一个绘图窗口上同时绘制多条曲线,命令格式为:
plot(x1,Y1,x2,Y2,…,xn,Yn) (40)
x一、Y一、x二、Y2等一系列变量是一些向量对,每个x-y对均可以用图解表示出来,于是能够在一幅图上画出多条曲线。多重变量的优势是它容许不一样长度的向量在同一幅图上显示出来。每一对向量采用不一样的线型以示区别。
另外,在一幅图上叠画一条以上的曲线时,也能够利用hold命令。hold命令能够保持当前的图形,而且防止删除和修改比例尺。所以,后来画出的那条曲线将会重叠在原曲线图上。当再次输入命令hold,会使当前的图形复原。也能够用带参数的hold命令――hold on 和hold off来启动或关闭图形保持。
(3)图形的线型和颜色
为了区分多幅图形的重叠表示,MATLAB提供了一些绘图选项,能够用不一样的线型或颜色来区分多条曲线,经常使用选项见下表1
表1MATLAB绘图命令的多种选项
选项 |
意义 |
选项 |
意义 |
′-′ |
实线 |
′--′ |
短划线 |
′: ′ |
虚线 |
′-.′ |
点划线 |
′r′ |
红色 |
′*′ |
用星号绘制各个数据点 |
′b′ |
蓝色 |
′o′ |
用圆圈绘制各个数据点 |
′g′ |
绿色 |
′.′ |
用圆点绘制各个数据点 |
′y′ |
黄色 |
′×′ |
用叉号绘制各个数据点 |
表1中绘出的各个选项有一些能够并列使用,可以对一条曲线的线型和颜色同时做出规定。例如′--g′表示绿色的短划线。带有选项的曲线绘制命令的调用格式为:
plot(X1,Y1,S1,X2,Y2,S2,…) (41)
(4)加进网络线、图形标题、x轴和y轴标记
一旦在屏幕上显示出图形,就能够依次输入如下相应的命令将网络格线、图形标题、x、y轴标记叠加在图形上。命令格式以下:
grid(网络线) (42)
title(′图形标题′) (43)
xlabel(′x轴标记′) (44)
ylabel(′y轴标记′) (45)
函数引号内的字符串将被写到图形的坐标轴上或标题位置。
(5)在图形屏幕上书写文本。
若是想在图形窗口中书写文字,能够单击按钮,选择屏幕上一点,点击鼠标,在光标处输入文字。另外一种输入文字的方法是用text()命令。它能够在屏幕上以(x,y)为坐标的某处书写文字,命令格式以下:
text(x,y,′text′) (46)
例如,利用语句
text(3,0.45,′sint′)
将从点(3,0.45)开始,水平的写出"sint"。
(6)自动绘图算法及手工坐标轴定标
在MATLAB图形窗口中,图形的横、纵坐标是自动标定的,在另外一幅图形画出以前,这幅图形做为现行图将保持不变,可是在另外一幅图形画出后,原图形将被删除,坐标轴自动地从新标定。关于瞬态响应曲线、根轨迹、伯德图、奈魁斯特图等的自动绘图算法已经设计出来,它们对于各种系统具备普遍的适用性,可是并不是老是理想的。所以,在某些状况下,可能须要放弃绘图命令中的坐标轴自动标定特性,由用户本身设定坐标范围,能够在程序中加入下列语句:
v=[x-min x-max y-min y-max] (47)
axis(v) (48)
式中v是一个四元向量。axis(v)把坐标轴定标创建在规定的范围内。对于对数坐标图,v的元素应为最小值和最大值的经常使用对数。
执行axis(v)会把当前的坐标轴标定范围保持到后面的图中,再次键入axis可恢复系统的自动标定特性。
Axis(′sguare′)可以把图形的范围设定在方形范围内。对于方形长宽比,其斜率为1的直线恰位于45o 上,它不会因屏幕的不规则形状而变形。Axis(′normal′)将使长宽比恢复到正常状态。
5、线性系统的频域分析。
(1)频率特性函数。
设线性系统传递函数为:
则频率特性函数为:
由下面的MATLAB语句可直接求出G(jw)。
i=sqrt(-1) % 求取-1的平方根 (49)
GW=polyval(num,i*w)./polyval(den,i*w) (50)
其中(num,den)为系统的传递函数模型。而w为频率点构成的向量,点右除(./)运算符表示操做元素点对点的运算。从数值运算的角度来看,上述算法在系统的极点附近精度不会很理想,甚至出现无穷大值,运算结果是一系列复数返回到变量GW中。
(2)用MATLAB做奈魁斯特图。
控制系统工具箱中提供了一个MATLAB函数nyquist( ),该函数能够用来直接求解Nyquist 阵列或绘制奈氏图。当命令中不包含左端返回变量时,nyquist()函数仅在屏幕上产生奈氏图,命令调用格式为:
nyquist(num,den) (51)
nyquist(num,den,w) (52)
或者
nyquist(G) (53)
nyquist(G,w) (54)
该命令将画出下列开环系统传递函数的奈氏曲线:
若是用户给出频率向量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 考虑二阶典型环节:
试利用MATLAB画出奈氏图。
利用下面的命令,能够得出系统的奈氏图,如图19所示。
>> num=[0,0,1];
den=[1,0.8,1];
nyquist(num,den)
% 设置坐标显示范围
v=[-2,2,-2,2];
axis(v)
grid
图19 二阶环节奈氏图
title(′Nyquist Plot of G(s)=1/(s^2+0.8s+1) ′)
(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等分的,即在两个十进制数和 之间产生一个由50个点组成的份量,向量中的点数50是一个默认值。例如要在弧度/秒与弧度/秒之间的频区画伯德图,则输入命令时,,在此频区自动按对数距离等分红50个频率点,返回到工做空间中,即
w=logspace(-1,2)
要对计算点数进行人工设定,则采用公式(63)。例如,要在与之间产生100个对数等分点,可输入如下命令:
w=logspace(0,3,100)
在画伯德图时,利用以上各式产生的频率向量w,能够很方便地画出但愿频率的伯德图。
因为伯德图是半对数坐标图且幅频图和相频图要同时在一个绘图窗口中绘制,所以,要用到半对数坐标绘图函数和子图命令。
利用工做空间中的向量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 轴均取对数变换。
MATLAB容许将一个图形窗口分红多个子窗口,分别显示多个图形,这就要用到subplot()函数,其调用格式为:
subplot(m,n,k)
该函数将把一个图形窗口分割成m×n个子绘图区域,m为行数,n为列数,用户能够经过参数k调用各子绘图区域进行操做,子图区域编号为按行从左至右编号。对一个子图进行的图形设置不会影响到其它子图,并且容许各子图具备不一样的坐标系。例如,subplot(4,3,6)则表示将窗口分割成4×3个部分。在第6部分上绘制图形。 MATLAB最多容许9×9的分割。
例21 给定单位负反馈系统的开环传递函数为:
试画出伯德图。
利用如下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个子图,均是自动完成的。
图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(?Frequency(rad/s)?)
ylabel(?Gain(dB)?)
% 第二个图画伯德图相频部分。
subplot(2,1,2);
semilogx(w,phase, ?-r?);
grid
xlabel(?Frequency(rad/s)?)
ylabel(′Phase(deg) ′)
图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 已知三阶系统开环传递函数为:
利用下面的MATLAB程序,画出系统的奈氏图,求出相应的幅值裕量和相位裕量,并求出闭环单位阶跃响应曲线。
>> 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,故奈氏曲线与实轴的交点离临界点(-1,j0)很近,且相位裕量也只有7.1578o,因此系统尽管稳定,但其性能不会太好。观察闭环阶跃响应图,能够看到波形有较强的振荡。
若是系统的相角裕量γ>45o,咱们通常称该系统有较好的相角裕量。
例23 考虑一个新的系统模型,开环传递函数为:
由下面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′)
图23 较理想的系统响应
能够看出,该系统有无穷大幅值裕量,且相角裕量高达85.4365o。因此系统的闭环响应是较理想的,闭环响应图如图23.
(1) 时间延迟系统的传递函数模型
带有延迟环节e-Ts的系统不具备有理函数的标准形式,在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]决定)
一样,要想对G句柄所对应模型的延迟时间'Td'进行修改,则可调用下面命令
Set(G, ′Td′,T)
其中T为延迟时间。由此修改后,模型G就已具备时间延迟特性。
(2) 时间延迟系统的频域响应
含有一个延迟环节的系统,其开环频域响应为
可见,该系统的幅频特性不变,只加大了相位滞后。
例24 考虑系统的开环模型为:
当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 所示,试设计一个串联校订装置,使系统知足幅值裕量大于10分贝,相位裕量≥45o
图26 校订前系统
解:为了知足上述要求,咱们试探地采用超前校订装置Gc(s),使系统变为图27的结构。
图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
能够看出,这个系统有无穷大的幅值裕量,而且其相位裕量=28o,幅值穿越频率Wcp=47rad/sec。
引入一个串联超前校订装置:
咱们能够经过下面的MATLAB语句得出校订先后系统的Bode图如图28,校订先后系统的阶跃响应图如图29。其中、、ts1分别为校订前系统的幅值穿越频率、相角裕量、调节时间,二、、ts2分别为校订后系统的幅值穿越频率、相角裕量、调节时间。
>> G1=tf(100,[0.04,1,0]); % 校订前模型
G2=tf(100*[0.025,1],conv([0.04,1,0],[0.01,1])) % 校订后模型
% 画伯德图,校订前用实线,校订后用短划线。
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 校订先后系统的阶跃响应图
能够看出,在这样的控制器下,校订后系统的相位裕量由增长到48o,调节时间由0.28s减小到0.08s。系统的性能有了明显的提升,知足了设计要求。
5 自动控制理论模拟实验
《自动控制理论》是一门理论性和实践性很强的专业基础课,前面咱们经过计算机仿真,能够方便地研究系统性能,验证理论的正确性,加深对理论知识的理解。本节咱们再经过电子模拟实验,学习和掌握系统模拟电路的构成和测试技术,进一步培养学生的实际动手能力和分析、研究问题的能力。
在控制理论课程中,大部分院校目前拥有的实验设备是电子模拟学习机。这种专为教学实验制造的电子模拟学习机,体积较小,使用方便,实验箱中备有多个运算放大器构成的独立单元,再加上经常使用的电阻、电容等器件,经过手工连线、能够构成多种特性的被控对象和控制器。
在基础训练阶段,实验手段采用模拟方法,除了灵活方便以外,还具备如下两个优势:
固然,对于未来从事实际工做的学生来讲,仅仅掌握模拟实验方法仍是不够的,应在此基础上进行一些以实际系统为主要设备的实验训练。
以自控理论电子模拟学习机为核心的一组基本实验设备和仪器,共同完成对各类实验对象的模拟和测试任务,传统的测试手段下,构成基本实验必备仪器有如下几种:
图30 传统仪器组合
按照被测系统的数学模型,在电子模拟学习机上用基本运放单元模拟出相应的电路模型,而后按图30所示的方法进行模拟实验测试。
随着计算机软、硬件的快速发展,人们愈来愈多地利用计算机实现的虚拟仪器代替传统仪器。目前,大多数实验室都是用计算机来实现信号的产生、测量与显示、系统的控制及数据处理,使实验过程更加方便,功能更强大。如今的模拟实验组件是按图31来实现的。
图31 计算机仿真模拟实验
A/D、D/A卡起模拟信号与数字信号的转换做用,还可产生不一样的输入信号(阶跃、三角、正弦等),供实验时选用。使用时用RS232串口电缆将A/D、D/A卡与计算机链接起来。若是配备打印机,则可在实验的同时将实验结果打印输出。因为计算机能够方便地输入数据、观察数据,初学者能够在屏幕的提示下进行实验过程,使学习变的更加轻松。
实验一. 典型环节及阶跃响应测试
控制系统的模拟实验是采用复合网络法来模拟各类典型环节,即利用运算放大器和RC组成的不一样输入网络和反馈网络组合,模拟出各类典型环节,而后按照给定系统的结构图将这些模拟环节链接起来,便获得了相应的模拟系统。而后将输入信号加到模拟系统的输入端,使系统产生动态响应。这时,可利用计算机或示波器等测试仪器,测量系统的输出,即可观测到系统的动态响应过程,并进行性能指标的测量。若改变系统的某一参数,还可进一步分析研究参数对系统性能的影响。
在如下的实验过程当中,为了更好地检验实验结果,避免过多地出现错误操做,咱们将每一环节的正确结果,经过Simulink仿真软件绘出正确的图形,以便于读者检验实验结果的正确性。
(1)最大超调量
利用示波器或计算机显示器上测到的输出波形,读出响应最大值和稳态值所具备的刻度值,代入下式算出超调量:
(60)
(2)峰值时间
根据示波器或显示器上输出的波形最大值,找出这一点在水平方向上所具备的刻度值,便可换算出或读出峰值时间。
(3)调节时间
一样,读出水平方向上对应输出从零到进入5%或2%偏差带时所占的刻度值,便可获得调节时间。
(1)比例环节
模拟线路如图32
图32 比例放大电路
,
通常可经过改变电阻来调整放大倍数。
因为输入信号是从运算放大器的反相输入端输入,因此输出信号和输入信号在相位上正好相反,传递函数中出现负号。为了观测方便,能够从输入端输入负阶跃信号。也能够在输出端链接一个反相器,如图33。
图33
取,,将模拟学习机上手动阶跃信号(或信号发生器置于"手动阶跃")引入环节输入端,观测输出波形,并做记录。(为便于比较,应将输入信号与输出信号同时送入双踪示波器或计算机,两路信号同时在一个坐标系下显示。绘制曲线时,也用这种形式)。
(a)仿真模块 (b) 仿真输出
图34 比例环节的Simulink仿真
图34为Simulink 的仿真模块。为便于观察,阶跃信号输入时间设置为1s(系统默认值),后面的各个例题也都适当调整输入时间。增益(Gain)模块的增益放大倍数设为2。另外,也能够用鼠标双击各模块,设置适合其它参数。
图35 积分电路
(2)积分环节
模拟线路如图35 所示。
,
积分时间常数可经过改变电阻或电容来选择。
取,,按上述一样方法观测阶跃响应波形。用Simulink仿真的环节模块图如图36 (a)。因为积分环节附带的增益比较大(积分时间常数T=0.1),Scope(示波器)绘出图形的辐值显示范围并非很理想。咱们能够在Scope的显示图中点击鼠标右键,选Axes properties菜单,在弹出的对话框中设置Y-max属性为100,则输出结果如图36 (b)所示。
(a)仿真模块 (b) 仿真输出
图36 积分环节的Simulink仿真
图37 微分电路
(3)微分环节
模拟线路图如图37 。
其中: ,
微分时间常数可经过改变和来选取。令,,,按上述一样步骤进行模拟和测试,观察微分环节的阶跃响应波形。用Simulink仿真的模块图为图38(a),在Scope绘出的图形中调整横纵坐标,得出的时域响应图如图38(b)所示。
(4)惯性环节
(a)仿真模块 (b)仿真输出
图38 微分环节的Simulink仿真
模拟线路如图39。
图39 惯性环节电路图
,,
取,,观测其阶跃响应输出,测出,并与理论值(或4T)相比较,用Simulink仿真结果如图40。
(a)仿真模块 (b)仿真输出
图40 惯性环节的Simulink仿真
(5)振荡环节
其模拟线路如图41 。
图41 振荡环节
对应的结构图如图42。
图42 振荡环节结构图
,,
改变可改变的大小,改变可改变的大小。按表2给出的参数测量阶跃响应,并记入表中。
表2 不一样的和所对应的时域性能指标
记 录 参 数 |
|
(ms) |
(ms) |
阶跃响应波形 |
|
( ) |
|
100% |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
( ) |
|
|
|
|
|
|
|
|
|
|
当,时用Simulink仿真结果如图43。
(a)仿真模块 (b)仿真输出
图43 振荡环节的Simulink 仿真
按实验给出的响应曲线,求出,,,与理论计算值比较,并得出由实验结果求惯性环节传递函数的方法。
按实验给出的欠阻尼下的响应曲线,求出和,与理论值相比较,并肯定参数,,最终可得出由实验结果求振荡环节传递函数的方法。
讨论振荡环节性能指标与,的关系。
实验二 系统频率特性测量
利用简单仪器测量频率特性,测量精度是较差的,但物理意义明显,波形直观是其特色。本实验经过"李沙育图形"法进行频率特性测试,可使学生经过实验观测到物理系统的频率响应,并根据测量值算出频率特性的幅值和相角,经过实验能够掌握测试频率特性的基本原理和方法。
一个稳定的线性系统,在正弦信号的做用下,它的稳态输出将是一个与输入信号同频率的正弦信号,但其振幅和相位却随输入信号的频率不一样而变化,测取不一样频率下系统的输出、输入信号的振幅比及相位差,便可求得这个系统的幅频特性和相频特性。
设线性系统输入和稳态响应分别为如下两个正弦信号:
幅频特性 (69)
相频特性 (70)
若以为横轴,以为纵轴,而以做为参变量,则随的变化,和所肯定的点的轨迹,将在x-y平面上描绘出一条封闭的曲线(一般是一个椭圆)。这就是所谓的"李沙育图形",如图44。
图44 "李沙育图形"原理图
不断的改变的频率,就能够得到一系列形状不一样的李沙育图形。由此求出各个频率所对应的相位差和幅值比,就可得到系统的频率特性。
幅值比由测量数据按式(69)直接求出;而相位差的具体求法以下:
令=0,则
即得 (71)
显然上式仅当时成立,"李沙育图形"在四个象限的形状如图45 所示,注意箭头方向。
实际的控制系统通常为相位滞后系统,即频率特性的相频是负的角度,相频特性滞后角按"李沙育图形"法,应按下式肯定:
第四象限: (72)
(逆时针)
第三象限: (73)
(逆时针)
第二象限: (74)
(顺时针)
第一象限: (75)
(顺时针)
图45 "李沙育图形"形状
(1)给出三阶系统模拟电路如图46 所示
图46 三阶系统模拟图
对应的系统结构图如图47 所示。
图47 系统结构图
选元件,,,
则开环传递函数为:
图48 测试电路
(2)断开闭环系统模拟电路图46中主反馈线路,按开环三阶系统在学习机上接好线路,并将有关测试仪器按图48链接。
将超低频正弦输入信号输入系统,调节输入信号幅度使被测对象在避免饱和的状况下,输出幅度尽量大,以便于测量。而后调节示波器Y轴增益(量程范围),使在所取信号幅度下,图像达到满刻度。
(3)在示波器上测量此时输入信号幅值(用表示),并记录在表3中,此后在输出幅度能有效测出时,通常再也不改变输入信号的幅度。
按表中给定的测点依次改变输入信号频率,测试并记录于表3中。
为了提升读数精度,对示波器的X,Y轴增益可随时调节,以得到较好的"李沙育图形"。注意在X轴与Y轴增益不一致时,"李沙育图形"的形状可能会有所变化。读数后按相应的增益正确折算出,值。另外,在转折频率附近以及穿越频率附近应多测几点。
表3 频率特性测试结果记录
|
8 |
10 |
15 |
20 |
30 |
40 |
60 |
80 |
100 |
200 |
f(Hz) |
1.27 |
1.6 |
2.4 |
3.2 |
4.8 |
6.4 |
9.6 |
12.7 |
16 |
32 |
2 |
|
|
|
|
|
|
|
|
|
|
2 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
李沙育图形 |
|
|
|
|
|
|
|
|
|
|
3. 实验报告
附:(1)用如下MATLAB命令绘制的开环系统伯德图见图49。
>> G=tf(10,[conv([0,1,1],[0.001,0.1,1])]);
Bode(G)
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表示系统输入在各个时刻的值。
若是想求取和,可输入以下程序
>> i=length(u);
while(u(i)>0&u(i)<0)
i=i-1;
end
y0=abs(y(i))
ym=max(y)
对频率为8时给定的u和y,可得以下结果:
y0 = ym =
8.9149 9.2467
(a)=8时的"李沙育图形" (b)=60时的"李沙育图形"
图50 MATLAB仿真"李沙育图形"
图51(a) Simulink仿真图
图51(b)响应波形
图51(c)Simulink绘出的"李沙育图形"
(3)应用Simulink仿真工具输出的响应曲线及画出的"李沙育图形"如图51所示。
在图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).
实验三 连续系统的频率法串联校订
频率法串联校订,主要是根据工程上提出的频率指标, ,等,在频率特性曲线上进行计算测量,最终得出欲增长的校订装置的功能,而后验证校订后的实际结果,通过反复调整,直至满意为止。
本实验的校订原理已在前面章节介绍,理论计算过程能够由读者本身完成。本节只经过对校订先后系统性能的测量,来进一步体会串联校订的实际效果。
(1)加深理解串联校订装置对系统动态性能的校订做用。
(2)对给定系统进行串联校订设计,并经过模拟实验检验设计的正确性。
对于给定的单位反馈闭环系统如图52 ,其对应的开环传递函数为:
图52 原闭环系统结构图
其中
系统的模拟电路图如图53 所示,图中开关S断开时对应校订前的系统,当S接通时串入超前校订环节。
图53 模拟电路图
取原系统电路元件数值为:,,,, ,有,,,
则当S断开时,原系统开环传递函数为:
由下列MATLAB命令给出校订前系统的频域指标。
>> G=tf (120, [0.03,1,0]) ;
[Gm,Pm,Wcg,Wcp]=margin(G)
结果显示
Gm = Pm = Wcg = Wcp =
Inf 29.4646 Inf 59.0015
由显示结果得,前校订系统的频域指标为:=,=59.0rad/sec。
图54 串联超前校订系统结构图
当S开关闭合时,串入有源超前校订装置,对应结构图如图54。
其中:,,
选取,
由如下MATLAB命令绘出校订先后开环系统Bode图如图55,在串联超前校订做用下,将超前角叠加在原系统附近,使校订后的系统=,=71.8 rad/sec 。
>> G=tf (120,[0.03,1,0]); %校订前模型。
Gc=tf ([0.02,1], [0.01,1]);
G_0=Gc*G; %校订后模型。
bode(G) %校订前伯德图采用默认画法。
grid
hold
bode(G_0, ':' ) %校订后伯德图用虚线画。
[Gm,Pm,Wcg,Wcp]=margin(G_0)
结果显示
Gm = Pm = Wcg = Wcp =
Inf 44.3717 Inf 71.8021
可见,系统的相角裕量由增长到,幅值穿越频率由59.0 rad/sec增长到71.8rad/sec,系统的性能有了明显的提升。
图55 串联超前校订先后系统的伯德图
(2)串联滞后校订
图56 串联滞后校订模拟电路图
仍对上述图53系统进行串联滞后校订,校订前原系统不变。当S开关闭合后,串入滞后校订装置Gc(s),对应模拟电路图如图56所示。
其中:,,
取 ,,
利用下列MATLAB命令将系统校订先后的Bode图绘于图57 。
>> G=tf (120, [0.03, 1, 0]); %校订前模型。
Gc=tf ([0.75, 1], [3, 1]);
G_0=Gc*G; %校订后模型。
bode (G) %校订前伯德图采用默认画法。
grid
hold
bode(G_0, ':' ) %校订后伯德图用虚线画。
[Gm,Pm,Wcg,Wcp]=margin(G_0)
结果显示
Gm = Pm = Wcg = Wcp =
Inf 51.5713 Inf 24.2824
可见,系统的相角裕度增长到,幅值穿越频率却减少到=24.3 rad/dec,这种技术致使系统响应速度有所下降。
图57 串联滞后校订先后系统的伯德图
(a) 仿真模块 (b) 阶跃响应曲线
图58 校订前实验系统
附: 采用Simulink仿真工具得出的校订先后的仿真模块图及输出波形分别如图5八、5九、60所示。
(a) 仿真模块 (b) 阶跃响应曲线
图59 超前校订实验图
(a) 仿真模块 (b)阶跃响应曲线
图60 滞后校订实验图
表4 串联校订测试记录
|
校订前 |
超前校订后 |
滞后校订后 |
阶跃响应曲线 |
|
|
|
|
|
|
|
|
|
|
|
|