diff --git a/docs/make.jl b/docs/make.jl
index b4eaaf9..1d37ff9 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -13,7 +13,7 @@ tutorial_files = readdir(tutorial_dir)
 md_tutorial_files = [split(file, ".")[1] * ".md" for file in tutorial_files]
 benchmark_files = [joinpath("benchmarks", e) for e in readdir(benchmarks_dir)]
 
-include_tutorial = false
+include_tutorial = true
 
 if include_tutorial
     for file in tutorial_files
diff --git a/src/DynamicVehicleScheduling/DynamicVehicleScheduling.jl b/src/DynamicVehicleScheduling/DynamicVehicleScheduling.jl
index 6b6c5b2..38bd82a 100644
--- a/src/DynamicVehicleScheduling/DynamicVehicleScheduling.jl
+++ b/src/DynamicVehicleScheduling/DynamicVehicleScheduling.jl
@@ -111,7 +111,10 @@ function Utils.generate_policies(b::DynamicVehicleSchedulingBenchmark)
     return (lazy, greedy)
 end
 
-function Utils.generate_statistical_model(b::DynamicVehicleSchedulingBenchmark)
+function Utils.generate_statistical_model(
+    b::DynamicVehicleSchedulingBenchmark; seed=nothing
+)
+    Random.seed!(seed)
     return Chain(Dense((b.two_dimensional_features ? 2 : 14) => 1), vec)
 end
 
diff --git a/src/DynamicVehicleScheduling/anticipative_solver.jl b/src/DynamicVehicleScheduling/anticipative_solver.jl
index 863517f..b2e2452 100644
--- a/src/DynamicVehicleScheduling/anticipative_solver.jl
+++ b/src/DynamicVehicleScheduling/anticipative_solver.jl
@@ -48,30 +48,44 @@ function anticipative_solver(
     model_builder=highs_model,
     two_dimensional_features=env.instance.two_dimensional_features,
     reset_env=true,
-    nb_epochs=typemax(Int),
+    nb_epochs=nothing,
     seed=get_seed(env),
+    verbose=false,
 )
     if reset_env
         reset!(env; reset_rng=true, seed)
     end
 
+    @assert !is_terminated(env)
+
     start_epoch = current_epoch(env)
-    end_epoch = min(last_epoch(env), start_epoch + nb_epochs - 1)
+    end_epoch = if isnothing(nb_epochs)
+        last_epoch(env)
+    else
+        min(last_epoch(env), start_epoch + nb_epochs - 1)
+    end
     T = start_epoch:end_epoch
+    TT = (start_epoch + 1):end_epoch # horizon without start epoch
+
+    starting_state = deepcopy(env.state)
 
     request_epoch = [0]
-    for t in T
+    request_epoch = vcat(request_epoch, fill(start_epoch, customer_count(starting_state)))
+    for t in TT
         request_epoch = vcat(request_epoch, fill(t, length(scenario.indices[t])))
     end
-    customer_index = vcat(1, scenario.indices[T]...)
-    service_time = vcat(0.0, scenario.service_time[T]...)
-    start_time = vcat(0.0, scenario.start_time[T]...)
+
+    customer_index = vcat(starting_state.location_indices, scenario.indices[TT]...)
+    service_time = vcat(
+        starting_state.state_instance.service_time, scenario.service_time[TT]...
+    )
+    start_time = vcat(starting_state.state_instance.start_time, scenario.start_time[TT]...)
 
     duration = env.instance.static_instance.duration[customer_index, customer_index]
     (; epoch_duration, Δ_dispatch) = env.instance
 
     model = model_builder()
-    set_silent(model)
+    verbose || set_silent(model)
 
     nb_nodes = length(customer_index)
     job_indices = 2:nb_nodes
@@ -136,29 +150,25 @@ function anticipative_solver(
         value.(y), env, customer_index, epoch_indices
     )
 
-    epoch_indices = Vector{Int}[]
-    N = 1
-    indices = [1]
     index = 1
