[PYTHON] Résolvez un simple problème de voyageur de commerce à l'aide d'une machine Boltzmann avec recuit simulé
: blue_car: Objectif de cet article
** Le recuit simulé ** est appliqué pour résoudre l'un des problèmes d'optimisation combinatoire typiques, le problème du voyageur de commerce. De plus, le modèle de machine de Boltzmann est adopté comme modèle de calcul et la programmation est effectuée en Python pour trouver la solution optimale.
: blue_car: À propos de la théorie de la méthode de cuisson
Les explications détaillées des formules mathématiques et des théories de base liées à la méthode de cuisson, qui sont reprises ici, sont données dans le livre de référence "Algorithme d'optimisation" (P120-) décrit dans la "Source" inférieure, donc je les omettrai ici. Je voudrais écrire principalement sur la programmation.
: blue_car: Problème de vendeur de patrouille simple
En tant que problème de voyageur de commerce, nous utiliserons le même problème que nous avons abordé ici.
■ Résoudre le problème du voyageur de commerce avec l'algorithme génétique (GA) et sa bibliothèque (vcopt) https://qiita.com/ttlabo/items/9b00a29702d1205f6996
Considérez le tableau 1, où chaque ligne montre le nom de la ville et chaque colonne montre l'ordre des visites, comme indiqué ci-dessous. Si vous visitez les villes dans un certain ordre, la valeur du tableau sera 1, sinon la valeur sera 0. Par exemple, dans le tableau 1, l'ordre de patrouille est ** Ville E-> Ville A-> Ville D-> Ville C-> Ville B-> Ville E ** (retourne finalement à la ville de départ). Le tableau 2 définit également la distance entre les villes.
Appliquez le modèle de machine Boltzmann au problème du voyageur de commerce (TSP) qui parcourt ces cinq villes pour résoudre le problème.
: blue_car: Formule de la formule énergétique (hamiltonien) de la fonction objectif
Citant le livre de référence, l'idée de base de la méthode de cuisson est définie comme suit.
** Concept de base de la méthode de cuisson **
Fonction objective $ f (x) $: Énergie du système (hamiltonien) $ F (x) $ variable $ x $: L'état du système qui change selon une certaine distribution de probabilité Paramètre de température $ T $: variable qui contrôle la probabilité de transition d'état
En ce qui concerne ce qui précède, changez progressivement le paramètre de température $ T $ d'une valeur élevée à une valeur faible. Recherchez le point stable (valeur minimale) de $ f (x) $.
Nous devons maintenant définir la fonction objectif hamiltonienne. L'hamiltonien conçoit et formule de manière à ce que la solution optimale puisse être trouvée lorsque la valeur devient le minimum. Par exemple, ajoutez un terme qui augmente la valeur hamiltonienne lorsque vous prenez une route ou une condition de violation qui ne convient pas pour résoudre le problème. Ce terme est appelé le terme de pénalité.
Dans ce numéro, considérez ce qui suit comme hamiltonien $ H $.
\begin{eqnarray}
Hamiltonien H&=& H_A + H_B\\
H_A &=&Terme pour trouver la distance entre les villes\\
H_B &=&Terme de pénalité
\end{eqnarray}
Tout d'abord, considérons $ H_A $. Supposons que la distance entre la ville $ i $ et la ville $ j $ soit $ d_ {i, j} $, et ajustez la variable $ x $ dans chaque cellule du tableau 1 pour faire $ x_ {i, k} $. Où $ x_ {i, k} $ est la variable qui détermine s'il faut visiter la ville $ i $ au $ k $ ème endroit. Une variable qui prend 1 si vous visitez et 0 si vous ne le faites pas, et s'appelle une variable binaire.
L'itinéraire total peut être calculé en ajoutant la distance $ d_ {i, j} $ lorsque la ville $ i $ est visitée au $ k $ ème endroit et la ville $ j $ est visitée dans l'ordre $ k + 1 $ ème. C'est très bien, alors définissez $ H_A $ comme suit:
H_A = K_1\sum_{i=1}^{5}\sum_{j=1}^{5}\sum_{k=1}^{5}d_{i,j}x_{i,k}x_{j,k+1}… La cérémonie(1)
Dans la formule ci-dessus, $ K_1 $ est une constante. Faites attention à la plage de $ k $ dans le troisième $ \ sum $ (symbole sigma). Lorsque $ k $ = 5, $ x $, qui est la cible du calcul sigma, devient $ x_ {i, 5} x_ {j, 6} $, et la sixième ville à visiter est requise. Dans ce numéro, nous reviendrons sur la première ville que nous avons visitée (la première).
Ici, comme vous ne pouvez pas visiter la même ville à la suite, définissez la valeur de $ d_ {i, i} $ sur une valeur élevée et pénalisez que la valeur de $ H_A $ ne soit pas sélectionnée comme solution. .. Autrement dit, définissez la valeur de $ d_ {1,1} $, $ d_ {2,2} $ ,,, $ d_ {5,5} $. Définissez également la distance de la zone vide en bas à gauche du tableau 2 en tenant compte de la combinaison entre les villes. Vous pouvez le définir comme une matrice symétrique. Par conséquent, nous redéfinissons le tableau 2 comme suit.
Ensuite, considérez $ H_B $. Pour le tableau 1, chaque ligne et colonne contient 1 nombre de 1 et 4 nombres de 0. Cela signifie que vous ne visiterez pas la même ville plusieurs fois ou visiterez plusieurs villes dans le même ordre en même temps. En d'autres termes, la somme de chaque ligne et de chaque colonne est toujours 1. S'il n'est pas 1, il peut être considéré comme un terme de pénalité qui augmente la valeur, et $ H_B $ est défini comme suit.
H_B = K_2\Biggl[\sum_{i=1}^{5}\bigg(\sum_{k=1}^{5}x_{i,k}-1\bigg)^2 + \sum_{k=1}^{5}\bigg(\sum_{i=1}^{5}x_{i,k}-1\bigg)^2\Biggl]… La cérémonie(2)
Dans la formule ci-dessus, soit $ K_2 $ une constante.
Maintenant que nous avons défini $ H_A $ et $ H_B $, hamiltonien $ H $ est:
H = K_1\sum_{i=1}^{5}\sum_{j=1}^{5}\sum_{k=1}^{5}d_{i,j}x_{i,k}x_{j,k+1} + K_2\Biggl[\sum_{i=1}^{5}\bigg(\sum_{k=1}^{5}x_{i,k}-1\bigg)^2 + \sum_{k=1}^{5}\bigg(\sum_{i=1}^{5}x_{i,k}-1\bigg)^2\Biggl]… La cérémonie(3)
: blue_car: Pensez aux algorithmes du programme
Nous allons créer un algorithme basé sur le livre de référence. Tout d'abord, considérons la formule logique suivante, qui est à la base du modèle de machine de Boltzmann. La formule suivante s'appelle la distribution de Gibbs-Boltzmann et est une formule qui décrit le comportement des systèmes de particules tels que la température et le gaz en physique (dynamique statistique). $ f (x) $ est l'énergie du système et $ T $ est le paramètre correspondant à la température.
g(x,T) = \frac{1}{\sum_{x}^{}exp\big(\frac{-f(x)}{T}\big)}exp\big(\frac{-f(x)}{T}\big)… La cérémonie(4)
Dans le calcul du modèle de machine de Boltzmann, la température $ T $ est rapprochée de 0. En approchant la température à 0, il est possible d'approcher une distribution d'état d'équilibre par le théorème d'Ergod de la chaîne de Markov (livre de référence P123), et la solution optimale peut être obtenue.
Ici, la température $ T $ ci-dessus est définie comme une fonction du temps $ t $ comme suit. Au lieu d'abaisser la température, modifiez $ t $ pour effectuer le calcul.
T(t) = \frac{k}{ln(t+2)}… La cérémonie(5)\\
(t=0,1,2,...)
Dans l'équation ci-dessus, $ k $ est une constante et $ ln $ est une logarithmique naturelle.
L'énergie du système (hamiltonien) change chaque fois que la température T est abaissée, donc lorsque cette énergie change, définissez un indice pour déterminer s'il faut incorporer ce changement. Cet indicateur est appelé la probabilité d'acceptation. La probabilité d'acceptation est dérivée de l'équation de distribution de Gibbs-Boltzmann (4) et de la fonction d'acceptation $ F_B $ ci-dessous (livre de référence P129). Ici, $ A (x, y) $ est la probabilité d'acceptation, $ E (x) $ ou $ E (y) $ est l'énergie du système (hamiltonien), et $ g (x, T) $ est la formule de distribution de Gibbs-Boltzmann. ça ira. $ X $ est l'état avant le changement d'énergie, et $ y $ est l'état après le changement d'énergie.
Probabilité d'acceptation A(x,y) = F_B\bigg(\frac{g(y,T)}{g(x,T)}\bigg)\\
= \frac{1}{1+exp\big(\frac{E(y)-E(x)}{T}\big)}… La cérémonie(6)\\
Fonction d'acceptation F_B(s) = \frac{s}{1+s}
Maintenant que nous avons les expressions nécessaires, nous allons construire l'algorithme.
Considérez les conditions initiales du programme. Comme mentionné ci-dessus, le temps $ t $ commence à 0, s'incrémente et augmente de 1. L'incrément est effectué dans le traitement en boucle et le nombre de boucles est en boucle. $ t $ passe de 0 à loop-1.
Considérez l'état initial de la variable de type tableau bidimensionnel $ x $ (c'est le type dict dans le programme Python). $ X $ est une variable qui détermine s'il faut ou non visiter chaque ville, et s'appelle un déterminant. C'est une variable binaire qui prend 0 ou 1, et cette fois, elle est considérée comme un tableau bidimensionnel 5x5 comme indiqué dans la figure ci-dessous (le numéro d'index commence à 0).
Initialisez ce tableau de variables de décision $ x $. La méthode d'initialisation doit être définie sur 0 ou 1 au hasard.
Après le processus d'initialisation, le processus de boucle démarre. Ce processus de boucle équivaut à abaisser la température.
Le nombre de boucles est détenu par la boucle variable, mais il n'y a pas de nombre spécifique. Déterminé en fonction de l'état de fonctionnement du programme.
Tout d'abord, obtenez la température $ T $. La température $ T $ est déterminée par l'équation (5).
Sélectionnez au hasard une des variables de décision $ x $. Puisque $ x $ est un tableau bidimensionnel de matrices 5x5, nous obtenons des coordonnées aléatoires $ (i, j) $ dans chacune des directions x et y.
Inverse la valeur $ x (i, j) $ obtenue aléatoirement. Inverser signifie changer la valeur de $ x (i, j) $ à 0 s'il est 1 et à 1 s'il est 0. Inverser ce $ x (i, j) $ s'appelle "changer l'état". Calculez comment l'énergie du système (hamiltonien) change avant et après le changement d'état, et décidez s'il faut incorporer le changement. Dans ce programme, un nouveau tableau bidimensionnel dont l'état a été changé est sécurisé comme $ y $.
La probabilité de décider d'intégrer ou non ce changement est appelée probabilité d'acceptation (livre de référence P128). La probabilité d'acceptation est calculée par l'équation (6). L'état avant le changement (énergie du système) est $ E (x) $, et l'état après le changement est $ E (y) $. L'énergie du système (hamiltonien) est calculée par l'équation (3).
Concernant les changements ci-dessus, si la probabilité d'acceptation dépasse un certain seuil, elle sera acceptée, et si elle ne dépasse pas, elle ne sera pas acceptée. Accepter signifie accepter l'inversion d'une valeur $ x (i, j) $ choisie au hasard, et ne pas accepter signifie laisser la valeur $ x (i, j) $ telle quelle. Le programme définit le seuil avec une variable appelée juri_hantei. De plus, s'il est accepté, transférez le tableau bidimensionnel $ y $ vers $ x $.
Après avoir défini l'état modifié en fonction de la probabilité d'acceptation, avancez le temps d'un pas en avant avec $ t = t + 1 $.
Répétez la boucle ci-dessus. La solution finale est $ x $ à la fin de la boucle.
: blue_car: Disposition des algorithmes
Les algorithmes ci-dessus sont organisés dans l'ordre procédural comme indiqué dans la figure ci-dessous.
sa_simulation.py
1.Définition variable
Nombre de boucles=Réglez sur 50000.
Distance entre les villes di,Définissez j.
Constante hamiltonienne K1,Définissez K2.
Constante k pour déterminer la température(Équation 5)Est défini.
k1,Les valeurs de k2 et k sont déterminées en observant le mouvement du programme.
Seuil de jugement d'acceptation basé sur la probabilité d'acceptation juri_Définissez hantei.
juri_La valeur de hantei est déterminée en regardant le mouvement du programme.
2.Initialisation
Heure t=Réglez sur 0.
Variable x(Un tableau bidimensionnel)Définissez au hasard 0 ou 1 pour chaque élément de.
Définissez les fonctions utilisées dans le programme.
3.Démarrer le traitement de la boucle
formule(5)Obtenez plus de température Tp.
Au hasard x[i,j]Sélectionner. 1 dimension(axe y)Direction i, 2D(axe x)Dans chaque direction j
Identifier aléatoirement x[i,j]Sélectionner.
Même nombre d'éléments que x(Nombre de dimensions)Avec x[i,j]Seule la valeur de(1 si 0, 0 si 1)Et
Créez un tableau à deux dimensions y avec d'autres valeurs identiques à x.
Probabilité d'acceptation(Équation 6)Est calculé.
Si la probabilité d'acceptation dépasse le seuil, y remplace x.(Sinon, laissez x tel quel et ne faites rien)
t=t+Que ce soit 1.
Revenez à la position de départ de la boucle.
4.Obtenir des résultats
Le x finalement obtenu est affiché à l'écran comme la solution optimale.
Calculez la distance totale et affichez-la à l'écran.
: blue_car: Mise en œuvre du programme
La programmation cette fois est la suivante.
sa_simulation.py
import random
import math
# =====================================================================
# 1.Définition constante
# =====================================================================
#Nombre de boucles
loop = 50000
#Définition de la distance
d = [
[100, 30, 30, 25, 10],
[ 30, 100, 30, 45, 20],
[ 30, 30, 100, 25, 20],
[ 25, 45, 25, 100, 30],
[ 10, 20, 20, 30, 100]
]
#Une dimension d'un tableau à deux dimensions x(direction y)Nombre d'éléments de
len_x_y = 5
#2D du tableau 2D x(direction x)Nombre d'éléments de
len_x_x = 5
#Constante hamiltonienne
k1 = 0.01
k2 = 1
#Constante pour déterminer la température
k = 5
#Seuil pour juger d'accepter ou non en fonction de la probabilité d'acceptation
juri_hantei = 0.05
# =====================================================================
# 2.Initialisation
# =====================================================================
#Initialisation de l'heure
t = 0
#Aléatoire 0 pour les éléments de tableau 2D,Définir sur 1
x = {}
for i in range(len_x_y):
for j in range(len_x_x):
x[i,j] = random.randint(0,1)
# =====================================================================
#Définition des fonctions
# =====================================================================
#Hamiltonien
def H(x):
HA = 0
sumA = 0
sumB = 0
for i in range(len_x_y):
for j in range(len_x_y):
for k in range(len_x_x):
k_from = k
if k == len_x_x - 1:
k_to = 0
else:
k_to = k + 1
sumA += d[i][j] * x[i,k_from] * x[j, k_to]
HA += k1 * sumA
sumB1 = 0
for i in range(len_x_y):
for k in range(len_x_x):
sumB1 += x[i,k]
sumB += sumB1 * sumB1
sumB1 = 0
sumB2 = 0
for i in range(len_x_y):
for k in range(len_x_x):
sumB2 += x[i,k]
sumB -= 2 * sumB2
sumB2 = 0
for i in range(len_x_y):
sumB += 1
sumB3 = 0
for k in range(len_x_x):
for i in range(len_x_y):
sumB3 += x[i,k]
sumB += sumB3 * sumB3
sumB3 = 0
sumB4 = 0
for k in range(len_x_x):
for i in range(len_x_y):
sumB4 += x[i,k]
sumB -= 2 * sumB4
sumB4 = 0
for k in range(len_x_x):
sumB += 1
HA += k2 * sumB
return HA
#Définition de la fonction de température
def T(t):
T = k / math.log(t + 2)
return T
#Calcul de la probabilité d'acceptation
def A(x,y,T):
tmp = H(y) - H(x)
A = 1 / (1 + math.exp(tmp / T))
return A
# =====================================================================
# 3.Traitement en boucle
# =====================================================================
while t < loop:
#Obtenez la température
Tp = T(t)
#Sélectionnez au hasard un dans le tableau bidimensionnel x(Coordonnées aléatoires i,Obtenez j)
r_i = random.randint(0,len_x_x - 1)
r_j = random.randint(0,len_x_x - 1)
#Créez un tableau bidimensionnel y dans l'état y avec une seule coordonnée sélectionnée au hasard inversée.
y = {}
for i in range(len_x_y):
for j in range(len_x_x):
if i == r_i and j == r_j:
#Inverser
y[i,j] = -1 * (x[i,j] - 1)
else:
y[i,j] = x[i,j]
#Calculer la probabilité d'acceptation
juri = round(A(x,y,Tp), 2)
#Afficher les variables de temps et la probabilité d'acceptation
print('t=',t,'juri=',juri)
#Déterminer s'il faut accepter en fonction de la probabilité d'acceptation
if juri >= juri_hantei:
x = y
print('▼ Accepté.')
#Temps d'avance t
t = t + 1
# =====================================================================
# 4.Affichage des résultats(Ville visitée)
# =====================================================================
print('-----Affichage des résultats-------')
for i in range(len_x_y):
line = ''
for j in range(len_x_x):
line += ' '
if x[i,j] == 1:
line += 'o'
else:
line += '-'
print(line)
# =====================================================================
# 5.Affichage des résultats(Distance totale)
# =====================================================================
#Liste des villes à visiter
city = []
for k in range(len_x_x):
for i in range(len_x_y):
if x[i,k] == 1:
city.append(i)
#Ajouter un itinéraire de retour de la dernière ville visitée à la première ville visitée
city.append(city[0])
#Calculez la distance totale
td = 0
for c in range(len(city) - 1):
td += d[city[c]][city[c+1]]
print('total distance:', td)
: blue_car: Vérification du résultat de l'exécution
La solution obtenue par la méthode de recuit simulé est une solution approximative, pas une solution exacte. Itérer les calculs pour enfin déterminer une meilleure solution. Cette solution est appelée la solution optimale.
Par conséquent, la solution optimale obtenue peut différer légèrement à chaque fois en fonction du nombre d'itérations. Répétez le calcul itératif plusieurs fois et concluez que la solution la plus fréquente est la solution optimale finale.
La solution optimale cette fois était la suivante.
Il s'est avéré que l'ordre de patrouille des villes devrait être ** E → A → D → C → B → E ** à partir de la ville E. À propos, si la ville A est le point de départ, ce sera A → D → C → B → E → A. La distance totale est de ** 110 **.
Les points les plus difficiles cette fois ont été ** la détermination des variables k1 et k2 pour l'hamiltonien, les variables k pour la température et le seuil juri_hantei **.
Si la relation de grandeur entre k1 et k2 n'est pas bien ajustée, le terme de pénalité ne fonctionnera pas efficacement, et la solution aura des routes de patrouille qui se chevauchent et un petit nombre de villes visitées. Si la variable de température k n'est pas fixée à la valeur correcte, un débordement peut se produire dans le traitement de la fonction exp, ou la probabilité d'acceptation peut être trop petite pour être déterminée. En définissant le seuil juri_hantei à la bonne valeur, la forme idéale du scénario de traitement est que la première moitié des 50 000 processus de boucle accepte plus de changements d'état, et la seconde moitié accepte moins parce que l'hamiltonien a convergé. J'ai pu l'intégrer.
la fin
: blue_car: Source
■ Livre de référence
"Algorithme d'optimisation" Tomoharu Nagao [Auteur] Shokodo
: blue_car: Informations connexes
■ Résolvez un simple problème de recouvrement de sommets avec Blueqat https://qiita.com/ttlabo/items/b4ab222ef9c5ee99233c
■ Un exemple de résolution du problème du sac à dos avec Blueqat https://qiita.com/ttlabo/items/537564c90056a6d56879
■ Résoudre le problème du voyageur de commerce avec l'algorithme génétique (GA) et sa bibliothèque (vcopt) https://qiita.com/ttlabo/items/9b00a29702d1205f6996
■ Résolution des problèmes de sac à dos avec les outils OR de Google - Pratiquez les bases des problèmes d'optimisation des combinaisons https://qiita.com/ttlabo/items/bcadc03004418f32f49d
■ [Partie 1] Qu'est-ce que l'optimisation? - Matériel d'étude pour l'apprentissage de l'optimisation mathématique https://qiita.com/ttlabo/items/e6970c6e85cce9ff4e34
: blue_car: Opinions etc.
Si vous avez des opinions ou des corrections, veuillez nous en informer.
Recommended Posts