粒子群算法及其python实现

import numpy as np
import matplotlib.pyplot as plt
import copy


def F(x, y):
    return 3 * (1 - x)**2 * np.exp(-(x**2) - (y + 1)**2) - 10 * ( x / 5 - x**3 - y**5) * np.exp(-x**2 - y**2) - 1 / 3**np.exp(-(x + 1)**2 -y**2)


class PSOIndividual:
    '''
    创建pso种群个体
    '''
    def __init__(self, vardim, bound, vmax):
        self.vardim = vardim
        self.bound = bound
        self.vmax = vmax
        self.fitness = 0.  # 个体适应度

    def generate_pos(self):
        # 随机初始化创建粒子位置
        len_p = self.vardim
        self.position = np.zeros(len_p)
        rnd = np.random.random(size=len_p)
        for i in range(len_p):
            self.position[i] = self.bound[0, i] + (self.bound[1, i] - self.bound[0, i]) * rnd[i]

    def calculateFitness(self):
        x, y = self.position[0], self.position[1]
        self.fitness = F(x, y)


class PSOAlgorithm:
    def __init__(self, sizepop, vardim, bound, vmax, maxiter, params):
        self.sizepop = sizepop      # 种群规模
        self.vardim = vardim        # 决策变量维度
        self.bound = bound          # 决策变量上下界
        self.vmax = vmax            # 粒子最大速度
        self.maxiter = maxiter      # 迭代次数
        self.params = params        # params[w,c1,c2,r]
        self.population = []        # 种群
        self.fitness = np.zeros((self.sizepop, 1))   # 适应度
        self.trace = np.zeros((self.maxiter, 2))     # 历史最优和平均适应度

    # 初始化种群
    def initialize(self):
        for i in range(0,self.sizepop):
            ind = PSOIndividual(self.vardim, self.bound, self.vmax)
            ind.generate_pos()
            self.population.append(ind)

    # 种群适应度
    def evaluate(self):
        for i in range(self.sizepop):
            self.population[i].calculateFitness()
            self.fitness[i] = self.population[i].fitness

    def solve(self):
        self.t = 0  # number of iter
        self.initialize()
        self.evaluate()
        gbest = np.max(self.fitness)                 # 全局最佳
        gbestIndex = np.argmax(self.fitness)
        self.gbest = copy.deepcopy(self.population[gbestIndex])
        self.avefit = np.mean(self.fitness)          # 平均适应度
        self.pbest = copy.deepcopy(self.population)  # 用来记录每个粒子的历史最佳位置
        self.pbestfit = copy.deepcopy(self.fitness)  # 用来记录每个粒子的历史最佳距离
        self.trace[self.t, 0] = self.gbest.fitness
        self.trace[self.t, 1] = self.avefit
        print("Generation %d: optimal function value is: %f; average function value is %f"
              % (self.t, self.trace[self.t, 0], self.trace[self.t, 1]))
        while (self.t < self.maxiter-1):
            self.t += 1
            # 初始化速度
            self.velocity = np.random.random(size=(self.sizepop, self.vardim))*self.vmax
            self.VelocityUpdate()    # 速度更新
            self.PositionUpdate()    # 位置更新
            self.evaluate()
            gbest = np.max(self.fitness)
            gbestIndex = np.argmax(self.fitness)
            if gbest > self.gbest.fitness:
                self.gbest = copy.deepcopy(self.population[gbestIndex])
            self.avefit = np.mean(self.fitness)
            self.trace[self.t, 0] = self.gbest.fitness
            self.trace[self.t, 1] = self.avefit
            print("Generation %d: optimal function value is: %f; average function value is %f"
                  % (self.t, self.trace[self.t, 0], self.trace[self.t, 1]))
        print("Optimal function value is: %f; " % self.trace[self.t, 0])
        print("Optimal solution is:", self.gbest.position)
        self.printResult()

    def VelocityUpdate(self):
        '''
        params[w,c1,c2,r]
        w   惯性权重
        c1  跟踪自身最佳位置的权重，认知
        c2  跟踪种群最佳位置的权重，社会
        r1,r2 [0,1]间的随机数
        r   位置更新时的速度系数，通常为1
        '''
        new_vel=[]
        r1 = np.random.random( size=self.sizepop)
        r2 = np.random.random( size=self.sizepop)
        for i in range(self.sizepop):
            if self.fitness[i] > self.pbestfit[i]:
                # 单个粒子历史最优位置更新
                self.pbest[i].position = self.population[i].position
                self.pbestfit[i] = self.fitness[i]
            # 速度更新公式
            vel = self.velocity[i]*self.params[0] + self.params[1]*r1[i]*(
                self.pbest[i].position-self.population[i].position) + self.params[2]*r2[i] * (self.gbest.position-self.population[i].position)
            # 判断是否超过最大速度
            for j in range(self.vardim):
                if vel[j] > self.vmax:
                    vel[j] = self.vmax
                elif vel[j] < -self.vmax:
                    vel[j] = -self.vmax
            new_vel.append(vel)
        self.velocity=new_vel


    def PositionUpdate(self):
        for i in range(self.sizepop):
            # 位置更新公式
            pos = self.population[i].position + self.params[3]*self.velocity[i]
            # 判断是否超过边界
            for j in range(self.vardim):
                if pos[j] > self.bound[1, j]:
                    pos[j] = self.bound[1, j]
                elif pos[j] < self.bound[0, j]:
                    pos[j] = self.bound[0, j]
            self.population[i].position=pos

    def printResult(self):       
        # plot the result of the pso        
        x = np.arange(0, self.maxiter)
        y1 = self.trace[:, 0]
        y2 = self.trace[:, 1]
        plt.plot(x, y1, 'r', label='optimal value')
        plt.plot(x, y2, 'g', label='average value')
        plt.xlabel("Iteration")
        plt.ylabel("function value")
        plt.title("Particle Swarm Optimization")
        plt.legend()
        plt.show()


if __name__ == "__main__":
    bound = np.tile([[-5], [0]], 2)
    # def __init__(self,sizepop,vardim,bound,vmax,maxiter,params):
    # params[w,c1,c2,r]
    pso = PSOAlgorithm(100, 2, bound, 1, 100, [0.1, 2, 2, 1])
    pso.solve()