智能优化算法---遗传算法(GA)

说明

基础内容源自 智能优化算法及其MATLAB实例(第二版)源程序

源程序为 matlab 代码,笔者在此基础上将 matlab 代码转化为 python 代码实现

一、标准遗传算法求函数极值

GA21.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
%%%%%%%%%%%%%%%%%%%%标准遗传算法求函数极值%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%初始化参数%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; %清除所有变量
close all; %清图
clc; %清屏
NP=50; %种群数量
L=20; %二进制数串长度
Pc=0.8; %交叉率
Pm=0.1; %变异率
G=100; %最大遗传代数
Xs=10; %上限
Xx=0; %下限
f=randint(NP,L); %随机获得初始种群
% 现在比较新的 matlab 版本中,需要将 randint 函数换成 randi 函数
%%%%%%%%%%%%%%%%%%%%%%%%%遗传算法循环%%%%%%%%%%%%%%%%%%%%%%%%
for k=1:G
%%%%%%%%%%%%将二进制解码为定义域范围内十进制%%%%%%%%%%%%%%
for i=1:NP
U=f(i,:);
m=0;
for j=1:L
m=U(j)*2^(j-1)+m;
end
x(i)=Xx+m*(Xs-Xx)/(2^L-1);
Fit(i)= func1(x(i));
end
maxFit=max(Fit); %最大值
minFit=min(Fit); %最小值
rr=find(Fit==maxFit);
fBest=f(rr(1,1),:); %历代最优个体
xBest=x(rr(1,1));
Fit=(Fit-minFit)/(maxFit-minFit); %归一化适应度值
%%%%%%%%%%%%%%%%%%基于轮盘赌的复制操作%%%%%%%%%%%%%%%%%%%
sum_Fit=sum(Fit);
fitvalue=Fit./sum_Fit;
fitvalue=cumsum(fitvalue);
ms=sort(rand(NP,1));
fiti=1;
newi=1;
while newi<=NP
if (ms(newi))<fitvalue(fiti)
nf(newi,:)=f(fiti,:);
newi=newi+1;
else
fiti=fiti+1;
end
end
%%%%%%%%%%%%%%%%%%%%%%基于概率的交叉操作%%%%%%%%%%%%%%%%%%
for i=1:2:NP
p=rand;
if p<Pc
q=randint(1,L);
for j=1:L
if q(j)==1;
temp=nf(i+1,j);
nf(i+1,j)=nf(i,j);
nf(i,j)=temp;
end
end
end
end
%%%%%%%%%%%%%%%%%%%基于概率的变异操作%%%%%%%%%%%%%%%%%%%%%%%
i=1;
while i<=round(NP*Pm)
h=randint(1,1,[1,NP]); %随机选取一个需要变异的染色体
for j=1:round(L*Pm)
g=randint(1,1,[1,L]); %随机需要变异的基因数
nf(h,g)=~nf(h,g);
end
i=i+1;
end
f=nf;
f(1,:)=fBest; %保留最优个体在新种群中
trace(k)=maxFit; %历代最优适应度
end
xBest; %最优个体
figure
plot(trace)
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')

func1.m

1
2
3
4
%%%%%%%%%%%%%%%%%%%%%%%%%适应度函数%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function result=func1(x)
fit= x+10*sin(5*x)+7*cos(4*x);
result=fit;

GA21.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
import numpy as np
import matplotlib.pyplot as plt
import math

# 用来正常显示中文标签及负号
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 定义目标函数 func1
def func1(x):
return x+10*sin(5*x)+7*cos(4*x) # func1.m 对应的适应度函数

# 初始化参数
NP = 50 # 种群数量
L = 20 # 二进制数串长度
Pc = 0.8 # 交叉率
Pm = 0.1 # 变异率
G = 100 # 最大遗传代数
Xs = 10 # 上限
Xx = 0 # 下限

# 随机获得初始种群
f = np.random.randint(2, size=(NP, L)) # 表示生成一个(NP, L)形状的二维数组,元素只能为0、1

