%初始化初始点X0、精度theta、惩罚因子M1和放大系数C
theta=0.00001;
r=3;
C=0.99;
k=1;
syms x1 x2
% x0=input('请输入初始点');
% alpha=input('请输入障碍函数');
% ff=input('请输入目标函数');
x0=[3,4];
alpha=1/(x1-1)+1/x2;
ff=1/3*((x1+1)^3)+x2;
% f=r*alpha+ff;%增广函数
% g=jacobian(f);%梯度
% Hesse=jacobian(jacobian(f));%Hesse矩阵
% H1=inv(Hesse);
% g0=subs(g,{x1, x2},X0);%将X0代入到梯度中
% f0=subs(f,{x1, x2},X0);%将XO代入到增广函数中
while 1
% time=0;
% while 1
% time=time+1;
% %牛顿法求出Xk
% P=-inv(subs(Hesse,{x1,x2},X0))*(subs(g,{x1, x2},X0))';
% a=double(P);
% % G=subs(H1,{x1, x2},X0);%将XO代入到Hesse矩阵的逆矩阵中
% % P=-G*g0';%迭代方向
% X=X0'+P;%迭代出新的点
% % g1= subs(g, {x1, x2},X');%新的点代入梯度中
% % a=double(norm(g1));%计算当前点的梯度的模
% b=double(norm(subs(g,{x1, x2},X')))
% if b<0.01 %符合精度条件,跳出循环
% break;
% else
% X0=X';%不符合精度条件,更新点和梯度
% % g0=g1;
% end
% end
eval(['F=@(x) 1/3*(x(1)+1)^3+x(2)+' num2str(r) '*(1/(x(1)-1))+' num2str(r) '*(1/x(2));']);
X=fminunc(F,x0);
alpha1=subs(alpha, {x1, x2},X);%将用牛顿法求出的Xk代入罚函数中
a=r*double(alpha1);
if a<theta %如果符合精度条件,则将该点代入目标函数中求值
% if k>=300
F=subs(ff, {x1, x2},X);
X=double(X);
F=double(F);
fprintf('迭代了%d步\n',k);
fprintf('最优解为X=[%f,%f]\n',X(1),X(2));
fprintf('最优值为%f',F);
break;
end
r=r*C; %不符合精度条件,则增大惩罚因子
k=k+1;
% f=r*alpha+ff;%增广函数
% g=jacobian(f);%梯度
% Hesse=jacobian(jacobian(f));%Hesse矩阵
% H1=inv(Hesse);
x0=X;
% g0=subs(g,{x1, x2},X0);%将X0代入到梯度中
% f0=subs(f,{x1, x2},X0);%将XO代入到增广函数中
end