禁忌搜索_连续最优化_python

参考论文：
【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个区中取（满足下式），这样可以保证每次迭代候选集合的解相对均匀。
  
  半径$h_i$需要根据分区方式不同，有不同的方式，如：

  • （1）几何分区
    

  • （2）线性分区
    

  • （3）等容分区
    
    

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

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

三、禁忌表

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

• （1）禁忌目标函数及目标函数$\delta$邻域
  具体操作为：若禁忌表里已有目标函数值f(x)，候选解$x^{\prime }$的目标函数值为f(x’)且满足$|f(x)-f(x^{\prime })|<=\delta$，则禁忌f(x’)，反之，则将f(x’)放入禁忌表

• （2）禁忌解及解的$\delta$邻域
  具体为若禁忌表中已有解x，候选解x’，判断若$||x-x^{\prime }||<=\delta$（及候选解在之前的解的邻域内），则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.py

import numpy as np
import random
import matplotlib.pyplot as plt
import time


class 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], cost

    def domain_search(self):
        '''
        完全随机不分区
        :return:
        '''
        x = random.uniform(-2, 2)
        y = random.uniform(-2, 2)
        cost = self.aim_func(x, y)
        return [x, y], cost

    def 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 False
        else:
            return True

    def 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 = 1
        while 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 = 0
                    break
            if 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:  # 判断该禁忌列表长度是否以达到限制，是的话移除最初始的move
                self.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]):  # 解在禁忌表的邻域范围内
                        continue
                if len(self.taboo_list) >= self.taboo_list_length:  # 判断该禁忌列表长度是否以达到限制，是的话移除最初始的move
                    self.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_optimum


def 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几何分区的方式,

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

# -*- coding: utf-8 -*- 
# @Time : 2022/3/4 13:51 
# @Author : Orange
# @File : part_domain_search.py


import numpy as np
import random
import matplotlib.pyplot as plt
import time


class 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_0
        self.h_k = h_k

    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], cost

    def domain_search(self):
        '''
        完全随机不分区
        :return:
        '''
        x = random.uniform(-2, 2)
        y = random.uniform(-2, 2)
        cost = self.aim_func(x, y)
        return [x, y], cost

    def geometrical_domain_search(self):
        '''
        几何分区
        :return:
        '''
        ## radii h_i
        k = self.candidate_count
        h = [self.delta] * (k + 1)
        for i in range(1, k + 1):
            h[k - i + 1] = self.h_k / (2 ** i)

        return h

    def 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 False
        else:
            return True

    def 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]) + s
            candidate_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:  # 判断该禁忌列表长度是否以达到限制，是的话移除最初始的move
                self.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]):  # 解在禁忌表的邻域范围内
                        continue
                if len(self.taboo_list) >= self.taboo_list_length:  # 判断该禁忌列表长度是否以达到限制，是的话移除最初始的move
                    self.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_optimum


def 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^{\prime }=M\frac{x}{||x||}+s$，其中$h_{i-1}<M<h_i$,为随机数，x为任意随机数