# 遗传算法循环
trace = np.zeros(G) # 创建一个全为0的数组
for k in range(G):
# 将二进制解码为定义域范围内十进制
x = np.zeros(NP)
for i in range(NP):
m = np.sum(f[i,:] * 2 ** np.arange(L)[::-1]) # 计算二进制数组f[i,:]的十进制等价值
x[i] = Xx + m * (Xs - Xx) / (2**L - 1) # 将十进制数映射到区间上下界
Fit = np.array([func1(xi) for xi in x]) # 计算适应度值

maxFit = np.max(Fit) # 最大值
minFit = np.min(Fit) # 最小值
rr = np.where(Fit == maxFit)[0] # 使用np.where函数来找出数组Fit中等于maxFit的元素的索引,np.where会返回一个元组,其中包含行号和列号,[0]表示行号
fBest = f[rr[0], :] # 历代最优个体
xBest = x[rr[0]] # 最大适应度个体对应的数值(自变量)
Fit = (Fit - minFit) / (maxFit - minFit) # 归一化适应度值

# 基于轮盘赌的复制操作
sum_Fit = np.sum(Fit)
fitvalue = Fit / sum_Fit
fitvalue = np.cumsum(fitvalue) # 计算归一化适应度值的累积和
ms = np.sort(np.random.rand(NP))
fiti = 0
newi = 0
nf = np.zeros((NP, L), dtype=int) # 创建一个大小为(NP, L)的零矩阵nf,用于存储新种群的个体
while newi < NP:
if ms[newi] < fitvalue[fiti]:
nf[newi, :] = f[fiti, :]
newi += 1
else:
fiti += 1

# 基于概率的交叉操作
for i in range(0, NP, 2):
p = np.random.rand()
if p < Pc:
q = np.random.randint(1, L + 1)
mask = np.random.randint(1, L + 1, size=L) == q
nf[i, mask], nf[i + 1, mask] = nf[i + 1, mask], nf[i, mask]

# 基于概率的变异操作
i = 0
while i <= round(NP * Pm):
h = np.random.randint(1, NP + 1) - 1
for _ in range(round(L * Pm)):
g = np.random.randint(1, L + 1) - 1
nf[h, g] = 1 - nf[h, g]
i += 1

f = nf
f[0, :] = fBest # 保留最优个体在新种群中
trace[k] = maxFit # 历代最优适应度

# 输出最优个体
print("xBest:", xBest)

# 绘图
plt.plot(trace)
plt.xlabel('迭代次数')
plt.ylabel('目标函数值')
plt.title('适应度进化曲线')
plt.show()
1
xBest: 7.8567865913263235
alt text

二、实值遗传算法求函数极值

