Programación Lineal Avanzada con Python: Branch & Bound, CVXPY y TSP

Hemos recorrido un largo camino: en el nivel básico aprendiste a formular y resolver PL con SciPy y PuLP, y en el nivel intermedio dominaste la Programación Entera Mixta, los problemas clásicos y el análisis de sensibilidad. En este artículo avanzado abriremos la caja negra de los solvers, trabajaremos con solvers de nivel industrial, resolveremos el Problema del Viajante (TSP), usaremos CVXPY para programación cuadrática y exploraremos cómo escalar la optimización a problemas reales con miles de variables.

Branch and Bound: La Mecánica Interna de los Solvers MIP

Antes de usar solvers avanzados como Gurobi o HiGHS, es fundamental entender el algoritmo que los impulsa: Branch and Bound (B&B). Comprender B&B te ayuda a modelar mejor, establecer cotas útiles y depurar modelos que no convergen.

¿Cómo funciona Branch and Bound?

El algoritmo trabaja en tres fases iterativas:

1. Relajación LP (Bounding): Se resuelve el problema ignorando la restricción de integralidad. Esto da una cota superior (para maximización) del valor óptimo entero posible.
2. Branching: Si la solución relajada no es entera, se elige una variable fraccionaria \( x_j = 3.7 \) y se crean dos subproblemas: uno con \( x_j \leq 3 \) y otro con \( x_j \geq 4 \).
3. Poda (Pruning): Se descartan ramas donde la cota superior es menor o igual a la mejor solución entera encontrada hasta el momento (cota inferior). Esto evita explorar millones de subproblemas innecesarios.

El proceso genera un árbol de búsqueda. La calidad de la formulación MIP determina qué tan compacto es ese árbol:

Un gap MIP pequeño indica que la relajación LP es una buena aproximación del problema entero, lo que acelera enormemente la resolución.

Tip: Si tu solver tarda mucho, el problema casi siempre es un gap MIP grande. Mejorar la formulación (añadir restricciones válidas, ajustar cotas Big-M, usar presolve) reduce el gap y puede hacer la diferencia entre segundos y horas de cómputo.

Visualizando Branch and Bound en Python

import pulp

def branch_and_bound_demo(c, A_ub, b_ub, nivel=0, mejor=[None]):
    """
    Demostracion simplificada del concepto Branch and Bound.
    NO usar en produccion: es solo ilustrativo.
    """
    indent = "  " * nivel
    prob = pulp.LpProblem(f"Nivel_{nivel}", pulp.LpMaximize)
    n = len(c)
    x = [pulp.LpVariable(f"x{i}", lowBound=b_ub[i][0], upBound=b_ub[i][1])
         for i in range(n)]

    prob += pulp.lpSum(c[i] * x[i] for i in range(n))
    for fila, rhs in A_ub:
        prob += pulp.lpSum(fila[i] * x[i] for i in range(n)) <= rhs

    prob.solve(pulp.PULP_CBC_CMD(msg=0))

    if prob.status != 1:
        print(f"{indent}→ INFACTIBLE, poda.")
        return

    z_relax = pulp.value(prob.objective)
    vals = [xi.value() for xi in x]
    print(f"{indent}Nivel {nivel}: z_relajado={z_relax:.2f}, x={[f'{v:.2f}' for v in vals]}")

    # Si la cota superior no mejora la mejor solucion entera, podar
    if mejor[0] is not None and z_relax <= mejor[0]:
        print(f"{indent}→ PODA: z_relajado={z_relax:.2f} <= mejor={mejor[0]:.2f}")
        return

    # Buscar variable fraccionaria
    frac_idx = next((i for i, v in enumerate(vals)
                     if abs(v - round(v)) > 1e-4), None)

    if frac_idx is None:
        # Solucion entera encontrada
        print(f"{indent}→ SOLUCION ENTERA: z={z_relax:.2f}")
        if mejor[0] is None or z_relax > mejor[0]:
            mejor[0] = z_relax
        return

    frac_val = vals[frac_idx]
    print(f"{indent}Branch en x{frac_idx}={frac_val:.2f}: "
          f"rama ≤{int(frac_val)} y rama ≥{int(frac_val)+1}")

