数学实验之7计算机模拟

1、计算机模拟计算机模拟概述 连续系统MATLAB的实现范例一范例二海港系统卸载货物的模拟 可靠性问题 离散系统下一页上一页主 页实验目的1、掌握连续系统模拟的基本方法；2、掌握离散系统模拟的基本方法；(包括 分布的产生、性能指标的设置和计算、 模拟的设计和运行及输出结果的分析等)3、掌握一般的编程要领；返 回下一页上一页主 页概述 计算机模拟是利用计算机对系统系统的结构和行为进行动态演示，以评价或预测系统的行为效果，为决策者提供信息的一种方法。它是解决较复杂的实际问题的一条有效途径。 计算机模拟也可以说是用计算机程序直接建立真实系统的模型，并通过计算了解系统随时间变化的行为或特征。 应用领域：航

2、空、机电、冶金、社会经济、交通运输、生态系统等。 计算机模拟分为连续系统连续系统模拟和离散系统离散系统模拟。返 回下一页上一页主 页 状态随时间连续变化的系统称为连续系统。通常该系统的模型一般可以用微分方程的形式表达，通过一些物理机理推导出来。模拟结果往往是近似的。例如，1)飞机运动轨迹（x(t),y(t),z(t)， 2) 浓度变化问题。连续系统(x1(t),xn(t)模拟方法：将时间变量t离散化，通过某些物理机理或特征，找出相邻时刻下状态之间的关系。即 xi(t+) = fi(x1(t), , xn(t)返 回下一页上一页主 页离散系统 离散系统是指系统状态只在有限的时间点或可数的时间点上

3、发生变化的系统。并假设离散系统状态的变化是在一个时间点上瞬间完成。模型一般用流程图或网络来表示。可能涉及到随机事件等。 关键：模拟步骤、数据收集、模拟时钟与连续系统的区别是时刻不等间距。返 回下一页上一页主 页排队问题 机械故障等候修理 飞机跑道分配日常生活中经常遇到的排队问题： 自选商场收款台收费 医院里病人等候就诊离散模拟的典型问题之一下一页上一页主 页输入情况： 顾客到达时间和服务时间。系统状态： 排队等候的顾客数目（队长）L(t) 服务员是否在工作或服务效率等；简图：,.2 , 1,isDxyiiii关系：系统在什么条件下处于忙碌(空闲)状态？当yixi+1时，置S=0，系统处于空闲状

4、态；第二顾客接受服务时间s2x50 x1x2x3x4 y1 y2 y3 y4 y5D2T?1nDDnii单服务台离散模拟的典型问题之一下一页上一页主 页 排队系统中，顾客到达时刻数据如何收集？对每个顾客的服务时间如何？X： x1，x2，xn第一个顾客到达的时刻第二个顾客到达的时刻计算机遵循某种规则进行随机抽样。S: s1,s2,sn第一个顾客服务时刻返 回下一页上一页主 页随机数的产生1、均匀随机数(均匀分布U0,1) rand()2、产生其他分布随机数的方法 逆变换法、舍选法、近似抽样法等。xttfxFd)()(其中是X的密度函数0F(x)1, 令r = F(x) xexFrFx1)(, )

5、(1例如，指数分布Matalb软件的实现下一页上一页主 页随机数的产生例如指数分布：解出逆变换： ( )1xrF xe Matalb软件的实现1ln(1)xr 在MATLAB下输入： r = rand (1,5)； x = (1/3)*log (1- r)输出结果为：x = 0.0171 0.4882 0.1665 0.2405 0.0384下一页上一页主 页3、(非)常见分布随机数如何产生？（离散）经验分布函数法: 设X服从如下分布律Xf 2 4 6 8 0.1 0.4 0.3 0.2累积分布函数为Matalb软件的实现0,20.1,24( )0.5,460.8,681,8xxF xxxx下

6、一页上一页主 页由反函数的定义知，设随机数r0，1，其转换关系如下 Matalb软件的实现0 2 4 6 810.80.60.40.2rx 0,0.1) 0.1,0.5) 0.5,0.8) 0.8, 1 2 4 6 8rx下一页上一页主 页r=rand(10); y=;for i=1:10 if 0=r(i)&r(i)0.1 y(i)=2; elseif 0.1=r(i)&r(i)0.5 y(i)=4; elseif 0.5=r(i)&r(i)=0 idle(i)=timediff;wait(i)=0; else idle(i)=0;wait(i)=-timediff;

