因为个人比较熟悉 PyTorch，所以在做优化问题的时候多是用梯度下降法 —— But，梯度下降法不是万能的。对于一些自变量与因变量无直接运算 (加减乘除......) 的时候，不会产生梯度信息，所以这个时候梯度下降法是没有办法优化的。所以迫不得已学了粒子群算法，在实现的过程中加了自己的一些 idea

大致的思路如下：

根据用户所设置的各个坐标的取值范围生成指定规模的粒子群 (位置向量的集合)
粒子群算法类似鱼群觅食，当有一条鱼找到食物的时候，一定范围内的鱼都会朝着食物移动。我在网上看的多是 2 个最优粒子引领群体移动，我这里做的是一定百分比的优粒子引领群体移动
适应度函数由用户定义。使用适应度函数计算出每个粒子的适应度，根据适应度选出一定百分比的优粒子，并根据适应度计算这些优粒子的权值；计算粒子间的欧式距离，结合优粒子的权值生成优粒子的“影响力”。再以影响力为权值，与方向向量 (趋势)、上一轮移动量 (惯性) 加权得到每个粒子的移动量
为保证粒子群搜索范围的广度，在求解的过程中对粒子群的移动量乘上随机比例，使得粒子群在朝着最优解搜索之外还具有一定的随机性

优化器对象

fitness

用户自定义的适应度函数，对于某粒子，适应度越大表示该粒子越好

_new_unit

coord_range 规定了每个坐标的取值上下限，将下限记为 bias，上限 - 下限记为 scale，生成 [0, 1] 的随机数 x，变换为坐标的公式为：

$coord=x\cdot scale+bias$

对该函数指定数量，即可按照该规则生成指定规模的粒子群，用于初始化 / 补全粒子群

_particle_filter

输入布尔矩阵 / 索引矩阵，保留满足条件的粒子，用于重叠检测 (空值检测已移除) 的筛选

_fitness_factor

使用用户自定义的适应度函数计算每一个粒子的适应度 fitness：

$fitness = \{f_0, f_1, ...f_n\}$

使用 np.percentile 函数计算优粒子 (一定百分比) 的适应度边界，若粒子的适应度低于该值则置 0

适应度中最大值对应的粒子为局部最优粒子，与全局最优粒子对比后，选择是否更新全局最优粒子

将全局最优粒子的适应度追加到 fitness，得：

$fitness = \{f_0, f_1, ...f_n, f_{best}\}$

将 fitness 归一化，作为适应度影响因子

考虑到全局最优粒子的影响因子之和可能出现过大的情况 (=2 或远大于其它粒子的影响因子)，所以：

$f_{best} = 1.2-max\{f_0, f_1, ...f_n\}$

该操作削弱了全局最优粒子的影响力，增强了粒子群的局部搜索能力

_move_pace

使用 _fitness_factor 计算适应度影响因子 (行向量)，记为 $\alpha_{fitness}$

计算粒子间的方向向量 direct，并计算欧式距离 dist (为减少计算量，未使用开方)，则距离影响因子 (二维矩阵) 为：

$\alpha_{dist}=\frac{\max\{dist\}-dist}{\max\{dist\}}$

粒子间的影响力矩阵 (相互作用力，为二维矩阵) 为：

$\beta = \alpha_{fitness} \cdot \alpha_{dist}$

因为 $\beta$ 中，有的粒子适应度低于优粒子的适应度边界，所以影响力矩阵会是一个稀疏矩阵

每个粒子的移动量 = $\sum$ 其它粒子对它的影响力 × 其它粒子与它的方向向量

fit

初始化粒子群，初始化粒子惯性
for 循环体：

如果要求坐标是整数，则对粒子群的坐标取整
将越界的坐标变换到距离最近的边界
使用 _move_pace 计算粒子互相影响下产生的移动量 (同时更新全局最优粒子)，在无进展的次数达到一定值时退出
如果在随机搜索的轮次内，则将粒子的移动量乘上 [-1, 1] 的随机比例
更新粒子的坐标：粒子真实移动量 = 粒子互相影响下产生的移动量 + 粒子惯性 (粒子在上一轮的真实移动量) × 惯性权值
筛除重叠粒子，补全粒子群

