tips:题目来源于书籍,出处未知,如有侵权私删
问题背景:
位于缅因州米利诺克特的大北方纸业公司(Great Northern Paper Company)在佩诺布斯科特河(Penobscot River)上经营着一座水力发电站。水从大坝通过管道输送到发电站。水流过管道的速度因外部条件而异。
该电站有三个不同的水力涡轮机,每个水力涡轮机都有一个已知的(和独特的)功率函数,该函数给出了作为到达涡轮机的水流函数的发电量。进水可以以不同的体积分配给每个涡轮机,因此目标是确定如何在涡轮机之间分配水,以便在任何流速下提供最大的总能量产生。
使用实验证据和 伯努利方程,为每台涡轮机的功率输出以及允许的运行流量确定了以下二次模型:
题目要求:
- 如果所有三个涡轮机都被使用,我们希望确定 每个涡轮机的Qi的flow,这将提供最大的总能量产生。我们的局限性在于,流必须与总传入流相加,并且必须遵守给定的域限制。最大化总能量生产 KW1 + KW2 + KW3,要求最佳流量分配以及最大输出功率。
算法分析:
模拟退火算法是来源于固体退火原理,将固体加温,固体内部粒子......吧啦吧啦一大堆。
理科生请看图(没有针对文科生的意思)
tips:图来源于作者 Zachary_coding
好复杂,不如先看这个,是不是很熟悉(手动狗头)(tips: 此图来源百度百科)
整个算法就是模拟这个过程,那么思路不就是:随机初始温度,进行活化分子降能,无论是哪个初始点,都要降到产物(C)的过程。
那么难点就在于怎么进行 ’‘ 活化分子降能 ’‘,在降能的同时还要兼顾随机到最优的位置,以及降能方向,这就还要考虑到精度的问题。
(1)假如说初始温度在100,那就是在随机点进行随机步骤数(-1,1)之间的前提下,进行探索最优结果,迭代次数也是在一个正常范围。
那么温度如何下降?可以设置一个百分比,也就是下降速率,在此基础上进行下降。同时次过程和(1)过程,也就是模拟退火算法的一个基本思路。那么再来看这个图
是不是就清晰思路就清晰很多。
在这里还要考虑的一个点就是,局部最优和全局最优的问题,那么如何确保,搜索到的最优点是全局最优呢?那么(-1,1)这个就很巧妙,然后就是随机初始点也很巧妙。那么如何清晰的观察全局的迭代过程呢?是不是就得要一个迭代次数与目标函数值之间的关系图。那么思路就很明确。
metropolis准则:
tips:词条来自于百度学术
说白了就是,结果带进p这个式子里面当它大于(0,1)这个区间时,我们认为,这个答案是可接受的
题目求解:
tips:这里我详解代码过程
库导入:
import numpy as np
from matplotlib import pyplot as plt # 没有代码提示,个人认为还是很难写的(手动狗头)
# 运行图之后中文不显示和负号的问题
plt.rcParams['font.sans-serif'] = ['SimHei'] # 或者 ['Microsoft YaHei'] plt.rcParams['axes.unicode_minus'] = False
模型参数定义:
这里的的重要参数实际只有三个,初始温度,终止温度,下降速率。
T = 100
t = 1e-6
DeclineRate = 0.99
interate = 0
目标函数定义:
这里的约束条件处理很巧妙(假如用scipy模块,或者拉格朗日乘子。这个约束条件又怎么处理呢?)
def objective_function(Q1,Q2,Q): Q3 = Q - Q1 -Q2 kw1 = (-19.89 + 0.1277 * Q1 - 4.08e-5 * Q1**2) * (170 - 1.6e-6 * Q**2) kw2 = (-24.52 + 0.1358 * Q2 - 4.69e-5 * Q2**2) * (170 - 1.6e-6 * Q**2) kw3 = (-27.02 + 0.1380 * Q3 - 3.84e-5 * Q3**2) * (170 - 1.6e-6 * Q**2) kw = kw1 + kw2 + kw3 return kw
定义变元范围:
这个初始点是随机的,那么在定义域内的随机点,是不是就刚好满足这两个条件。
Q1 = np.random.uniform(250,1110) Q2 = np.random.uniform(250,1110) Q = np.random.uniform(750,3445)
函数主体:
# 存储目标函数值
kws = [ ]
while T > t:
# 先计算初始随机点的函数值
f = objective_function(Q1,Q2,Q)
# 随机探索(假如初始温度是1000,那么还是(-1,1)吗?)
Q1_new = Q1 + np.random.uniform(-1,1)
Q2_new = Q2 + np.random.uniform(-1,1)
Q_new = Q + np.random.uniform(-1,1)
# 在范围内才有用
if 250<=Q1_new<=1110:
if 250<=Q2_new<=1110:
if 750<=Q_new<= 3445:
f_new = objective_function(Q1_new,Q2_new,Q3_new)
if f>f_new : # 想想这里为什么是 f > f_new , 反过来是不是就是爬坡?
Q1 = Q1_new
Q2 = Q2_new
Q = Q_new
else :
p = np.exp(- (f_new - f) / T) # metropolis 准则
r = np.random.uniform(0,1)
if p > r:
Q1 = Q1_new
Q2 = Q2_new
Q = Q_new
# 假如这里直接存放,那应该是f还是f_new?为什么?
kws.append(objective_function(Q1,Q2,Q))
interate = interate + 1
plt.plot(range(1,interate+1),kws,color = 'r',label = 'interate_curve')
plt.xlabel('迭代次数')
plt.ylabel('目标函数函数值')
plt.titile('迭代曲线')
plt.legend()
plt.show()
完整代码:
# 模拟退火 # 随机初始点 --- 根据条件摆动 --- 判断条件 from matplotlib import pyplot as plt import numpy as np plt.rcParams['font.sans-serif'] = ['SimHei'] # 或者 ['Microsoft YaHei'] plt.rcParams['axes.unicode_minus'] = False # 变元域 # Q1 = (250,1110) # Q2 = (250,1110) # Q3 = (250,1225) # Q = 750,3445) # 定义目标函数 def qfunction(Q1,Q2,Q): Q3 = Q - Q1 -Q2 kw1 = (-19.89 + 0.1277 * Q1 - 4.08e-5 * Q1**2) * (170 - 1.6e-6 * Q**2) kw2 = (-24.52 + 0.1358 * Q2 - 4.69e-5 * Q2**2) * (170 - 1.6e-6 * Q**2) kw3 = (-27.02 + 0.1380 * Q3 - 3.84e-5 * Q3**2) * (170 - 1.6e-6 * Q**2) kw = kw1 + kw2 + kw3 return kw # 定义模型参数 T_initial = 1000 T_end = 1e-6 DeclineRate = 0.999 # 0.999 # metropolis准则 iterate = 0 Q1 = np.random.uniform(250,1110) Q2 = np.random.uniform(250,1110) Q = np.random.uniform(750,3445) # 存储元组 Q1_range = [] Q2_range = [] Q_range = [] T_range = [] # kw_range = [] while T_initial >T_end : # 初始随机数下的功率 kw = qfunction(Q1,Q2,Q) # 决策 - 变元摆动 Q1_new = Q1 + np.random.uniform(-10,10) Q2_new = Q2 + np.random.uniform(-10,10) Q_new = Q + np.random.uniform(-10,10) if 250 <= Q1_new <=1110 : if 250 <= Q2_new <=1110 : if 250 <= Q_new <=1225 : kw_new = qfunction(Q1_new ,Q2_new ,Q_new ) if kw_new > kw : # 这里其实是爬坡了 Q1 = Q1_new Q2 = Q2_new Q = Q_new else: p = np.exp( -(kw_new - kw) / T_initial) r = np.random.uniform(0,1) if p > r : # 满足 metropolis 准则 接受当前解 Q1 = Q1_new Q2 = Q2_new Q = Q_new T_initial = T_initial * DeclineRate kw_range.append(qfunction(Q1,Q2,Q)) iterate = iterate + 1 print(Q1,Q2,Q,qfunction(Q1,Q2,Q)) print("迭代次数:",iterate) plt.plot(range(1,iterate+1),kw_range,color = 'r',label = 'iterate_curve') plt.title('迭代曲线') plt.xlabel('迭代次数/n') plt.ylabel('输出功率/kw') plt.legend() plt.show()
运行结果:
这里的运行图不止这一个,还有另一种图。个人感觉是函数关系的原因,这个目标函数关系整个是单调的。看下面这个图:
还有就是,这里的 DeclineRate = 0.999,但是在这个下降速率为0.99999的时候,整个迭代曲线是平滑的直线,如下图(还是没想清楚原因)
总结:
博主纯小白,很多问题或者见解不周的地方还请读者大大们指正出来,有疑惑的地方提出来大家共同学习,共同进步。