clear	
	clc
    x=xlsread('送貨.xlsx','C3:C32');
    y=xlsread('送貨.xlsx','D3:D32');
    num=xlsread('送貨.xlsx','A3:A32');
	a = 0.99;	% 溫度衰減函式的引數
	t0 = 97; tf = 3; t = t0;
	Markov_length = 10000;	% Markov鏈長度
	coordinates = [num,x,y];
	coordinates(:,1) = [];
	amount = size(coordinates,1); 	% 城市的數目
	% 通過向量化的方法計算距離矩陣
	dist_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';
	dist_matrix = abs(coor_x_tmp1-coor_x_tmp2)+ ...
					abs(coor_y_tmp1-coor_y_tmp2);

	sol_new = 1:amount;         % 產生初始解
% sol_new是每次產生的新解；sol_current是當前解；sol_best是冷卻中的最好解；
	E_current = inf;E_best = inf; 		% E_current是當前解對應的迴路距離；
% E_new是新解的迴路距離；
% E_best是最優解的
	sol_current = sol_new; sol_best = sol_new;          
	p = 1;

	while t>=tf
		for r=1:Markov_length		% Markov鏈長度
			% 產生隨機擾動
			if (rand < 0.5)	% 隨機決定是進行兩交換還是三交換
				% 兩交換
				ind1 = 0; ind2 = 0;
				while (ind1 == ind2)
					ind1 = ceil(rand.*amount);
					ind2 = ceil(rand.*amount);
				end
				tmp1 = sol_new(ind1);
				sol_new(ind1) = sol_new(ind2);
				sol_new(ind2) = tmp1;
			else
				% 三交換
				ind1 = 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);
				end
				tmp1 = ind1;tmp2 = ind2;tmp3 = ind3;
				% 確保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;
				end
				
				tmplist1 = 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

			%檢查是否滿足約束
			
			% 計算目標函式值（即內能）
			E_new = 0;
			for i = 1 : (amount-1)
				E_new = E_new + ...
					dist_matrix(sol_new(i),sol_new(i+1));
			end
			% 再算上從最後一個城市到第一個城市的距離
			E_new = E_new + ...
				dist_matrix(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
		t=t.*a;		% 控制引數t（溫度）減少為原來的a倍
	end

	disp('最優解為：')
	disp(sol_best)
	disp('最短距離：')
	disp(E_best)