最新matlab_课件第7章_控制系统的仿真(共36张PPT课件).pptx

《最新matlab_课件第7章_控制系统的仿真(共36张PPT课件).pptx》由会员分享，可在线阅读，更多相关《最新matlab_课件第7章_控制系统的仿真(共36张PPT课件).pptx（36页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第第7章章 控制系统仿真控制系统仿真常微分方程常微分方程(wi fn fn chn)的数值求解的数值求解微分方程模型微分方程模型传递函数模型传递函数模型第一页，共三十六页。7.1 微分方程的求解微分方程的求解 在现在数学研究和工程实践中，很多数学模型都是用微分在现在数学研究和工程实践中，很多数学模型都是用微分方程确定的，很多基本方程本身就是一个微分方程，因此求微方程确定的，很多基本方程本身就是一个微分方程，因此求微分方程非常重要，但是分方程非常重要，但是(dnsh)大部分的微分方程目前难以求得其大部分的微分方程目前难以求得其解析解，因此人们只有利用计算机强大的计算功能来求其数值解析解，因此人们

2、只有利用计算机强大的计算功能来求其数值解。解。MATLAB主要使用龙格主要使用龙格-库塔法求解微分方程。库塔法求解微分方程。 在控制系统仿真中，常用的求微分方程数值解的函数是在控制系统仿真中，常用的求微分方程数值解的函数是ode23和和ode45。第二页，共三十六页。1. ode23 在在MATLAB中，函数中，函数ode23采用采用2-3阶龙格阶龙格-库塔法求库塔法求解微分方程解微分方程(wi fn fn chn)。 t,y=ode23(odefun,tspan,y0) t,y=ode23(odefun,tspan,y0,options) odefun：定义微分方程的形式：定义微分方程的形式

3、y=f(t,y) tspan=t0,tfinal：表示微分方程的积分限从：表示微分方程的积分限从t0（始（始 值）到值）到tfinal（终值），该积分限也（终值），该积分限也 可以是一些离散的点。可以是一些离散的点。 y0:初始状态列向量初始状态列向量 options：积分参数，包括：积分参数，包括RelTol（相对误差）和（相对误差）和 AbsTol（绝对误差），可省略。（绝对误差），可省略。第三页，共三十六页。例：使用例：使用ode23函数求解常微分方程函数求解常微分方程y=-y+x2+4x+1，x=1时，时，y=1。解：首先创建解：首先创建(chungjin)函数函数fun1.mfunc

4、tion dy=fun1(x,y)dy=-y+x2+4*x+1;在命令窗口中输入在命令窗口中输入 x,y=ode23(fun1,1,4,1); dy=-y+x.2+4*x+1; plot(x,y,x,dy); legend(y,dy)第四页，共三十六页。2. ode45 在在MATLAB中，函数中，函数ode45采用普通采用普通(ptng)4-5阶龙格阶龙格-库塔法求解微分方程。其使用方法与库塔法求解微分方程。其使用方法与ode23函数的使用方函数的使用方法基本相同。法基本相同。 ode45函数是大部分场合的首选算法，函数是大部分场合的首选算法，ode23函数主函数主要适用于精度较低的情形。要

