蒙特卡羅方法教材課件

數(shù)學(xué)建模專題之

蒙特卡羅方法主講：李培巒時間：2014-08-032022/12/221數(shù)學(xué)建模專題之

蒙特卡羅方法主講：李培巒2022/12/19內(nèi)容提綱1.引言2.MonteCarlo模擬基本思想3.隨機(jī)數(shù)生成函數(shù)4.應(yīng)用實例舉例5.隨機(jī)游走模擬6.機(jī)械零件的可靠度計算7.排隊論模擬8.MonteCarlo方法求解規(guī)劃問題9.MonteCarlo方法預(yù)測搜救區(qū)域2022/12/222內(nèi)容提綱1.引言2022/12/192引言(Introduction)MonteCarlo方法：蒙特卡羅方法，又稱隨機(jī)模擬方法，屬于計算數(shù)學(xué)的一個分支，它是在上世紀(jì)四十年代中期為了適應(yīng)當(dāng)時原子能事業(yè)的發(fā)展而發(fā)展起來的。亦稱統(tǒng)計模擬方法，statisticalsimulationmethod利用隨機(jī)數(shù)進(jìn)行數(shù)值模擬的方法MonteCarlo名字的由來：MonteCarlo是摩納哥（monaco)最大的城市，該城以賭博聞名。Monte-Carlo,Monaco2022/12/223引言(Introduction)MonteCarlo方法：引言(Introduction)MonteCarlo方法的起源：二十世紀(jì)四十年代中期，由于科學(xué)技術(shù)的發(fā)展和電子計算機(jī)的發(fā)明，蒙特卡羅方法作為一種獨立的方法被提出來，并首先在核武器的試驗與研制中得到了應(yīng)用(中子的鏈鎖反應(yīng))。但其基本思想并非新穎，可追溯到18世紀(jì)后半葉的蒲豐（Buffon）隨機(jī)投針實驗（1777年）。NicholasMetropolis(1915-1999)JohnVonNeumann(1903-1957)

2022/12/224引言(Introduction)MonteCarlo方法的引言(Introduction)MonteCarlo方法的應(yīng)用：物理：核物理，熱力學(xué)與統(tǒng)計物理，粒子輸運問題等數(shù)學(xué)：多重積分、解微分方程、非線性方程組求解等工程領(lǐng)域：真空技術(shù)，水力學(xué)，激光技術(shù)等經(jīng)濟(jì)學(xué)領(lǐng)域：期權(quán)定價、項目管理、投資風(fēng)險決策等其他領(lǐng)域：化學(xué)、醫(yī)學(xué)，生物，生產(chǎn)管理、系統(tǒng)科學(xué)、公用事業(yè)等方面。隨著科學(xué)技術(shù)的發(fā)展，其應(yīng)用范圍將更加廣泛。2022/12/225引言(Introduction)MonteCarlo方法的MonteCarlo方法的基本思想Buffon投針實驗1768年，法國數(shù)學(xué)家ComtedeBuffon利用投針實驗估計的值dL2022/12/226MonteCarlo方法的基本思想Buffon投針實驗17SolutionThepositioningoftheneedlerelativetonearbylinescanbedescribedwitharandomvectorwhichhascomponents:Therandomvectorisuniformlydistributedontheregion[0,d)×[0,).Accordingly,ithasprobabilitydensityfunction1/d.Theprobabilitythattheneedlewillcrossoneofthelinesisgivenbytheintegral2022/12/227SolutionThepositioningofthe例1.蒲豐投針問題利用關(guān)系式：求出π值其中Ｎ為投計次數(shù)，n為針與平行線相交次數(shù)。這就是古典概率論中著名的蒲豐氏問題。2022/12/228例1.蒲豐投針問題利用關(guān)系式：2022/12/198一些人進(jìn)行了實驗，其結(jié)果列于下表：實驗者年份投計次數(shù)π的實驗值沃爾弗(Wolf)185050003.1596斯密思(Smith)185532043.1553福克斯(Fox)189411203.1419拉查里尼(Lazzarini)19013408312/229實驗者年份投計次數(shù)π的實驗值沃爾弗(Wolf)1850500基本思想

MS基本思想：為了求解數(shù)學(xué)、物理、工程技術(shù)或隨機(jī)服務(wù)系統(tǒng)等方面的問題，首先構(gòu)造一個模型（概率統(tǒng)計模型），使所求問題的解正好是該模型的參數(shù)或特征量或有關(guān)量，然后通過模擬（統(tǒng)計試驗），給出模型參數(shù)或特征量的估計值，最后得出所求問題的近似解。

MS特點：1.方法新穎、應(yīng)用面廣、實用性強(qiáng)；

2.隨機(jī)模擬方面的算法簡單，但計算量大；3.模擬結(jié)果具有隨機(jī)性，精度較低；4.模擬結(jié)果的收斂過程服從概率規(guī)律。2022/12/2210基本思想MS基本思想：為了求解數(shù)學(xué)、物理、工程技術(shù)或隨機(jī)服例1

在我方某前沿防守地域，敵人以一個炮排（含兩門火炮）為單位對我方進(jìn)行干擾和破壞．為躲避我方打擊，敵方對其陣地進(jìn)行了偽裝并經(jīng)常變換射擊地點．經(jīng)過長期觀察發(fā)現(xiàn)，我方指揮所對敵方目標(biāo)的指示有50％是準(zhǔn)確的，而我方火力單位，在指示正確時，有1/3的射擊效果能毀傷敵人一門火炮，有1/6的射擊效果能全部毀傷敵人火炮．現(xiàn)在希望能用某種方式把我方將要對敵人實施的20次打擊結(jié)果顯現(xiàn)出來，確定有效射擊的比率及毀傷敵方火炮的平均值。分析：這是一個概率問題，可以通過理論計算得到相應(yīng)的概率和期望值.但這樣只能給出作戰(zhàn)行動的最終靜態(tài)結(jié)果，而顯示不出作戰(zhàn)行動的動態(tài)過程.為了能顯示我方20次射擊的過程，現(xiàn)采用模擬的方式。舉例2022/12/2211例1在我方某前沿防守地域，敵人以一個炮排（含兩門火炮）為單需要模擬出以下兩件事：

