结构可靠度分析与设计的编程实践

结构可靠度分析与设计的编程实践

1033002 班 1103300223 江莹

摘 要：基于《荷载与结构设计方法》中讲授的结构可靠性分析与设计的基本原理，对课

程中所给出的例题利用Matlab 软件编制了相应的计算机程序。通过此次编程实践，加深了自己对结构可靠性分析与设计的认识了理解。

1.引言

不确定性是自然界中普遍存在的一种客观现象，工程设计中的不确定性有多种不同的形式，人们认识最早、目前得到广泛应用的是随机性，人们用概率的方法研究。结构可靠度方法则是结构可靠性设计方法中的重要一项。结构的可靠度是结构在规定时间内，规定条件下结构能够完成预定功能的概率。从简单到复杂或精确程度的不同，先后提出的可靠度计算方法有MVFOSM （一次二阶矩法）、AFOSM （改进的一次二阶矩法）、RF 、MCS （蒙特卡洛数值模拟法）等方法。本文根据《荷载与结构设计方法》中讲授的结构可靠度分析与设计的基本原理，对课程中所给出的例题利用Matlab 软件编制了相应的计算机程序。

2.结构可靠度分析的基本原理

当用概率描述结构的可靠性时，就需要根据结构中基本随机变形变量或综合随机变量的概率分布进行计算。在实际工程中，结构的功能函数往往是由多个随机变量组成的非线性函数，而且这些随机变量并不都服从正态分布或对数正态分布，因此不能直接采用相应的公式计算可靠指标，而需要作出某些近似简化后进行计算。下面本文将介绍分析结构可靠度的几种常用方法： 2.1.AFOSM 法（改进的一次二阶矩法）

1974年Hasofer 和Lind 科学地对可靠指标进行了定义，引入了验算点的概念，使得一阶二次矩模式有了进一步的发展，由于分析中要迭代求解验算点，验算点是可靠分析中的一个关键点，所以人们又将这种方法称为验算点法。

随机变量服从正态分布的情形。、 2.1.1.功能函数为线性函数

01n

i i

i Z a a X ==+∑

式中：

01

,a a (i=0,1,2,…n)为常数。

为进一步在标准正太坐标系中研究可靠指标的几何意义，按下式将随机变量

(1,2,3,,)

i X i n =???变换为标准正态随机变量

(1,2,3,,)

i Y i n =???

(1,2,,)

i

i

i X i X

X Y i n μσ-=

=???

则结构功能函数可由

(1,2,3,,)

i Y i n =???表示为：

123001

1

1

(,,,,)()i i i i n n n

n i X X i i X i X i

i i i Z g Y Y Y Y a a Y a a a Y μσμσ====???=++=++∑∑∑

从而功能函数的平均值和标准差为：

01

221

i

i

n

Z i X

i n

Z i X i a a a μμσσ

===+=

∑∑

可靠度指标为：

01

221

i

i

n

i X i Z n

Z

i X i a a a μμβσσ

==+==

∑∑

由上式计算的可靠指标与结构失效概率的对应关系：

()p

f

β=Φ-

2.1.2.功能函数为非线性函数

为计算可靠指标，将非线性功能函数在验算点处泰勒级数展开，则结构功能函数的一次展开式为：

*****1

2

3

1

(,,,,)()n

n

P

i i i i

g Z g x x x x X x X =?=???+-?∑

则Z 的平均值和标准差可分别近似表示为

*****1

2

3

1

2

1

(,,,,)()

'

(

)i i

n

Z n

P

X i i i n

Z P

X i i

g g x x x x x X g X μμσσ==?=???+-??=

?∑

∑

由此可求得可靠指标：

*****

1

2

3

1

2

1

(,,,,)()

(

)i i

n

n

P X i i i

Z

n

Z

P

X i i

g g x x x x x X g

X μμβσσ==????+-?=

=??∑

∑

2.2.RF 法

