fix bugs and add column generation

Author: BatyLeo

Project.toml

@@ -4,6 +4,7 @@ authors = ["Members of JuliaDecisionFocusedLearning"]
 version = "0.1.0"
 
 [deps]
+ConstrainedShortestPaths = "b3798467-87dc-4d99-943d-35a1bd39e395"
 DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
@@ -27,6 +28,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
+ConstrainedShortestPaths = "0.5.1"
 DataDeps = "0.7"
 Distributions = "0.25"
 DocStringExtensions = "0.9"

src/StochasticVehicleScheduling/StochasticVehicleScheduling.jl

@@ -2,6 +2,8 @@ module StochasticVehicleScheduling
 
 using ..Utils
 using DocStringExtensions: TYPEDEF, TYPEDFIELDS, TYPEDSIGNATURES
+using ConstrainedShortestPaths:
+    stochastic_routing_shortest_path, stochastic_routing_shortest_path_with_threshold
 using Distributions: Distribution, LogNormal, Uniform
 using Graphs:
     AbstractGraph,
@@ -16,7 +18,15 @@ using Graphs:
     inneighbors,
     outneighbors
 using JuMP:
-    Model, @variable, @objective, @constraint, optimize!, value, objective_value, set_silent
+    Model,
+    @variable,
+    @objective,
+    @constraint,
+    optimize!,
+    value,
+    objective_value,
+    set_silent,
+    dual
 using Printf: @printf
 using Random: Random, AbstractRNG, MersenneTwister
 using SparseArrays: sparse
@@ -32,6 +42,7 @@ include("instance/instance.jl")
 
 include("solution/solution.jl")
 include("solution/exact_algorithms/mip.jl")
+include("solution/exact_algorithms/column_generation.jl")
 
 """
 $TYPEDFIELDS
@@ -74,6 +85,7 @@ function Utils.generate_statistical_model(bench::StochasticVehicleSchedulingBenc
 
 export StochasticVehicleSchedulingBenchmark
 
-export compact_linearized_mip, compact_mip
+export compact_linearized_mip,
+    compact_mip, column_generation_algorithm, evaluate_solution, is_feasible
 
 end

src/StochasticVehicleScheduling/instance/instance.jl

@@ -87,7 +87,7 @@ $TYPEDSIGNATURES
 Returns the number of scenarios in instance.
 """
 function get_nb_scenarios(instance::Instance)
