基于MATLAB的可靠度指标与失效概率的计算

Tooto_ttt

已于 2022-04-20 09:54:39 修改

阅读量1.5w

点赞数 12

分类专栏：可靠性分析理论文章标签： matlab

于 2022-04-14 11:16:44 首次发布

本文链接：https://blog.csdn.net/m0_56146217/article/details/124143392

版权

可靠性分析理论专栏收录该内容

2 篇文章

订阅专栏

(5条消息) 基于Abaqus的边坡可靠度计算_步步为营！的博客-CSDN博客https://blog.csdn.net/qq_41710383/article/details/115006487?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522165026719616782089328511%2522%252C%2522scm%2522%253A%252220140713.130102334.pc%255Fblog.%2522%257D&request_id=165026719616782089328511&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~blog~first_rank_ecpm_v1~rank_v31_ecpm-6-115006487.nonecase&utm_term=%E5%8F%AF%E9%9D%A0%E5%BA%A6&spm=1018.2226.3001.4450可靠性理论是以概率论和数理统计为理论基础的，且概率论是理论分析和实践相结合的一种方法，其相比定值分析更有优势也更符合实际情况。

一、可靠性度量指标

可靠性度量指标有以下三种：（1）可靠概率（又称可靠度），用 $P_{r}$ 表示；结构在规定的时间内，在规定的条件下，完成预定功能的能力，称为结构可靠性。对于特定情况而言，可具体解释成在规定的参照时期内，结构将不会达到某一特定极限状态的概率，即所谓的结构可靠度。（2）失效概率，用 $P_{f}$ 表示；可靠概率和失效概率是互补的。（3）可靠指标，用 $\beta$ 表示，可靠指标越大，表示越安全。

失效概率 $P_{f}$ 与可靠指标 $\beta$ 的关系为：

二、极限状态功能函数（功能函数、极限状态方程）

不同的模型有不同的极限状态方程，在这一步神经网络能够更准确的对模型的极限状态功能函数进行拟合。如我们在考虑直立式防波堤时，主要考虑两个破坏形式，一是滑动破坏，一个是倾覆破坏。则其两种极限状态的极限状态方程可以确定下来。要注意的是，功能函数其单位都得是相同的（如都是长度、力、力矩等）。在一些常见工程的各类设计规范中，也会有提及。

三、基于蒙特卡罗抽样的Matlab程序（防波堤）

在实际工程的设计中，有很多影响结构可靠性的不确定因素，这些因素的不确定性最终是表现为随机变量的形式。例如在防波堤的设计中，水平波浪力和波浪浮托力实际上都是随机变量，其直接影响着结构的使用寿命。在可靠度计算过程中，蒙特卡罗（Monte Carlo，MC）法避开了结构可靠度中的数学困难，通过试验模型来找到最佳选择。而且更容易编程。蒙特卡罗法能够直接计算出失效概率，从而利用得出可靠指标。

第一步：通过Excel输入计算所需的数据，计算各参数的均值和标准差，可以先绘制出各随机参数的概率密度函数图形（绘制这一步可有可无）。

第二步：根据随机参数的分布类型产生相应类型的1组随机数矩阵（2个随机变量，即2行1列）（即水平波浪力和波浪浮托力各一个）（这里就体现了蒙特卡罗法）（利用MATLAB可以快速产生各种分布类型的随机数的功能，函数功能可以见Matlab教程。如对于正态分布函数来说，可以利用normrnd（mu,sigma,m,n）来产生指定均值和标准差的m*n阶正态分布随机数矩阵）

在产生第k组随机数的这一步，可以利用之前博文所提到的LHS方法先进行抽样，这样能够更好的去抽样excel中的数据，然后再将抽样好的数据中的第k组去代入功能函数。

第三步：再将产生的这一组随机数代入到结构已知特定的极限状态方程中，得出Z值。如Z值小于0，为失效区（Z=0为极限状态，Z>0为可靠区），此时计为一次失效。然后当抽样次数少于预定的值时，继续抽样计算Z值，根据循环去统计Z<0的次数，最后失效概率为 $P_{f}=\frac{n}{N}$ ，由于可靠度与其互补，可以得到可靠度 $P_{r}=1-P_{f}$ ，以及由上述表2-1关系得到可靠指标 $\beta$ 。

四、Gumbel分布

Gumbel分布能够预测极端情况。其CDF是可逆的，即能够从均匀分布的随机变量中生成Gumbel分布的随机变量。

(5条消息) Matlab中极值I型分布(Gumbel分布)的evcdf和evpdf及evinv函数_leo2351960的博客-CSDN博客_极值i型分布https://blog.csdn.net/leo2351960/article/details/39210807?utm_medium=distribute.pc_aggpage_search_result.none-task-blog-2~aggregatepage~first_rank_ecpm_v1~rank_v31_ecpm-5-39210807.pc_agg_new_rank&utm_term=%E6%9E%81%E5%80%BC1%E5%9E%8B%E5%88%86%E5%B8%83%E5%87%BD%E6%95%B0&spm=1000.2123.3001.4430

