多目标优化问题
问题描述
假设存在五类物品,每类物品中又包含四种具体物品,现要求从这五类物品中分别选择一种物品放入背包中,使得背包内物品的总价值最大,总体积最小,并且背包的总质量不超过92kg。背包问题的数学模型为。
m
a
x
P
X
=
∑
i
=
1
4
P
i
X
maxP_{X} =\sum_{i=1}^{4} P_{i}X
maxPX=i=1∑4PiX
m a x R X = ∑ i = 1 4 R i X maxR_{X} =\sum_{i=1}^{4} R_{i}X maxRX=i=1∑4RiX
s . t . C X X ≤ 92 s.t. C_{X}X\le 92 s.t.CXX≤92
其中
P
X
P_{X}
PX表示背包内物品价值;
R
X
R_{X}
RX表示背包内物品体积,
C
X
C_{X}
CX表示选择物品的质量信息;
X
X
X表示选择物品。
P
P
P为每个物品的价值,
R
R
R为每个物品的体积,
C
C
C为每个物品的质量,
P
,
R
,
C
P,R,C
P,R,C取值如下:
P
=
[
3
4
9
15
2
4
6
8
10
2.5
5
7
10
12
3
3
5
10
10
2
]
,
R
=
[
0.2
0.3
0.4
0.6
0.1
0.25
0.35
0.38
0.45
0.15
0.3
0.37
0.5
0.5
0.2
0.3
0.32
0.45
0.6
0.2
]
P= \begin{bmatrix} &3 &4 &9 &15 &2\\ &4 &6 &8 &10 &2.5\\ &5 &7 &10 &12 &3\\ &3 &5 &10 &10 &2 \end{bmatrix},R= \begin{bmatrix} &0.2 &0.3 &0.4 &0.6 &0.1\\ &0.25 &0.35 &0.38 &0.45 &0.15\\ &0.3 &0.37 &0.5 &0.5 &0.2\\ &0.3 &0.32 &0.45 &0.6 &0.2 \end{bmatrix}
P=
345346759810101510121022.532
,R=
0.20.250.30.30.30.350.370.320.40.380.50.450.60.450.50.60.10.150.20.2
C = [ 10 13 24 32 4 12 15 22 26 5.2 14 18 25 28 6.8 14 14 28 32 6.8 ] C= \begin{bmatrix} &10 &13 &24 &32 &4\\ &12 &15 &22 &26 &5.2\\ &14 &18 &25 &28 &6.8\\ &14 &14 &28 &32 &6.8 \end{bmatrix} C= 1012141413151814242225283226283245.26.86.8
算法流程
上述背包问题有两个优化目标,即背包内的物品的总价值和总体积,约束条件有两条:一是五类物品必须每一类选择一种放入背包;二是背包物品的总质量不超过92。在这种条件下,多个优化目标优化结果可能会是相冲突的,即如果其中一个优化目标好,另外一个优化目标可能会变差。在这里,就可能出现背包内总价值变大,总体积会变大的情况,与优化目标总体积越小越好违背。因此,最终的优化结果是一个非劣解集,可从非劣解集中随机选择一个作为最优解,或者使用其它方法选择最优解。算法流程图如下所示:
代码实现
PSO_MT类设计
#ifndef PSO_MT_H
#define PSO_MT_H
#include <array>
#include <tuple>
#include <vector>
#include <algorithm>
#include <random>
#include <ctime>
#include <iostream>
using namespace std;
using uint = unsigned int;
const double eps = 1e-10;
const double c1 = 0.8;
const double c2 = 0.8;
// 速度权重系数的上下限
const double wmax = 1.2;
const double wmin = 0.1;
// 粒子维度
const uint dim = 5;
struct Particle{
vector<int> solution;
vector<double> speed; // 粒子速度
double value; // 个体价值
double volume; // 个体体积
double weight; // 个体重量
};
class PSO_MT
{
private:
// 5类物品,注意由于初始化array,放在类里面是使用静态类常量,这样编译器编译时就知道
const static uint dim = 5;
// 每类物品可选4种
const static uint item_num = 4;
// 物品价值信息
array<array<double,item_num>,dim> item_value_info{
{
{3, 4, 5, 3},{4, 6, 7, 5},{9, 8, 10, 10},{15, 10, 12, 10},{2, 2.5, 3, 2}
}
};
// 物品体积信息
array<array<double,item_num>,dim> item_volume_info{
{
{0.2, 0.25, 0.3, 0.3},{0.3, 0.35, 0.37, 0.32},
{0.4, 0.38, 0.5, 0.45},{0.6, 0.45, 0.5, 0.6},{0.1, 0.15, 0.2, 0.2}
}
};
// 物品重量信息
array<array<double,item_num>,dim> item_weight_info{
{
{10, 12, 14, 14},{13, 15, 18, 14},{24, 22, 25, 28},
{32, 26, 28, 32},{4, 5.2, 6.8, 6.8}
}
};
// 最大背包容量
const double limit_weight = 92.0;
// 种群数量
uint pop_size;
// 种群
vector<Particle> groups;
// 最大迭代次数
uint maxItems;
// 个体最优集合
vector<Particle> pbest_groups;
// 非劣解集
vector<Particle> pareto_optimal_set;
// 群体最优粒子
Particle gbest_particle;
public:
PSO_MT();
PSO_MT(uint pop_size, uint maxItems);
// 初始化粒子群
void initPopulation();
// 初始化非劣解集,并随机选择一个粒子作为初始的全局最优
void initParetoOptimalSet();
// 算法入口
void pso();
// 更新粒子群速度及位置
void updatePosAndSpeed(const uint k);
// 判断粒子A是否支配粒子B
bool isDominateParticle(const Particle& pa, const Particle& pb);
// 合并非劣解集
void mergeParetoOptimalSet();
// 筛选非劣解集,并去除非劣解集中的重复粒子
void updateParetoOptimalSet(const vector<Particle>& merged_set);
// 从最终的非劣解集中随机选择一个粒子作为最优解
void selectBestByRandom();
// 想从非劣解集中选择价值最大的粒子作为最优粒子
void selectBestByMaxValue();
// 想从非劣解集中选择体积最小的粒子作为最优粒子
void selectBestByMinVolume();
void printBestParticle(string txt, const Particle p);
// 计算粒子适应值
tuple<double, double, double> fitness(vector<int> solution);
};
#endif // PSO_MT_H
初始化粒子群
一个粒子实际上就表示从五类物品中依次选择一个物品放入背包的方案,这个方案就是指粒子的“位置”,同时初始化粒子的“速度”,默认初始化为0,而每个粒子的"位置"实际上又具有相应的适应值,即价值、体积、质量。
void PSO_MT::initPopulation()
{
std::time_t now = std::time(nullptr);
std::tm* cur_tm = std::localtime(&now);
int seed = cur_tm->tm_hour * 3600 + cur_tm->tm_min * 60 + cur_tm->tm_sec;
std::default_random_engine dre(seed);
// 产生1~4的随机整数,表示从当前类种挑选一个物品
std::uniform_int_distribution<int> uid(1, 4);
for(uint i = 0; i < pop_size; i++)
{
vector<int> cur_solution;
for(uint j = 0; j < dim; j++)
{
cur_solution.push_back(uid(dre));
}
// 计算当前解的各个适应值,包括价值、体积、重量
auto fitness_res = fitness(cur_solution);
Particle pcle;
pcle.solution = cur_solution;
pcle.speed = vector<double>(dim, 0.0); // 初速度置为0
pcle.value = get<0>(fitness_res);
pcle.volume = get<1>(fitness_res);
pcle.weight = get<2>(fitness_res);
// 将当前粒子加入粒子群
groups.push_back(pcle);
}
// 初始化时的个体最优即为当前种群的所有个体
pbest_groups = groups;
}
计算粒子的适应值
粒子有三个指标,即价值、体积、质量
tuple<double, double, double> PSO_MT::fitness(vector<int> solution)
{
double value = 0.0;
double volume = 0.0;
double weight = 0.0;
for(uint i = 0; i < dim; i++)
{
int idx = solution.at(i) - 1;
uint uidx = static_cast<uint>(idx);
value += item_value_info[i][uidx];
volume += item_volume_info[i][uidx];
weight += item_weight_info[i][uidx];
}
return make_tuple(value, volume, weight);
}
粒子群移动
前面的文章中提到,每个粒子拥有个体最优粒子,另外种群中有群体最优粒子,于是乎种群中的每个粒子跟随这两个支配粒子移动,更新每个粒子的位置和速度。
void PSO_MT::updatePosAndSpeed(const uint k)
{
std::time_t now = std::time(nullptr);
std::tm* cur_tm = std::localtime(&now);
int seed = cur_tm->tm_hour * 3600 + cur_tm->tm_min * 60 + cur_tm->tm_sec;
std::default_random_engine dre(seed);
// 产生1~4的随机整数,表示从当前类种挑选一个物品
std::uniform_real_distribution<double> ure;
std::uniform_int_distribution<int> uid(1, 4);
double w_value = wmax - (wmax - wmin) * pow(k * 1.0 / maxItems, 2.0);
for(uint i = 0; i < pop_size; ++i)
{
bool is_beyond_bound = false;
// 取出当前粒子的位置和速度信息
vector<int> pcle_solution = groups.at(i).solution;
vector<double> p_speed = groups.at(i).speed;
// 粒子移动后的位置信息
vector<int> pcle_new_solution(dim, 0);
// 更新速度和位置
for(uint j = 0; j < dim; ++j)
{
double p_speed_j = p_speed.at(j);
double p_pos_j = pcle_solution.at(j);
double pbest_pos_j = pbest_groups.at(i).solution.at(j); // 个体最优粒子的位置信息
double gbest_pos_j = gbest_particle.solution.at(j);
double p_speed1_j = w_value * p_speed_j + c1 * ure(dre) * (pbest_pos_j - p_pos_j) +
c2 * ure(dre) * (gbest_pos_j - p_pos_j);
p_speed.at(j) = p_speed1_j;
double cur_pos = p_pos_j + p_speed1_j;
if(cur_pos <= eps || cur_pos > 4.0)
{
is_beyond_bound = true;
break;
}
pcle_new_solution.at(j) = static_cast<int>(ceil(cur_pos)) ;
}
// 如果存在越界情况,则新解随机生成
if(is_beyond_bound)
{
for(uint j = 0; j < dim; ++j)
{
pcle_new_solution.at(j) = uid(dre);
}
}
// 移动后的粒子
Particle moved_pcle;
auto moved_pcle_fitness = fitness(pcle_new_solution);
moved_pcle.value = get<0>(moved_pcle_fitness);
moved_pcle.volume = get<1>(moved_pcle_fitness);
moved_pcle.weight = get<2>(moved_pcle_fitness);
moved_pcle.solution = pcle_new_solution;
moved_pcle.speed = p_speed;
// 粒子i的历史个体最优粒子
Particle pbest_pcle = pbest_groups.at(i);
// 历史个体最优粒子支配新粒子,不作操作
if(isDominateParticle(pbest_pcle, moved_pcle))
{
}
// 新粒子支配历史个体最优粒子
else if(isDominateParticle(moved_pcle, pbest_pcle))
{
pbest_groups.at(i) = moved_pcle;
}
// 两者互不支配
else{
if(ure(dre) > 0.5)
{
pbest_groups.at(i) = moved_pcle;
}
}
groups.at(i) = moved_pcle;
}
}
合并非劣解集
将上一代的非劣解集和当前的个体最优粒子合并,合并为非劣解集。
void PSO_MT::mergeParetoOptimalSet()
{
vector<Particle> merged_set;
// merged_set.resize(pareto_optimal_set.size() + pbest_groups.size());
// error: no match for 'operator<' (operand types are 'Particle' and 'Particle')
// merge(pbest_groups.begin(), pbest_groups.end(),
// pareto_optimal_set.begin(), pareto_optimal_set.end(), merged_set.begin());
move(pbest_groups.cbegin(), pbest_groups.cend(), back_inserter(merged_set));
move(pareto_optimal_set.cbegin(), pareto_optimal_set.cend(), back_inserter(merged_set));
// 筛选并去重非劣解集
updateParetoOptimalSet(merged_set);
}
合并后的非劣解集,在合并集里重新筛选支配粒子,并对支配粒子进行去重操作
void PSO_MT::updateParetoOptimalSet(const vector<Particle> &merged_set)
{
uint set_size = merged_set.size();
vector<Particle> curParetoOptimalSet;
for(uint i = 0; i < set_size; i++)
{
bool flag = false;
Particle pi = merged_set.at(i);
for(uint j = 0; j < set_size; j++)
{
if(i == j)
{
continue;
}
Particle pj = merged_set.at(j);
if(isDominateParticle(pj, pi))
{
flag = true;
break;
}
}
if(flag)
{
continue;
}
curParetoOptimalSet.push_back(pi);
}
// 更新非劣解集,将curParetoOptimalSet中去重后的粒子集合更新为最新的非劣解集
// TODO 2024-05-11 理论上curParetoOptimalSet也不能为空
pareto_optimal_set.clear();
pareto_optimal_set.push_back(curParetoOptimalSet.at(0));
uint cur_set_size = curParetoOptimalSet.size();
if(cur_set_size > 0)
{
for(uint i = 1; i < cur_set_size; i++)
{
bool is_same = false;
Particle pi = curParetoOptimalSet.at(i);
vector<int> solution_i = pi.solution;
// 判断粒子pi是否和已有非劣解重复
for(uint j = 0; j < pareto_optimal_set.size(); j++)
{
vector<int> solution_j = pareto_optimal_set.at(j).solution;
if(solution_i == solution_j)
{
is_same = true;
break;
}
}
if(is_same)
{
continue;
}
pareto_optimal_set.push_back(pi);
}
}
}
判断一个粒子是否支配另一个粒子
- 判断粒子pa是否支配粒子pb,判断依据如下:
- 粒子pb较粒子pa价值小且体积大;
- 粒子pb较粒子pa价值相同但体积大;
- 粒子pb较粒子pa体积相同但价值小
- 粒子pb重量超标
bool PSO_MT::isDominateParticle(const Particle &pa, const Particle &pb)
{
// 满足下面条件,说明粒子pa支配粒子pb
if(pb.weight > limit_weight || (pa.value > pb.value && pa.volume < pb.volume) ||
(abs(pa.value - pb.value) < eps && pa.volume < pb.volume)
|| (abs(pa.volume - pb.volume) < eps && pa.value > pb.value))
{
return true;
}
return false;
}
程序运行结果
随机从非劣解集中选择一个粒子作为最优解(这里运行了多次,可以看到最优解是随机选择的):
best solution:
1 1 2 3 1
value:29,volume:1.48,weight:77
best solution:
1 3 1 1 1
value:36,volume:1.67,weight:88
best solution:
1 4 1 1 1
value:34,volume:1.62,weight:84
选择价值最大的粒子作为最优解:
In max value condition,plan:
3 2 3 1 2
value:38.5,volume:1.9,weight:91.2
选择体积最小的粒子作为最优解:
In min volume condition,plan:
1 4 2 2 1
value:28,volume:1.45,weight:76
从运行结果上来看,背包体积基本随着物品价值变大而变大。如果倾向于价值指标,那么体积指标就会变差;如果倾向于体积指标,那么价值指标就会变差。因此最优解的选择可以根据实际需要制定选择方案,如可以采用加权方法,给价值和体积指标分配权重,变相转化为一个单目标优化问题。