自动控制原理的MATLAB仿真与实践第3章-线性系统的时域分析课件

11第3章主要内容：

3.1时域分析的拉普拉斯变换法3.1.1连续时间函数的拉普拉斯变换3.1.2时域函数的拉氏反变换法3.1.3时域函数的部分分时展开法3.2时域分析的函数命令方法3.3作业与实验2第3章主要内容：23.1时域分析的拉普拉斯变换法3

利用拉普拉斯变换法对线性系统进行时域分析的一般步骤：（1）对系统的传递函数模型进行部分分式展开，将其变为简单传递函数之和；（2）利用拉普拉斯反变换，得到系统的输出时间响应函数；（3）绘制系统的响应曲线;（4）通过改变系统的参数，观察系统输出响应的变化情况，对系统的时域特性（如性能指标等）进行分析。本节介绍上述过程的MATLAB实现方法。3.1时域分析的拉普拉斯变换法3利用拉普拉斯变2.1.1生成模型的常用函数命令格式4连续时间信号f(t)的单边拉普拉斯变换F(s)的定义：记为：拉普拉斯反变换的定义：记为：2.1.1生成模型的常用函数命令格式4连续时间信号f(t)s=tf(‘s’)：生成以s为变量的递函数s。此时，s既是传递函数也是指定变量。

sys=tf(num,den)：生成传递函数模型，num,den分别为模型的分子和分母多项式系数向量。sys=tf(num,den,‘td’,v)：生成延迟时间td=v的传递函数模型。sys=zpk(z,p,k)：生成零极点增益模型，z,p,k分别为零点、极点和增益向量。[num,den]=rmodel(n,P)：随机生成一个n阶连续的传递函数模型，该系统具有P个输出。[num,den]=ord2(wn,z)：生成固有频率为wn，阻尼系数为z的连续二阶系统模型系统。[num,dent]=pade(L,n)：返回延迟环节G(s)=e-Ls近似为n阶多项式传递函数的num和den。5s=tf(‘s’)：生成以s为变量的递函数s。此时，s既是传2.1.2有理多项式分式传递函数模型的建立6

控制理论中，系统传递函数的有理多项式分式模型通常表示为：

在MATLAB中，多项式是用其降幂系数的排列构成的数组来表示的，称为多项式系数向量。

系统传递函数的分子和分母的多项式系数向量分别表示为：num=[bm,bm-1,…,b1,b0]den=[an,an-1,…,a1,a0]2.1.2有理多项式分式传递函数模型的建立67

在MATLAB语言中有专门对信号进行正、反拉氏变换的函数命令。常用拉氏变换的函数命令格式及其说明见表3-1。F=laplace(f)：f(t)的拉氏变换，结果为F(s)，默认变量为sF=laplace(f,v)：f(t)拉氏变换，并用v代替s，结果为F(v)F=laplace(f,v,u)：f(v)拉氏变换，并用u代替v，结果为F(u)f=ilaplace(F)

：F(s)的拉氏反变换，结果为f(t)，变量为tf=ilaplace(F,u)：F(s)拉氏反变换，用u代替s，结果为f(u)f=ilaplace(F,v,u)：

F(v)拉氏反变换，用u代v，结果为f(u)说明：默认情况下，MTLAB给出的是单边拉氏变换；使用函数laplace()及ilaplace()之前，需要用syms命令对所有要用到的变量进行符号变量的定义。7在MATLAB语言中有专门【例3-1】试求函数

的拉氏变换式，并用拉氏反变换观察变换结果。解：MATLAB程序如下：

>>clear;%清除所有变量symstAwbs

%定义符号变量t,A,w,b,sft=A*sin(w*t+b);%定义f(t)的符号函数ft的表达式Fs=laplace(ft)

%求ft的拉氏变换式Fs,即F(s)运行结果：Fs=A*(cos(b)*w/(s^2+w^2)+sin(b)*s/(s^2+w^2))8【例3-1】试求函数当取参数A=1、b=0、w=1时，带入Fs>>A=1;b=0;w=1;%设定参数值Fs=A*(cos(b)*w/(s^2+w^2)+sin(b)*s/(s^2+w^2))

%计算给定参数下的F(s)运行结果：Fs=1/(s^2+1)即得到

F(s)=L[f(t)]=L[sin(t)]=1/(s2+1)9当取参数A=1、b=0、w=1时，带入Fs9

可利用拉氏反变换对上述结果进行检验：>>ft=ilaplace(Fs)%求Fs的拉氏反变换式ft运行结果：ft=sin(t)即

f(t)=L-1[F(s)]=L-1[1/(s2+1)]=sin(t)说明:MATLAB中的拉氏变换的函数命令laplace()，在默认情况下是指拉氏右变换，其运行结果是单边函数。如L[f(t)]=laplace(sin(t))，实际上是指对ft=sin(t)u(t)的双边拉氏变换，

其中u(t)为单位阶跃函数。10可利用拉氏反变换对上述结果进行检验：10

利用MATLAB拉氏反变换可以求出传递函数的时间函数，或系统的时域输出响应。【例3-3】试求传递函数

的时间函数g(t)。解：求解G(s)的拉氏反变换程序如下：113.1.2时域函数的拉氏反变换法利用MATLAB拉氏反变换可以求出传递函数1>>clear;symss;Gs=(s^3+4*s^2+8*s+5)/(s^2+3*s+4);gt=ilaplace(Gs)运行结果：gt=Dirac(t)+Dirac(1,t)-1/7*exp(-3/2*t)*7^(1/2)*sin(1/2*7^(1/2)*t)+exp(-3/2*t)*cos(1/2*7^(1/2)*t)其中，即得：12>>clear;symss;12【例3-4】某质量-弹簧-阻尼(m-k-c)系统如图3-1所示。其中，f(t)为外加力，x(t)为m的位移。试求：1）系统的MATLAB传递函数；2）输入为单位阶跃（f(t)=u(t)）时的时间响应函数x(t)；3）绘制m=3;k=2;c=1时的时域响应曲线。13【例3-4】某质量-弹簧-阻尼(m-k-c)系统如图3-1所解：1）由牛顿第二定律可得到该系统模型的数学表示：上式的拉氏变换后的传递函数模型为：生成上述模型的MATLAB模型：>>symsmkcsGs=1/(m*s^2+c*s+k)%生成G(s)的符号传递函数Gs运行结果：Gs=1/(m*s^2+c*s+k)14解：1）由牛顿第二定律可得到该系统模型的数学表示：14152）f(t)=u(t)的拉氏变换为F(s)=1/s。所以输出传函为求上式的MATLAB拉氏反变换：>>Xs=Gs*(1/s)%建立输出传函X(s)xt=ilaplace(Xs)%求解输出X(s)的时域响应

152）f(t)=u(t)的拉氏变换为F(s)=1/s。所以16运行结果：xt=-4*exp(-1/2*c/m*t)/(4*k*m-c^2)*m*cos(1/2*((4*k*m-c^2)/m^2)^(1/2)*t)+1/k*exp(-1/2*c/m*t)/(4*k*m-c^2)*c^2*cos(1/2*((4*k*m-c^2)/m^2)^(1/2)*t)-1/k*exp(-1/2*c/m*t)/(4*k*m-c^2)*m*((4*k*m-c^2)/m^2)^(1/2)*c*sin(1/2*((4*k*m-c^2)/m^2)^(1/2)*t)+1/k3）绘制m=3;k=2;c=1时的响应曲线。程序如下：>>symstm=3;k=2;c=1;%参数赋值