Solvers de Nivel Industrial: HiGHS y Gurobi

En producción, PuLP es solo una interfaz: delega la resolución a un solver externo. La elección del solver impacta enormemente el rendimiento.

Usando HiGHS directamente con highspy

pip install highspy

import highspy

# Crear modelo HiGHS directamente (sin PuLP como intermediario)
h = highspy.Highs()
h.silent()  # suprime output

# Problema: max 300*x1 + 500*x2
# s.t. 2*x1 + 5*x2 <= 40
#      3*x1 + 2*x2 <= 30
#      x1, x2 >= 0

# Agregar variables: (costo, lower_bound, upper_bound)
h.addVar(0, 1e30)   # x1 con lb=0, ub=inf
h.addVar(0, 1e30)   # x2

# Definir sentido: maximizar
h.changeColCost(0, -300)  # negamos para maximizar
h.changeColCost(1, -500)

# Agregar restricciones
# Restriccion 1: 2*x1 + 5*x2 <= 40
h.addRow(-1e30, 40, 2, [0, 1], [2.0, 5.0])
# Restriccion 2: 3*x1 + 2*x2 <= 30
h.addRow(-1e30, 30, 2, [0, 1], [3.0, 2.0])

# Resolver
h.run()

# Obtener solucion
info = h.getInfoValue("primal_solution_status")[1]
sol  = h.getSolution()
print(f"x1 (Sillas): {sol.col_value[0]:.2f}")
print(f"x2 (Mesas):  {sol.col_value[1]:.2f}")
print(f"Ganancia:    ${-h.getInfoValue('objective_function_value')[1]:,.2f}")

Usando Gurobi con gurobipy

Gurobi es el solver más rápido del mercado para MIP. Ofrece licencias académicas gratuitas para estudiantes e investigadores.

pip install gurobipy  # requiere licencia activa

import gurobipy as gp
from gurobipy import GRB

# Crear modelo
modelo = gp.Model("Fabrica")
modelo.setParam("OutputFlag", 0)  # silenciar output

# Variables
x1 = modelo.addVar(lb=0, vtype=GRB.CONTINUOUS, name="Sillas")
x2 = modelo.addVar(lb=0, vtype=GRB.CONTINUOUS, name="Mesas")

# Funcion objetivo
modelo.setObjective(300 * x1 + 500 * x2, GRB.MAXIMIZE)

# Restricciones
modelo.addConstr(2 * x1 + 5 * x2 <= 40, "Madera")
modelo.addConstr(3 * x1 + 2 * x2 <= 30, "Trabajo")

# Resolver
modelo.optimize()

if modelo.Status == GRB.OPTIMAL:
    print(f"Sillas: {x1.X:.0f}")
    print(f"Mesas:  {x2.X:.0f}")
    print(f"Ganancia: ${modelo.ObjVal:,.0f}")

    # Precios sombra (valores duales)
    for constr in modelo.getConstrs():
        print(f"  {constr.ConstrName}: precio_sombra={constr.Pi:.2f}, holgura={constr.Slack:.2f}")

Programación Cuadrática con CVXPY: El Modelo de Markowitz

La programación cuadrática (QP) extiende la programación lineal permitiendo una función objetivo cuadrática. El ejemplo más importante en finanzas es el modelo de Markowitz para optimización de portafolios de inversión.

Teoría del Modelo de Markowitz

Dado un portafolio con \( n \) activos, queremos encontrar los pesos \( w_i \) (fracción del capital invertida en cada activo) que minimicen el riesgo (varianza) para un rendimiento esperado dado:

Donde \( \Sigma \) es la matriz de covarianzas y \( \boldsymbol{\mu} \) el vector de rendimientos esperados.

pip install cvxpy numpy pandas yfinance

import numpy as np
import cvxpy as cp
import pandas as pd

# ============================================================
# Modelo de Markowitz con CVXPY
# ============================================================

