Un Algoritmo de Generación de Columnas para el Capacitated Vehicle Routing Problem

El algoritmo de Generación de Columnas sirve para resolver problemas de programación lineal cuyo tamaño es demasiado grande, descomponiendo el problema original en partes tractables, siendo una herramienta eficaz en el arsenal de cualquier optimizador. Aquí enseñaré como implementar dicho esquema para el problema de Ruteo de Vehículos con Capacidad (CVRP en inglés), las ideas claves y su comparación contra una implementación más simple, en Julia.

Dicha implementación está fuertemente inspirada por la implementación original impartida por José Emmanuel Gómez Rocha en la III Escuela de Primavera de Investigación de Operaciones.

Formulación matemática original
Implementación del modelo original
Generación de Columnas
1. Modelo de GC del CVRP
2. Modelo Precios o pricing
Implementación de Generación de Columnas
Referencias

Formulación matemática original

El CVRP es un problema fundamental en la logística y optimización del transporte, representado con un grafo, que consiste en lo siguiente:

Un depósito central donde los vehículos inician y terminan sus rutas
Un conjunto de clientes, cada uno con:

Una ubicación específica (típicamente en coordenadas)
Una demanda conocida (cantidad de productos que necesitan)

Una flota de vehículos idénticos, cada uno con:

Una capacidad máxima fija
Los mismos costos de operación

El objetivo es encontrar el conjunto óptimo de rutas que:

Atienda a todos los clientes exactamente una vez
Asegure que la demanda total en cada ruta no exceda la capacidad del vehículo
Minimice la distancia/costo total de viaje
Cada vehículo debe comenzar y terminar en el depósito

Conjuntos

$k = \{1,...,p\}$ : vehículos.
$V = \{0, 1, \ldots, n\}$ : conjunto de vértices, donde el vértice 0 representa el depósito y los vértices 1 a $n$ representan a los clientes.
$E = \{(i, j) : i, j \in V, i \neq j\}$ : conjunto de aristas que conectan los vértices.
$S \subset V$ : subconjunto del conjunto de vértices.

Parámetros

$c_{ij}$ : costo de la distancia recorrida en una arista que conecta los vértices $i,j$ .
$Q$ : capacidad del vehículo, se considera homogénea.
$d_i$ : demanda cliente $i$

Variables de decisión

$X_{ijk} \in \{0,1\}$ : tiene un valor de 1 si el arco del nodo $i$ al nodo $j$ está en la ruta óptima y es conducido por vehículo $k$ , con $X_{ijk} = 0$ si $i=j$ .

Función objetivo

Busca minimizar los costos de recorrer los clíentes usando los vehículos.

\text{Min} \quad Z = \sum_{k=1}^{p}\sum_{i=0}^{n}\sum_{j=0}^{n}c_{ij}X_{ijk}

Restricciones

De los clientes

De flujo:

\sum_{i=0: i \neq l}^{n}X_{ilk} = \sum_{j=0: j \neq l}^{n}X_{ljk} \forall\ k = \{1,...,p\}, l = \{0, 1, \ldots, n\}

Un cliente es visitado alguna vez:

\sum_{k=1}^{p}\sum_{j=0}^{n}X_{ijk} =1\ \forall\ i\ \{0, 1, \ldots, n\}

De los vehículos

De capacidad:

\sum_{i=0}^{n}\sum_{j=0:j \neq i}^{n}d_iX_{ijk} \leq Q\ \forall\ k = \{1,...,p\}

Que salgan los vehículos del depot:

\sum_{j=1}^{n}X_{0jk} \leq 1\ \forall\ k = \{1,...,p\}

De los subtours

Eliminación de subtours:

\sum_{k=1}^{p}\sum_{i=0}^{S}\sum_{j=0: i \neq j}^{S}X_{ijk} \leq |S|-1 \quad \forall\ S \subset V\setminus\{0\}s

La restricción 6 incrementa de manera exponencial conforme aumenta el número de clientes, por lo cual, se utilizan las conocidas "Lazy constraints" mediante el uso de un callback.

Implementación del modelo original

Empezamos a definir la función principal que resuelve el problema:

function solve_vrp(data)
    # Calcular distancias
    dist = data["distance"]
    # Extraer parametros
    n_clients = data["n_clients"]
    demand = data["demands"]
    cap = data["capacity"]
    # Número de vehículos, se ajusta en base a la instancia
    n_vehicles = data["n_vehicles"]

Posteriormente definimos el modelo de JuMP. Al usar Lazy Constraints y Callbacks, debemos de definir un modelo directo, definimos la variable binaria de decisión $X_{ijk}$

# Inicializar el modelo
    model = direct_model(Gurobi.Optimizer())
    @variable(model, x[1:n_clients+1, 1:n_clients+1, 1:n_vehicles], Bin)

Definimos las restricciones de los clientes

# Restricciones de clientes
    for k in 1:n_clients+1
        for u in 1:n_vehicles
            # Restricción de flujo, cada vehículo que entra a un cliente debe de salir, equivalente a la restricción 2
            @constraint(model, sum(x[i, k, u] for i in 1:n_clients+1 if i != k) ==
                                sum(x[k, j, u] for j in 1:n_clients+1 if j != k))
        end
    end

    for l in 2:n_clients+1
        for u in 1:n_vehicles
            # Un vehículo no puede ir de un mismo cliente a si mismo
            @constraint(model, x[l, l, u] == 0)
        end
    end

    for i in 2:n_clients+1
        # Cada cliente debe ser visitado una vez, equivalente a la restricción 3
        @constraint(model, sum(x[i, j, u] for u in 1:n_vehicles, j in 1:n_clients+1 if j != i) == 1)
    end

