禁忌搜索_连续最优化

文章目录

一、引言
二、领域选择方式
- 从球内取解的方式
三、禁忌表
四、流程图
五、示例
- 5.1完全随机初始化
- 5.2 领域采用几何分区的方式

参考论文：
【1】Fitting_of_tabu_search_to_optimize_functions_of_continuous_variables
【2】基于连续函数优化的禁忌搜索算法，张晓菲,张火明
【3】一种改进的禁忌搜索算法及其在连续全局优化中的应用，郭崇慧,　岳晓晖

参考博客：TS求解二次函数的最优值

一、引言

禁忌搜索算法多用于离散变量的组合优化，如路径规划、生产调度等。但禁忌搜索也是可以用在连续变量的寻优的，其主要难点在于领域选择方式（TSP问题中是随机交换两个城市，构成候选解），第二个则在于禁忌表的生成。

二、领域选择方式

设置一个初始解s(s是多维向量，下一次的中心为下一个放入禁忌表的点)，以s为中心画个大“球”，半径为r，下一个候选解就从这个“球”内取。

从球内取解的方式

（1）完全随机生成解：由于候选集合长度不会太长（一般4——20），因此完全随机生成解，会使集合的解不够均匀，这样搜索到最优解的耗时较长。
（2）分散寻找解：对"球"分区，若候选集合长度为4，将球分为4个区，第i个候选解都从第i个区中取（满足下式），这样可以保证每次迭代候选集合的解相对均匀。

半径hih_ihi需要根据分区方式不同，有不同的方式，如：
- （1）几何分区
- （2）线性分区
- （3）等容分区

注意：梯度求解：可以根据梯度，每次求下降梯度对应的解的位置（候选解可以按上一轮最优解的梯度方向，按不同步长得到多个候选解，将适应度最高的且不在tabo_list的解放到tabo_list中），但是现实中的问题往往不可导，所以梯度算法有一定局限性。这里不做讲解。

若有约束，如一个分量x_l<x_i<x_u，则可以考虑:

三、禁忌表

对于连续变量不在禁忌单独的解或者目标函数值（数量太多了，禁忌它们无意义，总能找到不同的点，很难找到最优值，但迭代次数足够多，多，哪怕是禁忌解，就像是完全遍历解一样，也可以找到一个还可以的最优值，不建议禁忌解），因此有以下两种方式：

（1）禁忌目标函数及目标函数δ\deltaδ邻域
具体操作为：若禁忌表里已有目标函数值f(x)，候选解x′x'x′的目标函数值为f(x’)且满足∣f(x)−f(x′)∣<=δ|f(x)-f(x')|<=\delta∣f(x)−f(x′)∣<=δ，则禁忌f(x’)，反之，则将f(x’)放入禁忌表
（2）禁忌解及解的δ\deltaδ邻域
具体为若禁忌表中已有解x，候选解x’，判断若∣∣x−x′∣∣<=δ||x-x'||<=\delta∣∣x−x′∣∣<=δ（及候选解在之前的解的邻域内），则x’已经被禁忌，判断下一个点是否被禁忌。

四、流程图

五、示例

5.1完全随机初始化

不进行分区，在x[-2,2]和y=[-2,2]的范围内搜索函数最小值
禁忌表按**三、(2)**的方式禁忌目标解及其δ\deltaδ邻域，
注意这里邻域选择也不是按照纯随机，也是判断候选解是否在候选表里的其它解的δ\deltaδ邻域范围内，若在邻域范围内，重新寻找候选解，这也相当于一定意义上的分区。（若采用完全随机，达到同样的精度，搜索时间会变长）

