当量正态化的matlab程序

案例：已知非线性状态方程 567fr-0.5H2 =0, f 服从正态分布， =0.6 ，变异系数 =0.131 ， r 服从正态分布， =2.18 ， =0.03 ； H 服从对数正态分布， =32.8, =0.03. 用 JC 法计算可靠指标 及设计验算点坐标（ f* ， r* ， H* ）。 解：功能函数梯度为 g （ f ， r ， H ） = （ 567r ， 567f ， -H ） T 。关于如何确定对数正态分布的参数见附录。

Matlab 代码如下：

clear;clc;

muX=[0.6;2.18;32.8]; cvX=[13.1;3;3]/100; sigmaX=cvX.*muX; （初始均值，变异系数，标准差）

sLn=sqrt(log(1+(sigmaX(3)/muX(3))^2));mLn=log(muX(3))-sLn^2/2; （对数正态分布的初始正态化）

muX1=muX;sigmaX1=sigmaX;

x=muX; normX=eps;

while abs(norm(x)-normX)/normX>1e-6

normX=norm(x);

g=567*x(1)*x(2)-x(3)^2/2;

gX=[567*x(2);567*x(1);-x(3)];

cdfX=logncdf(x(3),mLn,sLn); （求当量正态分布函数）

pdfX=lognpdf(x(3),mLn,sLn); （求当量正态分布的概率密度函数） nc=norminv(cdfX); （求当量正态分布函数的反函数）

sigmaX1(3)=normpdf(nc)/pdfX; （求当量正态分布函数的标准差） muX1(3)=x(3)-nc*sigmaX1(3);

gs=gX.*sigmaX1;

alphaX=-gs/norm(gs); （cosx）

bbeta=(g+gX'*(muX1-x))/norm(gs) (可靠度计算)

x=muX1+bbeta*sigmaX1.*alphaX （新的样本点）

end

结果：设计验算点坐标（ f* ， r* ， H* ） = （ 0.4561 ， 2.1590 ， 33.4178 ），beta =1.9645 。