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

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