Learning Notes of Simulated Annealing Algorithms

Concept:

Internal energy degradation law of metallurgical metals.

Advantages: Avoid reaching local optimum in reasoning, and achieve global optimum through probability of random events.

Algorithmic principles:

Random heuristic search algorithm.

It consists of three parts: initial solution, solution space and objective function. Meaning is initial state, midway process, optimization goal.

(1) Initial solution:

The first solution is the starting point of iteration.

(2) Solution space:

The set of feasible solutions.

(3) The objective function:

Description of the optimal solution.

The flow chart of the algorithm is as follows:

Cooling probability:

When the system temperature is TTT, the probability of *** cooling with energy difference of dEdEdE *** is as follows:
p(dE)=edEkT p(dE)=e^{\frac{dE}{kT}} p(dE)=ekTdE​
When calculating the minimum value, the strategy is as follows:
Δf&lt;0→BS→1Δf≥0→TS→p(−Δf) \Delta f &lt; 0 \rightarrow BS \rightarrow 1 \\ \Delta f \geq 0 \rightarrow TS \rightarrow p(-\Delta f) Δf<0→BS→1Δf≥0→TS→p(−Δf)
When the maximum value is obtained, the positive and negative symbols can be changed.

C++ code

Solving Traveling Salesman Problem by Simulated Annealing

#include <iostream>
#include <string.h>
#include <stdlib.h>
#include <algorithm>
#include <stdio.h>
#include <time.h>
#include <math.h>

#define N     30      //City quantity
#define T     3000    //initial temperature
#define EPS   1e-8    //Termination temperature
#define DELTA 0.98    //Temperature attenuation rate

#define LIMIT 1000   //Upper limit of probability selection
#define OLOOP 20    //Number of external circulation
#define ILOOP 100   //Number of internal cycles

using namespace std;

//Define route structure
struct Path
{
    int citys[N];
    double len;
};

//Defining Urban Point Coordinates
struct Point
{
    double x, y;
};

Path bestPath;        //Recording Optimal Path
Point p[N];       //The coordinates of each city
double w[N][N];   //Path Length Between Two Cities
int nCase;        //Test times

double dist(Point A, Point B)
{
    return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));
}

void GetDist(Point p[], int n)
{
    for(int i = 0; i < n; i++)
        for(int j = i + 1; j < n; j++)
            w[i][j] = w[j][i] = dist(p[i], p[j]);
}

void Input(Point p[], int &n)
{
    scanf("%d", &n);
    for(int i = 0; i < n; i++)
        scanf("%lf %lf", &p[i].x, &p[i].y);
}

void Init(int n)
{
    nCase = 0;
    bestPath.len = 0;
    for(int i = 0; i < n; i++)
    {
        bestPath.citys[i] = i;
        if(i != n - 1)
        {
            printf("%d--->", i);
            bestPath.len += w[i][i + 1];
        }
        else
            printf("%d\n", i);
    }
    printf("\nInit path length is : %.3lf\n", bestPath.len);
    printf("-----------------------------------\n\n");
}

void Print(Path t, int n)
{
    printf("Path is : ");
    for(int i = 0; i < n; i++)
    {
        if(i != n - 1)
            printf("%d-->", t.citys[i]);
        else
            printf("%d\n", t.citys[i]);
    }
    printf("\nThe path length is : %.3lf\n", t.len);
    printf("-----------------------------------\n\n");
}

Path GetNext(Path p, int n)
{
    Path ans = p;
    int x = (int)(n * (rand() / (RAND_MAX + 1.0)));
    int y = (int)(n * (rand() / (RAND_MAX + 1.0)));
    while(x == y)
    {
        x = (int)(n * (rand() / (RAND_MAX + 1.0)));
        y = (int)(n * (rand() / (RAND_MAX + 1.0)));
    }
    swap(ans.citys[x], ans.citys[y]);
    ans.len = 0;
    for(int i = 0; i < n - 1; i++)
        ans.len += w[ans.citys[i]][ans.citys[i + 1]];
    cout << "nCase = " << nCase << endl;
    Print(ans, n);
    nCase++;
    return ans;
}

void SA(int n)
{
    double t = T;
    srand((unsigned)(time(NULL)));
    Path curPath = bestPath;
    Path newPath = bestPath;
    int P_L = 0;
    int P_F = 0;
    while(1)       //External circulation, main update parameter t, simulated annealing process
    {
        for(int i = 0; i < ILOOP; i++)    //Internal circulation to find the optimal value at a certain temperature
        {
            newPath = GetNext(curPath, n);
            double dE = newPath.len - curPath.len;
            if(dE < 0)   //If you find a better value, update it directly
            {
                curPath = newPath;
                P_L = 0;
                P_F = 0;
            }
            else
            {
                double rd = rand() / (RAND_MAX + 1.0);
                //If you find a worse solution than the current one, accept it with a certain probability, and the probability will become smaller and smaller.
                if(exp(dE / t) > rd && exp(dE / t) < 1)
                    curPath = newPath;
                P_L++;
            }
            if(P_L > LIMIT)
            {
                P_F++;
                break;
            }
        }
        if(curPath.len < bestPath.len)
            bestPath = curPath;
        if(P_F > OLOOP || t < EPS)
            break;
        t *= DELTA;
    }
}

