现代智能优化算法

1.基本概念

现代智能优化算法本质上是一类算法，是为了解决NP-hard问题而诞生的启发式算法，常见的启发式算法主要包括模拟退火、遗传算法等。

1.1什么是NP-hard问题？

在这之前，必须了解一个概念，叫做多项式时间：在计算复杂度理论中，指的是一个问题的计算时间不大于问题大小的多项式倍数。通俗来说就是一定规模的问题，总有一个时间范围内可以将它解决，这个范围就是多项式时间。
所谓的非确定性是指，可用一定数量的运算去解决多项式时间内可解决的问题。也就是说，我们可以很容易验证答案是否可行，但是无法快速求解，因为NP-hard问题所包含的解集实在是太多，常见的NP-hard问题有TSP旅行商问题。
如：假设一个推销员需要从香港出发，经过广州，北京，上海，…，等 n 个城市， 最后返回香港。 任意两个城市之间都有飞机直达，但票价不等。假设公司只给报销 C 元钱，问是否存在一个行程安排，使得他能遍历所有城市，而且总的路费小于 C？

1.2算法简介

根据目前常用的情况，本文将选择模拟退火和遗传算法进行讲解。为了更方便入门，本文选用实例十分简单。

2.模拟退火

2.1算法简介

模拟退火的基本思想来自于材料统计学，里面涉及到一些物理量，但我们并不需要理解其中物理量的含义。只需要知道如何求解即可。
在高温条件下,粒子的能量较高,可以自由运动和重新排列。在低温条件下,粒子能量较低。如果从高温开始,非常缓慢地降温(这个过程被称为退火),粒子就可以在每个温度下达到热平衡。当系统完全被冷却时,最终形成处于低能状态的晶体。

2.2简单实例理解思路

首先，我们求解一个函数最大值时，有什么手段？
1.求导
2.取间隔为0.01，使用蒙特卡洛模拟1遍历区间
3.启发式算法
这里我们举例如下：
求解函数z = sin((x²+y²)^0.5)在x、y在区间[-5,5]的最小值。这个函数的最大值很容易用肉眼看出来大概范围，这里我们主要是作为示例说明。
绘图有：

具体实现代码如下：

import math
import random
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号


def func(x,y):
    R = math.sqrt(x**2 + y**2)
    Z = math.sin(R)
    return Z#模拟算法是求原函数的最小值的，则将函数反向取负数即可求最大值
class SA:
    def __init__(self,func,n=100,t0=100,td=0.01,alpha=0.99):#alph要合理设置
        self.func=func#将函数放入类属性中，转化为类方法
        self.n=n#内循环次数
        self.t0=t0#初始温度
        self.td=td#最低温度
        self.t=t0#当前温度
        self.alpha=alpha#温度降低系数
        self.x = [random.uniform(-5,5) for i in range(n)]
        self.y = [random.uniform(-5,5) for i in range(n)]#初始化x,y,范围为-5~5
        self.best_f=[]#初始化存放每个温度下的最好值
        self.T=[]#初始化每一个温度
        self.best_xy=[]
    def new_xy(self,x,y):
        while 1:
            new_x=x+self.t*random.uniform(-1,1)#产生新的值，波动范围与温度有关，可适当设置
            new_y=y+self.t*random.uniform(-1,1)
            if (-5<=new_x<=5) and (-5<=new_y<=5):#在定义域中，即可
                break
        return new_x,new_y
    def isretain(self,f_new,f):#这里利用的是Metrospolis原则
        ist=0
        dc=f_new-f
        if dc<=0:#如果新的值小于等于原来值，保留
            ist=1
        else:#如果新的值大于原来的值，让概率进行选择,此时dc>0
            p=math.exp(-dc/self.t)#让值始终处于0-1之间，e的x次方，x小于0
            if p>=random.random():
                ist=1
        return ist
    def best(self):
        f_list=[]
        for i in range(self.n):
            f=self.func(self.x[i],self.y[i])#将该温度下值进行存储，找到最小值
            f_list.append(f)
        f_best=min(f_list)
        idx=f_list.index(f_best)#得到最佳值下标，以便后面找寻自变量
        return f_best,idx
    def run(self):#方法启动器，启动整个方法
        cnt=0#计数外循环次数
        while self.t > self.td:#外循环迭代
            for i in range(self.n):#内循环迭代
                f=self.func(self.x[i],self.y[i])#计算该次内循环的值
                new_x,new_y=self.new_xy(self.x[i],self.y[i])#生成新的xy，分别替代原来100个里面的x和y
                #再次计算新的xy所对应的新的值，以此来给定概率是否保留
                f_new=self.func(new_x,new_y)
                if self.isretain(f_new,f):#如果保留，替换对应位置的x[i],y[i]
                    self.x[i]=new_x
                    self.y[i]=new_y#就这样迭代，让n个样本点逐渐向着新的值靠近
            #记录该温度下的最好成绩，即最大值
            best_f,idx=(self.best())
            best_x=self.x[idx]
            best_y=self.y[idx]#保存最佳值对应的x，y
            self.best_xy.append((best_x,best_y))#以列表中套元祖保存
            self.best_f.append(best_f)
            self.T.append(self.t)#传递当前的最佳值和温度给类
            self.t=self.t*self.alpha#温度按照一定比率下降
            cnt+=1
            if cnt%10==0:
                print(f"已经迭代温度{cnt}次，此时最佳值为{best_f}")
        all_best_f=min(self.best_f)
        all_best_idx=self.best_f.index(all_best_f)
        print(f"总共迭代了{cnt}次,最佳值为{all_best_f},对应的自变量x为{self.best_xy[all_best_idx][0]},y为{self.best_xy[all_best_idx][1]}")

