第7章 控制系统的MATLAB仿真

1、 本章主要教学内容在MATLAB中描述控制系统的数学模型系统方框图模型(mxng)的简化及应用控制系统的时域分析控制系统的频域分析利用MATLAB语言进行仿真编程的具体应用第7章控制系统(kn zh x tn)的MATLAB仿真1共一百零二页本章教学目的(md)及要求熟悉MATLAB的基本应用掌握利用MATLAB建立数学模型的方法熟悉控制系统的时域和频域分析掌握MATLAB的仿真编程应用第7章控制系统(kn zh x tn)的MATLAB仿真2共一百零二页 7.1 控制系统的模型表示 MATLAB 提供了数学模型的建立函数和各模型之间的转换功能函数，可以分别采用传递函数、零极点增益、状态空间以

2、及动态结构图等4种数学模型来表示控制系统，前3种是用数学表达式描述的系统模型，每种模型都有连续(linx)系统及离散系统两种类别的表示；而动态结构图是基于传递函数的图形化形式，是采用MATLAB中提供的SIMULINK结构图来实现的。MATLAB中使用的数学模型之间的转换也很方便，使得采用MATLAB编制的程序更加简单、精炼而高效。第7章控制系统(kn zh x tn)的MATLAB仿真3共一百零二页 7.1.1 系统的传递函数模型表示 传递函数模型(mxng)通常表示线性定常时不变系统（LTI），可以是连续的时间系统，也可以是离散的时间系统。 对于离散的时间系统，其脉冲传递函数可表示为：第7

3、章控制系统(kn zh x tn)的MATLAB仿真 不论是连续的还是离散的时间系统，其传递函数的分子/分母多项式均按s或z的降幂来排列。在MATLAB中可直接采用分子/分母多项式系数构成的两个向量num与den来表示系统，即：4共一百零二页 在MATLAB中，可用函数命令tf()来建立控制系统的传递函数模型，其调用格式和功能(gngnng)分别为：（1）sys= tf(num, den)；（2）sys= tf(num, den, Ts) ；。（3）sys= tf(M) ；（4）tfsys= tf(sys)； 第7章控制系统(kn zh x tn)的MATLAB仿真5共一百零二页第7章7.1.

4、2 零极点增益模型 当连续系统的传递函数表达式采用系统增益、系统零点与系统极点来表示时，称之为系统零极点增益模型。系统零极点增益模型是传递函数模型的一种特殊(tsh)形式。 离散系统的传递函数零极点增益模型： 控制系统(kn zh x tn)的MATLAB仿真6共一百零二页第7章 在MATLAB里，连续与离散系统都可直接用向量z、p、k构成的矢量(shling)组z, p, k来表示系统，即： 控制系统(kn zh x tn)的MATLAB仿真7共一百零二页第7章 在MATLAB中，可用函数命令zpk()来建立控制系统的零极点增益模型，其调用格式(g shi)和功能分别为：（1）sys= zp

5、k(num, den) ；（2）sys= zpk(num,den, Ts)； （3）sys= zpk(M)； （4）tfsys= zpk(sys)；控制系统(kn zh x tn)的MATLAB仿真8共一百零二页 7.1.3 状态空间模型 控制系统在主要工作区域内的一定条件下可近似为线性时不变（LTI）模型，连续LTI对象系统总是能用一阶微分方程组来表示(biosh)，写成矩阵形式即为状态空间模型：第7章控制系统(kn zh x tn)的MATLAB仿真其中： 式（a）系统的状态方程，是由n个一阶微分方程组成的微分方程组；式（b）系统的输出方程，是由1个线性代数方程组成的；9共一百零二页第7章

6、离散系统的状态(zhungti)空间模型可表示为：控制系统(kn zh x tn)的MATLAB仿真式中：U系统的控制输入向量；X系统的状态向量；Y系统的输出向量；k特定时刻的采样点；A状态矩阵，由控制对象的参数决定；B控制矩阵；C输出矩阵；D直接传输矩阵。10共一百零二页 MATLAB中的函数ss()可用来建立控制系统(xtng)的状态空间模型，或者将传递函数模型与零极点增益模型转换为系统(xtng)状态空间模型。ss()函数的调用格式为：（1）sys= ss(a, b, c, d)；（2）sys= ss(a, b, c, d, Ts)；（3）sys= ss(d)；该函数等价于sys= ss

7、( , , ,d)。（4）sys_ss= ss(sys)；第7章控制系统(kn zh x tn)的MATLAB仿真11共一百零二页第7章【例7.1】已知某系统的状态(zhungti)空间表达式为：控制系统(kn zh x tn)的MATLAB仿真试采用MATLAB语言求出该系统的状态空间模型。12共一百零二页解：采用状态空间模型表示时，可在MATLAB命令窗口中输入以下命令： A=1 0 0 0;0 1 0 0;0 0 1 0;-1 -5 0 -2; B=0;0;0;1; C=3 2 1 0;D=0; sys=ss(A,B,C,D)上述指令(zhlng)执行后可得指定系统的状态空间模型为：a

