1、TSP 问 题 的 求 解摘 要旅 行 商 问 题 ( Traveling Salesman Problem, TSP) 代 表 一 类 组 合 优 化 问 题 ,在 计 算 机 网 络 、 公 路 交 通 分 布 等 多 种 实 际 问 题 中 都 有 重 要 意 义 。 “旅 行 商 问 题 ”也常 被 称 为 “旅 行 推 销 员 问 题 ”, 其 实 质 为 是 指 一 名 推 销 员 要 拜 访 多 个 地 点 时 , 如何 找 到 在 拜 访 每 个 地 点 一 次 后 再 回 到 起 点 的 最 短 路 径 。针 对 该 题 求 解 经 过 30 个 城 市 旅 行 的 最 短
2、路 径 问 题 , 我 们 采 用 三 种 方 法 解 决 :法 一 : ( 模 拟 退 火 算 法 ) 借 助 Matlab程 序 采 用 模 拟 退 火 算 法 进 行 分 析 解 决 。 即首 先 在 固 定 温 提 条 件 下 求 出 当 前 温 度 下 的 最 短 路 径 ( 局 部 最 优 解 ) , 然 后 改 变 温度 后 再 次 求 局 部 最 优 后 得 到 最 终 的 最 短 路 径 ( 全 局 最 优 )法 二 : ( 遗 传 算 法 ) 借 助 Matlab 程 序 采 用 遗 传 算 法 进 行 分 析 解 决 。 即法 三 : ( 线 性 规 划 ) 找 出 该
3、线 性 目 标 规 划 的 目 标 函 数 及 约 束 条 件 , 借 助 Lingo 软件 求 得 该 TSP 问 题 的 最 短 路 径关 键 词 : TSP问 题 模 拟 退 火 算 法 线 性 规 划 遗 传 算 法一 、 问 题 重 述1.1引 言TSP 是 典 型 的 组 合 优 化 问 题 , 并 且 是 一 个 NP-hard 问 题 , TSP 简 单 描 述 为 :一 名 商 人 欲 到 n个 不 同 的 城 市 去 推 销 商 品 , 每 2个 城 市 i和 j之 间 的 距 离 为 dij,如 何 选 择 一 条 路 径 使 得 商 人 每 个 城 市 走 一 遍 后
4、回 到 起 点 , 所 走 的 路 径 最 短 。 用 数学 符 号 表 示 为 :设 n 维 向 量 s =(c1 , c2 , , cn )表 示 一 条 路 经 , 目 标 函 数 为 :minC(c1,c2,cn)= 11 ),1()1,(ni cncdcicid 。1.2问 题 提 出题 目 给 出 30个 城 市 的 分 布 图 象 及 具 体 坐 标 , 要 求 出 经 过 30个 城 市 的 最 短 旅行 路 径 。 本 文 依 次 解 决 以 下 问 题 : 寻 找 目 标 函 数 ; 模 拟 退 火 算 法 中 初 始 温 度 T0的 确 定 ; 根 据 Metropoli
5、s 准 则 ( 等 温 ) 的 概 率 判 断 是 否 接 受 新 状 态 ; 通 过 循 环 先求 出 局 部 最 优 解 再 得 出 最 终 的 最 短 路 径 。二 、 问 题 分 析2.1城 市 编 号 及 坐 标为 方 便 问 题 说 明 , 我 们 先 给 30 个 城 市 编 号 并 给 出 坐 标 值 如 下 表 1 所 示 :城 市 编 号 X 坐 标 Y 坐 标 城 市 编 号 X 坐 标 Y 坐 标 城 市 编 号 X 坐 标 Y 坐 标1 41 94 11 64 60 21 87 762 37 84 12 18 54 22 18 403 54 67 13 22 60 2
6、3 13 404 25 62 14 83 46 24 82 75 7 64 15 91 38 25 62 326 2 99 16 25 38 26 58 357 68 58 17 24 42 27 45 218 71 44 18 58 69 28 41 269 54 62 19 71 71 19 44 3510 83 69 20 74 78 30 4 50表 12.2对 于 该 旅 行 中 的 最 短 路 径 问 题 , 利 用 模 拟 退 火 算 法 解 决 问 题 的 过 程 如 下 :( 外 循 环 ) ( 内 循 环 )即 温 度 作 为 外 循 环 利 用 Metropolis 抽
7、样 稳 定 准 则 求 出 当 前 温 度 下 的 最 短 路径 ( 局 部 最 优 解 ) , 然 后 改 变 温 度 后 再 次 求 局 部 最 优 后 得 到 最 终 的 最 短 路 径 ( 全局 最 优 )2.3对 于 该 旅 行 中 的 最 短 路 径 问 题 , 利 用 遗 传 算 法 解 决 问 题 的 流 程 如 下 :随 机 产 生 初 始 路 径S0, 令 Sbest=S0,并 计 算初 始 距 离 L0 设 置 初 始 温 度 Tk, 且 10,0,1 ktt kk Metropolis 抽 样 稳定 准 则 求 相 应 温度 下 最 优 路 径最 短 路 径三 、 符
8、号 说 明序 号 符 号 含 义1 dij i和 j之 间 的 距 离 为 dij2 Tk 初 始 温 度 为 Tk3 S0 初 始 路 径 为 S04 i,j 当 前 状 态 i , 新 状 态 j5 Sj,Si 当 前 状 态 路 径 Si , 新 状 态 路 径 Sj6 Li 路 径 S长 度四 、 模 型 的 建 立 与 求 解4.1模 型 基 础( 1) 模 拟 退 火 算 法模 拟 退 火 算 法 背 景 : 模 拟 退 火 算 法 来 源 于 固 体 物 质 降 温 原 理 , 在 高 温 条 件下 , 粒 子 将 自 由 运 动 , 在 达 到 稳 定 状 态 后 , 再 缓
9、慢 的 降 低 其 温 度 , 则 固 体 物 质最 终 会 形 成 最 低 能 量 的 状 态 。 可 以 将 这 种 思 想 应 用 于 求 解 组 合 优 化 问 题 。 将 组 合优 化 问 题 解 空 间 中 的 一 个 解 对 应 于 固 体 降 温 过 程 中 的 一 个 状 态 , 将 目 标 函 数 对 应于 该 状 态 下 的 能 量 。 以 控 制 参 数 T来 模 拟 固 体 的 温 度 。 对 于 每 一 个 T, 进 行 迭 代过 程 : “ 解 变 换 产 生 新 解 判 别 准 则 新 解 的 取 舍 ” , 采 用 Metropolis 准 则来 决 定 解
10、的 取 舍 , 随 着 温 度 的 降 低 , 该 算 法 有 可 能 从 局 部 极 值 区 域 跳 出 , 从 而 达到 全 局 最 优 解 。Boltzmann概 率 分 布在 温 度 T, 分 子 停 留 在 状 态 r 满 足 Boltzmann概 率 分 Ds BB BTk sETZ TZk rrEE Tk rETZrEEP )(exp)( )(Boltzmann0 )()(exp)(1)( 子 :为 概 率 分 布 的 标 准 化 因常 数 。为 的 能 量 ,表 示 状 态机 变 量 ,表 示 分 子 能 量 的 一 个 随则 在 同 一 个 温 度 , 分 子 停 留 在 能
11、 量 小 的 状 态 的 概 率 比 停 留 在 能 量 大 的 状 态 的概 率 要 大 。Metropolis 准 则 ( 等 温 ) 以 概 率 接 受 新 状 态若 在 温 度 T, 当 前 状 态 i 新 状 态 j1.若 Ej=tau%Len2-Len1Riffind(Len2-Len1R)=0)path(find(Len2-Len1R)=0), : )=path2(find(Len2-Len1R)=0),:);Len1(find(Len2-Len1R)=0)=Len2(find(Len2-Len1R)=0); TempMinD,TempIndex=min(Len1);%TempM
12、inDTracePath(N,:)=path(TempIndex,:);Distance(N,:)=TempMinD;N=N+1;%T=T*0.9m_num=0;else m_num=m_num+1;enditer_num=iter_num+1;endT=T*0.9;%m_num,iter_num,NendMinD,Index=min(Distance);BestPath=TracePath(Index,:);disp(MinD)%T1=clock%更 新 路 线 子 程 序functionp2=ChangePath2(p1,CityNum)%globalp2;while(1)R=unidrn
13、d(CityNum,1,2);ifabs(R(1)-R(2)1break;endendR=unidrnd(CityNum,1,2);I=R(1);J=R(2);%len1=D(p(I),p(J)+D(p(I+1),p(J+1);%len2=D(p(I),p(I+1)+D(p(J),p(J+1);ifI=rand %交 叉 概 率 PcSelCh(i,:),SelCh(i+1,:)=intercross(SelCh(i,:),SelCh(i+1,:);endend%输 入 :%a 和 b 为 两 个 待 交 叉 的 个 体%输 出 :%a 和 b 为 交 叉 后 得 到 的 两 个 个 体5 变
14、 异 操 作 函 数 Mutate的 代 码 :function a,b=intercross(a,b)L=length(a);r1=randsrc(1,1,1:L);r2=randsrc(1,1,1:L);if r1=r2a0=a;b0=b;s=min(r1,r2);e=max(r1,r2);for i=s:ea1=a;b1=b;a(i)=b0(i);b(i)=a0(i);x=find(a=a(i);y=find(b=b(i);i1=x(x=i);i2=y(y=i);if isempty(i1)a(i1)=a1(i);endif isempty(i2)b(i2)=b1(i);endenden
15、d6 进 化 逆 转 操 作 函 数 Reverse 代 码 :% SelCh 变 异 后 的 个 体function SelCh=Mutate(SelCh,Pm)NSel,L=size(SelCh);for i=1:NSelif Pm=randR=randperm(L);SelCh(i,R(1:2)=SelCh(i,R(2:-1:1);endend7 画 出 所 给 路 线 的 轨 迹 图 函 数 DrawPath 的 代 码 :%画 图CityPosition=4194;3784;5467;2562;764;299;6858;7144;5462;8369;6460;1854;2260;83
16、46;9138;2538;2442;5869;7171;7478;8776;1840;1340;827;6232;5835;4521;4126;4435;450;x=41,37,54,25,7,2,68,71,54,83,64,18,22,83,91,25,24,58,71,74,87,18,13,82,62,58,45,41,44,4;y=94,84,67,62,64,99,58,44,62,69,60,54,60,46,38,38,42,69,71,78,76,40,40,7,32,35,21,26,35,50;set(gcf,color,w);scatter(x,y,k,filled);
17、holdon;x1=x(30),x(12);y1=y(30),y(12);plot(x1,y1,k);holdon;x2=x(12),x(13);y2=y(12),y(13);plot(x2,y2,k);holdon;x3=x(13),x(4);y3=y(13),y(4);plot(x3,y3,k);holdon;x4=x(4),x(5);y4=y(4),y(5);plot(x4,y4,k);holdon;x5=x(5),x(6);y5=y(5),y(6);plot(x5,y5,k);holdon;x6=x(6),x(1);y6=y(6),y(1);plot(x6,y6,k);holdon;x
18、7=x(1),x(2);y7=y(1),y(2);plot(x7,y7,k);holdon;x8=x(2),x(3);y8=y(2),y(3);plot(x8,y8,k);holdon;x9=x(3),x(9);y9=y(3),y(9);plot(x9,y9,k);holdon;x10=x(9),x(18),;y10=y(9),y(18);plot(x10,y10,k);holdon;x11=x(18),x(19);y11=y(18),y(19);plot(x11,y11,k);holdon;x12=x(19),x(20);y12=y(19),y(20);plot(x12,y12,k);hol
19、don;x13=x(20),x(21);y13=y(20),y(21);plot(x13,y13,k);holdon;x14=x(21),x(10);y14=y(21),y(10);plot(x14,y14,k);holdon;x15=x(10),x(11);y15=y(10),y(11);plot(x15,y15,k);holdon;x16=x(11),x(7);y16=y(11),y(7);plot(x16,y16,k);holdon;x17=x(7),x(8);y17=y(7),y(8);plot(x17,y17,k);holdon;x18=x(8),x(14);y18=y(8),y(1
20、4);plot(x18,y18,k);holdon;x19=x(14),x(15);y19=y(14),y(15);plot(x19,y19,k);holdon;x20=x(15),x(24);y20=y(15),y(24);plot(x20,y20,k);holdon;x21=x(24),x(25);y21=y(24),y(25);plot(x21,y21,k);holdon;x22=x(25),x(26);y22=y(25),y(26);plot(x22,y22,k);holdon;x23=x(26),x(27);y23=y(26),y(27);plot(x23,y23,k);holdon
21、;x24=x(27),x(28);y24=y(27),y(28);plot(x24,y24,k);holdon;x25=x(28),x(29);y25=y(28),y(29);plot(x25,y25,k);holdon;x26=x(29),x(16);y26=y(29),y(16);plot(x26,y26,k);holdon;x27=x(16),x(17);y27=y(16),y(17);plot(x27,y27,k);holdon;x28=x(17),x(22);y28=y(17),y(22);plot(x28,y28,k);holdon;x29=x(22),x(23);y29=y(22
22、),y(23);plot(x29,y29,k);holdon;x30=x(23),x(30);y30=y(23),y(30);plot(x30,y30,k);8 %D 距 离 矩 阵%NIND 为 种 群 个 数%X 参 数 是 中 国 10 个 城 市 的 坐 标 (初 始 给 定 )%MAXGEN 为 停 止 代 数 , 遗 传 到 第 MAXGEN代 时 程 序 停 止 ,MAXGEN 的 具 体 取 值 视 问 题 的 规 模 和耗 费 的 时 间 而 定%m 为 适 值 淘 汰 加 速 指 数 ,最 好 取 为 1,2,3,4,不 宜 太 大%Pc 交 叉 概 率%Pm 变 异 概
23、率%输 出 :%R 为 最 短 路 径%Rlength 为 路 径 长 度clearclcclose allX=41 94; 37 84; 54 67; 25 6; 27 64; 2 99; 68 58; 71 44; 54 62; 83 69;64 60; 18 54; 22 60; 83 46; 91 38; 25 38; 24 42; 58 69; 71 71; 74 78;87 76; 18 40; 13 40; 82 7; 62 32; 58 35; 45 21; 41 26; 44 35; 4 50;D=Distanse(X); %生 成 距 离 矩 阵N=size(D,1); %
24、城 市 个 数%遗 传 参 数NIND=100; %种 群 大 小MAXGEN=200; %最 大 遗 传 代 数Pc=0.9; %交 叉 概 率Pm=0.05; %变 异 概 率GGAP=0.9; %代 沟% 初 始 化 种 群Chrom=InitPop(NIND,N);% 画 出 随 机 解 的 路 径 图%DrawPath(Chrom(1,:),X);%titlepause(0.0001)% 输 出 随 机 解 的 路 径 和 总 距 离disp(初 始 种 群 中 的 一 个 随 机 值 :)OutputPath(Chrom(1,:);Rlength=PathLength(D,Chro
25、m(1,:);disp(总 距 离 : ,num2str(Rlength);disp()% 优 化gen=0;figure;hold on;box onxlim(0,MAXGEN)title(优 化 过 程 )xlabel(代 数 )ylabel(最 优 值 )ObjV=PathLength(D,Chrom); %计 算 路 径 长 度preObjV=min(ObjV);while gen,num2str(R(i);enddisp(p)计 算 个 体 路 线 长 度 函 数 PathLength 代 码 :%D 两 两 城 市 之 间 的 距 离%Chrom 个 体 的 轨 迹function
26、len=PathLength(D,Chrom)row,col=size(D);NIND=size(Chrom,1);len=zeros(NIND,1);fori=1:NINDp=Chrom(i,:)Chrom(i,1);i1=p(1:end-1);i2=p(2:end);len(i,1)=sum(D(i1-1)*col+i2);end重 插 入 子 代 得 到 新 种 群 的 函 数 Reins代 码 :%Chrom 父 代 的 种 群%SelCh 子 代 种 群%ObjV 父 代 适 应 度%输 出%Chrom 组 合 父 代 与 子 代 后 得 到 的 新 种 群functionChrom
27、=Reins(Chrom,SelCh,ObjV)NIND=size(Chrom,1);NSel=size(SelCh,1);TobjV,index=sort(ObjV);Chrom=Chrom(index(1:NIND-NSel),:);SelCh;( 3) Lingo程 序MODEL:SETS:city/130/:U;!U(i)=cicyNo;links(city,city):distance,!thedistanceofthematrix;x;!x(i,j)=1ifweuselinki,j;ENDSETS!distancematrix;!10000,thedistance ofimposs
28、iblelink;data:distance=file(distance.txt);ENDDATAn=size(city);!minimizetotaldistanceofthelinks;MIN=SUM(links:distance*x);!theentrance;FOR(city(k):SUM(city(i)|i#ne#k:x(i,k)=1;!itmustbedeparted;SUM(city(j)|j#ne#k:x(k,j)=1;);FOR(city(a);FOR(city(j)|a#gt#1#and#j#ne#k:U(a)-U(j)+n*x(a,j)=n-1););FOR(links:BIN(x);!itmakesthexs0or1;FOR(city(a)|k#gt#1:U(a)=n-1);END