7、 endMatlab程序 start(i)=arrive(i)+wait(i); finish(i)=start(i)+unload(i); harbor(i)=wait(i)+unload(i); if harbor(i)MaxHar MaxHar=harbor(i); else end if wait(i)MaxWait MaxWait=wait(i); else endendF=finish(n);Hartime=sum(harbor)/nWaittime=sum(wait)/nIdletime=sum(idle)/FMaxHar,MaxWait,下一页上一页主 页海港系统卸载货物的模拟

8、海港系统卸载货物的模拟 结果实验次数HartimeMaxHarWaittimeMaxWaitIdletime1111254431680.12297238301710.16393194291420.17488202221470.22589202221190.22696295292320.197110284442040.16897217291430.16当输入船只数N=100，根据上述程序模拟8次，得到如下结果： 下一页上一页主 页海港系统卸载货物的模拟海港系统卸载货物的模拟 如果相邻两只船到达港口的间隔时间以及在港设备卸载货物的服务时间不服从15145和4590分钟的均匀分布。为了更好地模拟实际

9、海港系统，需要收集该系统的历史数据，比如人们对到达港口的1200艘船只进行了调查，得到如下数据：进一步说明 下一页上一页主 页海港系统卸载货物的模拟海港系统卸载货物的模拟 间隔时间 频率卸载时间 频率15-24 0.00925-34 0.02935-44 0.03545-54 0.05155-64 0.09065-74 0.16175-84 0.20085-94 0.17295-104 0.125105-114 0.071115-124 0.037125-134 0.017135-145 0.00345-49 0.01750-54 0.04555-59 0.09560-64 0.08665-6

10、9 0.13070-74 0.18575-79 0.20880-84 0.14385-90 0.091 下一页上一页主 页海港系统卸载货物的模拟海港系统卸载货物的模拟 .009.038.073.124.214.375.575.747.872.943.980.997.1.000.20.40.60.81.0时间频率152535 45 55 65 75 85 95 105 115125 135 145图3 相邻两艘船只到达港口的间隔时间xy下一页上一页主 页海港系统卸载货物的模拟海港系统卸载货物的模拟 利用累积直方图，作分段线性插值折线。例如在图8.3中，区间15，25和25，35的中点分别是20和

11、30，在直角坐标(时间，频率)下，用直线连接平面点(20,0.009), (30,0.038)，可以得到直线方程： 344.816.8966, 0.0090.038xyy200.00930200.0380.009xy下一页上一页主 页海港系统卸载货物的模拟海港系统卸载货物的模拟 产生均匀随机数的区间 对应时间范围 线性插值函数 00.009y1520 x555.615.0000 xy0.0090.038y2030 x344.816.8966xy下一页上一页主 页海港系统卸载货物的模拟海港系统卸载货物的模拟 MATLAB编程如下：(jge.m)a=rand;if 0=a&a0.009 b

12、=555.6*a+15;elseif 0.009=a&a0.038 b=334.8*a+16.8966;elseif 0.038=a&a0.073 b=285.7*a+19.1429; end主程序：n=N;wait=;idle=;harbor=;jge; %调用子程序fwu;between(1)=b;unload(1)=c; 下一页上一页主 页海港系统卸载货物的模拟海港系统卸载货物的模拟 1. 如果港口可以同时卸载多艘船只的货物，将是一个多窗口的服务问题，如何进行模拟？2. 如果海港系统不仅仅考虑卸载货物问题，如还考虑装载货物问题，或卸载货物以后再装载等问题，又应该如何分析？