# Rendimientos historicos simulados (en la practica usar yfinance)
np.random.seed(42)
n_activos  = 5
n_periodos = 252  # un año de datos diarios

tickers = ["AMZN", "MSFT", "GOOG", "TSLA", "NVDA"]

# Simular rendimientos diarios (log-returns)
rendimientos_sim = np.random.multivariate_normal(
    mean  = [0.0008, 0.0006, 0.0007, 0.0012, 0.0010],
    cov   = np.array([
        [0.0004, 0.0002, 0.0001, 0.0000, 0.0001],
        [0.0002, 0.0003, 0.0002, 0.0001, 0.0001],
        [0.0001, 0.0002, 0.0003, 0.0001, 0.0002],
        [0.0000, 0.0001, 0.0001, 0.0009, 0.0003],
        [0.0001, 0.0001, 0.0002, 0.0003, 0.0006],
    ]),
    size  = n_periodos
)

# Estadisticas del portafolio
mu    = rendimientos_sim.mean(axis=0) * 252      # rendimiento anual
Sigma = np.cov(rendimientos_sim.T)   * 252       # covarianza anual

print("Rendimientos anualizados esperados:")
for ticker, r in zip(tickers, mu):
    print(f"  {ticker}: {r*100:.1f}%")

# ============================================================
# Frontera eficiente: resolver para distintos niveles de riesgo
# ============================================================
w = cp.Variable(n_activos)
r_objetivo = cp.Parameter()

# Problema QP: minimizar varianza del portafolio
riesgo    = cp.quad_form(w, Sigma)
problema  = cp.Problem(
    cp.Minimize(riesgo),
    [
        w @ mu >= r_objetivo,   # rendimiento minimo deseado
        cp.sum(w) == 1,         # pesos suman 100%
        w >= 0,                  # sin ventas en corto
    ]
)

# Calcular la frontera eficiente
resultados = []
for r_min in np.linspace(mu.min(), mu.max() * 0.95, 30):
    r_objetivo.value = r_min
    try:
        problema.solve(solver=cp.CLARABEL, warm_start=True)
        if problema.status == "optimal":
            resultados.append({
                "rendimiento": r_min * 100,
                "riesgo_std":  np.sqrt(riesgo.value) * 100,
                "pesos":       w.value.round(4)
            })
    except Exception:
        continue

# Mostrar portafolio de minima varianza
min_var = min(resultados, key=lambda r: r["riesgo_std"])
print(f"\n--- Portafolio de Minima Varianza ---")
print(f"Rendimiento anual: {min_var['rendimiento']:.2f}%")
print(f"Riesgo (std):      {min_var['riesgo_std']:.2f}%")
print(f"Sharpe approx:     {min_var['rendimiento']/min_var['riesgo_std']:.2f}")
print("\nAsignacion:")
for ticker, peso in zip(tickers, min_var["pesos"]):
    if peso > 0.001:
        print(f"  {ticker}: {peso*100:.1f}%")

import matplotlib.pyplot as plt

# Graficar la frontera eficiente
riesgos     = [r["riesgo_std"]    for r in resultados]
rendims     = [r["rendimiento"]   for r in resultados]

fig, ax = plt.subplots(figsize=(9, 6))
ax.plot(riesgos, rendims, 'b-o', markersize=4, linewidth=2, label="Frontera Eficiente")

# Activos individuales
for i, ticker in enumerate(tickers):
    ax.scatter(np.sqrt(Sigma[i, i]) * 100, mu[i] * 100,
               s=100, zorder=5)
    ax.annotate(ticker, (np.sqrt(Sigma[i, i])*100, mu[i]*100),
                textcoords="offset points", xytext=(6, 4), fontsize=9)

ax.set_xlabel("Riesgo — Desviacion Estandar Anual (%)", fontsize=11)
ax.set_ylabel("Rendimiento Esperado Anual (%)", fontsize=11)
ax.set_title("Frontera Eficiente de Markowitz", fontsize=13)
ax.legend(); ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("frontera_eficiente.png", dpi=150)
plt.show()

El Problema del Viajante (TSP) como MIP

