PuLP

Introduction

In this example, we will explore how to submit a PuLP program to Multivac in order to solve a TSP-type problem.

PuLP is a Python optimization library that provides a high-level interface to solve linear programming, mixed-integer programming, and quadratic programming problems. It can use several solvers, such as GLPK, COIN, CPLEX, Gurobi, and others, to find optimal solutions for complex optimization problems.

The travelling salesman problem (TSP) is a classic combinatorial optimization problem where the goal is to find the shortest route that visits a given set of cities exactly once and returns to the starting point. It is an NP-hard problem, which means there are no known algorithms that can solve large instances efficiently in terms of execution time. However, approximate or heuristic solutions can be found for large instances. The TSP has many applications in logistics, transportation, circuit design, and more.

TSP Example Code

In this example we have 20 nodes with their corresponding coordinates:

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

coordinates = [[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]]

We want to obtain an optimal route that visits each node exactly once:

_images/TSP.jpg

We will create a file called tsp.py to compute the result:

#!/usr/bin/python3
import pulp as pl
import sys
# Ask GLPK which solvers are available
solvers = pl.listSolvers(onlyAvailable=True)
print('Solvers', solvers)
print('We will use 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']

# Create the problem
prob = pl.LpProblem("TSP", pl.LpMinimize)

# Parameters
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}

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

# Constraints
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

# Solve the problem
# You can use any solver here: ['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']
# You can also configure gapRel, time limit, and mip = 1 for mixed-integer programming

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

# Check problem status
if resultat == pl.LpStatusOptimal:
   print("An optimal solution was found.")
   # Print solution
   print("Optimal route:")
   for i in id_ciutat:
      for j in id_ciutat:
            if i != j and x[i, j].value() == 1:
               print(f"Point {i} -> Point {j}")
elif resultat == pl.LpStatusInfeasible:
   print("The problem is infeasible.")
elif resultat == pl.LpStatusUnbounded:
   print("The problem is unbounded.")
else:
   print("No optimal solution could be found.")

Sudoku Example Code

In this example we solve the well-known Sudoku puzzle, a classic mathematical challenge where we must fill a 9x9 grid with numbers from 1 to 9, without repeated numbers in any row, column, or 3x3 subgrid.

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

prob = pl.LpProblem("Sudoku")

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

# Constraints

# Constraint to ensure that each cell contains exactly one number
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

# Constraints to ensure each row contains different numbers
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

# Constraints to ensure each column contains different numbers
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

# Constraints to ensure each 3x3 subgrid contains different numbers
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

# Add initial Sudoku values
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

# Solve the problem
resultat = prob.solve(pl.PULP_CBC_CMD(mip=1, timeLimit=600, gapRel=0.001))


# Check problem status
if resultat == pl.LpStatusOptimal:
   print("An optimal solution was found.")
   # Print solution
   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("The problem is infeasible.")
elif resultat == pl.LpStatusUnbounded:
   print("The problem is unbounded.")
else:
   print("No optimal solution could be found.")

Execution

We will create a Python virtual environment named venv. This step is only required the first time you create the venv; afterwards you can reuse it. You need to install the PuLP package: pip install pulp

If you want to install all solvers, install the following packages: 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

We recommend using the mvac_crear_venv script if you want to use docplex.

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

We will create a configuration file, for example tsp.slurm, containing the execution settings for this PuLP script.

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" # To know which machine executed the job
   "source $ROUTE/venv/bin/activate"
   "python3 $ROUTE/tsp.py"
   "deactivate"

)

We will submit this script from iocex using the following command:

multivac tsp.slurm

Once execution is complete, the output of our program will be visible in the same directory where the job was launched.

Providing an initial solution for the solver

For more information, visit the PuLP website.