[2]當(dāng)指示正確時，我方火力單位的射擊結(jié)果情況[1]觀察所對目標(biāo)的指示正確與否模擬試驗有兩種結(jié)果，每一種結(jié)果出現(xiàn)的概率都是1/2．因此，可用投擲一枚硬幣的方式予以確定，當(dāng)硬幣出現(xiàn)正面時為指示正確，反之為不正確．模擬試驗有三種結(jié)果：毀傷一門火炮的可能性為1/3(即2/6)，毀傷兩門的可能性為1/6，沒能毀傷敵火炮的可能性為1/2(即3/6)．這時可用投擲骰子的方法來確定：如果出現(xiàn)的是１、２、３三個點：則認(rèn)為沒能擊中敵人；如果出現(xiàn)的是４、５點：則認(rèn)為毀傷敵人一門火炮；若出現(xiàn)的是６點：則認(rèn)為毀傷敵人兩門火炮．問題分析2022/12/2212需要模擬出以下兩件事：[2]i：要模擬的打擊次數(shù)；k1：沒擊中敵人火炮的射擊總數(shù)；k2：擊中敵人一門火炮的射擊總數(shù)；k3：擊中敵人兩門火炮的射擊總數(shù)；E：有效射擊比率；E1：20次射擊平均每次毀傷敵人的火炮數(shù)．符號說明2022/12/2213i：要模擬的打擊次數(shù)；符號說明2022模擬框圖初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子點數(shù)?k1=k1+1k2=k2+1k3=k3+1k1=k1+1i＜20?E=(k2+k3)/20E1=(0*k1+1*k2+2*k3)/20停止硬幣正面?YNNY1，2，34，562022/12/2214模擬框圖初始化:i=0,k1=0,k2=0,k3=0i=i+模擬結(jié)果2022/12/2215模擬結(jié)果2022/12/1915從以上模擬結(jié)果可計算出：2022/12/2216從以上模擬結(jié)果可計算出：2022/12/1916理論計算2022/12/2217理論計算2022/12/1917結(jié)果比較

雖然模擬結(jié)果與理論計算不完全一致，但它卻能更加真實地表達(dá)實際戰(zhàn)斗動態(tài)過程．要使結(jié)果接近理論計算值，必須加大模擬次數(shù)，這就要求使用計算機(jī)模擬了。

用蒙特卡洛方法進(jìn)行計算機(jī)模擬的步驟：[1]設(shè)計一個邏輯框圖，即建立模擬模型．[2]根據(jù)流程圖編寫程序，模擬隨機(jī)現(xiàn)象．可通過生成具有各種概率分布的隨機(jī)數(shù)來模擬隨機(jī)現(xiàn)象，進(jìn)行模擬試驗．[3]分析模擬結(jié)果，給出所求問題的近似解（求解）.2022/12/2218結(jié)果比較雖然模擬結(jié)果與理論計算不完全一致，但它卻能更注：rand(n)=rand(n,n)Matlab中的隨機(jī)數(shù)生成函數(shù)randperm(m)生成一個由1:m

組成的隨機(jī)排列randn(m,n)生成一個滿足正態(tài)m

n

的隨機(jī)矩陣rand(m,n)

生成一個滿足均勻分布的m

n

隨機(jī)矩陣，矩陣的每個元素都在(0,1)

之間。perms(1:n)生成由1:n

組成的全排列，共n!

個結(jié)果2022/12/2219注：rand(n)=rand(n,n)Matlab中的隨機(jī)name

的取值可以是'norm'or'Normal''unif'or'Uniform''poiss'or'Poisson''beta'or'Beta''exp'or'Exponential''gam'or'Gamma''geo'or'Geometric''unid'or'DiscreteUniform'......random('name',A1,A2,A3,M,N)Matlab中的隨機(jī)數(shù)生成函數(shù)2022/12/2220name的取值可以是'norm'or'Normalfix(x):

截尾取整，直接將小數(shù)部分舍去floor(x):

不超過x

的最大整數(shù)ceil(x):

不小于x

的最小整數(shù)round(x):

四舍五入取整Matlab中的取整函數(shù)2022/12/2221fix(x):截尾取整，直接將小數(shù)部分舍去Matla模擬隨機(jī)投擲均勻硬幣，驗證國徽朝上與朝下的概率是否都是1/2

