爬山算法的缺点

爬山算法是一种简单的贪心搜索算法，它的原理是，只向高出走，如果发现移动后的高度变小了，那么就不移动，如果高度变高了就移动。该算法每次从当前解的临近解空间中选择一个最优解作为当前解，直到达到一个局部最优解。

爬山算法实现很简单，其主要缺点是会陷入局部最优解，而不一定能搜索到全局最优解。假设B点为当前解，爬山算法搜索到A点这个局部最优解就可能会停止搜索，因为在A点无论向那个方向小幅度移动都不能得到更优的解。

模拟退火算法实现

模拟退火算法从某一较高初温出发，伴随温度参数的不断下降，结合 一定的概率 突跳特性 在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。

核心思想：服从目标就一定接受目标，如果不服从目标就按概率接受。关键点就是这个概率，模拟每一步走的步长也是关键点。

爬山算法和模拟退火的比喻

  • 爬山算法：兔子朝着比现在高的地方跳去。它找到了不远处的最高山峰。但是这座山不一定是珠穆朗玛峰。这就是爬山算法，它不能保证局部最优值就是全局最优值。
  • 模拟退火：兔子喝醉了（它本来知道自己要去哪里）。它随机地跳了很长时间。这期间，它可能走向高处，也可能踏入平地。但是，它渐渐清醒了并朝最高方向跳去。这就是模拟退火。

模拟退火算法的数学原理

核心算法：

背后的原理和物理上的一模一样。能量越低越有秩序，越靠近最优解。

  • 温度越高接受错误解的概率越高
  • 温度越低接受错误解的概率越低

错误解：目标是求最小值，正确解是解变小的方向，错误解就是变大的方向。在爬山算法中接受错误解的概率就是0。

问题案例

问题

求

Python代码

研究函数

1
2
3
def func(self, x, y):
    res = 4*x**2 - 2.1*x**4 + x**6/3 + 4*y**2 + 4*y**4
    return res

初始化

1
2
3
4
5
6
7
8
9
10
11
12
class SA:  # 接下来的方法都是这个对象的
    def __init__(self, func, N=100, T0=100, Tf=0.01, alpha=0.99):
        self.func = func
        self.N = N        		# 每一代点的个数
        self.alpha = alpha		# 降温系数
        self.T0 = T0			# 初始温度
        self.Tf = Tf			# 温度终值
        self.T = T0				# 当前温度
        self.x = [random()*11-5 for i in range(N)]   # 随机生成iter个x的值
        self.y = [random()*11-5 for i in range(N)]	# 随机生成iter个y的值
        self.most_best = []
        self.history = {'f': [], 'T': []}

生成点，由于 $x$ 的范围已经规定了，所以随机生成的点就是random()*11-5

参数解释：

  • self不解释
  • func：待研究函数
  • iter：迭代次数
  • T0：初始温度 一般是100左右，如果解的范围比较大，可以设置大一点
  • Tf：温度终值 一般接近于0，但不等于0（看下一个函数，就知道等于0时就没有意义了，不扰动了）
  • alpha：降温系数 一般取0.85～0.99

温度变化公式：$T = \alpha T$

画个温度曲线图演示一下，0.99不太明显，用0.9演示：

温度下降曲线

随机扰动产生新解

1
2
3
4
5
6
7
def generate_new(self, x, y):
    while True:
        x_new = x + self.T*(random() - random())
        y_new = y + self.T*(random() - random())
        if (-5 <= x_new <= 5) & (-5 <= y_new <= 5):
            break
    return x_new, y_new

random() - random()理解为就是一个随机数就可以了。

所以 $\Delta x$ 与当前温度和随机值有关。所以在温度比较高的时候 $\Delta x$ 比较大，寻找新解的步长比较大。也就是突跳特性。

  • 前期找到更多新区域
  • 后期找到局部最优解
  • 局部最优解的出全局最优解

Metroplis准则

1
2
3
4
5
6
7
8
9
def Metropolis(self, f, f_new):
    if f_new <= f:
        return 1
    else:
        p = math.exp((f - f_new) / self.T)
        if random() < p:
            return 1
        else:
            return 0

这个函数就是判断是否选取新解。

上面提到了如果是正确解就一定选择，如果是错误解就按概率选择。

2-3行体现了正确解一定选择。

第5行是计算当前温度下的选择概率。

第6行生成一个随机数，来模拟是否选择。

代码很简单。重点是概率的计算。这里是：

更规范的描述这一步的公式：

当然，要根据题目分析正确解的情况。

获取最优目标函数值

1
2
3
4
5
6
7
8
def best(self):
    f_list = []  # 保存每个点通过函数计算后的值
    for i in range(self.N): # 每一代有 N 个点
        f = self.func(self.x[i], self.y[i])
        f_list.append(f)
    f_best = min(f_list) # 求这一代中的最优解
    idx = f_list.index(f_best) # 得到最优解的索引
    return f_best, idx

