在上一篇博客 粒子群优化算法(1)中介绍了基本的粒子群算法,基本粒子群算法是基于连续空间(区间)进行搜索,然而在一些实际的工程应用中,我们的待求解的变量可能并不是历需的,而实一种离散型的变量。这就需要对基本的粒子群算法做出一些相应的改进。
在离散粒子群算法中,将离散问题空间映射到连续粒子运动空间,并做适当的修改。任然保留经典粒子群算法的速度-位置更新策略。粒子在状态空间的取值只限于0,1两个值,而速度的每一个位代表的是粒子位置所对应的位取值为0/1的可能性。因此在离散粒子群算法中,粒子速度的更新公式依然保持不变,但是个体最优位置和全局最优位置每一位的取值只能为0/1。
离散粒子群算法速度和位置的更新方式如下所示:
粒子速度的更新方式不变,但是其位置的更新方式如下所示:
(利用sigmoid函数将例子的速度映射到0-1之间)
则粒子的速度可以表示为:
例如,对于经典的0-1背包问题:如果采用粒子群优化算法求解,只能是采用离散的粒子群算法:
问题描述:
有一个容量为V的背包,有N件物品,第i件物品的体积是c(i),价值是w(i).那么将哪些物品放入背包中,可以使得价值最大。假设具体的数据如下:
假设N=10, V = 300, 每件物品的体积为:[95, 75, 23, 73, 50, 22, 6, 57, 89, 98], 物品的价值为[89, 59,19, 43, 100, 72, 44, 16, 7, 64]
具体的计算流程如下:
1. 对速度进行初始化,对粒子的位置进行初始化,这里粒子的位置采用的二进制的表示方法。1表示选择该物体,0表示不选择该物体。其他的参数的初始化的过程与前面的方法一样。
源代码:
clc;
clear all;
NP = 100; % 种群个数
G = 200; % 迭代次数
D = 10; % 决策变量的维度
c1 = 1.5; % 学习因子
c2 = 1.5;
w_max = 0.8; % 惯性权重
w_min = 0.4;
v_max = 5; % 粒子的速度限制
v_min = -5;
V = 300; % 背包容量
capacity = [95 75 23 73 50 22 6 57 89 98]; % 物品的体积
weight = [89 59 19 43 100 72 44 16 7 64]; % 物品的价值
penality = 2;
% 初始化种群个体
x = (rand(NP,D)>0.5); % 产生均匀分布的二进制串 randn产生的是符合正态分布的随机数
v = v_min + rand(NP,D)*(v_max - v_min); % 速度进行初始化
% 初始化个体最优
individual_best = x; % 每个个体的历史最优
pbest = zeros(NP, 1); % 个体最优位置对应的适应度值
for k=1:NP
pbest(k, 1) = func(individual_best(k, :), capacity, weight, V, penality);
end
% 初始化全局最优
global_best = zeros(1, D);
global_best_fit = eps;
for k=1:NP
temp = func(individual_best(k, :), capacity, weight, V, penality);
if temp > global_best_fit
global_best = individual_best(k, :);
global_best_fit = temp;
end
end
% 进行迭代
for gen = 1:G
% 计算动态惯性权值
w = w_max - (w_max-w_min) * gen / G;
for k=1:NP
% 更新速度
v(k, :) = w * v(k, :) + c1 * rand() * (individual_best(k, :) - x(k, :)) + c2 * rand() * (global_best - x(k, :));
% 边界条件处理 % 边界吸收
for t=1:D
if v(k, D) > v_max
v(k, D) = v_max;
end
if v(k, D) < v_min
v(k, D) = v_min;
end
end
% 使用sigmoid函数对速度进行映射
vs(k, :) = 1./(1+exp(-v(k, :)));
% 更新粒子的位置
for t=1:D
if vs(k, t)>rand()
x(k, t) = 1;
else
x(k, t) = 0;
end
end
end
% 计算个体历史最优与全局最优
% 个体历史最优
for k=1:NP
old_fitness = func(individual_best(k, :), capacity, weight, V, penality);
new_fitness = func(x(k, :), capacity, weight, V, penality);
if new_fitness > old_fitness
individual_best(k, :) = x(k, :);
pbest(k, 1) = new_fitness;
end
end
% 全局最优
for k=1:NP
temp = func(individual_best(k, :), capacity, weight, V, penality);
if temp > global_best_fit
global_best = individual_best(k, :);
global_best_fit = temp;
end
end
global_optimal(gen) = global_best_fit; % 记录每次迭代中
end
figure(1)
plot(global_optimal);
% 定义适应度函数
function res = func(x, capacity, weight, bag_volume, penality)
% 适应度函数的输入参数
% x: 可行解 二进制串
% capacity: 物品的体积
% bag_volume: 背包的容积
% penality: 惩罚系数
fitness = sum(x.*weight); % 总的价值
total_volume = sum(x.*capacity); % 总的体积
if total_volume <= bag_volume
res = fitness;
else
res = fitness - penality * (total_volume - bag_volume);
end
end
运行结果:
除此之外,还可以使用离散粒子群算法求解单变量的优化问题:
例如,求解函数的极值,函数的表达式如下所示:
其中:
可以画出函数的图像,如下所示:
源代码:
clc;
clear all;
xx = 0:0.05:9;
f_x = xx + 6*sin(4.*xx) + 9*cos(6.*xx);
figure(1)
plot(xx, f_x);
title('f(x)=x+6sin(4x)+9cos(6x)');
xlabel('x');
ylabel('f(x)');
NP = 100; % 种群个数
G = 200; % 迭代次数
D = 10; % 二进制编码的为位数
c1 = 1.5; % 学习因子
c2 = 1.5;
w_max = 0.8; % 惯性权重
w_min = 0.4;
v_max = 5; % 粒子的速度限制
v_min = -5;
x_min = 0;
x_max = 9;
mode = 'min'; % 模式选择,是求函数的最大值函数函数的最小值 'min' , 'max'
% 初始化种群个体
x = (rand(NP,D)>0.5); % 产生均匀分布的二进制串 randn产生的是符合正态分布的随机数
v = v_min + rand(NP,D)*(v_max - v_min); % 速度进行初始化
% 初始化个体最优
individual_best = x; % 每个个体的历史最优
pbest = zeros(NP, 1); % 个体最优位置对应的适应度值
for k=1:NP
pbest(k, 1) = func(individual_best(k, :), x_max, x_min, D);
end
% 初始化全局最优
if strcmp(mode, 'max') % 求函数的最大值
global_best = zeros(1, D);
global_best_fit = eps;
for k=1:NP
temp = func(individual_best(k, :), x_max, x_min, D);
if temp > global_best_fit
global_best = individual_best(k, :);
global_best_fit = temp;
end
end
else % 求函数的最小值
global_best = zeros(1, D);
global_best_fit = inf;
for k=1:NP
temp = func(individual_best(k, :), x_max, x_min, D);
if temp < global_best_fit
global_best = individual_best(k, :);
global_best_fit = temp;
end
end
end
% 开始迭代
for gen=1:G
% 计算动态惯性权值
w = w_max - (w_max-w_min) * gen / G;
for k=1:NP
% 更新速度
v(k, :) = w * v(k, :) + c1 * rand() * (individual_best(k, :) - x(k, :)) + c2 * rand() * (global_best - x(k, :));
% 边界条件处理 % 边界吸收
for t=1:D
if v(k, D) > v_max
v(k, D) = v_max;
end
if v(k, D) < v_min
v(k, D) = v_min;
end
end
% 使用sigmoid函数对速度进行映射
vs(k, :) = 1./(1+exp(-v(k, :)));
% 更新粒子的位置
for t=1:D
if vs(k, t)>rand()
x(k, t) = 1;
else
x(k, t) = 0;
end
end
end
% 更新粒子的历史最有位置
if strcmp(mode, 'min')
% 个体最优位置
for k=1:NP
old_fitness = func(individual_best(k, :), x_max, x_min, D);
new_fitness = func(x(k, :), x_max, x_min, D);
if new_fitness < old_fitness
individual_best(k, :) = x(k, :); % 个体最优位置更新
end
end
% 全局最优位置
for k=1:NP
temp = func(individual_best(k, :), x_max, x_min, D);
if temp < global_best_fit
global_best = individual_best(k, :);
global_best_fit = temp;
end
end
else
% 个体最优位置
for k=1:NP
old_fitness = func(individual_best(k, :), x_max, x_min, D);
new_fitness = func(x(k, :), x_max, x_min, D);
if new_fitness > old_fitness
individual_best(k, :) = x(k, :); % 个体最优位置更新
end
end
% 全局最优位置
for k=1:NP
temp = func(individual_best(k, :), x_max, x_min, D);
if temp > global_best_fit
global_best = individual_best(k, :);
global_best_fit = temp;
end
end
end
% 记录函数的适应度值
fitness_optimal(gen) = global_best_fit;
end
figure(2)
plot(fitness_optimal);
xlabel('x')
ylabel('f(x)');
title(['适应度值 ', num2str(fitness_optimal(end)), ' | ', num2str(x_min+(x_max-x_min)/(2^D-1)*binary2dec(global_best))]);
function res = func(x, xmax, xmin, D)
% 参数 x 二进制串
% vmax : 速度的最大值
% vmin : 速度的最小值
% D : 二进制的编码位数
temp_ = binary2dec(x);
temp = xmin + (xmax-xmin)/(2^D-1) * temp_;
res = temp + 6*sin(4*temp) + 9*cos(6*temp);
end
function res = binary2dec(x)
total_value = 0;
for k=length(x):-1:1
total_value = total_value + x(k)*2^(length(x)-k);
end
res = total_value;
end
求取最大值,将mode参数设置为‘max’即可:
求取最小值,将mode参数设置为‘min’即可: