metaheuristic 元启发式方法。一些随机搜索算法诸如进化算法、蚁群算法、粒子群算法这类具有启发式框架的智能算法称为元启发式算法。

Metaheuristic Network网站对于现代启发式算法给出的定义为：“Metaheuristic是一个用来定义启发式算法的概念集，这些启发式算法可以用来求解不同的优化问题。换句话说Metaheuristic可以被看成是一种算法框架，这种算法框架通过微小的改动可以运用到不同的优化问题。”

下图介绍了两种不同种类的Metaheuristics，我们主要用左边的，尤其是Iterated Greedy。

Metaheuristics

1 基本原则

• Intensivierung：在有可能存在最优解的值域中寻找尽可能接近最优解的方案。主要通过改变目标函数来实现。
• Diversifizierung：在尚未分析过的值域范围内进行搜索，通常我们会把随机状态等控制机制作为目标函数。

I 和D 的过程用图像表示如下：

1. 首先我们通过启发式算法可以得到一个初始方案S，这个方案通常不是最优解。

   
   
   初始方案

2. 随后我们执行local search（LS），找到一个局部最优解，这其实是I 过程：

   
   
   LS

3. 找到局部最优后，我们需要探索目前尚未到达的领域，这就是D 过程。通过加入Perturbation， P移动到值域的其他位置，并创建一个新的解决方案。然后以这个新的方案作为初始方案，重复执行前面的步骤，直到满足结束条件。

   
   
   Pertubation
   
   
   重复执行，直到满足结束条件

2 Iterated Greedy简介

我在这里只用Iterated Greedy算法。原因如下：
其他算法诸如蚁群算法等，可能能提供非常接近最优解的方案，有些算法的运行速度也很快，可以弥补python运行慢的缺陷。但是这些算法通常存在诸多不足，例如算法太复杂，适用面很窄，需要设置过多参数导致很难实现。

Iterated Greedy的优势在于，它由两个简单的阶段构成：

1. Destrucrtion解构：
   把当前的完整解决方案进行部分分解，例如原本11个订单，我从中去除两个订单进行重新安排（Construction）
2. Construction构建：
   把上一步解构的方案借助constructiveHeuristics重新处理成一个完整的方案。

acceptanceKriterium 接受条件

D和I
更好的解决方案总会被接受
更坏的方案以特定的可能性接受，接受的概率如下图
类似模拟退火法：

模拟退火概率计算

3 Iterated Greedy过程

1. 借助NEH算法创建初始方案p
2. 对p执行本地查询LS
3. 在停止条件满足之前，执行：
   • 解构： 从p中随机删去d个订单，得到p'（是个不完全排列，缺了d个元素），然后把这d个订单增加到pr中
   • 构建：再一次利用NEH，重新构建一个完整的方案p''，pr中的订单被插入到p'中（在前面Neighborhood类中，我们定义了一种NB，叫做Insertion，就是用于执行最优插入的）
   • 对重构的完整方案p''执行LS
   • 当达到结束要求或达到可接受的条件（acceptanceKriterium）时，把p''作为当前方案p（currentSolution）
   • 如果p'' 比目前为止找到的最优解pb 还要好，就把p'' 作为最优解pb存储
4. 返回找到的最优解pb

4 Iterated Greedy：初始方案的创建

from Solver import *

data = InputData("../TestInstancesJson/Small/VFR10_5_5_SIST.json")

solver = Solver(data, 1060) #1014, 754 |  1060, 736     736 已经是优化的结果了，第五中，会出现两个一样的结果

startSolution = solver.ConstructionPhase("NEH")  # 用neh创建初始方案

# 用deepcopy这样被修改的就不是start，而是current。这一步的目的是把初始方案作为当前方案，
# 后面如果发现更好的方案，则把这个更好的方案作为current；经过了deep copy后，我们最初建立的方案数据不会消失
# 后面可以用初始方案作比较
currentSolution = deepcopy(startSolution) 

5 Destruction 解构

从当前解决方案中随机删除numberJobsToRemove个订单。这里numberJobsToRemove是Iterated Greedy的一个参数。
输出的结果是被移除的订单集合removedJobs和一个不完整的解决方案partialPermutation。

numberJobsToRemove = 2 # 1st parameter

