启发式算法代码

一，模拟退火算法

参考文献：
（1）胡山鹰，陈丙珍，非线性规划问题全局优化的模拟退火法，清华大学学报，1997,37(6):5-9
（2）李歧强，具有约束指导的模拟退火算法，系统工程，2001,19(3):49-55

import random
import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")
plt.rcParams['axes.unicode_minus'] = False #显示负号
plt.rcParams['font.family'] = ['sans-serif']
plt.rcParams['font.sans-serif'] = ['SimHei']  # 散点图标签可以显示中文

# 定义优化问题的目标函数，最小化问题
def Objective_function(X):  # 目标函数和约束条件    最小化
    X = X.flatten() #将X变为一维数组
    x1 = X[0]
    x2 = X[1]
    x3 = X[2]
    p1 = (max(0, 6 * x1 + 5 * x2 - 60)) ** 2     #表达的含义是函数表达式小于0
    p2 = (max(0, 10 * x1 + 20 * x2 - 150)) ** 2
    #如果是等式约束，可以转化成表达式=0，然后目标函数-10000*表达式
    fx = 100.0 * (x2 - x1) ** 2.0 + (1 - x1) ** 2.0 * x3
    return fx + 10000 * (p1 + p2)     #施加惩罚项
def main_work():       #有理数，整数，0-1变量
    sa = SA(x_min=[[-30,-30],[-1],[]],x_max=[[30,30],[1],[]],t1=1000,t0=1e-4,k=0.98,Markov=500,step=0.1)
    fv = sa.optimize()
    plt.title('模拟退火')
    plt.plot(range(1,len(fv)+1), fv, color='r')
    plt.show()
