Latin超立方采樣技術(shù)及其在結(jié)構(gòu)可靠性分析中的應(yīng)用二、 拉丁超立方采樣通過第一部分前言對(duì)拉丁超立方抽樣概念的理解，結(jié)合第二部分對(duì)其原理的介紹，用matlab編輯拉丁超立方抽樣的程序，具體如下：functions=lhsamp(n,k)%產(chǎn)生一個(gè)n行k列的拉丁超立方抽樣矩陣%s:元素介于(0.0,1.0)之間n*k的拉丁超立方抽樣矩陣%k:輸入變量數(shù)(維數(shù))%n:每個(gè)輸入變量抽取的樣本數(shù)(每一維的采樣數(shù))s=zeros(n,k);fori=1:k;s(:,i)=rand(1,n)'/n+(randperm(n)-1)'/n;%rand(1,n):產(chǎn)生n個(gè)值介于0.0到1.0的隨機(jī)數(shù)%randperm(n):產(chǎn)生正整數(shù)1,2,3，...n的隨機(jī)排列end%得到拉丁超立方抽樣矩陣后，根據(jù)每個(gè)變量的分布函數(shù)，根據(jù)Xnk=f-1(Un)之間的關(guān)系，由每個(gè)變量對(duì)應(yīng)的抽樣結(jié)果Un反算出對(duì)該變量的真實(shí)抽樣點(diǎn)Xk三、 統(tǒng)計(jì)相關(guān)的減小方程1、Latin超立方抽樣可能隨機(jī)的引進(jìn)了一定的統(tǒng)計(jì)相關(guān)，所以用文中采用的Spearman系數(shù)法來減小統(tǒng)計(jì)相關(guān)性。首先，根據(jù)Spearman相關(guān)系數(shù)的計(jì)算公式，用matlab編輯程序如下:functioncoeff=Spearman(X,Y)%本函數(shù)用于實(shí)現(xiàn)Spearman相關(guān)系數(shù)的計(jì)算操作%X輸入的數(shù)值序列%Y輸入的數(shù)值序列%coeff兩個(gè)輸入數(shù)值序列X，Y的相關(guān)系數(shù)iflength(X)~=length(Y)error'兩個(gè)數(shù)值數(shù)列的維數(shù)不相等');endN=length(X);%得到序列的長度Xrank=zeros(1,N);%存儲(chǔ)X中各元素的排行Yrank=zeros(1,N);%存儲(chǔ)Y中各元素的排行%計(jì)算Xrank中的各個(gè)值fori=1:Ncount=1;forj=1:NifX(i)<X(j)count=count+1;endendXrank(i)=count;end%計(jì)算Yrank中的各個(gè)值fori=1:Ncount=1;forj=1:NifY(i)<Y(j)count=count+1;endendYrank(i)=count;end%利用X,Y的序數(shù)排列計(jì)算Spearman相關(guān)系數(shù)A=6*sum((Xrank-Yrank).入2);B=N*(N—1)*(N+1);coeff=1-A/B;end2、 仿照文中例子做一個(gè)K=5個(gè)輸入變量和N=10個(gè)模擬的算例，以驗(yàn)證上述程序的正確性及Spearman系數(shù)對(duì)統(tǒng)計(jì)相關(guān)性的減小作用。3、 首先用sample.m函數(shù)生成一個(gè)K=5個(gè)輸入變量和N=10個(gè)模擬的秩數(shù)隨機(jī)排列表，見表1：functionR=sample(n,k)%產(chǎn)生一個(gè)n行k列的隨機(jī)抽樣矩陣%k:輸入變量數(shù)(維數(shù))%n:每個(gè)輸入變量抽取的樣本數(shù)(每一維的采樣數(shù))R=zeros(n,k);fori=1:k;R(:,i)=randperm(n)';%randperm(n):產(chǎn)生正整數(shù)1,2,3，...n的隨機(jī)排列

end表1 K=5個(gè)輸入變量和N=10個(gè)模擬的秩數(shù)隨機(jī)排列表模擬1未修正表變量2345166225231101023778414849575591366155887287798436110991036310102494修正表模擬變量12345166123231782377103148494105592546154967288698435189910675101023107矩陣各列間的統(tǒng)計(jì)相關(guān)由序相關(guān)矩陣T描述，其元素Tij是R的i列和j列間的Spearman系數(shù)。在matlab命令窗口執(zhí)行T.m腳本文件，調(diào)用前述計(jì)算Spearman系數(shù)的Spearman.m程序得到秩相關(guān)矩陣T：腳本文件T.m：forj=1:5;fori=1:5;T(i,j)=Spearman(R(:,i),R(:,j));endend

