智能优化算法（Ga,PSO,SA）高度模块化（可直接调用）python实现

智能优化算法（Ga,PSO,SA）高度模块化（可直接调用）python实现

为啥做这篇文章

此篇文章基于本人数学建模实验课程的智能算法研究，老师要求用matlab实现一种优化算法分析，并且实现函数封装模块化。本人由于数模需要且比较擅长python，于是在这里记录下算法实现，方便日后比赛利用。

三大智能优化算法本质探究

不要被智能蒙蔽了双眼，其实本质就是思维的转变。
我们之前的传统算法基本是遍历穷举进行可行范围内寻求目标函数的最优解。但在这里有如下几个基本思维步骤（不涉及具体实现）。
1.全局可行解随机选取
2.可行解带入目标函数求值
3.值筛选，选取当前最优值
4.当前最优值对应解进行有目标随机浮动，实现有目标随机遍历
5.2，3，4步骤进行重复。
我觉得这些步骤都看起来很简单，了解了这些基本步骤后对智能算法的原理以及实现也就更加清楚了。各种算法的差异主要体现在可行解的收敛算法。

遗传算法（ga）

遗传算法有个特殊的地方，就是自变量编码问题。我们习惯性将对应维度的值转化成二进制编码然后串联组成DNA序列，然后具体计算再涉及编码解码的问题（其实本人觉得有点多此一举了）。
这里需要将可行解理解为DNA表现型。值筛选就按适应度（与目标函数有关）计算。值得关注的就是可行解的收敛算法了。
GA运用了我们大自然优胜劣汰（就是寻当前最优解），然后逐步扩大当前最优解的涉及影响范围（就是使尽量多的后代与其保持相同的DNA序列）。但这有个通病，就是容易导致可行解逐渐陷入局部最优解，而不是找到我们最需要的解（可以理解局部最优就是陷入极小值点，而全局最优是最小值点，高中数学就告诉我们这两者不一样。当然这是为了方便理解，对于离散非连续的函数可不能这么理解）。
为了解决这个问题，我们就利用变异，染色体交叉的操作来实现可行解跳出舒适圈（实现寻找其他最优解）。提高了遗传算法的感受野，更为容易寻找到目标函数的最优解。

粒子群算法（PSO）

粒子群我个人觉得就可行解的收敛算法与GA不一致罢了。但我更觉得粒子群符合我们的社会行为。
粒子群随机分布，会互相联系，告知对方自己寻找到的最优解（群最优解），以及心知肚明自己碰到的最优解，这些会影响自己的加速度。然后每次迭代会改变自己的速度以及位置（高中物理就告诉我们速度是多维的，速度会改变多个维度变量的数值），最终迭代结束，找到最优解。
该算法可通过画粒子pos图（最好通过PCA将变量加权后进行降维到三维进行绘制）得出收敛的结论（现象是粒子群逐渐汇聚，但偶尔有几个跑去找其他地方的最优解了）。

模拟退火算法（SA）

个人觉得退火挺有意思的。这里没有种群的概念了哈。纯物理现象。
首先温度很高，增加了自变量的感受野（也就是刚开始不会收敛），然后逐步迭代后温度降低，对非最优解的容忍度降低，逐步进入收敛过程。

本文不介绍各个算法具体原理，只是进行收敛算法的探究。

python代码实现

首先给出目标函数

其中n=30(也就是维度)。各个变量的边界条件为[-5.12,5.12]
带入代码中。

首先给出GA代码

import numpy as np
import matplotlib.pyplot as plt

def decode(x, a, b):
    """解码，即基因型到表现型"""
    xt = 0
    for i in range(len(x)):
        xt = xt + x[i] * np.power(2, i)
    return a + xt * (b - a) / (np.power(2, len(x)) - 1)

def decode_X(X: np.array,lb,ub,dna_num):
    """对整个种群的基因解码，上面的decode是对某个染色体的某个变量进行解码"""
    X2 = np.zeros((X.shape[0],dna_num))
    for i in range(X.shape[0]):
        x=[]
        for j in range(dna_num):
            m=decode(X[i,10*j:10*(j+1)],lb,ub)#进行解码
        x.append(m)
        X2[i, :] = np.array(x)
    return X2

def select(X, fitness):
    """根据轮盘赌法选择优秀个体"""
    fitness = 1 / fitness  # fitness越小表示越优秀，被选中的概率越大，做 1/fitness 处理
    fitness = fitness / fitness.sum()  # 归一化    
    idx = np.array(list(range(X.shape[0])))
    X2_idx = np.random.choice(idx, size=X.shape[0], p=fitness)  # 根据概率选择
    X2 = X[X2_idx, :]
    return X2

