數(shù)學(xué)建模協(xié)會(huì)第二講(模擬)

基于matlab的計(jì)算機(jī)模擬舒興明TelQ：117562750一、模擬是什么模擬(也稱仿真)：就是用計(jì)算機(jī)程序在計(jì)算機(jī)上模仿各種實(shí)際系統(tǒng)的運(yùn)行過程，并通過計(jì)算了解系統(tǒng)隨時(shí)間變化的行為或特性。計(jì)算機(jī)模擬：是在已經(jīng)建立的數(shù)學(xué)、邏輯模型之上，通過計(jì)算機(jī)實(shí)驗(yàn)，對(duì)一個(gè)系統(tǒng)按照一定的決策原則或作業(yè)規(guī)則，由一個(gè)狀態(tài)變換為另外一個(gè)狀態(tài)的行為進(jìn)行描述和分析。實(shí)際問題數(shù)學(xué)、邏輯模型計(jì)算機(jī)模型數(shù)學(xué)、計(jì)算機(jī)解實(shí)際解二、為何要模擬（1）實(shí)際系統(tǒng)建立之前，要對(duì)系統(tǒng)的行為或結(jié)果進(jìn)行分析研究；（2）有些真實(shí)系統(tǒng)做實(shí)驗(yàn)會(huì)影響系統(tǒng)運(yùn)行，例如，在生產(chǎn)中任意改變工藝系數(shù)可能導(dǎo)致廢品，在經(jīng)濟(jì)活動(dòng)中隨意將一個(gè)決策付諸行動(dòng)會(huì)導(dǎo)致經(jīng)濟(jì)混亂；（3）在系統(tǒng)上做多次試驗(yàn)，很難保證每次試驗(yàn)的操作條件相同，因而對(duì)實(shí)驗(yàn)結(jié)果好壞很難作出正確的判斷；（4）當(dāng)人是系統(tǒng)的一部分時(shí)，他的行為往往實(shí)驗(yàn)結(jié)果有所影響，這時(shí)，最好進(jìn)行模擬研究；（5）實(shí)驗(yàn)時(shí)間太長，費(fèi)用太大，或者有危險(xiǎn)，使得試驗(yàn)不容易進(jìn)行；（6）有些系統(tǒng)一旦建立起來后無法復(fù)原，例如，建立大型企業(yè)，要分析社會(huì)和經(jīng)濟(jì)效益，不能用建立起來試試看的辦法。三、哪些問題適合模擬研究（1）難以用數(shù)學(xué)公式表示的系統(tǒng)，或者沒有求解數(shù)學(xué)模型的有效方法；（2）雖然可以用解析的方法解決問題，但是問題的分析與計(jì)算過于復(fù)雜，這時(shí)計(jì)算機(jī)仿真可能提供簡(jiǎn)單可行的求解方法；（3）希望在較短時(shí)間內(nèi)觀察到系統(tǒng)的發(fā)展全過程，以估計(jì)某參數(shù)對(duì)系統(tǒng)行為的影響；（4）難以在實(shí)際環(huán)境中進(jìn)行試驗(yàn)和觀察，計(jì)算機(jī)仿真是唯一的方法；（5）需要對(duì)系統(tǒng)或過程進(jìn)行長期運(yùn)行的比較，從大量方案中尋找最優(yōu)。四、模擬的分類模擬是系統(tǒng)狀態(tài)隨時(shí)間而變化的動(dòng)態(tài)寫照，因此，通常時(shí)間是模擬的主要自變量，其它的變量為因變量。（1）按照模擬過程中因變量的變化情況，可以將模擬分為離散、連續(xù)、混合3種類型；（2）如果采用模擬計(jì)算機(jī)、采用數(shù)字計(jì)算機(jī)以及聯(lián)合使用則分為模擬仿真、數(shù)字仿真以及混合仿真；（3）根據(jù)仿真變量的特征分為隨機(jī)模擬仿真和模糊模擬仿真。五、隨機(jī)模擬的基礎(chǔ)基礎(chǔ)知識(shí)：線性代數(shù)部分向量矩陣的表述；微積分的基本積分和微分表述；概率論的基本分布；基本統(tǒng)計(jì)知識(shí)等等。模擬（蒙特卡羅方法）基本原理：概率論的大數(shù)定律獨(dú)立隨機(jī)序列的和按概率收斂于這些隨機(jī)序列的數(shù)學(xué)期望的和（常數(shù)），即六、matlab統(tǒng)計(jì)工具箱中產(chǎn)生隨機(jī)數(shù)的函數(shù)1、概率密度函數(shù)pdf分布類型調(diào)用格式備注二項(xiàng)分布binoy=binopdf(x,n,p)返回參數(shù)為p，n的二項(xiàng)分布在x處的分布律值指數(shù)分布expy=exppdf(x,mu)返回參數(shù)為mu的指數(shù)分布在x處的密度函數(shù)值Χ2分布chi2y=chi2pdf(x,v)返回參數(shù)為v的卡方分布在x處的密度函數(shù)值F分布fy=fpdf(x,v1,v2)返回參數(shù)為v1,v2的F分布在x處密度函數(shù)值幾何分布geoy=geopdf(x,p)返回參數(shù)為p的幾何分布在x處的分布律值超幾何分布hygey=hygepdf(x,M,k,n)返回參數(shù)為M,k,n的超幾何分在x處的分布律值常見隨機(jī)變量的密度函數(shù)表分布類型調(diào)用格式備注正態(tài)分布normnormpdf(x,mu,sig)返回參數(shù)為mu和sig的正態(tài)分布在x處的密度值泊松分布poisspoisspdf(x,lamda)返回參數(shù)為lamda的泊松分布在x處的分布律值T分布ttpdf(x,v)返回參數(shù)為v的t分布在x處的密度函數(shù)值連續(xù)均勻分布unifunifpdf(x,a,b)返回在[a,b]上服從均勻分布的密度函數(shù)值離散均勻分布unidunidpdf(x,N)返回1~N的均勻分布概率分布律值密度函數(shù)續(xù)注意：二項(xiàng)分布(bino)，幾何分布(geo)，超幾何分布(hyge),泊松分布(poiss)，指數(shù)分布(exp)，正態(tài)分布(norm)，卡方分布(chi2)，t分布(t)，F(xiàn)分布(f)，均勻連續(xù)分布(unif)，均勻離散分布(unid)。2、常見的分布函數(shù)cdf分布類型調(diào)用格式備注二項(xiàng)分布y=binocdf(x,n,p)返回參數(shù)為p，n的二項(xiàng)分布在x處的分布值指數(shù)分布y=expcdf(x,mu)返回參數(shù)為mu的指數(shù)分布在x處的分布函數(shù)值Χ2分布y=chi2cdf(x,v)返回參數(shù)為v的卡方分布在x處的分布函數(shù)值F分布y=fcdf(x,v1,v2)返回參數(shù)為v1,v2的F分布在x處分布函數(shù)值幾何分布y=geocdf(x,p)返回參數(shù)為p的幾何分布在x處的分布值超幾何分布y=hygecdf(x,M,k,n)返回參數(shù)為M,k,n的超幾何分在x處的分布值常見隨機(jī)變量的分布函數(shù)數(shù)表分布類型調(diào)用格式備注正態(tài)分布normcdf(x,mu,sig)返回參數(shù)為mu和sig的正態(tài)分布在x處的分布值泊松分布poisscdf(x,lamda)返回參數(shù)為lamda的泊松分布在x處的分布值T分布tcdf(x,v)返回參數(shù)為v的t分布在x處的分布函數(shù)值連續(xù)均勻分布unifcdf(x,a,b)返回在[a,b]上服從均勻分布的分布函數(shù)值離散均勻分布unidcdf(x,N)返回1~N的均勻分布概率分布值分布函數(shù)續(xù)將分布函數(shù)的概率調(diào)用p=cdf(‘name’,x,a1,a2)調(diào)過來用，就是逆分布函數(shù)x=icdf(‘name’,p,a1,a2)。>>x1=icdf('norm',0.95,0,1)x1=1.644853626951473>>x2=icdf('t',0.95,10)x2=1.812461122811677>>x3=icdf('exp',0.95,3)x3=8.987196820661970>>x4=icdf('chi2',0.95,10)x4=18.3070380532751433、模擬隨機(jī)數(shù)的產(chǎn)生格式1：x=random(‘name’,A1,A2,A3,m,n)產(chǎn)生參數(shù)為A1，A2，A3的name分布的m行n列的矩陣。>>x=random('norm',0,1,3,4)x=0.53770.8622-0.43362.76941.83390.31880.3426-1.3499-2.2588-1.30773.57843.0349>>y=random('unif',1.2,2.2,3,4)y=2.15721.34191.99221.23571.68541.62182.15952.04912.00032.11571.85572.1340產(chǎn)生3行4列的參數(shù)為0和1的正態(tài)分布隨機(jī)數(shù)產(chǎn)生在[1.2,2.2]上服從均勻分布的隨機(jī)數(shù)格式2：x=namernd(A1,A2,m,n)產(chǎn)生m行，n列的服從參數(shù)為A1，A2的name分布的隨機(jī)數(shù)。>>x=trnd(3,4)x=0.6880-1.9018-0.6416-1.14291.0967-1.64313.3100-2.16210.9627-1.1476-0.7321-1.18490.07281.10551.62900.6073>>y=trnd(3,3,4)y=-4.3874-1.88741.0409-1.5729-1.08610.8758-2.57105.65671.75710.1385-0.25441.2256產(chǎn)生自由度為3的服從t分布的隨機(jī)數(shù)（4階方陣）。產(chǎn)生3行4列自由度為3的服從t分布的隨機(jī)數(shù)。七、隨機(jī)樣本的簡(jiǎn)單描述性統(tǒng)計(jì)量1、中心趨勢(shì)（位置）度量樣本的幾何平均m=geomean(X)樣本的調(diào)和平均m=harmmean(X)樣本的算術(shù)平均m=mean(X)樣本的中位數(shù)m=median(X)樣本的剔出極端數(shù)據(jù)的平均m=trimmean(X,percent)>>X=random('unif',1,8,10,1);>>[geomean(X),harmmean(X),mean(X),median(X),trimmean(X,0.2)]ans=4.59033.74005.24016.06785.24012、分散度量樣本四分位數(shù)間距m=iqr(X)樣本數(shù)據(jù)的平均絕對(duì)偏差m=mad(X)即mean(abs(X-mean(X)))樣本極差m=range(X)樣本方差m=var(X)樣本標(biāo)準(zhǔn)差m=std(X)>>X=random('unid',20,10,1);>>m=[iqr(X),mad(X),range(X),var(X),std(X)]m=9.00005.160018.000035.73335.9777八、隨機(jī)模擬計(jì)算積分例1如下圖所示，在正方形內(nèi)有1/4單位圓。向正方形內(nèi)投小石頭，假設(shè)每次都能夠投進(jìn)正方形內(nèi)且等可能落在正方形內(nèi)任何一點(diǎn)。問，小石頭落在1/4單位圓內(nèi)(包含邊界)的概率多大？01x1y分析：假設(shè)投入正方形內(nèi)的石頭有n塊，有k塊落入了1/4單位圓內(nèi)。P為小石頭落入1/4單位圓內(nèi)的概率。那么根據(jù)Bernoull(伯努利)大數(shù)定律，有即，當(dāng)實(shí)驗(yàn)次數(shù)n充分大時(shí)，頻率和概率之差小于任意數(shù)ε的概率趨于1。而另外一方面由幾何概型有在實(shí)際操作中，實(shí)驗(yàn)次數(shù)n不可能趨于無窮大，所以有(數(shù)學(xué)模型)對(duì)于估計(jì)只有不斷重復(fù)做實(shí)驗(yàn)，這種試驗(yàn)可以具體去操作，(均勻投石塊，然后數(shù)數(shù)，這樣需要較高成本)。也可以讓計(jì)算機(jī)去重復(fù)試驗(yàn)，但是需要將數(shù)學(xué)模型轉(zhuǎn)化為計(jì)算機(jī)模擬模型(讓計(jì)算機(jī)完成均勻投石塊，自動(dòng)計(jì)數(shù)，也需要成本，但成本比起前者，幾乎可以忽略。)。用計(jì)算機(jī)模擬投石塊過程和步驟如下：1、自動(dòng)生成隨機(jī)點(diǎn)[0,1]x[0,1]，模擬石塊在正方形內(nèi)的任意位置，用(xi,yi)表示，共n個(gè)點(diǎn)；2、判別(xi,yi)是否滿足xi2+yi2≤1，即判別石塊是否落在1/4單位圓內(nèi)，共k個(gè)點(diǎn)滿足；3、整理、統(tǒng)計(jì)模擬結(jié)果，用4k/n估計(jì)π。符號(hào)設(shè)置：符號(hào)說明n試驗(yàn)次數(shù)（投石塊次數(shù)）k石塊落入1/4圓內(nèi)次數(shù)(xi,yi)第i塊石頭落點(diǎn)坐標(biāo)pai圓周率計(jì)算機(jī)模擬的流程圖如下初始化：k=0,i=1i>n是Pai=4k/n否在[0,1]x[0,1]上產(chǎn)生隨機(jī)點(diǎn)(xi,yi)r=xi^2+yi^2r<=1是k=k+1i=i+1否設(shè)計(jì)一個(gè)m文件，專門針對(duì)給定n，模擬計(jì)算pai的估計(jì)值functionpai=fpai(n)ifn<=0error('nhastobeapositiveinteger');endk=0;i=1;whilei<=nx=random('unif',0,1,1,2);r=sum(x.^2);ifr<=1k=k+1;endi=i+1;endpai=4*k/n;模擬的穩(wěn)定性描述>>clear>>tic>>n=50:10:5000;>>fork=1:length(n)pai(k)=fpai(n(k));end>>plot(n,pai,'*')TocElapsedtimeis223.847533seconds.>>mean(pai)ans=3.1405??！值得注意的是，每次模擬即或程序相同，不同人不同計(jì)算時(shí)間，結(jié)果一般不同(因?yàn)殡S機(jī))。functionpai=fpai1(n)ifn<=0error('nhastobeapositiveinteger');endk=0;x=random('unif',0,1,n,2);fork1=1:nx1=x(k1,:);r=sum(x1.^2).^0.5;ifr<=1k=k+1;endendpai=4*k/n;把while語句改為for語句，同時(shí)不用每次都調(diào)用random函數(shù)，m文件如下>>clearticn=50:10:5000;fork=1:length(n)pai(k)=fpai1(n(k));endplot(n,pai,'*')tocElapsedtimeis7.372741seconds.時(shí)間大大縮短了??！計(jì)算定積分分析：若隨機(jī)變量X的概率分布密度是P(x)(a≤x≤b)，則隨機(jī)變量Y=f(X)的數(shù)學(xué)期望為當(dāng)X是[a,b]上的均勻分布時(shí)，有將p(x)代入上面的期望算子，有例2而另外一方面，根據(jù)大數(shù)定理，設(shè)y1,y2,…,yn為來自總體Y的一組樣本，有其中所以有設(shè)x1,x2,…,xn為來自服從[a,b]上均勻分布總體X的一組樣本。且yi=f(xi)，i=1,2,…,n。那么計(jì)算的步驟如下：1、在[2，3]上抽樣x1,x2,…,xn；2、計(jì)算3、整理估計(jì)定積分；為了對(duì)比模擬結(jié)果，先把這個(gè)積分求出來：>>int('(1+x^2)^0.5',2,3)ans=2.6947540047794075420843777595616然后根據(jù)上面的模擬步驟設(shè)計(jì)一個(gè)m文件functionI=jifen1(n)x=random('unif',2,3,1,n);y=(1+x.^2).^0.5;I=(3-2)/n*sum(y);>>clear>>tic;n=50:10:5000;fork=1:length(n)I(k)=jifen1(n(k));endplot(n,I,'*');tocElapsedtimeis0.334083seconds.例3求在空間坐標(biāo)系下，兩個(gè)圓柱體相交部分的體積。兩個(gè)圓柱體的柱面方程分別為根據(jù)二重積分的知識(shí)，相交部分的體積為xyz=16/3計(jì)算的模擬步驟如下（1）在空間[0,1]x[0,1]x[0,1]上產(chǎn)生均勻隨機(jī)點(diǎn)(xi,yi,zi)，i=1,2,…,n；（2）判斷是否成立，若有k個(gè)；（3）用8k/n來估計(jì)V的值。把以上步驟編寫成m文件functionV=erchjfen(n)x=random('unif',0,1,n,3);k=0;fork1=1:nr1=x(k1,1)^2+x(k1,2)^2;r2=x(k1,1)^2+x(k1,3)^2;if(r1<=1)&(r2<=1)k=k+1;endendV=8*k/n;>>clear>>tic>>n=20:20:4000;>>fork=1:length(n)v(k)=erchjfen(n(k));end;plot(n,v,'*');tocElapsedtimeis72.712497seconds.九、評(píng)估事件發(fā)生的概率有些隨機(jī)變量，由于不容易找到（或者不存在）準(zhǔn)確的概率分布律或者密度函數(shù)，不能算出此類隨機(jī)變量表達(dá)的隨機(jī)的事件的概率，需要用到模擬。其中X1服從均勻分布U[2,5]，X2服從指數(shù)分布exp(3)，X3和X4分別服從正態(tài)分布N(3,2)和N(1,1)。例4求下面事件表達(dá)的概率1、在總體X=(X1,X2,X3,X4)中抽樣，得到樣本(x1k,x2k,x3k,x4k)，k=1,2,…,N；2、判斷隨機(jī)樣本是否滿足不等式組，假設(shè)滿足不等式組的有n個(gè)；3、用n/N估計(jì)p。（大數(shù)定律作保障，或者點(diǎn)估計(jì)原理）步驟：按照上述步驟編寫m文件functionp=gailv1(n)x1=random('unif',2,5,n,1);x2=random('exp',3,n,1);x3=random('norm',3,2,n,1);x4=random('norm',1,1,n,1);k1=0;fork=1:nr1=x1(k)+x2(k)^2;r2=x3(k)+x4(k)^2;if(r1>=3)&(r2<=9)k1=k1+1;endendp=k1/n;>>n=20:20:4000;>>tic>>fork=1:length(n)p(k)=gailv1(n(k));end;plot(n,p,'*');tocElapsedtimeis38.580145seconds.>>p(length(n))ans=0.8375十、求隨機(jī)變量的臨界值（上下分為數(shù)）求使下式成立的最大f值步驟：1、產(chǎn)生N個(gè)獨(dú)立隨機(jī)序列(X1k,X2k,X3k),k=1,2,…,N；2、計(jì)算fk=f(X1k,X2k,X3k),k=1,2,…,N;3、從從大到小的序列f1,f2,…,fN中取f0.8N作為f的估計(jì)值，或從小到大的f0.2N作為f的估計(jì)值。例5其中，X1服從均勻分布U[1,3]，X2服從指數(shù)分布exp(1)，X3服從正態(tài)分布N(2,1)。根據(jù)上面的流程，編寫m文件如下functionf=linjie(n)x1=random('unif',1,3,n,1);x2=random('exp',1,n,1);x3=random('norm',2,1,n,1);fork=1:ny(k)=x1(k)+x2(k)^2+x3(k)^3;endz=sort(y);k1=round(0.2*n);f=z(k1);>>clearn=50:50:10000;ticfork=1:length(n)f(k)=linjie(n(k));end;plot(n,f,'*');tocElapsedtimeis11.335281seconds.十一、報(bào)童問題有個(gè)報(bào)童每天從報(bào)社批發(fā)報(bào)紙，每份報(bào)紙1元，然后再以零售價(jià)1.5元每份售出，如果賣不出，則退回報(bào)社，退回價(jià)格為每份1.2元。某個(gè)社區(qū)某個(gè)月記錄在案的每天報(bào)紙的需求量如下表需求量50454038353028天數(shù)2367732頻率0.06670.10.20.23330.23330.10.0667問：報(bào)童每天該批發(fā)多少分報(bào)紙才是最佳策略？合理假設(shè)（1）報(bào)童兜售報(bào)紙的社區(qū)比較穩(wěn)定，即每個(gè)月的報(bào)紙的需求量的分布抽樣差不多；（2）把上表給的結(jié)果看成社區(qū)每天報(bào)紙的需求分布律。例6符號(hào)假設(shè)a訂購價(jià)（元/份）b零售價(jià)（元/份）c退回價(jià)（元/份）r每天報(bào)紙的需求量P(r)每天報(bào)紙的需求量的抽樣分布律x報(bào)童報(bào)紙的批發(fā)量（份）y報(bào)童每天的獲利（元）建立模型每天報(bào)紙的需求量r是個(gè)隨機(jī)變量，報(bào)童每天的盈利由r和批發(fā)量x決定，也是隨機(jī)變量，即目標(biāo)1每天盈利期往最大maxE(y)數(shù)學(xué)模型1（數(shù)學(xué)期望模型）r~P(r)模型1求解將報(bào)童每天盈利函數(shù)編寫成m文件functionEy=baotong(x)globalabcPrifb<a|c<aerror('bandcmustbelargerthana');endn=size(Pr,2);r=Pr(1,:);p=Pr(2,:);fork=1:nifx<=r(k)y(k)=(b-a)*x;elsey(k)=(b-a)*r(k)-(c-a)*(x-r(k));endendEy=sum(y.*p);在命令窗口調(diào)用上述函數(shù)，計(jì)算得>>clear>>globalabcPr>>a=1;b=1.5;c=1.2;>>Pr=[50454038353028;0.06670.10.20.23330.23330.10.0667];>>x=1:1:60;>>fork=1:length(x)Ey(k)=baotong(x(k));end>>plot(x,Ey,'*')>>Ey1=max(Ey);>>m=find(Ey==Ey1)m=40最佳批發(fā)為40份報(bào)紙。數(shù)學(xué)模型2（最低盈利機(jī)會(huì)最大模型）目標(biāo)2r~P(r)每天至少賺16元的概率最大（機(jī)會(huì)最大）將每天賺16元的機(jī)會(huì)計(jì)算程序編制為m文件：functionp=baotongjh(x,ybar)globalabcPrifb<a|c<aerror('bandcmustbelargerthana');endn=size(Pr,2);r=Pr(1,:);pr=Pr(2,:);fork=1:nifx<=r(k)Ey(k)=(b-a)*x;elseEy(k)=(b-a)*r(k)-(c-a)*(x-r(k));endendpy=Ey>ybar;p=sum(py.*pr);>>clear>>globalabcPr>>a=1;b=1.5;c=1.2;>>Pr=[50454038353028;0.06670.10.20.23330.23330.10.0667];>>ybar=15;>>x=1:60;>>fork=1:60p(k)=baotongjh(x(k),ybar);end>>plot(x,p,'*')由此可見，當(dāng)報(bào)紙批發(fā)量在30~45之間時(shí)，報(bào)童每天可以賺15元錢的概率最大有83%。十二、利用模擬法求解非線性規(guī)劃例7用隨機(jī)模擬方法求解給定的規(guī)劃問題求解此問題的流程圖如下初始化n=0,A=[]n>1000YN在[0,10]上產(chǎn)生均勻隨機(jī)數(shù)x1，x2=(10-x1)/3，y=x1+2x2y>2NYn=n+1,A=[A;[x1,x2]]z=0,k=1:1000,xx=[0,0]x1=A(k,1),x2=A(k,2)z1=f(x1,x2)z1>zYz=z1Nk=k+1xx=[x1,x2]依照上述流程編寫一個(gè)求此非線性規(guī)劃最大值的m文件：function[z,xx]=nlpf(N)n=0;A=[];whilen<=Nx1=random('unif',0,10);x2=(10-x1)/3;y=x1+2*x2;ify>2n=n+1;A=[A;[x1,x2]];endendz=0;xx=[0,0];fork=1:Nx11=A(k,1);x22=A(k,2);z1=-2*x11^2-x22^2+x11*x22+8*x11+3*x22;ifz1>zz=z1;xx=[x11,x22];endend>>clear>>A=[];>>N=20:5:400;>>fork=1:length(N)[z,xx]=nlpf(N(k));A=[A;z,xx];end>>plot(N,A(:,1),'*')十三、一個(gè)排隊(duì)系統(tǒng)的模擬例7一個(gè)排隊(duì)系統(tǒng)的仿真某店只有一個(gè)收款臺(tái)，顧客到達(dá)收款臺(tái)的時(shí)間間隔服從均值為4.5的負(fù)指數(shù)分布，每個(gè)顧客的服務(wù)時(shí)間服從均值為3.2，標(biāo)準(zhǔn)差為0.6的正態(tài)分布，對(duì)100位顧客去收款臺(tái)繳款的排隊(duì)過程進(jìn)行仿真。分析：排隊(duì)系統(tǒng)涉及的參考變量為：顧客的等待時(shí)間Wq，顧客的逗留時(shí)間W,排隊(duì)長Lq,隊(duì)長L，服務(wù)臺(tái)的空閑時(shí)間I，服務(wù)臺(tái)的繁忙時(shí)間B，排隊(duì)系統(tǒng)的狀態(tài)概率P。最后需要統(tǒng)計(jì)出有關(guān)指標(biāo)的平均值均方差。合理假設(shè)1、進(jìn)入系統(tǒng)的顧客沒有人因?yàn)椴辉敢獾榷x開系統(tǒng)；2、先到先服務(wù)；模擬步驟：1、產(chǎn)生初始事件表，并從初始事件表中找出最靠前的事件；2、啟動(dòng)仿真時(shí)鐘；3、判別當(dāng)前事件是哪一個(gè)事件，如果是顧客到達(dá)事件，則啟動(dòng)顧客到達(dá)子程序；如果是服務(wù)事件，則啟動(dòng)服務(wù)結(jié)束子程序；4、重復(fù)2-3，直到仿真完畢；5、整理、分析并輸出結(jié)果。顧客到達(dá)的子程序1、產(chǎn)生下一個(gè)顧客到來的時(shí)刻，記入事件表；2、判別服務(wù)臺(tái)是否忙s=1?，若是，則顧客則排隊(duì)長L=L+1;否則，服務(wù)臺(tái)空閑,L=0;3、產(chǎn)生一個(gè)服務(wù)結(jié)束時(shí)刻，記入事件表；4、統(tǒng)計(jì)所需的數(shù)據(jù)；服務(wù)結(jié)束子程序1、結(jié)束服務(wù)，產(chǎn)生結(jié)束時(shí)間，記入事件表；2、判別是否有顧客排隊(duì)，L=0?，若L=0，則服務(wù)臺(tái)空閑；否則，排隊(duì)長L=L-1;3、產(chǎn)生服務(wù)結(jié)束時(shí)刻，記入事件表；4、統(tǒng)計(jì)所需的數(shù)據(jù)；產(chǎn)生初始事件表從事件表中找出最靠前的事件仿真時(shí)鐘步進(jìn)是哪一類事件顧客到達(dá)服務(wù)結(jié)束產(chǎn)生下一個(gè)顧客到達(dá)紀(jì)錄，記入事件表結(jié)束服務(wù)，產(chǎn)生結(jié)束時(shí)間，記入事件表出納是否忙？出納空閑，s=0顧客排隊(duì)，Lq=Lq+1否是產(chǎn)生服務(wù)結(jié)束時(shí)刻，記入事件表是否顧客排隊(duì)隊(duì)長減1，Lq-1出納空閑，s=0否產(chǎn)生服務(wù)結(jié)束時(shí)刻，記入事件表是統(tǒng)計(jì)所需數(shù)據(jù)仿真是否完畢是結(jié)束程序否按照流程，編寫模擬程序function[ndaoda,nlikai,t,LLq,LL,BI,tdaoda,tfuwu]=paiduisimu(n,a,b,c)t=0；LL=[];LLq=[];BI=[];tdaoda=[];tfuwu=[];ndaoda=0;nlikai=0;L=0;Lq=0;t1=exprnd(a);t2=+inf;s=0;whilenlikai<n;ift1<t2t=t1;r1=exprnd(a);t1=t+r1;ifs==1L=L+1;Lq=max(L-1,0);elseL=1;s=1;Lq=max(0,L-1);r2=normrnd(b,c);t2=t+r2;endLL=[LL,L];LLq=[LLq,Lq];tdaoda=[tdaoda,t];