n=10000;%給定試驗次數(shù)m=0;%m表示試驗成功（國徽朝上）的次數(shù)fori=1:nx=randperm(2)-1;%randperm(2）生成1和2的隨機(jī)排列y=x(1);ify==0%0表示國徽朝上，1表示國徽朝下m=m+1;endendfprintf('國徽朝上的頻率為：%f\n',m/n);小實例一：投擲硬幣2022/12/2222模擬隨機(jī)投擲均勻硬幣，驗證國徽朝上與朝下的概率是否都是1隨機(jī)投擲骰子，驗證各點出現(xiàn)的概率是否為1/6

n=10000;m1=0;m2=0;m3=0;m4=0;m5=0;m6=0;fori=1:nx=randperm(6);y=x(1);switchycase1,m1=m1+1;case2,m2=m2+1;case3,m3=m3+1;case4,m4=m4+1;case5,m5=m5+1;otherwise,m6=m6+1;endend...%

輸出結(jié)果小實例二：投擲骰子2022/12/2223隨機(jī)投擲骰子，驗證各點出現(xiàn)的概率是否為1/6n=100

用蒙特卡羅(MonteCarlo)投點法計算

的值n=10000;a=2;m=0;%m表示落入圓內(nèi)的次數(shù)fori=1:nx=rand(1)*a/2;y=rand(1)*a/2;if(x^2+y^2<=(a/2)^2)m=m+1;endendfprintf('計算出來的pi為：%f\n',4*m/n);小實例三：蒙特卡羅投點法2022/12/2224用蒙特卡羅(MonteCarlo)投點法計算小實例三：蒙特卡羅投點法ezplot('x^2+y^2-1',[-1.11.1]);holdonaxisequalplot([-1-111-1],[-111-1-1]);N=0;fork=1:100N_point=10000;xy=(rand(2,N_point)-0.5)*2;d=sqrt(xy(1,:).^2+xy(2,:).^2);N(k)=length(find(d<=1));endp1=sum(N)/(N_point*k)p2=pi*1^2/4pi0=p1*4p_xy=(rand(500,2)-0.5).*2;plot(p_xy(:,1),p_xy(:,2),'.');2022/12/2225小實例三：蒙特卡羅投點法ezplot('x^2+y^2-1'

在畫有許多間距為d

的等距平行線的白紙上，隨機(jī)投擲一根長為l(l

d)的均勻直針，求針與平行線相交的概率，并計算的值小實例四：蒲豐投針實驗2022/12/2226在畫有許多間距為d的等距平行線的白紙上，隨機(jī)投擲一根長n=100000;L=0.5;d=1;m=0;fori=1:nalpha=rand(1)*pi;y=rand(1)*d/2;ify<=L/2*sin(alpha)m=m+1;endendfprintf('針與平行線相交的頻率為：%f\n',m/n);fprintf('計算出來的pi為：%f\n',2*n*L/(m*d));小實例四源程序2022/12/2227n=100000;小實例四源程序2022/12/1927

已知某變量滿足如圖所示的直方圖分布，求該變量的一個序列，要求該序列滿足該分布。例：病人視網(wǎng)膜疾病術(shù)后觀察時間一般為[5,15]天，其人數(shù)分布如直方圖所示，求滿足該分布的術(shù)后觀察天數(shù)序列。小實例五：直方圖采樣2022/12/2228已知某變量滿足如圖所示的直方圖分布，求該變量的一個序列，要小實例五：直方圖采樣天數(shù)概率（頻率）分布函數(shù)值區(qū)間53/1013/101162/1015/101277/10112/1013816/10128/1014911/10139/10151018/10157/10161115/10172/10171211/10183/1018139/10192/1019145/10197/10110154/1011112022/12/2229小實例五：直方圖采樣天數(shù)概率（頻率）分布函數(shù)值區(qū)間53/10小實例五：直方圖采樣x=[5:15];X=[];%

X模擬所得各觀察天數(shù)的患者y=[3,2,7,16,11,18,15,11,9,5,4];yy=y(1);%yy表示各對應(yīng)天數(shù)的累積人數(shù)值fori=2:length(y)yy=[yy,y(i)+yy(i-1)];endyy=yy/yy(end);%yy表示各對應(yīng)天數(shù)的累積分布函數(shù)值fori=1:10000r=rand(1);d=size(find(r>=yy),2)+1;X(i)=x(d);endhist(X,x)%思考：如何生產(chǎn)指定分布隨機(jī)數(shù)?2022/12/2230小實例五：直方圖采樣x=[5:15];X=[];%隨機(jī)游走模擬隨機(jī)游走是一種基于運用[0，1]區(qū)間的均勻分布隨機(jī)數(shù)序列來進(jìn)行的計算。利用蒙特卡羅方法多次模擬醉漢行走，統(tǒng)計在X處占總次數(shù)的比值，即為概率PN(x)。2022/12/2231隨機(jī)游走模擬隨機(jī)游走是一種基于運用[0，1]區(qū)間的均勻分布隨隨機(jī)游走模擬設(shè)定朝右走的概率P,總步數(shù)N,模擬次數(shù)M,X,j=1,S=0;i=1;x=0;i<=Nj<=M產(chǎn)生隨機(jī)數(shù)確定行走方向計算所在位置xi=i+1j=j+1X==xS=S+1YNYNYN計算S/Mi步數(shù)，j模擬的次數(shù)，S為M次中成功次數(shù)x表示N步后實際位置,X表示指定的位置2022/12/2232隨機(jī)游走模擬設(shè)定朝右走的概率P,總步數(shù)N,模擬次數(shù)M,程序P=0.5;N=10;M=100;X=0;S=0;forj=1:Mx=0;fori=1:Nif(rand(1)<=P)x=x+1;elsex=x-1;endendif(X==x)S=S+1;endendPP=S/M2022/12/2233程序P=0.5;N=10;M=100;X=0;2022機(jī)械零件的可靠度計算可靠度計算主要研究機(jī)械零件在一定分布的隨機(jī)應(yīng)力S（即零件的正常工作應(yīng)力）作用下的可靠程度。應(yīng)力和強(qiáng)度都是隨機(jī)變量，影響應(yīng)力和強(qiáng)度的各種因素也是隨機(jī)變量。因此，機(jī)械零件的可靠度計算過程中，利用蒙特卡羅方法來處理對隨機(jī)量的計算問題，具有特有的優(yōu)勢。問題：設(shè)計一個拉桿，若受外力P均值為20000N，標(biāo)準(zhǔn)差2000N，拉桿的半徑D均值為20mm，標(biāo)準(zhǔn)差0.5mm，材料強(qiáng)度S均值412MPa，標(biāo)準(zhǔn)差15.6MPa，計算其可靠度。2022/12/2234機(jī)械零件的可靠度計算可靠度計算主要研究機(jī)械零件在一定分布的隨機(jī)械零件的可靠度計算問題分析：2022/12/2235機(jī)械零件的可靠度計算問題分析：2022/12/1935機(jī)械零件的可靠度計算理論推導(dǎo)：2022/12/2236機(jī)械零件的可靠度計算理論推導(dǎo)：2022/12/1936機(jī)械零件的可靠度計算理論推導(dǎo)：2022/12/2237機(jī)械零件的可靠度計算理論推導(dǎo)：2022/12/1937機(jī)械零件的可靠度計算蒙特卡羅模擬：通過多次試驗，統(tǒng)計可靠零件個數(shù)，求出其占總零件數(shù)的比值，即為可靠度。設(shè)定總零件數(shù)N=1000，可靠零件數(shù)St=0;令P_mean=20000;P_deta=2000;D_mean=20;D_deta=0.5;Q_mean=412*10^6；Q_deta=15.6*10^6;i=1;i<=N根據(jù)已知分布產(chǎn)生隨機(jī)數(shù)p,d,q;計算應(yīng)力s=(4*p)/(pi*d^2)YNs<=qSt=St+1,i=i+1;計算可靠度：St/N2022/12/2238機(jī)械零件的可靠度計算蒙特卡羅模擬：通過多次試驗，統(tǒng)計可靠零件機(jī)械零件的可靠度計算程序：N=100000;St=0;P_mean=20000;P_deta=2000;D_mean=20*10^-3;D_deta=0.5*10^-3;Q_mean=412*10^6;Q_deta=150.6*10^6;fori=1:Np=P_mean+P_deta*randn(1);%p為正態(tài)分布隨機(jī)數(shù)d=D_mean+D_deta*randn(1);q=Q_mean+Q_deta*randn(1);s=(4*p)/(pi*d^2);ifs<=qSt=St+1;endendSt_p=St/N;fprintf('可信度為：%f\n',St_p);2022/12/2239機(jī)械零件的可靠度計算程序：N=100000;St=0;20排隊問題隨機(jī)模擬排隊論主要研究隨機(jī)服務(wù)系統(tǒng)的工作過程。在排隊系統(tǒng)中，服務(wù)對象的到達(dá)時間和服務(wù)時間都是隨機(jī)的。排隊論通過對隨機(jī)服務(wù)現(xiàn)象的統(tǒng)計研究，找出反映這些隨機(jī)現(xiàn)象平均特性的規(guī)律指標(biāo)，如排隊的長度、等待的時間及服務(wù)利用率。2022/12/2240排隊問題隨機(jī)模擬排隊論主要研究隨機(jī)服務(wù)系統(tǒng)的工作過程。在排隊系統(tǒng)的假設(shè)：（1）顧客源是無窮的；（2）排隊的長度沒有限制；（3）到達(dá)系統(tǒng)的顧客按先后順序依次進(jìn)入服務(wù)，“先到先服務(wù)”。

在某商店有一個售貨員，顧客陸續(xù)來到，售貨員逐個地接待顧客．當(dāng)?shù)絹淼念櫩洼^多時，一部分顧客便須排隊等待，被接待后的顧客便離開商店．設(shè)：1．顧客到來間隔時間服從參數(shù)為0.1

的指數(shù)分布．2．對顧客的服務(wù)時間服從［４，15］上的均勻分布．3．排隊按先到先服務(wù)規(guī)則，隊長無限制．

假定一個工作日為8小時，時間以分鐘為單位。[1]模擬一個工作日內(nèi)完成服務(wù)的個數(shù)及顧客平均等待時間t．[2]模擬100個工作日，求出平均每日完成服務(wù)的個數(shù)及每日顧客的平均等待時間。單服務(wù)員的排隊模型模擬2022/12/2241系統(tǒng)的假設(shè)：在某商店有一個售貨員，顧客陸續(xù)來到，售貨w：總等待時間；ci：第i個顧客的到達(dá)時刻；bi：第i個顧客開始服務(wù)時刻；ei：第i個顧客服務(wù)結(jié)束時刻；xi：第i-1個顧客與第i個顧客之間到達(dá)的間隔時間；yi：對第i個顧客的服務(wù)時間。符號說明2022/12/2242w：總等待時間；符號說明2022/12/1942c1b1c3c4c5c2e1b2e2b3e3b4e4b5ci=ci-1+xi%到達(dá)時刻、時間間隔ei=bi+yi%結(jié)束服務(wù)時刻、開始服務(wù)時刻、服務(wù)時間bi=max(ci,ei-1)t思路分析2022/12/2243c1b1c3c4c5c2e1b2e2b3e3b4e4b5ci初始化：令i=1，ei-1=0，w=0產(chǎn)生間隔時間隨機(jī)數(shù)xi~參數(shù)為0.1的指數(shù)分布ci=xi,bi=xi

產(chǎn)生服務(wù)時間隨機(jī)數(shù)yi~[4,15]的均勻分布ei=bi+yi累計等待時間：w=w+bi-ci準(zhǔn)備下一次服務(wù)：i=i+1產(chǎn)生間隔時間隨機(jī)數(shù)xi~參數(shù)為0.1的指數(shù)分布ci=ci-1+xi

確定開始服務(wù)時間：bi=max(ci,ei-1)bi＞480?YNi=i-1,t=w/i輸出結(jié)果：完成服務(wù)個數(shù)：m=i

平均等待時間：t停止[1]模擬一日ToMatlab(simu1)[2]模擬100日ToMatlab(simu2)b,服務(wù)時刻；c,到達(dá)時刻；e,結(jié)束時刻；2022/12/2244初始化：令i=1，ei-1=0，w=0產(chǎn)生間隔時間隨機(jī)數(shù)xi源代碼simu1function[i,t]=simu1i=1;e=0;w=0;x(i)=random('exp',10);%lambda=1/10c(i)=x(i);b(i)=x(i);while(b(i)<=480)y(i)=11*rand(1)+4;e(i)=b(i)+y(i);w=w+b(i)-c(i);i=i+1;x(i)=random('exp',10);c(i)=c(i-1)+x(i);b(i)=max(c(i),e(i-1));endi=i-1;t=w/i;fprintf('一天完成服務(wù)人數(shù)：%f位\n',i);fprintf('每位顧客平均等待時間：%f分鐘\n',t);2022/12/2245源代碼simu1function[i,t]=simu1源代碼simu2functionsimu2M(1)=0;T(1)=0;fori=1:100[M(i),T(i)]=simu1;endmean_M=mean(M);mean_T=mean(T);fprintf(‘一百天完成服務(wù)人數(shù)：%f位\n',mean_M);fprintf(‘一百天每位顧客平均等待時間：%f分鐘\n',mean_T);2022/12/2246源代碼simu2functionsimu22022/12/用蒙特卡洛法解非線性規(guī)劃問題2022/12/2247用蒙特卡洛法解非線性規(guī)劃問題2022/12/1947基本假設(shè)試驗點的第j個分量xj服從[aj,bj]內(nèi)的均勻分布．符號假設(shè)求解過程先產(chǎn)生一個隨機(jī)數(shù)作為初始試驗點,以后則將上一個試驗點的第j個分量隨機(jī)產(chǎn)生，其它分量不變而產(chǎn)生一新的試驗點．這樣，每產(chǎn)生一個新試驗點只需一個新的隨機(jī)數(shù)分量．當(dāng)K>MAXK或P>MAXP時停止迭代．2022/12/2248基本假設(shè)試驗點的第j個分量xj服從[aj,框圖初始化:給定MAXK,MAXP;k=0,p=0,Q:大整數(shù)xj=aj+R(bj-aj)j=1,2,…,nj=0j=j+1,p=p+1P>MAXP?YNxj=aj+R(bj-aj)gi(X)≥0?i=1,2…nYNj<n?f(X)≥Q?YNX*=X,Q=f(X)k=k+1k>MAXK?YN輸出X,Q,停止YN2022/12/2249框圖初始化:給定MAXK,MAXP;k=0,p=0,Q:大

在Matlab中編程，共需三個Ｍ－文件:randlp.m,mylp.m和

lpconst.m.主程序為randlp.m.%mylp.mfunctionz=mylp(x)%目標(biāo)函數(shù)z=2*x(1)^2+x(2)^2-x(1)*x(2)-8*x(1)-3*x(2);%轉(zhuǎn)化為求最小值問題%lpconst.mfunction

lpc=lpconst(x)%約束條件if3*x(1)+x(2)^2-10<0.5&3*x(1)+x(2)^2-10>-0.5;

%約束條件的誤差為

lpc=1;elselpc=0;end2022/12/2250在Matlab中編程，共需三個Ｍ－文件:randl%randlp.m

function[sol,r1,r2]=randlpdebug=1;a=0;

%試驗點下界b=10;%試驗點上界n=1000;%試驗點個數(shù)r1=unifrnd(a,b,n,1);%[a,b]均勻分布隨機(jī)數(shù)矩陣r2=unifrnd(a,b,n,1);sol=[r1(1)r2(1)];z0=inf;fori=1:nx1=r1(i);x2=r2(i);lpc=lpconst([x1x2]);

iflpc==1z=mylp([x1x2]);

ifz<z0z0=z;

sol=[x1x2];

end

endend2022/12/2251%randlp.m2022/12/1951蒙特卡洛方法預(yù)測最佳搜救區(qū)域問題介紹

根據(jù)搜救目標(biāo)最后已知位置和時間,依據(jù)目標(biāo)海上漂流風(fēng)壓特性,綜合海區(qū)不同時段風(fēng)速、流速等環(huán)境信息及其不確定性影響,預(yù)測目標(biāo)漂移后的搜索區(qū)域,隨時間推移計算目標(biāo)漂流軌跡,同時求出某一具體時刻目標(biāo)在某可能區(qū)域的概率。準(zhǔn)確地預(yù)測漂浮物的漂浮區(qū)域?qū)τ诩霸绱_定遇險對象搜救區(qū)域、成功實施援救具有決定意義。準(zhǔn)確的搜索區(qū)域劃定包含兩個要求:①搜救區(qū)域要能以最大概率包含搜救對象,不至于遺漏搜救對象;②搜救區(qū)域確定盡可能細(xì)致,盡可能小,使搜救力量可以集中在最短時間內(nèi)搜索可能性最高的區(qū)域。2022/12/2252蒙特卡洛方法預(yù)測最佳搜救區(qū)域問題介紹2022/12/1952蒙特卡洛方法預(yù)測最佳搜救區(qū)域影響漂流的因素漂浮物漂流模型:漂浮物漂流位置受事發(fā)位置、總流壓速度、風(fēng)壓、搜索目標(biāo)自身風(fēng)壓特性等因素的影響。事發(fā)位置可由報告得知,包含一定半徑誤差,總流速可以通過海上流速測量設(shè)備獲得,不同海區(qū)海流速度隨時間而變化,包含一定的測量誤差。波浪的影響,對于搜救對象外形大小遠(yuǎn)小于波浪波長的搜救對象可以忽略不計。而在一些特殊區(qū)域,如沿岸潮流、島嶼、潮汐、內(nèi)河水流等因素則需要特殊考慮。

2022/12/2253蒙特卡洛方法預(yù)測最佳搜救區(qū)域影響漂流的因素2022/12/1蒙特卡洛方法預(yù)測最佳搜救區(qū)域遇險事故位置的概率分布遇險位置因事發(fā)條件不同通常可以看作滿足一定分布率的位置分布。當(dāng)遇險對象最后時間已知,失事位置未定,事故位置分布在一個圓形區(qū)域內(nèi),符合最后已知點為圓心,事故位置偏差距離為方差的基點分布;當(dāng)搜尋目標(biāo)沿一定航線航行,最后時間未知,搜尋目標(biāo)位置分布可以看作以航線為基準(zhǔn)線,線的兩側(cè)滿足正態(tài)分布的基線分布。當(dāng)失事位置可能是捕魚區(qū)或其他作業(yè)區(qū),事故的位置一般認(rèn)為是區(qū)域平均分布。對于一些特殊情況,根據(jù)經(jīng)驗劃分區(qū)域遇險位置概率分布。2022/12/2254蒙特卡洛方法預(yù)測最佳搜救區(qū)域遇險事故位置的概率分布2022/蒙特卡洛方法預(yù)測最佳搜救區(qū)域風(fēng)壓作用模型

漂浮物在水面上受風(fēng)力作用和水面下受到海流作用影響而使漂浮物偏離風(fēng)向漂移,形成風(fēng)壓角。風(fēng)壓角可能在風(fēng)速方向順時針右偏,也可能以風(fēng)速方向逆時針左偏,以左偏矢量稱作橫風(fēng)向為負(fù),右偏為正。將風(fēng)壓矢量依照風(fēng)速方向和風(fēng)速垂直方向分解,風(fēng)速方向分解得到矢量稱作順風(fēng)向矢量,垂直風(fēng)速分解得到的矢量稱作橫風(fēng)向矢量DWL=|L|sin(90°-α)(3)CWL=|L|cos(90°-α)(4)式中:L———總風(fēng)壓矢量;α———風(fēng)壓角;DWL———順向風(fēng)壓方向矢量;CWL———橫向風(fēng)壓方向矢量;2022/12/2255蒙特卡洛方法預(yù)測最佳搜救區(qū)域風(fēng)壓作用模型DWL=|L|蒙特卡洛方法預(yù)測最佳搜救區(qū)域漂流位置確定

漂浮物的漂流速度Vdrift滿足:Vdrift=Vcurr+Vleeway(5)式中:Vcurr———海流對地總流速;Vleeway———漂浮物相對水面環(huán)境的速度,即風(fēng)壓速度。

漂流當(dāng)前位置POScurr滿足:POScurr=POSlkp+Vdrift×(tcurr–tlkp)(6)式中:POSlkp———事發(fā)位置;tcurr———當(dāng)前預(yù)測時間;tlkp———事發(fā)時間。2022/12/2256蒙特卡洛方法預(yù)測最佳搜救區(qū)域漂流位置確定2022/12/19流程圖2022/12/2257流程圖2022/12/1957模擬結(jié)果2022/12/2258模擬結(jié)果2022/12/1958模擬結(jié)果2022/12/2259模擬結(jié)果2022/12/1959總結(jié)應(yīng)用范圍：隨機(jī)性問題非隨機(jī)性問題一般步驟：建立模型，畫出流程圖根據(jù)隨機(jī)變量的分布進(jìn)行抽樣模擬統(tǒng)計出概率或求解。2022/12/2260總結(jié)應(yīng)用范圍：2022/12/1960謝謝大家!2022/12/2261謝謝大家!2022/12/1961數(shù)學(xué)建模專題之

蒙特卡羅方法主講：李培巒時間：2014-08-032022/12/2262數(shù)學(xué)建模專題之

蒙特卡羅方法主講：李培巒2022/12/19內(nèi)容提綱1.引言2.MonteCarlo模擬基本思想3.隨機(jī)數(shù)生成函數(shù)4.應(yīng)用實例舉例5.隨機(jī)游走模擬6.機(jī)械零件的可靠度計算7.排隊論模擬8.MonteCarlo方法求解規(guī)劃問題9.MonteCarlo方法預(yù)測搜救區(qū)域2022/12/2263內(nèi)容提綱1.引言2022/12/192引言(Introduction)MonteCarlo方法：蒙特卡羅方法，又稱隨機(jī)模擬方法，屬于計算數(shù)學(xué)的一個分支，它是在上世紀(jì)四十年代中期為了適應(yīng)當(dāng)時原子能事業(yè)的發(fā)展而發(fā)展起來的。亦稱統(tǒng)計模擬方法，statisticalsimulationmethod利用隨機(jī)數(shù)進(jìn)行數(shù)值模擬的方法MonteCarlo名字的由來：MonteCarlo是摩納哥（monaco)最大的城市，該城以賭博聞名。Monte-Carlo,Monaco2022/12/2264引言(Introduction)MonteCarlo方法：引言(Introduction)MonteCarlo方法的起源：二十世紀(jì)四十年代中期，由于科學(xué)技術(shù)的發(fā)展和電子計算機(jī)的發(fā)明，蒙特卡羅方法作為一種獨立的方法被提出來，并首先在核武器的試驗與研制中得到了應(yīng)用(中子的鏈鎖反應(yīng))。但其基本思想并非新穎，可追溯到18世紀(jì)后半葉的蒲豐（Buffon）隨機(jī)投針實驗（1777年）。NicholasMetropolis(1915-1999)JohnVonNeumann(1903-1957)