8、= x1 x2 x3 x4 x1 1 0 0 0 x2 0 1 0 0 x3 0 0 1 0 x4 -1 -5 0 -2第7章控制系统(kn zh x tn)的MATLAB仿真13共一百零二页第7章 b = u1 x1 0 x2 0 x3 0 x4 1 c = x1 x2 x3 x4 y1 3 2 1 0 d = u1 y1 0 Continuous-time model. 控制系统(kn zh x tn)的MATLAB仿真14共一百零二页 7.1.4 系统不同模型(mxng)间的相互转换 在MATLAB 6.1的信号处理工具箱与控制系统工具箱中，提供了传递函数模型、零极点增益模型与状态空间模

9、型之间转换的函数：ss2tf()，ss2zp()，tf2ss()，tf2zp()，zp2ss()，zp2tf()。 这些函数之间的转换功能如表7-1所示。 第7章控制系统(kn zh x tn)的MATLAB仿真15共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真表7-1 数学模型之间的转换函数(hnsh)及其功能函 数 名函 数 功 能ss2tf将系统状态空间模型转换为传递函数模型ss2zp将系统状态空间模型转换为零极点增益模型tf2ss将系统传递函数模型转换为状态空间模型tf2zp将系统传递函数模型转换为零极点增益模型zp2ss将系统零极点增益模型转换为状态空间模型zp

10、2tf将系统零极点增益模型转换为传递函数模型16共一百零二页【例7.2】 已知某系统(xtng)的传递函数为：第7章控制系统(kn zh x tn)的MATLAB仿真 试用MATLAB语言求出该系统的传递函数模型、状态空间模型和零极点增益模型。解：（1）求系统的传递函数模型在MATLAB命令窗口输入以下命令：num=12 24 12 20;den=2 4 6 2 2;sys=tf(num, den)17共一百零二页第7章执行以上语句后可得系统的传递函数模型为：12 s3 + 24 s2 + 12 s + 20-2 s4 + 4 s3 + 6 s2 + 2 s + 2（2）求系统的状态空间模型该

11、系统的状态空间模型可以通过MATLAB的模型转换函数来完成。在MATLAB命令窗口输入(shr)以下命令：a,b,c,d=tf2ss(num,den);sys=ss(a,b,c,d) 控制系统(kn zh x tn)的MATLAB仿真18共一百零二页执行完上述语句后，可得系统的状态空间模型的状态矩阵a，系统控制矩阵b，系统输出(shch)矩阵c，系统直接传输矩阵d分别为：a = x1 x2 x3 x4 x1 -2 -3 -1 -1 x2 1 0 0 0 x3 0 1 0 0 x4 0 0 1 0第7章控制系统(kn zh x tn)的MATLAB仿真19共一百零二页第7章 b = u1 x1

12、1 x2 0 x3 0 x4 0 c = x1 x2 x3 x4 y1 6 12 6 10 d = u1 y1 0 Continuous-time model.控制系统(kn zh x tn)的MATLAB仿真20共一百零二页由以上(yshng)数据可写出系统的状态空间模型为：第7章控制系统(kn zh x tn)的MATLAB仿真21共一百零二页第7章（3）求系统的零极点增益模型该系统的零极点增益模型也可以通过MATLAB的模型转换函数来完成。在MATLAB命令窗口(chungku)输入以下命令： z, p, k=tf2zp(num, den); sys= zpk(z, p, k)执行以上语

13、句后可得系统的零极点增益模型为：Zero/pole/gain: 6 (s+1.929) (s2 + 0.07058s + 0.8638)-(s2 + 0.08663s + 0.413) (s2 + 1.913 s + 2.421)控制系统(kn zh x tn)的MATLAB仿真22共一百零二页7.2 环节方框图模型的化简 7.2.1 环节串联连接的化简 多个环节串联的连接形式是控制系统(kn zh x tn)最基本的组成结构形式之一。控制系统(kn zh x tn)的环节串联及其化简就是模块方框图模型的串联及其化简。可以用MATLAB的函数命令series()将串联模块进行等效变换。 使用s

14、eries()函数命令不必做多项式的乘除运算即可实现两个环节传递函数的串联连接。如果令sys1= tf(num1,den1)，sys2= tf(num2,den2)，其命令格式为：sys= series(sysl,sys2)第7章控制系统(kn zh x tn)的MATLAB仿真23共一百零二页如果已知两个环节的状态空间模型矩阵组分别为：（a1,b1,c1,d1）与（a2,b2,c2,d2），则求两个环节串联连接等效系统状态空间模型a,b,c,d矩阵组的命令格式为：a,b,c,d= series(a1,bl,cl,dl,a2,b2,c2,d2) series()函数命令还可以将多个环节按两两串

