Cet article est le 13e jour du Furukawa Lab Advent_calendar. Cet article a été rédigé par un étudiant du Furukawa Lab dans le cadre de leur apprentissage. Le contenu peut être ambigu ou l'expression peut être légèrement différente.

introduction

Bonjour à tous Connaissez-vous le problème du voyageur de commerce? C'est une sorte de problème d'optimisation de combinaison courant, mais dans cet article, j'aimerais parler d'un cas triste où j'ai essayé de résoudre correctement le problème du voyageur de commerce. Pour conclure d'abord, l'erreur d'arrondi fait peur. Comme c'était une histoire quand j'étais un débutant dans le programme, j'espère que vous pourrez la lire avec un faible obstacle dans tout.

Problème de vendeur de patrouille

Tout d'abord, je voudrais expliquer le problème du voyageur de commerce. Le problème du voyageur de commerce est "le problème de trouver la distance la plus courte pour revenir au point de départ après avoir traversé toutes les villes une fois". Par exemple, étant donné les quatre villes A, B, C et D comme indiqué dans la figure suivante, la distance la plus courte est la suivante. Dans ce cas, le nombre de villes est petit et la disposition est facile, donc la distance la plus courte peut être trouvée immédiatement, mais comme le nombre de villes augmente un peu, la distance la plus courte ne peut pas être vue d'un coup d'œil.

En fait, le nombre de combinaisons d'ordre de visite des villes est de (nombre de villes-1)!, Donc s'il s'agit de 30 villes ou 40 villes, il y aura un grand nombre de commandes de visites, il est donc très difficile de vérifier toutes les commandes de visites. C'est difficile, n'est-ce pas? Par conséquent, j'ai essayé de résoudre ce problème d'optimisation de combinaison de manière appropriée afin d'obtenir une bonne solution.

Méthode SA

La méthode la plus courante pour résoudre les problèmes d'optimisation de combinaison est la méthode appelée méthode SA (Simulated-an earing). Traduit en japonais, c'est une méthode grillée. La méthode de recuit définit la fonction d'énergie $ E $. La fonction énergétique dans ce cas est la distance totale qui fait le tour de toutes les villes et revient au point de départ. Nous modifierons l'ordre des visites dans les villes afin que cette distance totale diminue.

Dans ce cas, la distance totale sera plus longue qu'avant le remplacement. Par conséquent, je voudrais rejeter ce remplacement, mais dans la méthode SA, le remplacement est adopté même si la distance devient longue avec la probabilité indiquée dans la formule suivante.

p=\exp \left\{-\frac{1}{T}\left(E_{t+1}-E_{t}\right)\right\}

Ici, $ T $ dans la formule est un paramètre appelé température, qui agit comme du bruit. Plus le $ T $ est élevé, plus il est facile d'adopter le remplacement. De plus, si la distance diminue à la suite du remplacement, elle sera adoptée avec une probabilité de 1. Cette méthode s'appelle la méthode de la métropole.

Pour résumer le flux jusqu'à présent

Changer l'ordre des visites
Décidez s'il convient d'adopter le remplacement par la méthode Metropolis Après cela, répétez ceci pour le nombre de villes et faites un pas. Ce sera.

De plus, dans la méthode SA, la température $ T $ est progressivement réduite en fonction du nombre d'étapes par la formule suivante.

T_{t+1}=\alpha T_{t}

De plus, $ \ alpha $ prend la valeur de $ 0 <\ alpha <1 $ comme indicateur de la vitesse à laquelle la température baisse.

En abaissant la température de cette manière, le remplacement a plus de chances de se produire lorsque le nombre d'étapes est petit, et le remplacement a moins de chances de se produire lorsque le nombre d'étapes est grand, c'est-à-dire à la fin du jeu. Ensuite, la première personne peut sortir même si elle tombe dans une solution locale, et la dernière personne se mettra à jour pour qu'elle puisse décider exactement. Cela vient du fait que la résistance de l'acier est augmentée en partant d'une température élevée lors du forgeage et en le refroidissant progressivement, c'est pourquoi on l'appelle la méthode de cuisson. Cette méthode est très polyvalente dans les problèmes d'optimisation combinée.

Rébellion d'algorithme

Eh bien, c'est là que je voulais parler. J'essayais de résoudre le problème du voyageur de commerce en utilisant la méthode de bronzage ci-dessus. L'algorithme lui-même est simple, j'ai donc pu le construire immédiatement. C'est pourquoi j'ai essayé de le tester avec un plan de ville très simple. La disposition de la ville la plus simple est probablement circulaire.

Cette fois, je voudrais parler des prémisses de ce plan de ville. À propos, à titre de revue, la procédure de la méthode de cuisson est

Changer l'ordre des visites
Décidez s'il convient d'adopter le remplacement par la méthode Metropolis Après cela, répétez ceci pour le nombre de villes et faites un pas. C'était. Commençons par changer l'ordre des visites de la ville. La valeur initiale est donnée sous la forme d'un nombre aléatoire approprié.

