粒子群PSO优化算法学习笔记及其python实现（附讲解如何使用python语言sko.PSO工具包）

本文首先讲解粒子群算法原理，接着讲解python粒子群工具包，最后通过一个函数极大值来手写实现粒子群算法（不调用工具包）。

总裁余(余登武)

27189人浏览 · 2020-07-02 21:13:51

总裁余(余登武) · 2020-07-02 21:13:51 发布

算法描述

粒子群算法思想来源于实际生活中鸟捕食的过程。假设在一个n维的空间中，有一群鸟（m只）在捕食，食物位于n维空间的某个点上，对于第i只鸟某一时刻来说，有两个向量描述，一个是鸟的位置向量，第二个是鸟的速度。假设鸟能够判断一个位置的好坏，所谓“好坏”，就是离食物更近了还是更远了。鸟在捕食的过程中会根据自己的经验以及鸟群中的其他鸟的位置决定自己的速度，根据当前的位置和速度，可以得到下一刻的位置，这样每只鸟通过向自己和鸟群学习不断的更新自己的速度位置，最终找到食物，或者离食物足够近的点。更新速度和位置的表达式如下。

在这里插入图片描述
其效果图见下面：

sko.PSO 工具包讲解

python语言
首先要下载这个工具包。
这个anaconda下载大家都会
参数详解
如代码：
pso = PSO(func=demo_func, dim=3, pop=40, max_iter=150, lb=[0, -1, 0.5], ub=[1, 1, 1], w=0.8, c1=0.5, c2=0.5)

参数	说明
func	类型function，所要求得优化函数
dim	类型int，维数，即函数的参数数
pop	类型int,种群的大小，也就是粒子的数量。我们使用“pop”来与GA保持一致。默认40
max_iter	类型int,iter迭代的最大值默认150
lb	类型列表,下限。每个参数的下限
ub	类型列表,上限。每个参数的上限
W	对应公式里的惯性权重，默认0.8
C1	学习因子1，默认0.5
C2	学习因子2，默认0.5

属性	说明
pbest_x	array_like, shape is (pop,dim)历史上每个粒子的最佳位置
pbest_y	array_like, shape is (pop,1)历史上最好的粒子图像
gbest_x	array_like, shape is (1,dim)general best location for all particles in history
gbest_y	float历史上所有粒子的最佳图像
gbest_y_hist	list每个迭代的gbest_y

算例：有限制的粒子群

来源于官方文档例子

第一步，定义问题

def demo_func(x):
    x1, x2, x3 = x
    return x1 ** 2 + (x2 - 0.05) ** 2 + x3 ** 2

第二步，做粒子群算法

from sko.PSO import PSO

pso = PSO(func=demo_func, dim=3, pop=40, max_iter=150, lb=[0, -1, 0.5], ub=[1, 1, 1], w=0.8, c1=0.5, c2=0.5)
pso.run()
print('best_x is ', pso.gbest_x, 'best_y is', pso.gbest_y)

结果：
在这里插入图片描述
可以发现x的值落在啦限制区间内。
第三步，画出结果

import matplotlib.pyplot as plt

plt.plot(pso.gbest_y_hist)
plt.show()

在这里插入图片描述

算例：无限制的粒子群

def demo_func(x):
    x1, x2, x3 = x
    return x1 ** 2 + (x2 - 0.05) ** 2 + x3 ** 2


# %% Do PSO
from sko.PSO import PSO

pso = PSO(func=demo_func, dim=3)
pso.run()
print('best_x is ', pso.gbest_x, 'best_y is', pso.gbest_y)

# %% Plot the result
import matplotlib.pyplot as plt

plt.plot(pso.gbest_y_hist)
plt.show()

结果：
在这里插入图片描述

可以发现x1的值无限接近0，x2无限接近0.05,x3的值无限接近0.

所谓优化，我的理解是对一个问题求出它足够好的解，即使这个解不是最优解。如题中解接近0而不是0。在现实生活中，这个解满足要求。

