模拟退火算法的通俗理解及代码实现

本篇文章部分内容引自百科

问题来源:

旅行商问题，即TSP问题（Traveling Salesman Problem）又译为旅行推销员问题、货郎担问题，是数学领域中著名问题之一。假设有一个旅行商人要拜访n个城市，他必须选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。路径的选择目标是要求得的路径路程为所有路径之中的最小值。

作为图论问题

可以用无向加权图来对TSP建模，则城市是图的顶点，道路是图的边，道路的距离就是该边的长度。它是起点和终点都在一个特定顶点，访问每个顶点恰好一次的最小化问题。通常，该模型是一个完全图（即每对顶点由一条边连接）。如果两个城市之间不存在路径，则增加一条非常长的边就可以完成图，而不影响计算最优回路。

非对称和对称

在对称TSP问题中，两座城市之间来回的距离是相等的，形成一个无向图。这种对称性将解的数量减少了一半。在非对称TSP问题中，可能不是双向的路径都存在，或是来回的距离不同，形成了有向图。交通事故、单行道和出发与到达某些城市机票价格不同等都是打破这种对称性的例子。
相关问题.

图论中的一个等价形式是：给定一个加权完全图（顶点表示城市，边表示道路，权重就会是道路的成本或距离）, 求一权值最小的哈密尔顿回路。

返回到起始城市的要求不会改变问题的计算复杂度，见哈密顿路径问题。

另一个相关问题是瓶颈旅行商问题（bottleneck TSP）：求加权图中权重最大的边最小的哈密尔顿回路。问题在运输和物流之外都有相当广泛的实际意义。一个典型的例子是在印刷电路板制造中：规划打孔机在PCB版上钻孔的路线。在机械加工或钻孔应用中，“城市”是需要加工的部分或需要钻的（不同大小）的孔，而“旅行成本”包括更换机具所用的时间（单机作业排序问题）。

广义旅行商问题，又称“旅行政客问题”，处理“国家”中有（一个或多个）“城市”，而旅行商需要在每个“国家”访问恰好一座“城市”。其中一种应用是在求解下料问题时，想要最小化刀具改变次数中。另一种应用与半导体制造业中的打孔有关。令人惊喜的是，Behzad与Modarres证明了广义旅行商问题可以转化为相同城市数量的标准旅行商问题 ，只是要改变距离矩阵。

优先顺序旅行推销员问题处理城市之间存在访问次序的问题。

旅行购买者问题涉及购买一系列产品的购买者。他可以在若干城市购买这些产品，但价格会有不同，也不是所有城市都有售相同的商品。目标是在城市的子集中间找到一条路径，使得总成本（旅行成本 + 购买成本）最小，并且能够买到所有需求的商品。

要解决的问题:
求出某个目标函数的最值（某一给定的函数，旅行的路程或费用，买书的花费，利润）
其中，最大值问题可以通过给目标函数添加负号转换为求最小值问题

求解的通用步骤
1，蒙特卡洛模拟(这里用其求解最值问题)
随机生成多组解，验证这些解是否满足题目的约束条件
若满足将其保存到一个可行集中，然后计算这个可行集内每个解对应的目标函数值，在里面找到最值即可。
对比穷举法:按照某一特定规则来穷举整个空间，而蒙特卡洛模拟是随机搜索，这是两者的区别
共同缺点:如果求解的空间非常非常大的时候，这里两种方法所耗费的计算资源将难以估计

即np-hard问题
其中，NP是指非确定性多项式（non-deterministic polynomial，缩写NP）。所谓的非确定性是指，可用一定数量的运算去解决多项式时间内可解决的问题。
那么后面加个hard，意思就是用代数问题难以求解

盲目搜索
盲目搜索方法又叫非启发式搜索，是一种无信息搜索，一般只适用于求解比较简单的问题，盲目搜索通常是按预定的搜索策略进行搜索，而不会考虑到问题本身的特性。常用的盲目搜索有宽度优先搜索和深度优先搜索两种。