def Destruction(solution):
    # 思考：怎么随机删除？？？？
    # choice() 方法返回一个列表，元组或字符串的随机项。注意：choice()是不能直接访问的，需要导入 random 模块，然后通过 random 静态对象调用该方法。
    removedJobs = solver.RNG.choice(solver.InputData.n, size = numberJobsToRemove, replace=False).tolist() # 有放回和无放回。InputData.n是jobs的数量

    partialPermutation = [ i for i in solution.Permutation if i not in removedJobs]

    return removedJobs, partialPermutation # 两个返回值，也就是说，调用该函数的时候需要两个变量

# 使用解构函数把currentSolution分解成两部分
removedJobs, partialPermutation = Destruction(currentSolution)

print(f'removed jobs: {removedJobs}')
print(f'partial permutation: {partialPermutation}')

解构Destruction

• 详解numpy.random.choice方法的使用
  random.choice(a, size=None, replace=True, p=None)
  参数表：
  • a: 一维数组或int型变量
  • size： 可选参数，整型或由整型组成的元组，表示返回多少个随机样本。默认值为none，即返回一个样本。
  • replace：可选参数，bool型，表示抽样是否有放回。默认值是 true，表示有放回，即同一个元素可以多次抽样。在我们的例子中需要无放回抽样，因此设这个参数的值为false

*注意：solver.RNG.choice的结果每次都一样，是因为我们设置了随机数种子，目的就是只要是用同一个种子作为参数构造出的solver的属性RNG都是同一个。

6 Construction构建

将removedJobs重新加回partialPermutation，并插入到最佳位置（NEH）的顺序，并返回新的完整解决方案。这个插入过程是通过排列Permutation实现的，看一下Construction函数的参数表可知，需要两个列表。
*注：前面讲算法的时候说过，重构函数执行之后会生成一个新的方案newSolution，新方案就是通过把removedJobs插入到最佳位置得到的。最优位置的选择需要在Construction函数中借助 solver.EvaluationLogic.DetermineBestInsertion(completeSolution, i)来实现

def Construction(removedJobs, partialPermutation):
    completeSolution = Solution(solver.InputData.InputJobs, partialPermutation)  # completeSolution在这里还没有complete
    
    # 注意： 重构时要把删除的元素逐一插入回最优位置，因此要用for循环
    for i in removedJobs:
        solver.EvaluationLogic.DetermineBestInsertion(completeSolution, i) # best insertion in neh

    return completeSolution

# newSolution现在是完整的
newSolution = Construction(removedJobs, partialPermutation)
print(f'new complete solution after construction: \n{newSolution}')

重构Construction

7 接受新方案newSolution的准则

我们通过简单的解构和重构，必然会得出一个新方案newSolution，那么我们接受新方案newSolution为当前方案currentSolution的前提是，如果

• 新的解决方案比目前的解决方案更好（总会接受新方案作为当前方案）。
• 达到了接受标准（随机决定是否接受新方案）。

- 新的解决方案（newSolution），
- 当前解决方案（currentSolution）。

接受更坏方案的概率

幂函数表达式用于计算概率probability，表示接受坏结果的可能性。概率跟一个随机数比较，概率大于随机数，则接受坏结果。

*公式中，T叫做baseTemperature，是Iterated Greedy的第二个参数，T越大，则接受更坏方案的概率也越大。新的最优方案存储在SolutionPool中（numberJobsToRemove是第一个参数）

下面开始对newSolution进行评估和比较：

print(f'SOlution before checking for acceptance: \n{newSolution}')

baseTemperature = 1 # Iterated Greedy第二个参数
currentBestObj = solver.SolutionPool.GetLowestMakespanSolution().Makespan 
# 上面这行表示，目前所有解决方案（solutionPool）中用时最短的方案的总时长。把这个最短的时间作为Obj

#"""新方案和当前方案的用时进行比较，若新方案更短，则接受这个方案"""
if newSolution.Makespan < currentSolution.Makespan:
    currentSolution = newSolution 

    #"""接受新方案只是第一步，因为新方案目前仅仅优于currentSolution，现在要判断newSolution是否为最优"""
    if newSolution.Makespan < currentBestObj:
        solver.SolutionPool.AddSolution(newSolution)
        currentBestObj = newSolution.Makespan

