1、1,第三章 Matlab编程(数值积分法仿真),3.1 连续系统数值积分法仿真编程思路,目的：针对下面的系统编写通用的数值积分法仿真程序,(3.1-1),对这样的系统进行仿真，实际上涉及到控制的计算、状态的数值积分计算和输出的计算。3个函数g,f,h确定后，就可以完整地描述一个系统。,给定初值：u0，y0,系统中的向量都采用列向量表示,2,function t,y=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,InteMethod,StateModel,OutputFile, ControlFile) % 函数功能：用数值积分法仿真 % 输入参数：tstar
2、t, tstop,h 分别是起始时间、结束时间和仿真步长，是标量 % x0,u0是状态、输入的初值，都是列向量 % cnty是输出变量的个数 % InteMethod时数值积分方法，可以是EUL ,RK2,RK4 % StateModel是状态模型的文件名 % ControlFile是控制作用的文件名 % OutputFile是系统输出的文件名 % 输出参数：t是仿真结果的时间序列 % y是仿真结果系统的输出序列,(1)数值积分法仿真框架函数,3,function t,y=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,InteMethod,StateMode
3、l,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; %当前一步的x curu=u0; %当前一步的u cury=y0; %当前一步的y for i=1:1:cntt-1curu=eval(ControlFile,(t(i),h,curx,curu,cury);%计算控制时传递的参数：
4、当前时间，步长，当前状态和输出curx=w_StepIntegral(t(i),h,curx,curu,InteMethod,StateModel);%单步积分运算cury=eval(OutputFile,(t(i),curx,curu);%计算输出y(:,i+1)=cury;%将输出加入到输出序列里 end,4,function 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。,5,function NewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel); if InteMethod = RK4 k1=eval(StateModel,(curt,curx,curu)
6、;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 = RK2k1=eval(StateModel,(curt,curx,curu);k2=eval(StateModel,(curt+h,curx+h*k1,curu);NewX=curx+0.5*h*(k1
7、+k2); else %欧拉法EULQk=eval(StateModel,(curt,curx,curu);%eval用来执行一个函数，传递函数名和函数的输入参数NewX=curx+h*Qk; end,6,函数w_DigiInteSimu和w_StepIntegral构造了一个数值积分法仿真的框架，并不涉及具体的系统。 具体的系统由StateModel,ControlFile,OutputFile参个参数决定，实际上就是三个函数文件名，这三个函数输入输出参数必须遵循特定的格式，在准备好由这3个函数描述的系统后，调用w_DigiInteSimu即可进行仿真。 还需要一个调用w_DigiInteS
8、imu进行仿真的程序，它指定模型文件，指定初始参数，并且对仿真结果绘图。,7,3.2 算法稳定性分析仿真编程,针对下面的系统，求用解析法、欧拉法和RK4分别求解，计算欧拉法最大允许步长，将步长从0.1逐渐增大，比较三种解的效果。,解：1）该系统是稳定的，解析解为,2）用欧拉法计算本例时，其步长应该满足,3）RK4步长条件式是一个高阶不等式，无法直接求解，只能用试探法确定RK4的步长上限。,(3.2-1),8,我们使用3.1介绍的两个基本函数，同时要针对系统模型编写3个函数来描述该系统，最后编写一个实现仿真初值设置、仿真求解和仿真结果显示的函数。,function derX=w_SysState
9、Equ(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)式的状态微分方程,9,function OutputY=w_SysOutput(curt,curx,curu) % function OutputY=w_SysOutput(curt,curx,curu) %
10、 函数功能：计算系统的输出向量 % 输入参数：curt是当前时间，curx,curu分别是当前的状态和控制向量 % 输出参数：OutputY是计算出来的输出向量 OutputY(1)=curx(1);%根据具体的模型写输出方程,（2）系统输出函数w_SysOutput,function ControlU=w_SysControl(curt,h,curx,curu,cury) % ControlU=w_sysControl(curt,h,curx,curu,cury) % 函数功能： 计算系统的控制，如u=PID(t,curx,cury) % 输入参数：curt是当前时间，h是数值积分步长 cu
11、rx,curu,cury分别是当前的状态、控制和输出向量 % 输出参数：ControlU是计算出来的控制向量 ControlU=curu; %根据实际系统写控制,（3）系统控制计算函数w_SysControl,10,（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(tstar
12、t,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、和u 2)设置模型的输出方程、状态方程和控制方程所在函数文件名 3)调用w_DigiInteSimu进行仿真计算 4)将计算结果绘图进行比较。,11,修改simu_Stability函数里的参数h，就可以作出不同步长时的仿真结果。下面是步长从0.1不断增大时的仿真结果,h=0.1时的仿真曲线。 1】RK4法的曲线与解析解重合，说明RK4法非常准确。 2】欧拉法有较大误差。,12,h=0.4时的仿真曲线。 1】RK4法还是比较准确。 2】欧拉法误差很大，出现衰减震荡。 3】两种数值算法的步长条件仍然满足，算法仍然稳定。,13,h=0.6时的仿真曲线。欧拉法的步长条件已经不满足，仿真解发散,h=0
14、.6时的仿真曲线。RK4法的步长条件仍然满足。但误差增大,14,h=0.7时的仿真曲线。 RK4法的步长条件已经不满足，仿真结果发散。,15,3.3 卫星发射仿真,卫星发射运动方程,G为重力系数,系统是高阶微分方程，不能直接用于仿真计算，需要转化为状态方程。定义状态变量：,运动方程用极坐标表示，同时还需要用直角坐标输出。,1. 模型转换,16,得到系统仿真模型,取X-Y平面坐标作为输出,状态初值为：,v是初始发射速度，可分别取,该系统没有控制，主要是研究初始发射速度对轨道的影响。,17,2. 模型描述函数构造,function derX=sat_StateEqu(curt,curx,curu)
15、 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)输出方程
16、,(3)控制方程 由于该系统没有控制，因而采用与3.2节同样的控制函数。,18,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, OutputF
17、ile, 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求解
18、plot(y10(1,:),y10(2,:),-,y9(1,:),y9(2,:),-,y8(1,:),y8(2,:),:);%绘制叠加曲线 xlabel(x); ylabel(y); title(卫星发射轨道);,19,直角坐标显示的卫星发射轨道：初始速度越大，规大越大。 V不同时，需要的步长和计算步数不一样。,20,3.4 线性系统仿真编程,前面所介绍的仿真程序是针对(3.1-1)所表示的通用的非线性系统，针对如下所表示的线性状态空间模型,(3.4-1),我们虽然也可以用前面介绍的程序仿真，但比较麻烦，实际上系统由矩阵A、B、C、D就可以指定微分方程和输出方程，而不必单独用函数文件来表示。
19、为此，我们单独编写针对系统(3.4-1)的线性系统数值积分法仿真程序。程序仍然由仿真框架和单步积分构成。,21,(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是线性状态空间模型的系数矩阵 % In
20、teMethod是数值积分方法，可以是EUL ,RK2,RK4 % ControlFile是控制作用函数，没有控制算法时，可以省略 % 输出参数：t是仿真结果的时间序列 % y是仿真结果系统的输出序列,22,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;
21、% eval(OutputFile,(tstart,x0,u0);%计算初始输出 y(:,1)=y0 ;%将y0作为输出的第1列 curx=x0; %当前一步的x curu=u0; %当前一步的u cury=y0; %当前一步的y for i=1:1:cntt-1curu=eval(ControlFile,(t(i),h,curx,curu,cury,A,B,C,D);%计算控制时传递的参数：当前时间，步长，当前状态和输出,模型系数矩阵curx=w_LinearStepIntegral(h,curx,curu,InteMethod,A,B);%单步积分运算,针对线性模型cury=C*curx+
22、D*curu;%计算输出y(:,i+1)=cury;%将输出加入到输出序列里 end,23,(2)线性系统单步积分函数,function NewX=w_LinearStepIntegral(h,curx,curu,InteMethod,A,B); % 函数功能：单步积分运算,与模型方程有关 % 输入参数：h是数值积分步长，curx,curu分别是当前的状态和控制向量 % InteMethod是积分算法，字符串类型， % 可以取EUL,RK2,RK4,由于要用于判断，字符串必须等长 % A,B是线性状态模型的微分方程的系数矩阵 % 输出参数：NewX是这一步计算的新的状态向量 Bu=B*curu
23、; 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 = RK2k1=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
24、 %欧拉法EULNewX=curx+h*(A*curx+Bu); end,24,3.5二阶系统传递函数仿真,1. 模型转换,写为状态模型就是,基本参数取值为,针对本例，线性系统仿真程序进行仿真。,25,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,
25、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,y1
26、0,*y,t,y15,.b);%绘制叠加曲线 xlabel(x); ylabel(y); title(二阶系统阶跃输入响应);,26,二阶系统在阶跃输入作用下，不同阻尼比时的输出曲线。,27,3.6面向结构图仿真,对下图的系统进行面向结构图的变换，并进行仿真。,其中,28,(1) 模型处理,每个环节用3个参数来描述，得到的矩阵,29,(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,-
27、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(面向结构图的仿真);,30,G=2,1,1;0,1,3;5,1,6;0,1,1;,只绘制y1和y3，可以看出两个输出的区别，它们都趋向于稳定。,31,同时绘制4个输出时，由于y2,y4变化很大，发散，所以y1和y3看起来好像重叠了。,32,附录：plot函数参数取值意义 Plot函数的一般调用是：plot(x1,y1,option1,x2,y2,option2,.),
