第3-1章-连续系统数值积分法仿真Matlab编程课件.ppt

1、1第三章第三章 Matlab编程编程(数值积分法仿真数值积分法仿真)3.1 连续系统数值积分法仿真编程思路连续系统数值积分法仿真编程思路目的：针对下面的系统编写通用的数值积分法仿真程序),(),(),(uxthyuxtfxyuxtgu(3.1-1)对这样的系统进行仿真，实际上涉及到控制的计算、状态的数值积分计算和输出的计算。3个函数g,f,h确定后，就可以完整地描述一个系统。给定初值：u0，y0,系统中的向量都采用列向量表示2function t,y=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,InteMethod,StateModel,OutputFil

2、e,ControlFile)%函数功能：用数值积分法仿真%输入参数：tstart,tstop,h 分别是起始时间、结束时间和仿真步长，是标量%x0,u0是状态、输入的初值，都是列向量%cnty是输出变量的个数%InteMethod时数值积分方法，可以是EUL,RK2,RK4%StateModel是状态模型的文件名%ControlFile是控制作用的文件名%OutputFile是系统输出的文件名%输出参数：t是仿真结果的时间序列%y是仿真结果系统的输出序列(1)数值积分法仿真框架函数数值积分法仿真框架函数3function t,y=w_DigiInteSimu(tstart,tstop,h,x0

3、,u0,cnty,InteMethod,StateModel,OutputFile,ControlFile)t=tstart:h:tstop;%t数一个行序列cntt=size(t,2);%返回列数y=zeros(cnty,cntt);%构造一个空矩阵，用来存储结果y0=eval(OutputFile,(tstart,x0,u0);%计算初始输出y(:,1)=y0;%将cury作为输出的第1列curx=x0;%当前一步的xcuru=u0;%当前一步的ucury=y0;%当前一步的yfor i=1:1:cntt-1 curu=eval(ControlFile,(t(i),h,curx,curu,

4、cury);%计算控制时传递的参数：当前时间，步长，当前状态和输出 curx=w_StepIntegral(t(i),h,curx,curu,InteMethod,StateModel);%单步积分运算 cury=eval(OutputFile,(t(i),curx,curu);%计算输出 y(:,i+1)=cury;%将输出加入到输出序列里end4function NewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel)%函数功能：单步积分运算,与模型方程有关%输入参数：curt是当前时间，h是数值积分步长%curx,curu分别

5、是当前的状态和控制向量%InteMethod是积分算法，字符串类型，可以取EUL,RK2,RK4%StateModel是状态模型文件名称，字符串类型%输出参数：NewX是这一步计算的新的状态向量(2)单步数值积分法函数单步数值积分法函数单步数值积分函数只是对微分方程组StateModel进行一步的计算，计算法由InteMethod参数指定，可以上欧拉法，RK2或RK4。5function NewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel);if InteMethod=RK4 k1=eval(StateModel,(curt,c

6、urx,curu);k2=eval(StateModel,(curt+0.5*h,curx+0.5*h*k1,curu);k3=eval(StateModel,(curt+0.5*h,curx+0.5*h*k2,curu);k4=eval(StateModel,(curt+h,curx+h*k3,curu);NewX=curx+h*(k1+2*k2+2*k3+k4)/6;elseif InteMethod=RK2 k1=eval(StateModel,(curt,curx,curu);k2=eval(StateModel,(curt+h,curx+h*k1,curu);NewX=curx+0.

7、5*h*(k1+k2);else%欧拉法EUL Qk=eval(StateModel,(curt,curx,curu);%eval用来执行一个函数，传递函数名和函数的输入参数 NewX=curx+h*Qk;end6函数w_DigiInteSimu和w_StepIntegral构造了一个数值积分法仿真的框架，并不涉及具体的系统。具体的系统由StateModel,ControlFile,OutputFile参个参数决定，实际上就是三个函数文件名，这三个函数输入输出参数必须遵循特定的格式，在准备好由这3个函数描述的系统后，调用w_DigiInteSimu即可进行仿真。还需要一个调用w_DigiInt

8、eSimu进行仿真的程序，它指定模型文件，指定初始参数，并且对仿真结果绘图。73.2 算法稳定性分析仿真编程算法稳定性分析仿真编程针对下面的系统，求用解析法、欧拉法和RK4分别求解，计算欧拉法最大允许步长，将步长从0.1逐渐增大，比较三种解的效果。1)0(4xxx 解：1）该系统是稳定的，解析解为tetx4)(2）用欧拉法计算本例时，其步长应该满足5.0141hh3）RK4步长条件式是一个高阶不等式，无法直接求解，只能用试探法确定RK4的步长上限。(3.2-1)8我们使用3.1介绍的两个基本函数，同时要针对系统模型编写3个函数来描述该系统，最后编写一个实现仿真初值设置、仿真求解和仿真结果显示的

9、函数。function derX=w_SysStateEqu(curt,curx,curu)%function derX=w_stateEqu(curt,curx,curu)%函数功能：在此函数里写系统的状态方程%输入参数：curt当前时间，curx和curu是当前状态和控制%输出参数：derX是状态的X的导数值derX(1)=-4*curx(1);（1）状态模型函数）状态模型函数w_SysStateEqu在该函数里描述(3.2-1)式的状态微分方程9function OutputY=w_SysOutput(curt,curx,curu)%function OutputY=w_SysOutpu

10、t(curt,curx,curu)%函数功能：计算系统的输出向量%输入参数：curt是当前时间，curx,curu分别是当前的状态和控制向量%输出参数：OutputY是计算出来的输出向量OutputY(1)=curx(1);%根据具体的模型写输出方程（2）系统输出函数）系统输出函数w_SysOutputfunction ControlU=w_SysControl(curt,h,curx,curu,cury)%ControlU=w_sysControl(curt,h,curx,curu,cury)%函数功能：计算系统的控制，如u=PID(t,curx,cury)%输入参数：curt是当前时间，h

11、是数值积分步长 curx,curu,cury分别是当前的状态、控制和输出向量%输出参数：ControlU是计算出来的控制向量ControlU=curu;%根据实际系统写控制（3）系统控制计算函数）系统控制计算函数w_SysControl10（4）仿真组织函数）仿真组织函数function simu_Stability%函数功能：稳定性分析仿真tstart=0;tstop=7;h=0.7;StateModel=w_SysStateEqu;ControlFile=w_SysControl;OutputFile=w_SysOutput;x0=1;u0=0;cnty=1;t,y1=w_DigiInte

12、Simu(tstart,tstop,h,x0,u0,cnty,EUL,StateModel,OutputFile,ControlFile);%用欧拉法求解t,y4=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,RK4,StateModel,OutputFile,ControlFile);%用RK4求解y5=exp(-4*t);%计算解析解plot(t,y1,-m,t,y2,-.r,t,y5,-k);%绘制叠加曲线xlabel(时间);ylabel(y);title(算法稳定性分析仿真);该函数的主要任务是：1)设置仿真起始时间、结束时间、仿真步长、初始的x和

13、u2)设置模型的输出方程、状态方程和控制方程所在函数文件名3)调用w_DigiInteSimu进行仿真计算4)将计算结果绘图进行比较。11修改simu_Stability函数里的参数h，就可以作出不同步长时的仿真结果。下面是步长从0.1不断增大时的仿真结果h=0.1时的仿真曲线。1】RK4法的曲线与解析解重合，说明RK4法非常准确。2】欧拉法有较大误差。00.20.40.60.811.21.41.61.8200.10.20.30.40.50.60.70.80.91时 间y算 法 稳 定 性 分 析 仿 真EulerRK4解 析12h=0.4时的仿真曲线。1】RK4法还是比较准确。2】欧拉法误差

14、很大，出现衰减震荡。3】两种数值算法的步长条件仍然满足，算法仍然稳定。00.511.522.533.54-0.6-0.4-0.200.20.40.60.81时 间y算 法 稳 定 性 分 析 仿 真EulerRK4解 析13h=0.6时的仿真曲线。欧拉法的步长条件已经不满足，仿真解发散h=0.6时的仿真曲线。RK4法的步长条件仍然满足。但误差增大0123456-30-20-100102030时 间y算 法 稳 定 性 分 析 仿 真Euler解 析012345600.20.40.60.81时 间y算 法 稳 定 性 分 析 仿 真RK4解 析14h=0.7时的仿真曲线。RK4法的步长条件已经不

15、满足，仿真结果发散。0123456700.511.52时 间y算 法 稳 定 性 分 析 仿 真Euler解 析153.3 卫星发射仿真卫星发射仿真卫星发射运动方程 dtddtdrrdtddtdrGdtrdr2222222sec/4014083kmG G为重力系数系统是高阶微分方程，不能直接用于仿真计算，需要转化为状态方程。定义状态变量：4321,xrxxrx运动方程用极坐标表示，同时还需要用直角坐标输出。1.模型转换模型转换16143424121342312xxxxxxxGxxxxx得到系统仿真模型取X-Y平面坐标作为输出)sin()cos(212211xxyxxy状态初值为：)0()0(0

16、)0(0)0(6400)0(14321xvxxxkmxv是初始发射速度，可分别取sec/12sec,/10sec,/8kmkmkmv 该系统没有控制，主要是研究初始发射速度对轨道的影响。172.模型描述函数构造模型描述函数构造function derX=sat_StateEqu(curt,curx,curu)G=401408;derX(1)=curx(3);derX(2)=curx(4);derX(3)=-G/(curx(1)*curx(1)+curx(1)*curx(4)*curx(4);derX(4)=-2*curx(3)*curx(4)/curx(1);derX=derX;%转换为列向量

17、(1)状态方程状态方程function OutputY=sat_Output(curt,curx,curu)OutputY(1)=curx(1)*cos(curx(2);OutputY(2)=curx(1)*sin(curx(2);OutputY=OutputY;%转换为列向量(2)输出方程输出方程(3)控制方程控制方程由于该系统没有控制，因而采用与3.2节同样的控制函数。183.仿真组织函数仿真组织函数function simu_Satelite%函数功能：对卫星发射轨道仿真tstart=0;StateModel=sat_StateEqu;ControlFile=w_SysControl;O

18、utputFile=sat_Output;u0=0;cnty=2;h=200;tstop=100*h;v=10;x0=6400,0,0,v/6400;t,y10=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,RK4,StateModel,OutputFile,ControlFile);%用RK4求解h=200;tstop=50*h;v=9;x0=6400,0,0,v/6400;t,y9=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,RK4,StateModel,OutputFile,ControlFile);%用RK4求解

19、h=100;tstop=60*h;v=8;x0=6400,0,0,v/6400;t,y8=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,RK4,StateModel,OutputFile,ControlFile);%用RK4求解plot(y10(1,:),y10(2,:),-,y9(1,:),y9(2,:),-,y8(1,:),y8(2,:),:);%绘制叠加曲线xlabel(x);ylabel(y);title(卫星发射轨道);19直角坐标显示的卫星发射轨道：初始速度越大，规大越大。V不同时，需要的步长和计算步数不一样。-3-2.5-2-1.5-1-0.5

20、00.51x 104-1.5-1-0.500.511.5x 104xy卫 星 发 射 轨 道v=10v=9v=8203.4 线性系统仿真编程线性系统仿真编程前面所介绍的仿真程序是针对(3.1-1)所表示的通用的非线性系统，针对如下所表示的线性状态空间模型DuCxyBuAxx(3.4-1)我们虽然也可以用前面介绍的程序仿真，但比较麻烦，实际上系统由矩阵A、B、C、D就可以指定微分方程和输出方程，而不必单独用函数文件来表示。为此，我们单独编写针对系统(3.4-1)的线性系统数值积分法仿真程序。程序仍然由仿真框架和单步积分构成。21(1)线性系统数值积分仿真框架线性系统数值积分仿真框架functio

21、n t,y=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,InteMethod,ControlFile)%函数功能：用数值积分法仿真对线性系统dx=AX+Bu,y=Cx+Du进行仿真%输入参数：tstart,tstop,h 分别是起始时间、结束时间和仿真步长，是标量%x0,u0是状态、输入的初值，都是列向量%A,B,C,D是线性状态空间模型的系数矩阵%InteMethod是数值积分方法，可以是EUL,RK2,RK4%ControlFile是控制作用函数，没有控制算法时，可以省略%输出参数：t是仿真结果的时间序列%y是仿真结果系统的输出序列22functio

22、n t,y=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,InteMethod,ControlFile)t=tstart:h:tstop;%t数一个行序列cntt=size(t,2);%返回列数cnty=size(C,1);%返回y的维数y=zeros(cnty,cntt);%构造一个空矩阵，用来存储结果y0=C*x0+D*u0;%eval(OutputFile,(tstart,x0,u0);%计算初始输出y(:,1)=y0 ;%将y0作为输出的第1列curx=x0;%当前一步的x curu=u0;%当前一步的u cury=y0;%当前一步的yfor i=

23、1:1:cntt-1 curu=eval(ControlFile,(t(i),h,curx,curu,cury,A,B,C,D);%计算控制时传递的参数：当前时间，步长，当前状态和输出,模型系数矩阵 curx=w_LinearStepIntegral(h,curx,curu,InteMethod,A,B);%单步积分运算,针对线性模型 cury=C*curx+D*curu;%计算输出 y(:,i+1)=cury;%将输出加入到输出序列里end23(2)线性系统单步积分函数线性系统单步积分函数function NewX=w_LinearStepIntegral(h,curx,curu,InteM

24、ethod,A,B);%函数功能：单步积分运算,与模型方程有关%输入参数：h是数值积分步长，curx,curu分别是当前的状态和控制向量%InteMethod是积分算法，字符串类型，%可以取EUL,RK2,RK4,由于要用于判断，字符串必须等长%A,B是线性状态模型的微分方程的系数矩阵%输出参数：NewX是这一步计算的新的状态向量Bu=B*curu;if InteMethod=RK4 k1=A*curx+Bu;k2=A*(curx+0.5*h*k1)+Bu;k3=A*(curx+0.5*h*k2)+Bu;k4=A*(curx+h*k3)+Bu;NewX=curx+h*(k1+2*k2+2*k3

25、+k4)/6;elseif InteMethod=RK2 k1=A*curx+Bu;%eval(StateModel,(curt,curx,curu);k2=A*(curx+h*k1)+Bu;%eval(StateModel,(curt+h,curx+h*k1,curu);NewX=curx+0.5*h*(k1+k2);else%欧拉法EUL NewX=curx+h*(A*curx+Bu);end243.5二阶系统传递函数仿真二阶系统传递函数仿真2222)()()(nnnsssusysG1.模型转换模型转换写为状态模型就是Tnnnxxyuxxxx212212210102102.0,1,5.0,

26、1hun基本参数取值为针对本例，线性系统仿真程序进行仿真。252.仿真组织程序仿真组织程序function simu_StateSpace%函数功能：对一个状态空间模型进行仿真tstart=0;tstop=20;h=0.1;B=0;1;C=1,0;D=0;x0=0;0;u0=1;%单位阶跃输入dampratio=0.1;A=0,1;-1,-2*dampratio;t,y1=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,RK4);dampratio=0.5;A(2,2)=-2*dampratio;t,y5=w_LinearSimu(tstart,tstop,

27、h,x0,u0,A,B,C,D,RK4);dampratio=1;A(2,2)=-2*dampratio;t,y10=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,RK4);dampratio=1.5;A(2,2)=-2*dampratio;t,y15=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,RK4);plot(t,y1,-r,t,y5,:g,t,y10,*y,t,y15,.b);%绘制叠加曲线xlabel(x);ylabel(y);title(二阶系统阶跃输入响应);260246810121416182000.

28、20.40.60.811.21.41.61.8xy二 阶 系 统 阶 跃 输 入 响 应=0.1=0.5=1.0=1.5二阶系统在阶跃输入作用下，不同阻尼比时的输出曲线。273.6面向结构图仿真面向结构图仿真对下图的系统进行面向结构图的变换，并进行仿真。G1(s)G2(s)G3(s)G4(s)0y1u1y2u2y3u3y4u4y其中2,1)(,56)(,3)(,21)(4321ssGssGssGssG28(1)模型处理模型处理0004321432104423312011001000010100010000yWWYyyyyyuuuuyuyyuyyuyu每个环节用3个参数来描述，得到的矩阵1106

29、15310112G29(2)仿真组织程序仿真组织程序利用第2章编写的结构图到状态空间模型转换的函数w_bd2ss进行模型处理function simu_BlockDiagram%函数功能：面向结构图的仿真G=2,1,1;0,1,3;5,1,6;0,1,1;W=0,0,0,0;1,0,-2,0;0,1,0,1;0,0,0,0;W0=1;0;0;1;P,Q=w_bd2ss(G,W,W0);%先通过面向结构图的模型变换得到状态模型tstart=0;tstop=5;h=0.1;C=1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1;%44个状态都输出D=0;x0=0;0;0;0;u0=1;

30、%单位阶跃输入t,y=w_LinearSimu(tstart,tstop,h,x0,u0,P,Q,C,D,RK4);%size(y)%plot(t,y(1,:),-r,t,y(2,:),:g,t,y(3,:),*m,t,y(4,:),.b);%绘制叠加曲线plot(t,y(1,:),-r,t,y(3,:),*m);%绘制叠加曲线xlabel(time);ylabel(y);title(面向结构图的仿真);30G=2,1,1;0,1,3;5,1,6;0,1,1;00.511.522.533.544.5500.10.20.30.40.5timey面 向 结 构 图 的 仿 真y1y3只绘制y1和y3，可以看出两个输出的区别，它们都趋向于稳定。3100.511.522.533.544.55-505timey面 向 结 构 图 的 仿 真y1y2y3y4同时绘制4个输出时，由于y2,y4变化很大，发散，所以y1和y3看起来好像重叠了。32附录：plot函数参数取值意义Plot函数的一般调用是：plot(x1,y1,option1,x2,y2,option2,.)