Example of Simulated Annealing Algorithm and Its matlab implementation

The longitude and latitude of the enemy's 100 targets are known, as shown in table 1.
Table 1 longitude and latitude data table
Longitude Latitude
53.7121 15.3046 51.1758 0.0322 46.3253 28.2753 30.3313
56.5432 21.4188 10.8198 16.2529 22.7891 23.1045 10.1584
20.1050 15.4562 1.9451 0.2057 26.4951 22.1221 31.4847
26.2418 18.1760 44.0356 13.5401 28.9836 25.9879 38.4722
28.2694 29.0011 32.1910 5.8699 36.4863 29.7284 0.9718
8.9586 24.6635 16.5618 23.6143 10.5597 15.1178 50.2111
8.1519 9.5325 22.1075 18.5569 0.1215 18.8726 48.2077
31.9499 17.6309 0.7732 0.4656 47.4134 23.7783 41.8671
43.5474 3.9061 53.3524 26.7256 30.8165 13.4595 27.7133
23.9222 7.6306 51.9612 22.8511 12.7938 15.7307 4.9568
21.5051 24.0909 15.2548 27.2111 6.2070 5.1442 49.2430
17.1168 20.0354 34.1688 22.7571 9.4402 3.9200 11.5812
52.1181 0.4088 9.5559 11.4219 24.4509 6.5634 26.7213
37.5848 16.8474 35.6619 9.9333 24.4654 3.1644 0.7775
14.4703 13.6368 19.8660 15.1224 3.1616 4.2428 18.5245
58.6849 27.1485 39.5168 16.9371 56.5089 13.7090 52.5211
38.4300 8.4648 51.8181 23.0159 8.9983 23.6440 50.1156
13.7909 1.9510 34.0574 23.3960 23.0624 8.4319 19.9857
-274-40.8801 14.2978 58.8289 14.5229 18.6635 6.7436 52.8423
39.9494 29.5114 47.5099 24.0664 10.1121 27.2662 28.7812
8.0831 27.6705 9.1556 14.1304 53.7989 0.2199 33.6490
1.3496 16.8359 49.9816 6.0828 19.3635 17.6622 36.9545
15.7320 19.5697 11.5118 17.3884 44.0398 16.2635 39.7139
6.9909 23.1804 38.3392 19.9950 24.6543 19.6057 36.9980
4.1591 3.1853 40.1400 20.3030 23.9876 9.4030 41.1084
We have a base with longitude and latitude ). Assume that the speed of our plane is 1000 km/h.
We will send an airplane to depart from the base, review all enemy targets, and return to the original base. At each target of the enemy
The reconnaissance time of punctuation is not counted, and the time spent on the plane is calculated (assuming that our aircraft's cruise time can be sufficient ).
This is a traveling salesman problem. We sequentially numbered the base as 1, and the enemy's targets as 2, 3 ,...,
101. At last, we repeat the base number 102 (so that it is easy to calculate in the program ). Distance matrix 102 102
)(
× = Ij d D,
Where ij d represents j I, the distance between two points, 102, 2, 1, L = j I, where D is a real symmetric matrix. The problem is
Find a shortest path starting from and traveling through all the intermediate points to reach.
Given the geographical coordinates (longitude and latitude), we must calculate the actual distance between two points. Set
B A, the geographical coordinates of the two points are respectively ),(
1 1
Y x ,),(
2 2
Y x, over B A, the arc length of the big circle of two points is two points
. Take the earth center as the coordinate origin O, take the equator plane as the XOY plane, and take the level where the coil is located 0 degrees
The plane is the XOZ plane to create a three-dimensional Cartesian coordinate system. Then the Cartesian coordinates of B A and two points are:
) Sin, cos sin, cos (
1 1 1 1
Y R y x R
) Sin, cos sin, cos (
2 2 2 2
Y R y x R B
6370 = R indicates the Earth's radius.
B A, the actual distance between two points
Bytes
Bytes
Bytes
Bytes
Bytes
Bytes
Bytes
Bytes
Bytes
Bytes
=
OB
OB
R d
OA
OA
Arccos,
Simplified
] Sin cos) (arccos [cos
2 1 2 1 2 1
Y x R d + −=.
The simulated annealing algorithm is described as follows:
(1) solution space
A list of all the fixed start points and endpoints in the space S table: {102,101, 2, 1 L }.
} 102,} 101, 3, 2 {), (, 1 | ),,{(
102 101 2 1 102 1 = π the cyclic arrangement is L S
Each cyclic arrangement represents a loop of reconnaissance 100 targets.
I = π indicates the j point in the I-th reconnaissance,
The initial solution can be) 102, 2, 1 (L. In this paper, we use the Monte Carlo method to obtain a better initial solution.
(2) target functions
The target function is the path length or cost function of all targets. Our requirements
Σ
=
+
=
101
1
2 1 102
), (Min
I
I
D f
π L
An Iteration consists of the following three steps:
(3) generation of new solutions
① 2 Transformation Method
-275-any sequence number v u, (v u <) switches the order between u and v. The new path is:
102 1 1 1 1 1 π L ++ −v u v u
② 3 Transformation Method
Select the serial number v u, and w. After the path between u and v is inserted to w, the corresponding new path is
W v u <)
102 1 1 1 1 π L ++ −w v u w v u
(4) cost function difference
For the 2 transformation method, the path difference can be expressed
)()(
1 1 1 1 + −
+ −+ = △
V u
D f
π
(5) Acceptance Criteria
Bytes
Bytes
Bytes
≥△△−
<△
=
0)/exp (
0 1
F T f
F
P
If 0 <△f, the new path is accepted. Otherwise, use probability)/exp (T f △−accept the new path, that is, if
)/Exp (T f △−random number greater than 0 to 1 is accepted.
(6) Cooling
Use the selected cool-down factor α to reduce the temperature, that is, T α temperature, to get a new temperature. Here we take
999. 0 = α.
(7) conditions for ending
Use the selected termination Temperature
30
10
−
= E: Determine whether the annealing process is over. If e T <, the algorithm ends and the output is
Current status.
We compile the following matlab program as follows:
Clc, clear
Load sj.txt % load the data of the enemy's 100 targets. The data is stored in plain text according to the position in the table.
File sj.txt
X = sj (:,); x = x (:);
Y = sj (:,); y = y (:);
Sj = [x y];
D1 = [70, 40];
Sj = [d1; sj; d1];
Sj = sj * pi/180;
% Distance matrix d
D = zeros (102 );
For I = 1:101
For j = I + :102
Temp = cos (sj (I, 1)-sj (j, 1) * cos (sj (I, 2) * cos (sj (j, 2 )) + sin (sj (I, 2) * sin (sj (j, 2 ));
D (I, j) = 6370 * acos (temp );
End
End
D = d + d ';
S0 = []; Sum = inf;
Rand ('state', sum (clock ));
For j = 1:1000
S = [1 1 + randperm (100), 102];
Temp = 0;
-276-for I =
Temp = temp + d (S (I), S (I + 1 ));
End
If temp <Sum
S0 = S; Sum = temp;
End
End
E = 0.1 ^ 30; L = 20000; at = 0.999; T = 1;
% Annealing Process
For k = 1: L
% Generate new solution
C = 2 + floor (100 * rand (1, 2 ));
C = sort (c );
C1 = c (1); c2 = c (2 );
% Calculate the cost function value
Df = d (S0 (c1-1), S0 (c2) + d (S0 (c1), S0 (c2 + 1)-d (S0 (c1-1 ), s0 (c1)-d (S0 (c2), S0 (c2 + 1 ));
% Acceptance Criteria
If df <0
S0 = [S0 (1: c1-1), S0 (c2:-1: c1), S0 (c2 +)];
Sum = Sum + df;
Elseif exp (-df/T)> rand (1)
S0 = [S0 (1: c1-1), S0 (c2:-1: c1), S0 (c2 +)];
Sum = Sum + df;
End
T = T *;
If T <e
Break;
End
End
% Output cruise path and path length
S0, Sum