基于MATLAB的控制工程

MATLABstandsfor“MatrixLaboratory”

MATLAB功能:

数学计算数据分析&可视化软件开发仿真开始界面Command-Window(命令窗口)

Workspace&Directory(工作空间&当前路径窗口)Command-History(命令历史窗)基本Matlab命令1、窗口命令close：该命令用于关闭当前的图形窗口。closeall：该命令用于关闭所有的图形窗口。clc：该命令用于清除命令窗口中的内容，光标将回到窗口的左上角。holdon:不刷新当前窗口。holdoff:刷新当前窗口。2、工作空间管理命令who:该命令用于列出当前工作区间的所有变量。whos：该命令用于列出当前工作区间的所有变量，并显示变量的大小、类型及其占用的磁盘空间。clear：该命令用于从工作区间清除所有的变量。注意：1、MATLAB中所有命令及表达式必须在英文状态下输入（汉字除外），而且MATLAB严格区分字母的大小写。2、所有命令都可通过help来显示该命令的帮助信息，如helpsin(显示正弦函数sin的帮助信息)。3、所有命令都必须以小写字母形式输入才能正确执行，否则出错。变量1、变量的命名变量名必须以字母开头（不能超过19个字符），之后可以是任意字母、数字或下划线；变量名称区分字母的大小写；变量中不能包含有标点符号。变量的表达格式如下：变量=表达式

在命令窗口状态下，输入的命令均存在于workspace工作空间中，只要输入who或者whos即可看到工作空间中的变量。如果要清除变量，可用clear命令。MATLAB所定义的特殊变量及其意义

）基本虚数单位（即函数的输出变量个数nargout函数的输入变量个数nargin非数（NotANumber）,如0/0NaN或nan∞值，无限大，如1/0inf圆周率π值3.14159265...piMATLAB定义的正的极小值=2.2204e-16eps列出所有定义过的变量名称whoi或j用于结果的缺省变量名,最近的计算结果的变量名(即answer的缩写)ans意义变量名基本算术运算符符号意义符号意义符号意义+相加^矩阵乘方\左除-相减.^数组乘方/右除*矩阵相乘./数组右除.*数组点乘.\数组左除关系运算符关系运算符用于比较数、矩阵、字符串之间的大小或不等关系，其返回值为0或1。<：小于>：大于<=：小于等于>=：大于等于==：等于~=：不等于逻辑运算符基本语句（1）if…else…end语句语法： if条件式1

语句段1 elseif条件式2

语句段2 ... else

语句段n+1 end基本语句（2）for...end循环结构语法：

for循环变量=array

循环体 end说明：循环体被循环执行，执行的次数就是array的列数，array可以是向量也可以是矩阵，循环变量依次取array的各列，每取一次循环体执行一次。基本语句（3）while...end循环结构语法： while表达式

循环体 end说明：只要表达式为逻辑真，就执行循环体；一旦表达式为假，就结束循环。表达式可以是向量也可以是矩阵，如果表达式为矩阵则当所有的元素都为真才执行循环体，如果表达式为nan，MATLAB认为是假，不执行循环体。基本语句（4）switch…case开关结构语法： switch开关表达式 case表达式1

语句段1 case表达式2

语句段2 ... otherwise

语句段n end说明：(1)将开关表达式依次与case后面的表达式进行比较，如果表达式1不满足，则与下一个表达式2比较，如果都不满足则执行otherwise后面的语句段n；一旦开关表达式与某个表达式相等，则执行其后面的语句段。(2)开关表达式只能是标量或字符串。(3)case后面的表达式可以是标量、字符串或元胞数组，如果是元胞数组则将开关表达式与元胞数组的所有元素进行比较，只要某个元素与开关表达式相等，就执行其后的语句段。基本语句（5）流程转向控制语句break命令break命令可以使包含break的最内层的for或while语句强制终止，立即跳出该结构，执行end后面的命令，break命令一般和If结构结合使用。continue命令continue命令用于结束本次for或while循环，只结束本次循环而继续进行下次循环。向量的建立线性等间距生成法A=[start：step：end]start为起始值,step为步长,end为终止值。当步长为1时可省略step参数。step取正数，起始值start<终止值end；step取负数，起始值start>终止值end。》a=[1:2:10]a=13579定数线性采样法a=linspace(n1,n2,n)在线性空间上，行矢量的值从n1到n2，数据个数为n，缺省n为100。》a=linspace(1,10,10)a=12345678910定数对数采样法a=logspace(n1,n2,n)

