模拟退火算法

模拟退火算法

文章主要参考：

• 模拟退火算法（Python）by 学习啊ZzZ
• 模拟退火算法（Simulated Annealing，SA） by 卡卡西~
  
• 【数之道17】金属冷却处理中隐藏的智慧-模拟退火优化算法 by FunInCode
  

概念阐述

​ 简述：模拟算法属于一种元启发式算法，所谓元启发式算法，是一类基于生物、物理、化学等自然现象和社会现象启发而构造的只能优化算法。主要结合了随机算法和局部搜索算法的启发优化式算法，它不是直接寻求问题的全局最优解，而是在允许的时间和空间内，提供一个接近最优的解。其设计通常是基于直观或经验。

​ 退火源自冶金学，是指将材料加热后再经过特定速率冷却。加热后粒子能量增大，在随机移动中找到使自己能量达到较小值的位置。退火冷却速度慢，使得原子可以有较多可能找到内能比原先更低的位置。

​ 模拟退火的原理也同金属退火原理近似，将搜寻空间的每一点当作晶体分子，在特定的温度下不断搜寻使自身能量达到较低值的位置。

结构简述

实际问题

以实际问题为例，求函数:
$$f(x)=x^3\ast cos(x)$$
在定义域内的最小值。
$$-3.1<=x<=6.16$$
对原函数进行求导后：
$$f^{\prime }(x)=2\ast x^2\ast cos(x)-3\ast x^3\ast sin(x)$$
令导函数为零并于原函数图象重合有：
定义域内极值点有：

易得区间内最小值的自变量为3.809，函数值为-43.41。

算法结构

模拟粒子在降温过程中随机搜索最优解的过程，我们的算法主结构将包含以下几个部分。

• 能量函数：主要是用于反应当前状态下的能量状态。
• 领域内新解生成函数：用于随机产生当前解附近的邻域解，以此模拟粒子状态变化。
• 新解接受准则：主要是使用Metropolis准则既判断能量的局部最低值，并且由当前温度决定粒子活跃程度，既接受较差新解的可能性。新解接受准则主要包含以下三条原则：
  • 首先，在固定的温度下，接受使得目标函数下降的候选解的概率要大于是目标函数上升的候选解的概率。
  • 二是随着温度的下降，接受使得目标函数上升的解的概率要逐渐减小。
  • 三是当温度区域零时，只能接受使目标函数下降的解。

这些函数将包括在模拟退火的主要过程中：

模拟退火的主要框架有

• ​ 由候选点的位置以及当前目标函数的值来模拟当前温度下粒子整体状态与能量。
  • ​ 由经由新解接受准则保留后的新的候选领域点来模拟当前温度下一定数量的粒子在随机跃迁后的结束状态。
    • 随机跃迁由领域内新解生成函数生成下一个领域点
    • 再经过能量函数计算当前目标函数值
    • 最后再根据新解接受准则决定是否保留当前候选点
  • 实现温度下降，来作为模拟退火算法的结束条件。

算法实现

初始状态：

def init(start, end):
    return random.uniform(start, end)

能量函数：

def energy(x):
    return pow(x, 3) * math.cos(x)

领域内新解生成函数：

def generate_x(start, end, x):
    while True:
        x += random.unifom(start, end)
        if start <= x <= end
        	return x

新解接受准则：

def metropolis(e, new_e, t):
    if new_e < e :
        return True
    else:
        return random.random() < math.exp((e - new_e) / t)

搜寻解集内最优解：

def search(all_x):
    best_X = all_x[0]
    best_e = energy(best_x)
    for x in all_x:
        e = energy(x)
        if e < best_e:
            best_x = x
    return best_x

主函数：

def SimulatedAnealing(Temperature, final_t, turns, alpha, a, b):
    all_x = []
    all_x.append(inite(a, b))
    curr_tem = Temperature
    while curr_tem > final_t:
    	x = search(all_x)
        e = energy(x)
        for int i in range(turns):
            new_x = generate(x)
            new_e = energy(new_x)
            if metropolis(e, new_e):
                x = new_x
                e = new e
        all_x.append(x)
        curr_tem = curr_tem * alpha
    return all_x

运行：

