【優化算法】Greedy Randomized Adaptive Search算法超詳細解析，附代碼實現TSP問題求解

時間 2019-12-09

標籤優化算法 greedy randomized adaptive search 算法詳細解析代碼實現 tsp 問題求解简体版

原文原文鏈接

01 概述

Greedy Randomized Adaptive Search，貪婪隨機自適應搜索（GRAS），是組合優化問題中的多起點元啓發式算法，在算法的每次迭代中，主要由兩個階段組成：構造（construction）和局部搜索（ local search）。構造（construction）階段主要用於生成一個可行解，然後該初始可行解會被放進局部搜索進行鄰域搜索，直到找到一個局部最優解爲止。python

02 總體框架

如上面所說，其實整一個算法的框架相對於其餘算法來講還算比較簡單明瞭，你們能夠先看如下總體的僞代碼：
git

GRAS主要由兩部分組成：github

Greedy_Randomized_Construction:在貪心的基礎上，加入必定的隨機因素，構造初始可行解。
Local Search：對上面構造的初始可行解進行鄰域搜索，直到找到一個局部最優解。

而後再多說兩句：算法

Repair是什麼鬼？
有時候因爲隨機因素的加入，Greedy_Randomized_Construction階段生成的解不必定都是可行解，因此爲了保證下一步的Local Search能繼續進行，加入repair算子，對解進行修復，保證其可行。架構
不是說自適應（Adaptive）嗎？我怎麼沒看到Adaptive 的過程？
別急，這個後面具體舉例的時候會詳細講到。app

03 舉個例子說明

爲了你們能更加深刻理解該算法，咱們舉一個簡單的例子來爲你們詳細講解算法的流程。框架

好了，相信你們都看懂上面的問題了（看不懂也別問我，攤手）。對於上述問題，咱們來一步一個腳印用GRAS來求解之，來，跟緊小編的腳步……dom

強調了不少次，GRAS由兩部分組成：Greedy_Randomized_Construction和Local Search，因此，在求解具體問題的時候，完成這兩部分的設計，而後按照第二節所示的框架搭起來就能夠。ide

3.1 Greedy_Randomized_Construction

這裏仍是老規矩，先上僞代碼給你們看看，而後咱們再進行講解，畢竟對於算法來講，僞代碼的做用不言而喻。
函數

第1行，一開始解是一個空集。
第2行，初始化候選元素的集合，這裏候選元素是指能放進Solution的元素（也就是目前Solution裏面沒有的），好比１，２，３……。
第3行，對候選集合的每一個元素進行評估，計算將元素x放入Solution會致使目標函數f改變量delta_x。
第5行，根據delta_x對各個元素排序，選取部分較好的候選元素組成RCL表（貪心性體如今這裏）。
第6行，隨機在RCL中選取一個元素放進Solution。（算法的隨機性）
第八、9行，更新候選元素集合，而後對每一個元素進行從新評估計算delta值。（算法的自適應性體如今這裏）

相信通過上面如此詳細的介紹，你們都懂了吧！

3.2 Local Search

關於Local Search方面的內容，相信你們學習heuristic這麼久了，就不用我多說什麼了吧：

簡單看一下僞代碼便可，主要是鄰域算子的設計，而後就是在鄰域裏面進行搜索，找到一個局部最優解爲止。而後關於鄰域搜索，有best-improving or first-improving strategy 兩種策略，這個下次有時間出個專題給你們講明白一些相關概念吧。

04 再論Greedy_Randomized_Construction

前面咱們說了，Greedy_Randomized_Construction用於生成初始解，既然是Greedy_Randomized兩個結合體，那麼確定就有一個權重分配的問題，即，是Greedy成分多一點呢？仍是Randomized成分多一點好呢？所以，爲了控制這兩個小老弟的權重，防止某個傢伙在該過程當中用力過猛致使解不那麼好的狀況，咱們引入一個參數α：

其餘部分就再也不多說，能夠看到，上面的α參數主要是控制RCL的長度：

當α=0時，純貪心，只能選取最優的候選元素。
當α=1時，純隨機，全部候選元素均可隨機選。

05 代碼實現

因爲小編精力有限，就不從頭寫一遍了，從GitHub上找了一個感受還不錯的算法給你們，也是求解TSP問題的。不過說實在的，python寫算法的速度是很慢的，不管是速度仍是算法架構等方面都不推薦你們用matlab或者python寫大型優化算法。

運行結果以下：

代碼算例以及相關運行結果請關注公衆號【程序猿聲】，後臺回覆：GRAS，便可下載

############################################################################

