续篇：继续对例题的讨论
第二问：模拟退火算法
(2) 在第一问的基础上，求四条回路，它们加起来能遍历每一个节点，并且令这四条回路的总长度尽可能小。
第二问用数学语言来描述比较繁琐，字面含义已经足够清晰。它就是想要这种感觉的一个解：
Figure_4.png
我试着在上一问已经写好的模拟退火算法的基础上做一些改动。
首先，将回路集合表示为：
x = [[i_0, ..., i_n_0], [j_0, ..., j_n_1], [k_0, ..., k_n_2], [l_0, ..., l_n_3]；将“在附近邻域随机扰动”，设为令
x中多个元素随机交换位置（可以跨组别），以及一定概率地令一些元素跳到另一个组别中去（会改变不同组别的顶点数），这样的取值范围就覆盖了整个解空间。（只要覆盖了整个解空间，剩下的就交给上天安排了）
# 核心改动在于扰动函数
def get_next_x(x):
x_next = deepcopy(x)
# 多次交换
for i in range(int(T / T_0 * 8 + 2)):
# 任选两点，交换
g_1, g_2 = randint(0, num_group - 1), randint(0, num_group - 1)
i_1, i_2 = randint(0, len(x_next[g_1]) - 1), randint(0, len(x_next[g_2]) - 1)
x_next[g_1][i_1], x_next[g_2][i_2] = x_next[g_2][i_2], x_next[g_1][i_1]
# 多次跳动
for i in range(int(T / T_0 * 8) + 1):
# 任选一点和一边，将点移至该边上，
g_1, g_2 = randint(0, num_group - 1), randint(0, num_group - 1)
# 不能令某一组变空，边需在另一组
if len(x_next[g_1]) == 1 or g_1 == g_2:
continue
i_1, i_2 = randint(0, len(x_next[g_1]) - 1), randint(0, len(x_next[g_2]))
x_next[g_2].insert(i_2, x_next[g_1][i_1])
del x_next[g_1][i_1]
return x_next
此外，像这种从一维加到二维时时，还要注意浅拷贝和深拷贝的坑，因为以前被坑过所以很怕。
完整代码如下：
# q_2_1.py
# 模拟退火算法求解最短哈密顿回路
from basic import *
from random import randint, uniform, shuffle
from math import exp, log
from copy import deepcopy
num_group = 4
# 计算回路长度
def get_path_length(x):
length = 0
for path in x:
for i in range(len(path) - 1):
length += dist[path[i]][path[i + 1]]
length += dist[len(path) - 1][0]
return length
# 给予x一个随机位移
def get_next_x(x):
x_next = deepcopy(x)
# 多次交换
for i in range(int(T / T_0 * 8 + 2)):
# 任选两点，交换
g_1, g_2 = randint(0, num_group - 1), randint(0, num_group - 1)
i_1, i_2 = randint(0, len(x_next[g_1]) - 1), randint(0, len(x_next[g_2]) - 1)
x_next[g_1][i_1], x_next[g_2][i_2] = x_next[g_2][i_2], x_next[g_1][i_1]
# 多次跳动
for i in range(int(T / T_0 * 8) + 1):
# 任选一点和一边，将点移至该边上，
g_1, g_2 = randint(0, num_group - 1), randint(0, num_group - 1)
# 不能令某一组变空，边需在另一组
if len(x_next[g_1]) == 1 or g_1 == g_2:
continue
i_1, i_2 = randint(0, len(x_next[g_1]) - 1), randint(0, len(x_next[g_2]))
x_next[g_2].insert(i_2, x_next[g_1][i_1])
del x_next[g_1][i_1]
return x_next
# 采纳概率
# 自己瞎捣鼓出来的
def get_accpectable_probability(y, y_next, T):
return exp(-T_0 / T * (y_next - y) / y)
# # 温度初值
T_0 = 100000
T = T_0
# 温度终止值
T_min = 25
# x: 回路
# 例: v_0-v_1-v_2-v_0: x = [0, 1, 2]
# x初值随机
x_0 = list(range(num_pos))
shuffle(x_0)
x = [x_0[0: 7], x_0[7: 15], x_0[15: 22], x_0[22: 30]]
# y: 回路长度
y = get_path_length(x)
# 内循环次数
k = 200
# 计数器
t = 0
# 存储求解过程
x_list = [x]
y_list = [y]
x_best = deepcopy(x)
y_best = get_path_length(x)
# 开始模拟退火
while T > T_min:
# 内循环
for i in range(k):
y = get_path_length(x)
# 给予x随机位移
x_next = get_next_x(x)
# 试探新解
y_next = get_path_length(x_next)
if y_next < y:
# 更优解，一定接受新解
x = x_next
# 更新最优记录
if y_next < y_best:
y_best = y_next
x_best = deepcopy(x_next)
else:
# 更劣解，一定概率接受新解
p = get_accpectable_probability(y, y_next, T)
r = uniform(0, 1)
if r < p:
x = x_next
# 做记录
x_list.append(x)
y_list.append(y_next)
t += 1
if t % 1000 == 0:
print(T)
# 降温
T = T_0 / (1 + t) # 快速降温
# T = T_0 / log(1 + t) # 经典降温
for path in x_best:
path += [path[0]]
# 输出解
print('min length:', y_best)
print('path:', x_best)
plt.figure(1)
ax1 = plt.subplot(1, 2, 1)
ax2 = plt.subplot(1, 2, 2)
# 绘制连线图
plt.sca(ax1)
for path in x_best:
show_path(path)
# 绘制退火图
plt.sca(ax2)
x = np.linspace(1, k * len(x_list), len(x_list))
y = np.linspace(1, len(x_list), len(x_list))
for i in range(len(x_list)):
y[i] = y_list[i]
plt.plot(x, y)
plt.title('模拟退火进程', fontproperties = 'SimHei')
plt.xlabel('遍历次数', fontproperties = 'SimHei')
plt.ylabel('回路长度', fontproperties = 'SimHei')
plt.show()
多次运行，得到以下结果：
| 运行次数 | 运行时间 | 回路总长度 |
|---|---|---|
| 1 | 44.7s | 447.0 |
| 2 | 43.5s | 438.3 |
| 3 | 42.9s | 503.0 |
| 4 | 43.4s | 466.6 |
| 5 | 43.8s | 515.8 |
| 6 | 44.1s | 459.1 |
| … | … | … |
时间上可以接受，以下为一些解的形状：
Figure_5_1.png
Figure_5_2.png
Figure_5_3.png
跑了大概二十次不到，得到的最优解总长度为444.7，形状如下：（回路的长度按边长两次算，单个点的回路长度记为零）
Figure_5_4.png
大体上，这个解勉强可以令人满意（感觉最优解应该要把点7孤立，图中没有做到）。但我们也发现一个问题，就是MTSP问题中没有其他条件加以限定的话，很容易往孤立点、孤立边这种奇怪的方向发展。换言之，这种问题通常缺乏现实意义。
第二问：遗传算法[1]
为了找到更好看的解，我觉得尝试用过都说好的遗传算法。
遗传算法（Genetic Algorithm, GA）是模拟达尔文生物进化论的自然选择和遗传学机理的生物进化过程的计算模型，是一种通过模拟自然进化过程搜索最优解的方法。
其主要特点是直接对结构对象进行操作，不存在求导和函数连续性的限定；具有内在的隐并行性和更好的全局寻优能力；采用概率化的寻优方法，不需要确定的规则就能自动获取和指导优化的搜索空间，自适应地调整搜索方向。
遗传算法以一种群体中的所有个体为对象，并利用随机化技术指导对一个被编码的参数空间进行高效搜索。其中，选择、交叉和变异构成了遗传算法的遗传操作；参数编码、初始群体的设定、适应度函数的设计、遗传操作设计、控制参数设定五个要素组成了遗传算法的核心内容。
遗传算法在我看来，是一个远比模拟退火复杂的算法，一个是因为它是一个种群（一群粒子）在同时随机，而退火是一个粒子在随机。另一个是因为它的效果好坏很大程度上取决于建模的结构（尤其是编码规则）是否与其机理（交叉和变异）相吻合。相比之下，退火只要能取值范围能覆盖解空间就基本完事了（没错我就是过河拆桥）。
在这道题中，我认为最困难的地方就在于设定遗传算法的编码模式，使之能适应交叉互换这一过程。为了简化思路，我决定让每一组回路所包含的顶点数固定。显然，这样的处理比较差劲，这样直接就ban掉了很大一部分解空间，但我一时半会儿也没有好的办法，只好先做着试试。
设置染色体编码和解读规则
设染色体编码gene = [i_0, i_1, i_2, ..., i_29]为一个长度为30的整数序列，再设一个长度为30的整数序列group_index = [g_0, g_1, g_2, ..., g_29]表示一个每个位置的属于哪个组别。
# 将整数基因序列转换为回路列表
def get_path_list(gene):
path_list = [[] for i in range(num_group)]
for i in range(len(gene)):
path_list[group_index[i]].append(gene[i])
return path_list
令种群内所有染色体共用同一个group_index来解读，这样它们所表示的回路总长度便不会在微小的编码距离内发生巨大变化。
计算基因的适应度
每个基因（个体）都会表现出对环境不同的适应度，适应度低的基因更容易被淘汰，适应度高的基因更容易保留到下一代。在这里，适应度可取4条回路的总长度的倒数。
# 计算4条回路的总长度
def get_length(path_list):
length = 0
for path in path_list:
for i in range(len(path) - 1):
length += dist[path[i]][path[i + 1]]
length += dist[path[-1]][path[0]]
return length
# 获取该基因的适应度，适应度越高越容易存活
def get_gene_fitness(gene):
# 返回对应4条回路的总长度的倒数
path_list = get_path_list(gene)
length = get_length(path_list)
return 1 / length
自然选择过程
自然选择下，适应度低的个体所占比例会下降，适应度高的个体所占比例会上升（有些优势基因会重复出现），总种群数保持不变。
其中，选择方式为轮盘选择，每个基因（个体）被选中的概率等同于其适应度在适应度总和中所占比例。因此一般来说，优势个体和劣势个体的概率差距并不大，需要很多代才能看出一个大致的进化方向，这也因而保留了基因多样性。
# 根据适应度选择，随机淘汰低适应度基因，复制高适应度基因
# 返回gene_lib_next
def selection(gene_lib, gene_lib_fitness):
# 轮盘选择，根据适应度分配被选中的概率
sum_gene_lib_fitness = sum(gene_lib_fitness)
gene_lib_next = []
for i in range(len(gene_lib)):
# 向新基因库添加基因，种群总数不变
r = uniform(0, sum_gene_lib_fitness)
for i in range(len(gene_lib)):
r -= gene_lib_fitness[i]
if r <= 0:
# 选中i基因，加入新基因库
gene_lib_next.append(gene_lib[i].copy())
break
return gene_lib_next
交叉互换
交叉互换是遗传算法中独一无二的部分。两个基因（个体）通过交叉互换彼此的一部分的方式产生新的基因进入后续的演化进程。这个操作的效果需要从编码上给予配合才能发挥100%的潜力，在当前的这个编码中其实效果等同于彼此随机交换。
这一部分的核心思想（对不合法情况进行修正）是从网上借鉴的[2]，我觉得非常巧妙。
# 对基因交叉互换后产生的不合法情况进行修正
def gene_amend(gene, low, high):
# 已知gene = gene_dad[0: low] + gene_mon[low: high] + gene_dad[high:]
# 其中中段可能出现与前后段重复的基因，不符合约束条件
# 保证中段不变，对前后段不合法情况进行修复
cross_gene = gene[low: high]
raw_gene = gene[0: low] + gene[high:]
# raw_gene中重复的下标
repeated_index = [
i for i in range(len(raw_gene))
if raw_gene[i] in cross_gene
]
# 缺失的基因值
missed_g = [
g for g in range(1, 30)
if g not in gene
]
shuffle(missed_g)
# 一一填补
for j in range(len(missed_g)):
raw_gene[repeated_index[j]] = missed_g[j]
return raw_gene[0: low] + cross_gene + raw_gene[low:]
# 以gene_dad, gene_mom做染色体交叉，返回得到的2个新子代基因
def gene_crossover(gene_dad, gene_mom):
cut_1, cut_2 = randint(0, len(gene_dad) - 1), randint(0, len(gene_mom) - 1)
if cut_1 > cut_2:
cut_1, cut_2 = cut_2, cut_1
# 交叉互换
if cut_1 == cut_2:
gene_son, gene_daughter = gene_dad.copy(), gene_mom.copy()
else:
gene_son = gene_dad[0: cut_1] + gene_mom[cut_1: cut_2] + gene_dad[cut_2:]
gene_daughter = gene_dad[0: cut_1] + gene_mom[cut_1: cut_2] + gene_dad[cut_2:]
# 修复不合法情况
gene_son = gene_amend(gene_son, cut_1, cut_2)
gene_daughter = gene_amend(gene_daughter, cut_1, cut_2)
return gene_son, gene_daughter
突变
突变较交叉来说比较简单，但是是基因多样性的根本来源。如果突变概率或突变能力设置过低，则演化一定时间后所有基因将会全部趋同甚至相同，效果极差。
# 对gene做基因突变处理，返回突变得到的新基因
def mutation(gene):
# 重复多次，强化突变的影响力
for j in range(i_mutation):
# 基因突变触发概率
if uniform(0, 1) > p_mutation:
return gene.copy()
i_1, i_2 = randint(0, len(gene) - 1), randint(0, len(gene) - 1)
# 突变：两个点位位置互换
gene_next = gene.copy()
gene[i_1], gene[i_2] = gene[i_2], gene[i_1]
return gene_next
繁衍得到子代
一一配对，得到子代，简单粗暴。
# 获取该群体繁衍得到的子代
def get_offspring(gene_lib):
half = int(len(gene_lib) / 2)
gene_dad = gene_lib[0: half]
gene_mom = gene_lib[half:]
shuffle(gene_dad)
offspring = []
for j in range(half):
# 交叉互换概率
if uniform(0, 1) > p_crossover:
gene_son, gene_daughter = gene_dad[j].copy(), gene_mom[j].copy()
else:
gene_son, gene_daughter = gene_crossover(gene_dad[j], gene_mom[j])
offspring.append(gene_son)
offspring.append(gene_daughter)
return offspring
完整代码
# q_3_1.py
# 模拟退火算法求解最短哈密顿回路
from basic import *
from random import randint, uniform, shuffle
from copy import deepcopy
'''
1. 随机产生种群。
2. 根据策略判断个体的适应度，是否符合优化准则，若符合，输出最佳个体及其最优解，结束。否则，进行下一步。
3. 依据适应度选择父母，适应度高的个体被选中的概率高，适应度低的个体被淘汰。
4. 用父母的染色体按照一定的方法进行交叉，生成子代。
5. 对子代染色体进行变异。
'''
# 将整数基因序列转换为回路列表
def get_path_list(gene):
path_list = [[] for i in range(num_group)]
for i in range(len(gene)):
path_list[group_index[i]].append(gene[i])
return path_list
# 计算4条回路的总长度
def get_length(path_list):
length = 0
for path in path_list:
for i in range(len(path) - 1):
length += dist[path[i]][path[i + 1]]
length += dist[path[-1]][path[0]]
return length
# 获取该基因的适应度，适应度越高越容易存活
def get_gene_fitness(gene):
# 返回对应4条回路的总长度的倒数
path_list = get_path_list(gene)
length = get_length(path_list)
return 1 / length
# 根据适应度选择，随机淘汰低适应度基因，复制高适应度基因
# 返回gene_lib_next
def selection(gene_lib, gene_lib_fitness):
# 轮盘选择，根据适应度分配被选中的概率
sum_gene_lib_fitness = sum(gene_lib_fitness)
gene_lib_next = []
for i in range(len(gene_lib)):
# 向新基因库添加基因，种群总数不变
r = uniform(0, sum_gene_lib_fitness)
for i in range(len(gene_lib)):
r -= gene_lib_fitness[i]
if r <= 0:
# 选中i基因，加入新基因库
gene_lib_next.append(gene_lib[i].copy())
break
return gene_lib_next
# 对基因交叉互换后产生的不合法情况进行修正
def gene_amend(gene, low, high):
# 已知gene = gene_dad[0: low] + gene_mon[low: high] + gene_dad[high:]
# 其中中段可能出现与前后段重复的基因，不符合约束条件
# 保证中段不变，对前后段不合法情况进行修复
cross_gene = gene[low: high]
raw_gene = gene[0: low] + gene[high:]
# raw_gene中重复的下标
repeated_index = [
i for i in range(len(raw_gene))
if raw_gene[i] in cross_gene
]
# 缺失的基因值
missed_g = [
g for g in range(1, 30)
if g not in gene
]
shuffle(missed_g)
# 一一填补
for j in range(len(missed_g)):
raw_gene[repeated_index[j]] = missed_g[j]
return raw_gene[0: low] + cross_gene + raw_gene[low:]
# 以gene_dad, gene_mom做染色体交叉，返回得到的2个新基因
def gene_crossover(gene_dad, gene_mom):
cut_1, cut_2 = randint(0, len(gene_dad) - 1), randint(0, len(gene_mom) - 1)
if cut_1 > cut_2:
cut_1, cut_2 = cut_2, cut_1
# 交叉互换
if cut_1 == cut_2:
gene_son, gene_daughter = gene_dad.copy(), gene_mom.copy()
else:
gene_son = gene_dad[0: cut_1] + gene_mom[cut_1: cut_2] + gene_dad[cut_2:]
gene_daughter = gene_dad[0: cut_1] + gene_mom[cut_1: cut_2] + gene_dad[cut_2:]
# 修复不合法情况
gene_son = gene_amend(gene_son, cut_1, cut_2)
gene_daughter = gene_amend(gene_daughter, cut_1, cut_2)
return gene_son, gene_daughter
# 对gene做基因突变处理，返回突变得到的新基因
def mutation(gene):
# 重复多次，强化突变的影响力
for j in range(i_mutation):
# 基因突变触发概率
if uniform(0, 1) > p_mutation:
return gene.copy()
i_1, i_2 = randint(0, len(gene) - 1), randint(0, len(gene) - 1)
# 突变：两个点位位置互换
gene_next = gene.copy()
gene[i_1], gene[i_2] = gene[i_2], gene[i_1]
return gene_next
# 获取该群体繁衍得到的子代
def get_offspring(gene_lib):
half = int(len(gene_lib) / 2)
gene_dad = gene_lib[0: half]
gene_mom = gene_lib[half:]
shuffle(gene_dad)
offspring = []
for j in range(half):
# 交叉互换概率
if uniform(0, 1) > p_crossover:
gene_son, gene_daughter = gene_dad[j].copy(), gene_mom[j].copy()
else:
gene_son, gene_daughter = gene_crossover(gene_dad[j], gene_mom[j])
offspring.append(gene_son)
offspring.append(gene_daughter)
return offspring
# 分组数量
num_group = 4
# 交叉互换触发概率
p_crossover = 0.4
# 基因突变触发概率
p_mutation = 0.1
# 基因突变的影响能力
i_mutation = 5
# 种群大小
population = 1000
# 初始化
# 基因库（种群）
gene_lib = [list(range(0, num_pos)) for i in range(population)]
for gene in gene_lib:
shuffle(gene)
# 基因解读：各点所在组
group_index = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3]
shuffle(group_index)
# 遗传迭代次数
iter_max = 1000
# 记录历史最优
fitness_best = 0
gene_best = []
for iter_i in range(iter_max):
fitness_val = [get_gene_fitness(gene) for gene in gene_lib]
# 记录最优解
fitness_max = max(fitness_val)
best_index = fitness_val.index(fitness_max)
if fitness_max > fitness_best:
fitness_best = fitness_max
gene_best = gene_lib[best_index]
# 输出信息
if iter_i % 50 == 0:
print(min(fitness_val), max(fitness_val))
# 自然选择
gene_lib_next = selection(gene_lib, fitness_val)
# 交叉互换，产生后代
offspring = get_offspring(gene_lib_next)
offspring = deepcopy(gene_lib_next)
# 基因突变
offspring = [mutation(gene) for gene in offspring]
gene_lib = offspring.copy()
# 输出结果
path_list = get_path_list(gene_best)
length = get_length(path_list)
print('min length:', length)
print('path', path_list)
可视化
for path in path_list:
path += [path[0]]
show_path(path)
plt.show()
运行结果
| 运行次数 | 运行时间 | 回路总长度 |
|---|---|---|
| 1 | 80.2s | 578.8 |
| 2 | 78.1s | 664.0 |
| 3 | 79.8s | 563.8 |
| 4 | 79.6s | 528.4 |
| 5 | 78.3s | 620.92 |
| … | … | … |
结果非常令人失望，部分解如下：
Figure_6_1.png
Figure_6_2.png
Figure_6_3.png
可以看出，与模拟退火中缩环为点的情况有着质的差别，最小长度的差距稳定在100以上。这相当于给（在这种编码下的）遗传算法判了死刑。
检讨：为什么不如人
明明是希望遗传算法能凭借着更复杂的规则，更多的搜索粒子（基因群），来对模拟退火的结果进行更上一层楼的优化，然而效果却十分让人失望，这是为什么呢？
主要原因有二：
一是，把每个回路所包含的点数写死了， 而且与模拟退火所呈现的最优解方向大相径庭，因此最优解和较优解根本就不在搜索范围之内。这是编码规则一开始就决定的。
二是：交叉互换的算法在这个编码规则下，效果完全等同于不同组间的随机点位互换，并没有发挥出“交叉互换”真正的优势。在该算法中，交叉互换是不同组间随机点位互换，突变是同组内随机点位互换，它们加起来才跟模拟退火中的随机互换等效。然而遗传算法又做不到模拟退火中的跨组跳动，不能该边每一组的包含点数，因而只会表现得更差而不是更好。
当然，这不是说遗传算法很差劲，而只是我的这个编码很差劲。这也体现了遗传算法中编码的重要性。
待续
参考资料
-
超详细的遗传算法(Genetic Algorithm)解析 番茄鸡蛋炒饭被抢注啦 https://www.jianshu.com/p/ae5157c26af9 ↩
-
遗传算法——解决M-TSP多旅行商问题（基于python基本实现） 一只黍离 https://blog.csdn.net/weixin_42201701/article/details/103746595 ↩
网友评论