初探 模拟退火算法、物理思想、代码实现（simulated annealing algorithm）

模拟退火算法概述

模拟退火的书籍和贴子非常多，这个贴子讲的就简单易懂
https://zhuanlan.zhihu.com/p/404204692。

算法流程简述

我简单总结一下模拟退火的大致流程，

• 做一个初始解，
• 外循环会控制温度参数T，
• 外循环的每一个step包含一个内循环，
• 内循环是遵循一些控制条件的局部搜索，
• 局部搜索会逐个得到的解，可能是更优的解，也可能是更差的解。
  • 模拟退火的精髓之一，以一定几率接受这个差的解。
  • 模拟退火的精髓之二，接受差解的几率和T相关。
• 经过上述的由外至内的多次循环，在提前结束条件和迭代次数等退出条件控制下，结束搜索，同时给出所有搜索过的相对最优解。

物理中的退火

我这里主要讲一下我对模拟退火这个名字的理解。为什么这个算法要被叫做退火？退火在局部搜索中，到底起到了什么作用？
（ps，纯属我个人的解读，欢迎大家批评指正。）

在物理中，固体物质中的分子在不断的振动，这种振动指的是，在分子附近的位置无规律的移动（特别地，在气体相中，就不仅仅是简单的振动，而是大范围的自由运行）。当温度比较高的时候，这种振动就越发强烈，当温度比较低的时候，这种振动就变得很微弱。

对比算法中的模拟退火

高温的情况。如果在一个邻域动作中，如果解从一个点移动到另一个点，可以看做是这种振动，当温度高的时候，很剧烈的振动也是合理的，也就是说，解的移动可以导致目标值很大的变化。也就是说，高温的情况，保证了搜索的广度。

低温的情况。和高温是相对的，在温度低的时候，目标值很剧烈的变化，被认为是几率很小的。因此，降低温度，就形成了搜索的广度到搜索的精细度的变化。

这种低温，大距离振动几率降低的情况，在物理上是通过玻尔兹曼分布形成的。这个玻尔兹曼分布，可以理解为，在随机振动的过程中，从一个目标值，随机跳到另一个目标值的几率。如果简单来说，比如说屋子里很热，真的是每个空气分子的运动速度都很快，而导致这些分子都高速撞到我们的脸上让我们热起来吗？并不是的，而是在高温的时候，分子能量高速度快的几率更大了。对于算法来说，也就是，在温度高的时候，大跨度、更剧烈的目标值跃迁的几率更大了。

所以模拟退火算法借用了玻尔兹曼分布的概率密度函数，
f = exp(E_j - E_i) / T，
当随机值大于这个值的时候，小于这个值的时候接受解的劣化。

举个例子，温度高的时候，解在local search中，可以有更多机会保留剧烈变化的差解。

对于算法来说，其中E是目标函数，也就是说，在算法流程中，外循环的T，逐渐降低。进而导致，内循环中，剧烈劣化的解被保留的几率越来越低。

简而言之，算法刚开始搜的时候，差解差很多，也可以保留，说不定后面能找到更好的解呢；而搜的时间久了，差很多的解就不那么想要了，变成了精细地在局部搜。而这些精确还是粗糙，都是以保留的概率的变化的形式实现的。

而恰好，上面这种算法先粗后细的思路，和退火过程中分子的振动由剧烈变微弱，很类似，还提供了现成的概率密度函数。所有就可以用模拟退火来表示这个算法。

代码实践

背包问题，容量8的背包，5个货物，大小为2, 3, 5, 1, 4，价值为2, 5, 8, 3, 6，追求总价值最大。算法流程图，参考https://zhuanlan.zhihu.com/p/404204692

import copy
import random

import numpy as np


