蒙特卡洛模擬

《蒙特卡洛模擬》由會員上傳分享，免費(fèi)在線閱讀，更多相關(guān)內(nèi)容在教育資源-天天文庫。

1、是一個(gè)論文里看到的，內(nèi)容如下：兩元件A,B組成的串聯(lián)系統(tǒng),已知兩元件A,B的壽命分布均為威布爾分布,則其概率密度函數(shù)和分布函數(shù)分別為:式中,m為形狀參數(shù);η為尺度參數(shù);δ為位置參數(shù)。由上式可得威布爾分布的直接抽樣方法如下:假設(shè)元件A的參數(shù)γ=0小時(shí),β=5,η=80小時(shí);元件B的參數(shù)γ=0小時(shí),β=4,η=150小時(shí)。則故障率的計(jì)算及曲線的繪制步驟為:(1)根據(jù)式(9)對元件進(jìn)行抽樣;(2)對系統(tǒng)進(jìn)行50000次仿真;(3)統(tǒng)計(jì)壽命區(qū)間Δt=20小時(shí)時(shí),每個(gè)區(qū)間的故障次數(shù);(4)計(jì)算各個(gè)元件的故障率及系統(tǒng)的故障率,對于串聯(lián)系統(tǒng)來說系統(tǒng)的故障率為各元

2、件故障率的乘積;(5)繪制故障率曲線,如圖2所示,其中A表示元件A的故障率,B表示元件B的故障率,S表示系統(tǒng)的故障率。結(jié)果是這樣的：表格1元件和系統(tǒng)的故障率計(jì)算統(tǒng)計(jì)間隔△tA的故障率λaB的故障率λb系統(tǒng)的故障率λ0000200.0012511e-0060.001252400.0102823e-0050.010312600.0225630.0020720.024635800.0172140.005940.0231541000.0070120.0083090.0153211200.0128450.0106960.0235411400.0176870.

3、0096380.0273251600.0145680.0089770.0235451800.0115460.0056680.0172142000.0136980.0055150.0192132200.0156250.0058760.0215012400.0135460.0078770.0214232600.0128540.0084480.0213022800.0136250.0090070.0226323000.0140320.0088130.0228453200.0131420.0082120.0213543400.0130210.0069330

4、.0199543600.0139850.0065550.020543800.0140320.0067530.0207854000.0137210.0073110.0210324200.0136540.0073210.0209754400.0136140.0073550.0209694600.0136320.0073350.0209674800.0136050.0073680.0209735000.0136120.0073570.0209695200.0136150.0073570.0209725400.0136140.0073490.0209635

5、600.0136130.0073540.0209675800.0136140.0073550.020969主程序functionmain()n=50000;%仿真次數(shù)%產(chǎn)生均勻分布隨機(jī)數(shù)z1=rand(1,n);%元件A隨機(jī)m1=5;eta1=80;sgma1=0;lambda1=webull(z1,m1,eta1,sgma1,n);z2=rand(1,n);%元件B隨機(jī)m2=4;eta2=150;sgma2=0;figurelambda1=webull(z2,m2,eta2,sgma2,n);disp('系統(tǒng)的故障率')lambda=lambda

6、1+lambda2;調(diào)用的子程序functionlambda=webull(z,m,eta,sgma,n)E=eta*(-log(1-z)).^(1/m)+sgma;%轉(zhuǎn)換成威布爾分布的隨機(jī)數(shù)y=zeros(1,30);fori=1:nfort=1:30if(E(i)>20*t)&&(E(i)<=20*(t+1))%對于每一個(gè)隨機(jī)數(shù)判斷在哪個(gè)微區(qū)間內(nèi)，則該區(qū)間的故障數(shù)加1y(t)=y(t)+1;endendendy=y/n;%故障次數(shù)除以抽樣次數(shù)為故障率;t=1:30;plot(t,y,'.-')gridonxlabel('時(shí)間t（間隔為20小時(shí)）

7、')ylabel('故障率')title('故障率曲線')legend('元件A','元件B','系統(tǒng)')gridofflambda=y;