Latin超立方抽样学习报告牛亚运.docx
- 文档编号:27942323
- 上传时间:2023-07-06
- 格式:DOCX
- 页数:16
- 大小:26.54KB
Latin超立方抽样学习报告牛亚运.docx
《Latin超立方抽样学习报告牛亚运.docx》由会员分享,可在线阅读,更多相关《Latin超立方抽样学习报告牛亚运.docx(16页珍藏版)》请在冰豆网上搜索。
Latin超立方抽样学习报告牛亚运
Latin超立方采样技术及其在结构可靠性分析中的应用
2、拉丁超立方采样
通过第一部分前言对拉丁超立方抽样概念的理解,结合第二部分对其原理的介绍,用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)'/n;
%rand(1,n):
产生n个值介于0.0到1.0的随机数
%randperm(n):
产生正整数1,2,3,...n的随机排列
end
%得到拉丁超立方抽样矩阵后,根据每个变量的分布函数,根据Xnk=f-1(Un)之间的关系,由每个变量对应的抽样结果Un反算出对该变量的真实抽样点Xk
3、统计相关的减小方程
1、Latin超立方抽样可能随机的引进了一定的统计相关,所以用文中采用的Spearman系数法来减小统计相关性。
首先,根据Spearman相关系数的计算公式,用matlab编辑程序如下:
functioncoeff=Spearman(X,Y)
%本函数用于实现Spearman相关系数的计算操作
%X:
输入的数值序列
%Y:
输入的数值序列
%coeff:
两个输入数值序列X,Y的相关系数
iflength(X)~=length(Y)
error('两个数值数列的维数不相等');
end
N=length(X);%得到序列的长度
Xrank=zeros(1,N);%存储X中各元素的排行
Yrank=zeros(1,N);%存储Y中各元素的排行
%计算Xrank中的各个值
fori=1:
N
count=1;
forj=1:
N
ifX(i) count=count+1; end end Xrank(i)=count; end %计算Yrank中的各个值 fori=1: N count=1; forj=1: N ifY(i) count=count+1; end end Yrank(i)=count; end %利用X,Y的序数排列计算Spearman相关系数 A=6*sum((Xrank-Yrank).^2); B=N*(N-1)*(N+1); coeff=1-A/B; end 2、仿照文中例子做一个K=5个输入变量和N=10个模拟的算例,以验证上述程序的正确性及Spearman系数对统计相关性的减小作用。 3、首先用sample.m函数生成一个K=5个输入变量和N=10个模拟的秩数随机排列表,见表1: functionR=sample(n,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个模拟的秩数随机排列表 未修正表 模拟变量 12345 166225 23110102 377841 484957 559136 615588 728779 8436110 9910363 10102494 修正表 模拟变量 12345 166123 231782 3771031 4849410 559254 615496 728869 843518 9910675 101023107 矩阵各列间的统计相关由序相关矩阵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)); end end 表2表1秩数随机排列表的秩相关矩阵 未修正表 变量变量 12345 11.00000.0667-0.2121-0.0909-0.4909 20.06671.0000-0.5152-0.3455-0.0182 3-0.2121-0.51521.00000.3697-0.0909 4-0.0909-0.34550.36971.0000-0.2970 5-0.4909-0.0182-0.0909-0.29701.0000 修正表 变量变量 12345 11.00000.06670.0061-0.0061-0.0061 20.06671.0000-0.0061-0.1758-0.1152 30.0061-0.00611.0000-0.09090.1152 4-0.0061-0.1758-0.09091.00000.0424 5-0.0061-0.11520.11520.04241.0000 T是正定的对称矩阵,可用Choiesky分解将T分解为T=Q*QT 在matlab命令运行窗口运行Q=chol(T),得到Q: Q= 1.00000.0667-0.2121-0.0909-0.4909 00.9978-0.5021-0.34020.0146 000.83840.2142-0.2239 0000.9111-0.3168 00000.7799 修正后的RB=R*Q-1,在matlab中运行RB=R*inv(Q),得到RB: RB= 6.00005.61257.26513.180813.4605 3.00000.801813.16718.478111.6619 7.00006.547915.23513.950811.5447 8.00003.474414.84014.093619.8692 5.00008.68607.66015.233115.0029 1.00004.94439.17828.567916.9100 2.00007.884213.57767.633219.6500 4.00002.73949.80950.212818.1910 9.00009.421011.49808.296816.0068 10.00001.33638.10169.469617.5709 按矩阵RB列中次序重新排列输入矩阵R中的值,则序数随机排列表各列间统计相关性便减小了,得到修正后的R,见表1。 然后根据修正后的R,验证相关性是否减小,及求出修正后的T,见表2. 结论: 由表2易得,修正后的各列间的统计相关性明显减小,用Spearman系数法可有效减小拉丁超立方抽样随机引进的统计相关。 4、在结构可靠性分析中的应用 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、结构模型集: ζ αg αs γy 1 2% 0.1 0.01 0.002 2 2% 0.1 0.01 0.0025 3 2% 0.1 0.01 0.003 4 2% 0.1 0.03 0.002 5 2% 0.1 0.03 0.0025 6 2% 0.1 0.03 0.003 7 2% 0.1 0.05 0.002 8 2% 0.1 0.05 0.0025 9 2% 0.1 0.05 0.003 10 2% 0.17 0.01 0.002 11 2% 0.17 0.01 0.0025 12 2% 0.17 0.01 0.003 13 2% 0.17 0.03 0.002 14 2% 0.17 0.03 0.0025 15 2% 0.17 0.03 0.003 16 2% 0.17 0.05 0.002 17 2% 0.17 0.05 0.0025 18 2% 0.17 0.05 0.003 19 2% 0.25 0.01 0.002 20 2% 0.25 0.01 0.0025 21 2% 0.25 0.01 0.003 22 2% 0.25 0.03 0.002 23 2% 0.25 0.03 0.0025 24 2% 0.25 0.03 0.003 25 2% 0.25 0.05 0.002 26 2% 0.25 0.05 0.0025 27 2% 0.25 0.05 0.003 28 4% 0.1 0.01 0.002 29 4% 0.1 0.01 0.0025 30 4% 0.1 0.01 0.003 31 4% 0.1 0.03 0.002 32 4% 0.1 0.03 0.0025 33 4% 0.1 0.03 0.003 34 4% 0.1 0.05 0.002 35 4% 0.1 0.05 0.0025 36 4% 0.1 0.05 0.003 37 4% 0.17 0.01 0.002 38 4% 0.17 0.01 0.0025 39 4% 0.17 0.01 0.003 40 4% 0.17 0.03 0.002 41 4% 0.17 0.03 0.0025 42 4% 0.17 0.03 0.003 43 4% 0.17 0.05 0.002 44 4% 0.17 0.05 0.0025 45 4% 0.17 0.05 0.003 46 4% 0.25 0.01 0.002 47 4% 0.25 0.01 0.0025 48 4% 0.25 0.01 0.003 49 4% 0.25 0.03 0.002 50 4% 0.25 0.03 0.0025 51 4% 0.25 0.03 0.003 52 4% 0.25 0.05 0.002 53 4% 0.25 0.05 0.0025 54 4% 0.25 0.05 0.003 55 6% 0.1 0.01 0.002 56 6% 0.1 0.01 0.0025 57 6% 0.1 0.01 0.003 58 6% 0.1 0.03 0.002 59 6% 0.1 0.03 0.0025 60 6% 0.1 0.03 0.003 61 6% 0.1 0.05 0.002 62 6% 0.1 0.05 0.0025 63 6% 0.1 0.05 0.003 64 6% 0.17 0.01 0.002 65 6% 0.17 0.01 0.0025 66 6% 0.17 0.01 0.003 67 6% 0.17 0.03 0.002 68 6% 0.17 0.03 0.0025 69 6% 0.17 0.03 0.003 70 6% 0.17 0.05 0.002 71 6% 0.17 0.05 0.0025 72 6% 0.17 0.05 0.003 73 6% 0.25 0.01 0.002 74 6% 0.25 0.01 0.0025 75 6% 0.25 0.01 0.003 76 6% 0.25 0.03 0.002 77 6% 0.25 0.03 0.0025 78 6% 0.25 0.03 0.003 79 6% 0.25 0.05 0.002 80 6% 0.25 0.05 0.0025 81 6% 0.25 0.05 0.003 2、地震时程集: ωg ζg t f(t) 1 3 0.6 5 f(t)1 2 3 0.6 5 f(t)2 3 3 0.6 5 f(t)3 4 3 0.6 10 f(t)1 5 3 0.6 10 f(t)2 6 3 0.6 10 f(t)3 7 3 0.6 1.5 f(t)1 8 3 0.6 1.5 f(t)2 9 3 0.6 1.5 f(t)3 10 3 0.36 5 f(t)1 11 3 0.36 5 f(t)2 12 3 0.36 5 f(t)3 13 3 0.36 10 f(t)1 14 3 0.36 10 f(t)2 15 3 0.36 10 f(t)3 16 3 0.36 1.5 f(t)1 17 3 0.36 1.5 f(t)2 18 3 0.36 1.5 f(t)3 19 3 0.84 5 f(t)1 20 3 0.84 5 f(t)2 21 3 0.84 5 f(t)3 22 3 0.84 10 f(t)1 23 3 0.84 10 f(t)2 24 3 0.84 10 f(t)3 25 3 0.84 1.5 f(t)1 26 3 0.84 1.5 f(t)2 27 3 0.84 1.5 f(t)3 28 6 0.6 5 f(t)1 29 6 0.6 5 f(t)2 30 6 0.6 5 f(t)3 31 6 0.6 10 f(t)1 32 6 0.6 10 f(t)2 33 6 0.6 10 f(t)3 34 6 0.6 1.5 f(t)1 35 6 0.6 1.5 f(t)2 36 6 0.6 1.5 f(t)3 37 6 0.36 5 f(t)1 38 6 0.36 5 f(t)2 39 6 0.36 5 f(t)3 40 6 0.36 10 f(t)1 41 6 0.36 10 f(t)2 42 6 0.36 10 f(t)3 43 6 0.36 1.5 f(t)1 44 6 0.36 1.5 f(t)2 45 6 0.36 1.5 f(t)3 46 6 0.84 5 f(t)1 47 6 0.84 5 f(t)2 48 6 0.84 5 f(t)3 49 6 0.84 10 f(t)1 50 6 0.84 10 f(t)2 51 6 0.84 10 f(t)3 52 6 0.84 1.5 f(t)1 53 6 0.84 1.5 f(t)2 54 6 0.84 1.5 f(t)3 55 9 0.6 5 f(t)1 56 9 0.6 5 f(t)2 57 9 0.6 5 f(t)3 58 9 0.6 10 f(t)1 59 9 0.6 10 f(t)2 60 9 0.6 10 f(t)3 61 9 0.6 1.5 f(t)1 62 9 0.6 1.5 f(t)2 63 9 0.6 1.5 f(t)3 64 9 0.36 5 f(t)1 65 9 0.36 5 f(t)2 66 9 0.36 5 f(t)3 67 9 0.36 10 f(t)1 68 9 0.36 10 f(t)2 69 9 0.36 10 f(t)3 70 9 0.36 1.5 f(t)1 71 9 0.36 1.5 f(t)2 72 9 0.36 1.5 f(t)3 73 9 0.84 5 f(t)1 74 9 0.84 5 f(t)2 75 9 0.84 5 f(t)3 76 9 0.84 10 f(t)1 77 9 0.84 10 f(t)2 78 9 0.84 10 f(t)3 79 9 0.84 1.5 f(t)1 80 9 0.84 1.5 f(t)2 81 9 0.84 1.5 f(t)3 在matlab中调用sample.m函数, functionR=sample(n,k) %产生一个n行k列的随机抽样矩阵 %k: 输入变量数(维数) %n: 每个输入变量抽取的样本数(每一维的采样数) R=zeros(n,k); fori=1: k; R(: i)=randperm(n)'; %randperm(n): 产生正整数1,2,3,...n的随机排列 end 令n=81,k=2生成随机结构模型集与地震时程集的随机组合,见下表: 组合编号 结构模型 地震时程 组合编号 结构模型 地震时程 组合编号 结构模型 地震时程 1 75 42 28 69 31 55 55 74 2 37 4 29 10 23 56 54 34 3 39 39 30 11 6 57 56 27 4 30 33 31 1 66 58 35 30 5 76 14 32 49 45 59 46 22 6 25 59 33 38 29 60 53 18 7 21 41 34 48 20 61 43 44 8 66 79 35 79 32 62 4 5 9 17 58 36 78 69 63 40 64 10 62 36 37 20 10 64 33 75 11 50 15 38 8 57 65 6 78 12 64 72 39 67 76 66 65 11 13 3 65 40 59 3 67 22 8 14 5 49 41 15 68 68 71 7 15 52 40 42 51 24 69 77 54 16 16 37 43 7 9 70 18 17 17 63 25 44 27 61 71 58 67 18 12 13 45 36 53 72 74 62 19 23 51 46 44 16 73 72 73 20 45 26 47 68 50 74 24 63 21 31 46 48 26 35 75 32 52 22 80 47 49 57 71 76 2 80 23 34 70 50 14 43 77 81 2 24 29 6 51 28 81 78 73 56 25 60 38 52 70 19 79 61 21 26 9 28 53 42 1 80 13 77 27 19 12 54 47 48 81 41 55
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Latin 立方 抽样 学习 报告 亚运
![提示](https://static.bdocx.com/images/bang_tan.gif)