当随机变量不服从正态分布时，则需要将其进行当量正态化将随机变量i X 等效为正态随机变量'i X ，当量正态化的条件是，在验算点处使非正态随机变量i X 的概率分布函数值与当量正态随机变量'i X 的概率分布函数值相等，i X 的概率密度函数值与当量正态随机变

量'i X 的概率密度函数值相等，即

*'

**'

*'

**'

'

()()()

1

()(

)()

i i i i

i i i i

i

i X X i

X i X i

X X i X i X X x F x F x x f x f x μσμ?σσ-=Φ=-=

=

则：

*1*''

1*'*[()]{[()]}

()

i

i

i

i

i

i X i X i X

X i X X i x F x F x f x μσ?σ--=-ΦΦ=

当量正态化后：

*****1

2

3

'1

2

'1

'

'2

'1

*''(,,,,)()

()cos (1,2,3,,)(

)(1,2,3,,)

i i

i

i

i

j

i i i n

n

P

X i i i

Z

n

Z

P X i i

P

X

i

X X n

P

X j j

i X X X g g x x x x x X g

X

g X i n g X x i n μμβσσσαθσμασβ===????+-?=

=????==-

=?????=+=???∑

∑∑

上述三式构成非正态随机变量情况下可靠指标的迭代计算公式。

3.结构可靠性设计的基本原理

结构设计首先要保证结构的安全性，其次是保证结构的适用性、耐久性和偶然作用下的整体稳定性。为实现这些目的，就要对结构进行合理的设计。

3.1结构可靠度校准

如果按照结构设计规范规定的设计表达式进行可靠度计算，按规范设计的结构或结构构件的可靠度指标，在这种情况下反映的是一种结构或一类结构构件的可靠指标，代表了设计规范的可靠水平，这一过程称为结构可靠度校准。 结构可靠度校准的步骤： (1).假设一个荷载效应比ρ； (2).确定构件抗力的特征值k R ：

(1)k Gk R K S ρ=+

(3).确定基本变量的均值和标准差：

均值： ,,G G Q Q R R k S S Sk S S Qk R S S μλμλμλ=== 标准差： ,,G G G Q Q Q R R R S S S S S S σδμσδμσδμ===

(4). 确定极限状态方程：

0G Q R S S --=

(5).利用RF 法求解可靠度指标β；

(6).变化荷载效应比，计算不同可靠指标的均值

3.2分项系数的校准

如果给定了结构或构件的设计表达式且已知变量的概率分布和统计参数，即可求得按设计表达式设计的结构或构件的可靠指标，可靠指标至于作用效应标准值的比值有关，而与作用效应标准值本身的大小无关。这样可以通过先设定作用分项系数G γ、Q γ抗力分项系数

R γ和作用组合系数ψ，然后计算相应的可靠指标，并与目标可靠指标进行比较，从而确定

一组使计算的可靠指标与目标可靠指标最接近的系数值，即确定的G γ、Q γ、R γ和ψ使下式最小

2

()i i T H w ββ=-∑

下面是用迭代算法求分项系数的过程：

（1）确定极限状态方程和设计方程的表达式。

①确定基本变量的概率分布类型和近似参数

②最多有两个未知的均值需要求解。一个是抗力的均值R μ，另一个是相应的荷载效

应均值i S μ，荷载效应比率用于建立荷载与抗力均值之间的联系

③第一次迭代，利用极限状态方程 获得两个未知变量均值的关系。 （2）利用上述均值获得初始的设计点*

i x 。

（3）对于非正态分布对应的设计点取值*

i x ，利用当量正态化公式确定等效的正态均值i e

X μ和正态标准差i e

X σ。

（4）计算n 个变量的方向余弦i α

*

*

2

1()i

i

X

P i

i n

X P

i i

g

X g X σασ=?-

?=??∑

（5）计算n 个变量设计点的坐标值*

i x

*[](1,2,3,,)i i i X i X x i n μαβσ=+=???