def crossover(X, c):
    """按顺序选择2个个体以概率c进行交叉操作"""
    for i in range(0, X.shape[0], 2):
        xa = X[i, :]
        xb = X[i + 1, :]
        for j in range(X.shape[1]):
            # 产生0-1区间的均匀分布随机数，判断是否需要进行交叉替换
            if np.random.rand() <= c:
                xa[j], xb[j] = xb[j], xa[j]
        X[i, :] = xa
        X[i + 1, :] = xb
    return X    

def mutation(X, m):
    """变异操作"""
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            if np.random.rand() <= m:
                X[i, j] = (X[i, j] + 1) % 2
    return X
def fitness_func1(X):
    return (X**2-10*np.cos(2*np.pi*X)+10).sum(axis=1)   

def ga(fun_c,c,m,lb,ub,num,dna_num,iter_num):
    """遗传算法主函数"""
    #fun_c:传入的优化函数（传入参数应该为矩阵）
    #c:交叉概率
    #m：变异概率
    #num:种群数量
    #dna_num:dna的数量（表现型）,也就是变量个数
    #lb,ub:变量的上下边界（为矩阵的形式）
    best_fitness = []  # 记录每次迭代的效果
    best_xlist = []
    X0 = np.random.randint(0, 2, (num, 10*dna_num))  # 随机初始化种群，为50*500的0-1矩阵
    for i in range(iter_num):
        X1 = decode_X(X0,lb,ub,dna_num)  # 染色体解码
        fitness = fun_c(X1)  # 计算个体适应度
        X2 = select(X0, fitness)  # 选择操作
        X3 = crossover(X2, c)  # 交叉操作
        X4 = mutation(X3, m)  # 变异操作
        # 计算一轮迭代的效果
        X5 = decode_X(X4,lb,ub,dna_num)
        fitness = fun_c(X5)
        best_fitness.append(fitness.min())
        x= X5[fitness.argmin()]
        best_xlist.append(x)
        X0 = X4
    # 多次迭代后的最终效果
    print("最优值是：%.5f" % best_fitness[-1])
    print("最优解是：",best_xlist[-1])
    plt.plot(best_fitness, color='r')
    plt.show()
ga(fitness_func1,0.1,0.1,-5.12,5.12,50,30,1000)

        

得到的结果：

接下来给出PSO代码：

import numpy as np
import matplotlib.pyplot as plt

def fit_fun(X):  # 适应函数(根据自身情况进行修改)
    return (X**2-10*np.cos(2*np.pi*X)+10).sum(axis=1)


class Particle:
    # 初始化
    def __init__(self, x_max, max_vel, dim):
        self.__pos = np.random.uniform(-x_max, x_max, (1, dim))  # 粒子的位置
        self.__vel = np.random.uniform(-max_vel, max_vel, (1, dim))  # 粒子的速度
        self.__bestPos = np.random.random((1, dim))  # 粒子最好的位置
        self.__fitnessValue = fit_fun(self.__pos)  # 适应度函数值

    def set_pos(self, value):
        self.__pos = value

    def get_pos(self):
        return self.__pos

    def set_best_pos(self, value):
        self.__bestPos = value

    def get_best_pos(self):
        return self.__bestPos

    def set_vel(self, value):
        self.__vel = value

    def get_vel(self):
        return self.__vel

    def set_fitness_value(self, value):
        self.__fitnessValue = value

    def get_fitness_value(self):
        return self.__fitnessValue


class PSO:
    def __init__(self, dim, size, iter_num, x_max, max_vel, tol,  C1, C2, W,best_fitness_value=float('Inf')):
        self.C1 = C1
        self.C2 = C2
        self.W = W      
        self.dim = dim  # 粒子的维度
        self.size = size  # 粒子个数
        self.iter_num = iter_num  # 迭代次数
        self.x_max = x_max
        self.max_vel = max_vel  # 粒子最大速度
        self.tol = tol  # 截至条件
        self.best_fitness_value = best_fitness_value
        self.best_position = np.zeros((1, dim))  # 种群最优位置
        self.fitness_val_list = []  # 每次迭代最优适应值

        # 对种群进行初始化
        self.Particle_list = [Particle(self.x_max, self.max_vel, self.dim) for i in range(self.size)]

    def set_bestFitnessValue(self, value):
        self.best_fitness_value = value

    def get_bestFitnessValue(self):
        return self.best_fitness_value

    def set_bestPosition(self, value):
        self.best_position = value

    def get_bestPosition(self):
        return self.best_position

    # 更新速度
    def update_vel(self, part):
        vel_value = self.W * part.get_vel() + self.C1 * np.random.uniform(-1,1) * (part.get_best_pos() - part.get_pos()) \
                    + self.C2 * np.random.uniform(-1,1) * (self.get_bestPosition() - part.get_pos())
        vel_value[vel_value > self.max_vel] = self.max_vel
        vel_value[vel_value < -self.max_vel] = -self.max_vel
        part.set_vel(vel_value)
    # 更新位置
    def update_pos(self, part):
        pos_value = part.get_pos() + part.get_vel()
        pos_value[pos_value>self.x_max]=self.x_max
        pos_value[pos_value<-self.x_max]=-self.x_max
        part.set_pos(pos_value)
        value = fit_fun(part.get_pos())
        if value < part.get_fitness_value():
            part.set_fitness_value(value)
            part.set_best_pos(pos_value)
        if value < self.get_bestFitnessValue():
            self.set_bestFitnessValue(value)
            self.set_bestPosition(pos_value)

    def update_ndim(self):
        for i in range(self.iter_num):
            for part in self.Particle_list:
                self.update_vel(part)  # 更新速度
                self.update_pos(part)  # 更新位置
            self.fitness_val_list.append(self.get_bestFitnessValue())  # 每次迭代完把当前的最优适应度存到列表
            if self.get_bestFitnessValue() < self.tol:
                break
        return self.fitness_val_list, self.get_bestPosition()

