基于模拟退火算法的多元极值问题

花生酥干

于 2024-04-23 09:25:38 发布

阅读量965

点赞数 31

文章标签：模拟退火算法算法机器学习

本文链接：https://blog.csdn.net/OAOOzz/article/details/138095744

版权

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的时候，整个迭代曲线是平滑的直线，如下图（还是没想清楚原因）

总结：

博主纯小白，很多问题或者见解不周的地方还请读者大大们指正出来，有疑惑的地方提出来大家共同学习，共同进步。

花生酥干

关注

31
点赞
踩
18

收藏

觉得还不错? 一键收藏
2
评论
基于模拟退火算法的多元极值问题

这里的的重要参数实际只有三个，初始温度，终止温度，下降速率。T = 100t = 1e-6这里的约束条件处理很巧妙（假如用scipy模块，或者拉格朗日乘子。这个约束条件又怎么处理呢？return kw这个初始点是随机的，那么在定义域内的随机点，是不是就刚好满足这两个条件。博主纯小白，很多问题或者见解不周的地方还请读者大大们指正出来，有疑惑的地方提出来大家共同学习，共同进步。
复制链接

扫一扫