（6）通过求解极限状态方程确定两个未知均值新的关系

****123(,,,)0n g x x x x ???=

（7）重复3-6步直到{}i α收敛 （8）当收敛条件满足，计算分项系数

*/i

X i ri x X γ=

4.例题的程序运行结果展示

4.1Homework2.1的结果展示

The example 5.4

迭代次数 E *

I * W * β

1 2e7 8e-4 10

2 9.690e+06

5.680e-04 10.0256 2.5782 3 5.369e+06 5.754e-04 10.0512 3.2896 4 4.541e+06

6.702e-04 10.0534 3.2135 5 4.419e+06

7.050e-04 10.0462 3.1821 6 4.380e+06 7.119e-04 10.0440 3.1806 7 4.371e+06

7.134e-04

10.0436

3.1805

The problem 5.3

迭代次数

1 2 β 3.1537 3.1537 1α

-0.9997 -0.9997 2α

0.0232

0.0232

4.2.Homework2.2的结果展示

The example 3.5

迭代次数

R * Q * G * β

1 250 85

50 2 167.4692 116.6522 50.8169 3.9046 3 204.2590 153.5989 50.6601 3.4566 4 213.8648 163.4033 50.4615 3.3373 5 214.6858 164.2532 50.4326 3.3344 6 214.7308

164.3004 50.4304

3.3344

The example 5.9

迭代次数R*Q*β

1 200 100

2 123.0790 123.0790 4.3903

3 156.9136 156.9136 3.9476

4 167.5417 167.5417 3.7674

5 168.4619 168.4619 3.7633

6 168.499

7 168.4997 3.7633

The example 5.10

迭代次数R*Q*β

1 150 120

2 147.5024 139.3502 3.3600

3 147.5582 139.3925 3.3640

4 147.5574 139.3918 3.3640

The example 5.11

迭代次数Z*Fy*M*β

1 100 40 2000

2 93.9069 24.6440 2314.2419 4.4794

3 95.711

4 30.2853 2898.6496 4.3254

4 96.6872 32.5887 3150.9067 4.0376

5 96.8361 32.904

6 3186.3504 4.0223

6 96.8516 32.9273 3189.0605 4.0221

7 96.8529 32.9288 3189.2505 4.0221

The problem 5.4

迭代次数Y*X*β

1 12.5 24

2 9.4185 28.2554 2.3823

3 10.5011 31.5033 2.2468

4 10.7061 32.1183 2.2139

5 10.725

6 32.1769 2.2132

6 10.7272 32.1816 2.2132

4.3Homework2.3的结果展示

The example 5.11

p=3.0000e-005 β= 4.0128

f

The problem 5.4

p=0.0139 β=2.1996

f

4.4 Homework3的结果展示

The example 6.3

k R =225.4844 gk A =278.3758

The example 6.7

oR γ=1.2746 oQ γ=1.2478

5.结论与讨论

结构可靠度设计与分析方法，从简单到复杂或精确程度的不同，先后提出的可靠度计算方法有MVFOSM 、AFOSM 、RF 、MCS 等方法。通过此次编程实践，加深了自己对结构可靠度分析与设计的认识了理解。

6.参考文献

1. 中华人民共和国国家标准. 工程结构可靠性设计统一标准（GB50153—2008），北京：中

国建筑工业出版社，2009. 2. 中华人民共和国国家标准. 建筑结构可靠度设计统一标准（GB50068—2001），北京：中

国建筑工业出版社，2001. 3. 张明著. 结构可靠度分析 — 方法与程序. 北京：科学出版社，2009年.

4. 李国强等编著. 工程结构荷载与可靠度设计原理（第三版）. 北京：中国建筑工业出版

社，2005. 5. 贡金鑫，魏巍巍. 工程结构可靠性设计原理. 北京：机械工业出版社，2007.

6. Nowak, A. S., and K. R. Collins (2000). Reliability of Structures. McGraw-Hill, New York,