2022/12/2265引言(Introduction)MonteCarlo方法的引言(Introduction)MonteCarlo方法的應(yīng)用：物理：核物理，熱力學(xué)與統(tǒng)計物理，粒子輸運問題等數(shù)學(xué)：多重積分、解微分方程、非線性方程組求解等工程領(lǐng)域：真空技術(shù)，水力學(xué)，激光技術(shù)等經(jīng)濟(jì)學(xué)領(lǐng)域：期權(quán)定價、項目管理、投資風(fēng)險決策等其他領(lǐng)域：化學(xué)、醫(yī)學(xué)，生物，生產(chǎn)管理、系統(tǒng)科學(xué)、公用事業(yè)等方面。隨著科學(xué)技術(shù)的發(fā)展，其應(yīng)用范圍將更加廣泛。2022/12/2266引言(Introduction)MonteCarlo方法的MonteCarlo方法的基本思想Buffon投針實驗1768年，法國數(shù)學(xué)家ComtedeBuffon利用投針實驗估計的值dL2022/12/2267MonteCarlo方法的基本思想Buffon投針實驗17SolutionThepositioningoftheneedlerelativetonearbylinescanbedescribedwitharandomvectorwhichhascomponents:Therandomvectorisuniformlydistributedontheregion[0,d)×[0,).Accordingly,ithasprobabilitydensityfunction1/d.Theprobabilitythattheneedlewillcrossoneofthelinesisgivenbytheintegral2022/12/2268SolutionThepositioningofthe例1.蒲豐投針問題利用關(guān)系式：求出π值其中Ｎ為投計次數(shù)，n為針與平行線相交次數(shù)。這就是古典概率論中著名的蒲豐氏問題。2022/12/2269例1.蒲豐投針問題利用關(guān)系式：2022/12/198一些人進(jìn)行了實驗，其結(jié)果列于下表：實驗者年份投計次數(shù)π的實驗值沃爾弗(Wolf)185050003.1596斯密思(Smith)185532043.1553?？怂?Fox)189411203.1419拉查里尼(Lazzarini)19013408312/2270實驗者年份投計次數(shù)π的實驗值沃爾弗(Wolf)1850500基本思想