if __name__ == '__main__':
    pso = PSO(30, 20, 100000, 5.12, 1, 1e-12, 2, 5,8)
    fit_var_list, best_pos = pso.update_ndim()
    print("最优位置:" + str(best_pos))
    print("最优解:" + str(min(fit_var_list)))
    plt.plot(range(len(fit_var_list)), fit_var_list, alpha=0.5)
    plt.show()

结果如下：
接下来给出SA代码：

import math
import numpy as np
from random import random
from numpy.random import rand
import matplotlib.pyplot as plt
def func(X):                  #函数优化问题
    X=np.array(X)
    return (X**2-10*np.cos(2*np.pi*X)+10).sum()
class SA:
    def __init__(self, func,lb,ub,iter,dim, T0, Tf, alpha):
        self.func = func         #定义优化函数
        self.iter = iter         #内循环迭代次数
        self.alpha = alpha       #降温系数
        self.T0 = T0             #初始温度T0
        self.Tf = Tf             #温度终值Tf
        self.T = T0              #当前温度(这里只是初始化)
        self.ub=ub               #自变量最大值
        self.lb=lb               #自变量最小值
        self.dim=dim             #自变量维数（也就是数量)
        self.X =np.array([lb+(ub-lb)*rand()  for i in range(dim)
                ])#限定取值区间
        self.most_best =[]
        self.history = {'f': [], 'T': []}

    def generate_new(self, X):   #扰动产生新解的过程
            X=np.array(X)
            X_new=X+self.T*np.random.uniform(-1,1,X.shape)/100
            X_new[X_new<self.lb]=self.lb
            X_new[X_new>self.ub]=self.ub
            return X_new 

    def Metrospolis(self, f, f_new):   #Metropolis准则
        if f_new<= f:
            return 1
        else:
            p = math.exp((f - f_new) / self.T)
            if random() < p:
                return 1
            else:
                return 0

    def best(self):    #获取最优目标函数值
        f = self.func(self.X)
        return f   #f_best,idx分别为在该温度下，迭代L次之后目标函数的最优解和最优解的下标

    def run(self):
        count = 0
        #外循环迭代，当前温度小于终止温度的阈值
        while self.T > self.Tf:       
            for i in range(self.iter): 
                f = self.func(self.X)                    #f为迭代一次后的值
                X_new = self.generate_new(self.X) #产生新解
                f_new = self.func(X_new)                        #产生新值
                if self.Metrospolis(f, f_new):                         #判断是否接受新值
                    self.X=X_new
            # 迭代L次记录在该温度下最优解
            ft = self.best()
            self.history['f'].append(ft)
            self.history['T'].append(self.T)
            #温度按照一定的比例下降（冷却）
            self.T = self.T * self.alpha        
            count += 1
            
            # 得到最优解
        f_best = self.best()
        print(f"F={f_best}, x={self.X}")
sa = SA(func,-5.12,5.12,1000,30,100,0.01,0.98)
sa.run()
plt.plot(sa.history['T'], sa.history['f'])
plt.title('SA')
plt.xlabel('T')
plt.ylabel('f')
plt.gca().invert_xaxis()
plt.show()

结果如下：

不知道为啥哈，SA算法不进人意，大家可自行修改超参数进行完善。我累了。

算法注意事项

1.首先我们需要注意对于变量范围的划分，我们最好首先均匀计算选点绘出离散图然后选择一个有效区间。
2.超参数一定要认真调试！我的最后一个就是超参数调不好了，靠！！！