ndaoda=ndaoda+1;nlikai=nlikai;BI=[BI,s];elseift1>t2t=t2;tfuwu=[tfuwu,t];ifLq>0s=1;L=L-1;Lq=Lq-1;r2=normrnd(b,c);t2=t+r2;elses=0;L=0;Ls=max(L-1,0);t2=+inf;endLL=[LL,L];LLq=[LLq,Lq];BI=[BI,s];nlikai=nlikai+1;ndaoda=ndaoda;

elseift1==t2t=t1;tdaoda=[tdaoda,t];tfuwu=[tfuwu,t];ndaoda=ndaoda+1;nlikai=nlikai+1;L=L;Lq=Lq;s=s;LL=[LL,L];LLq=[LLq,Lq];BI=[BI,s];t1=t+exprnd(a);t2=t+normrnd(b,c);endend>>[ndaoda,nlikai,t,LLq,LL,BI,tdaoda,tfuwu]=paiduisimu(n,a,b,c);>>[mean(LLq),mean(LL),mean(BI)]ans=0.96021.76620.8060由此可見，這個(gè)服務(wù)窗口一直處以繁忙狀態(tài)（空閑率)的概率為1-0.8060;平均排隊(duì)長（mean(L)=0.9602人，不用再增加一個(gè)窗口來緩解壓力。十四、航班延誤問題某航空公司經(jīng)營A，B，C三個(gè)城市的航線，這些航線每天班次起飛與到達(dá)時(shí)間如下表所示。設(shè)飛機(jī)在機(jī)場(chǎng)停留的損失費(fèi)大致與停留時(shí)間成正比，又每架飛機(jī)從降落到下班起飛至少需2小時(shí)準(zhǔn)備時(shí)間，試決定一個(gè)使停留費(fèi)用損失為最小的分派飛行方案。另外，根據(jù)抽樣資料，三個(gè)城市A，B，C之間的航班延誤時(shí)間分別服從負(fù)指數(shù)分布：城市A到城市B的航班延誤時(shí)間服從均值為1的指數(shù)分布；城市A到城市C的航班延誤時(shí)間服從均值為1.5的指數(shù)分布；城市B到城市C的航班延誤時(shí)間服從均值為0.8的指數(shù)分布；城市B到城市A的航班延誤時(shí)間服從均值為1.2的指數(shù)分布；城市C到城市A的航班延誤時(shí)間服從均值為0.8的指數(shù)分布；城市C到城市B的航班延誤時(shí)間服從均值為0.5的指數(shù)分布；航班號(hào)

