數(shù)學(xué)實(shí)驗(yàn)之7計(jì)算機(jī)模擬

1、計(jì)算機(jī)模擬計(jì)算機(jī)模擬概述 連續(xù)系統(tǒng)MATLAB的實(shí)現(xiàn)范例一范例二海港系統(tǒng)卸載貨物的模擬 可靠性問(wèn)題 離散系統(tǒng)下一頁(yè)上一頁(yè)主 頁(yè)實(shí)驗(yàn)?zāi)康?、掌握連續(xù)系統(tǒng)模擬的基本方法；2、掌握離散系統(tǒng)模擬的基本方法；(包括 分布的產(chǎn)生、性能指標(biāo)的設(shè)置和計(jì)算、 模擬的設(shè)計(jì)和運(yùn)行及輸出結(jié)果的分析等)3、掌握一般的編程要領(lǐng)；返 回下一頁(yè)上一頁(yè)主 頁(yè)概述 計(jì)算機(jī)模擬是利用計(jì)算機(jī)對(duì)系統(tǒng)系統(tǒng)的結(jié)構(gòu)和行為進(jìn)行動(dòng)態(tài)演示，以評(píng)價(jià)或預(yù)測(cè)系統(tǒng)的行為效果，為決策者提供信息的一種方法。它是解決較復(fù)雜的實(shí)際問(wèn)題的一條有效途徑。 計(jì)算機(jī)模擬也可以說(shuō)是用計(jì)算機(jī)程序直接建立真實(shí)系統(tǒng)的模型，并通過(guò)計(jì)算了解系統(tǒng)隨時(shí)間變化的行為或特征。 應(yīng)用領(lǐng)域：航

2、空、機(jī)電、冶金、社會(huì)經(jīng)濟(jì)、交通運(yùn)輸、生態(tài)系統(tǒng)等。 計(jì)算機(jī)模擬分為連續(xù)系統(tǒng)連續(xù)系統(tǒng)模擬和離散系統(tǒng)離散系統(tǒng)模擬。返 回下一頁(yè)上一頁(yè)主 頁(yè) 狀態(tài)隨時(shí)間連續(xù)變化的系統(tǒng)稱為連續(xù)系統(tǒng)。通常該系統(tǒng)的模型一般可以用微分方程的形式表達(dá)，通過(guò)一些物理機(jī)理推導(dǎo)出來(lái)。模擬結(jié)果往往是近似的。例如，1)飛機(jī)運(yùn)動(dòng)軌跡（x(t),y(t),z(t)， 2) 濃度變化問(wèn)題。連續(xù)系統(tǒng)(x1(t),xn(t)模擬方法：將時(shí)間變量t離散化，通過(guò)某些物理機(jī)理或特征，找出相鄰時(shí)刻下?tīng)顟B(tài)之間的關(guān)系。即 xi(t+) = fi(x1(t), , xn(t)返 回下一頁(yè)上一頁(yè)主 頁(yè)離散系統(tǒng) 離散系統(tǒng)是指系統(tǒng)狀態(tài)只在有限的時(shí)間點(diǎn)或可數(shù)的時(shí)間點(diǎn)上

3、發(fā)生變化的系統(tǒng)。并假設(shè)離散系統(tǒng)狀態(tài)的變化是在一個(gè)時(shí)間點(diǎn)上瞬間完成。模型一般用流程圖或網(wǎng)絡(luò)來(lái)表示。可能涉及到隨機(jī)事件等。 關(guān)鍵：模擬步驟、數(shù)據(jù)收集、模擬時(shí)鐘與連續(xù)系統(tǒng)的區(qū)別是時(shí)刻不等間距。返 回下一頁(yè)上一頁(yè)主 頁(yè)排隊(duì)問(wèn)題 機(jī)械故障等候修理 飛機(jī)跑道分配日常生活中經(jīng)常遇到的排隊(duì)問(wèn)題： 自選商場(chǎng)收款臺(tái)收費(fèi) 醫(yī)院里病人等候就診離散模擬的典型問(wèn)題之一下一頁(yè)上一頁(yè)主 頁(yè)輸入情況： 顧客到達(dá)時(shí)間和服務(wù)時(shí)間。系統(tǒng)狀態(tài)： 排隊(duì)等候的顧客數(shù)目（隊(duì)長(zhǎng)）L(t) 服務(wù)員是否在工作或服務(wù)效率等；簡(jiǎn)圖：,.2 , 1,isDxyiiii關(guān)系：系統(tǒng)在什么條件下處于忙碌(空閑)狀態(tài)？當(dāng)yixi+1時(shí)，置S=0，系統(tǒng)處于空閑狀

4、態(tài)；第二顧客接受服務(wù)時(shí)間s2x50 x1x2x3x4 y1 y2 y3 y4 y5D2T?1nDDnii單服務(wù)臺(tái)離散模擬的典型問(wèn)題之一下一頁(yè)上一頁(yè)主 頁(yè) 排隊(duì)系統(tǒng)中，顧客到達(dá)時(shí)刻數(shù)據(jù)如何收集？對(duì)每個(gè)顧客的服務(wù)時(shí)間如何？X： x1，x2，xn第一個(gè)顧客到達(dá)的時(shí)刻第二個(gè)顧客到達(dá)的時(shí)刻計(jì)算機(jī)遵循某種規(guī)則進(jìn)行隨機(jī)抽樣。S: s1,s2,sn第一個(gè)顧客服務(wù)時(shí)刻返 回下一頁(yè)上一頁(yè)主 頁(yè)隨機(jī)數(shù)的產(chǎn)生1、均勻隨機(jī)數(shù)(均勻分布U0,1) rand()2、產(chǎn)生其他分布隨機(jī)數(shù)的方法 逆變換法、舍選法、近似抽樣法等。xttfxFd)()(其中是X的密度函數(shù)0F(x)1, 令r = F(x) xexFrFx1)(, )

