通过python中的模拟退火找到噪声函数的全局最小值

Question

我试图从百位数百美元挑战中找到函数的全局最小值，问题 #4 作为模拟退火的练习。

作为我理解和编写代码方法的基础，我参考了在网上免费找到的全局优化算法第 3 版书。

因此，我最初想出了以下代码：

嘈杂的功能：

def noisy_func(x, y):
    return (math.exp(math.sin(50*x)) +
            math.sin(60*math.exp(y)) +
            math.sin(70*math.sin(x)) +
            math.sin(math.sin(80*y)) -
            math.sin(10*(x + y)) +
            0.25*(math.pow(x, 2) +
            math.pow(y, 2)))

用于改变值的函数：

def mutate(X_Value, Y_Value):

    mutationResult_X = X_Value + randomNumForInput()
    mutationResult_Y = Y_Value + randomNumForInput()

    while mutationResult_X > 4 or mutationResult_X < -4:
        mutationResult_X = X_Value + randomNumForInput()

    while mutationResult_Y > 4 or mutationResult_Y < -4:
        mutationResult_Y = Y_Value + randomNumForInput()

    mutationResults = [mutationResult_X, mutationResult_Y]
    return mutationResults

randomNumForInput 只是返回一个介于 4 和 -4 之间的随机数。 （搜索的间隔限制。）因此它相当于random.uniform(-4, 4)。

这是程序的核心功能。

def simulated_annealing(f):
    """Peforms simulated annealing to find a solution"""
    #Start by initializing the current state with the initial state
    #acquired by a random generation of a number and then using it
    #in the noisy func, also set solution(best_state) as current_state
    #for a start
    pCurSelect  = [randomNumForInput(),randomNumForInput()]
    current_state = f(pCurSelect[0],pCurSelect[1])
    best_state = current_state
    #Begin time monitoring, this will represent the
    #Number of steps over time
    TimeStamp = 1

    #Init current temp via the func, using such values as to get the initial temp
    initial_temp = 100
    final_temp = .1
    alpha = 0.001
    num_of_steps = 1000000
    #calculates by how much the temperature should be tweaked
    #each iteration
    #suppose the number of steps is linear, we'll send in 100
    temp_Delta = calcTempDelta(initial_temp, final_temp, num_of_steps)
    #set current_temp via initial temp
    current_temp = getTemperature(initial_temp, temp_Delta)

    #max_iterations = 100
    #initial_temp = get_Temperature_Poly(TimeStamp)

    #current_temp > final_temp
    while current_temp > final_temp:
        #get a mutated value from the current value
        #hence being a 'neighbour' value
        #with it, acquire the neighbouring state
        #to the current state
        neighbour_values = mutate(pCurSelect[0], pCurSelect[1])
        neighbour_state = f(neighbour_values[0], neighbour_values[1])

        #calculate the difference between the newly mutated
        #neighbour state and the current state
        delta_E_Of_States = neighbour_state - current_state

        # Check if neighbor_state is the best state so far

        # if the new solution is better (lower), accept it
        if delta_E_Of_States <= 0:
            pCurSelect = neighbour_values
            current_state = neighbour_state
            if current_state < best_state:
                best_state = current_state

        # if the new solution is not better, accept it with a probability of e^(-cost/temp)
        else:
            if random.uniform(0, 1) < math.exp(-(delta_E_Of_States) / current_temp):
                pCurSelect = neighbour_values
                current_state = neighbour_state
        # Here, we'd decrement the temperature or increase the timestamp, normally
        """current_temp -= alpha"""

        #print("Run number: " + str(TimeStamp) + " current_state = " + str(current_state) )
        #increment TimeStamp
        TimeStamp = TimeStamp + 1

        # calc temp for next iteration
        current_temp = getTemperature(current_temp, temp_Delta)

    #print("Iteration Count: " + str(TimeStamp))
    return best_state

alpha 不用于此实现，但是使用以下函数线性调节温度：

def calcTempDelta(T_Initial, T_Final, N):
    return((T_Initial-T_Final)/N)

def getTemperature(T_old, T_new):
    return (T_old - T_new)

这就是我实现本书第 245 页中描述的解决方案的方式。然而，这个实现并没有返回给我噪声函数的全局最小值，而是它附近的局部最小值之一。

我以这种方式实施解决方案的原因有两个：

1. 它作为线性温度调节的工作示例提供给我，因此是一个工作模板。