NY . 张川导读，重庆：重庆大学出版社，2005年3月

附录

附件1：Homework 2

%AFOSM homework 2_1_1 % limit state g=EI-310.5w

% E ,I and w are all normal variables

clear;clc

mu_w=10;deta_w=0.04;mu_E=2e7;deta_E=0.25;mu_I=8e-4;deta_I=1.5/8; sigma_E=mu_E*deta_E;sigma_I=mu_I*deta_I;sigma_w=mu_w*deta_w;

Estar=mu_E;Istar=mu_I;wstar=mu_w;

i=2;A(1:100,1:4)=zeros(1:100,1:4);A(1,4)=1; while abs(A(i,4)-A(i-1,4))>=1e-4

a_E=Istar*sigma_E; a_I=Estar*sigma_I; a_w=-310.5*sigma_w; sqrt1=sqrt(a_E^2+a_I^2+a_w^2);

beta=(Estar*Istar-310.5*wstar+Istar*(mu_E-Estar)+Estar*(mu_I-Istar)-310.5*(mu_w-wstar))/sqrt 1;

alpha_E=-a_E/sqrt1; alpha_I=-a_I/sqrt1;alpha_w=-a_w/sqrt1;

Estar=mu_E+alpha_E*beta*sigma_E; Istar=mu_I+alpha_I*beta*sigma_I; wstar=mu_w+alpha_w*beta*sigma_w;

c=[Estar Istar wstar beta]; A(i+1,1:4)=c;

i=i+1

end

s=xlswrite('2_1_1',A(2:i,:));

beta

Estar

Istar

wstar

*******************************************************************************

******************************************************************************* %RF 2_2_1

%limit state g=R-G-Q

%R lognorm G normal Q extrem I

clear;clc

mu_G=50;mu_Q=85;mu_R=250;sigma_G=2.5;sigma_Q=17;sigma_R=25;

Gstar=mu_G;Rstar=mu_R;Qstar=mu_Q;

i=2;A(1:100,1:4)=zeros(1:100,1:4);A(1,4)=1;

while abs(A(i,4)-A(i-1,4))>=1e-4

%R服从对数正态分布

mu_lnR=log(mu_R/sqrt(1+(sigma_R/mu_R)^2)); mu_R1=Rstar*(1-log(Rstar)+mu_lnR);

sigma_lnR=sqrt(log(1+(sigma_R/mu_R)^2)); sigma_R1=Rstar*sigma_lnR;

%Q服从极值Ⅰ型分布

alpha=pi/(sqrt(6)*sigma_Q); u=mu_Q-0.5772/alpha;

F_Q=exp(-exp(-alpha*(Qstar-u)));f_Q=alpha*exp(-alpha*(Qstar-u)-exp(-alpha*(Qstar-u)));

sigma_Q1=normpdf(norminv(F_Q))/f_Q;

mu_Q1=Qstar-norminv(F_Q)*sigma_Q1;

%G服从正态分布

mu_G1=mu_G;sigma_G1=sigma_G;

%正态化后极限状态方程变为g=R1-G1-Q1

beta=(Rstar-Gstar-Qstar+((mu_R1-Rstar)-(mu_G1-Gstar)-(mu_Q1-Qstar)))/sqrt(sigma_R1^2+sig

ma_G1^2+sigma_Q1^2);

alpha_R1=-sigma_R1/sqrt(sigma_R1^2+sigma_G1^2+sigma_Q1^2);

alpha_G1=sigma_G1/sqrt(sigma_R1^2+sigma_G1^2+sigma_Q1^2);

alpha_Q1=sigma_Q1/sqrt(sigma_R1^2+sigma_G1^2+sigma_Q1^2);

Rstar=mu_R1+alpha_R1*sigma_R1*beta;

Qstar=mu_Q1+alpha_Q1*sigma_Q1*beta;

Gstar=mu_G1+alpha_G1*sigma_G1*beta;

