計(jì)算機(jī)模擬版本3

1、計(jì)算機(jī)模擬計(jì)算機(jī)模擬主講： 何仁斌概述 計(jì)算機(jī)模擬是用計(jì)算機(jī)對(duì)系統(tǒng)系統(tǒng)的結(jié)構(gòu)和行為進(jìn)行動(dòng)態(tài)演示，以評(píng)價(jià)或預(yù)測(cè)系統(tǒng)的行為效果，為決策者提供信息的一種方法。它是解決較復(fù)雜的實(shí)際問(wèn)題的一條有效途徑。 計(jì)算機(jī)模擬也可以說(shuō)是用計(jì)算機(jī)程序計(jì)算機(jī)程序直接建立真實(shí)系統(tǒng)的模型，并通過(guò)計(jì)算了解系統(tǒng)隨時(shí)間變化的行為或特征。 應(yīng)用領(lǐng)域：航空、機(jī)電、冶金、社會(huì)經(jīng)濟(jì)、交通運(yùn)輸、生態(tài)系統(tǒng)等。 計(jì)算機(jī)模擬分為連續(xù)系統(tǒng)連續(xù)系統(tǒng)模擬和離散系統(tǒng)離散系統(tǒng)模擬。 狀態(tài)隨時(shí)間連續(xù)變化的系統(tǒng)稱(chēng)為連續(xù)系統(tǒng)。通常該系統(tǒng)的模型一般可以用微分方程的形式表達(dá)，通過(guò)一些物理機(jī)理推導(dǎo)出來(lái)。模擬結(jié)果往往是近似的。例如，追逐問(wèn)題，濃度問(wèn)題。一、連續(xù)系統(tǒng)(x

2、1(t),xn(t)1、追逐問(wèn)題、追逐問(wèn)題abbccdda狀態(tài)變量a(x(t), y(t) 假設(shè)：v=1 m/s, 追逐的方向是對(duì)準(zhǔn)追逐目標(biāo)a(0,10)b(10,10)d(0,0)c(0,10) 試確定四人追逐的運(yùn)動(dòng)軌跡。a(x(t+t), y(t+ t)模擬方法：模擬方法：1. 建立平面直角坐標(biāo)系；建立平面直角坐標(biāo)系；2. 以時(shí)間間隔以時(shí)間間隔t 進(jìn)行采樣，并且計(jì)算各個(gè)進(jìn)行采樣，并且計(jì)算各個(gè) 時(shí)刻下的狀態(tài)：時(shí)刻下的狀態(tài)： a: (x1(t), y1(t) (x1(t+ t ), y1(t+ t ) b: (x2(t), y2(t) (x2(t+ t ), y2(t+ t ) (x1(t+

3、t ), y1(t+ t ) （由幾何原理）（由幾何原理）(x1(t)+v t cos(), y1(t) +v t sin()2122121212)()(,sin,cosyyxxddyydxxx1,y1x2,y2在 t時(shí)刻下3. 選取足夠小的選取足夠小的t ，模擬到任意兩人的距，模擬到任意兩人的距 離小于離小于v t 為止。為止。4. 連接四人在各時(shí)刻下的位置，就得到所求的連接四人在各時(shí)刻下的位置，就得到所求的 運(yùn)動(dòng)軌跡。運(yùn)動(dòng)軌跡。 matlab (simu2.m)v=1;dt=0.05;d=20; % x=x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8)x=0

4、0 0 10 10 10 10 0;x(9)=x(1);x(10)=x(2);hold onaxis(equal);axis(0 10 0 10);while(d0.1) for i=1:2:7 d=sqrt(x(i)-x(i+2)2+(x(i+1)-x(i+3)2); x(i)=x(i)+v*dt*(x(i+2)-x(i)/d; x(i+1)=x(i+1)+v*dt*(x(i+3)-x(i+1)/d; plot(x(i),x(i+1),.) end x(9)=x(1);x(10)=x(2);endhold onabcd 某軍的一導(dǎo)彈基地發(fā)現(xiàn)正北方向某軍的一導(dǎo)彈基地發(fā)現(xiàn)正北方向120km處處海

5、面上有敵艇一艘以海面上有敵艇一艘以90km/h的速度向正東方向的速度向正東方向行駛。該基地立即發(fā)射導(dǎo)彈跟蹤追擊敵艇，導(dǎo)彈行駛。該基地立即發(fā)射導(dǎo)彈跟蹤追擊敵艇，導(dǎo)彈速度為速度為450km/h，自動(dòng)導(dǎo)航系統(tǒng)使導(dǎo)彈在任一時(shí)，自動(dòng)導(dǎo)航系統(tǒng)使導(dǎo)彈在任一時(shí)刻都能對(duì)準(zhǔn)敵艇。試問(wèn)導(dǎo)彈在何時(shí)何處擊中敵艇？刻都能對(duì)準(zhǔn)敵艇。試問(wèn)導(dǎo)彈在何時(shí)何處擊中敵艇？2、導(dǎo)彈跟蹤問(wèn)題、導(dǎo)彈跟蹤問(wèn)題p(x(t), y(t)oxm(vet, h)ya(0,h)正東方向正北方向h =120（千米）90ev450wvkm/hkm/h敵導(dǎo)(xk(t), yk(t)simu1.mstest1.mhuatu.m 微分方程建模計(jì)算機(jī)模擬消去t：微

6、分方程建模xtyxy90120)tan(dd1）2）450)()(45022dtdydtdxsdtds表示弧長(zhǎng)0)0( , 0)0(450901)()120(222xxdydxydyxd x=dsolve(d2x*(120-y)/sqrt(dx)2+1)=90/450,dx(0)=0,x(0)=0,y)微分方程求解simplify(x) ans=(45*120(1/5)*(y - 120)(4/5) - 1800*(-1)(4/5) + (-1)(3/5)*202500(1/5)*(y - 120)(6/5)/(72*(-(-1)(3/5)(1/2) -(45*120(1/5)*(y - 12

7、0)(4/5) - 1800*(-1)(4/5) + (-1)(3/5)*202500(1/5)*(y - 120)(6/5)/(72*(-(-1)(3/5)(1/2)y=120eval(x)ans = 25.0000 + 0.0000i -25.0000 - 0.0000i如何決定導(dǎo)彈位置 p2(x2, y2)?計(jì)算機(jī)模擬1)當(dāng)t = 0時(shí)，導(dǎo)彈位置o(0, 0); 敵艇位置a(0, 120)；2)當(dāng)t= 時(shí)，導(dǎo)彈位置p1(x1, y1); 敵艇位置m1(90 , 120)；3)當(dāng)t=2時(shí)，導(dǎo)彈位置p2(x2, y2); 敵艇位置m2(902, 120)；11190120)tan(xy112

8、112sin450cos450yyxx212111212111)120()90(120sin,)120()90(90cosyxyyxxp(x(t), y(t)oxm(vet, h)ya(0,h)正東方向正北方向h =120（千米）90ev450wvkm/hkm/h敵導(dǎo)(xk(t), yk(t) function num,y_j,l,t=simu1(x,y,t,eps)k=0;while k1000 (simu1.m) p=90*k*t-x; q=120-y; d1=p/(p2+q2)0.5; d2=q/(p2+q2)0.5; x=x+450*t*d1; y=y+450*t*d2; if (ab

9、s(q)0.1) % 終止條件終止條件() for k=1:280 p=90*k*t-x(1); q=120-x(2); d1=p/(p2+q2)0.5; d2=q/(p2+q2)0.5; x(1)=x(1)+450*t*d1; x(2)=x(2)+450*t*d2; x1(1)=90*k*t; x1(2)=120; h1=line(color,0,0.2,0.4,linewidth,2); h2=line(color,0,0.6,0.9,linewidth,3); set(h1,xdata,x1(1),ydata,x1(2); set(h2,xdata,x(1),ydata,x(2); en

10、dendhold on（2 2）如果當(dāng)基地發(fā)射導(dǎo)彈的同時(shí)，敵艇立即）如果當(dāng)基地發(fā)射導(dǎo)彈的同時(shí)，敵艇立即由儀器發(fā)覺(jué)。假定敵艇為高速快艇，它即刻由儀器發(fā)覺(jué)。假定敵艇為高速快艇，它即刻以以135135千米千米/ /小時(shí)的速度與導(dǎo)彈方向成垂直的小時(shí)的速度與導(dǎo)彈方向成垂直的方向逃逸，問(wèn)導(dǎo)彈何時(shí)何地?fù)糁袛惩Х较蛱右?，?wèn)導(dǎo)彈何時(shí)何地?fù)糁袛惩?？不停地改變逃跑策略，運(yùn)動(dòng)軌跡如何？不停地改變逃跑策略，運(yùn)動(dòng)軌跡如何？ 修正敵艇速度：修正敵艇速度：90 13590 135千米千米/ /小時(shí)小時(shí)思考：二、離散系統(tǒng) 離散系統(tǒng)離散系統(tǒng)是指系統(tǒng)狀態(tài)只在有限的時(shí)間點(diǎn)或可數(shù)的時(shí)間點(diǎn)上發(fā)生變化的系統(tǒng)。并假設(shè)離散系統(tǒng)狀態(tài)的變化是在一

11、個(gè)時(shí)間點(diǎn)上瞬間完成。模型一般用流程圖或網(wǎng)絡(luò)來(lái)表示。其間可能涉及到隨機(jī)事件等。 關(guān)鍵：模擬步驟、數(shù)據(jù)收集、模擬時(shí)鐘在我方某前沿防守地域，敵人以一個(gè)炮排（含兩門(mén)火炮）為單位對(duì)我方進(jìn)行干擾和破壞為躲避我方打擊，敵方對(duì)其陣地進(jìn)行了偽裝并經(jīng)常變換射擊地點(diǎn) 經(jīng)過(guò)長(zhǎng)期觀(guān)察發(fā)現(xiàn)，我方指揮所對(duì)敵方目標(biāo)的指示有50是準(zhǔn)確的，而我方火力單位，在指示正確時(shí)，有1/3的射擊效果能毀傷敵人一門(mén)火炮，有1/6的射擊效果能全部消滅敵人 現(xiàn)在希望能用某種方式把我方將要對(duì)敵人實(shí)施的20次打擊結(jié)果顯現(xiàn)出來(lái)，確定有效射擊的比率及毀傷敵方火炮的平均值。分析分析： 這是一個(gè)概率問(wèn)題，可以通過(guò)理論計(jì)算得到相應(yīng)的概率和期望值.但這樣只能給出

12、作戰(zhàn)行動(dòng)的最終靜態(tài)結(jié)果，而顯示不出作戰(zhàn)行動(dòng)的動(dòng)態(tài)過(guò)程. 為了能顯示我方20次射擊的過(guò)程，現(xiàn)采用模擬的方式。射擊問(wèn)題 需要模擬出以下兩件事： 1. 問(wèn)題分析問(wèn)題分析2 2 當(dāng)指示正確時(shí)，我方火力單位的射擊結(jié)果情況當(dāng)指示正確時(shí)，我方火力單位的射擊結(jié)果情況1 1 觀(guān)察所對(duì)目標(biāo)的指示正確與否觀(guān)察所對(duì)目標(biāo)的指示正確與否模擬試驗(yàn)有兩種結(jié)果，每一種結(jié)果出現(xiàn)的概率都是1/2 因此，可用投擲一枚硬幣的方式予以確定可用投擲一枚硬幣的方式予以確定，當(dāng)硬幣出現(xiàn)正面時(shí)為指示正確，反之為不正確 模擬試驗(yàn)有三種結(jié)果：毀傷一門(mén)火炮的可能性為1/3(即2/6)，毀傷兩門(mén)的可能性為1/6，沒(méi)能毀傷敵火炮的可能性為1/2(即3/6

13、) 這時(shí)可用投擲骰子的方法來(lái)確定可用投擲骰子的方法來(lái)確定：如果出現(xiàn)的是、三個(gè)點(diǎn)：則認(rèn)為沒(méi)能擊中敵人；如果出現(xiàn)的是、點(diǎn)：則認(rèn)為毀傷敵人一門(mén)火炮；若出現(xiàn)的是點(diǎn)：則認(rèn)為毀傷敵人兩門(mén)火炮2. 符號(hào)假設(shè)符號(hào)假設(shè)i：要模擬的打擊次數(shù)； k1：沒(méi)擊中敵人火炮的射擊總數(shù)； k2：擊中敵人一門(mén)火炮的射擊總數(shù)；k3：擊中敵人兩門(mén)火炮的射擊總數(shù)e：有效射擊比率； e1：20次射擊平均每次毀傷敵人的火炮數(shù)3. 模擬框圖模擬框圖初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子點(diǎn)數(shù)?k1=k1+1k2=k2+1k3=k3+1k1=k1+1i20?e=(k2+k3)/20 e1=0*k1/20+1*k2/20+2

14、*k3/20停止硬幣正面?ynny1，2，34，564. 模擬結(jié)果模擬結(jié)果消滅敵人火炮數(shù)消滅敵人火炮數(shù)試驗(yàn)試驗(yàn)序號(hào)序號(hào)投硬幣投硬幣結(jié)果結(jié)果指示指示正確正確指示指示不正確不正確擲骰子擲骰子結(jié)果結(jié)果正正反正正反正正反反消滅敵人火炮數(shù)消滅敵人火炮數(shù)試驗(yàn)試驗(yàn)序號(hào)序號(hào)投硬幣投硬幣結(jié)果結(jié)果指示指示正確正確指示指示不正確不正確擲骰子擲骰子結(jié)果結(jié)果正反正反正正正正反正從以上模擬結(jié)果可計(jì)算出： e=7/20=0.35 20322041201301e=0.55. 理論計(jì)算理論計(jì)算設(shè)：觀(guān)察所對(duì)目標(biāo)指示正確確觀(guān)察所對(duì)目標(biāo)指示不正10ja0：射中敵方火炮的事件；a1：射中敵方一門(mén)火炮的事件；a2：射中敵方兩門(mén)火炮的事件

15、則由全概率公式：e = p(a0) = p(j=0)p(a0j=0) + p(j=1)p(a0j=1) = 25. 02121021p(a1) = p(j=0)p(a1j=0) + p(j=1)p(a1j=1) = 613121021p(a2) = p(j=0)p(a2j=0) + p(j=1)p(a2j=1) = 1216121021e1 = 33. 012126116. 結(jié)果比較結(jié)果比較 理論計(jì)算和模擬結(jié)果的比較分類(lèi)項(xiàng)目無(wú)效射擊有效射擊平均值模擬理論 雖然模擬結(jié)果與理論計(jì)算不完全一致，但它卻能更加真實(shí)地表達(dá)實(shí)雖然模擬結(jié)果與理論計(jì)算不完全一致，但它卻能更加真實(shí)地表達(dá)實(shí)際戰(zhàn)斗動(dòng)態(tài)過(guò)程際戰(zhàn)斗動(dòng)態(tài)

16、過(guò)程 蒙特卡羅方法蒙特卡羅方法(monte carlo method) 蒙特卡羅方法蒙特卡羅方法(monte carlo method)，或稱(chēng)計(jì)，或稱(chēng)計(jì)算機(jī)隨機(jī)模擬方法，是一種基于算機(jī)隨機(jī)模擬方法，是一種基于“隨機(jī)數(shù)隨機(jī)數(shù)”的計(jì)的計(jì)算方法。源于美國(guó)在第二次世界大戰(zhàn)研制原子彈算方法。源于美國(guó)在第二次世界大戰(zhàn)研制原子彈的的“曼哈頓計(jì)劃曼哈頓計(jì)劃”，該計(jì)劃的主持人之一數(shù)學(xué)家，該計(jì)劃的主持人之一數(shù)學(xué)家馮馮諾伊曼諾伊曼用馳名世界的用馳名世界的賭城賭城摩納哥的摩納哥的monte carlo來(lái)命名這種方法，為它蒙上了一層神秘色來(lái)命名這種方法，為它蒙上了一層神秘色彩。彩。 蒙特卡羅方法的蒙特卡羅方法的基本思想

17、基本思想很早以前就被人們所很早以前就被人們所發(fā)現(xiàn)和利用。早在發(fā)現(xiàn)和利用。早在17世紀(jì)，人們就知道世紀(jì)，人們就知道用事件發(fā)用事件發(fā)生的生的“頻率頻率”來(lái)決定事件的來(lái)決定事件的“概率概率”。19世紀(jì)人世紀(jì)人們用們用蒲豐投針蒲豐投針的方法來(lái)計(jì)算圓周率的方法來(lái)計(jì)算圓周率，上世紀(jì)，上世紀(jì)40年年代電子計(jì)算機(jī)的出現(xiàn)，特別是近年來(lái)代電子計(jì)算機(jī)的出現(xiàn)，特別是近年來(lái)高速電子計(jì)高速電子計(jì)算機(jī)算機(jī)的出現(xiàn)，使得用數(shù)學(xué)方法在計(jì)算機(jī)上大量、的出現(xiàn)，使得用數(shù)學(xué)方法在計(jì)算機(jī)上大量、快速地模擬這樣的試驗(yàn)成為可能?？焖俚啬M這樣的試驗(yàn)成為可能。蒲豐投針實(shí)驗(yàn)近似計(jì)算圓周率蒲豐投針實(shí)驗(yàn)：蒲豐投針實(shí)驗(yàn)： 法國(guó)科學(xué)家蒲豐法國(guó)科學(xué)家蒲豐(

18、buffon)在在1777年提年提出的蒲豐投針實(shí)驗(yàn)是早期幾何概率一個(gè)出的蒲豐投針實(shí)驗(yàn)是早期幾何概率一個(gè)非常著名的例子。蒲豐投針實(shí)驗(yàn)的重要非常著名的例子。蒲豐投針實(shí)驗(yàn)的重要性并非是為了求得比其它方法更精確的性并非是為了求得比其它方法更精確的值，而是它值，而是它開(kāi)創(chuàng)了開(kāi)創(chuàng)了使用使用隨機(jī)數(shù)處理確定性數(shù)學(xué)問(wèn)隨機(jī)數(shù)處理確定性數(shù)學(xué)問(wèn)題題的先河，是用的先河，是用偶然性方法去解決確定性計(jì)算偶然性方法去解決確定性計(jì)算的的前導(dǎo)，由此可以領(lǐng)略到從前導(dǎo)，由此可以領(lǐng)略到從“概率土壤概率土壤”上開(kāi)出的上開(kāi)出的一朵瑰麗的鮮花一朵瑰麗的鮮花蒙特卡羅方法蒙特卡羅方法(mc) 蒲豐投針實(shí)驗(yàn)蒲豐投針實(shí)驗(yàn)可歸結(jié)為下面的數(shù)學(xué)問(wèn)題可歸結(jié)

19、為下面的數(shù)學(xué)問(wèn)題：平面：平面上畫(huà)有距離為上畫(huà)有距離為a的一些平行線(xiàn)，向平面上任意投一的一些平行線(xiàn)，向平面上任意投一根長(zhǎng)為根長(zhǎng)為l(la)的針，假設(shè)針落在任意位置的可能性的針，假設(shè)針落在任意位置的可能性相同，試求相同，試求針與平行線(xiàn)相交的概率針與平行線(xiàn)相交的概率p(從而求從而求)蒲豐投針實(shí)驗(yàn)近似計(jì)算圓周率蒲豐投針實(shí)驗(yàn)：蒲豐投針實(shí)驗(yàn)： 如右圖所示，以如右圖所示，以m表示針落下后的中點(diǎn)，表示針落下后的中點(diǎn)，以以x表示表示m到最近一條到最近一條平行線(xiàn)的距離，以平行線(xiàn)的距離，以表示針與此線(xiàn)的交角：表示針與此線(xiàn)的交角：針落地的所有可能結(jié)果滿(mǎn)足：針落地的所有可能結(jié)果滿(mǎn)足：其樣本空間視作矩形區(qū)域其樣本空間視作

20、矩形區(qū)域, 面積是面積是:針與平行線(xiàn)相交的條件：針與平行線(xiàn)相交的條件：它是樣本空間它是樣本空間子集子集a，面積是：，面積是：syms l phi; int(l/2*sin(phi),phi,0,pi); %ans=l因此，針與平行線(xiàn)相交的概率為：因此，針與平行線(xiàn)相交的概率為：從而有：從而有： 特別當(dāng)特別當(dāng) 時(shí)時(shí)02,0 xa( )2sa 0sin2,0 xl0( )sin2s aldl( )/ ( )2 /ps asl a 2l ap2al1 p蒲豐投針實(shí)驗(yàn)近似計(jì)算圓周率蒲豐投針實(shí)驗(yàn)近似計(jì)算圓周率蒲豐投針實(shí)驗(yàn)的計(jì)算機(jī)模擬：蒲豐投針實(shí)驗(yàn)的計(jì)算機(jī)模擬：format long; %設(shè)置設(shè)置15位顯示

21、精度位顯示精度a=1; l=0.6; %兩平行線(xiàn)間的寬度和針長(zhǎng)兩平行線(xiàn)間的寬度和針長(zhǎng)figure; axis(0,pi,0,a/2); %初始化繪圖板初始化繪圖板set(gca,nextplot,add); %初始化繪圖方式為疊加初始化繪圖方式為疊加counter=0; n=2010; %初始化計(jì)數(shù)器和設(shè)定投針次數(shù)初始化計(jì)數(shù)器和設(shè)定投針次數(shù)x=unifrnd(0,a/2,1,n); phi=unifrnd(0,pi,1,n); %樣本空間樣本空間for i=1:n if x(i)l*sin(phi(i)/2 %滿(mǎn)足此條件表示針與線(xiàn)的相交滿(mǎn)足此條件表示針與線(xiàn)的相交 plot(phi(i),x(i

22、),r.); frame(i)=getframe; %描點(diǎn)并取幀描點(diǎn)并取幀 counter=counter+1; %統(tǒng)計(jì)針與線(xiàn)相交的次數(shù)統(tǒng)計(jì)針與線(xiàn)相交的次數(shù) endendfren=counter/n; pihat=2*l/(a*fren) %用頻率近似計(jì)算用頻率近似計(jì)算蒲豐投針實(shí)驗(yàn)近似計(jì)算圓周率蒲豐投針實(shí)驗(yàn)的計(jì)算機(jī)模擬：蒲豐投針實(shí)驗(yàn)的計(jì)算機(jī)模擬： 意大利數(shù)學(xué)家拉澤里尼得到了準(zhǔn)確到意大利數(shù)學(xué)家拉澤里尼得到了準(zhǔn)確到6位小數(shù)的位小數(shù)的值，不過(guò)他的實(shí)驗(yàn)因?yàn)橹?，不過(guò)他的實(shí)驗(yàn)因?yàn)樘珳?zhǔn)確太準(zhǔn)確而受到了質(zhì)疑而受到了質(zhì)疑蒲豐投針實(shí)驗(yàn)計(jì)算圓周率蒙特卡羅投點(diǎn)法是蒲豐投針實(shí)驗(yàn)的推廣：蒙特卡羅投點(diǎn)法是蒲豐投針實(shí)驗(yàn)的推廣

23、： 在一個(gè)邊長(zhǎng)為在一個(gè)邊長(zhǎng)為a的正方形內(nèi)隨機(jī)投點(diǎn)，該點(diǎn)落在此正方形的的正方形內(nèi)隨機(jī)投點(diǎn)，該點(diǎn)落在此正方形的內(nèi)切圓中的概率應(yīng)為該內(nèi)切圓與正方形的面積比值，即內(nèi)切圓中的概率應(yīng)為該內(nèi)切圓與正方形的面積比值，即n=10000; a=2; m=0; for i=1:n x=rand(1)*a; y=rand(1)*a; if ( (x-a/2)2+(y-a/2)2 = (a/2)2 ) m=m+1; endenddisp(近似為: ,num2str(4*m/n);%以下為向量化程序以下為向量化程序x=rand(1,10000);y=rand(1,10000);num2str(length(find(x-

24、1).2+(y-1).2=1)*4/10000)22/ 2:/ 4aaxyo(a/2,a/2)一些monte carlo思想思想的例子例1(相遇問(wèn)題) 甲、乙兩船在24小時(shí)內(nèi)獨(dú)立地隨機(jī)到達(dá)碼頭. 設(shè)兩船到達(dá)碼頭時(shí)刻都服從0,24上的均勻分布，如果甲船到達(dá)碼頭后停留2小時(shí),乙船到達(dá)碼頭后停留1小時(shí).問(wèn)兩船相遇的概率有多大？（可用幾何概率，此處略） 分析：如果知道甲、乙兩船到達(dá)的時(shí)刻x和y，兩船能相遇的條件就是：兩船到達(dá)的時(shí)刻x和y可以隨機(jī)生成，生成一組到達(dá)時(shí)刻，可以確定是否能相遇。如果重復(fù)很多次，統(tǒng)計(jì)相遇的比例就可近似為相遇的概率。21xyxyxy且這也是這也是monte carlo思想思想.兩

25、船到達(dá)碼頭時(shí)刻服從0,24上的均勻分布，甲船停留2小時(shí),乙船停留1小時(shí),相遇概率？方法一 利用軟件產(chǎn)生一組0，24上的隨機(jī)數(shù)xi:10.4 5.3 0.3 9.0 20.1 6.8 7.8 5.4分別代表甲船到達(dá)碼頭時(shí)間； 再產(chǎn)生一組0，24上的隨機(jī)數(shù)yi:6.0 3.4 8.1 17.5 0.8 13.4 14.0分別代表乙船到達(dá)碼頭時(shí)間； 把（xi，yi）(i1，2)當(dāng)作一天中甲船乙船到達(dá)的時(shí)間，統(tǒng)計(jì)出能相遇的頻數(shù)，計(jì)算出相遇的頻率做為相遇的概率。可以使用軟件（excel、c等實(shí)現(xiàn)）兩船到達(dá)碼頭時(shí)刻服從0,24上的均勻分布，甲船停留2小時(shí),乙船停留1小時(shí),相遇概率？方法二clc; n=10

26、00; c=0; %模擬次數(shù)n，相遇次數(shù)c清零for i=1:n %重復(fù)n次到達(dá)時(shí)間 x=unifrnd(0,24,1,1); %甲船到達(dá)時(shí)間x(隨機(jī)數(shù)) y=unifrnd(0,24,1,1); %乙船到達(dá)時(shí)間y(隨機(jī)數(shù)) if(x=y & y=x+2) |(y=x & x=y+1) c=c+1; %如果能相遇，則計(jì)數(shù)器加1 endendp=c/n %顯示相遇的概率近似值注意：每次運(yùn)行的結(jié)果一般都不一樣%計(jì)算機(jī)模擬程序 （報(bào)童訣竅的簡(jiǎn)化版）報(bào)童每天清晨從報(bào)社購(gòu)進(jìn)報(bào)紙零售，晚上退回沒(méi)有賣(mài)掉的報(bào)紙.若每份報(bào)紙的購(gòu)進(jìn)價(jià)為b=0.75元，售出價(jià)為a=1元，退回價(jià)為c=0.6元.每天需

27、求量x是離散型隨機(jī)變量，其分布為例2x500510520p0.340.360.30問(wèn)：如果報(bào)童每天購(gòu)進(jìn)報(bào)紙為n=510份，每天的平均利潤(rùn)是多少？方法二：如果我們知道每天的需求量，可直接計(jì)算利潤(rùn)。而每天需求量可以按分布生成（隨機(jī)模擬思想）方法一：概率方法（略）售出價(jià)a=1，購(gòu)進(jìn)價(jià)b=0.75，退回價(jià)c=0.6，購(gòu)進(jìn)數(shù)量n=510份需求量x500510520概率p0.34 0.36 0.30按照需求量的分布規(guī)律，隨機(jī)生成n=20個(gè)數(shù)據(jù):510 520 500 510 520 500 500 500 500510 500 500 500 520 510 520 510 510代表20天的需求量，計(jì)算

28、出報(bào)童在這20天的總利潤(rùn)和平均利潤(rùn)，用平均利潤(rùn)來(lái)近似報(bào)童的平均收入。這也是這也是monte carlo方法方法.問(wèn)問(wèn)題題如何按分布規(guī)律產(chǎn)生隨機(jī)數(shù)據(jù)？隨機(jī)數(shù)據(jù)很多時(shí)，如何編程？報(bào)童訣竅的簡(jiǎn)化版需求量x500510520概率p0.340.360.30問(wèn)題1：如何產(chǎn)生以下分布規(guī)律的隨機(jī)數(shù)據(jù)？注:rand(m,n)可以生成0,1上均勻分布隨機(jī)數(shù) 把 0,1分成長(zhǎng)度為0.34、0.36、0.30的三個(gè)區(qū)間 0,0.34、(0.34,0.70、(0.70,1 用rand(1,1)產(chǎn)生1個(gè)0,1上均勻分布隨機(jī)數(shù)，如該數(shù)在0,0.34、 (0.34,0.70或(0.70,1內(nèi)，相當(dāng)于該天的需求量相應(yīng)為500、

29、510和520 重復(fù)多次就可以若干天的需求量了生成n=20天的需求量的matlab代碼可以為報(bào)童訣竅的簡(jiǎn)化版生成n=20天的需求量的matlab代碼（定義函數(shù)）function y=randfun1(n)y=zeros(1,n); for i=1:n t=rand(1,1); if t=0.34 x=500; elseif tn) y=n*(a-b);else y=x*(a-b)-(n-x)*(b-c)end()()()()n abxnyx abnx bcxn模擬程序代碼n=1000;x=randfun1(n);y=0;for i=1:n y=y+fun2(x(i);endy/n報(bào)童訣竅的簡(jiǎn)化

30、版a=1;b=0.75;c=0.6;n=510;n=510; n=1000; y=0;for i=1:n t=rand(1,1); if (t=0.34) x=500; elseif (tn) y=y+n*(a-b); else y=y+x*(a-b)-(n-x)*(b-c); endendfprintf(1,平均利潤(rùn)=%.3f,y/n);也可完整模擬程序：根據(jù)隨機(jī)數(shù)t，計(jì)算需求量x值根據(jù)需求量x，計(jì)算利潤(rùn)并累加到y(tǒng)中。顯示平均利潤(rùn)售出價(jià)a 購(gòu)進(jìn)價(jià)b 退回價(jià)c購(gòu)進(jìn)數(shù)n 總利潤(rùn)y 模擬天數(shù)n報(bào)童訣竅的簡(jiǎn)化版 粒子分離器的某關(guān)鍵參數(shù)（記作y）由7個(gè)零件的參數(shù)（記作x1，x2，x7）決定，12345

31、670 1.5160.75xxxxxxxbbbccbb標(biāo) 定 值.容 差 等 級(jí)320.561.1644220.85315216712.62 10.36()()xy=174.42()()xxxxxxxxx x均值=標(biāo)定值 標(biāo)準(zhǔn)差=容差等級(jí)*標(biāo)定值3 關(guān)鍵參數(shù)y的目標(biāo)值是1.50，當(dāng)偏離為0.1但未達(dá)到0.3時(shí)，產(chǎn)品為次品，損失為1000元；當(dāng)偏離達(dá)到0.3時(shí)，產(chǎn)品為廢品，損失為9000元。由于工藝原因，7個(gè)零件參數(shù)可以看著是正態(tài)隨機(jī)變量，在后面的標(biāo)定值及容差等級(jí)情況下，求產(chǎn)品的平均損失？容差等級(jí)a=1% b=5%c=15%零件參數(shù)設(shè)計(jì)（1997a）例312345670 20

32、.160.75xxxxxxxbbbccbb標(biāo) 定 值.容 差 等 級(jí)均值=標(biāo)定值 標(biāo)準(zhǔn)差=容差等級(jí)*標(biāo)定值3容差等級(jí)a=1% b=5%c=15%2(,)iiixn ab零件損失q是一個(gè)隨機(jī)變量，平均損失就是期望e(q)()0(0)1000(1000)9000(9000)e qp qp qp q1000(0 1|1.5 |0.3)9000(|1.5 |0.3)pypy. 由于y的表達(dá)式很復(fù)雜，要想計(jì)算y的分布和上述概率很困難，我們必須尋找較為有效的近似方法。模擬：通過(guò)產(chǎn)生指定分布的隨機(jī)數(shù)，來(lái)代表7個(gè)零件的參數(shù)值，計(jì)算y值，確定損失大小。多做幾次可得均值零件參數(shù)%編寫(xiě)計(jì)算y值的

33、函數(shù)：y=funy(x1,x2,x3,x4,x5,x6,x7)%編寫(xiě)計(jì)算損失費(fèi)函數(shù):q=funq(y) d1=0.1,0.3,0.1,0.1,1.5,16,0.75 %標(biāo)定值d2=5,5,5,15,15,5,5./100 %容差等級(jí)a=d1b=d1.*d2./3x1=normrnd(a(1),b(1),1,1);x7=normrnd(a(7),b(7),1,1)y=funy(x1,x2,x3,x4,x5,x6,x7)q=funq(y）產(chǎn)生一組零件參數(shù)，相當(dāng)制作一個(gè)產(chǎn)品。重復(fù)n1000次（即生產(chǎn)n個(gè)產(chǎn)品），求出損失費(fèi)用的平均.(代碼略)按指定分布產(chǎn)生隨機(jī)數(shù)，作為零件參數(shù)計(jì)算該產(chǎn)品的關(guān)鍵參數(shù)y值和

34、損失的大小零件期望a零件標(biāo)準(zhǔn)差b零件參數(shù)用蒙特卡洛法解非線(xiàn)性規(guī)劃問(wèn)題用蒙特卡洛法解非線(xiàn)性規(guī)劃問(wèn)題對(duì)于非非線(xiàn)線(xiàn)性性規(guī)規(guī)劃劃問(wèn)問(wèn)題題： min f(x) xne s.t. 0)(xgi i=1,2, m jjjbxa j=1,2, n用蒙特卡洛法求求解解的的基基本本思思想想是：在估計(jì)的區(qū)域 (x1,x2,xn)| xj,jjba,j=1,2,n內(nèi)隨機(jī)取若干實(shí)驗(yàn)點(diǎn)，然后從試驗(yàn)點(diǎn)中找出可行點(diǎn)，再?gòu)目尚悬c(diǎn)中選擇最小點(diǎn)基本假設(shè)基本假設(shè) 試驗(yàn)點(diǎn)的第j個(gè)分量xj服從aj ,bj內(nèi)的均勻分布符號(hào)假設(shè)符號(hào)假設(shè) p: 試驗(yàn)點(diǎn)總數(shù)；maxp：最大試驗(yàn)點(diǎn)總數(shù)；k: 可行點(diǎn)總數(shù)； maxk：最大可行點(diǎn)數(shù);x*：迭代產(chǎn)生的

35、最優(yōu)點(diǎn)； q：迭代產(chǎn)生的最小值 f(x*)，其初始值為計(jì)算機(jī)所能表示的最大數(shù)求解過(guò)程求解過(guò)程 先產(chǎn)生一個(gè)隨機(jī)數(shù)作為初始試驗(yàn)點(diǎn),以后則將上一個(gè)試驗(yàn)點(diǎn)的第j個(gè)分量隨機(jī)產(chǎn)生，其它分量不變而產(chǎn)生一新的試驗(yàn)點(diǎn)這樣，每產(chǎn)生一個(gè)新試驗(yàn)點(diǎn)只需一個(gè)新的隨機(jī)數(shù)分量當(dāng)kmaxk或pmaxp時(shí)停止迭代框框 圖圖初始化:給定maxk,maxp;k=0,p=0,q:大整數(shù)xj=aj+r(bj-aj) j=1,2,nj=0j=j+1,p=p+1pmaxp?ynxj=aj+r(bj-aj)gi(x)0?i=1,2nynjmaxk?yn輸出x,q,停止yn 例例 max 21212221382xxxxxxz s.t 10321

36、 xx 01x 02x 在matlab軟件包中編程，共需三個(gè)文件:randlp.m, mylp.m, lpconst.m.主程序?yàn)閞andlp.m.% mylp.m% mylp.mfunction z=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)化為求最小值問(wèn)題% % l lp pc co on ns st t. .m mfunction lpc=lpconst(x)%約束條件if 3*x(1)+x(2)-10=-0.5%約束條件的誤差為5 . 0 lpc=1;else lpc=0;end% randlp.m% randlp.

37、m function sol,r1,r2=randlp(a,b,n) %隨機(jī)模擬解非線(xiàn)性規(guī)劃debug=1;a=0; %試驗(yàn)點(diǎn)下界b=10; %試驗(yàn)點(diǎn)上界n=1000; %試驗(yàn)點(diǎn)個(gè)數(shù)r1=unifrnd(a,b,n,1);%n1階的a,b均勻分布隨機(jī)數(shù)矩陣r2=unifrnd(a,b,n,1);sol=r1(1) r2(1);z0=inf;for i=1:n x1=r1(i); x2=r2(i); lpc=lpconst(x1 x2)； if lpc=1 z=mylp(x1 x2)； if zz0 z0=z； sol=x1 x2； end endendto matlab(randlp)返回隨機(jī)

38、數(shù)的產(chǎn)生1、均勻隨機(jī)數(shù)(均勻分布u0,1) rand()2、產(chǎn)生其他隨機(jī)數(shù)的方法 逆變換法、舍選法、近似抽樣法等。 設(shè)概率分布函數(shù)設(shè)概率分布函數(shù)f(x)是嚴(yán)格單調(diào)增的，是嚴(yán)格單調(diào)增的，f的反函數(shù)記作的反函數(shù)記作f-1。先產(chǎn)生。先產(chǎn)生uu(0,1),再取再取xf-1即為所求。即為所求。如指數(shù)分布，分布函數(shù)如指數(shù)分布，分布函數(shù)f(x)=1-exp(-),可得可得)1ln()(1uufx3、(非)常見(jiàn)分布隨機(jī)數(shù)如何產(chǎn)生？（離散）經(jīng)驗(yàn)分布函數(shù)法: 設(shè)一組實(shí)際數(shù)據(jù)，將它們分組整理形成頻率圖或表格。x0xf 3.5 5.5 8 0.3 0.5 0.2 3.5 3.5, 5.5） 5.5, 8) 8構(gòu)造經(jīng)驗(yàn)

39、分布函數(shù)xf 0 0.3 0.8 1由均勻隨機(jī)數(shù)r0,1，決定x的抽樣值。 當(dāng)0r0.3，抽樣值xi（0, 3.5； xi=3.5 當(dāng)0.3r0.8，抽樣值xi（3.5, 5.5； xi=5.5 當(dāng)r0.8，抽樣值xi（5.5, 8； xi=8r=rand;if 0r & r0.3 y=3.5;elseif 0.3=r & r0.8 y=5.5;else y=8;endy % 產(chǎn)生一個(gè)隨機(jī)數(shù)產(chǎn)生一個(gè)隨機(jī)數(shù)matlab程序：r=rand(10); lisanrnd1.my=;for i=1:10 if 0r(i)&r(i)0.3 y(i)=3.5; elseif 0.3=

40、r(i)&r(i)0為常數(shù)，則稱(chēng)x服從參數(shù)為 的指數(shù)分布指數(shù)分布。000)(xxexft指數(shù)分布的期望值為 1排隊(duì)服務(wù)系統(tǒng)中顧客到達(dá)率為常數(shù)時(shí)的到達(dá)間隔、故障率為常數(shù)時(shí)零件的壽命都服從指數(shù)分布。指數(shù)分布在排隊(duì)論、可靠性分析中有廣泛應(yīng)用。注意：注意：matlab中，產(chǎn)生參數(shù)為 的指數(shù)分布的命令為exprnd( )1例例 顧客到達(dá)某商店的間隔時(shí)間服從參數(shù)為顧客到達(dá)某商店的間隔時(shí)間服從參數(shù)為0.10.1的指數(shù)分布的指數(shù)分布 指數(shù)分布的均值為指數(shù)分布的均值為1/0.1=101/0.1=10。 指兩個(gè)顧客到達(dá)商店的平均間隔時(shí)間是指兩個(gè)顧客到達(dá)商店的平均間隔時(shí)間是1010個(gè)單位時(shí)間個(gè)單位時(shí)間. .

41、即平均即平均1010個(gè)單位個(gè)單位時(shí)間到達(dá)時(shí)間到達(dá)1 1個(gè)顧客個(gè)顧客. . 顧客到達(dá)的間隔時(shí)間可用顧客到達(dá)的間隔時(shí)間可用exprnd(10)exprnd(10)模擬。模擬。設(shè)離散型隨機(jī)變量x的所有可能取值為0,1,2,且取各個(gè)值的概率為其中 0為常數(shù)，則稱(chēng)x服從參數(shù)為 的帕松分布帕松分布。, 2 , 1 , 0,!)(kkekxpk5產(chǎn)生 mn 階參數(shù)為的帕松分布的隨機(jī)數(shù)矩陣：poissrnd poissrnd (,m, n)帕松分布在排隊(duì)系統(tǒng)、產(chǎn)品檢驗(yàn)、天文、物理等領(lǐng)域有廣泛應(yīng)用。帕松分布的期望值為如相繼兩個(gè)事件出現(xiàn)的間隔時(shí)間服從參數(shù)為 的指數(shù)分布，則在單位時(shí)間間隔內(nèi)事件出現(xiàn)的次數(shù)服從參數(shù)為

42、的泊松分布即單位時(shí)間內(nèi)該事件出現(xiàn)k次的概率為：, 2 , 1 , 0,!)(kkekxpk反之亦然。指數(shù)分布與帕松分布的關(guān)系： (1) (1)指兩個(gè)顧客到達(dá)商店的平均間隔時(shí)間是指兩個(gè)顧客到達(dá)商店的平均間隔時(shí)間是1010個(gè)單位時(shí)間個(gè)單位時(shí)間. .即平均即平均1010個(gè)單個(gè)單位時(shí)間到達(dá)位時(shí)間到達(dá)1 1個(gè)顧客個(gè)顧客. . (2) (2)指一個(gè)單位時(shí)間內(nèi)平均到達(dá)指一個(gè)單位時(shí)間內(nèi)平均到達(dá)0.10.1個(gè)顧客個(gè)顧客例例 (1)(1)顧客到達(dá)某商店的間隔時(shí)間服從參數(shù)為顧客到達(dá)某商店的間隔時(shí)間服從參數(shù)為0.10.1的指數(shù)分布的指數(shù)分布 (2)(2)該商店在單位時(shí)間內(nèi)到達(dá)的顧客數(shù)服從參數(shù)為該商店在單位時(shí)間內(nèi)到達(dá)的

43、顧客數(shù)服從參數(shù)為0.10.1的帕松分布的帕松分布 排隊(duì)問(wèn)題 機(jī)械故障等候修理 飛機(jī)跑道日常生活中經(jīng)常遇到的排隊(duì)問(wèn)題： 自選商場(chǎng)收款臺(tái) 醫(yī)院里病人等候就診輸入情況：輸入情況： 顧客到達(dá)時(shí)間和服務(wù)時(shí)間。顧客到達(dá)時(shí)間和服務(wù)時(shí)間。系統(tǒng)狀態(tài)：系統(tǒng)狀態(tài)： 排隊(duì)等候的顧客數(shù)目（隊(duì)長(zhǎng)）排隊(duì)等候的顧客數(shù)目（隊(duì)長(zhǎng)）l(t)l(t) 服務(wù)員是否在工作或服務(wù)效率等；服務(wù)員是否在工作或服務(wù)效率等；簡(jiǎn)圖：第二顧客接受服務(wù)時(shí)間s2x50 x1x2x3x4 y1 y2 y3 y4 y5d2,.2 , 1,isdxyiiii關(guān)系：系統(tǒng)在什么條件下處于空閑狀態(tài)？(yixi+1) 排隊(duì)系統(tǒng)中，顧客到達(dá)時(shí)刻數(shù)據(jù)如何收集？對(duì)每個(gè)顧客的

44、服務(wù)時(shí)間如何？x： x1，x2，xn第一個(gè)顧客到達(dá)的時(shí)刻第二個(gè)顧客到達(dá)的時(shí)刻計(jì)算機(jī)遵循某種規(guī)則進(jìn)行隨機(jī)抽樣。s: s1,x2,xn排隊(duì)服務(wù)系統(tǒng)的模擬1、模擬目的2、系統(tǒng)假設(shè)與輸入3、系統(tǒng)的狀態(tài)4、初始條件和終止條件5、模擬過(guò)程的運(yùn)行6、系統(tǒng)的性能指標(biāo)7、與分析模型的比較8、其它排隊(duì)系統(tǒng)1、模擬目的如果知道了顧客到達(dá)的時(shí)間和服務(wù)時(shí)間的統(tǒng)計(jì)規(guī)律（大量實(shí)際數(shù)據(jù)或概率分布），怎樣通過(guò)模擬得到衡量系統(tǒng)的性能指標(biāo)，并與分析模型的結(jié)果進(jìn)行比較。2、系統(tǒng)假設(shè)與輸入首行作如下基本假設(shè)：（1）顧客源是無(wú)窮的（2）排隊(duì)的長(zhǎng)度沒(méi)有限制（3）到達(dá)系統(tǒng)的顧客按先后順序進(jìn)入服務(wù)假設(shè)由概率分布產(chǎn)生了如下的數(shù)據(jù)：k123456

45、789ik013413182sk231341341其中ik表示第k位顧客與第k-1位顧客到達(dá)的時(shí)間間隔，sk是第k位顧客的服務(wù)時(shí)間，下面就以這組數(shù)據(jù)作為系統(tǒng)的輸入。3、系統(tǒng)的狀態(tài)n隊(duì)長(zhǎng)：排隊(duì)等候的顧客數(shù)目。記作l(t).n服務(wù)員的狀態(tài)用s(t)表示，當(dāng)服務(wù)員工作時(shí)，令s(t)=1，空閑時(shí)令s(t)=0.n事件：引起系統(tǒng)狀態(tài)l(t)和s(t)改變的行為。 “顧客到達(dá)”和“顧客離開(kāi)”設(shè)第k位顧客到達(dá)和離開(kāi)時(shí)刻分別為ak和dk，則有, 2 , 1,),max( 11ksdadiaakkkkkkk4、初始條件和終止條件不妨假設(shè)t=0是第1位顧客到達(dá)的時(shí)刻，此時(shí)服務(wù)員處于空閑狀態(tài)。模擬的終止時(shí)刻是t（給

46、定值），或模擬直到第n個(gè)（給定值）顧客進(jìn)入服務(wù)的時(shí)間。5、模擬過(guò)程 進(jìn)行模擬時(shí)，可以事先產(chǎn)生一批數(shù)據(jù)ik,sk,再計(jì)算ak,dk，然后讓時(shí)鐘t按ak,dk從小到大的順序跳動(dòng)；也可在時(shí)鐘t跳到某一事件發(fā)生時(shí)，才產(chǎn)生一個(gè)所需要的i或s。這里采用后一種方法。 設(shè)當(dāng)前時(shí)鐘為t，隊(duì)長(zhǎng)l(t)記作wl，服務(wù)員狀態(tài)s(t)記作ss，t以后下一個(gè)“顧客到達(dá)”事件的發(fā)生時(shí)刻記作at, t以后下一個(gè)“顧客離開(kāi)”事件的發(fā)生時(shí)刻記作dt.at0wl=wl+1令ss1產(chǎn)生s置dt=t+s產(chǎn)生i置att+it=t令ss0wlwl1產(chǎn)生s置dtt+s停止置dt999是否是否是否否是nst=6;it=5;wl=0;ss=0;

47、at=0;dt=99999;t=10;t=0;ttt=0;www=0;sss=0;nwhile (t=t)n if at0n wl=wl-1; s=exprnd(1/st,1,1); dt=t+s;n elsen ss=0; dt=99999;n endn endn ttt=ttt,t; www=www,wl; sss=sss,ss;nendn ttt(2:length(ttt)n www(2:length(ttt)n sss(2:length(ttt)6、系統(tǒng)的性能指標(biāo)n平均隊(duì)長(zhǎng)n平均待時(shí)間n服務(wù)利用率dttltlt0)(1tnltdtlnwt)()(10)()(10tdtstut7、與分析

48、模型的比較對(duì)于一般的排隊(duì)服務(wù)系統(tǒng)，很難建立分析模型并求出結(jié)果，即使是單服務(wù)員系統(tǒng)，也只有在顧客到達(dá)間隔和服務(wù)時(shí)間服從某種特殊的概率分布時(shí)，才可以在穩(wěn)態(tài)情況下由分析模型得到上述三個(gè)指標(biāo)的解析形式的結(jié)果。分析模型基本假設(shè)：1、到達(dá)間隔服從參數(shù)的指數(shù)分布；2、服務(wù)時(shí)間服從參數(shù)為的指數(shù)分布；3、排隊(duì)規(guī)則是“先到先服務(wù)”，顧客數(shù)目和排隊(duì)長(zhǎng)度都無(wú)限制。 排隊(duì)論中將這樣的單服務(wù)員系統(tǒng)記作m/m/1。服務(wù)強(qiáng)度：穩(wěn)態(tài)（1）穩(wěn)態(tài)平均隊(duì)長(zhǎng)）( 12l穩(wěn)態(tài)平均等待時(shí)間）( l8、其它的排隊(duì)服務(wù)系統(tǒng)n有兩個(gè)（或多個(gè)）服務(wù)員n串聯(lián)服務(wù)系統(tǒng)，每個(gè)顧客接受多次服務(wù)n從排隊(duì)顧客中選服務(wù)時(shí)間最短的最先服務(wù)n隊(duì)的長(zhǎng)度有限制n帶有優(yōu)

49、化目標(biāo)的排隊(duì)服務(wù)系統(tǒng) 庫(kù)存太多造成浪費(fèi)或資金積壓，庫(kù)存太少不能滿(mǎn)足需求也造成損失，需要進(jìn)行決策：何時(shí)進(jìn)貨，進(jìn)多少，使得平均費(fèi)用最少，而收益最大。庫(kù)存系統(tǒng)庫(kù)存系統(tǒng)方案甲：按前一天的銷(xiāo)售量作為當(dāng)天的生產(chǎn)量；方案乙：按前二天的平均銷(xiāo)售量作為當(dāng)天的生產(chǎn)量； 假定市場(chǎng)對(duì)該產(chǎn)品的每天需要量需要量是一個(gè)隨機(jī)變量，但從以往的統(tǒng)計(jì)分析得知它服從正態(tài)分布n(135,22.42）.例例1 某企業(yè)生產(chǎn)易變質(zhì)的產(chǎn)品。（如蛋糕）當(dāng)天生產(chǎn)的產(chǎn)品必須當(dāng)天售出，否則就會(huì)變質(zhì)。該產(chǎn)品單位成本為2.5元，單位產(chǎn)品售價(jià)為5元。企業(yè)為避免存貨過(guò)多而造成損失，擬從以下兩種貨存方案中選出一個(gè)較優(yōu)的方案.(如何決定當(dāng)天的生產(chǎn)數(shù)量？)模擬的基

50、本思路： 需要獲得市場(chǎng)對(duì)該產(chǎn)品的需要量的樣本 值；(產(chǎn)生抽樣值) 尋找貨存量與需要量之間的關(guān)系； 按照兩種不同方案（貨存量需要量）計(jì)算出經(jīng) t 天后企業(yè)的利潤(rùn)值(累計(jì))； 比較大小，從中選出一個(gè)較優(yōu)的方案；方案甲：貨存量方案甲：貨存量y1分析兩種狀態(tài)下分析兩種狀態(tài)下xy1和和xy1的利潤(rùn)值的利潤(rùn)值l1置初始狀態(tài)置初始狀態(tài)獲取產(chǎn)品的需求量獲取產(chǎn)品的需求量xx= normrnd（135，22.5）n = 1方案乙：貨存量方案乙：貨存量y2分析兩種狀態(tài)下分析兩種狀態(tài)下xy2和和xy2的利潤(rùn)值的利潤(rùn)值l2累計(jì)利潤(rùn)值累計(jì)利潤(rùn)值tl1累計(jì)利潤(rùn)值累計(jì)利潤(rùn)值tl2n=n+1判斷：判斷：nt ( t天天)比較計(jì)

51、算：比較計(jì)算：maxtl1，tl2輸出最大利潤(rùn)值及方案輸出最大利潤(rùn)值及方案matlab編制的程序： kucun.m。運(yùn)行程序：kutest.m初始狀態(tài)：t=100;s1=136;s21=126;s22=148; tl1,tl2=kucun(t,s1,s21,s22)觀(guān)察計(jì)算結(jié)果 如果多數(shù)情況有：tl1tl2, 則認(rèn)為方案乙較優(yōu).思考：1）一般模擬結(jié)果波動(dòng)性較大，如何減少這種波動(dòng)？ 2）修改上述程序?？煽啃詥?wèn)題可靠性問(wèn)題 一設(shè)備上有三個(gè)相同的軸承，每個(gè)軸承正常工作壽命為隨機(jī)變量，概率分布如下：壽命/h1000110012001300140015001600170018001900概率0.100.

52、90.120.020.060.050.05有軸承損壞設(shè)備停止工作檢修工準(zhǔn)備開(kāi)始更換部件，稱(chēng)為一個(gè)延遲時(shí)間，它也是隨機(jī)變量，分布如下：延遲時(shí)間 /min51015概 率主要費(fèi)用：1、設(shè)備停工損失費(fèi)： 5 元/分鐘；2、檢修工人的工時(shí)費(fèi)：12 元/小時(shí)；3、軸承的成本費(fèi)： 16元/個(gè)更換軸承所需要的時(shí)間： 一個(gè) 兩個(gè) 三個(gè) 20 30 40 （min）問(wèn)題：現(xiàn)在有兩種方案：方案一：損壞一個(gè)軸承只更換一個(gè)軸承；方案二：一旦有軸承損壞就全部更換；試通過(guò)計(jì)算機(jī)模擬對(duì)以上兩種方案做出評(píng)價(jià)。 隨機(jī)數(shù)怎樣產(chǎn)生？ 模擬時(shí)選用時(shí)間步長(zhǎng)法還是 事件步長(zhǎng)法？關(guān)于隨機(jī)數(shù)的產(chǎn)生

53、見(jiàn) lrnd.m（零件壽命） yrnd.m（延遲時(shí)間） 方案一的數(shù)學(xué)模型：kekao1.m目標(biāo)函數(shù) min c =ui / t（三個(gè)軸承同時(shí)損壞）（兩個(gè)軸承同時(shí)損壞）只有一個(gè)軸承損壞）3166040125)40(2166030125)30(166020125)20(iiiitttu其中 損失費(fèi)工時(shí)費(fèi)成本費(fèi)ti 表示延遲時(shí)間 設(shè)事件發(fā)生時(shí)刻 ti，它是由軸承工作壽命li、延遲時(shí)間和修理時(shí)間迭加產(chǎn)生?？偟倪\(yùn)行時(shí)間為t（10000）。例： （手工模擬）事件類(lèi)型發(fā)生時(shí)刻延遲時(shí)間a1400 h5 minb1500 h15 minc1500 h15 min方案一的初始事件表 t=0t=1400； cost= (5+20)5+12 1/3+16=145； 產(chǎn)生下一個(gè)a事件發(fā)生時(shí)刻1400+25/60+1000（li）=2400小時(shí)25分鐘損失費(fèi)工時(shí)費(fèi)成本費(fèi)6020事件類(lèi)型發(fā)生時(shí)刻延遲時(shí)間b1500 h15 minc1500 h15 mina2400 h25min5 min第一次刷新后的事件表 t=1400 t=1500 ； cost=145+ (15+30)5+12 1/2+16 2 = 408 ； 產(chǎn)生下一個(gè)b事件發(fā)生時(shí)刻15