數(shù)模競(jìng)賽培訓(xùn)之計(jì)算機(jī)模擬

數(shù)學(xué)建模培訓(xùn)計(jì)算機(jī)模擬主講：何仁斌概述

計(jì)算機(jī)模擬是用計(jì)算機(jī)對(duì)系統(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ī)程序直接建立真實(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)模擬和離散系統(tǒng)模擬。

狀態(tài)隨時(shí)間連續(xù)變化的系統(tǒng)稱為連續(xù)系統(tǒng)。通常該系統(tǒng)的模型一般可以用微分方程的形式表達(dá)，通過(guò)一些物理機(jī)理推導(dǎo)出來(lái)。模擬結(jié)果往往是近似的。例如，追逐問(wèn)題，濃度問(wèn)題。一、連續(xù)系統(tǒng)(x1(t),…,xn(t))1、追逐問(wèn)題A→BB→CC→DD→A狀態(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)系；2.以時(shí)間間隔△t進(jìn)行采樣，并且計(jì)算各個(gè)時(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+△t),y1(t+△t))（由幾何原理）≈(x1(t)+v△tcos(θ),y1(t)+v△tsin(θ)θx1,y1x2,y2在t時(shí)刻下3.選取足夠小的△t，模擬到任意兩人的距離小于v△t為止。4.連接四人在各時(shí)刻下的位置，就得到所求的運(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=[0001010

10

100];x(9)=x(1);x(10)=x(2);holdonaxis('equal');axis([010010]);while(d>0.1)

fori=1:2:7d=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);endholdonABCD

某軍的一導(dǎo)彈基地發(fā)現(xiàn)正北方向120km處海面上有敵艇一艘以90km/h的速度向正東方向行駛。該基地立即發(fā)射導(dǎo)彈跟蹤追擊敵艇，導(dǎo)彈速度為450km/h，自動(dòng)導(dǎo)航系統(tǒng)使導(dǎo)彈在任一時(shí)刻都能對(duì)準(zhǔn)敵艇。試問(wèn)導(dǎo)彈在何時(shí)何處擊中敵艇？2、導(dǎo)彈跟蹤問(wèn)題P(x(t),y(t))OxM(vet,H)yA(0,H)正東方向正北方向H=120（千米）km/hkm/h敵導(dǎo)(xk(t),yk(t))simu1.mstest1.mhuatu.m

微分方程建模計(jì)算機(jī)模擬α消去t：微分方程建模1）2）>>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-120)^(4/5)-1800*(-1)^(4/5)+(-1)^(3/5)*202500^(1/5)*(y-120)^(6/5))/(72*(-(-1)^(3/5))^(1/2))>>y=120>>eval(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(90×2τ,120)；function[num,y_j,L,T]=simu1(x,y,t,eps)k=0;whilek<1000(simu1.m)p=90*k*t-x;q=120-y;d1=p/(p^2+q^2)^0.5;d2=q/(p^2+q^2)^0.5;x=x+450*t*d1;y=y+450*t*d2;

if(abs(q)<eps)

break;

end

k=k+1;endnum=k;L=x;y_j=y;T=L/90;(stest1.m)t=0.1;x=0;y=0;[num,y_j,L,T]=simu1(x,y,t,eps)計(jì)算結(jié)果：Num=2;y_j=131.2155;L=22.675（km）;T=0.2519（h）導(dǎo)彈的位置(x(t),y(t))能否用圖形描述它們的追逐運(yùn)動(dòng)過(guò)程？敵艇的位置(水平方向(90*k*t,120)x=[0,0];(huatu.m)axis('equal');axis([0,140,0,140]);gridonholdont=0.001;%步長(zhǎng)while(abs(x(2)-120)>0.1)%終止條件(<)

fork=1:280p=90*k*t-x(1);q=120-x(2);d1=p/(p^2+q^2)^0.5;d2=q/(p^2+q^2)^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));

endendholdon（2）如果當(dāng)基地發(fā)射導(dǎo)彈的同時(shí)，敵艇立即由儀器發(fā)覺(jué)。假定敵艇為高速快艇，它即刻以135千米/小時(shí)的速度與導(dǎo)彈方向成垂直的方向逃逸，問(wèn)導(dǎo)彈何時(shí)何地?fù)糁袛惩?？不停地改變逃跑策略，運(yùn)動(dòng)軌跡如何？

