前一段师姐帮学姐用 粒子群算法计算带线性等式约束的最优问题。matlab自带的粒子群算法并没有这个功能,所以对原有的粒子群增加了两个步骤来进行线性等式约束。
没理解基础粒子群算法的童鞋请先去别的文章处弄懂,我是在这里看的 推荐一下 https://blog.csdn.net/nightmare_dimple/article/details/74331679
增加的两个步骤。
1. 引入惩罚因子,在离线性等式较远的地方让他适应度(函数值)无穷大(我找的是最小值)。
2.速度分情况给值,如果不在那条线附近(线性等式约束),那么他的速度只是群体学习速度。
我觉得这种思想(惩罚因子和速度选择)也是可以解决非线性约束的,有需要的童鞋可以验证一下。
clear
% 初始化种群
f_1 = @(x)0.01011.*x.^2-1.49559.*x+230.38805;
f_2 = @(x)0.0103.*x.^2-1.52432.*x+233.45944;
f_min = @(x,y)2*f_1(x/16800)*x+2*f_2(y/9600)*y;
tic;
N = 5000; % 初始种群个数 粒子数为1000个
d = 2; % 空间维数
threshold =45800; % x1+x2+x3+x4≥45800 的阈值(线性不等式约束)
ger = 180; % 最大迭代次数 (我这个点比较多,迭代次数不用那么多)
limit = [14000, 15220;7700,8700]; % X 的上下限
vlimit = [- 7.5, 7.5;-7.5, 7.5]; % 种群移动速度的上下限
c_1 = 0.8; % 惯性权重
c_2 = 0.5; % 自我学习因子
c_3 = 0.5; % 群体学习因子
for i = 1:d
x(:,i) = limit(i, 1) + (limit(i, 2) - limit(i, 1)) * rand(N, 1);%初始种群的位置
end
v = rand(N, d); % 初始种群的速度
xm = x; % 每个个体的历史最佳位置
ym = zeros(1, d); % 种群的历史最佳位置
fxm = zeros(N, 1)+inf; % 每个个体的历史最佳适应度
fym = inf; % 种群历史最佳适应度
iter = 1; % 初始迭代次数为一
times = 1; % 进度条计数用的
bar = waitbar(0,'读取数据中...'); % 查看进度用的进度面板
ezplot('(96/168)*x',[10000,16000]);
hold on
while iter <= ger
str=['计算中...',num2str(100*iter/ger),'%']; % 百分比形式显示处理进程
waitbar(iter/ger,bar,str) % 更新进度条bar,配合bar使用
for i = 1:N
fx(i) = f_min(x(i,1),x(i,2)) ; % 个体当前适应度
if ((2*x(i,1)+2*x(i,2))<threshold) % 若不满足x1+x2+x3+x4≥45800要求则惩罚因子无穷大
fx(i) = inf;
end
if((abs(9600*x(i,1)-16800*x(i,2)))>=500) % 若不满足|9600*x1-16800*x2|>=500 则惩罚因子无穷大
fx(i) = inf;
end
end
for i = 1:N
if (fxm(i) > fx(i))
fxm(i) = fx(i); % 更新个体历史最佳适应度(即:看谁的f_min更小,并把每个粒子历史上最小的f_min存起来)
xm(i,:) = x(i,:); % 更新个体历史最佳位置
end
end
if fym > min(fxm)
[fym, nmin] = min(fxm); % 更新群体历史最佳适应度 找最小值
ym = xm(nmin, :); % 更新群体历史最佳位置
end
for i = 1:N
if(((2*x(i,1)+2*x(i,2))<threshold))
v(i,:) = c_3 * rand *(ym - x(i,:)); %如果不满足x1+x2+x3+x4≥45800,则速度方向取决于群体学习
elseif((9600*x(i,1)~=16800*x(i,2))) %如果不满足|9600*x1-16800*x2|>=500,则速度方向取决于群体学习
v(i,:) = c_3 * rand *(ym - x(i,:));
else
v(i,:) = v(i,:) * c_1 + c_2 * rand *(xm(i,:) - x(i,:)) + c_3 * rand *(ym - x(i,:));% 速度更新
end
end
% 边界速度处理
for i=1:d
for j=1:N
if v(j,i)>vlimit(i,2)
v(j,i)=vlimit(i,2);
end
if v(j,i) < vlimit(i,1)
v(j,i)=vlimit(i,1);
end
end
end
if iter % 100 == 0
ezplot('(96/168)*x',[10000,16000]);
hold on
end
x = x + v;% 位置更新
for i=1:d % 边界位置处理
for j=1:N
if x(j,i)>limit(i,2)
x(j,i)=limit(i,2);
end
if x(j,i) < limit(i,1)
x(j,i)=limit(i,1);
end
end
end
if times >= 10 %这是画动态图用的
cla;
scatter(x(:,1),x(:,2));title('状态位置变化');
xlabel(" X1 ");
ylabel(" X3 ");
pause(0.5); % 暂停0.5秒
times=0;
end
iter = iter+1;
times=times+1;
end
toc;