La programación de restricciones es una técnica de elección para resolver un problema de satisfacción de restricciones. En este artículo veremos que también se adapta bien a problemas de optimización pequeños y medianos. Usando el conocido problema del vendedor ambulante (TSP) como ejemplo, detallaremos todos los pasos que conducen a un modelo eficiente.
Por simplicidad, consideraremos el caso simétrico del TSP (la distancia entre dos ciudades es la misma en cada dirección opuesta).
Todos los ejemplos de código de este artículo utilizan Nucsun solucionador rápido de restricciones escrito 100% en Python que actualmente estoy desarrollando como proyecto paralelo. NuCS se publica bajo el licencia MIT.
Citando Wikipedia:
Dada una lista de ciudades y las distancias entre cada par de ciudades, ¿cuál es la ruta más corta posible que visita cada ciudad exactamente una vez y regresa a la ciudad de origen?
Este es un problema NP-difícil. De ahora en adelante, consideremos que existen norte ciudades.
La formulación más ingenua de este problema es decidir, para cada posible borde entre ciudades, si pertenece a la solución óptima. El tamaño del espacio de búsqueda es 2ⁿ⁽ⁿ⁻¹⁾ᐟ² que es aproximadamente 8.8e130 para n=30 (mucho mayor que el número de átomos en el universo).
Es mucho mejor encontrar, para cada ciudad, su sucesora. La complejidad se vuelve ¡norte! que es aproximadamente 2.6e32 para n=30 (mucho más pequeño pero aún muy grande).
A continuación, compararemos nuestros modelos con las siguientes instancias pequeñas de TSP: GR17, GR21 y GR24.
GR17 es un 17 nodos TSP simétricos, sus costos están definidos por 17 x 17 matriz simétrica de costos sucesores:
[
[0, 633, 257, 91, 412, 150, 80, 134, 259, 505, 353, 324, 70, 211, 268, 246, 121],
[633, 0, 390, 661, 227, 488, 572, 530, 555, 289, 282, 638, 567, 466, 420, 745, 518],
[257, 390, 0, 228, 169, 112, 196, 154, 372, 262, 110, 437, 191, 74, 53, 472, 142],
[91, 661, 228, 0, 383, 120, 77, 105, 175, 476, 324, 240, 27, 182, 239, 237, 84],
[412, 227, 169, 383, 0, 267, 351, 309, 338, 196, 61, 421, 346, 243, 199, 528, 297],
[150, 488, 112, 120, 267, 0, 63, 34, 264, 360, 208, 329, 83, 105, 123, 364, 35],
[80, 572, 196, 77, 351, 63, 0, 29, 232, 444, 292, 297, 47, 150, 207, 332, 29],
[134, 530, 154, 105, 309, 34, 29, 0, 249, 402, 250, 314, 68, 108, 165, 349, 36],
[259, 555, 372, 175, 338, 264, 232, 249, 0, 495, 352, 95, 189, 326, 383, 202, 236],
[505, 289, 262, 476, 196, 360, 444, 402, 495, 0, 154, 578, 439, 336, 240, 685, 390],
[353, 282, 110, 324, 61, 208, 292, 250, 352, 154, 0, 435, 287, 184, 140, 542, 238],
[324, 638, 437, 240, 421, 329, 297, 314, 95, 578, 435, 0, 254, 391, 448, 157, 301],
[70, 567, 191, 27, 346, 83, 47, 68, 189, 439, 287, 254, 0, 145, 202, 289, 55],
[211, 466, 74, 182, 243, 105, 150, 108, 326, 336, 184, 391, 145, 0, 57, 426, 96],
[268, 420, 53, 239, 199, 123, 207, 165, 383, 240, 140, 448, 202, 57, 0, 483, 153],
[246, 745, 472, 237, 528, 364, 332, 349, 202, 685, 542, 157, 289, 426, 483, 0, 336],
[121, 518, 142, 84, 297, 35, 29, 36, 236, 390, 238, 301, 55, 96, 153, 336, 0],
]
Echemos un vistazo a la primera fila:
[0, 633, 257, 91, 412, 150, 80, 134, 259, 505, 353, 324, 70, 211, 268, 246, 121]
Estos son los costos para los posibles sucesores de nodo 0 en el circuito. Si exceptuamos el primer valor 0 (no queremos el sucesor del nodo 0 ser nodo 0) entonces el valor mínimo es 70 (cuando el nodo 12 es el sucesor del nodo 0) y el máximo es 633 (cuando el nodo 1 es el sucesor del nodo 0). Esto significa que el costo asociado al sucesor del nodo 0 en el circuito oscila entre 70 y 633.
Vamos a modelar nuestro problema reutilizando el Problema de circuito disponible en el mercado en NuCS. Pero primero comprendamos qué sucede detrás de escena. El Problema de circuito es en sí misma una subclase de la Permutación problema, otro modelo disponible ofrecido por NuCS.
El problema de la permutación
El problema de permutación define dos modelos redundantes: los modelos sucesores y predecesores.
def __init__(self, n: int):
"""
Inits the permutation problem.
:param n: the number variables/values
"""
self.n = n
shr_domains = [(0, n - 1)] * 2 * n
super().__init__(shr_domains)
self.add_propagator((list(range(n)), ALG_ALLDIFFERENT, []))
self.add_propagator((list(range(n, 2 * n)), ALG_ALLDIFFERENT, []))
for i in range(n):
self.add_propagator((list(range(n)) + [n + i], ALG_PERMUTATION_AUX, [i]))
self.add_propagator((list(range(n, 2 * n)) + [i], ALG_PERMUTATION_AUX, [i]))
El modelo de sucesores (las primeras n variables) define, para cada nodo, su sucesor en el circuito. Los sucesores tienen que ser diferentes. El modelo de predecesores (las últimas n variables) define, para cada nodo, su predecesor en el circuito. Los predecesores tienen que ser diferentes.
Ambos modelos están relacionados con las reglas (ver la ALG_PERMUTATION_AUX restricciones):
- si éxito[i] =j entonces pred[j] = yo
- si pred[j] = yo entonces éxito[i] =j
- si pred[j] ≠ yo entonces éxito[i] ≠j
- si éxito[i] ≠j entonces pred[j] ≠ yo
El problema del circuito
El problema del circuito refina los dominios de los sucesores y predecesores y agrega restricciones adicionales para prohibir subciclos (no entraremos en ellos aquí por razones de brevedad).
def __init__(self, n: int):
"""
Inits the circuit problem.
:param n: the number of vertices
"""
self.n = n
super().__init__(n)
self.shr_domains_lst[0] = [1, n - 1]
self.shr_domains_lst[n - 1] = [0, n - 2]
self.shr_domains_lst[n] = [1, n - 1]
self.shr_domains_lst[2 * n - 1] = [0, n - 2]
self.add_propagator((list(range(n)), ALG_NO_SUB_CYCLE, []))
self.add_propagator((list(range(n, 2 * n)), ALG_NO_SUB_CYCLE, []))
El modelo TSP
Con la ayuda del problema del circuito, modelar el TSP es una tarea sencilla.
Consideremos un nodo icomo se vio antes costos[i] es la lista de posibles costos para los sucesores de i. Si j es el sucesor de i entonces el costo asociado es costos[i]ⱼ. Esto se implementa mediante la siguiente línea donde costos_succ si el índice inicial de los sucesores cuesta:
self.add_propagators([([i, self.succ_costs + i], ALG_ELEMENT_IV, costs[i]) for i in range(n)])
Simétricamente, para los costos predecesores obtenemos:
self.add_propagators([([n + i, self.pred_costs + i], ALG_ELEMENT_IV, costs[i]) for i in range(n)])
Finalmente, podemos definir el costo total sumando los costos intermedios y obtenemos:
def __init__(self, costs: List[List[int]]) -> None:
"""
Inits the problem.
:param costs: the costs between vertices as a list of lists of integers
"""
n = len(costs)
super().__init__(n)
max_costs = [max(cost_row) for cost_row in costs]
min_costs = [min([cost for cost in cost_row if cost > 0]) for cost_row in costs]
self.succ_costs = self.add_variables([(min_costs[i], max_costs[i]) for i in range(n)])
self.pred_costs = self.add_variables([(min_costs[i], max_costs[i]) for i in range(n)])
self.total_cost = self.add_variable((sum(min_costs), sum(max_costs))) # the total cost
self.add_propagators([([i, self.succ_costs + i], ALG_ELEMENT_IV, costs[i]) for i in range(n)])
self.add_propagators([([n + i, self.pred_costs + i], ALG_ELEMENT_IV, costs[i]) for i in range(n)])
self.add_propagator(
(list(range(self.succ_costs, self.succ_costs + n)) + [self.total_cost], ALG_AFFINE_EQ, [1] * n + [-1, 0])
)
self.add_propagator(
(list(range(self.pred_costs, self.pred_costs + n)) + [self.total_cost], ALG_AFFINE_EQ, [1] * n + [-1, 0])
)
Tenga en cuenta que no es necesario tener modelos sucesores y predecesores (uno sería suficiente), pero es más eficiente.
Usemos la estrategia de ramificación predeterminada del RetrocederSolvernuestras variables de decisión serán los sucesores y predecesores.
solver = BacktrackSolver(problem, decision_domains=decision_domains)
solution = solver.minimize(problem.total_cost)
La solución óptima se encuentra en 248 en una MacBook Pro M2 con Python 3.12, Numpy 2.0.1, Numba 0.60.0 y NuCS 4.2.0. Las estadísticas detalladas proporcionadas por NuCS son:
{
'ALG_BC_NB': 16141979,
'ALG_BC_WITH_SHAVING_NB': 0,
'ALG_SHAVING_NB': 0,
'ALG_SHAVING_CHANGE_NB': 0,
'ALG_SHAVING_NO_CHANGE_NB': 0,
'PROPAGATOR_ENTAILMENT_NB': 136986225,
'PROPAGATOR_FILTER_NB': 913725313,
'PROPAGATOR_FILTER_NO_CHANGE_NB': 510038945,
'PROPAGATOR_INCONSISTENCY_NB': 8070394,
'SOLVER_BACKTRACK_NB': 8070393,
'SOLVER_CHOICE_NB': 8071487,
'SOLVER_CHOICE_DEPTH': 15,
'SOLVER_SOLUTION_NB': 98
}
En particular, hay 8.070.393 retrocesos. Intentemos mejorar esto.
NuCS ofrece una heurística basada en el arrepentimiento (diferencia entre el mejor y el segundo mejor costo) para seleccionar la variable. Luego elegiremos el valor que minimice el coste.
solver = BacktrackSolver(
problem,
decision_domains=decision_domains,
var_heuristic_idx=VAR_HEURISTIC_MAX_REGRET,
var_heuristic_params=costs,
dom_heuristic_idx=DOM_HEURISTIC_MIN_COST,
dom_heuristic_params=costs
)
solution = solver.minimize(problem.total_cost)
Con estas nuevas heurísticas, la solución óptima se encuentra en 38s y las estadisticas son:
{
'ALG_BC_NB': 2673045,
'ALG_BC_WITH_SHAVING_NB': 0,
'ALG_SHAVING_NB': 0,
'ALG_SHAVING_CHANGE_NB': 0,
'ALG_SHAVING_NO_CHANGE_NB': 0,
'PROPAGATOR_ENTAILMENT_NB': 12295905,
'PROPAGATOR_FILTER_NB': 125363225,
'PROPAGATOR_FILTER_NO_CHANGE_NB': 69928021,
'PROPAGATOR_INCONSISTENCY_NB': 1647125,
'SOLVER_BACKTRACK_NB': 1647124,
'SOLVER_CHOICE_NB': 1025875,
'SOLVER_CHOICE_DEPTH': 36,
'SOLVER_SOLUTION_NB': 45
}
En particular, hay 1.647.124 retrocesos.
Podemos seguir mejorando diseñando un heurística personalizada que combina el arrepentimiento máximo y el dominio más pequeño para la selección de variables.
tsp_var_heuristic_idx = register_var_heuristic(tsp_var_heuristic)
solver = BacktrackSolver(
problem,
decision_domains=decision_domains,
var_heuristic_idx=tsp_var_heuristic_idx,
var_heuristic_params=costs,
dom_heuristic_idx=DOM_HEURISTIC_MIN_COST,
dom_heuristic_params=costs
)
solution = solver.minimize(problem.total_cost)
La solución óptima ahora se encuentra en 11 y las estadisticas son:
{
'ALG_BC_NB': 660718,
'ALG_BC_WITH_SHAVING_NB': 0,
'ALG_SHAVING_NB': 0,
'ALG_SHAVING_CHANGE_NB': 0,
'ALG_SHAVING_NO_CHANGE_NB': 0,
'PROPAGATOR_ENTAILMENT_NB': 3596146,
'PROPAGATOR_FILTER_NB': 36847171,
'PROPAGATOR_FILTER_NO_CHANGE_NB': 20776276,
'PROPAGATOR_INCONSISTENCY_NB': 403024,
'SOLVER_BACKTRACK_NB': 403023,
'SOLVER_CHOICE_NB': 257642,
'SOLVER_CHOICE_DEPTH': 33,
'SOLVER_SOLUTION_NB': 52
}
En particular, hay 403.023 retrocesos.
La minimización (y más generalmente la optimización) se basa en una ramificarse y encuadernarse algoritmo. El mecanismo de retroceso permite explorar el espacio de búsqueda tomando decisiones (derivación). Partes del espacio de búsqueda son podadas por delimitando la variable objetivo.
Al minimizar una variable tse puede agregar la restricción adicional t < s siempre que se encuentre una solución intermedia s.
NuCS ofrece dos modos de optimización correspondientes a dos formas de aprovechar t < s:
- el REINICIAR El modo reinicia la búsqueda desde cero y actualiza los límites de la variable de destino.
- el CIRUELA PASA El modo modifica los puntos de elección para tener en cuenta los nuevos límites de la variable objetivo.
Probemos ahora el CIRUELA PASA modo:
solution = solver.minimize(problem.total_cost, mode=PRUNE)
La solución óptima se encuentra en 5,4s y las estadisticas son:
{
'ALG_BC_NB': 255824,
'ALG_BC_WITH_SHAVING_NB': 0,
'ALG_SHAVING_NB': 0,
'ALG_SHAVING_CHANGE_NB': 0,
'ALG_SHAVING_NO_CHANGE_NB': 0,
'PROPAGATOR_ENTAILMENT_NB': 1435607,
'PROPAGATOR_FILTER_NB': 14624422,
'PROPAGATOR_FILTER_NO_CHANGE_NB': 8236378,
'PROPAGATOR_INCONSISTENCY_NB': 156628,
'SOLVER_BACKTRACK_NB': 156627,
'SOLVER_CHOICE_NB': 99143,
'SOLVER_CHOICE_DEPTH': 34,
'SOLVER_SOLUTION_NB': 53
}
En particular, sólo hay 156.627 retrocesos.
La siguiente tabla resume nuestros experimentos:
Puedes encontrar todo el código correspondiente. aquí.
Por supuesto, hay muchas otras vías que podríamos explorar para mejorar estos resultados:
- diseñar un restricción redundante por el costo total
- mejorar la ramificación explorando nuevas heurísticas
- use un algoritmo de consistencia diferente (NuCS viene con afeitado)
- calcular los límites inferior y superior utilizando otras técnicas
El problema del viajante ha sido objeto de extensos estudios y abundante literatura. En este artículo esperamos haber convencido al lector de que es posible encontrar soluciones óptimas a problemas de tamaño mediano en muy poco tiempo, sin tener muchos conocimientos sobre el problema del viajante.