2. 虽然我试图理解本书第 248-249 页中列出的其他形式的温度调节，但我并不完全清楚变量“Ts”是如何计算的，即使在尝试查看了一些在本书引用的引用来源中，它对我来说仍然是深奥的。因此我想，我宁愿先尝试使这个“简单”的解决方案正常工作，然后再尝试其他温度淬火方法（对数、指数等）。

从那时起，我尝试了多种方法，通过代码的各种不同迭代来获取噪声函数的全局最小值，这对于一次发布在这里来说太过分了。我尝试过对这段代码进行不同的重写：

1. 减少每次迭代的随机滚动数，以便每次在更小的范围内搜索，这会导致更一致但仍然不正确的结果。

2. 以不同的增量变异，比如说，在 -1 和 1 之间，等等。效果相同。

3. 重写 mutate 以便通过某个步长检查与当前点的相邻点，并通过从当前点的 x/y 值添加/减少所述步长来检查相邻点，检查新生成的点之间的差异和当前点（基本上是 E 的增量），并返回适当的值，其中任何一个产生与当前函数的最短距离，因此是其最近的邻近邻居。

4. 减少搜索发生的间隔限制。

正是在这些解决方案中，涉及步长/减少限制/按象限检查邻居，我使用了由一些恒定 alpha 乘以 time_stamp 组成的运动。

我尝试过的这些和其他解决方案都没有奏效，要么产生更不准确的结果（尽管在某些情况下结果更一致），要么在一种情况下根本不起作用。

因此，我一定遗漏了一些东西，无论是与温度调节有关，还是我应该在算法中进行下一步（变异）的精确方式（公式）。

我知道有很多内容需要吸收和研究，但如果您能提供任何建设性的批评/帮助/建议，我将不胜感激。 如果展示其他解决方案尝试的代码位有任何帮助，我会在需要时发布它们。

Answer 1

重要的是你要跟踪你在做什么。 我在frigidum上放了一些重要的提示

alpha 冷却通常运作良好，它确保您不会加速通过有趣的最佳位置，大约 0.1 的提案被接受。

确保你的建议不要太粗略，我举了一个例子，我只改变 x 或 y，但从不改变两者。这个想法是退火将采取最好的方式，或者进行巡回演出，让方案决定。

我为算法使用了 frigidum 包，但它与您的代码几乎相同。另请注意，我有 2 个提案，一个大变化和一个小变化，组合通常效果很好。

最后，我注意到它跳得很厉害。一个小的变化是在你进入最后 5% 的冷却之前选择迄今为止最好的。

我使用/安装了 frigidum

!pip install frigidum

并进行了小改动以使用 numpy 数组；

import math

def noisy_func(X):
    x, y = X
    return (math.exp(math.sin(50*x)) +
            math.sin(60*math.exp(y)) +
            math.sin(70*math.sin(x)) +
            math.sin(math.sin(80*y)) -
            math.sin(10*(x + y)) +
            0.25*(math.pow(x, 2) +
            math.pow(y, 2)))


import frigidum
import numpy as np
import random

def random_start():
    return np.random.random( 2 ) * 4

def random_small_step(x):
    if np.random.random() < .5:
        return np.clip( x + np.array( [0, 0.02 * (random.random() - .5)] ), -4,4)
    else:
        return np.clip( x + np.array( [0.02 * (random.random() - .5), 0] ), -4,4)


def random_big_step(x):
    if np.random.random() < .5:
        return np.clip( x + np.array( [0, 0.5 * (random.random() - .5)] ), -4,4)
    else:
        return np.clip( x + np.array( [0.5 * (random.random() - .5), 0] ), -4,4)

local_opt = frigidum.sa(random_start=random_start, 
                        neighbours=[random_small_step, random_big_step], 
                        objective_function=noisy_func, 
                        T_start=10**2, 
                        T_stop=0.00001, 
                        repeats=10**4, 
                        copy_state=frigidum.annealing.copy)

---

上面的输出是

---
Neighbour Statistics: 
(proportion of proposals which got accepted *and* changed the objective function)
   random_small_step                : 0.451045
   random_big_step                  : 0.268002
---
(Local) Minimum Objective Value Found: 
   -3.30669277

使用上面的代码，有时我会低于 -3，但我也注意到有时它会在 -2 附近找到一些东西，而不是卡在最后一个阶段。

因此，一个小的调整是重新退火最后阶段的退火，使用迄今为止最好的。

希望对您有所帮助，如有任何问题，请告诉我。