-    return size(instance.delays, 2)
+    return size(instance.intrinsic_delays, 2)
 end
 
 """

src/StochasticVehicleScheduling/solution/exact_algorithms/column_generation.jl

@@ -0,0 +1,141 @@
+"""
+    delay_sum(path, slacks, delays)
+
+Evaluate the total delay along path.
+"""
+function delay_sum(path, slacks, delays)
+    nb_scenarios = size(delays, 2)
+    old_v = path[1]
+    R = delays[old_v, :]
+    C = 0.0
+    for v in path[2:(end - 1)]
+        @. R = max(R - slacks[old_v, v], 0) + delays[v, :]
+        C += sum(R) / nb_scenarios
+        old_v = v
+    end
+    return C
+end
+
+"""
+    column_generation(instance::Instance)
+
+Note: If you have Gurobi, use `grb_model` as `model_builder` instead of `glpk_model`.
+"""
+function column_generation(
+    instance::Instance; only_relaxation=false, model_builder=highs_model
+)
+    (; graph, slacks, intrinsic_delays, vehicle_cost, delay_cost) = instance
+
+    nb_nodes = nv(graph)
+    job_indices = 2:(nb_nodes - 1)
+
+    model = model_builder()
+    set_silent(model)
+
+    @variable(model, λ[v in 1:nb_nodes])
+
+    @objective(model, Max, sum(λ[v] for v in job_indices))
+
+    initial_paths = [[1, v, nb_nodes] for v in job_indices]
+    @constraint(
+        model,
+        con[p in initial_paths],
+        delay_cost * delay_sum(p, slacks, intrinsic_delays) + vehicle_cost -
+        sum(λ[v] for v in job_indices if v in p) >= 0
+    )
+    @constraint(model, λ[1] == 0)
+    @constraint(model, λ[nb_nodes] == 0)
+
+    new_paths = Vector{Int}[]
+    cons = []
+
+    while true
+        optimize!(model)
+        λ_val = value.(λ)
+        (; c_star, p_star) = stochastic_routing_shortest_path(
+            graph, slacks, intrinsic_delays, λ_val ./ delay_cost
+        )
+        λ_sum = sum(λ_val[v] for v in job_indices if v in p_star)
+        path_cost = delay_cost * c_star + λ_sum + vehicle_cost
+        if path_cost - λ_sum > -1e-10
+            break
+        end
+        push!(new_paths, p_star)
+        push!(
+            cons,
+            @constraint(
+                model, path_cost - sum(λ[v] for v in job_indices if v in p_star) >= 0
+            )
+        )
+    end
+
+    c_low = objective_value(model)
+    paths = cat(initial_paths, new_paths; dims=1)
+    c_upp, y, _ = compute_solution_from_selected_columns(instance, paths)
+
+    # If relaxation requested or solution is optimal, return
+    if c_upp ≈ c_low || only_relaxation
+        return value.(λ),
+        objective_value(model), cat(initial_paths, new_paths; dims=1), dual.(con),
+        dual.(cons)
+    end
+
+    # else, try to close the gap
+    threshold = (c_upp - c_low - vehicle_cost) / delay_cost
+    λ_val = value.(λ)
+    additional_paths, costs = stochastic_routing_shortest_path_with_threshold(
+        graph, slacks, intrinsic_delays, λ_val ./ delay_cost; threshold
+    )
+
+    return value.(λ),
+    objective_value(model),
+    unique(cat(initial_paths, new_paths, additional_paths; dims=1)),
+    dual.(con),
+    dual.(cons)
+end
+
+"""
+    compute_solution_from_selected_columns(instance::Instance, paths[; bin=true])
+
+Note: If you have Gurobi, use `grb_model` as `model_builder` instead od `glpk_model`.
+"""
+function compute_solution_from_selected_columns(
+    instance::Instance, paths; bin=true, model_builder=highs_model
+)
+    (; graph, slacks, intrinsic_delays, vehicle_cost, delay_cost) = instance
+
+    nb_nodes = nv(graph)
+    job_indices = 2:(nb_nodes - 1)
+
+    model = model_builder()
+    set_silent(model)
+
+    if bin
+        @variable(model, y[p in paths], Bin)
+    else
+        @variable(model, y[p in paths] >= 0)
+    end
+
+    @objective(
+        model,
+        Min,
+        sum(
+            (delay_cost * delay_sum(p, slacks, intrinsic_delays) + vehicle_cost) * y[p] for
+            p in paths
+        )
+    )
+
+    @constraint(model, con[v in job_indices], sum(y[p] for p in paths if v in p) == 1)
+
+    optimize!(model)
+
+    sol = value.(y)
+    return objective_value(model), sol, paths[[sol[p] for p in paths] .== 1.0]
+end
+
+function column_generation_algorithm(instance::Instance)
+    _, _, columns, _, _ = column_generation(instance)
+    _, _, sol = compute_solution_from_selected_columns(instance, columns)
+    col_solution = solution_from_paths(sol, instance)
+    return col_solution
+end

src/StochasticVehicleScheduling/solution/solution.jl

@@ -1,16 +1,17 @@
 # TODO: only keep value field ?
 """
-    Solution
+$TYPEDEF
 
 Should always be associated with an `Instance`.
 
 # Fields
-- `value::BitVector`: for each graph edge of instance, 1 if selected, else 0
-- `path_value::BitMatrix`: each row represents a vehicle, each column a task.
-    1 if task is done by the vehicle, else 0
+$TYPEDFIELDS
 """
 struct Solution
+    "for each graph edge of instance, 1 if selected, else 0"
     value::BitVector       # for every edge of graph 1 if selected, else 0
+    "each row represents a vehicle, each column a task.
+    1 if task is done by the vehicle, else 0"
     path_value::BitMatrix  # list of vehicles paths
 end
 
@@ -23,7 +24,7 @@ function get_nb_vehicles(solution::Solution)
 end
 
 """