，在对数空间上，行矢量的值从10

n1

到10n2，数据个数为n，缺省n为50。这个指令为建立对数频域轴坐标提供了方便。》a=logspace(1,3,3)a=101001000矩阵运算1、矩阵的要素整个矩阵必须用方括号“[]”括起来；矩阵A(m,n)表示矩阵变量A是一个m行n列的矩阵。矩阵行与行之间必须用分号“；”或在输入完一行之后直接用回车键【Enter】隔离；矩阵元素之间必须以逗号“，”或空格号分离。矩阵在MATLAB中是按先列后行的方式储存。绘制平面曲线(绘制单条曲线)plot(x,y)：以x为横坐标、y为纵坐标绘制二维图形x,y是同维数的向量；plot(y)：相当于x=[1,2,…,length(y)]时情形。绘制平面曲线（绘制多条曲线)1.plot(x,[y1;y2;…]),x是横坐标向量，[y1;y2;…]是由若干函数的纵坐标拼成的矩阵2.plot(x,y1),holdon,plot(x,y2),holdoff3.plot(x,y1,x,y2,…)4.plotyy两个坐标系，用于绘制不同尺度的函数。xlabel(‘’):x轴标注，ylabel(‘‘):y轴标注title(‘‘):图标名subplot(n,m,p):窗口分成n行m列块，在第p块绘制图形ezplot(fun)ezplot(fun,[min,max])eaplot是一个易用的一元函数绘图函数。特别是在绘制含有符号变量的函数的图像时，ezplot要比plot更方便。因为plot绘制图形时要指定自变量的范围，而ezplot无需数据准备，直接绘出图形常用波形函数cos():余弦波sin():正弦波square():方波square(,duty):特殊周期的方波sawtooth():锯齿波sawtooth(t,0.5):三角波

MATLAB的内部函数是有限的，有时为了研究某一个函数的各种性态，需要为MATLAB定义新函数，为此必须编写函数文件.函数文件是文件名后缀为M的文件，这类文件的第一行必须是一特殊字符function开始，格式为：

function因变量名=函数名（自变量名）函数值的获得必须通过具体的运算实现，并赋给因变量.

4、M文件M文件建立方法：1.在Matlab中，点:File->New->M-file2.在编辑窗口中输入程序内容

3.点：File->Save，存盘，M文件名必须与函数名一致。Matlab的应用程序也以M文件保存。M文件：纯文本文件，扩展名为M，用任何字处理软件都可以编辑修改。注意：M文件的文件名必须是字母开头，不能用系统已有函数名。练习

先建立M文件fun1.m定义函数f（x），再在Matlab命令窗口输入fun1(2),fun1(-1)即可。矩阵运算1、矩阵的要素整个矩阵必须用方括号“[]”括起来；矩阵A(m,n)表示矩阵变量A是一个m行n列的矩阵。矩阵行与行之间必须用分号“；”或在输入完一行之后直接用回车键【Enter】隔离；矩阵元素之间必须以逗号“，”或空格号分离。矩阵在MATLAB中是按先列后行的方式储存。一些常用的特殊矩阵单位矩阵：eye(m,n)；eye(m)零矩阵：zeros(m,n)；zeros(m)纯1矩阵：ones(m,n)；ones(m)对角矩阵：对角元素向量V=[a1,a2,…,an]，A=diag(V)随机矩阵：rand(m,n)产生一个m×n的均匀分别的随机矩阵矩阵的修改X=A;Then,todeletethesecondcolumnofX,useX(:,2)=[]ThischangesXtoX=16213511897124141A=16321351011896712415141X=A;X=[A;[1,3,5,7]]ThischangesXtoX=16221351111897712414411357A=16321351011896712415141增加行取子阵>>A([1,3,4],[1,3])ThischangesAtoA=16297414A=16321351011896712415141矩阵运算求和-----sum(A)转置-------‘对角矩阵-----diag(A)矩阵求逆-----inv(A)行列式的值——det(A)矩阵的特征值——eig(A)PolynomialFunction