# Created by: Prof. Valdecy Pereira, D.Sc.
# UFF - Universidade Federal Fluminense (Brazil)
# email:  valdecy.pereira@gmail.com
# Course: Metaheuristics
# Lesson: Local Search-GRASP

# Citation: 
# PEREIRA, V. (2018). Project: Metaheuristic-Local_Search-GRASP, File: Python-MH-Local Search-GRASP.py, GitHub repository: <https://github.com/Valdecy/Metaheuristic-Local_Search-GRASP>

############################################################################

# Required Libraries
import pandas as pd
import random
import numpy  as np
import copy
import os
from matplotlib import pyplot as plt 

# Function: Tour Distance
def distance_calc(Xdata, city_tour):
    distance = 0
    for k in range(0, len(city_tour[0])-1):
        m = k + 1
        distance = distance + Xdata.iloc[city_tour[0][k]-1, city_tour[0][m]-1]            
    return distance

# Function: Euclidean Distance 
def euclidean_distance(x, y):       
    distance = 0
    for j in range(0, len(x)):
        distance = (x.iloc[j] - y.iloc[j])**2 + distance   
    return distance**(1/2) 

# Function: Initial Seed
def seed_function(Xdata):
    seed = [[],float("inf")]
    sequence = random.sample(list(range(1,Xdata.shape[0]+1)), Xdata.shape[0])
    sequence.append(sequence[0])
    seed[0] = sequence
    seed[1] = distance_calc(Xdata, seed)
    return seed

# Function: Build Distance Matrix
def buid_distance_matrix(coordinates):
    Xdata = pd.DataFrame(np.zeros((coordinates.shape[0], coordinates.shape[0])))
    for i in range(0, Xdata.shape[0]):
        for j in range(0, Xdata.shape[1]):
            if (i != j):
                x = coordinates.iloc[i,:]
                y = coordinates.iloc[j,:]
                Xdata.iloc[i,j] = euclidean_distance(x, y)        
    return Xdata

# Function: Tour Plot
def plot_tour_distance_matrix (Xdata, city_tour):
    m = Xdata.copy(deep = True)
    for i in range(0, Xdata.shape[0]):
        for j in range(0, Xdata.shape[1]):
            m.iloc[i,j] = (1/2)*(Xdata.iloc[0,j]**2 + Xdata.iloc[i,0]**2 - Xdata.iloc[i,j]**2)    
    m = m.values
    w, u = np.linalg.eig(np.matmul(m.T, m))
    s = (np.diag(np.sort(w)[::-1]))**(1/2) 
    coordinates = np.matmul(u, s**(1/2))
    coordinates = coordinates.real[:,0:2]
    xy = pd.DataFrame(np.zeros((len(city_tour[0]), 2)))
    for i in range(0, len(city_tour[0])):
        if (i < len(city_tour[0])):
            xy.iloc[i, 0] = coordinates[city_tour[0][i]-1, 0]
            xy.iloc[i, 1] = coordinates[city_tour[0][i]-1, 1]
        else:
            xy.iloc[i, 0] = coordinates[city_tour[0][0]-1, 0]
            xy.iloc[i, 1] = coordinates[city_tour[0][0]-1, 1]
    plt.plot(xy.iloc[:,0], xy.iloc[:,1], marker = 's', alpha = 1, markersize = 7, color = 'black')
    plt.plot(xy.iloc[0,0], xy.iloc[0,1], marker = 's', alpha = 1, markersize = 7, color = 'red')
    plt.plot(xy.iloc[1,0], xy.iloc[1,1], marker = 's', alpha = 1, markersize = 7, color = 'orange')
    return

# Function: Tour Plot
def plot_tour_coordinates (coordinates, city_tour):
    coordinates = coordinates.values
    xy = pd.DataFrame(np.zeros((len(city_tour[0]), 2)))
    for i in range(0, len(city_tour[0])):
        if (i < len(city_tour[0])):
            xy.iloc[i, 0] = coordinates[city_tour[0][i]-1, 0]
            xy.iloc[i, 1] = coordinates[city_tour[0][i]-1, 1]
        else:
            xy.iloc[i, 0] = coordinates[city_tour[0][0]-1, 0]
            xy.iloc[i, 1] = coordinates[city_tour[0][0]-1, 1]
    plt.plot(xy.iloc[:,0], xy.iloc[:,1], marker = 's', alpha = 1, markersize = 7, color = 'black')
    plt.plot(xy.iloc[0,0], xy.iloc[0,1], marker = 's', alpha = 1, markersize = 7, color = 'red')
    plt.plot(xy.iloc[1,0], xy.iloc[1,1], marker = 's', alpha = 1, markersize = 7, color = 'orange')
    return

