基于坐标轮换法的罚函数求解有约束二维问题

Matlab实现程序如下：

function fa
x0=[3 3];                                           %初始点
r0=1.7;                                             %惩罚因子增长速度
c=1;                                                 %初始惩罚因子值
deta=0.001;                                     %收敛条件
%  坐标轮换法求极值点
x_temp=x0;
x_last=x0;
fla=1;

fanwei=0:0.01:10;                               %画图范围
[X,Y] = meshgrid(fanwei,fanwei);   
Z=(X-2).^2+(Y-1).^2; 
contour(X,Y,Z,20);
hold on;
plot(fanwei,fanwei.^2);
plot(fanwei,2-fanwei);
packer(:,1)=x0;
counter=2;

while           fla>0.0001
c=c*r0;
first_series=test(c,2,x_temp(2));
x_temp(1)=first_series;

sec_series=test(c,1,first_series);
x_temp(2)=sec_series;

fla=rms(x_temp-x_last);
x_last=x_temp;
packer(:,counter)=x_temp;
counter=counter+1;
end
plot(packer(1,:),packer(2,:),'ro');
plot(x_temp(1),x_temp(2),'r*');
x_optimzation=x_temp

end




function h=test(alpha,which,one_x)   %which为不变的量
i=which;
if i==1
    j=2;
else
    j=1;
end                        %j为寻优的变量
step=0.1;
n=0;
x=zeros(2,3);
x(i,:)=one_x;
% interval=[-10:0.01:20];
% plot(interval,f(interval));
% hold on
while f(x(:,1))<=f(x(:,2))|f(x(:,2))>=f(x(:,3))
   
    x(j,1)=x(j,2);
     x(j,2)=x(j,3);
    x(j,3)= x(j,3)+step*2^n;
    n=n+1;
   if x(j,3)>10^5

   step=-0.1;
   n=0;
   end
  
end
%  plot(x,f(x),'o');
 a=x(j,1);b=x(j,3);
x1=0.382*(b-a)+a;
x2=0.618*(b-a)+a;
xx1=zeros(2,1);
xx1(j)=x1;xx1(i)=one_x;
xx2=zeros(2,1);
xx2(j)=x2;xx2(i)=one_x;
aa=zeros(2,1);
aa(j)=a;aa(i)=one_x;
bb=zeros(2,1);
bb(j)=b;bb(i)=one_x;


while (b-a)/b>0.0001
if f(xx1)<f(xx2)
    b=xx2(j);
    xx2(j)=xx1(j);
    xx1(j)=0.382*(b-a)+a;
else
    a=xx1(j);
    xx1(j)=x2;
    xx2(j)=0.681*(b-a)+a;
    
end

end
xo=(a+b)/2;
% plot([x1,x2],f([x1,x2]),'r*');
h=xo;



function a=f(x)
a=(x(1)-2)^2+(x(2)-1)^2+alpha*max(x(1)^2-x(2),0)+alpha*max(x(1)+x(2)-2,0);   %%%%惩罚函数
end



end