Ahora las restricciones de los vehículos

# Restricciones de vehículos
    for u in 1:n_vehicles
        # Restricción de capacidad, equivalente a la Restricción 4
        @constraint(model, sum(demand[i] * x[i, j, u] for i in 1:n_clients+1, j in 1:n_clients+1 if i != j) <= cap)

        # Restricción de salida del depósito, equivalente a la Restricción 5
        @constraint(model, sum(x[1, j, u] for j in 2:n_clients+1) <= 1)
    end

Finalmente, definimos una Lazy Constraint como Callback para las restricciones de subtours. Una Lazy Constraint es una clase de restricciones que permite agregar de manera dinámica al modelo durante la optimización una restricción que está siendo violada, ya que no todas las restricciones se añaden al inicio modelo, solo conforme sea necesario.

Cada vez que llegamos a una nueva solución (definida por el código de callback MIPSOL), obtenemos el valor de nuestra variable de decisión $X$ y le damos forma de un grafo dirigido, y revisamos si existen ciclos que representan la ruta de cada vehículo que no empiecen en el depósito, indicando un subtour infactible. Arreglamos dicha ruta y agregamos la restricción explícita de dicho subtour y actualizamos el modelo.

# Definir el callback de eliminación de subtours
    function subtour_elim(cb_data, cb_where::Cint)
        if cb_where == GRB_CB_MIPSOL
            Gurobi.load_callback_variable_primal(cb_data, cb_where)
            x_sol = callback_value.(cb_data, x)
            x_sol = reshape(x_sol, n_clients + 1, n_clients + 1, n_vehicles)
            x_sol .= x_sol .> 0.5

            G_sol = [SimpleDiGraph(n_clients + 1) for _ in 1:n_vehicles]

            for u in 1:n_vehicles
                for i in 1:n_clients+1
                    for j in 1:n_clients+1
                        if x_sol[i, j, u] == 1
                            add_edge!(G_sol[u], i, j)
                        end
                    end
                end
            end

            for u in 1:n_vehicles
                for cycle in simplecycles(G_sol[u])
                    if cycle[1] != 1  # Si el ciclo no empieza con el depósito algo está mal
                        subtour = cycle
                        append!(subtour, cycle[1])  # Prepara la constraint para eliminar este ciclo que representa un subtour
                        constraint = @build_constraint(sum(x[subtour[i], subtour[i+1], u] for i in 1:(length(subtour)-1)) <= length(subtour) - 2)
                        MOI.submit(model, MOI.LazyConstraint(cb_data), constraint)
                    end
                end
            end
        end
    end

Posteriormente, definimos la función objetivo y ajustamos parámetros del modelo y optimizamos.

# Función objetivo: minimizar la distancia total
    @objective(model, Min, sum(dist[i, j] * x[i, j, u] for i in 1:n_clients+1, j in 1:n_clients+1, u in 1:n_vehicles))
    MOI.set(model, MOI.RawOptimizerAttribute("LazyConstraints"), 1)
    MOI.set(model, Gurobi.CallbackFunction(), subtour_elim)
    set_optimizer_attribute(model, "TimeLimit", 600)
    optimize!(model)

Finalmente, extraemos la solución del modelo y la interpretamos y graficamos

# Extraer valores de x
    x_sol = value.(x)

    routes = []
    for u in 1:n_vehicles
        G_sol = SimpleDiGraph(n_clients + 1)
        for i in 1:n_clients+1
            for j in 1:n_clients+1
                if x_sol[i, j, u] > 0.5
                    add_edge!(G_sol, i, j)
                end
            end
        end

        # Encuentra las rutas
        for cycle in simplecycles(G_sol)
            if 1 in cycle  # Considera las rutas que empiezan en el depósito
                # La ruta debe de empezar y acabar en el depósito
                depot_index = findfirst(x -> x == 1, cycle)
                route = circshift(cycle, 1 - depot_index)

                # Añade el depósito al final para completar la ruta para su impresión
                push!(route, 1)

                # Convertir a indices empezando en 0 para paridad
                zero_indexed_route = route .- 1

                println("Ruta para el vehículo $u: ", zero_indexed_route)
                push!(routes, zero_indexed_route)
            end
        end
    end
end

Todo esto se llama desde la función principal del programa