这样做就与我们上面提到的计算出全部求解方案一样，于是我们就引出了新概念，启发式搜索

启发式搜索
启发式搜索(Heuristically Search)又称为有信息搜索(Informed Search)，它是利用问题拥有的启发信息来引导搜索，达到减少搜索范围、降低问题复杂度的目的，这种利用启发信息的搜索过程称为启发式搜索。即加速找到最优解方法与寻找到最优解的方法都是启发式搜索

模拟退火算法就是启发式算法的一种
类似的还有蚁群(模拟蚂蚁觅食)算法，粒子群(模拟鸟类飞行)算法，遗传(模拟基因遗传)算法

我们先引入一个概念，爬山法

爬山法是指经过评价当前的问题状态后，限于条件去增加这一状态与目标状态的差异，经过迂回前进，最终达到解决问题的总目标。就如同爬山一样，为了到达山顶，有时不得不先上矮山顶，然后再下来，这样翻越一个个的小山头，直到最终达到山顶。可以说，爬山法是一种"以退为进"的方法，往往具有"退一步进两步"的作用，后退乃是为了更有效地前进。爬山法也叫逐个修改法、瞎子摸象法。

解多变量无约束最优化问题的一类方法。有的书上称直接法或直接搜索法，是通过点的直接移动产生的目标值有所改善的点，经过这样的移动，逐步到达使目标函数最优的点。如果我们把目标函数的几何图形看成一个山峰，那么点的直接移动就像人在爬山，选择方向，逐步向山顶移动。可分为轴向搜索法、单纯形调优法、Powell法等。轴向搜索法是以沿坐标轴方向移动为基础的搜索方法，在进行每一轮沿坐标轴方向搜索时，是从一参考点出发，依次沿平行于各个坐标轴方向连续作对应的目标函数值改进的搜索移动，并以最后获得的点作为下一轮迭代点。同时，为提高求解的效率．还要采取某些加快收敛的措施。 [1]
问题求解的过程就是努力沟通问题的起始状态和目标状态之间的联系链条，由起始状态出发，逐步向目标推移、逼近的过程。在思维课题的求解活动中，人们几乎总是一直关注着所要达到的最终目标，试图不断地向目标逼近。这就像在登山活动中运动员时刻把顶峰放在心目之中，力图接近它、占领它一样。在解决思维课题时，我们经常自觉或不自觉地运用着能够尽量向目标靠拢的方法。这也就是通常所说的爬山法。 [2]
简单地说，爬山法就是按照下述原则进行试探的方法。这就是在每一状态下都选取这样的步骤进行试探：由这个步骤所达到的终结状态是所有可容许步骤所能达到的终结状态中，最接近最终目标的一个。这样的步骤称为最优步骤。按照这样的原则依次选取步骤，顺序试探下去，直到最终目标为止。如果达不到最终目标，可以回过头来，从已经经历过的某一中间状态开始，改用直接效果稍差一点的次优步骤，沿着另一条分支途径再行试探下去。当然，也可以一下子回转到整个问题的起始状态，沿着另一条全新的途径进行试探。到底应当返回到什么状态开始另行试探，这要由解题者根据具体情况作出自己的判断。

步骤:
1在解空间中随机生成一个初始解
2根据初始解的位置，我们在下一步有两种走法，向左邻域一步或向右邻域一步，步长越小越好但会增大计算量
3比较不同走法下目标函数的大小，并决定下一步往哪个放向走，比如向右走目标函数较大，因此我们往右走一小步
4不断重复这个步骤，直到找到一个极大值点(或定义域边缘处)，然后结束搜索

但是这样做很容易陷入局部最优解

于是我们引入了向左走或者向右走的概率，就有了模拟退火

