python实现粒子群优化算法

粒子群优化算法（Particle Swarm Optimization，PSO）属于进化算法的一种，通过迭代寻找最优解。

PSO算法包含以下部分：
价值函数y = f ( x , a ) y=f(x,a)y=f(x,a)，其中a aa为模型参数，x xx为待求变量，我们希望求得令y yy最小的x xx
N NN个粒子，每个粒子为一个矢量，其维数与待求变量一致
每个粒子有两个属性，速度v vv和位置x xx
粒子最优位置l b e s t lbestlbest，每个粒子有一个粒子最优位置，为其历史上到达过的使价值函数最低的位置
全局最优位置g b e s t gbestgbest，为所有粒子到达过的使价值函数最低的位置

PSO算法的基本流程如下：

随机初始化N个粒子的位置、速度
对每次迭代：
	更新各粒子最优位置、最优值
	更新全局最优位置、最优值
	检查是否满足结束条件
	对每个粒子：
		v = C * v + M_1 * R_1 * (gbest - x) + M_2 * R_2 * (lbest - x) 
		x = x + v		

python实现PSO粒子群算法的完整程序：

import numpy as np
from inspect import isfunction
import csv


class Particle:
    def __init__(self, dimension, velocity_limit, value_limit):
        self.position = np.array(value_limit[0]) + \
                        np.array((value_limit[1] - value_limit[0]) * np.random.random(size=dimension))
        self.velocity = 2 * velocity_limit * np.random.random(size=dimension) - velocity_limit
        self.local_best_position = np.array([0.0 for i in range(dimension)])
        self.local_best_value = np.array(10000)
        self.dimension = dimension
        self.velocity_limit = velocity_limit
        self.value_limit = value_limit

    def update_position(self):
        assert len(self.position) == len(self.velocity)
        self.position = np.array(self.position) + np.array(self.velocity)
        for i in range(self.dimension):
            if self.position[i] > self.value_limit[1]:
                self.position[i] = self.value_limit[1]
            elif self.position[i] < self.value_limit[0]:
                self.position[i] = self.value_limit[0]
        return np.array(self.position)

    def update_best(self, value_function, function_parameter):
        value = value_function(self.position, function_parameter)
        if value < self.local_best_value:
            self.local_best_position = self.position
            self.local_best_value = value

    def set_velocity(self, velocity):
        assert len(velocity) == self.dimension
        for i in range(self.dimension):
            if velocity[i] > self.velocity_limit:
                velocity[i] = self.velocity_limit
            elif velocity[i] < -self.velocity_limit:
                velocity[i] = -self.velocity_limit
        self.velocity = velocity


'''
update_parameter=[C,w1,w2], dimension=每个粒子维度
particle_num=粒子数, value_function=价值函数，形式为f=(x, a)，a为价值函数参数，价值越低越好
function_parameter=价值函数的参数, velocity_limit=粒子速度限制，为标量，整形或浮点数，默认0.1
position_limit=粒子位置限制，为二维矢量[xmin, xmax]，默认[0, 1], print_interval=打印中间结果的间隔
evaluate_limit=最大迭代次数，默认10000, value_limit=最低价值，得到使f低于此值得x时停止迭代，默认1.0
'''
class PSO:
    def __init__(self, update_parameter, dimension, particle_num, value_function, function_parameter, velocity_limit=None,
                 position_limit=None, print_interval=100, evaluate_limit=10000, value_limit=1.0):
        assert len(update_parameter) == 3
        assert isinstance(dimension, int)
        assert isinstance(particle_num, int)
        assert isfunction(value_function)
        if velocity_limit is None:
            velocity_limit = 0.1
        assert isinstance(velocity_limit, int) or isinstance(velocity_limit, float)
        if position_limit is None:
            position_limit = [0, 1]
        assert len(position_limit) == 2
        assert position_limit[0] < position_limit[1]
        assert isinstance(print_interval, int)
        assert isinstance(evaluate_limit, int)
        assert isinstance(value_limit, int) or isinstance(value_limit, float)
        self.update_parameter = np.array(update_parameter)
        self.dimension = np.array(dimension)
        self.particle_num = particle_num
        self.function = value_function
        self.print_interval = print_interval
        self.function_parameter = function_parameter
        self.velocity_limit = np.array(velocity_limit)
        self.position_limit = np.array(position_limit)
        self.global_best_position = np.array([0.0 for i in range(dimension)])
        self.global_best_value = np.array(10000)
        self.particle_list = [Particle(self.dimension, self.velocity_limit, self.position_limit)
                              for i in range(self.particle_num)]
        self.evaluate_limit = evaluate_limit
        self.value_limit = value_limit
        print('pso initialized')

    def update_particle(self, function_parameter=None):
        if function_parameter is None:
            function_parameter = self.function_parameter
        for particle in self.particle_list:
            particle.set_velocity(self.update_parameter[0] * particle.velocity +
                                  self.update_parameter[1] * np.random.random() *
                                  (self.global_best_position - particle.position) +
                                  self.update_parameter[2] * np.random.random() *
                                  (particle.local_best_position - particle.position))

            particle.update_position()
            particle.update_best(self.function, function_parameter)
        for particle in self.particle_list:
            if particle.local_best_value < self.global_best_value:
                self.global_best_position = particle.local_best_position
                self.global_best_value = particle.local_best_value

    def train(self):
        for i in range(self.evaluate_limit):
            if self.random_parameter:
                num = np.random.randint(0, len(self.function_parameter) - 1)
                self.update_particle(self.function_parameter[num])
            else:
                self.update_particle()
            if self.global_best_value < self.value_limit:
                print('iteration %d, global best value = %f' % (i, float(self.global_best_value)))
                break
            if i % self.print_interval == 0:
                print('iteration %d, global best value = %f' % (i, float(self.global_best_value)))
        print('pso training has ended')
        return self.global_best_position