Vous pouvez voir en un coup d'œil que ce n'est pas l'itinéraire le plus court. Ensuite, décidons de la ville à remplacer. Ceci est également décidé par des nombres aléatoires, mais à ce moment-là, j'utilisais le langage C, donc en multipliant la distribution uniforme qui produit la valeur de $ 0 <x <1 $ par le nombre de villes $ N $ (type int), $ 0 Je générais un nombre aléatoire de type int de <x <N-1 $. C'était la cause de la rébellion.

Une simulation a été réalisée en utilisant un nombre aléatoire qui produit cet entier auto-créé. Au début, le nombre d'étapes était d'environ 100, mais cela fonctionne assez bien. Il a été remplacé 600 fois. Mais un jour, un incident s'est produit. J'ai décidé d'essayer 10000 étapes (remplacement 60000 fois) par hasard et j'ai tourné la simulation. J'ai eu ce résultat. J'ai pensé que c'était étrange. Je pensais que quelque chose n'allait pas lorsque la distance totale inférieure à la valeur théorique obtenue en calculant manuellement les coordonnées était sortie. C'est naturel, car ** la ville a disparu **. De plus, ce résultat n'est pas reproductible, et une ville disparaît, deux villes disparaissent et rien ne disparaît. Si vous êtes habitué à la programmation, vous pouvez tout de suite penser que "le résultat n'est pas reproductible = le nombre aléatoire est mauvais", mais à ce moment-là j'ai commencé et je suis resté coincé dans ma tête. À la fin de l'histoire, j'ai même eu une illusion comme un certain Skynet, en disant: "L'algorithme ne juge-t-il pas que la ville devrait disparaître pour obtenir la distance la plus courte?"

Ochi

Comme certains d'entre vous l'ont peut-être remarqué, la cause était que 1 lui-même était sorti en raison de l'erreur d'arrondi dans la partie qui produit un nombre aléatoire uniforme de $ 0 <x <1 $. Si 1 est sorti, il sera remplacé par la Nième ville. Puisque le tableau de langage c est compté à partir de 0, le tableau est préparé de 0 à N-1, donc le Nième tableau n'existe pas. N'ayant rien acquis de là et le remplaçant, la ville a disparu.

Problème de SOM et de voyageur de commerce

À partir de ce moment, ce sera une histoire complètement différente, mais récemment, j'ai à nouveau essayé le problème des vendeurs itinérants, alors j'aimerais y publier les résultats. Cette fois, j'ai essayé de résoudre le problème du voyageur de commerce en utilisant SOM. SOM est très simplement des données d'observation $ \ mathbf {Y} = (\ mathbf {y} _1, \ mathbf {y} _2, ... \ mathbf {y} _N) \ in \ mathbb {R} ^ {D × N} $ et variables latentes

$ \ mathbf {X} = (\ mathbf {x} _1, \ mathbf {x} _2, ... \ mathbf {x} _K) \ in \ mathbb {R} ^ {M × K} $ Pour estimer le mappage $ \ mathbf {f}: \ mathbb {R} ^ {M} → \ mathbb {R} ^ {D} $

Pour l'algorithme détaillé, il y a Document rédigé par un professeur de notre laboratoire, donc veuillez vous y référer. Je suis heureux. Par exemple, si vous donnez des données 3D et un espace latent 2D, vous obtiendrez un tel résultat.

La figure de gauche montre les données d'observation et la cartographie, et la figure de droite montre l'espace latent.

De plus, étant donné les données bidimensionnelles et l'espace latent unidimensionnel, un tel résultat est produit. Cliquez ici pour une image montrant le résultat final ああああああ.gif

À propos, les données fournies sont les coordonnées de 48 villes des États-Unis appelées att48. C'est un peu regrettable. Je sens que je peux le faire. Le goulot d'étranglement ici est que les points de départ et d'arrivée ne sont pas connectés. Si vous pouvez le faire, vous avez terminé. Par conséquent, je vais donner une variable latente qui relie le point de départ et le point final. Oui, nous plaçons les variables latentes dans un cercle sur l'espace latent. Voici le code source et le résultat de l'exécution de SOM avec des variables latentes circulaires.

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
from numpy.random import*
import matplotlib.animation as anm
import math
import random