5、(1例如，指數(shù)分布Matalb軟件的實(shí)現(xiàn)下一頁(yè)上一頁(yè)主 頁(yè)隨機(jī)數(shù)的產(chǎn)生例如指數(shù)分布：解出逆變換： ( )1xrF xe Matalb軟件的實(shí)現(xiàn)1ln(1)xr 在MATLAB下輸入： r = rand (1,5)； x = (1/3)*log (1- r)輸出結(jié)果為：x = 0.0171 0.4882 0.1665 0.2405 0.0384下一頁(yè)上一頁(yè)主 頁(yè)3、(非)常見(jiàn)分布隨機(jī)數(shù)如何產(chǎn)生？（離散）經(jīng)驗(yàn)分布函數(shù)法: 設(shè)X服從如下分布律Xf 2 4 6 8 0.1 0.4 0.3 0.2累積分布函數(shù)為Matalb軟件的實(shí)現(xiàn)0,20.1,24( )0.5,460.8,681,8xxF xxxx下

6、一頁(yè)上一頁(yè)主 頁(yè)由反函數(shù)的定義知，設(shè)隨機(jī)數(shù)r0，1，其轉(zhuǎn)換關(guān)系如下 Matalb軟件的實(shí)現(xiàn)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下一頁(yè)上一頁(yè)主 頁(yè)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,下一頁(yè)上一頁(yè)主 頁(yè)海港系統(tǒng)卸載貨物的模擬

8、海港系統(tǒng)卸載貨物的模擬 結(jié)果實(shí)驗(yàn)次數(shù)HartimeMaxHarWaittimeMaxWaitIdletime1111254431680.12297238301710.16393194291420.17488202221470.22589202221190.22696295292320.197110284442040.16897217291430.16當(dāng)輸入船只數(shù)N=100，根據(jù)上述程序模擬8次，得到如下結(jié)果： 下一頁(yè)上一頁(yè)主 頁(yè)海港系統(tǒng)卸載貨物的模擬海港系統(tǒng)卸載貨物的模擬 如果相鄰兩只船到達(dá)港口的間隔時(shí)間以及在港設(shè)備卸載貨物的服務(wù)時(shí)間不服從15145和4590分鐘的均勻分布。為了更好地模擬實(shí)際

9、海港系統(tǒng)，需要收集該系統(tǒng)的歷史數(shù)據(jù)，比如人們對(duì)到達(dá)港口的1200艘船只進(jìn)行了調(diào)查，得到如下數(shù)據(jù)：進(jìn)一步說(shuō)明 下一頁(yè)上一頁(yè)主 頁(yè)海港系統(tǒng)卸載貨物的模擬海港系統(tǒng)卸載貨物的模擬 間隔時(shí)間 頻率卸載時(shí)間 頻率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 下一頁(yè)上一頁(yè)主 頁(yè)海港系統(tǒng)卸載貨物的模擬海港系統(tǒng)卸載貨物的模擬 .009.038.073.124.214.375.575.747.872.943.980.997.1.000.20.40.60.81.0時(shí)間頻率152535 45 55 65 75 85 95 105 115125 135 145圖3 相鄰兩艘船只到達(dá)港口的間隔時(shí)間xy下一頁(yè)上一頁(yè)主 頁(yè)海港系統(tǒng)卸載貨物的模擬海港系統(tǒng)卸載貨物的模擬 利用累積直方圖，作分段線性插值折線。例如在圖8.3中，區(qū)間15，25和25，35的中點(diǎn)分別是20和

11、30，在直角坐標(biāo)(時(shí)間，頻率)下，用直線連接平面點(diǎn)(20,0.009), (30,0.038)，可以得到直線方程： 344.816.8966, 0.0090.038xyy200.00930200.0380.009xy下一頁(yè)上一頁(yè)主 頁(yè)海港系統(tǒng)卸載貨物的模擬海港系統(tǒng)卸載貨物的模擬 產(chǎn)生均勻隨機(jī)數(shù)的區(qū)間 對(duì)應(yīng)時(shí)間范圍 線性插值函數(shù) 00.009y1520 x555.615.0000 xy0.0090.038y2030 x344.816.8966xy下一頁(yè)上一頁(yè)主 頁(yè)海港系統(tǒng)卸載貨物的模擬海港系統(tǒng)卸載貨物的模擬 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; %調(diào)用子程序fwu;between(1)=b;unload(1)=c; 下一頁(yè)上一頁(yè)主 頁(yè)海港系統(tǒng)卸載貨物的模擬海港系統(tǒng)卸載貨物的模擬 1. 如果港口可以同時(shí)卸載多艘船只的貨物，將是一個(gè)多窗口的服務(wù)問(wèn)題，如何進(jìn)行模擬？2. 如果海港系統(tǒng)不僅僅考慮卸載貨物問(wèn)題，如還考慮裝載貨物問(wèn)題，或卸載貨物以后再裝載等問(wèn)題，又應(yīng)該如何分析？

13、思考：下一頁(yè)上一頁(yè)主 頁(yè)可靠性問(wèn)題可靠性問(wèn)題 一設(shè)備上有三個(gè)相同的軸承，每個(gè)軸承正常工作壽命為隨機(jī)變量，概率分布如下：壽命/h1000110012001300140015001600170018001900概率0.100.130.250.130.090.120.020.060.050.05有軸承損壞設(shè)備停止工作檢修工準(zhǔn)備開(kāi)始更換部件，稱為一個(gè)延遲時(shí)間，它也是隨機(jī)變量，分布如下：延遲時(shí)間 /min 5 1015概 率0.60.30.1范例二下一頁(yè)上一頁(yè)主 頁(yè)主要費(fèi)用：1、設(shè)備停工損失費(fèi)： 5 元/分鐘；2、檢修工人的工時(shí)費(fèi)：12 元/小時(shí)；3、軸承的成本費(fèi)： 16元/個(gè)更換軸承所需要的時(shí)間： 一

14、個(gè) 兩個(gè) 三個(gè) 20 30 40 （min）范例二下一頁(yè)上一頁(yè)主 頁(yè)問(wèn)題：現(xiàn)在有兩種方案：方案一：損壞一個(gè)軸承只更換一個(gè)軸承；方案二：一旦有軸承損壞就全部更換；試通過(guò)計(jì)算機(jī)模擬對(duì)以上兩種方案做出評(píng)價(jià)。 隨機(jī)數(shù)怎樣產(chǎn)生？ 模擬時(shí)選用時(shí)間步長(zhǎng)法還是 事件步長(zhǎng)法？范例二下一頁(yè)上一頁(yè)主 頁(yè)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下一頁(yè)上一頁(yè)主 頁(yè)關(guān)于隨機(jī)數(shù)的產(chǎn)生見(jiàn) Lrnd.m（零件壽命） Yrnd.m（延遲時(shí)間） 方案一的數(shù)學(xué)模型：kekao1.m目標(biāo)函數(shù) min c =Ui / T（三個(gè)軸承同時(shí)損壞）（兩個(gè)軸承同時(shí)損壞）只有一個(gè)軸承損壞）3166040125)40(2166030125)30(166020125)20(iiiitttU

16、其中 損失費(fèi)工時(shí)費(fèi)成本費(fèi)ti 表示延遲時(shí)間 下一頁(yè)上一頁(yè)主 頁(yè)方案二的數(shù)學(xué)模型：目標(biāo)函數(shù) min c =Ui / T其中31612325)40(iitU10000)(1iniiLt方案二的情況比較單一。Kekao3.m模擬框圖如下：范例二下一頁(yè)上一頁(yè)主 頁(yè)數(shù)據(jù)初始化T100000產(chǎn)生下一個(gè)事件發(fā)生時(shí)刻根據(jù)不同方案確定更換策略產(chǎn)生延遲時(shí)間計(jì)算當(dāng)前更換費(fèi)用累積費(fèi)用產(chǎn)生新軸承的壽命時(shí)間確定下一個(gè)事件發(fā)生的時(shí)刻系統(tǒng)時(shí)刻跳轉(zhuǎn)到下一個(gè)事件發(fā)生時(shí)刻輸出總費(fèi)用 停止否是范例二下一頁(yè)上一頁(yè)主 頁(yè)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)計(jì)算結(jié)果：ZU=7850（元），C=0.3354（元/小時(shí)）。下一頁(yè)上一頁(yè)主 頁(yè)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)計(jì)算結(jié)果：ZU=14800（元），C=0.6271。（元/小時(shí)）范例二下一頁(yè)上一頁(yè)主 頁(yè)計(jì)算結(jié)果：方案一：c = 0.3307 元/小時(shí)；方案二：c = 0.5486（0.7234）元/小時(shí)；顯然方案一優(yōu)于方案二。范例二返 回下一頁(yè)上一頁(yè)主 頁(yè)練習(xí)