目录
Preface
一.极值问题Extremum problem
1.参数初始化 Parameter initialization
2.解码 decode
3.轮盘赌,选择淘汰Roulette, choose to eliminate
4.交叉 cross
5.变异 variation
6.完成迭代 Iteration completed
二.旅行商问题TSP
Reference article
Preface
关于遗传算法的知识点,可以参考其他博主的文章,写的都很详尽易懂。作为数模小白,作者在这里分享关于代码的理解和每一步操作的意义。
一.极值问题Extremum problem
1.参数初始化 Parameter initialization
function result = func1(x) %待求的函数
fit = x+10*sin(5*x)+7*cos(4*x);
result = fit;
end
NP = 100; %种群数量
L = 20; %二进制位串长度
Pc = 0.8; %交叉率
Pm = 0.01; %变异率
G = 100; %最大遗传代数
Xs = 10; %上限
Xx = 0; %下限
f = rand(NP,L); %随机获得初始种群
trace=zeros(1,G); %用来记录每次迭代的最优解
2.解码 decode
for k = 1:G %这是外层循环,对应的end在后面
for i = 1:NP
U = f(i,:); %第i行,解码
m = 0;
for j = 1:L
m = U(j)*2^(j-1)+m;%二进制转化为10进制
end
x(i) = Xx+m*(Xs-Xx)/(2^L-1);%代回定义域
Fit(i) = func1(x(i));
end
maxFit = max(Fit); %最大值,寻求的最优个体
minFit = min(Fit); %最小值
trace(k) = maxFit; %历代最优适应度,记录下来
rr = find(Fit==maxFit);
fBest = f(rr(1,1),:); %历代最优个体
xBest = x(rr(1,1));
Fit = (Fit-minFit)/(maxFit-minFit); %归一化适应度值
3.轮盘赌,选择淘汰Roulette, choose to eliminate
sum_Fit = sum(Fit);
fitvalue = Fit./sum_Fit; %返回概率矩阵
fitvalue = cumsum(fitvalue);%累计和,a(j)=a(1)+a(2)+...a(j)
ms = sort(rand(NP,1));%轮盘赌
fiti = 1;
newi = 1;
while (newi <=NP) && (fiti<=NP)
if (ms(newi)) < fitvalue(fiti)
nf(newi,:) = f(fiti,:);%f是初始种群,储存被选中的个体
newi = newi+1; %对于优秀个体(占比率较多的)多复制的可能更大
else
fiti = fiti+1;
end
end
4.交叉 cross
for i = 1:2:NP
p = rand;
if p < Pc %进行交叉
q = rand(1,L);
for j = 1:L
if q(j)<Pc
temp = nf(i+1,j);
nf(i+1,j) = nf(i,j);
nf(i,j) = temp;
end
end
end
end
5.变异 variation
i = 1;
while i <= round(NP*Pm)
h = randi([1,NP],1,1); %随机选取一个需要变异的染色体,返回1*1介于[1,NP]的伪整数矩阵
for j = 1:round(L*Pm)
g = randi([1,L],1,1); %随机选取需要变异的基因数
nf(h,g) =~ nf(h,g);
end
i = i+1;
end
6.完成迭代 Iteration completed
f = nf;
f(1,:) = fBest; %保留最优个体在新种群中
end %与前面的for对应
x=linspace(Xx,Xs,10000);
plot(x,func1(x));legend('x+10*sin(5*x)+7*cos(4*x)','location','northwest')
title('function');
figure
plot(trace)
xlabel('Iterations')
ylabel('Objective function value')
title('Fitness evolution curve')
从上面可以发现,GA很快就收敛到了目标值,没有陷入局部最优值,效率还是很快的
二.旅行商问题TSP
关于代码的解释上面已给出,下面分享TSP问题的代码,稍难的已给出注释。另可参考阿白数模笔记之蚁群算法(1):TSP问题,遍历全球244个国家和地区的最短路径和
function [out]=distance(d)
m=size(d,1);
out=zeros(m,m);
for i=1:m
for j=i:m
a=d(i,:);
b=d(j,:);
out(i,j)= ((a-b)*(a-b)')^0.5; % 计算最短直线距离
out(j,i)=out(i,j);
end
end
end
tic
clc;clear
load('data2.mat')%C 31个地区的坐标
N = size(C,1); %TSP 问题的规模,即城市数目
D =distance(C); %任意两个城市距离间隔矩阵
NP = 200; %种群规模
G = 5000; %最大遗传代数
% f = zeros(NP,N); %用于存储种群
F = []; %种群更新中间存储
for i = 1:NP
f(i,:) = randperm(N); %p = randperm(n) 返回行向量,其中包含从 1 到 n 没有重复元素的整数随机排列。
end
R = f(1,:); %存储最优种群
len = zeros(NP,1); %存储路径长度
fitness = zeros(NP,1); %存储归一化适应值
gen = 0;
while gen < G
%%%%%%%%%%%%%%%计算路径长度%%%%%%%%%%%%%%%%
for i = 1:NP
len(i,1) = D(f(i,N),f(i,1));
for j = 1:(N-1)
len(i,1) = len(i,1)+D(f(i,j),f(i,j+1));
end
end
maxlen = max(len); %最长路径
minlen = min(len); %最短路径
%%%%%%%%%%%%%%%更新最短路径%%%%%%%%%%%%%%%
rr = find(len==minlen);
R = f(rr(1,1),:);
%%%%%%%%%%%%%%计算归一化适应值%%%%%%%%%%%%%%
for i = 1:length(len)
fitness(i,1) = (1-((len(i,1)-minlen)/(maxlen-minlen+0.001)));
end
%%%%%%%%%%%%%%%%%选择操作%%%%%%%%%%%%%%%%
nn = 0;
for i = 1:NP
if fitness(i,1) >= rand
nn = nn+1;
F(nn,:) = f(i,:);
end
end
[aa,bb] = size(F);
while aa < NP
nnper = randperm(nn);
A = F(nnper(1),:);
B = F(nnper(2),:);
%%%%%%%%%%%%%%%交叉操作%%%%%%%%%%%%%%%
W = ceil(N/10); %交叉点个数
p = unidrnd(N-W+1); %随机选择交叉范围,从 p 到 p+W,
%r = unidrnd(n) 从由最大值 n 指定的离散均匀分布中生成随机数。
for i = 1:W
x = find(A==B(1,p+i-1));
y = find(B==A(1,p+i-1));
temp = A(1,p+i-1);
A(1,p+i-1) = B(1,p+i-1);
B(1,p+i-1) = temp;
temp = A(1,x);
A(1,x) = B(1,y);
B(1,y) = temp;
end
%%%%%%%%%%%%%%%%%变异操作%%%%%%%%%%%%%
%但迭代到后面为防止重要基因丢失,变异概率应随着迭代进行而降低
for k=1:floor((G-gen)/G*0.2*N)+1
p1 = floor(1+N*rand());
%Y = floor(X) 将 X 的每个元素四舍五入到小于或等于该元素的最接近整数。
p2 = floor(1+N*rand());
while p1==p2
p1 = floor(1+N*rand());
p2 = floor(1+N*rand());
end
tmp = A(p1);
A(p1) = A(p2);
A(p2) = tmp;
tmp = B(p1);
B(p1) = B(p2);
B(p2) = tmp;
F = [F;A;B];
[aa,bb] = size(F);
end
end
if aa > NP
F = F(1:NP,:); %保持种群规模为 NP
end
f = F; %更新种群
f(1,:) = R; %保留每代最优个体
clear F;
gen = gen+1;
Rlength(gen) = minlen;
for i = 1:N-1
plot([C(R(i),1),C(R(i+1),1)],[C(R(i),2),C(R(i+1),2)],'bo-');
hold on
end
plot([C(R(N),1),C(R(1),1)],[C(R(N),2),C(R(1),2)],'ro-');
title(['Minimum distance:',num2str(minlen)]);
hold off
end
figure
plot(Rlength)
xlabel('Iterations')
ylabel('Objective function value')
title('Fitness evolution curve')
toc
历时 233.815222 秒。
Reference article
[1]遗传算法小结及算法实例(附Matlab代码)