-    get_routes(solution::Solution)
+$TYPEDSIGNATURES
 
 Compute routes of solution.
 """
@@ -42,7 +43,7 @@ function get_routes(solution::Solution)
 end
 
 """
-    Solution(value::BitVector, instance::Instance)
+$TYPEDSIGNATURES
 
 Create a Solution from a BitVector value.
 """
@@ -61,7 +62,7 @@ function Solution(value::BitVector, instance::Instance)
 end
 
 """
-    Solution(value::BitVector, instance::Instance)
+$TYPEDSIGNATURES
 
 Create a Solution from a BitMatrix path value.
 """
@@ -81,7 +82,7 @@ function Solution(path_value::BitMatrix, instance::Instance)
 end
 
 """
-    solution_from_paths(paths, instance::Instance)
+$TYPEDSIGNATURES
 
 Create a Solution from routes.
 """
@@ -103,7 +104,7 @@ function solution_from_paths(paths, instance::Instance)
 end
 
 """
-    solution_from_paths(paths, instance::Instance)
+$TYPEDSIGNATURES
 
 Create a Solution from a JuMP array.
 """
@@ -161,7 +162,7 @@ function basic_path_solution(graph::AbstractGraph)
 end
 
 """
-    basic_solution(graph::AbstractGraph)
+$TYPEDSIGNATURES
 
 Create a solution with one vehicle per task.
 """
@@ -179,13 +180,7 @@ function basic_solution(instance::Instance)
 end
 
 """
-    evaluate_task(
-        i_task::Integer,
-        instance::Instance,
-        old_task_index::Integer,
-        old_delay::Real,
-        scenario::Int,
-    )
+$TYPEDSIGNATURES
 
 Evaluate the total delay of task `i_task` in `scenario`, knowing that current delay from task
 `old_task_index` is `old_delay`.
@@ -197,9 +192,9 @@ function evaluate_task(
     old_delay::Real,
     scenario::Int,
 )
-    (; slacks, delays) = instance
+    (; slacks, intrinsic_delays) = instance
 
-    delay = delays[i_task, scenario]
+    delay = intrinsic_delays[i_task, scenario]
     slack = slacks[old_task_index, i_task][scenario]
 
     return delay + max(old_delay - slack, 0)
@@ -240,7 +235,7 @@ function evaluate_scenario(path_value::BitMatrix, instance::Instance, scenario_i
 end
 
 """
-    evaluate_scenario(solution::Solution, instance::Instance, scenario_index::Int)
+$TYPEDSIGNATURES
 
 Compute total delay of scenario.
 """
@@ -249,7 +244,7 @@ function evaluate_scenario(solution::Solution, instance::Instance, scenario_inde
 end
 
 """
-    evaluate_scenario(path_value::BitMatrix, instance::Instance, scenario_index::Int)
+$TYPEDSIGNATURES
 
 Compute total weighted objective of solution.
 """
@@ -267,7 +262,7 @@ function evaluate_solution(path_value::BitMatrix, instance::Instance)
 end
 
 """
-    evaluate_scenario(path_value::BitMatrix, instance::Instance, scenario_index::Int)
+$TYPEDSIGNATURES
 
 Compute total weighted objective of solution.
 """
@@ -314,7 +309,7 @@ function to_array(path_value::BitMatrix, instance::Instance)
 end
 
 """
-    to_array(solution::Solution, instance::Instance)
+$TYPEDSIGNATURES
 
 Returns a BitMatrix, with value true at each index (i, j) if corresponding edge of graph
 is selected in the solution
@@ -324,11 +319,11 @@ function to_array(solution::Solution, instance::Instance)
 end
 
 """
-    is_admissible(solution::Solution, instance::Instance)
+$TYPEDSIGNATURES
 
 Check if `solution` is an admissible solution of `instance`.
 """
-function is_admissible(solution::Solution, instance::Instance)
+function is_feasible(solution::Solution, instance::Instance)
     graph = instance.graph
     nb_nodes = nv(graph)
     nb_tasks = nb_nodes - 2