-    for epoch in 1:last_epoch(env)
-        M = length(scenario.indices[epoch])
-        indices = vcat(indices, (N + 1):(N + M))
-        push!(epoch_indices, copy(indices))
+    indices = collect(1:(customer_count(starting_state) + 1)) # current known indices in global indexing
+    epoch_indices = [indices] # store global indices present at each epoch
+    N = length(indices) # current last index known in global indexing
+    for epoch in TT # 1:last_epoch(env)
+        # remove dispatched customers from indices
+        dispatched = vcat(epoch_routes[index]...)
+        indices = setdiff(indices, dispatched)
+
+        M = length(scenario.indices[epoch]) # number of new customers in epoch
+        indices = vcat(indices, (N + 1):(N + M))  # add global indices of customers in epoch
+        push!(epoch_indices, copy(indices)) # store global indices present at each epoch
         N = N + M
-        if epoch in T
-            dispatched = vcat(epoch_routes[index]...)
-            index += 1
-            indices = setdiff(indices, dispatched)
-        end
+        index += 1
     end
 
-    indices = vcat(1, scenario.indices...)
-    start_time = vcat(0.0, scenario.start_time...)
-    service_time = vcat(0.0, scenario.service_time...)
-
     dataset = map(enumerate(T)) do (i, epoch)
         routes = epoch_routes[i]
