一种混合均值中心反向学习粒子群优化算法的python程序实现

一种混合均值中心反向学习粒子群优化算法的python程序实现

最优控制作业，看了一篇优化粒子群优化算法的论文，用python跑了一下效果还可以。

主体程序程序

import numpy as np
import math
import random
import matplotlib.pyplot as plt
import function


class Particles:
    def __init__(self, num, d, iterMax, c1, c2, wMax, wMin, xMax, xMin, vMax, vMin, x=0, v=0):
        self.num = num  # 粒子个数
        self.d = d  # 粒子维数
        self.iterMax = iterMax  # 最大迭代次数
        self.c1 = c1  # 学习因子
        self.c2 = c2  # 学习因子
        self.wMax = wMax  # 惯性权重最大值
        self.wMin = wMin  # 惯性权重最小值
        self.xMax = xMax  # 位置最大值
        self.xMin = xMin  # 位置最小值
        self.vMax = vMax  # 速度最大值
        self.vMin = vMin  # 速度最小值
        self.x = x  # 位置
        self.v = v  # 速度
        self.init()

    def init(self):
        self.x = np.random.rand(self.num, self.d) * (self.xMax - self.xMin) + self.xMin  # 初始化位置
        print(type(self.x))
        self.v = np.random.rand(self.num, self.d) * (self.vMax - self.vMin) + self.vMin  # 初始化速度


particles = Particles(20, 2, 200, 2, 2, 1, 0.4, 20, -20, 1, -1)
p_best = np.ones((particles.num, 1))  # 定义历史最优值
p_Position = np.ones((particles.num, particles.d))  # 定义历史最优位置
g_best = float('inf')  # 定义整个族群的历史最优值
g_Position = np.ones((1, particles.d))  # 定义整个族群的历史最优位置
g_best_curve = np.ones((1, particles.iterMax))  # 定义“Fitness evolution curve”点
max_HMC = np.zeros((1, particles.d))  # 混合均值中心最大范围
min_HMC = np.zeros((1, particles.d))  # 混合均值中心最小范围


def g_best_find():
    """
    族群最优获取
    """
    global g_best
    global g_Position
    for i in range(particles.num):
        if p_best[i] < g_best:
            g_best = p_best[i]
            g_Position = p_Position[i]


def p_best_find():
    """
    最优获取
    """
    global p_best
    global p_Position
    for i in range(particles.num):
        p = function.function0(particles.x[i][0], particles.x[i][1])
        if p < p_best[i]:
            p_best[i] = p
            p_Position[i] = particles.x[i]


def MC():
    """
    均值计算
    @return: 均值
    """
    x_average = np.zeros((1, particles.d))
    for i in range(particles.num):
        x_average = x_average + particles.x[i]
    x_average = x_average / particles.num
    return x_average


def MVF():
    """
    适应值均值计算
    @return: 适应值均值
    """
    fitness_average = 0
    for i in range(particles.num):
        fitness_average = fitness_average + function.function0(particles.x[i][0], particles.x[i][1])
    fitness_average = fitness_average / particles.num
    return fitness_average


def PMC(fitness_average):
    """
    偏均值中心计算
    @param fitness_average: 适应值均值
    @return: 偏均值中心
    """
    x_P_average = np.zeros((1, particles.d))
    n = 0
    for i in range(particles.num):
        if function.function0(particles.x[i][0], particles.x[i][1]) < fitness_average:
            x_P_average = x_P_average + particles.x[i]
            n += 1
    x_P_average = x_P_average / n
    return x_P_average


def HMC(x_average, x_P_average):
    """
    混合均值中心计算
    @param x_average: 均值中心
    @param x_P_average: 偏均值中心
    @return: 混合均值中心
    """
    if function.function0(x_average[0][0], x_average[0][1]) \
            < function.function0(x_P_average[0][0], x_P_average[0][1]):
        return x_P_average
    else:
        return x_average


def OBL_HMC(a, b, HMC):
    """
    反向混合均值中心计算
    @param a: 历史最小混合均值中心
    @param b: 历史最大混合均值中心
    @param HMC: 混合均值中心
    @return: 反向混合均值中心
    """
    return random.random() * (a + b) - HMC