function main()
    filename = ARGS[1]
    dist = [
        0 548 776 696 582 274 502 194 308 194 536 502 388 354 468 776 662;
        548 0 684 308 194 502 730 354 696 742 1084 594 480 674 1016 868 1210;
        776 684 0 992 878 502 274 810 468 742 400 1278 1164 1130 788 1552 754;
        696 308 992 0 114 650 878 502 844 890 1232 514 628 822 1164 560 1358;
        582 194 878 114 0 536 764 388 730 776 1118 400 514 708 1050 674 1244;
        274 502 502 650 536 0 228 308 194 240 582 776 662 628 514 1050 708;
        502 730 274 878 764 228 0 536 194 468 354 1004 890 856 514 1278 480;
        194 354 810 502 388 308 536 0 342 388 730 468 354 320 662 742 856;
        308 696 468 844 730 194 194 342 0 274 388 810 696 662 320 1084 514;
        194 742 742 890 776 240 468 388 274 0 342 536 422 388 274 810 468;
        536 1084 400 1232 1118 582 354 730 388 342 0 878 764 730 388 1152 354;
        502 594 1278 514 400 776 1004 468 810 536 878 0 114 308 650 274 844;
        388 480 1164 628 514 662 890 354 696 422 764 114 0 194 536 388 730;
        354 674 1130 822 708 628 856 320 662 388 730 308 194 0 342 422 536;
        468 1016 788 1164 1050 514 514 662 320 274 388 650 536 342 0 764 194;
        776 868 1552 560 674 1050 1278 742 1084 810 1152 274 388 422 764 0 798;
        662 1210 754 1358 1244 708 480 856 514 468 354 844 730 536 194 798 0
    ]

    dist = [dist[i+1, j+1] for i in 0:size(dist, 1)-1, j in 0:size(dist, 2)-1] # se arregla para el índices basados en 1

    # Demandas
    demand = [0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]
    demand = [demand[i+1] for i in 0:length(demand)-1] # se arregla para el índices basados en 1

    # Número de vehículos
    n_vehicles = 5

    # Capacidad homogénea de los vehículos
    cap = 15
    data = Dict()
    data["demands"] = demand
    data["distance"] = dist
    data["n_vehicles"] = n_vehicles
    data["capacity"] = cap
    data["n_clients"] = length(demand) - 1
    solve_vrp(data)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end

Resolviendo para la instancia del código obtenemos una solución entera con valor de función objetivo de 6208, con una mejor cota con valor de función objetivo de 5330.49, obteniendo un gap de optimalidad del 14.1 % en 5 minutos, obteniendo las siguientes rutas:

Ruta para vehículo 1: [0, 5, 6, 2, 8, 0]
Ruta para vehículo 2: [0, 13, 15, 11, 12, 0]
Ruta para vehículo 3: [0, 1, 3, 4, 7, 0]
Ruta para vehículo 4: [0, 14, 16, 10, 9, 0]

No sabemos si 6208 representa la solución óptima, pero podemos mejorar esto con Generación de Columnas

Generación de Columnas

En optimización lineal, podemos representar a nuestras variables como las columnas de nuestra matriz de coeficientes y a las restricciones como nuestras filas. A pesar de que la programación lineal se resuelve en tiempo polinomial, podemos tener problemas con un número muy grande de variables (columnas). La intuición de la generación de columnas es resolver un problema restringido que no incluya explícitamente todas las variables, empezando con un subconjunto de dichas variables, agregando de manera iterativa nuevas variables con potencial para mejorar la función objetivo, generando columnas en la matriz. Una vez que ya ninguna otra columna mejore la función objetivo, se detiene el algoritmo.

El algoritmo considera dos problemas: el problema maestro restringido y el subproblema. El problema maestro es el problema original en el que sólo se considera un subconjunto de variables. El subproblema es un nuevo problema creado para identificar una variable mejoradora (es decir, que puede mejorar la función objetivo del problema maestro). La gran dificultad a la hora de implementar el algoritmo está en definir una variable que mejore al problema maestro, esto se consigue al encontrar la variable con el mayor costo reducido negativo, tratando con problemas de minimización. El costo reducido nos dice cuánto cambiaría la función objetivo si forzamos que una unidad de la variable no básica entre a la solución. Si el costo reducido es positivo (en un problema de minimización), la variable debe permanecer en cero pues aumentarla empeoraría la función objetivo, si es negativo, podría mejorar la solución al introducir esa variable y si es cero, hay soluciones alternativas óptimas.

Determinar cual variable tiene el mayor costo reducido negativo es la responsabilidad del subproblema, y formular el subproblema es el arte detrás de la generación de columnas. El subproblema también es conocido como problema de pricing o precios, depende de la estructura del problema original, siendo la función objetivo del problema de precios el costo reducido de las variables duales de las restricciones del problema maestro.

La generación de columnas está íntimamente ligada a la teoría de la dualidad en programación lineal. Cuando resolvemos el problema maestro restringido, obtenemos valores duales óptimos que podemos interpretar como "precios sombra" o el valor marginal de nuestros recursos. Estos valores duales son cruciales porque nos ayudan a determinar si vale la pena añadir nuevas variables (columnas) a nuestro problema maestro, evaluando si el costo reducido (calculado como el coeficiente original menos el producto de las duales por los coeficientes de la columna) es favorable para nuestra función objetivo.

El proceso es iterativo: resolvemos el maestro restringido, usamos sus duales en el subproblema para buscar columnas prometedoras (aquellas con costo reducido negativo en minimización), y si las encontramos, las añadimos al maestro. La teoría de dualidad nos garantiza que cuando no existan más columnas con costo reducido favorable, hemos alcanzado el óptimo global. Esto es particularmente poderoso porque nos permite manejar problemas con un número muy grande de variables potenciales sin tener que considerarlas todas explícitamente, generando solo aquellas que son verdaderamente útiles para mejorar nuestra solución, alcanzando el óptimo global de forma rigurosa.

Modelo de GC del CVRP

Sea $\Omega$ el conjunto de rutas de vehículos factibles. Sea $c_k$ el costo de la ruta $r_k \in \Omega$ . Sea $a_{jk} = 1$ si la ruta $r_k$ visita al cliente $j$ y 0 en caso contrario. Sea $b_{ijk} = 1$ si la ruta $r_k$ viaja directamente del cliente $i$ al cliente $j$ y 0 en caso contrario. Note las siguientes relaciones:

c_k = \sum_{i=0}^n \sum_{j=0}^n c_{ij}b_{ijk} \qquad \qquad a_{jk} = \sum_{i=0}^n b_{ijk}

Sea $z_k = 1$ si se selecciona la ruta $r_k$ y 0 en caso contrario. El modelo de Cobertura de Conjuntos para el Problema de Ruteo de Vehículos tiene el siguiente objetivo:

\min \sum_{k \in \Omega}c_k y_k

y las siguientes restricciones:

\sum_{k=1}^{p} y_{k} \leq p

\sum_{r_k \in \Omega} a_{jk} y_k \geq 1 \qquad \forall j \in \{1,...,n\}

Asegurar que todos los clientes sean visitados.

Para claridad, la relajación lineal de cualquier formulación MIP es donde todas las variables enteras se permiten ser continuas. Dado que el problema maestro requieren que cada $z_k \in \{0, 1\}$ , su relajación lineal se reducen a:

0 \leq y_k \leq 1 \qquad \forall r_k \in \Omega'

Modelo Precios o pricing

El Problema de Precios para el problema maestro (set covering )difiere del CVRP en las siguientes maneras:

El índice $k$ se elimina de las variables de arco, $X$ , ya que ahora se resuelve para una sola ruta por vehículo.
El objetivo es minimizar el costo reducido de la ruta.

Por la definición de $b_{ijk}$ , tenemos que:

b_{ijk} = X_{ijk} \qquad \forall i, j \in \{0, ..., n\}, \; k \in \Omega

Luego, por la definición de $a_{jk}$ :

a_{jk} = \sum_{i=0}^n X_{ijk} \qquad \forall j \in \{0, ..., n\}, \; k \in \Omega

Por lo tanto, nuestra expresión para los costos reducidos se puede reescribir de la siguiente manera:

c_k - A_k^T y = {\sum_{i = 0}^{n}{\sum_{j = 0}^{n}{c_{ij}x_{ijk}}}} - {\sum_{i = 1}^{n} \Bigg( W_i \sum_{j = 0}^{n} x_{ijk} \Bigg)}

Por lo tanto, la formulación es la siguiente:

\min {\sum_{i = 0}^{n}{\sum_{j = 0}^{n}\Bigg(c_{ij}-W_i}}\Bigg)X_{ij}

Sujeto a:

Retricción de flujo:

\sum_{i=0: i \neq l}^{n}X_{il} = \sum_{j=0: j \neq l}^{n}X_{lj} \forall\ l = \{0, 1, \ldots, n\} \tag{15}

La ruta sale del depósito como máximo una vez:

\sum_{j = 1}^{n}{X_{0j}}\ \leq 1

Restricción de capacidad:

\sum_{i=0}^{n}\sum_{j=0:j \neq i}^{n}d_iX_{ij} \leq Q\

Implementación de Generación de Columnas

Primero, necesitamos un subconjunto de columnas factibles para empezar el problema restringido maestro, podría verse como un warm start del modelo de optimización, pero no es del todo correcto. Para hacerlo, resolvemos el problema con un heurístico, el algoritmo de savings o de Clarke & Wright.

El algoritmo de Clarke & Wright (también conocido como algoritmo de ahorros) funciona de la siguiente manera:

Inicialización:

Se comienza con una solución donde cada cliente tiene su propia ruta (depósito → cliente → depósito)
Esto significa que inicialmente tenemos n rutas si hay n clientes

Cálculo de Ahorros:

Para cada par de clientes i,j se calcula el ahorro:
savings(i,j) = c(0,i) + c(0,j) - c(i,j)
Donde c(i,j) es el costo/distancia entre i y j, y 0 representa el depósito

La intuición es: ¿cuánto nos ahorramos si en lugar de hacer dos viajes separados al depósito, conectamos estos clientes?

Proceso de Unión:

Se ordenan los ahorros de mayor a menor
Para cada par de clientes (i,j) en orden de ahorro descendente:

Se verifica si se pueden unir sus rutas sin violar restricciones:

Los nodos no deben estar ya en la misma ruta
Los nodos deben estar conectados directamente al depósito
La unión debe respetar la capacidad del vehículo

Si todas las condiciones se cumplen, se unen las rutas

Criterio de Parada:

El algoritmo termina cuando ya no hay más uniones factibles, es decir, cuando no quedan pares de rutas que se puedan unir sin violar restricciones

