运筹系列39：ALNS使用

1. ALNS介绍

ALNS(Adaptive Large Neighborhood Search)是现在routing和scheduling里面用的很多的启发式算法，发表于2010年。其基本思路是不断destroying和repairing问题。
定义$X$为问题$I$的可行解集，c(X)表示要优化的目标。我们从一个初始解x1出发，搜索N(x1)领域范围内的最优值x2，然后再搜索N(x2)范围内的最优值x3……，这是最优梯度下降法。
第一个需要关注的问题就是领域N如何定义。在VRP问题中，2-opt算子和relocate算子可以到达的新解称为其领域。领域的大小定义了算法复杂度，例如刚刚两个算法都是n方复杂度。
ALNS的核心是destroy和repair构成的领域，如下图，很好理解。
destroy阶段，我们可以按照一定的规则优先destroy有问题的局部点，例如上述CVRP问题中的交错的路线，再比如特别长的路线，这是worst destroy或者叫做critical destroy，完全随机的叫做random destroy，根据历史信息来的叫做history based destroy。
repair阶段或者是使用启发式算法，也可以使用精确求解算法。
ALNS与一般LNS不同的地方在于，会使用多种destroy和search方法，每种方法以一定的概率出现。

这里介绍github上star最多的python框架。使用pip install alns 安装此框架。重点如下：

• 两个基本类，ALNS用于运行程序，State用于存储解
• objective用于定义目标函数
• alns.criteria来判断每次的解是否接受，已实现的包括①HillClimbing：myopic的算法；RecordToRecordTravel：设置了update的条件；SimulatedAnnealing：概率大于随机产生数时进行更新。

介绍两个例子，TSP问题和CSP问题，来说明使用方法。

2. TSP问题

2.1 tsplib95库

文档见https://tsplib95.readthedocs.io/en/stable/，2020-5-8更新了0.7.1版本，API有一些变化。
tsp文件见TSPLIB problem files，tsplib95的几个常用方法：

• load加载问题数据为data，data.node_coords是一个字典：
  
  data.trace_tours接受求解结果列表(二维)，返回总距离列表(一维)。
• load也可以加载结果，使用tours返回结果列表
  下面是xqf131.tsp问题示例：

from alns import ALNS, State
from alns.criteria import HillClimbing
import copy
import itertools
import numpy.random as rnd
import networkx as nx
import tsplib95
import tsplib95.distances as distances
import matplotlib.pyplot as plt
data = tsplib95.load('xqf131.tsp')

# These we will use in our representation of a TSP problem: a list of
# (city, coord)-tusples.
cities = [(city, tuple(coord)) for city, coord in data.node_coords.items()]
solution = tsplib95.load('xqf131.opt.tour')
optimal = data.trace_tours(solution)[0]

2.2 绘图展示

tsplib95会把数据转成networkx的格式，使用networkx的绘图功能：

def draw_graph(graph, only_nodes=False):
    """
    Helper method for drawing TSP (tour) graphs.
    """
    fig, ax = plt.subplots(figsize=(12, 6))
    func = nx.draw_networkx
    if only_nodes:
        func = nx.draw_networkx_nodes
    func(graph, data.node_coords, node_size=25, with_labels=False, ax=ax)

2.3 定义state

定义state，里面保存了当前解（点和边），此外定义了计算费用的objective函数和辅助绘图的to_graph函数。

class TspState(State):
    def __init__(self, nodes, edges):
        self.nodes = nodes
        self.edges = edges

    def copy(self):
        return copy.deepcopy(self)

    def objective(self):
        return sum(distances.euclidean(node[1], self.edges[node][1]) for node in self.nodes)
    
    def to_graph(self):
        graph = nx.Graph()

        for node, coord in self.nodes:
            graph.add_node(node, pos=coord)
 
        for node_from, node_to in self.edges.items():
            graph.add_edge(node_from[0], node_to[0])

        return graph

数据结构是这样的：

2.4 destroy算子

下面定义destroy算子，注意输入数据，第一项是当前解（完整），第二个是随机数发生器。

degree_of_destruction = 0.25
# 每次要destroy的算子个数
def edges_to_remove(state):
    return int(len(state.edges) * degree_of_destruction)
# 随机移除边
def random_removal(current, random_state):
    destroyed = current.copy()    
    for idx in random_state.choice(len(destroyed.nodes),
                                   edges_to_remove(current),
                                   replace=False):
        del destroyed.edges[destroyed.nodes[idx]]

    return destroyed

2.5 repair算子

下面定义repair算子。greedy_repair输入第一个是当前解（残缺），第二个是随机数发生器。
would_form_subcycle是为了阻止形成环。

# 不断追寻to_node下去，如果中间有断，则返回false退出；如果遇到from_node，则返回true退出
def would_form_subcycle(from_node, to_node, state):
    for step in range(1, len(state.nodes)):
        if to_node not in state.edges:
            return False
        to_node = state.edges[to_node]        
        if from_node == to_node and step != len(state.nodes) - 1:
            return True

    return False
    
