问题描述
在无约束条件下，此函数本应该在（0，0）处有最小值0，但在上面的约束条件下，最后的结果应该是（1，1）处，最小值为2。
引入罚函数方法
通过罚函数，我们可以将粒子的活动范围限制在规定的约束区间内。
实施方法
文章来源：https://uudwc.com/A/DN3Er
代码（有详细注释）
// An highlighted block
clc;
clear;
close;
f=@(x,y) 20+x.^2+y.^2-10*cos(2*pi.*x)-10*cos(2*pi.*y);
x0=[-5.12:0.05:5.12];
y0=x0;
[X,Y]=meshgrid(x0,y0);
Z=f(X,Y);
figure(1)
mesh(X,Y,Z);
title(['初始粒子状态']);
hold on;
%设置PSO的相关参数
N=50; %种群大小
w_max=0.8; %速度惯性因子
c1=2; %自我学习因子
c2=2; %群体学习因子
dim=2; %可行解的维数
max_iter=100; %最大迭代次数
xlimit=[-5,5];%粒子位置限制范围
vlimit=[-0.5,0.5]; %粒子速度限制范围
%初始化种群参数
x=xlimit(1)+(xlimit(2)-xlimit(1))*rand(N,dim); %随机产生一个位置的二维的数组，范围在xlimit内
v=vlimit(1)+(vlimit(2)-vlimit(1))*rand(N,dim); %随机产生一个速度的二维的数组，范围在vlimit内
x_min=x; %每个个体的历史最佳位置，初始化为原位置
y_min=zeros(1,dim); %种群的历史最佳位置，先初始化为0
fx_min=ones(N,1)*inf; %每个个体的历史最佳适应度，初始化为无穷
fy_min=inf; %种群的历史最佳适应度，初始化为无穷
x_min_e=ones(N,1)*inf; %个体历史最优位置的罚函数
y_min_e=inf; %全局历史最优位置的罚函数
fitness_value_list=[]; %记录最佳历史适应度
scatter3(x(:,1),x(:,2),f(x(:,1),x(:,2)),'*r'); %初始点的散点图
%开始迭代
iter=1;
while iter<=max_iter %等于迭代次数时结束迭代
fx=f(x(:,1),x(:,2)); %计算x的适应度值
e1=max(0,-x(:,1)-x(:,2)+1.9); %约束条件1情况下的罚函数,约束条件:x+y>1.9
e2=max(0,x(:,1).*x(:,2)-1.2); %约束条件2情况下的罚函数,约束条件:x*y<1.2
e=e1+e2; %总的罚函数
%取个体最优适应度及其位置
for i=1:N
if x_min_e(i,:) <= 0.00001 && e(i,:)<=0.00001 %情况1，历史和当前，粒子都可行。方法：比较适应度大小，选取适应度小的粒子
if fx(i)<fx_min(i)
fx_min(i)=fx(i);
x_min(i,:)=x(i,:);
end
end
if x_min_e(i,:)<= 0.00001 && e(i,:) > 0.00001 %情况2，历史可行，当前不可行。方法：选取历史可行
fx_min(i)=fx_min(i);
x_min(i,:)=x_min(i,:);
x_min_e(i,:)=x_min_e(i,:);
end
if x_min_e(i,:) > 0.00001 && e(i,:) <= 0.00001 %情况3，历史不可行，当前可行。方法：选取当前可行
fx_min(i)=fx(i);
x_min(i,:)=x(i,:);
x_min_e(i,:)=e(i,:);
end
if x_min_e(i,:) > 0.00001 && e(i,:)>0.00001 %情况4，历史和当前都不可行。方法：选取罚函数小的粒子位置和适应度
if e(i,:)<x_min_e(i,:)
x_min_e(i,:)=e(i,:);
fx_min(i)=fx(i);
x_min(i,:)=x(i,:);
end
end
end
%取群体最优适应度及其位置
if y_min_e<=0.00001 && min(x_min_e)<=0.00001
index=find(x_min_e==0);
[a,n_min]=min(f(x(index,1),x(index,2)));
index_min=index(n_min); %首先找到符合约束条件的点
if f(x(index_min,1),x(index_min,2))<fy_min %如果满足约束条件点的适应度小于全局最优
fy_min=f(x(index_min,1),x(index_min,2)); %更新全局适应度
y_min=x_min(index_min,:); %更新最优位置
end
end
if y_min_e>0.0001 && min(x_min_e)<= 0.00001
index=find(x_min_e==0);
[a,n_min]=min(f(x(index,1),x(index,2)));
index_min=index(n_min); %首先找到符合约束条件的点
fy_min=f(x(index_min,1),x(index_min,2)); %更新全局最优适应度
y_min=x_min(index_min,:); %更新最优位置
y_min_e=min(x_min_e); %更新最小罚函数
end
if y_min_e<=0.0001 && min(x_min_e)>0.00001 %说明全局位置符合，当前位置不符合，选取全局位置
end
if y_min_e>0.00001 && min(x_min_e)>0.00001
if min(x_min_e)<y_min_e
[a,n_min]=min(x_min_e);
fy_min=f(x(n_min,1),x(n_min,2));
y_min=x_min(n_min,:);
y_min_e=min(x_min_e);
end
end
%更新速度
w=w_max-w_max*(iter/max_iter); %使w的值随迭代次数增大不断递减，提高其局部搜索能力，加快收敛
v=w*v+c1*rand*(x_min-x)+c2*rand*(repmat(y_min,N,1)-x);
v(v>vlimit(2))=vlimit(2);
v(v<vlimit(1))=vlimit(1);
%更新位置
x=x+v;
x(x>xlimit(2))=xlimit(2);
x(x<xlimit(1))=xlimit(1);
fitness_value_list(iter)=fy_min;
figure(2)
subplot(1,2,1)
mesh(X,Y,Z);
hold on;
scatter3(x(:,1),x(:,2),f(x(:,1),x(:,2)),'*r');
pause(0.1);
hold off;
title(['迭代次数为：',num2str(iter), '粒子状态']);
subplot(1,2,2)
plot(fitness_value_list);
title(['适应度值']);
iter=iter+1;
end
figure(3)
mesh(X,Y,Z);
hold on;
scatter3(x(:,1),x(:,2),f(x(:,1),x(:,2)),'*r');
title(['迭代次数为：',num2str(iter-1), '时的粒子状态']);
disp(['最优值为',num2str(fy_min)])
disp(['位置为',num2str(y_min)])
运行结果
文章来源地址https://uudwc.com/A/DN3Er