# -*- coding: utf-8 -*-
# @Time : 2022/3/1 14:43
# @Author : Orange
# @File : algorithm.py.pyimport numpy as np
import random
import matplotlib.pyplot as plt
import timeclass Taboo_search:def __init__(self, delata, candidate_count, taboo_list_length, iteration_count, is_random=True):self.delta = delata  # 城市列表self.candidate_count = candidate_count  # 候选集合长度self.taboo_list_length = taboo_list_length  # 禁忌长度self.iteration_count = iteration_count  # 迭代次数self.global_optimal_solution, self.global_optimum = self.random_point()  # 生成初始解并计算适应度函数值self.taboo_list = []  # 禁忌表def aim_func(self, x, y):return (1 + (x + y + 1) ** 2 * (19 - 14 * x + 3 * x ** 2 - 14 * y + 6 * x * y + 3 * y ** 2)) * (30 + (2 * x - 3 * y) ** 2 * (18 - 32 * x + 12 * x ** 2 + 48 * y - 36 * x * y + 27 * y ** 2))# 随机生成初始线路def random_point(self):x = random.uniform(-2, 2)y = random.uniform(-2, 2)cost = self.aim_func(x, y)return [x, y], costdef domain_search(self):'''完全随机不分区:return:'''x = random.uniform(-2, 2)y = random.uniform(-2, 2)cost = self.aim_func(x, y)return [x, y], costdef is_in_near_domain(self, delta, before_point, now_point):"""利用二范数判断:param delta:邻域范围:return:"""diff_point = np.array(before_point) - np.array(now_point)norm_2 = np.linalg.norm(diff_point, keepdims=True)if norm_2[0] <= delta:  # 候选解在禁忌表里的解的邻域范围内print("二范数：", norm_2)return Falseelse:return Truedef single_search(self):# 生成候选集合列表init_candidate_point, init_cadidate_point_cost = self.domain_search()candidate_list = [init_candidate_point]candidate_cost_list = [init_cadidate_point_cost]point_not_in_cadidate_list_domain_flag = 1while len(candidate_list) <= self.candidate_count:  # 在候选集合里放candidate_count条不重复路径candidate_point, candidate_point_cost = self.domain_search()"""注：注释掉下面的for...break 就是完全随机"""for i in candidate_list:if not self.is_in_near_domain(self.delta, i, candidate_point):  # 判断候选解是否在候选解集合的邻域内point_not_in_cadidate_list_domain_flag = 0breakif point_not_in_cadidate_list_domain_flag:candidate_list.append(candidate_point)candidate_cost_list.append(candidate_point_cost)# 将 candidate_list按candidate_cost_list降序排序temp_candidate = zip(candidate_cost_list, candidate_list)sort_temp_candidate = sorted(temp_candidate, key=lambda x: x[0], reverse=False)candidate_cost_list = [i[0] for i in sort_temp_candidate]candidate_list = [i[1] for i in sort_temp_candidate]if candidate_cost_list[0] < self.global_optimum:# 候选集的最好解是否优于全局最优解self.global_optimum = candidate_cost_list[0]self.global_optimal_solution = candidate_list[0]if len(self.taboo_list) >= self.taboo_list_length:  # 判断该禁忌列表长度是否以达到限制，是的话移除最初始的moveself.taboo_list.remove(self.taboo_list[0])self.taboo_list.append(candidate_list[0])  # 将该候选点加入到禁忌列表else:for i in range(1, len(candidate_cost_list)):for j in self.taboo_list:if not self.is_in_near_domain(self.delta, j, candidate_list[i]):  # 解在禁忌表的邻域范围内continueif len(self.taboo_list) >= self.taboo_list_length:  # 判断该禁忌列表长度是否以达到限制，是的话移除最初始的moveself.taboo_list.remove(self.taboo_list[0])self.taboo_list.append(candidate_list[i])return# 进行taboo_search直到达到终止条件:循环100次def taboo_search(self):for i in range(self.iteration_count):self.single_search()return self.global_optimal_solution, self.global_optimumdef plot_fig(aim_func):x = np.linspace(-2, 2, 50)  # x轴范围y = np.linspace(-2, 2, 50)  # y轴范围点a, b = np.meshgrid(x, y)  # 生成X-O-Y平面z = aim_func(a, b)  # 在每一个网格点上的对应的z轴的取值ax = plt.subplot(111, projection='3d')ax.set_title('tabo_search')ax.plot_surface(a, b, z, rstride=2, cstride=1, cmap=plt.cm.Spectral)# 设置坐标轴标签ax.set_xlabel('X')ax.set_ylabel('Y')ax.set_zlabel('Z')plt.show()if __name__ == '__main__':ts_random = Taboo_search(delata=0.001, candidate_count=5, taboo_list_length=5, iteration_count=30000)# plot_fig(ts_random.aim_func)start_time1 = time.time()global_optimal_solution, global_optimum = ts_random.taboo_search()end_time1 = time.time()duration1 = (end_time1 - start_time1)# route_greedy, cost_greedy = ts_greedy.taboo_search()print("最优解：", global_optimal_solution)print("最小目标函数值：", global_optimum)print("期望最优解：(0,-1)，期望最小目标函数值：3")print("随机TS耗时：", end_time1 - start_time1)# draw_line_pic(route_greedy, cost_greedy, None, "greedy")