function savings_algorithm(data)
    # Crear grafo dirigido con metadatos
    dist = data["distance"]
    n_clients = data["n_clients"]
    demand = data["demands"]
    cap = data["capacity"]
    n_vehicles = data["n_vehicles"]
    G = MetaDiGraph(n_clients + 1)  # +1 para el depósito (Julia usa índices base 1)
    n_routes = n_clients
    

    # Añadir arcos desde el depósito a todos los clientes y viceversa
    for i in 2:(n_clients+1)  # el depósito es 1 en Julia
        add_edge!(G, 1, i)
        set_prop!(G, 1, i, :weight, dist[1, i-1])
        add_edge!(G, i, 1)
        set_prop!(G, i, 1, :weight, dist[i-1, 1])
    end

    # Calcular ahorros
    savings = Tuple{Tuple{Int,Int},Float64}[]
    for i in 2:(n_clients+1)
        for j in (i+1):(n_clients+1)
            if i != j
                # Ajustar índices para la matriz de distancias que está en base 0
                temp = dist[i, 1] + dist[1, j] - dist[i, j]
                if temp > 0
                    push!(savings, ((i, j), temp))
                end
            end
        end
    end

    # Ordenar ahorros en orden descendente
    sort!(savings, by=x -> x[2], rev=true)

    # Fusionar rutas
    for (sav, _) in savings
        # Condición 1: Verificar si los nodos no están en la misma ruta
        condition_1 = true
        for cycle in simplecycles(G)
            push!(cycle, 1)  # Añadir depósito al final del ciclo
            if (sav[1] in cycle) && (sav[2] in cycle)
                condition_1 = false
                break
            end
        end

        # Condición 2: Verificar si los nodos están directamente conectados al depósito
        condition_2a = false
        condition_2b = false
        for cycle in simplecycles(G)
            push!(cycle, 1)
            if cycle[end-1] != sav[2]  # Cambiado de == a !=
                condition_2b = true
                break
            end
        end
        condition_2a = false
        for cycle in simplecycles(G)
            push!(cycle, 1)
            if cycle[2] == sav[1]  # Verificar que esta lógica coincida con Python
                condition_2a = true
                break
            end
        end

        # Condición 3: Verificar restricción de capacidad
        condition_3 = true
        sav0_ruta_dem = 0.0
        sav1_ruta_dem = 0.0
        condition_3 = true
        sav0_ruta_dem = 0
        sav1_ruta_dem = 0

        for cycle in simplecycles(G)
            if sav[1] in cycle
                # Cambiado para sumar correctamente las demandas de todos los nodos no depósito
                sav0_ruta_dem = sum(demand[i] for i in cycle if i != 1)
            end
            if sav[2] in cycle
                # Cambiado para sumar correctamente las demandas de todos los nodos no depósito
                sav1_ruta_dem = sum(demand[i] for i in cycle if i != 1)
            end
        end

        if sav0_ruta_dem + sav1_ruta_dem > cap
            condition_3 = false
        end

        # Fusionar rutas si se cumplen todas las condiciones
        if condition_1 && condition_2a && condition_2b && condition_3
            if has_edge(G, sav[1], 1) && has_edge(G, 1, sav[2])
                rem_edge!(G, sav[1], 1)
                rem_edge!(G, 1, sav[2])
                add_edge!(G, sav[1], sav[2])
                set_prop!(G, sav[1], sav[2], :weight, dist[sav[1], sav[2]])
                n_routes -= 1
            elseif has_edge(G, 1, sav[1]) && has_edge(G, sav[2], 1)
                rem_edge!(G, 1, sav[1])
                rem_edge!(G, sav[2], 1)
                add_edge!(G, sav[2], sav[1])
                set_prop!(G, sav[2], sav[1], :weight, dist[sav[2], sav[1]])  # Índices corregidos
                n_routes -= 1
            end
        end
    end

    # Imprimir rutas y sus demandas
    for cycle in simplecycles(G)
        push!(cycle, 1)  # Añadir depósito al final
        # Calcular demanda total para la ruta excluyendo el depósito
        total_demand = sum(demand[i] for i in cycle[2:end-1])
        println("Demanda: $total_demand, Ruta: $cycle")
    end

    return G, n_routes
end

Ya que obtuvimos rutas iniciales, debemos de convertirlas a la forma de las matrices $a, b$ y el vector $c$ que acepta el MIP, de las ecuaciones 12, 13 y 14.

function create_solution_matrices(G::MetaDiGraph, n_clients::Int, n_routes::Int, dist)
    # Inicializar matrices
    a = zeros(Int, n_clients, n_routes)
    rutas = Vector{Vector{Int}}()

    # Crear matriz de asignación cliente-ruta
    for (k, ruta) in enumerate(simplecycles(G))
        push!(rutas, ruta)
        for i in 2:n_clients+1
            if i in ruta
                a[i-1, k] = 1
            end
        end
    end

    # Crear matriz de incidencia arco-ruta
    b = zeros(Int, n_clients + 1, n_clients + 1, n_routes)
    for (k, ruta) in enumerate(rutas)
        # Crear una copia para evitar modificar la original
        ruta_temp = copy(ruta)
        # Añadir depósito al inicio y final
        pushfirst!(ruta_temp, 1)  # Añadir depósito al inicio
        push!(ruta_temp, 1)      # Añadir depósito al final
        # Establecer los arcos en la matriz b
        for r in 1:length(ruta_temp)-1
            from_node = ruta_temp[r]
            to_node = ruta_temp[r+1]
            b[from_node, to_node, k] = 1
        end
    end

    # Calcular costos de las rutas
    c = zeros(Float64, n_routes)
    for (k, ruta) in enumerate(rutas)
        # Obtener ruta completa con depósito
        ruta_temp = copy(ruta)
        pushfirst!(ruta_temp, 1)  # Añadir depósito al inicio
        push!(ruta_temp, 1)      # Añadir depósito al final
        # Calcular costo
        for r in 1:length(ruta_temp)-1
            from_node = ruta_temp[r]
            to_node = ruta_temp[r+1]
            c[k] += dist[from_node, to_node]
        end
    end
    return a, b, c, rutas
end

Ya que definimos esas funciones auxiliares y podemos generar un subconjunto de columnas factibles, empezamos a construir nuestra función principal de Generación de Columnas

