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

2、分布規(guī)律出現(xiàn)。模擬時(shí)產(chǎn)生的隨機(jī)數(shù)與事件的對應(yīng)關(guān)系如表4。表 4 到達(dá)車數(shù)的概率及其對應(yīng)的隨機(jī)數(shù)到達(dá)車數(shù)概率累積概率對應(yīng)的隨機(jī)數(shù)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 表示產(chǎn)生的隨機(jī)數(shù)，a2 表示到達(dá)的車數(shù)，a3 表示需要卸貨車數(shù)，a4 表 示實(shí)際卸貨車數(shù)，a5 表示推遲卸貨車數(shù)。編寫程序如下：clear rand(state,sum(100*clock);n=50000;m=2 a1=ran

3、d(n,1);a2=a1; %a2初始化a2(find(a10.23)=0; a2(find(0.23=a1&a10.53)=1; a2(find(0.53=a1&a10.83)=2; a2(find(0.83=a1&a10.93),1)=3; a2(find(0.93=a1&a1=0.98)=5;a3=zeros(n,1);a4=zeros(n,1);a5=zeros(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);

4、 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銀行計(jì)劃安置自動取款機(jī)，已知A型機(jī)的價(jià)格是B型機(jī)的2倍，而A型機(jī)的性能平均服務(wù)率也是B型機(jī)的2倍，問應(yīng)該購置1臺 A 型機(jī)還是2臺 B 型機(jī)。為了通過模擬回答這類問題，作如下具體假設(shè)，顧客平均每分鐘到達(dá)1位， A 型 機(jī)的平均服務(wù)時(shí)間為0.9分鐘，B 型機(jī)為1.8分鐘，顧客到達(dá)間隔和服務(wù)時(shí)間都服從指數(shù)分布，2臺B型機(jī)

5、采取M/M/2模型（排一隊(duì)），用前100名顧客（第 1 位顧客到達(dá)時(shí)取款機(jī)前為空）的平均等待時(shí)間為指標(biāo)，對A型機(jī)和B型機(jī)分別作1000次模擬，進(jìn)行比較。理論上已經(jīng)得到，A型機(jī)和B型機(jī)前100名顧客的平均等待時(shí)間分別為1(100)=4.13，2(100)=3.70，即 B 型機(jī)優(yōu)。對于M/M/1模型，記第k位顧客的到達(dá)時(shí)刻為ck，離開時(shí)刻為gk，等待時(shí)間為wk，它們很容易根據(jù)已有的到達(dá)間隔ik和服務(wù)時(shí)間sk按照以下的遞推關(guān)系得到（w1 = 0，設(shè)c1,g1已知）：ck=ck1+ik，gk=max(ck,gk1)+ sk，wk=max(0, gk1 ck)， k=2,3,L。在模擬A型機(jī)時(shí)，用cs

6、pan表示到達(dá)間隔時(shí)間，sspan表示服務(wù)時(shí)間，ctime表示到達(dá)時(shí)間，gtime表示離開時(shí)間，wtime表示等待時(shí)間。我們總共模擬了m次，每次n個(gè)顧客。程序如下：ticrand(state,sum(100*clock);n=100;m=1000;mu1=1;mu2=0.9; for j=1:m cspan=exprnd(mu1,1,n); sspan=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(

7、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 seconds. 類似地，模擬B型機(jī)的程序如下：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(m

8、u2,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(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

9、_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=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),

10、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 seconds. 黃河小浪底調(diào)水調(diào)沙問題5.1問題的提出2004年6月至7月黃河進(jìn)行了第三次調(diào)水調(diào)沙試驗(yàn)，特別是首次由小浪底、三門峽和萬家寨三大水庫聯(lián)合調(diào)度，采用接力式防洪預(yù)泄放水，形成人造洪峰進(jìn)行調(diào)沙試驗(yàn)獲得成功。整個(gè)試驗(yàn)期為20多天，小浪底從6月19日開始預(yù)泄放水，直到7月13日恢復(fù)正

11、常供水結(jié)束。小浪底水利工程按設(shè)計(jì)攔沙量為75.5 億m3，在這之前，小浪底共積泥沙達(dá)14.15億t。這次調(diào)水調(diào)沙試驗(yàn)一個(gè)重要目的就是由小浪底上游的三門峽和萬家寨水庫泄洪，在小浪底形成人造洪峰，沖刷小浪底庫區(qū)沉積的泥沙，在小浪底水庫開閘泄洪以后，從6月27日開始三門峽水庫和萬家寨水庫陸續(xù)開閘放水，人造洪峰于29日先后到達(dá)小浪底，7月3日達(dá)到最大流量2700m3/s，使小浪底水庫的排沙量也不斷地增加。表7是由小浪底觀測站從6月29日到7月10日檢測到的試驗(yàn)數(shù)據(jù)。表 7 觀測數(shù)據(jù)日期6.296.307.17.27.37.4時(shí)間8:0020:008:0020:008:0020:008:0020:008

12、:0020:008:0020:00水流量180019002100220023002400250026002650270027202650含沙量326075859098100102108112115116日期7.57.67.77.87.97.10時(shí)間8:0020:008:0020:008:0020:008:0020:008:0020:008:0020:00水流量26002500230022002000185018201800175015001000900含沙量11812011810580605030262085現(xiàn)在，根據(jù)試驗(yàn)數(shù)據(jù)建立數(shù)學(xué)模型研究下面的問題：（1）給出估計(jì)任意時(shí)刻的排沙量及總排沙量

13、的方法；（2）確定排沙量與水流量的關(guān)系。5.2模型的建立與求解已知給定的觀測時(shí)刻是等間距的，以6月29日零時(shí)刻開始計(jì)時(shí)，則各次觀測時(shí)刻（離開始時(shí)刻6月29日零時(shí)刻的時(shí)間）分別為ti =3600(12i4)，i=1,2,L,24，其中計(jì)時(shí)單位為秒 。第1次觀測的時(shí)刻t1=28800，最后一次觀測的時(shí)刻t24=1022400 。記第i（i= 1,2,L,24）次觀測時(shí)水流量為vi，含沙量為ci，則第i次觀測時(shí)的排沙量為yi=ci*vi 。有關(guān)的數(shù)據(jù)見表8。表 8 插值數(shù)據(jù)對應(yīng)關(guān)系 單位：排沙量為 kg節(jié)點(diǎn)12345678時(shí)刻2880072000115200158400201600244800288

14、000331200排沙量57600114000157500187000207000235200250000265200節(jié)點(diǎn)910111213141516時(shí)刻374400417600460800504000547200590400633600676800排沙量286200302400312800307400306800300000271400231000節(jié)點(diǎn)1718192021222324時(shí)刻7200007632008064008496008928009360009792001022400排沙量1600001110009100054000455003000080004500對于問題（1），根據(jù)所

15、給問題的試驗(yàn)數(shù)據(jù)，要計(jì)算任意時(shí)刻的排沙量，就要確定出排沙量隨時(shí)間變化的規(guī)律，可以通過插值來實(shí)現(xiàn)?？紤]到實(shí)際中的排沙量應(yīng)該是時(shí)間的連續(xù)函數(shù)，為了提高模型的精度，我們采用三次樣條函數(shù)進(jìn)行插值。利用 MATLAB 函數(shù)，求出三次樣條函數(shù)，得到排沙量y=y(t)與時(shí)間的關(guān)系，然后進(jìn)行積分，就可以得到總的排沙量最后求得總的排沙量為1.844 109 t，計(jì)算的 Matlab 程序如下：clc,clearload data.txt%data.txt 按照原始數(shù)據(jù)格式把水流量和排沙量排成4行12列l(wèi)iu=data(1,3,:);liu=liu;liu=liu(:);sha=data(2,4,:);sha=s

16、ha;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%求得插值多項(xiàng)式的系數(shù)矩陣，每一行是一個(gè)區(qū)間上多項(xiàng)式的系數(shù)。TL=quadl(tt)ppval(pp,tt),t1,t2)也可以利用 3 次 B 樣條函數(shù)進(jìn)行插值，求得總的排沙量也為1.844 109 t，計(jì)算 的 Matlab 程序如下：clc,clearload data.txt%data.txt 按照原始數(shù)據(jù)格式把水流量和排沙量排成4行12列l(wèi)iu=data(1,3,:);liu=liu;liu=l

17、iu(:);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(pp,pp)%把 B 樣條函數(shù)轉(zhuǎn)化為 pp 格式TL=quadl(tt)fnval(pp,tt),t1,t2)對于問題（2），研究排沙量與水量的關(guān)系，從試驗(yàn)數(shù)據(jù)可以看出，開始排沙量是隨著水流量的增加而增長，而后是隨著水流量的減少而減少。顯然，變化規(guī)律并非是線性 的關(guān)系，為此，把問題分為兩部分，從開始水流量增加到最大值 2720m3/s（即

18、增長的過 程）為第一階段，從水流量的最大值到結(jié)束為第二階段，分別來研究水流量與排沙量的關(guān)系。畫出排沙量與水流量的散點(diǎn)圖。畫散點(diǎn)圖的程序如下：load data.txt liu=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(12:24),*)從散點(diǎn)圖可以看出，第一階段基本上是線性關(guān)系，第二階段準(zhǔn)備依次用二次、三次、 四次曲線來擬合，看

19、哪一個(gè)模型的剩余標(biāo)準(zhǔn)差小就選取哪一個(gè)模型。最后求得第一階段排沙量y與水流量v之間的預(yù)測模型為y = 250.5655v 373384.4661第二階段的預(yù)測模型為一個(gè)四次多項(xiàng)式。y =2.7693107v4+0.0018v34.092v2+3891.0441v1.32262749668106計(jì)算的 Matlab 程序如下：clc, clearload data.txt%data.txt 按照原始數(shù)據(jù)格式把水流量和排沙量排成4行12列l(wèi)iu=data(1,3,:); liu=liu; liu=liu(:); sha=data(2,4,:); sha=sha; sha=sha(:); y=sha.

20、*liu;%以下是第一階段的擬合format long enihe1_1=polyfit(liu(1:11),y(1:11),1)%擬合一次多項(xiàng)式，系數(shù)排列從高次冪到低次冪nihe1_2=polyfit(liu(1:11),y(1:11),2)yhat1_1=polyval(nihe1_1,liu(1:11);%求預(yù)測值yhat1_2=polyval(nihe1_2,liu(1:11);%以下求誤差平凡和與剩余標(biāo)準(zhǔn)差cha1_1=sum(y(1:11)-yhat1_1).2); rmse1_1=sqrt(cha1_1/9)cha1_2=sum(y(1:11)-yhat1_2).2); rmse

21、1_2=sqrt(cha1_2/8)%以下是第二階段的擬合for j=1:3 str1=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);); 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 例:當(dāng)實(shí)際問題可以用馬爾可夫鏈來描述時(shí)，首先要確定它的狀態(tài)空間及參數(shù)集合，然后確定它的一步轉(zhuǎn)移概率。關(guān)于這一概率的確定，可以由問題的內(nèi)在規(guī)律得到，也可以由過去經(jīng)驗(yàn)給出，還可以根據(jù)觀測數(shù)據(jù)來估計(jì)。例4 某計(jì)算機(jī)機(jī)房的一臺計(jì)算機(jī)經(jīng)常出故障，研究者每隔15分鐘觀察一次計(jì)算機(jī)的運(yùn)行狀態(tài)，收集了24小時(shí)的數(shù)據(jù)（共作97次觀察）。用1表示正常狀態(tài)，用