模擬退火演算法解100城市以上的TSP問題
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 11 100 56 9 101 56 17 102 56 25 103 56 33 104 56 41 105 64 41 106 72 41 107 72 49 108 56 49 109 48 51 110 56 57 111 56 65 112 48 63 113 48 73 114 56 73 115 56 81 116 48 83 117 56 89 118 56 97 119 104 97 120 104 105 121 104 113 122 104 121 123 104 129 124 104 137 125 104 145 126 116 145 127 124 145 128 132 145 129 132 137 130 140 137 131 148 137 132 156 137 133 164 137 134 172 125 135 172 117 136 172 109 137 172 101 138 172 93 139 172 85 140 180 85 141 180 77 142 180 69 143 180 61 144 180 53 145 172 53 146 172 61 147 172 69 148 172 77 149 164 81 150 148 85 151 124 85 152 124 93 153 124 109 154 124 125 155 124 117 156 124 101 157 104 89 158 104 81 159 104 73 160 104 65 161 104 49 162 104 41 163 104 33 164 104 25 165 104 17 166 92 9 167 80 9 168 72 9 169 64 21 170 72 25 171 80 25 172 80 25 173 80 41 174 88 49 175 104 57 176 124 69 177 124 77 178 132 81 179 140 65 180 132 61 181 124 61 182 124 53 183 124 45 184 124 37 185 124 29 186 132 21 187 124 21 188 120 9 189 128 9 190 136 9 191 148 9 192 162 9 193 156 25 194 172 21 195 180 21 196 180 29 197 172 29 198 172 37 199 172 45 200 180 45 201 180 37 202 188 41 203 196 49 204 204 57 205 212 65 206 220 73 207 228 69 208 228 77 209 236 77 210 236 69 211 236 61 212 228 61 213 228 53 214 236 53 215 236 45 216 228 45 217 228 37 218 236 37 219 236 29 220 228 29 221 228 21 222 236 21 223 252 21 224 260 29 225 260 37 226 260 45 227 260 53 228 260 61 229 260 69 230 260 77 231 276 77 232 276 69 233 276 61 234 276 53 235 284 53 236 284 61 237 284 69 238 284 77 239 284 85 240 284 93 241 284 101 242 288 109 243 280 109 244 276 101 245 276 93 246 276 85 247 268 97 248 260 109 249 252 101 250 260 93 251 260 85 252 236 85 253 228 85 254 228 93 255 236 93 256 236 101 257 228 101 258 228 109 259 228 117 260 228 125 261 220 125 262 212 117 263 204 109 264 196 101 265 188 93 266 180 93 267 180 101 268 180 109 269 180 117 270 180 125 271 196 145 272 204 145 273 212 145 274 220 145 275 228 145 276 236 145 277 246 141 278 252 125 279 260 129 280 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
當然，還可以繼續修改引數進行嘗試，進而得到更優的解