FunctionDescriptionconv：Multiplypolynomials.（乘）deconv：Dividepolynomials.（除）poly：Polynomialwithspecifiedroots（求特殊多项式系数）polyder：Polynomialderivative.（求导）polyfit：Polynomialcurvefitting.（曲线拟合）polyval：Polynomialevaluation.（求值）polyvalm：Matrixpolynomialevaluation.（矩阵求值）residue：Partial-fractionexpansion(residues).（部分分式展开）roots：Findpolynomialroots.（求特征多项式根）RepresentingPolynomialsP=[1，0，-2，-5];PolynomialRootsPartialFractionExpansionb=[-4,8]；a=[1,6,8]；[r,p]=residue(b,a)MATLAB在机械工程控制中的应用控制系统模型控制系统的时域分析控制系统的频域分析系统稳定性分析控制系统校正及PID控制器仿真MATLAB在机械工程控制中的应用控制系统模型控制系统的时域分析控制系统的频域分析系统稳定性分析控制系统校正及PID控制器仿真

控制系统模型传递函数模型状态方程模型零极点增益模型

对线性定常系统，式中s的系数均为常数，且a1不等于零，

分子系数:num=[b1,b2,…,bm,bm+1]

分母系数:den=[a1,a2,…,an,an+1]

注意：按s降幂进行排列。连续系统的传递函数模型连续系统的传递函数如下：用函数tf(num,den)就可得到系统的传递函数。

G=tf(num,den)在MATLAB中，用[z,p,K]矢量组表示零极点增益模型。

z=[z1,z2,…,zm]，p=[p1,p2,...,pn]，K=[k]. G=zpk(z,p,k)零极点增益模型K为系统增益，zi为零点，pj为极点函数tf2zp()可以用来求传递函数的零极点和增益。[z,p,k]=tf2zp(num,den)模型之间的转换

模型转换的函数包括：tf2zp：传递函数模型转换为零极点增益模型.zp2tf：零极点增益模型转换为传递函数模型.

同一系统可用不同的模型表示，它们各有特点，相互之间可以进行转换。传递函数模型零极点增益模型zp2tftf2zp例子>>a1=[521];b1=[47132];>>G=tf(a1,b1)Transferfunction:5s^2+2s+1-----------------------------4s^4+7s^3+s^2+3s+2当传递函数的分子和分母由若干个多项式的乘积表示时，它可由MATLAB提供的多项式乘法运算函数conv()来处理，以获得分子多项式式向量。此函数的调用格式为:

c=conv(a，b)>>num=4*conv([12]，[166])；>>den=conv([10]，conv([11]，conv([11]，[1325])))；>>G1=tf(num，den)系统串联sys1sys2系统的串联连接Uy[num,den]=series(num1,den1,num2,den2)并联系统并联连接sys1sys2UU1YY1[num,den]=parallel(num1,den1,num2,den2)反馈Sys1Sys2Y2+UY系统的反馈连接[num,den]=feedback(num1,den1,num2,den2，sign)5400.1++——G3G1G2G4>>num1=[10];den1=[11];num2=[10];den2=[20.5];>>num3=[540];den3=[1];num4=[0.1];den4=[1];>>[na,da]=series(num1,den1,num2,den2);>>[nb,db]=feedback(na,da,num4,den4,-1);>>[nc,dc]=series(num3,den3,nb,db);>>[num,den]=cloop(nc,dc,-1);>>printsys(num,den)

回车，结果为：num/den= 54000--------------------------------2s^2+2.5s+5401.5MATLAB在机械工程控制中的应用控制系统模型控制系统的时域分析控制系统的频域分析系统稳定性分析控制系统校正及PID控制器仿真

控制系统的时域分析求取系统单位阶跃响应：step()求取系统的脉冲响应：impulse()step(num,den)：可得到由传递函数表示的系统的单位阶跃响应。时间向量t由系统模型的特性自动生成。

step(num,den,t)，可利用用户指定的时间矢量t来绘制响应曲线。t为指定的仿真时间向量，一般可以由t=0:step:end等步长地产生出来。

求取系统单位阶跃响应：step()[y,x,t]=step(num,den)；[y,x,t]=step(num,den,t)；[y,t]=step(G,t)；

求取脉冲（激励）响应的调用方法系统的单位脉冲响应：impulse()

