最近有一个作业是用遗传算法求解旅行商问题,在网上找了博主的代码,但结果方面还是存在不小的误差,希望有大神能解答一下
主函数:
clc
clear
load("berlin52.mat"); % 加在数据
location = load("berlin52.mat"); % 开始对数据进行可视化
x = berlin52(:,1);
y = berlin52(:,2);
plot(x,y,'bo','MarkerSize', 4);
xlabel('城市横坐标');
ylabel('城市纵坐标');
grid on;
NIND = 300; %种群大小
MAXGEN = 1200; %最大迭代次数
Pc = 0.9; %交叉概率,相当于基因遗传的时候染色体交叉
Pm = 0.05; %染色体变异
GGAP = 0.9; %这个是代沟,通过遗传方式得到的子代数为父代数*GGAP
D = Distance(berlin52); %通过这个函数可以计算i,j两点之间的距离
N = size(D,1); %计算有多少个坐标点
Chrom = InitPop(NIND,N); %%初始化种群
gen = 0;
trace = zeros(1,MAXGEN);
ObjV = PathLength(D,Chrom);%计算当前路线长度,即上面随机产生的那些个体路径
preObjV = min(ObjV);%当前最优解
while gen<MAXGEN
%%计算适应度
ObjV = PathLength(D,Chrom); %计算路线长度
pause(0.0001);
preObjV = min(ObjV);
trace(1,gen+1)=preObjV;
FitnV = Fitness(ObjV);
%选择
SelCh = Select(Chrom,FitnV,GGAP);
SelCh = Reverse(SelCh,D);
%交叉操作
SelCh = Recombin(SelCh,Pc);
SelCh = Reverse(SelCh,D);
%变异
SelCh = Mutate(SelCh,D,Pm);
Chrom = Reins(Chrom,SelCh,ObjV);
gen = gen + 1;
SelCh = Reverse(SelCh,D);
end
%画出最优解的路线图
ObjV = PathLength(D,Chrom); %计算路线长度
[minObjV,minInd] = min(ObjV);
DrawPath(Chrom(minInd(1),:),berlin52)
figure();
plot([1:MAXGEN],trace(1,:));
title('优化过程');
xlabel('迭代次数');
ylabel('当前最优解');
%%输出最优解的距离
disp(['旅行商走过的总距离:',num2str(ObjV(minInd(1)))]);
disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
Distance函数:
function D = Distance(x)
row = size(x,1);
D = zeros(row,row);
for i = 1:row
for j = i+1:row
D(i,j) = ((x(i,1)-x(j,1))^2+(x(i,2)-x(j,2))^2)^0.5;
D(j,i) = D(i,j);
end
end
DrawPath函数:
function DrawPath(Chrom,X)
R = [Chrom(1,:) Chrom(1,1)]; %第一个随机解(个体)【Chrom(1,:)取第一行数据】,一共有n个城市,但是这里R有n+1个值,因为后面又补了一个Chrom(1,1),“是为了让路径最后再回到起点”
figure;
hold on
plot(X(:,1),X(:,2),'bo','MarkerSize', 4);%X(:,1),X(:,2)分别代表的X轴坐标和Y轴坐标
plot(X(Chrom(1,1),1),X(Chrom(1,1),2),'rv','MarkerSize',9)%标记起点
A = X(R,:); %A是将之前的坐标顺序用R打乱后重新存入A中
row = size(A,1); %row为坐标数+1
for i = 2:row
% 直接使用 plot 函数连接两个点
plot(A(i-1:i,1), A(i-1:i,2), 'g-', 'LineWidth', 1);
end
hold off
xlabel('横坐标x')
ylabel('纵坐标y')
title('旅行商轨迹图')
box on
初始化函数:
function Chrom = InitPop(NIND,N)
Chrom = zeros(NIND,N); % 定义存储种群的变量
for i = 1:NIND
Chrom(i,:) = randperm(N);
end
Fitness函数
function FitnV = Fitness(len) %适应度函数
%输入:
%len 个体的长度(TSP的距离)
%输出:
%FitnV 个体的适应度值
FitnV = 1./len;
Intercross:
function [a,b] = intercross(a,b)
%输入:
%a和b为两个待交叉的个体
%输出:
%a和b为交叉后得到的两个个体
L = length(a);
%随机产生交叉区段
r1 = randsrc(1,1,[1:L]); % 这里先随机选出两个位置,位置不能超过
r2 = randsrc(1,1,[1:L]); %
if r1~=r2
a0 = a;
b0 = b;
s = min([r1,r2]);
e = max([r1,r2]);
for i = s:e
a1 = a;
b1 = b;
%第一次互换
a(i) = b0(i);
b(i) = a0(i);
%寻找相同的城市
x = find(a==a(i)); % 如果有相同的城市,x会得到包含i的两个值,y同理
y = find(b==b(i));
%第二次互换产生新的解
i1 = x(x~=i); % 将位置不是i但相同的城市标记出来
i2 = y(y~=i);
if ~isempty(i1)
a(i1)=a1(i);
end
if ~isempty(i2)
b(i2)=b1(i); % 注意,原文这里有误,应该是b(i2)
end
end
end
% 这里增加一段代码,r1=r2时,两个个体只在一点交叉
if r1 == r2
a0 = a;
b0 = b;
a(r1) = b0(r1);
b(r1) = a0(r1);
x = find(a==a(r1));
y = find(b==b(r1));
i1 = x(x~=r1);
i2 = y(y~=r1);
if ~isempty(i1)
a(i1) = a0(r1);
end
if ~isempty(i2)
b(i2) = b0(r1);
end
end
mutate
function SelCh = Mutate(SelCh,D,Pm)
%变异操作
%SelCh 被选择的个体
%Pm 变异概率
%输出:
%SelCh 变异后的个体
index = SelCh;
col = size(SelCh,2); %返回SelCh的列数
lenSelCh = PathLength(D,SelCh);
[NSel,L] = size(SelCh);
for i = 1:NSel
if Pm >= rand
R = randperm(L);
index(i,R(1:2)) = index(i,R(2:-1:1)); % 将个体i中R(1)和R(2)这两个位置的城市互换
p = [index(i,:) index(i,1)];
i1 = p(1:end-1);
i2 = p(2:end);
DIndexi = sum(D((i1-1)*col+i2)); % 计算出变异后个体的路径距离
if DIndexi < lenSelCh(i) % 如果变异后路径距离比变异前更小,则保留变异
SelCh(i,:) = index(i,:);
end
end
end
Pathlength
function len = PathLength(D,Chrom)
%计算所有个体的路线长度
%输入
%D两两城市之间的距离
%Chrom个体的轨迹
[~,col] = size(D); %返回D的列数
NIND = size(Chrom,1);%NIND等于初始种群
len = zeros(NIND,1);%初始化一个大小等于NIND的len来记录长度
for i = 1:NIND
p = [Chrom(i,:) Chrom(i,1)];%构造p矩阵保存路线图 将第一行路线提出 再加上第一个构成回路
i1 = p(1:end-1);%i1从第一个开始遍历到倒数第二个
i2 = p(2:end);%i2从第二个开始遍历到倒数第一个
len(i,1) = sum(D((i1-1)*col+i2));%计算出每种路线(种群的个体)的长度,这里相当于把D矩阵沿行展开
end
recombin:
function SelCh = Recombin(SelCh,Pc)
%交叉操作
%输入:
%SelCh 被选择的个体
%Pc 交叉概率
%输出:
%SelCh 交叉后的个体
NSel = size(SelCh,1);
for i = 1:2:NSel - mod(NSel,2) % 若不减去余数且NSel如果是奇数,最后一个i会等于NSel,i+1报错
if Pc>=rand %交叉概率PC
[SelCh(i,:),SelCh(i+1,:)] = intercross(SelCh(i,:),SelCh(i+1,:));
end
end
Reins
function Chrom = Reins(Chrom,SelCh,ObjV)
%重插入子代的种群
%输入:
%Chrom 父代的种群
%SelCh 子代的种群
%ObjV 父代适应度
%输出:
%Chrom 组合父代与子代后得到的新种群
NIND = size(Chrom,1);
NSel = size(SelCh,1);
[TobjV,index] = sort(ObjV);
Chrom = [Chrom(index(1:NIND-NSel),:);SelCh];
Reverse:
function SelCh = Reverse(SelCh,D)
%进化逆转函数
%输入:
%SelCh 被选择的个体
%D 各城市的距离矩阵
%输出:
%SelCh 进化逆转后的个体
[row,col] = size(SelCh);
ObjV = PathLength(D,SelCh);
SelCh1 = SelCh;
for i = 1:row
r1 = randsrc(1,1,[1:col]);
r2 = randsrc(1,1,[1:col]);
mininverse = min([r1 r2]);
maxinverse = max([r1 r2]);
SelCh1(i,mininverse:maxinverse) = SelCh1(i,maxinverse:-1:mininverse);
end
ObjV1 = PathLength(D,SelCh1);%计算路线长度
index = ObjV1<ObjV;
SelCh(index,:)=SelCh1(index,:);
select:
function SelCh = Select(Chrom,FitnV,GGAP)
%输入:
%Chrom 种群
%FitnV 适应度值
%GGAP 选择概率
%输出:
%SelCh 被选择的个体
NIND = size(Chrom,1);%种群个体总数
NSel = max(floor(NIND * GGAP+.5),2);%确定下一代种群的个体数,如果不足二个就计为二个
ChrIx = Sus(FitnV,NSel);
SelCh = Chrom(ChrIx,:);
Sus
function NewChrIx = Sus(FitnV,Nsel)
%输入:
%FitnV 个体的是适应度值
%Nsel 被选个体的数目
%输出:
%NewChrIx 被选择个体的索引号
[Nind,ans] = size(FitnV);%Nind为FitnV的行数也就是个体数 ans为列数1
cumfit = cumsum(FitnV);%对适应度累加 例如 1 2 3 累加后 1 3 6 用来计算累积概率
trials = cumfit(Nind)/Nsel * (rand + (0:Nsel-1)');%cumfit(Nind)表示的是矩阵cumfit的第Nind个元素 A.'是一般转置,A'是共轭转置 rand返回一个在区间 (0,1) 内均匀分布的随机数
%cumfit(Nind)/Nsel平均适应度 * (rand +(0:Nsel-1)')从0开始到89的转置矩阵(行矩阵变列矩阵)加上每一项加上一个0-1的随机数
%生成随机数矩阵 用来筛选
Mf = cumfit(:,ones(1,Nsel));%将生成的累积概率
Mt = trials(:,ones(1,Nind))';
[NewChrIx,ans] = find(Mt<Mf & [zeros(1,Nsel);Mf(1:Nind-1,:)]<= Mt);%寻找满足条件的个体 返回返回数组 X 中每个元素的行和列下标
[ans,shuf] = sort(rand(Nsel,1));%生成Nsel*1的随机数矩阵 按升序对 A 的元素进行排序 返回选择的shuf 随机打乱个体顺序 保证后续遗传算子操作的随机性
NewChrIx = NewChrIx(shuf);%返回shuf索引的矩阵元素
数据集为berlin52,
NAME: berlin52
TYPE: TSP
COMMENT: 52 locations in Berlin (Groetschel)
DIMENSION: 52
EDGE_WEIGHT_TYPE: EUC_2D
NODE_COORD_SECTION
1 565.0 575.0
2 25.0 185.0
3 345.0 750.0
4 945.0 685.0
5 845.0 655.0
6 880.0 660.0
7 25.0 230.0
8 525.0 1000.0
9 580.0 1175.0
10 650.0 1130.0
11 1605.0 620.0
12 1220.0 580.0
13 1465.0 200.0
14 1530.0 5.0
15 845.0 680.0
16 725.0 370.0
17 145.0 665.0
18 415.0 635.0
19 510.0 875.0
20 560.0 365.0
21 300.0 465.0
22 520.0 585.0
23 480.0 415.0
24 835.0 625.0
25 975.0 580.0
26 1215.0 245.0
27 1320.0 315.0
28 1250.0 400.0
29 660.0 180.0
30 410.0 250.0
31 420.0 555.0
32 575.0 665.0
33 1150.0 1160.0
34 700.0 580.0
35 685.0 595.0
36 685.0 610.0
37 770.0 610.0
38 795.0 645.0
39 720.0 635.0
40 760.0 650.0
41 475.0 960.0
42 95.0 260.0
43 875.0 920.0
44 700.0 500.0
45 555.0 815.0
46 830.0 485.0
47 1170.0 65.0
48 830.0 610.0
49 605.0 625.0
50 595.0 360.0
51 1340.0 725.0
52 1740.0 245.0
EOF