El Problema del Viajante (TSP, Travelling Salesman Problem) es uno de los problemas NP-hard más famosos: dado un conjunto de ciudades, encontrar la ruta más corta que visite todas exactamente una vez y regrese al origen.

La formulación MIP más compacta usa las restricciones de Miller-Tucker-Zemlin (MTZ) para eliminar subtours:

import pulp
import numpy as np
import itertools

def resolver_tsp_mtz(ciudades, distancias):
    """
    Resuelve el TSP con formulacion MTZ usando PuLP.
    ciudades: lista de nombres
    distancias: dict {(i,j): distancia}
    """
    n = len(ciudades)
    idx = list(range(n))

    tsp = pulp.LpProblem("TSP_MTZ", pulp.LpMinimize)

    # Variables binarias x[i][j] = 1 si vamos de ciudad i a ciudad j
    x = {(i, j): pulp.LpVariable(f"x_{i}_{j}", cat='Binary')
         for i in idx for j in idx if i != j}

    # Variables auxiliares MTZ para eliminar subtours
    u = {i: pulp.LpVariable(f"u_{i}", lowBound=0, upBound=n - 1)
         for i in idx}

    # Funcion objetivo: minimizar distancia total
    tsp += pulp.lpSum(distancias[i, j] * x[i, j]
                      for i in idx for j in idx if i != j)

    # Restricciones de grado: cada ciudad se visita exactamente una vez
    for i in idx:
        tsp += pulp.lpSum(x[i, j] for j in idx if j != i) == 1  # salida
        tsp += pulp.lpSum(x[j, i] for j in idx if j != i) == 1  # entrada

    # Restricciones MTZ (eliminacion de subtours)
    for i in idx[1:]:
        for j in idx[1:]:
            if i != j:
                tsp += u[i] - u[j] + n * x[i, j] <= n - 1

    # Resolver con tiempo limite
    tsp.solve(pulp.PULP_CBC_CMD(msg=0, timeLimit=60))

    if tsp.status != 1:
        return None, None

    # Reconstruir la ruta
    ruta = [0]
    actual = 0
    while len(ruta) < n:
        siguiente = next(j for j in idx if j != actual
                         and x[actual, j].value() > 0.5)
        ruta.append(siguiente)
        actual = siguiente
    ruta.append(0)  # regresar al origen

    distancia_total = pulp.value(tsp.objective)
    return ruta, distancia_total

# Ejemplo: 8 ciudades con coordenadas
np.random.seed(7)
n_ciudades  = 8
coordenadas = {i: (np.random.uniform(0, 100), np.random.uniform(0, 100))
               for i in range(n_ciudades)}
nombres     = [f"Ciudad_{i}" for i in range(n_ciudades)]

# Calcular distancias euclidianas
distancias = {}
for i in range(n_ciudades):
    for j in range(n_ciudades):
        if i != j:
            xi, yi = coordenadas[i]
            xj, yj = coordenadas[j]
            distancias[i, j] = np.sqrt((xi - xj)**2 + (yi - yj)**2)

ruta, dist_total = resolver_tsp_mtz(nombres, distancias)

if ruta:
    print(f"Distancia total optima: {dist_total:.2f} km")
    print("Ruta optima:")
    for paso, ciudad_idx in enumerate(ruta):
        print(f"  Paso {paso+1}: {nombres[ciudad_idx]}")

Tip: Las restricciones MTZ hacen que el modelo sea polinomial en variables, pero el TSP sigue siendo NP-hard. Para instancias grandes (>100 ciudades) se usan heurísticas como Lin-Kernighan o metaheurísticas (algoritmos genéticos, colonia de hormigas). Para <20 ciudades, el MIP con Gurobi lo resuelve en segundos.

Optimización a Gran Escala con Pandas

En problemas industriales reales, los datos viven en DataFrames. Integrar pandas con PuLP o CVXPY es una habilidad esencial para construir modelos mantenibles.

import pandas as pd
import pulp

# Ejemplo: planificacion de produccion multi-producto multi-periodo
# Los datos vienen de un DataFrame (simulando una hoja de calculo o BD)