c=[Rstar Qstar Gstar beta]; A(i+1,1:4)=c;

i=i+1;

end

s=xlswrite('2_2_1',A（2：i,:）);

beta

Rstar

Gstar

Qstar

*******************************************************************************

%RF 2_2_21

%limit state g=R-Q

%R lognorm Q extrem I

clear;clc

mu_Q=100;mu_R=200;sigma_Q=12;sigma_R=20;

Rstar=mu_R;Qstar=mu_Q;

i=2;A(1:100,1:3)=zeros(1:100,1:3);A(1,3)=1;

while abs(A(i,3)-A(i-1,3))>=1e-4

%R服从对数正态分布

mu_lnR=log(mu_R/sqrt(1+(sigma_R/mu_R)^2)); mu_R1=Rstar*(1-log(Rstar)+mu_lnR);

sigma_lnR=sqrt(log(1+(sigma_R/mu_R)^2)); sigma_R1=Rstar*sigma_lnR;

%Q服从极值Ⅰ型分布

alpha=pi/(sqrt(6)*sigma_Q); u=mu_Q-0.5772/alpha;

F_Q=exp(-exp(-alpha*(Qstar-u)));f_Q=alpha*exp(-alpha*(Qstar-u)-exp(-alpha*(Qstar-u)));

sigma_Q1=normpdf(norminv(F_Q))/f_Q;

mu_Q1=Qstar-norminv(F_Q)*sigma_Q1;

%正态化后极限状态方程变为g=R1-Q1

beta=(Rstar-Qstar+((mu_R1-Rstar)-(mu_Q1-Qstar)))/sqrt(sigma_R1^2+sigma_Q1^2);

alpha_R1=-sigma_R1/sqrt(sigma_R1^2+sigma_Q1^2);

alpha_Q1=sigma_Q1/sqrt(sigma_R1^2+sigma_Q1^2);

Rstar=mu_R1+alpha_R1*sigma_R1*beta;

Qstar=mu_Q1+alpha_Q1*sigma_Q1*beta;

c=[Rstar Qstar beta]; A(i+1,1:3)=c;

i=i+1;

end

s=xlswrite('2_2_21',A(2:i,:));

Rstar

Qstar

beta

******************************************************************************* %RF 2_2_22

%limit state g=R-Q

%R lognorm Q normal

clear;clc

mu_R=input('mu_R=');deta_R=input('deta_R=');

mu_Q=input('mu_Q=');deta_Q=input('deta_Q=');

sigma_R=mu_R*deta_R;sigma_Q=mu_Q*deta_Q;

Rstar=mu_R;Qstar=mu_Q;

i=2;A(1:100,1:3)=zeros(1:100,1:3);A(1,3)=1;

while abs(A(i,3)-A(i-1,3))>=1e-4

%R服从对数正态分布

mu_lnR=log(mu_R/sqrt(1+(sigma_R/mu_R)^2)); mu_R1=Rstar*(1-log(Rstar)+mu_lnR);

sigma_lnR=sqrt(log(1+(sigma_R/mu_R)^2)); sigma_R1=Rstar*sigma_lnR;

%Q服从正态分布

mu_Q1=mu_Q;sigma_Q1=sigma_Q;

%正态化后极限状态方程变为g=R1-Q1

beta=(Rstar-Qstar+((mu_R1-Rstar)-(mu_Q1-Qstar)))/sqrt(sigma_R1^2+sigma_Q1^2);

alpha_R1=-sigma_R1/sqrt(sigma_R1^2+sigma_Q1^2+sigma_Q1^2);

alpha_Q1=sigma_Q1/sqrt(sigma_R1^2+sigma_Q1^2+sigma_Q1^2);

Rstar=mu_R1+alpha_R1*sigma_R1*beta;

Qstar=mu_Q1+alpha_Q1*sigma_Q1*beta;

c=[Rstar Qstar beta]; A(i+1,1:3)=c;

