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