# Datos de productos
productos_df = pd.DataFrame({
    "producto":      ["P1", "P2", "P3", "P4", "P5"],
    "ganancia":      [120,  200,   80,  150,  180],
    "mat_prima_kg":  [1.5,  2.0,  0.8,  1.2,  1.8],
    "horas_maquina": [2.0,  3.5,  1.0,  2.5,  3.0],
    "demanda_max":   [200,  150,  300,  180,  120],
})

# Recursos disponibles
recursos = {"mat_prima_kg": 500, "horas_maquina": 800}

# Construir modelo dinamicamente desde el DataFrame
modelo = pulp.LpProblem("Produccion_Multi_Producto", pulp.LpMaximize)

# Variables: una por cada fila del DataFrame
variables = {
    row["producto"]: pulp.LpVariable(
        f"x_{row['producto']}",
        lowBound=0,
        upBound=row["demanda_max"]
    )
    for _, row in productos_df.iterrows()
}

# Funcion objetivo desde el DataFrame
modelo += pulp.lpSum(
    row["ganancia"] * variables[row["producto"]]
    for _, row in productos_df.iterrows()
)

# Restricciones de recursos desde el DataFrame
for recurso, disponible in recursos.items():
    modelo += (
        pulp.lpSum(
            row[recurso] * variables[row["producto"]]
            for _, row in productos_df.iterrows()
        ) <= disponible,
        f"Recurso_{recurso}"
    )

modelo.solve(pulp.PULP_CBC_CMD(msg=0))

# Resultados en un DataFrame
productos_df["produccion_optima"] = productos_df["producto"].map(
    lambda p: variables[p].value()
)
productos_df["ganancia_generada"] = (
    productos_df["produccion_optima"] * productos_df["ganancia"]
)

print(f"Ganancia total: ${productos_df['ganancia_generada'].sum():,.0f}\n")
print(productos_df[["producto", "produccion_optima", "ganancia_generada"]]
      .to_string(index=False))

# Verificar uso de recursos
print("\n--- Uso de Recursos ---")
for recurso, disponible in recursos.items():
    usado = (productos_df["produccion_optima"] * productos_df[recurso]).sum()
    print(f"  {recurso}: {usado:.1f} / {disponible} ({usado/disponible*100:.1f}%)")

Column Generation: Cuando Hay Demasiadas Variables

En algunos problemas, el número de variables es exponencial (problema de corte de barras, generación de horarios, rutas de vehículos). La técnica de Column Generation resuelve estos problemas sin enumerar todas las variables:

1. Problema Maestro Restringido (RMP): Empieza con un subconjunto pequeño de columnas (variables) y resuelve el LP.
2. Subproblema de Pricing: Usa los precios duales del RMP para encontrar una nueva columna con costo reducido negativo (columna que mejora la solución).
3. Iteración: Agrega la nueva columna al RMP y resuelve de nuevo. Repite hasta que no haya columnas con costo reducido negativo.

Si \( \bar{c}_j \geq 0 \) para todas las columnas posibles, la solución actual es óptima.

import pulp
import numpy as np