function solve_cvrp_cg(data)
    # Obtener datos del problema
    demand = data["demands"]
    cap = data["capacity"]
    n_clients = data["n_clients"]
    dist = data["distance"]
    n_vehicles = data["n_vehicles"]

    # Registrar tiempo de inicio
    time_start = Dates.now()

    # Generar solución inicial usando algoritmo de ahorros
    G, n_routes = savings_algorithm(data)

    # Crear matrices necesarias para el modelo
    a, b, c, rutas = create_solution_matrices(G, n_clients, n_routes, dist)

    # Inicializar problema maestro
    master = Model(Gurobi.Optimizer)
    set_silent(master)

    # Añadir variables de decisión
    @variable(master, y[k=1:n_routes] >= 0)

    # Establecer función objetivo
    @objective(master, Min, sum(c[k] * y[k] for k = 1:n_routes))

    # Añadir restricciones
    # Restricción de número máximo de vehículos
    @constraint(master, fleet_size, sum(y[k] for k = 1:n_routes) <= n_vehicles)
    
    # Restricción de cobertura: cada cliente debe ser visitado al menos una vez
    @constraint(master, covering[i=1:n_clients],
        sum(a[i, k] * y[k] for k = 1:n_routes) >= 1)

    # Optimización inicial
    optimize!(master)
    println("Iteracion 0: \nValor Funcion Objetivo: ", round(objective_value(master), digits=2))

Definimos el modelo del problema maestro con variables de decisión $y_k$ de la ecuación 8 y lo optimizamos para obtener la primera iteración. Ahora tenemos que definir el subproblema. Primero, preparamos datos para el problema de pricing

# Crear matriz de distancias expandida para el pricing
    dist_2 = fill(999999.0, n_clients + 2, n_clients + 2)
    dist_2[1:n_clients+1, 1:n_clients+1] = dist
    for i in 2:n_clients+1
        dist_2[i, n_clients+2] = dist[i, 1]
        dist_2[i, 1] = 999999.0
    end
    dist_2[n_clients+2, n_clients+2] = 0

    # Vector de demanda expandido
    dem_2 = vcat(demand, demand[1])

Definimos el bucle principal de iteraciones de la generación de columnas, obtenemos los valores duales del problema maestro, configuramos y resolvemos el problema de pricing con dichos valores duales y buscamos rutas nuevas que puedan mejorar la función objetivo, reutilizando el mismo callback de eliminación de subtours.

# Column generation iterations
    for iter in 1:150
        # Obtener valores duales del problema maestro de la restriccion de cobertura y de los vehiculos usados
        w = [dual.(covering)..., dual(fleet_size)]

        # Inicializamos problema de pricing
        pricing = direct_model(Gurobi.Optimizer())
        set_silent(pricing)

        # Variable de decision Xij
        @variable(pricing, x[i=1:n_clients+2, j=1:n_clients+2], Bin)

        # Restriccion de capacidad
        @constraint(pricing, capacity,
            sum(dem_2[i] * sum(x[i, j] for j in 1:n_clients+2)
                for i in 1:n_clients+2) <= cap)

        # Restriccion de entrada
        @constraint(pricing, entry_flow,
            sum(x[i, n_clients+2] for i in 1:n_clients+1) == 1)

        # Restriccion de salida
        @constraint(pricing, exit_flow,
            sum(x[1, j] for j in 2:n_clients+2) == 1)

        # Restriccion de conservacion de flujo
        @constraint(pricing, flow_conservation[k in 2:n_clients+1],
            sum(x[i, k] for i in 1:n_clients+2) ==
            sum(x[k, j] for j in 1:n_clients+2))

        # Restriccion de no ciclos en un mismo cliente
        @constraint(pricing, no_loops[i in 1:n_clients+2], x[i, i] == 0)

        # Objetivo de costo reducido
        @objective(pricing, Min,
            sum(dist_2[i, j] * x[i, j] for i in 1:n_clients+2, j in 1:n_clients+2) - # minimizamos la matriz de distancias ajustada
            sum(w[i-1] * sum(x[i, j] for j in 1:n_clients+2) for i in 2:n_clients+1) -  # Ajustando indices i-1, minizamos los valores duales
            w[end] * sum(x[1, j] for j in 2:n_clients+2)
        )

        # Callback de subtours
        function subtour_elim(cb_data, cb_where::Cint)
            if cb_where == GRB_CB_MIPSOL
                # Cargamos los valores actuales
                Gurobi.load_callback_variable_primal(cb_data, cb_where)
                x_sol = callback_value.(cb_data, x)
                x_sol = reshape(x_sol, n_clients + 2, n_clients + 2)
                x_sol .= x_sol .> 0.5

                # Creamos grafo dirigido de la solucion actual
                G_sol = SimpleDiGraph(n_clients + 2)
                for i in 1:n_clients+2
                    for j in 1:n_clients+2
                        if x_sol[i, j] == 1
                            add_edge!(G_sol, i, j)
                        end
                    end
                end

                # Revisamos los ciclos del grafo de solucion
                for cycle in simplecycles(G_sol)
                    if cycle[1] != 1 
                        subtour = copy(cycle)
                        push!(subtour, cycle[1])  # Cierra el ciclo si no empieza en el depot

                        # Añade la restriccion
                        constraint = @build_constraint(
                            sum(x[subtour[i], subtour[i+1]]
                                for i in 1:(length(subtour)-1)) <= length(subtour) - 2
                        )
                        MOI.submit(pricing, MOI.LazyConstraint(cb_data), constraint)
                    end
                end
            end
        end

        # Ajusta parametros y añade el callback al modelo
        MOI.set(pricing, MOI.RawOptimizerAttribute("LazyConstraints"), 1)
        MOI.set(pricing, Gurobi.CallbackFunction(), subtour_elim)

        # Resuelve el problema de pricing
        optimize!(pricing)

        # Revisa terminacion
        if objective_value(pricing) >= -0.0001
            println("\nGeneracion de columnas completada!")
            break
        end

