系统建模与仿真-第6次课-第三章(2).ppt

1、1,回顾,2,其中：,四阶龙格库塔法(简称RK4)公式,3,（3.11）,实时二阶龙格库塔法：,4,四阶五级的实时RK计算公式如下：,（3.12）,5,所有龙格库塔公式都有以下特点：,(1)在计算 时只用到 ，而不直接用 ， 等项。换句话说，在后一步的计算中，仅仅利用前一步的计算结果，所以称为单步法。显然它不仅能使存储量减小，而且此法可以自启动，即已知初值后，不必用别的方法来帮助，就能由初值逐步计算得到后续各时间点上的仿真值。,(2)步长 在整个计算中并不要求固定，可以根据精度要求改变，但是在一步中算若干个系数 (俗称龙格库塔系数)，则必须用同一个步长 。,6,(3)龙格库塔法的精度取决于步长

2、 的大小及方法的 阶次。许多计算实例表明：为达到相同的精度，四阶方法的 可以比二阶方法的 大10倍，而四阶方法的每步计算量仅比二阶方法大1倍，所以总的计算量仍比二阶方法小。正是由于上述原因，一般系统进行数字仿真常用四阶龙格库塔公式。值得指出的是：高于四阶的方法由于每步计算量将增加较多，而精度提高不快，因此使用也比较少。,7,3.2.4 面向方程的龙格-库塔法仿真举例,为了加深对数值积分法的理解，以便能使用仿真程序来解决实际问题，本小节将给出面向方程的龙格库塔法仿真实例。“面向方程”是指仿真对象是以一阶微分方程组形式描述的系统，如下所示：,（3.13）,8,其中， 是n维状态向量； 是k维参数向