[y,x,t]=impulse(num,den)；[y,x,t]=impulse(num,den,t)；格式：impulse函数用于计算系统的单位脉冲（冲激）响应，当不带输出变量时，可在当前图形窗口中绘出系统的响应曲线。典型二阶系统求当=8,=0.2时,在时间向量t=[0:0.01:3]下的单位阶跃响应。wn=8;%频率kosi=0.2;%阻尼比num=[wn^2];%传递函数分子den=[1,2*kosi*wn,wn^2];%传递函数分母t=[0:0.01:3];%时间向量G=tf(num,den);%得到传递函数[y,t]=step(G,t)%也可以用step(num,den,t),得到响应曲线例子典型二阶系统求当=8,=0.2时,在时间向量t=[0:0.01:3]下的单位脉冲响应。wn=8;%频率kosi=0.2;%阻尼比num=[wn^2];%传递函数分子den=[12*kosi*wnwn^2];%传递函数分母t=[0:0.01:3];%时间向量G=tf(num,den);%得到传递函数impulse(num,den,t)%得到响应曲线例子其他常用时域分析函数

initial：连续系统的零输入响应lsim：连续系统对任意输入的响应格式：[y,x,t]=initial(num,den,x0)[y,x,t]=lsim(num,den,u,t)x0为初始状态,t为时间向量，u为输入clearall;t=[0:0.01:0.8];nG=[50];tao=0;dG=[0.051+50*tao50];G1=tf(nG,dG);tao=0.0125;dG=[0.051+50*tao50];G2=tf(nG,dG);tao=0.025;dG=[0.051+50*tao50];G3=tf(nG,dG)[y1,T]=impulse(G1,t);[y1a,T]=step(G1,t);[y2,T]=impulse(G2,t);[y2a,T]=step(G2,t);[y3,T]=impulse(G3,t);[y3a,T]=step(G3,t);subplot(121),plot(T,y1,'--',T,y2,'-.',T,y3,'-')legend('tao=0','tao=0.0125','tao=0.025')xlabel('t(sec)'),ylabel('x(t)');gridon;subplot(122),plot(T,y1a,'--',T,y2a,'-.',T,y3a,'-')legend('tao=0','tao=0.0125','tao=0.025')gridon;xlabel('t(sec)'),ylabel('x(t)');例3-8-1MATLAB在机械工程控制中的应用控制系统模型控制系统的时域分析控制系统的频域分析系统稳定性分析控制系统校正及PID控制器仿真

控制系统的动态特性分析

