Python Mathematical Modeling Learning Simulated Annealing Algorithm Multivariable Function Optimization Example Analysis
- 2021-12-09 09:35:18
- OfStack
1. Simulated annealing algorithm
Annealing is a process in which metal cools slowly from molten state and finally reaches the equilibrium state with the lowest energy. Based on the similarity between the optimization process and metal annealing process, the simulated annealing algorithm takes the optimization objective as the energy function, the solution space as the state space, and simulates the thermal motion of particles with random disturbance to solve the optimization problem ([1] KIRKPATRICK, 1988).
Simulated annealing algorithm is simple in structure and consists of temperature update function state generation function state acceptance function and termination criteria of inner cycle and outer cycle.
Temperature update function refers to the realization scheme of slow reduction of annealing temperature, also known as cooling schedule;
State generation function refers to the method of randomly generating new candidate solutions from current solutions.
State acceptance function refers to the mechanism of accepting candidate solutions, which usually adopts Metropolis criterion.
The external cycle is a temperature cycle controlled by the cooling schedule;
Inner cycle is the number of times a new solution is produced by cyclic iteration at every 1 temperature, which is also called Markov chain length.
The basic flow of simulated annealing algorithm is as follows:
(1) Initialization: initial temperature T, initial solution state s, iteration times L;
(2) For each temperature state, L cycles are repeated and new solutions are accepted probabilistically:
(3) Generating a new solution s 'from the current solution s by a transformation operation;
(4) Calculate the energy difference E, that is, the difference between the objective function of the new solution and the objective function of the original solution;
(5) If E
<
0 accepts s 'as the new current solution, otherwise accepts s' as the new current solution with probability exp (-® E/T);
(6) After L internal cycles are completed in each temperature state, the temperature T is lowered until the end temperature is reached.
2. Multivariable function optimization problem
Classical function optimization problems and combinatorial optimization problems are selected as test cases.
Question 1: The Schwefel test function is a complex multimodal function with a large number of local extremum regions.
In this paper, we take d=10, x= [-500,500], and the function is the global minimum f (X) = 0.0 at X= (420.9687, … 420.9687).
F (X) = 418.9829 × n-Σ (i = 1, n) [xi* sin (tick (xi))
The basic scheme of simulated annealing algorithm is as follows: the temperature is controlled to decay exponentially according to T (k) = a * T (k-1), and the attenuation coefficient is a; For example, formula (1) accepts a new solution according to Metropolis criterion. For problem 1 (Schwefel function), a new solution is generated by applying a random perturbation of normal distribution to one independent variable of the current solution.
3. Simulated annealing algorithm Python program
# 模拟退火算法 程序：多变量连续函数优化
# Program: SimulatedAnnealing_v1.py
# Purpose: Simulated annealing algorithm for function optimization
# Copyright 2021 YouCans, XUPT
# Crated：2021-04-30
# = 关注 Youcans，分享原创系列 https://blog.csdn.net/youcans =
# -*- coding: utf-8 -*-
import math # 导入模块
import random # 导入模块
import pandas as pd # 导入模块
import numpy as np # 导入模块 numpy，并简写成 np
import matplotlib.pyplot as plt # 导入模块 matplotlib.pyplot，并简写成 plt
from datetime import datetime
# 子程序：定义优化问题的目标函数
def cal_Energy(X, nVar):
# 测试函数 1： Schwefel 测试函数
# -500 <= Xi <= 500
# 全局极值：(420.9687,420.9687,...),f(x)=0.0
sum = 0.0
for i in range(nVar):
sum += X[i] * np.sin(np.sqrt(abs(X[i])))
fx = 418.9829 * nVar - sum
return fx
# 子程序：模拟退火算法的参数设置
def ParameterSetting():
cName = "funcOpt" # 定义问题名称
nVar = 2 # 给定自变量数量，y=f(x1,..xn)
xMin = [-500, -500] # 给定搜索空间的下限，x1_min,..xn_min
xMax = [500, 500] # 给定搜索空间的上限，x1_max,..xn_max
tInitial = 100.0 # 设定初始退火温度(initial temperature)
tFinal = 1 # 设定终止退火温度(stop temperature)
alfa = 0.98 # 设定降温参数，T(k)=alfa*T(k-1)
meanMarkov = 100 # Markov链长度，也即内循环运行次数
scale = 0.5 # 定义搜索步长，可以设为固定值或逐渐缩小
return cName, nVar, xMin, xMax, tInitial, tFinal, alfa, meanMarkov, scale
# 模拟退火算法
def OptimizationSSA(nVar,xMin,xMax,tInitial,tFinal,alfa,meanMarkov,scale):
# ====== 初始化随机数发生器 ======
randseed = random.randint(1, 100)
random.seed(randseed) # 随机数发生器设置种子，也可以设为指定整数
# ====== 随机产生优化问题的初始解 ======
xInitial = np.zeros((nVar)) # 初始化，创建数组
for v in range(nVar):
# random.uniform(min,max) 在 [min,max] 范围内随机生成1个实数
xInitial[v] = random.uniform(xMin[v], xMax[v])
# 调用子函数 cal_Energy 计算当前解的目标函数值
fxInitial = cal_Energy(xInitial, nVar)
# ====== 模拟退火算法初始化 ======
xNew = np.zeros((nVar)) # 初始化，创建数组
xNow = np.zeros((nVar)) # 初始化，创建数组
xBest = np.zeros((nVar)) # 初始化，创建数组
xNow[:] = xInitial[:] # 初始化当前解，将初始解置为当前解
xBest[:] = xInitial[:] # 初始化最优解，将当前解置为最优解
fxNow = fxInitial # 将初始解的目标函数置为当前值
fxBest = fxInitial # 将当前解的目标函数置为最优值
print('x_Initial:{:.6f},{:.6f},\tf(x_Initial):{:.6f}'.format(xInitial[0], xInitial[1], fxInitial))
recordIter = [] # 初始化，外循环次数
recordFxNow = [] # 初始化，当前解的目标函数值
recordFxBest = [] # 初始化，最佳解的目标函数值
recordPBad = [] # 初始化，劣质解的接受概率
kIter = 0 # 外循环迭代次数，温度状态数
totalMar = 0 # 总计 Markov 链长度
totalImprove = 0 # fxBest 改善次数
nMarkov = meanMarkov # 固定长度 Markov链
# ====== 开始模拟退火优化 ======
# 外循环，直到当前温度达到终止温度时结束
tNow = tInitial # 初始化当前温度(current temperature)
while tNow >= tFinal: # 外循环，直到当前温度达到终止温度时结束
# 在当前温度下，进行充分次数(nMarkov)的状态转移以达到热平衡
kBetter = 0 # 获得优质解的次数
kBadAccept = 0 # 接受劣质解的次数
kBadRefuse = 0 # 拒绝劣质解的次数
# ---内循环，循环次数为Markov链长度
for k in range(nMarkov): # 内循环，循环次数为Markov链长度
totalMar += 1 # 总 Markov链长度计数器
# ---产生新解
# 产生新解：通过在当前解附近随机扰动而产生新解，新解必须在 [min,max] 范围内
# 方案 1：只对 n元变量中的1个进行扰动，其它 n-1个变量保持不变
xNew[:] = xNow[:]
v = random.randint(0, nVar-1) # 产生 [0,nVar-1]之间的随机数
xNew[v] = xNow[v] + scale * (xMax[v]-xMin[v]) * random.normalvariate(0, 1)
# random.normalvariate(0, 1)：产生服从均值为0、标准差为 1 的正态分布随机实数
xNew[v] = max(min(xNew[v], xMax[v]), xMin[v]) # 保证新解在 [min,max] 范围内
# ---计算目标函数和能量差
# 调用子函数 cal_Energy 计算新解的目标函数值
fxNew = cal_Energy(xNew, nVar)
deltaE = fxNew - fxNow
# ---按 Metropolis 准则接受新解
# 接受判别：按照 Metropolis 准则决定是否接受新解
if fxNew < fxNow: # 更优解：如果新解的目标函数好于当前解，则接受新解
accept = True
kBetter += 1
else: # 容忍解：如果新解的目标函数比当前解差，则以1定概率接受新解
pAccept = math.exp(-deltaE / tNow) # 计算容忍解的状态迁移概率
if pAccept > random.random():
accept = True # 接受劣质解
kBadAccept += 1
else:
accept = False # 拒绝劣质解
kBadRefuse += 1
# 保存新解
if accept == True: # 如果接受新解，则将新解保存为当前解
xNow[:] = xNew[:]
fxNow = fxNew
if fxNew < fxBest: # 如果新解的目标函数好于最优解，则将新解保存为最优解
fxBest = fxNew
xBest[:] = xNew[:]
totalImprove += 1
scale = scale*0.99 # 可变搜索步长，逐步减小搜索范围，提高搜索精度
# ---内循环结束后的数据整理
# 完成当前温度的搜索，保存数据和输出
pBadAccept = kBadAccept / (kBadAccept + kBadRefuse) # 劣质解的接受概率
recordIter.append(kIter) # 当前外循环次数
recordFxNow.append(round(fxNow, 4)) # 当前解的目标函数值
recordFxBest.append(round(fxBest, 4)) # 最佳解的目标函数值
recordPBad.append(round(pBadAccept, 4)) # 最佳解的目标函数值
if kIter%10 == 0: # 模运算，商的余数
print('i:{},t(i):{:.2f}, badAccept:{:.6f}, f(x)_best:{:.6f}'.\
format(kIter, tNow, pBadAccept, fxBest))
# 缓慢降温至新的温度，降温曲线：T(k)=alfa*T(k-1)
tNow = tNow * alfa
kIter = kIter + 1
# ====== 结束模拟退火过程 ======
print('improve:{:d}'.format(totalImprove))
return kIter,xBest,fxBest,fxNow,recordIter,recordFxNow,recordFxBest,recordPBad
# 结果校验与输出
def ResultOutput(cName,nVar,xBest,fxBest,kIter,recordFxNow,recordFxBest,recordPBad,recordIter):
# ====== 优化结果校验与输出 ======
fxCheck = cal_Energy(xBest,nVar)
if abs(fxBest - fxCheck)>1e-3: # 检验目标函数
print("Error 2: Wrong total millage!")
return
else:
print("\nOptimization by simulated annealing algorithm:")
for i in range(nVar):
print('\tx[{}] = {:.6f}'.format(i,xBest[i]))
print('\n\tf(x):{:.6f}'.format(fxBest))
return
# 加粗样式
def main():
# 参数设置，优化问题参数定义，模拟退火算法参数设置
[cName, nVar, xMin, xMax, tInitial, tFinal, alfa, meanMarkov, scale] = ParameterSetting()
# print([nVar, xMin, xMax, tInitial, tFinal, alfa, meanMarkov, scale])
# 模拟退火算法
[kIter,xBest,fxBest,fxNow,recordIter,recordFxNow,recordFxBest,recordPBad] \
= OptimizationSSA(nVar,xMin,xMax,tInitial,tFinal,alfa,meanMarkov,scale)
# print(kIter, fxNow, fxBest, pBadAccept)
# 结果校验与输出
ResultOutput(cName, nVar,xBest,fxBest,kIter,recordFxNow,recordFxBest,recordPBad,recordIter)
if __name__ == '__main__':
main()
4. Program running results
x_Initial:-143.601793,331.160277, f(x_Initial):959.785447
i:0,t(i):100.00, badAccept:0.469136, f(x)_best:300.099320
i:10,t(i):81.71, badAccept:0.333333, f(x)_best:12.935760
i:20,t(i):66.76, badAccept:0.086022, f(x)_best:2.752498
...
i:200,t(i):1.76, badAccept:0.000000, f(x)_best:0.052055
i:210,t(i):1.44, badAccept:0.000000, f(x)_best:0.009448
i:220,t(i):1.17, badAccept:0.000000, f(x)_best:0.009448
improve:18
Optimization by simulated annealing algorithm:
x[0] = 420.807471
x[1] = 420.950005
f(x):0.003352
The above is Python mathematical modeling learning simulated annealing algorithm multivariable function optimization example analysis of the details, more information about mathematical modeling simulated annealing algorithm please pay attention to other related articles on this site!