数据批量做粒子群优化

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# @Author: yudengwu
# @Date  : 2020/8/18
def demo_func(x):
    x1, x2, x3 = x
    return x1 ** 2 + (x2 - 0.05) ** 2 + x3 ** 2

from sko.PSO import PSO
import numpy as np
import pandas as pd



data={'lb':[[0,-1,0.5],[1,1,1],[2,3,4]],'ub':[[1,1,1],[2,2,2],[4,5,6]]}
data=pd.DataFrame(data)

print(data.shape[0])

def pso(lb,ub):
    pso = PSO(func=demo_func, dim=3, pop=40, max_iter=150, lb=lb, ub=ub, w=0.8, c1=0.5, c2=0.5)
    pso.run()
    print('best_x is ', pso.gbest_x, 'best_y is', pso.gbest_y)


for i in range(data.shape[0]):
    pso(data['lb'][i], data['ub'][i])

结果

python 实现粒子群

这是一个简单算例，通过该算例简单感受下粒子群
求取函数x + 16 * np.sin(5 * x) + 10 * np.cos(4 * x) 的最大值


import numpy as np
import matplotlib.pyplot as plt
# 粒子（鸟）
class particle:
    def __init__(self):
        self.pos = 0  # 粒子当前位置
        self.speed = 0
        self.pbest = 0  # 粒子历史最好位置


class PSO:
    def __init__(self):
        self.w = 0.5  # 惯性因子
        self.c1 = 1  # 自我认知学习因子
        self.c2 = 1  # 社会认知学习因子
        self.gbest = 0  # 种群当前最好位置
        self.N = 20  # 种群中粒子数量
        self.POP = []  # 种群
        self.iter_N = 100  # 迭代次数

    # 适应度值计算函数
    def fitness(self, x):
        return x + 16 * np.sin(5 * x) + 10 * np.cos(4 * x)

    # 找到全局最优解
    def g_best(self, pop):
        for bird in pop:
            if bird.fitness > self.fitness(self.gbest):
                self.gbest = bird.pos

    # 初始化种群
    def initPopulation(self, pop, N):
        for i in range(N):
            bird = particle()#初始化鸟
            bird.pos = np.random.uniform(-10, 10)#均匀分布
            bird.fitness = self.fitness(bird.pos)
            bird.pbest = bird.fitness
            pop.append(bird)

        # 找到种群中的最优位置
        self.g_best(pop)

    # 更新速度和位置
    def update(self, pop):
        for bird in pop:
            # 速度更新
            speed = self.w * bird.speed + self.c1 * np.random.random() * (
                bird.pbest - bird.pos) + self.c2 * np.random.random() * (
                self.gbest - bird.pos)

            # 位置更新
            pos = bird.pos + speed

            if -10 < pos < 10: # 必须在搜索空间内
                bird.pos = pos
                bird.speed = speed
                # 更新适应度
                bird.fitness = self.fitness(bird.pos)

                # 是否需要更新本粒子历史最好位置
                if bird.fitness > self.fitness(bird.pbest):
                    bird.pbest = bird.pos

    # 最终执行
    def implement(self):
        # 初始化种群
        self.initPopulation(self.POP, self.N)

        # 迭代
        for i in range(self.iter_N):
            # 更新速度和位置
            self.update(self.POP)
            # 更新种群中最好位置
            self.g_best(self.POP)


pso = PSO()
pso.implement()

best_x=0
best_y=0
for ind in pso.POP:
    #print("x=", ind.pos, "f(x)=", ind.fitness)
    if ind.fitness>best_y:
        best_y=ind.fitness
        best_x=ind.pos
print(best_y)
print(best_x)




x = np.linspace(-10, 10, 100000)


def fun(x):
    return x + 16 * np.sin(5 * x) + 10 * np.cos(4 * x)
y=fun(x)
plt.plot(x, y)

plt.scatter(best_x,best_y,c='r',label='best point')
plt.legend()
plt.show()