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

. . Latin超立方采样技术与其在结构可靠性分析中的应用二、 拉丁超立方采样通过第一部分前言对拉丁超立方抽样概念的理解，结合第二部分对其原理的介绍，用matlab编辑拉丁超立方抽样的程序，具体如下：function s=lhsamp(n,k)%产生一个n行k列的拉丁超立方抽样矩阵%s:元素介于（0.0,1.0）之间n*k的拉丁超立方抽样矩阵%k:输入变量数(维数）%n:每个输入变量抽取的样本数（每一维的采样数）s=zeros(n,k);for i=1:k; s(:,i)=rand(1,n)/n+(randperm(n)-1)/n;%rand(1,n)：产生n个值介于0.0到1.0的随机数%randperm(n):产生正整数1,2,3，.n的随机排列end%得到拉丁超立方抽样矩阵后，根据每个变量的分布函数，根据Xnk=f-1(Un)之间的关系,由每个变量对应的抽样结果Un反算出对该变量的真实抽样点Xk三、 统计相关的减小方程1、Latin超立方抽样可能随机的引进了一定的统计相关，所以用文中采用的Spearman系数法来减小统计相关性。首先，根据Spearman相关系数的计算公式，用matlab编辑程序如下：function coeff = Spearman(X , Y) % 本函数用于实现Spearman相关系数的计算操作 % X：输入的数值序列 % Y：输入的数值序列 % coeff：两个输入数值序列X，Y的相关系数 if length(X) = length(Y) error(两个数值数列的维数不相等); endN = length(X); %得到序列的长度 Xrank = zeros(1 , N); %存储X中各元素的排行 Yrank = zeros(1 , N); %存储Y中各元素的排行 %计算Xrank中的各个值 for i = 1 : N count = 1;for j = 1 : N if X(i) X(j) count = count + 1; endend Xrank(i) = count ;end%计算Yrank中的各个值 for i = 1 : N count = 1; for j = 1 : N if Y(i) Y(j) count = count + 1; endend Yrank(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:function R=sample(n,k)%产生一个n行k列的随机抽样矩阵%k:输入变量数(维数）%n:每个输入变量抽取的样本数（每一维的采样数）R=zeros(n,k);for i=1:k; R(:,i)=randperm(n);%randperm(n):产生正整数1,2,3，.n的随机排列end表1 K=5个输入变量和N=10个模拟的秩数随机排列表 未修正表 模拟 变量 1 2 3 4 5 1 6 6 2 2 5 2 3 1 10 10 2 3 7 7 8 4 1 4 8 4 9 5 7 5 5 9 1 3 6 6 1 5 5 8 8 7 2 8 7 7 9 8 4 3 6 1 10 9 9 10 3 6 3 10 10 2 4 9 4 修正表 模拟 变量 1 2 3 4 5 1 6 6 1 2 3 2 3 1 7 8 2 3 7 7 10 3 1 4 8 4 9 4 10 5 5 9 2 5 4 6 1 5 4 9 6 7 2 8 8 6 9 8 4 3 5 1 8 9 9 10 6 7 5 10 10 2 3 10 7 矩阵各列间的统计相关由序相关矩阵T描述，其元素Tij是R的i列和j列间的Spearman系数。在matlab命令窗口执行T.m脚本文件，调用前述计算Spearman系数的Spearman.m程序得到秩相关矩阵T：脚本文件T.m：for j=1:5;for i=1:5; T(i,j)=Spearman(R(:,i),R(:,j);endend 表2 表1秩数随机排列表的秩相关矩阵 未修正表 变量 变量 1 2 3 4 5 1 1.0000 0.0667 -0.2121 -0.0909 -0.4909 2 0.0667 1.0000 -0.5152 -0.3455 -0.0182 3 -0.2121 -0.5152 1.0000 0.3697 -0.0909 4 -0.0909 -0.3455 0.3697 1.0000 -0.2970 5 -0.4909 -0.0182 -0.0909 -0.2970 1.0000 修正表变量 变量 1 2 3 4 5 1 1.0000 0.0667 0.0061 -0.0061 -0.0061 2 0.0667 1.0000 -0.0061 -0.1758 -0.1152 3 0.0061 -0.0061 1.0000 -0.0909 0.1152 4 -0.0061 -0.1758 -0.0909 1.0000 0.0424 5 -0.0061 -0.1152 0.1152 0.0424 1.0000T是正定的对称矩阵，可用Choiesky分解将T分解为T=Q*QT 在matlab命令运行窗口运行Q=chol(T)，得到Q：Q = 1.0000 0.0667 -0.2121 -0.0909 -0.4909 0 0.9978 -0.5021 -0.3402 0.0146 0 0 0.8384 0.2142 -0.2239 0 0 0 0.9111 -0.3168 0 0 0 0 0.7799修正后的RB=R*Q-1，在matlab中运行RB=R*inv(Q)，得到RB：RB= 6.0000 5.6125 7.2651 3.1808 13.4605 3.0000 0.8018 13.1671 8.4781 11.6619 7.0000 6.5479 15.2351 3.9508 11.5447 8.0000 3.4744 14.8401 4.0936 19.8692 5.0000 8.6860 7.6601 5.2331 15.0029 1.0000 4.9443 9.1782 8.5679 16.9100 2.0000 7.8842 13.5776 7.6332 19.6500 4.0000 2.7394 9.8095 0.2128 18.1910 9.0000 9.4210 11.4980 8.2968 16.0068 10.0000 1.3363 8.1016 9.4696 17.5709按矩阵RB列中次序重新排列输入矩阵R中的值，则序数随机排列表各列间统计相关性便减小了，得到修正后的R，见表1。然后根据修正后的R，验证相关性是否减小，与求出修正后的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);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，得到相应的不用模拟次数的失效概率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.018,0.024);xlabel(Latin超立方采样模拟数);ylabel(失效概率Pf);title(可靠性计算结果)2） 例题2由题意知，若对非线性结构进行地震可靠性分析，首先需分别建立结构模型集和地震时程集。然后可利用Latin超立方采样技术将这些地震时程和结构模型匹配成81个地震结构系统。1、结构模型集：gsy12%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.002112%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.0025756%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.0032、 地震时程集：wggtf(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.3610f(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)3346p0.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)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)3709p0.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函数，function R=sample(n,k)%产生一个n行k列的随机抽样矩阵%k:输入变量数(维数）%n:每个输入变量抽取的样本数（每一维的采样数）R=zeros(n,k);for i=1:k; R(:,i)=randperm(n);%randperm(n):产生正整数1,2,3，.n的随机排列end令n=81,k=2生成随机结构模型集与地震时程集的随机组合，见下表：组合编号结构模型地震时程组合编号结构模型地震时程组合编号结构模型地震时程1754228693155557423742910235654343393930116575627430333116658353057614324945594622625593338296053187214134482061434486679357932624591758367869634064106236372010643375115015388576567812647239677666651113365405936722814549411568687171552404251246977541616374379701817176325442761715867181213453653727462192351464416737273204526476850742463213146482635753252228047495771762802334705014437781224296512881787356256038527019796121269285342180137727191254474881415512 / 12