#"""若新方案并未没有提供任何优化，那么这里计算公式的结果，试图获取有多大概率接受这个新方案，即使新方案并不好"""
else:
    newObjValue = newSolution.Makespan
    currentObjValue = currentSolution.Makespan

    # 计算所有订单加工时间的总和-->公式中的分子
    totalProcessingTime = sum(x.ProcessingTime(i) for x in solver.InputData.InputJobs for i in range(len(x.Operations)))
    temperature = baseTemperature * totalProcessingTime / (solver.InputData.n * solver.InputData.m * 10)
    probability = math.exp(-(newObjValue - currentObjValue) / temperature)

    randomNumber = solver.RNG.random()
    if randomNumber < probability:
        currentSolution = newSolution

print(f'SOlution after checking for acceptance: \n {currentSolution}')

分两种情况检查新方案是否可接受

*注意：判断一个新方案是否可接受取决于两方面：1. 如果新方案更好，那么无论如何都会接受新方案；2. 如果新方案并没有优化，则视作WorseSolution，即使是worseSolution也是要按照公式计算出的概率来衡量是否要接受这个不好的方案

为了看起来更加直观，我把接受差方案的过程写成了函数AcceptWorseSolution。注意看参数表，以便于确定何时调用这个函数。

def AcceptWorseSolution(currentObjectiveValue, newObjectiveValue):  # 这个函数的参数表分别是 Cp''和Cp
    temperature = baseTemperature * solver.InputData.TotalProcessingTime / (solver.InputData.n * solver.InputData.m * 10)
    probability = math.exp(-(newObjectiveValue - currentObjectiveValue) / temperature)

    return solver.RNG.random() <= probability  # 返回一个逻辑值，为真时，正如函数名那样，我们要接受新方案

8 局部搜索Local Search

局部搜索可以被用于优化currentSOlution，但不是必须的。通常会使用IterativeImprovement结合Insertion邻域一起使用。

print(f'Solution before LS: \n{currentSolution}')

localSearch = IterativeImprovement(solver.InputData, 'BestImprovement', ['Insertion'])
localSearch.Initialize(solver.EvaluationLogic, solver.SolutionPool)

currentSolution = localSearch.Run(currentSolution)

print(f'Solution after LS:\n{currentSolution}')

LS

9 Iterated Greedy作为各部分的控制机制

Iterated Greedy是一个迭代的解构D和构建C组成的序列。
在一个循环中被反复执行，直到达到停止标准。

在这里的例子中，停止标准是迭代次数maxIterations，但时间限制或没有优化的迭代次数也是可以的。

后面可能还会讲到怎么设计一个没有优化的迭代次数，这里可以先思考一下。

# 迭代的停止标准：总的迭代次数，即只要执行过IterGreedy，计数器就要加1。
#对于没有优化的迭代次数，可以理解为执行过AcceptWorseSolution的次数，即计数新方案不如当前方案的次数。
# AcceptWorseSolution返回一个bool值，当elif执行时，考虑到概率问题，这条语句执行的次数可能回很少，
# 因此无优化迭代的最大值通常较小
maxIterations = 10 

def Run(currentSolution):
    currentSolution = localSearch.Run(currentSolution)

    currentBestObj = solver.SolutionPool.GetLowestMakespanSolution().Makespan
    iteration = 0

    while iteration < maxIterations:
        removedJobs, partialPermutation = Destruction(currentSolution)
        newSolution = Construction(removedJobs, partialPermutation)

        newSolution = localSearch.Run(newSolution)

        if newSolution.Makespan < currentSolution.Makespan:
            currentSolution = newSolution

            if newSolution.Makespan < currentBestObj:
                print(f'new best solution found in iteration {iteration}: {currentSolution}')
                solver.SolutionPool.AddSolution(newSolution)
                currentBestObj = newSolution.Makespan

        elif AcceptWorseSolution(currentSolution.Makespan, newSolution.Makespan):
            currentSolution = newSolution

        iteration += 1

    return solver.SolutionPool.GetLowestMakespanSolution()

现在尝试运行一下：

solver = Solver(data, 1060)

startSolution = solver.ConstructionPhase("NEH")

bestSolution = Run(startSolution)