i=i+1;

end

s=xlswrite('2_2_22',A(2:i,:));

beta

Rstar

Qstar

******************************************************************************* %RF 2_2_31

%g=Z*Fy-M

% Z are normal variables; Fy are lognormal variables; M are extrem type I variables

clear;clc

mu_Z=100; deta_Z=0.04;

mu_Fy=40;deta_Fy=0.1;

mu_M=2000;deta_M=0.1;

sigma_Z=mu_Z*deta_Z;sigma_Fy=mu_Fy*deta_Fy;sigma_M=mu_M*deta_M;

sigma_Z1=sigma_Z;sigma_Fy1=sigma_Fy;sigma_M1=sigma_M;

beta=0; alpha_Z1=0; alpha_Fy1=0; alpha_M1=0;

Zstar=mu_Z+alpha_Z1*beta*sigma_Z1;

Fystar=mu_Fy+alpha_Fy1*beta*sigma_Fy1;

Mstar=mu_M+alpha_M1*beta*sigma_M1;

i=2;A(1:100,1:4)=zeros(1:100,1:4);A(1,4)=1;

while abs(A(i,4)-A(i-1,4))>1e-4

%Z服从正态分布

mu_Z1=mu_Z;sigma_Z1=sigma_Z;

%Fy服从对数正态分布

mu_lnFy=log(mu_Fy/sqrt(1+(sigma_Fy/mu_Fy)^2));

mu_Fy1=Fystar*(1-log(Fystar)+mu_lnFy);

sigma_lnFy=sqrt(log(1+(sigma_Fy/mu_Fy)^2)); sigma_Fy1=Fystar*sigma_lnFy;

%M服从极值Ⅰ型分布

alpha=pi/(sqrt(6)*sigma_M); u=mu_M-0.5772/alpha;

F_M=exp(-exp(-alpha*(Mstar-u)));f_M=alpha*exp(-alpha*(Mstar-u)-exp(-alpha*(Mstar-u)));

sigma_M1=normpdf(norminv(F_M))/f_M;

mu_M1=Mstar-norminv(F_M)*sigma_M1;

sqrt1=sqrt((Fystar*sigma_Z1)^2+(Zstar*sigma_Fy1)^2+sigma_M1^2);

alpha_Z1=-sigma_Z1*Fystar/sqrt1;

alpha_Fy1=-sigma_Fy1*Zstar/sqrt1;

alpha_M1=sigma_M1/sqrt1;

%正态化后极限状态方程变为g=Z1*Fy1*M1

b=solve('(mu_Z1+alpha_Z1*x*sigma_Z1)*(mu_Fy1+alpha_Fy1*x*sigma_Fy1)-(mu_M1+alpha_ M1*x*sigma_M1)=0');

c=eval(b);

beta=c(1);

Zstar=mu_Z1+alpha_Z1*beta*sigma_Z1;

Fystar=mu_Fy1+alpha_Fy1*beta*sigma_Fy1;

Mstar=mu_M1+alpha_M1*beta*sigma_M1;

c=[Zstar Fystar Mstar beta]; A(i+1,1:4)=c;

i=i+1;

end

s=xlswrite('2_2_31',A(1:i,:));

beta

Zstar

Fystar

Mstar

******************************************************************************* %RF 2_2_32

%limit state g=3Y-X

%Y lognormal; X extrem type I

clear;clc

mu_Y=12.5;deta_Y=0.125;

mu_X=24;deta_X=0.15;

sigma_Y=mu_Y*deta_Y;sigma_X=mu_X*deta_X;

sigma_Y1=sigma_Y;sigma_X1=sigma_X;

beta=0; alpha_Y1=0; alpha_X1=0;

Ystar=mu_Y+alpha_Y1*beta*sigma_Y1;

Xstar=mu_X+alpha_X1*beta*sigma_X1;

i=2;A(1:100,1:3)=zeros(1:100,1:3);A(1,3)=1;

