python数学建模--模拟退火算法求解一元函数极值

算法模型

解空间bound、目标函数func、初始解s

基本思想

1. 设置参数：初始温度T、初始解$s_0$、降温系数$\delta$，
2. 定义目标函数func，生成初始函数值$f_0=func(s_0)$
3. 迭代过程开始
4. 产生新解$s_i$，计算新解对应的函数值$f_i=func(s_i)$
5. 计算增量$\Delta f=f_i-f_{i-1}$，若$\Delta f<0$表明新解对应的函数值比较小，此时接受新解和新的函数值，若$\Delta >0$则表明新解对应的函数值比较大，此时按照$Metropolis$准则接受新解
6. 检查是否满足终止条件，若满足则终止程序返回当前最优函数值，若不满足则循环步骤4~5

带约束条件的一元函数

函数表达式及图像

求解函数$y=sin(x^2-1)+2cos(2x),x\in [-8,8]$的极大值

# 本段代码用于绘制函数图像
import numpy as np
import matplotlib.pyplot as plt

bound = [-8, 8]
func = lambda x: np.sin(x**2-1)+2*np.cos(2*x)
plt.ylim([-4,4])
inputs=np.arange(bound[0],bound[1],0.1)
plt.title(r'$y=sin(x^2-1)+2cos(2x)$')
plt.xlabel('x')
plt.ylabel('f')
# 绘制两条正交直线
plt.axvline(x=[0.0],ls='--',color='g')
plt.axhline(y=[0.0],ls='--',color='g')
y=[func(x) for x in inputs]
plt.plot(inputs,y,'--')

xmax=inputs[y.index(max(y))]
# print(xmax,max(y))
plt.scatter(xmax,max(y),color='r')
plt.annotate(f'{round(xmax,2),round(max(y),2)}',xy=(xmax,max(y)),xytext=(round(xmax-1.5,2),round(max(y)+0.5,2)))
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

函数图像如下

退火算法实现

import numpy as np
import random
import matplotlib.pyplot as plt

# 设置参数：温度、最低温度、降温系数、解空间
T = 20000
Tmin = 1e-4
delta = 0.95
bound = [-8, 8]
# 定义目标函数,注意我们要求的是极大值，函数值取反即可
func = lambda x: -(np.sin(x ** 2 - 1) + 2 * np.cos(2 * x))
# 初始值定义（也表示最优值）
x = random.uniform(bound[0], bound[1])
f = func(x)

# 开启matplotlib交互模式，准备绘制动图
plt.ion()
# 绘制目标函数图像
inputs = np.arange(bound[0], bound[1], 0.1)
plt.plot(inputs, [func(x) for x in inputs], '--')
plt.ylim([-4, 4])

epoch = 0
while T > Tmin:
    is_change = False
    # 产生新解及新函数值
    x_new = x + (random.random() * 2 - 1) * T
    while not -8 < x_new < 8:# 保证生成的解在解空间中，使得每次温度降低是有意义的
        x_new = x + (random.random() * 2 - 1) * T
    f_new = func(x_new)
    # 计算增量
    if f_new < f:  # 以1的概率接受新解
        x, f = x_new, f_new
        is_change = True
    else:
        if pow(np.e, (f - f_new) / T) > random.random():  # 根据Metropolis准则接受新解
            x, f = x_new, f_new
            is_change = True
    plt.title(f'epoch:{epoch},ischange:{is_change}')
    plt.plot(x_new, f_new, 'go')

    plt.pause(0.5)
    T *= delta
    epoch += 1
plt.ioff()
plt.show()
print(f'函数的最大值为：{f},此时变量的值为：{x}')
>>> 函数的最值为：-2.9098918846391335,此时变量的值为：-3.0008540944272024
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

求解过程可视化

运行结果

加速（设置plt.pause()即可）

最终效果

上图中红色的点表明新解没有被Metropolis准则接受
遗留的问题：在程序运行的过程中笔者观察到下面的规律

1. 红色的点在开始的时候不会生成，也就是说温度较高的时候，产生的新解几乎全部被接受了，这应该是模拟退火算法的一个特点，但是笔者可能还不太明白
2. 在程序即将结束时，红色的点和绿色的点在极值点附近反复生成（如图中点稠密的位置），这个原因笔者暂时也没有弄明白

求解过程分析

首先观察程序的运行结果，与原图中我们得到的极值点几乎是完全一样的，也就是说我们花费了这么久研究算法得到的结果是有意义的！

得到了满意的结果之后，我们再来分析一下过程：实际上，笔者在研究退火算法时发现，它与爬山算法的代码十分相近，但不同的是，退火算法的Metropolis准则十分有效，他可以保证算法不拘泥于局部最优解，而是以一定的概率探索新的解

在上面的程序中，当新解对应的函数值f_new

然后我们再来观察一下求解过程的可视化图像，其中红点表示的是：在此次循环中产生的新解与旧解相比较次，而且没有被Metropolis准则接受，我们再来看下面两个图

# 本段程序用于生成图一
plt.ylim([-0.1,1.1])
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
plt.title(r'$f_{new} \geq f$情况下:$e^{\frac{f-f_{new}}{T}}$与P的比较')
plt.scatter(range(len(list1)),[float(i) for i in list1.keys()],color='g',label=r'$e^{\frac{f-f_{new}}{T}}$')
plt.scatter(range(len(list1)),list(list1.values()),color='b',label='P')
plt.xlabel('次数')
plt.ylabel('概率值大小')
plt.legend(loc='upper right')
plt.show()
1
2
3
4
5
6
7
8
9
10
11

运行结果

# 本段程序用于生成图二
plt.subplots_adjust(hspace=1)
plt.subplot(211)
plt.title('温度T的变化')
plt.plot(range(len(t_list)),t_list,color='g')
plt.subplot(212)
plt.title('是否更新')
plt.scatter(change.keys(),change.values(),color='b')
plt.show()
1
2
3
4
5
6
7
8
9

运行结果

我们可以从上面的图中发现在温度较高的时候，$e^{\frac{f-f_{new}}{T}}$的值会比较大。此时在Metropolis准则的作用下新解被接受，也就是说，即使现在产生的新解并不是很理想，算法还是决定接受，这样做有利于在解空间中寻找另外的优解而不拘泥于局部极值，实际上从上面的动态图中我们也可以看出来这一过程：在开始的时候即使下一个点对应的函数值大于前一个点对应的函数值，但是它仍然是绿色的，而且点的出现不是单区域爆发性的，而是整个区域内都有生成