修正敵艇速度：90→135千米/小時(shí)思考：二、離散系統(tǒng)

離散系統(tǒng)是指系統(tǒng)狀態(tài)只在有限的時(shí)間點(diǎn)或可數(shù)的時(shí)間點(diǎn)上發(fā)生變化的系統(tǒng)。并假設(shè)離散系統(tǒng)狀態(tài)的變化是在一個(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)期觀察發(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)的概率和期望值.但這樣只能給出作戰(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)題分析[2]當(dāng)指示正確時(shí)，我方火力單位的射擊結(jié)果情況[1]觀察所對(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)．

這時(shí)可用投擲骰子的方法來(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è)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+1i＜20?E=(k2+k3)/20E1=0*k1/20+1*k2/20+2*k3/20停止硬幣正面?YNNY1，2，34，564.模擬結(jié)果5.理論計(jì)算6.結(jié)果比較

雖然模擬結(jié)果與理論計(jì)算不完全一致，但它卻能更加真實(shí)地表達(dá)實(shí)際戰(zhàn)斗動(dòng)態(tài)過(guò)程．蒙特卡羅方法(MonteCarlomethod)

蒙特卡羅方法(MonteCarlomethod)，或稱計(jì)算機(jī)隨機(jī)模擬方法，是一種基于“隨機(jī)數(shù)”的計(jì)算方法。源于美國(guó)在第二次世界大戰(zhàn)研制原子彈的“曼哈頓計(jì)劃”，該計(jì)劃的主持人之一數(shù)學(xué)家馮·諾伊曼用馳名世界的賭城—摩納哥的MonteCarlo—來(lái)命名這種方法，為它蒙上了一層神秘色彩。蒙特卡羅方法的基本思想很早以前就被人們所發(fā)現(xiàn)和利用。早在17世紀(jì)，人們就知道用事件發(fā)生的“頻率”來(lái)決定事件的“概率”。19世紀(jì)人們用蒲豐投針的方法來(lái)計(jì)算圓周率π，上世紀(jì)40年代電子計(jì)算機(jī)的出現(xiàn)，特別是近年來(lái)高速電子計(jì)算機(jī)的出現(xiàn)，使得用數(shù)學(xué)方法在計(jì)算機(jī)上大量、快速地模擬這樣的試驗(yàn)成為可能。蒲豐投針實(shí)驗(yàn)近似計(jì)算圓周率π蒲豐投針實(shí)驗(yàn)：法國(guó)科學(xué)家蒲豐(Buffon)在1777年提出的蒲豐投針實(shí)驗(yàn)是早期幾何概率一個(gè)非常著名的例子。蒲豐投針實(shí)驗(yàn)的重要性并非是為了求得比其它方法更精確的π值，而是它開(kāi)創(chuàng)了使用隨機(jī)數(shù)處理確定性數(shù)學(xué)問(wèn)題的先河，是用偶然性方法去解決確定性計(jì)算的前導(dǎo)，由此可以領(lǐng)略到從“概率土壤”上開(kāi)出的一朵瑰麗的鮮花——蒙特卡羅方法(MC)蒲豐投針實(shí)驗(yàn)可歸結(jié)為下面的數(shù)學(xué)問(wèn)題：平面上畫(huà)有距離為a的一些平行線，向平面上任意投一根長(zhǎng)為l(l<a)的針，假設(shè)針落在任意位置的可能性相同，試求針與平行線相交的概率P(從而求π)蒲豐投針實(shí)驗(yàn)近似計(jì)算圓周率π蒲豐投針實(shí)驗(yàn)：如右圖所示，以M表示針落下后的中點(diǎn)，以x表示M到最近一條平行線的距離，以φ表示針與此線的交角：針落地的所有可能結(jié)果滿足：其樣本空間視作矩形區(qū)域Ω,面積是:針與平行線相交的條件：它是樣本空間Ω子集A，面積是：syms