def greedy_repair(current, random_state):
    visited = set(current.edges.values())
    shuffled_idcs = random_state.permutation(len(current.nodes))
    nodes = [current.nodes[idx] for idx in shuffled_idcs]

    while len(current.edges) != len(current.nodes):
        node = next(node for node in nodes 
                    if node not in current.edges)
        unvisited = {other for other in current.nodes
                     if other != node
                     if other not in visited
                     if not would_form_subcycle(node, other, current)}

        # Closest visitable node.
        nearest = min(unvisited,
                      key=lambda other: distances.euclidean(node[1], other[1]))

        current.edges[node] = nearest
        visited.add(nearest)
    return current

2.6 初始化

这里使用greedy repair进行初始化（当然也可以用其他的初始化方法），打印中间过程如下，边是不断生成的。注意这里的randomState是numpy.random里面的方法。

最后结果如下：

2.7 定义ALNS并启动iterate

实例化，注意alns的使用方法，随机数发生器时刻备着，定义完之后紧接着加入destroy operator和repair operator：

alns = ALNS(random_state)
alns.add_destroy_operator(random_removal)
alns.add_destroy_operator(path_removal)
alns.add_destroy_operator(worst_removal)
alns.add_repair_operator(greedy_repair)

启动alns的iterate函数，其中第二个参数是更新策略的权值，4个数分别表示获得destroy算子里获得全局最优、获得一步最优，以及repair算子里接受、拒绝的概率。

criterion = HillClimbing()
result = alns.iterate(initial_solution, [3, 2, 1, 0.5], 0.8, criterion,iterations=5000, collect_stats=True)
solution = result.best_state
objective = solution.objective()
print('Best heuristic objective is {0}.'.format(objective))
print('This is {0:.1f}% worse than the optimal solution, which is {1}.'.format(100 * (objective - optimal) / optimal, optimal))

2.8 绘图

_, ax = plt.subplots(figsize=(12, 6))
result.plot_objectives(ax=ax, lw=2)

3. cutting stock problem

import copy

from functools import partial

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

from alns import ALNS, State
from alns.criteria import HillClimbing
SEED = 5432
with open('640.csp') as file:
    data = file.readlines()

NUM_LINES = int(data[0])
BEAM_LENGTH = int(data[1])

# Beams to be cut from the available beams
BEAMS = [int(length)
         for datum in data[-NUM_LINES:]
         for length, amount in [datum.strip().split()]
         for _ in range(int(amount))]

print("Each available beam is of length:", BEAM_LENGTH)
print("Number of beams to be cut (orders):", len(BEAMS))
class CspState(State):
    """
    Solution state for the CSP problem. It has two data members, assignments
    and unassigned. Assignments is a list of lists, one for each beam in use.
    Each entry is another list, containing the ordered beams cut from this 
    beam. Each such sublist must sum to at most BEAM_LENGTH. Unassigned is a
    list of ordered beams that are not currently assigned to one of the
    available beams.
    """

    def __init__(self, assignments, unassigned=None):
        self.assignments = assignments
        self.unassigned = []
        
        if unassigned is not None:
            self.unassigned = unassigned

    def copy(self):
        """
        Helper method to ensure each solution state is immutable.
        """
        return CspState(copy.deepcopy(self.assignments),
                        self.unassigned.copy())

    def objective(self):
        """
        Computes the total number of beams in use.
        """
        return len(self.assignments)

    def plot(self):
        """
        Helper method to plot a solution.
        """
        _, ax = plt.subplots(figsize=(12, 6))

        ax.barh(np.arange(len(self.assignments)), 
                [sum(assignment) for assignment in self.assignments], 
                height=1)

        ax.set_xlim(right=BEAM_LENGTH)
        ax.set_yticks(np.arange(len(self.assignments), step=10))

        ax.margins(x=0, y=0)

        ax.set_xlabel('Usage')
        ax.set_ylabel('Beam (#)')

        plt.draw_if_interactive()

def wastage(assignment):
    """
    Helper method that computes the wastage on a given beam assignment.
    """
    return BEAM_LENGTH - sum(assignment)

degree_of_destruction = 0.25

def beams_to_remove(num_beams):
    return int(num_beams * degree_of_destruction)

def random_removal(state, random_state):
    """
    Iteratively removes randomly chosen beam assignments.
    """
    state = state.copy()

    for _ in range(beams_to_remove(state.objective())):
        idx = random_state.randint(state.objective())
        state.unassigned.extend(state.assignments.pop(idx))

    return state

def greedy_insert(state, random_state):
    """
    Inserts the unassigned beams greedily into the first fitting
    beam. Shuffles the unassigned ordered beams before inserting.
    """
    random_state.shuffle(state.unassigned)

    while len(state.unassigned) != 0:
        beam = state.unassigned.pop(0)

        for assignment in state.assignments:
            if beam <= wastage(assignment):
                assignment.append(beam)
                break
        else:
            state.assignments.append([beam])

    return state

rnd_state = rnd.RandomState(SEED)
state = CspState([], BEAMS.copy())
initial_solution = greedy_insert(state, rnd_state)
print("Initial solution has objective value:", initial_solution.objective())
alns = ALNS(rnd_state)
alns.add_destroy_operator(random_removal)
alns.add_repair_operator(greedy_insert)
criterion = HillClimbing()
result = alns.iterate(initial_solution, [3, 2, 1, 0.5], 0.8, criterion,
                      iterations=5000, collect_stats=True)
solution = result.best_state
objective = solution.objective()