§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
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;
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)，则模型为
变形为：
以小步长将时间离散化，时间每前进一步，就按照“前进尤拉方法”计算，当时就认为击中。得计算公式：
程序（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,0.2,0.25(小时)时，导弹及敌舰的位置：
t=0.075(小时)时
t=0.15(小时)时
t=0.2(小时)时
t=0.25(小时)时
（2）任意t时刻，导弹位置(x,y)，敌舰位置(u,v)，则模型是什么？
t=0时：x,y,u,v,=?
t时刻：
模型为
变形为：
以小步长将时间离散化，时间每前进一步，就按照“前进尤拉方法”计算，当都很小时就认为击中。得计算公式：
只需对前面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，输出
函数文件一：
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+ycsj时刻 开始更换三个轴承
t=t+40时刻 修好了！ 产生本次停工损失
再循环（有限次循环）到“sm=函数文件一”
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
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;
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)，则模型为
变形为：
以小步长将时间离散化，时间每前进一步，就按照“前进尤拉方法”计算，当时就认为击中。得计算公式：
程序（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,0.2,0.25(小时)时，导弹及敌舰的位置：
t=0.075(小时)时
t=0.15(小时)时
t=0.2(小时)时
t=0.25(小时)时
（2）任意t时刻，导弹位置(x,y)，敌舰位置(u,v)，则模型是什么？
t=0时：x,y,u,v,=?
t时刻：
模型为
变形为：
以小步长将时间离散化，时间每前进一步，就按照“前进尤拉方法”计算，当都很小时就认为击中。得计算公式：
只需对前面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，输出
函数文件一：
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+ycsj时刻 开始更换三个轴承
t=t+40时刻 修好了！ 产生本次停工损失
再循环（有限次循环）到“sm=函数文件一”