5、适用于精度较低的情形。例：解经典非线性方程，范得波（例：解经典非线性方程，范得波（Van der Pol）微分方程）微分方程 （w=2）。）。0)1(222 xdtdxxdtxd 当当t=0时，时，dx/dt=1，d2x/dt2=0。第五页，共三十六页。0)1(222 xdtdxxdtxd 解解：（：（1）将高阶微分方程式等价变换成一阶微分方程组。）将高阶微分方程式等价变换成一阶微分方程组。令令y1=x且且y2=dx/dt则则dy1/dt=y2，dy2/dt=w(1-y12)y2-y1（2）编写）编写M文件表示文件表示(biosh)该微分方程，给定当前时间及该微分方程，给定当前时间及y1和和y

6、2的当前值，返回上述的导数值，并将导数值以列向量的形的当前值，返回上述的导数值，并将导数值以列向量的形式给出。式给出。function fun2=vdpol(t,y)fun2=y(2) 2*(1-y(1)2)*y(2)-y(1)第六页，共三十六页。（3）计算结果如下：）计算结果如下： t,y=ode45(vdpol,0 30,1;0); y1=y(:,1); y2=y(:,2); plot(t,y1,:b,t,y2,-r) legend(速度速度(sd),位移位移)第七页，共三十六页。3. 定积分的数值解法定积分的数值解法 MATLAB软件使用软件使用quad函数进行定积分的数值解法。函数进行

7、定积分的数值解法。使用格式为：使用格式为： q = quad(fun,a,b) fun：被积分函数：被积分函数 a、b：积分上下限：积分上下限(xixin)例：计算下列定积分例：计算下列定积分xxxd521203 function y = myfun(x) y = 1./(x.3-2*x-5); Q = quad(myfun,0,2)Q = -0.4605F = (x)1./(x.3-2*x-5);Q = quad(F,0,2);第八页，共三十六页。7.2 微分方程模型微分方程模型7.2.1 方法描述方法描述 微分方程模型是数学模型的一种主要形式。当采用一微分方程模型是数学模型的一种主要形式。

8、当采用一阶微分方程的数值积分法进行数值计算时，阶微分方程的数值积分法进行数值计算时，应把高阶微分方应把高阶微分方程变换成程变换成n个一阶微分方程形式个一阶微分方程形式。对于微分方程而言，除了少。对于微分方程而言，除了少数可以得到解析解外，大多数只能采用数值解法数可以得到解析解外，大多数只能采用数值解法(ji f)。 在在MATLAB中，使用中，使用ode函数建立微分方程模型。函数建立微分方程模型。7.2.2 范例分析范例分析例：在例：在RC低通滤波器电路中，电阻低通滤波器电路中，电阻R=5，理想电压源为，理想电压源为Vi=20V，电容，电容C=70F。分析电容元件的时域特性。分析电容元件的时域

9、特性。第九页，共三十六页。（1）分析：电容电压）分析：电容电压(diny)和电流的关和电流的关系系dttdVCtIdttICtVCC)()()(1)( 根据根据(gnj)基尔霍夫定律，可得出微分方程基尔霍夫定律，可得出微分方程dttdVCtIVtVtRICiC)()()()( 使用使用(shyng)ode函数时，对微分方程进行如下假函数时，对微分方程进行如下假设设)(tVyC RCyVdtdyi 第十页，共三十六页。（2）建立状态导数函数）建立状态导数函数function dy=cap(t,y)Vi=20;R=5;C=70e-6;dy=(Vi-y)/(R*C);（3）使用）使用(shyng)o

10、de函数进行仿真，仿真时间函数进行仿真，仿真时间00.006s，Vc初始值为初始值为0V。 t,y=ode45(cap,0,0.006,0); plot(t,y) axis(0 0.006 0 25) title(Vc-Time) xlabel(Time/sec) ylabel(Vc/V)第十一页，共三十六页。 当电压源为直流电压源时，加载在电容上的电压随时当电压源为直流电压源时，加载在电容上的电压随时间呈抛物线增大，稳态值为电源电压。电容电压在间呈抛物线增大，稳态值为电源电压。电容电压在t=0时时取得取得(qd)最小值，最小值为最小值，最小值为0；电容电压在；电容电压在t=0.0023时达到

11、时达到最大值，为最大值，为20V。第十二页，共三十六页。例：用简单的例：用简单的LC谐振电路组成滤波器电路，其电路方程是二阶谐振电路组成滤波器电路，其电路方程是二阶微分方程。观察该微分方程。观察该RLC电路中，相关电气元件的时域变化情况电路中，相关电气元件的时域变化情况(qngkung)。假设电源为直流电压源。假设电源为直流电压源Vi=20V，电阻，电阻R=5，电容，电容C=70F，电感，电感L=70mH。第十三页，共三十六页。（1）分析）分析(fnx)：根据电路分析：根据电路分析(fnx)，由基尔霍夫定律可以得出，由基尔霍夫定律可以得出微分方程微分方程CtIdttdVttILtVtVtVtV

12、tRICLiCL)()(d)(d)()()()()( 在利用在利用(lyng)ode函数时，对微分方程作出如下假设：函数时，对微分方程作出如下假设：)2(*)1()(1)2()2(1)1()()2()()1(yRytVLdtdyyCdtdytIytVyiC 第十四页，共三十六页。（2）创建）创建(chungjin)状态导数函数状态导数函数 function dy=RLC(t,y) Vi=20; R=5; C=70e-06; L=70e-03; dy=zeros(2,1); dy(1)=y(2)/C; dy(2)=(Vi-y(1)-R*y(2)/L;（3）使用）使用ode函数进行仿真，仿真时间函

13、数进行仿真，仿真时间00.12s，Vc初始值初始值为为0V，I初始值为初始值为0A。第十五页，共三十六页。t,y=ode45(RLC,0 0.12,0;0);figure(1) subplot(2,1,1); plot(t,y(:,1); title(Vc-Time) xlabel(Time/sec) ylabel(Vc/V) subplot(2,1,2); plot(t,y(:,2); title(I-Time) xlabel(Time/sec) ylabel(I/A) figure(2) plot(y(:,2),y(:,1) title(Vc-I) xlabel(I/A) ylabel(V

14、c/V)第十六页，共三十六页。 在电容在电容(dinrng)电压与时间曲线中，可以得出如下结电压与时间曲线中，可以得出如下结论：当电压源为直流电压源时，加载在电容论：当电压源为直流电压源时，加载在电容(dinrng)上的上的电压随时间呈震荡衰减；在电路电流与时间曲线中，可电压随时间呈震荡衰减；在电路电流与时间曲线中，可以得出结论，当电压源为直流电压源时，流经电路的电以得出结论，当电压源为直流电压源时，流经电路的电流随时间呈震荡衰减；在电容流随时间呈震荡衰减；在电容(dinrng)电压与电路电流的电压与电路电流的曲线图中可以得出如下结论：当处于稳态时，电容曲线图中可以得出如下结论：当处于稳态时，

15、电容(dinrng)电压为电压为20V，电路电流为，电路电流为0A。第十七页，共三十六页。7.3 传递函数模型传递函数模型 传递函数模型在一般控制系统中运用十分广泛，比传递函数模型在一般控制系统中运用十分广泛，比微分方程更加微分方程更加(gnji)方便使用。方便使用。传递函数是输出值拉普拉斯变传递函数是输出值拉普拉斯变换后的函数与输入值拉普拉斯变换后的函数之间的比值换后的函数与输入值拉普拉斯变换后的函数之间的比值。这主。这主要介绍拉普拉斯变换、传递函数的零点、极点和增益、以及传要介绍拉普拉斯变换、传递函数的零点、极点和增益、以及传递函数的部分分式展开。递函数的部分分式展开。7.3.1 拉普拉斯

16、变换拉普拉斯变换 jj -0d)(j21)(d)()( sesFtftetfsFstst F(s)是时域函数是时域函数(hnsh)f(t)在在s域的映像，因此称域的映像，因此称F(s)为为f(t)的象函数，函数的象函数，函数f(t)为为F(s)的原函数。的原函数。第十八页，共三十六页。1. 线性电路的线性电路的s域解法域解法(ji f)（1）电阻方程）电阻方程)()()()(sGVsIsRIsV （2）电容）电容(dinrng)方程方程)()0()()(1)0()(ssCVCVsIsIsCsVsV （3）电感）电感(din n)方程方程)(1)0()()()0()(sVsLsIsIssLILI

17、sV 第十九页，共三十六页。2. s域的传递函数域的传递函数 设电路的输入信号为设电路的输入信号为f(t)，对应的，对应的s域函数域函数(hnsh)F(s)。设。设电路的输出信号为电路的输出信号为g(t)，对应的，对应的s域函数域函数(hnsh)为为G(s)。则电。则电路在路在s域的传递函数域的传递函数(hnsh)为：为：nnnnmmmmasasasabsbsbsbsFsGsH 11101110)()()( 系统系统(xtng)在在MATLAB中可以方便地中可以方便地由分子和分母系数构由分子和分母系数构成的两个向量唯一地确定出来，成的两个向量唯一地确定出来，这两个向量分别用这两个向量分别用nu

18、m和和den表示。表示。 num=b0,b1,bm-1,bm den=a0,a1,an-1,an注意：注意：它们都是按它们都是按s的降幂进行排列的。的降幂进行排列的。第二十页，共三十六页。例：对下列例：对下列(xili)传递函数作幅频曲线和相频曲线。传递函数作幅频曲线和相频曲线。65431)(23 sssssH建立建立(jinl)文件文件control1.mnum=1 1;den=3 4 5 6;h,w=freqs(num,den);amp=abs(h);ang=angle(h);subplot(2,1,1);semilogx(w,amp);title(Amplitude-Frequency)

19、;xlabel(Frequency/rad);ylabel(Amplitude);subplot(2,1,2);semilogx(w,ang);title(Angle-Frequency);xlabel(Frequency/rad);ylabel(Angle);第二十一页，共三十六页。7.3.2 传递函数的零点传递函数的零点(ln din)、极点和增益、极点和增益1. 利用传递函数求解零点、极点和增益利用传递函数求解零点、极点和增益 零点、极点和增益是传递函数的三个特性。利用这零点、极点和增益是传递函数的三个特性。利用这些特性可以判断系统的稳定性。些特性可以判断系统的稳定性。MATLAB可以将

20、传递函可以将传递函数形式数形式nnnnmmmmasasasabsbsbsbsFsGsH 11101110)()()()()()()()()()(2121nmpspspszszszsksFsGsH 转换为零点、极点转换为零点、极点(jdin)和增益的形式：和增益的形式： 如果系统的所有极点都位于复平面的左半平面，则系统如果系统的所有极点都位于复平面的左半平面，则系统稳定；相反稳定；相反(xingfn)，系统是不稳定的。如果系统为稳定系统，系统是不稳定的。如果系统为稳定系统，且系统的所有零点都位于复平面的左半平面，则系统为最小且系统的所有零点都位于复平面的左半平面，则系统为最小系统；相反系统；相反

21、(xingfn)，系统为非最小系统。，系统为非最小系统。第二十二页，共三十六页。 在在MATLAB中，利用中，利用(lyng)tf2zp函数求解零点、极点和函数求解零点、极点和增益。增益。z,p,k=tf2zp(num,den) z：传递函数的零点矢量：传递函数的零点矢量 p：传递函数的极点矢量：传递函数的极点矢量 k：传递函数的增益：传递函数的增益 num：传递函数的分子系数矢量：传递函数的分子系数矢量 den：传递函数的分母系数矢量：传递函数的分母系数矢量例：根据传递函数求解零点、极点和增益，并判断该系统是例：根据传递函数求解零点、极点和增益，并判断该系统是否稳定。否稳定。98765434

22、32)(234562 sssssssssH第二十三页，共三十六页。建立建立(jinl)文件文件control2.mnum=2 3 4;den=3 4 5 6 7 8 9;z,p,k=tf2zp(num,den);disp(zeros:);zdisp(poles:);pdisp(gain3:);k control2zeros:z = -0.7500 + 1.1990i -0.7500 - 1.1990ipoles:p = 0.7099 + 0.9561i 0.7099 - 0.9561i -1.0903 + 0.5189i -1.0903 - 0.5189i -0.2863 + 1.1701i

23、-0.2863 - 1.1701igain3:k = 0.6667结论：从计算结果可以得知，结论：从计算结果可以得知，因为传递函数包括因为传递函数包括(boku)实部实部为正的极点，所以系统是不为正的极点，所以系统是不稳定的。稳定的。第二十四页，共三十六页。2. 根据极点、零点和增益根据极点、零点和增益(zngy)求解传递函数求解传递函数 该部分介绍将传递函数的零点、极点和增益的形式转换为该部分介绍将传递函数的零点、极点和增益的形式转换为传递函数形式。传递函数形式。 在在MATLAB软件中，利用软件中，利用zp2tf函数求解传递函数，函数求解传递函数，zp2tf函数的表达式为：函数的表达式为：

24、 num,den=zp2tf(z,p,k) num：传递函数的分子系数矢量：传递函数的分子系数矢量 den：传递函数的分母系数矢量：传递函数的分母系数矢量 z：传递函数的零点矢量（列向量）：传递函数的零点矢量（列向量） p：传递函数的极点矢量（列向量）：传递函数的极点矢量（列向量） k：传递函数的增益：传递函数的增益第二十五页，共三十六页。例：根据零点、极点和增益例：根据零点、极点和增益(zngy)求解传递函数求解传递函数)5)(4)(32)(32()1)(2(2)( ssisissssH编写编写(binxi)control3.m文件文件z=-1;-2;p=-2+3i;-2-3i;-4;-5;

25、k=2;num den=zp2tf(z,p,k);disp(num:);numdisp(den:);den control3num:num = 0 0 2 6 4den:den = 1 13 69 197 2602601976913462)(2342 sssssssH第二十六页，共三十六页。3. 绘制零点极点图绘制零点极点图 在在MATLAB软件中，可以根据传递函数的零、极点软件中，可以根据传递函数的零、极点和增益绘制零点极点图；也可以根据传递函数绘制零点和增益绘制零点极点图；也可以根据传递函数绘制零点极点图。所使用函数为极点图。所使用函数为zplane。（1）zplane(z,p) z：传递

26、函数的零点矢量：传递函数的零点矢量 p：传递函数的极点矢量：传递函数的极点矢量（2）zplane(num,den) num：传递函数的分子：传递函数的分子(fnz)系数矢量系数矢量 den：传递函数的分母系数矢量：传递函数的分母系数矢量 在零点极点图中，极点使用在零点极点图中，极点使用“”表示，零点使用表示，零点使用“o”表示。表示。第二十七页，共三十六页。例：根据传递函数绘制零点例：根据传递函数绘制零点(ln din)极点图极点图7809714981321834)(23452 ssssssssH编写编写(binxi)文件文件control4.mnum=1 4 3;den=1 18 132 4

27、98 971 780;figure(1)zplane(num,den);title(Zero-Pole Plane);xlabel(Real Part);ylabel(Imaginary Part);z,p,k=tf2zpk(num,den);disp(zeros:);zdisp(poles:);pdisp(gain3:);kfigure(2)zplane(z,p);title(Zero-Pole Plane);xlabel(Real Part);ylabel(Imaginary Part);第二十八页，共三十六页。 control4zeros:z = 0 0 0 -3 -1由此可以得到传递函

28、数零点、极点增益由此可以得到传递函数零点、极点增益(zngy)形式为：形式为：)4)(3)(23)(23)(5()1)(3()( ssisisssssHpoles:p = -5.0000 -3.0000 + 2.0000i -3.0000 - 2.0000i -4.0000 -3.0000 gain3:k = 1第二十九页，共三十六页。 从上可以看出，利用传递函数和极点零点可以画出从上可以看出，利用传递函数和极点零点可以画出相同的零点极点图。从图中可以看出，系统相同的零点极点图。从图中可以看出，系统(xtng)的极点的极点全部位于系统全部位于系统(xtng)的左半平面，因此该系统的左半平面，因

29、此该系统(xtng)稳定。稳定。第三十页，共三十六页。7.3.3 传递函数的部分分式展开传递函数的部分分式展开1. 根据传递函数求解部分分式展开根据传递函数求解部分分式展开 部分分式展开是传递函数的一种形式部分分式展开是传递函数的一种形式(xngsh)，可以将，可以将传递函数形式传递函数形式(xngsh)变换为部分分式展开形式变换为部分分式展开形式(xngsh)。)()()()(2211skpsrpsrpsrsFsGsHnn 式中，式中，r为系统的留数矢量，为系统的留数矢量，p为系统的极点矢量，为系统的极点矢量，k(s)为系统的常数矢量。在为系统的常数矢量。在MATLAB软件中，利用软件中，利

30、用(lyng)residue函函数求解部分分式展开。数求解部分分式展开。 r,p,k=residue(num,den) r：部分分式的留数：部分分式的留数 k：部分分式的常数矢量：部分分式的常数矢量第三十一页，共三十六页。例：根据传递函数求解例：根据传递函数求解(qi ji)部分分式展开部分分式展开9876543432)(234562 sssssssssH建立建立(jinl)文件文件control5.mnum=2 3 4;den=3 4 5 6 7 8 9;r,p,k=residue(num,den);disp(residue:);rdisp(poles:);pdisp(direct term

31、s:)k第三十二页，共三十六页。 control5residue:r = 0.0600 - 0.0831i 0.0600 + 0.0831i 0.0520 - 0.0454i 0.0520 + 0.0454i -0.1120 - 0.0665i -0.1120 + 0.0665i根据所得结果，便可以根据所得结果，便可以(ky)写出部分分式展开。写出部分分式展开。poles:p = -1.0903 + 0.5189i -1.0903 - 0.5189i -0.2863 + 1.1701i -0.2863 - 1.1701i 0.7099 + 0.9561i 0.7099 - 0.9561idir

32、ect terms:k = )9561070990(0665011200)9561070990(0665011200)1701128630(0454005200)1701128630(0454005200)5189009031(0831006000)5189009031(0831006000)(i. - .s-i. .-i. .s-i. - .-i. - .-s-i. .i. .-s-i. - .i. - .-s-i. .i.si.sH 第三十三页，共三十六页。2. 根据部分根据部分(b fen)分式求解传递函数分式求解传递函数 MATLAB软件可以将部分分式展开的形式转换为传软件可以将部分分

33、式展开的形式转换为传递函数形式，使用的函数为递函数形式，使用的函数为residue，表达式为：，表达式为： num,den=residue(r,p,k) r：部分分式的留数矢量（列向量）：部分分式的留数矢量（列向量） p：部分分式的极点矢量（列向量）：部分分式的极点矢量（列向量） k：部分分式的常数矢量（行向量）：部分分式的常数矢量（行向量） num：传递函数分子矢量：传递函数分子矢量 den：传递函数分母矢量：传递函数分母矢量例：根据部分分式展开求解传递函数例：根据部分分式展开求解传递函数1)1(2)1(2332211)( sisiisissssH第三十四页，共三十六页。编写编写(binxi

34、)文件文件control6.mr=1;2;3;2i;-2i;p=1;2;3;1+i;1-i;k=1 1;num,den=residue(r,p,k);disp(num:);numdisp(den:);den在命令窗口在命令窗口(chungku)中执行该文件中执行该文件 control6num:num = 1 -7 23 -53 92 -102 48den:den = 1 -8 25 -40 34 -12 根据以上输出的传递函数的分子系数矢量和分母根据以上输出的传递函数的分子系数矢量和分母(fnm)系数矢量，得出以下传递函数：系数矢量，得出以下传递函数：123440258481029253237)(234523456 ssssssssssssH第三十五页，共三十六页。内容(nirng)总结第7章 控制系统仿真。odefun：定义微分方程的形式y=f(t,y)。值）到tfinal（终值），该积分限也。 y1=y(:,1)。 y2=y(:,2)。电容电压在t=0时取得最小值，最小值为0。（1）分析：根据电路分析，由基尔霍夫定律可以得出(d ch)微分方程。系统在MATLAB中可以方便地由分子和分母系数构成的两个向量唯一地确定出来，这两个向量分别用num和den表示。den =第三十六页，共三十六页。