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

本文主要是介绍模拟退火算法介绍以及一个最简单模型的Python实现方式，对大家解决编程问题具有一定的参考价值，需要的程序猿们随着小编来一起学习吧！

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

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

        直到找到极大值爬山法就能停止搜索。但是有一个问题，爬山法很有可能会停留在局部最优解里，那么如何克服呢？假设我们求一个函数的最大值，爬山法告诉我们假如先确定了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越小，那么一旦我们达到局部最优解，岂不是越难出去了？于是，我们有以下结论：

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

        Ct就是C的一个有正向关系的函数，我们的搜索过程就（以最大值为例）应该如下所示：

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

        那么就有：

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

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

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

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

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

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数据（此时解的位置在图上发生了变化）

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

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

这篇关于模拟退火算法介绍以及一个最简单模型的Python实现方式的文章就介绍到这儿，希望我们推荐的文章对大家有所帮助，也希望大家多多支持为之网！