模拟退火算法解100城市以上的TSP问题

模拟退火能有助于避免陷入局部最优解，理论上来说初始温度足够高，热平衡时间足够长即可搜索到全局最优解，用matlab编程实现如下。

clc;clear all;close all;coordinates=[ 1 288 149  2 288 129  3 270 133  4 256 141  5 256 157  6 246 157  7 236 169  8 228 169  9 228 161 10 220 169 11 212 169 12 204 169 13 196 169 14 188 169 15 196 161 16 188 145 17 172 145 18 164 145 19 156 145 20 148 145 21 140 145 22 148 169 23 164 169 24 172 169 25 156 169 26 140 169 27 132 169 28 124 169 29 116 161 30 104 153 31 104 161 32 104 169 33  90 165 34  80 157 35  64 157 36  64 165 37  56 169 38  56 161 39  56 153 40  56 145 41  56 137 42  56 129 43  56 121 44  40 121 45  40 129 46  40 137 47  40 145 48  40 153 49  40 161 50  40 169 51  32 169 52  32 161 53  32 153 54  32 145 55  32 137 56  32 129 57  32 121 58  32 113 59  40 113 60  56 113 61  56 105 62  48  99 63  40  99 64  32  97 65  32  89 66  24  89 67  16  97 68  16 109 69   8 109 70   8  97 71   8  89 72   8  81 73   8  73 74   8  65 75   8  57 76  16  57 77   8  49 78   8  41 79  24  45 80  32  41 81  32  49 82  32  57 83  32  65 84  32  73 85  32  81 86  40  83 87  40  73 88  40  63 89  40  51 90  44  43 91  44  35 92  44  27 93  32  25 94  24  25 95  16  25 96  16  17 97  24  17 98  32  17 99  44  11100  56   9101  56  17102  56  25103  56  33104  56  41105  64  41106  72  41107  72  49108  56  49109  48  51110  56  57111  56  65112  48  63113  48  73114  56  73115  56  81116  48  83117  56  89118  56  97119 104  97120 104 105121 104 113122 104 121123 104 129124 104 137125 104 145126 116 145127 124 145128 132 145129 132 137130 140 137131 148 137132 156 137133 164 137134 172 125135 172 117136 172 109137 172 101138 172  93139 172  85140 180  85141 180  77142 180  69143 180  61144 180  53145 172  53146 172  61147 172  69148 172  77149 164  81150 148  85151 124  85152 124  93153 124 109154 124 125155 124 117156 124 101157 104  89158 104  81159 104  73160 104  65161 104  49162 104  41163 104  33164 104  25165 104  17166  92   9167  80   9168  72   9169  64  21170  72  25171  80  25172  80  25173  80  41174  88  49175 104  57176 124  69177 124  77178 132  81179 140  65180 132  61181 124  61182 124  53183 124  45184 124  37185 124  29186 132  21187 124  21188 120   9189 128   9190 136   9191 148   9192 162   9193 156  25194 172  21195 180  21196 180  29197 172  29198 172  37199 172  45200 180  45201 180  37202 188  41203 196  49204 204  57205 212  65206 220  73207 228  69208 228  77209 236  77210 236  69211 236  61212 228  61213 228  53214 236  53215 236  45216 228  45217 228  37218 236  37219 236  29220 228  29221 228  21222 236  21223 252  21224 260  29225 260  37226 260  45227 260  53228 260  61229 260  69230 260  77231 276  77232 276  69233 276  61234 276  53235 284  53236 284  61237 284  69238 284  77239 284  85240 284  93241 284 101242 288 109243 280 109244 276 101245 276  93246 276  85247 268  97248 260 109249 252 101250 260  93251 260  85252 236  85253 228  85254 228  93255 236  93256 236 101257 228 101258 228 109259 228 117260 228 125261 220 125262 212 117263 204 109264 196 101265 188  93266 180  93267 180 101268 180 109269 180 117270 180 125271 196 145272 204 145273 212 145274 220 145275 228 145276 236 145277 246 141278 252 125279 260 129280 280 133];%数据集为280个城市的X,Y坐标a=0.99;t0=120;tf=3;t=t0;%设定初始温度和冷却温度Markov_length=50000;%MarKov链长度coordinates(:,1)=[];amount=size(coordinates,1);%计算距离矩阵dist_marix=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';dist_marix=sqrt((coor_x_tmp1-coor_x_tmp2).^2+(coor_y_tmp1-coor_y_tmp2).^2);sol_new=1:amount;%随机设定一个初始解E_current=inf;E_best=inf;sol_current=sol_new;sol_best=sol_new;p=1;k=1;while t>=tf    for i=1:Markov_length        %产生随机扰动        if(rand<0.5)          ind1=0;ind2=0;          while(ind1==ind2)              ind1=ceil(rand.*amount);              ind2=ceil(rand.*amount);          end          temp1=sol_new(ind1);          sol_new(ind1)=sol_new(ind2);          sol_new(ind2)=temp1;        else            ind1=0;ind2=0;ind3=0;            while(ind1==ind2)||(ind1==ind3)||(ind2==ind3)||(abs(ind2-ind1)==1)                ind1=ceil(rand.*amount);                ind2=ceil(rand.*amount);                ind3=ceil(rand.*amount);            end            temp1=ind1;temp2=ind2;temp3=ind3;            if (ind1<ind2)&&(ind2<ind3);                elseif (ind2<ind1)&&(ind3>ind1)                    ind2=temp1;ind1=temp2;                elseif (ind3<ind1)&&(ind2<ind3)                    ind1=temp2;ind2=temp3;ind3=temp1;                elseif (ind1<ind3)&&(ind3<ind2)                    ind3=temp2;ind2=temp3;                elseif (ind3<ind2)&&(ind2<ind1)                    ind1=temp3;ind3=temp1;                elseif (ind3<ind1)&&(ind1<ind2)                    ind1=temp3;ind2=temp1;ind3=temp2;            end            templist=sol_new((ind1+1):(ind2-1));            templist2=sol_new(ind3+1:end);            len1=size(templist,2);            sol_new(ind1+1:ind2-1)=[];            sol_new(ind3+1-len1:end)=[];            sol_new=[sol_new templist templist2];        end        %若含约束的规划问题这里需检查是否满足约束条件        %计算内能，即计算目标函数值        E_new=0;        for j=1:(amount-1)            E_new=E_new+dist_marix(sol_new(j),sol_new(j+1));        end        E_new= E_new+dist_marix(sol_new(amount),sol_new(1));        if E_new<E_current            E_current=E_new;            sol_current=sol_new;            if(E_new<E_best)                E_best=E_new;                sol_best=sol_new;            end        else            if(rand<exp(-(E_new-E_current)./t))                E_current=E_new;                sol_current=sol_new;            else                sol_new=sol_current;            end        end    end    Y(k)=E_new;    X(k)=k;    k=k+1;    t=t.*a;end    

