Navigator
TSP问题-GA求解
TSP问题已经被证明是NP-hard
问题,近似求解TSP问题是算法设计的一个重要研究课题,对该问题的求解算法主要分为两类:一种是与问题特征相关的启发式搜索算法,主要包括动态规划法,分支定界法等,另一类是独立于问题的智能优化算法,如模拟退火法,禁忌搜索算法,蚁群算法,遗传算法和粒子群算法等。
选择算子(selection
)
从群体中选择优秀的个体,淘汰劣质的个体的操作叫做选择,选择的目标是将优化的个体直接遗传到下一代或通过配对交叉产生新的个体再遗传到下一代。目前常用的选择算子有适应度比例方法,随机遍历抽样法和局部选择法。其中轮盘赌(Roulette Wheel Selection
) 是最简单也是最常用的选择方法,在该方法中,各个个体的选择概率和其适应度成比例,设群体大小为 n n n,其中个体 i i i的适应度为 f i f_i fi,则 i i i被选择的概率为:
P i = f i ∑ j = 1 n f j P_i=\frac{f_i}{\sum_{j=1}^n f_j} Pi=∑j=1nfjfi
概率反映了个体i
的适应度在整个群体的个体适应度总和中所占的比例,个体适应度越大,其被选择的概率越高,反之亦然。计算出群体中各个个体的选择概率后,为了选择交配个体,需要进行多轮选择。每一轮产生一个 [ 0 , 1 ] [0, 1] [0,1]之间的均匀随机数,将该随机数作为选择指针来确定被选择个体,当个体被选择后,可以随机组成交配对,以供后面的交叉操作。
交叉算子(crossover
)
交叉就是将两个父代个体部分结构加以替换重组而生成新个体的操作,通过交叉,可以提高GA的搜索能力。根据编码表示方法的不同,分为实值重组和二进制交叉两类算法:
实值重组(real valued recombination
):离散重组,中间重组,线性重组和扩展线性重组
二进制交叉(binary valued crossover
):单点交叉,多点交叉,均匀交叉,洗牌交叉和缩小代理交叉
其中最常用的交叉算子是单点交叉(one-point crossover
),在个体串中随机设定一个交叉点,实行交叉时,该点前后两部分结构进行互换,并生成两个新个体。
变异算子(mutation
)
变异算子的基本内容是对群体中个体串的某些基因座上的基因值作变动,根据个体编码表示方法的不同,分为实值变异和二进制变异,一般变异算子操作分两步完成:
- 对群体中所有个体以事先设置的变异概率判断是否进行变异
- 对进行变异的个体随机变异位进行变异
变异算子的引入主要目的有两个:
- 遗传算法具有局部随机搜索能力,当GA通过交叉算子已经接近最优解邻域时,利用变异算子这种局部随机搜索能力可以加速向最优解收敛,在这种情况下需要设置较小的变异概率,防止破坏接近最优解状态下的基因块。
- GA为了维持群体多样性,防止出现未成熟收敛现象,此时收敛概率需要取较大值。
在GA中,交叉算子为主要算子,变异算子为辅助算子,通过两种算子的相互配合与竞争操作,使得GA同时具有兼顾全局和局部的均衡搜索能力。
TSP问题建模
给定 n n n个城市和各城市之间的距离,寻找一条遍历所有城市且每个城市只被访问过一次的路径,并保证总路径距离最短。设 G = ( V , E ) G=(V, E) G=(V,E)为赋权图, V = { 1 , 2 , … , n } V=\{1, 2, \dots, n\} V={
1,2,…,n}为顶点集, E E E为边集,各顶点之间距离为 C i j ( C i j > 0 , i , j ∈ V ) C_{ij} (C_{ij}>0, i,j \in V) Cij(Cij>0,i,j∈V)并设置指示变量
x i j = { 1 On the optimal path 0 O.W. x_{ij}= \begin{cases} 1 & \text{On the optimal path}\\ 0 & \text{O.W.} \end{cases} xij={
10On the optimal pathO.W.
那么整个TSP问题的数学模型表示如下:
min Z = ∑ i ≠ j c i j x i j s . t . { ∑ j ≠ i x i j = 1 j ∈ v ∑ i ≠ j x j i = 1 j ∈ v ∑ i , j ∈ s x i j ≤ ∣ K ∣ − 1 k ⊂ v x i j ∈ { 0 , 1 } i , j ∈ v \min Z=\sum_{i\neq j} c_{ij}x_{ij}\\ s.t. \begin{cases} \sum_{j\neq i} x_{ij} = 1 & j\in v\\ \sum_{i\neq j} x_{ji} = 1 & j\in v\\ \sum_{i, j\in s}x_{ij}\leq |K|-1 & k\subset v\\ x_{ij}\in \{0, 1\} & i,j \in v \end{cases} minZ=i=j∑cijxijs.t.⎩⎪⎪⎪⎨⎪⎪⎪⎧∑j=ixij=1∑i=jxji=1∑i,j∈sxij≤∣K∣−1xij∈{
0,1}j∈vj∈vk⊂vi,j∈v
Algorithmic Framework
- 种群初始化,在TSP问题中,使用实数编码,为 1 ∼ n 1\sim n 1∼n的实数排列,初始化的参数:种群个数 M M M,染色体基因个数 N N N(表示城市的个数),迭代次数 C C C,交叉概率 P c Pc Pc和变异概率 P m Pm Pm.
- 适应度函数,将一个随机全排列的总距离的倒数作为适应度函数,即距离越短,适应度函数越好
clc;
clear all;
close all;
n = 50; % number of cities
xy = 10*rand(n, 2); % coordination of cities
popSize = 60; % number of population
numIter = 1e4; % number of iterations
showProg = 1;
showResult = 1;
a = meshgrid(1:n);
dmat = reshape(sqrt(sum((xy(a, :)-xy(a',:)).^2, 2)), n, n);
[optRoute, minDist] = tsp_ga(xy, dmat, popSize, numIter, showProg, showResult);
TSP-GA求解
function varargout = tsp_ga(xy, dmat, popSize, numIter, showProg, showResult)
nargs = 6;
for k=nargin: nargs-1
switch k
case 0
xy = 10*rand(50, 2);
case 1
N = size(xy, 1);
a = meshgrid(1:N);
dmat = reshape(sqrt(sum((xy(a, :)-xy(a', :)).^2, 2)), N, N);
case 2
popSize = 100;
case 3
numIter = 1e4;
case 4
showProg = 1;
case 5
showResult =1;
otherwise
end
end
% verify inputs
[N, dims] = size(xy);
[nr, nc] = size(dmat);
if N~=nr || N~=nc
error('Invalid XY or DMAT inputs!');
end
n = N;
% sanity checks
popSize = 4*ceil(popSize/4);
numIter = max(1, round(real(numIter(1))));
showProg = logical(showProg(1));
showResult = logical(showResult(1));
% initialize the population (real-coded)
pop = zeros(popSize, n);
pop(1, :) = (1:n);
for k = 2:popSize
pop(k, :) = randperm(n);
end
% RUN GA
globalMin = Inf;
totalDist = zeros(1, popSize);
distHistory = zeros(1, numIter);
tmpPop = zeros(4, n);
newPop = zeros(popSize, n);
if showProg
pfig = figure('Name', 'TSP_GA|Current Best Solution', 'Numbertitle', 'off');
end
for iter = 1:numIter
% Evaluate each population member (calculate total distance)
for p=1:popSize
d = dmat(pop(p, n), pop(p, 1));
for k=2:n
d = d+dmat(pop(p, k-1), pop(p, k)); % crossover, mutation
end
totalDist(p) = d;
end
% Find the best route in the population
[minDist, index] = min(totalDist);
distHistory(iter) = minDist;
% sub
if minDist < globalMin
globalMin = minDist;
optRoute = pop(index, :);
% Plot the best route
if showProg
figure(pfig);
rte = optRoute([1:n 1]);
if dims>2
plot3(xy(rte, 1), xy(rte, 2), xy(rte, 3), 'r.-');
else
plot(xy(rte, 1), xy(rte, 2), 'r.-');
end
title(sprintf('Total Distance=%1.4f, Iteration=%d', minDist, iter));
end
end
% GA
randomOrder = randperm(popSize);
for p = 4:4:popSize
rtes = pop(randomOrder(p-3: p), :);
dists = totalDist(randomOrder(p-3: p));
[~, idx] = min(dists);
bestOf4Route = rtes(idx, :);
routeInsertionPoints = sort(ceil(n*rand(1, 2)));
I = routeInsertionPoints(1);
J = routeInsertionPoints(2);
for k = 1:4
tmpPop(k, :) = bestOf4Route;
switch k
case 2 % Flip
tmpPop(k, I:J) = tmpPop(k, J:-1:I);
case 3 % Swap
tmpPop(k, [I J]) = tmpPop(k, [J I]);
case 4 % Slide
tmpPop(k, I:J) = tmpPop(k, [I+1:J I]);
otherwise
end
end
newPop(p-3:p, :) = tmpPop;
end
pop = newPop;
end
% figure
if showResult
% Plot GA Results
figure('Name', 'TSP_GA | Results', 'Numbertitle', 'off');
subplot(2, 2, 1);
% figure(1);
pclr = ~get(0, 'DefaultAxesColor');
if dims > 2
plot3(xy(:, 1), xy(:, 2), xy(:, 3), '.', 'Color', pclr);
else
plot(xy(:, 1), xy(:, 2), '.', 'Color', pclr);
end
title('Position of Cities');
grid on;
subplot(2, 2, 2);
imagesc(dmat(optRoute, optRoute));
title('Distance Matrix---imagesc');
subplot(2, 2, 3);
rte = optRoute([1:n 1]);
if dims > 2
plot3(xy(rte, 1), xy(rte, 2), xy(rte, 3), 'r.-');
else
plot(xy(rte, 1), xy(rte, 2), 'r.-');
end
title(sprintf('min distance = 1.4%f', minDist));
grid on;
subplot(2, 2, 4);
plot(distHistory, 'b', 'LineWidth', 2);
title('The best fitness curve');
grid on
set(gca, 'XLim', [0 numIter+1], 'YLim', [0 1.1*max([1 distHistory])]);
end
% Return Outputs
if nargout
varargout{1} = optRoute;
varargout{2} = minDist;
end
end
%% plot
figure;
hold on;
for i=1:n
plot(xy(:, 1), xy(:, 2), 'k.');
text(xy(i, 1), xy(i, 2)+0.08, num2str(i));
end
for i=1:n-1
plot([xy(optRoute(1, i), 1), xy(optRoute(1, i+1), 1)], [xy(optRoute(1, i), 2), xy(optRoute(1, i+1), 2)]);
end
for i=n
plot([xy(optRoute(1, i), 1), xy(optRoute(1, 1), 1)], [xy(optRoute(1, i), 2), xy(optRoute(1, 1), 2)]);
end
hold off;
Reference
MATLAB优化算法案例分析与应用 余胜威