最优解： [-0.0009558686926132154, -1.0005154455994396]
最小目标函数值： 3.000238791878117
期望最优解：(0,-1)，期望最小目标函数值：3
随机TS耗时： 2.1790268421173096

5.2 领域采用几何分区的方式

按照论文Fitting_of_tabu_search_to_optimize_functions_of_continuous_variables几何分区的方式,

根据“球”分区的思路，我们算法的未知数有hk,h0h_k,h_0hk,h0，候选集合长度k，禁忌表长度m

# -*- coding: utf-8 -*-
# @Time : 2022/3/4 13:51
# @Author : Orange
# @File : part_domain_search.pyimport numpy as np
import random
import matplotlib.pyplot as plt
import timeclass Taboo_search:def __init__(self, delata, candidate_count, taboo_list_length, iteration_count, h_0, h_k):self.delta = delata  # 城市列表self.candidate_count = candidate_count  # 候选集合长度self.taboo_list_length = taboo_list_length  # 禁忌长度self.iteration_count = iteration_count  # 迭代次数self.global_optimal_solution, self.global_optimum = self.random_point()  # 生成初始解并计算适应度函数值self.taboo_list = []  # 禁忌表self.h_0 = h_0self.h_k = h_kdef aim_func(self, x, y):return (1 + (x + y + 1) ** 2 * (19 - 14 * x + 3 * x ** 2 - 14 * y + 6 * x * y + 3 * y ** 2)) * (30 + (2 * x - 3 * y) ** 2 * (18 - 32 * x + 12 * x ** 2 + 48 * y - 36 * x * y + 27 * y ** 2))# 随机生成初始线路def random_point(self):x = random.uniform(-2, 2)y = random.uniform(-2, 2)cost = self.aim_func(x, y)return [x, y], costdef domain_search(self):'''完全随机不分区:return:'''x = random.uniform(-2, 2)y = random.uniform(-2, 2)cost = self.aim_func(x, y)return [x, y], costdef geometrical_domain_search(self):'''几何分区:return:'''## radii h_ik = self.candidate_counth = [self.delta] * (k + 1)for i in range(1, k + 1):h[k - i + 1] = self.h_k / (2 ** i)return hdef is_in_near_domain(self, delta, before_point, now_point):"""利用二范数判断:param delta:邻域范围:return:"""diff_point = np.array(before_point) - np.array(now_point)norm_2 = np.linalg.norm(diff_point, keepdims=True)if norm_2[0] <= delta:  # 候选解在禁忌表里的解的邻域范围内print("二范数：", norm_2)return Falseelse:return Truedef single_search(self):# 生成候选集合列表candidate_list = []candidate_cost_list = []h = self.geometrical_domain_search()s = self.taboo_list[-1]  #for i in range(self.candidate_count-1):  # 在候选集合里放candidate_count条不重复路径"""采用几何分区的方式生成候选解"""X= [random.uniform(-2, 2),random.uniform(-2, 2)]candidate_point = np.array(X) / (np.linalg.norm(np.array(X), keepdims=True)[0]) * random.uniform(h[i], h[i+1]) + scandidate_point_cost = self.aim_func(candidate_point[0], candidate_point[1])candidate_list.append(candidate_point)candidate_cost_list.append(candidate_point_cost)# 将 candidate_list按candidate_cost_list降序排序temp_candidate = zip(candidate_cost_list, candidate_list)sort_temp_candidate = sorted(temp_candidate, key=lambda x: x[0], reverse=False)candidate_cost_list = [i[0] for i in sort_temp_candidate]candidate_list = [i[1] for i in sort_temp_candidate]if candidate_cost_list[0] < self.global_optimum:# 候选集的最好解是否优于全局最优解self.global_optimum = candidate_cost_list[0]self.global_optimal_solution = candidate_list[0]if len(self.taboo_list) >= self.taboo_list_length:  # 判断该禁忌列表长度是否以达到限制，是的话移除最初始的moveself.taboo_list.remove(self.taboo_list[0])self.taboo_list.append(candidate_list[0])  # 将该候选点加入到禁忌列表else:for i in range(1, len(candidate_cost_list)):for j in self.taboo_list:if not self.is_in_near_domain(self.delta, j, candidate_list[i]):  # 解在禁忌表的邻域范围内continueif len(self.taboo_list) >= self.taboo_list_length:  # 判断该禁忌列表长度是否以达到限制，是的话移除最初始的moveself.taboo_list.remove(self.taboo_list[0])self.taboo_list.append(candidate_list[i])return# 进行taboo_search直到达到终止条件:循环100次def taboo_search(self):self.taboo_list.append(self.global_optimal_solution)for i in range(self.iteration_count):self.single_search()return self.global_optimal_solution, self.global_optimumdef plot_fig(aim_func):x = np.linspace(-2, 2, 50)  # x轴范围y = np.linspace(-2, 2, 50)  # y轴范围点a, b = np.meshgrid(x, y)  # 生成X-O-Y平面z = aim_func(a, b)  # 在每一个网格点上的对应的z轴的取值ax = plt.subplot(111, projection='3d')ax.set_title('tabo_search')ax.plot_surface(a, b, z, rstride=2, cstride=1, cmap=plt.cm.Spectral)# 设置坐标轴标签ax.set_xlabel('X')ax.set_ylabel('Y')ax.set_zlabel('Z')plt.show()if __name__ == '__main__':ts_random = Taboo_search(delata=0.01, candidate_count=5, taboo_list_length=5, iteration_count=6000, h_0=0.01,h_k=1)# plot_fig(ts_random.aim_func)start_time1 = time.time()global_optimal_solution, global_optimum = ts_random.taboo_search()end_time1 = time.time()duration1 = (end_time1 - start_time1)# route_greedy, cost_greedy = ts_greedy.taboo_search()print("最优解：", global_optimal_solution)print("最小目标函数值：", global_optimum)print("期望最优解：(0,-1)，期望最小目标函数值：3")print("随机TS耗时：", end_time1 - start_time1)# draw_line_pic(route_greedy, cost_greedy, None, "greedy")

最优解： [-6.48899113e-04 -1.00113449e+00]
最小目标函数值： 3.000503656979898
期望最优解：(0,-1)，期望最小目标函数值：3
随机TS耗时： 0.8189139366149902s

可以看出达到同样的精度比随机的方式要快很多。

注：

生成s’的方式
令s′=Mx∣∣x∣∣+ss'=M \frac{x}{||x||}+ss′=M∣∣x∣∣x+s，其中hi−1<M<hih_{i-1}<M<h_{i}hi−1<M<hi,为随机数，x为任意随机数