from typing import Sequence, Optionalimport numpy as np
from tqdm import trangeDTYPE = np.float16def normalize(data, axis=None):''' 归一化axis=0: 各列归一化axis=1: 各行归一化'''data = np.array(data)min_ = data.min(axis=axis, keepdims=True)max_ = data.max(axis=axis, keepdims=True)if min_ == max_:min_ = 0return (data - min_) / (max_ - min_)class Particle_Swarm_Optimizer:''' 粒子群优化器particle_size: 粒子群规模coord_range: 每个坐标的取值范围integer: 坐标是否整数well_percent: 优粒子百分比learn_rate: 学习率best_unit: 已知最优个体'''def __init__(self, particle_size: int,coord_range: Sequence[Sequence],integer: bool = False,well_percent: float = 0.05,learn_rate: float = 1.,best_unit: Optional[Sequence] = None):assert particle_size * well_percent > 2, '优粒子百分比 / 群体规模 过小'# 记录系统参数self._particle_size = particle_sizeself._coord_range = np.array(coord_range, dtype=DTYPE)self._coord_scale = self._coord_range[:, 1] - self._coord_range[:, 0]self._integer = integerself._well_percent = (1 - well_percent) * 100self._learn_rate = learn_rate# 记录最优个体if best_unit:self.best_unit = np.array(best_unit, DTYPE)self.best_fitness = self.fitness(self.best_unit)else:self.best_unit = Noneself.best_fitness = - np.inf# 记录不产生最优个体的次数self._angry = 0def fitness(self, particle):''' 适应度函数 (max -> best)'''raise NotImplementedErrordef fit(self, epochs: int,patience: int = 50,inertia_weight: float = 0.2,random_epochs_percent: float = 0.2,prefix='PSO_fit'):''' epochs: 训练轮次patience: 允许搜索无进展的次数inertia_weight: 惯性权值random_epochs_percent: 随机搜索轮次百分比'''# 生成粒子群self.particle = self._new_unit(self._particle_size)self.inertia = np.zeros_like(self.particle)# 随机搜索轮次数random_epochs = int(random_epochs_percent * epochs)pbar = trange(epochs)for epoch in pbar:# 取整操作self.particle = np.round(self.particle, 0) if self._integer else self.particle# 越界处理for coord_idx, (amin, amax) in enumerate(self._coord_range):self.particle[:, coord_idx] = np.clip(self.particle[:, coord_idx],a_min=amin, a_max=amax)# 粒子互相影响下产生的移动量move_pace = self._move_pace()# 收敛检测if self._angry == patience: break# 根据轮次生成随机移动量if epoch < random_epochs:# 随机比例: [-1, 1]move_pace *= 2 * np.random.random(self.particle.shape) - 1self.particle += (move_pace + inertia_weight * self.inertia) * self._learn_rateself.inertia = move_pace# 重叠检测_, unique = np.unique(self.particle, return_index=True, axis=0)self._particle_filter(unique)# 群体补全need = self._particle_size - len(self.particle)if need:self.particle = np.concatenate([self.particle, self._new_unit(need)], axis=0)self.inertia = np.concatenate([self.inertia, np.zeros([need, len(self._coord_range)])])# 展示进度pbar.set_description((f'%-10s' + '%-10.4g') % (prefix, self.best_fitness))pbar.close()return self.best_unitdef _new_unit(self, num: int):''' 产生指定规模的群体'''x = np.random.random([num, len(self._coord_range)])group = x * self._coord_scale + self._coord_range[:, 0]return groupdef _particle_filter(self, cond: np.array):''' 粒子筛选'''self.particle = self.particle[cond]self.inertia = self.inertia[cond]def _fitness_factor(self):''' 适应度因子'''fitness = np.array(list(map(self.fitness, self.particle)), dtype=DTYPE)# 只保留优粒子的适应度well_bound = np.percentile(fitness, self._well_percent)fitness *= fitness >= well_bound# 局部最优的个体cur_best_index = fitness.argmax()cur_best_fitness = fitness[cur_best_index]# 更新全局最优的个体if cur_best_fitness > self.best_fitness:self._angry = 0self.best_fitness = cur_best_fitnessself.best_unit = self.particle[cur_best_index].copy()else:self._angry += 1# 转化为与最优适应度的比例fitness = normalize(np.append(fitness, self.best_fitness))# 削弱最优适应度的影响力fitness[-1] = 1.1 - fitness[:-1].max()return fitnessdef _move_pace(self):''' 根据 距离、适应度 产生的移动量'''# 适应度因子fitness_factor = self._fitness_factor()# 粒子间的距离refer = np.append(self.particle, self.best_unit[None], axis=0)direct = refer[None] - self.particle[:, None]dist = ((direct / self._coord_scale) ** 2).sum(axis=-1)# 归一化: 距离因子dist_max = dist.max(axis=1, keepdims=True)dist_factor = (dist_max - dist) / dist_max# 粒子间的影响力influence = fitness_factor * dist_factormove_pace = (direct * influence[..., None]).sum(axis=1)return move_pace

