計算機模擬-排隊論

1、word可編輯當排隊系統(tǒng)的到達間隔時間和效勞時間的概率分布很復雜時，或不能用公式給出 時，那么就不能用解析法求解。這就需用隨機模擬法求解，現舉例說明。例 1 設某倉庫前有一卸貨場，貨車一般是夜間到達，白天卸貨，每天只能卸貨2車，假設一天內到達數超過2車，那么就推遲到次日卸貨。根據表3所示的數據，貨車到達數的概率分布相對頻率平均為1.5車/天，求每天推遲卸貨的平均車數。表 3 到達車數的概率到達車數012345 6概率0.230.300.300.10.050.020.00解這是單效勞臺的排隊系統(tǒng)，可驗證到達車數不服從泊松分布，效勞時間也不服從指數分布這是定長效勞時間。隨機模擬法首先要求事件能按歷

2、史的概率分布規(guī)律出現。模擬時產生的隨機數與事件的對應關系如表4。表 4 到達車數的概率及其對應的隨機數到達車數概率累積概率對應的隨機數0123450.230.300.300.10.050.020.230.530.830.930.981.000 x < 0.230.23 x < 0.530.53 x < 0.830.83 x < 0.930.93 x < 0.980.98 x 1.00我們用 a1 表示產生的隨機數，a2 表示到達的車數，a3 表示需要卸貨車數，a4 表 示實際卸貨車數，a5 表示推遲卸貨車數。編寫程序如下：clear rand('state

3、',sum(100*clock);n=50000;m=2 a1=rand(n,1);a2=a1; %a2初始化a2(find(a1<0.23)=0; a2(find(0.23<=a1&a1<0.53)=1; a2(find(0.53<=a1&a1<0.83)=2; a2(find(0.83<=a1&a1<0.93),1)=3; a2(find(0.93<=a1&a1<0.98),1)=4; a2(find(a1>=0.98)=5;a3=zeros(n,1);a4=zeros(n,1);a5=ze

4、ros(n,1); %a2初始化a3(1)=a2(1);if a3(1)<=m a4(1)=a3(1);a5(1)=0;elsea4(1)=m;a5(1)=a2(1)-m;endfor i=2:na3(i)=a2(i)+a5(i-1); if a3(i)<=m a4(i)=a3(i);a5(i)=0; else a4(i)=m;a5(i)=a3(i)-m;endend a=a1,a2,a3,a4,a5;sum(a)/n m = 2ans = 0.4985 1.4909 2.3782 1.4909 0.8874 例2銀行方案安置自動取款機，A型機的價格是B型機的2倍，而A型機的性能平

5、均效勞率也是B型機的2倍，問應該購置1臺 A 型機還是2臺 B 型機。為了通過模擬答復這類問題，作如下具體假設，顧客平均每分鐘到達1位， A 型 機的平均效勞時間為0.9分鐘，B 型機為1.8分鐘，顧客到達間隔和效勞時間都服從指數分布，2臺B型機采取M/M/2模型排一隊，用前100名顧客第 1 位顧客到達時取款機前為空的平均等待時間為指標，對A型機和B型機分別作1000次模擬，進行比擬。理論上已經得到，A型機和B型機前100名顧客的平均等待時間分別為1(100)=4.13，2(100)=3.70，即 B 型機優(yōu)。對于M/M/1模型，記第k位顧客的到達時刻為ck，離開時刻為gk，等待時間為wk，

6、它們很容易根據已有的到達間隔ik和效勞時間sk按照以下的遞推關系得到w1 = 0，設c1,g1：ck=ck1+ik，gk=max(ck,gk1)+ sk，wk=max(0, gk1 ck)， k=2,3,L。在模擬A型機時，用cspan表示到達間隔時間，sspan表示效勞時間，ctime表示到達時間，gtime表示離開時間，wtime表示等待時間。我們總共模擬了m次，每次n個顧客。程序如下：ticrand('state',sum(100*clock);n=100;m=1000;mu1=1;mu2=0.9; for j=1:m cspan=exprnd(mu1,1,n); ssp

7、an=exprnd(mu2,1,n); ctime(1)=cspan(1); gtime(1)=ctime(1)+sspan(1); wtime(1)=0; for i=2:n ctime(i)=ctime(i-1)+cspan(i); gtime(i)=max(ctime(i),gtime(i-1)+sspan(i); wtime(i)=max(0,gtime(i-1)-ctime(i); end result1(j)=sum(wtime)/n;endresult_1=sum(result1)/mtoc result_1 = 4.0467Elapsed time is 0.445770 se

8、conds. 類似地，模擬B型機的程序如下：ticrand('state',sum(100*clock);n=100;m=1000;mu1=1;mu2=1.8;for j=1:m cspan=exprnd(mu1,1,n); sspan=exprnd(mu2,1,n); ctime(1)=cspan(1); ctime(2)=ctime(1)+cspan(2); gtime(1:2)=ctime(1:2)+sspan(1:2); wtime(1:2)=0;flag=gtime(1:2); for i=3:n ctime(i)=ctime(i-1)+cspan(i); gtime

9、(i)=max(ctime(i),min(flag)+sspan(i); wtime(i)=max(0,min(flag)-ctime(i); flag=max(flag),gtime(i); end result2(j)=sum(wtime)/n;end result_2=sum(result2)/m toc result_2 = 3.7368Elapsed time is 1.453880 seconds. 可以用下面的程序與上面的程序比擬了解編程的效率問題。ticclearrand('state',sum(100*clock);n=100;m=1000;mu1=1;mu2

10、=0.9;for j=1:m ctime(1)=exprnd(mu1); gtime(1)=ctime(1)+exprnd(mu2); wtime(1)=0; for i=2:n ctime(i)=ctime(i-1)+exprnd(mu1); gtime(i)=max(ctime(i),gtime(i-1)+exprnd(mu2); wtime(i)=max(0,gtime(i-1)-ctime(i); end result(j)=sum(wtime)/n;end result=sum(result)/mtoc result = 4.2162Elapsed time is 3.854620

11、seconds. 黃河小浪底調水調沙問題5.1問題的提出2004年6月至7月黃河進行了第三次調水調沙試驗，特別是首次由小浪底、三門峽和萬家寨三大水庫聯(lián)合調度，采用接力式防洪預泄放水，形成人造洪峰進行調沙試驗獲得成功。整個試驗期為20多天，小浪底從6月19日開始預泄放水，直到7月13日恢復正常供水結束。小浪底水利工程按設計攔沙量為75.5 億m3，在這之前，小浪底共積泥沙達14.15億t。這次調水調沙試驗一個重要目的就是由小浪底上游的三門峽和萬家寨水庫泄洪，在小浪底形成人造洪峰，沖刷小浪底庫區(qū)沉積的泥沙，在小浪底水庫開閘泄洪以后，從6月27日開始三門峽水庫和萬家寨水庫陸續(xù)開閘放水，人造洪峰于29

12、日先后到達小浪底，7月3日到達最大流量2700m3/s，使小浪底水庫的排沙量也不斷地增加。表7是由小浪底觀測站從6月29日到7月10日檢測到的試驗數據。表 7 觀測數據日期6.296.307.17.27.37.4時間8:0020:008:0020:008:0020:008:0020:008:0020:008:0020:00水流量180019002100220023002400250026002650270027202650含沙量326075859098100102108112115116日期7.57.67.77.87.97.10時間8:0020:008:0020:008:0020:008:00

13、20:008:0020:008:0020:00水流量26002500230022002000185018201800175015001000900含沙量11812011810580605030262085現在，根據試驗數據建立數學模型研究下面的問題：1給出估計任意時刻的排沙量及總排沙量的方法；2確定排沙量與水流量的關系。5.2模型的建立與求解給定的觀測時刻是等間距的，以6月29日零時刻開始計時，那么各次觀測時刻離開始時刻6月29日零時刻的時間分別為ti =3600(12i4)，i=1,2,L,24，其中計時單位為秒 。第1次觀測的時刻t1=28800，最后一次觀測的時刻t24=1022400

14、。記第ii= 1,2,L,24次觀測時水流量為vi，含沙量為ci，那么第i次觀測時的排沙量為yi=ci*vi 。有關的數據見表8。表 8 插值數據對應關系 單位：排沙量為 kg節(jié)點12345678時刻2880072000115200158400202200244800288000331200排沙量57600114000157500187000207000235200250000265200節(jié)點910111213141516時刻374400417600460800504000547200590400633600676800排沙量2862003024003128003074003068003000

15、00271400231000節(jié)點1718192021222324時刻7200007632008064008496008928009360009792001022400排沙量1600001110009100054000455003000080004500對于問題1，根據所給問題的試驗數據，要計算任意時刻的排沙量，就要確定出排沙量隨時間變化的規(guī)律，可以通過插值來實現?？紤]到實際中的排沙量應該是時間的連續(xù)函數，為了提高模型的精度，我們采用三次樣條函數進行插值。利用 MATLAB 函數，求出三次樣條函數，得到排沙量y=y(t)與時間的關系，然后進行積分，就可以得到總的排沙量最后求得總的排沙量為1.84

16、4 ×109 t，計算的 Matlab 程序如下：clc,clearload data.txt%data.txt 按照原始數據格式把水流量和排沙量排成4行12列l(wèi)iu=data(1,3,:);liu=liu'liu=liu(:);sha=data(2,4,:);sha=sha'sha=sha(:);y=sha.*liu;y=y'i=1:24;t=(12*i-4)*3600;t1=t(1);t2=t(end);pp=csape(t,y);xsh=pp.coefs%求得插值多項式的系數矩陣，每一行是一個區(qū)間上多項式的系數。TL=quadl(tt)ppval(pp,

17、tt),t1,t2)也可以利用 3 次 B 樣條函數進行插值，求得總的排沙量也為1.844 ×109 t，計算 的 Matlab 程序如下：clc,clearload data.txt%data.txt 按照原始數據格式把水流量和排沙量排成4行12列l(wèi)iu=data(1,3,:);liu=liu'liu=liu(:);sha=data(2,4,:);sha=sha'sha=sha(:);y=sha.*liu;y=y'i=1:24;t=(12*i-4)*3600;t1=t(1);t2=t(end);pp=spapi(4,t,y)%三次 B 樣條pp2=fn2fm

18、(pp,'pp')%把 B 樣條函數轉化為 pp 格式TL=quadl(tt)fnval(pp,tt),t1,t2)對于問題2，研究排沙量與水量的關系，從試驗數據可以看出，開始排沙量是隨著水流量的增加而增長，而后是隨著水流量的減少而減少。顯然，變化規(guī)律并非是線性 的關系，為此，把問題分為兩局部，從開始水流量增加到最大值 2720m3/s即增長的過 程為第一階段，從水流量的最大值到結束為第二階段，分別來研究水流量與排沙量的關系。畫出排沙量與水流量的散點圖。畫散點圖的程序如下：load data.txt liu=data(1,3,:); liu=liu'liu=liu(:)

19、; sha=data(2,4,:); sha=sha'sha=sha(:); y=sha.*liu;subplot(1,2,1), plot(liu(1:11),y(1:11),'*')subplot(1,2,2), plot(liu(12:24),y(12:24),'*')從散點圖可以看出，第一階段根本上是線性關系，第二階段準備依次用二次、三次、 四次曲線來擬合，看哪一個模型的剩余標準差小就選取哪一個模型。最后求得第一階段排沙量y與水流量v之間的預測模型為y = 250.5655v 373384.4661第二階段的預測模型為一個四次多項式。y =2.7

20、693×107v4+0.0018v34.092v2+3891.0441v1.32262749668×106計算的 Matlab 程序如下：clc, clearload data.txt%data.txt 按照原始數據格式把水流量和排沙量排成4行12列l(wèi)iu=data(1,3,:); liu=liu' liu=liu(:); sha=data(2,4,:); sha=sha' sha=sha(:); y=sha.*liu;%以下是第一階段的擬合format long enihe1_1=polyfit(liu(1:11),y(1:11),1)%擬合一次多項式，系

21、數排列從高次冪到低次冪nihe1_2=polyfit(liu(1:11),y(1:11),2)yhat1_1=polyval(nihe1_1,liu(1:11);%求預測值yhat1_2=polyval(nihe1_2,liu(1:11);%以下求誤差平凡和與剩余標準差cha1_1=sum(y(1:11)-yhat1_1).2); rmse1_1=sqrt(cha1_1/9)cha1_2=sum(y(1:11)-yhat1_2).2); rmse1_2=sqrt(cha1_2/8)%以下是第二階段的擬合for j=1:3 str1=char('nihe2_' int2str(j

22、) '=polyfit(liu(12:24),y(12:24),' int2str(j+1) ')'); eval(str1) str2=char('yhat2_' int2str(j) '=polyval(nihe2_' int2str(j) ',liu(12:24);'); eval(str2) str3=char('cha2_' int2str(j) '=sum(y(12:24)-yhat2_' int2str(j) ').2);' 'rmse2_' int2str(j) '=sqrt(cha2_' int2str(j) '/(11-j)'); eval(str3)end 例:當實際問題可以用馬爾