一. 内容简介
多目标粒子群算法python实现(双目标)
二. 软件环境
2.1vsCode
2.2Anaconda
version: conda 22.9.0
三. 代码
3.1 代码
import numpy as np
import matplotlib.pyplot as plt
from itertools import combinations
import warnings
# 主函数
def MOPSO(params, MultiObj):
Np = params['Np']
Nr = params['Nr']
maxgen = params['maxgen']
W = params['W']
C1 = params['C1']
C2 = params['C2']
ngrid = params['ngrid']
maxvel = params['maxvel']
u_mut = params['u_mut']
fun = MultiObj['fun']
# 表示有几个变量，
nVar = MultiObj['nVar']
# 变量的最小值，
var_min = MultiObj['var_min']
# 变量的最大值
var_max = MultiObj['var_max']
# 初始化例子的取值，位置代表变量取值
# 200x3
POS = np.random.rand(Np, nVar) * (var_max - var_min) + var_min
# 每个粒子的每个方向速度
# 200x3
# 前面行，后面列
VEL = np.zeros((Np, nVar))
# 代表每个粒子的适应度，有两个目标的适应度
# 200x2
POS_fit = fun(POS)
# 第一轮都是最优
# 用于记录每个粒子的最优参数
PBEST = POS.copy()
# 用于记录每个粒子的最优适应度
# 适应度要尽可能的小
PBEST_fit = POS_fit.copy()
# 检查支配关系
# DOMINATED里面=1的点，是要被淘汰的，
DOMINATED = checkDomination(POS_fit)
# 精英库，将最优节点添加进精英库
REP = {}
# 精英库参数添加，~DOMINATED = 1 时候，会被添加进去
REP['pos'] = POS[~DOMINATED, :]
# 最优适应度添加
REP['pos_fit'] = POS_fit[~DOMINATED, :]
# REP['pos'] = POS[np.logical_not(DOMINATED), :]
# REP['pos_fit'] = POS_fit[np.logical_not(DOMINATED), :]
# 更新了网格质量
REP = updateGrid(REP, ngrid)
# 求出最大速度
maxvel = (var_max - var_min) * maxvel / 100
# 代数
gen = 1
# 打印代数，以及精英粒子数
print(f"Generation #0 - Repository size: {REP['pos'].shape[0]}")
stopCondition = False
while not stopCondition:
# 选取领头节点
h = selectLeader(REP)
# 更新速度
VEL = W * VEL + C1 * np.random.rand(Np, nVar) * (PBEST - POS) + C2 * np.random.rand(Np, nVar) * (np.tile(REP['pos'][h, :], (Np, 1)) - POS)
# 更新位置
POS = POS + VEL
# 粒子变异
POS = mutation(POS, gen, maxgen, Np, var_max, var_min, nVar, u_mut)
# 超出边界的，对速度和边界做一个修正
POS, VEL = checkBoundaries(POS, VEL, maxvel, var_max, var_min)
# 计算适应度
POS_fit = fun(POS)
# 更新精英库
REP = updateRepository(REP, POS, POS_fit, ngrid)
REP['reduce'] = 0
if REP['pos'].shape[0] > Nr:
# 精英库满了，在比较，删除相对不好的
REP = deleteFromRepository(REP, REP['pos'].shape[0] - Nr, ngrid)
# 返回1说明，说明前面那个点要更优秀
# 也就是说 pos_best 中 1 的参数，是要比PBEST_fit更好的
pos_best = dominates(POS_fit, PBEST_fit)
# dominates返回值中1是可以明确比出好坏的，而0是不可以的
# ~就相当于给1变成0，0变成1，也就是说，best_pos中1是说明比不出来好坏的参数
best_pos = ~dominates(PBEST_fit, POS_fit)
# 那么就让其中接近一半的点都为0，因为比不出来好坏，就部分替换，本轮最优解部分替换之前的最优解
best_pos[np.random.rand(Np) >= 0.5] = False
# pos_best中为1的点是更好的点
if np.sum(pos_best) > 1:
# 先给表现好的点进行更新
PBEST_fit[pos_best, :] = POS_fit[pos_best, :]
PBEST[pos_best, :] = POS[pos_best, :]
if np.sum(best_pos) > 1:
# 给比不出来好坏的点更新
PBEST_fit[best_pos, :] = POS_fit[best_pos, :]
PBEST[best_pos, :] = POS[best_pos, :]
if REP['reduce'] > 0:
# 都是从0开始的
print(f"Generation #{gen} - Repository size: {REP['pos'].shape[0]} - Eliminate size: {REP['eliminate']} - Add size: {REP['add']} - Delete size: {REP['reduce']}")
else:
print(f"Generation #{gen} - Repository size: {REP['pos'].shape[0]} - Eliminate size: {REP['eliminate']} - Add size: {REP['add']} - Delete size: {0}")
gen += 1
if gen > maxgen:
stopCondition = True
h_rep = plt.plot(REP['pos_fit'][:, 0], REP['pos_fit'][:, 1], 'ok')
def updateRepository(REP, POS, POS_fit, ngrid):
# DOMINATED里面=1的点，是要被淘汰的，也就是说1是不好的点
DOMINATED = checkDomination(POS_fit)
# 获取原来精英库的粒子库
# REP['pos'].shape[0]看的这个数组的行数，从1开始
countElite = REP['pos'].shape[0]
countCivilian = POS.shape[0] - sum(DOMINATED)
# 将这些优秀的点添加进入精英库中，axis表示拼接的方向，axis = 0 表示要从行方向拼接
REP['pos'] = np.concatenate((REP['pos'], POS[~DOMINATED, :]), axis=0)
REP['pos_fit'] = np.concatenate((REP['pos_fit'], POS_fit[~DOMINATED, :]), axis=0)
# 在对精英库，进行支配排序
DOMINATED = checkDomination(REP['pos_fit'])
# 求原来库中被淘汰的个数，从0到countElite，但是并不包含countElite，刚好可以给原来的精英库覆盖，因为数组从0开始
REP['eliminate'] = sum(DOMINATED[:countElite])
REP['add'] = countCivilian - sum(DOMINATED[countElite:])
# 用优秀点覆盖原来的点
REP['pos_fit'] = REP['pos_fit'][~DOMINATED, :]
REP['pos'] = REP['pos'][~DOMINATED, :]
# 更新网格质量
REP = updateGrid(REP, ngrid)
# 返回精英库
return REP
def checkBoundaries(POS, VEL, maxvel, var_max, var_min):
# 和malba 中 repmat 类似，功能一样
# np.tile(array, (m, n))
# 相当于array在行方向变为m个
# 相当于array在列方向变成n个
MAXLIM = np.tile(var_max, (POS.shape[0], 1))
MINLIM = np.tile(var_min, (POS.shape[0], 1))
MAXVEL = np.tile(maxvel, (POS.shape[0], 1))
MINVEL = -MAXVEL
# 如果速度过大过小的的，都会拿一个数值限制
VEL[VEL > MAXVEL] = MAXVEL[VEL > MAXVEL]
VEL[VEL < MINVEL] = MINVEL[VEL < MINVEL]
# 碰壁以后，速度变反
VEL[POS > MAXLIM] = -VEL[POS > MAXLIM]
# 修正为边界
POS[POS > MAXLIM] = MAXLIM[POS > MAXLIM]
# 碰壁以后，速度变反
VEL[POS < MINLIM] = -VEL[POS < MINLIM]
# 修正为边界
POS[POS < MINLIM] = MINLIM[POS < MINLIM]
# 返回位置，速度
return POS, VEL
def checkDomination(fitness):
# 粒子个数
Np = fitness.shape[0]
# 用于存放支配关系
dom_vector = np.zeros(Np)
# 排列组合算是
# all_perm = list(combinations(range(1, Np+1), 2))
# matlab从1开始，python从0开始
all_perm = list(combinations(range(0, Np+0), 2))
# 给上下拼接到一起
all_perm = np.concatenate((all_perm, np.flip(all_perm, axis=1)), axis=0)
c = all_perm[:, 0]
# 返回1说明，说明前面那个点要更优秀
d = dominates(fitness[all_perm[:, 0], :], fitness[all_perm[:, 1], :])
# 并不能根据直接的支配关系来找比较好的点，A>B，C>A，这中A就会被淘汰的，所以只能通过d==1找出一定会被淘汰的点
dominated_particles = np.unique(all_perm[d == 1, 1])
# 把要淘汰的点标记为1
dom_vector[dominated_particles] = 1
# 转化为布尔值，不然后边没办法按位取反
dom_vector = dom_vector.astype(bool)
# 返回要淘汰的点
return dom_vector
def dominates(x, y):
# all axis = 1 说明要从行方向来看，1行都是1才返回1
# any axis = 1 说明要从行方向来看，1行有1就返回1
# 这块是要x全部都要小于y,而且x和y还不能说两个维度都相等
# 返回1说明，x这个点要更优秀
result = np.all(x <= y, axis=1) & np.any(x < y, axis=1)
return result
def updateGrid(REP, ngrid):
# 目标适应度个数
ndim = REP['pos_fit'].shape[1]
# 20个网格，21和节点，REP里面的点就是帕利托前沿的点了
REP['hypercube_limits'] = np.zeros((ngrid + 1, ndim))
# 把两个适应度均分ngrid份
for dim in range(ndim):
REP['hypercube_limits'][:, dim] = np.linspace(np.min(REP['pos_fit'][:, dim]), np.max(REP['pos_fit'][:, dim]), ngrid + 1)
# 查看精英库里面的节点个数，0代表行数,1代表列数
npar = REP['pos_fit'].shape[0]
# 记录索引，将下标转换为矩阵的线性索引，
REP['grid_idx'] = np.zeros(npar, dtype=int) # 先行后列
# 网格位置
REP['grid_subidx'] = np.zeros((npar, ndim), dtype=int)
for n in range(npar):
idnames = ""
# 这块算是计算了两个方向的拥挤度了
for d in range(ndim):
# 这是找这个适应度在第几格，这个相比于matlab代码去除了减一，因为python从0开始的
REP['grid_subidx'][n, d] = np.where(REP['pos_fit'][n, d] <= REP['hypercube_limits'][:, d])[0][0]
# 这块是有一个边缘处理的
if REP['grid_subidx'][n, d] == -1:
REP['grid_subidx'][n, d] = 0
# 这个就不用了
# idnames += "," + str(REP['grid_subidx'][n, d])
# REP['grid_idx'][n] = eval(f"np.ravel_multi_index([{idnames[1:]}], [ngrid] * ndim)")
row = REP['grid_subidx'][n, 0]
col = REP['grid_subidx'][n, 1]
if ndim == 2:
array_shape = (ngrid,ngrid)
if ndim == 3:
array_shape = (ngrid,ngrid,ngrid)
# 获得索引位置
REP['grid_idx'][n] = sub2ind(array_shape,row,col)
# 网格质量这块，20不够的，原来的代码可能是只考虑一个维度的网格，后来考虑；两个维度的网格，可能忘记改了
REP['quality'] = np.zeros((ngrid*ngrid, 2))
# 去除重复的，而且还会唯一元素出现的个数，
ids, counts = np.unique(REP['grid_idx'], return_counts=True)
for i in range(len(ids)):
REP['quality'][i, 0] = ids[i]
# 这个是用来统计同一个网格出现的次数，这个可以用其他来替代，
# REP['quality'][i, 1] = 10 / np.sum(REP['grid_idx'] == ids[i])
# 这是我自己修改的
REP['quality'][i, 1] = 10 / counts[i]
return REP
def selectLeader(REP):
# 按概率走的，都能选到，只是说，点少的，更容易选中，有点轮盘赌的味道
prob = np.cumsum(REP['quality'][:, 1])
# 按概率选一个网格
sel_hyp = REP['quality'][np.where(np.random.rand() * np.max(prob) <= prob)[0][0], 0]
# idx是点的顺序
idx = np.arange(len(REP['grid_idx']))
# 把这个网格里面的点选出来
selected = idx[REP['grid_idx'] == sel_hyp]
# 再从选出来的点中，随机取一个
selected = selected[np.random.randint(len(selected))]
# 返回选的点
return selected
def deleteFromRepository(REP, n_extra, ngrid):
REP['reduce'] = n_extra
# 先确定有多少个点
crowding = np.zeros(REP['pos'].shape[0])
for m in range(REP['pos_fit'].shape[1]):
# 对适应度排序，排序是为了算距离
m_fit = REP['pos_fit'][:, m]
# id记录的是原来的位置
idx = np.argsort(m_fit)
# 排序完的，1，2，3，4，这种是对应的距离
m_fit_sorted = m_fit[idx]
# 计算点和相邻两个点的距离,np.inf表示无限大
m_up = np.concatenate((m_fit_sorted[1:], [np.inf]))
m_down = np.concatenate(([np.inf], m_fit_sorted[:-1]))
distance = (m_up - m_down) / (max(m_fit_sorted) - min(m_fit_sorted))
# 再把顺序换回去，idx_back就是每个粒子距离再distance中的位置
idx_back = np.argsort(idx)
# 因为 inf - inf 会变成 NAN 会报错，所以要忽略这个，后边会对这个进行处理
# 忽略 'invalid value encountered in add' 类型的 RuntimeWarning
warnings.filterwarnings('ignore', category=RuntimeWarning, message="invalid value encountered in add")
crowding = crowding + distance[idx_back]
# 恢复 RuntimeWarning
warnings.resetwarnings()
# 给nan变成inf，方便排序
crowding[np.isnan(crowding)] = np.inf
# del_idx = np.argsort(crowding)
# del_idx = del_idx[:n_extra]
# 排序，选出距离最小的点
del_idx = np.argsort(crowding)[:n_extra]
# axis = 0 表示沿行进行，np.delete 是 NumPy 库中的一个函数，用于删除数组中指定的元素。
REP['pos'] = np.delete(REP['pos'], del_idx, axis=0)
REP['pos_fit'] = np.delete(REP['pos_fit'], del_idx, axis=0)
# 更新网格
REP = updateGrid(REP, ngrid)
return REP
def mutation(POS, gen, maxgen, Np, var_max, var_min, nVar, u_mut):
# 给种群数分三份，np.floor 用于取整
fract = Np / 3 - np.floor(Np / 3)
if fract < 0.5:
# np.round 四舍五入 np.ceil向上取整
sub_sizes = [np.ceil(Np / 3), np.round(Np / 3), np.round(Np / 3)]
else:
sub_sizes = [np.round(Np / 3), np.round(Np / 3), np.floor(Np / 3)]
# 用来确定粒子
cum_sizes = np.cumsum(sub_sizes)
# 第二个种群全部变异
nmut = np.round(u_mut * sub_sizes[1])
if nmut > 0:
# 按概率抽选出要变异的粒子，replace=false表示是不放回抽样
idx = cum_sizes[0] + np.random.choice(np.arange(sub_sizes[1]), size=int(nmut), replace=False)
# 给取整，map可以对数组批量处理，list转换结果
idx = list(map(int, idx))
# 对位置进行随机
POS[idx, :] = np.random.rand(int(nmut), nVar) * (var_max - var_min) + var_min
# 第三个种群部分变异，**是指数
per_mut = (1 - gen / maxgen) ** (5 * nVar)
nmut = np.round(per_mut * sub_sizes[2])
if nmut > 0:
idx = cum_sizes[1] + np.random.choice(np.arange(sub_sizes[2]), size=int(nmut), replace=False)
# 给取整，map可以对数组批量处理，list转换结果
idx = list(map(int, idx))
# 对位置进行随机
POS[idx, :] = np.random.rand(int(nmut), nVar) * (var_max - var_min) + var_min
# 返回位置
return POS
def sub2ind(shape, row, col):
# 将行索引和列索引转换为线性索引
return row * shape[1] + col
# 设置参数
params = {
# 粒子数
'Np': 500,
# 精英库
'Nr': 1000,
'maxgen': 200,
'W': 0.4,
'C1': 2,
'C2': 2,
# 网格个数
'ngrid': 20,
'maxvel': 5,
'u_mut': 0.5
}
# 设置多目标函数
MultiObj = {
'fun': lambda x: np.column_stack(
(-10 * (np.exp(-0.2 * np.sqrt(x[:, 0] ** 2 + x[:, 1] ** 2)) + np.exp(-0.2 * np.sqrt(x[:, 1] ** 2 + x[:, 2] ** 2))),
np.sum(np.abs(x) ** 0.8 + 5 * np.sin(x ** 3), axis=1))
),
'nVar': 3,
# 这块是直接处理了
'var_min': np.array([-5, -5, -5]),
'var_max': np.array([5, 5, 5])
}
REP = MOPSO(params, MultiObj)
3.2 效果
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