-        epoch_customers = epoch_indices[epoch]
+        epoch_customers = epoch_indices[i]
 
         y_true = VSPSolution(
             Vector{Int}[
@@ -167,7 +177,7 @@ function anticipative_solver(
             max_index=length(epoch_customers),
         ).edge_matrix
 
-        location_indices = indices[epoch_customers]
+        location_indices = customer_index[epoch_customers]
         new_coordinates = env.instance.static_instance.coordinate[location_indices]
         new_start_time = start_time[epoch_customers]
         new_service_time = service_time[epoch_customers]
@@ -184,7 +194,7 @@ function anticipative_solver(
         epoch_duration = env.instance.epoch_duration
         Δ_dispatch = env.instance.Δ_dispatch
         planning_start_time = (epoch - 1) * epoch_duration + Δ_dispatch
-        if epoch == last_epoch
+        if epoch == end_epoch
             # If we are in the last epoch, all requests must be dispatched
             is_must_dispatch[2:end] .= true
         else
@@ -193,6 +203,7 @@ function anticipative_solver(
                 new_start_time[2:end]
         end
         is_postponable[2:end] .= .!is_must_dispatch[2:end]
+        # TODO: avoid code duplication with add_new_customers!
 
         state = DVSPState(;
             state_instance=static_instance,
diff --git a/src/DynamicVehicleScheduling/plot.jl b/src/DynamicVehicleScheduling/plot.jl
index 98defa2..93343b6 100644
--- a/src/DynamicVehicleScheduling/plot.jl
+++ b/src/DynamicVehicleScheduling/plot.jl
@@ -192,49 +192,7 @@ Plot a given DVSPState with routes overlaid. This version accepts routes as a Bi
 where entry (i,j) = true indicates an edge from location i to location j.
 """
 function plot_routes(state::DVSPState, routes::BitMatrix; kwargs...)
-    # Convert BitMatrix to vector of route vectors
-    n_locations = size(routes, 1)
-    route_vectors = Vector{Int}[]
-
-    # Find all outgoing edges from depot (location 1)
-    depot_destinations = findall(routes[1, :])
-
-    # For each destination from depot, reconstruct the route
-    for dest in depot_destinations
-        if dest != 1  # Skip self-loops at depot
-            route = Int[]
-            current = dest
-            push!(route, current)
-
-            # Follow the route until we return to depot
-            while true
-                # Find next location (should be unique for valid routes)
-                next_locations = findall(routes[current, :])
-
-                # Filter out the depot for intermediate steps
-                non_depot_next = filter(x -> x != 1, next_locations)
-
-                if isempty(non_depot_next)
-                    # Must return to depot, route is complete
-                    break
-                elseif length(non_depot_next) == 1
-                    # Continue to next location
-                    current = non_depot_next[1]
-                    push!(route, current)
-                else
-                    # Multiple outgoing edges - this shouldn't happen in valid routes
-                    # but we'll take the first one
-                    current = non_depot_next[1]
-                    push!(route, current)
-                end
-            end
-
-            if !isempty(route)
-                push!(route_vectors, route)
-            end
-        end
-    end
-
+    route_vectors = decode_bitmatrix_to_routes(routes)
     return plot_routes(state, route_vectors; kwargs...)
 end
 
diff --git a/src/DynamicVehicleScheduling/state.jl b/src/DynamicVehicleScheduling/state.jl
index 0d0a177..5b6e37d 100644
--- a/src/DynamicVehicleScheduling/state.jl
+++ b/src/DynamicVehicleScheduling/state.jl
@@ -149,7 +149,8 @@ function is_feasible(state::DVSPState, routes::Vector{Vector{Int}}; verbose::Boo
     if all(is_dispatched[is_must_dispatch])
         return true
     else
-        verbose && @warn "Not all must-dispatch requests are dispatched"
+        verbose &&
+            @warn "Not all must-dispatch requests are dispatched $(is_dispatched[is_must_dispatch])"
         return false
     end
 end
@@ -180,6 +181,58 @@ function apply_routes!(
     return c
 end
 
+function decode_bitmatrix_to_routes(routes::BitMatrix)
+    # Convert BitMatrix to vector of route vectors
+    n_locations = size(routes, 1)
+    route_vectors = Vector{Int}[]
+
+    # Find all outgoing edges from depot (location 1)
+    depot_destinations = findall(routes[1, :])
+
+    # For each destination from depot, reconstruct the route
+    for dest in depot_destinations
+        if dest != 1  # Skip self-loops at depot
+            route = Int[]
+            current = dest
+            push!(route, current)
+
+            # Follow the route until we return to depot
+            while true
+                # Find next location (should be unique for valid routes)
+                next_locations = findall(routes[current, :])
+
+                # Filter out the depot for intermediate steps
+                non_depot_next = filter(x -> x != 1, next_locations)
+
+                if isempty(non_depot_next)
+                    # Must return to depot, route is complete
+                    break
+                elseif length(non_depot_next) == 1
+                    # Continue to next location
+                    current = non_depot_next[1]
+                    push!(route, current)
+                else
+                    throw(
+                        ErrorException(
+                            "Invalid route: multiple outgoing edges from location $current"
+                        ),
+                    )
+                end
+            end
+
+            if !isempty(route)
+                push!(route_vectors, route)
+            end
+        end
+    end
+    return route_vectors
+end
+
+function apply_routes!(state::DVSPState, routes::BitMatrix; check_feasibility::Bool=true)
+    route_vectors = decode_bitmatrix_to_routes(routes)
+    return apply_routes!(state, route_vectors; check_feasibility)
+end
+
 function cost(state::DVSPState, routes::Vector{Vector{Int}})
     return cost(routes, duration(state.state_instance))
 end
diff --git a/src/Utils/interface.jl b/src/Utils/interface.jl
index 1fa1f65..1aabda6 100644
--- a/src/Utils/interface.jl
+++ b/src/Utils/interface.jl
@@ -232,7 +232,7 @@ Generate a vector of environments for the given dynamic benchmark and dataset.
 """
 function generate_environments(
     bench::AbstractDynamicBenchmark,
-    dataset::Vector{<:DataSample};
+    dataset::AbstractArray{<:DataSample};
     seed=nothing,
     rng=MersenneTwister(seed),
     kwargs...,
diff --git a/test/dynamic_vsp.jl b/test/dynamic_vsp.jl
index 0f890c0..c2ea1f3 100644
--- a/test/dynamic_vsp.jl
+++ b/test/dynamic_vsp.jl
@@ -46,4 +46,9 @@
     y2 = maximizer(θ2; instance=instance2)
     @test size(x, 1) == 2
     @test size(x2, 1) == 14
+
+    anticipative_value, solution = generate_anticipative_solution(b, env; reset_env=true)
+    reset!(env; reset_rng=true)
+    cost = sum(step!(env, sample.y_true) for sample in solution)
+    @test isapprox(cost, anticipative_value; atol=1e-5)
 end