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

1、基于基于matlab的計(jì)算機(jī)模擬的計(jì)算機(jī)模擬舒興明TelQ： 117562750一、模擬是什么一、模擬是什么模擬模擬(也稱仿真仿真)：就是用計(jì)算機(jī)程序在計(jì)算機(jī)上模仿各種實(shí)際系統(tǒng)的運(yùn)行過程，并通過計(jì)算了解系統(tǒng)隨時(shí)間變化的行為或特性。計(jì)算機(jī)模擬計(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í)際問題數(shù)學(xué)、邏輯模型數(shù)學(xué)、邏輯模型計(jì)算機(jī)模型計(jì)算機(jī)模型數(shù)學(xué)、計(jì)算機(jī)解數(shù)學(xué)、計(jì)算機(jī)解實(shí)際解實(shí)際解二、為何要模擬二、為何要模擬（1）實(shí)際系統(tǒng)建立之前，要對(duì)系統(tǒng)的行為或結(jié)果進(jìn)行分

2、析研究；）實(shí)際系統(tǒng)建立之前，要對(duì)系統(tǒng)的行為或結(jié)果進(jìn)行分析研究；（2）有些真實(shí)系統(tǒng)做實(shí)驗(yàn)會(huì)影響系統(tǒng)運(yùn)行，例如，在生產(chǎn)中任）有些真實(shí)系統(tǒng)做實(shí)驗(yàn)會(huì)影響系統(tǒng)運(yùn)行，例如，在生產(chǎn)中任意改變工藝系數(shù)可能導(dǎo)致廢品，在經(jīng)濟(jì)活動(dòng)中隨意將一個(gè)決策付意改變工藝系數(shù)可能導(dǎo)致廢品，在經(jīng)濟(jì)活動(dòng)中隨意將一個(gè)決策付諸行動(dòng)會(huì)導(dǎo)致經(jīng)濟(jì)混亂；諸行動(dòng)會(huì)導(dǎo)致經(jīng)濟(jì)混亂；（3）在系統(tǒng)上做多次試驗(yàn)，很難保證每次試驗(yàn)的操作條件相同，）在系統(tǒng)上做多次試驗(yàn)，很難保證每次試驗(yàn)的操作條件相同，因而對(duì)實(shí)驗(yàn)結(jié)果好壞很難作出正確的判斷；因而對(duì)實(shí)驗(yàn)結(jié)果好壞很難作出正確的判斷；（4）當(dāng)人是系統(tǒng)的一部分時(shí)，他的行為往往實(shí)驗(yàn)結(jié)果有所影響，）當(dāng)人是系統(tǒng)的一部分時(shí)，他的

3、行為往往實(shí)驗(yàn)結(jié)果有所影響，這時(shí)，最好進(jìn)行模擬研究；這時(shí)，最好進(jìn)行模擬研究；（5）實(shí)驗(yàn)時(shí)間太長(zhǎng)，費(fèi)用太大，或者有危險(xiǎn)，使得試驗(yàn)不容易）實(shí)驗(yàn)時(shí)間太長(zhǎng)，費(fèi)用太大，或者有危險(xiǎn)，使得試驗(yàn)不容易進(jìn)行；進(jìn)行；（6）有些系統(tǒng)一旦建立起來(lái)后無(wú)法復(fù)原，例如，建立大型企業(yè)，）有些系統(tǒng)一旦建立起來(lái)后無(wú)法復(fù)原，例如，建立大型企業(yè)，要分析社會(huì)和經(jīng)濟(jì)效益，不能用建立起來(lái)試試看的辦法。要分析社會(huì)和經(jīng)濟(jì)效益，不能用建立起來(lái)試試看的辦法。三、哪些問題適合模擬研究三、哪些問題適合模擬研究（1） 難以用數(shù)學(xué)公式表示的系統(tǒng)，或者沒有求解難以用數(shù)學(xué)公式表示的系統(tǒng)，或者沒有求解數(shù)學(xué)模型的有效方法；數(shù)學(xué)模型的有效方法；（2）雖然可以用解析的

4、方法解決問題，但是問題的）雖然可以用解析的方法解決問題，但是問題的分析與計(jì)算過于復(fù)雜，這時(shí)計(jì)算機(jī)仿真可能提供簡(jiǎn)分析與計(jì)算過于復(fù)雜，這時(shí)計(jì)算機(jī)仿真可能提供簡(jiǎn)單可行的求解方法；單可行的求解方法；（3）希望在較短時(shí)間內(nèi)觀察到系統(tǒng)的發(fā)展全過程，）希望在較短時(shí)間內(nèi)觀察到系統(tǒng)的發(fā)展全過程，以估計(jì)某參數(shù)對(duì)系統(tǒng)行為的影響；以估計(jì)某參數(shù)對(duì)系統(tǒng)行為的影響；（4）難以在實(shí)際環(huán)境中進(jìn)行試驗(yàn)和觀察，計(jì)算機(jī)仿）難以在實(shí)際環(huán)境中進(jìn)行試驗(yàn)和觀察，計(jì)算機(jī)仿真是唯一的方法；真是唯一的方法；（5）需要對(duì)系統(tǒng)或過程進(jìn)行長(zhǎng)期運(yùn)行的比較，從大）需要對(duì)系統(tǒng)或過程進(jìn)行長(zhǎng)期運(yùn)行的比較，從大量方案中尋找最優(yōu)。量方案中尋找最優(yōu)。四、模擬的分類四、

5、模擬的分類 模擬是系統(tǒng)狀態(tài)隨時(shí)間而變化的動(dòng)態(tài)寫照，因此，模擬是系統(tǒng)狀態(tài)隨時(shí)間而變化的動(dòng)態(tài)寫照，因此，通常時(shí)間是模擬的主要自變量，其它的變量為因變量。通常時(shí)間是模擬的主要自變量，其它的變量為因變量。（1）按照模擬過程中因變量的變化情況，可以將模擬）按照模擬過程中因變量的變化情況，可以將模擬分為離散、連續(xù)、混合分為離散、連續(xù)、混合3種類型；種類型；（2）如果采用模擬計(jì)算機(jī)、采用數(shù)字計(jì)算機(jī)以及聯(lián)合）如果采用模擬計(jì)算機(jī)、采用數(shù)字計(jì)算機(jī)以及聯(lián)合使用則分為模擬仿真、數(shù)字仿真以及混合仿真；使用則分為模擬仿真、數(shù)字仿真以及混合仿真；（3）根據(jù)仿真變量的特征分為隨機(jī)模擬仿真和模糊模）根據(jù)仿真變量的特征分為隨機(jī)模