MS基本思想：為了求解數(shù)學(xué)、物理、工程技術(shù)或隨機(jī)服務(wù)系統(tǒng)等方面的問題，首先構(gòu)造一個模型（概率統(tǒng)計模型），使所求問題的解正好是該模型的參數(shù)或特征量或有關(guān)量，然后通過模擬（統(tǒng)計試驗），給出模型參數(shù)或特征量的估計值，最后得出所求問題的近似解。

MS特點：1.方法新穎、應(yīng)用面廣、實用性強(qiáng)；

2.隨機(jī)模擬方面的算法簡單，但計算量大；3.模擬結(jié)果具有隨機(jī)性，精度較低；4.模擬結(jié)果的收斂過程服從概率規(guī)律。2022/12/2271基本思想MS基本思想：為了求解數(shù)學(xué)、物理、工程技術(shù)或隨機(jī)服例1

在我方某前沿防守地域，敵人以一個炮排（含兩門火炮）為單位對我方進(jìn)行干擾和破壞．為躲避我方打擊，敵方對其陣地進(jìn)行了偽裝并經(jīng)常變換射擊地點．經(jīng)過長期觀察發(fā)現(xiàn)，我方指揮所對敵方目標(biāo)的指示有50％是準(zhǔn)確的，而我方火力單位，在指示正確時，有1/3的射擊效果能毀傷敵人一門火炮，有1/6的射擊效果能全部毀傷敵人火炮．現(xiàn)在希望能用某種方式把我方將要對敵人實施的20次打擊結(jié)果顯現(xiàn)出來，確定有效射擊的比率及毀傷敵方火炮的平均值。分析：這是一個概率問題，可以通過理論計算得到相應(yīng)的概率和期望值.但這樣只能給出作戰(zhàn)行動的最終靜態(tài)結(jié)果，而顯示不出作戰(zhàn)行動的動態(tài)過程.為了能顯示我方20次射擊的過程，現(xiàn)采用模擬的方式。舉例2022/12/2272例1在我方某前沿防守地域，敵人以一個炮排（含兩門火炮）為單需要模擬出以下兩件事：