在MATLAB中,可以先利用step求出系统的单位阶跃响应,再通过编程求出系统的动态特性.num=[5];%分子den=[125];%分母wendingzhi=polyval(num,0)/polyval(den,0);%求出稳态值[y,x,t]=step(num,den);%阶跃响应L=length(t);%响应时间长度[Y,K]=max(y);%峰值Y以及Y所在的序号fengzhishijian=t(K);%求峰值时间chaotiaoliang=100*(Y-wendingzhi)/wendingzhi;%求超调量求系统的调节时间、稳态值、峰值、峰值时间、超调量、上升时间。%求上升时间fori1=1:Lify(i1)<0.1*wendingzhit1=t(i1);elsebreakendendforj1=1:Lify(j1)<1*wendingzhit2=t(j1);elsebreakendend%上升时间shangshengshijian=t2-t1；%求调节时间while(y(L)>0.98*wendingzhi)&(y(L)<1.02*wendingzhi)L=L-1;endtiaojieshijian=t(L);disp('调节时间'),tiaojieshijiandisp(‘稳态值'),wendingzhidisp('峰值'),Ydisp('峰值时间'),fengzhishijiandisp('超调量'),chaotiaoliangdisp('上升时间'),shangshengshijian调节时间:tiaojieshijian=3.6994稳定值:wendingzhi=1峰值:Y=1.2076峰值时间:fengzhishijian=1.5460超调量:chaotiaoliang=20.7555上升时间:shangshengshijian=0.8282例3-8-3closeall;%t=0:0.001:1;%yss=1;dta=0.02;%nG=[50];tao=0;dG=[0.051+50*tao50];G1=tf(nG,dG);tao=0.0125;dG=[0.051+50*tao50];G2=tf(nG,dG);tao=0.025;dG=[0.051+50*tao50];G3=tf(nG,dG);y1=step(G1,t);y2=step(G2,t);y3=step(G3,t);%r=1;whiley1(r)<yss;r=r+1;endtr1=(r-1)*0.001;%[ymax,tp]=max(y1);tp1=(tp-1)*0.001;%

mp1=(ymax-yss)/yss;%s=1001;whiley1(s)>1-dta&y1(s)<1+dta;s=s-1;endts1=(s-1)*0.001;%r=1;whiley2(r)<yss;r=r+1;endtr2=(r-1)*0.001;[ymax,tp]=max(y2);tp2=(tp-1)*0.001;mp2=(ymax-yss)/yss;s=1001;whiley2(s)>1-dta&y3(s)<1+dta;s=s-1;endts2=(s-1)*0.001;%r=1;whiley3(r)<yss;r=r+1;endtr3=(r-1)*0.001;[ymax,tp]=max(y3);tp3=(tp-1)*0.001;mp3=(ymax-yss)/yss;s=1001;whiley3(s)>1-dta&y3(s)<1+dta;s=s-1;endts3=(s-1)*0.001%[tr1tp1mp1ts1;tr2tp2mp2ts2;tr3tp3mp3ts3]利用MATLAB进行频率特性分析求取系统对数频率特性:bode()[mag,phase,w]=bode(num,den)[mag,phase,w]=bode(num,den,w)

格式：幅值相角频率求取系统对数频率特性:bode()bode(num,den)：可绘制出以连续时间多项式传递函数表示的系统的波特图。

bode(num,den,w)：利用指定的角频率矢量w绘制出系统的波特图。当输出变量[mag,pha,w]或[mag,pha]引用函数时，可得到系统波特图相应的幅值mag、相角pha及角频率点矢量w或只是返回幅值与相角。相角以度为单位，幅值可转换为分贝单位：magdb=20×log10(mag)

二阶系统，已知=2，阻尼比=0.7。试求系统的传递函数，并绘制系统的Bode图。[num,den]=ord2(2,0.7);G=tf(num,den);bode(G);gridon;%给Bode图加上网格

先利用ord2函数求出系统的传递函数模型，再利用bode函数绘制系统的Bode图。例子

求取系统奈奎斯特图：nyquist()[re,im,w]=nyquist(num,den)[re,im,w]=nyquist(num,den,w)

格式：nyquist(num,den)可绘制出以连续时间多项式传递函数表示的系统的极坐标图。nyquist(num,den,w)：可利用指定的角频率矢量w绘制出系统的极坐标图。当不带返回参数时，直接在屏幕上绘制出系统的极坐标图（图上用箭头表示w的变化方向，负无穷到正无穷）。当带输出变量[re,im,w]引用函数时，可得到系统频率特性函数的实部re和虚部im及角频率点w矢量（为正的部分）。可以用plot(re,im)绘制出对应w从负无穷到零变化的部分。nyquist(num,den)%注意不能直接用nyquist(G)，这样会造成输入向量不够，gridon求系统的nyquist图.num=[203];den=[5231];G=tf(num,den)例子s作的系统的nyquist和Bode图%建模k=24;nunG1=k*[0.251];denG1=conv([52],[0.052]);%求实频和虚频函数[re,im]=nyquist(nunG1,denG1);%画nyquist图plot(re,im);gridonpause%画Bode图w=logspace(-2,3,100);bode(nunG1,denG1,w);gridon判断系统的稳定性已知系统的开环传递函数为绘制系统nyquist图并判断闭环系统稳定性先用nyquist函数绘制出系统的nyquist图;再利用nyquist判据判断闭环系统的稳定性.思路：例子从上图可看出nyquist曲线逆时针包围(-1,j0)点1圈，而从开环系统的零极点图可知开环系统在右平面有1极点，所以闭环系统稳定。求相对稳定性第一步:先求出系统传递函数G=tf(num,den)第二步:利用margin函数求出幅值裕量,相位裕度以及相应的频率.[Gm,Pm,Wcp,Wcg]=margin(G)第三步:将幅值裕量转化为幅值裕度GmdB=20*log10(Gm)利用MATLAB直接求相对稳定性例子1:num=conv([2.5],[0.11]);%求取传函分子2.5*(0.1*S+1)的系数den=conv([100],[0.21]);%求取传函分母的系数G=tf(num,den)%求取传函[G