lphi;int('l/2*sin(phi)',phi,0,pi);%ans=l因此，針與平行線相交的概率為：從而有：特別當(dāng)時(shí)蒲豐投針實(shí)驗(yàn)近似計(jì)算圓周率π蒲豐投針實(shí)驗(yàn)近似計(jì)算圓周率π蒲豐投針實(shí)驗(yàn)的計(jì)算機(jī)模擬：formatlong;%設(shè)置15位顯示精度a=1;

l=0.6;

%兩平行線間的寬度和針長(zhǎng)figure;axis([0,pi,0,a/2]);%初始化繪圖板set(gca,'nextplot','add');%初始化繪圖方式為疊加counter=0;

n=2010;

%初始化計(jì)數(shù)器和設(shè)定投針次數(shù)x=unifrnd(0,a/2,1,n);phi=unifrnd(0,pi,1,n);

%樣本空間Ωfori=1:nifx(i)<l*sin(phi(i))/2

%滿足此條件表示針與線的相交

plot(phi(i),x(i),'r.');frame(i)=getframe;%描點(diǎn)并取幀

counter=counter+1;%統(tǒng)計(jì)針與線相交的次數(shù)endendfren=counter/n;

pihat=2*l/(a*fren)

%用頻率近似計(jì)算π%movie(frame,1)%播放幀動(dòng)畫(huà)1次蒲豐投針實(shí)驗(yàn)近似計(jì)算圓周率π蒲豐投針實(shí)驗(yàn)的計(jì)算機(jī)模擬：意大利數(shù)學(xué)家拉澤里尼得到了準(zhǔn)確到6位小數(shù)的π值，不過(guò)他的實(shí)驗(yàn)因?yàn)樘珳?zhǔn)確而受到了質(zhì)疑蒲豐投針實(shí)驗(yàn)計(jì)算圓周率π蒙特卡羅投點(diǎn)法是蒲豐投針實(shí)驗(yàn)的推廣：在一個(gè)邊長(zhǎng)為a的正方形內(nèi)隨機(jī)投點(diǎn)，該點(diǎn)落在此正方形的內(nèi)切圓中的概率應(yīng)為該內(nèi)切圓與正方形的面積比值，即n=10000;a=2;m=0;fori=1:nx=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,n);y=rand(1,n);num2str(length(find((x-1).^2+(y-1).^2<=1))*4/n)xyo(a/2,a/2)一些MonteCarlo思想的例子[例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ì)相遇的比例就可近似為相遇的概率。這也是MonteCarlo思想.兩船到達(dá)碼頭時(shí)刻服從[0,24]上的均勻分布，甲船停留2小時(shí),乙船停留1小時(shí),相遇概率？方法一利用軟件產(chǎn)生一組[0，24]上的隨機(jī)數(shù)xi:

10.45.30.39.020.16.87.8……5.4

分別代表甲船到達(dá)碼頭時(shí)間；再產(chǎn)生一組[0，24]上的隨機(jī)數(shù)yi:

6.03.48.117.50.813.4……14.0

