[PYTHON] À propos du problème du voyageur de commerce
Récemment, j'ai eu la chance d'entrer en contact avec «l'optimisation mathématique» grâce à mon travail. C'est une bonne idée, alors j'aimerais publier les résultats d'une petite touche privée.
Ce que j'ai fait était un exemple d'optimisation mathématique, le «problème du vendeur circulaire» Comme cela peut être mis en œuvre avec la pulpe d'environnement libre, j'ai essayé de l'implémenter tout en préparant l'environnement.
Les problèmes que je veux essayer sont les suivants.
◆ 2.12 Problème de vendeur de patrouille https://www.msi.co.jp/nuopt/docs/v19/examples/html/02-12-00.html
J'ai dessiné une image simple ci-dessous. A, B, C, D sont les points sur la carte.
Quel itinéraire a la plus petite distance totale en quittant A et en revenant d'un seul coup? Il semble que le problème d'optimisation s'appelle le «problème du vendeur de patrouille». Cela s'appelle la difficulté NP, et cela semble assez difficile.
Cependant, l'approche avec le modeleur appelé pulp qui est touché dans l'environnement python est gratuite !!! C'est facile à faire, je vais donc me référer à l'article suivant.
◆ Optimisation mathématique à partir du problème du voyageur de commerce https://qiita.com/panchovie/items/6509fb54e3d53f4766aa
La grande chose est que nous avons tout, de l'installation de pulp à l'échantillon de code. Mettons immédiatement le premier paramètre dans l'exemple de code. Qu'en est-il du code suivant qui semble compiler pour le moment?
import pulp as pp
#Définition du modèle d'optimisation
mip_model = pp.LpProblem("tsp_mip", pp.LpMinimize)
N = 4
BigM = 10000
c_ = [
[0,6,5,5], # pointA > [pointA,pointB,pointC,pointD]
[6,0,7,4], # pointB > [pointA,pointB,pointC,pointD]
[5,7,0,3], # pointC > [pointA,pointB,pointC,pointD]
[5,4,3,0], # pointD > [pointA,pointB,pointC,pointD]
]
#Définition variable
x = [[pp.LpVariable("x(%s,%s)"%(i, j), cat="Binary") for i in range(N)] for j in range(N)]
u = [pp.LpVariable("u(%s)"%(i), cat="Continuous", lowBound=1.0, upBound=(N)) for i in range(N)]
#Définition et enregistrement de l'indice d'évaluation (équation (1))
objective = pp.lpSum(c_[i][j] * x[i][j] for i in range(N) for j in range(N) if i != j)
mip_model += objective
#Expression conditionnelle(2)Enregistrement de
for i in range(N):
mip_model += pp.lpSum(x[i][j] for j in range(N) if i != j) == 1
#Expression conditionnelle(3)Enregistrement de
for i in range(N):
mip_model += pp.lpSum(x[j][i] for j in range(N) if i != j) == 1
for i in range(N):
mip_model += x[i][i] == 0
#Expression conditionnelle(4) (Contrainte MTZ)
for i in range(N):
for j in range(N):
if i != j:
mip_model += u[i] + 1.0 - BigM * (1.0 - x[i][j]) <= u[j]
#Effectuer l'optimisation
status = mip_model.solve()
#Comprendre les résultats
print("Status: {}".format(pp.LpStatus[status]))
print("Optimal Value [a.u.]: {}".format(objective.value()))
for i in range(N):
for j in range(N):
if i != j:
print("x[%d][%d]:%f" % (i,j,x[i][j].value()))
for i in range(len(u)):
print("u[%d] %f" % (i,u[i].value()))
Si vous l'exécutez dans cet état ...
Status: Infeasible
Optimal Value [a.u.]: 18.0
x[0][1]:0.000000
x[0][2]:1.000000
x[0][3]:0.000000
x[1][0]:0.999800
x[1][2]:0.000000
x[1][3]:0.000200
x[2][0]:0.000200
x[2][1]:0.000000
x[2][3]:0.999800
x[3][0]:0.000000
x[3][1]:1.000000
x[3][2]:0.000000
Il semble qu'il n'ait pas été bien calculé. Il semble que les conditions de contrainte ne soient pas bien satisfaites, donc après divers examens, il semble que la formule de contrainte MTZ soit étrange en l'état. Si vous écrivez les quatre formules pour x à u dans l'état final, ce sera comme suit.
u[0] < u[2]\\
u[1] < u[0]\\
u[2] < u[3]\\
u[3] < u[1]
Quand il s'agit de satisfaire tout cela,
u[0] < u[2] < u[3] < u[1] < u[0] < ...
Il augmentera à l'infini, mais comme la plage effective de u est comprise entre 1 et 4, il semble qu'il y ait une contradiction. Par conséquent, en supprimant la restriction sur la variable qui retourne au point A à la fin, reconstruisons la formule comme suit.
u[0] < u[2]\\
u[2] < u[3]\\
u[3] < u[1]
Cela modifiera l'implémentation et changera la formule de condition MTZ en:
#Expression conditionnelle(4) (Contrainte MTZ)
for i in range(N):
for j in range(1,N):
if i != j:
mip_model += u[i] + 1.0 - BigM * (1.0 - x[i][j]) <= u[j]
Le code source est résumé ci-dessous.
sample_route.py
import pulp as pp
#Définition du modèle d'optimisation
mip_model = pp.LpProblem("tsp_mip", pp.LpMinimize)
N = 4
BigM = 10000
c_ = [
[0,6,5,5], # pointA > [pointA,pointB,pointC,pointD]
[6,0,7,4], # pointB > [pointA,pointB,pointC,pointD]
[5,7,0,3], # pointC > [pointA,pointB,pointC,pointD]
[5,4,3,0], # pointD > [pointA,pointB,pointC,pointD]
]
#Définition variable
x = [[pp.LpVariable("x(%s,%s)"%(i, j), cat="Binary") for i in range(N)] for j in range(N)]
u = [pp.LpVariable("u(%s)"%(i), cat="Continuous", lowBound=1.0, upBound=(N)) for i in range(N)]
#Définition et enregistrement de l'indice d'évaluation (équation (1))
objective = pp.lpSum(c_[i][j] * x[i][j] for i in range(N) for j in range(N) if i != j)
mip_model += objective
#Expression conditionnelle(2)Enregistrement de
for i in range(N):
mip_model += pp.lpSum(x[i][j] for j in range(N) if i != j) == 1
#Expression conditionnelle(3)Enregistrement de
for i in range(N):
mip_model += pp.lpSum(x[j][i] for j in range(N) if i != j) == 1
for i in range(N):
mip_model += x[i][i] == 0
#Expression conditionnelle(4) (Contrainte MTZ)
for i in range(N):
for j in range(1,N):
if i != j:
mip_model += u[i] + 1.0 - BigM * (1.0 - x[i][j]) <= u[j]
#Effectuer l'optimisation
status = mip_model.solve()
#Comprendre les résultats
print("Status: {}".format(pp.LpStatus[status]))
print("Optimal Value [a.u.]: {}".format(objective.value()))
for i in range(N):
for j in range(N):
if i != j:
print("x[%d][%d]:%f" % (i,j,x[i][j].value()))
for i in range(len(u)):
print("u[%d] %f" % (i,u[i].value()))
Le résultat est le suivant.
Status: Optimal
Optimal Value [a.u.]: 18.0
x[0][1]:0.000000
x[0][2]:1.000000
x[0][3]:0.000000
x[1][0]:1.000000
x[1][2]:0.000000
x[1][3]:0.000000
x[2][0]:0.000000
x[2][1]:0.000000
x[2][3]:1.000000
x[3][0]:0.000000
x[3][1]:1.000000
x[3][2]:0.000000
u[0] 1.000000
u[1] 4.000000
u[2] 2.000000
u[3] 3.000000
Il a convergé en toute sécurité. Les résultats ne sont pas graphiques,
Point A → Point C → Point D → Point B → Point A
Semble être le plus court et la distance est de 18. Où le problème a été répertorié
Point A → Point B → Point D → Point C → Point A
Donc, après tout, la distance est de 18 et elle semble être la plus courte, donc il y a quelques réponses.
L'optimisation mathématique n'est très récente que dans des domaines que nous n'avons pas encore abordés. Il existe un environnement où vous pouvez vous sentir libre de l'essayer, alors je vais profiter de cette occasion pour étudier diverses choses.
Recommended Posts