if __name__=='__main__':
    #构建GA类的一个实例对象
    sa=SA(func)#将待求函数发送
    sa.run()
    #画图看看迭代的结果
    plt.plot(sa.T, sa.best_f,)
    plt.xlabel('温度')
    plt.ylabel('函数值')
    plt.gca().invert_xaxis()#反向x轴
    plt.show()

运行结果如下：
已经迭代温度890次，此时最佳值为-0.9999983220396803
已经迭代温度900次，此时最佳值为-0.9999839439032638
已经迭代温度910次，此时最佳值为-0.9999999909611703
总共迭代了917次,最佳值为-0.9999999999928478,对应的自变量x为3.280114346731557,y为3.383414757360287

2.3核心思想及代码解释

核心实现在于初始化温度为100，在不断迭代中使得温度变为0，在这之间不断迭代最优解。
代码步骤：
1.首先利用随机值生成初始解，大小为n自由设置
2.当当前温度大于要求温度0时，循环遍历每个解（x，y）
3.对于每个解，先计算得分（一般为每个解打分，这个打分是自己设置的，如旅游商问题所花总费用，我们这里即为函数值）
4.随后将每个解生成新的解

def new_xy(self,x,y):
        while 1:
            new_x=x+self.t*random.uniform(-1,1)#产生新的值，波动范围与温度有关，可适当设置
            new_y=y+self.t*random.uniform(-1,1)
            if (-5<=new_x<=5) and (-5<=new_y<=5):#在定义域中，即可
                break
        return new_x,new_y

生成条件为当前温度*随机值+原来值，确保在可行域下。这里可以明确，当温度越来越低后，新的值和原来的值的差距会逐渐缩小，这也是启发式算法通用的核心思想：在初期广度搜索，后期在找到可能的最优值后再最优值左右微微的波动。
5.随后再次计算新的得分（此题为函数值），并利用函数判断是否保留新的值，如果保留则将原来值替换，代码如下：

def isretain(self,f_new,f):#这里利用的是Metrospolis原则
        ist=0
        dc=f_new-f
        if dc<=0:#如果新的值小于等于原来值，保留
            ist=1
        else:#如果新的值大于原来的值，让概率进行选择,此时dc>0
            p=math.exp(-dc/self.t)#让值始终处于0-1之间，e的x次方，x小于0
            if p>=random.random():
                ist=1
        return ist

这里的核心思想在于，得分越高的值，保留下来的几率比原来大，这和遗传算法中的遗传实现也很类似，核心在于将更高的值更好的可能下继承下来。
6.在将解集（n大小）的所有值都更新完毕后，开始温度的更新，这里温度的更新格式为self.t=self.t*self.alpha#温度按照一定比率下降，即按照一定比例下降，alpha为自由设定值。

