1. Simulated annealing algorithm
Simulated annealing algorithm draws lessons from the idea of statistical physics. It is a simple and general heuristic optimization algorithm, and has probabilistic global optimization performance in theory. Therefore, it has been widely used in scientific research and engineering.
Annealing is a process in which the metal cools slowly from the molten state and finally reaches the equilibrium state with the lowest energy. Based on the similarity between the solution process of the optimization problem and the metal annealing process, the simulated annealing algorithm takes the optimization objective as the energy function, the solution space as the state space, and the random disturbance simulates the thermal motion of particles to solve the optimization problem ([1] KIRKPATRICK,1988).
The structure of simulated annealing algorithm is simple, which is composed of temperature update function, state generation function, state acceptance function and inner loop and outer loop termination criteria.
Temperature update function refers to the realization scheme of slow reduction of annealing temperature, also known as cooling schedule;
State generating function refers to the method of randomly generating new candidate solutions from the current solution;
State acceptance function refers to the mechanism of accepting candidate solutions, which usually adopts Metropolis criterion;
External circulation is a temperature cycle controlled by the cooling schedule;
The inner loop is the number of times that the loop iteration produces a new solution at each temperature, also known as the length of Markov chain.
The basic flow of simulated annealing algorithm is as follows:
(1) Initialization: initial temperature T, initial solution state s, number of iterations L;
(2) For each temperature state, repeat L cycles to generate and probabilistic accept a new solution:
(3) The new solution s' is generated from the current solution s by the 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, accept s' as the new current solution, otherwise accept s' as the new current solution with probability exp(- ∆ E/T);
(6) After completing L internal cycles in each temperature state, reduce the temperature T until the termination temperature is reached.
2. Multivariable function optimization problem
The classical function optimization problem and combinatorial optimization problem are selected as test cases.
Question 1: Schwefel test function is a complex multimodal function with a large number of local extremum regions.
F(X)=418.9829×n-∑(i=1,n)〖xi* sin(√(|xi|)) 〗
In this paper, we take d = 10, x = [- 500500], and the function is the global minimum f(X)=0.0 at X=(420.9687,... 420.9687).
The basic scheme of using simulated annealing algorithm: control the temperature to decay exponentially according to T(k) = a * T(k-1), and the attenuation coefficient is a; As shown in equation (1), accept the new solution according to Metropolis criteria. For problem 1 (Schwefel function), a new solution is generated by applying a random disturbance of normal distribution to an independent variable of the current solution.
3. Simulated annealing algorithm Python program
# Simulated annealing algorithm program: Multivariable continuous function optimization
# Program: SimulatedAnnealing_v1.py
# Purpose: Simulated annealing algorithm for function optimization
# Copyright 2021 YouCans, XUPT
# Crated: 2021-04-30
# -*- coding: utf-8 -*-
import math # Import module
import random # Import module
import pandas as pd # Import module
import numpy as np # Import the module numpy, which is abbreviated as np
import matplotlib.pyplot as plt # Import the module Matplotlib Pyplot and abbreviated as plot
from datetime import datetime
# Subroutine: defines the objective function of the optimization problem
def cal_Energy(X, nVar):
# Test function 1: Schwefel test function
# -500 <= Xi <= 500
# Global extremum: (420.9687420.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
# Subroutine: parameter setting of simulated annealing algorithm
def ParameterSetting():
cName = "funcOpt" # Define question name
nVar = 2 # Given the number of arguments, y = f (x1,... Xn)
xMin = [-500, -500] # Lower limit of given search space, x1_min,..xn_min
xMax = [500, 500] # Upper limit of given search space, x1_max,..xn_max
tInitial = 100.0 # Set initial annealing temperature
tFinal = 1 # Set stop temperature
alfa = 0.98 # Set cooling parameters, T(k)=alfa*T(k-1)
meanMarkov = 100 # Markov chain length, that is, the number of inner loop operations
scale = 0.5 # Define the search step size, which can be set to a fixed value or gradually reduced
return cName, nVar, xMin, xMax, tInitial, tFinal, alfa, meanMarkov, scale
# Simulated annealing algorithm
def OptimizationSSA(nVar,xMin,xMax,tInitial,tFinal,alfa,meanMarkov,scale):
# ======Initialize random number generator======
randseed = random.randint(1, 100)
random.seed(randseed) # The random number generator can set the seed or set it to a specified integer
# ======Initial solution of stochastic optimization problem======
xInitial = np.zeros((nVar)) # Initialize, create array
for v in range(nVar):
# random.uniform(min,max) randomly generates a real number in the range of [min,max]
xInitial[v] = random.uniform(xMin[v], xMax[v])
# Call sub function cal_Energy calculates the objective function value of the current solution
fxInitial = cal_Energy(xInitial, nVar)
# ======Simulated annealing algorithm initialization======
xNew = np.zeros((nVar)) # Initialize, create array
xNow = np.zeros((nVar)) # Initialize, create array
xBest = np.zeros((nVar)) # Initialize, create array
xNow[:] = xInitial[:] # Initialize the current solution and set the initial solution as the current solution
xBest[:] = xInitial[:] # Initialize the optimal solution and set the current solution as the optimal solution
fxNow = fxInitial # Set the objective function of the initial solution to the current value
fxBest = fxInitial # Set the objective function of the current solution as the optimal value
print('x_Initial:{:.6f},{:.6f},\tf(x_Initial):{:.6f}'.format(xInitial[0], xInitial[1], fxInitial))
recordIter = [] # Initialization, number of external cycles
recordFxNow = [] # Initialization, the objective function value of the current solution
recordFxBest = [] # Initialization, the objective function value of the optimal solution
recordPBad = [] # Initialization, acceptance probability of inferior solution
kIter = 0 # Number of external loop iterations, number of temperature states
totalMar = 0 # Total Markov chain length
totalImprove = 0 # fxBest improvement times
nMarkov = meanMarkov # Fixed length Markov chain
# ======Start simulated annealing optimization======
# The external circulation ends when the current temperature reaches the termination temperature
tNow = tInitial # Initialize current temperature
while tNow >= tFinal: # The external circulation ends when the current temperature reaches the termination temperature
# At the current temperature, a sufficient number of state transitions (nMarkov) are performed to achieve thermal equilibrium
kBetter = 0 # Number of times to obtain high-quality solutions
kBadAccept = 0 # Number of times to accept inferior solutions
kBadRefuse = 0 # Number of rejections of inferior solutions
# ---Inner loop, the number of cycles is the length of Markov chain
for k in range(nMarkov): # Inner loop, the number of cycles is the length of Markov chain
totalMar += 1 # Total Markov chain length counter
# ---Generate new solutions
# Generate a new solution: generate a new solution by randomly disturbing near the current solution. The new solution must be within the range of [min,max]
# Scheme 1: only one of the n variables is disturbed, and the other n-1 variables remain unchanged
xNew[:] = xNow[:]
v = random.randint(0, nVar-1) # Generate random numbers between [0,nVar-1]
xNew[v] = xNow[v] + scale * (xMax[v]-xMin[v]) * random.normalvariate(0, 1)
# random. Normal variable (0, 1): generate a random real number that obeys a normal distribution with a mean of 0 and a standard deviation of 1
xNew[v] = max(min(xNew[v], xMax[v]), xMin[v]) # Ensure that the new solution is within the range of [min,max]
# ---Calculate the objective function and energy difference
# Call sub function cal_Energy calculates the objective function value of the new solution
fxNew = cal_Energy(xNew, nVar)
deltaE = fxNew - fxNow
# ---Accept the new solution according to Metropolis guidelines
# Acceptance judgment: decide whether to accept the new solution according to Metropolis criteria
if fxNew < fxNow: # Better solution: if the objective function of the new solution is better than the current solution, the new solution is accepted
accept = True
kBetter += 1
else: # Tolerant solution: if the objective function of the new solution is worse than the current solution, the new solution will be accepted with a certain probability
pAccept = math.exp(-deltaE / tNow) # Calculate the state transition probability of the tolerant solution
if pAccept > random.random():
accept = True # Accept inferior solutions
kBadAccept += 1
else:
accept = False # Reject inferior solutions
kBadRefuse += 1
# Save new solution
if accept == True: # If the new solution is accepted, the new solution is saved as the current solution
xNow[:] = xNew[:]
fxNow = fxNew
if fxNew < fxBest: # If the objective function of the new solution is better than the optimal solution, the new solution is saved as the optimal solution
fxBest = fxNew
xBest[:] = xNew[:]
totalImprove += 1
scale = scale*0.99 # Variable search step size, gradually reduce the search range and improve the search accuracy
# ---Data sorting after internal circulation
# Complete the search of current temperature, save data and output
pBadAccept = kBadAccept / (kBadAccept + kBadRefuse) # Acceptance probability of inferior solution
recordIter.append(kIter) # Current external circulation times
recordFxNow.append(round(fxNow, 4)) # Objective function value of current solution
recordFxBest.append(round(fxBest, 4)) # Objective function value of optimal solution
recordPBad.append(round(pBadAccept, 4)) # Objective function value of optimal solution
if kIter%10 == 0: # Modular operation, remainder of quotient
print('i:{},t(i):{:.2f}, badAccept:{:.6f}, f(x)_best:{:.6f}'.\
format(kIter, tNow, pBadAccept, fxBest))
# Slowly cool down to a new temperature. Cooling curve: T(k)=alfa*T(k-1)
tNow = tNow * alfa
kIter = kIter + 1
# ======End the simulated annealing process======
print('improve:{:d}'.format(totalImprove))
return kIter,xBest,fxBest,fxNow,recordIter,recordFxNow,recordFxBest,recordPBad
# Result verification and output
def ResultOutput(cName,nVar,xBest,fxBest,kIter,recordFxNow,recordFxBest,recordPBad,recordIter):
# ======Verification and output of optimization results======
fxCheck = cal_Energy(xBest,nVar)
if abs(fxBest - fxCheck)>1e-3: # Objective function test
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
# main program
def main():
# Parameter setting, optimization problem parameter definition, simulated annealing algorithm parameter setting
[cName, nVar, xMin, xMax, tInitial, tFinal, alfa, meanMarkov, scale] = ParameterSetting()
# print([nVar, xMin, xMax, tInitial, tFinal, alfa, meanMarkov, scale])
# Simulated annealing algorithm
[kIter,xBest,fxBest,fxNow,recordIter,recordFxNow,recordFxBest,recordPBad] \
= OptimizationSSA(nVar,xMin,xMax,tInitial,tFinal,alfa,meanMarkov,scale)
# print(kIter, fxNow, fxBest, pBadAccept)
# Result verification and output
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
Copyright description:
Original works
Copyright 2021 YouCans, XUPT
Crated: 2021-05-01