while abs(A(i,3)-A(i-1,3))>=1e-4

%Y服从对数正态分布

mu_lnY=log(mu_Y/sqrt(1+(sigma_Y/mu_Y)^2)); mu_Y1=Ystar*(1-log(Ystar)+mu_lnY);

sigma_lnY=sqrt(log(1+(sigma_Y/mu_Y)^2)); sigma_Y1=Ystar*sigma_lnY;

%X服从极值Ⅰ型分布

alpha=pi/(sqrt(6)*sigma_X); u=mu_X-0.5772/alpha;

F_M=exp(-exp(-alpha*(Xstar-u)));f_M=alpha*exp(-alpha*(Xstar-u)-exp(-alpha*(Xstar-u)));

sigma_X1=normpdf(norminv(F_M))/f_M;

mu_M1=Xstar-norminv(F_M)*sigma_X1;

sqrt1=sqrt((3*sigma_Y1)^2+sigma_X1^2);

alpha_Y1=-sigma_Y1*3/sqrt1;

alpha_X1=sigma_X1/sqrt1;

%正态化后极限状态方程变为g=3Y1-X1

b=solve('3*(mu_Y1+alpha_Y1*x*sigma_Y1)-(mu_M1+alpha_X1*x*sigma_X1)=0');

c=eval(b);

beta=c(1);

Ystar=mu_Y1+alpha_Y1*beta*sigma_Y1;

Xstar=mu_M1+alpha_X1*beta*sigma_X1;

c=[Ystar Xstar beta]; A(i+1,1:3)=c;

i=i+1;

end

s=xlswrite('2_2_32',A(1:i,:));

beta

Ystar

Xstar

******************************************************************************* %2_31

Monte Carlo Simulation

%judgement:Fy*Z-M<0 (limit state g=Fy*Z-M)

%Fy are lognormal variables Z are normal variables M are extreme i variables

clear;clc

mu_Fy=40;deta_Fy=0.1;sigma_Fy=mu_Fy*deta_Fy;

mu_lnFy=log(mu_Fy/(sqrt(1+deta_Fy^2))); sigma_lnFy=sqrt(log(1+deta_Fy^2));

mu_Z=100;deta_Z=0.04;sigma_Z=mu_Z*deta_Z;

mu_M=2000;deta_M=0.1;sigma_M=mu_M*deta_M;

alpha=pi/(sqrt(6)*sigma_M); u=mu_M-0.5772/alpha;

Fy=lognrnd(mu_lnFy,sigma_lnFy,10e5,1);

Z=normrnd(mu_Z,sigma_Z,10e5,1);

p=rand(10e5,1);

M=u-log(-log(p))/alpha;

a=find(Z.*Fy-M<=0);

N=length(a);

pf=N/10e5

beta=-norminv(pf)

******************************************************************************* %2_32

Monte Carlo Simulation homework 2_32

%judgement3Y-X<0 (limit state g=3Y-X)

%Y are lognormal variables X are extreme i variables

clear; clc;

mu_Y=12.5;deta_Y=0.125;

mu_X=24;deta_X=0.15;

sigma_Y=mu_Y*deta_Y;sigma_X=mu_X*deta_X;

alpha=pi/(sqrt(6)*sigma_X); u=mu_X-0.5772/alpha;

p=rand(10e5,1);

X=u-log(-log(p))/alpha;

mu_lnY=log(mu_Y/(sqrt(1+deta_Y^2))); sigma_lnY=sqrt(log(1+deta_Y^2));

Y=lognrnd(mu_lnY,sigma_lnY,10e5,1);

a=find(3*Y-X<=0);

N=length(a);

pf=N/10e5

beta=-norminv(pf)

*******************************************************************************

附录2 Homework 3

%homework 3_1

%reliability-based partial safety factor design

%R are lognorm variables Q are extreme I viariables G are normal variables

clear; clc;

