Simulated annealing algorithm solves the traveling quotient problem (MATLAB)

The simulated annealing algorithm solves the traveling quotient problem.
The new solution based on probability is mainly composed of two ways: two Exchange and three Exchange
The second exchange is to select two cities in the TSP Circuit to exchange directly
The three-switch is to select three points in the TSP Circuit, p1,p2,p3, and then exchange the city between P1,P2 directly with the corresponding length of the city before P3
The method for generating new solutions is not unique, as long as the new solution can be guaranteed to contain the solution space where the optimal solution resides.
There are two main cases of accepting a new solution:
The new solution is better than the historical optimal solution, then the new solution is accepted.
The new solution is better than the current solution, and the new solution is accepted with a certain probability, and with the decrease of temperature. The probability of acceptance will also decrease.
The following is the TSP code.

Clearclca =0.99;% parameters of the temperature decay functionT0 = the; tf =3; t = t0; Markov_length =10000;% Markov chain lengthcoordinates =[1    565.0   575.0;2     25.0   185.0;3    345.0   750.0;4    945.0   685.0;5    845.0   655.0;6    880.0   660.0;7     25.0   230.0;8    525.0  1000.0;9    580.0  1175.0;Ten   650.0  1130.0; One  1605.0   620.0; A  1220.0   580.0; -  1465.0   200.0; -  1530.0     5.0; the   845.0   680.0; -   725.0   370.0; -   145.0   665.0; -   415.0   635.0; +   510.0   875.0; -   560.0   365.0; +   300.0   465.0; A   520.0   585.0; at   480.0   415.0; -   835.0   625.0; -   975.0   580.0; -  1215.0   245.0; -  1320.0   315.0; -  1250.0   400.0; in   660.0   180.0; -   410.0   250.0; to   420.0   555.0; +   575.0   665.0; -  1150.0  1160.0; the   700.0   580.0; *   685.0   595.0; $   685.0   610.0;Panax Notoginseng   770.0   610.0; -   795.0   645.0; the   720.0   635.0; +   760.0   650.0; A   475.0   960.0; the    95.0   260.0; +   875.0   920.0; -   700.0   500.0; $   555.0   815.0; $   830.0   485.0; -  1170.0    65.0; -   830.0   610.0; the   605.0   625.0; -   595.0   360.0;Wuyi  1340.0   725.0; the  1740.0   245.0; ]; Coordinates (:,1) =[]; amount =size(Coordinates,1);% of Cities% calculate the distance matrix by vectorization methodDist_matrix =Zeros(Amount, amount); coor_x_tmp1 = coordinates (:,1) *ones(1, amount); coor_x_tmp2 =Coor_x_tmp1 '; coor_y_tmp1 = coordinates (:,2) *ones(1, amount); coor_y_tmp2 =Coor_y_tmp1 ';d Ist_matrix =sqrt((COOR_X_TMP1-COOR_X_TMP2). ^2+ ... (COOR_Y_TMP1-COOR_Y_TMP2). ^2); sol_new =1: Amount;% produces initial solutionThe sol_new is the new solution of each generation, and the sol_current is the current solution; Sol_best is the best solution in cooling ;E_current =inf; E_best =inf;% e_current is the circuit distance corresponding to the current solution;% e_new is the new solution of the loop distance;% e_best is the optimal solutionSol_current = sol_new; Sol_best = Sol_new;p =1; whileT>=tf forR=1: Markov_length% Markov chain length        % generates random disturbances        if(Rand<0.5)% randomly decides whether to make two or three exchanges            % of two exchangesInd1 =0; Ind2 =0; while(Ind1 = = ind2) Ind1 =Ceil(Rand.*amount); Ind2 =Ceil(Rand.*amount);EndTMP1 = Sol_new (IND1);            Sol_new (IND1) = Sol_new (Ind2); Sol_new (Ind2) = TMP1;Else            % of three exchangesInd1 =0; Ind2 =0; Ind3 =0; while(Ind1 = = Ind2) | | (Ind1 = = Ind3) ... | | (Ind2 = = Ind3) | | (ABS(ind1-ind2) = =1) Ind1 =Ceil(Rand.*amount); Ind2 =Ceil(Rand.*amount); Ind3 =Ceil(Rand.*amount);EndTMP1 = IND1;TMP2 = Ind2;tmp3 = Ind3;% ensure Ind1 < Ind2 < Ind3            if(Ind1 < Ind2) && (Ind2 < IND3);ElseIf(Ind1 < Ind3) && (Ind3 < ind2) Ind2 = Tmp3;ind3 = TMP2;ElseIf(Ind2 < Ind1) && (Ind1 < Ind3) Ind1 = Tmp2;ind2 = TMP1;ElseIf(Ind2 < Ind3) && (Ind3 < ind1) Ind1 = Tmp2;ind2 = Tmp3; Ind3 = TMP1;ElseIf(Ind3 < Ind1) && (Ind1 < ind2) Ind1 = Tmp3;ind2 = TMP1; Ind3 = TMP2;ElseIf(Ind3 < Ind2) && (Ind2 < ind1) Ind1 = Tmp3;ind2 = TMP2; Ind3 = TMP1;EndTmplist1 = Sol_new ((ind1+1):(ind2-1)); Sol_new ((ind1+1):(ind1+ind3-ind2+1) = ... sol_new ((ind2):(ind3)); Sol_new ((ind1+ind3-ind2+2): Ind3) = ... tmplist1;End        % Check to see if constraints are met        % calculated target function value (i.e. internal energy)E_new =0; for I=1: (amount-1) E_new = E_new + ... dist_matrix (sol_new (I), Sol_new (I+1));End        % again counting the distance from the last city to the first cityE_new = e_new + ... dist_matrix (sol_new (amount), Sol_new (1));ifE_new < E_current e_current = e_new; Sol_current = sol_new;ifE_new < E_best% Save the best solution in the cooling processE_best = e_new; Sol_best = sol_new;End        Else            % if the target function value of the new solution is less than the current solution,            % will only accept the new solution at a certain probability            if Rand<Exp(-(e_new-e_current)./t)                E_current = e_new; Sol_current = sol_new;ElseSol_new = sol_current;End        End    Endt=T.*a;% control parameter T (temperature) reduced to original a TimesEnd for I=1:length(coordinates) plot (coordinates (I,1), Coordinates (I,2),' r* '); Hold on;End; X=coordinates ([Sol_best sol_best (1)],1); Y=coordinates ([Sol_best sol_best (1)],2);p lot (x, y);disp(' optimal solution is: ')disp(sol_best)disp(' Shortest distance: ')disp(e_best)

The final image obtained is as follows, under visual conditions, should be satisfied with the condition (not necessarily the optimal, but also almost)

Simulated annealing algorithm solves the traveling quotient problem (MATLAB)