# Function: Rank Cities by Distance
def ranking(Xdata, city = 0):
    rank = pd.DataFrame(np.zeros((Xdata.shape[0], 2)), columns = ['Distance', 'City'])
    for i in range(0, rank.shape[0]):
        rank.iloc[i,0] = Xdata.iloc[i,city]
        rank.iloc[i,1] = i + 1
    rank = rank.sort_values(by = 'Distance')
    return rank

# Function: RCL
def restricted_candidate_list(Xdata, greediness_value = 0.5):
    seed = [[],float("inf")]
    sequence = []
    sequence.append(random.sample(list(range(1,Xdata.shape[0]+1)), 1)[0])
    for i in range(0, Xdata.shape[0]):
        count = 1
        rand = int.from_bytes(os.urandom(8), byteorder = "big") / ((1 << 64) - 1)
        if (rand > greediness_value and len(sequence) < Xdata.shape[0]):
            next_city = int(ranking(Xdata, city = sequence[-1] - 1).iloc[count,1])
            while next_city in sequence:
                count = np.clip(count+1,1,Xdata.shape[0]-1)
                next_city = int(ranking(Xdata, city = sequence[-1] - 1).iloc[count,1])
            sequence.append(next_city)
        elif (rand <= greediness_value and len(sequence) < Xdata.shape[0]):
            next_city = random.sample(list(range(1,Xdata.shape[0]+1)), 1)[0]
            while next_city in sequence:
                next_city = int(random.sample(list(range(1,Xdata.shape[0]+1)), 1)[0])
            sequence.append(next_city)
    sequence.append(sequence[0])
    seed[0] = sequence
    seed[1] = distance_calc(Xdata, seed)
    return seed

# Function: 2_opt
def local_search_2_opt(Xdata, city_tour):
    tour = copy.deepcopy(city_tour)
    best_route = copy.deepcopy(tour)
    seed = copy.deepcopy(tour)  
    for i in range(0, len(tour[0]) - 2):
        for j in range(i+1, len(tour[0]) - 1):
            best_route[0][i:j+1] = list(reversed(best_route[0][i:j+1]))           
            best_route[0][-1]  = best_route[0][0]                          
            best_route[1] = distance_calc(Xdata, best_route)           
            if (best_route[1] < tour[1]):
                tour[1] = copy.deepcopy(best_route[1])
                for n in range(0, len(tour[0])): 
                    tour[0][n] = best_route[0][n]          
            best_route = copy.deepcopy(seed) 
    return tour

# Function: GRASP
def greedy_randomized_adaptive_search_procedure(Xdata, city_tour, iterations = 50, rcl = 25, greediness_value = 0.5):
    count = 0
    best_solution = copy.deepcopy(city_tour)
    while (count < iterations):
        rcl_list = []
        for i in range(0, rcl):
            rcl_list.append(restricted_candidate_list(Xdata, greediness_value = greediness_value))
        candidate = int(random.sample(list(range(0,rcl)), 1)[0])
        city_tour = local_search_2_opt(Xdata, city_tour = rcl_list[candidate])
        while (city_tour[0] != rcl_list[candidate][0]):
            rcl_list[candidate] = copy.deepcopy(city_tour)
            city_tour = local_search_2_opt(Xdata, city_tour = rcl_list[candidate])
        if (city_tour[1] < best_solution[1]):
            best_solution = copy.deepcopy(city_tour) 
        count = count + 1
        print("Iteration =", count, "-> Distance =", best_solution[1])
    print("Best Solution =", best_solution)
    return best_solution

######################## Part 1 - Usage ####################################

X = pd.read_csv('Python-MH-Local Search-GRASP-Dataset-01.txt', sep = '\t') #17 cities = 1922.33
seed = seed_function(X)
lsgrasp = greedy_randomized_adaptive_search_procedure(X, city_tour = seed, iterations = 5, rcl = 5, greediness_value = 0.5)
plot_tour_distance_matrix(X, lsgrasp) # Red Point = Initial city; Orange Point = Second City # The generated coordinates (2D projection) are aproximated, depending on the data, the optimum tour may present crosses.

Y = pd.read_csv('Python-MH-Local Search-GRASP-Dataset-02.txt', sep = '\t') # Berlin 52 = 7544.37
X = buid_distance_matrix(Y)
seed = seed_function(X)
lsgrasp = greedy_randomized_adaptive_search_procedure(X, city_tour = seed, iterations = 10, rcl = 15, greediness_value = 0.5)
plot_tour_coordinates (Y, lsgrasp) # Red Point = Initial city; Orange Point = Second City