分別代表乙船到達(dá)碼頭時(shí)間；把（xi，yi）(i＝1，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=1000;c=0;

%模擬次數(shù)N，相遇次數(shù)c清零fori=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元.每天需求量X是離散型隨機(jī)變量，其分布為[例2]X500510520P0.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.340.360.30按照需求量的分布規(guī)律，隨機(jī)生成N=20個(gè)數(shù)據(jù):510520500510520500500

500

500510500500

500520510520510510代表20天的需求量，計(jì)算出報(bào)童在這20天的總利潤(rùn)和平均利潤(rùn)，用平均利潤(rùn)來(lái)近似報(bào)童的平均收入。這也是MonteCarlo方法.問(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、510和520重復(fù)多次就可以若干天的需求量了生成N=20天的需求量的matlab代碼可以為報(bào)童訣竅的簡(jiǎn)化版生成N=20天的需求量的matlab代碼（定義函數(shù)）functiony=randfun1(N)y=zeros(1,N);fori=1:Nt=rand(1,1);

ift<=0.34

X=500;

elseift<=0.70

X=510;

else

X=520;

end

y(i)=X;end根據(jù)隨機(jī)數(shù)t的范圍，確定需求量X值，并保存到數(shù)組的相應(yīng)位置中（關(guān)鍵部分）產(chǎn)生1行N列的全0矩陣，目的分配好矩陣大小，可以省略文件名為randfun1.m報(bào)童訣竅的簡(jiǎn)化版問(wèn)題2：如何從需求量計(jì)算利潤(rùn)？售出價(jià)a=1，購(gòu)進(jìn)價(jià)b=0.75，退回價(jià)c=0.6，購(gòu)進(jìn)數(shù)n=510函數(shù)文件名fun2.m

functiony=fun2(x)a=1;%售出價(jià)

b=0.75;%購(gòu)進(jìn)價(jià)

c=0.6;%退回價(jià)

n=510;%購(gòu)進(jìn)數(shù)量if(x>n)

y=n*(a-b);

else

y=x*(a-b)-(n-x)*(b-c)

end模擬程序代碼N=1000;x=randfun1(N);y=0;fori=1:Ny=y+fun2(x(i));endy/N報(bào)童訣竅的簡(jiǎn)化版a=1;b=0.75;c=0.6;n=510;

n=510;N=1000;y=0;fori=1:Nt=rand(1,1);

if(t<=0.34)x=500;

elseif(t<=0.70)x=510;

elsex=520;

end

if(x>n)y=y+n*(a-b);

elsey=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）決定，均值=標(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）[例3]均值=標(biāo)定值標(biāo)準(zhǔn)差=容差等級(jí)*標(biāo)定值÷3容差等級(jí)A=1%B=5%C=15%零件損失Q是一個(gè)隨機(jī)變量，平均損失就是期望E(Q)由于y的表達(dá)式很復(fù)雜，要想計(jì)算y的分布和上述概率很困難，我們必須尋找較為有效的近似方法。模擬：通過(guò)產(chǎn)生指定分布的隨機(jī)數(shù)，來(lái)代表7個(gè)零件的參數(shù)值，計(jì)算y值，確定損失大小。多做幾次可得均值零件參數(shù)%編寫(xiě)計(jì)算y值的函數(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=D1

b=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ù)N＝1000次（即生產(chǎn)N個(gè)產(chǎn)品），求出損失費(fèi)用的平均.(代碼略)按指定分布產(chǎn)生隨機(jī)數(shù)，作為零件參數(shù)計(jì)算該產(chǎn)品的關(guān)鍵參數(shù)y值和損失的大小零件期望a零件標(biāo)準(zhǔn)差b零件參數(shù)用蒙特卡洛法解非線性規(guī)劃問(wèn)題基本假設(shè)試驗(yàn)點(diǎn)的第j個(gè)分量xj服從[aj,bj]內(nèi)的均勻分布．符號(hào)假設(shè)求解過(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)K>MAXK或P>MAXP時(shí)停止迭代．框圖初始化:給定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,停止YN

在Matlab軟件包中編程，共需三個(gè)Ｍ－文件:randlp.m,

mylp.m,

lpconst.m.主程序?yàn)閞andlp.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)化為求最小值問(wèn)題%randlp.m

function[sol,r1,r2]=randlp(a,b,n)

%隨機(jī)模擬解非線性規(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);

%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])；

if

lpc==1z=mylp([x1x2])；

ifz<z0z0=z；

sol=[x1x2]；

end

endendToMatlab(randlp)返回隨機(jī)數(shù)的產(chǎn)生1、均勻隨機(jī)數(shù)(均勻分布U[0,1])rand()2、產(chǎn)生其他隨機(jī)數(shù)的方法逆變換法、舍選法、近似抽樣法等。

設(shè)概率分布函數(shù)F(x)是嚴(yán)格單調(diào)增的，F(xiàn)的反函數(shù)記作F-1。先產(chǎn)生U~U(0,1),再取X＝F-1即為所求。如指數(shù)分布，分布函數(shù)F(x)=1-exp(-λ),可得