15、联的形式多次递归调用加以连接，进行等效化简。 sys= series(sysl,sys2)命令可以用命令sys=sys1*sys2*sysn取代，不仅省掉“series()”字符，且可以实现(shxin)多个环节的串联等效传递函数的求取。第7章控制系统(kn zh x tn)的MATLAB仿真24共一百零二页第7章【例7.3】 已知双闭环调速系统电流(dinli)环内的前向通道3个模块传递函数分别为：控制系统(kn zh x tn)的MATLAB仿真试求串联连接的等效传递函数及其等效状态空间模型。25共一百零二页解：（1）根据(gnj)MATLAB程序设计的基本方法和函数命令series()，

16、可以编写出MATLAB程序如下：n1=0.0128 1;d1=0.04 0;sys1=tf(n1,d1);n2=30;d2=0.00167 1;sys2=tf(n2,d2);n3=2.5;d3=0.0128 1;sys3=tf(n3,d3);sys=sys1*sys2*sys3s1=ss(sys1);s2=ss(sys2);s3=ss(sys3);sys12=series(s1,s2);sys123=series(sys12,s3)第7章控制系统(kn zh x tn)的MATLAB仿真26共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真（2）在MATLAB命令窗口输入程序

17、名，程序运行后得到如下(rxi)电流环内前向通道的等效传递函数及等效状态空间模型：Transfer function: 0.96 s + 75-8.55e-007 s3 + 0.0005788 s2 + 0.04 s a = x1 x2 x3 x1 -78.13 2246 0 x2 0 -598.8 800 x3 0 0 0 27共一百零二页 b = u1 x1 0 x2 40.96 x3 4 c = x1 x2 x3 y1 12.21 0 0 d = u1 y1 0Continuous-time model.第7章控制系统(kn zh x tn)的MATLAB仿真28共一百零二页第7章 7.

18、2.2 环节并联连结的化简 环节并联是指多个(du )环节的输入信号相同，所有环节输出的代数和为其总输出。 采用parallel()函数命令可以等效化简两个环节的并联连接。parallel()函数命令调用格式为：num，den= parallel（numl,denl,num2,den2）该命令已由命令sys=sysl+sys2+sysn所取代，省掉了“parallel()”字符，且可以实现多个环节的并联等效处理。parallel()函数命令调用格式还有：a,b,c,d= parallel(a1,bl,cl,dl,a2,b2,c2,d2) 控制系统(kn zh x tn)的MATLAB仿真29共

19、一百零二页第7章【例7.4】 已知两个环节的传递函数分别为： 试求两环节并联连接等效传递函数的num与den向量(xingling)及等效的状态空间模型。控制系统(kn zh x tn)的MATLAB仿真30共一百零二页第7章解：（1）采用函数命令parallel()，给出MATLAB程序如下：num1=1;den1=1 2;sys1=tf(num1,den1);num2=2 1;den2=1 1 2;sys2=tf(num2,den2);s1=ss(sys1);s2=ss(sys2);sys=sys1+sys2sys12=parallel(s1,s2)（2）在MATLAB命令窗口运行程序后得

20、到等效传递函数及等效状态空间(kngjin)模型为：Transfer function: 3 s2 + 6 s + 4-s3 + 3 s2 + 4 s + 4 控制系统(kn zh x tn)的MATLAB仿真31共一百零二页第7章a = x1 x2 x3 x1 -2 0 0 x2 0 -1 -1 x3 0 2 0b = u1 x1 1 x2 2 x3 0c = x1 x2 x3 y1 1 1 0.25d = u1 y1 0控制系统(kn zh x tn)的MATLAB仿真32共一百零二页 7.2.3 环节反馈连结的化简 利用MATLAB中的feedback()函数命令可将两个(lin )环节

21、按反馈形式进行连接后求其等效传递函数。feedback()函数命令形式：（1）G=feedback(G1,G2,sign)； （2）a,b,c,d=feedback(a1,b1,c1,d1,a2,b2,c2,d2,sign)（3）num,den=feedback(num1,den1,num2,den2,sign)；第7章控制系统(kn zh x tn)的MATLAB仿真33共一百零二页第7章 【例7.5】已知晶闸管-直流电机单闭环调速系统(xtng)的动态结构图如图7-4所示，求该闭环系统(xtng)的传递函数。控制系统(kn zh x tn)的MATLAB仿真图7-4 直流单闭环调速系统34

22、共一百零二页第7章解：（1）给出MATLAB程序(chngx)如下：n1=1;d1=0.017 1;s1=tf(n1,d1);n2=1;d2=0.075 0;s2=tf(n2,d2);sys1=feedback(s1*s2,1)n3=0.049 1;d3=0.088 0;s3=tf(n3,d3);n4=44;d4=0.00167 1;s4=tf(n4,d4);n5=1;d5=0.1925;s5=tf(n5,d5);n6=0.01178;d6=1;s6=tf(n6,d6);sysq=sys1*s3*s4*s5;sys=feedback(sysq,s6)控制系统(kn zh x tn)的MATLA