Al resolver el subproblema, debemos de extraer la ruta obtenida y verificar si es una ruta nueva, para generar y agregar la columna correspondiente al problema maestro.

# Extraer la ruta nueva obtenida del grafo
        sol = value.(x)
        sol = reshape(sol, n_clients + 2, n_clients + 2)
        sol = round.(Int, sol)
        # Encontrar la secuencia de nodos siguiendo los arcos activos
        function extract_path(sol, n_clients)
            path = [1]  # Empezar en el depósito
            current = 1
            while current != n_clients + 2  # Hasta llegar al depósito artificial
                for j in 1:n_clients+2
                    if sol[current, j] > 0.5
                        push!(path, j)
                        current = j
                        break
                    end
                end
            end
            pop!(path)  # Remover el depósito artificial
            return path
        end

        new_route = extract_path(sol, n_clients)

        # Remueve el depot artificial si esta en la ruta
        if !isempty(new_route) && new_route[end] == n_clients + 2
            pop!(new_route)
        end

        # Revisa si la ruta es nueva
        if !any(route == new_route for route in rutas)
            # Actualiza el problema maestro con dicha ruta nueva
            master, rutas, n_routes, a, b, c = update_master_problem(
                master, new_route, dist, n_clients, n_vehicles,
                rutas, n_routes
            )

            optimize!(master)
            println("Iteracion: ", iter)
            println("Funcion objetivo del maestro: ", objective_value(master))
            println("Funcion objetivo del pricing: ", objective_value(pricing))
        end
    end

Defino la función de actualziación del problema maestro y $a, b, c$

# Actualizar master
    function update_master_problem(master::Model, new_route::Vector{Int},
        dist,
        n_clientes::Int, n_vehiculos::Int,
        rutas::Vector{Vector{Int}}, n_rutas::Int)
        # Añadir la nueva ruta al conjunto de rutas
        push!(rutas, new_route)
        n_rutas += 1

        # Reconstruir matrices a, b, y c con la nueva ruta
        # Matriz de asignación cliente-ruta (a[i,k] = 1 si el cliente i está en la ruta k)
        a = zeros(Int, n_clientes, n_rutas)
        for (k, ruta) in enumerate(rutas)
            for i in 2:n_clientes+1  # Comenzar desde 2 ya que 1 es el depósito
                if i in ruta
                    a[i-1, k] = 1
                end
            end
        end

        # Matriz de incidencia arco-ruta (b[i,j,k] = 1 si la ruta k usa el arco (i,j))
        b = zeros(Int, n_clientes + 1, n_clientes + 1, n_rutas)
        for (k, ruta) in enumerate(rutas)
            ruta_with_depot = [ruta..., 1]  # Añadir depósito al final
            for r in 1:length(ruta_with_depot)-1
                b[ruta_with_depot[r], ruta_with_depot[r+1], k] = 1
            end
        end

        # Calcular costos para cada ruta sumando los costos de los arcos utilizados
        c = zeros(Float64, n_rutas)
        for (k, ruta) in enumerate(rutas)
            ruta_with_depot = [ruta..., 1]
            for r in 1:length(ruta_with_depot)-1
                c[k] += b[ruta_with_depot[r], ruta_with_depot[r+1], k] *
                        dist[ruta_with_depot[r], ruta_with_depot[r+1]]
            end
        end

        # Resetear y recrear el problema maestro
        empty!(master)

        # Añadir variables de decisión y[k] que indican si se usa la ruta k
        @variable(master, y[k=1:n_rutas] >= 0)

        # Establecer función objetivo: minimizar el costo total de las rutas seleccionadas
        @objective(master, Min, sum(c[k] * y[k] for k = 1:n_rutas))

        # Añadir restricciones
        # Límite en el número de vehículos disponibles
        @constraint(master, fleet_size,
            sum(y[k] for k = 1:n_rutas) <= n_vehiculos)
        
        # Cada cliente debe ser visitado al menos una vez
        @constraint(master, covering[i=1:n_clientes],
            sum(a[i, k] * y[k] for k = 1:n_rutas) >= 1)

        # Configurar parámetros del solver para modo silencioso
        set_silent(master)

        return master, rutas, n_rutas, a, b, c
    end

Finalmente, defino una función que nos ayuda a extraer las rutas finales del problema maestro y las rutas generadas:

function extract_final_routes(y_val::Vector{Float64}, rutas::Vector{Vector{Int}},
        n_clientes::Int, dist)
        # Crear un grafo dirigido para almacenar todas las rutas
        G_final = MetaDiGraph(n_clientes + 1)  # +1 para el depósito
        routes = Vector{Vector{Int}}()

        # Procesar cada ruta seleccionada
        for (k, route) in enumerate(rutas)
            if y_val[k] > 0.5  # Si la ruta está seleccionada en la solución
                # Añadir ruta al grafo
                route_with_depot = [route..., 1]  # Añadir depósito al final
                for i in 1:(length(route_with_depot)-1)
                    from_node = route_with_depot[i]
                    to_node = route_with_depot[i+1]
                    # Añadir arco si no existe
                    if !has_edge(G_final, from_node, to_node)
                        add_edge!(G_final, from_node, to_node)
                        # Añadir distancia como peso del arco
                        set_prop!(G_final, from_node, to_node, :weight, dist[from_node, to_node])
                    end
                end
                
                # Almacenar la ruta (convertir a indexación base-0)
                zero_indexed_route = route_with_depot .- 1  # Convertir a indexación base 0
                push!(routes, zero_indexed_route)
                println("Ruta $k: ", zero_indexed_route)
            end
        end

        # Calcular costo total sumando los pesos de todos los arcos
        total_cost = sum(get_prop(G_final, e.src, e.dst, :weight) for e in edges(G_final))
        println("Costo total: ", round(total_cost, digits=2))
        
        return G_final, routes
    end

Integramos dicha función al finalizar nuestro ciclo de iteraciones de pricing

# Obtener solucion final
    y_val = value.(master[:y])
    selected_routes = [rutas[k] for k in 1:n_routes if y_val[k] > 0.5]

    y_val = value.(master[:y])
    G_final, routes = extract_final_routes(y_val, rutas, n_clients, dist)
    time_end = Dates.now()

    delta = time_end - time_start

    println(delta)

    return routes, objective_value(master)
end

function main()
    dist = [
        0 548 776 696 582 274 502 194 308 194 536 502 388 354 468 776 662;
        548 0 684 308 194 502 730 354 696 742 1084 594 480 674 1016 868 1210;
        776 684 0 992 878 502 274 810 468 742 400 1278 1164 1130 788 1552 754;
        696 308 992 0 114 650 878 502 844 890 1232 514 628 822 1164 560 1358;
        582 194 878 114 0 536 764 388 730 776 1118 400 514 708 1050 674 1244;
        274 502 502 650 536 0 228 308 194 240 582 776 662 628 514 1050 708;
        502 730 274 878 764 228 0 536 194 468 354 1004 890 856 514 1278 480;
        194 354 810 502 388 308 536 0 342 388 730 468 354 320 662 742 856;
        308 696 468 844 730 194 194 342 0 274 388 810 696 662 320 1084 514;
        194 742 742 890 776 240 468 388 274 0 342 536 422 388 274 810 468;
        536 1084 400 1232 1118 582 354 730 388 342 0 878 764 730 388 1152 354;
        502 594 1278 514 400 776 1004 468 810 536 878 0 114 308 650 274 844;
        388 480 1164 628 514 662 890 354 696 422 764 114 0 194 536 388 730;
        354 674 1130 822 708 628 856 320 662 388 730 308 194 0 342 422 536;
        468 1016 788 1164 1050 514 514 662 320 274 388 650 536 342 0 764 194;
        776 868 1552 560 674 1050 1278 742 1084 810 1152 274 388 422 764 0 798;
        662 1210 754 1358 1244 708 480 856 514 468 354 844 730 536 194 798 0
    ]

    dist = [dist[i+1, j+1] for i in 0:size(dist, 1)-1, j in 0:size(dist, 2)-1]
    demand = [0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]
    demand = [demand[i+1] for i in 0:length(demand)-1]
    n_vehicles = 5
    cap = 15
    data = Dict()
    data["demands"] = demand
    data["distance"] = dist
    data["n_vehicles"] = n_vehicles
    data["capacity"] = cap
    data["n_clients"] = length(demand) - 1
    solve_cvrp_cg(data)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end

Ejecutamos todo el código y medimos el tiempo, incluyendo el heurístico de Clarke & Wright y obtenemos en 16 segundos los siguientes resultados con 48 iteraciones del subproblema, con un valor de función objetivo de 6208, el mismo que el valor obtenido con el modelo original, y ya que la generación de columnas paró al no encontrar un costo reducido negativo, significa que es el óptimo global del problema:

Ruta para vehículo 1: [0, 7, 1, 3, 4, 0]
Ruta para vehículo 2: [0, 12, 11, 15, 13, 0]
Ruta para vehículo 3: [0, 9, 10, 16, 14, 0]
Ruta para vehículo 4: [0, 8, 6, 2, 5, 0]

Podemos ver entonces como la generación de columnas supera ampliamente a un modelo ingenuo, garantizando optimalidad en una fracción del tiempo asignado al modelo original. Si bien las rutas están al reves, esto es un artefacto de la interpretación de los ciclos del grafo de solución, ya que la matriz de distancias es simétrica.

Trabajo a futuro incluye una formulación mejor del problema de pricing, ya que ese es el principal problema que se resuelve y que la mayor parte del tiempo del algoritmo se emplea en resolver iterativamente dichos problemas, agregar cortes o esquemas más complejos de Branch-and.

Referencias

Uchoa, Eduardo and Pessoa, Artur and Moreno, Lorenza. "Optimizing with Column Generation: Advanced Branch-Cut-and-Price Algorithms (Part I)" (2024)

A tutorial on using Column Generation to solve MIP problems

"Técnicas de descomposición para problemas de optimización de gran escala en la administración de la cadena de suministro", José Emmanuel Gómez Rocha, III Escuela de Primavera de Investigación de Operaciones