(5条消息) MATLAB实现各种概率密度函数（概率密度/分布/逆概率分布函数）_WW、forever的博客-CSDN博客_matlab概率密度分布图函数https://blog.csdn.net/qq_44246618/article/details/123301478?utm_medium=distribute.pc_aggpage_search_result.none-task-blog-2~aggregatepage~first_rank_ecpm_v1~rank_v31_ecpm-17-123301478.pc_agg_new_rank&utm_term=%E6%9E%81%E5%80%BC1%E5%9E%8B%E5%88%86%E5%B8%83%E5%87%BD%E6%95%B0&spm=1000.2123.3001.4430 里面包含了详细的关于各种分布的基础知识

五、基于Matlab的串并联可靠度计算

(5条消息) 2.18 串联系统与并联系统可靠度计算_锋戈的博客-CSDN博客_串并联系统可靠度计算例题https://blog.csdn.net/qq_30353463/article/details/117842059?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522165026719616782089328511%2522%252C%2522scm%2522%253A%252220140713.130102334.pc%255Fblog.%2522%257D&request_id=165026719616782089328511&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~blog~first_rank_ecpm_v1~rank_v31_ecpm-5-117842059.nonecase&utm_term=%E5%8F%AF%E9%9D%A0%E5%BA%A6&spm=1018.2226.3001.4450

六、响应面法

(5条消息) 结构可靠性分析中响应面方法的基本思想_朽木可燃的博客-CSDN博客_响应面法https://blog.csdn.net/SKY_CM/article/details/90611823?ops_request_misc=&request_id=&biz_id=102&utm_term=%E5%A4%B1%E6%95%88%E6%A6%82%E7%8E%87%E8%AE%A1%E7%AE%97&utm_medium=distribute.pc_search_result.none-task-blog-2~blog~sobaiduweb~default-6-90611823.nonecase&spm=1018.2226.3001.4450响应面法程序转载自

muX = [2;4]; 
sigmaX=[0.4;0.8]; 
rhoX =[1,0;0,1];
sLn = sqrt(log(1+(sigmaX./muX).^2));
mLn = log(muX)-sLn.^2/2;x = muX;
normX = eps;
n = length(muX); 
f = 2;d = diag(f*sigmaX);
gM = 1+(x(1)-1)*x(2);
while abs(norm(x)-normX)/normX > 1e-6
normX = norm(x);
a1 = x'; a2 = a1.*a1;
a3 = repmat(a1,n,1);
a3 = a3-d; 
a4 =a3.*a3;
a5 =a3+2*d;
a6 = a5.*a5;
a = [ones(2*n+1,1),[a1,a2;a3,a4;a5,a6]];
gv = 1+(a(:,2)-1).*a(:,3);
lambda = a\gv;
x = muX; 
normX1 = eps;
while abs(norm(x)-normX1)/normX1 > 1e-6
normX1 = norm(x);
g = lambda'*[1;x;x.*x];
gX = lambda(2:n+1)+2*lambda(n+2:2*n+1).*x;
cdfX = logncdf(x,mLn,sLn);
pdfX = lognpdf(x,mLn,sLn);
nc = norminv(cdfX);
sigmaX1 = normpdf(nc)./pdfX;muX1 = x-nc.*sigmaX1;
gs =gX.*sigmaX1;
alphaX = -rhoX*gs/sqrt(gs'*rhoX*gs);
bbeta = (g+gX'*(muX1-x))/sqrt(gs'*rhoX*gs);
x = muX1+bbeta*sigmaX1.*alphaX;
end
g = 1+(x(1)-1)*x(2);
x = muX+gM/(gM-g)*(x-muX)
end

(5条消息) 随机响应面法_步步为营！的博客-CSDN博客_随机响应面法https://blog.csdn.net/qq_41710383/article/details/115406601?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522165027187416780269851640%2522%252C%2522scm%2522%253A%252220140713.130102334.pc%255Fblog.%2522%257D&request_id=165027187416780269851640&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~blog~first_rank_ecpm_v1~rank_v31_ecpm-2-115406601.nonecase&utm_term=%E8%B6%85%E6%8B%89%E4%B8%81&spm=1018.2226.3001.4450

七、JC法

JC法程序转载于

muX = [2;4]; 
sigmaX = [0.4;0.8]; 
rhoX = [1,0;0,1];
sLn = sqrt(log(1+(sigmaX./muX).^2));
mLn = log(muX)-sLn.^2/2;
muX1 = muX;
sigmaX1 = sigmaX;
x = muX
normX = eps;
while abs(norm(x)-normX)/normX > 1e-6
normX = norm(x);
g = 1+prod(x)-x(2);
gX = [x(2);x(1)-1];
cdfX = logncdf(x,mLn,sLn); 
pdfX = lognpdf(x,mLn,sLn);
nc = norminv(cdfX);
sigmaX1 = normpdf(nc)./pdfX;
muX1 = x-nc.*sigmaX1;
gs = gX.*sigmaX1;
alphaX = -rhoX*gs/sqrt(gs'*rhoX*gs);
bbeta = (g+gX'*(muX1-x))/sqrt(gs'*rhoX*gs);
x = muX1+bbeta*sigmaX1.*alphaX;
end