起飛城市

起飛時(shí)間

到達(dá)城市

到達(dá)時(shí)間101A9：00B

12：00102

A10：00

B

13：00103

A

22：00

C

2：00（次日)104

B4：00A7：00

105

B11：00A

14：00

106

C7：00

A

11：00

107

B13：00

C

18：00

108C

15：00B

20：00各個(gè)航班往返城市及其時(shí)間如下表：合理假設(shè)（1）城市A，B，C的到達(dá)與起飛航班數(shù)目一樣，每個(gè)城市都用到達(dá)航班的飛機(jī)作為起飛航班的飛機(jī);（2）沒有多余的飛機(jī)可用；（3）只算一天內(nèi)的延誤時(shí)間；（3）飛機(jī)到達(dá)后2個(gè)小時(shí)任何時(shí)刻都可以起飛；（4）先到先服務(wù)原則；（5）概率小于0.0001為小概率事件，不考慮。符號(hào)設(shè)置xijxij=1,城市A的第i到達(dá)作為第j離開，否則xij=0yijyij=1,城市B的第i到達(dá)作為第j離開，否則yij=0zijzij=1,城市C的第i到達(dá)作為第j離開，否則zij=0t1ijA城市第i到達(dá)作為第j離開的等待時(shí)間t2ijB城市第i到達(dá)作為第j離開的等待時(shí)間t3ijC城市第i到達(dá)作為第j離開的等待時(shí)間TEF第i城市到第j城市的延誤時(shí)間E=A,B,C,F=A,B,C,E≠F同時(shí)根據(jù)假設(shè)（5），所有航班延誤超過某個(gè)時(shí)刻后是小概率事件，認(rèn)為不發(fā)生。1、可行解的產(chǎn)生方法k<=nNstop初始化：y=[],k=1Y產(chǎn)生隨機(jī)數(shù)a1~U(0,1)a1<=1/3x=[100]YNa1>2/3Yx=[001]Nx=[010]產(chǎn)生隨機(jī)數(shù)a2~U(0,1)a2<=0.5確定x的0所在列第二行1的位置確定剩1行1列1的位置y=[y;x]k=k+1functiony=hbzhipai3(n)y=[];k=1;whilek<=na1=random('unif',0,1);ifa1<=1/3x=[1,0,0];a2=random('unif',0,1);ifa2<=0.5x=[x;[0,1,0]];x=[x;[001]];elsex=[x;[001]];x=[x;[010]];endelseifa1>2/3x=[001];a2=random('unif',0,1);ifa2<=0.5x=[x;[100]];x=[x;[010]];elsex=[x;[010]];x=[x;[100]];endelsex=[010];a2=random('unif',0,1);ifa2<=0.5x=[x;[100]];x=[x;[001]];elsex=[x;[001]];x=[x;[100]];endendy=[y;x];k=k+1;end但是，根據(jù)這個(gè)問題的特殊性，城市A的航班安排只有廈門6個(gè)可行解：x1=[100;010;001];x2=[100;001;010];x3=[010;001;100];x4=[010;100;001];x5=[001;010;100];x6=[001;100;010];2、模擬在給定延誤情況下的等待時(shí)間航班號(hào)\航班號(hào)101(9)102(10)103(22)B：104(7)2315B：105(14)19208C：106(11)222311城市A的航班間的不延誤等待時(shí)間t1ij2.1航班延誤時(shí)A城市的各個(gè)等待時(shí)間表由于B和C到A城市的航班延誤時(shí)間服從的分布分別為當(dāng)TBA和TCA取定時(shí)，航班飛機(jī)的等待時(shí)間就跟著確定，如下流程初始化：tba1~e(1.2)tba1<=1Yt1(104,102)=3-tba1t1(104,103)=15-tba1t1(104,101)=25-tba1Ntba1>1,tba1<=13Yt1(104,103)=15-tba1t1(104,101)=24-tba1t1(104,102)=25-tba1Nt1(104,101)=26-tba1t1(104,102)=17-tba1t1(104,103)=39-tba1初始化：tba2~e(1.2)tba2<=6Yt1(105,102)=20-tba2t1(105,103)=8-tba2t1(105,101)=19-tba2Ntba2>6,tba2<=17Yt1(105,103)=24-tba2t1(105,101)=19-tba2t1(105,102)=20-tba2Nt1(105,101)=43-tba2t1(105,102)=20-tba2t1(105,103)=32-tba2tba2>17,tba2<=18YNt1(105,101)=43-tba2t1(105,102)=44-tba2t1(105,103)=32-tba2初始化：tba3~e(0.8)tba3<=9Yt1(106,102)=23-tba3t1(106,103)=11-tba3t1(106,101)=22-tba3Ntba3>9,tba3<=20Yt1(106,103)=33-tba3t1(106,101)=22-tba3t1(106,102)=23-tba3Nt1(106,101)=46-tba3t1(106,102)=23-tba3t1(106,103)=35-tba3tba3>20,tba3<=21YNt1(106,103)=35-tba3t1(106,101)=46-tba3t1(106,102)=47-tba3按照這三個(gè)流程，編寫M文件elseif(tba2>17)&(tba2<=18)t12=[432032]-tba2;elset12=[434432]-tba2;endtba3=random('exp',0.8);iftba3<=9t13=[222311]-tba3;elseif(tba3>9)&(tba3<=20)t13=[222333]-tba3;elseif(tba3>20)&(tba3<=21)t13=[462335]-tba3;elset13=[464735]-tba3;endtr1=[t11;t12;t13];t1=[t1;tr1];k=k+1;endfunctiont1=hbwaita(m)t1=[];k=1;whilek<=mtba1=random('exp',1.2);iftba1<=1t11=[25315]-tba1;elseif(tba1>1)&(tba1<=13)t11=[242515]-tba1;elset11=[262739]-tba1;endtba2=random('exp',1.2);iftba2<=6t12=[19208]-tba2;elseif(tba2>6)&(tba2<=17)t12=[192024]-tba2;2.2尋找城市A的最小費(fèi)用航班安排模型按照普通指派問題，數(shù)學(xué)模型為步驟如下：（1）針對(duì)每一個(gè)可行解x；（2）產(chǎn)生m個(gè)航班延誤時(shí)間t1；（3）計(jì)算可行解x在各種延誤時(shí)間t1的總的安排時(shí)間，并把這m個(gè)時(shí)間平均，記錄其為x對(duì)應(yīng)的目標(biāo)；（4）比較每個(gè)解所對(duì)應(yīng)的目標(biāo)，選最小的最為參考最優(yōu)解。2.3尋找最優(yōu)解按照上述流程編寫m文件：function[tt]=zuiyoua1(x,m)tt1=[];fork1=1:mt1=hbwaita(1);yt=x.*t1;tt11=sum(sum(yt));tt1=[tt1;tt11];endtt=mean(tt1);在命令窗口計(jì)算，結(jié)果顯示如下>clearm=10000;x1=[100;010;001];x2=[100;001;010];x3=[010;001;100];x4=[010;100;001];x5=[001;010;100];x6=[001;100;010];tt1=zuiyoua1(x1,m);tt2=zuiyoua1(x2,m);tt3=zuiyoua1(x3,m);tt4=zuiyoua1(x4,m);tt5=zuiyoua1(x5,m);tt6=zuiyoua1(x6,m);[tt1,tt2,tt3,tt4,tt5,tt6]ans=52.366652.480439.511739.058953.812553.8071計(jì)算結(jié)果顯示，可行解x4比較好！X3其實(shí)也不錯(cuò)。對(duì)于城市B和城市C，可以用相同的方法模擬計(jì)算結(jié)果。DVD在線租賃