收敛过程如图。

得到E_best =


   2.9724e+03                  

sol_best =


  1 至 20 列


   232   229   230   251   250   247   244   241   243   242     2     1   280     3   277   276   275   260   259   258


  21 至 40 列


   257   256   255   254   253   208   209   252   249   248   278   279     4     5     6     7     9     8    10    11


  41 至 60 列


    12    13    14    24    23    25    22    26    27    28    29    30   125   124   123   154   155   153   156   152


  61 至 80 列


   151   177   176   181   183   184   185   163   162   161   174   107   106   173   172   170   101   100    99    98


  81 至 100 列


    97    96    95    94    93    92    91    90    82    83    84    85    65    64    63    62    61    59    44    58


  101 至 120 列


    57    56    45    46    47    55    54    48    49    53    52    51    50    37    36    34    33    32    31    35


  121 至 140 列


    38    40    39    41    42    43    60   118   117   115   111   110   112    88   113    87   114   116    86    66


  141 至 160 列


    68    69    67    70    71    72    73    74    75    77    78    76    79    80    81    89   109   108   104   105


  161 至 180 列


   103   102   169   171   168   167   166   164   165   188   189   187   186   190   191   192   193   197   194   195


  181 至 200 列


   196   202   200   201   198   199   143   142   141   140   266   138   137   136   135   270   134   269   268   267


  201 至 220 列


   265   264   263   262   261   274   273   272    15   271    16    17    18   133    20    19   132   131    21   130


  221 至 240 列


   129   128   127   126   122   121   120   119   157   158   159   160   175   182   180   179   178   150   149   139


  241 至 260 列


   148   147   146   145   144   203   204   205   206   207   212   210   211   214   215   213   216   218   217   220


  261 至 280 列


   221   219   222   223   224   225   226   227   228   233   234   235   236   237   238   240   245   239   246   231

当然，还可以继续修改参数进行尝试，进而得到更优的解