class SA:
    def __init__(self, x_min,x_max,t1,t0,k,Markov,step):
        # ====== 初始化随机数发生器 ======
        randseed = random.randint(1, 100)
        random.seed(randseed)  # 随机数发生器设置种子，也可以设为指定整数

        self.n = len(np.array(x_min).ravel())  # 自变量数量
        self.x_min = x_min  # 给定搜索空间的下限
        self.x_max = x_max  # 给定搜索空间的上限
        self.t1 = t1  # 初始退火温度
        self.t0 = t0  # 终止退火温度，不能是0，因为只会无限趋近于0
        self.k = k  # 降温参数，T(i)=k*T(i-1)
        self.Markov = Markov  # Markov链长度，内循环运行次数
        self.step = step  # 搜索步长

    def optimize(self):
        # ====== 随机产生优化问题的初始解 ======
        xInitial = np.zeros((self.n))  # 初始化，创建数组
        l1 = len(self.x_min[0])
        l2 = len(self.x_min[0]+self.x_min[1])
        for v in range(self.n):
            if v < l1:
                xInitial[v] = np.random.uniform(self.x_min[0][v], self.x_max[0][v])     #有理数
            elif v >= l1 and v < l2:
                xInitial[v] = np.random.randint(self.x_min[1][v-l1], self.x_max[1][v-l1]+1)   #整数
            else:
                xInitial[v] = np.random.randint(0,2)    #0-1变量
        fxInitial = Objective_function(xInitial)

        ##################################      模拟退火算法初始化      ##############################
        xNew = np.zeros((self.n))  # 初始化，创建数组
        xNow = np.zeros((self.n))  # 初始化，创建数组
        xBest = np.zeros((self.n))  # 初始化，创建数组
        xNow[:] = xInitial[:]  # 初始化当前解，将初始解置为当前解
        xBest[:] = xInitial[:]  # 初始化最优解，将当前解置为最优解
        fxNow = fxInitial  # 将初始解的目标函数置为当前值
        fxBest = fxInitial  # 将当前解的目标函数置为最优值
        #print('初始解:{:.6f},{:.6f},\t初始值:{:.6f}'.format(xInitial[0], xInitial[1], fxInitial))

        recordIter = []  # 初始化，外循环次数
        recordFxNow = []  # 初始化，当前解的目标函数值
        recordFxBest = []  # 初始化，最佳解的目标函数值
        recordPBad = []  # 初始化，劣质解的接受概率
        n_Iter = 0  # 外循环迭代次数
        totalMar = 0  # 总计 Markov 链长度
        totalImprove = 0  # fxBest 改善次数
        nMarkov = self.Markov  # 固定长度 Markov链

        ###################################       开始模拟退火优化     ####################################
        # 外循环
        fv = []
        tNow = self.t1  # 当前温度
        while tNow > self.t0:  # 外循环，直到当前温度达到终止温度时结束
            kBetter = 0  # 获得优质解的次数
            kBadAccept = 0  # 接受劣质解的次数
            kBadRefuse = 0  # 拒绝劣质解的次数

            # 内循环
            for k in range(nMarkov):  # 内循环，循环次数为Markov链长度
                totalMar += 1  # 总 Markov链长度计数器

                # ---产生新解
                # 产生新解：通过在当前解附近随机扰动而产生新解，新解必须在 [min,max] 范围内
                # 方案 1：只对 n元变量中的一个进行扰动，其它 n-1个变量保持不变
                xNew[:] = xNow[:]
                v = random.randint(0, self.n - 1)  # 产生 [0,n-1]之间的随机数
                if v < l1:
                    xNew[v] = xNow[v] + self.step * (self.x_max[0][v] - self.x_min[0][v]) * random.normalvariate(0, 1)
                    # random.normalvariate(0, 1)：产生服从均值为0、标准差为 1 的正态分布随机实数
                    xNew[v] = max(min(xNew[v], self.x_max[0][v]), self.x_min[0][v])  # 保证新解在 [min,max] 范围内
                elif v >= l1 and v < l2:
                    xNew[v] = xNow[v] + int(self.step * (self.x_max[1][v-l1] - self.x_min[1][v-l1]) * random.normalvariate(0, 1))
                    # random.normalvariate(0, 1)：产生服从均值为0、标准差为 1 的正态分布随机实数
                    xNew[v] = max(min(xNew[v], self.x_max[1][v-l1]), self.x_min[1][v-l1])  # 保证新解在 [min,max] 范围内
                else:
                    xNew[v] = np.random.randint(0,2)    #0-1变量
                    xNew[v] = max(min(xNew[v], self.x_max[2][v-l2]), self.x_min[2][v-l2])  # 保证新解在 [min,max] 范围内

                # 计算目标函数和能量差
                fxNew = Objective_function(xNew)
                deltaE = fxNew - fxNow

                # 按 Metropolis 准则接受新解
                if fxNew < fxNow:  # 更优解：如果新解的目标函数好于当前解，则接受新解
                    accept = True
                    kBetter += 1
                else:  # 容忍解：如果新解的目标函数比当前解差，则以一定概率接受新解
                    pAccept = np.exp(-np.float64(deltaE) / np.float64(tNow))        # 计算容忍解的状态迁移概率
                    if pAccept > random.random():
                        accept = True  # 接受劣质解
                        kBadAccept += 1
                    else:
                        accept = False  # 拒绝劣质解
                        kBadRefuse += 1

                # 保存新解
                if accept == True:  # 如果接受新解，则将新解保存为当前解
                    xNow[:] = xNew[:]
                    fxNow = fxNew
                    if fxNew < fxBest:  # 如果新解的目标函数好于最优解，则将新解保存为最优解
                        fxBest = fxNew
                        xBest[:] = xNew[:]
                        totalImprove += 1
                        self.step = self.step * 0.99  # 可变搜索步长，逐步减小搜索范围，提高搜索精度

            # 完成当前温度的搜索，保存数据和输出
            pBadAccept = kBadAccept / (kBadAccept + kBadRefuse)  # 劣质解的接受概率
            recordIter.append(n_Iter)  # 当前外循环次数
            recordFxNow.append(round(fxNow, 4))  # 当前解的目标函数值
            recordFxBest.append(round(fxBest, 4))  # 最佳解的目标函数值
            recordPBad.append(round(pBadAccept, 4))  # 最佳解的目标函数值

            #if n_Iter % 10 == 0:  # 模运算，商的余数
                #print('迭代次数:{},温度:{:.3f},最优值:{:.6f}'.format(n_Iter,tNow,fxBest))

            # 缓慢降温至新的温度，降温曲线：T(k)=k*T(k-1)
            tNow = tNow * self.k
            n_Iter = n_Iter + 1
            fxBest = Objective_function(xBest)
            ###############################    结束模拟退火过程    #######################################
        print('--------------模拟退火-------------')
        print('提升次数:{:d}'.format(totalImprove))
        print("求解结果:")
        for i in range(self.n):
            print('\tx[{}] = {:.9f}'.format(i, xBest[i]))
        print('\tf(x):{:.9f}'.format(Objective_function(xBest)))
        return recordFxBest#,fxBest,n_Iter, xBest,  fxNow, recordIter, recordFxNow,  recordPBad