%matplotlib nbagg
fig = plt.figure(1,figsize=(6, 4))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
class som():
    def __init__(self,N,K,sigmin,sigmax,r,w,zeta):
        self.N = N
        self.K = K
        self.sigmin = sigmin
        self.sigmax = sigmax
        self.r = r
        self.w = w
        self.zeta = zeta
        self.hist =[]
        self.hist2 = []
        self.hist3 = []
        self.sigma=0
    def fit(self,T):
        for i in range(T):
            #print(i)
            if i==0 :
                self.kn=np.random.randint(self.K,size=(self.N))
            #if i%10 == 0 :
            #    self.sigma = self.sigma+0.5
            self.sigma=max(self.sigmax-(i/self.r)*(self.sigmax-self.sigmin),self.sigmin)
            #print(self.kn.shape)
            self.test = (self.zeta[self.kn][None,:,:]-self.zeta[:,None,:])**2
            #print("self,test=",self.test.shape)
            #print("self,zeta=",self.zeta[self.kn][None,:,:].shape)
            #print("self,zeta=",self.zeta[:,None,:].shape)
            self.h_kn= np.exp(-1/(2*self.sigma**2)*np.sum((self.zeta[self.kn][None,:,:]-self.zeta[:,None,:])**2,axis=2))
            #print(self.h_kn)
            self.g_k = np.sum(self.h_kn,axis=1)
            self.y_test =(1/(self.g_k[:,None]))*self.h_kn @ self.w 
            #print("y_test:{}",self.y_test.shape)
            self.kn= np.argmin(np.sum((self.w[:,None,:]-self.y_test[None,:,:])**2,axis=2),axis=1)
            self.X1=np.reshape(self.y_test,(K,2))
            self.hist.append(self.X1)
            self.hist2.append(self.kn)
            self.hist3.append(self.sigma)
            print(self.sigma)
        print(np.array(self.hist).shape)
        return self.hist,self.hist2,self.hist3
            #self.hist.setdefault('y',[]).append(self.y_test) 
            #self.hist.setdefault('k',[]).append(self.kn) 
    #def history(self):
        #return self.hist
a=0.1
c=-0.1
N=48#Le nombre de données
D=2#Nombre de dimensions des données
L=1#Nombre de dimensions de l'espace latent
K=300#Nombre de nœuds
T=100
zeta=np.zeros((K,2))
distance = 0
sigmax=5.0#Valeur initiale de sigma
sigmin=0.01#Valeur minimale de sigma
latent_space=np.random.normal
r=T-5#Inclinaison
rad = np.linspace(-np.pi, np.pi, K)
zetax = np.sin(rad)+1
zetay = np.cos(rad)+1
zeta[:,0] = zetax
zeta[:,1] =zetay
train_data = np.loadtxt('att48.txt')
w = np.reshape(train_data,(48,D))
i=0
SOM = som(N,K,sigmin,sigmax,r,w,zeta)
kekka,k,sigma = SOM.fit(T)
kekka = np.array(kekka)
k = np.array(k)
sigma = np.array(sigma)
#print(k.shape)
for i in range(K):
    #print(("x1",kekka[T-1,i,0]))
    #print(("x2",kekka[T-1,i+1,0]))
    #print(("y1",kekka[T-1,i,1]))
    #print(("y2",kekka[T-1,i+1,1]))
    #print(("x_delta",kekka[T-1,i,0]-kekka[T-1,i+1,0]))
    #print(("y_delta",kekka[T-1,i,1]-kekka[T-1,i+1,1]))
    if(i==K-1):
        distance += np.sqrt(abs(kekka[T-1,i,0]-kekka[T-1,0,0])**2+abs(kekka[T-1,i,1]-kekka[T-1,0,1])**2)
    else:
        distance += np.sqrt(abs(kekka[T-1,i,0]-kekka[T-1,i+1,0])**2+abs(kekka[T-1,i,1]-kekka[T-1,i+1,1])**2)
    #print("delta",np.sqrt(abs(kekka[T-1,i,0]-kekka[T-1,i+1,0])**2+abs(kekka[T-1,i,1]-kekka[T-1,i+1,1])**2))
    print(distance)
#print(kekka.shape)
def update(i, zeta,w):
#for i in range(T):
    print(i)
    ax1.cla()
    ax2.cla()
    ax1.scatter(w[:,0],w[:,1],s=100,color='r')
    ax1.scatter(kekka[i,:,0],kekka[i,:,1],color = 'k',marker = '*')
    ax1.plot(kekka[i,:,0],kekka[i,:,1],color='k')
    ax1.plot(kekka[i,:,0].T,kekka[i,:,1].T,color='k')
    ax2.scatter(zeta[k[i,:]][:,0],zeta[k[i,:]][:,1],c=w[:,0])
    ax1.set_xlabel("X_axis")
    ax1.set_ylabel("Y_axis")
    ax2.set_title('latentspace')
    ax2.set_xlabel("X_axis")
    ax2.set_ylabel("Y_axis")

ani = anm.FuncAnimation(fig, update, fargs = (zeta,w), interval = 100, frames = T-1,repeat = False)


ani.save("RSom_TSP_Node300.gif", writer = 'imagemagick')

S'il s'agit d'un gif, le résultat final disparaîtra en un instant, donc je prépare également une image. c'est ici

Le temps de calcul est rapide, donc je pense que c'est une assez bonne ligne. Après cela, je pense que ce sera encore mieux si nous pouvons sortir de la solution locale en ajoutant du bruit.