print(f'Best solution found.\n{bestSolution}') 

10 Iterated Greedy 类

与IterativeImprovement算法类似，现在要为IteratedGreedy创建一个单独的类，以便该类的一个实例可以传递给求解器Solver。

像IterativeImprovement一样，Iterated Greedy应该继承自ImprovementAlgorithm。

必要的参数是以下属性：

• NumberJobsToRemove
• BaseTemperature
• MaxIterations
• Local search算法的实例

EvaluationLogic, SolutionPool以及随机数生成器RNG都由求解器solver传递给算法。

添加IteratedGreedy类，以便它可以作为一种算法传递给求解器。

成员函数：Konstruktor, Initialize, Destruction, Construction, AcceptWorseSolution, Run

class IteratedGreedy(ImprovementAlgorithm):
    def __init__(self, inputData, numberJobsToRemove, baseTemperature, maxIterations, localSearchAlgorithm=None):
        super().__init__(inputData)

        self.NumberJobsToRemove = numberJobsToRemove
        self.BaseTemperature = baseTemperature
        self.MaxIterations = maxIterations
        
    def Destruction(self, solution): 
        # 思考：怎么随机删除？？？？
        removedJobs = self.RNG.choice(solver.InputData.n, size = self.NumberJobsToRemove, replace=False).tolist() # 哟放回和无放回

        partialPermutation = [ i for i in solution.Permutation if i not in removedJobs]

        return removedJobs, partialPermutation

    def Construction(self, removedJobs, partialPermutation):
        completeSolution = Solution(solver.InputData.InputJobs, partialPermutation)
        for i in removedJobs:
            solver.EvaluationLogic.DetermineBestInsertion(completeSolution, i) # best insertion in neh

        return completeSolution


    def AcceptWorseSolution(self, currentObjectiveValue, newObjectiveValue):  
        temperature = self.BaseTemperature * solver.InputData.TotalProcessingTime / (solver.InputData.n * solver.InputData.m * 10)
        probability = math.exp(-(newObjectiveValue - currentObjectiveValue) / temperature)

        return self.RNG.random() <= probability 

    def Run(self, currentSolution):
        currentSolution = localSearch.Run(currentSolution)

        currentBestObj = solver.SolutionPool.GetLowestMakespanSolution().Makespan
        iteration = 0

        while iteration < self.MaxIterations:
            removedJobs, partialPermutation = Destruction(currentSolution)
            newSolution = Construction(removedJobs, partialPermutation)

            newSolution = localSearch.Run(newSolution)

            if newSolution.Makespan < currentSolution.Makespan:
                currentSolution = newSolution

                if newSolution.Makespan < currentBestObj:
                    print(f'new best solution found in iteration {iteration}: {currentSolution}')
                    solver.SolutionPool.AddSolution(newSolution)
                    currentBestObj = newSolution.Makespan

            elif AcceptWorseSolution(currentSolution.Makespan, newSolution.Makespan):
                currentSolution = newSolution

            iteration += 1

        return solver.SolutionPool.GetLowestMakespanSolution()

    data = InputData("../TestInstancesJson/Small/VFR10_5_5_SIST.json") #"../TestInstancesJson/Small/VFR40_10_3_SIST.json"
# 有些很大的数据，python处理起来会很慢，可以用Java

    # 这就是Iterated Greedy的localSearchAlgorithm参数
    insertion = IterativeImprovement(data, 'BestImprovement', ['Insertion'])
    # taillardInsertion = IterativeImprovement(data, 'BestImprovement', ['TaillardInsertion'])
    iteratedGreedy = IteratedGreedy( # IterateGreedy是一种优化算法，而solver的Run local search函数需要两个参数，第一个是启发式算法，第二个是优化算法
        data, 
        numberJobsToRemove=2, 
        baseTemperature=1, 
        maxIterations=100,
        localSearchAlgorithm=insertion)

    solver = Solver(data, 1010)

    # 在进行本地搜索前需要准备好用于优化启发式算法结果的优化算法。上一章讲的优化算法是创建邻域，通过插入和换位来获得优化。这一张讲的Iterated Greedy则是另一种思路，解构和重构
    solver.RunLocalSearch(constructiveSolutionMethod='NEH', algorithm=iteratedGreedy)