引言
TSP(旅行商问题)，一个旅行商需要途径14个城市，不可以返回已经走过的城市。每个城市有对于的坐标，求旅行商走完14个城市的最短路径。
1. 贪婪算法
贪婪算法即每一步都采取最优解，即每一次都选择与所在城市距离最短的城市。初始城市随机设定。该方法简洁，但是寻找最优解的精度不够。
当城市数量达到一定数量后，容易陷入局部最优。
%% 贪心算法解决TSP问题%% 坐标数据初始化clear; clc; close allload CityPosition1.mat%% 建立距离矩阵amount = size(X,1); % 城市的数目% 通过向量化的方法计算距离矩阵% 计算矩阵中两两点之间的距离，可以通过转置相减的方法coor_x_tmp1 = repmat(X(:,1),1,amount);coor_x_tmp2 = coor_x_tmp1';coor_y_tmp1 = repmat(X(:,2),1,amount);coor_y_tmp2 = coor_y_tmp1';D = sqrt((coor_x_tmp1-coor_x_tmp2).^2 + (coor_y_tmp1-coor_y_tmp2).^2);%% 贪婪算法求解最短距离city_order = linspace(1,14,14); % 城市标号stp = randperm(14,1); % 随机选择一个点作为初始点city_order(stp) = []; % 从城市标号中删去该点visit = [stp]; % 用于存放已经访问的城市obj = 0; % 用于保存最短路径for i = 1:13dis = []; % 用来存放将要访问城市的距离for k = [city_order]distance = D(stp,k);dis = [dis,distance];enddis_min = min(dis);obj = obj + dis_min;m = find(dis == dis_min);visit = [visit,city_order(m)];stp = city_order(m);city_order(m) = [];end
2. 模拟退火算法
模拟退火算法(Simulated Annealing，SA)最早的思想是由N. Metropolis等人于1953年提出。1983 年,S. Kirkpatrick 等成功地将退火思想引入到组合优化领域。它是基于Monte-Carlo迭代求解策略的一种随机寻优算法，其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法从某一较高初温出发，伴随温度参数的不断下降,结合概率突跳特性解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。
求解TSP问题中，初始解为1-14的序列。新解的产生有两种方式，一种为二交换，即交换其中两个数字的位置；另一种为三交换。其他算法部分与常规的模拟退火相同。
在编写模拟退火算法时，最好加入CV函数，本题没用到约束条件，因此CV（罚函数）=0。遇到有约束条件的问题，采用罚函数处理时最简洁的。
%% 模拟退火解决旅行商问题%% 坐标数据初始化clear; clc; close allload CityPosition1.mat%% 建立距离矩阵city_order = linspace(1,14,14); % 城市标号amount = size(X,1); % 城市的数目% 通过向量化的方法计算距离矩阵% 计算矩阵中两两点之间的距离，可以通过转置相减的方法coor_x_tmp1 = repmat(X(:,1),1,amount);coor_x_tmp2 = coor_x_tmp1';coor_y_tmp1 = repmat(X(:,2),1,amount);coor_y_tmp2 = coor_y_tmp1';D = sqrt((coor_x_tmp1-coor_x_tmp2).^2 + (coor_y_tmp1-coor_y_tmp2).^2);%% 模拟退火算法确定最短路径%% 设置求解问题的参数problem.numVar = 14; %变量个数problem.fun = @(x,D)obj_fun(x,D); %优化目标函数名称problem.fun_CV = @(x)obj_fun_CV(x); %约束条件%% 模拟退火的参数SAParameters.temperature = 1000; % 初始温度 设置的足够大的话，可以在初始拥有更好的性能SAParameters.kb = 0.3; % 温度系数SAParameters.alpha = 0.99; % 降温系数SAParameters.penalty = 1.5; % 惩罚系数SAParameters.num = 1000; % 马尔可夫链长度SAParameters.Tmin = 1; % 结束温度%% 初始化解 产生一个可行解 实际上因为加入了罚函数可以考虑直接生成一个范围内的解即可while(1)variables = 1:amount; variables = variables'; % 初始化一条可行的路径CV = problem.fun_CV(variables);if(CV==0)break;endendvar_final = variables; % 初始化最终最优解T = SAParameters.temperature; % 初始化温度E0_OBJ = problem.fun(variables,D); % 初始化目标函数值E0_CV = problem.fun_CV(variables); % 初始化CV值E0 = E0_OBJ+SAParameters.penalty*E0_CV; % 最终目标值E_OBJ_f = E0; % 初始化最佳温度while (T>=SAParameters.Tmin) % 开始降温for i = 1:SAParameters.num % 马尔科夫链%% 通过二交换或三交换产生一个新解if (rand < 0.5)% 随机决定是进行两交换还是三交换% 两交换法产生新的解ind1 = 0; ind2 = 0;while (ind1 == ind2)ind1 = ceil(rand.*amount);ind2 = ceil(rand.*amount);endtmp1 = variables(ind1);variables(ind1) = variables(ind2);variables(ind2) = tmp1; % 随机交换2个变量的位置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);endtmp1 = ind1;tmp2 = ind2;tmp3 = ind3;% 确保ind1 < ind2 < ind3if (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 = variables((ind1+1):(ind2-1));variables((ind1+1):(ind1+ind3-ind2+1)) = variables((ind2):(ind3));variables((ind1+ind3-ind2+2):ind3) = tmplist1;end%% 移动后的目标值计算E_OBJ = problem.fun(variables,D); % 移动后的目标函数值E_CV = problem.fun_CV(variables); % 移动后的CV值E = E_OBJ+SAParameters.penalty*E_CV;dE = E-E0;if (E_OBJ<=E_OBJ_f && E_CV==0)var_final = variables; % 保留精英E_OBJ_f=E_OBJ;endprob=exp(-dE/SAParameters.kb/T);if(dE>0 && rand()>prob)var_final = variables; % 还原endE0_OBJ=problem.fun(variables,D); %初始目标函数值E0_CV=problem.fun_CV(variables); %初始CV值E0=E0_OBJ+SAParameters.penalty*E0_CV;endT = T*SAParameters.alpha; % 降温endfunction f = obj_fun(var,D) % 目标函数f = 0;for i=1:13f = f + D(var(i),var(i+1));endendfunction CV = obj_fun_CV(x) % 约束条件函数CV = 0;end
3. 遗传算法
遗传算法的起源可追溯到20世纪60年代初期。1967年,美国密歇根大学J. Holland教授的学生 Bagley在他的博士论文中首次提出了遗传算法这一术语,并讨论了遗传算法在博弈中的应用,但早期研究缺乏带有指导性的理论和计算工具的开拓。1975年, J. Holland等提出了对遗传算法理论研究极为重要的模式理论,出版了专著《自然系统和人工系统的适配》,在书中系统阐述了遗传算法的基本理论和方法,推动了遗传算法的发展。20世纪80年代后,遗传算法进入兴盛发展时期，被广泛应用于自动控制、生产计划、图像处理、机器人等研究领域。
本题中，由于初始解为序列形式，本质上可以视为染色体编码，因此交叉操作可以采用二交换法，变异操作采用三交换法。选择与竞标赛操作与常规遗传算法相同。当迭代次数达到1000次以上后，遗传算法求解的效果为三种方法中最佳的。
主函数
%% 遗传算法解决TSP问题clc;clear;%% 数据初始化addpath('obj_function');%将目标函数文件夹添加到搜索路径中 包含两个文件 一个存放目标函数 一个存放约束条件load CityPosition1.matcity_order = linspace(1,14,14); % 城市标号amount = size(X,1); % 城市的数目% 通过向量化的方法计算距离矩阵% 计算矩阵中两两点之间的距离，可以通过转置相减的方法coor_x_tmp1 = repmat(X(:,1),1,amount);coor_x_tmp2 = coor_x_tmp1';coor_y_tmp1 = repmat(X(:,2),1,amount);coor_y_tmp2 = coor_y_tmp1';D = sqrt((coor_x_tmp1-coor_x_tmp2).^2 + (coor_y_tmp1-coor_y_tmp2).^2);%% 设置问题的参数problem.name = 'OBJ_Fitness';problem.numVar = 14; %变量个数problem.fun = @(x,D)OBJ_Fitness(x,D); % 优化目标函数名称problem.fun_CV = @(x)OBJ_Fitness_CV(x); % 约束条件%% 遗传算法参数 结构体变量GAParameters.popsize = 40; % 种群大小GAParameters.maxgen = 50; % 最大进化代数GAParameters.pc = 0.7; % 交叉概率GAParameters.pm = 0.4; % 变异概率GAParameters.psm = 0.1; % 进行模拟退火运算的概率GAParameters.alpha = 0.9; % 模拟退火的降温系数GAParameters.temperature = 100; % 模拟退火的温度GAParameters.step = 3; % 模拟退火的基本步数 为整数%% 种群初始化P = initialize_pop(problem,GAParameters,D);%% 主程序k = 0;% P表示父种群 Q表示子种群while k < GAParameters.maxgenP = sort_pop(P); % 按照适应度和CV值排序大小对种群进行排序Q = tournament_select(P,GAParameters); % 锦标赛选择%% 遗传操作Q=simulated_binary_cross(problem,GAParameters,Q,D);%模拟二进制交叉Q=polynomal_mutate(Q,problem,GAParameters,D);%多项式变异R=[P;Q];%P和Q混合R=sort_pop(R);P = R(1:GAParameters.popsize);%产生新一代群体k = k+1;%更新迭代次数end% 最优值obj = min([P.obj]);vmin = P.var;
种群初始化函数
遗传算法基本参数设置，每个种群初始染色体编码的话就全部设置为1-14的序列。
%% 初始化种群function pop = initialize_pop(problem,GAParameters,D)individual.var = 0; % 个体的变量值individual.obj = 0; % 个体的目标值individual.CV = 0; % 个体的CV值pop = repmat(individual,GAParameters.popsize,1); % 产生种群大小的个体numVar = problem.numVar; % 变量的个数variables = linspace(1,numVar,numVar); % 变量元素初始化for k = 1:GAParameters.popsizepop(k).var = variables; %给每一个个体（染色体）赋值pop(k).obj = problem.fun(variables,D); %每一个染色体的目标函数值适应度pop(k).CV = problem.fun_CV(variables); %每一个个体的约束情况 CV值endend
交叉函数
交叉函数采用二变换，即任意交换序列中两个元素的位置。
%% 交叉操作function value_new = simulated_binary_cross(problem,GAParameters,P,D)amount = 14;P1 = P;if (rand < GAParameters.pc)% 随机决定是进行两交换还是三交换% 两交换法产生新的解for k = 1:GAParameters.popsizeind1 = 0; ind2 = 0;while (ind1 == ind2)ind1 = ceil(rand.*amount);ind2 = ceil(rand.*amount);endtmp1 = P1(k).var(ind1);P1(k).var(ind1) = P1(k).var(ind2);P1(k).var(ind2) = tmp1; % 随机交换2个变量的位置endendfor k = 1:GAParameters.popsizeif problem.fun(P1(k).var,D) < problem.fun(P(k).var,D)P(k).var = P1(k).var; P(k).obj = problem.fun(P1(k).var,D);endendvalue_new = P;end
变异函数
由于直接将访问城市顺序作为序列，二进制中的变异没有意义，因此采用三变换来代替变异的作用。
%% 变异操作function value_new = polynomal_mutate(P,problem,GAParameters,D)amount = 14;P1 = P;if (rand < GAParameters.pm)for k = 1:GAParameters.popsize 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);endtmp1 = ind1;tmp2 = ind2;tmp3 = ind3;% 确保ind1 < ind2 < ind3if (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 = P(k).var((ind1+1):(ind2-1));P1(k).var((ind1+1):(ind1+ind3-ind2+1)) = P1(k).var((ind2):(ind3));P1(k).var((ind1+ind3-ind2+2):ind3) = tmplist1;endendfor k = 1:GAParameters.popsizeif problem.fun(P1(k).var,D) < problem.fun(P(k).var,D)P(k).var = P1(k).var; P(k).obj = problem.fun(P1(k).var,D);endendvalue_new = P;end
锦标赛选择函数
函数的适应度就是距离的总值。采用竞标赛规则，两两比较大小，选出距离更小的那一个，模拟遗传中选择的过程。
function newpop = tournament_select(P,GAParameters) % 锦标赛法newpop = P;for k = 1:GAParameters.popsizecandidate = randperm(GAParameters.popsize,2); %随机产生两个序号candidate_win = min(candidate); %序号小的为胜者 因为之前已经排序过了newpop(k).var = P(candidate_win).var; %结构体数组定义直接加维度endend
排序函数
对函数适应度进行排序，由于本题没有罚函数的影响，该函数可以直接并入竞标赛选择函数，但是当存在罚函数时，需要先进行排序，然后进行竞标赛选择。
function pop=sort_pop(pop)%第二步：然后是适应度值排序[~,I]=sort([pop.obj]);pop=pop(I);%第三步：CV值排序[~,I]=sort([pop.CV]);pop=pop(I);end
总结
以上就是三种启发式算法对TSP问题的求解，最终答案徘徊在26左右。当城市数量增加时，贪婪算法便很难求出最优解了，最好还是采用遗传算法和模拟退火。考虑时间复杂度和空间复杂度，当城市数量增加到一定量级后，我认为遗传算法是最佳的寻优算法。当然，模拟退火也可以结合粒子群，本质上可以与遗传算法类似。变化很多，但内核不变，更好的寻优算法，在于找到一个更高效的梯度寻优模式。