3、(非)常見(jiàn)分布隨機(jī)數(shù)如何產(chǎn)生？（離散）經(jīng)驗(yàn)分布函數(shù)法:①設(shè)一組實(shí)際數(shù)據(jù)，將它們分組整理形成頻率圖或表格。x≥0Xf3.55.580.30.50.2

＜3.5[3.5,5.5）[5.5,8)≥8構(gòu)造經(jīng)驗(yàn)分布函數(shù)XF00.30.81②由均勻隨機(jī)數(shù)R∈[0,1]，決定X的抽樣值。當(dāng)0<R＜0.3，抽樣值xi∈（0,3.5]；xi=3.5

當(dāng)0.3≤R＜0.8，抽樣值xi∈（3.5,5.5]；xi=5.5

當(dāng)R≥0.8，抽樣值xi∈（5.5,8]；xi=8r=rand;if0<r&r<0.3y=3.5;elseif

0.3<=r&r<0.8y=5.5;else

y=8;endy%產(chǎn)生一個(gè)隨機(jī)數(shù)Matlab程序：r=rand(10);lisanrnd1.my=[];fori=1:10if0<r(i)&r(i)<0.3

y(i)=3.5;

elseif

0.3<=r(i)&r(i)<0.8

y(i)=5.5;else

y(i)=8;endendy%產(chǎn)生隨機(jī)變量數(shù)組。

計(jì)算結(jié)果y=Columns1through75.55.53.55.53.55.55.5Columns8through108.05.55.5MATLAB的實(shí)現(xiàn)常見(jiàn)分布的隨機(jī)變量獲取抽樣值1、U[0,1]R=rand(m,n)

2、U[a,b]R=unifrnd(a,b,m,n)3、E()

R=exprnd(1/

,m,n)4、N(0,1)R=randn(m,n)5、N(mu,sigma)R=normrnd(mu,sigma,m,n)6、B(N,p)R=binornd(N,p,m,n)7、P(

)R=poissrnd(

,m,n)以上語(yǔ)句都是產(chǎn)生m

n

的隨機(jī)矩陣.(hist(x))產(chǎn)生模擬隨機(jī)數(shù)的計(jì)算機(jī)命令2．產(chǎn)生m

n階［０，１］均勻分布的隨機(jī)數(shù)矩陣：rand(m,n)

產(chǎn)生一個(gè)［０，１］均勻分布的隨機(jī)數(shù)：rand1．產(chǎn)生m

n階［a，b］均勻分布U（a，b）的隨機(jī)數(shù)矩陣：

unifrnd(a,b,m,n)

產(chǎn)生一個(gè)［a，b］均勻分布的隨機(jī)數(shù)：unifrnd(a,b)

當(dāng)只知道一個(gè)隨機(jī)變量取值在（a，b）內(nèi)，但不知道（也沒(méi)理由假設(shè)）它在何處取值的概率大，在何處取值的概率小，就只好用U（a，b）來(lái)模擬它。當(dāng)研究對(duì)象視為大量相互獨(dú)立的隨機(jī)變量之和，且其中每一種變量對(duì)總和的影響都很小時(shí)，可以認(rèn)為該對(duì)象服從正態(tài)分布。機(jī)械加工得到的零件尺寸的偏差、射擊命中點(diǎn)與目標(biāo)的偏差、各種測(cè)量誤差、人的身高、體重等，都可近似看成服從正態(tài)分布。若連續(xù)型隨機(jī)變量X的概率密度函數(shù)為其中>0為常數(shù)，則稱X服從參數(shù)為的指數(shù)分布。指數(shù)分布的期望值為

排隊(duì)服務(wù)系統(tǒng)中顧客到達(dá)率為常數(shù)時(shí)的到達(dá)間隔、故障率為常數(shù)時(shí)零件的壽命都服從指數(shù)分布。指數(shù)分布在排隊(duì)論、可靠性分析中有廣泛應(yīng)用。注意：Matlab中，產(chǎn)生參數(shù)為的指數(shù)分布的命令為exprnd()例顧客到達(dá)某商店的間隔時(shí)間服從參數(shù)為0.1的指數(shù)分布