all_x = SimulatedAnealing(3000, 10, 200, 0.98, -3.1, 6.16)
#四舍五入保留三位小数
print(round(search(all_x), 3), round(energy(search(all_X)), 3))
#每次迭代产生的最优解
interation = len(all_x)
#最优解下的目标函数值,利用矩阵计算得出一维矩阵
all_energy = all_x ** 3 * np.cos(all_x)
#设置x轴名称
plt.xlabel('Interation')
plt.ylabel('F(x)')
#绘制数据
plt.plot(range(iteration), all_energy)
plt.show()

输出：

可以看到算法给出的最小值是比较精确的。

实际运用——TCP问题

算法实现

产生初始解：

#numbers为所有城市的数量
def init(numbers):
    # 随机选择给的范围内不重复的数构成np数组
    path_random = random.sample(range(numbers), k=numbers)
    # 将np数组转化为python列表
    return path_random

领域内新解生成：

这里主要选择的方法是随机选择两个点，将其内部路径翻转后获得一个新的随机路径。

def generate(pth):
    # 获得路径长度
    length = len(pth)
    # 随机产生首位两点,由于是作为列表下标的选择，所以可以从零开始.
    choice = random.choices(list(range(length)), k=2)
    # 获得子串
    lower = min(choice[0], choice[1])
    upper = max(choice[0], choice[1])
    mid_path = pth[lower:upper + 1]
    mid_path.reverse()
    new_path = pth[:lower]
    new_path.extend(mid_path)
    new_path.extend(pth[upper + 1:])
    return new_path

能量函数:

# length_mat为各个点之间的最短距离矩阵
def energy(pth, lenth_mat):
    length = len(pth)
    # 首先加上从终点返回值起始点的距离长度
    tol_len = length_mat[pth[0]][pth[length - 1]]
    # 迭代加上每个城市间的最短距离
    for city_st in range(length - 1):
        tol_len += length_mat[pth[city_st]][pth[city_st + 1]]
    return tol_len

metropolis准则：

def metropolis(e, new_e, t):
    if new_e < e:
        return True
    else:
        return random.random() < math.exp((e - new_e) / t)

搜索最优路径

def search(all_path, length_mat):
    best_e = energy(all_path[0], length_mat)
    best_pth = all_path[0]
    for pth in all_path:
        e = energy(pth, length_mat)
        if e < best_e:
            best_pth = pth
            best_e = e
    return best_pth

模拟退火主函数

def SimulatedAnealing(temperature, final_tem, interation, alpha, citys, length_mat):
    all_pth = []
    all_length = []
    # 初始化解
    all_pth.append(init(citys))
    all_length.append(energy(all_pth[0], length_mat))
    curr_tem = temperature
    while curr_tem > final_tem:
        pth = search(all_pth, length_mat)
        length = energy(pth, length_mat)
        for i in range(interation):
            new_pth = generate(pth)
            new_length = energy(new_pth, length_mat)
            if metropolis(length, new_length, curr_tem):
                pth = new_pth
                length = new_length
        all_pth.append(pth)
        all_length.append(length)
        curr_tem = curr_tem * alpha
    return all_pth, all_length

length_mat函数

#此函数用于创造最短路径矩阵，用于方便计算
#首先创建一个图对象
graph = nx.Graph()
#加入节点
grapg.add_nodes_from(range(1, 7))
edges = [
    (1, 2, 4), (1, 4, 2), (1, 5, 7),
    (2, 3, 3), (2, 4, 5), (2, 6, 8),
    (3, 4, 6), (3, 6, 7),
    (4, 5, 2), (4, 6, 3),
    (5, 6, 8), 
]
graph.add_weightd_edges_from(edges)
shortest_length = dict(nx.shortest_path_length(graph, weight = 'weight'))
length_mat = []
for i in range(1, 7):
    i_to_j = []
    for j in range(1, 7):
        i_to_j.append(shortest_length[i][j])
    length_mat.append(i_to_j)

运行：

all_path, all_length= SimulatedAnealing(3000, pow(10, -1), 200, 0.98, 6, length_mat)
print(search(all_path, length_mat), round(e(search(all_path, length_mat), length_mat)))
interation = len(all_path)
all_path = np.array(all_path)
all_length = np.array(all_length)
plt.xlabel("Iteration")
plt.ylabel("Length")
plt.plot(range(iteration), all_length)
plt.show

运行结果：