16运行结果：xt=-4*exp(-1/2*c/m*t)/(4*k*m-c^2)*m*cos(1/2*((4*k*m-c^2)/m^2)^(1/2)*t)+1/k*exp(-1/2*c/m*t)/(4*k*m-c^2)*c^2*cos(1/2*((4*k*m-c^2)/m^2)^(1/2)*t)-1/k*exp(-1/2*c/m*t)/(4*k*m-c^2)*m*((4*k*m-c^2)/m^2)^(1/2)*c*sin(1/2*((4*k*m-c^2)/m^2)^(1/2)*t)+1/k运行结果：xt=-1/2*exp(-1/6*t)*cos(1/6*23^(1/2)*t)-1/46*exp(-1/6*t)*23^(1/2)*sin(1/6*23^(1/2)*t)+1/217xt=-4*exp(-1/2*c/m*t)/(4*k*m-c>>t=[0:0.1:10]*pi;%取绘图区间采样点xt=-1/2*exp(-1/6*t).*cos(1/6*23^(1/2)*t)-...1/46*exp(-6*t).*23^(1/2).*sin(1/6*23^(1/2)*t)+1/2;%计算所有采样点上的函数值plot(t,xt);grid%绘图，加栅格结果如图3-2所示。18>>t=[0:0.1:10]*pi;1919

对于复杂的分式多项式传函模型，首先，采用部分分式法将其分解成如下部分分式展开的形式：其中，pi为极点，

ri为有理函数Y(s)的留数，即k为余项；其次，对其求拉氏反变换；最后，由拉氏变换简表得到时间响应函数。203.1.3时域函数的部分分时展开法对于复杂的分式多项式传函模型，首先，采用部20321通常有下列几中情况及其相应的处理方法：1）单极点时，即pi为不相同极点。则Y(s)的拉普拉斯逆变换为：2）有m重极点时。假设有m个重极点

，Y(s)的拉普拉斯逆变换为：3）有共轭极点时，逆变换较为复杂，需具体分析处理（参考教材），

21通常有下列几中情况及其相应的处理方法：，22

在MATLAB中，对于上述三种情况均可利用如下函数命令格式[r,p,k]=residue(num,den)对分式多项式进行部分分式展开，然后再根据极点的不同情况进行相应的处理。其中，num,den为Y(s)的分子和分母多项式系数向量，[r,p,k]为返回的留数，极点，

和余项。22在MATLAB中，对于上述三种情况均可利用如下【例3-5】试求输出传递函数的时间响应。解：MATLAB程序如下：>>clear;num=[1-2];den=[13310];%系数向量[r,p,k]=residue(num,den)%求部分分式参数运行结果：r=2.00002.00003.0000-2.000023【例3-5】试求输出传递函数23p=-1.0000-1.0000-1.00000k=[]可见，系统有3个重根p1=p2=p3=-1和p4=0，则由式（3-3）和式（3-4）得到由式（3-4）知，m=3,对应的逆变换为y(t)=L-1[Y(s)]。可以直接利用拉氏变换表得出结果。也可以通过MATLAB求解，程序如下：24p=24>>symssft1=ilaplace(3/(s+1)^3),ft2=ilaplace(2/(s+1)^2),ft3=ilaplace(2/(s+1)),ft4=ilaplace(-2/s)%求各分式逆变换得结果：ft1=3/2*t^2*exp(-t)ft2=2*t*exp(-t)ft3=2*exp(-t)ft4=-2>>ft=ft1+ft2+ft3+ft4%求f(t)的结果：ft=3/2*t^2*exp(-t)+2*t*exp(-t)+2*exp(-t)-2即25>>symss25【例3-7】设一个RLC网络结构如图所示。其中ui为输入信号，uc为输出信号，网络各变量的初始状态均为零。试求：1）建立RLC网络的传递函数

；2）当ui=Eu(t)，用拉氏逆变换函数命令求uc；3）绘制不同参数时的单位阶跃响应uc(t)。解：1）建立RLC网络的传递函数

(1)建立时域模型。根据R、L、C电压定律及回路电压方程

，可得RLC网络的微分方程模型：其中，

26【例3-7】设一个RLC网络结构如图所示。26（2）建立传函模型

对上述微分模型的两边取拉氏变换，经整理得到RLC网络拉氏变换的标准二阶系统传递函数：2）当ui=Eu(t)时，由于L[ui]=L[E]=E/s，所以输出的拉氏变换为：利用MATLAB拉氏逆变换函数命令求uc。程序如下：27（2）建立传函模型27>>clear

;symstsEwnetwd%定义变量t,s,E,ωn,ζUc=E*wn^2/((s^2+2*et*wn*s+wn^2)*s);%建立函数uc=ilaplace(Uc,s,t)%拉氏逆变换求输出uc(t)simple(uc)%简化结果结果：ans=E*wn^2*(1/wn^2+1/(4*et^2*wn^2-4*wn^2)^(1/2)*(1/(-et*wn+1/2*(4*et^2*wn^2-4*wn^2)^(1/2))*exp((-et*wn+1/2*(4*et^2*wn^2-4*wn^2)^(1/2))*t)-1/(-et*wn-1/2*(4*et^2*wn^2-4*wn^2)^(1/2))*exp((-et*wn-1/2*(4*et^2*wn^2-4*wn^2)^(1/2))*t)))28>>clear

;283）绘制不同参数时的单位阶跃响应uc(t)。>>clear

;E=1;wn=3;t=0:0.02:6;%设置参数E和wn,时间采样点tet=[0.1,0.2,0.5,0.95];%设置4个不同的阻尼向量etfori=1:4;%对et循环计算uc=E*wn^2*(1/wn^2+1/(4*et(i)^2*wn^2-4*wn^2)^(1/2)*(1/(-et(i)*wn+1/2*(4*et(i)^2*wn^2-4*wn^2)^(1/2))*exp((-et(i)*wn+1/2*(4*et(i)^2*wn^2-4*wn^2)^(1/2)).*t)-1/(-et(i)*wn-1/2*(4*et(i)^2*wn^2-4*wn^2)^(1/2)).*exp((-et(i)*wn-1/2*(4*et(i)^2*wn^2-4*wn^2)^(1/2)).*t)));plot(t,uc,'k-');holdon;%绘制响应曲线并保持图形end;293）绘制不同参数时的单位阶跃响应uc(t)。29title('uc(t)的二阶欠阻尼响应曲线')gtext('et=0.1')%在图中对应的响应曲线的合适位置%点击鼠标左键添加阻尼参数标示'et=0.1'gtext('et=0.2')gtext('et=0.5')gtext('et=0.95')xlabel('t');ylabel('uc(t)');%添加x和y轴标示't'和'uc(t)'grid;%添加网格运行得到响应曲线见图3-4。30title('uc(t)的二阶欠阻尼响应曲线')3031313.2时域分析的函数命令方法

上述介绍的MATLAB拉氏变换法，不但可以用来对通常情况下的微分方程、分式多项式、传递函数及系统结构等函数进行数学上的处理与分析，还可用于控制系统的时域分析。MATLAB还提供了特殊的函数命令，专用于对系统的分析。本节介绍常用MATLAB线性系统时域分析的函数命令格式及其使用方法。表3-2列出了常用时域响应函数命令格式及说明。323.2时域分析的函数命令方法上述介绍的MA【例3-8】已知系统传递函数为

试求其单位脉冲响应。解：使用专用函数命令来实现，程序如下：>>clear;n=[3];d=[1,1,3];%传递函数的系数向量Gs=tf(n,d);%建立模型G(s)impulse(Gs)%G(s)的单位脉冲响应运行结果如图3-4所示。33【例3-8】已知系统传递函数为333434【例3-9】已知系统传递函数为求其单位阶跃响应。解：程序如下：>>clear;n=[1,7];d=[1,1.5,8,7];Gs=tf(n,d);step(Gs)运行结果如图3-5所示。35【例3-9】已知系统传递函数为3536【例3-14】已知某二阶系统的超调量Mp=16.3%。试根据下列条件确定二阶系统的开环传递函数G(s)：1）上升时间tr=0.1s；2）峰值时间tp=0.114s；3）调节时间ts=2s。

解：由二阶系统模型结构可知，其传递函数由两个参数决定：阻尼比ζ和自然振荡频率ωn。可根据二阶系统的性能指标中的任意两个进行求解即可确定模型参数。所以可以利用MATLAB的方程求解函数命令实现。具体格式：

X=solve('eq1','eq2',...,'eqn',val1,val2,...,valn)其中，eqi和vali，i=1,2,…,n,分别为方程和待求变量。36【例3-14】已知某二阶系统的超调量Mp=16.3%。371）求解同时满足超调量Mp=16.3%和上升时间tr=0.1s时的阻尼比ζ和自然振荡频率ωn。即求解下列方程组中的两个变量ζ和ωn。设ζ=x，ωn=y，则上述方程改写为371）求解同时满足超调量Mp=16.3%和上升时间tr=038计算程序如下：>>clear;symsxyMptrx=solve('Mp-exp(-x*pi/sqrt(1-x^2))=0',x);%求x=ζx=subs(x,Mp,0.163)y=solve('tr-(pi-atan(sqrt(1-x^2)/x))/(y*sqrt(1-x^2))=0',y);%求y=ωny=subs(y,tr,0.1);y=subs(y,x)得到的结果为：x=0.5000；y=24.1852。即得二阶系统模型为：38计算程序如下：392）求解同时满足超调量Mp=16.3%和峰值时间tp=0.114s时的阻尼比ζ和自然振荡频率ωn。即求解下列方程组计算程序如下：>>clear;symsxyMptpx=solve('Mp-exp(-x*pi/sqrt(1-x^2))=0',x);%求阻尼比x=subs(x,Mp,0.163)y=solve('tp-pi/(y*sqrt(1-x^2))=0',y);%求y=ωny=subs(y,tp,0.114);y=subs(y,x)所得结果：x=ζ=0.5000；y=ωn=31.8219。将参数ζ和ωn代入G(s)即可。392）求解同时满足超调量Mp=16.3%和峰值时间tp=040【例2-12】已知传递函已知传递函数G1和G2分别为试以G1作为前向通道，以G2作为反馈通道时的负反馈等效闭环模型G。解：实现反馈的等效传函程序如下：方法1：利用多项式系数求解。>>clear;num1=[1,11];den1=[1,3,21];num2=[2,1];den2=[3,5,