连续系统数值积分法仿真Matlab编程公开课一等奖优质课大赛微课获奖课件.pptx

第三章第三章 Matlab编程编程(数值积分法仿真数值积分法仿真)3.1 连续系统数值积分法仿真编程思绪连续系统数值积分法仿真编程思绪目的：针对下面系统编写通用数值积分法仿真程序(3.1-1)对这样系统进行仿真，事实上涉及到控制计算、状态数值积分计算和输出计算。3个函数g,f,h拟定后，就能够完整地描述一个系统。给定初值：u0，y0,系统中向量都采用列向量表示1第1页第1页function t,y=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,InteMethod,StateModel,OutputFile,ControlFile)%函数功效：用数值积分法仿真%输入参数：tstart,tstop,h 分别是起始时间、结束时间和仿真步长，是标量%x0,u0是状态、输入初值，都是列向量%cnty是输出变量个数%InteMethod时数值积分办法，能够是EUL,RK2,RK4%StateModel是状态模型文献名%ControlFile是控制作用文献名%OutputFile是系统输出文献名%输出参数：t是仿真结果时间序列%y是仿真结果系统输出序列(1)数值积分法仿真框架函数数值积分法仿真框架函数2第2页第2页function t,y=w_DigiInteSimu(tstart,tstop,h,x0,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,cury);%计算控制时传递参数：当前时间，步长，当前状态和输出 curx=w_StepIntegral(t(i),h,curx,curu,InteMethod,StateModel);%单步积分运算 cury=eval(OutputFile,(t(i),curx,curu);%计算输出 y(:,i+1)=cury;%将输出加入到输出序列里end3第3页第3页function NewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel)%函数功效：单步积分运算,与模型方程相关%输入参数：curt是当前时间，h是数值积分步长%curx,curu分别是当前状态和控制向量%InteMethod是积分算法，字符串类型，能够取EUL,RK2,RK4%StateModel是状态模型文献名称，字符串类型%输出参数：NewX是这一步计算新状态向量(2)单步数值积分法函数单步数值积分法函数单步数值积分函数只是对微分方程组StateModel进行一步计算，计算法由InteMethod参数指定，能够上欧拉法，RK2或RK4。4第4页第4页function NewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel);if InteMethod=RK4 k1=eval(StateModel,(curt,curx,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.5*h*(k1+k2);else%欧拉法EUL Qk=eval(StateModel,(curt,curx,curu);%eval用来执行一个函数，传递函数名和函数输入参数 NewX=curx+h*Qk;end5第5页第5页函数w_DigiInteSimu和w_StepIntegral结构了一个数值积分法仿真框架，并不涉及详细系统。详细系统由StateModel,ControlFile,OutputFile参个参数决定，事实上就是三个函数文献名，这三个函数输入输出参数必须遵循特定格式，在准备好由这3个函数描述系统后，调用w_DigiInteSimu即可进行仿真。还需要一个调用w_DigiInteSimu进行仿真程序，它指定模型文献，指定初始参数，并且对仿真结果绘图。6第6页第6页3.2 算法稳定性分析仿真编程算法稳定性分析仿真编程针对下面系统，求用解析法、欧拉法和RK4分别求解，计算欧拉法最大允许步长，将步长从0.1逐步增大，比较三种解效果。解：1）该系统是稳定，解析解为2）用欧拉法计算本例时，其步长应当满足3）RK4步长条件式是一个高阶不等式，无法直接求解，只能用试探法拟定RK4步长上限。(3.2-1)7第7页第7页我们使用3.1简介两个基本函数，同时要针对系统模型编写3个函数来描述该系统，最后编写一个实现仿真初值设置、仿真求解和仿真结果显示函数。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)式状态微分方程8第8页第8页function OutputY=w_SysOutput(curt,curx,curu)%function OutputY=w_SysOutput(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是数值积分步长 curx,curu,cury分别是当前状态、控制和输出向量%输出参数：ControlU是计算出来控制向量ControlU=curu;%依据实际系统写控制（3）系统控制计算函数）系统控制计算函数w_SysControl9第9页第9页（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_DigiInteSimu(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和u2)设置模型输出方程、状态方程和控制方程所在函数文献名3)调用w_DigiInteSimu进行仿真计算4)将计算结果绘图进行比较。10第10页第10页修改simu_Stability函数里参数h，就能够作出不同时长时仿真结果。下面是步长从0.1不停增大时仿真结果h=0.1时仿真曲线。1】RK4法曲线与解析解重叠，阐明RK4法非常准确。2】欧拉法有较大误差。11第11页第11页h=0.4时仿真曲线。1】RK4法还是比较准确。2】欧拉法误差很大，出现衰减震荡。3】两种数值算法步长条件仍然满足，算法仍然稳定。12第12页第12页h=0.6时仿真曲线。欧拉法步长条件已经不满足，仿真解发散h=0.6时仿真曲线。RK4法步长条件仍然满足。但误差增大13第13页第13页h=0.7时仿真曲线。RK4法步长条件已经不满足，仿真结果发散。14第14页第14页3.3 卫星发射仿真卫星发射仿真卫星发射运动方程G为重力系数系统是高阶微分方程，不能直接用于仿真计算，需要转化为状态方程。定义状态变量：运动方程用极坐标表示，同时还需要用直角坐标输出。1.模型转换模型转换15第15页第15页得到系统仿真模型取X-Y平面坐标作为输出状态初值为：v是初始发射速度，可分别取该系统没有控制，主要是研究初始发射速度对轨道影响。16第16页第16页2.模型描述函数结构模型描述函数结构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;%转换为列向量(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节同样控制函数。17第17页第17页3.仿真组织函数仿真组织函数function simu_Satelite%函数功效：对卫星发射轨道仿真tstart=0;StateModel=sat_StateEqu;ControlFile=w_SysControl;OutputFile=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求解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(卫星发射轨道);18第18页第18页直角坐标显示卫星发射轨道：初始速度越大，规大越大。V不同时，需要步长和计算步数不同。19第19页第19页3.4 线性系统仿真编程线性系统仿真编程前面所简介仿真程序是针对(3.1-1)所表示通用非线性系统，针对下列所表示线性状态空间模型(3.4-1)我们即使也能够用前面简介程序仿真，但比较麻烦，事实上系统由矩阵A、B、C、D就能够指定微分方程和输出方程，而不必单独用函数文献来表示。为此，我们单独编写针对系统(3.4-1)线性系统数值积分法仿真程序。程序仍然由仿真框架和单步积分构成。20第20页第20页(1)线性系统数值积分仿真框架线性系统数值积分仿真框架function 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是仿真结果系统输出序列21第21页第21页function 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;%当前一步xcuru=u0;%当前一步ucury=y0;%当前一步yfor i=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;%将输出加入到输出序列里end22第22页第22页(2)线性系统单步积分函数线性系统单步积分函数function NewX=w_LinearStepIntegral(h,curx,curu,InteMethod,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+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);end23第23页第23页3.5二阶系统传递函数仿真二阶系统传递函数仿真1.模型转换模型转换写为状态模型就是基本参数取值为针对本例，线性系统仿真程序进行仿真。24第24页第24页2.仿真组织程序仿真组织程序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,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(二阶系统阶跃输入响应);25第25页第25页二阶系统在阶跃输入作用下，不同阻尼比时输出曲线。26第26页第26页3.6面向结构图仿真面向结构图仿真对下图系统进行面向结构图变换，并进行仿真。G1(s)G2(s)G3(s)G4(s)其中27第27页第27页(1)模型处理模型处理每个环节用3个参数来描述，得到矩阵28第28页第28页(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;%单位阶跃输入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(面向结构图仿真);29第29页第29页G=2,1,1;0,1,3;5,1,6;0,1,1;只绘制y1和y3，能够看出两个输出区别，它们都趋向于稳定。30第30页第30页同时绘制4个输出时，由于y2,y4改变很大，发散，因此y1和y3看起来仿佛重叠了。31第31页第31页附录：plot函数参数取值意义Plot函数普通调用是：plot(x1,y1,option1,x2,y2,option2,.)选项意义选项意义-实线.点号数据点-虚线*星号数据点-.点划线+加号数据点:点线y黄色g绿色m洋红色b蓝色c青色w白色r红色k黑色32第32页第32页