diff --git a/.gitignore b/.gitignore
index f65c259..2b3da69 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,6 +2,7 @@
 /docs/src/index.md
 data
 scripts
+*heuristic_algorithms
 .DS_Store
 
 # Files generated by invoking Julia with --code-coverage
diff --git a/Project.toml b/Project.toml
index c6588f0..bc396a9 100644
--- a/Project.toml
+++ b/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"
@@ -11,35 +12,43 @@ Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
-InferOpt = "4846b161-c94e-4150-8dac-c7ae193c601f"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Metalhead = "dbeba491-748d-5e0e-a39e-b530a07fa0cc"
 NPZ = "15e1cf62-19b3-5cfa-8e77-841668bca605"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+SCIP = "82193955-e24f-5292-bf16-6f2c5261a85f"
 SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [compat]
+ConstrainedShortestPaths = "0.6.0"
 DataDeps = "0.7"
 Distributions = "0.25"
 DocStringExtensions = "0.9"
-Flux = "0.14"
+Flux = "0.16"
 Graphs = "1.11"
 HiGHS = "1.9"
-Images = "0.26"
-InferOpt = "0.6"
+Images = "0.26.1"
 Ipopt = "1.6"
 JuMP = "1.22"
 LinearAlgebra = "1"
-Metalhead = "0.9"
+Metalhead = "0.9.4"
 NPZ = "0.4"
 Plots = "1"
+Printf = "1.11.0"
 Random = "1"
+Requires = "1.3.0"
 SimpleWeightedGraphs = "1.4"
 SparseArrays = "1"
+Statistics = "1.11.1"
+StatsBase = "0.34.4"
 julia = "1.6"
 
 [extras]