6、擬仿真和模糊模擬仿真。擬仿真。五、隨機(jī)模擬的基礎(chǔ)五、隨機(jī)模擬的基礎(chǔ)基礎(chǔ)知識(shí)基礎(chǔ)知識(shí)：線性代數(shù)部分向量矩陣的表述；微積分的：線性代數(shù)部分向量矩陣的表述；微積分的基本積分和微分表述；概率論的基本分布；基本統(tǒng)計(jì)基本積分和微分表述；概率論的基本分布；基本統(tǒng)計(jì)知識(shí)等等。知識(shí)等等。模擬（蒙特卡羅方法）基本原理模擬（蒙特卡羅方法）基本原理：概率論的大數(shù)定律概率論的大數(shù)定律 獨(dú)立隨機(jī)序列的和按概率收斂于這些隨機(jī)序列獨(dú)立隨機(jī)序列的和按概率收斂于這些隨機(jī)序列的數(shù)學(xué)期望的和（常數(shù)），即的數(shù)學(xué)期望的和（常數(shù)），即0, 1| )(|lim1niiniinXEXp六、六、matlab統(tǒng)計(jì)工具箱中產(chǎn)生隨機(jī)數(shù)的函數(shù)統(tǒng)計(jì)工具箱

7、中產(chǎn)生隨機(jī)數(shù)的函數(shù)1、概率密度函數(shù)、概率密度函數(shù)pdf分布類型分布類型調(diào)用格式調(diào)用格式備注備注 二項(xiàng)分布二項(xiàng)分布binoy=binopdf(x,n,p)返回參數(shù)為返回參數(shù)為p，n的二項(xiàng)分的二項(xiàng)分布在布在x處的分布律值處的分布律值指數(shù)分布指數(shù)分布expy=exppdf(x,mu)返回參數(shù)為返回參數(shù)為mu的指數(shù)分布的指數(shù)分布在在x處的密度函數(shù)值處的密度函數(shù)值2分布分布chi2y=chi2pdf(x,v)返回參數(shù)為返回參數(shù)為v的卡方分布的卡方分布在在x處的密度函數(shù)值處的密度函數(shù)值F分布分布fy=fpdf(x,v1,v2)返回參數(shù)為返回參數(shù)為v1,v2的的F分布分布在在x處密度函數(shù)值處密度函數(shù)值幾何分

8、布幾何分布geoy=geopdf(x,p)返回參數(shù)為返回參數(shù)為p的幾何分布的幾何分布在在x處的分布律值處的分布律值超幾何分布超幾何分布hygey=hygepdf(x,M,k,n)返回參數(shù)為返回參數(shù)為M,k,n的超幾的超幾何分在何分在x處的分布律值處的分布律值常見隨機(jī)變量的密度函數(shù)表常見隨機(jī)變量的密度函數(shù)表分布類型分布類型調(diào)用格式調(diào)用格式備注備注 正態(tài)分布正態(tài)分布normnormpdf(x,mu,sig)返回參數(shù)為返回參數(shù)為mu和和sig的正的正態(tài)分布在態(tài)分布在x處的密度值處的密度值泊松分布泊松分布poisspoisspdf(x,lamda)返回參數(shù)為返回參數(shù)為lamda的泊松的泊松分布在分布在

9、x處的分布律值處的分布律值T分布分布ttpdf(x,v)返回參數(shù)為返回參數(shù)為v的的t分布在分布在x處處的密度函數(shù)值的密度函數(shù)值連續(xù)均勻分布連續(xù)均勻分布unifunifpdf(x,a,b)返回在返回在a,b上服從均勻分上服從均勻分布的密度函數(shù)值布的密度函數(shù)值離散均勻分布離散均勻分布unidunidpdf(x,N)返回返回1N的均勻分布概率的均勻分布概率分布律值分布律值密度函數(shù)續(xù)密度函數(shù)續(xù)注意：二項(xiàng)分布注意：二項(xiàng)分布(bino)，幾何分布，幾何分布(geo)，超幾何分布，超幾何分布(hyge),泊松分布泊松分布(poiss)，指數(shù)分布，指數(shù)分布(exp)，正態(tài)分布，正態(tài)分布(norm)，卡方分，卡

10、方分布布(chi2)，t分布分布(t)，F(xiàn)分布分布(f)，均勻連續(xù)分布，均勻連續(xù)分布(unif)，均勻離，均勻離散分布散分布(unid)。2、常見的分布函數(shù)、常見的分布函數(shù)cdf分布類型分布類型調(diào)用格式調(diào)用格式備注備注 二項(xiàng)分布二項(xiàng)分布y=binocdf(x,n,p)返回參數(shù)為返回參數(shù)為p，n的二項(xiàng)分的二項(xiàng)分布在布在x處的分布值處的分布值指數(shù)分布指數(shù)分布y=expcdf(x,mu)返回參數(shù)為返回參數(shù)為mu的指數(shù)分布的指數(shù)分布在在x處的分布函數(shù)值處的分布函數(shù)值2分布分布y=chi2cdf(x,v)返回參數(shù)為返回參數(shù)為v的卡方分布的卡方分布在在x處的分布函數(shù)值處的分布函數(shù)值F分布分布y=fcdf(

11、x,v1,v2)返回參數(shù)為返回參數(shù)為v1,v2的的F分布分布在在x處分布函數(shù)值處分布函數(shù)值幾何分布幾何分布y=geocdf(x,p)返回參數(shù)為返回參數(shù)為p的幾何分布的幾何分布在在x處的分布值處的分布值超幾何分布超幾何分布y=hygecdf(x,M,k,n)返回參數(shù)為返回參數(shù)為M,k,n的超幾的超幾何分在何分在x處的分布值處的分布值常見隨機(jī)變量的分布函數(shù)常見隨機(jī)變量的分布函數(shù) 數(shù)表數(shù)表分布類型分布類型調(diào)用格式調(diào)用格式備注備注 正態(tài)分布正態(tài)分布normcdf(x,mu,sig)返回參數(shù)為返回參數(shù)為mu和和sig的正的正態(tài)分布在態(tài)分布在x處的分布值處的分布值泊松分布泊松分布poisscdf(x,la

12、mda)返回參數(shù)為返回參數(shù)為lamda的泊松的泊松分布在分布在x處的分布值處的分布值T分布分布tcdf(x,v)返回參數(shù)為返回參數(shù)為v的的t分布在分布在x處處的分布函數(shù)值的分布函數(shù)值連續(xù)均勻分布連續(xù)均勻分布unifcdf(x,a,b)返回在返回在a,b上服從均勻分上服從均勻分布的分布函數(shù)值布的分布函數(shù)值離散均勻分布離散均勻分布unidcdf(x,N)返回返回1N的均勻分布概率的均勻分布概率分布值分布值分布函數(shù)續(xù)分布函數(shù)續(xù) 將分布函數(shù)的概率調(diào)用將分布函數(shù)的概率調(diào)用p=cdf(name,x,a1,a2)調(diào)調(diào)過來(lái)用，就是逆分布函數(shù)過來(lái)用，就是逆分布函數(shù)x=icdf(name,p,a1,a2)。 x1

13、=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)生、模擬隨機(jī)數(shù)的產(chǎn)生格式格式1：x=random(name,A1,A2,A3,m,n)產(chǎn)生參數(shù)為產(chǎn)生參數(shù)為A1，A2，A3的的name分布的分布的m行行n列的矩陣。列的矩陣。 x=random(norm,0,1,3,4)x = 0.5377

