5—2[解] 用应力表示的极限状态方程为()223Z ,,,04u u ql g f q b d f bd ==-=利用书本上导出的梯度优化法计算公式,利用matlab 编写了电算程序(见附录1),计算结失效概率为:1()0.927f P β=-Φ=5-7[解] 用应力表示的极限状态方程为223332()16(,,)()044Y Y M f M d Z g M T f d ππ==+-= 利用书本上导出的梯度优化法计算公式,利用matlab 编写了电算程序(见附录2),计算结果见下表。
结果分析:因为题目中存在两个非正态变量,需要将其转化为正态变量。
而利用梯度最优化方法,是利用标准正态变量进行迭代且迭代过程中用到了转为了正态分布时的均值和方差。
所以,此时就需要利用每次迭代出的标准正态变量,求解非正态变量及其对应的正态分布时的均值与方差,同样需要迭代。
怀疑是在求解非正态变量及其对应的正态分布时的均值与方差时,迭代出现了问题。
所以用网格搜索法(附录3)求解,结果仍不理想。
β=,以-0.01为步为进一步找到问题所在,尝试利用改进一次二阶矩法。
取初值 2.0长,发现直到β为负值时,结果仍未收敛。
由于时间有限,未能找到错误原因,希望老师能给予点拨。
%例题5-2(梯度优化法)syms Yf Yq Yb YdG=1377.88*Yf+3936.8-100390.2*(8.928*Yq+111.6)/((0.6225*Yb+12.45)*(1.5 39*Yd+25.65)^2);gradient=jacobian(G,[Yf,Yq,Yb,Yd]);Y=zeros(10,4);for i=1:10beta(i)=sqrt(Y(i,:)*Y(i,:)')Pf=1-normpdf(beta(i))Yf=Y(i,1);Yq=Y(i,2);Yb=Y(i,3);Yd=Y(i,4);tiduzhi=eval(gradient);shuzhi(i)=eval(G)a(i,:)=-tiduzhi./norm(tiduzhi) %normΪÇóÏòÁ¿µÄÄ£Y(i+1,:)=(Y(i,:)*a(i,:)'+shuzhi(i)/norm(tiduzhi))*a(i,:);endX(:,1)=Y(:,1)*1377.88+3936.8;X(:,2)=Y(:,2)*8.928+111.6;X(:,3)=Y(:,3)*0.6225+12.45;X(:,4)=Y(:,4)*1.539+25.65;图一程序运行图%例题5-7(梯度优化法——迭代)X=[26 17];Sigma=[4.68 2.38];M=[26 17];i=1;Y=zeros(5,3)for p=1:2x(p)=X(i,p);sigma=Sigma(p);m=M(p);a(p)=1.2825/sigma , k(p)=m-0.5772/a(p)F=exp(-exp(-a(p)*(x(p)-k(p))));f=a(p)*exp(-a(p)*(x(p)-k(p))-exp(-a(p) *(x(p)-k(p))));sigma_tran(i,p)=normcdf(norminv(F))/f;m_tran(i,p)=x(p)-sigma_tran(i,p)*norminv(F);Y(i,p)=(x(p)-m_tran(i,p))/sigma_tran(i,p);endfor i=1:5syms YM YT YfG=(159.155*(YM*sigma_tran(i,1)+m_tran(i,1)))^2/4+(79.577*(YT*sigm a_tran(i,2)+m_tran(i,2)))^2-(Yf*840+7*10^3)^2/4;gradient=jacobian(G,[YM,YT,Yf]);beta(i)=sqrt(Y(i,:)*Y(i,:)');Pf=1-normpdf(beta(i));YM=Y(i,1);YT=Y(i,2);Yf=Y(i,3);tiduzhi(i,:)=eval(gradient);shuzhi(i)=eval(G);b(i,:)=-tiduzhi(i,:)./norm(tiduzhi(i,:));Y(i+1,:)=(Y(i,:)*b(i,:)'+shuzhi(i)/norm(tiduzhi(i,:)))*b(i,:);Xattemp(1,1)=Y(i+1,1)*sigma_tran(i,1)+m_tran(i,1);Xattemp(1,2)=Y(i+1, 2)*sigma_tran(i,2)+m_tran(i,2);for p=1:2for M=1:10syms x mm sigmaaQ=Y(i+1,p)-(x-mm)/sigmaa;Gradient=jacobian(Q,[x]);x=Xattemp(M,p)F=exp(-exp(-a(p)*(x-k(p))));f=a(p)*exp(-a(p)*(x-k(p))-exp(-a(p)*(x-k( p))));sigma_tra(M,p)=normcdf(norminv(F))/f;m_tra(M,p)=x-sigma_tra(M,p)*norminv(F);mm=m_tra(M,p);sigmaa=sigma_tra(M,p);Tiduzhi=eval(Gradient);Shuzhi(i)=eval(Q);c(M)=-Tiduzhi/norm(Tiduzhi);Xattemp(M+1,p)=(Xattemp(M,p)*c(M)'+Shuzhi(i)/norm(Tiduzhi))*c(M);if Shuzhi(i)<=1breakendendX(i+1,p)=x;sigma_tran(i+1,p)=sigmaa;m_tran(i+1,p)=mm; endendX(:,3)=Y(:,3)*840+7000;图二程序运行图附录三%例题5-7(梯度优化法——网格搜索)X=[26 17];Sigma=[4.68 2.38];M=[26 17]; %µÚÒ»¸öÊÇÍ侨£¬µÚ¶þ¸öÊÇŤ¾Øi=1;Y=zeros(10,3)for p=1:2x(p)=X(i,p);sigma=Sigma(p);m=M(p);a(p)=1.2825/sigma , k(p)=m-0.5772/a(p)F=exp(-exp(-a(p)*(x(p)-k(p))));f=a(p)*exp(-a(p)*(x(p)-k(p))-exp(-a(p) *(x(p)-k(p))));sigma_tran(i,p)=normcdf(norminv(F))/f;m_tran(i,p)=x(p)-sigma_tran(i,p)*norminv(F);Y(i,p)=(x(p)-m_tran(i,p))/sigma_tran(i,p);clear sigma x mendfor i=1:100syms YM YT YfG=(159.155*(YM*sigma_tran(i,1)+m_tran(i,1)))^2/4+(79.577*(YT*sigma_tr an(i,2)+m_tran(i,2)))^2-(Yf*840+7*10^3)^2/4;gradient=jacobian(G,[YM,YT,Yf]);beta(i)=sqrt(Y(i,:)*Y(i,:)');Pf=1-normpdf(beta(i));YM=Y(i,1);YT=Y(i,2);Yf=Y(i,3);tiduzhi=eval(gradient);shuzhi(i)=eval(G);b(i,:)=-tiduzhi./norm(tiduzhi); %normΪÇóÏòÁ¿µÄÄ£Y(i+1,:)=(Y(i,:)*b(i,:)'+shuzhi(i)/norm(tiduzhi))*b(i,:);Xattemp(i,1)=Y(i+1,1)*sigma_tran(i,1)+m_tran(i,1);Xattemp(i,2)=Y(i+1, 2)*sigma_tran(i,2)+m_tran(i,2);for p=1:2for x=Xattemp(i,p)-20:0.01:Xattemp(i,p)+20F=exp(-exp(-a(p)*(x-k(p))));f=a(p)*exp(-a(p)*(x-k(p))-exp(-a(p)*(x-k( p))));sigma_tra(p)=normcdf(norminv(F))/f;m_tra(p)=x-sigma_tra(p)*norminv(F);mm=m_tra(p);sigmaa=sigma_tra(p);Q(i)=Y(i+1,p)-(x-mm)/sigmaa;if Q(i)<=0.01breakendendX(i+1,p)=x;sigma_tran(i+1,p)=sigmaa;m_tran(i+1,p)=mm; endend图三程序运行图附录四%例题5-7(改进一次二阶矩法)Sigma0=[4.68 2.38 840];M0=[26 17 7*10^3];%第一次迭代beta(1)=2.5;X_design(1,1)=26;X_design(1,2)=17;X_design(1,3)=7*10^3;y(1)=1;y(2)=1;y(3)=1;while y(1)>0.01&y(2)>0.01&y(3)>0.01for p=1:2x(p)=X_design(1,p);sigma=Sigma0(p);m=M0(p);a(p)=1.2825/sigma , k(p)=m-0.5772/a(p)F=exp(-exp(-a(p)*(x(p)-k(p))));f=a(p)*exp(-a(p)*(x(p)-k(p))-exp(-a(p) *(x(p)-k(p))));sigma_tran(1,p)=normcdf(norminv(F))/f;m_tran(1,p)=x(p)-sigma_tran(1,p)*norminv(F);clear sigma x m f Fendm_tran(1,3)=M0(3);sigma_tran(1,3)=Sigma0(3);syms M T fuG=(159.155*M)^2/4+(79.577*T)^2-(fu)^2/4;gradient=jacobian(G,[M,T,fu]);M=X_design(1,1);T=X_design(1,2);fu=X_design(1,3);tiduzhi=eval(gradient);alpha(1,:)=sigma_tran(1,:).*tiduzhi/norm(sigma_tran(1,:).*tiduzhi);X_circle=m_tran(1,:)-beta(1)*alpha(1,:).*sigma_tran(1,:);y(1)=abs((X_circle(1)-X_design(1,1))/X_design(1,1));y(2)=abs((X_circle(2)-X_design(1,2))/X_design(1,2));y(3)=abs((X_circle(3)-X_design(1,3))/X_design(1,3));X_design(1,:)=X_circle;endM=X_design(1,1);T=X_design(1,2);fu=X_design(1,3);shuzhi(1)=eval(G);clear y X_circle M T fu G tiduzhi gradient%第二次迭代beta(2)=2;X_design(2,:)=X_design(1,:);y(1)=1;y(2)=1;y(3)=1;while y(1)>0.01&y(2)>0.01&y(3)>0.01for p=1:2x(p)=X_design(2,p);sigma=Sigma0(p);m=M0(p);a(p)=1.2825/sigma , k(p)=m-0.5772/a(p)F=exp(-exp(-a(p)*(x(p)-k(p))));f=a(p)*exp(-a(p)*(x(p)-k(p))-exp(-a(p) *(x(p)-k(p))));sigma_tran(2,p)=normcdf(norminv(F))/f;m_tran(2,p)=x(p)-sigma_tran(2,p)*norminv(F);clear sigma x m f Fendm_tran(2,3)=M0(3);sigma_tran(2,3)=Sigma0(3);syms M T fuG=(159.155*M)^2/4+(79.577*T)^2-(fu)^2/4;gradient=jacobian(G,[M,T,fu]);M=X_design(2,1);T=X_design(2,2);fu=X_design(2,3);tiduzhi=eval(gradient);alpha(2,:)=sigma_tran(2,:).*tiduzhi/norm(sigma_tran(2,:).*tiduzhi);X_circle=m_tran(2,:)-beta(2)*alpha(2,:).*sigma_tran(2,:);y(1)=abs((X_circle(1)-X_design(2,1))/X_design(2,1));y(2)=abs((X_circle(2)-X_design(2,2))/X_design(2,2));y(3)=abs((X_circle(3)-X_design(2,3))/X_design(2,3));X_design(2,:)=X_circle;endM=X_design(2,1);T=X_design(2,2);fu=X_design(2,3);shuzhi(2)=eval(G);clear y X_circle M T fu G tiduzhi gradient%之后迭代i=2;while abs(shuzhi(i))>0.01i=i+1;beta(i)=beta(i-1)-0.01%shuzhi(i-1)*(beta(i-1)-beta(i-2))/(shuzhi(i-1) -shuzhi(i-2));X_design(i,:)=X_design(i-1,:);y(1)=1;y(2)=1;y(3)=1;while y(1)>0.01&y(2)>0.01&y(3)>0.01for p=1:2x(p)=X_design(i,p);sigma=Sigma0(p);m=M0(p);a(p)=1.2825/sigma , k(p)=m-0.5772/a(p)F=exp(-exp(-a(p)*(x(p)-k(p))));f=a(p)*exp(-a(p)*(x(p)-k(p))-exp(-a(p) *(x(p)-k(p))));sigma_tran(i,p)=normcdf(norminv(F))/f;m_tran(i,p)=x(p)-sigma_tran(i,p)*norminv(F);clear sigma x m f Fendm_tran(i,3)=M0(3);sigma_tran(i,3)=Sigma0(3);syms M T fuG=(159.155*M)^2/4+(79.577*T)^2-(fu)^2/4;gradient=jacobian(G,[M,T,fu]);M=X_design(i,1);T=X_design(i,2);fu=X_design(i,3);tiduzhi=eval(gradient);alpha(i,:)=sigma_tran(i,:).*tiduzhi/norm(sigma_tran(i,:).*tiduzhi); X_circle=m_tran(i,:)-beta(i)*alpha(i,:).*sigma_tran(i,:);y(1)=abs((X_circle(1)-X_design(i,1))/X_design(i,1));y(2)=abs((X_circle(2)-X_design(i,2))/X_design(i,2));y(3)=abs((X_circle(3)-X_design(i,3))/X_design(i,3));X_design(i,:)=X_circle;endM=X_design(i,1);T=X_design(i,2);fu=X_design(i,3);shuzhi(i)=eval(G);clear y X_circle M T fu G gradientend图四部分beta值图五部分设计验算点。