表2表1秩數(shù)隨機(jī)排列表的秩相關(guān)矩陣未修正表變量1變量523411.00000.0667-0.2121-0.0909-0.490920.06671.0000-0.5152-0.3455-0.01823-0.2121-0.51521.00000.3697-0.09094-0.0909-0.34550.36971.0000-0.29705-0.4909-0.0182-0.0909-0.29701.0000變量12修正表變量34511.00000.06670.0061-0.0061-0.006120.06671.0000-0.0061-0.1758-0.115230.0061-0.00611.0000-0.09090.11524-0.0061-0.1758-0.09091.00000.04245-0.0061-0.11520.11520.04241.0000T是正定的對(duì)稱矩陣，可用Choiesky分解將T分解為T=Q*Qt在matlab命令運(yùn)行窗口運(yùn)行Q=chol(T)，得到Q：q=1.00000.0667-0.2121-0.0909-0.490900.9978-0.5021-0.34020.0146000.83840.2142-0.22390000.9111-0.316800000.7799修正后的RB=R*Q-i,在matlab中運(yùn)行RB=R*inv(Q),得到RB：B B BRb=6.0000 5.6125 7.26513.1808 13.46053.0000 0.8018 13.16718.4781 11.66197.0000 6.5479 15.23513.9508 11.54478.0000 3.4744 14.84014.0936 19.86925.0000 8.6860 7.66015.2331 15.00291.0000 4.94439.1782 8.5679 16.91002.0000 7.8842 13.57767.6332 19.65006.0000 5.6125 7.26513.1808 13.46053.0000 0.8018 13.16718.4781 11.66197.0000 6.5479 15.23513.9508 11.54478.0000 3.4744 14.84014.0936 19.86925.0000 8.6860 7.66015.2331 15.00291.0000 4.94439.1782 8.5679 16.91002.0000 7.8842 13.57767.6332 19.65004.00002.73949.80950.212818.19109.00009.421011.49808.296816.006810.00001.33638.10169.469617.5709按矩陣RB列中次序重新排列輸入矩陣R中的值，則序數(shù)隨機(jī)排列表各列間統(tǒng)計(jì)相關(guān)性便減小了，得到修正后的R,見表1。然后根據(jù)修正后的R,驗(yàn)證相關(guān)性是否減小，及求出修正后的T,見表2.結(jié)論：由表2易得，修正后的各列間的統(tǒng)計(jì)相關(guān)性明顯減小，用Spearman系數(shù)法可有效減小拉丁超立方抽樣隨機(jī)引進(jìn)的統(tǒng)計(jì)相關(guān)。四、在結(jié)構(gòu)可靠性分析中的應(yīng)用1)例題1由題意知，F(xiàn)y、S分布類型和參數(shù)已知，隨機(jī)變量Fy和S用前述的Latin超立方采樣和逆變換隨機(jī)產(chǎn)生，求得分布函數(shù)值后，通過逆分布函數(shù)變換產(chǎn)生隨機(jī)變量。在matlab命令運(yùn)行窗口中運(yùn)行Pf.m腳本文件：p=lhsamp(N,2);Fy=norminv(p(:,1),262000000,26200000);S=norminv(p(:,2),0.00082,4.1e-5);Pfi=1-exp(-(97722./(Fy.*S)).入5.18);Pf=sum(Pfi)/N分別取N=10,20,30,50,70,100,150,200,250,300,350,400,得到相應(yīng)的不用模擬次數(shù)的失效概率Pf。用matlab繪圖如下：在matlab命令運(yùn)行窗口輸入下列指令運(yùn)行x=[10,20,30,50,70,100,150,200,250,300,350,400];y=[0.0186,0.0212,0.0207,0.0207,0.0203,0.0203,0.0207,0.0204,0.0203,0.0205,0.0208,0.0208];plot(x,y)axis([0,400,0.018,0.024]);xlabel('Latin超立方采樣模擬數(shù)')；ylabel('失效概率Pf);title('可靠性計(jì)算結(jié)果')

bb率概效失2）例題2由題意知，若對(duì)非線性結(jié)構(gòu)進(jìn)行地震可靠性分析，首先需分別建立結(jié)構(gòu)模型集和地震時(shí)程集。然后可利用Latin超立方采樣技術(shù)將這些地震時(shí)程和結(jié)構(gòu)模型匹配成81個(gè)地震一結(jié)構(gòu)系統(tǒng)。1、結(jié)構(gòu)模型集：ZagasYy12%0.10.010.00222%0.10.010.002532%0.10.010.00342%0.10.030.00252%0.10.030.002562%0.10.030.00372%0.10.050.00282%0.10.050.002592%0.10.050.003102%0.170.010.002

112%0.170.010.0025122%0.170.010.003132%0.170.030.002142%0.170.030.0025152%0.170.030.003162%0.170.050.002172%0.170.050.0025182%0.170.050.003192%0.250.010.002202%0.250.010.0025212%0.250.010.003222%0.250.030.002232%0.250.030.0025242%0.250.030.003252%0.250.050.002262%0.250.050.0025272%0.250.050.003284%0.10.010.002294%0.10.010.0025304%0.10.010.003314%0.10.030.002324%0.10.030.0025334%0.10.030.003344%0.10.050.002354%0.10.050.0025364%0.10.050.003374%0.170.010.002384%0.170.010.0025394%0.170.010.003404%0.170.030.002414%0.170.030.0025424%0.170.030.003434%0.170.050.002444%0.170.050.0025454%0.170.050.003464%0.250.010.002474%0.250.010.0025484%0.250.010.003494%0.250.030.002504%0.250.030.0025514%0.250.030.003524%0.250.050.002534%0.250.050.0025544%0.250.050.003556%0.10.010.002566%0.10.010.0025576%0.10.010.003586%0.10.030.002596%0.10.030.0025606%0.10.030.003616%0.10.050.002626%0.10.050.0025636%0.10.050.003646%0.170.010.002656%0.170.010.0025666%0.170.010.003676%0.170.030.002686%0.170.030.0025696%0.170.030.003706%0.170.050.002716%0.170.050.0025726%0.170.050.003736%0.250.010.002746%0.250.010.0025花6%0.250.010.003766%0.250.030.002776%0.250.030.0025786%0.250.030.003796%0.250.050.002806%0.250.050.0025816%0.250.050.003

2、地震時(shí)程集:wgZtf(t)13兀0.65f(t)123兀0.65f(t)233兀0.65f(t)343兀0.610f(t)153兀0.610f(t)263兀0.610f(t)373兀0.61.5f(t)183兀0.61.5f(t)293兀0.61.5f(t)3103兀0.365f(t)1113兀0.365f(t)2123兀0.365f(t)3133兀0.3610f(t)1143兀0.3610f(t)2153兀0.3610f(t)3163兀0.361.5f(t)1173兀0.361.5f(t)2183兀0.361.5f(t)3193兀0.845f(t)1203兀0.845f(t)2213兀0.845f(t)3223兀0.8410f(t)1233兀0.8410f(t)2243兀0.8410f(t)3253兀0.841.5f(t)1263兀0.841.5f(t)2273兀0.841.5f(t)3286兀0.65f(t)1296兀0.65f(t)2306兀0.65f(t)3316兀0.610f(t)1326兀0.610f(t)2336兀0.610f(t)3

346兀0.61.5f(t)1356兀0.61.5f(t)2366兀0.61.5f(t)3376兀0.365f(t)1386兀0.365f(t)2396兀0.365f(t)3406兀0.3610f(t)1416兀0.3610f(t)2426兀0.3610f(t)3436兀0.361.5f(t)1446兀0.361.5f(t)2456兀0.361.5f(t)3466兀0.845f(t)1476兀0.845f(t)2486兀0.845f(t)3496兀0.8410f(t)1506兀0.8410f(t)2516兀0.8410f(t)3526兀0.841.5f(t)1536兀0.841.5f(t)2546兀0.841.5f(t)3559兀0.65f(t)1569兀0.65f(t)2579兀0.65f(t)3589兀0.610f(t)1599兀0.610f(t)2609兀0.610f(t)3619兀0.61.5f(t)1629兀0.61.5f(t)2639兀0.61.5f(t)3649兀0.365f(t)1659兀0.365f(t)2669兀0.365f(t)3679兀0.3610f(t)1689兀0.3610f(t)2699兀0.3610f(t)3709兀0.361.5f(t)1719兀0.361.5f(t)2729兀0.361.5f(t)3739兀0.845f(t)1749兀0.845f(t)2花9兀0.845f(t)3769兀0.8410f(t)1779兀0.8410f(t)2789兀0.8410f(t)3799兀0.841.5f(t)1809兀0.841.5f(t)2819兀0.841.5f(t)3在matla