2024美赛数学建模思路 – 案例：退火算法

本文介绍: # 0 赛题思路（赛题出来以后第一时间在CSDN分享）

## 0 赛题思路

（赛题出来以后第一时间在CSDN分享）

https://blog.csdn.net/dc_sinor?type=blog

在热力学上，退火（annealing）现象指物体逐渐降温的物理现象，温度愈低，物体的能量状态会低；够低后，液体开始冷凝与结晶，在结晶状态时，系统的能量状态最低。大自然在缓慢降温（亦即，退火）时，可“找到”最低能量状态：结晶。但是，如果过程过急过快，快速降温（亦称「淬炼」，quenching）时，会导致不是最低能态的非晶形。

如下图所示，首先（左图）物体处于非晶体状态。我们将固体加温至充分高（中图），再让其徐徐冷却，也就退火（右图）。加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小（此时物体以晶体形态呈现）。

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

class SA(object):

    def __init__(self, interval, tab='min', T_max=10000, T_min=1, iterMax=1000, rate=0.95):
        self.interval = interval                                    # 给定状态空间 - 即待求解空间
        self.T_max = T_max                                          # 初始退火温度 - 温度上限
        self.T_min = T_min                                          # 截止退火温度 - 温度下限
        self.iterMax = iterMax                                      # 定温内部迭代次数
        self.rate = rate                                            # 退火降温速度
        #############################################################
        self.x_seed = random.uniform(interval[0], interval[1])      # 解空间内的种子
        self.tab = tab.strip()                                      # 求解最大值还是最小值的标签: 'min' - 最小值；'max' - 最大值
        #############################################################
        self.solve()                                                # 完成主体的求解过程
        self.display()                                              # 数据可视化展示

    def solve(self):
        temp = 'deal_' + self.tab                                   # 采用反射方法提取对应的函数
        if hasattr(self, temp):
            deal = getattr(self, temp)
        else:
            exit('>>>tab标签传参有误："min"|"max"<<<')
        x1 = self.x_seed
        T = self.T_max
        while T >= self.T_min:
            for i in range(self.iterMax):
                f1 = self.func(x1)
                delta_x = random.random() * 2 - 1
                if x1 + delta_x >= self.interval[0] and x1 + delta_x <= self.interval[1]:   # 将随机解束缚在给定状态空间内
                    x2 = x1 + delta_x
                else:
                    x2 = x1 - delta_x
                f2 = self.func(x2)
                delta_f = f2 - f1
                x1 = deal(x1, x2, delta_f, T)
            T *= self.rate
        self.x_solu = x1                                            # 提取最终退火解

    def func(self, x):                                              # 状态产生函数 - 即待求解函数
        value = np.sin(x**2) * (x**2 - 5*x)
        return value

    def p_min(self, delta, T):                                      # 计算最小值时，容忍解的状态迁移概率
        probability = np.exp(-delta/T)
        return probability

    def p_max(self, delta, T):
        probability = np.exp(delta/T)                               # 计算最大值时，容忍解的状态迁移概率
        return probability

    def deal_min(self, x1, x2, delta, T):
        if delta < 0:                                               # 更优解
            return x2
        else:                                                       # 容忍解
            P = self.p_min(delta, T)
            if P > random.random(): return x2
            else: return x1

    def deal_max(self, x1, x2, delta, T):
        if delta > 0:                                               # 更优解
            return x2
        else:                                                       # 容忍解
            P = self.p_max(delta, T)
            if P > random.random(): return x2
            else: return x1

    def display(self):
        print('seed: {}nsolution: {}'.format(self.x_seed, self.x_solu))
        plt.figure(figsize=(6, 4))
        x = np.linspace(self.interval[0], self.interval[1], 300)
        y = self.func(x)
        plt.plot(x, y, 'g-', label='function')
        plt.plot(self.x_seed, self.func(self.x_seed), 'bo', label='seed')
        plt.plot(self.x_solu, self.func(self.x_solu), 'r*', label='solution')
        plt.title('solution = {}'.format(self.x_solu))
        plt.xlabel('x')
        plt.ylabel('y')
        plt.legend()
        plt.savefig('SA.png', dpi=500)
        plt.show()
        plt.close()


if __name__ == '__main__':
    SA([-5, 5], 'max')