def column_generation_cutting_stock(longitud_barra, pedidos):
    """
    Problema de corte de barras (Cutting Stock Problem) con Column Generation.
    longitud_barra: longitud de cada barra disponible
    pedidos: dict {longitud_pieza: cantidad_demandada}
    """
    piezas   = list(pedidos.keys())
    demanda  = list(pedidos.values())
    m        = len(piezas)

    # Patrones iniciales: uno por pieza (corte trivial)
    # patron[j] = vector con cuantas piezas de cada tipo corta el patron j
    patrones = []
    for i, l in enumerate(piezas):
        patron = [0] * m
        patron[i] = int(longitud_barra // l)
        patrones.append(patron)

    iteracion = 0
    while True:
        iteracion += 1

        # ---- Problema Maestro Restringido ----
        rmp = pulp.LpProblem("RMP", pulp.LpMinimize)
        x   = [pulp.LpVariable(f"x_{j}", lowBound=0) for j in range(len(patrones))]

        # Minimizar numero de barras usadas
        rmp += pulp.lpSum(x)

        # Satisfacer demanda de cada pieza
        for i in range(m):
            rmp += pulp.lpSum(patrones[j][i] * x[j]
                              for j in range(len(patrones))) >= demanda[i]

        rmp.solve(pulp.PULP_CBC_CMD(msg=0))

        # Precios duales (pi)
        pi = [c.pi for c in rmp.constraints.values()]

        # ---- Subproblema de Pricing (Knapsack) ----
        # Encontrar patron que maximiza valor dual - 1
        # max sum(pi[i] * y[i]) s.t. sum(piezas[i]*y[i]) <= longitud_barra
        pricing = pulp.LpProblem("Pricing", pulp.LpMaximize)
        y = [pulp.LpVariable(f"y_{i}", lowBound=0, cat='Integer') for i in range(m)]

        pricing += pulp.lpSum(pi[i] * y[i] for i in range(m))
        pricing += pulp.lpSum(piezas[i] * y[i] for i in range(m)) <= longitud_barra

        pricing.solve(pulp.PULP_CBC_CMD(msg=0))

        costo_reducido = 1 - pulp.value(pricing.objective)
        print(f"  Iteracion {iteracion}: barras={pulp.value(rmp.objective):.1f}, "
              f"costo_reducido={costo_reducido:.4f}")

        # Si el costo reducido >= 0, no hay columna que mejore: optimo
        if costo_reducido >= -1e-6:
            break

        # Agregar nuevo patron al RMP
        nuevo_patron = [int(round(y[i].value())) for i in range(m)]
        patrones.append(nuevo_patron)

    return pulp.value(rmp.objective), patrones

# Ejemplo: barras de 100 cm, necesitamos piezas de distintas longitudes
longitud = 100
pedidos  = {25: 80, 35: 50, 45: 30, 60: 20}

print("Resolviendo Cutting Stock con Column Generation...")
barras_necesarias, patrones_finales = column_generation_cutting_stock(longitud, pedidos)
print(f"\nBarras totales necesarias: {barras_necesarias:.0f}")
print(f"Patrones de corte generados: {len(patrones_finales)}")

Checklist para Modelos MIP en Producción

Antes de desplegar un modelo de optimización en un sistema productivo, verifica estos puntos:

• Validación: ¿El modelo produce soluciones factibles y sensatas en datos conocidos?
• Tiempo de resolución: ¿Termina dentro del tiempo disponible para el proceso de negocio?
• Robustez: ¿Qué pasa si faltan datos o cambian los parámetros? ¿El modelo sigue siendo factible?
• Interpretabilidad: ¿Los tomadores de decisión entienden la solución y confían en ella?
• Warm start: ¿Puedes proporcionar una solución inicial de buena calidad para acelerar el solver?
• Logging: ¿Guardas el gap MIP, el tiempo de resolución y la solución en cada ejecución?

Conclusión

Con este artículo completas la trilogía de programación lineal con Python a nivel avanzado. Has aprendido:

• La mecánica interna de Branch and Bound y por qué la formulación del modelo importa tanto.
• Las diferencias entre solvers: HiGHS (open-source, producción), Gurobi (industrial, más rápido) y CPLEX.
• Programación Cuadrática con CVXPY: el modelo de Markowitz para optimización de portafolios de inversión y la frontera eficiente.
• El Problema del Viajante (TSP) formulado con restricciones MTZ y resuelto como MIP.
• Cómo integrar pandas con PuLP para construir modelos escalables desde datos reales.
• Los fundamentos de Column Generation para problemas con espacio de variables exponencial.

La programación lineal y sus extensiones son la columna vertebral de la optimización de decisiones en logística, finanzas, manufactura y telecomunicaciones. Dominarlas te posiciona para resolver problemas de alto impacto que ningún modelo de machine learning puede abordar directamente. El siguiente paso es explorar la optimización estocástica (cuando los parámetros son inciertos) y la optimización robusta para decisiones bajo incertidumbre.