[2]當(dāng)指示正確時，我方火力單位的射擊結(jié)果情況[1]觀察所對目標(biāo)的指示正確與否模擬試驗有兩種結(jié)果，每一種結(jié)果出現(xiàn)的概率都是1/2．因此，可用投擲一枚硬幣的方式予以確定，當(dāng)硬幣出現(xiàn)正面時為指示正確，反之為不正確．模擬試驗有三種結(jié)果：毀傷一門火炮的可能性為1/3(即2/6)，毀傷兩門的可能性為1/6，沒能毀傷敵火炮的可能性為1/2(即3/6)．這時可用投擲骰子的方法來確定：如果出現(xiàn)的是１、２、３三個點：則認(rèn)為沒能擊中敵人；如果出現(xiàn)的是４、５點：則認(rèn)為毀傷敵人一門火炮；若出現(xiàn)的是６點：則認(rèn)為毀傷敵人兩門火炮．問題分析2022/12/2273需要模擬出以下兩件事：[2]i：要模擬的打擊次數(shù)；k1：沒擊中敵人火炮的射擊總數(shù)；k2：擊中敵人一門火炮的射擊總數(shù)；k3：擊中敵人兩門火炮的射擊總數(shù)；E：有效射擊比率；E1：20次射擊平均每次毀傷敵人的火炮數(shù)．符號說明2022/12/2274i：要模擬的打擊次數(shù)；符號說明2022模擬框圖初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子點數(shù)?k1=k1+1k2=k2+1k3=k3+1k1=k1+1i＜20?E=(k2+k3)/20E1=(0*k1+1*k2+2*k3)/20停止硬幣正面?YNNY1，2，34，562022/12/2275模擬框圖初始化:i=0,k1=0,k2=0,k3=0i=i+模擬結(jié)果2022/12/2276模擬結(jié)果2022/12/1915從以上模擬結(jié)果可計算出：2022/12/2277從以上模擬結(jié)果可計算出：2022/12/1916理論計算2022/12/2278理論計算2022/12/1917結(jié)果比較

雖然模擬結(jié)果與理論計算不完全一致，但它卻能更加真實地表達(dá)實際戰(zhàn)斗動態(tài)過程．要使結(jié)果接近理論計算值，必須加大模擬次數(shù)，這就要求使用計算機(jī)模擬了。

用蒙特卡洛方法進(jìn)行計算機(jī)模擬的步驟：[1]設(shè)計一個邏輯框圖，即建立模擬模型．[2]根據(jù)流程圖編寫程序，模擬隨機(jī)現(xiàn)象．可通過生成具有各種概率分布的隨機(jī)數(shù)來模擬隨機(jī)現(xiàn)象，進(jìn)行模擬試驗．[3]分析模擬結(jié)果，給出所求問題的近似解（求解）.2022/12/2279結(jié)果比較雖然模擬結(jié)果與理論計算不完全一致，但它卻能更注：rand(n)=rand(n,n)Matlab中的隨機(jī)數(shù)生成函數(shù)randperm(m)生成一個由1:m

組成的隨機(jī)排列randn(m,n)生成一個滿足正態(tài)m

n

的隨機(jī)矩陣rand(m,n)

生成一個滿足均勻分布的m

n

隨機(jī)矩陣，矩陣的每個元素都在(0,1)

之間。perms(1:n)生成由1:n

組成的全排列，共n!

個結(jié)果2022/12/2280注：rand(n)=rand(n,n)Matlab中的隨機(jī)name

的取值可以是'norm'or'Normal''unif'or'Uniform''poiss'or'Poisson''beta'or'Beta''exp'or'Exponential''gam'or'Gamma''geo'or'Geometric''unid'or'DiscreteUniform'......random('name',A1,A2,A3,M,N)Matlab中的隨機(jī)數(shù)生成函數(shù)2022/12/2281name的取值可以是'norm'or'Normalfix(x):

截尾取整，直接將小數(shù)部分舍去floor(x):

不超過x

的最大整數(shù)ceil(x):

不小于x

的最小整數(shù)round(x):

四舍五入取整Matlab中的取整函數(shù)2022/12/2282fix(x):截尾取整，直接將小數(shù)部分舍去Matla模擬隨機(jī)投擲均勻硬幣，驗證國徽朝上與朝下的概率是否都是1/2