class SimulatedAnnealingForKnapsackProblem:
    def __init__(self):
        self.knapsack_size = 8
        self.cargo_num = 5
        self.cargo_sizes = np.array([2, 3, 5, 1, 4])
        self.cargo_values = np.array([2, 5, 8, 3, 6])

        self.T = 10
        self.T_min = 1
        self.sum_value = 0
        self.sum_value_best = 0
        self.early_stop_no_better_step_count = 4
        self.seed = 1
        self.solution = np.zeros(self.cargo_num).astype(int)
        self.absolute_best_solution = np.zeros(self.cargo_num).astype(int)
        self.relative_best_solution = np.zeros(self.cargo_num).astype(int)

    def solve_absolute_best_solution(self):
        print("solve_absolute_best_solution() start...")
        pow_2_n_list_small = [np.power(2, n) for n in range(self.cargo_num)]
        pow_2_n_list_small.reverse()
        pow_2_n_list_big = [np.power(2, self.cargo_num)] + pow_2_n_list_small[:-1]

        sum_value_max = 0
        for o in range(0, np.power(2, self.cargo_num)):
            self.solution = np.array([o % pow_2_n_list_big[i] // pow_2_n_list_small[i] for i in range(self.cargo_num)])
            sum_size = np.dot(self.solution, self.cargo_sizes)
            if sum_size > self.knapsack_size:
                continue
            sum_value = np.dot(self.solution, self.cargo_values)
            if sum_value >= sum_value_max:
                sum_value_max = sum_value
                self.absolute_best_solution = copy.deepcopy(self.solution)
        self.sum_value_best = sum_value_max
        print(f"\tabsolute_best_solution: {self.absolute_best_solution}, sum_value_max: {self.sum_value_best}")

    def solve_simulated_annealing(self):
        print("solve_simulated_annealing() start...")
        random.seed(self.seed)
        self._initialize_solution()
        old_best_sum_value = self.sum_value_best
        no_better_step_count = 0
        while self.T >= self.T_min:
            self._local_search()
            if self.sum_value_best > old_best_sum_value:
                old_best_sum_value = self.sum_value_best
                no_better_step_count = 0
            else:
                no_better_step_count += 1
            if no_better_step_count >= self.early_stop_no_better_step_count:
                break
            self.T -= 1
        print(f"\trelative_best_solution: {self.relative_best_solution}, sum_value_max: {self.sum_value_best}")

    def _initialize_solution(self):
        self.solution = np.zeros(self.cargo_num).astype(int)
        self.sum_value = np.dot(self.solution, self.cargo_values)
        self.sum_value_best = self.sum_value
        self.T = 10

    def _local_search(self):
        n_step = 3
        i_step = 0
        while i_step < n_step:
            tmp_solution = copy.deepcopy(self.solution)
            rand_neigh = random.randint(0, 1)
            neigh_str = "neigh flip" if rand_neigh == 0 else "neigh swap"
            if rand_neigh == 0:
                rand_idx = random.randint(0, self.cargo_num - 1)
                tmp_solution[rand_idx] = int(not tmp_solution[rand_idx])
            else:
                rand_idx1 = random.randint(0, self.cargo_num - 1)
                rand_idx2 = random.randint(0, self.cargo_num - 2)
                if rand_idx2 == rand_idx1:
                    rand_idx2 = self.cargo_num - 1
                if tmp_solution[rand_idx1] == tmp_solution[rand_idx2]:
                    continue
                tmp_solution[rand_idx1], tmp_solution[rand_idx2] = tmp_solution[rand_idx2], tmp_solution[rand_idx1]
            sum_weight = np.dot(tmp_solution, self.cargo_sizes)
            if sum_weight > self.knapsack_size:
                continue
            sum_value = np.dot(tmp_solution, self.cargo_values)
            if sum_value < self.sum_value:
                p = np.exp((sum_value - self.sum_value) / self.T)
                if random.random() > p:
                    continue
            self.solution = copy.deepcopy(tmp_solution)
            self.sum_value = sum_value
            if self.sum_value > self.sum_value_best:
                self.relative_best_solution = copy.deepcopy(self.solution)
                self.sum_value_best = self.sum_value
                print(f"\tT: {self.T}, i_step: {i_step}, neigh: {neigh_str}, "
                      f"solution: {self.relative_best_solution}, sum_value: {self.sum_value_best}, better")
            else:
                print(f"\tT: {self.T}, i_step: {i_step}, neigh: {neigh_str}, "
                      f"solution: {self.solution}, sum_value: {self.sum_value}")
            i_step += 1


if __name__ == '__main__':
    a = SimulatedAnnealingForKnapsackProblem()
    a.solve_absolute_best_solution()
    a.solve_simulated_annealing()

运行结果：

solve_absolute_best_solution() start...
	absolute_best_solution: [0 1 0 1 1], sum_value_max: 14
solve_simulated_annealing() start...
	T: 10, i_step: 0, neigh: neigh flip, solution: [0 0 0 0 1], sum_value: 6, better
	T: 10, i_step: 1, neigh: neigh flip, solution: [0 0 0 1 1], sum_value: 9, better
	T: 10, i_step: 2, neigh: neigh flip, solution: [1 0 0 1 1], sum_value: 11, better
	T: 9, i_step: 0, neigh: neigh flip, solution: [0 0 0 1 1], sum_value: 9
	T: 9, i_step: 1, neigh: neigh flip, solution: [0 0 0 0 1], sum_value: 6
	T: 9, i_step: 2, neigh: neigh swap, solution: [0 1 0 0 0], sum_value: 5
	T: 8, i_step: 0, neigh: neigh swap, solution: [0 0 0 0 1], sum_value: 6
	T: 8, i_step: 1, neigh: neigh swap, solution: [1 0 0 0 0], sum_value: 2
	T: 8, i_step: 2, neigh: neigh flip, solution: [1 0 0 1 0], sum_value: 5
	T: 7, i_step: 0, neigh: neigh flip, solution: [1 1 0 1 0], sum_value: 10
	T: 7, i_step: 1, neigh: neigh flip, solution: [1 1 0 0 0], sum_value: 7
	T: 7, i_step: 2, neigh: neigh swap, solution: [1 0 0 0 1], sum_value: 8
	T: 6, i_step: 0, neigh: neigh flip, solution: [1 0 0 1 1], sum_value: 11
	T: 6, i_step: 1, neigh: neigh swap, solution: [1 0 1 1 0], sum_value: 13, better
	T: 6, i_step: 2, neigh: neigh flip, solution: [0 0 1 1 0], sum_value: 11
	T: 5, i_step: 0, neigh: neigh swap, solution: [0 1 0 1 0], sum_value: 8
	T: 5, i_step: 1, neigh: neigh swap, solution: [1 0 0 1 0], sum_value: 5
	T: 5, i_step: 2, neigh: neigh flip, solution: [1 0 1 1 0], sum_value: 13
	T: 4, i_step: 0, neigh: neigh swap, solution: [1 1 0 1 0], sum_value: 10
	T: 4, i_step: 1, neigh: neigh flip, solution: [1 1 0 0 0], sum_value: 7
	T: 4, i_step: 2, neigh: neigh flip, solution: [1 1 0 1 0], sum_value: 10
	T: 3, i_step: 0, neigh: neigh swap, solution: [1 0 0 1 1], sum_value: 11
	T: 3, i_step: 1, neigh: neigh swap, solution: [0 1 0 1 1], sum_value: 14, better
	T: 3, i_step: 2, neigh: neigh flip, solution: [0 0 0 1 1], sum_value: 9
	T: 2, i_step: 0, neigh: neigh swap, solution: [0 0 1 1 0], sum_value: 11
	T: 2, i_step: 1, neigh: neigh swap, solution: [0 0 0 1 1], sum_value: 9
	T: 2, i_step: 2, neigh: neigh swap, solution: [1 0 0 0 1], sum_value: 8
	T: 1, i_step: 0, neigh: neigh swap, solution: [0 0 0 1 1], sum_value: 9
	T: 1, i_step: 1, neigh: neigh flip, solution: [0 1 0 1 1], sum_value: 14
	T: 1, i_step: 2, neigh: neigh flip, solution: [0 1 0 0 1], sum_value: 11
	relative_best_solution: [0 1 0 1 1], sum_value_max: 14