if __name__ == '__main__':
    main_work()

二，粒子群算法

import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")
plt.rcParams['axes.unicode_minus'] = False #显示负号
plt.rcParams['font.family'] = ['sans-serif']
plt.rcParams['font.sans-serif'] = ['SimHei']  # 散点图标签可以显示中文

def Objective_function(X):  # 目标函数和约束条件    最小化
    X = X.flatten() #将X变为一维数组
    x1 = X[0]
    x2 = X[1]
    x3 = X[2]
    p1 = (max(0, 6 * x1 + 5 * x2 - 60)) ** 2     #表达的含义是函数表达式小于0
    p2 = (max(0, 10 * x1 + 20 * x2 - 150)) ** 2
    #如果是等式约束，可以转化成表达式=0，然后目标函数-10000*表达式
    fx = 100.0 * (x2 - x1) ** 2.0 + (1 - x1) ** 2.0 * x3
    return fx + 10000 * (p1 + p2)     #施加惩罚项

def main_work():
    # ( 粒子个数, 最大迭代次数,x_min,x_max, max_vel, 阈值, C1=2, C2=2, W=1)   范围都是左闭右开
    pso = PSO(10, 10000,[[-30,-30],[-1],[]],[[30,30],[1],[]], 30, -1000, C1=1, C2=2, W=1)
    fit_var_list, best_pos = pso.update_ndim()
    print("最优位置:" + str(best_pos))
    print(f"最优解为：{fit_var_list[-1]:.9f}")
    plt.title('粒子群')
    plt.plot([i for i in range(1,len(fit_var_list)+1)],fit_var_list, color='r')
    plt.show()

class particle:
    # 初始化
    def __init__(self, x_min, x_max, max_vel, dim):
        pos = np.zeros((dim))
        self.l1 = len(x_min[0])
        self.l2 = len(x_min[0]+x_min[1])
        for v in range(dim):           #生成初始解
            if v < self.l1:
                pos[v] = np.random.uniform(x_min[0][v], x_max[0][v])     #有理数
            elif v >= self.l1 and v < self.l2:
                pos[v] = np.random.randint(x_min[1][v-self.l1], x_max[1][v-self.l1]+1)   #整数
            else:
                pos[v] = np.random.randint(0,2)    #0-1变量
        self.__pos = np.array(pos)  # 粒子的位置
        self.__vel = np.random.uniform(-max_vel, max_vel, (1, dim))  # 粒子的速度
        self.__bestPos = np.zeros((1, dim))  # 粒子最好的位置
        self.__fitnessValue = Objective_function(self.__pos)  # 适应度函数值
