§8-4海港系统卸载货物的模拟 84.1关于原问题的模拟 读书P19-P122之“流程框图” 关于下面主程序(syp123)的说明: g是100个数,服从“15到145的均匀分布”,代表100个 船到来的间隔时间; G是100个数,代表100个船的到来时刻 f是100个数,服从“45到90的均匀分布”,代表100个 船的卸载时间; 是100个数,代表100个船的离开时刻 g=15+130*rand(1,100) G(1)=g(1) fori=2:100 G(=G(-l)+g(, f=45+45*rand(1,100); F(1)=g(1)+f(1) fori=2:100 F(=max(F(i-1)G()+f(); d ZXZSJ=sum(f); ZKXSJ=F(100)-ZXZSJ KXL=ZKXSJ/F(100) dd(1)0tl()=d(1)+f(1) dd()=max(F(-)G(0), tl(Fdd(o+f(o n(dd),pjtl-mean(tl) zedd=max(dd), zctl-=max(tl) 由于随机因素,每次运行,机器会产生不同的随 机数,故计算结果也会不同。下表给出了十次计算得 到的结果 次平均停留最长停留「平均等待|最长等待空闲率 112 275 189 0.16 4134 0.1l 123 218 0.13 202 21 146 0.2 93 127 0.19 0.13 书P124之“想(1)”时间间隔改为“10到90的均 匀分布”,则结果如下: 次平均停留最长停留「平均等待|最长等待空闲 781 1389 713 1323 0.0180 752 1358 1307 0.0053
§8—4 海港系统卸载货物的模拟 8.4.1 关于原问题的模拟 读书 P119-----P122 之“流程框图” 关于下面主程序（syp123）的说明： g 是 100 个数，服从“15 到 145 的均匀分布”，代表 100 个 船到来的间隔时间； G 是 100 个数，代表 100 个船的到来时刻； f 是 100 个数，服从“45 到 90 的均匀分布”，代表 100 个 船的卸载时间； F 是 100 个数，代表 100 个船的离开时刻 clear g=15+130*rand(1,100); G(1)=g(1); for i=2:100 G(i)=G(i-1)+g(i); end f=45+45*rand(1,100); F(1)=g(1)+f(1); for i=2:100 F(i)=max(F(i-1),G(i))+f(i); end ZXZSJ=sum(f);ZKXSJ=F(100)-ZXZSJ; KXL=ZKXSJ/F(100) dd(1)=0;tl(1)=dd(1)+f(1); for i=2:100 dd(i)=max(F(i-1)-G(i),0); tl(i)=dd(i)+f(i); end pjdd=mean(dd),pjtl=mean(tl) zcdd=max(dd),zctl=max(tl) 由于随机因素，每次运行，机器会产生不同的随 机数，故计算结果也会不同。下表给出了十次计算得 到的结果： 次 平均停留 最长停留 平均等待 最长等待 空闲率 1 101 322 36 267 0.22 2 88 200 20 146 0.20 3 112 275 44 189 0.16 4 134 296 65 246 0.11 5 123 269 56 181 0.16 6 117 271 50 218 0.13 7 88 202 21 146 0.22 8 92 193 25 127 0.19 9 103 256 37 189 0.15 10 106 291 37 223 0.13 书 P124 之“想(1)”：时间间隔改为“10 到 90 的均 匀分布”，则结果如下： 次 平均停留 最长停留 平均等待 最长等待 空闲率 1 781 1389 713 1323 0.0180 2 752 1358 686 1307 0.0053
822 1524 753 1448|0.0079 8.4.2根据统计数据确定间隔时间、卸载时间 看书P124之“表89”,如何根据统计规律来产生 间隔时间、卸栽时间随机数呢? 法一:书P125-P127之方法。 法二:命令g=mand(1,100)产生0到1的100个均 匀数,逐个检验每个g(i),若它介于0--0.009则令 g(i)=20就是第i个间隔时间,若它介于0009-0038 则令g(i)=30就是第ⅰ个间隔时间,若它介于 0.038-0073则令g(}=40就是第ⅰ个间隔时间,等等。 以同样方法产生卸载时间。 下面两个函数文件就是用“法二”的方法,根据 表89”的数据,产生100个船的到来间隔g、卸载 时间f function g=syp125hswjl g=rand(1,100) fori=1:100 ifg(i)<0.0009 g()=20 elseif g(i<0.038 g()=30 elseif g(i<0.073 g()=40 elseif g(i<0. 124 g()=50 elseif g(<0. 214 g() elseif g(i)<0.375 g(1)=70 elseif g(<0.575 g()=80 elseif g(i<0.747 g(1)= elseif g(i)<0.872 g()=100 elseif g(<0.943 g()=110 elseif g(<0.980 g()=120 elseif g(i)×0.997 g(1)=130; else g()= end function f=syp125hswj2 f=rand(1,100) fori=1:100 if()<0.017 f(i)=47;
3 822 1524 753 1448 0.0079 8.4.2 根据统计数据确定间隔时间、卸载时间 看书 P124 之“表 8.9”，如何根据统计规律来产生 间隔时间、卸栽时间随机数呢？ 法一：书 P125----P127 之方法。 法二：命令 g=rand(1,100)产生 0 到 1 的 100 个均 匀数，逐个检验每个 g(i)，若它介于 0---0.009 则令 g(i)=20 就是第 i 个间隔时间，若它介于 0.009---0.038 则令 g(i)=30 就是第 i 个间隔时间，若它介于 0.038---0.073 则令 g(i)=40 就是第 i 个间隔时间，等等。 以同样方法产生卸载时间。 下面两个函数文件就是用“法二”的方法，根据 “表 8.9”的数据，产生 100 个船的到来间隔 g、卸载 时间 f ： function g=syp125hswj1 g=rand(1,100); for i=1:100 if g(i)<0.0009 g(i)=20; elseif g(i)<0.038 g(i)=30; elseif g(i)<0.073 g(i)=40; elseif g(i)<0.124 g(i)=50; elseif g(i)<0.214 g(i)=60; elseif g(i)<0.375 g(i)=70; elseif g(i)<0.575 g(i)=80; elseif g(i)<0.747 g(i)=90; elseif g(i)<0.872 g(i)=100; elseif g(i)<0.943 g(i)=110; elseif g(i)<0.980 g(i)=120; elseif g(i)<0.997 g(i)=130; else g(i)=140; end end 第二个： function f=syp125hswj2 f=rand(1,100); for i=1:100 if f(i)<0.017 f(i)=47;
seif f(i<0. 062 elseif f(i)×0.157 elseif f(i)×0.243 f(i)=62; elseif f(i)×0.373 f()=67 elseif f(<0.558 f(i)= iff(i)<0.7 f(i)=77; f(i)<0.909 else g()=87 d d 再将主程序(sypl23)中的两句修改 g=15+130*rand(1,100),改为g=syp125 howl f=45+45*rand(1,100)改为f=syp12hsw2; 执行结果 次平均停留|最长停留平均等待|最长等待空闲率 147 70 0.26 81 18 0.23 0.21 119 0.23 12 0.18 §8-5连续系统的计算机模拟 当系统的状态是随时间连续变化时,就称为连 续系统。如处在飞行中的导弹,其位置坐标就是连续 变化的。对连续系统的计算机模拟,就是近似地获取 在一些离散的时刻点上的数值,利用这些数值模拟系 统的运行过程 例:t=0时刻,在导弹基地正北120km处,敌舰向东 行驶速度90km/h,同时导弹发出,速度450kmh。飞 行中,任意时刻导弹的方向总是指向目标 (1)若敌舰不变速、不变向,何时何地击中? (2)从0时刻起,敌舰发现导弹,任意时刻速度大 小为135,方向总与导弹垂直,何时何地击中? 解:(1)任意t时刻导弹位置(xy),则模型为 dx 90t-x 450 dt 0 变形为
elseif f(i)<0.062 f(i)=52; elseif f(i)<0.157 f(i)=57; elseif f(i)<0.243 f(i)=62; elseif f(i)<0.373 f(i)=67; elseif f(i)<0.558 f(i)=72; elseif f(i)<0.766 f(i)=77; elseif f(i)<0.909 f(i)=82; else g(i)=87; end end 再将主程序（syp123）中的两句修改： g=15+130*rand(1,100); 改为 g=syp125hswj1; f=45+45*rand(1,100); 改为 f=syp125hswj2; 执行结果： 次 平均停留 最长停留 平均等待 最长等待 空闲率 1 68 147 7 70 0.26 2 81 194 18 120 0.23 3 81 210 17 157 0.22 4 72 133 6 61 0.21 5 71 119 7 49 0.23 6 83 182 17 112 0.18 §8—5 连续系统的计算机模拟 当 系统的状态 是随时间连续变化时，就称为连 续系统。如处在飞行中的导弹，其位置坐标就是连续 变化的。对连续系统的计算机模拟，就是近似地获取 在一些离散的时刻点上的数值，利用这些数值模拟系 统的运行过程。 例：t=0 时刻，在导弹基地正北 120km 处，敌舰向东 行驶速度 90km/h，同时导弹发出，速度 450 km/h。飞 行中，任意时刻导弹的方向总是指向目标。 （1）若敌舰不变速、不变向，何时何地击中？ （2）从 0 时刻起，敌舰发现导弹，任意时刻速度大 小为 135，方向总与导弹垂直，何时何地击中？ 解：（1）任意 t 时刻导弹位置(x,y)，则模型为 = = = + − − = = = | 0, | 0. 450, ( ) ( ) , 90 120 0 0 2 2 t t x y dt dx dy t x y dx dy 变形为：
dx 450.90t 450·120 dt√90-x)2+(120-y)2dt√(90n-x)2+(120-y)2 以小步长τ将时间离散化,时间每前进一步,就 按照“前进尤拉方法”计算xy,当20-y≈0时就 认为击中。得计算公式: x1=0,y1=0,k=1,2,3 000025, +T 90kr-x4)2+( 450·120 Vk=V+T (90kr-x)2+(120-y) 程序(sp13la) x=Oy=0t=0.00025; fork=1:10000 120- X=x+450*t*abs(p)w,y=y+450*t*abs(q)/ if abs((120-y)=0.01 k, jzsjk*t,jZwz=X 结果:k=l11js=0.2778jzz=249963 图形演示: x=Oy=0t=0.00025; hold on plot(0,OJ.150J)plot(40.40]0.150]) plot([O,40][0,0)plot(0,40]150,150grid fork=1:10000 p=90*k*t-x q=120-y; w=sqrt(p 2+q2); X-X+450 t*abs(p)w, y=y+450*t*abs(q plot(x, y, . ), plot(90*kt, 120, .) if abs((120-y)=0.01 break end k, jzsk*t, jZwzX hold off
2 2 2 2 (90 ) (120 ) 450 120 , (90 ) (120 ) 450 90 t x y y dt dy t x y t x dt dx − + − − = − + − − = 以小步长 将时间离散化，时间每前进一步，就 按照“前进尤拉方法”计算 x, y ，当 120 − y 0 时就 认为击中。得计算公式： . (90 ) (120 ) 450 120 , (90 ) (120 ) 450 90 0, 0, 1,2,3,......, 0.00025, 2 2 1 2 2 1 1 1 k k k k k k k k k k k x y y y y k x y k x x x x y k − + − − = + − + − − = + = = = = + + 程序（syp131a） clear x=0;y=0;t=0.00025; for k=1:10000 p=90*k*t-x;q=120-y;w=sqrt(p^2+q^2); x=x+450*t*abs(p)/w;y=y+450*t*abs(q)/w; if abs(120-y)<=0.01 break end end k,jzsj=k*t,jzwz=x 结果： k =1111 jzsj =0.2778 jzwz =24.9963 图形演示： clear x=0;y=0;t=0.00025; hold on plot([0,0],[0.150]),plot([40,40],[0.150]) plot([0,40],[0,0]),plot([0,40],[150,150]),grid for k=1:10000 p=90*k*t-x;q=120-y;w=sqrt(p^2+q^2); x=x+450*t*abs(p)/w;y=y+450*t*abs(q)/w; plot(x,y,'.'),plot(90*k*t,120,'.') if abs(120-y)<=0.01 break end end k,jzsj=k*t,jzwz=x hold off
下面几个图,分别显示时间为0.075,0.15,02 025(小时)时,导弹及敌舰的位置 t=0075(小时) 152025 t=0.15(小时)时 t=0.2(小时)时
下面几个图，分别显示时间为 0.075, 0.15, 0.2, 0.25(小时)时，导弹及敌舰的位置： t=0.075(小时)时 t=0.15(小时)时 t=0.2(小时)时
t=0.25(小时)时 (2)任意t时刻,导弹位置(xy),敌舰位置(u,v), 则模型是什么? t=0时:xy,uv,=? 时刻:=?中 dy dv dx dt 模型为(ay+时=40.c)+(d)=135 xl=0=0,yl=0=0.,ul=0=0,v1-=120 变形为 450.l-x 450· dt√(-x2+(v-y)2d√-x)2+(v-y)2 135-ydh135y-x
t=0.25(小时)时 （2）任意 t 时刻，导弹位置(x,y)，敌舰位置(u,v)， 则模型是什么？ t=0 时：x,y,u,v,=? t 时刻： = ? = ? du dv dx dy 模型为 = = = = = + = + − − = − − − = = = = = | 0, | 0, | 0, | 120 . 135, ( ) ( ) 450, ( ) ( ) , , 0 0 0 0 2 2 2 2 t t t t x y u v dt du dv dt dx dy v y u x du dv u x v y dx dy 变形为： . ( ) ( ) 135 , ( ) ( ) 135 , ( ) ( ) 450 , ( ) ( ) 450 2 2 2 2 2 2 2 2 u x v y u x dt dv u x v y v y dt du u x v y v y dt dy u x v y u x dt dx − + − − = − − + − − = − + − − = − + − − =
以小步长τ将时间离散化,时间每前进一步,就按照 “前进尤拉方法”计算xyu,y,当-y-x都很小 时就认为击中。得计算公式 x1=0,y1=0,1=0,v1=120,k=1,2,3…,r=000025 450 yu )2+(v-y)2 k+1 √a4-x)+(v4-y)2 只需对前面sypl3la略做修改补充,得程序(sypl31b) x=Oy=0;u=0V=120;t=0.00025; hold on fork=1:10000 p=u-X, q=v-y w=sqrt(p 2+q2); X=X+450*t*abs(p)w;y=y+450*t*abs(q)w u=u+135*t*abs(qw, v=v-135*t*abs(p)w; if abs(v-y )=0.3&abs(u-x =0.3 end k,jzs=*t,jzwz=x,y] 结果 k=1065jzsj=0.266 zWz=(328595,10.3626 图形演示: =0y=0u=0V=120;t=0.00025 hold plo(.oJ[0.150plot(40.40][0.150 plo([0.400,0 D).plot(040]150,150]grid fork=1:10000 p=u-x q=V-y w=sqrt(p 2+q2); x=x+450*t*abs(p)/w,y=y+450*t*abs(q)w u=u+135*t*abs(q w, v=v-135* t*abs(p)/w; plot(x, y, . ),plot(u, v,!) if abs(v-y)=0.3&abs(u-xk=0.3 break d d k,jzsjk*,jZwz=x, yI
以小步长 将时间离散化，时间每前进一步，就按照 “前进尤拉方法”计算 x, y,u,v ，当 v − y , u − x 都很小 时就认为击中。得计算公式： . ( ) ( ) 135 ...... , , ...... , ( ) ( ) 450 0, 0, 0, 120, 1,2,3,......, 0.00025, 2 2 1 1 1 2 2 1 1 1 1 1 k k k k k k k k k k k k k k k k k k u x v y u x u v v y u x v y u x x x x y u v k − + − − = = − = − + − − = + = = = = = = + + + + 只需对前面 syp131a 略做修改补充，得程序（syp131b） clear x=0;y=0;u=0;v=120;t=0.00025; hold on for k=1:10000 p=u-x;q=v-y;w=sqrt(p^2+q^2); x=x+450*t*abs(p)/w;y=y+450*t*abs(q)/w; u=u+135*t*abs(q)/w;v=v-135*t*abs(p)/w; if abs(v-y)<=0.3&abs(u-x)<=0.3 break end end k,jzsj=k*t,jzwz=[x,y] 结果： k=1065 jzsj=0.2662 jzwz =(32.8595,110.3626) 图形演示： clear x=0;y=0;u=0;v=120;t=0.00025; hold on plot([0,0],[0.150]),plot([40,40],[0.150]) plot([0,40],[0,0]),plot([0,40],[150,150]),grid for k=1:10000 p=u-x;q=v-y;w=sqrt(p^2+q^2); x=x+450*t*abs(p)/w;y=y+450*t*abs(q)/w; u=u+135*t*abs(q)/w;v=v-135*t*abs(p)/w; plot(x,y,'.'),plot(u,v,'.') if abs(v-y)<=0.3&abs(u-x)<=0.3 break end end k,jzsj=k*t,jzwz=[x,y]
§8-6操练 1.补充讲“函数文件可以被反复调用” 如:对于多个自然数k,输出∑2,R 函数文件 function pfh=hswjl(k) for F =l-k pfh=pfh+r2 end 函数文件 function kjc=hswj2(k) kjc=l for Flk e小c=kj 主程序 for k1: 5 pfh=hswjI(k) kjc=hswj2(k) Ik, pfh, kjc] end 2.关于操练一之第二方案的提示 函数文件一:产生3个寿命数(服从“表8.13”) 函数文件二:产生1个延迟时间(服从“表8.14” 主程序框图: t=0 函数文件 停工时刻t=tmin(sm)ycsj=函数文件二
§8—6 操 练 1．补充讲“函数文件可以被反复调用” 如：对于多个自然数 k，输出 , !. 1 2 i k k i = 函数文件一： function pfh=hswj1(k) pfh=0; for i=1:k pfh=pfh+i^2; end 函数文件二： function kjc=hswj2(k) kjc=1; for i=1:k kjc=kjc*i; end 主程序： clear for k=1:5 pfh=hswj1(k); kjc=hswj2(k); [k,pfh,kjc] end 2．关于操练一之第二方案的提示 函数文件一：产生 3 个寿命数（服从“表 8.13”）. 函数文件二：产生 1 个延迟时间（服从“表 8.14”）. 主程序框图： t=0 sm=函数文件一 停工时刻 t=t+min(sm) ycsj=函数文件二
t= t+ycs]时刻开始更换三个轴承 t=t+40时刻修好了!产生本次停工损失 再循环(有限次循环)到“sm=函数文件一
t=t+ycsj 时刻 开始更换三个轴承 t=t+40 时刻 修好了！ 产生本次停工损失 再循环（有限次循环）到“sm=函数文件一