7.反复迭代1-6步，一直到温度最终低于阈值。

2.4模拟退火总结

可见，整个遗传算法的思路并不难，但是这个问题很简单，当遇到一个复杂问题的时候，我们在构建产生新的可行解时会更加困难，可能需要考虑各种各样的限制。同时，得分也会更加难以确定。但但我们理解到算法思想后，便能够更加只有的去设计具体细节，以便更好的收敛。
旅行推销员问题（英语：Travelling salesman problem, TSP）是这样一个问题：给定一系列城市和每对城市之间的距离，求解访问每一座城市一次并回到起始城市的最短回路。对于旅行商问题，我们应该怎么改进我们的算法来求解？
这里我提一下思路，由于篇幅限制暂且不进行求解。
在生成解集时，利用随机数生成满足题意的多条路线，可以使用深度优先+随机生成。在更新解集时，随机更改一条路线下的的某段城市点，但是确保修改后的路线是可行的。在衡量得分的时候，旅行商问题即为直接求和取负即可。

3.遗传算法

3.1算法简介

遗传算法是用于解决最优化问题的一种搜索算法。从名字来看，遗传算法借用了生物学里达尔文的进化理论：”适者生存，不适者淘汰“，将该理论以算法的形式表现出来就是遗传算法的过程。
遗传算法的整体算法思想和模拟退火类似，但是他采用了一个更奇妙的表现方式表达了整个迭代更新。

3.2简单实例理解思路

求函数y = x + 16 * sin(5x) + 10 * cos(4x)的最大值，设置x区间范围为：[-10,10]。
具体实现代码如下：

import random
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号

a=-10
b=10#设置定义域

def func(x):
    y=x + 16 * np.sin(5 * x) + 10 * np.cos(4 * x)
    return y
def ori_popular(num):
    global a,b
    ori_popular=[]
    for i in range(num):
        x = random.uniform(a, b)#在此范围内生成一个随机浮点数
        ori_popular.append(x)
    return ori_popular

def encode(popular):#将十进制小数转化为二进制
    global a,b
    popular_gene = []
    for i in range(0, len(popular)):
        data = int((popular[i]-(a)) / b * 2**18)  # 染色体序列为18bit
        bin_data = bin(data)  # 整形转换成二进制是以字符串的形式存在的
        for j in range(len(bin_data)-2, 18):  # 序列长度不足补0
            bin_data = bin_data[0:2] + '0' + bin_data[2:]
        popular_gene.append(bin_data)
    return popular_gene

def decode(popular_gene):#解码，将二进制转化为十进制
    fitness,xs = [],[]
    global a,b
    for i in range(len(popular_gene)):
        x = (int(popular_gene[i], 2) / 2**18) * b + a
        value = x + 16 * np.sin(5 * x) + 10 * np.cos(4 * x) 
        xs.append(x)
        fitness.append(value)
    return fitness,xs

#这里采用轮盘赌选择方法
def nature_choice(popular_gene):
    scores,xs=decode(popular_gene)#这里将函数值作为得分即可
    for i in range(len(scores)):
        if scores[i]<=0:
            scores[i]=0#剔除负数概率
    sum_score=np.sum(scores)#计算得到100个个体的总分，方面之后进行概率抽取
    #计算各个个体被选中的概率
    prob=[]
    for i in range(len(scores)):
        prob.append(scores[i]/sum_score)
    popular_new=[]#初始化一个新的已进化种群
    prob=np.array(prob)#将列表转化为数组，便于后面的概率选择
    for j in range(int(len(scores)/2)):#进行种群的一半次交配
        #随机按照概率选择一堆夫妻繁衍后代，概率越大，得分越大，越可能被选择
        for k in range(2):
            father=np.random.choice(popular_gene, p=prob.ravel())
            mother=np.random.choice(popular_gene, p=prob.ravel())
        #进行交配，产生下一代
        # 设置交叉率为0.66。
        is_mix = random.randint(0, 2)
        temp=[father,mother]
        x1=random.randint(1,10)
        x2=random.randint(11,17)
        if is_mix:#假设一对父母生两个，保持种群数量不变
            temp[0] = temp[0][0:x1] + temp[1][x1:x2] + temp[0][x2:]
            temp[0] = temp[1][0:x1] + temp[0][x1:x2] + temp[1][x2:]
        popular_new.append(temp[0])
        popular_new.append(temp[1])#将子女放入下一个进化时代
    return popular_new