diff --git a/README.md b/README.md
index 6eb7f89..af77a5b 100644
--- a/README.md
+++ b/README.md
@@ -15,5 +15,6 @@ Currently, this package provides the following benchmark problems (many more to
 - `FixedSizeShortestPathBenchmark`: shortest path problem with on a graph with fixed size.
 - `WarcraftBenchmark`: shortest path problem on image maps
 - `PortfolioOptimizationBenchmark`: portfolio optimization problem.
+- `StochasticVehicleSchedulingBenchmark`: stochastic vehicle scheduling problem.
 
 See the [documentation](https://JuliaDecisionFocusedLearning.github.io/DecisionFocusedLearningBenchmarks.jl/stable/) for more details.
diff --git a/docs/make.jl b/docs/make.jl
index 16461f0..e95ff1e 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -7,13 +7,15 @@ cp(joinpath(@__DIR__, "..", "README.md"), joinpath(@__DIR__, "src", "index.md");
 md_dir = joinpath(@__DIR__, "src")
 tutorial_dir = joinpath(@__DIR__, "src", "tutorials")
 benchmarks_dir = joinpath(@__DIR__, "src", "benchmarks")
+api_dir = joinpath(@__DIR__, "src", "api")
 
+api_files = map(x -> joinpath("api", x), readdir(api_dir))
 tutorial_files = readdir(tutorial_dir)
 md_tutorial_files = [split(file, ".")[1] * ".md" for file in tutorial_files]
 benchmark_files = readdir(benchmarks_dir)
 md_benchmark_files = [split(file, ".")[1] * ".md" for file in benchmark_files]
 
-include_tutorial = true
+include_tutorial = false
 
 if include_tutorial
     for file in tutorial_files
@@ -35,9 +37,9 @@ makedocs(;
             "benchmarks/fixed_size_shortest_path.md",
             "benchmarks/warcraft.md",
             "benchmarks/portfolio_optimization.md",
+            "benchmarks/vsp.md",
         ],
-        "API reference" =>
-            ["api/interface.md", "api/decision_focused.md", "api/warcraft.md"],
+        "API reference" => api_files,
     ],
 )
 
diff --git a/docs/src/api/interface.md b/docs/src/api/0_interface.md
similarity index 100%
rename from docs/src/api/interface.md
rename to docs/src/api/0_interface.md
diff --git a/docs/src/api/decision_focused.md b/docs/src/api/decision_focused.md
deleted file mode 100644
index a41f439..0000000
--- a/docs/src/api/decision_focused.md
+++ /dev/null
@@ -1,15 +0,0 @@
-# Decisions-focused learning paper
-
-## Public
-
-```@autodocs
-Modules = [DecisionFocusedLearningBenchmarks.FixedSizeShortestPath, DecisionFocusedLearningBenchmarks.PortfolioOptimization, DecisionFocusedLearningBenchmarks.SubsetSelection]
-Private = false
-```
-
-## Private
-
-```@autodocs
-Modules = [DecisionFocusedLearningBenchmarks.FixedSizeShortestPath, DecisionFocusedLearningBenchmarks.PortfolioOptimization, DecisionFocusedLearningBenchmarks.SubsetSelection]
-Public = false
-```
diff --git a/docs/src/api/fixed_shortest_path.md b/docs/src/api/fixed_shortest_path.md
new file mode 100644
index 0000000..df50a9f
--- /dev/null
+++ b/docs/src/api/fixed_shortest_path.md
@@ -0,0 +1,15 @@
+# Fixed-size shortest path
+
+## Public
+
+```@autodocs
+Modules = [DecisionFocusedLearningBenchmarks.FixedSizeShortestPath]
+Private = false
+```
+
+## Private
+
+```@autodocs
+Modules = [DecisionFocusedLearningBenchmarks.FixedSizeShortestPath]
+Public = false
+```
diff --git a/docs/src/api/portfolio_optimization.md b/docs/src/api/portfolio_optimization.md
new file mode 100644
index 0000000..5b0102b
--- /dev/null
+++ b/docs/src/api/portfolio_optimization.md
@@ -0,0 +1,15 @@
+# Subset selection
+
+## Public
+
+```@autodocs
+Modules = [DecisionFocusedLearningBenchmarks.PortfolioOptimization]
+Private = false
+```
+
+## Private
+
+```@autodocs
+Modules = [DecisionFocusedLearningBenchmarks.PortfolioOptimization]
+Public = false
+```
diff --git a/docs/src/api/subset_selection.md b/docs/src/api/subset_selection.md
new file mode 100644
index 0000000..76b686d
--- /dev/null
+++ b/docs/src/api/subset_selection.md
@@ -0,0 +1,15 @@
+# Subset selection
+
+## Public
+
+```@autodocs
+Modules = [DecisionFocusedLearningBenchmarks.SubsetSelection]
+Private = false
+```
+
+## Private
+
+```@autodocs
+Modules = [DecisionFocusedLearningBenchmarks.SubsetSelection]
+Public = false
+```
diff --git a/docs/src/api/vsp.md b/docs/src/api/vsp.md
new file mode 100644
index 0000000..96e4cdb
--- /dev/null
+++ b/docs/src/api/vsp.md
@@ -0,0 +1,15 @@
+# Stochastic Vehicle Scheduling
+
+## Public
+
+```@autodocs
+Modules = [DecisionFocusedLearningBenchmarks.StochasticVehicleScheduling]
+Private = false
+```
+
+## Private
+
+```@autodocs
+Modules = [DecisionFocusedLearningBenchmarks.StochasticVehicleScheduling]
+Public = false
+```
diff --git a/docs/src/benchmarks/vsp.md b/docs/src/benchmarks/vsp.md
new file mode 100644
index 0000000..9c45ff2
--- /dev/null
+++ b/docs/src/benchmarks/vsp.md
@@ -0,0 +1,3 @@
+# Stochastic Vehicle Scheduling
+
+[`StochasticVehicleSchedulingBenchmark`](@ref).
diff --git a/src/DecisionFocusedLearningBenchmarks.jl b/src/DecisionFocusedLearningBenchmarks.jl
index 21265e3..2122788 100644
--- a/src/DecisionFocusedLearningBenchmarks.jl
+++ b/src/DecisionFocusedLearningBenchmarks.jl
@@ -1,10 +1,10 @@
 module DecisionFocusedLearningBenchmarks
 
 using DataDeps
-using HiGHS
-using InferOpt
+using Requires: @require
 
 function __init__()
+    # Register the Warcraft dataset
     ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"
     register(
         DataDep(
@@ -14,6 +14,10 @@ function __init__()
             post_fetch_method=unpack,
         ),
     )
+
+    # Gurobi setup
+    @info "If you have Gurobi installed and want to use it, make sure to `using Gurobi` in order to enable it."
+    @require Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b" include("gurobi_setup.jl")
     return nothing
 end
 
@@ -25,6 +29,7 @@ include("SubsetSelection/SubsetSelection.jl")
 include("Warcraft/Warcraft.jl")
 include("FixedSizeShortestPath/FixedSizeShortestPath.jl")
 include("PortfolioOptimization/PortfolioOptimization.jl")
+include("StochasticVehicleScheduling/StochasticVehicleScheduling.jl")
 
 using .Utils
 using .Argmax
@@ -33,13 +38,15 @@ using .SubsetSelection
 using .Warcraft
 using .FixedSizeShortestPath
 using .PortfolioOptimization
+using .StochasticVehicleScheduling
 
 # Interface
 export AbstractBenchmark, DataSample
 export generate_dataset
 export generate_statistical_model
-export generate_maximizer
-export plot_data
+export generate_maximizer, maximizer_kwargs
+export objective_value
+export plot_data, plot_instance, plot_solution
 export compute_gap
 
 # Export all benchmarks
@@ -49,5 +56,6 @@ export SubsetSelectionBenchmark
 export WarcraftBenchmark
 export FixedSizeShortestPathBenchmark
 export PortfolioOptimizationBenchmark
+export StochasticVehicleSchedulingBenchmark
 
 end # module DecisionFocusedLearningBenchmarks
diff --git a/src/FixedSizeShortestPath/FixedSizeShortestPath.jl b/src/FixedSizeShortestPath/FixedSizeShortestPath.jl
index 7eea5a3..fd60de2 100644
--- a/src/FixedSizeShortestPath/FixedSizeShortestPath.jl
+++ b/src/FixedSizeShortestPath/FixedSizeShortestPath.jl
@@ -54,6 +54,12 @@ function FixedSizeShortestPathBenchmark(;
     return FixedSizeShortestPathBenchmark(g, grid_size, p, deg, ν)
 end
 
+function Utils.objective_value(
+    ::FixedSizeShortestPathBenchmark, θ::AbstractArray, y::AbstractArray
+)
+    return -dot(θ, y)
+end
+
 """
 $TYPEDSIGNATURES
 
@@ -149,5 +155,6 @@ function Utils.generate_statistical_model(bench::FixedSizeShortestPathBenchmark)
 end
 
 export FixedSizeShortestPathBenchmark
+export generate_dataset, generate_maximizer, generate_statistical_model
 
 end
diff --git a/src/PortfolioOptimization/PortfolioOptimization.jl b/src/PortfolioOptimization/PortfolioOptimization.jl
index b881909..308770a 100644
--- a/src/PortfolioOptimization/PortfolioOptimization.jl
+++ b/src/PortfolioOptimization/PortfolioOptimization.jl
@@ -2,12 +2,12 @@ module PortfolioOptimization
 
 using ..Utils
 using DocStringExtensions: TYPEDEF, TYPEDFIELDS, TYPEDSIGNATURES
-using Distributions
+using Distributions: Uniform, Bernoulli
 using Flux: Chain, Dense
-using Ipopt
-using JuMP
-using LinearAlgebra
-using Random
+using Ipopt: Ipopt
+using JuMP: @variable, @objective, @constraint, optimize!, value, Model, set_silent
+using LinearAlgebra: I
+using Random: MersenneTwister
 
 """
 $TYPEDEF
@@ -127,5 +127,6 @@ function Utils.generate_statistical_model(bench::PortfolioOptimizationBenchmark)
 end
 
 export PortfolioOptimizationBenchmark
+export generate_dataset, generate_maximizer, generate_statistical_model
 
 end
diff --git a/src/StochasticVehicleScheduling/StochasticVehicleScheduling.jl b/src/StochasticVehicleScheduling/StochasticVehicleScheduling.jl
new file mode 100644
index 0000000..c864225
--- /dev/null
+++ b/src/StochasticVehicleScheduling/StochasticVehicleScheduling.jl
@@ -0,0 +1,234 @@
+module StochasticVehicleScheduling
+
+export StochasticVehicleSchedulingBenchmark
+export generate_dataset, generate_maximizer, generate_statistical_model
+export plot_instance, plot_solution
+export compact_linearized_mip, compact_mip, column_generation_algorithm, local_search
+export evaluate_solution, is_feasible
+
+using ..Utils
+using DocStringExtensions: TYPEDEF, TYPEDFIELDS, TYPEDSIGNATURES
+using ConstrainedShortestPaths:
+    stochastic_routing_shortest_path, stochastic_routing_shortest_path_with_threshold
+using Distributions: Distribution, LogNormal, Uniform, DiscreteUniform
+using Flux: Chain, Dense
+using Graphs:
+    AbstractGraph,
+    SimpleDiGraph,
+    add_edge!,
+    nv,
+    ne,
+    edges,
+    src,
+    dst,
+    has_edge,
+    inneighbors,
+    outneighbors
+using JuMP:
+    JuMP, Model, @variable, @objective, @constraint, optimize!, value, set_silent, dual
+using Plots: Plots, plot, plot!, scatter!, annotate!, text
+using Printf: @printf
+using Random: Random, AbstractRNG, MersenneTwister
+using SparseArrays: sparse
+using Statistics: quantile, mean
+
+include("utils.jl")
+include("instance/constants.jl")
+include("instance/task.jl")
+include("instance/district.jl")
+include("instance/city.jl")
+include("instance/features.jl")
+include("instance/instance.jl")
+
+include("solution/solution.jl")
+include("solution/algorithms/mip.jl")
+include("solution/algorithms/column_generation.jl")
+include("solution/algorithms/local_search.jl")
+
+include("maximizer.jl")
+
+"""
+$TYPEDEF
+
+Data structure for a stochastic vehicle scheduling benchmark.
+
+# Fields
+$TYPEDFIELDS
+"""
+@kwdef struct StochasticVehicleSchedulingBenchmark <: AbstractBenchmark
+    "number of tasks in each instance"
+    nb_tasks::Int = 25
+    "number of scenarios in each instance"
+    nb_scenarios::Int = 10
+end
+
+function Utils.objective_value(
+    ::StochasticVehicleSchedulingBenchmark, sample::DataSample, y::BitVector
+)
+    return evaluate_solution(y, sample.instance)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Create a dataset of `dataset_size` instances for the given `StochasticVehicleSchedulingBenchmark`.
+If you want to not add label solutions in the dataset, set `compute_solutions=false`.
+By default, they will be computed using column generation.
+Note that computing solutions can be time-consuming, especially for large instances.
+You can also use instead `compact_mip` or `compact_linearized_mip` as the algorithm to compute solutions.
+If you want to provide a custom algorithm to compute solutions, you can pass it as the `algorithm` keyword argument.
+If `algorithm` takes keyword arguments, you can pass them as well directly in `kwargs...`.
+If `store_city=false`, the coordinates and unnecessary information about instances will not be stored in the dataset.
+"""
+function Utils.generate_dataset(
+    benchmark::StochasticVehicleSchedulingBenchmark,
+    dataset_size::Int;
+    compute_solutions=true,
+    seed=nothing,
+    rng=MersenneTwister(0),
+    algorithm=column_generation_algorithm,
+    store_city=true,
+    kwargs...,
+)
+    (; nb_tasks, nb_scenarios) = benchmark
+    Random.seed!(rng, seed)
+    instances = [
+        Instance(; nb_tasks, nb_scenarios, rng, store_city) for _ in 1:dataset_size
+    ]
+    features = get_features.(instances)
+    if compute_solutions
+        solutions = [algorithm(instance; kwargs...).value for instance in instances]
+        return [
+            DataSample(; x=feature, instance, y_true=solution) for
+            (instance, feature, solution) in zip(instances, features, solutions)
+        ]
+    end
+    # else
+    return [
+        DataSample(; x=feature, instance) for
+        (instance, feature) in zip(instances, features)
+    ]
+end
+
+"""
+$TYPEDSIGNATURES
+"""
+function Utils.generate_maximizer(bench::StochasticVehicleSchedulingBenchmark)
+    return vsp_maximizer
+end
+
+"""
+$TYPEDSIGNATURES
+"""
+function Utils.generate_statistical_model(
+    ::StochasticVehicleSchedulingBenchmark; seed=nothing
+)
+    Random.seed!(seed)
+    return Chain(Dense(20 => 1; bias=false), vec)
+end
+
+"""
+$TYPEDSIGNATURES
+"""
+function plot_instance(
+    ::StochasticVehicleSchedulingBenchmark,
+    sample::DataSample{<:Instance{City}};
+    color_scheme=:lightrainbow,
+    kwargs...,
+)
+    (; tasks, district_width, width) = sample.instance.city
+    ticks = 0:district_width:width
+    max_time = maximum(t.end_time for t in sample.instance.city.tasks[1:(end - 1)])
+    fig = plot(;
+        xlabel="x",
+        ylabel="y",
+        gridlinewidth=3,
+        aspect_ratio=:equal,
+        size=(500, 500),
+        xticks=ticks,
+        yticks=ticks,
+        xlims=(-1, width + 1),
+        ylims=(-1, width + 1),
+        clim=(0.0, max_time),
+        label=nothing,
+        colorbar_title="Time",
+    )
+    scatter!(
+        fig,
+        [tasks[1].start_point.x],
+        [tasks[1].start_point.y];
+        label=nothing,
+        marker=:rect,
+        markersize=10,
+    )
+    annotate!(fig, (tasks[1].start_point.x, tasks[1].start_point.y, text("0", 10)))
+    for (i_task, task) in enumerate(tasks[2:(end - 1)])
+        (; start_point, end_point) = task
+        points = [(start_point.x, start_point.y), (end_point.x, end_point.y)]
+        plot!(fig, points; color=:black, label=nothing)
+        scatter!(
+            fig,
+            points[1];
+            markersize=10,
+            marker=:rect,
+            marker_z=task.start_time,
+            colormap=:turbo,
+            label=nothing,
+        )
+        scatter!(
+            fig,
+            points[2];
+            markersize=10,
+            marker=:rect,
+            marker_z=task.end_time,
+            colormap=:turbo,
+            label=nothing,
+            # color=palette[max(floor(Int, task.end_time), 1)],
+        )
+        annotate!(fig, (points[1]..., text("$(i_task)", 10)))
+    end
+    return fig
+end
+
+"""
+$TYPEDSIGNATURES
+"""
+function plot_solution(
+    ::StochasticVehicleSchedulingBenchmark, sample::DataSample{<:Instance{City}}; kwargs...
+)
+    (; tasks, district_width, width) = sample.instance.city
+    ticks = 0:district_width:width
+    solution = Solution(sample.y_true, sample.instance)
+    path_list = compute_path_list(solution)
+    fig = plot(;
+        xlabel="x",
+        ylabel="y",
+        legend=false,
+        gridlinewidth=3,
+        aspect_ratio=:equal,
+        size=(500, 500),
+        xticks=ticks,
+        yticks=ticks,
+        xlims=(-1, width + 1),
+        ylims=(-1, width + 1),
+    )
+    for path in path_list
+        X = Float64[]
+        Y = Float64[]
+        (; start_point, end_point) = tasks[path[1]]
+        (; x, y) = end_point
+        push!(X, x)
+        push!(Y, y)
+        for task in path[2:end]
+            (; start_point, end_point) = tasks[task]
+            push!(X, start_point.x)
+            push!(Y, start_point.y)
+            push!(X, end_point.x)
+            push!(Y, end_point.y)
+        end
+        plot!(fig, X, Y; marker=:circle)
+    end
+    return fig
+end
+
+end
diff --git a/src/StochasticVehicleScheduling/instance/city.jl b/src/StochasticVehicleScheduling/instance/city.jl
new file mode 100644
index 0000000..5f962b0
--- /dev/null
+++ b/src/StochasticVehicleScheduling/instance/city.jl
@@ -0,0 +1,263 @@
+"""
+$TYPEDEF
+
+Data structure for a city in the vehicle scheduling problem.
+Contains all the relevant information to build an instance of the problem.
+
+# Fields
+$TYPEDFIELDS
+"""
+struct City
+    "city width (in minutes)"
+    width::Int
+    # Objectives ponderation
+    "cost of a vehicle in the objective function"
+    vehicle_cost::Float64
+    "cost of one minute delay in the objective function"
+    delay_cost::Float64
+    # Tasks
+    "number of tasks to fulfill"
+    nb_tasks::Int
+    "tasks list (see [`Task`](@ref)), that should be ordered by start time"
+    tasks::Vector{Task}
+    # Stochastic specific stuff
+    "idth (in minutes) of each district"
+    district_width::Int
+    "districts matrix (see [`District`](@ref)), indices corresponding to their relative positions"
+    districts::Matrix{District}
+    "a log-normal distribution modeling delay between districts"
+    random_inter_area_factor::LogNormal{Float64}
+    "size (nb_scenarios, 24), each row correspond to one scenario, each column to one hour of the day"
+    scenario_inter_area_factor::Matrix{Float64}
+end
+
+"""
+$TYPEDSIGNATURES
+
+Constructor for [`City`](@ref).
+"""
+function City(;
+    nb_scenarios=default_nb_scenarios,
+    width=default_width,
+    vehicle_cost=default_vehicle_cost,
+    nb_tasks=default_nb_tasks,
+    tasks=Vector{Task}(undef, nb_tasks + 2),
+    district_width=default_district_width,
+    districts=Matrix{District}(undef, width ÷ district_width, width ÷ district_width),
+    delay_cost=default_delay_cost,
+    random_inter_area_factor=default_random_inter_area_factor,
+    scenario_inter_area_factor=zeros(nb_scenarios, 24),
+)
+    return City(
+        width,
+        vehicle_cost,
+        delay_cost,
+        nb_tasks,
+        tasks,
+        district_width,
+        districts,
+        random_inter_area_factor,
+        scenario_inter_area_factor,
+    )
+end
+
+"""
+$TYPEDSIGNATURES
+
+- Creates a city from `city_kwargs`
+- Depot location at city center
+- Randomize tasks, and add two dummy tasks : one `source` task at time=0 from the depot,
+    and one `destination` task ending at time=end at depot
+- Roll every scenario.
+"""
+function create_random_city(;
+    αᵥ_low=default_αᵥ_low,
+    αᵥ_high=default_αᵥ_high,
+    first_begin_time=default_first_begin_time,
+    last_begin_time=default_last_begin_time,
+    district_μ=default_district_μ,
+    district_σ=default_district_σ,
+    task_μ=default_task_μ,
+    task_σ=default_task_σ,
+    seed=nothing,
+    rng=Random.default_rng(),
+    city_kwargs...,
+)
+    isnothing(seed) || Random.seed!(rng, seed)
+    city = City(; city_kwargs...)
+    init_districts!(city, district_μ, district_σ; rng=rng)
+    init_tasks!(
+        city, αᵥ_low, αᵥ_high, first_begin_time, last_begin_time, task_μ, task_σ; rng=rng
+    )
+    generate_scenarios!(city; rng=rng)
+    compute_perturbed_end_times!(city)
+    return city
+end
+
+"""
+$TYPEDSIGNATURES
+
+Initialize the districts of the city.
+"""
+function init_districts!(
+    city::City, district_μ::Distribution, district_σ::Distribution; rng::AbstractRNG
+)
+    nb_scenarios = size(city.scenario_inter_area_factor, 1)
+    nb_district_per_edge = city.width ÷ city.district_width
+    for x in 1:nb_district_per_edge
+        for y in 1:nb_district_per_edge
+            μ = rand(rng, district_μ)
+            σ = rand(rng, district_σ)
+            city.districts[x, y] = District(;
+                random_delay=LogNormal(μ, σ), nb_scenarios=nb_scenarios
+            )
+        end
+    end
+    return nothing
+end
+
+"""
+$TYPEDSIGNATURES
+
+Draw the tasks of the city.
+"""
+function init_tasks!(
+    city::City,
+    αᵥ_low::Real,
+    αᵥ_high::Real,
+    first_begin_time::Real,
+    last_begin_time::Real,
+    task_μ::Distribution,
+    task_σ::Distribution;
+    rng::AbstractRNG,
+)
+    nb_scenarios = size(city.scenario_inter_area_factor, 1)
+
+    point_distribution = Uniform(0, city.width)
+    start_time_distribution = Uniform(first_begin_time, last_begin_time)
+    travel_time_multiplier_distribution = Uniform(αᵥ_low, αᵥ_high)
+
+    for i_task in 1:(city.nb_tasks)
+        start_point = draw_random_point(point_distribution; rng=rng)
+        end_point = draw_random_point(point_distribution; rng=rng)
+
+        start_time = rand(rng, start_time_distribution)
+        end_time =
+            start_time +
+            rand(rng, travel_time_multiplier_distribution) *
+            distance(start_point, end_point)
+
+        μ = rand(rng, task_μ)
+        σ = rand(rng, task_σ)
+        random_delay = LogNormal(μ, σ)
+
+        city.tasks[i_task + 1] = Task(;
+            type=job::TaskType,
+            start_point=start_point,
+            end_point=end_point,
+            start_time=start_time,
+            end_time=end_time,
+            random_delay=random_delay,
+            nb_scenarios=nb_scenarios,
+        )
+    end
+
+    # add start and final "artificial" tasks
+    city_center = Point(city.width / 2, city.width / 2)  # ? hard coded ?
+    city.tasks[1] = Task(;
+        type=depot_start::TaskType,
+        start_point=city_center,
+        end_point=city_center,
+        start_time=0.0,
+        end_time=0.0,
+        random_delay=ZERO_UNIFORM,
+        nb_scenarios=nb_scenarios,
+    )
+    final_task_time = 24 * 60.0 # ? hard coded ?
+    city.tasks[end] = Task(;
+        type=depot_end::TaskType,
+        start_point=city_center,
+        end_point=city_center,
+        start_time=final_task_time,
+        end_time=final_task_time,
+        random_delay=ZERO_UNIFORM,
+        nb_scenarios=nb_scenarios,
+    )
+
+    # sort tasks by start time
+    sort!(city.tasks; by=task -> task.start_time, rev=false)
+    return nothing
+end
+
+"""
+    get_district(point::Point, city::City)
+
+Return indices of the `city` district containing `point`.
+"""
+function get_district(point::Point, city::City)
+    return trunc(Int, point.x / city.district_width) + 1,
+    trunc(Int, point.y / city.district_width) + 1
+end
+
+"""
+$TYPEDSIGNATURES
+
+Draw all delay scenarios for the city.
+"""
+function generate_scenarios!(city::City; rng::AbstractRNG)
+    # roll all tasks
+    for task in city.tasks
+        roll(task, rng)
+    end
+
+    # roll all districts
+    for district in city.districts
+        roll(district, rng)
+    end
+
+    # roll inter-district
+    nb_scenarios, nb_hours = size(city.scenario_inter_area_factor)
+    for s in 1:nb_scenarios
+        previous_delay = 0.0
+        for h in 1:nb_hours
+            previous_delay =
+                (previous_delay + 0.1) * rand(rng, city.random_inter_area_factor)
+            city.scenario_inter_area_factor[s, h] = previous_delay
+        end
+    end
+    return nothing
+end
+
+"""
+$TYPEDSIGNATURES
+
+Compute the end times of the tasks for each scenario.
+"""
+function compute_perturbed_end_times!(city::City)
+    nb_scenarios = size(city.scenario_inter_area_factor, 1)
+
+    for task in city.tasks[2:(end - 1)]
+        start_time = task.start_time
+        end_time = task.end_time
+        start_point = task.start_point
+        end_point = task.end_point
+
+        origin_x, origin_y = get_district(start_point, city)
+        destination_x, destination_y = get_district(end_point, city)
+        origin_district = city.districts[origin_x, origin_y]
+        destination_district = city.districts[destination_x, destination_y]
+
+        scenario_start_time = task.scenario_start_time
+        origin_delay = origin_district.scenario_delay
+        destination_delay = destination_district.scenario_delay
+        inter_area_delay = city.scenario_inter_area_factor
+
+        for s in 1:nb_scenarios
+            ξ₁ = scenario_start_time[s]
+            ξ₂ = ξ₁ + origin_delay[s, hour_of(ξ₁)]
+            ξ₃ = ξ₂ + end_time - start_time + inter_area_delay[s, hour_of(ξ₂)]
+            task.scenario_end_time[s] = ξ₃ + destination_delay[s, hour_of(ξ₃)]
+        end
+    end
+    return nothing
+end
diff --git a/src/StochasticVehicleScheduling/instance/constants.jl b/src/StochasticVehicleScheduling/instance/constants.jl
new file mode 100644
index 0000000..276c868
--- /dev/null
+++ b/src/StochasticVehicleScheduling/instance/constants.jl
@@ -0,0 +1,23 @@
+# --- Default values for structure attributes ---
+
+const ZERO_UNIFORM = LogNormal(-Inf, 1.0)  # always returns 0 (useful for type stability reasons)
+
+const default_delay_cost = 2.0             # cost of one minute of delay
+const default_vehicle_cost = 1000.0        # cost of one vehicle
+const default_width = 50                   # width (in minutes) of the squared city
+const default_αᵥ_low = 1.2                 # used for drawing random tasks
+const default_αᵥ_high = 1.6                # used for drawing random tasks
+const default_first_begin_time = 60.0 * 6  # Start of time window at 6AM
+const default_last_begin_time = 60.0 * 20  # End of time window at 8PM
+
+const default_district_width = 10          # width (in minutes) of each squared district
+const default_random_inter_area_factor = LogNormal(0.02, 0.05)
+
+const default_district_μ = Uniform(0.8, 1.2)
+const default_district_σ = Uniform(0.4, 0.6)
+
+const default_task_μ = Uniform(1.0, 3.0)
+const default_task_σ = Uniform(0.8, 1.2)
+
+const default_nb_tasks = 10               # Number of tasks in an instnace
+const default_nb_scenarios = 1            # Number of scenrios in an instance
diff --git a/src/StochasticVehicleScheduling/instance/district.jl b/src/StochasticVehicleScheduling/instance/district.jl
new file mode 100644
index 0000000..864a4c3
--- /dev/null
+++ b/src/StochasticVehicleScheduling/instance/district.jl
@@ -0,0 +1,55 @@
+"""
+$TYPEDEF
+
+Data structure for a district in the vehicle scheduling problem.
+
+# Fields
+$TYPEDFIELDS
+"""
+struct District
+    "log-normal distribution modeling the district delay"
+    random_delay::LogNormal{Float64}
+    "size (nb_scenarios, 24), observed delays for each scenario and hour of the day"
+    scenario_delay::Matrix{Float64}
+end
+
+"""
+$TYPEDSIGNATURES
+
+Constructor for [`District`](@ref).
+Initialize a district with a given number of scenarios, with zeros in `scenario_delay`.
+"""
+function District(; random_delay::LogNormal{Float64}, nb_scenarios::Int)
+    return District(random_delay, zeros(nb_scenarios, 24))
+end
+
+"""
+$TYPEDSIGNATURES
+
+Return one scenario of future delay given current delay and delay distribution.
+"""
+function scenario_next_delay(
+    previous_delay::Real, random_delay::Distribution, rng::AbstractRNG
+)
+    return previous_delay / 2.0 + rand(rng, random_delay)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Populate `scenario_delay` with delays drawn from `random_delay` distribution
+for each (scenario, hour) pair.
+"""
+function roll(district::District, rng::AbstractRNG)
+    nb_scenarios, nb_hours = size(district.scenario_delay)
+    # Loop on scenarios
+    for s in 1:nb_scenarios
+        previous_delay = 0.0
+        # Loop on hours
+        for h in 1:nb_hours
+            previous_delay = scenario_next_delay(previous_delay, district.random_delay, rng)
+            district.scenario_delay[s, h] = previous_delay
+        end
+    end
+    return nothing
+end
diff --git a/src/StochasticVehicleScheduling/instance/features.jl b/src/StochasticVehicleScheduling/instance/features.jl
new file mode 100644
index 0000000..81a807d
--- /dev/null
+++ b/src/StochasticVehicleScheduling/instance/features.jl
@@ -0,0 +1,137 @@
+"""
+$TYPEDSIGNATURES
+
+Compute the achieved travel time of scenario `scenario` from `old_task_index` to `new_task_index`.
+"""
+function get_perturbed_travel_time(
+    city::City, old_task_index::Int, new_task_index::Int, scenario::Int
+)
+    old_task = city.tasks[old_task_index]
+    new_task = city.tasks[new_task_index]
+
+    origin_x, origin_y = get_district(old_task.end_point, city)
+    destination_x, destination_y = get_district(new_task.start_point, city)
+
+    ξ₁ = old_task.scenario_end_time[scenario]
+    ξ₂ = ξ₁ + city.districts[origin_x, origin_y].scenario_delay[scenario, hour_of(ξ₁)]
+    ξ₃ =
+        ξ₂ +
+        distance(old_task.end_point, new_task.start_point) +
+        city.scenario_inter_area_factor[scenario, hour_of(ξ₂)]
+    return ξ₃ + city.districts[destination_x, destination_y].scenario_delay[
+        scenario, hour_of(ξ₃)
+    ] - ξ₁
+end
+
+"""
+$TYPEDSIGNATURES
+
+Compute slack for features.
+"""
+function compute_slacks(city::City, old_task_index::Int, new_task_index::Int)
+    old_task = city.tasks[old_task_index]
+    new_task = city.tasks[new_task_index]
+
+    travel_time = distance(old_task.end_point, new_task.start_point)
+    perturbed_end_times = old_task.scenario_end_time
+    perturbed_start_times = new_task.scenario_start_time
+
+    return perturbed_start_times .- (perturbed_end_times .+ travel_time)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Returns a matrix of features of size (20, nb_edges).
+For each edge, compute the following features (in the same order):
+- travel time
+- vehicle_cost if edge is connected to source, else 0
+- 9 deciles of the slack
+- cumulative probability distribution of the slack evaluated in [-100, -50, -20, -10, 0, 10, 50, 200, 500]
+"""
+function compute_features(city::City)
+    graph = create_VSP_graph(city)
+
+    cumul = [-100, -50, -20, -10, 0, 10, 50, 200, 500]
+    nb_features = 2 + 9 + length(cumul)
+    features = zeros(Float32, nb_features, ne(graph))
+
+    # features indices
+    travel_time_index = 1
+    connected_to_source_index = 2
+    slack_deciles_indices = 3:11
+    slack_cumulative_distribution_indices = 12:nb_features
+
+    for (i, edge) in enumerate(edges(graph))
+        # compute travel time
+        features[travel_time_index, i] = distance(
+            city.tasks[src(edge)].end_point, city.tasks[dst(edge)].start_point
+        )
+        # if edge connected to source node
+        features[connected_to_source_index, i] = src(edge) == 1 ? city.vehicle_cost : 0.0
+
+        # slack related features
+        slacks = compute_slacks(city, src(edge), dst(edge))
+        # compute deciles
+        features[slack_deciles_indices, i] = quantile(slacks, [0.1 * i for i in 1:9])
+        # compute cumulative distribution
+        features[slack_cumulative_distribution_indices, i] = [
+            mean(slacks .<= x) for x in cumul
+        ]
+    end
+    return features
+end
+
+"""
+$TYPEDSIGNATURES
+
+Compute slack for instance.
+TODO: differentiate from other method
+"""
+function compute_slacks(city::City, graph::AbstractGraph)
+    (; tasks) = city
+    N = nv(graph)
+    slack_list = [
+        [
+            (dst(e) < N ? tasks[dst(e)].scenario_start_time[ω] : Inf) -
+            (tasks[src(e)].end_time + get_perturbed_travel_time(city, src(e), dst(e), ω))
+            for ω in 1:get_nb_scenarios(city)
+        ] for e in edges(graph)
+    ]
+    I = [src(e) for e in edges(graph)]
+    J = [dst(e) for e in edges(graph)]
+    return sparse(I, J, slack_list)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Compute delays for instance.
+"""
+function compute_delays(city::City)
+    nb_tasks = get_nb_tasks(city)
+    nb_scenarios = get_nb_scenarios(city)
+    ε = zeros(nb_tasks, nb_scenarios)
+    for (index, task) in enumerate(city.tasks)
+        ε[index, :] .= task.scenario_end_time .- task.end_time
+    end
+    return ε
+end
+
+"""
+$TYPEDSIGNATURES
+
+Returns the number of tasks in city.
+"""
+function get_nb_tasks(city::City)
+    return length(city.tasks)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Returns the number of scenarios in city.
+"""
+function get_nb_scenarios(city::City)
+    return size(city.scenario_inter_area_factor, 1)
+end
diff --git a/src/StochasticVehicleScheduling/instance/instance.jl b/src/StochasticVehicleScheduling/instance/instance.jl
new file mode 100644
index 0000000..090c69a
--- /dev/null
+++ b/src/StochasticVehicleScheduling/instance/instance.jl
@@ -0,0 +1,117 @@
+"""
+$TYPEDEF
+
+Instance of the stochastic VSP problem.
+
+# Fields
+$TYPEDFIELDS
+"""
+struct Instance{CC,G<:AbstractGraph,M1<:AbstractMatrix,M2<:AbstractMatrix,F,C}
+    "graph computed from `city` with the `create_VSP_graph(city::City)` method"
+    graph::G
+    "features matrix computed from `city`"
+    features::Matrix{F}
+    "slack matrix"
+    slacks::M1
+    "intrinsic delays scenario matrix"
+    intrinsic_delays::M2
+    "cost of a vehicle"
+    vehicle_cost::C
+    "cost of one minute delay"
+    delay_cost::C
+    "associated city"
+    city::CC
+end
+
+"""
+$TYPEDSIGNATURES
+
+Return the acyclic directed graph corresponding to `city`.
+Each vertex represents a task. Vertices are ordered by start time of corresponding task.
+There is an edge from task u to task v the (end time of u + tie distance between u and v <= start time of v).
+"""
+function create_VSP_graph(city::City)
+    # Initialize directed graph
+    nb_vertices = city.nb_tasks + 2
+    graph = SimpleDiGraph(nb_vertices)
+    starting_task = 1
+    end_task = nb_vertices
+    job_tasks = 2:(city.nb_tasks + 1)
+
+    travel_times = [
+        distance(task1.end_point, task2.start_point) for task1 in city.tasks,
+        task2 in city.tasks
+    ]
+
+    # Create existing edges
+    for iorigin in job_tasks
+        # link every task to base
+        add_edge!(graph, starting_task, iorigin)
+        add_edge!(graph, iorigin, end_task)
+
+        for idestination in (iorigin + 1):(city.nb_tasks + 1)
+            travel_time = travel_times[iorigin, idestination]
+            origin_end_time = city.tasks[iorigin].end_time
+            destination_begin_time = city.tasks[idestination].start_time # get_prop(graph, idestination, :task).start_time
+
+            # there is an edge only if we can reach destination from origin before start of task
+            if origin_end_time + travel_time <= destination_begin_time
+                add_edge!(graph, iorigin, idestination)
+            end
+        end
+    end
+
+    return graph
+end
+
+"""
+$TYPEDSIGNATURES
+
+Constructor for [`Instance`](@ref).
+Build an `Instance` for the stochatsic vehicle scheduling problem, with `nb_tasks` tasks and `nb_scenarios` scenarios.
+"""
+function Instance(;
+    nb_tasks::Int,
+    nb_scenarios::Int,
+    rng::AbstractRNG=Random.default_rng(),
+    store_city=true,
+    kwargs...,
+)
+    city = create_random_city(; rng=rng, nb_tasks, nb_scenarios, kwargs...)
+    graph = create_VSP_graph(city)
+    features = compute_features(city)
+    slacks = compute_slacks(city, graph)
+    intrinsic_delays = compute_delays(city)
+    return Instance(
+        graph,
+        features,
+        slacks,
+        intrinsic_delays,
+        city.vehicle_cost,
+        city.delay_cost,
+        store_city ? city : nothing,
+    )
+end
+
+"""
+$TYPEDSIGNATURES
+
+Returns the number of scenarios in instance.
+"""
+function get_nb_scenarios(instance::Instance)
+    return size(instance.intrinsic_delays, 2)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Returns the number of tasks in `instance`.
+"""
+get_nb_tasks(instance::Instance) = nv(instance.graph) - 2
+
+"""
+$TYPEDSIGNATURES
+
+Returns the feature matrix associated to `instance`.
+"""
+get_features(instance::Instance) = instance.features
diff --git a/src/StochasticVehicleScheduling/instance/task.jl b/src/StochasticVehicleScheduling/instance/task.jl
new file mode 100644
index 0000000..0479c61
--- /dev/null
+++ b/src/StochasticVehicleScheduling/instance/task.jl
@@ -0,0 +1,90 @@
+@enum TaskType depot_start job depot_end
+
+"""
+$TYPEDEF
+
+Data structure for a task in the vehicle scheduling problem.
+
+# Fields
+$TYPEDFIELDS
+"""
+struct Task
+    "type of the task (depot start, job, or depot end)"
+    type::TaskType
+    "starting location of the task"
+    start_point::Point
+    "end location of the task"
+    end_point::Point
+    "start time (in minutes) of the task"
+    start_time::Float64
+    "end time (in minutes) of the task"
+    end_time::Float64
+    "lognormal distribution modeling the task start delay"
+    random_delay::LogNormal{Float64}
+    "size (nb_scenarios), observed delayed start times for each scenario"
+    scenario_start_time::Vector{Float64}
+    "size (nb_scenarios), observed delayed end times for each scenario"
+    scenario_end_time::Vector{Float64}
+end
+
+"""
+$TYPEDSIGNATURES
+
+Constructor for [`Task`](@ref).
+"""
+function Task(;
+    type::TaskType=job,
+    start_point::Point,
+    end_point::Point,
+    start_time::Float64,
+    end_time::Float64,
+    nb_scenarios::Int,
+    random_delay::LogNormal{Float64}=ZERO_UNIFORM,
+)
+    return Task(
+        type,
+        start_point,
+        end_point,
+        start_time,
+        end_time,
+        random_delay,
+        zeros(nb_scenarios) .+ start_time,
+        zeros(nb_scenarios) .+ end_time,
+    )
+end
+
+"""
+$TYPEDSIGNATURES
+
+Return the number of scenarios for the given task.
+"""
+function nb_scenarios(task::Task)
+    return length(task.scenario_start_time)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Populate `scenario_start_time` with delays drawn from the `random_delay` distribution of
+the given task for each scenario.
+"""
+function roll(task::Task, rng::AbstractRNG)
+    S = nb_scenarios(task)
+    task.scenario_start_time .= task.start_time .+ rand(rng, task.random_delay, S)
+    return nothing
+end
+
+function Base.show(io::IO, task::Task)
+    @printf(
+        "(%.2f, %.2f) -> (%.2f, %.2f), [%.2f, %.2f], %s, %s, %s",
+        task.start_point.x,
+        task.start_point.y,
+        task.end_point.x,
+        task.end_point.y,
+        task.start_time,
+        task.end_time,
+        task.type,
+        task.scenario_start_time,
+        task.scenario_end_time,
+    )
+end
diff --git a/src/StochasticVehicleScheduling/maximizer.jl b/src/StochasticVehicleScheduling/maximizer.jl
new file mode 100644
index 0000000..66e9eb3
--- /dev/null
+++ b/src/StochasticVehicleScheduling/maximizer.jl
@@ -0,0 +1,45 @@
+"""
+$TYPEDSIGNATURES
+
+Given arcs weights θ, solve the deterministic VSP problem associated to `instance`.
+"""
+function vsp_maximizer(
+    θ::AbstractVector; instance::Instance, model_builder=highs_model, silent=true
+)
+    (; graph) = instance
+
+    model = model_builder()
+    silent && set_silent(model)
+
+    nb_nodes = nv(graph)
+    job_indices = 2:(nb_nodes - 1)
+
+    @variable(model, y[i=1:nb_nodes, j=1:nb_nodes; has_edge(graph, i, j)], Bin)
+
+    @objective(
+        model,
+        Max,
+        sum(θ[i] * y[src(edge), dst(edge)] for (i, edge) in enumerate(edges(graph)))
+    )
+
+    @constraint(
+        model,
+        flow[i in job_indices],
+        sum(y[j, i] for j in inneighbors(graph, i)) ==
+            sum(y[i, j] for j in outneighbors(graph, i))
+    )
+    @constraint(
+        model, demand[i in job_indices], sum(y[j, i] for j in inneighbors(graph, i)) == 1
+    )
+
+    optimize!(model)
+
+    solution = falses(ne(graph))
+    for (i, edge) in enumerate(edges(graph))
+        if isapprox(value(y[edge.src, edge.dst]), 1; atol=1e-3)
+            solution[i] = true
+        end
+    end
+
+    return solution
+end
diff --git a/src/StochasticVehicleScheduling/solution/algorithms/column_generation.jl b/src/StochasticVehicleScheduling/solution/algorithms/column_generation.jl
new file mode 100644
index 0000000..dbd2fd6
--- /dev/null
+++ b/src/StochasticVehicleScheduling/solution/algorithms/column_generation.jl
@@ -0,0 +1,193 @@
+"""
+    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;
+    model_builder=highs_model,
+    bounding,
+    use_convex_resources,
+    scenario_range=nothing,
+    silent=true,
+)
+    (; graph, vehicle_cost, delay_cost) = instance
+
+    Ω = isnothing(scenario_range) ? (1:get_nb_scenarios(instance)) : scenario_range
+    intrinsic_delays = instance.intrinsic_delays[:, Ω]
+    slacks = deepcopy(instance.slacks)
+    for edge in edges(graph)
+        slacks[src(edge), dst(edge)] = slacks[src(edge), dst(edge)][Ω]
+    end
+
+    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;
+            bounding,
+            use_convex_resources,
+        )
+        λ_sum = sum(λ_val[v] for v in job_indices if v in p_star)
+        path_cost = delay_cost * c_star + λ_sum + vehicle_cost
+        silent || @info path_cost - λ_sum
+        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 = JuMP.objective_value(model)
+    columns = unique(cat(initial_paths, new_paths; dims=1))
+    return columns, c_low::Float64, value.(λ)
+end
+
+"""
+$TYPEDSIGNATURES
+
+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,
+    scenario_range=nothing,
+    silent=true,
+)
+    (; graph, vehicle_cost, delay_cost) = instance
+
+    nb_nodes = nv(graph)
+    job_indices = 2:(nb_nodes - 1)
+
+    Ω = isnothing(scenario_range) ? (1:get_nb_scenarios(instance)) : scenario_range
+    intrinsic_delays = instance.intrinsic_delays[:, Ω]
+    slacks = deepcopy(instance.slacks)
+    for edge in edges(graph)
+        slacks[src(edge), dst(edge)] = slacks[src(edge), dst(edge)][Ω]
+    end
+
+    model = model_builder()
+    silent && 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 JuMP.objective_value(model)::Float64,
+    sol,
+    paths[isapprox.([sol[p] for p in paths], 1.0)]
+end
+
+"""
+$TYPEDSIGNATURES
+
+Solve input instance using column generation.
+"""
+function column_generation_algorithm(
+    instance::Instance;
+    scenario_range=nothing,
+    bounding=true,
+    use_convex_resources=true,
+    silent=true,
+    close_gap=false,
+)
+    Ω = isnothing(scenario_range) ? (1:get_nb_scenarios(instance)) : scenario_range
+    columns, c_low, λ_val = column_generation(
+        instance;
+        bounding=bounding,
+        use_convex_resources=use_convex_resources,
+        scenario_range=Ω,
+        silent=silent,
+    )
+    c_upp, _, sol = compute_solution_from_selected_columns(
+        instance, columns; scenario_range=Ω, silent=silent
+    )
+
+    if close_gap && abs(c_upp - c_low) > 1e-6
+        (; vehicle_cost, delay_cost, graph, slacks) = instance
+        intrinsic_delays = instance.intrinsic_delays[:, Ω]
+        slacks = deepcopy(instance.slacks)
+        for edge in edges(graph)
+            slacks[src(edge), dst(edge)] = slacks[src(edge), dst(edge)][Ω]
+        end
+        threshold = (c_upp - c_low - vehicle_cost) / delay_cost
+        additional_paths, _ = stochastic_routing_shortest_path_with_threshold(
+            graph, slacks, intrinsic_delays, λ_val ./ delay_cost; threshold
+        )
+        columns = unique(cat(columns, additional_paths; dims=1))
+
+        _, _, sol = compute_solution_from_selected_columns(
+            instance, columns; scenario_range=Ω, silent=silent
+        )
+    end
+
+    col_solution = solution_from_paths(sol, instance)
+    return col_solution
+end
diff --git a/src/StochasticVehicleScheduling/solution/algorithms/local_search.jl b/src/StochasticVehicleScheduling/solution/algorithms/local_search.jl
new file mode 100644
index 0000000..b4f0f0f
--- /dev/null
+++ b/src/StochasticVehicleScheduling/solution/algorithms/local_search.jl
@@ -0,0 +1,157 @@
+"""
+$TYPEDSIGNATURES
+
+Return the optimal solution of the deterministic VSP problem associated to `instance`.
+The objective function is `vehicle_cost * nb_vehicles + include_delays * delay_cost * sum_of_travel_times`
+Note: If you have Gurobi, use `grb_model` as `model_builder` instead od `highs_model`.
+"""
+function solve_deterministic_VSP(
+    instance::Instance; include_delays=true, model_builder=highs_model, verbose=false
+)
+    (; city, graph) = instance
+
+    travel_times = [
+        distance(task1.end_point, task2.start_point) for task1 in city.tasks,
+        task2 in city.tasks
+    ]
+
+    model = model_builder()
+    verbose || set_silent(model)
+
+    nb_nodes = nv(graph)
+    job_indices = 2:(nb_nodes - 1)
+
+    @variable(model, x[i=1:nb_nodes, j=1:nb_nodes; has_edge(graph, i, j)], Bin)
+
+    @objective(
+        model,
+        Min,
+        instance.city.vehicle_cost * sum(x[1, j] for j in job_indices) +
+            include_delays *
+        instance.city.delay_cost *
+        sum(
+            travel_times[i, j] * x[i, j] for i in 1:nb_nodes for
+            j in 1:nb_nodes if has_edge(graph, i, j)
+        )
+    )
+
+    @constraint(
+        model,
+        flow[i in job_indices],
+        sum(x[j, i] for j in inneighbors(graph, i)) ==
+            sum(x[i, j] for j in outneighbors(graph, i))
+    )
+    @constraint(
+        model, demand[i in job_indices], sum(x[j, i] for j in inneighbors(graph, i)) == 1
+    )
+
+    optimize!(model)
+
+    solution = solution_from_JuMP_array(value.(x), graph)
+
+    return JuMP.objective_value(model), solution
+end
+
+"""
+$TYPEDSIGNATURES
+
+Select one random (uniform) task and move it to another random (uniform) feasible vehicle
+"""
+function move_one_random_task!(path_value::BitMatrix, graph::AbstractGraph)
+    nb_tasks = size(path_value, 2)
+    selected_task = rand(DiscreteUniform(1, nb_tasks))
+    selected_vehicle = find_first_one(@view path_value[:, selected_task])
+
+    can_be_inserted = Int[]
+    # do not empty if already empty
+    empty_encountered = false  #sum(@view path_value[selected_vehicle, :]) == 1 ? true : false
+    for i in 1:nb_tasks
+        if i == selected_vehicle
+            continue
+        end
+        # else
+        is_empty = false
+        if selected_task > 1
+            before = @view path_value[i, 1:(selected_task - 1)]
+            if any(before)
+                aaa = find_first_one(reverse(before))
+                @assert aaa >= 0
+                precedent_task = selected_task - aaa
+                if !has_edge(graph, precedent_task + 1, selected_task + 1)
+                    continue
+                end
+            elseif empty_encountered
+                continue
+            else # if !empty_encountered
+                is_empty = true
+            end
+        end
+
+        if selected_task < nb_tasks
+            after = @view path_value[i, (selected_task + 1):end]
+            if any(after)
+                bbb = find_first_one(@view path_value[i, (selected_task + 1):end])
+                @assert bbb >= 0
+                next_task = selected_task + bbb
+
+                if !has_edge(graph, selected_task + 1, next_task + 1)
+                    continue
+                end
+            elseif empty_encountered
+                continue
+            elseif !empty_encountered && is_empty
+                empty_encountered = true
+            end
+        end
+
+        push!(can_be_inserted, i)
+    end
+    if length(can_be_inserted) == 0
+        @warn "No space to be inserted" selected_task path_value
+        return nothing
+    end
+    new_vehicle = rand(can_be_inserted)
+    path_value[selected_vehicle, selected_task] = false
+    path_value[new_vehicle, selected_task] = true
+    return nothing
+end
+
+"""
+$TYPEDSIGNATURES
+
+Very simple local search heuristic, using the neighborhood defined by `move_one_random_task`
+"""
+function _local_search(solution::Solution, instance::Instance; nb_it::Integer=100)
+    best_solution = copy(solution.path_value)
+    best_value = evaluate_solution(solution, instance)
+    history_x = [0]
+    history_y = [best_value]
+
+    candidate_solution = copy(solution.path_value)
+    for it in 1:nb_it
+        move_one_random_task!(candidate_solution, instance.graph)
+
+        value = evaluate_solution(candidate_solution, instance)
+        if value <= best_value # keep changes
+            best_solution = copy(candidate_solution)
+            best_value = value
+            push!(history_x, it)
+            push!(history_y, best_value)
+        else # revert changes
+            candidate_solution = copy(best_solution)
+        end
+    end
+    return Solution(best_solution, instance), best_value, history_x, history_y
+end
+
+"""
+$TYPEDSIGNATURES
+
+Very simple heuristic, using [`local_search`](@ref)
+    initialised with the solution of the deterministic Linear program
+"""
+function local_search(instance::Instance; num_iterations=1000)
+    _, initial_solution = solve_deterministic_VSP(instance)
+    sol, _, _, _ = _local_search(initial_solution, instance; nb_it=num_iterations)
+    return sol
+end
diff --git a/src/StochasticVehicleScheduling/solution/algorithms/mip.jl b/src/StochasticVehicleScheduling/solution/algorithms/mip.jl
new file mode 100644
index 0000000..e202569
--- /dev/null
+++ b/src/StochasticVehicleScheduling/solution/algorithms/mip.jl
@@ -0,0 +1,153 @@
+"""
+$TYPEDSIGNATURES
+
+Returns the optimal solution of the Stochastic VSP instance, by solving the associated compact MIP.
+Quadratic constraints are linearized using Mc Cormick linearization.
+Note: If you have Gurobi, use `grb_model` as `model_builder` instead of `highs_model`.
+"""
+function compact_linearized_mip(
+    instance::Instance; scenario_range=nothing, model_builder=scip_model, silent=true
+)
+    (; graph, slacks, intrinsic_delays, vehicle_cost, delay_cost) = instance
+    nb_nodes = nv(graph)
+    job_indices = 2:(nb_nodes - 1)
+    nodes = 1:nb_nodes
+
+    # Pre-processing
+    ε = intrinsic_delays
+    Rmax = maximum(sum(ε; dims=1))
+    nb_scenarios = size(ε, 2)
+    Ω = isnothing(scenario_range) ? (1:nb_scenarios) : scenario_range
+
+    # Model definition
+    model = model_builder()
+    silent && set_silent(model)
+
+    # Variables and objective function
+    @variable(model, y[u in nodes, v in nodes; has_edge(graph, u, v)], Bin)
+    @variable(model, R[v in nodes, ω in Ω] >= 0) # propagated delay of job v
+    @variable(model, yR[u in nodes, v in nodes, ω in Ω; has_edge(graph, u, v)] >= 0) # yR[u, v] = y[u, v] * R[u, ω]
+
+    @objective(
+        model,
+        Min,
+        delay_cost * sum(sum(R[v, ω] for v in job_indices) for ω in Ω) / nb_scenarios # average total delay
+            +
+            vehicle_cost * sum(y[1, v] for v in job_indices) # nb_vehicles
+    )
+
+    # Flow contraints
+    @constraint(
+        model,
+        flow[i in job_indices],
+        sum(y[j, i] for j in inneighbors(graph, i)) ==
+            sum(y[i, j] for j in outneighbors(graph, i))
+    )
+    @constraint(
+        model,
+        unit_demand[i in job_indices],
+        sum(y[j, i] for j in inneighbors(graph, i)) == 1
+    )
+
+    # Delay propagation constraints
+    @constraint(model, [ω in Ω], R[1, ω] == ε[1, ω])
+    @constraint(model, R_delay_1[v in job_indices, ω in Ω], R[v, ω] >= ε[v, ω])
+    @constraint(
+        model,
+        R_delay_2[v in job_indices, ω in Ω],
+        R[v, ω] >=
+            ε[v, ω] + sum(
+            yR[u, v, ω] - y[u, v] * slacks[u, v][ω] for u in nodes if has_edge(graph, u, v)
+        )
+    )
+
+    # Mc Cormick linearization constraints
+    @constraint(
+        model,
+        R_McCormick_1[u in nodes, v in nodes, ω in Ω; has_edge(graph, u, v)],
+        yR[u, v, ω] >= R[u, ω] + Rmax * (y[u, v] - 1)
+    )
+    @constraint(
+        model,
+        R_McCormick_2[u in nodes, v in nodes, ω in Ω; has_edge(graph, u, v)],
+        yR[u, v, ω] <= Rmax * y[u, v]
+    )
+
+    # Solve model
+    optimize!(model)
+    solution = value.(y)
+
+    sol = solution_from_JuMP_array(solution, graph)
+    return sol
+end
+
+"""
+$TYPEDSIGNATURES
+
+Returns the optimal solution of the Stochastic VSP instance, by solving the associated compact quadratic MIP.
+Note: If you have Gurobi, use `grb_model` as `model_builder` instead of `highs_model`.
+
+!!! warning
+    You need to use a solver that supports quadratic constraints to use this method.
+"""
+function compact_mip(
+    instance::Instance; scenario_range=nothing, model_builder=scip_model, silent=true
+)
+    (; graph, slacks, intrinsic_delays, vehicle_cost, delay_cost) = instance
+    nb_nodes = nv(graph)
+    job_indices = 2:(nb_nodes - 1)
+    nodes = 1:nb_nodes
+
+    # Pre-processing
+    ε = intrinsic_delays
+    nb_scenarios = size(ε, 2)
+    Ω = isnothing(scenario_range) ? (1:nb_scenarios) : scenario_range
+
+    # Model definition
+    model = model_builder()
+    silent && set_silent(model)
+
+    # Variables and objective function
+    @variable(model, y[u in nodes, v in nodes; has_edge(graph, u, v)], Bin)
+    @variable(model, R[v in nodes, ω in Ω] >= 0) # propagated delay of job v
+    @variable(model, yR[u in nodes, v in nodes, ω in Ω; has_edge(graph, u, v)] >= 0) # yR[u, v] = y[u, v] * R[u, ω]
+
+    @objective(
+        model,
+        Min,
+        delay_cost * sum(sum(R[v, ω] for v in job_indices) for ω in Ω) / nb_scenarios # average total delay
+            +
+            vehicle_cost * sum(y[1, v] for v in job_indices) # nb_vehicles
+    )
+
+    # Flow contraints
+    @constraint(
+        model,
+        flow[i in job_indices],
+        sum(y[j, i] for j in inneighbors(graph, i)) ==
+            sum(y[i, j] for j in outneighbors(graph, i))
+    )
+    @constraint(
+        model,
+        unit_demand[i in job_indices],
+        sum(y[j, i] for j in inneighbors(graph, i)) == 1
+    )
+
+    # Delay propagation constraints
+    @constraint(model, [ω in Ω], R[1, ω] == ε[1, ω])
+    @constraint(model, R_delay_1[v in job_indices, ω in Ω], R[v, ω] >= ε[v, ω])
+    @constraint(
+        model,
+        R_delay_2[v in job_indices, ω in Ω],
+        R[v, ω] >=
+            ε[v, ω] +
+        sum(y[u, v] * (R[u, ω] - slacks[u, v][ω]) for u in nodes if has_edge(graph, u, v))
+    )
+
+    # Solve model
+    optimize!(model)
+    solution = value.(y)
+
+    sol = solution_from_JuMP_array(solution, graph)
+    return sol
+end
diff --git a/src/StochasticVehicleScheduling/solution/solution.jl b/src/StochasticVehicleScheduling/solution/solution.jl
new file mode 100644
index 0000000..109d712
--- /dev/null
+++ b/src/StochasticVehicleScheduling/solution/solution.jl
@@ -0,0 +1,386 @@
+# TODO: only keep value field ?
+"""
+$TYPEDEF
+
+Should always be associated with an `Instance`.
+
+# Fields
+$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
+
+# function get_nb_vehicles(path_value::BitMatrix)
+#     return sum(any(path_value; dims=2))
+# end
+
+# function get_nb_vehicles(solution::Solution)
+#     return get_nb_vehicles(solution.path_value)
+# end
+
+# """
+# $TYPEDSIGNATURES
+
+# Compute routes of solution.
+# """
+# function get_routes(solution::Solution)
+#     res = Vector{Int}[]
+#     for vehicle in 1:get_nb_vehicles(solution)
+#         resv = Int[]
+#         for (index, value) in enumerate(solution.path_value[vehicle, :])
+#             if value
+#                 push!(resv, index + 1)
+#             end
+#         end
+#         push!(res, resv)
+#     end
+#     return res
+# end
+
+"""
+$TYPEDSIGNATURES
+
+Create a Solution from a BitVector value.
+"""
+function Solution(value::BitVector, instance::Instance)
+    (; graph) = instance
+    nb_tasks = nv(graph)
+    is_selected = falses(nb_tasks, nb_tasks)
+    for (i, edge) in enumerate(edges(graph))
+        if value[i]
+            is_selected[src(edge), dst(edge)] = true
+        end
+    end
+
+    path_sol = path_solution_from_JuMP_array(is_selected, graph)
+    return Solution(value, path_sol)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Create a Solution from a BitMatrix path value.
+"""
+function Solution(path_value::BitMatrix, instance::Instance)
+    graph = instance.graph
+    solution = falses(ne(graph))
+
+    mat = to_array(path_value, instance)
+
+    for (edge_index, edge) in enumerate(edges(graph))
+        if mat[edge.src, edge.dst]
+            solution[edge_index] = true
+        end
+    end
+
+    return Solution(solution, path_value)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Create a Solution from routes.
+"""
+function solution_from_paths(paths, instance::Instance)
+    (; graph) = instance
+    mat = falses(nv(graph), nv(graph))
+    for p in paths
+        for i in 1:(length(p) - 1)
+            mat[p[i], p[i + 1]] = true
+        end
+    end
+
+    res = falses(ne(graph))
+    for (i, edge) in enumerate(edges(graph))
+        res[i] = mat[src(edge), dst(edge)]
+    end
+
+    return Solution(res, instance)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Create a Solution from a JuMP array.
+"""
+function solution_from_JuMP_array(x::AbstractArray, graph::AbstractGraph)
+    sol = falses(ne(graph)) # init
+
+    for (a, edge) in enumerate(edges(graph))
+        if x[src(edge), dst(edge)] >= 0.5
+            sol[a] = true
+        end
+    end
+
+    return Solution(sol, path_solution_from_JuMP_array(x, graph))
+end
+
+function path_solution_from_JuMP_array(x::AbstractArray, graph::AbstractGraph)
+    nb_tasks = nv(graph)
+    sol = falses(nb_tasks - 2, nb_tasks - 2) # init
+    job_indices = 2:(nb_tasks - 1)
+
+    start = [i for i in job_indices if x[1, i] ≈ 1]
+    for (v_index, task) in enumerate(start)
+        current_task = task
+        while current_task < nb_tasks
+            sol[v_index, current_task - 1] = true
+            next_tasks = [
+                i for i in outneighbors(graph, current_task) if x[current_task, i] >= 0.5
+            ]
+            # TODO : there is a more efficient way to search for next task (but more dangerous)
+            if length(next_tasks) == 1
+                current_task = next_tasks[1]
+            elseif length(next_tasks) == 0
+                @warn "No new task :("
+                @info current_task
+                @info outneighbors(graph, current_task)
+                @info [x[current_task, i] for i in outneighbors(graph, current_task)]
+                current_task = nb_tasks
+            elseif length(next_tasks) > 1
+                @warn "Flow constraint is broken..."
+                current_task = next_tasks[1]
+            end
+        end
+    end
+    return sol
+end
+
+# function basic_path_solution(graph::AbstractGraph)
+#     nb_tasks = nv(graph) - 2
+#     sol = falses(nb_tasks, nb_tasks)
+#     for i_task in 1:nb_tasks
+#         sol[i_task, i_task] = true
+#     end
+#     return sol
+# end
+
+# """
+# $TYPEDSIGNATURES
+
+# Create a solution with one vehicle per task.
+# """
+# function basic_solution(instance::Instance)
+#     graph = instance.graph
+#     value = falses(ne(graph))
+
+#     for (a, edge) in enumerate(edges(graph))
+#         if edge.src == 1 || edge.dst == nv(graph)
+#             value[a] = true
+#         end
+#     end
+
+#     return Solution(value, basic_path_solution(graph))
+# end
+
+"""
+$TYPEDSIGNATURES
+
+Evaluate the total delay of task `i_task` in `scenario`, knowing that current delay from task
+`old_task_index` is `old_delay`.
+"""
+function evaluate_task(
+    i_task::Integer,
+    instance::Instance,
+    old_task_index::Integer,
+    old_delay::Real,
+    scenario::Int,
+)
+    (; slacks, intrinsic_delays) = instance
+
+    delay = intrinsic_delays[i_task, scenario]
+    slack = slacks[old_task_index, i_task][scenario]
+
+    return delay + max(old_delay - slack, 0)
+end
+
+"""
+    evaluate_scenario(path_value::BitMatrix, instance::Instance, scenario_index::Int)
+
+Compute total delay of scenario.
+"""
+function evaluate_scenario(path_value::BitMatrix, instance::Instance, scenario_index::Int)
+    total_delay = 0.0
+    nb_tasks = get_nb_tasks(instance)
+
+    for i_vehicle in 1:nb_tasks
+        # no delay if no tasks
+        if !any(@view path_value[i_vehicle, :])
+            continue
+        end
+
+        task_delay = 0.0
+        old_task_index = 1 # always start at depot
+
+        for i_task in 1:nb_tasks
+            # check if task is done by this vehicle
+            if !path_value[i_vehicle, i_task]
+                continue
+            end
+            task_delay = evaluate_task(
+                i_task + 1, instance, old_task_index, task_delay, scenario_index
+            )
+            old_task_index = i_task + 1
+
+            total_delay += task_delay
+        end
+    end
+    return total_delay
+end
+
+"""
+$TYPEDSIGNATURES
+
+Compute total delay of scenario.
+"""
+function evaluate_scenario(solution::Solution, instance::Instance, scenario_index::Int)
+    return evaluate_scenario(solution.path_value, instance, scenario_index)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Compute total weighted objective of solution.
+"""
+function evaluate_solution(path_value::BitMatrix, instance::Instance)
+    nb_scenarios = get_nb_scenarios(instance)
+
+    average_delay = 0.0
+    for s in 1:nb_scenarios
+        average_delay += evaluate_scenario(path_value, instance, s)
+    end
+    average_delay /= nb_scenarios
+
+    nb_vehicles = sum(any(path_value; dims=2))
+    return instance.vehicle_cost * nb_vehicles + instance.delay_cost * average_delay
+end
+
+"""
+$TYPEDSIGNATURES
+
+Compute total weighted objective of solution.
+"""
+function evaluate_solution(solution::Solution, instance::Instance)
+    return evaluate_solution(solution.path_value, instance)
+end
+
+function evaluate_solution(value::BitVector, instance::Instance)
+    return evaluate_solution(Solution(value, instance), instance)
+end
+
+function to_array(path_value::BitMatrix, instance::Instance)
+    graph = instance.graph
+    nb_nodes = nv(graph)
+    nb_tasks = nb_nodes - 2
+    # job_indices = 2:(nb_nodes - 1)
+    mat = falses(nb_nodes, nb_nodes)
+
+    # check each task used once and only once
+    for i in 1:nb_tasks
+        if !any(@view path_value[i, :])
+            continue
+        end
+        # else
+        current_task = 1
+        while true
+            index_shift = find_first_one(@view path_value[i, current_task:end])
+            if index_shift === -1
+                mat[current_task, nb_nodes] = true
+                break
+            end
+            next_task = current_task + index_shift
+            if !has_edge(graph, current_task, next_task)
+                @warn "Flow not respected" current_task next_task
+                @warn "" outneighbors(graph, current_task)
+                return false
+            end
+            mat[current_task, next_task] = true
+            current_task = next_task
+        end
+    end
+
+    return mat
+end
+
+"""
+$TYPEDSIGNATURES
+
+Returns a BitMatrix, with value true at each index (i, j) if corresponding edge of graph
+is selected in the solution
+"""
+function to_array(solution::Solution, instance::Instance)
+    return to_array(solution.path_value, instance)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Check if `solution` is an admissible solution of `instance`.
+"""
+function is_feasible(solution::Solution, instance::Instance; verbose=true)
+    graph = instance.graph
+    nb_nodes = nv(graph)
+    nb_tasks = nb_nodes - 2
+    job_indices = 2:(nb_nodes - 1)
+    mat = falses(nb_nodes, nb_nodes)
+
+    # check each task used once and only once
+    for i in 1:nb_tasks
+        if !any(@view solution.path_value[i, :])
+            continue
+        end
+        # else
+        current_task = 1
+        while true
+            index_shift = find_first_one(@view solution.path_value[i, current_task:end])
+            if index_shift == -1
+                mat[current_task, nb_nodes] = true
+                break
+            end
+            next_task = current_task + index_shift
+            if !has_edge(graph, current_task, next_task)
+                verbose && @warn "Flow not respected" current_task next_task
+                return false
+            end
+            mat[current_task, next_task] = true
+            current_task = next_task
+        end
+    end
+
+    if !all(sum(solution.path_value; dims=1) .== 1)
+        verbose && @warn "One task done by more than one vehicle (or less than once)"
+        return false
+    end
+
+    for i in job_indices
+        s1 = sum(mat[j, i] for j in inneighbors(graph, i))
+        s2 = sum(mat[i, j] for j in outneighbors(graph, i))
+        if s1 != s2 || s1 != 1
+            verbose && @warn "Flow is broken" i s1 s2
+            return false
+        end
+    end
+
+    return true
+end
+
+function compute_path_list(solution::Solution)
+    (; path_value) = solution
+    paths = Vector{Int64}[]
+    for v in axes(path_value, 1)
+        path = [1]
+        for (i, elem) in enumerate(path_value[v, :])
+            if elem == 1
+                push!(path, i + 1)
+            end
+        end
+        push!(path, size(path_value, 2) + 2)
+        push!(paths, path)
+    end
+    return paths
+end
diff --git a/src/StochasticVehicleScheduling/utils.jl b/src/StochasticVehicleScheduling/utils.jl
new file mode 100644
index 0000000..bca1b3c
--- /dev/null
+++ b/src/StochasticVehicleScheduling/utils.jl
@@ -0,0 +1,99 @@
+"""
+$TYPEDEF
+
+2D point data structure.
+"""
+struct Point
+    x::Float64
+    y::Float64
+end
+
+"""
+$TYPEDSIGNATURES
+
+Returns a Point with random x and y, drawn from distrib.
+"""
+function draw_random_point(distrib::Distribution; rng)
+    return Point(rand(rng, distrib), rand(rng, distrib))
+end
+
+"""
+$TYPEDSIGNATURES
+
+Returns euclidean distance between p₁ and p₂.
+"""
+function distance(p₁::Point, p₂::Point)
+    return sqrt((p₁.x - p₂.x) * (p₁.x - p₂.x) + (p₁.y - p₂.y) * (p₁.y - p₂.y))
+end
+
+"""
+$TYPEDSIGNATURES
+
+Returns hour of the day corresponding to minutes amount.
+"""
+function hour_of(minutes::Real)::Int
+    return min(24, trunc(Int, minutes / 60) + 1)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Returns index of first non zero element of A.
+"""
+function find_first_one(A::AbstractVector)
+    for (i, elem) in enumerate(A)
+        if elem
+            return i
+        end
+    end
+    return -1
+end
+
+# Config stuff, probably not needed in this package
+# """
+#     recursive_namedtuple(x)
+
+# Convert recursively a Dict to a NamedTuple.
+# """
+# recursive_namedtuple(x::Any) = x
+# function recursive_namedtuple(d::Dict)
+#     return namedtuple(Dict(k => recursive_namedtuple(v) for (k, v) in d))
+# end
+
+# """
+#     recursive_convert(x)
+
+# Convert recursively a NamedTuple to a Dict.
+# """
+# recursive_convert(x::Any) = x
+# function recursive_convert(x::NamedTuple)
+#     nt = NamedTuple((k, recursive_convert(v)) for (k, v) in zip(keys(x), x))
+#     return convert(Dict, nt)
+# end
+
+# """
+#     read_config(config_file::String)
+
+# Read a Yaml config into a NamedTuple.
+# """
+# function read_config(config_file::String)
+#     return recursive_namedtuple(YAML.load_file(config_file; dicttype=Dict{Symbol,Any}))
+# end
+
+# """
+#     save_config(config::NamedTuple, save_path::String)
+
+# Save a NamedTuple config to yaml file.
+# """
+# function save_config(config::NamedTuple, save_path::String)
+#     return YAML.write_file(save_path, recursive_convert(config))
+# end
+
+# """
+#     save_config(config::Dict, save_path::String)
+
+# Save Dict config to yaml file.
+# """
+# function save_config(config::Dict, save_path::String)
+#     return YAML.write_file(save_path, config)
+# end
diff --git a/src/SubsetSelection/SubsetSelection.jl b/src/SubsetSelection/SubsetSelection.jl
index 7d05135..0e738a5 100644
--- a/src/SubsetSelection/SubsetSelection.jl
+++ b/src/SubsetSelection/SubsetSelection.jl
@@ -91,5 +91,6 @@ function Utils.generate_statistical_model(bench::SubsetSelectionBenchmark; seed=
 end
 
 export SubsetSelectionBenchmark
+export generate_dataset, generate_maximizer, generate_statistical_model
 
 end
diff --git a/src/Utils/Utils.jl b/src/Utils/Utils.jl
index 65f17d8..60b5b92 100644
--- a/src/Utils/Utils.jl
+++ b/src/Utils/Utils.jl
@@ -2,20 +2,29 @@ module Utils
 
 using DocStringExtensions: TYPEDEF, TYPEDFIELDS, TYPEDSIGNATURES
 using Flux: softplus
+using HiGHS: HiGHS
+using JuMP: Model
 using LinearAlgebra: dot
+using SCIP: SCIP
 using SimpleWeightedGraphs: SimpleWeightedDiGraph
+using StatsBase: StatsBase
+using Statistics: mean
 
 include("data_sample.jl")
 include("interface.jl")
 include("grid_graph.jl")
 include("misc.jl")
+include("model_builders.jl")
 
 export DataSample
 
 export AbstractBenchmark
 export generate_dataset,
     generate_statistical_model, generate_maximizer, plot_data, compute_gap
+export maximizer_kwargs
 export grid_graph, get_path, path_to_matrix
 export neg_tensor, squeeze_last_dims, average_tensor
+export scip_model, highs_model
+export objective_value
 
 end
diff --git a/src/Utils/data_sample.jl b/src/Utils/data_sample.jl
index 33b7bb4..e9a8a3c 100644
--- a/src/Utils/data_sample.jl
+++ b/src/Utils/data_sample.jl
@@ -6,7 +6,9 @@ Data sample data structure.
 # Fields
 $TYPEDFIELDS
 """
-@kwdef struct DataSample{F,S,C,I}
+@kwdef struct DataSample{
+    I,F<:AbstractArray,S<:Union{AbstractArray,Nothing},C<:Union{AbstractArray,Nothing}
+}
     "features"
     x::F
     "target cost parameters (optional)"
@@ -16,3 +18,59 @@ $TYPEDFIELDS
     "instance object (optional)"
     instance::I = nothing
 end
+
+"""
+$TYPEDSIGNATURES
+
+Fit the given transform type (`ZScoreTransform` or `UnitRangeTransform`) on the dataset.
+"""
+function StatsBase.fit(transform_type, dataset::AbstractVector{<:DataSample}; kwargs...)
+    x = hcat([d.x for d in dataset]...)
+    return StatsBase.fit(transform_type, x; kwargs...)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Transform the features in the dataset.
+"""
+function StatsBase.transform(t, dataset::AbstractVector{<:DataSample})
+    return map(dataset) do d
+        (; instance, x, θ_true, y_true) = d
+        DataSample(; instance, x=StatsBase.transform(t, x), θ_true, y_true)
+    end
+end
+
+"""
+$TYPEDSIGNATURES
+
+Transform the features in the dataset in place.
+"""
+function StatsBase.transform!(t, dataset::AbstractVector{<:DataSample})
+    for d in dataset
+        StatsBase.transform!(t, d.x)
+    end
+end
+
+"""
+$TYPEDSIGNATURES
+
+Reconstruct the features in the dataset.
+"""
+function StatsBase.reconstruct(t, dataset::AbstractVector{<:DataSample})
+    return map(dataset) do d
+        (; instance, x, θ_true, y_true) = d
+        DataSample(; instance, x=StatsBase.reconstruct(t, x), θ_true, y_true)
+    end
+end
+
+"""
+$TYPEDSIGNATURES
+
+Reconstruct the features in the dataset in place.
+"""
+function StatsBase.reconstruct!(t, dataset::AbstractVector{<:DataSample})
+    for d in dataset
+        StatsBase.reconstruct!(t, d.x)
+    end
+end
diff --git a/src/Utils/interface.jl b/src/Utils/interface.jl
index 44b21fa..96a2a39 100644
--- a/src/Utils/interface.jl
+++ b/src/Utils/interface.jl
@@ -3,10 +3,12 @@ $TYPEDEF
 
 Abstract type interface for a benchmark problem.
 
-The following methods exist for benchmarks:
+The following methods are mandatory for benchmarks:
 - [`generate_dataset`](@ref)
 - [`generate_statistical_model`](@ref)
 - [`generate_maximizer`](@ref)
+
+The following methods are optional:
 - [`plot_data`](@ref)
 - [`objective_value`](@ref)
 - [`compute_gap`](@ref)
@@ -14,7 +16,7 @@ The following methods exist for benchmarks:
 abstract type AbstractBenchmark end
 
 """
-    generate_dataset(::AbstractBenchmark, dataset_size::Int) -> Vector{<:DataSample}
+    generate_dataset(::AbstractBenchmark, dataset_size::Int; kwargs...) -> Vector{<:DataSample}
 
 Generate a `Vector` of [`DataSample`](@ref)  of length `dataset_size` for given benchmark.
 Content of the dataset can be visualized using [`plot_data`](@ref), when it applies.
@@ -22,7 +24,7 @@ Content of the dataset can be visualized using [`plot_data`](@ref), when it appl
 function generate_dataset end
 
 """
-    generate_maximizer(::AbstractBenchmark)
+    generate_maximizer(::AbstractBenchmark; kwargs...)
 
 Generates a maximizer function.
 Returns a callable f: (θ; kwargs...) -> y, where θ is a cost array and y is a solution.
@@ -30,21 +32,35 @@ Returns a callable f: (θ; kwargs...) -> y, where θ is a cost array and y is a
 function generate_maximizer end
 
 """
-    generate_statistical_model(::AbstractBenchmark)
+    generate_statistical_model(::AbstractBenchmark; kwargs...)
 
 Initializes and return an untrained statistical model of the CO-ML pipeline.
-It's usually a Flux model, that takes a feature matrix x as iinput, and returns a cost array θ as output.
+It's usually a Flux model, that takes a feature matrix x as input, and returns a cost array θ as output.
 """
 function generate_statistical_model end
 
 """
-    plot_data(::AbstractBenchmark, args...)
+    plot_data(::AbstractBenchmark, ::DataSample; kwargs...)
 
 Plot a data sample from the dataset created by [`generate_dataset`](@ref).
 Check the specific benchmark documentation of `plot_data` for more details on the arguments.
 """
 function plot_data end
 
+"""
+    plot_instance(::AbstractBenchmark, instance; kwargs...)
+
+Plot the instance object of the sample.
+"""
+function plot_instance end
+
+"""
+    plot_solution(::AbstractBenchmark, sample::DataSample, [solution]; kwargs...)
+
+Plot `solution` if given, else plot the target solution in the sample.
+"""
+function plot_solution end
+
 """
     compute_gap(::AbstractBenchmark, dataset::Vector{<:DataSample}, statistical_model, maximizer) -> Float64
 
@@ -55,10 +71,61 @@ function compute_gap end
 """
 $TYPEDSIGNATURES
 
+For simple benchmarks where there is no instance object, maximizer does not need any keyword arguments.
+"""
+function maximizer_kwargs(
+    ::AbstractBenchmark, sample::DataSample{Nothing,F,S,C}
+) where {F,S,C}
+    return NamedTuple()
+end
+
+"""
+$TYPEDSIGNATURES
+
+For benchmarks where there is an instance object, maximizer needs the instance object as a keyword argument.
+"""
+function maximizer_kwargs(::AbstractBenchmark, sample::DataSample)
+    return (; instance=sample.instance)
+end
+
+"""
+$TYPEDSIGNATURES
+
 Default behaviour of `objective_value`.
 """
-function objective_value(::AbstractBenchmark, θ̄::AbstractArray, y::AbstractArray)
-    return dot(θ̄, y)
+function objective_value(::AbstractBenchmark, θ::AbstractArray, y::AbstractArray)
+    return dot(θ, y)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Compute the objective value of given solution `y`.
+"""
+function objective_value(
+    bench::AbstractBenchmark, sample::DataSample{I,F,S,C}, y::AbstractArray
+) where {I,F,S,C<:AbstractArray}
+    return objective_value(bench, sample.θ_true, y)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Compute the objective value of the target in the sample (needs to exist).
+"""
+function objective_value(
+    bench::AbstractBenchmark, sample::DataSample{I,F,S,C}
+) where {I,F,S<:AbstractArray,C}
+    return objective_value(bench, sample, sample.y_true)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Check if the benchmark is a minimization problem.
+"""
+function is_minimization_problem(::AbstractBenchmark)
+    return true
 end
 
 """
@@ -67,19 +134,23 @@ $TYPEDSIGNATURES
 Default behaviour of `compute_gap` for a benchmark problem where `features`, `solutions` and `costs` are all defined.
 """
 function compute_gap(
-    bench::AbstractBenchmark, dataset::Vector{<:DataSample}, statistical_model, maximizer
+    bench::AbstractBenchmark,
+    dataset::AbstractVector{<:DataSample},
+    statistical_model,
+    maximizer,
+    op=mean,
 )
-    res = 0.0
-
-    for sample in dataset
-        x = sample.x
-        θ̄ = sample.θ_true
-        ȳ = sample.y_true
-        θ = statistical_model(x)
-        y = maximizer(θ)
-        target_obj = objective_value(bench, θ̄, ȳ)
-        obj = objective_value(bench, θ̄, y)
-        res += (target_obj - obj) / abs(target_obj)
-    end
-    return res / length(dataset)
+    check = is_minimization_problem(bench)
+
+    return op(
+        map(dataset) do sample
+            target_obj = objective_value(bench, sample)
+            x = sample.x
+            θ = statistical_model(x)
+            y = maximizer(θ; maximizer_kwargs(bench, sample)...)
+            obj = objective_value(bench, sample, y)
+            Δ = check ? obj - target_obj : target_obj - obj
+            return Δ / abs(target_obj)
+        end,
+    )
 end
diff --git a/src/Utils/model_builders.jl b/src/Utils/model_builders.jl
new file mode 100644
index 0000000..95df58b
--- /dev/null
+++ b/src/Utils/model_builders.jl
@@ -0,0 +1,20 @@
+"""
+$TYPEDSIGNATURES
+
+Initialize a HiGHS model (with disabled logging).
+"""
+function highs_model()
+    model = Model(HiGHS.Optimizer)
+    # set_attribute(model, "log_to_console", false)
+    return model
+end
+
+"""
+$TYPEDSIGNATURES
+
+Initialize a SCIP model (with disabled logging).
+"""
+function scip_model()
+    model = Model(SCIP.Optimizer)
+    return model
+end
diff --git a/src/Warcraft/Warcraft.jl b/src/Warcraft/Warcraft.jl
index b577f54..669a828 100644
--- a/src/Warcraft/Warcraft.jl
+++ b/src/Warcraft/Warcraft.jl
@@ -25,6 +25,10 @@ Does not have any field.
 """
 struct WarcraftBenchmark <: AbstractBenchmark end
 
+function Utils.objective_value(::WarcraftBenchmark, θ::AbstractArray, y::AbstractArray)
+    return -dot(θ, y)
+end
+
 """
 $TYPEDSIGNATURES
 
diff --git a/src/gurobi_setup.jl b/src/gurobi_setup.jl
new file mode 100644
index 0000000..a5ff979
--- /dev/null
+++ b/src/gurobi_setup.jl
@@ -0,0 +1,20 @@
+using DocStringExtensions: TYPEDSIGNATURES
+using JuMP: Model
+
+@info "Creating a GRB_ENV const for AircraftRoutingBase..."
+# Gurobi package setup (see https://github.com/jump-dev/Gurobi.jl/issues/424)
+const GRB_ENV = Ref{Gurobi.Env}()
+GRB_ENV[] = Gurobi.Env()
+export GRB_ENV
+
+"""
+$TYPEDSIGNATURES
+
+Create an empty Gurobi model.
+"""
+function grb_model()
+    model = Model(() -> Gurobi.Optimizer(GRB_ENV[]))
+    return model
+end
+
+export grb_model
diff --git a/test/vsp.jl b/test/vsp.jl
new file mode 100644
index 0000000..595b38c
--- /dev/null
+++ b/test/vsp.jl
@@ -0,0 +1,46 @@
+@testitem "Stochastic VSP" begin
+    using DecisionFocusedLearningBenchmarks
+    using DecisionFocusedLearningBenchmarks.StochasticVehicleScheduling
+    using Graphs
+    using Plots
+
+    b = StochasticVehicleSchedulingBenchmark(; nb_tasks=25, nb_scenarios=10)
+
+    N = 5
+    dataset = generate_dataset(b, N; seed=0)
+    mip_dataset = generate_dataset(b, N; seed=0, algorithm=compact_mip)
+    mipl_dataset = generate_dataset(b, N; seed=0, algorithm=compact_linearized_mip)
+    local_search_dataset = generate_dataset(b, N; seed=0, algorithm=local_search)
+    @test length(dataset) == N
+
+    figure_1 = plot_instance(b, dataset[1])
+    @test figure_1 isa Plots.Plot
+    figure_2 = plot_solution(b, dataset[1])
+    @test figure_2 isa Plots.Plot
+
+    maximizer = generate_maximizer(b)
+    model = generate_statistical_model(b)
+
+    gap = compute_gap(b, dataset, model, maximizer)
+    gap_mip = compute_gap(b, mip_dataset, model, maximizer)
+    gap_mipl = compute_gap(b, mipl_dataset, model, maximizer)
+    gap_local_search = compute_gap(b, local_search_dataset, model, maximizer)
+
+    @test gap >= 0 && gap_mip >= 0 && gap_mipl >= 0 && gap_local_search >= 0
+    @test gap_mip ≈ gap_mipl rtol = 1e-2
+    @test gap_mip >= gap_local_search
+    @test gap_mip >= gap
+
+    for sample in dataset
+        x = sample.x
+        instance = sample.instance
+        E = ne(instance.graph)
+        @test size(x) == (20, E)
+        θ = model(x)
+        @test length(θ) == E
+        y = maximizer(θ; instance=instance)
+        @test length(y) == E
+        solution = StochasticVehicleScheduling.Solution(y, instance)
+        @test is_feasible(solution, instance)
+    end
+end