指數(shù)分布的均值為1/0.1=10。指兩個(gè)顧客到達(dá)商店的平均間隔時(shí)間是10個(gè)單位時(shí)間.即平均10個(gè)單位時(shí)間到達(dá)1個(gè)顧客.顧客到達(dá)的間隔時(shí)間可用exprnd(10)模擬。設(shè)離散型隨機(jī)變量X的所有可能取值為0,1,2,…,且取各個(gè)值的概率為其中>0為常數(shù)，則稱X服從參數(shù)為的帕松分布。帕松分布在排隊(duì)系統(tǒng)、產(chǎn)品檢驗(yàn)、天文、物理等領(lǐng)域有廣泛應(yīng)用。帕松分布的期望值為如相繼兩個(gè)事件出現(xiàn)的間隔時(shí)間服從參數(shù)為的指數(shù)分布，則在單位時(shí)間間隔內(nèi)事件出現(xiàn)的次數(shù)服從參數(shù)為的泊松分布．即單位時(shí)間內(nèi)該事件出現(xiàn)k次的概率為：反之亦然。指數(shù)分布與帕松分布的關(guān)系：(1)指兩個(gè)顧客到達(dá)商店的平均間隔時(shí)間是10個(gè)單位時(shí)間.即平均10個(gè)單位時(shí)間到達(dá)1個(gè)顧客.(2)指一個(gè)單位時(shí)間內(nèi)平均到達(dá)0.1個(gè)顧客例(1)顧客到達(dá)某商店的間隔時(shí)間服從參數(shù)為0.1的指數(shù)分布

(2)該商店在單位時(shí)間內(nèi)到達(dá)的顧客數(shù)服從參數(shù)為0.1的帕松分布排隊(duì)問(wèn)題③機(jī)械故障等候修理④飛機(jī)跑道日常生活中經(jīng)常遇到的排隊(duì)問(wèn)題：①自選商場(chǎng)收款臺(tái)②醫(yī)院里病人等候就診輸入情況：顧客到達(dá)時(shí)間和服務(wù)時(shí)間。系統(tǒng)狀態(tài)：排隊(duì)等候的顧客數(shù)目（隊(duì)長(zhǎng)）L(t)

服務(wù)員是否在工作或服務(wù)效率等；簡(jiǎn)圖：第二顧客接受服務(wù)時(shí)間s2x50x1x2x3x4y1y2y3y4y5D2關(guān)系：系統(tǒng)在什么條件下處于空閑狀態(tài)？(yi<xi+1)

排隊(duì)系統(tǒng)中，顧客到達(dá)時(shí)刻數(shù)據(jù)如何收集？對(duì)每個(gè)顧客的服務(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ù)：k123456789…ik013413182…sk231341341…其中ik表示第k位顧客與第k-1位顧客到達(dá)的時(shí)間間隔，sk是第k位顧客的服務(wù)時(shí)間，下面就以這組數(shù)據(jù)作為系統(tǒng)的輸入。3、系統(tǒng)的狀態(tài)隊(duì)長(zhǎng)：排隊(duì)等候的顧客數(shù)目。記作L(t).服務(wù)員的狀態(tài)用S(t)表示，當(dāng)服務(wù)員工作時(shí)，令S(t)=1，空閑時(shí)令S(t)=0.事件：引起系統(tǒng)狀態(tài)L(t)和S(t)改變的行為。

“顧客到達(dá)”和“顧客離開(kāi)”設(shè)第k位顧客到達(dá)和離開(kāi)時(shí)刻分別為ak和dk，則有

