Latin超立方抽样学习报告-牛亚运

1、页脚Latin超立方采样技术及其在结构可靠性分析中的应用超立方采样技术及其在结构可靠性分析中的应用二、拉丁超立方采样二、拉丁超立方采样通过第一部分前言对拉丁超立方抽样概念的理解，结合第二部分对其原理的介绍，用 matlab 编辑拉丁超立方抽样的程序，具体如下：functions=lhsamp(n,k)%产生一个 n 行 k 列的拉丁超立方抽样矩阵%s：元素介于(0.0,1.0)之间 n*k 的拉丁超立方抽样矩阵%k：输入变量数(维数)%n:每个输入变量抽取的样本数(每一维的采样数)s=zeros(n,k);fori=1：k;s(：,i)=rand(1,n)/n+(randperm(n)-1)/

2、n;%rand(l,n):产生 n 个值介于 0.0 到 1.0 的随机数%randperm(n):产生正整数 1,2,3，.n 的随机排列end%得到拉丁超立方抽样矩阵后，根据每个变量的分布函数，根据 Xnk=f-1(Un)之间的关系，由每个变量对应的抽样结果 Un 反算出对该变量的真实抽样点 Xk三、统计相关的减小方程三、统计相关的减小方程1、Latin 超立方抽样可能随机的引进了一定的统计相关，所以用文中采用的 Spearman 系数法来减小统计相关性。首先，根据 Spearman 相关系数的计算公式，用 matlab 编辑程序如下：functioncoeff=Spearman(X,Y)

3、%本函数用于实现 Spearman 相关系数的计算操作%X：输入的数值序列%Y：输入的数值序列%coeff：两个输入数值序列 X,Y 的相关系数iflength(X)=length(Y)error(两个数值数列的维数不相等)；endN=length(X);%得到序列的长度Xrank=zeros(1,N);%存储 X 中各元素的排行Yrank=zeros(1,N);%存储 Y 中各元素的排行%计算 Xrank 中的各个值页脚fori=1:Ncount=1;forj=1:NifX(i)X(j)count=count+1;endendXrank(i)=count;end%计算 Yrank 中的各个值

4、fori=1:Ncount=1;forj=1:NifY(i)Y(j)count=count+1;endendYrank(i)=count;end%利用 X,Y 的序数排列计算 Spearman 相关系数A=6*sum(Xrank-Yrank).2)；B=N*(N-1)*(N+1);coeff=1-A/B；end2、仿照文中例子做一个 K=5 个输入变量和 N=10 个模拟的算例，以验证上述程序的正确性及 Spearman 系数对统计相关性的减小作用。3、首先用 sample.m 函数生成一个 K=5 个输入变量和 N=10 个模拟的秩数随机排列表，见表 1:functionR=sample(n

5、,k)%产生一个 n 行 k 列的随机抽样矩阵%k：输入变量数(维数)%n：每个输入变量抽取的样本数(每一维的采样数)R=zeros(n,k)；fori=1：k；R(：,i)=randperm(n)；%randperm(n)：产生正整数 1,2,3,.n 的随机排列页脚end表表1K=5个输入变量和个输入变量和N=10个模拟的秩数随机排列表个模拟的秩数随机排列表模拟1未修正表变量5234166225231101023778414849575591366155887287798436110991036310102494模拟12修正表变量345166123231782377103148494105

6、592546154967288698435189910675101023107矩阵各列间的统计相关由序相关矩阵T描述，其元素Tij是R的i列和j列间的Spearman系数。 在matlab命令窗口执行T.m脚本文件， 调用前述计算Spearman系数的Spearman.m程序得到秩相关矩阵T：脚本文件 T.m：forj=1:5;fori=1:5;T(i,j)=Spearman(R(:,i),R(:,j);endend页脚表表2表表1秩数随机排列表的秩相关矩阵秩数随机排列表的秩相关矩阵变量未修正表变量5123411.00000.0667-0.2121-0.0909-0.490920.06671.

7、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.000修正表变量12变量34511.00000.06670.0061-0.0061-0.006120.06671.0000-0.0061-0.1758-0.115230.0061-0.00611.0000-0.09090.1154-0.0061-0.1758-0.09091.00000.0425-0.0061-0.11520.11520.04241.

8、000T是正定的对称矩阵，可用Choiesky分解将T分解为T=Q*QT在matlab命令运行窗口运行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修正后的R=R*Q-1，在matlab中运行R=R*inv(Q)，得到R：BBBR=B6.00005.61257.26513.180813.46053.00000.801813.16718.478111.6617.00006.547915.23513.

9、950811.5448.00003.474414.84014.093619.8695.00008.68607.66015.233115.0021.00004.94439.17828.567916.9102.00007.884213.57767.633219.6504.00002.73949.80950.212818.1919.00009.421011.49808.296816.00610.00001.33638.10169.469617.570页脚按矩阵RB列中次序重新排列输入矩阵R中的值，则序数随机排列表各列间统计相关性便减小了，得到修正后的R,见表1。然后根据修正后的R,验证相关性是否减小

10、，及求出修正后的T,见表2.结论结论：由表2易得，修正后的各列间的统计相关性明显减小，用Spearman系数法可有效减小拉丁超立方抽样随机引进的统计相关。四、在结构可靠性分析中的应用四、在结构可靠性分析中的应用1) )例题例题1由题意知，Fy、S分布类型和参数已知，随机变量Fy和S用前述的Latin超立方采样和逆变换随机产生，求得分布函数值后，通过逆分布函数变换产生随机变量。在matlab命令运行窗口中运行Pf.m脚本文件：p=lhsamp(N,2);Fy=norminv(p(:,1),262000000,26200000);S=norminv(p(:,2),0.00082,4.1e-5);P