int main(int argc, const char * argv[]) {
    freopen("TSP.data", "r", stdin);
    int n;
    Input(p, n);
    GetDist(p, n);
    Init(n);
    SA(n);
    Print(bestPath, n);
    printf("Total test times is : %d\n", nCase);
    return 0;
}

The key to this code is GetNext(Path p, int n),SA(int n) function, GetNext(Path p, int n) for random selection of solutions, SA(int n) for selection of whether to accept points, how to accept points, and the number of iterations and termination conditions.

Python code

import numpy as np
import matplotlib.pyplot as plt
import random

class SA(object):

    def __init__(self, interval, tab='min', T_max=10000, T_min=1, iterMax=1000, rate=0.95):
        self.interval = interval                                    # Given state space - that is, the space to be solved
        self.T_max = T_max                                          # Initial Annealing Temperature-Temperature Upper Limit
        self.T_min = T_min                                          # Cut-off Annealing Temperature-Lower Temperature Limit
        self.iterMax = iterMax                                      # Number of internal iterations at constant temperature
        self.rate = rate                                            # Annealing cooling rate
        #############################################################
        self.x_seed = random.uniform(interval[0], interval[1])      # Seeds in Solution Space
        self.tab = tab.strip()                                      # Label for maximum or minimum:'min'- minimum;'max' - Maximum
        #############################################################
        self.solve()                                                # Complete the Solving Process of the Subject
        self.display()                                              # Visual display of data
        
    def solve(self):
        temp = 'deal_' + self.tab                                   # Extraction of corresponding functions by reflection method
        if hasattr(self, temp):
            deal = getattr(self, temp)
        else:
            exit('>>>tab The label is incorrect:"min"|"max"<<<')  
        x1 = self.x_seed
        T = self.T_max
        while T >= self.T_min:
            for i in range(self.iterMax):
                f1 = self.func(x1)
                delta_x = random.random() * 2 - 1
                if x1 + delta_x >= self.interval[0] and x1 + delta_x <= self.interval[1]:   # Binding random solutions to a given state space
                    x2 = x1 + delta_x
                else:
                    x2 = x1 - delta_x
                f2 = self.func(x2)
                delta_f = f2 - f1
                x1 = deal(x1, x2, delta_f, T)
            T *= self.rate
        self.x_solu = x1                                            # Extraction of final annealing solution       
        
    def func(self, x):                                              # STATE GENERATING FUNCTION-THE FUNCTION TO BE SOLUTED
        value = np.sin(x**2) * (x**2 - 5*x)
        return value
        
    def p_min(self, delta, T):                                      # When calculating the minimum, the state transition probability of the tolerant solution
        probability = np.exp(-delta/T)
        return probability
        
    def p_max(self, delta, T):
        probability = np.exp(delta/T)                               # The state transition probability of the tolerant solution when calculating the maximum value
        return probability
        
    def deal_min(self, x1, x2, delta, T):
        if delta < 0:                                               # Better solution
            return x2
        else:                                                       # Tolerance solution
            P = self.p_min(delta, T)
            if P > random.random(): return x2
            else: return x1
            
    def deal_max(self, x1, x2, delta, T):
        if delta > 0:                                               # Better solution
            return x2
        else:                                                       # Tolerance solution
            P = self.p_max(delta, T)
            if P > random.random(): return x2
            else: return x1
        
    def display(self):
        print('seed: {}\nsolution: {}'.format(self.x_seed, self.x_solu))
        plt.figure(figsize=(6, 4))
        x = np.linspace(self.interval[0], self.interval[1], 300)
        y = self.func(x)
        plt.plot(x, y, 'g-', label='function')
        plt.plot(self.x_seed, self.func(self.x_seed), 'bo', label='seed')
        plt.plot(self.x_solu, self.func(self.x_solu), 'r*', label='solution')
        plt.title('solution = {}'.format(self.x_solu))
        plt.xlabel('x')
        plt.ylabel('y')
        plt.legend()
        plt.savefig('SA.png', dpi=500)
        plt.show()
        plt.close()

        
if __name__ == '__main__':
    SA([-5, 5], 'max')

Code excerpts
https://www.cnblogs.com/xxhbdk/p/9192750.html
https://www.cnblogs.com/ranjiewen/p/6084052.html