數(shù)學(xué)實驗之7計算機模擬

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

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

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

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

5、(1例如，指數(shù)分布Matalb軟件的實現(xiàn)下一頁上一頁主 頁隨機數(shù)的產(chǎn)生例如指數(shù)分布：解出逆變換： ( )1xrF xe Matalb軟件的實現(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下一頁上一頁主 頁3、(非)常見分布隨機數(shù)如何產(chǎn)生？（離散）經(jīng)驗分布函數(shù)法: 設(shè)X服從如下分布律Xf 2 4 6 8 0.1 0.4 0.3 0.2累積分布函數(shù)為Matalb軟件的實現(xiàn)0,20.1,24( )0.5,460.8,681,8xxF xxxx下

6、一頁上一頁主 頁由反函數(shù)的定義知，設(shè)隨機數(shù)r0，1，其轉(zhuǎn)換關(guān)系如下 Matalb軟件的實現(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下一頁上一頁主 頁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; endMatlab程序

7、 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,下一頁上一頁主 頁海港系統(tǒng)卸載貨物的模擬海港系統(tǒng)卸載貨物的模擬

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

9、的歷史數(shù)據(jù)，比如人們對到達港口的1200艘船只進行了調(diào)查，得到如下數(shù)據(jù)：進一步說明 下一頁上一頁主 頁海港系統(tǒng)卸載貨物的模擬海港系統(tǒng)卸載貨物的模擬 間隔時間 頻率卸載時間 頻率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-69 0.13070-74

10、 0.18575-79 0.20880-84 0.14385-90 0.091 下一頁上一頁主 頁海港系統(tǒng)卸載貨物的模擬海港系統(tǒng)卸載貨物的模擬 .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下一頁上一頁主 頁海港系統(tǒng)卸載貨物的模擬海港系統(tǒng)卸載貨物的模擬 利用累積直方圖，作分段線性插值折線。例如在圖8.3中，區(qū)間15，25和25，35的中點分別是20和30，在直角坐標(時間，

11、頻率)下，用直線連接平面點(20,0.009), (30,0.038)，可以得到直線方程： 344.816.8966, 0.0090.038xyy200.00930200.0380.009xy下一頁上一頁主 頁海港系統(tǒng)卸載貨物的模擬海港系統(tǒng)卸載貨物的模擬 產(chǎn)生均勻隨機數(shù)的區(qū)間 對應(yīng)時間范圍 線性插值函數(shù) 00.009y1520 x555.615.0000 xy0.0090.038y2030 x344.816.8966xy下一頁上一頁主 頁海港系統(tǒng)卸載貨物的模擬海港系統(tǒng)卸載貨物的模擬 MATLAB編程如下：(jge.m)a=rand;if 0=a&a0.009 b=555.6*a+15;else

12、if 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; 下一頁上一頁主 頁海港系統(tǒng)卸載貨物的模擬海港系統(tǒng)卸載貨物的模擬 1. 如果港口可以同時卸載多艘船只的貨物，將是一個多窗口的服務(wù)問題，如何進行模擬？2. 如果海港系統(tǒng)不僅僅考慮卸載貨物問題，如還考慮裝載貨物問題，或卸載貨物以后再裝載等問題，又應(yīng)該如何分析？思考：下一頁上一頁主 頁可靠性問題可靠性問題 一

13、設(shè)備上有三個相同的軸承，每個軸承正常工作壽命為隨機變量，概率分布如下：壽命/h1000110012001300140015001600170018001900概率0.100.130.250.130.090.120.020.060.050.05有軸承損壞設(shè)備停止工作檢修工準備開始更換部件，稱為一個延遲時間，它也是隨機變量，分布如下：延遲時間 /min 5 1015概 率0.60.30.1范例二下一頁上一頁主 頁主要費用：1、設(shè)備停工損失費： 5 元/分鐘；2、檢修工人的工時費：12 元/小時；3、軸承的成本費： 16元/個更換軸承所需要的時間： 一個 兩個 三個 20 30 40 （min）范例

14、二下一頁上一頁主 頁問題：現(xiàn)在有兩種方案：方案一：損壞一個軸承只更換一個軸承；方案二：一旦有軸承損壞就全部更換；試通過計算機模擬對以上兩種方案做出評價。 隨機數(shù)怎樣產(chǎn)生？ 模擬時選用時間步長法還是 事件步長法？范例二下一頁上一頁主 頁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; elseif 0.7a&a=0.82 l=1500; elseif 0.82a&

15、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下一頁上一頁主 頁關(guān)于隨機數(shù)的產(chǎn)生見 Lrnd.m（零件壽命） Yrnd.m（延遲時間） 方案一的數(shù)學(xué)模型：kekao1.m目標函數(shù) min c =Ui / T（三個軸承同時損壞）（兩個軸承同時損壞）只有一個軸承損壞）3166040125)40(2166030125)30(166020125)20(iiiitttU其中 損失費工時費成本費ti 表示延遲時間 下一頁上一頁主 頁方案二的數(shù)學(xué)模型：目標函數(shù) min c =Ui /

16、 T其中31612325)40(iitU10000)(1iniiLt方案二的情況比較單一。Kekao3.m模擬框圖如下：范例二下一頁上一頁主 頁數(shù)據(jù)初始化T100000產(chǎn)生下一個事件發(fā)生時刻根據(jù)不同方案確定更換策略產(chǎn)生延遲時間計算當前更換費用累積費用產(chǎn)生新軸承的壽命時間確定下一個事件發(fā)生的時刻系統(tǒng)時刻跳轉(zhuǎn)到下一個事件發(fā)生時刻輸出總費用 停止否是范例二下一頁上一頁主 頁c1=5;c2=12;c3=16; (kekao1.m)g1=20/60;g2=30/60;g3=40/60;Lrnd; % Random function Lrnd gets a random life(three)yrnd;

17、% 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),lm(3),t(j);lm=sort(L);yrnd;endU;t; zU=sum(U) C=zU/t(50)

18、計算結(jié)果：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)計算結(jié)果：ZU=14800（元），C=0.6271。（元/小時）范例二下一頁上一頁主 頁計算結(jié)果：方案一：c = 0.3307 元/小時；方案二：c = 0.5486（0.7234）元/小時；顯然方案一優(yōu)于方案二。范例二返 回下一頁上一頁主 頁練習 某港口有