#__开头的为私有属性，只在类内存在
    def set_pos(self, value):
        pos = np.array(value).ravel()
        self.__pos = pos
    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,  size, iter_num, x_min,x_max, max_vel, tol, best_fitness_value=float('Inf'), C1=2, C2=2, W=1):
        self.C1 = C1      #加速常数1，控制局部最优解
        self.C2 = C2      #加速常数2，控制全局最优解
        self.W = W        #惯性因子
        self.dim = len(np.array(x_min).ravel())  # 粒子的维度，变量个数
        self.size = size  # 粒子个数
        self.iter_num = iter_num  # 迭代次数
        self.x_min = x_min    #x 的下限
        self.x_max = x_max     # x 的上限
        self.max_vel = max_vel  # 粒子最大速度
        self.tol = tol  # 截止条件
        self.l1 = len(x_min[0])
        self.l2 = len(x_min[0]+x_min[1])
        self.best_fitness_value = best_fitness_value
        self.best_position = np.zeros((1, self.dim))  # 种群最优位置
        self.fitness_val_list = []  # 每次迭代最优适应值
        # 对种群进行初始化
        self.Particle_list = [particle(self.x_min,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.rand() * (part.get_best_pos() - part.get_pos()) \
                    + self.C2 * np.random.rand() * (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 = np.array(pos_value).ravel()
        for v in range(len(pos_value)):           #生成初始解
            if v < self.l1:
                pos_value[v] = max(min(pos_value[v], self.x_max[0][v]), self.x_min[0][v])  # 保证新解在 [min,max] 范围内
            elif v >= self.l1 and v < self.l2:
                pos_value[v] = max(min(pos_value[v], self.x_max[1][v-self.l1]), self.x_min[1][v-self.l1])  # 保证新解在 [min,max] 范围内
            else:
                pos_value[v] = max(min(pos_value[v], self.x_max[1][v-self.l2]), self.x_min[1][v-self.l2])  # 保证新解在 [min,max] 范围内
        part.set_pos(pos_value)
        value = Objective_function(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())  # 每次迭代完把当前的最优适应度存到列表
            print('第{}次最佳适应值为{}'.format(i, self.get_bestFitnessValue()))#################################################
            if self.get_bestFitnessValue() < self.tol:
                break
        print('--------------粒子群--------------')
        return self.fitness_val_list, self.get_bestPosition()
if __name__ == '__main__':
    main_work()

三，改进的鲸鱼算法

import numpy as np
from numpy import random
from copy import deepcopy
from tqdm import tqdm
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")
np.set_printoptions(threshold=np.inf) # threshold 指定超过多少使用省略号，np.inf代表无限大
np.set_printoptions(suppress=True) #不以科学计数法输出
plt.rcParams['axes.unicode_minus'] = False #显示负号
plt.rcParams['font.family'] = ['sans-serif']
plt.rcParams['font.sans-serif'] = ['SimHei']  # 散点图标签可以显示中文
def fun(X):  # 目标函数和约束条件
    x = X.flatten() #将X变为一维数组
    fx = 0
    for i in range(len(x)-1):
        a = x[i]**2 - 10*np.cos(2*np.pi*x[i]) + 10
        fx += a
    fx += -7*x[-1]
    return fx      #施加惩罚项

s = np.zeros((1,30))
sub = np.array(s-10).ravel()  # 自变量下限
up = np.array(s+10).ravel()  # 自变量上限
type = np.array(s).ravel()    #-1是有理数，0是整数，1是0-1变量

def dd2(best_x, x):  #欧氏距离
    best_x = np.array(best_x)   #转化成numpy数组
    x = np.array(x)          #转化成numpy数组
    c = np.sum(pow(x - best_x, 2), axis=1)    #求方差，在行上的标准差
    d = pow(c, 0.5)   #标准差
    return d
def new_min(arr):  #求最小
    min_data = min(arr)   #找到最小值
    key = np.argmin(arr)  #找到最小值的索引
    return min_data, key
def type_x(xx,type,n):  #变量范围约束
    for v in range(n):
        if type[v] == -1:
            xx[v] = np.maximum(sub[v], xx[v])
            xx[v] = np.minimum(up[v], xx[v])
        elif type[v] == 0:
            xx[v] = np.maximum(sub[v], int(xx[v]))
            xx[v] = np.minimum(up[v], int(xx[v]))
        else:
            xx[v] = np.maximum(sub[v], random.randint(0,2))
            xx[v] = np.minimum(up[v], random.randint(0,2))
    return xx
def woa(sub,up,type,nums,det):
    n = len(sub)  # 自变量个数
    num = nums * n  # 种群大小
    x = np.zeros([num, n])  #生成保存解的矩阵
    f = np.zeros(num)   #生成保存值的矩阵
    for s in range(num):      #随机生成初始解
        for v in range(n):
            rand_data = np.random.uniform(0,1)
            x[s, v] = sub[v] + (up[v] - sub[v]) * rand_data
        x[s, :] = type_x(x[s, :],type,n)
        f[s] = fun(x[s, :])
    best_f, a = new_min(f)  # 记录历史最优值
    best_x = x[a, :]  # 记录历史最优解
    trace = np.array([deepcopy(best_f)]) #记录初始最优值,以便后期添加最优值画图
    ############################ 改进的鲸鱼算法 ################################
    xx = np.zeros([num, n])
    ff = np.zeros(num)
    Mc = (up - sub) * 0.1  # 猎物行动最大范围
    for ii in tqdm(range(det)):      #设置迭代次数，进入迭代过程
        # 猎物躲避,蒙特卡洛模拟，并选择最佳的点作为下一逃跑点 #########！！！创新点
        d = dd2(best_x, x)  #记录当前解与最优解的距离
        d.sort()  #从小到大排序,d[0]恒为0
        z = np.exp(-d[1] / np.mean(Mc))  # 猎物急躁系数
        z = max(z, 0.1)     #决定最终系数
        yx = []  #初始化存储函数值
        dx = []  #初始化存储解
        random_rand = random.random(n) #0-1的随机数
        for i in range(50):    #蒙特卡洛模拟的次数
            m = [random.choice([-1, 1]) for _ in range(n)] #随机的-1和1
            asd = best_x + Mc * z * ((det-ii )/det) * random_rand * m   #最优解更新公式
            xd = type_x(asd,type,n)  #对自变量进行限制
            if i < 1:
                dx = deepcopy(xd)
            else:
                dx = np.vstack((dx,xd))   #存储每一次的解
            yx=np.hstack((yx,fun(xd)))    #存储每一次的值
        best_t, a = new_min(yx)  # 选择最佳逃跑点
        best_c = dx[a, :]   #最佳逃跑点
        if best_t < best_f:   #与鲸鱼算法得到的最优值对比
            best_f = best_t   #更新最优值
            best_x = best_c   #更新最优解
        ############################# 鲸鱼追捕 #################################
        w = (ii / det)**3   #自适应惯性权重!!!创新点
        a = (2 - 2*ii/det)*(1- w)  #a随迭代次数从2非线性下降至0！！！创新点
        pp=0.7 if ii <= 0.5*det else 0.4
        for i in range(num):
            r1 = np.random.rand()  # r1为[0,1]之间的随机数
            r2 = np.random.rand()  # r2为[0,1]之间的随机数
            A = 2 * a * r1 - a
            C = 2 * r2
            b = 1     #螺旋形状系数
            l = np.random.uniform(-1,1)  #参数l
            p = np.random.rand()
            if p < pp:
                if abs(A) >= 1:
                    rand_leader = np.random.randint(0, num)
                    X_rand = x[rand_leader, :]
                    D_X_rand = abs(C * X_rand - x[i, :])
                    xx[i, :] = w*X_rand - A * D_X_rand
                    xx[i, :] = type_x(xx[i, :],type,n) #对自变量进行限制
                elif abs(A) < 1:
                    D_Leader = abs(C * best_x - x[i, :])
                    xx[i, :] = w*best_x - A * D_Leader
                    xx[i, :] = type_x(xx[i, :],type,n) #对自变量进行限制
            elif p >= pp:
                D = abs(best_x - x[i, :])
                xx[i, :] = D*np.exp(b*l)*np.cos(2*np.pi*l) + (1-w)*best_x   #完整的气泡网捕食公式
                xx[i, :] = type_x(xx[i, :],type,n) #对自变量进行限制
            ff[i] = fun(xx[i, :])
            if len(np.unique(ff[:i]))/(i+1) <= 0.1:     #limit阈值 + 随机差分变异！！！创新点
                xx[i,:] = (r1*(best_x-xx[i,:]) +
                           r2*(x[np.random.randint(0,num),:] - xx[i,:]))
                xx[i, :] = type_x(xx[i, :],type,n) #对自变量进行限制
                ff[i] = fun(xx[i, :])
        #将上一代种群与这一代种群以及最优种群结合，选取排名靠前的个体组成新的种群
        F = np.hstack((np.array([best_f]), f, ff))
        F, b = np.sort(F,axis=-1,kind='stable'), np.argsort(F)#按小到大排序,获得靠前的位置
        X = np.vstack(([best_x], x, xx))[b, :]
        f = F[:num]  #新种群的位置
        x = X[:num, :]  #新种群的位置
        best_f, a = new_min(f)  # 记录历史最优值
        best_x = x[a , :]  # 记录历史最优解
        trace = np.hstack((trace, [best_f]))
    return best_x,best_f,trace

best_x,best_f,trace = woa(sub,up,type,20,60)     #种群大小，迭代次数
#种群大小可以为自变量个数，迭代次数看情况
print('最优解为：')
print(best_x)
print('最优值为：')
print(float(best_f))

plt.title('鲸鱼算法')
plt.plot(range(1,len(trace)+1),trace, color='r')
plt.show()