模拟退火算法改进

import numpy as np
import matplotlib.pyplot as plt
import math
import random
from scipy.stats import norm
from mpl_toolkits.mplot3d import Axes3D

# 目标函数
def Function(x, y):
    return -20 * np.exp(-0.2*np.sqrt(0.5*(x*x+y*y)))\
           -np.exp(0.5*(np.cos(2*np.pi*x)+np.cos(2*np.pi*y)))+20+np.e

# 初始化状态
def initN(low, high):
    '''
    随机生成一组状态取能量最低的状态
    :param low:
    :param high:
    :return:
    '''
    x = random.uniform(low, high)
    y = random.uniform(low, high)
    z = Function(x, y)
    for i in range(20):
        x1 = random.uniform(low, high)
        y1 = random.uniform(low, high)
        z1 = Function(x1, y1)
        if z1 < z:
            x = x1
            y = y1
            z = z1
    return x, y

# 绘制目标函数
def figureFuc(low, high):
    X = np.linspace(low, high, 500)
    Y = np.linspace(low, high, 500)
    XX, YY = np.meshgrid(X, Y)
    Z = -20 * np.exp(-0.2*np.sqrt(0.5*(XX*XX+YY*YY)))\
           -np.exp(0.5*(np.cos(2*np.pi*XX)+np.cos(2*np.pi*YY)))+20+np.e
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.plot_surface(XX, YY, Z, cmap=plt.get_cmap('rainbow'))
    plt.show()

# Metropolis准则接受新解
def Metropolis(detaF, T):  # detaF为f(n+1) - f(x)
    if detaF < 0:
        return 1
    else:
        pTk = math.exp(-detaF/T)
        if pTk > random.random():
            return 1
        else:
            return 0

# 利用正态分布产生新解
def get_normal_random_number(x,y,scale):  # 正态分布
    '''
    :param x: 均值
    :param y: 均值
    :param scale: 方差
    :return:
    '''
    fx = np.random.normal(x, scale)
    x = norm.ppf(fx)
    fy = np.random.normal(y, scale)
    y = norm.ppf(fy)
    return x, y

# 利用均匀分布产生新解
def get_uniform_random_number(low, high):
    '''
    :param low:
    :param high:
    :return:
    '''
    x = np.random.uniform(low, high)
    y = np.random.uniform(low, high)
    return x, y

# 冷却函数
def descT(T, k):
    # return T/np.log(1 + k)
    return 0.9*T

# 主函数
def startMain():
    # 初始化
    low = -5
    high = 5
    T = 10000
    Tmin = 10
    k = 1
    # figureFuc(low, high)  # 画图

    #x = random.uniform(low, high)
    #y = random.uniform(low, high)
    x, y = initN(low, high)
    z = Function(x, y)
    min_value = z
    record_value = []  # 用数组记录被接受的新解并绘图，方便分析
    while(T > Tmin and k <= 1000):
        x1, y1 = get_normal_random_number(x, y, 2)  # 利用正态分布产生新解
        # x1, y1 = get_uniform_random_number(low, high)  # 利用随机分布产生新解
        if x1 < low or x1 > high or y1 < low or y1 > high:   # 新解不在定义域内时跳出本次循环
            break
        z1 = Function(x1, y1)  # 计算新解的目标函数值
        deltaE = z1 - z
        min_value = min(min_value, z1)
        if Metropolis(deltaE, T) == 1:  # 接受按照Metropolis准则接受新解
            x = x1
            y = y1
            z = z1
            record_value.append(z)
        if deltaE > 0:
            T = descT(T, k)
        else:
            k += 1
    print('迭代到组后的解：', z)
    print('记录下的最优解：', min_value)

    #  打印解的变化曲线
    x=[i+1 for i in range(len(record_value))]
    plt.plot(x, record_value)
    plt.show()

if __name__ == "__main__":
    startMain()