14、0.8622 -0.4336 2.7694 1.8339 0.3188 0.3426 -1.3499 -2.2588 -1.3077 3.5784 3.0349 y=random(unif,1.2,2.2,3,4)y = 2.1572 1.3419 1.9922 1.2357 1.6854 1.6218 2.1595 2.0491 2.0003 2.1157 1.8557 2.1340產(chǎn)生產(chǎn)生3行行4列的參數(shù)為列的參數(shù)為0和和1的正態(tài)分布隨機(jī)數(shù)的正態(tài)分布隨機(jī)數(shù)產(chǎn)生在產(chǎn)生在1.2,2.2上服上服從均勻分布的隨機(jī)數(shù)從均勻分布的隨機(jī)數(shù)格式格式2：x=namernd(A1,A2,m,n) 產(chǎn)生產(chǎn)生m行

15、，行，n列的服從參數(shù)為列的服從參數(shù)為A1，A2的的name分布的隨機(jī)分布的隨機(jī)數(shù)。數(shù)。 x=trnd(3,4)x = 0.6880 -1.9018 -0.6416 -1.1429 1.0967 -1.6431 3.3100 -2.1621 0.9627 -1.1476 -0.7321 -1.1849 0.0728 1.1055 1.6290 0.6073 y=trnd(3,3,4)y = -4.3874 -1.8874 1.0409 -1.5729 -1.0861 0.8758 -2.5710 5.6567 1.7571 0.1385 -0.2544 1.2256產(chǎn)生自由度為產(chǎn)生自由度為3的的

16、服從服從t分布的隨機(jī)分布的隨機(jī)數(shù)（數(shù)（4階方陣）。階方陣）。產(chǎn)生產(chǎn)生3行行4列自由度列自由度為為3的服從的服從t分布的分布的隨機(jī)數(shù)。隨機(jī)數(shù)。七、隨機(jī)樣本的簡(jiǎn)單描述性統(tǒng)計(jì)量七、隨機(jī)樣本的簡(jiǎn)單描述性統(tǒng)計(jì)量1、中心趨勢(shì)（位置）度量、中心趨勢(shì)（位置）度量樣本的幾何平均樣本的幾何平均m=geomean(X)樣本的調(diào)和平均樣本的調(diào)和平均m=harmmean(X)樣本的算術(shù)平均樣本的算術(shù)平均m=mean(X)樣本的中位數(shù)樣本的中位數(shù)m=median(X)樣本的剔出極端數(shù)據(jù)的平均樣本的剔出極端數(shù)據(jù)的平均 m=trimmean(X,percent) X=random(unif,1,8,10,1); geomea

