大家好，又见面了，我是你们的朋友全栈君。

Ｍatlab 非线性有约束规划的粒子群算法

粒子群算法的基本认识

简单介绍：通过群体中个体之间的协作和信息共享来寻找最优解。

适用于连续函数极值问题，对于非线性，多峰问题均有较强的全局搜索能力。

主要掌握两点

1.粒子的速度和位置

速度代表移动的快慢，位置代表移动的方向。
位置对应每个自变量，速度一般设置为变量范围的10%~20%。

2.粒子的更新规则

具体实例

ｍatlab代码

clear;close;clc
%% 约束条件和目标函数构建
fun = @(x) x(1)^2 + x(2)^2 + x(3)^2 + 8;
bind1 = @(x) x(1)^2 - x(2) + x(3)^2 >= 0;
bind2 = @(x) x(1) + x(2)^2 + x(3)^2 <= 20;
% 不太适合等式约束
ekc = 1e-10;
bind3 = @(x) abs(-x(1) - x(2)^2 + 2) <= ekc;
bind4 = @(x) abs(x(2) + 2*x(3)^2 - 3) <= ekc;
%% 初始化
popsize = 500; % 粒子个数
dim = 3; % 维度
max_iter = 100; % 最大迭代次数
xlimit_max = [2 3 20]'; % 由等式约束推出位置边界
xlimit_min = zeros(dim,1); 
vlimit_max = 1*ones(dim,1);
vlimit_min = -1*ones(dim,1);
w = 0.6; % 惯性权重
c1 = 0.5;c2 = 1.5;
pr = 0.4; % 变异率
pop_x = zeros(dim,popsize);  % 当前粒子位置
pop_v = zeros(dim,popsize); % 当前粒子速度
fitness_pop = zeros(1,popsize); % 粒子群当前位置适应度函数
fitness_lbest = zeros(1,popsize); % 个体粒子的历史最优极值
rand('state',sum(clock));
for j = 1:popsize 
% 位置初始化
pop_x(1,j) = xlimit_min(1) + rand*(xlimit_max(1) - xlimit_min(1));
pop_x(2,j) = sqrt(2-pop_x(1,j));
pop_x(3,j) = sqrt((3 - pop_x(2,j))/2);
% 速度初始化
for  i = 1:dim
pop_v(i,j) = vlimit_min(i) + rand*(vlimit_max(i) - vlimit_min(i));
end
end
%% 初始化个体极值
lbest = pop_x; % 个体历史最佳极值记录
for j =1: popsize 
if bind1(pop_x(:,j))
if bind2(pop_x(:,j))
fitness_lbest(j) = fun(pop_x(:,j));
else fitness_lbest(j) = 500;
end
else fitness_lbest(j) = 500;
end
end
%% 初始化全局极值
popbest = pop_x(:,1);
fitness_popbest = fitness_lbest(1);
for j = 2:popsize 
if fitness_lbest(j) < fitness_popbest
fitness_popbest = fitness_lbest(j);
popbest = pop_x(:,j);
end
end
tic
%% 粒子群迭代
iter = 1; % 当前迭代次数
record = zeros(max_iter,1); % 记录每次迭代的全局极小值
format long;
while iter <= max_iter
for j = 1:popsize 
% 更新速度 边界处理
pop_v(:,j) = w*pop_v(:,j) + c1*rand*(lbest(:,j) - pop_x(:,j)) +...
c2*rand*(popbest - pop_x(:,j));
for i = 1:dim 
if pop_v(i,j) > vlimit_max(i)
pop_v(i,j) = vlimit_max(i);
elseif pop_v(i,j) < vlimit_min(i) 
pop_v(i,j) = vlimit_min(i);
end
end
% 更新位置 边界处理 修正位置 （等式约束）
pop_x(:,j) = pop_x(:,j) + pop_v(:,j);
for i = 1:dim 
if pop_x(i,j) > xlimit_max(i)
pop_x(i,j)  = xlimit_max(i);
elseif pop_x(i,j) < xlimit_min(i)
pop_x(i,j) = xlimit_min(i);
end
end
% 进行自适应变异
if rand < pr 
i = ceil(dim*rand);
pop_x(i,j) = xlimit_min(i) + rand*(xlimit_max(i) - xlimit_min(i));
end
% 约束条件限制 类似罚函数法
if bind1(pop_x(:,j))
if bind2(pop_x(:,j))
if bind3(pop_x(:,j))
if bind4(pop_x(:,j))
fitness_pop(j) = fun(pop_x(:,j));
else fitness_pop(j) = 500;
end
else fitness_pop(j) = 500;
end
else fitness_pop(j) = 500;
end
else fitness_pop(j) = 500;
end
% 当前适应度与个体历史最佳适应度作比较
if fitness_pop(j) < fitness_lbest(j)
lbest(:,j) = pop_x(:,j);
fitness_lbest(j) = fitness_pop(j);
end
% 个体历史最佳适应度与种群历史最佳适应度作比较
if fitness_popbest > fitness_lbest(j)
fitness_popbest = fitness_lbest(j);
popbest = lbest(:,j);
end
end
record(iter) = fitness_popbest;
iter = iter + 1;
end
toc
%% 输出解
minx = popbest
miny = fitness_popbest
plot(record,'r-');
title('粒子群算法迭代过程');
xlabel('迭代次数');
ylabel('当前迭代最佳函数值');

结果：

从以上结果可以看出，粒子群算法几乎一开始就保持收敛，说明对于小规模的粒子群，它收敛的迅速

Lingo求解全局最小值比较：