模拟退火算法求函数最大、小值——python实现

模拟退火算法（Simulate Anneal，SA）是一种通用概率演算法，用来在一个大的搜寻空间内找寻命题的最优解。模拟退火是由S.Kirkpatrick, C.D.Gelatt和M.P.Vecchi在1983年所发明的。V.Černý在1985年也独立发明此演算法。模拟退火算法是解决TSP问题的有效方法之一。

模拟退火的出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法是一种通用的优化算法，其物理退火过程由加温过程、等温过程、冷却过程这三部分组成。（by 百度百科）

本文讨论用python实现简单模拟退火法，来求函数在某个区间范围内的最大/最小值。

---

1. 模拟退火算法

基本流程图：

下面是模拟退火算法的主程序simAnneal.py，实现了算一维，二维及多维函数在给定区间内的最大/最小值。

# -*- coding: utf-8 -*-

'''
=========================================
| kiterun |
| 2017/08/11 |
| [e-mail] |
=========================================
'''
from random import random
import math
import sys

class SimAnneal(object):
'''
 Simulated annealing algorithm 
 '''
def __init__(self, target_text='min'):
        self.target_text = target_text

def newVar(self, oldList, T):
'''
 :old : list
 :return : list, new solutions based on old solutions
 :T : current temperature
 '''
        newList = [i + (random()*2-1) for i in oldList]
return newList

def juge(self, func, new, old, T):
'''
 matropolise conditions: to get the maximun or minmun
 :new : new solution data from self.newX
 :old : old solution data
 :T : current temperature

 '''
        dE = func(new) - func(old) if self.target_text == 'max' else func(old) - func(new)
if dE >= 0:
            x, ans = new, func(new)
else:
if math.exp(dE/T) > random():
                x, ans = new,func(new)
else:
                x, ans = old, func(old)
return [x, ans]

class OptSolution(object):
'''
 find the optimal solution.

 '''
def __init__(self, temperature0=100, temDelta=0.98,
 temFinal=1e-8, Markov_chain=2000, 
 result=0, val_nd=[0]):
# initial temperature
        self.temperature0 = temperature0
# step factor for decreasing temperature
        self.temDelta = temDelta
# the final temperature
        self.temFinal = temFinal
# the Markov_chain length (inner loops numbers)
        self.Markov_chain = Markov_chain
# the final result
        self.result = result
# the initial coordidate values: 1D [0], 2D [0,0] ...
        self.val_nd = val_nd

def mapRange(self, oneDrange):
return (oneDrange[1]-oneDrange[0])*random() + oneDrange[0]


def soulution(self, SA_newV, SA_juge, juge_text, ValueRange, func):
'''
 calculate the extreme value: max or min value
 :SA_newV : function from class SimAnneal().newVar
 :SA_juge : function from class SimAnneal().juge_*
 :ValueRange : [], range of variables, 1D or 2D or 3D...
 :func : target function obtained from user

 '''
        Ti = self.temperature0
        ndim = len(ValueRange)
        f = max if juge_text=='max' else min
        loops =0

while Ti > self.temFinal:
            res_temp, val_temp = [], []
            preV = [[self.mapRange(ValueRange[j]) for i in range(self.Markov_chain)] for j in range(ndim)]
            newV = [ SA_newV(preV[j], T=Ti) for j in range(ndim)]

for i in range(self.Markov_chain):
                boolV = True
for j in range(ndim):
                    boolV &= (ValueRange[j][0]<= newV[j][i] <= ValueRange[j][1])
if boolV == True:
                    res_temp.append(SA_juge(new=[newV[k][i] for k in range(ndim)], func=func, old=[preV[k][i] for k in range(ndim)], T=Ti)[-1])
                    val_temp.append(SA_juge(new=[newV[k][i] for k in range(ndim)], func=func, old=[preV[k][i] for k in range(ndim)], T=Ti)[0])
else:
continue
                loops += 1

# get the index of extreme value
            idex = res_temp.index(f(res_temp))
            result_temp = f(self.result, f(res_temp))
# update the cooordidate of current extrema value 
            self.val_nd = self.val_nd if result_temp == self.result else val_temp[idex]
# update the extreme value
            self.result = result_temp
# update the current temperature
            Ti *= self.temDelta
            print(self.val_nd, self.result)

        print(loops)

#print('The extreme value = %f' %self.result[-1])

---

2. 利用这个程序我们来找一下一个二维函数的最大/最小值问题。

二维函数：
f(x,y)=y sin(2πx)+xcos(2πy)
example.py

# -*- coding: utf-8 -*-

'''
======================================
| kiterun |
| 2017/08/11 |
| [e-mail] |
======================================
'''

from random import random
import math
import sys
from time import time
from simAnneal_dev import SimAnneal
from simAnneal_dev import OptSolution
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

def func(w):
    x, = w
    fx = x + 10*math.sin(5*x) + 7*math.cos(4*x)
return fx

def func2(w):
    x, y = w
    fxy = y*np.sin(2*np.pi*x) + x*np.cos(2*np.pi*y)