2005年國賽B題舒興明TeQ117562750原問題考慮如下的在線DVD租賃問題。顧客繳納一定數(shù)量的月費(fèi)成為會(huì)員，訂購DVD租賃服務(wù)。會(huì)員對(duì)哪些DVD有興趣，只要在線提交訂單，網(wǎng)站就會(huì)通過快遞的方式盡可能滿足要求。會(huì)員提交的訂單包括多張DVD，這些DVD是基于其偏愛程度排序的。網(wǎng)站會(huì)根據(jù)手頭現(xiàn)有的DVD數(shù)量和會(huì)員的訂單進(jìn)行分發(fā)。每個(gè)會(huì)員每個(gè)月租賃次數(shù)不得超過2次，每次獲得3張DVD。會(huì)員看完3張DVD之后，只需要將DVD放進(jìn)網(wǎng)站提供的信封里寄回（郵費(fèi)由網(wǎng)站承擔(dān)），就可以繼續(xù)下次租賃。請(qǐng)考慮以下問題：（1）網(wǎng)站正準(zhǔn)備購買一些新的DVD，通過問卷調(diào)查1000個(gè)會(huì)員，得到了愿意觀看這些DVD的人數(shù)（表1給出了其中5種DVD的數(shù)據(jù)）。此外，歷史數(shù)據(jù)顯示，60%的會(huì)員每月租賃DVD兩次，而另外的40%只租一次。假設(shè)網(wǎng)站現(xiàn)有10萬個(gè)會(huì)員，對(duì)表1中的每種DVD來說，應(yīng)該至少準(zhǔn)備多少張，才能保證希望看到該DVD的會(huì)員中至少50%在一個(gè)月內(nèi)能夠看到該DVD？如果要求保證在三個(gè)月內(nèi)至少95%的會(huì)員能夠看到該DVD呢？（2）表2中列出了網(wǎng)站手上100種DVD的現(xiàn)有張數(shù)和當(dāng)前需要處理的1000位會(huì)員的在線訂單（表2的數(shù)據(jù)格式示例如下表2，具體數(shù)據(jù)請(qǐng)從/mcm05/problems2005c.asp下載），如何對(duì)這些DVD進(jìn)行分配，才能使會(huì)員獲得最大的滿意度？請(qǐng)具體列出前30位會(huì)員（即C0001~C0030）分別獲得哪些DVD。（3）繼續(xù)考慮表2，并假設(shè)表2中DVD的現(xiàn)有數(shù)量全部為0。如果你是網(wǎng)站經(jīng)營管理人員，你如何決定每種DVD的購買量，以及如何對(duì)這些DVD進(jìn)行分配，才能使一個(gè)月內(nèi)95%的會(huì)員得到他想看的DVD，并且滿意度最大？（4）如果你是網(wǎng)站經(jīng)營管理人員，你覺得在DVD的需求預(yù)測(cè)、購買和分配中還有哪些重要問題值得研究？請(qǐng)明確提出你的問題，并嘗試建立相應(yīng)的數(shù)學(xué)模型。表1對(duì)1000個(gè)會(huì)員調(diào)查的部分結(jié)果DVD名稱DVD1DVD2DVD3DVD4DVD5愿意觀看的人數(shù)200100502510表2現(xiàn)有DVD張數(shù)和當(dāng)前需要處理的會(huì)員的在線訂單DVD編號(hào)D001D002D003D004…DVD現(xiàn)有數(shù)量10401520…會(huì)員在線訂單C00016000…C00020000…C00030003…C00040000…………………注：D001~D100表示100種DVD,C0001~C1000表示1000個(gè)會(huì)員,會(huì)員的在線訂單用數(shù)字1,2,…表示，數(shù)字越小表示會(huì)員的偏愛程度越高，數(shù)字0表示對(duì)應(yīng)的DVD當(dāng)前不在會(huì)員的在線訂單中。合理假設(shè)（1）抽樣調(diào)查的1000人具有代表性，即頻率可以代替概率；（2）DVD租賃業(yè)務(wù)按規(guī)則進(jìn)行，客戶在規(guī)定時(shí)間內(nèi)歸還DVD；（3）每次租賃都是3張DVD；（4）每個(gè)會(huì)員每月的租賃次數(shù)與對(duì)DVD的偏好是獨(dú)立的；（5）各個(gè)會(huì)員之間的選擇和偏好也是相互獨(dú)立的；符號(hào)設(shè)置Xi會(huì)員i每月租賃的次數(shù)，Xi=1,2,i=1,2,…,100000;Yij會(huì)員i對(duì)第j類DVD的選擇，Yij=3,0,j=1,2,3,4,5;i=1,2,…,100000Rj第j類DVD的最少定購張數(shù)，j=1,2,3,4,5;問題（1）的模型根據(jù)條件和假設(shè)（5）會(huì)員i的每月定購次數(shù)X1,X2,…,X100000是服從同分布的隨機(jī)變量，且分布規(guī)律為Xi21P(Xi)0.60.4且X1,…,X100000相互獨(dú)立。在100000個(gè)會(huì)員中抽取1000個(gè)會(huì)員，樣本容量n=1000足夠大，會(huì)員對(duì)DVD的選擇頻率可以直接作為選擇概率的點(diǎn)估計(jì)