4、初始條件和終止條件不妨假設(shè)t=0是第1位顧客到達(dá)的時(shí)刻，此時(shí)服務(wù)員處于空閑狀態(tài)。模擬的終止時(shí)刻是T（給定值），或模擬直到第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.AT<DTWL=0,SS=0,AT=0,DT=999置t=ATSS=0置t=DTWL>0WL=WL+1令SS＝1產(chǎn)生s置DT=t+s產(chǎn)生i置AT＝t+it>=T令SS＝0WL＝WL－1產(chǎn)生s置DT＝t+s停止置DT＝999是否是否是否否是6、系統(tǒng)的性能指標(biāo)平均隊(duì)長(zhǎng)平均待時(shí)間服務(wù)利用率7、與分析模型的比較對(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）8、其它的排隊(duì)服務(wù)系統(tǒng)有兩個(gè)（或多個(gè)）服務(wù)員串聯(lián)服務(wù)系統(tǒng)，每個(gè)顧客接受多次服務(wù)從排隊(duì)顧客中選服務(wù)時(shí)間最短的最先服務(wù)隊(duì)的長(zhǎng)度有限制帶有優(yōu)化目標(biāo)的排隊(duì)服務(wù)系統(tǒng)

庫(kù)存太多造成浪費(fèi)或資金積壓，庫(kù)存太少不能滿足需求也造成損失，需要進(jìn)行決策：何時(shí)進(jìn)貨，進(jìn)多少，使得平均費(fèi)用最少，而收益最大。庫(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ù)量？)模擬的基本思路：需要獲得市場(chǎng)對(duì)該產(chǎn)品的需要量的樣本值；(產(chǎn)生抽樣值)②尋找貨存量與需要量之間的關(guān)系；③按照兩種不同方案（貨存量≥{≤}需要量）計(jì)算出經(jīng)T天后企業(yè)的利潤(rùn)值(累計(jì))；④比較大小，從中選出一個(gè)較優(yōu)的方案；方案甲：貨存量Y1分析兩種狀態(tài)下X＞Y1和X≤Y1的利潤(rùn)值L1置初始狀態(tài)獲取產(chǎn)品的需求量XX=normrnd（135，22.5）N=1方案乙：貨存量Y2分析兩種狀態(tài)下X＞Y2和X≤Y2的利潤(rùn)值L2累計(jì)利潤(rùn)值TL1累計(jì)利潤(rùn)值TL2N=N+1判斷：N≤T(T天)比較計(jì)算：max{TL1，TL2}輸出最大利潤(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)觀察計(jì)算結(jié)果如果多數(shù)情況有：TL1<TL2,則認(rèn)為方案乙較優(yōu).思考：1）一般模擬結(jié)果波動(dòng)性較大，如何減少這種波動(dòng)？2）修改上述程序?？煽啃詥?wèn)題

一設(shè)備上有三個(gè)相同的軸承，每個(gè)軸承正常工作壽命為隨機(jī)變量，概率分布如下：壽命/h1000110012001300140015001600170018001900概率0.100.130.250.130.090.120.020.060.050.05有軸承損壞→設(shè)備停止工作→檢修工準(zhǔn)備開(kāi)始更換部件，稱為一個(gè)延遲時(shí)間，它也是隨機(jī)變量，分布如下：延遲時(shí)間/min51015概率0.60.30.1主要費(fèi)用：1、設(shè)備停工損失費(fèi)：5元/分鐘；2、檢修工人的工時(shí)費(fèi)：12元/小時(shí)；3、軸承的成本費(fèi)：16元/個(gè)更換軸承所需要的時(shí)間：一個(gè)兩個(gè)三個(gè)

203040（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)生見(jiàn)Lrnd.m（零件壽命）

Yrnd.m（延遲時(shí)間）

方案一的數(shù)學(xué)模型：kekao1.m目標(biāo)函數(shù)minc=∑Ui

/T其中損失費(fèi)工時(shí)費(fèi)成本費(fèi)ti

表示延遲時(shí)間

設(shè)事件發(fā)生時(shí)刻Ti，它是由軸承工作壽命Li、延遲時(shí)間和修理時(shí)間迭加產(chǎn)生?？偟倪\(yùn)行時(shí)間為T(mén)（10000）。例：（手工模擬）事件類型發(fā)生時(shí)刻延遲時(shí)間A1400h5minB1500h15minC1500h15min方案一的初始事件表t=0→t=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)事件類型發(fā)生時(shí)刻延遲時(shí)間B1500h15minC1500h15minA2400h25min5min第一次刷新后的事件表t=1400→t=1500；

cost=145+(15+30)×5+12×1/2+16×2=408；③產(chǎn)生下一個(gè)B事件發(fā)生時(shí)刻1500+45/60+1200（Li）=2700小時(shí)45分鐘c1=5;c2=12;c3=16;(kekao1