GA22.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
%%%%%%%%%%%%%%%%%%%%实值遗传算法求函数极值%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; %清除所有变量
close all; %清图
clc; %清屏
D=10; %基因数目
NP=100; %染色体数目
Xs=20; %上限
Xx=-20; %下限
G=1000; %最大遗传代数
f=zeros(D,NP); %初始种群赋空间
nf=zeros(D,NP); %子种群赋空间
Pc=0.8; %交叉概率
Pm=0.1; %变异概率
f=rand(D,NP)*(Xs-Xx)+Xx; %随机获得初始种群
%%%%%%%%%%%%%%%%%%%%%%按适应度升序排列%%%%%%%%%%%%%%%%%%%%%%%
for np=1:NP
MSLL(np)=func2(f(:,np));
end
[SortMSLL,Index]=sort(MSLL);
Sortf=f(:,Index);
%%%%%%%%%%%%%%%%%%%%%%%遗传算法循环%%%%%%%%%%%%%%%%%%%%%%%%%%
for gen=1:G
%%%%%%%%%%%%%%采用君主方案进行选择交叉操作%%%%%%%%%%%%%%%%
Emper=Sortf(:,1); %君主染色体
NoPoint=round(D*Pc); %每次交叉点的个数
PoPoint=randint(NoPoint,NP/2,[1 D]); %交叉基因的位置
nf=Sortf;
for i=1:NP/2
nf(:,2*i-1)=Emper;
nf(:,2*i)=Sortf(:,2*i);
for k=1:NoPoint
nf(PoPoint(k,i),2*i-1)=nf(PoPoint(k,i),2*i);
nf(PoPoint(k,i),2*i)=Emper(PoPoint(k,i));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%变异操作%%%%%%%%%%%%%%%%%%%%%%%%%
for m=1:NP
for n=1:D
r=rand(1,1);
if r<Pm
nf(n,m)=rand(1,1)*(Xs-Xx)+Xx;
end
end
end
%%%%%%%%%%%%%%%%%%%%%子种群按适应度升序排列%%%%%%%%%%%%%%%%%%
for np=1:NP
NMSLL(np)=func2(nf(:,np));
end
[NSortMSLL,Index]=sort(NMSLL);
NSortf=nf(:,Index);
%%%%%%%%%%%%%%%%%%%%%%%%%产生新种群%%%%%%%%%%%%%%%%%%%%%%%%%%
f1=[Sortf,NSortf]; %子代和父代合并
MSLL1=[SortMSLL,NSortMSLL]; %子代和父代的适应度值合并
[SortMSLL1,Index]=sort(MSLL1); %适应度按升序排列
Sortf1=f1(:,Index); %按适应度排列个体
SortMSLL=SortMSLL1(1:NP); %取前NP个适应度值
Sortf=Sortf1(:,1:NP); %取前NP个个体
trace(gen)=SortMSLL(1); %历代最优适应度值
end
Bestf=Sortf(:,1); %最优个体
trace(end) %最优值
figure
plot(trace)
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')

GA22.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
# 这里需要补一个对应的python代码
import numpy as np
import matplotlib.pyplot as plt

# 定义目标函数
def func2(x):
return np.sum(x**2)

# 初始化参数
D = 10 # 基因数目
NP = 100 # 染色体数目
Xs = 20 # 上限
Xx = -20 # 下限
G = 1000 # 最大遗传代数
Pc = 0.8 # 交叉概率
Pm = 0.1 # 变异概率

# 随机获得初始种群
f = np.random.rand(D, NP) * (Xs - Xx) + Xx

# 按适应度升序排列
MSLL = np.array([func2(f[:, np]) for np in range(NP)]) # 对每一列应用函数func2,接收一列数据并返回一个值,所以MSLL是一个一维数组
Index = np.argsort(MSLL) # 对一维数组进行排序,得到一个一维整数数组
Sortf = f[:, Index] # Sortf是一个与f形状相同的数组,但它的列按照MSLL中值的升序排序