def nature_variation(popular_gene):#对种群整体进行变异可能性操作，减小收敛于局部最优解的可能
    for i in range(len(popular_gene)):
        is_variation = random.uniform(0, 1)#随机生成0-1的任意数
        if is_variation < 0.05:#设置变异率
            rand = random.randint(2, 19)
            #假设变异，修改基因
            if popular_gene[i][rand] == '0':
                popular_gene[i] = popular_gene[i][0:rand] + '1' + popular_gene[i][rand+1:]
            else:
                popular_gene[i] = popular_gene[i][0:rand] + '0' + popular_gene[i][rand+1:]
    return popular_gene

if __name__=='__main__':
    num=100#定义种群的数目
    n=300#定义进化迭代次数
    popular = ori_popular(num)#得到一个随机的初始种群
    #1.对初始种群进行编码，模拟得到基因
    ori_popular_gene=encode(popular)#这里设置为18位基因
    new_popular_gene=ori_popular_gene#只是为了形式仪式感
    y=[]
    all_max_score=[]
    cnt=0
    fig = plt.figure()
    x=np.linspace(a, b, 100)
    for i in range(n):
        new_popular_gene=nature_choice(new_popular_gene)
        new_popular_gene=nature_variation(new_popular_gene)
        new_scores,xs=decode(new_popular_gene)#解码，得到新的种群的各个个体得分
        sum_score=np.sum(new_scores)#得到新种群的总得分，便于后面观察种群进化水平
        y.append(sum_score/num)#返回种群平均值
        max_score=np.max(new_scores)#得到该时代优秀的个体，即函数值最大的个体
        all_max_score.append(max_score)
        cnt+=1#技术种群的迭代次数
        if cnt%20==0:
            print(f"已经迭代进化{cnt}次，此时种群平均值为{sum_score/num}，最佳值为{max_score}")
        #画实观图
        plt.clf()
        plt.xlim((-10,10))
        plt.ylim((-40,40))
        plt.plot(x,func(x))
        plt.scatter(xs,new_scores, s=10, color='r')
        t=(n-1)*0.001
        plt.pause(t)#设置停留时间与迭代次数有关
    # 画进化图
    fig=plt.figure()
    nx = np.linspace(0, n, n)
    ax1 = fig.add_subplot(111)#坐标轴
    ax1.plot(nx, y,'r')#画出种群目前的情况
    ax1.plot(nx, all_max_score,'b')#画出种群目前的情况
    plt.legend(["平均值","最佳值"])
    plt.show()

运行结果如下：

已经迭代进化220次，此时种群平均值为32.05702212366586，最佳值为33.31557993182841
已经迭代进化240次，此时种群平均值为31.837646996627015，最佳值为33.27762194440547
已经迭代进化260次，此时种群平均值为31.45958626656901，最佳值为33.19787983857166
已经迭代进化280次，此时种群平均值为32.23062520430858，最佳值为33.19787983857166
已经迭代进化300次，此时种群平均值为32.087583602871064，最佳值为33.11304968140086

3.3核心思想即代码解释

和模拟退火算法一样，本质思想就是随机生成后产生新的解，并使得更好的解更多的留下来，最终迭代多次后求得最优解。
在遗传算法中我们可以将整个流程分为如下的
1.首先生成初始可行解，也就是初始种群（n个），比如在此题中即为区间内生成随机数，在旅行商问题中即为生成初始可行路径。
2.接下来将初始可行解进行编码，将其变为基因格式，常见的基因格式有二进制编码，onehot编码等等。
3.接下来的进行自然选择过程