简单地理解就是
物理的退火降温的过程，
一个处于很高温度的物体，要给它降温，使物体的内能降到最低。
一般来讲，越快越好，让它的温度迅速地降低。
然而，过快地降温使得物体来不及有序地收缩，难以形成结晶。结晶态才是物体的内能最低的形态。
正确的做法，是慢慢降温，也就是退火，才能使得物体的每一个粒子都有足够的时间找到自己的最佳位置并紧密有序地排列。开始温度高的时候，粒子活跃地运动并逐渐找到一个合适的状态。在这过程中温度也会越降越低，温度低下来以后，粒子也渐渐稳定下来，相较于以前不那么活跃了。这时候就可以慢慢形成最终稳定的结晶态了。

简单来说就是新的结果与旧的结果越接近，就越愿意接受新解

那么，要如何计算这个解的概率呢

P 正比于 1/|(旧解-新解)|

那么，我们现在需要一个函数将P与旧解和新解的关系连接起来，同时该函数要在0~1之间，
于是我们想到了

故我们假设

这样同时也解决了1/|(旧解-新解)|中分母为0的问题

接受新解的概率P越大就表明算法的搜索更容易接受不好的解，这时候就相当于增加了算法的搜索范围

现在我们来理解下模拟退火的参数与过程

1 T,温度(注意，退火算法里面最重要的就是温度)

2 ΔT,温度变化率，每次温度等于上一次温度 * ΔT，实际应用时一般取0.95−0.99，模拟慢慢降温

3 x,当前选择的解

4 Δx,旧解与新解的差值

5 x1,当前的目标解，等于x+Δx

6 Δf,当前解的函数值与目标解函数值之差，等于f(x)−f(x1)

给一个初始解x，让它不断变动。模拟变动的大小随温度的降低而降低，每次的Δx应该在一个大小与T成正比的值域内随机取值。

这时候就要考虑是否将当前解变为目标解。因为我们还是需要贪心，所以如果f(x1)<f(x)，那么接受目标解，x=x1

那如果f(x1)>f(x)呢？当然要以一定概率接受它啦，这样才能跳出局部最优解，去搜寻更优的解，弥补贪心的局限性。那么这个概率应该是多少呢？同样要模拟变动的大小随温度的降低而降低。这个概率就是

将概率引入时间变量

这里取倒数是为了是Ct关于t递增(因为0.95是小数)

如此反复选择直到T趋近于0（可以设一个EPS）这时候我们就认为当前的x就是最优解

#本功能实现最小值的求解#

from matplotlib import pyplot as plt
import numpy as np
import random
import math
plt.ion()#把matplotlib改为交互状态

#初始值设定
hi=3
lo=-3
alf=0.95
T=100

#目标函数
def f(x):
    return 11*np.sin(x)+7*np.cos(5*x)##注意这里要是np.sin

#可视化
def visual(x):
    plt.cla()
    plt.axis([lo-1,hi+1,-20,20])
    m=np.arange(lo,hi,0.0001)
    plt.plot(m,f(m))
    plt.plot(x,f(x),marker='o',color='black',markersize='4')
    plt.title('temperature={}'.format(T))
    plt.pause(0.1)#如果不停啥都看不见

#随机产生初始值
def init():
    return random.uniform(lo,hi)

#新解的随机产生
def new(x):
    x1=x+T*random.uniform(-1,1)
    if (x1<=hi)&(x1>=lo):
        return x1
    elif x1<lo:
        rand=random.uniform(-1,1)
        return rand*lo+(1-rand)*x
    else:
        rand=random.uniform(-1,1)
        return rand*hi+(1-rand)*x

#p函数
def p(x,x1):
    return math.exp(-abs(f(x)-f(x1))/T)

def main():
    global x
    global T
    x=init()
    while T>0.0001:
        visual(x)
        for i in range(500):
            x1=new(x)
            if f(x1)<=f(x):
                x=x1
            else:
                if random.random()<=p(x,x1):
                    x=x1
                else:
                    continue
        T=T*alf
    print('最小值为：{}'.format(f(x)))

main()