23、B仿真35共一百零二页第7章 （2）程序运行结果(ji gu)为：Transfer function: 1-0.001275 s2 + 0.075 s + 1 Transfer function: 2.156 s + 44-3.607e-008 s4 + 2.372e-005 s3 + 0.001299 s2+ 0.04234 s + 0.5183 控制系统(kn zh x tn)的MATLAB仿真36共一百零二页第7章由以上运算数据(shj)可得单闭环系统的小闭环的传递函数为：单闭环系统的闭环的传递函数（略去分母的项）为：控制系统(kn zh x tn)的MATLAB仿真37共一百零二页第7

24、章7.3 控制系统时域分析的MATLAB实现7.3.1 概述1系统时域响应仿真(fn zhn)的主要问题时域分析法是根据系统微分方程用拉普拉斯变换求解系统的动态响应。典型的动态过程响应有单位阶跃响应、单位斜坡响应、单位加速度响应与单位脉冲响应等。时域分析的另外一个目的是求解响应的性能指标。通常将控制系统跟踪或复现阶跃输入信号响应的指标作为系统控制性能的指标。阶跃响应的一般性能指标有：峰值时间tp、超调量%、调节时间ts、及稳态误差ess。 控制系统(kn zh x tn)的MATLAB仿真38共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真 2时域分析的MATLAB实现 时

25、域响应的MATLAB仿真方法有两种，一种是在MATLAB的函数指令(zhlng)方式下进行时域仿真；另一种是在SIMULINK窗口菜单操作方式下进行时域仿真。 对于连续系统，可利用MATLAB提供的求取单位阶跃响应函数step()、单位脉冲响应函数impulse()、零输入响应函数initial()等处理系统的响应。 对于离散系统，可利用MATLAB提供的单位阶跃响应函数dstep()、单位脉冲响应函数dimpulse()、零输入响应函数dinitial()等来处理。 39共一百零二页第7章 7.3.2 连续系统的单位阶跃响应 利用step()函数可求出连续系统的单位阶跃响应，该命令的调用格式

26、和功能如下：（1）step(sys)； （2）step(sys，t)；（3）step(sysl,sys2,，sysn)；（4）y，t，x=step(sys)； 求出了系统的阶跃响应曲线，要分析系统的稳定性能，还需要计算系统的性能指标。包括峰值时间tp、超调量%、调节时间ts、及稳态误差ess等。这些参数可以(ky)直接从图上读取，也可通过MABLAB的M函数来求取。 控制系统(kn zh x tn)的MATLAB仿真40共一百零二页第7章下面介绍利用M函数(hnsh)来计算性能指标。定义该指标函数为terg函数，其程序语句如下：function mp,tp,ess,b1,b2,sigma,n,