def mutation():
    """
    混合均值中心进化&反向混合均值中心进化
    """
    global g_best
    global g_Position
    global max_HMC
    global min_HMC
    mc = MC()
    mvf = MVF()
    pmc = PMC(mvf)
    hmc = HMC(mc, pmc)
    for i in range(particles.d):
        if max_HMC[0][i] < hmc[0][i]:
            max_HMC[0][i] = hmc[0][i]
        if min_HMC[0][i] > hmc[0][i]:
            min_HMC[0][i] = hmc[0][i]
    g = function.function0(hmc[0][0], hmc[0][1])
    if g < g_best:
        g_best = g
        g_Position = hmc
    obl_hmc = OBL_HMC(min_HMC, max_HMC, hmc)
    g = function.function0(obl_hmc[0][0], obl_hmc[0][1])
    if g < g_best:
        g_best = g
        g_Position = hmc


def main():
    """
    主函数
    @return: None
    """
    global g_best
    global max_HMC
    global min_HMC
    """ 最优初始化"""
    for i in range(particles.num):
        p_best[i] = function.function0(particles.x[i][0], particles.x[i][1])
        p_Position[i] = particles.x[i]
    """  迭代开始  """
    """  迭代初始化  """
    p_best_find()
    g_best_find()
    """ HMC范围初始化 """
    mc = MC()
    mvf = MVF()
    pmc = PMC(mvf)
    hmc = HMC(mc, pmc)
    max_HMC = min_HMC = hmc
    """
    """
    w = 1
    for i in range(particles.iterMax):
        """
        """
        w = 0.25 * math.exp(-0.5 * i / particles.iterMax)
        for j in range(particles.num):
            particles.v[j] = w * particles.v[j] + particles.c1 * random.random() \
                             * (p_Position[j] - particles.x[j]) + particles.c2 \
                             * random.random() * (g_Position - particles.x[j])
            for n in range(particles.d):  # 速度界限判断
                if particles.v[j][n] > particles.vMax:
                    particles.v[j][n] = particles.vMax
                if particles.v[j][n] < particles.vMin:
                    particles.v[j][n] = particles.vMin
            particles.x[j] = particles.x[j] + particles.v[j]  # 粒子位置更新
            """
            """
            for n in range(particles.d):  # 位置界限判断
                if particles.x[j][n] > particles.xMax:
                    particles.x[j][n] = particles.xMax
                if particles.x[j][n] < particles.xMin:
                    particles.x[j][n] = particles.xMin
        p_best_find()
        g_best_find()
        mutation()  # 混合均值中心进化&反向混合均值中心进化
        g_best_curve[0][i] = g_best
    """                                           
    """
    print("函数最优值：" + str(g_best))
    print(g_Position)
    c = np.arange(0, particles.iterMax, 1).tolist()
    gb_curve_list = g_best_curve.tolist()
    plt.plot(c, gb_curve_list[0], "r")
    plt.title("Fitness evolution curve")
    #plt.savefig('./H图/HCOPSO_Himm10.jpg')
    plt.show()


if __name__ == "__main__":
    main()
    # print("hello")

测试函数（以下程序文件取名function.py,把要测试的函数改名为function0()）

import math


def function1(x1, x2):
    """
    Ackley function
    @param x1:
    @param x2:
    @return: y值
    """
    y = (-20 * math.exp(-0.2 * math.sqrt((x1 ** 2 + x2 ** 2) / 2)) -
         math.exp((math.cos(2 * math.pi * x1) + math.cos(2 * math.pi * x2)) / 2)
         + 20 + math.e)
    return y


def function2(x1, x2):
    """
    Himmelblau function
    @param x1:
    @param x2:
    @return: y值
    """
    y = (x1 ** 2 + x2 - 11) ** 2 + (x1 + x2 ** 2 - 7) ** 2
    return y


def function0(x1, x2):
    """
    Rastrigin function
    @param x1:
    @param x2:
    @return: y值
    """
    y = 2 * 10 + x1 ** 2 - 10 * math.cos(2 * math.pi * x1) + x2 ** 2 - 10 * math.cos(2 * math.pi * x2)
    return y