# 遗传算法循环
trace = []
for gen in range(G):
# 采用君主方案进行选择交叉操作
Emper = Sortf[:, 0] # 君主染色体
NoPoint = round(D * Pc) # 每次交叉点的个数
PoPoint = np.random.randint(0, D, (NoPoint, NP // 2)) # 交叉基因的位置
nf = Sortf.copy()
for i in range(NP // 2):
nf[:, 2*i] = Emper
nf[:, 2*i + 1] = Sortf[:, 2*i + 1]
for k in range(NoPoint):
nf[PoPoint[k, i], 2*i] = nf[PoPoint[k, i], 2*i + 1]
nf[PoPoint[k, i], 2*i + 1] = Emper[PoPoint[k, i]]

# 变异操作
for m in range(NP):
for n in range(D):
if np.random.rand() < Pm:
nf[n, m] = np.random.rand() * (Xs - Xx) + Xx

# 子种群按适应度升序排列
NMSLL = np.array([func2(nf[:, np]) for np in range(NP)])
Index = np.argsort(NMSLL)
NSortf = nf[:, Index]

# 产生新种群
f1 = np.hstack((Sortf, NSortf)) # 子代和父代合并
MSLL1 = np.hstack((MSLL, NMSLL)) # 子代和父代的适应度值合并
Index = np.argsort(MSLL1) # 适应度按升序排列
Sortf1 = f1[:, Index] # 按适应度排列个体
SortMSLL = MSLL1[Index][:NP] # 取前NP个适应度值
Sortf = Sortf1[:, :NP] # 取前NP个个体
trace.append(SortMSLL[0]) # 历代最优适应度值

Bestf = Sortf[:, 0] # 最优个体
print("最优值:", trace[-1])

# 绘制适应度进化曲线
plt.plot(trace)
plt.xlabel('迭代次数')
plt.ylabel('目标函数值')
plt.title('适应度进化曲线')
plt.show()

三、遗传算法解决TSP问题

GA23.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
%%%%%%%%%%%%%%%%%%%%%%%%%遗传算法解决TSP问题%%%%%%%%%%%%%%%%%%%%%%%
clear all; %清除所有变量
close all; %清图
clc; %清屏
C=[1304 2312;3639 1315;4177 2244;3712 1399;3488 1535;3326 1556;...
3238 1229;4196 1044;4312 790;4386 570;3007 1970;2562 1756;...
2788 1491;2381 1676;1332 695;3715 1678;3918 2179;4061 2370;...
3780 2212;3676 2578;4029 2838;4263 2931;3429 1908;3507 2376;...
3394 2643;3439 3201;2935 3240;3140 3550;2545 2357;2778 2826;...
2370 2975]; %31个省会城市坐标
N=size(C,1); %TSP问题的规模,即城市数目
D=zeros(N); %任意两个城市距离间隔矩阵
%%%%%%%%%%%%%%%%%%%%%求任意两个城市距离间隔矩阵%%%%%%%%%%%%%%%%%%%%%
for i=1:N
for j=1:N
D(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5;
end
end
NP=200; %种群规模
G=1000; %最大遗传代数
f=zeros(NP,N); %用于存储种群
F=[]; %种群更新中间存储
for i=1:NP
f(i,:)=randperm(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
for i=1:W
x=find(A==B(p+i-1));
y=find(B==A(p+i-1));
temp=A(p+i-1);
A(p+i-1)=B(p+i-1);
B(p+i-1)=temp;
temp=A(x);
A(x)=B(y);
B(y)=temp;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%变异操作%%%%%%%%%%%%%%%%%%%%%%%%%
p1=floor(1+N*rand());
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
if aa>NP
F=F(1:NP,:); %保持种群规模为n
end
f=F; %更新种群
f(1,:)=R; %保留每代最优个体
clear F;
gen=gen+1
Rlength(gen)=minlen;
end
figure
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(['优化最短距离:',num2str(minlen)]);
figure
plot(Rlength)
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')

alt text alt text

alt text alt text

GA23.py

对应的python脚本如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
import numpy as np
import matplotlib.pyplot as plt

# 数据初始化
C = np.array([
[1304, 2312], [3639, 1315], [4177, 2244], [3712, 1399], [3488, 1535],
[3326, 1556], [3238, 1229], [4196, 1044], [4312, 790], [4386, 570],
[3007, 1970], [2562, 1756], [2788, 1491], [2381, 1676], [1332, 695],
[3715, 1678], [3918, 2179], [4061, 2370], [3780, 2212], [3676, 2578],
[4029, 2838], [4263, 2931], [3429, 1908], [3507, 2376], [3394, 2643],
[3439, 3201], [2935, 3240], [3140, 3550], [2545, 2357], [2778, 2826],
[2370, 2975]
]) # 31个省会城市坐标

N = C.shape[0] # 城市数目
D = np.zeros((N, N)) # 距离矩阵

# 计算距离矩阵
for i in range(N):
for j in range(N):
D[i, j] = np.sqrt((C[i, 0] - C[j, 0])**2 + (C[i, 1] - C[j, 1])**2)

NP = 200 # 种群规模
G = 1000 # 最大遗传代数

# 初始化种群
f = np.array([np.random.permutation(N) for _ in range(NP)])
R = f[0] # 存储最优种群
len_paths = np.zeros(NP) # 路径长度
fitness = np.zeros(NP) # 适应值

gen = 0
# 遗传算法主循环
Rlength = []

while gen < G:
# 计算路径长度
for i in range(NP):
path = f[i]
len_paths[i] = D[path[-1], path[0]]
for j in range(N - 1):
len_paths[i] += D[path[j], path[j + 1]]

maxlen = np.max(len_paths) # 最长路径
minlen = np.min(len_paths) # 最短路径

# 更新最短路径
best_index = np.argmin(len_paths)
R = f[best_index]

# 计算归一化适应值
fitness = 1 - ((len_paths - minlen) / (maxlen - minlen + 0.001))

# 选择操作
F = f[fitness >= np.random.rand(NP)]

while len(F) < NP:
nnper = np.random.permutation(len(F))
A = F[nnper[0]]
B = F[nnper[1]]

# 交叉操作
W = np.ceil(N / 10).astype(int) # 交叉点个数
p = np.random.randint(0, N - W + 1) # 选择交叉点
for i in range(W):
x = np.where(A == B[p + i])[0][0]
y = np.where(B == A[p + i])[0][0]
A[p + i], B[p + i] = B[p + i], A[p + i]
A[x], B[y] = B[y], A[x]

# 变异操作
p1, p2 = np.random.choice(N, 2, replace=False)
A[p1], A[p2] = A[p2], A[p1]
B[p1], B[p2] = B[p2], B[p1]

F = np.vstack([F, A, B])

if len(F) > NP:
F = F[:NP] # 保持种群规模

f = F # 更新种群
f[0] = R # 保留每代最优个体

gen += 1
Rlength.append(minlen)

# 绘图
plt.figure()
for i in range(N - 1):
plt.plot([C[R[i], 0], C[R[i + 1], 0]], [C[R[i], 1], C[R[i + 1], 1]], 'bo-')
plt.plot([C[R[-1], 0], C[R[0], 0]], [C[R[-1], 1], C[R[0], 1]], 'ro-')
plt.title(f'优化最短距离: {minlen:.2f}')
plt.show()

plt.figure()
plt.plot(Rlength)
plt.xlabel('迭代次数')
plt.ylabel('目标函数值')
plt.title('适应度进化曲线')
plt.show()

alt text

发现写的 python 脚本运行效果不好,不如 matlab 脚本

从网上找了一段python代码来比较下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
# 遗传算法求解TSP问题完整代码:
# 参考链接:http://t.csdnimg.cn/J1jSj
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
import math
import random

coord = np.array([
[1304, 2312], [3639, 1315], [4177, 2244], [3712, 1399], [3488, 1535],
[3326, 1556], [3238, 1229], [4196, 1044], [4312, 790], [4386, 570],
[3007, 1970], [2562, 1756], [2788, 1491], [2381, 1676], [1332, 695],
[3715, 1678], [3918, 2179], [4061, 2370], [3780, 2212], [3676, 2578],
[4029, 2838], [4263, 2931], [3429, 1908], [3507, 2376], [3394, 2643],
[3439, 3201], [2935, 3240], [3140, 3550], [2545, 2357], [2778, 2826],
[2370, 2975]
]) # 31个省会城市坐标
w, h = coord.shape
coordinates = np.zeros((w, h), float)
for i in range(w):
for j in range(h):
coordinates[i, j] = float(coord[i, j])

# 得到距离矩阵
distance = np.zeros((w, w))
for i in range(w):
for j in range(w):
distance[i, j] = distance[j, i] = np.linalg.norm(coordinates[i] - coordinates[j])

# 种群数
count = 300
# 进化次数
iter_time = 1000
# 最优选择概率
retain_rate = 0.3 # 适应度前30%可以活下来
# 弱者生存概率
random_select_rate = 0.5
# 变异
mutation_rate = 0.1
# 改良
gailiang_N = 3000
# 适应度
def get_total_distance(x):
dista = 0
for i in range(len(x)):
if i == len(x) - 1:
dista += distance[x[i]][x[0]]
else:
dista += distance[x[i]][x[i + 1]]
return dista

# 初始种群的改良
def gailiang(x):
distance = get_total_distance(x)
gailiang_num = 0
while gailiang_num < gailiang_N:
while True:
a = random.randint(0, len(x) - 1)
b = random.randint(0, len(x) - 1)
if a != b:
break
new_x = x.copy()
temp_a = new_x[a]
new_x[a] = new_x[b]
new_x[b] = temp_a
if get_total_distance(new_x) < distance:
x = new_x.copy()
gailiang_num += 1

# 自然选择
def nature_select(population):
grad = [[x, get_total_distance(x)] for x in population]
grad = [x[0] for x in sorted(grad, key=lambda x: x[1])]
# 强者
retain_length = int(retain_rate * len(grad))
parents = grad[: retain_length]
# 生存下来的弱者
for ruozhe in grad[retain_length:]:
if random.random() < random_select_rate:
parents.append(ruozhe)
return parents

# 交叉繁殖
def crossover(parents):
target_count = count - len(parents)
children = []
while len(children) < target_count:
while True:
male_index = random.randint(0, len(parents)-1)
female_index = random.randint(0, len(parents)-1)
if male_index != female_index:
break
male = parents[male_index]
female = parents[female_index]
left = random.randint(0, len(male) - 2)
right = random.randint(left, len(male) - 1)
gen_male = male[left:right]
gen_female = female[left:right]
child_a = []
child_b = []

len_ca = 0
for g in male:
if len_ca == left:
child_a.extend(gen_female)
len_ca += len(gen_female)
if g not in gen_female:
child_a.append(g)
len_ca += 1

len_cb = 0
for g in female:
if len_cb == left:
child_b.extend(gen_male)
len_cb += len(gen_male)
if g not in gen_male:
child_b.append(g)
len_cb += 1

children.append(child_a)
children.append(child_b)
return children

# 变异操作
def mutation(children):
for i in range(len(children)):
if random.random() < mutation_rate:
while True:
u = random.randint(0, len(children[i]) - 1)
v = random.randint(0, len(children[i]) - 1)
if u != v:
break
temp_a = children[i][u]
children[i][u] = children[i][v]
children[i][v] = temp_a

def get_result(population):
grad = [[x, get_total_distance(x)] for x in population]
grad = sorted(grad, key=lambda x: x[1])
return grad[0][0], grad[0][1]

population = []
# 初始化种群
index = [i for i in range(w)]
for i in range(count):
x = index.copy()
random.shuffle(x)
gailiang(x)
population.append(x)

distance_list = []
result_cur_best, dist_cur_best = get_result(population)
distance_list.append(dist_cur_best)

i = 0
while i < iter_time:
# 自然选择
parents = nature_select(population)

# 繁殖
children = crossover(parents)

# 变异
mutation(children)

# 更新
population = parents + children

result_cur_best, dist_cur_best = get_result(population)
distance_list.append(dist_cur_best)
i = i + 1
print(result_cur_best)
print(dist_cur_best)

for i in range(len(result_cur_best)):
result_cur_best[i] += 1

result_path = result_cur_best
result_path.append(result_path[0])

print(result_path)

# 画图

X = []
Y = []
for index in result_path:
X.append(coordinates[index-1, 0])
Y.append(coordinates[index-1, 1])

plt.rcParams['font.sans-serif'] = 'SimHei' # 设置中文显示
plt.rcParams['axes.unicode_minus'] = False
plt.figure(1)
plt.plot(X, Y, '-o')
for i in range(len(X)):
plt.text(X[i] + 0.05, Y[i] + 0.05, str(result_path[i]), color='red')
plt.xlabel('横坐标')
plt.ylabel('纵坐标')
plt.title('轨迹图')

plt.figure(2)
plt.plot(np.array(distance_list))
plt.title('优化过程')
plt.ylabel('最优值')
plt.xlabel('代数({}->{})'.format(0, iter_time))
plt.show()
alt text

这个明显效果好一些

对于上述tsp问题,还可以使用蚁群算法解决,参考: http://t.csdnimg.cn/C5NqG


智能优化算法---遗传算法(GA)
http://jrhu0048.github.io/2024/07/31/zhi-neng-you-hua-suan-fa-yi-chuan-suan-fa-ga/
作者
JR.HU
发布于
2024年7月31日
更新于
2024年10月15日
许可协议