27、pusi,T,f=targ(y,t) mp,tf=max(y); ct=length(t);mp; tp=t(tf);tm=max(tf);yss=(ct);q=1;m=q-1;while m3,控制系统(kn zh x tn)的MATLAB仿真41共一百零二页第7章 for a=(tm/100):0.01:tm j=0:a:tm; for i=1:length(j); if (y(i+1)-y(i)0 m=m+1; pm(m)=y(i); tp(m)=t(i); end end endendyss=y(ct);ess=1-yss;控制系统(kn zh x tn)的MATLAB仿真42共一百零

28、二页第7章 b1=pm(1)-yss;b2=pm(2)-yss;sigma=100*b1/yss;n=b1/b2;pusi=(b1-b2)/b1;T=(tp(2)-tp(1);f=1/Ttp=tp(1); 上述程序中：mp为曲线的峰值；tp为峰值时间；ess 为阶跃响应的余差；b1为第一峰值的超调量；b2为第二峰值的超调量；sigma为阶跃响应的超调量；n为阶跃响应的衰减比；pusi为阶跃响应的衰减率；T为衰减振荡(zhndng)周期；f振荡(zhndng)频率。控制系统(kn zh x tn)的MATLAB仿真43共一百零二页第7章7.3.3 连续系统的单位(dnwi)脉冲响应求连续系统单位

29、脉冲响应的函数为impulse()，该函数命令的调用格式和功能如下： （1）impulse(sys)；（2）impulse(sys，t)；（3）impulse(sysl,sys2,，sysN)； （4）y，t，x=impulse(sys)； 【例7.6】已知单位负反馈系统前向通道的传递函数为： ，试绘出系统的单位阶跃响应和脉冲响应曲线。并计算系统的阶跃响应性能指标。控制系统(kn zh x tn)的MATLAB仿真44共一百零二页第7章解：（1）求阶跃响应曲线(qxin)在命令窗口输入：num=4; den=1 1 0;sys=tf(num,den);closys=feedback(sys,1

30、);% 求闭环传递函数step(closys)可得如图7-5所示的阶跃响应曲线。控制系统(kn zh x tn)的MATLAB仿真45共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真图7-5 系统(xtng)的阶跃响应曲线46共一百零二页第7章 （2）求脉冲响应曲线在命令(mng lng)窗口继续输入：impulse(closys) ；可得如图7-6所示的脉冲响应曲线。控制系统(kn zh x tn)的MATLAB仿真图7-6 系统的脉冲响应曲线47共一百零二页第7章 （3）求系统的阶跃响应性能指标 首先在MATLAB程序编辑窗口编辑前面给出的targ函数，并保存(bocn

31、)函数名为targ.m。然后继续在MATLAB命令窗口输入以下语句：y,t=step(slosys);mp,tp,ess,b1,b2,sigma,n,pusi,T,f=targ(y,t)运行完上述语句后，可得以下性能指标：mp = 1.4441 ；tp = 1.6062ess = -0.0017 ；b1 = 0.4425b2 = 0.0857 ；sigma = 44.1728n = 5.1651 ；pusi = 0.8064T = 3.2125 ；f = 0.3113控制系统(kn zh x tn)的MATLAB仿真48共一百零二页第7章7.3.4 连续系统的零输入响应求连续系统零输入响应的函

32、数为initial()，该函数命令的格式(g shi)如下：（1）y，t，x= initial(sys)；（2）initial(sys,x0)； （3）initial(sys,x0, t)； （4）initial(sysl,sys2,，sysN,x0)。控制系统(kn zh x tn)的MATLAB仿真49共一百零二页第7章 7.4 控制系统频域分析的MATLAB实现7.4.1 概述1控制系统频域响应仿真的主要问题 对于频域分析，MATLAB仿真就是用先进的计算机仿真技术来精确绘制Bode图、Nyquist曲线图、Nichols曲线图等三种曲线，并计算系统的频域性能指标：剪切频率(pnl)（也

33、叫截止频率(pnl)）、- 穿越频率、相角稳定裕度、幅值稳定裕度，以便研究系统控制过程的稳定性、快速性及稳态精度等性能。控制系统(kn zh x tn)的MATLAB仿真50共一百零二页第7章 2MATLAB中频域分析的实现方法 频域分析在经典控制里占有重要的地位。借助于MATLAB软件，运行它所提供的频率分析函数，能够方便、简单、快捷地绘制Bode图、Nyquist曲线和Nichols图，并计算出频域性能指标，还可以借助于这些(zhxi)曲线对系统进行分析。 应用MATLAB提供的上述频域仿真函数与其他函数命令、语句，编制成MATLAB程序，这种在MATLAB的指令方式下进行频域仿真同样是最

34、常用的实现方法。控制系统(kn zh x tn)的MATLAB仿真51共一百零二页第7章 7.4.2 求连续系统的Bode图1求连续系统Bode图的函数 在MATLAB中求连续系统Bode图的函数为bode()，该函数命令的调用(dioyng)格式和功能为：（1）mag,phase,w=bode(sys)； （2）bode(sys)；（3）bode(sys,w)； （4）bode(sysl,sys2,sysn)或bode(sysl,sys2,，sysn,w)； 控制系统(kn zh x tn)的MATLAB仿真52共一百零二页第7章 2求离散系统Bode图的函数(hnsh)dbode()求离散

35、系统Bode图的函数为dbode()，该函数命令的调用格式为：（1）mag,phase,w=dbode(a,b,c,d,Ts)（2）mag,phase,w=dbode(a,b,c,d,Ts,iu) （3）mag,phase,w=dbode(a,b,c,d,Ts,iu,w)或 mag,phase,w=dbode(num,den,Ts,w) （4）mag,phase,w=dbode(num,den,Ts) 控制系统(kn zh x tn)的MATLAB仿真53共一百零二页第7章 7.4.3 绘制系统Nyquist曲线的函数 计算(j sun)或绘制系统Nyquist曲线的函数为nyquist()，

36、该函数命令的调用格式为：（1）re,im,w=nyquist(sys)（2）nyquist(sys)（3）nyquist(sys,w)（4）nyquist(sysl,sys2,sysn) （5）nyquist(sysl,sys2,sysn,w)控制系统(kn zh x tn)的MATLAB仿真54共一百零二页第7章 如果系统为离散系统，则相应的nyquist函数为 dnyquist()，函数命令(mng lng)调用为：（1）re,im,w= dnyquist(a,b,c,d,Ts)（2）re,im,w= dnyquist(a,b,c,d,Ts,iu)（3）re,im,w= dnyquist(

37、a,b,c,d,Ts,iu,w)（4）re,im,w= dnyquist(num,den,Ts)（5）re,im,w= dnyquist(num,den,Ts,w)控制系统(kn zh x tn)的MATLAB仿真55共一百零二页第7章 7.4.4 求连续系统(xtng)Nichols曲线的函数 求连续系统Nichols曲线的函数为nichols()，该函数命令的调用格式为：（1）mag,phase,w= nichols(sys)（2）nichols(sys)（3）nichols(sys,w)（4）nichols(sysl,sys2,sysn)（5）nichols(sysl,sys2,，sys

38、n,w)控制系统(kn zh x tn)的MATLAB仿真56共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真 7.4.5 求系统幅值裕度与相位裕度 求系统幅值裕度与相位裕度的函数为margin()，该函数命令(mng lng)的调用格式为：（1）Gm,Pm,Wcp,Wcg= margin(sys)； （2）Gm,Pm,Wcp,Wcg= margin(mag,phase,w)（3）margin(sys)； （4）margin(mag,phase,w)； 对于离散时间系统，可以先用函数dbode()计算频率响应，然后再调用margin()函数，即用margin(mag,pha

39、se,w)函数命令，再求系统性能指标或绘制Bode图。 57共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真 7.5 利用MATLAB语言进行编程仿真的应用7.5.1 面向传递函数的线性系统仿真 对线性连续系统进行单输入(shr)、单输出的仿真，可以将用户输入(shr)的系统传递函数模型转化为仿真计算模型，再应用数值积分法进行处理。采用四阶龙格库塔法保证系统仿真过程中具备一定的精度和性能指标要求，并可实现重复运行，便于研究参数的变化对系统动态性能的影响，且运行过程直观、形象，修改参数容易。58共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真 1仿真系统的

40、典型(dinxng)结构框图如图7-11所示的典型闭环系统结构，有：图7-11 典型闭环系统结构图 59共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真 2仿真程序框图与实现（1）程序框图的描述 如图7-12所示。（2）程序设计的实现 用MATLAB语言(yyn)编写的面向传递函数的线性系统仿真程序定义为CSS1(Control System Simulation 1) 。60共一百零二页第7章图7-12 面向(min xin)传递函数的线性系统仿真框图61共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真设源程序名为CSS1.m，保存为M文件。1）输入数

41、据a=;% n+1维分母系数向量c=;% m+1维分子系数向量 X0=;% 状态(zhungti)向量初值V=;% 反馈系数n=;% 系统阶次T0=;% 起始时间Tf=;% 终止时间h=;% 计算步长R=r;% 阶跃输入函数幅值 62共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真 2）形成开、闭环系数(xsh)阵c=c/a(1); % 首一化处理a=a/a(1);A=a(2:n+1); A=rot90(rot90(eye(n-1,n);-fliplr(A);% 形成能控标准形B=zeros(1,n-1),1;% 形成输入阵Bm1=length(c);% 分子系数向量维数m

42、+1C=fliplr(c),zeros(1,n-m1);% 形成输出阵CAb=A-B*C*V;% 形成闭环系数阵AbX=X0;y=0;t=T0;% 设初值，准备递推运算63共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真 3）运算求解N=round(Tf-T0)/h;% 确定输出点数for i=1:N% 四阶龙格-库塔法K1=Ab*X+B*R;K2=Ab*(X+h*K1/2)+B*R;K3=Ab*(X+h*K2/2)+B*R;K4=Ab*(X+h*K3)+B*R;% 求各次斜率KX=X+h*(K1+2*K2+2*K3+K4)/6;% 求状态(zhungti)y=y,C*X;

43、% 求输出并以向量形式保存 t=t,t(i)+h;% 输出对应时刻以向量形式保存end64共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真 4）输出结果t,y % 输出数据形式结果plot(t,y)% 输出曲线形式结果 以上即为采用MATLAB语言实现面向典型闭环系统数字仿真的程序，可编辑为CSS1.m文件存储起来。使用时，只要进入到MATLAB命令窗口，按要求输入系统参数和运行(ynxng)参数，再调用该文件，即可得到系统的仿真结果。65共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真 该程序中用到了几个MATLAB特殊功能函数，简要说明(shumng

44、)如下：eye(m,n)产生mn型I阵，即主对角元素均为1，其余元素为零的长方阵length( )求取括号内向量的维数zeros(m,n)产生mn型零阵rot90( )将括号内矩阵左旋90fliplr( )将括号内矩阵或向量左右翻转round( )对括号内运算取整plot(t,y)打印以t为横轴，y为纵轴的y(t)曲线 66共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真 7.5.2 面向结构图的线性系统仿真1仿真程序框图与程序实现（1）程序框图的描述 当控制系统的数学模型以动态结构图的形式给出时，可以用以下的方法处理，如图7-15所示。（2）程序设计的实现 采用MATLA

45、B语言编程，程序(chngx)名为CSS2.M，保存为M文件。67共一百零二页图7-15 面向(min xin)结构图的线性系统的仿真框图第7章68共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真输入系统原始数据：P % 输入各环节参数WIJ=% 输入连接阵非零元素(yun s)n=;% 环节个数（系统阶次）Y0=;% 阶跃输入幅值Yt0=% 各环节输出初值h=;% 计算步长L1=;%每隔L1个点打印一次 T0=;% 起始时间 Tf=;% 终止时间nout=Nout;% 输出环节编号 69共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真程序处理开始：A=

46、diag(P(:,1);B=diag(P(:,2);% 求开环系统系数(xsh)阵C=diag(P(:,3);D=diag(P(:,4);m=length(WIJ(:,1);% 求非零元素个数W0=zeros(n,1);W=zeros(n,n);% 建立初始W、W0阵for k=1:mif (WIJ(k,2)=0);W0(WIJ(k,1)=WIJ(k,3);% 求W0阵else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3);% 求W阵endend70共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真 Q=B-D*W;Qn=inv(Q);% 求Q和Q逆阵R=C*W-

47、A;V1=C*W0;% 求R,V1阵Ab=Qn*R;b1=Qn*V1;% 形成(xngchng)闭环系数阵Y=Yt0; y=Y(nout);t=T0;% 置初值N=round(Tf-T0)/(h*L1);% 总输出点数for i=1:N; % 每循环一次，输出一点数据for j=1:L1;% 每输出点之间计算L1次K1=Ab*Y+b1*Y0; K2=Ab*(Y+h*K1/2)+b1*Y0;K3=Ab*(Y+h*K2/2)+b1*Y0;K4=Ab*(Y+h*K3)+b1*Y0;Y=Y+h*(K1+2*K2+2*K3+K4)/6;% 龙格库塔法求解方程end71共一百零二页第7章控制系统(kn z

48、h x tn)的MATLAB仿真y=y,Y(nout);% 保存输出环节动态响应值t=t,t(i)+h*L1;% 保存时间向量endt,y% 输出响应数据plot(t,y)% 输出响应图形 程序中用到MATLAB特殊功能函数diag( )，其作用是产生以括号内向量元素(yun s)为主对角元素(yun s)的对角阵；P(:,I)用于取P阵中第I列的所有元素(yun s)作为单独列向量。72共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真7.5.3 线性系统的离散相似法仿真1仿真框图描述 离散相似法仿真线性系统的程序框图如图7-18所示。2仿真程序设计 本程序定义为CSS3.

49、M，输入(shr)系统各环节参数阵P，连接阵WIJ，运行参数n,h,L1,T0,Tf，阶跃输入(shr)幅值Y0，输出环节编号nout等，方法同CSS2程序。73共一百零二页图7-18 离散相似(xin s)法仿真线性系统程序框图第7章74共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真CSS3仿真程序设计如下：A=diag(P(:,1);B=diag(P(:,2);% 求开环系统(xtng)系数阵C=diag(P(:,3);D=diag(P(:,4);m=length(WIJ(:,1);% 求非零元素个数W0=zeros(n,1);W=zeros(n,n);% 建立初始W

50、、W0阵for k=1:mif (WIJ(k,2)=0);W0(WIJ(k,1)=WIJ(k,3);% 求W0阵else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3);% 求W阵endend75共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真for i=1:nif (A(i,i)=0);% 求积分或积分比例环节各系数FI(i)=1;FIM(i)=h*C(i,i)/B(i,i);FIJ(i)=h*h*C(i,i)/B(i,i)/2;FIC(i)=1;FID(i)=0;% 求积分环节系数if (D(i,i)=0);% 若为积分比例环节，修正(xizhng)系数FID

51、(i)=D(i,i)/B(i,i);elseendelse% 求惯性或惯性比例环节各系数FI(i)=exp(-h*A(i,i)/B(i,i);FIM(i)=(1-FI(i)*C(i,i)/A(i,i);FIJ(i)=h*C(i,i)/A(i,i)-FIM(i)*B(i,i)/A(i,i);76共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真FIC(i)=1;FID(i)=0;% 求惯性环节系数if (D(i,i)=0);% 若为惯性比例环节，修正(xizhng)系数FIC(i)=C(i,i)/D(i,i)-A(i,i)/B(i,i);FID(i)=D(i,i)/B(i,i)

52、;elseendendendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ub=Uk;t=T0;% 置初值 N=round(Tf-T0)/(h*L1);% 计算总输出点数 77共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真for k=1:Nfor l=1:L1Ub=Uk; % 保存前一次输入值Uk=W*Y+W0*Y0; % 求当前UkUdot=(Uk-Ub)/h; % 求当前Uk的导数(do sh)Uf=2*Uk-Ub; % 求下一步的输入UfX=FI.*X+FIM.*Uk+FIJ.*Udot;% 求状态向量Y=FIC.*X+FID.*Uf;% 求

53、输出endy=y,Y(nout);t=t,t(k)+h*L1;endt,yplot(t,y)78共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真该程序用到了几个MATLAB语言中的特殊运算符： 冒号符，用于生成行向量；= = 逻辑等号符，MATLAB中的规定写法； = 逻辑不等号符，MATLAB中的规定写法； * 点乘运算符，MATLAB中的规定写法，表示两个同维矩阵中各相应(xingyng)位置上的元素相乘运算。79共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真7.5.4 非线性系统离散相似法仿真1典型非线性环节的求解（1）饱和非线性环节：采用(ci

54、yng)MATLAB语言来编制M函数如下，函数名为satu( )：function Uc=satu(Ur,S1)if(abs(Ur)=S1)if(Ur0)Uc=S1;else Uc=-S1;endelse Uc=Ur;end 80共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真（2）死区非线性：采用MATLAB语言(yyn)来编制M函数如下，函数名为dead( )：function Uc=dead(Ur,S1)if (abs(Ur)=S1)if(Ur0)Uc=Ur-S1;else Uc=Ur+S1;endelse Uc=0;end81共一百零二页第7章控制系统(kn zh

55、x tn)的MATLAB仿真（3）滞环非线性：M函数(hnsh)名为backlash ( )：function Uc,Ubb=backlash(Urb,Ur,Ucb,S1)if (UrUrb)if(Ur-S1)=Ucb)Uc=Ur-S1；else Uc=Ucb；else if(UrUrb) if(Ur-S1)=Ucb)Uc=Ur+S1；else Uc=Ucb；else Uc=Ucb；endendUbb=Ur； 82共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真非线性系统的离散相似法仿真程序命名为css4.m。根据上述分析可以画出CSS4程序框图如7-21所示。CSS4程序

56、设计：输入数据(shj)：输入系统各环节参数阵P，连接阵WIJ，运行参数n,h,L1,T0,Tf,阶跃输入幅值Y0，输出环节编号nout(方法同CSS3)，以及非线性标志向量Z，参数向量S。程序开始：求开环各系数阵A,B,C,D 及连接关系矩阵 W,W0 (方法同CSS3)；求取各环节离散系数FI(i),FIM(i),FIJ(i),FIC(i),FID(i)。83共一百零二页第7章图7-21 离散相似(xin s)法仿真非线性系统程序框图84共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真求解运行(ynxng)部分如下：Y=zeros(n,1);X=Y;y=0; Uk=ze

57、ros(n,1);Ubb=Uk;t=T0;N=round(Tf-T0)/(h*L1);for k=1:Nfor l=1:L1Ub=Uk; Uk=W*Y+W0*Y0;for i=1:nif(Z(i)=0)if (Z(i)=1)Uk(i)=satu(Uk(i),S(i);end 85共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真if (Z(i)=2) Uk(i)=dead(Uk(i),S(i);endif (Z(i)=3) Uk(i),Ubb(i)=backlash(Ubb(i),Uk(i),Ub(i),S(i);endend end Udot=(Uk-Ub)/h;Uf=2*

58、Uk-Ub;X=FI.*X+FIM.*Uk+FIJ.*Udot; Yb=Y;Y=FIC.*X+FID.*Uf;86共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真for i=1:nif(Z(i)=0)if(Z(i)=4)Y(i)=satu(Y(i),S(i);endif(Z(i)=5)Y(i)=dead(Y(i),S(i);endif (Z(i)=6)Y(i),Ubb(i)=backlash(Ubb(i),Y(i),Yb(i),S(i);87共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真endendendendy=y,Y(nout);t=t,t(k)+

59、h*L1; endt,yplot(t,y) 88共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真7.5.5 采样控制系统仿真根据图7-25中所示的典型(dinxng)数字采样控制系统结构来分析采样控制系统的仿真方法。图7-25 采样控制系统结构图89共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真【例7.13】采样控制系统如图7-31所示，数字控制器D(Z)按有纹波最小拍系统设计，取采样周期T=1秒，外部输入信号(xnho)为单位阶跃函数r(t)=1(t)，仿真该系统的动态响应过程。有纹波最小拍系统： 90共一百零二页第7章控制系统(kn zh x tn

60、)的MATLAB仿真图7-31 采样(ci yn)控制系统结构图解：仿真该采样控制系统可以采用连续部分按环节离散化处理的方法，对数字控制器的仿真仍然采用离散化处理，后面连续的被控对象按连续系统的离散相似法仿真。91共一百零二页第7章控制系统(kn zh x tn)的MATLAB仿真（1）根据框图编写如下(rxi)的MATLAB程序，保存为CSS5.M。A=diag(P(:,1);B=diag(P(:,2); C=diag(P(:,3);D=diag(P(:,4);m=length(WIJ(:,1);W0=zeros(n,1);W=zeros(n,n);for k=1:mif (WIJ(k,2)