Pulp

Introducció

En aquest exemple, explorarem com enviar a Multivac un programa usant Pulp per a resoldre un problema de tipus TSP.

Pulp és una biblioteca d’optimització de Python que proporciona una interfície d’alt nivell per resoldre problemes de programació lineal, mixta-entera i de programació quadràtica. Utilitza diversos solucionadors, com ara GLPK, COIN, CPLEX, Gurobi, etc., per trobar solucions òptimes a problemes complexos d’optimització.

The travelling salesman problem (TSP), problema del Venedor Viatjant és un problema clàssic d’optimització combinatorial en què es busca trobar la ruta més curta que passi per un conjunt donat de ciutats exactament una vegada i torni al punt de partida. És un problema NP-hard, el que significa que no hi ha algoritmes coneguts que puguin resoldre instàncies grans de manera eficient en termes de temps d’execució. No obstant això, es poden trobar solucions aproximades o heurístiques per a instàncies grans. El TSP té moltes aplicacions en logística, transport, disseny de circuits, entre d’altres.

Codi d’exemple TSP

En aquest exemple tenim 20 nodes amb les seves coordenades corresponents:

id_ciutat = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20}

coordenades = [[2.3, 886.2],[832.3, 491.9],[646.3, 731.1],[323.9, 109.9],[934.2, 38.3],[605.2, 166.4],[401.7, 176.9],[710.3, 112.9],[390.3, 499.4],[565.9, 474.9],[43.7, 758.4],[603.3, 46.4],[745.6, 758.9],[336.4, 478.1],[150.8, 878.1],[831.0, 993.9],[701.8, 796.7],[239.2, 821.5],[534.2, 27.8],[608.4, 809.0]]

Tal que ves vulgui obtenir una ruta òptima passant només una vegada per cada node:

Crearem un arxiu per exemple anomenat tsp.py per a calcular el resultat:

#!/usr/bin/python3
import pulp as pl
import sys
# Preguntem a GPLK quins solvers té disponibles
solvers = pl.listSolvers(onlyAvailable=True)
print('Solvers', solvers)
print('Usarem CPLEX_PY')
# ['GLPK_CMD', 'PYGLPK', 'CPLEX_CMD', 'CPLEX_PY', 'CPLEX_DLL', 'GUROBI', 'GUROBI_CMD', 'MOSEK', 'XPRESS', 'PULP_CBC_CMD', 'COIN_CMD', 'COINMP_DLL', 'CHOCO_CMD', 'MIPCL_CMD', 'SCIP_CMD']

# Crear el problema
prob = pl.LpProblem("TSP", pl.LpMinimize)

# Paràmetres
id_ciutat = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20}
n_ciutats = len(id_ciutat)
coordenades = [
   [2.3, 886.2], [832.3, 491.9], [646.3, 731.1], [323.9, 109.9],
   [934.2, 38.3], [605.2, 166.4], [401.7, 176.9], [710.3, 112.9],
   [390.3, 499.4], [565.9, 474.9], [43.7, 758.4], [603.3, 46.4],
   [745.6, 758.9], [336.4, 478.1], [150.8, 878.1], [831.0, 993.9],
   [701.8, 796.7], [239.2, 821.5], [534.2, 27.8], [608.4, 809.0],
]
distancies = [[0.0] * n_ciutats for _ in range(n_ciutats)]
for i in id_ciutat:
   for j in id_ciutat:
      distancies[i - 1][j - 1] = ((coordenades[i - 1][0] - coordenades[j - 1][0]) ** 2 + (coordenades[i - 1][1] - coordenades[j - 1][1]) ** 2) ** 0.5

# Variables
x = {(i, j): pl.LpVariable(f'x_{i}_{j}', cat=pl.LpBinary) for i in id_ciutat for j in id_ciutat}
u = {i: pl.LpVariable(f'u_{i}', lowBound=0) for i in id_ciutat}

# Funció objectiu
prob += pl.lpSum(distancies[i - 1][j - 1] * x[i, j] for i in id_ciutat for j in id_ciutat if i != j)

# Restriccions
for j in id_ciutat:
   prob += pl.lpSum(x[i, j] for i in id_ciutat if i != j) == 1

for i in id_ciutat:
   prob += pl.lpSum(x[i, j] for j in id_ciutat if i != j) == 1

for i in range(2, n_ciutats + 1):
   for j in range(2, n_ciutats + 1):
      if i != j:
            prob += u[i] - u[j] + n_ciutats * x[i, j] <= n_ciutats - 1

# Resoldre el problema
# Aquí podem posar qualsevol solver: ['GLPK_CMD', 'PYGLPK', 'CPLEX_CMD',
# 'CPLEX_PY', 'CPLEX_DLL', 'GUROBI', 'GUROBI_CMD', 'MOSEK', 'XPRESS',
# 'PULP_CBC_CMD', 'COIN_CMD', 'COINMP_DLL', 'CHOCO_CMD', 'MIPCL_CMD', 'SCIP_CMD']
# A més podem configurar el gapRel, el límit en temps i mip = 1, que és de programació mixta entera

resultat = prob.solve(pl.CPLEX_PY(mip=1, timeLimit=600, gapRel=0.001)) # Més info a: https://coin-or.github.io/pulp/guides/how_to_configure_solvers.html