求解示例

自变量： $[x_1, x_2, x_3]$

适应度函数： $\sin{x_1}+\sin{x_2}+\sin{x_3}$

其中： $0\leq x_1\leq 2, 2\leq x_2\leq 4, 4\leq x_3\leq 7$

if __name__ == '__main__':import matplotlib.pyplot as plt# 定义 3 个自变量的范围COORD_RANGE = [0, 2], [2, 4], [4, 7]# 绘制正弦函数t = np.linspace(0, 7, 100)plt.plot(t, np.sin(t), color='deepskyblue')# 绘制自变量边界for bound in set(sum(COORD_RANGE, [])):plt.plot([bound] * 2, [-1, 1], color='aqua', linestyle='--')class My_PSO(Particle_Swarm_Optimizer):def fitness(self, particle):return np.sin(particle).sum()# 重写粒子群优化器, 并初始化pso = My_PSO(50, coord_range=COORD_RANGE, integer=False)best = pso.fit(100)print(best)# 绘制最优解plt.scatter(best, np.sin(best), marker='p', c='orange')plt.show()

最优解：[1.570796, 2.000000, 7.000000]

模型求解：[1.570619, 2.000000, 7.000000]

Python 粒子群算法 PSO相关推荐

粒子群优化算法和python代码_Python编程实现粒子群算法(PSO)详解
1 原理粒子群算法是群智能一种,是基于对鸟群觅食行为的研究和模拟而来的.假设在鸟群觅食范围,只在一个地方有食物,所有鸟儿看不到食物(不知道食物的具体位置),但是能闻到食物的味道(能知道食物距离自己位 ...
粒子群算法(PSO)Matlab实现(两种解法)
粒子群算法(PSO) 用途:可以用于寻求最优解问题生物机理:鸟群寻找湖泊在函数中,有很多是无法求出最优解的在这时,我们会采用软计算方法,而PSO算法,在软计算算法中有重要的地位: 好吧,这个仁者 ...
C语言实现粒子群算法(PSO)一
C语言实现粒子群算法(PSO)一最近在温习C语言,看的书是<C primer Plus>,忽然想起来以前在参加数学建模的时候,用过的一些智能算法,比如遗传算法.粒子群算法.蚁群算法等等. ...
【老生谈算法】标准粒子群算法(PSO)及其Matlab程序和常见改进算法——粒子群算法
1.算法详解: 1.原文下载: 本算法原文如下,有需要的朋友可以点击进行下载序号原文(点击下载) 本项目原文 [老生谈算法]标准粒子群算法(PSO)及其Matlab程序和常见改进算法.docx 2 ...
粒子群算法(PSO)以及Matlab实现
粒子群算法(PSO)以及Matlab实现算法背景粒子群优化算法(PSO)是一种进化计算技术(evolutionary computation),1995 年由Eberhart 博士和kennedy ...
粒子群算法(PSO)的C++实现
粒子群算法(PSO)的C++实现粒子群算法(PSO----Particle Swarm Optimization)是常用的智能算法之一,它模拟了鸟群觅食行为,是一种具有随机性的仿生算法 .PS ...
粒子群算法(PSO)初识
粒子群算法PSO是模拟群体智能所建立起来的一种优化算法,用于解决各种优化问题. 抽象问题实例化: 假设一群鸟在觅食,只有一个地方有食物,所有鸟儿都看不见食物(不知道食物的具体位置,知道了就不无需觅 ...
Python编程实现粒子群算法(PSO)详解
1 原理粒子群算法是群智能一种,是基于对鸟群觅食行为的研究和模拟而来的.假设在鸟群觅食范围,只在一个地方有食物,所有鸟儿看不到食物(不知道食物的具体位置),但是能闻到食物的味道(能知道食物距离自己位 ...
【ELM预测】基于粒子群算法PSO优化极限学习机预测含Matlab源码
1 模型为了提高空气质量预测精度,提出一种基于粒子群算法优化极限学习机的空气质量预测模型.运用粒子群算法优化极限学习机的初始权值和偏置,在保证预测误差最小的情况下实现空气质量最优预测.选择平均绝对百 ...
【回归预测-ELM预测】基于粒子群算法PSO优化极限学习机预测附matlab代码
1 内容介绍风电功率预测为电网规划提供重要的依据,研究风电功率预测方法对确保电网在安全稳定运行下接纳更多的风电具有重要的意义.针对极限学习机(ELM)回归模型预测结果受输入参数影响的问题,现将粒子群 ...

Python 粒子群算法 PSO