3、量； 是m维控制向量； 是n维函数向量； 是给定的初始状态。,9,10,令： ， ，而其中的,11,首先，令,计算,12,即,13,计算,14,即,15,计算,16,即,17,计算,18,即,19,即,20,例3.3下面对仿真对象作一说明。设卫星在空中运行的运动方程为：,（3.14）,其中 是重力系数( 401408km3s2)。卫星轨道采用极坐标表示，通过仿真研究卫星轨道。,21,（3.15）,这是一个二阶微分方程组，首先要将其转换成一阶微分方程组。若设 ， ， , , 则有：,22,23,平面极坐标与直角坐标的几何关系,24,建立方程组(3.15)的初始条件：,(卫星到地心的距离，即地球半

4、径),仿真研究的目的是要计算出卫星轨道参量 和 。,25,目标：用MATLAB语言编写4阶龙格-库塔法子程序，解决上述实际仿真问题。,其中,26,27,function KK=ff(t,X) k=401408; f1=X(3); f2=X(4); f3=-k/X(1)2+X(1)*X(4)2; f4=-2*X(3)*X(4)/X(1); KK=f1,f2,f3,f4;,28,29,30,31,h=0.01; x(1,:)=6400,0.7,8,0.32; for k=1:1:101 t(k)=(k-1)*h; xx=x(k,:); K1=ff(t(k),xx); K2=ff(t(k)+h/2,

5、xx+(h/2)*K1); K3=ff(t(k)+h/2,xx+(h/2)*K2); K4=ff(t(k)+h,xx+h*K3); x(k+1,:)=x(k,:)+(h*(K1+2*K2+2*K3+K4)/6); end for k=1:1:101 y(k,1)=x(k,1)*cos(x(k,2); y(k,2)=x(k,1)*sin(x(k,2); end plot(y(:,1),y(:,2),r-);,function KK=ff(t,X) k=401408; f1=X(3); f2=X(4); f3=-k/X(1)2+X(1)*X(4)2; f4=-2*X(3)*X(4)/X(1); K

6、K=f1,f2,f3,f4;,32,在MATLAB的工作空间中键入：y，得到,y = 1.0e+003 * 4.8950 4.1230 4.8819 4.1387 4.8687 4.1544 4.8556 4.1701 4.8425 4.1859 4.8293 4.2016 4.8162 4.2173 4.8031 4.2330 4.7899 4.2487 4.7768 4.2644 4.7637 4.2801 4.7505 4.2959 4.7374 4.3116 4.7243 4.3273 4.7111 4.3430,33,34,另外一种MATLAB程序设计方法,35,h=0.01; x=

7、6400,0.7,8,0.32; fid=fopen(zxc.dat,a+); for k=0:1:100 t=k*h; y1=x(1)*cos(x(2); y2=x(1)*sin(x(2); fprintf(fid,n%10.4f,%10.4f,y1,y2); xx=x; K1=ff(t,xx); K2=ff(t+h/2,xx+(h/2)*K1); K3=ff(t+h/2,xx+(h/2)*K2); K4=ff(t+h,xx+h*K3); x=x+(h*(K1+2*K2+2*K3+K4)/6); end fclose(fid); load zxc.dat; zxc; plot(zxc(:,1

8、),zxc(:,2),b-);,36,例3.4,考虑下述系统的输入输出仿真问题。,单输入单输出3阶非线性系统，构造状态变量,得,37,构造非线性状态反馈控制律,得到闭环系统状态方程,38,令,则,闭环系统的特征多项式为,39,由于,所以，上述特征多项式有一个零根和两个负实根。,闭环系统为临界稳定，现在用数字仿真来验证我们的结论。,首先，给定初始状态,40,将闭环系统的状态方程写成矩阵形式,其中,41,用MATLAB编写右端子函数,42,function K=F(t,X) f1=X(2); f2=X(3); f3=-2*X(3)+(sin(X(3)/3-1)*X(2); K=f1,f2,f3;,

9、设计RK-4的MATLAB程序,43,h=0.01; x=0.6,-0.2,0.26; fid1=fopen(zxc1.dat,a+); for k=0:1:20 t=k*h; u=-x(1)2+x(2); fprintf(fid1,%10.4f,%10.4f,%10.4f,%10.4f,%10.4fn,t,x(1),x(2),x(3),u); xx=x; K1=F(t,xx); K2=F(t+h/2,xx+(h/2)*K1); K3=F(t+h/2,xx+(h/2)*K2); K4=F(t+h,xx+h*K3); x=x+(h*(K1+2*K2+2*K3+K4)/6); end,44,fcl

10、ose(fid1); load zxc1.dat; zxc1; figure; plot(zxc1(:,1),zxc1(:,2),b-); hold on; plot(zxc1(:,1),zxc1(:,3),r-); plot(zxc1(:,1),zxc1(:,4),k-); legend(x_1,x_2,x_3); figure; plot(zxc1(:,1),zxc1(:,5),b-);,45,输出程序运行结果,46,可见数字仿真的结果证实了我们的理论分析结果。,47,例3.5,试用RK-4方法对下述系统的状态响应进行数字仿真。,系统模型,48,将闭环系统的状态方程写成矩阵形式,其中,49

11、,用MATLAB编写右端子函数,50,function K=F_1(t,X) f1=-3*X(1)+0.1*sin(X(2); f2=0.2*X(1)-5*X(2)+X(1)*X(3); f3=X(1)2/3-2*X(3); K=f1,f2,f3;,51,给定初始状态,52,主程序,h=0.01; x=0.3,-0.8,0.45; fid2=fopen(zxc2.dat,a+); for k=0:1:200 t=k*h; fprintf(fid2,%10.4f,%10.4f,%10.4f,%10.4fn,t,x(1),x(2),x(3); xx=x; K1=F_1(t,xx); K2=F_1(

12、t+h/2,xx+(h/2)*K1); K3=F_1(t+h/2,xx+(h/2)*K2); K4=F_1(t+h,xx+h*K3); x=x+(h*(K1+2*K2+2*K3+K4)/6); end,53,fclose(fid2); load zxc2.dat; zxc2; figure; plot(zxc2(:,1),zxc2(:,2),b-); hold on; plot(zxc2(:,1),zxc2(:,3),r-); plot(zxc2(:,1),zxc2(:,4),k-); legend(x_1,x_2,x_3); figure; plot3(zxc2(:,2),zxc2(:,3)

13、,zxc2(:,4),b-); grid;,54,仿真结果,55,56,给定带有时滞的微分方程组,及其初始条件,试分析4阶龙格-库塔法的算法结构。,例3.6,57,这是一个带有状态时滞的微分方程组，通常的4阶龙格-库塔法是无法处理的。,必须对4阶龙格-库塔法进行必要的改进。,具体方法如下：,58,令,59,对本例而言，4阶龙格-库塔法公式具体为,其中,不妨称为改进的龙格-库塔法,60,具体为,61,62,63,64,65,66,67,除了对上述过程进行改进以外，其余的都和常规的4阶龙格-库塔法一样。,主要改进措施是在微分方程右端函数的自变量中增加一项带有时滞的状态向量 ，即：,68,编写右端函

14、数的M文件：,function dx=ff_4(t,x,x_tao) dx(1)=-t*x(1)+3*x(2); dx(2)=x(1)-2*x_tao(2); dx=dx;,69,编写主程序的M文件：,h=0.001; tao=0.3; % 滞后时间 kp=tao/h; % 滞后步数 x(1,:)=0.5,-0.8; for k=1:1:2000 t(k)=(k-1)*h; xx=x(k,:); if k=kp xx_tao=0,0; else xx_tao=x(k-kp,:); end,70,K1=ff_4(t(k), xx, xx_tao); K2=ff_4(t(k)+h/2, xx+(h

15、/2)*K1, xx_tao+(h/2)*K1); K3=ff_4(t(k)+h/2, xx+(h/2)*K2, xx_tao+(h/2)*K2); K4=ff_4(t(k)+h, xx+h*K3, xx_tao+h*K3); x(k+1,:)=x(k,:)+(h*(K1+2*K2+2*K3+K4)/6); End t(2001)=2; plot(t,x(:,1),r-); hold on; plot(t,x(:,2),b-); legend(x_1,x_2);,71,运行结果如下：,72,如果计算到t=6s，则结果如下：,73,可见该时滞系统是渐近稳定的。,如果把滞后时间改为 ，,仿真结果如

16、下：,74,如果把滞后时间改为 ，,仿真结果如下：,75,可见当时滞达到一定程度时，该系统变为不稳定。,76,例3.8,给定一个带有时变参数的微分方程组,初始状态为 ，,时变参数 满足变化规律,77,试用4阶龙格-库塔法对该微分方程组进行数字仿真，时间范围为 。,编写右端函数的M文件,function dx=ff_5(t,x) t0= 0, 1, 2, 3, 4, 5, 6; p0=1, 0.8, 0.2, 0.1, 0.7, 0.4, 0.5; p=interp1(t0,p0,t,linear); A=-2, 0.3; 1, -13*p; dx=A*x;,78,编写龙格-库塔法主程序的M文件,h=0.001; x(1,:)=3,-4; for k=1:1:6000 t(k)=(k-1)*h; xx=x(k,:); K1=ff_5(t(k),xx); K2=ff_5(t(k)+h/2, xx+(h/2)*K1); K3=ff_5(t(k)+h/2, xx+(h/2)*K2); K4=ff_5(t(k)+h, xx+h*K3); x(k+1,:)=x(k,:)+(h*(K1+2*K2+2*K3+K4)/6); end,79