計算機模擬_排隊論

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

2、布規(guī)律出現(xiàn).模擬時產生的隨機數(shù)與事件的對應關系如表4.表4到達車數(shù)的概率及其對應的隨機數(shù)到達車數(shù)概率累積概率對應的隨機數(shù)00.23P0.230<x<0.2310.300.530.23<x<0.5320.300.830.53<x<0.8330.10.930.83<x<0.9340.050.980.93<x<0.9850.021.000.98<x<1.00我們用a1表示產生的隨機數(shù),a2表示到達的車數(shù),a3表示需要卸貨車數(shù),a4表示實際卸貨車數(shù),a5表示推遲卸貨車數(shù).編寫程序如下：clearrand('state

3、9;,sum(100*clock);n=50000;m=2a1=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=zeros(n,1);%

4、a2初始化a3(1)=a2(1);ifa3(1)<=ma4(1)=a3(1);a5(1)=0;elsea4(1)=m;a5(1)=a2(1)-m;endfori=2:na3(i)=a2(i)+a5(i-1);ifa3(i)<=ma4(i)=a3(i);a5(i)=0;elsea4(i)=m;a5(i)=a3(i)-m;endenda=a1,a2,a3,a4,a5;sum(a)/nm=2ans=0.49851.49092.37821.49090.8874例2銀行方案安置自動取款機,A型機的價格是B型機的2倍,而A型機的性能-平均效勞率也是B型機的2倍,問應該購置1臺A型機還是2臺B型

5、機.為了通過模擬答復這類問題,作如下具體假設,顧客平均每分鐘到達1位,A型機的平均效勞時間為0.9分鐘,B型機為1.8分鐘,顧客到達間隔和效勞時間都服從指數(shù)分布,2臺B型機采取M/M/2模型(排一隊),用前100名顧客(第1位顧客到達時取款機前為空)的平均等待時間為指標,對A型機和B型機分別作1000次模擬,進行比擬.理論上已經得到,A型機和B型機前100名顧客的平均等待時間分別為(100)=4.13,必100)=3.70,即B型機優(yōu).對于M/M/1模型,記第k位顧客的到達時刻為ck,離開時刻為gk,等待時間為wk,它們很容易根據已有的到達間隔ik和效勞時間sk根據以下的遞推關系得到(w1=0

6、,設c1,g1)：ck=ck-1+ik,gk=max(ck,gk-1)+sk,wk=max(0,gk-1-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;forj=1:mcspan=exprnd(mu1,1,n);sspan=exprnd(mu2,1,n);ctime(1)=cspan(1);gtime(

7、1)=ctime(1)+sspan(1);wtime(1)=0;fori=2:nctime(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);endresult1(j)=sum(wtime)/n;endresult_1=sum(result1)/mtocresult_1=4.0467Elapsedtimeis0.445770seconds.類似地,模擬B型機的程序如下：ticrand('state',sum(100*clock);

8、n=100;m=1000;mu1=1;mu2=1.8;forj=1:mcspan=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);fori=3:nctime(i)=ctime(i-1)+cspan(i);gtime(i)=max(ctime(i),min(flag)+sspan(i);wtime(i)=max(0,min(flag)-ctime(i);fl

9、ag=max(flag),gtime(i);endresult2(j)=sum(wtime)/n;endresult_2=sum(result2)/mtocresult_2=3.7368Elapsedtimeis1.453880seconds.可以用下面的程序與上面的程序比擬了解編程的效率問題.ticclearrand('state',sum(100*clock);n=100;m=1000;mu1=1;mu2=0.9;forj=1:mctime(1)=exprnd(mu1);gtime(1)=ctime(1)+exprnd(mu2);wtime(1)=0;fori=2:ncti

10、me(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);endresult(j)=sum(wtime)/n;endresult=sum(result)/mtocresult=4.2162Elapsedtimeis3.854620seconds.黃河小浪底調水調沙問題5.1問題的提出2004年6月至7月黃河進行了第三次調水調沙試驗,特別是首次由小浪底、三門峽和萬家寨三大水庫聯(lián)合調度,采用接力式防洪預泄放水,形成人造洪峰進行調沙試驗獲得成

11、功.整個試驗期為20多天,小浪底從6月19日開始預泄放水,直到7月13日恢復正常供水結束.小浪底水利工程按設計攔沙量為75.5億m3,在這之前,小浪底共積泥沙達14.15億t.這次調水調沙試驗一個重要目的就是由小浪底上游的三門峽和萬家寨水庫泄洪,在小浪底形成人造洪峰,沖刷小浪底庫區(qū)沉積的泥沙,在小浪底水庫開閘泄洪以后,從6月27日開始三門峽水庫和萬家寨水庫陸續(xù)開閘放水,人造洪峰于29日先后到達小浪底,7月3日到達最大流量2700m3/s,使小浪底水庫的排沙量也不斷地增加.表7是由小浪底觀測站從6月29日至I7月10日檢測到的試驗數(shù)據.表7觀測數(shù)據日期6.296.307.17.27.37.4時間

12、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:0020:008:0020:008:0020:00水流量26002500230022002000185018201800175015001000900含沙量11812011810580605030262085現(xiàn)在,根據試

13、驗數(shù)據建立數(shù)學模型研究下面的問題：(1)給出估計任意時刻的排沙量及總排沙量的方法；(2)確定排沙量與水流量的關系.5.2模型的建立與求解給定的觀測時刻是等間距的,以6月29日零時刻開始計時,那么各次觀測時刻(離開始時刻6月29日零時刻的時間)分別為ti=3600(12i-4),i=1,2,L,24,其中計時單位為秒.第1次觀測的時刻t1=28800,最后一次觀測的時刻t24=1022400.記第i(i=1,2,L,24)次觀測時水流量為vi,含沙量為ci,那么第i次觀測時的排沙量為yi=ci*vi.有關的數(shù)據見表8.表8插值數(shù)據對應關系單位：排沙量為kg節(jié)點12345678討亥U2880072

14、000115200158400202100244800288000331200排沙量57600114000157500187000207000235200250000265200節(jié)點910111213141516時刻374400417600460800504000547200590400633600676800排沙量286200302400312800307400306800300000271400231000節(jié)點1718192021222324時刻720000763200806400849600892800936000979200102240C琲沙量160000111000910005400

15、0455003000080004500對于問題(1),根據所給問題的試驗數(shù)據,要計算任意時刻的排沙量,就要確定出排沙量隨時間變化的規(guī)律,可以通過插值來實現(xiàn).考慮到實際中的排沙量應該是時間的連續(xù)函數(shù),為了提升模型的精度,我們采用三次樣條函數(shù)進行插值.利用MATLAB函數(shù),求出三次樣條函數(shù),得到排沙量y=y(t)與時間的關系,然后進行積分,就可以得到總的排沙量Z=p4y(tyit%最后求得總的排沙量為1.844X109t,計算的Matlab程序如下：clc,clearloaddata.txt%data.txt根據原始數(shù)據格式把水流量和排沙量排成4行12列l(wèi)iu=data(1,3,:);liu=li

16、u'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%充得插值多項式的系數(shù)矩陣,每一行是一個區(qū)間上多項式的系數(shù).TL=quadl(tt)ppval(pp,tt),t1,t2)也可以利用3次B樣條函數(shù)進行插值,求得總的排沙量也為1.844X109t,計算的Matlab程序如下：clc,clearloaddata.txt%data.txt根據原始數(shù)據格式把水流量和排沙量

17、排成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=H1);t2=t(end);pp=spapi(4,t,y)味次B樣條pp2=fn2fm(pp,'pp')%EB樣條函數(shù)轉化為pp格式TL=quadl(tt)fnval(pp,tt),t1,t2)對于問題(2),研究排沙量與水量的關系,從試驗數(shù)據可以看出,開始排沙量是隨著水流量的增加而增長,而后是隨著水流量的減少而減少.顯然

18、,變化規(guī)律并非是線性的關系,為此,把問題分為兩局部,從開始水流量增加到最大值2720m3/s(即增長的過程)為第一階段,從水流量的最大值到結束為第二階段,分別來研究水流量與排沙量的關系.畫出排沙量與水流量的散點圖.畫散點圖的程序如下：loaddata.txtliu=data(1,3,:);liu=liu'liu=liu(:);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

19、(12:24),'*')從散點圖可以看出,第一階段根本上是線性關系,第二階段準備依次用二次、三次、四次曲線來擬合,看哪一個模型的剩余標準差小就選取哪一個模型.最后求得第一階段排沙量y與水流量v之間的預測模型為y=250.5655v-373384.4661第二階段的預測模型為一個四次多項式.y=-2.7693X0-7v4+0.0018v3-4.092v2+3891.0441v-1.32262749668X106計算的Matlab程序如下：clc,clearloaddata.txt%data.txt根據原始數(shù)據格式把水流量和排沙量排成4行12列l(wèi)iu=data(1,3,:);liu

20、=liu'liu=liu(:);sha=data(2,4,:);sha=sha'sha=sha(:);y=sha.*liu;%以下是第一階段的擬合formatlongenihe1_1=polyfit(liu(1:11),y(1:11),1)%以合二次多項式,系數(shù)排列從高次幕到低次幕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).

21、A2);rmse1_1=sqrt(cha1_1/9)cha1_2=sum(y(1:11)-yhat1_2).A2);rmse1_2=sqrt(cha1_2/8)蛆刀是第二階段的擬合一一一forj=1:3str1=char('nihe2_'int2str(j)'=polyfit(liu(12:24),y(12:24),'int2str(j+1)')');eval(str1)str2=char('yhat2_'int2str(j)'=polyval(nihe2_'int2str(j)',liu(12:24);&

22、#39;);eval(str2)str3=char('cha2_'int2str(j)'=sum(y(12:24)-yhat2_'int2str(j)').A2);''rmse2_'int2str(j)'=sqrt(cha2_'int2str(j)'/(11-j)');eval(str3)end例:當實際問題可以用馬爾可夫鏈來描述時,首先要確定它的狀態(tài)空間及參數(shù)集合,然后確定它的一步轉移概率.關于這一概率確實定,可以由問題的內在規(guī)律得到,也可以由過去經驗給出,還可以根據觀測數(shù)據來估計.例4某計算機機房的一臺計算機經常出故障,研究者每隔15分鐘觀察一次計算機的運行狀態(tài),收集了24小時的數(shù)據(共作97次觀察).用1表示正常狀態(tài),用0表示不正常狀態(tài),所得的數(shù)據序列如下：1110