11、fi=1-exp(-(97722./(Fy.*S).5.18)；Pf=sum(Pfi)/N分别取N=10,20,30,50,70,100,150,200,250,300,350,400,得到相应的不用模拟次数的失效概率Pf。用matlab绘图如下：在matlab命令运行窗口输入下列指令运行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.01

12、8,0.024)；xlabel(Latin超立方采样模拟数);ylabel(失效概率Pf);title(可靠性计算结果)页脚50100150200250300350400Latin超立方采样模拟数N2）例题例题2由题意知，若对非线性结构进行地震可靠性分析，首先需分别建立结构模型集和地震时程集。 然后可利用Latin超立方采样技术将这些地震时程和结构模型匹配成81个地震结构系统。1、结构模型集：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%

13、0.10.050.00282%0.10.050.002592%0.10.050.003102%0.170.010.0020.0230.0220.0210.020.019可靠性计算结果0.0180页脚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.

14、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.1

15、70.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.002页脚474%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

16、.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.0025756%0.2

17、50.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、地震时程集：wgZgtf(t)13p0.65f(t)123p0.65f(t)233p0.65f(t)343p0.610f(t)153p0.610f(t)263p0.610f(t)373p0.61.5f(t)183p0.61.5f(t)293p0.61.5f(t)3103p0.365f(t)1113p0.365f(t)2123p0.365f(t)3133p0.36

18、10f(t)1143p0.3610f(t)2153p0.3610f(t)3163p0.361.5f(t)1173p0.361.5f(t)2183p0.361.5f(t)3193p0.845f(t)1203p0.845f(t)2213p0.845f(t)3223p0.8410f(t)1233p0.8410f(t)2243p0.8410f(t)3253p0.841.5f(t)1263p0.841.5f(t)2273p0.841.5f(t)3286p0.65f(t)1296p0.65f(t)2306p0.65f(t)3316p0.610f(t)1326p0.610f(t)2336p0.610f(t)

19、3页脚346p0.61.5f(t)1356p0.61.5f(t)2366p0.61.5f(t)3376p0.365f(t)1386p0.365f(t)2396p0.365f(t)3406p0.3610f(t)1416p0.3610f(t)2426p0.3610f(t)3436p0.361.5f(t)1446p0.361.5f(t)2456p0.361.5f(t)3466p0.845f(t)1476p0.845f(t)2486p0.845f(t)3496p0.8410f(t)1506p0.8410f(t)2516p0.8410f(t)3526p0.841.5f(t)1536p0.841.5f(t

20、)2546p0.841.5f(t)3559p0.65f(t)1569p0.65f(t)2579p0.65f(t)3589p0.610f(t)1599p0.610f(t)2609p0.610f(t)3619p0.61.5f(t)1629p0.61.5f(t)2639p0.61.5f(t)3649p0.365f(t)1659p0.365f(t)2669p0.365f(t)3679p0.3610f(t)1689p0.3610f(t)2699p0.3610f(t)3页脚709p0.361.5f(t)1719p0.361.5f(t)2729p0.361.5f(t)3739p0.845f(t)1749p0.845f(t)2759p0.845f(t)3769p0.8410f(t)1779p0.8410f(t)2789p0.8410f(t)3799p0.841.5f(t)1809p0.841.5f(t)2819p0.841.5f(t)3在matlab中调用sample.m函数，functionR=s