[ncgcake]

目标函数：122.5221*sqrt(0.9025+x1^2)*x2 约束：(1.2732e+007*sqrt(0.9025+x1^2)/(x1*x2)-4.2e+008)<=0,(25908*(0.0025^2+x2^2)*x1*x2-1.2732*(0.9025+x1^2)^1.5)+1/(x2-0.005)+1/(x1)>=0
clc
m=zeros(1,100);a=zeros(1,100);b=zeros(1,100);f0=zeros(1,100);%a b为最优点坐标，f0为最优点函数值，f1 f2最优点梯度。
syms x1 x2 e; %e为罚因子。
m(1)=1;c=0.2;a(1)=0.75;b(1)=0.08; %c为递增系数,赋初值。
f=122.5221*sqrt(0.9025+x1^2)*x2+e*(1/(42*(x1*x2)-1.2732*sqrt(0.9025+x1^2))+1/(25908*(0.0025^2+x2^2)*x1*x2-1.2732*(0.9025+x1^2)^1.5)+1/(x2-0.005)+1/(x1));
%x1为高度h，x2为管直径D。
fx1=diff(f,'x1');fx2=diff(f,'x2');fx1x1=diff(fx1,'x1');fx1x2=diff(fx1,'x2');fx2x1=diff(fx2,'x1');fx2x2=diff(fx2,'x2');%求偏导、海森元素。
for k=1:5 %外点法e迭代循环.
x1=a(k);x2=b(k);e=m(k);
for n=1:1000 %牛顿法求最优值。
f1=subs(fx1); %求解梯度值和海森矩阵
f2=subs(fx2);
f11=subs(fx1x1);
f12=subs(fx1x2);
f21=subs(fx2x1);
f22=subs(fx2x2);
if(double(sqrt(f1^2+f2^2))<0.01) %最优值收敛条件
a(k+1)=double(x1);b(k+1)=double(x2);f0(k+1)=double(subs(f));
break;
else %如未最优值收敛条件，对x1,x2进行迭代,直到寻找到合适的迭代点
X=[x1 x2]-(([f1 f2]*[f1 f2]')/([f1 f2]*[f11 f12;f21 f22]*[f1 f2]'))*[f1 f2]';
x1=X(1,1);x2=X(2,1);
end
end
if(double(sqrt((a(k+1)-a(k))^2+(b(k+1)-b(k))^2))<=0.01)&&(double(abs((f0(k+1)-f0(k))/f0(k)))<=0.01) %罚因子迭代收敛条件
h=a(k+1) %输出最优点h,D值，罚因子迭代次数，最优值
D=b(k+1)
times=k
zhiliang=f0(k+1)
break;
else
m(k+1)=c*m(k);

h=a(k+1) %输出最优点h,D值，罚因子迭代次数，最优值
D=b(k+1)
times=k
zhiliang=f0(k+1)
end
end
有什么错误，请详细指出来，谢谢了