return fxy

if __name__ == '__main__':

    targ = SimAnneal(target_text='max')
    init = -sys.maxsize # for maximun case
#init = sys.maxsize # for minimun case
    xyRange = [[-2, 2], [-2, 2]]
    xRange = [[0, 10]]
    t_start = time()

    calculate = OptSolution(Markov_chain=1000, result=init, val_nd=[0,0])
    output = calculate.soulution(SA_newV=targ.newVar, SA_juge=targ.juge, 
                        juge_text='max',ValueRange=xyRange, func=func2)

'''
 with open('out.dat', 'w') as f:
 for i in range(len(output)):
 f.write(str(output[i]) + '\n')
 '''
    t_end = time()
    print('Running %.4f seconds' %(t_end-t_start))

下图红点是找到的-2 <= x,y <= 2区间内最大值位置，程序运行结果在x,y=[1.7573972092698966, -1.9991309314219978]找到最大值3.7543430598423946：

下图红点是这该区间的最小值位置，运行结果在x,y = [-1.7612279505916202, -1.9998457808015955]找到最小值：-3.7560984312466141

---

3. 程序优化

上面的程序直接采用python的list做数据处理，众所周知numpy处理数据的速度快很多，所以对上面的程序进行小修改,主程序更改后为simAnneal_dev.py

# -*- coding: utf-8 -*-

'''
=========================================
| kiterun |
| 2017/08/11 |
| [e-mail] |
=========================================
'''
from random import random
import math
import sys
import numpy as np

class SimAnneal(object):
'''
 Simulated annealing algorithm 
 '''
def __init__(self, target_text='min'):
        self.target_text = target_text

def newVar(self, oldList, T):
'''
 :old : list
 :return : list, new solutions based on old solutions
 :T : current temperature
 '''
        newList = [i + (random()*2-1) for i in oldList]
return newList

def juge(self, func, new, old, T):
'''
 matropolise conditions: to get the maximun or minmun
 :new : new solution data from self.newX
 :old : old solution data
 :T : current temperature

 '''
        dE = func(new) - func(old) if self.target_text == 'max' else func(old) - func(new)
if dE >= 0:
            x, ans = new, func(new)
else:
if math.exp(dE/T) > random():
                x, ans = new, func(new)
else:
                x, ans = old, func(old)
return [x, ans]

class OptSolution(object):
'''
 find the optimal solution.

 '''
def __init__(self, temperature0=100, temDelta=0.98,
 temFinal=1e-8, Markov_chain=2000, 
 result=0, val_nd=[0]):
# initial temperature
        self.temperature0 = temperature0
# step factor for decreasing temperature
        self.temDelta = temDelta
# the final temperature
        self.temFinal = temFinal
# the Markov_chain length (inner loops numbers)
        self.Markov_chain = Markov_chain
# the final result
        self.result = result
# the initial coordidate values: 1D [0], 2D [0,0] ...
        self.val_nd = val_nd

# create unifrom distributed x,y ..., depend on value range
def mapRange(self, oneDrange):
return (oneDrange[1]-oneDrange[0])*random() + oneDrange[0]

def soulution(self, SA_newV, SA_juge, juge_text, ValueRange, func):
'''
 calculate the extreme value: max or min value
 :SA_newV : function from class SimAnneal().newVar
 :SA_juge : function from class SimAnneal().juge_*
 :ValueRange : [[],], range of variables, 1D or 2D or 3D...
 :func : target function obtained from user

 '''
        Ti = self.temperature0
        ndim = len(ValueRange)
        f = max if juge_text=='max' else min
        nf = np.amax if juge_text=='max' else np.amin
        loops = 0

while Ti > self.temFinal:
            res_temp = []
            preV = [[self.mapRange(ValueRange[j]) for i in range(self.Markov_chain)] for j in range(ndim)]
            newV = [ SA_newV(preV[j], T=Ti) for j in range(ndim)]

for i in range(self.Markov_chain):
                boolV = True
for j in range(ndim):
                    boolV &= (ValueRange[j][0]<= newV[j][i] <= ValueRange[j][1])
if boolV == True:
                    res_temp.append(SA_juge(new=[newV[k][i] for k in range(ndim)], 
                                    func=func, old=[preV[k][i] for k in range(ndim)], 
                                    T=Ti))
else:
continue
                loops += 1

# change list to numpy.array
            sol_temp = np.array(res_temp)
# find the extreme value
            extreme_temp = nf(sol_temp[:, 1])
# find the row No. and column No. of the extreme value
            re = np.where(sol_temp == extreme_temp)

            result_temp = f(self.result, extreme_temp)
# update the cooordidate of current extrema value 
            self.val_nd = self.val_nd if result_temp == self.result else sol_temp[re[0][0], 0]
# update the extreme value
            self.result = result_temp
# update the current temperature
            Ti *= self.temDelta

        output = [self.val_nd, self.result]
        print(output)
        print('Total loops = %d' %loops)
return output

运行的速度明显提升，快了将近一倍（没有做统计平均，只跑了单次结果）：