13、思考：下一页上一页主 页可靠性问题可靠性问题 一设备上有三个相同的轴承，每个轴承正常工作寿命为随机变量，概率分布如下：寿命/h1000110012001300140015001600170018001900概率0.100.130.250.130.090.120.020.060.050.05有轴承损坏设备停止工作检修工准备开始更换部件，称为一个延迟时间，它也是随机变量，分布如下：延迟时间 /min 5 1015概 率0.60.30.1范例二下一页上一页主 页主要费用：1、设备停工损失费： 5 元/分钟；2、检修工人的工时费：12 元/小时；3、轴承的成本费： 16元/个更换轴承所需要的时间： 一

14、个 两个 三个 20 30 40 （min）范例二下一页上一页主 页问题：现在有两种方案：方案一：损坏一个轴承只更换一个轴承；方案二：一旦有轴承损坏就全部更换；试通过计算机模拟对以上两种方案做出评价。 随机数怎样产生？ 模拟时选用时间步长法还是 事件步长法？范例二下一页上一页主 页a=rand; (Lrnd.m)if a=0.1 l=1000; elseif 0.1a&a=0.23 l=1100; elseif 0.23a&a=0.48 l=1200; elseif 0.48a&a=0.61 l=1300; elseif 0.61a&a=0.7 l=1400;

15、elseif 0.7a&a=0.82 l=1500; elseif 0.82a&a=0.84 l=1600; elseif 0.84a&a=0.90 l=1700; elseif 0.90a&a=0.95 l=1800; else l=1900;end下一页上一页主 页关于随机数的产生见 Lrnd.m（零件寿命） Yrnd.m（延迟时间） 方案一的数学模型：kekao1.m目标函数 min c =Ui / T（三个轴承同时损坏）（两个轴承同时损坏）只有一个轴承损坏）3166040125)40(2166030125)30(166020125)20(iiiitttU

16、其中 损失费工时费成本费ti 表示延迟时间 下一页上一页主 页方案二的数学模型：目标函数 min c =Ui / T其中31612325)40(iitU10000)(1iniiLt方案二的情况比较单一。Kekao3.m模拟框图如下：范例二下一页上一页主 页数据初始化T100000产生下一个事件发生时刻根据不同方案确定更换策略产生延迟时间计算当前更换费用累积费用产生新轴承的寿命时间确定下一个事件发生的时刻系统时刻跳转到下一个事件发生时刻输出总费用 停止否是范例二下一页上一页主 页c1=5;c2=12;c3=16; (kekao1.m)g1=20/60;g2=30/60;g3=40/60;Lrnd

17、; % Random function Lrnd gets a random life(three)yrnd; % Random function yrnd gets a random delate(one)lm=sort(l);U=;t=;for j=1:50if lm(1)lm(2) U(j)=c1*(y+g1*60)+c2*g1+c3;elseif lm(1)=lm(2)lm(3) U(j)=c1*(y+g2*60)+c2*g2+2*c3;else U(j)=c1*(y+g3*60)+c2*g3+3*c3;endLrnd1;%(one)t(j)=lm(1)+y/60+l1;L=lm(2)

18、,lm(3),t(j);lm=sort(L);yrnd;endU;t; zU=sum(U) C=zU/t(50)计算结果：ZU=7850（元），C=0.3354（元/小时）。下一页上一页主 页c1=5;c2=12;c3=16; (kekao3.m)g1=20/60;g2=30/60;g3=40/60;Lrnd;% (three)yrnd; % yLrnd1;lm=sort(l);U=;t=;for j=1:50 U(j)=c1*(y+g3*60)+c2*g3+3*c3; t(j)=lm(1)+y/60+l1;L=lm(2),lm(3),t(j);lm=sort(L);yrnd; endU;t;zU=sum(U)C=zU/t(50)计算结果：ZU=14800（元），C=0.6271。（元/小时）范例二下一页上一页主 页计算结果：方案一：c = 0.3307 元/小时；方案二：c = 0.5486（0.7234）元/小时；显然方案一优于方案二。范例二返 回下一页上一页主 页练习