主体函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
def run(self):
    while self.T > self.Tf: # 外循环，降温过程，判断温度是否降到期望温度
        for i in range(self.N): # 内循环遍历每一个解
            f = self.func(self.x[i], self.y[i]) # 计算初始解
            x_new, y_new = self.generate_new(self.x[i], self.y[i]) # 随机扰动，产生新解
            f_new = self.func(self.x[i], self.y[i]) # 计算扰动后的新解
            if self.Metropolis(f, f_new):
                self.x[i] = x_new
                self.y[i] = y_new
        ft, idx = self.best()
        self.history['f'].append(ft)
        self.history['T'].append(self.T)
        
        # 降温，别写成死循环了
        self.T = self.T*self.alpha
    # 输出最优解
    f_best, idx = self.best()
    print(f"F={f_best}, x={self.x[idx]}, y={self.y[idx]}")

整体代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
import math
from random import random
import matplotlib.pyplot as plt


def func(x, y):
    res = 4 * x ** 2 - 2.1 * x ** 4 + x ** 6 / 3 + 4 * y ** 2 + 4 * y ** 4
    return res


class SA:  # 接下来的方法都是这个对象的
    def __init__(self, func, N=100, T0=100, Tf=0.01, alpha=0.99):
        self.func = func
        self.N = N  # 每一代点的个数
        self.alpha = alpha  # 降温系数
        self.T0 = T0  # 初始温度
        self.Tf = Tf  # 温度终值
        self.T = T0  # 当前温度
        self.x = [random() * 11 - 5 for i in range(N)]  # 随机生成iter个x的值
        self.y = [random() * 11 - 5 for i in range(N)]  # 随机生成iter个y的值
        self.history = {'f': [], 'T': []}

    def generate_new(self, x, y):
        while True:
            x_new = x + self.T * (random() - random())
            y_new = y + self.T * (random() - random())
            if (-5 <= x_new <= 5) & (-5 <= y_new <= 5):
                break
        return x_new, y_new

    def Metropolis(self, f, f_new):
        if f_new <= f:
            return 1
        else:
            p = math.exp((f - f_new) / self.T)
            if random() < p:
                return 1
            else:
                return 0

    def best(self):
        f_list = []  # 保存每个点通过函数计算后的值
        for i in range(self.N):  # 每一代有 N 个点
            f = self.func(self.x[i], self.y[i])
            f_list.append(f)
        f_best = min(f_list)  # 求这一代中的最优解
        idx = f_list.index(f_best)  # 得到最优解的索引
        return f_best, idx

    def run(self):
        while self.T > self.Tf:  # 外循环，降温过程，判断温度是否降到期望温度
            for i in range(self.N):  # 内循环遍历每一个解
                f = self.func(self.x[i], self.y[i])  # 计算初始解
                x_new, y_new = self.generate_new(self.x[i], self.y[i])  # 随机扰动，产生新解
                f_new = self.func(self.x[i], self.y[i])  # 计算扰动后的新解
                if self.Metropolis(f, f_new):
                    self.x[i] = x_new
                    self.y[i] = y_new
            ft, idx = self.best()
            self.history['f'].append(ft)
            self.history['T'].append(self.T)

            # 降温，别写成死循环了
            self.T = self.T * self.alpha
        # 输出最优解
        f_best, idx = self.best()
        print(f"F={f_best}, x={self.x[idx]}, y={self.y[idx]}")


sa = SA(func)
sa.run()

# 画图
plt.plot(sa.history['T'], sa.history['f'])
plt.title('SA')
plt.xlabel('T')
plt.ylabel('f')
plt.gca().invert_xaxis()
plt.show()

最终结果：

F=1.4410512133214726, x=-0.2964965160195926, y=0.4748706225379135

image-20230214115435954

Matlab 代码

这里的Matlab代码实际上是根据上面的python代码改编的，写得不太严谨。

Matlab中有模拟退火算法的工具箱，所以也不用自己写模拟退火算法。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
clc, clear;

% 研究函数
fun = @(x, y) 4*x^2 - 2.1*x^4 + x^6/3 + 4*y^2 + 4*y^4;

% 相关参数
N = 100;
T0 = 100;
Tf = 0.01;
alpha = 0.99;

% 产生初始解
x = rand(2, 100)*10-5;
T = T0;

% 存储最优解
history = [];

while Tf < T
    for i = 1 : N
        f = fun(x(1, i), x(2, i));
        x_new = generate_new(x, i, T);
        f_new = fun(x_new(1, i), x_new(2, i));
        flag = Metropolis(f, f_new, T);
        if flag
            x(1, i) = x_new(1, i);
            x(2, i) = x_new(2, i);
        end
    end
    [ft, idx] = best(fun, x, N);
    history = [history, ft];
    T = T*alpha;
end
hs = size(history);
plot(1:hs(2), history);

function x_new = generate_new(x, i, T)
    x_new = x;
    while true
        x_new(1, i) = x(1, i) + T * (rand - rand);
        x_new(2, i) = x(2, i) + T * (rand - rand);
        if -5 <= x_new(1, i) <= 5 && -5 <= x_new(2, i) <= 5
            break
        end
    end
end

function flag = Metropolis(f, f_new, T)
    if f >= f_new
        flag = 1;
    else
        p = exp((f - f_new) / T);
            if rand() < p
                flag = 1;
            else
                flag = 0;
            end
    end
end

function [f_best, idx] = best(fun, x, N)
    f_list = [];
    for i = 1 : N
        fb = fun(x(1, i), x(2, i));
        f_list = [f_list, fb];
    end
    [f_best, idx] = min(f_list);
end

应用场景

优化类问题

  • 图论中的最短路径问题
  • 线性、非线性规划问题