#这里采用轮盘赌选择方法
def nature_choice(popular_gene):
    scores,xs=decode(popular_gene)#这里将函数值作为得分即可
    for i in range(len(scores)):
        if scores[i]<=0:
            scores[i]=0#剔除负数概率
    sum_score=np.sum(scores)#计算得到100个个体的总分，方面之后进行概率抽取
    #计算各个个体被选中的概率
    prob=[]
    for i in range(len(scores)):
        prob.append(scores[i]/sum_score)
    popular_new=[]#初始化一个新的已进化种群
    prob=np.array(prob)#将列表转化为数组，便于后面的概率选择
    for j in range(int(len(scores)/2)):#进行种群的一半次交配
        #随机按照概率选择一堆夫妻繁衍后代，概率越大，得分越大，越可能被选择
        for k in range(2):
            father=np.random.choice(popular_gene, p=prob.ravel())
            mother=np.random.choice(popular_gene, p=prob.ravel())
        #进行交配，产生下一代
        # 设置交叉率为0.66。
        is_mix = random.randint(0, 2)
        temp=[father,mother]
        x1=random.randint(1,10)
        x2=random.randint(11,17)
        if is_mix:#假设一对父母生两个，保持种群数量不变
            temp[0] = temp[0][0:x1] + temp[1][x1:x2] + temp[0][x2:]
            temp[0] = temp[1][0:x1] + temp[0][x1:x2] + temp[1][x2:]
        popular_new.append(temp[0])
        popular_new.append(temp[1])#将子女放入下一个进化时代
    return popular_new

具体来说，就是对种群（n个）进行解码并获取到种群得分，将得分对应权重的转为概率矩阵，pro = scores[i]/sum_score，也就是得分越高的越容易被选择。
4.接下来就是交配过程，我们可以假设每次交配产生两个孩子，保持种群总数不变，因此需要进行种群大小/2次交配，每次交配2只，每次配对的交配概率自由设置，可设置为逐渐减少的，也可以设置为固定概率的。随后依照得分概率在种群中随机选取两个个体，并将编码基因基因进行随机组合（也就是生物学上的基因的交叉结合），常用的手段有简单拼接等（即直接将父母之间的基因进行交叉组合，如0b11011和0b10121，直接将两者之间进行随机选取中间段拼接），还可以采用其他的如固定拼接，段落交叉等等，交叉的方式将决定算法的收敛情况。
5.得到新的孩子的基因后，依然需要判断，是否可行以及是否能够是否优于原来，并按照自己设置的方式选择是否保留，常用的手段有择优选取、有偏差的择优选取（在容忍程度下，尽管后代更差了，依然保留，因为可能距离最优解更近）。
6.在做完自然选择后，需要依照概率进行基因图片

def nature_variation(popular_gene):#对种群整体进行变异可能性操作，减小收敛于局部最优解的可能
    for i in range(len(popular_gene)):
        is_variation = random.uniform(0, 1)#随机生成0-1的任意数
        if is_variation < 0.05:#设置变异率
            rand = random.randint(2, 19)
            #假设变异，修改基因
            if popular_gene[i][rand] == '0':
                popular_gene[i] = popular_gene[i][0:rand] + '1' + popular_gene[i][rand+1:]
            else:
                popular_gene[i] = popular_gene[i][0:rand] + '0' + popular_gene[i][rand+1:]
    return popular_gene

也就是说，基因中的编码，有概率随机变化，如0突变为1，点（x，y）变为（x+1,y+1）等，其中的概率自由设置。
7.随后计算出整个种群的得分，保证种群稳步前进。并保存下最优的种群，并不断重复2-6过程，最终迭代获取最优解。

3.4遗传算法总结

相对模拟退火来说，遗传算法更加复杂一些，整体算法趣味也更大，其中有很多是我们可以自由设置的，如交叉结合方式，编码方式，突变概率和突变方式。
两个算法的思想非常类似，如果说非要做一点区分，那就是模拟退火更在意个人的优秀，而遗传算法则是站在整个种群的角度下，不断迭代使得整个种群向前发展。

4总结

启发式算法本质上思想都很类似，但是实现思想不太一致。
模拟退火的思想是，在不断搜索最优解的过程中不断降低波动范围（辐射范围），最终收敛于最优解附近。
遗传算法的思想在于不断搜索的过程中保证更高质量个体之间更有可能产生后代，最终整个种群将会趋于同质化，并且收敛于最优解。
而粒子群算法的思想在于，种群中的精英作为目标，整个种群朝着精英前进，最终精英之间交互产生，相互追逐，最终整个种群将会收敛于最优值附近。

参考资料

1.什么是NP-Hard
2.Python三维绘图–Matplotlib