17、n(X),harmmean(X),mean(X),median(X),trimmean(X,0.2)ans = 4.5903 3.7400 5.2401 6.0678 5.24012、分散度量、分散度量樣本四分位數(shù)間距樣本四分位數(shù)間距m=iqr(X)樣本數(shù)據(jù)的平均絕對(duì)偏差樣本數(shù)據(jù)的平均絕對(duì)偏差m=mad(X) 即即mean(abs(X-mean(X)樣本極差樣本極差m=range(X)樣本方差樣本方差m=var(X)樣本標(biāo)準(zhǔn)差樣本標(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.000

18、0 5.1600 18.0000 35.7333 5.9777八、隨機(jī)模擬計(jì)算積分八、隨機(jī)模擬計(jì)算積分例例1 如下圖所示，在正方形內(nèi)有如下圖所示，在正方形內(nèi)有1/4單位圓。向正方形內(nèi)投小單位圓。向正方形內(nèi)投小石頭，假設(shè)每次都能夠投進(jìn)正方形內(nèi)且等可能落在正方形內(nèi)任石頭，假設(shè)每次都能夠投進(jìn)正方形內(nèi)且等可能落在正方形內(nèi)任何一點(diǎn)。問，小石頭落在何一點(diǎn)。問，小石頭落在1/4單位圓內(nèi)單位圓內(nèi)(包含邊界包含邊界)的概率多大？的概率多大？01x1y分析：假設(shè)投入正方形內(nèi)的分析：假設(shè)投入正方形內(nèi)的石頭有石頭有n塊，有塊，有k塊落入了塊落入了1/4單位圓內(nèi)。單位圓內(nèi)。P為小石頭落為小石頭落入入1/4單位圓內(nèi)的概率

19、。那么單位圓內(nèi)的概率。那么根據(jù)根據(jù)Bernoull(伯努利伯努利)大數(shù)大數(shù)定律，定律，有有1|pnk|plimn 即，當(dāng)實(shí)驗(yàn)次數(shù)即，當(dāng)實(shí)驗(yàn)次數(shù)n充分大時(shí)，頻率和概率之差小于任意數(shù)充分大時(shí)，頻率和概率之差小于任意數(shù)的概率趨于的概率趨于1。而另外一方面由而另外一方面由幾何概型幾何概型有有4Ap樣本空間的測(cè)度的測(cè)度事件在實(shí)際操作中，實(shí)驗(yàn)次數(shù)在實(shí)際操作中，實(shí)驗(yàn)次數(shù)n不可能趨于無(wú)窮大，所以有不可能趨于無(wú)窮大，所以有nk4p44p(數(shù)學(xué)模型數(shù)學(xué)模型)對(duì)于估計(jì)對(duì)于估計(jì)nk4只有不斷重復(fù)做實(shí)驗(yàn)，這種試驗(yàn)可以只有不斷重復(fù)做實(shí)驗(yàn)，這種試驗(yàn)可以具體去操作，具體去操作，(均勻投石塊，然后數(shù)數(shù)，這樣需要較高成本均勻投石

20、塊，然后數(shù)數(shù)，這樣需要較高成本)。也。也可以讓計(jì)算機(jī)去重復(fù)試驗(yàn)，但是需要將數(shù)學(xué)模型轉(zhuǎn)化為計(jì)算可以讓計(jì)算機(jī)去重復(fù)試驗(yàn)，但是需要將數(shù)學(xué)模型轉(zhuǎn)化為計(jì)算機(jī)模擬模型機(jī)模擬模型(讓計(jì)算機(jī)完成均勻投石塊，自動(dòng)計(jì)數(shù)，也需要成本，讓計(jì)算機(jī)完成均勻投石塊，自動(dòng)計(jì)數(shù)，也需要成本，但成本比起前者，幾乎可以忽略。但成本比起前者，幾乎可以忽略。)。用計(jì)算機(jī)模擬投石塊過程和步驟如下：用計(jì)算機(jī)模擬投石塊過程和步驟如下：1、自動(dòng)生成隨機(jī)點(diǎn)、自動(dòng)生成隨機(jī)點(diǎn)0,1x0,1，模擬石塊在正方形內(nèi)的任，模擬石塊在正方形內(nèi)的任意位置，用意位置，用(xi,yi)表示，共表示，共n個(gè)點(diǎn)；個(gè)點(diǎn)；2、判別、判別(xi,yi)是否滿足是否滿足xi2

21、+yi211，即判別石塊是否落在，即判別石塊是否落在1/41/4單位圓內(nèi)，共單位圓內(nèi)，共k k個(gè)點(diǎn)滿足；個(gè)點(diǎn)滿足；3、整理、統(tǒng)計(jì)模擬結(jié)果，用、整理、統(tǒng)計(jì)模擬結(jié)果，用4k/n估計(jì)估計(jì) 。符號(hào)設(shè)置：符號(hào)設(shè)置：符號(hào)說(shuō)明n試驗(yàn)次數(shù)（投石塊次數(shù)）k石塊落入1/4圓內(nèi)次數(shù)(xi,yi)第i塊石頭落點(diǎn)坐標(biāo)pai圓周率計(jì)算機(jī)模擬的流程圖如下計(jì)算機(jī)模擬的流程圖如下初始化：初始化：k=0,i=1in是是Pai=4k/n否否在在0,1x0,1上產(chǎn)生隨機(jī)點(diǎn)上產(chǎn)生隨機(jī)點(diǎn)(xi,yi)r=xi2+yi2r=1是是k=k+1i=i+1否否設(shè)計(jì)一個(gè)設(shè)計(jì)一個(gè)m文件，專門針對(duì)給定文件，專門針對(duì)給定n，模擬計(jì)算，模擬計(jì)算pai的估

22、計(jì)值的估計(jì)值function pai=fpai(n)if n=0 error(n has to be a positive integer );endk=0;i=1;while i=n x=random(unif,0,1,1,2); r=sum(x.2); if r clear tic n=50:10:5000; for k=1:length(n) pai(k)=fpai(n(k); end plot(n,pai,*)TocElapsed time is 223.847533 seconds. mean(pai)ans = 3.1405！值得注意的是，每次模擬即或！值得注意的是，每次模擬即或程

23、序相同，不同人不同計(jì)算時(shí)間，程序相同，不同人不同計(jì)算時(shí)間，結(jié)果一般不同結(jié)果一般不同(因?yàn)殡S機(jī)因?yàn)殡S機(jī))。05001000150020002500300035004000450050002.82.933.13.23.33.43.53.63.7function pai=fpai1(n)if n=0 error(n has to be a positive integer );endk=0;x=random(unif,0,1,n,2);for k1=1:n x1=x(k1,:); r=sum(x1.2).0.5; if r clearticn=50:10:5000;for k=1:length(n)

24、 pai(k)=fpai1(n(k); endplot(n,pai,*)tocElapsed time is 7.372741 seconds.05001000150020002500300035004000450050002.9533.053.13.153.23.253.33.35時(shí)間大大縮短了！時(shí)間大大縮短了！計(jì)算定積分計(jì)算定積分dxx1322分析：若隨機(jī)變量分析：若隨機(jī)變量X的概率分布密度是的概率分布密度是P(x)(axb)xb)，則隨，則隨機(jī)變量機(jī)變量Y=f(X)Y=f(X)的數(shù)學(xué)期望為的數(shù)學(xué)期望為badx)x(p)x(f)X(f E)Y(E當(dāng)當(dāng)X是是a,b上的均勻分布時(shí)，有上的均勻分

25、布時(shí)，有bxa,ab1)x(p將將p(x)代入上面的期望算子，有代入上面的期望算子，有)X( f E)ab(dx)x( fdxab)x( f)X( f Ebaba例例2 而另外一方面，根據(jù)大數(shù)定理，設(shè)而另外一方面，根據(jù)大數(shù)定理，設(shè)y1,y2,yn為來(lái)自總體為來(lái)自總體Y的一組樣本，有的一組樣本，有1| )y(Ey|limn其中其中)X(f E)Y(E)y(E,yn1yn1ii所以有所以有niiniibaxfnabnyabyabXfEabYEabdxxf11)()()()()()()()( 設(shè)設(shè)x1,x2,xn為來(lái)自服從為來(lái)自服從a,b上均勻分布總體上均勻分布總體X的一組樣的一組樣本。且本。且yi

26、=f(xi)，i=1,2,n。dxx1322那么計(jì)算那么計(jì)算的步驟如下：的步驟如下：1、在、在2，3上抽樣上抽樣x1,x2,xn；2、計(jì)算、計(jì)算n,.,2, 1i ,x1y2ii3、整理估計(jì)定積分；、整理估計(jì)定積分；niiyndxx1322231為了對(duì)比模擬結(jié)果，先把這個(gè)積分求出來(lái)：為了對(duì)比模擬結(jié)果，先把這個(gè)積分求出來(lái)： int(1+x2)0.5,2,3)ans =2.6947540047794075420843777595616然后根據(jù)上面的模擬步驟設(shè)計(jì)一個(gè)然后根據(jù)上面的模擬步驟設(shè)計(jì)一個(gè)m文件文件function I=jifen1(n)x=random(unif,2,3,1,n);y=(1+

27、x.2).0.5;I=(3-2)/n*sum(y); clear tic;n=50:10:5000;for k=1:length(n) I(k)=jifen1(n(k); endplot(n,I,*);tocElapsed time is 0.334083 seconds.05001000150020002500300035004000450050002.662.682.72.722.742.762.78例例3求在空間坐標(biāo)系下，兩個(gè)圓柱體相交部分的體積。兩個(gè)求在空間坐標(biāo)系下，兩個(gè)圓柱體相交部分的體積。兩個(gè)1, 12222zxyx圓柱體的柱面方程分別為圓柱體的柱面方程分別為 根據(jù)二重積分的知識(shí)，

28、相交部分根據(jù)二重積分的知識(shí)，相交部分的體積為的體積為122212yxdxdyxVxyz0, 122218yxyxdxdyx=16/30, 122218yxyxdxdyxV計(jì)算計(jì)算的模擬步驟如下的模擬步驟如下（1）在空間）在空間0,1x0,1x0,1上產(chǎn)生均勻隨機(jī)點(diǎn)上產(chǎn)生均勻隨機(jī)點(diǎn)(xi,yi,zi)，i=1,2,n；（2）判斷）判斷1, 12222iiiizxyx是否成立，若有是否成立，若有k個(gè)；個(gè)；（3）用）用8k/n來(lái)估計(jì)來(lái)估計(jì)V的值。的值。把以上步驟編寫成把以上步驟編寫成m文件文件function V=erchjfen(n)x=random(unif,0,1,n,3);k=0;for k

29、1=1:n r1=x(k1,1)2+x(k1,2)2; r2=x(k1,1)2+x(k1,3)2; if (r1=1)&(r2 clear tic n=20:20:4000; for k=1:length(n)v(k)=erchjfen(n(k);end;plot(n,v,*);tocElapsed time is 72.712497 seconds.0500100015002000250030003500400055.25.45.65.866.26.46.6九、評(píng)估事件發(fā)生的概率九、評(píng)估事件發(fā)生的概率 有些隨機(jī)變量，由于不容易找到（或者不存在）準(zhǔn)確的概率有些隨機(jī)變量，由于不容易找到（

30、或者不存在）準(zhǔn)確的概率分布律或者密度函數(shù)，不能算出此類隨機(jī)變量表達(dá)的隨機(jī)的事件分布律或者密度函數(shù)，不能算出此類隨機(jī)變量表達(dá)的隨機(jī)的事件的概率，需要用到模擬。的概率，需要用到模擬。9, 3Prp243221XXXX 其中其中X1服從均勻分布服從均勻分布U2,5，X2服從指數(shù)分布服從指數(shù)分布exp(3)，X3和和X4分別服從正態(tài)分布分別服從正態(tài)分布 N(3,2)和和N(1,1)。例例4求下面事件表達(dá)的概率求下面事件表達(dá)的概率1、在總體、在總體X=(X1,X2,X3,X4)中抽樣，得到樣本中抽樣，得到樣本(x1k,x2k,x3k,x4k)， k=1,2,N；2、判斷隨機(jī)樣本是否滿足不等式組，假設(shè)滿足

31、不等式組的有、判斷隨機(jī)樣本是否滿足不等式組，假設(shè)滿足不等式組的有n個(gè)；個(gè)；3、用、用n/N估計(jì)估計(jì)p。（大數(shù)定律作保障，或者點(diǎn)估計(jì)原理）（大數(shù)定律作保障，或者點(diǎn)估計(jì)原理）步驟：步驟：按照上述步驟編寫按照上述步驟編寫m文件文件function p=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;for k=1:n r1=x1(k)+x2(k)2; r2=x3(k)+x4(k)2; if (r1=3)&(r2 n=20:20:

32、4000; tic for k=1:length(n)p(k)=gailv1(n(k);end;plot(n,p,*);tocElapsed time is 38.580145 seconds. p(length(n)ans =0.8375050010001500200025003000350040000.650.70.750.80.850.90.951十、求隨機(jī)變量的臨界值（上下分為數(shù)）十、求隨機(jī)變量的臨界值（上下分為數(shù)）求使下式成立的最大求使下式成立的最大f值值8 .0fXXXpr33221步驟：步驟：1、產(chǎn)生、產(chǎn)生N個(gè)獨(dú)立隨機(jī)序列個(gè)獨(dú)立隨機(jī)序列(X1k,X2k,X3k),k=1,2,N；

33、2、計(jì)算、計(jì)算fk=f(X1k,X2k,X3k),k=1,2,N;3、從從大到小的序列、從從大到小的序列f1,f2,fN中中 取取f0.8N作為作為f的估計(jì)值，的估計(jì)值，或從小到大的或從小到大的f0.2N作為作為f的估計(jì)值。的估計(jì)值。例例5其中，其中，X1服從均勻分布服從均勻分布U1,3，X2服從指數(shù)分布服從指數(shù)分布exp(1)，X3服從正態(tài)分布服從正態(tài)分布N(2,1)。根據(jù)上面的流程，編寫根據(jù)上面的流程，編寫m文件如下文件如下function f=linjie(n)x1=random(unif,1,3,n,1);x2=random(exp,1,n,1);x3=random(norm,2,1,

34、n,1);for k=1:n y(k)=x1(k)+x2(k)2+x3(k)3;endz=sort(y);k1=round(0.2*n);f=z(k1); clearn=50:50:10000;ticfor k=1:length(n)f(k)=linjie(n(k);end;plot(n,f,*);tocElapsed time is 11.335281 seconds.0100020003000400050006000700080009000100004.24.44.64.855.25.45.65.8十一、報(bào)童問題十一、報(bào)童問題 有個(gè)報(bào)童每天從報(bào)社批發(fā)報(bào)紙，每份報(bào)紙有個(gè)報(bào)童每天從報(bào)社批發(fā)報(bào)紙

35、，每份報(bào)紙1元，然后再以零元，然后再以零售價(jià)售價(jià)1.5元每份售出，如果賣不出，則退回報(bào)社，退回價(jià)格為每元每份售出，如果賣不出，則退回報(bào)社，退回價(jià)格為每份份1.2元。某個(gè)社區(qū)某個(gè)月記錄在案的每天報(bào)紙的需求量如下表元。某個(gè)社區(qū)某個(gè)月記錄在案的每天報(bào)紙的需求量如下表需求量需求量50454038353028天數(shù)天數(shù)2367732頻率頻率0.06670.10.20.23330.23330.10.0667問：報(bào)童每天該批發(fā)多少分報(bào)紙才是最佳策略？問：報(bào)童每天該批發(fā)多少分報(bào)紙才是最佳策略？合理假設(shè)合理假設(shè)（1）報(bào)童兜售報(bào)紙的社區(qū)比較穩(wěn)定，即每個(gè)月的報(bào)紙的需）報(bào)童兜售報(bào)紙的社區(qū)比較穩(wěn)定，即每個(gè)月的報(bào)紙的需求量

36、的分布抽樣差不多；求量的分布抽樣差不多；（2）把上表給的結(jié)果看成社區(qū)每天報(bào)紙的需求分布律。）把上表給的結(jié)果看成社區(qū)每天報(bào)紙的需求分布律。 例例6符號(hào)假設(shè)符號(hào)假設(shè)a訂購(gòu)價(jià)（元訂購(gòu)價(jià)（元/份）份）b零售價(jià)（元零售價(jià)（元/份）份）c退回價(jià)（元退回價(jià)（元/份）份）r每天報(bào)紙的需求量每天報(bào)紙的需求量P(r)每天報(bào)紙的需求量的抽樣分布律每天報(bào)紙的需求量的抽樣分布律x報(bào)童報(bào)紙的批發(fā)量（份）報(bào)童報(bào)紙的批發(fā)量（份）y報(bào)童每天的獲利（元）報(bào)童每天的獲利（元）建立模型建立模型 每天報(bào)紙的需求量每天報(bào)紙的需求量r是個(gè)隨機(jī)變量，報(bào)童每天的盈利由是個(gè)隨機(jī)變量，報(bào)童每天的盈利由r和和批發(fā)量批發(fā)量x決定，也是隨機(jī)變量，即決定

37、，也是隨機(jī)變量，即rxrxacrabrxxaby),)()(,)(目標(biāo)目標(biāo)1 每天盈利期往最大每天盈利期往最大 max E(y)數(shù)學(xué)模型數(shù)學(xué)模型1 （數(shù)學(xué)期望模型）（數(shù)學(xué)期望模型）rxrxacrabrxxabyxtsyE),)()(,)(0.)(maxr P(r)模型模型1求解求解 將報(bào)童每將報(bào)童每天盈利函數(shù)編天盈利函數(shù)編寫成寫成m文件文件function Ey=baotong(x)global a b c Prif ba|ca error(b and c must be larger than a);endn=size(Pr,2);r=Pr(1,:);p=Pr(2,:);for k=1:n

38、if x clear global a b c Pr a=1;b=1.5;c=1.2; Pr=50 45 40 38 35 30 28;0.0667 0.1 0.2 0.2333 0.2333 0.1 0.0667; x=1:1:60; for k=1:length(x)Ey(k)=baotong(x(k);end plot(x,Ey,*)0102030405060024681012141618 Ey1=max(Ey); m=find(Ey=Ey1)m = 40最佳批發(fā)為最佳批發(fā)為40份報(bào)紙。份報(bào)紙。數(shù)學(xué)模型數(shù)學(xué)模型2 （最低盈利機(jī)會(huì)最大模型）（最低盈利機(jī)會(huì)最大模型）目標(biāo)目標(biāo)2rxrxacra

39、brxxabyxtsyp),)()(,)(0.16maxr P(r)每天至少賺每天至少賺16元的概率最大元的概率最大（機(jī)會(huì)最大）（機(jī)會(huì)最大）將每天賺將每天賺16元的機(jī)會(huì)計(jì)算程序編制為元的機(jī)會(huì)計(jì)算程序編制為m文件：文件：function p=baotongjh(x,ybar)global a b c Prif ba|ca error(b and c must be larger than a);endn=size(Pr,2);r=Pr(1,:);pr=Pr(2,:);for k=1:n if xybar;p=sum(py.*pr); clear global a b c Pr a=1;b=1.5

40、;c=1.2; Pr=50 45 40 38 35 30 28;0.0667 0.1 0.2 0.2333 0.2333 0.1 0.0667; ybar=15; x=1:60; for k=1:60p(k)=baotongjh(x(k),ybar);end plot(x,p,*)010203040506000.10.20.30.40.50.60.70.80.9 由此可見，當(dāng)報(bào)紙批發(fā)由此可見，當(dāng)報(bào)紙批發(fā)量在量在3045之間時(shí)，報(bào)童每之間時(shí)，報(bào)童每天可以賺天可以賺15元錢的概率最元錢的概率最大有大有83%。十二、利用模擬法求解非線性規(guī)劃十二、利用模擬法求解非線性規(guī)劃例例7用隨機(jī)模擬方法求解給定的

41、規(guī)劃問題用隨機(jī)模擬方法求解給定的規(guī)劃問題0,22103. .382max21212121212221xxxxxxtsxxxxxxz求解此問題的流程圖如下求解此問題的流程圖如下初始化初始化n=0,A= n1000YN在在0,10上產(chǎn)生均勻隨機(jī)數(shù)上產(chǎn)生均勻隨機(jī)數(shù)x1，x2=(10-x1)/3，y=x1+2x2y2NYn=n+1,A=A;x1,x2z=0,k=1:1000,xx=0,0 x1=A(k,1),x2=A(k,2)z1=f(x1,x2)z1zYz=z1Nk=k+1xx=x1,x2依照上述流程編寫一個(gè)求此非線性規(guī)劃最大值的依照上述流程編寫一個(gè)求此非線性規(guī)劃最大值的m文件：文件：functio

42、n z,xx=nlpf(N)n=0;A=;while n2 n=n+1; A=A;x1,x2; endendz=0;xx=0,0;for k=1:N x11=A(k,1);x22=A(k,2); z1=-2*x112-x222+x11*x22+8*x11+3*x22; if z1z z=z1; xx=x11,x22; endend clear A=; N=20:5:400; for k=1:length(N)z,xx=nlpf(N(k);A=A;z,xx;end plot(N,A(:,1),*)05010015020025030035040014.6514.714.7514.814.8514.

43、914.951515.0515.1十三、一個(gè)排隊(duì)系統(tǒng)的模擬十三、一個(gè)排隊(duì)系統(tǒng)的模擬例例7一個(gè)排隊(duì)系統(tǒng)的仿真一個(gè)排隊(duì)系統(tǒng)的仿真 某店只有一個(gè)收款臺(tái)，顧客到達(dá)收款臺(tái)的時(shí)間間隔服從均某店只有一個(gè)收款臺(tái)，顧客到達(dá)收款臺(tái)的時(shí)間間隔服從均值為值為4.54.5的負(fù)指數(shù)分布，每個(gè)顧客的服務(wù)時(shí)間服從均值為的負(fù)指數(shù)分布，每個(gè)顧客的服務(wù)時(shí)間服從均值為3.23.2，標(biāo)準(zhǔn)差為標(biāo)準(zhǔn)差為0.60.6的正態(tài)分布，對(duì)的正態(tài)分布，對(duì)100100位顧客去收款臺(tái)繳款的排隊(duì)過位顧客去收款臺(tái)繳款的排隊(duì)過程進(jìn)行仿真。程進(jìn)行仿真。分析：分析： 排隊(duì)系統(tǒng)涉及的參考變量為：顧客的等待時(shí)間排隊(duì)系統(tǒng)涉及的參考變量為：顧客的等待時(shí)間WqWq，顧客的，

44、顧客的逗留時(shí)間逗留時(shí)間W,W,排隊(duì)長(zhǎng)排隊(duì)長(zhǎng)Lq,Lq,隊(duì)長(zhǎng)隊(duì)長(zhǎng)L L，服務(wù)臺(tái)的空閑時(shí)間，服務(wù)臺(tái)的空閑時(shí)間I I，服務(wù)臺(tái)的，服務(wù)臺(tái)的繁忙時(shí)間繁忙時(shí)間B B，排隊(duì)系統(tǒng)的狀態(tài)概率，排隊(duì)系統(tǒng)的狀態(tài)概率P P。最后需要統(tǒng)計(jì)出有關(guān)指標(biāo)。最后需要統(tǒng)計(jì)出有關(guān)指標(biāo)的平均值均方差。的平均值均方差。合理假設(shè)合理假設(shè)1 1、進(jìn)入系統(tǒng)的顧客沒有人因?yàn)椴辉敢獾榷x開系統(tǒng)；、進(jìn)入系統(tǒng)的顧客沒有人因?yàn)椴辉敢獾榷x開系統(tǒng)；2 2、先到先服務(wù)；、先到先服務(wù)；模擬步驟：模擬步驟：1 1、產(chǎn)生初始事件表，并從初始事件表中找出最靠前的、產(chǎn)生初始事件表，并從初始事件表中找出最靠前的 事件；事件；2 2、啟動(dòng)仿真時(shí)鐘；、啟動(dòng)仿真時(shí)鐘；3

45、3、判別當(dāng)前事件是哪一個(gè)事件，如果是顧客到達(dá)事件，則啟動(dòng)、判別當(dāng)前事件是哪一個(gè)事件，如果是顧客到達(dá)事件，則啟動(dòng)顧客到達(dá)子程序；如果是服務(wù)事件，則啟動(dòng)服務(wù)結(jié)束子程序；顧客到達(dá)子程序；如果是服務(wù)事件，則啟動(dòng)服務(wù)結(jié)束子程序；4 4、重復(fù)、重復(fù)2-32-3，直到仿真完畢；，直到仿真完畢；5 5、整理、分析并輸出結(jié)果。、整理、分析并輸出結(jié)果。顧客到達(dá)的子程序顧客到達(dá)的子程序1 1、產(chǎn)生下一個(gè)顧客到來(lái)的時(shí)刻，記入事件表；、產(chǎn)生下一個(gè)顧客到來(lái)的時(shí)刻，記入事件表；2 2、判別服務(wù)臺(tái)是否忙、判別服務(wù)臺(tái)是否忙s=1?s=1?，若是，則顧客則排隊(duì)長(zhǎng)，若是，則顧客則排隊(duì)長(zhǎng)L=L+1;L=L+1;否否則，服務(wù)臺(tái)空閑則，

46、服務(wù)臺(tái)空閑,L=0;,L=0;3 3、產(chǎn)生一個(gè)服務(wù)結(jié)束時(shí)刻，記入事件表；、產(chǎn)生一個(gè)服務(wù)結(jié)束時(shí)刻，記入事件表；4 4、統(tǒng)計(jì)所需的數(shù)據(jù)；、統(tǒng)計(jì)所需的數(shù)據(jù)；服務(wù)結(jié)束子程序服務(wù)結(jié)束子程序1 1、結(jié)束服務(wù)，產(chǎn)生結(jié)束時(shí)間，記入事件表；、結(jié)束服務(wù)，產(chǎn)生結(jié)束時(shí)間，記入事件表；2 2、判別是否有顧客排隊(duì)，、判別是否有顧客排隊(duì)，L=0?L=0?，若，若L=0L=0，則服務(wù)臺(tái)空閑；否則，則服務(wù)臺(tái)空閑；否則，排隊(duì)長(zhǎng)排隊(duì)長(zhǎng)L=L-1;L=L-1;3 3、產(chǎn)生服務(wù)結(jié)束時(shí)刻，記入事件表；、產(chǎn)生服務(wù)結(jié)束時(shí)刻，記入事件表；4 4、統(tǒng)計(jì)所需的數(shù)據(jù)；、統(tǒng)計(jì)所需的數(shù)據(jù)；產(chǎn)生初始事件表產(chǎn)生初始事件表從事件表中找出最靠前的事件從事件表中

47、找出最靠前的事件仿真時(shí)鐘步進(jìn)仿真時(shí)鐘步進(jìn)是哪一類事件是哪一類事件顧客顧客到達(dá)到達(dá)服務(wù)服務(wù)結(jié)束結(jié)束產(chǎn)生下一個(gè)顧客到達(dá)紀(jì)錄，記入事件表產(chǎn)生下一個(gè)顧客到達(dá)紀(jì)錄，記入事件表結(jié)束服務(wù)，產(chǎn)生結(jié)束時(shí)間，記入事件表結(jié)束服務(wù)，產(chǎn)生結(jié)束時(shí)間，記入事件表出納是否忙？出納是否忙？出納空閑，出納空閑，s=0s=0顧客排隊(duì)，顧客排隊(duì)，Lq=Lq+1Lq=Lq+1否否是是產(chǎn)生服務(wù)結(jié)束時(shí)刻，記入事件表產(chǎn)生服務(wù)結(jié)束時(shí)刻，記入事件表是否顧客排隊(duì)是否顧客排隊(duì)隊(duì)長(zhǎng)減隊(duì)長(zhǎng)減1 1，Lq-1Lq-1出納空閑，出納空閑，s=0 s=0 否否產(chǎn)生服務(wù)結(jié)束時(shí)刻，記入事件表產(chǎn)生服務(wù)結(jié)束時(shí)刻，記入事件表是是統(tǒng)計(jì)所需數(shù)據(jù)統(tǒng)計(jì)所需數(shù)據(jù)仿真是否完畢仿真

48、是否完畢是是結(jié)束程序結(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; while nlikain; if t1t2 t=t2; tfuwu=tfuwu,t; if Lq0 s=1; L=L-1; Lq=Lq-1; r2=normrnd(b,c); t2=t+r2; else s=0; L

49、=0; Ls=max(L-1,0); t2=+inf; end LL=LL,L; LLq=LLq,Lq; BI=BI,s; nlikai=nlikai+1; ndaoda=ndaoda; elseif t1=t2 t=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

50、,tfuwu=paiduisimu(n,a,b,c); mean(LLq),mean(LL),mean(BI)ans = 0.9602 1.7662 0.8060 由此可見，這個(gè)服務(wù)窗口一直處以繁忙狀態(tài)（空閑率由此可見，這個(gè)服務(wù)窗口一直處以繁忙狀態(tài)（空閑率)的的概率為概率為1-0.8060;平均排隊(duì)長(zhǎng)（平均排隊(duì)長(zhǎng)（mean(L)=0.9602人，不用再增人，不用再增加一個(gè)窗口來(lái)緩解壓力。加一個(gè)窗口來(lái)緩解壓力。十四、航班延誤問題十四、航班延誤問題 某航空公司經(jīng)營(yíng)某航空公司經(jīng)營(yíng)A A，B B，C C三個(gè)城市的航線，這些航線每天班三個(gè)城市的航線，這些航線每天班次起飛與到達(dá)時(shí)間如下表所示。次起飛與到達(dá)

51、時(shí)間如下表所示。 設(shè)飛機(jī)在機(jī)場(chǎng)停留的損失費(fèi)大致與停留時(shí)間成正比，又每設(shè)飛機(jī)在機(jī)場(chǎng)停留的損失費(fèi)大致與停留時(shí)間成正比，又每架飛機(jī)從降落到下班起飛至少需架飛機(jī)從降落到下班起飛至少需2 2小時(shí)準(zhǔn)備時(shí)間，試決定一個(gè)使小時(shí)準(zhǔn)備時(shí)間，試決定一個(gè)使停留費(fèi)用損失為最小的分派飛行方案。停留費(fèi)用損失為最小的分派飛行方案。 另外，根據(jù)抽樣資料，三個(gè)城市另外，根據(jù)抽樣資料，三個(gè)城市A A，B B，C C之間的航班延誤時(shí)之間的航班延誤時(shí)間分別服從負(fù)指數(shù)分布：間分別服從負(fù)指數(shù)分布：城市城市A到城市到城市B的航班延誤時(shí)間服從均值為的航班延誤時(shí)間服從均值為1的指數(shù)分布；的指數(shù)分布；城市城市A到城市到城市C的航班延誤時(shí)間服從均值

52、為的航班延誤時(shí)間服從均值為1.5的指數(shù)分布；的指數(shù)分布；城市城市B到城市到城市C的航班延誤時(shí)間服從均值為的航班延誤時(shí)間服從均值為0.8的指數(shù)分布；的指數(shù)分布；城市城市B到城市到城市A的航班延誤時(shí)間服從均值為的航班延誤時(shí)間服從均值為1.2的指數(shù)分布；的指數(shù)分布；城市城市C到城市到城市A的航班延誤時(shí)間服從均值為的航班延誤時(shí)間服從均值為0.8的指數(shù)分布；的指數(shù)分布；城市城市C到城市到城市B的航班延誤時(shí)間服從均值為的航班延誤時(shí)間服從均值為0.5的指數(shù)分布；的指數(shù)分布；航班號(hào)航班號(hào) 起飛城市起飛城市 起飛時(shí)間起飛時(shí)間 到達(dá)城市到達(dá)城市 到達(dá)時(shí)間到達(dá)時(shí)間101 A 9：00 B 12：00102 A 10

53、：00 B 13：00103 A 22：00 C 2：00（次日（次日)104 B 4：00 A 7：00105 B 11：00 A 14：00106 C 7：00 A 11：00107 B 13：00 C 18：00108 C 15：00 B 20：00各個(gè)航班往返城市及其時(shí)間如下表：各個(gè)航班往返城市及其時(shí)間如下表：合理假設(shè)合理假設(shè)（1）城市）城市A，B，C的到達(dá)與起飛航班數(shù)目一樣，每個(gè)城市都的到達(dá)與起飛航班數(shù)目一樣，每個(gè)城市都用到達(dá)航班的飛機(jī)作為起飛航班的飛機(jī)用到達(dá)航班的飛機(jī)作為起飛航班的飛機(jī);（2）沒有多余的飛機(jī)可用；）沒有多余的飛機(jī)可用；（3）只算一天內(nèi)的延誤時(shí)間；）只算一天內(nèi)的延誤時(shí)

54、間；（3）飛機(jī)到達(dá)后）飛機(jī)到達(dá)后2個(gè)小時(shí)任何時(shí)刻都可以起飛；個(gè)小時(shí)任何時(shí)刻都可以起飛；（4）先到先服務(wù)原則；）先到先服務(wù)原則；（5）概率小于）概率小于0.0001為小概率事件，不考慮。為小概率事件，不考慮。符號(hào)設(shè)置符號(hào)設(shè)置xijxij=1,城市城市A的第的第i到達(dá)作為第到達(dá)作為第j離開，否則離開，否則xij=0yijyij=1,城市城市B的第的第i到達(dá)作為第到達(dá)作為第j離開，否則離開，否則yij=0zijzij=1,城市城市C的第的第i到達(dá)作為第到達(dá)作為第j離開，否則離開，否則zij=0t1ijA城市第城市第i到達(dá)作為第到達(dá)作為第j離開的等待時(shí)間離開的等待時(shí)間t2ijB城市第城市第i到達(dá)作為第

55、到達(dá)作為第j離開的等待時(shí)間離開的等待時(shí)間t3ijC城市第城市第i到達(dá)作為第到達(dá)作為第j離開的等待時(shí)間離開的等待時(shí)間T TEFEF第第i城市到第城市到第j城市的延誤時(shí)間城市的延誤時(shí)間E=A,B,C,F=A,B,C,EF 同時(shí)根據(jù)假設(shè)（同時(shí)根據(jù)假設(shè)（5），所有航班延誤超過某個(gè)時(shí)刻后是），所有航班延誤超過某個(gè)時(shí)刻后是小概率事件，認(rèn)為不發(fā)生。小概率事件，認(rèn)為不發(fā)生。1、可行解的產(chǎn)生方法、可行解的產(chǎn)生方法k=nNstop初始化：初始化：y=,k=1Y產(chǎn)生隨機(jī)數(shù)產(chǎn)生隨機(jī)數(shù)a1U(0,1)a12/3Yx= 0 0 1Nx= 0 1 0產(chǎn)生隨機(jī)數(shù)產(chǎn)生隨機(jī)數(shù)a2U(0,1)a2=0.5確定確定x的的0所在所在列

56、第二行列第二行1的位置的位置確定剩確定剩1行行1列列1的位置的位置y=y;xk=k+1function y=hbzhipai3(n)y=;k=1;while k=n a1=random(unif,0,1); if a1=1/3 x=1,0,0; a2=random(unif,0,1); if a22/3 x=0 0 1; a2=random(unif,0,1); if a2=0.5 x=x;1 0 0;x=x;0 1 0; elsex=x;0 1 0;x=x;1 0 0; end else x=0 1 0; a2=random(unif,0,1); if a2=0.5 x=x;1 0 0; x

57、=x;0 0 1; else x=x;0 0 1;x=x;1 0 0; end end y=y;x;k=k+1;end 但是，根據(jù)這個(gè)問題的特殊性，城市但是，根據(jù)這個(gè)問題的特殊性，城市A的航班安排只有的航班安排只有廈門廈門6個(gè)可行解：個(gè)可行解：x1=1 0 0;0 1 0;0 0 1;x2=1 0 0;0 0 1;0 1 0;x3=0 1 0;0 0 1;1 0 0;x4=0 1 0;1 0 0;0 0 1;x5=0 0 1;0 1 0;1 0 0;x6=0 0 1;1 0 0;0 1 0;2、模擬在給定延誤情況下的等待時(shí)間、模擬在給定延誤情況下的等待時(shí)間航班號(hào)航班號(hào)101(9)102(10)

58、103(22)B：104(7)2315B：105(14)19208C：106(11)222311城市A的航班間的不延誤等待時(shí)間t1ij2.1 航班延誤時(shí)航班延誤時(shí)A城市的各個(gè)等待時(shí)間表城市的各個(gè)等待時(shí)間表由于由于B和和C到到A城市的航班延誤時(shí)間服從的分布分別為城市的航班延誤時(shí)間服從的分布分別為)8.0(),2.1(eTeTCABA當(dāng)當(dāng)TBA和和TCA取定時(shí)，航班飛機(jī)的等待時(shí)間就跟著確定，如下取定時(shí)，航班飛機(jī)的等待時(shí)間就跟著確定，如下流程流程初始化：初始化：tba1e(1.2)tba11,tba1=13Yt1(104,103)=15-tba1t1(104,101)=24-tba1t1(104,1

59、02)=25-tba1Nt1(104,101)=26-tba1t1(104,102)=17-tba1t1(104,103)=39-tba1初始化：初始化：tba2e(1.2)tba26,tba217,tba2=18YNt1(105,101)=43-tba2t1(105,102)=44-tba2t1(105,103)=32-tba2初始化：初始化：tba3e(0.8)tba39,tba320,tba317)&(tba2=18) t12=43 20 32-tba2; else t12=43 44 32-tba2; end tba3=random(exp,0.8); if tba39)&am

60、p;(tba320)&(tba3=21) t13=46 23 35-tba3; else t13=46 47 35-tba3; end tr1=t11;t12;t13; t1=t1;tr1; k=k+1;endfunction t1=hbwaita(m)t1=;k=1;while k=m tba1=random(exp,1.2); if tba11)&(tba1=13) t11=24 25 15-tba1; else t11=26 27 39-tba1; end tba2=random(exp,1.2); if tba26)&(tba2 clearm=10000;x1=1 0 0;0 1 0;0 0 1;x2=1 0 0