模拟退火算法介绍以及一个最简单模型的Python实现方式

        我们经常会看到以下四类问题,给定一个函数求极值、旅行商问题(TSP)、书店买书问题、背包问题。通常我们的解法是运用蒙特卡洛模拟or穷举法,但是当函数中自变量特别多时,这些方法的计算复杂度将非常非常大,显然不是我们在数模比赛中可以应用的。

        以上的两种方法运用的都是盲目搜索的理念,即每次获取新解的过程都是相互独立的,并没有运用在搜索过程中产生的中间信息;反之,如果利用了中间信息,我们称之为启发式搜索,例如机器学习中的梯度下降法等。我们看一下“爬山法”这个简单的启发式搜索算法的过程:

模拟退火算法介绍以及一个最简单模型的Python实现方式

        直到找到极大值爬山法就能停止搜索。但是有一个问题,爬山法很有可能会停留在局部最优解里,那么如何克服呢?假设我们求一个函数的最大值,爬山法告诉我们假如先确定了x(i),然后我们找到新点x(j),若后者函数值大于前者函数值,那么就取代前者,反之就不取代,为了避免就停留在局部最优点,我们引入概率P,表示在f(x(j))<f(x(i))的情况下我们依旧有概率接受x(j),但是如果前者大于后者我们还是选择直接取代。那么如何确定P呢?

        首先P一定介于[0,1]之间,可以设想假如P=1,就是我们的蒙特卡洛模拟or枚举法,也就是无论大于小于我们都接受,而P=0,就是我们的爬山法,一旦小于我们一定不接受。且,如果f(x(j))和f(x(i))越接近,我们就应该更乐于接受新解,P应该越大。为什么呢?直观理解,在越小的区间内,我们应该更多次的进行搜索,提高精准度,否则,越近P越小,那么一旦我们达到局部最优解,岂不是越难出去了?于是,我们有以下结论:

模拟退火算法介绍以及一个最简单模型的Python实现方式

        接受新解的P越大,我们搜索范围越广,搜索前期我们范围要广,后期要窄,那么P值后期就应该越小,所以P又和时间有关,且时间越长,P越小。我们引入以下式子:

模拟退火算法介绍以及一个最简单模型的Python实现方式        Ct就是C的一个有正向关系的函数,我们的搜索过程就(以最大值为例)应该如下所示:

模拟退火算法介绍以及一个最简单模型的Python实现方式

        但是,我们在实际建模中又会遇到其他问题,题目有约束条件怎么办?时间函数怎么定?搜索过程什么时候结束?怎么在A旁边随机生成B?因此,不同题目对应不同算法。而针对第二个问题,不同题目大体一致,这里就引入了模拟退火思想。

模拟退火算法介绍以及一个最简单模型的Python实现方式        那么就有:

模拟退火算法介绍以及一个最简单模型的Python实现方式

        怎么理解这两个结论呢?温度我们就对应时间,温度越高,时间就是越短的。那么,温度一定,也就是时间点一定,我们新旧值取值越接近,则概率越大;间距一定,如果温度越高,那么时间越短,时间越短,我们所要的是能够搜索的范围大,那么P值就应该大。

        在具体编程中,我们没法模拟时间的进行,但是我们可以模拟温度的慢慢变化。于是我们用高温代替时间刚开始(即我们刚开始循环),用低温表示时间进行了很久(即我们循环了很久),而每个温度下(时间点下)我们又进行多次搜索,然后降低温度,进行新的一轮,也就是2层嵌套的循环语句。

        那么什么时候停止?那有很多方法,温度足够低(时间够长),迭代次数够多(和前面一个意思),找到最优解M多次没有变化。

        那么最后一个问题,怎么从A更新到B?这个是具体题目不同的方法。我们可以自创、也可以查阅文献。下面以一个求给定函数最大值or最小值的更新方法:

模拟退火算法介绍以及一个最简单模型的Python实现方式

        那么,看一下这道题函数求最大值的代码吧(只展示主体部分,剩余部分可以私下与博主交流,包括动态可视化展示):

for iter in range(maxgen):                       #外循环, 采用的是指定最大迭代次数
    for i in range(1,(Lk+1)):                        #内循环,在每个温度下开始迭代
        y = np.random.randn(1,narvs)               #生成1行narvs列的N(0,1)随机数
        z = y / np.sqrt(np.sum(pow(y,2)))       #根据新解的产生规则计算z
        x_new = x0 + z*T[iter]                           #根据新解的产生规则计算x_new的值
        #如果这个新解的位置超出了定义域,就对其进行调整
        for j in range(narvs):
            if x_new[j] < x_lb[j]:
                r = np.random.rand(1)
                x_new[j] = r*x_lb[j]+(1-r)*x0[j];
            elif x_new[j] > x_ub[j]:
                r = np.random.rand(1)
                x_new[j] = r*x_ub[j]+(1-r)*x0[j]
        x1 = x_new                                 #将调整后的x_new赋值给新解x1
        y1 = Obj_fun1(x1)                          #计算新解的函数值
        if y1 > y0:                            #如果新解函数值大于当前解的函数值,更新当前解为新解
            x0 = x1
            y0 = y1
        else:
            p = np.exp(-(y0 - y1)/T[iter])               #根据Metropolis准则计算一个概率
            if np.random.rand(1) < p:               #生成一个随机数和这个概率比较,如果该随机数小于这个概率,更新当前解为新解
                x0 = x1 
                y0 = y1
        #判断是否要更新找到的最佳的解
        if y0 > max_y:                            #如果当前解更好,则对其进行更新
            max_y = y0                            #更新最大的y
            best_x = x0                           #更新找到的最好的x
    MAXY[iter] = max_y                           #保存本轮外循环结束后找到的最大的y
    T.append(alfa*T[iter])                                   #温度下降
    datax.append(x0)                  #更新散点图的x数据(此时解的位置在图上发生了变化)

        不足:当然模拟退火还有好多问题没有讨论,尽请关注下几期论文

本文主体思路借鉴某视频,侵权即删。         

上一篇:矩阵树定理的式子


下一篇:数据预处理ETL