题目如下:
通过模拟退火算法用python求解这个问题,有以下几个步骤:
- 搭建模拟退火算法的框架
- 将数学公式转写为python函数
- 设计自变量的获取方式和更新方式
- 带入算法进行计算
搭建模拟退火算法的框架
import random
import numpy as np
import math
t = 2000 #t是当前温度
T = 2000 #T是初始温度
dt = 0.993 #dt是温度的变化率
eps = 1e-14
n = 0 #n是外循环的运行次数
x = [0,0,0] #初始化解,用列表存储
while n < 200 : #这里的数字可以设置为希望外循环运行的次数
#中间是退火代码的实现方式,先留空
t = t * dt
n += 1
将数学公式转写为python函数
定义目标函数:
def target_fun(x:list):
return x[0]**2 + x[1] ** 2 + x[2]**2 + 8
定义约束条件:
def st_fun1(x:list):
if x[0]**2 - x[1] + x[2]**2 >= 0:
return True
else:
return False
这里可以将函数定义成lambda表达式,可以让我们的代码代码更简洁,这样写的效果和上面写的效果和用法是一样的:
target_fun = lambda x: x[0]**2 + x[1] ** 2 + x[2]**2 + 8
st_fun1 = lambda x: x[0]**2 - x[1] + x[2]**2 >= 0
st_fun2 = lambda x: x[0] + x[1]**2 +x[2]**2 <= 20
st_fun3 = lambda x: math.isclose(-x[0] - x[1]**2 + 2 , 0,rel_tol=0.1, abs_tol=0.001)
st_fun4 = lambda x: math.isclose(x[1] + 2*x[2]**2 , 3,rel_tol=1e-1, abs_tol=0.01)
st_fun5 = lambda x: x[0] >= 0 and x[1] >= 0 and x[2] >= 0
judgy_st = lambda x: st_fun1(x) and st_fun2(x) and st_fun3(x) and st_fun4(x) and st_fun5(x)
注意第三个约束条件和第四个约束条件是一个等式,但在计算过程中是很难严格相等的,因此要用math.isclose()在误差允许的范围内近似相等,就相当于约等于。math.isclose()函数的功能如下:
math.isclose(a, b, rel_tol=1e-9, abs_tol=0.0)
其中,a和b是要比较的两个数值,rel_tol是相对容差(默认值为1e-9),abs_tol是绝对容差(默认值为0)。如果两个数值的差小于等于这两个容差的任意一个,则认为它们近似相等。
当 a-b<rel_tol 或者 |a-b|<abs_tol 时就可认为他们相等
最后还要添加一个judge_st()函数,判断我们的解是否满足所有的约束条件。
设计自变量的获取方式和更新方式
代码如下:
move = 5 #x的定义域
t_move = t/T #当前温度除以初始温度的比率,随着温度降低而减小
def get_x(x:list,lb,ub,t_move): #随机获取一个浮点数作为解
i = random.randint(0,2)
lb = x[i]-t_move*(x[i]-lb)
ub = x[i]+t_move*(ub-x[1])
# print(t_move)
x[i] = random.uniform(lb,ub)
return x
jud_n = 0
best_x = get_x(x,0,move,t/T) #这里先初始化一个最优解
# breakpoint()
while judgy_st(best_x) != True: #这里要判读得到解是否符合约束条件,如果不符合要重新取一个解
best_x = get_x(best_x,0,move,t/T)
jud_n += 1
best_y = target_fun(best_x)
print(f"f = {best_y},x = {best_x},第{jud_n}次循环")
y = best_y
这里定义了一个get_x(x:list,lb,ub,t_move)函数作为获取解的方式,它接收一个当前解x,并对其中的随机一项进行更改,得到新解然后返回。lb和ub是x的变化范围,lb是下限,ub是上限,也就是说x会在lb和ub之间进行变化。不过有时候,我们希望随着温度降低,x的变化范围越来越小,那就可以用t_move影响lb和ub,随着温度的降低t_move可以缩短lb和ub的范围。反之,如果不想x的变化范围受到温度的影响,可以把t_move设为1。
接着重点来了
我们这里需要大概猜一下,x的定义域,也就是x的变化范围。定义域猜的好的话,算法计算的就比较快,如果定义域猜的不好,算法计算的时候就会很慢,很难计算出结果。
看到这两个约束条件,如果x=[0,0,4]是符合条件的,如果x=[0,0,5]就不符合条件了,那就猜个x的定义域是0到5吧,设move等于5,也就是x会在0到5之间随机取一个浮点数作为解。这里只是猜个大概就行。
带入算法进行计算
直接放上完整代码:
import random
import numpy as np
import math
t = 2000
T = 2000
dt = 0.993
eps = 1e-14
n = 0
x = [0,0,0]
y = 100
target_fun = lambda x: x[0]**2 + x[1] ** 2 + x[2]**2 + 8
st_fun1 = lambda x: x[0]**2 - x[1] + x[2]**2 >= 0
st_fun2 = lambda x: x[0] + x[1]**2 +x[2]**2 <= 20
st_fun3 = lambda x: math.isclose(-x[0] - x[1]**2 + 2 , 0,rel_tol=0.1, abs_tol=0.001)
st_fun4 = lambda x: math.isclose(x[1] + 2*x[2]**2 , 3,rel_tol=1e-1, abs_tol=0.01)
st_fun5 = lambda x: x[0] >= 0 and x[1] >= 0 and x[2] >= 0
judgy_st = lambda x: st_fun1(x) and st_fun2(x) and st_fun3(x) and st_fun4(x) and st_fun5(x)
move = 4
t_move = t/T
def get_x(x:list,lb,ub,t_move):
i = random.randint(0,2)
lb = x[i]-t_move*(x[i]-lb)
ub = x[i]+t_move*(ub-x[1])
# print(t_move)
x[i] = random.uniform(lb,ub)
return x
jud_n = 0
best_x = get_x(x,0,move,t/T)
# breakpoint()
while judgy_st(best_x) != True:
best_x = get_x(best_x,0,move,t/T)
jud_n += 1
best_y = target_fun(best_x)
print(f"f = {best_y},x = {best_x},第{jud_n}次循环")
y = best_y
while n!= 200 :
# print(t/T)
dx = get_x(x,0,move,t/T)
judgy_n = 0 #这里用一个judgy_n判断内循环的运行次数,如果内循环运行到一定次数还没有计算出解就跳出循环,防止出现死循环
while judgy_st(dx) != True:
dx = get_x(dx,0,move,1)
judgy_n += 1
if judgy_n == 1000000:#这里设置内循环达到1000000次时跳出循环
break
if judgy_n == 1000000:
continue
dy = target_fun(dx)
if dy < y:
y = dy
x = dx
if dy < best_y:
best_y = dy
best_x = dx
print(f"f = {best_y},x = {best_x},第{n}次循环")
elif np.exp((y - dy) / t) > random.random():
y = dy
x = dx
# print(f"f = {y},x = {x},第{n}次循环")
t = t * dt
n += 1
print(n)
计算结果如下:
这是视频中给出的解:
可以看到用模拟退火计算的结果与视频的结果还是比较相近的,而且退火算法计算的最小值更小。
总结
这个算法是具有通用性的,只要把目标函数和约束条件改一下,就可以计算类似这样的非线性规划问题,当然也可以解决线性规划问题。