Mathematic Modeling Method-Simulated annealing algorithm

First, Introduction

Hello everyone, see the title people should know what I am going to say today. "Simulated annealing algorithm", how to sound very burning feeling, haha and no, not at all, but very useful!! After reading this article you will know what I mean.

second, annealing phenomenon

First, let's look at what "annealing" is. It means the Demons are leaving, and the fire is going away. Almost haha, because the fire left the temperature is low. Okay, no kidding, seriously, annealing is just a matter of slowly heating the metal to a certain temperature, and then letting it cool slowly to room temperature. For example, Bo Master has a bowl of noodles, but the temperature is too high good to eat, then put on the table, and so it slowly cooled to the appropriate temperature can be eaten, the process is called annealing.

What does the simulated annealing algorithm have to do with this annealing phenomenon? In fact, there is a little relationship, the idea of simulated annealing algorithm is actually related to the phenomenon of annealing. Specifically how to listen to me slowly describe ha.

three, simulated annealing algorithm

First look at the figure below, if I have a function $y = f (x) $, the picture is the same as that. Now, I want to find the maximum value of this function . So what should we do? The simulated annealing algorithm is so thought of:

we first take a starting point $x = i$, red dashed line in the definition field of $x$, and get $y = f (i) $. Next, we can take the following strategy:

Let the red Line move one unit to the right, then make the following decision:

(1) If f (i + 1) > F (i), receive the move

It's good to understand that if you find that the new point you're moving to has a larger value, you're definitely receiving the move.

Now our red line is in the process of repeating (1) operation, and comes to the first maximum value of the function. Such as:

Now, as I proceed to the (1) step, I will find:

F (i + 1) < F (i), so here we are not going to accept the move. But here because we know what the figure looks like, we know that as long as the red Line passes through the "trough", we can find a real maximum value. But the computer will not know the appearance of the picture, then how to do it?

The core of the simulated annealing algorithm is here, it does not give up this new F (i + 1), it says, F (i + 1) you can be smaller than f (i), I give you the opportunity to give you a time to see if this f (i + 1) can be transformed into a "maximum value". So for f (i + 1) < F (i), do the following:

(2) if f (i + 1) < F (i), receive the movement at a certain probability and decrease the probability over time.

This probability is referred to the annealing process of metal smelting, which is the reason why this algorithm is called simulated annealing algorithm.

Iv. Metropolis Guidelines

In the annealing process of metal smelting, when the temperature is T, the probability that the energy difference for the DE particles is cooled is P (DE), which is expressed as:

$P \left (dE \right) = Exp\left ({\frac{de}{kt}}\right) $

where k is a constant, exp represents a natural index, and de < 0. The meaning of this formula is: the higher the temperature, the greater the probability of the cooling of the de by a single energy difference; the lower the temperature, the smaller the probability of cooling. And since DE is always less than 0 (otherwise it is not called annealing), so de/kt < 0, so the function value range of P (DE) is (0,1).

In the simulated annealing algorithm, we have a certain probability is defined as:

DC = f (i + 1)-F (i) < 0;

If exp (dc/t) >= Rand, the move is received, otherwise it is not received

the definition above is called　　 Metropolis Guidelines.

So. The process of the entire algorithm is this (pseudo code)