# Comprovar l'estat del problema
if resultat == pl.LpStatusOptimal:
   print("S'ha trobat una solució òptima.")
   # Imprimir la solució
   print("Ruta òptima:")
   for i in id_ciutat:
      for j in id_ciutat:
            if i != j and x[i, j].value() == 1:
               print(f"Punt {i} -> Punt {j}")
elif resultat == pl.LpStatusInfeasible:
   print("El problema és infactible.")
elif resultat == pl.LpStatusUnbounded:
   print("El problema no està acotat.")
else:
   print("No s'ha pogut trobar una solució òptima.")

Codi d’exemple Sudoku

En aquest exemple tenim el conegut Sudoku, clàssic repte matemàtic en el què hem d’emplenar una graella de 9x9 amb nombres de l’1 al 9 en la què no hi poden haver nombres repetits a la fila, columna o subgrups de 3x3.

#!/usr/bin/python3
import pulp as pl

prob = pl.LpProblem("Sudoku")

# Definir les variables
choices = pl.LpVariable.dicts("Choice", (range(1, 10), range(1, 10), range(1, 10)), cat='Binary')

# Restriccions

# Restricció per assegurar que cada cel·la contingui exactament un número
for i in range(1, 10):
   for j in range(1, 10):
      prob += pl.lpSum(choices[i][j][k] for k in range(1, 10)) == 1

# Restriccions per assegurar que cada fila contingui números diferents
for i in range(1, 10):
   for k in range(1, 10):
      prob += pl.lpSum(choices[i][j][k] for j in range(1, 10)) == 1

# Restriccions per assegurar que cada columna contingui números diferents
for j in range(1, 10):
   for k in range(1, 10):
      prob += pl.lpSum(choices[i][j][k] for i in range(1, 10)) == 1

# Restriccions per assegurar que cada subgraella 3x3 contingui números diferents
for k in range(1, 10):
   for i in range(1, 10, 3):
      for j in range(1, 10, 3):
            prob += pl.lpSum(choices[i+di][j+dj][k] for di in range(3) for dj in range(3)) == 1

# Afegir les posicions inicials del Sudoku
initial_values = {
   (1, 1, 5), (1, 2, 3), (1, 5, 7), (2, 1, 6),
   (2, 4, 1), (2, 5, 9), (2, 6, 5), (3, 2, 9),
   (3, 3, 8), (3, 8, 6), (4, 1, 8), (4, 5, 6),
   (4, 9, 3), (5, 1, 4), (5, 4, 8), (5, 6, 3),
   (5, 9, 1), (6, 1, 7), (6, 5, 2), (6, 9, 6),
   (7, 2, 6), (7, 7, 2), (7, 8, 8), (8, 4, 4),
   (8, 5, 1), (8, 6, 9), (8, 9, 5), (9, 5, 8),
   (9, 8, 7), (9, 9, 9)
}

for (i, j, k) in initial_values:
   prob += choices[i][j][k] == 1

# Resoldre el problema
resultat = prob.solve(pl.PULP_CBC_CMD(mip=1, timeLimit=600, gapRel=0.001))


# Comprovar l'estat del problema
if resultat == pl.LpStatusOptimal:
   print("S'ha trobat una solució òptima.")
   # Imprimir la solució
   for i in range(1, 10):
      for j in range(1, 10):
            for k in range(1, 10):
               if pl.value(choices[i][j][k]) == 1:
                  print(k, end=" ")
                  break
      print()

elif resultat == pl.LpStatusInfeasible:
   print("El problema és infactible.")
elif resultat == pl.LpStatusUnbounded:
   print("El problema no està acotat.")
else:
   print("No s'ha pogut trobar una solució òptima.")

Execució

Crearem un virtual environment de python anomenat venv. Aquest pas només és necessari fer la primera vegada que es vol crear el venv, les altres vegades es pot reusar el mateix. Es requereix instal·lar el paquet pulp: pip install pulp

Si volem instal·lar tots els solvers hem d’instal·lar els següents paquets: pip install docplex pip install cylp pip install glpk pip install gurobipy pip install mosek pip install xpress pip install choco-solver pip install mipcl-py

Recomanem usar l’script mvac_crear_venv en cas de voler usar docplex.

$ python3 -m venv venv
$ source venv/bin/activate
$ pip install pulp
$ deactivate

Crearem un document amb els parametres de configuració, per exemple: tsp.slurm` que conté la configuració d’execució per a aquest script de pulp.

VERSION=1.3
JOB_NAME=pulp_example
NAME_OUTPUT=out
PARTITION=all
N_TASKS=1
CPUS_PER_TASK=1
MAIL_TYPE=END,FAIL
MAIL_USER=nom.usuari@upc.edu
MEMORY=1G
BEGIN=now
TIME_LIMIT=00:05:00
LOG_OUTPUT=log
FORCED_NODES=
EXCLUDED_NODES=
ROUTE=~/
COMMANDS=(
   "hostname" #Per a saber a quina màquina s'ha executat
   "source $ROUTE/venv/bin/activate"
   "python3 $ROUTE/tsp.py"
   "deactivate"

)

Aquest script el llançarem des de iocex usant la següent comanda:

multivac tsp.slurm

El resultat del nostre programa un cop finalitzada l’execució serà visible al mateix directori a on hem fet l’execució.

Donar solució inicial al solver d’on començar

Per a més informació, visita el lloc web de Pulp.