计算机过程控制实验指导书

实验指导书工程学院与技术系实验1控制系统动态特性测试在设定值和扰动信号的作用下，过程控制系统的输出可由两条通道来产生：值附近。干扰通道：干扰信号对被控变量影响的通道。能的影响。(1)增益对控制系统的影响假设被控对象的传递函数为G0(s)=2e-10s/(35s+1),考核系统在单位阶跃信号作用下，运行下列matlab语句，观察响应曲线，得出有关结论。G0=tf(2,[351]);%2/35s+1[np,dp]=pade(10,2);Gp=tf(np,dp);%e-10sG1=G0*Gp;Kc=1:0.5:2.5;%Kc=11.52.02.5holdonfori=1:length(Kc)Gc=feedback(Kc(i)*G1,1);%阶跃响应%等键盘step(Gc);pause的灵敏度，有利于提高控制质量)Gses5s+1),干扰的传递函数为Gd(s)=2e-5s/(7s+1),G0=tf(1,[71]);[np,dp]=pade(5,2);Gp1=tf(np,dp);Gd=G0*Gp1;G1=tf(2,[351]);[[np,dp]=pade(10,2);Gp2=tf(np,dp);Go=G1*Gp2;Kc=1.5;G3=feedback(Go,Kc);G4=Gd*G3;holdonKd=1:4fori=1:length(Kd)G=Kd(i)*G4;tepG%等键盘pause从图形可以看出，随着扰动通道增益Kd的增加，系统的稳态误差增加，并且扰动作用小愈好)。(2)时间常数对控制系统的影响T一定的阻力而引起的，同T值(25，35，45)下系统的响应，运行下列matlab语句，观察响应曲线，得出有关结Kc=1.5;T=[253545];holdonfori=1:length(T)G0=tf(2,[T(i)1]);[np,dp]=pade(10*T(i)/35,2);Gp=tf(np,dp);G1=G0*Gp;Gc=feedback(Kc*G1,1);stepGc);pause;设T0=35，改变扰动通道的时间常数Td(10，35，50)，观察扰动阶跃响应曲线。44%研究扰动通道不同时间常数对系统的影响G0=tf(2,[351]);Kc=1.5;G1=feedback(G0,Kc);Td=[103550];holdonfori=1:length(Td)Gd=tf(5,[Td(i)1]);G=Gd*G1;epGpause系统扰动通道的时间常数Td越大，扰动对输出的影响越缓慢，有利于系统克服干d(3)时滞对控制系统的影响时滞t0指输出变量的变化落后于输入变量变化的时间。滞后的产生是由于介质的输送时，简单的控制系统已很难满足要求，要考虑负责方案。扰动通道的时滞扰动通道的td不会对系统的稳定性产生影响，仅仅表示扰动进入系统的时间先后对系有影响。实验2比例积分微分控制规律特性分析系统的余差及稳定性的影响。假设被控系统为Gp(s)=e-50s/(36s+1)rte-50s/rte-50s/(36s+1)─%%tf函数：传递函数定义参数：分子、分母G0=tf(1,[36,1]);[np,dp]=pade(50,2);G1=tf(np,dp);%GpGp=G0*G1;P=[0.5,0.7,0.9,1,1.5];oldonholdon;fori=1:length(P)Gc=feedback(P(i)*Gp,1,-1);%定义反馈结构step(Gc);%求阶跃响应pause;假设被控系统为Gp(s)=e-50s/(36s+1)，只采用积分策略，研究不同的Ki值下，%研究积分速度对系统调节的影%研究积分速度对系统调节的影响rG0=tf(1,[36,1]);[np,dp]=pade(50,2);G1=tf(np,dp);Gp=G0*G1;Ki=[0.005,0.01,0.015,0.02];holdon;fori=1:length(Ki)Gc=tf(Ki(i),[1,0])G=feedback(Gc*Gp,1,-1);tepGpause;结论：积分作用可以消除余差，但增大Ki将会降低系统的稳定性，甚至会导致系统不由于微分作用不单独采用，所以研究比例微分作用，改变微分时间常数Td，观察系统rG0=tf(1,[36,1]);[np,dp]=pade(50,2);G1=tf(np,dp);Gp=G0*G1;Kp=0.8;Td=20:5:35;holdon;fori=1:length(Td)Gc=tf(Kp*[Td(i),1],1)G=feedback(Gc*Gp,1);epGpause;1、零极点增益模型形式G(S)=k[(S-Z)(S-Z)…(S-Z)]/(S-P)(S-P)…(S-P)12m12nz,z…,z:系统零点；12,m12n统的传函模型给出以后，可以立即得出系统的零极点模型。k=[k];12mppp…;p];12n状态空间模型转换为零极点增益模型。ksys=zpk(z,p,k)G(S)=[6(S+1.9294)(S+0.0353±0.9287i)]/[(S+0.9567±1.2272i)(S-0.0433±0.6412i)]k=6;z=[-1.9294;-0.0353+0.9287*i;-0.0353-0.9287*i];p=[-0.9567+1.2272*i;-0.9567-1.2272*i;+0.0433+0.6412*i;+0.0433-0.6412*i];G=zpk(z,p,k)假设系统的模型为Gp=10/(s+1)(s+2)(s+3)(s+4),研究不同调节器下闭环系统阶跃响%%研究不同调节器下闭环系统阶跃响应。Gp=zpk([],[-1;-2;-3;-4],10);holdon;Gc=Kp;G1=feedback(Gp*Gc,1);stepG);pauseGc=tf(Ki,[1,0]);G2=feedback(Gp*Gc,1);stepG);pausePIGc=tf(Kp*[1,1/Ti],[1,0]);G3=feedback(Gp*Gc,1);stepG);pauseGc=tf(Kp*[Td,1],1);G4=feedback(Gp*Gc,1);stepG);pauseKp=7.4;Ti=1.5;Td=0.38;%PID调节Gc=tf(Kp*[Td,1,1/Ti],[1,0]);G5=feedback(Gp*Gc,1);stepG);得出各种控制策略的特点。画出各段语句的图形，得出相关结论。实验3数字PID控制算法的实现学习传递函数到离散方程的转换方法。Gp(s)=e-80s/(60s+1)，系统的采样周期为20s，延迟时间为4个采样周期，则被控对muk%%积分分离式PID%采样时间ts;%被控对象离散化sys=tf([1],[60,1],'inputdelay',80);%描述被控对象u_1=0;u_2=0;u_3=0;u_4=0;u_5=0;y_1=0;y_2=0;y_3=0;error_1=0;error_2=0;fork=1:1:200;time(k)=k*ts;%离散化对象y%Iseparationrin(k)=40;error(k)=rin(k)-yout(k);ktseiifM==1%采用分段积分分离方式ifabs(error(k))>=30&abs(error(k))<=40;beta=0.3;elseifabs(error(k))>=20&abs(error(k))<=30;beta=0.6;elseifabs(error(k))>=10&abs(error(k))<=20;beta=0.9;beta=1.0;elseifM==2%不采用积分分离方式beta=1.0;kp=0.80;ki=0.005;kd=3.0;u(k)=kp*error(k)+kd*(error(k)-error_1)/ts+beta*ki*ei;%控制器的输出限制ifuk)>=110u(k)=110;ifu(k)<=-110u(k)=-110;u_5=u_4;u_4=u_3;u_3=u_2;u_2=u_1;u_1=u(k);y_3=y_2;y_2=y_1;y_1=yout(k);error_2=error_1;error_1=error(k);plot(time,rin,'b',time,yout,'r');xlabel('time(s)');ylabel('rin,yout');画出各段语句的图形，得出相关结论。实验4PID调节器参数的工程整定如下图所示。该广义对象可以用如下数学模型来近似：K=y(∞)-y(0)/△u如果在曲线上取y(t1)=0.39y(t2)=0.63T=2(t2-t1)t=2t1-t2%采用曲线拟合法，可近似得到广义对象的一阶惯性加纯滞后对象的参数%假设被控对象的传递函数为1/(s+1)(2s+1)(5s+1)(10s+1)num=1;den=conv(conv([1,1],[2,1]),conv([5,1],[10,1]));Gp=tf(num,den);%原系统figure(1),holdonpGppausek=0.993;T=14.4;L=6.6;%请从曲线用鼠标读取数组，并计算相应参数值G0=tf(k,[T1]);[np,dp]=pade(L,2);G1=tf(np,dp);G=G0*G1;%近似模型pauseKp1=T/(1.1*k*L);Ti1=3.3*L;Gc1=tf(Kp1*[1,1/Ti1],[1,0]);G_c1=feedback(G*Gc1,1);figure(2),step(G_c1)观察曲线形状，求取有关数值，了解整定公式画出各段语句的图形，描述整定过程。MATLAB是Mathworks公司开发的一种集数值计算、符号计算和图形可视化三大基本功能于一体的功能强大、操作简单的优秀工程计算应用软件。MATLAB不仅可以处理代数问题师及科学家解决实际的技术问题。MATLAB的含义是矩阵实验室(MatrixLaboratory)，最初主要用于方便矩阵的存取，工具箱(Toolbox)的综合应用软件，不仅成为线性代数课程的标准工具，而且适合具有不同专业研究方向及工程应用需求的用户使用。数(称为M文件)，从而构成适合于其它领域的工具箱，大大扩展了MATLAB的应用范围。目前，MATLAB已成为国际控制界最流行的软件，控制界很多学者将自己擅长的CAD方法用(1)MATLAB的安装ndowsMATLAB以执行安装盘内的setup.exe文件启动MATLAB的安装程序。启动安装程序后，屏幕将显示程如下：码；收软件公司的协议；入用户名和公司名；选择MATLAB组件(Toolbox)；件安装路径和目录；Next界面如图1所示。安装完毕后，选择[Restartmycomputernow]选项以重新启动计算机。了。MATLAB启动过程界面如(2)MATLAB桌面系统BhPad)④路径浏览器(CurrentDirectoryBrowser)⑤帮助浏览器(HelpBrowser)⑥工作空间浏览器(WorkspaceBrowser)⑦数组编辑器(ArrayEditor)⑧M文件编辑调试器(Editor-Debugger)。MATLAB桌面平台(3)MATLAB命令窗口这种方式简单易用，但在编程过程中要修改整个程序比较困难，并且用户编写的程序不F5)运行整个程序。运行过程中的错误信息和运行结果显示在命令窗口中。整个程序的源代在介绍MATLAB的强大计算和图象处理功能前，我们可以先运行一个简单的程序。设系统的闭环传递函数为:求系统的时域响应图，可输入下面的命令：>>num=[1,4];den=[1,2,8];step(num,den)程序运行后会在一个新的窗口中显示出系统的时域动态响应曲线，如图4。用鼠标左键numdenGs两个完全不同的变量。应该注意所有的函数名均由小写字母构成。MATLAB是一个功能强大的工程应用软件，它提供了相当丰富的帮助信息，同时也提供了多种获得帮助的方法。如果用户第一次使用MATLAB，则建议首先在>>提示符下键入DEMO2MATLAB基本操作命令本节简单介绍与本书内容相关的一些基本知识和操作命令。(1)简单矩阵的输入这就是说，MATLAB只有一种数据形式，那就是矩阵，或者数的矩形阵列。标量可看作为1输入矩阵最便捷的方式为直接输入矩阵的元素，其定义如下：(1)元素之间用空格或逗号间隔；(2)用中括号([])把所有元素括起来；(3)用分号(；)指定行结束。>>a=[23563649MATLAB的矩阵输入方式很灵活，大矩阵可以分成n行输入，用回车符代替分号或用续行符号(…)将元素续写到下一行。例如：a=[123456789]在MATLAB中，矩阵元素不限于常量，可以采用任意形式的表达式。同时，除了直接输入方式之外，还可以采用其它方式输入矩阵，如：(1)利用内部语句或函数产生矩阵；(3)利用外部数据文件装入到指定矩阵。(2)复数矩阵输入MATLAB允许在计算或函数中使用复数。输入复数矩阵有两种方法：(1)a=[12;34]+i*[56;78](2)a=[1+5i2+6i;3+7i4+8i](3)MATLAB语句和变量执行语句一样。变量＝表达式[;]表达式[;]表达式可以由操作符、特殊符号、函数、变量名等组成。表达式的结果为一矩阵，它赋给左边的变量，同时显示在屏幕上。如果省略变量名和“＝”号，则MATLAB自动产生一个ans是MATLAB提供的固定变量，具有特定的功能，是不能由用户清除的。常用的固定MATAB允许在函数调用时同时返回多个变量,而一个函数又可以由多种格式进行调用，[返回变量列表]＝fun-name(输入变量列表)例如用bode()函数来求取或绘制系统的Bode图，可由下面的格式调用：和以分号“；”结束的特殊效用(5)获取工作空间信息统会自动将其保存在工作空间里。在退出程序之前，这些变量将被保留在存储器中。为了得到工作空间中的变量清单，可以在命令提示符>>后输入who或whos命令，当前存放在工作空间的所有变量便会显示在屏幕上。(6)常数与算术运算符其中i=1。epsIEEE其范围为：应该注意：(∕)右除法和(﹨)左除法这两种符号对数值操作时，其结果相同，其斜(7)选择输出格式不影响结果矩阵的计算和存储。所有计算都是以双精度方式完成的。(1)如果矩阵的所有元素都是整数，则矩阵以不带小数点的格式显示。x=[101]四①formatshort(短格式,也是系统默认格式)②formatshorte(短格式科学表示)③formatlong(长格式)④formatlonge(长格式科学表示)x短格式科学表示长格式科学表示一旦调用了某种格式，则这种被选用的格式将保持，直到对格式进行了改变为止。(8)MATLAB图形窗口当调用了一个产生图形的函数时，MATLAB会自动建立一个图形窗口。这个窗口还可分mpjpg作空间的内容。(9)剪切板的使用WindowsMATLAB间交换信息。(1)要将MATLAB的图形移到其它应用程序，首先按Alt-PrintScreen键，将图形复制到剪切板中，然后激活其它应用程序，选择edit(编辑)中的paste(粘贴)，就可以在应用程序中得到MATLAB中的图形。当然还可以借助于copytoBitmap或copytoMetafile选项来传递图形信息。(2)要将其它应用程序中的数据传递到MATLAB，应先将数据放入剪切板，然后在然后选择Edit中的paste，最后加上“]”,这样可将应用程序中的数据送入MATLAB(10)MATLAB编程指南M文件是用MATLAB语言编写的可在MATLAB环境中运行的磁盘文件。它为脚本文件(ScriptFile)和函数文件(FunctionFile),这两种文件的扩展名都是.m。(1)脚本文件是将一组相关命令编辑在一个文件中，也称命令文件。脚本文件的语句可以访问MATLAB工作空间中的所有数据，运行过程中产生的所有变量都是全局变量。例如下述语句如果以.m为扩展名存盘，就构成了M脚本文件,我们不妨将其文件名取为“Step_Response”。%用于求取一阶跃响应。num=[14];den=[128];step(num,den)用于求取一阶跃响应。容易地对想查询的各个函数的有关信息进行查询。该命令使用格式为：help命令或函数名注意：若用户把文件存放在自己的工作目录上，在运行之前应该使该目录处在MATLAB(2)函数文件是用于定义专用函数的，文件的第一行是以function作为关键字引导的，后面为注释和函数体语句。函数就像一个黑箱，把一些数据送进去，经加工处理，再把结果送出来。在函数体内使都是局部变量，执行完后，这些内部变量就被清除了。函数文件的文件名与函数名相同(文件名后缀为.m),它的执行与命令文件不同，不能例1求一系列数的平均数，该函数的文件名为“mean.m”functiony=mean(x)%这是一个用于求平均数的函数w=length(x);y=sum(x)/w;length函数表示取向量x的长度n量，程序执行完后，便不存在了。运行结果显示3MATLAB在控制系统中的应用MATLAB是国际控制界目前使用最广的工具软件，几乎所有的控制理论与应用分支中都SystemsToolbox)和仿真环境(Simulink)，学习MATLAB的应用。(1)用MATLAB建立传递函数模型线性系统的传递函数模型可一般地表示为：Gsmmnm(1)Gsmmnm(1)nn将系统的分子和分母多项式的系数按降幂的方式以向量的形式输入给两个变量num和 量构造出单个的传递函数对象。从而使得系统模型的输入和处理更加方便。G＝tf(num，den);(4)模型:>>num=[1，5];den=[1，2，3，4，5];G=tf(num，den)Transferfunction:s+5-----------------------------s^4+2s^3+3s^2+4s+5G(s)=6(s+5)(s2+3s+1)2(s+6)>>num=6*[1，5];den=conv(conv([1，3，1]，[1，3，1])，[1，6]);tf(num,den)运行结果Transferfunction:6s+30-----------------------------------------s^5+12s^4+47s^3+72s^2+37s+6其中conv()函数(标准的MATLAB函数)用来计算两个向量的卷积，多项式乘法当然也可以用这个函数来计算。该函数允许任意地多层嵌套，从而表示复杂的计算。2.零极点模型线性系统的传递函数还可以写成极点的形式：G(s)=K(5)KGainK(6)12m12n零极点对象，用于简单地表述零极点模型。该函数的调用格式为：G=zpkG=zpk(Z,P,KGain)G(s)=6(s+1.9294)(s+0.0353士0.9287j)(s+0.9567士1.2272j)(s0.0433士0.6412j)>>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.反馈系统结构图模型控控制系统工具箱中提供G=feedback(G1,G2,sign);取反馈连接下总的系统模(10)G(s)=，1(s+1)2G(s)=s>>G1=tf(1,[1,2,1]);G2=tf(1,[1,1]);G=feedback(G1,G2)运行结果:Transferfunction:s+1---------------------s^3+3s^2+3s+2若采用正反馈连接结构输入命令>>G=feedback(G1,G2,1)则得出如下结果:Transferfunction:s+1-----------------s^3+3s^2+3s>>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)Transferfunction:0.1s^5+10.75s^4+77.75s^3+278.6s^2+361.2s+120--------------------------------------------------------------------0.01s^6+1.1s^5+20.35s^4+110.5s^3+325.2s^2+384s+12044.有理分式模型与零极点模型的转换有了传递函数的G1=zpk(G)(11)对应的零极点格式可由下面的命令得出>>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)G1=tf(G)(12)ssj(s+1.5)((1)SIMULINK的启动：在MATLAB可以用下面的MATLAB命令立即得出其等效的传递函数模型。输入程序的过程中要注意>>Z=[-2,-7];P=[0,-3-2j,-3+2j,-1.5];K=6.8;G=zpk(Z,P,K);G1=tf(G)Transferfunction:6.8s^2+61.2s+95.2-------------------------------s^4+7.5s^3+22s^2+19.5s功能：Simul(仿真)与Link(连接)，亦即可以利用鼠标在模型窗口上“画”出所需的控制系输入命令相比，这样输入更容易，分析更直观。下面简单介绍SIMULINK建立系统模型的基本步骤：命令窗口的工具栏中单击按钮命令窗口的工具栏中单击按钮或者在命令提示符>>所示。这一模型库中含有许多子模型库，如Sources(输入源模块库)、(2)画出系统的各个模块：打开相应的子模块库，选择所需要的元素，用鼠标左键点中后拖图7simulink模型库图8模型编辑窗口到模型编辑窗口的合适位置。型为1/(s+1)的简单格式，必须通过修改得到实际的模块参数。要修改模块的参数，可对话框，提示用户修改模块参数。的连线。若需要从连线中引出节点，可在鼠标点击起始节点时按住Ctrl键，再将鼠标接模块库)图标，从中选取相应的输入输出端子。(1)画出所需模块，并给出正确的参数：左键双击该图标，打开参数设定的对话框，将参数steptime(阶跃时刻)设为0。在Math(数学)子模块库中选中加法器(sum)图标，拖到编辑窗口中，并双击该图igns在continuous(连续)子模块库中、选积分器(Integrator)和传递函数(TransferFcn)图标拖到编辑窗口中，并将传递函数分子(Numerator)改为〔900〕，分母 sinks(输出)子模块库中选择scope(示波器)和Out1(输出端口模块)图标并将之(3)将画出的所有模块按图9用鼠标连接起来，构成一个原系统的框图描述如图10所示。(4)选择仿真算法和仿真控制参数，启动仿真过程。在solver模板中设置响应的仿真范围StartTime(开始时间)和StopTime(终止时间)，仿真步长范围Maxinumstepsize(最大步长)和Mininumstepsize(最小步长)。对于上会“实时地”显示出仿真结果。输出结果如图11B(2)利用MATLAB进行时域分析1.线性系统稳定性分析线性系统稳定的充要条件是系统的特征根均位于S平面的左半部分。系统的零极点模B也可以用于系统的分析和计算。(1)直接求特征多项式的根vvrootsp(13)项式为：x5+3x3+2x2+x+1>>p=[1,0,3,2,1,1];v=roots(p)v=v0.3202+1.7042i0.3202-1.7042i-0.72090.0402+0.6780i0.0402-0.6780i利用多项式求根函数roots(),可以很方便的求出系统的零点和极点，然后根据零极点分析系统稳定性和其它性能。(2)由根创建多项式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=p1.0000-0.00003.00002.00001.00001.0000(3)多项式求值polyval(ppolyval(p,v)>>p=[1,0,3,2,1,1];polyval(p,x)结果显示ans=8==01nsdenasnasna1n式中a0，但是a和b中某些量可能为零。0ij(16)M(s)=r(1)+r(2)+...+r(n)+k(s)N(s)sp(1)sp(2)sp(n)12nr(1)=r,r(2)=r,…,r(n)=r为各极点的留数，k(s)为余项。12n该传递函数的部分分式展开由以下命令获得：>>num=[2,5,3,6];den=[1,6,11,6];r,p,k]=residue(num,den)命令窗口中显示如下结果-6.0000-4.0000000-3.0000-2.0000-1.00002G(s)=6+4+3+2numdenresiduerpk(17)>>[num,den]=residue(r,p,k)结果显示num=2.00005.00003.00006.00001.00006.000011.00006.0000开式将包括下列诸项： >>v=[-1,-1,-1]num=[0,1,2,3];den=poly(v);[r,p,k]=residue(num,den)结果显示0.00002.0000-1.0000-1.0000-1.0000polydenk=[]为空矩阵。G(s)=+++0(5)由传递函数求零点和极点。Z=tzero(G)(18)P=G.P{1}(19){}表示的矩阵元素，详细内容参阅后面章节。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=Z-7-2P=P0-3.0000+2.0000i-3.0000-2.0000i-1.5000(6)零极点分布图。在MATLAB中，可利用pzmap()函数绘制连续系统的零、极点图，从而分析系统的稳pzmap(num,den)(20)例14给定传递函数:>>num=[3,2,5,4,6];den=[1,3,4,2,7,2];pzmap(num,den)title(¹Pole-ZeroMap¹)%图形标题。动态特性分析。(1)时域响应解析算法――部分分式展开法。用拉氏变换法求系统的单位阶跃响应，可直接得出输出c(t)随时间t变化的规律，对于高阶系统，输出的拉氏变换象函数为：(21)对函数c(s)进行部分分式展开，我们可以用num,[den,0]来表示c(s)的分子和分母。递函数：>>num=[1,7,24,24]den=[1,10,35,50,24][r,p,k]=residue(num,[den,0])输出结果为r=p=k=-1.0000.0000-1.0000-1.0000-4.0000-3.0000-2.0000-1.00000C(s)=1+211+1+0(2)单位阶跃响应的求法：控制系统工具箱中给出了一个函数step()来直接求取线性系统的阶跃响应，如果已知式：stepnumden(22)stepnumdent(23)或stepG(24)stepGt(25)该函数将绘制出系统在单位阶跃输入条件下的动态响应图，同时给出稳态值。对于式统根据输出曲线的形状自行设定。cstepG(26)ot个或多个给定的向量绘制二维图形，详细介绍可以查阅后面的章节。>>num=[0,0,25];den=[1,4,25];step(num,den)grid绘制网格线。title(¹Unit-StepResponseofG(s)=25/(s^2+4s+25)¹)%图像标题我们还可以用下面的语句来得出阶跃响应曲线>>G=tf([0,0,25],[1,4,25]);c=step(G,t);%动态响应的幅值赋给变量cCss=dcgain(G)%求取稳态值。系统显示的图形类似于上一个例子，在命令窗口中显示了如下结果Css=1(3)求阶跃响应的性能指标MATLAB提供了强大的绘图计算功能，可以用多种方法求取系统的动态响应指标。我们首先介绍一种最简单的方法――游动鼠标法。对于例16，在程序运行完毕后，用鼠标左键间)和纵坐标(幅值)。按住鼠标左键在曲线上移动，可以找到曲线幅值最大的一点――即简单易用，但同时应注意它不适用于用plot()命令画出的图形。另一种比较常用的方法就是用编程方式求取时域响应的各项性能指标。与上一段介绍的游动鼠标法相比，编程方法稍微复杂，但通过下面的学习，读者可以掌握一定的编程技巧，通过前面的学习，我们已经可以用阶跃响应函数step()获得系统输出量，若将输出量ytstepG(27)到时域性能指标。①峰值时间(timetopeak)可由以下命令获得：Ykmaxy(28)timetopeak=t(k)(29)②最大(百分比)超调量(percentovershoot)可由以下命令得到：C=dcgain(G);Ykmaxy(30)percentovershoot=100*(Y-C)/C(31)dcgain()函数用于求取系统的终值，将终值赋给变量C，然后依据超调量的定义，一下循环语句while的使用。while式为：while<循环判断语句>循环体其中,循环判断语句为某种形式的逻辑判断表达式。要求出上升时间，可以用while语句编写以下程序得到：C=dcgain(G);n=1;whiley(n)<Cn=n+1;risetime=t(n)应的时刻，即为上升时间。则计算程序如下：C=dcgain(G);n=1;whiley(n)<0.1*Cn=n+1;m=1;whiley(n)<0.9*Cm=m+1;risetime=t(m)-t(n)④调节时间(setllingtime)可由while语句编程得到：C=dcgain(G);i=length(t);while(y(i)>0.98*C)&(y(i)<1.02*C)ii1;setllingtime=t(i)例17已知二阶系统传递函数为:Gssm>>G=zpk([],[-1+3*i,-1-3*i],3);%计算最大峰值时间和它对应的超调量。C=dcgain(G)plot(t,y)timetopeak=t(k)percentovershoot=100*(Y-C)/C升时间。n=1;whiley(n)<Cn=n+1;risetime=t(n)%计算稳态响应时间。ilengtht);while(y(i)>0.98*C)&(y(i)<1.02*C)setllingtime=t(i)运行后的响应图如图15,命令窗口中显示的结果为C=timetopeak=0.30001.0491percentovershoot=35.0914setllingtime=5337risetime=.6626图15二阶系统阶跃响应有兴趣的读者可以用本节介绍的游动鼠标法求取此二阶系统的各项性能指标。将它们与假设闭环系统中的开环传递函数可以表示为：G(s)=Ksm+b1sm1+...+bm1+bm=Knum=K(s+z1)(s+z2)...(s+zm)=KG(s)ksn+asn1+...+as+aden(s+p)(s+p)...(s+p)01n1n12n特征方程的根随参数K的变化而变化，即为闭环根轨迹。控制系统工具箱中提供了rlocus()函数，可以用来绘制给定系统的根轨迹，它的调用格式有以下几种特征方程的根随参数K的变化而变化，即为闭环根轨迹。控制系统工具箱中提供了rlocus()函数，可以用来绘制给定系统的根轨迹，它的调用格式有以下几种：rlocus(num，den)rlocus(num，den，K)rlocus(G)rlocus(G，K)(32)(33)以上给定命令可以在屏幕上画出根轨迹图，其中G为开环系统G(s)的对象模型,K为用0函数调用中需要返回参数，则调用格式将引入左端变量。如locusG此时屏幕上不显示图形，而生成变量R和K。行。若需要画出根轨迹，则需要采用以下命令：plot(R，¹¹)(37)plot()函数里引号内的部分用于选择所绘制曲线的类型，详细内容见表1。控制系统工具箱中还有一个rlocfind()函数，该函数允许用户求取根轨迹上指定点处的开环增益值，并将该增益下所有的闭环极点显示出来。这个函数的调用格式为：cfindG所有的闭环极点直接在根轨迹曲线上显示出来。：G(s)=K=KG(s)ks(s+1)(s+2)0>>G=tf(1,[conv([1,1],[1,2]),0]);rlocus(G);title(¹Root_LocusPlotofG(s)=K/[s(s+1)(s+2)]¹)xlabel(¹RealAxis¹)%给图形中的横坐标命名。ylabel(¹ImagAxis¹)%给图形中的纵坐标命名。用鼠标点击根轨迹上与虚轴相交的点，在命令窗口中可发现如下结果select_point=0.0000+1.3921iK.8142-2.29830-0.0085+1.3961i-0.0085-1.3961i化，观察对应根处阶跃响应的变化。考虑K=0.1,0.2,…,1,2,…,5,这些增益下闭环系统的>>holdoff;%擦掉图形窗口中原有的曲线。t=0:0.2:15;Y=[];forK=[0.1:0.1:1,2:5]GK=feedback(K*G,1);y=step(GK,t);Y=[Y,y];plot(t,Y)方式都是相似的，它们只是在如何给坐标轴进行分度和如何显示数据上有所差别。(1)二维图形绘制plotxy(39)例19如果想绘制出一个周期内的正弦曲线，则首先应该用t=0:0.01:2*pi(pi是系统自定义的常数,可用help命令显示其定义)命令来产生自变量t；然后由命令y=sin(t)(2)一幅图上画多条曲线。利用具有多个输入变量的plot()命令，可以在一个绘图窗口上同时绘制多条曲线，YnxYxYxy都可以用图解表示出来，因而可以在一幅图上画出多条曲线。多重变量的优点是它允许不同长度的向量在同一幅图上显示出来。每一对向量采用不同的线型以示区别。(3)图形的线型和颜色为了区分多幅图形的重叠表示，MATLAB提供了一些绘图选项，可以用不同的线型或颜色来区分多条曲线，常用选项见下表1表1MATLAB绘图命令的多种选项选项选项选项短划线线用星号绘制各个数据点色用圆圈绘制各个数据点用圆点绘制各个数据点黄色用叉号绘制各个数据点表1中绘出的各个选项有一些可以并列使用，能够对一条曲线的线型和颜色同时作出规定。例如′--g′表示绿色的短划线。带有选项的曲线绘制命令的调用格式为：plot(X1，Y1,S1,X2,Y2,S2,…)(41)grid络线)(42)title(′图形标题′)(43)函数引号内的字符串将被写到图形的坐标轴上或标题位置。(5)在图形屏幕上书写文本。如果想在如果想在图形窗口中书写文字，可以单击按钮，选择屏幕上一点，点击鼠标，在光标处输入文字。另一种输入文字的方法是用text()命令。它可以在屏幕上以(x,y)为坐标将从点(3，0.45)开始，水平的写出“sint”。(6)自动绘图算法及手工坐标轴定标在MATLAB图形窗口中，图形的横、纵坐标是自动标定的，在另一幅图形画出之前，这幅图形作为现行图将保持不变，但是在另一幅图形画出后，原图形将被删除，坐标轴自动地能需要放弃绘图命令中的坐标轴自动标定特性，由用户自己设定坐标范围，可以在程序中加入下列语句：vxminx-maxy-miny-max]()47axisv(48)式中v是一个四元向量。axis(v)把坐标轴定标建立在规定的范围内。对于对数坐标图，v的元素应为最小值和最大值的常用对数。统的自动标定特性。Axis(′sguare′)能够把图形的范围设定在方形范围内。对于方形长宽比，其斜率为1复到正常状态。5、线性系统的频域分析。(1)频率特性函数G(j)。nnG(jw)=a(j)n+a(j)n+...+a(j)+an1nisqrtGW=polyval(num，i*w)./polyval(den，i*w)(49)(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)该命令将画出下列开环系统传递函数的奈氏曲线：G(s)=num(s)den(s)计算这些点的频率响应。reimwnyquistG(55)或reimwnyquistGw(56)函数运行后不在屏幕上产生图形，而是将计算结果返回到矩阵re、im和w中。矩阵reim别表示频率响应的实部和虚部，它们都是由向量w中指定的频率点计算得到的。>>num=[0,0,1];den=[1,0.8,1];nyquist(num,den)%设置坐标显示范围v=[-2,2,-2,2];axis(v)title(′NyquistPlotofG(s)=1/(s^2+0.8s+1)′)(3)用MATLAB作伯德图控制系统工具箱里提供的bode()函数可以直接求取、绘制给定线性系统的伯德图。e上不显示频率响应图。命令的调用格式为：[mag,phase,w]=bode(num,den)[mag,phase,w]=bode(num,den,w)或magphasewbodeGmagphasew=bode(G,w)(59)矩阵mag、phase包含系统频率响应的幅值和相角，这些幅值和相角是在用户指定的分贝值，因为分贝是作幅频图时常用单位。可以由以下命令把幅值转变成分贝：magdb=20*log10(mag)(61)绘图时的横坐标是以对数分度的。为了指定频率的范围，可采用以下命令格式：或(62)(63)121212102空间中，即w=logspace(-1，2)要对计算点数进行人工设定，则采用公式(63)。例如，要在=1与=1000之间2w=logspace(0，3，100)用到半对数坐标绘图函数和子图命令。(1)对数坐标绘图函数semilogxxys(64)semilogyxys(65)(2)子图命令subplot(m，n，k)G(s)=10(s+1)s(s+7)>>num=10*[1,1];den=[1,7,0];bode(num,den)title(′BodeDiagramofG(s)=10*(s+1)/[s(s+7)]′)如果希望显示的频率范围窄一点，则程序修改为：>>num=10*[1,1];den=[1,7,0];[mag,phase,w]=bode(num,den,w);magdB=20*log10(mag)%增益值转化为分贝值。%第一个图画伯德图幅频部分。subplot(2,1,1);semilogx(w,magdB,′-r′)%用红线画title(′BodeDiagramofG(s)=10*(s+1)/[s(s+7)]′)xlabel(¹Frequency(rad/s)¹)ylabel(¹Gain(dB)¹)%第二个图画伯德图相频部分。subplot(2,1,2);semilogx(w,phase,¹-r¹);xlabel(¹Frequency(rad/s)¹)ylabel(′Phase(deg)′)采用游动鼠标法求取系统的幅值裕量和相位裕量。例如，我们可以在图20的幅频曲线上按住鼠标左键游动鼠标，找出纵坐标(Magnitude)趋近于零的点，从提示框图中读出其频率约为7.25dB。然后在相频曲线上用同样的方法找到横坐标(Frequence)最接近7.25dB的此外，控制系统工具箱中提供了margin()函数来求取给定线性系统幅值裕量和相位裕，该函数可以由下面格式来调用：[Gm,Pm,Wcg,Wcp]=margi