Objective: To find the minimum value of the system S (i): the       evaluation function value of the system at State y I: the current state of the system i + 1: The new state of the system       rate: Control the speed of cooling T0: The         initial temperature of the system (high temperature state) T_min:      The lower limit of temperature whi Le (T0 > t_min) {DE = s (i + 1)-S (i); if DE < 0//Receive move from S (i) to S (i+1) else if (exp (-de/t) > Random (0,1))//Receive from S (i　　　　　　　　) to S (i+1) Move else//do not receive movement from S (i) to S (i+1) T0 = R * T0; Cooling annealing (0 < R < 1; r larger, cooling slower; R smaller, cooling faster) i++;}

application of simulated annealing algorithm

There is such an ancient problem, known as the Traveling salesman problem (TSP, traveling salesman problem), is the famous problem in the field of mathematics, the problem is this: there is such a traveling businessman, he will visit N cities, but the path of each walk is limited, That is, each city can only visit once, and finally return to the original city of departure. How to select the path makes the total distance the minimum value.

It is famous because the problem has not yet found an effective algorithm. If the exact solution is to be solved only by the poor lifting of all the path combinations but as the number of cities increases its complexity will be high.

But fortunately, we have a simulated annealing algorithm, of course, it does not solve the problem accurately, but he can approximate the solution of the problem. The following are the specific practices:

1. Set the initial temperature T0, and randomly select a traverse path P (i) as the initial path, calculate its length L (P (i));

2. Randomly generate a new traverse path P (i+1) to calculate its length L (p (i + 1));

3. If L (p (i + 1)) < L (P (i)), the receiving P (i + 1) is a new path, otherwise the probability of the simulated annealing is received P (i + 1), and then cooled

4. Repeat steps 1 and 2 until the temperature reaches the minimum value of tmin.

　　There are a number of ways to generate a new traversal path, with 3 of them listed below:

1. Randomly select 2 nodes, exchanging the order of the 2 nodes in the path.

2. Randomly select 2 nodes to reverse the sequence of nodes between the 2 nodes in the path.

3. Randomly select 3 node m,n,k, then shift the node between node m and N to after node K.

six, simulated annealing algorithm matlab code

The following code is a simulated annealing algorithm for the TSP problem. The main functions are as follows:

Percent I Clear environment variable clear allclc%% ii.     Import city location data x = [16.4700 96.1000 16.4700 94.4400 20.0900 92.5400 22.3900 93.3700 25.2300 97.2400 22.0000 96.0500 20.4700 97.0200 17.2000 96.2900 16.3000 97.3800 14.0500 98.1200 16.5300 9 7.3800 21.5200 95.5900 19.4100 97.1300 20.0900 92.5500];%% iii.  Calculate the distance matrix D = Distance (X);    % calculated distance Matrix n = size (d,1); % of the city is the number of percent IV.                                                                  Initialization parameter T0 = 1e10;                                                                 % Initial temperature tf = 1e-30;                                                                     % terminating temperature%l = 2;                                                                    % of the number of iterations at each temperature q = 0.9;       % Cooling rate time = ceil (double (Solve ([Num2str (T0) ' * (0.9) ^x = ', Num2str (Tf)]));                                                                  % Count Iterations T0 * (0.9) ^x = Tfcount = 0;                                  % initialization iteration Count obj = zeros (time, 1);                     The% target value Matrix initializes path = zeros (time, n);                                                             % the optimal route matrix for each generation is initialized with a percent of V. Randomly generates an initial route S1 = Randperm (n);D Rawpath (S1, X)                                                             % draws the initial route disp (' a random value in the initial population: ') OutputPath (S1);                                                % The initial route is expressed in the form of an arrow Rlength = Pathlength (D, S1); % calculates the total length of the initial route disp ([' Total Distance: ', Num2str (rlength)]);                                                      Iterative optimizations while T0 > Tf count = count + 1;    % Update iteration Count%temp = zeros (l,n+1); % 1.    Generate new Solution S2 = Newanswer (S1); % 2.                                    The Metropolis rule determines whether to accept the new solution [S1, R] = Metropolis (S1, S2, D, T0); % Metropolis sampling algorithm% 3. Record the optimal route for each iteration if count = = 1 | |                                                         R < obj (count-1)% if the distance below the current temperature is less than the distance from the previous temperature, record the current distance of obj (count) = R; else obj (count) = obj(count-1);                                                     End Path (count,:) = S1;                                                            % record each iteration of the route T0 = q * T0; % at q rate of cooling end%% VII. Optimization process Iteration Diagram Figureplot (1:count, OBJ) xlabel (' Iteration count ') ylabel (' Distance ') title (' Optimization process ') grid on%% VIII. Draw the best path graph DrawPath (end,:), X)% The last line of the path matrix must be the optimal route, percent IX. The route and total distance of the output optimal solution Disp (' Optimal solution: ') S = Path (end,:);p = OutputPath (S);d ISP ([' Total Distance: ', Num2str (Pathlength (D, S)]);

in order to better maintain the code and make the code easier to read, some functions are functional packaging, as follows:

(1) DISTANCE.M: Calculate the distance between two cities

function D = Distance (citys) percent calculation 22 distance between cities% input: location coordinates of each city (citys)% output: 22 distance between cities (D) n = Size (citys, 1);D = zeros (n, n); for i = 1:n for    j = i + 1:n        D (i, j) = sqrt ((Citys (I, 1)-Citys (J, 1)) ^2 + (Citys (i, 2)-Citys (J, 2)) ^2);        D (j, i) = d (i, j);    EndEnd

(2) DRAWPATH.M: Show the path in the form of a graph

function DrawPath (route, citys) percent draw path function% input% route  to draw path   % citys  each city coordinate position% draw initial route figureplot ([Citys (Route, 1) ; Citys (Route (1), 1)], [Citys (Route, 2), Citys (Route (1), 2)], ' O '), grid on% give each location a sequence number for i = 1:size (citys, 1)    text (city S (i, 1), Citys (I, 2), ['   num2str (i)]), Endtext (Citys (Route (1), 1), Citys (Route (1), 2), '       starting point '), Text (Citys (route (end), 1), Citys (Route (end), 2), '       end Point ');

(3) Metropolis.m:metropolis criterion judgment

function [S,r] = Metropolis (S1, S2, D, T) percent input% S1:  Current solution% S2:   New Solution% D:    distance matrix (distance between 22 cities)% T:    Current temperature% output% s:
   
     Next current solution% R:   The route distance of the next current solution R1 = Pathlength (D, S1);         % calculation S1 route length n = length (S1);                 % City Number R2 = Pathlength (D,S2);          % calculates S2 route length dc = r2-r1;                   % difference in computing power if DC < 0                       If ability to lower accept new route    S = S2;                     % accepts new solution    R = R2;elseif exp (-dc/t) >= rand       % with exp (-dc/t) probability to accept the new route    S = S2;    R = R2;else                            % does not accept new route    S = S1;    R = R1;endend
   

(4) NEWANSWER.M: How to generate new solutions randomly

function S2 = Newanswer (S1) percent input% S1:    Current solution percent output% S2:   new Solution n = length (S1); S2 = S1;a = Round (rand (1, 2) * (n-1) + 1);    % produces two random positions for switching w = S2 (A (1)); S2 (A (1)) = S1 (A (2)); S2 (A (2)) = W;

(5) OUTPUTPATH.M: Describe the path direction in the form of arrows

function P = OutputPath (route) percent output path function% Input: path (r)% add starting point to R end (1) route = [Route, Route (1)];n = length (route);p = Num2str (Ro Ute (1)); For i = 2:n    p = [P, '-a ', Num2str (Route (i))];enddisp (p)

(6) PATHLENGTH.M: Calculates the length of the total path

function length = pathlength (D, route) percent calculates the path length between the start and end point% input:% D     22 distance between cities% route length = 0;n = Length (route); for I = 1: (n-1)    length = length + D (Route (i), route (i + 1)), Endlength = length + d (Route (N), Route (1));

after running the above code, the results are as follows:

A random value in the initial population:11-> 10-> 14-> 9-> 5-> 3-> 7-> 4-> 8-> 1-> 2-> 13-> 12-> 6 11 Total Distance: 62.7207 optimal solution:9-> 11-> 8-> 13-> 7-> 12-> 6-> 5-> 4-> 3-> 14-> 2-> ; 1-> 10-> 9 Total Distance: 29.3405

That is, by this algorithm, it gives us a shortest distance of 29.3405. We can not guarantee that this answer must be the shortest, this is the disadvantage of this algorithm. If you iterate enough times, the results get closer to the right results. However, an excessive number of iterations will take a long time to compute, so choose the number of iterations appropriately.

The path of the first random selection and the calculated optimal solution are as follows:

In addition, as the number of iterations increases, the total length of the change graph is as follows:

As you can see, the distance is constantly declining as we iterate. Therefore, the simulated annealing algorithm, in a sense, in the approximate solution of some problems, still played a great role.

Mathematic Modeling Method-Simulated annealing algorithm