n=10000;%給定試驗次數(shù)m=0;%m表示試驗成功（國徽朝上）的次數(shù)fori=1:nx=randperm(2)-1;%randperm(2）生成1和2的隨機(jī)排列y=x(1);ify==0%0表示國徽朝上，1表示國徽朝下m=m+1;endendfprintf('國徽朝上的頻率為：%f\n',m/n);小實例一：投擲硬幣2022/12/2283模擬隨機(jī)投擲均勻硬幣，驗證國徽朝上與朝下的概率是否都是1隨機(jī)投擲骰子，驗證各點出現(xiàn)的概率是否為1/6

n=10000;m1=0;m2=0;m3=0;m4=0;m5=0;m6=0;fori=1:nx=randperm(6);y=x(1);switchycase1,m1=m1+1;case2,m2=m2+1;case3,m3=m3+1;case4,m4=m4+1;case5,m5=m5+1;otherwise,m6=m6+1;endend...%

輸出結(jié)果小實例二：投擲骰子2022/12/2284隨機(jī)投擲骰子，驗證各點出現(xiàn)的概率是否為1/6n=100

用蒙特卡羅(MonteCarlo)投點法計算

的值n=10000;a=2;m=0;%m表示落入圓內(nèi)的次數(shù)fori=1:nx=rand(1)*a/2;y=rand(1)*a/2;if(x^2+y^2<=(a/2)^2)m=m+1;endendfprintf('計算出來的pi為：%f\n',4*m/n);小實例三：蒙特卡羅投點法2022/12/2285用蒙特卡羅(MonteCarlo)投點法計算小實例三：蒙特卡羅投點法ezplot('x^2+y^2-1',[-1.11.1]);holdonaxisequalplot([-1-111-1],[-111-1-1]);N=0;fork=1:100N_point=10000;xy=(rand(2,N_point)-0.5)*2;d=sqrt(xy(1,:).^2+xy(2,:).^2);N(k)=length(find(d<=1));endp1=sum(N)/(N_point*k)p2=pi*1^2/4pi0=p1*4p_xy=(rand(500,2)-0.5).*2;plot(p_xy(:,1),p_xy(:,2),'.');2022/12/2286小實例三：蒙特卡羅投點法ezplot('x^2+y^2-1'

在畫有許多間距為d

的等距平行線的白紙上，隨機(jī)投擲一根長為l(l

d)的均勻直針，求針與平行線相交的概率，并計算的值小實例四：蒲豐投針實驗2022/12/2287在畫有許多間距為d的等距平行線的白紙上，隨機(jī)投擲一根長n=100000;L=0.5;d=1;m=0;fori=1:nalpha=rand(1)*pi;y=rand(1)*d/2;ify<=L/2*sin(alpha)m=m+1;endendfprintf('針與平行線相交的頻率為：%f\n',m/n);fprintf('計算出來的pi為：%f\n',2*n*L/(m*d));小實例四源程序2022/12/2288n=100000;小實例四源程序2022/12/1927

已知某變量滿足如圖所示的直方圖分布，求該變量的一個序列，要求該序列滿足該分布。例：病人視網(wǎng)膜疾病術(shù)后觀察時間一般為[5,15]天，其人數(shù)分布如直方圖所示，求滿足該分布的術(shù)后觀察天數(shù)序列。小實例五：直方圖采樣2022/12/2289已知某變量滿足如圖所示的直方圖分布，求該變量的一個序列，要小實例五：直方圖采樣天數(shù)概率（頻率）分布函數(shù)值區(qū)間53/1013/101162/1015/101277/10112/1013816/10128/1014911/10139/10151018/10157/10161115/10172/10171211/10183/1018139/10192/1019145/10197/10110154/1011112022/12/2290小實例五：直方圖采樣天數(shù)概率（頻率）分布函數(shù)值區(qū)間53/10小實例五：直方圖采樣x=[5:15];X=[];%

X模擬所得各觀察天數(shù)的患者y=[3,2,7,16,11,18,15,11,9,5,4];yy=y(1);%yy表示各對應(yīng)天數(shù)的累積人數(shù)值fori=2:length(y)yy=[yy,y(i)+yy(i-1)];endyy=yy/yy(end);%yy表示各對應(yīng)天數(shù)的累積分布函數(shù)值fori=1:10000r=rand(1);d=size(find(r>=yy),2)+1;X(i)=x(d);endhist(X,x)%思考：如何生產(chǎn)指定分布隨機(jī)數(shù)?2022/12/2291小實例五：直方圖采樣x=[5:15];X=[];%隨機(jī)游走模擬隨機(jī)游走是一種基于運用[0，1]區(qū)間的均勻分布隨機(jī)數(shù)序列來進(jìn)行的計算。利用蒙特卡羅方法多次模擬醉漢行走，統(tǒng)計在X處占總次數(shù)的比值，即為概率PN(x)。2022/12/2292隨機(jī)游走模擬隨機(jī)游走是一種基于運用[0，1]區(qū)間的均勻分布隨隨機(jī)游走模擬設(shè)定朝右走的概率P,總步數(shù)N,模擬次數(shù)M,X,j=1,S=0;i=1;x=0;i<=Nj<=M產(chǎn)生隨機(jī)數(shù)確定行走方向計算所在位置xi=i+1j=j+1X==xS=S+1YNYNYN計算S/Mi步數(shù)，j模擬的次數(shù)，S為M次中成功次數(shù)x表示N步后實際位置,X表示指定的位置2022/12/2293隨機(jī)游走模擬設(shè)定朝右走的概率P,總步數(shù)N,模擬次數(shù)M,程序P=0.5;N=10;M=100;X=0;S=0;forj=1:Mx=0;fori=1:Nif(rand(1)<=P)x=x+1;elsex=x-1;endendif(X==x)S=S+1;endendPP=S/M2022/12/2294程序P=0.5;N=10;M=100;X=0;2022機(jī)械零件的可靠度計算可靠度計算主要研究機(jī)械零件在一定分布的隨機(jī)應(yīng)力S（即零件的正常工作應(yīng)力）作用下的可靠程度。應(yīng)力和強(qiáng)度都是隨機(jī)變量，影響應(yīng)力和強(qiáng)度的各種因素也是隨機(jī)變量。因此，機(jī)械零件的可靠度計算過程中，利用蒙特卡羅方法來處理對隨機(jī)量的計算問題，具有特有的優(yōu)勢。問題：設(shè)計一個拉桿，若受外力P均值為20000N，標(biāo)準(zhǔn)差2000N，拉桿的半徑D均值為20mm，標(biāo)準(zhǔn)差0.5mm，材料強(qiáng)度S均值412MPa，標(biāo)準(zhǔn)差15.6MPa，計算其可靠度。2022/12/2295機(jī)械零件的可靠度計算可靠度計算主要研究機(jī)械零件在一定分布的隨機(jī)械零件的可靠度計算問題分析：2022/12/2296機(jī)械零件的可靠度計算問題分析：2022/12/1935機(jī)械零件的可靠度計算理論推導(dǎo)：2022/12/2297機(jī)械零件的可靠度計算理論推導(dǎo)：2022/12/1936機(jī)械零件的可靠度計算理論推導(dǎo)：2022/12/2298機(jī)械零件的可靠度計算理論推導(dǎo)：2022/12/1937機(jī)械零件的可靠度計算蒙特卡羅模擬：通過多次試驗，統(tǒng)計可靠零件個數(shù)，求出其占總零件數(shù)的比值，即為可靠度。設(shè)定總零件數(shù)N=1000，可靠零件數(shù)St=0;令P_mean=20000;P_deta=2000;D_mean=20;D_deta=0.5;Q_mean=412*10^6；Q_deta=15.6*10^6;i=1;i<=N根據(jù)已知分布產(chǎn)生隨機(jī)數(shù)p,d,q;計算應(yīng)力s=(4*p)/(pi*d^2)YNs<=qSt=St+1,i=i+1;計算可靠度：St/N2022/12/2299機(jī)械零件的可靠度計算蒙特卡羅模擬：通過多次試驗，統(tǒng)計可靠零件機(jī)械零件的可靠度計算程序：N=100000;St=0;P_mean=20000;P_deta=2000;D_mean=20*10^-3;D_deta=0.5*10^-3;Q_mean=412*10^6;Q_deta=150.6*10^6;fori=1:Np=P_mean+P_deta*randn(1);%p為正態(tài)分布隨機(jī)數(shù)d=D_mean+D_deta*randn(1);q=Q_mean+Q_deta*randn(1);s=(4*p)/(pi*d^2);ifs<=qSt=St+1;endendSt_p=St/N;fprintf('可信度為：%f\n',St_p);2022/12/22100機(jī)械零件的可靠度計算程序：N=100000;St=0;20排隊問題隨機(jī)模擬排隊論主要研究隨機(jī)服務(wù)系統(tǒng)的工作過程。在排隊系統(tǒng)中，服務(wù)對象的到達(dá)時間和服務(wù)時間都是隨機(jī)的。排隊論通過對隨機(jī)服務(wù)現(xiàn)象的統(tǒng)計研究，找出反映這些隨機(jī)現(xiàn)象平均特性的規(guī)律指標(biāo)，如排隊的長度、等待的時間及服務(wù)利用率。2022/12/22101排隊問題隨機(jī)模擬排隊論主要研究隨機(jī)服務(wù)系統(tǒng)的工作過程。在排隊系統(tǒng)的假設(shè)：（1）顧客源是無窮的；（2）排隊的長度沒有限制；（3）到達(dá)系統(tǒng)的顧客按先后順序依次進(jìn)入服務(wù)，“先到先服務(wù)”。

在某商店有一個售貨員，顧客陸續(xù)來到，售貨員逐個地接待顧客．當(dāng)?shù)絹淼念櫩洼^多時，一部分顧客便須排隊等待，被接待后的顧客便離開商店．設(shè)：1．顧客到來間隔時間服從參數(shù)為0.1

的指數(shù)分布．2．對顧客的服務(wù)時間服從［４，15］上的均勻分布．3．排隊按先到先服務(wù)規(guī)則，隊長無限制．

假定一個工作日為8小時，時間以分鐘為單位。[1]模擬一個工作日內(nèi)完成服務(wù)的個數(shù)及顧客平均等待時間t．[2]模擬100個工作日，求出平均每日完成服務(wù)的個數(shù)及每日顧客的平均等待時間。單服務(wù)員的排隊模型模擬2022/12/22102系統(tǒng)的假設(shè)：在某商店有一個售貨員，顧客陸續(xù)來到，售貨w：總等待時間；ci：第i個顧客的到達(dá)時刻；bi：第i個顧客開始服務(wù)時刻；ei：第i個顧客服務(wù)結(jié)束時刻；xi：第i-1個顧客與第i個顧客之間到達(dá)的間隔時間；yi：對第i個顧客的服務(wù)時間。符號說明2022/12/22103w：總等待時間；符號說明2022/12/1942c1b1c3c4c5c2e1b2e2b3e3b4e4b5ci=ci-1+xi%到達(dá)時刻、時間間隔ei=bi+yi%結(jié)束服務(wù)時刻、開始服務(wù)時刻、服務(wù)時間bi=max(ci,ei-1)t思路分析2022/12/22104c1b1c3c4c5c2e1b2e2b3e3b4e4b5ci初始化：令i=1，ei-1=0，w=0產(chǎn)生間隔時間隨機(jī)數(shù)xi~參數(shù)為0.1的指數(shù)分布ci=xi,bi=xi

產(chǎn)生服務(wù)時間隨機(jī)數(shù)yi~[4,15]的均勻分布ei=bi+yi累計等待時間：w=w+bi-ci準(zhǔn)備下一次服務(wù)：i=i+1產(chǎn)生間隔時間隨機(jī)數(shù)xi~參數(shù)為0.1的指數(shù)分布ci=ci-1+xi

確定開始服務(wù)時間：bi=max(ci,ei-1)bi＞480?YNi=i-1,t=w/i輸出結(jié)果：完成服務(wù)個數(shù)：m=i

平均等待時間：t停止[1]模擬一日ToMatlab(simu1)[2]模擬100日ToMatlab(simu2)b,服務(wù)時刻；c,到達(dá)時刻；e,結(jié)束時刻；2022/12/22105初始化：令i=1，ei-1=0，w=0產(chǎn)生間隔時間隨機(jī)數(shù)xi源代碼simu1function[i,t]=simu1i=1;e=0;w=0;x(i)=random('exp',10);%lambda=1/10c(i)=x(i);b(i)=x(i);while(b(i)<=480)y(i)=11*rand(1)+4;e(i)=b(i)+y(i);w=w+b(i)-c(i);i=i+1;x(i)=random('exp',10);c(i)=c(i-1)+x(i);b(i)=max(c(i),e(i-1));endi=i-1;t=w/i;fprintf('一天完成服務(wù)人數(shù)：%f位\n',i);fprintf('每位顧客平均等待時間：%f分鐘\n',t);2022/12/22106源代碼simu1function[i,t]=simu1源代碼simu2functionsimu2M(1)=0;T(1)=0;fori=1:100[M(i),T(i)]=simu1;endmean_M=mean(M);mean_T=mean(T);fprintf(‘一百天完成服務(wù)人數(shù)：%f位\n',mean_M);fprintf(‘一百天每位顧客平均等待時間：%f分鐘\n',mean_T);2022/12/22107源代碼simu2functionsimu22022/12/用蒙特卡洛法解非線性規(guī)劃問題2022/12/22108用蒙特卡洛法解非線性規(guī)劃問題2022/12/1947基本假設(shè)試驗點的第j個分量xj服從[aj,bj]內(nèi)的均勻分布．符號假設(shè)求解過程先產(chǎn)生一個隨機(jī)數(shù)作為初始試驗點,以后則將上一個試驗點的第j個分量隨機(jī)產(chǎn)生，其它分量不變而產(chǎn)生一新的試驗點．這樣，每產(chǎn)生一個新試驗點只需一個新的隨機(jī)數(shù)分量．當(dāng)K>MAXK或P>MAXP時停止迭代．2022/12/22109基本假設(shè)試驗點的第j個分量xj服從[aj,框圖初始化:給定MAXK,MAXP;k=0,p=0,Q:大整數(shù)xj=aj+R(bj-aj)j=1,2,…,nj=0j=j+1,p=p+1P>MAXP?YNxj=aj+R(bj-aj)gi(X)≥0?i=1,2…nYNj<n?f(X)≥Q?YNX*=X,Q=f(X)k=k+1k>MAXK?YN輸出X,Q,停止YN2022/12/22110框圖初始化:給定MAXK,MAXP;k=0,p=0,Q:大

在Matlab中編程，共需三個Ｍ－文件:randlp.m,mylp.m和

lpconst.m.主程序為randlp.m.%mylp.mfunctionz=mylp(x)%目標(biāo)函數(shù)z=2*x(1)^2+x(2)^2-x(1)*x(2)-8*x