mu_G=50;mu_Q=70;mu_R=120;sigma_G=3.5;sigma_Q=20;sigma_R=120*0.17;

Gstar1=mu_G;Rstar1=mu_R;Qstar1=mu_Q;Qstar=0;Rstar=0;Gstar=0;

beta=3.2;

while abs([Gstar,Qstar,Rstar]-[Gstar1,Qstar1,Rstar1])>=1e-3

Gstar=Gstar1;Qstar=Qstar1;Rstar=Rstar1;

%R服从对数正态分布

mu_lnR=log(mu_R/sqrt(1+(sigma_R/mu_R)^2)); mu_R1=Rstar*(1-log(Rstar)+mu_lnR);

sigma_lnR=sqrt(log(1+(sigma_R/mu_R)^2)); sigma_R1=Rstar*sigma_lnR;

%Q服从极值Ⅰ型分布

alpha=pi/(sqrt(6)*sigma_Q); u=mu_Q-0.5772/alpha;

F_Q=exp(-exp(-alpha*(Qstar-u)));f_Q=alpha*exp(-alpha*(Qstar-u)-exp(-alpha*(Qstar-u)));

sigma_Q1=normpdf(norminv(F_Q))/f_Q;

mu_Q1=Qstar-norminv(F_Q)*sigma_Q1;

%G服从正态分布

mu_G1=mu_G;sigma_G1=sigma_G;

%正态化后极限状态方程变为g=R1-G1-Q1

alpha_R1=-sigma_R1/sqrt(sigma_R1^2+sigma_G1^2+sigma_Q1^2);

alpha_G1=sigma_G1/sqrt(sigma_R1^2+sigma_G1^2+sigma_Q1^2);

alpha_Q1=sigma_Q1/sqrt(sigma_R1^2+sigma_G1^2+sigma_Q1^2);

Qstar1=mu_Q1+alpha_Q1*sigma_Q1*beta;

Gstar1=mu_G1+alpha_G1*sigma_G1*beta;

Rstar1=mu_R1+alpha_R1*sigma_R1*beta;

end

Rstar=Gstar+Qstar;

mu_lnR=log(mu_R/sqrt(1+(sigma_R/mu_R)^2)); mu_R1=Rstar*(1-log(Rstar)+mu_lnR);

sigma_lnR=sqrt(log(1+(sigma_R/mu_R)^2)); sigma_R1=Rstar*sigma_lnR;

alpha_R1=-sigma_R1/sqrt(sigma_R1^2+sigma_G1^2+sigma_Q1^2);

mu_R1=Rstar-sigma_R1*beta*alpha_R1;

mu_R=sqrt(1+(sigma_R/mu_R)^2)*exp(mu_R1/Rstar-1+log(Rstar));

R_K=mu_R/1.33

A_gK=R_K/((0.01/0.31*1/0.02+1)*0.31)

******************************************************************************* %RF3_2

%g=R-Q

%R、Q值未给出，不妨令R=Q=1

clear;clc

mu_Q=1;mu_R=1;deta_Q=0.12;deta_R=0.1; beta=3.0;

alpha_R=0;alpha_Q=0;alpha_R1=1;alpha_Q1=1;

while abs([alpha_R1 alpha_Q1]-[alpha_R alpha_Q])>=1e-3

alpha_R=alpha_R1;alpha_Q=alpha_Q1;

sigma_R=deta_R*mu_R;

sigma_Q=mu_Q*deta_Q;

alpha_R1=-sigma_R/sqrt(sigma_R^2+sigma_Q^2);

alpha_Q1=sigma_Q/sqrt(sigma_R^2+sigma_Q^2);

Qstar=mu_Q+alpha_Q1*sigma_Q*beta;

Rstar=mu_R+alpha_R1*sigma_R*beta;

mu_R=Qstar/Rstar*mu_Q;

end

Rstar;

Qstar;

gamma_0R=mu_R/Rstar

gamma_0Q=Qstar/mu_Q