Implement anticipative solvers for both stochastic benchmarks

Author: BatyLeo

src/ContextualStochasticArgmax/ContextualStochasticArgmax.jl

@@ -118,6 +118,26 @@ end
 
 include("policies.jl")
 
+"""
+$TYPEDSIGNATURES
+
+Generates the anticipative solver for the benchmark.
+"""
+function Utils.generate_anticipative_solver(::ContextualStochasticArgmaxBenchmark)
+    return AnticipativeSolver()
+end
+
+"""
+$TYPEDSIGNATURES
+
+Generates the parametric anticipative solver for the benchmark.
+"""
+function Utils.generate_parametric_anticipative_solver(
+    ::ContextualStochasticArgmaxBenchmark
+)
+    return AnticipativeSolver()
+end
+
 export ContextualStochasticArgmaxBenchmark
 
 end

src/ContextualStochasticArgmax/policies.jl

@@ -1,5 +1,3 @@
-using Statistics: mean
-
 """
 $TYPEDSIGNATURES
 
@@ -28,3 +26,35 @@ Each policy has signature `(ctx_sample, scenarios) -> Vector{DataSample}`.
 function Utils.generate_baseline_policies(::ContextualStochasticArgmaxBenchmark)
     return (; saa=Policy("SAA", "argmax of mean scenarios", csa_saa_policy))
 end
+
+"""
+$TYPEDEF
+
+A policy that acts with perfect information about the future scenario.
+"""
+struct AnticipativeSolver end
+
+function Base.show(io::IO, ::AnticipativeSolver)
+    return print(io, "Anticipative solver for ContextualStochasticArgmaxBenchmark")
+end
+
+"""
+$TYPEDSIGNATURES
+
+Evaluate the anticipative policy for a given `scenario`. 
+Returns the optimal action `one_hot_argmax(scenario)`.
+"""
+function (::AnticipativeSolver)(scenario; context...)
+    return one_hot_argmax(scenario)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Evaluate the anticipative policy with a parametric prediction `θ` and a `scenario`.
+Returns the optimal action for the combined signal `one_hot_argmax(scenario + θ)`.
+"""
+function (::AnticipativeSolver)(θ, scenario; context...)
+    ξ = scenario + θ
+    return one_hot_argmax(ξ)
+end

src/StochasticVehicleScheduling/StochasticVehicleScheduling.jl

@@ -47,6 +47,7 @@ include("solution/algorithms/mip.jl")
 include("solution/algorithms/column_generation.jl")
 include("solution/algorithms/local_search.jl")
 include("solution/algorithms/deterministic_mip.jl")
+include("solution/algorithms/anticipative_solver.jl")
 
 include("maximizer.jl")
 
@@ -113,13 +114,21 @@ end
 $TYPEDSIGNATURES
 
 Return the anticipative solver: a callable `(scenario::VSPScenario; instance, kwargs...) -> y`
-that solves the 1-scenario stochastic VSP via column generation.
+that solves the 1-scenario stochastic VSP.
 """
 function Utils.generate_anticipative_solver(::StochasticVehicleSchedulingBenchmark)
-    return (scenario::VSPScenario; instance::Instance, kwargs...) -> begin
-        stochastic_inst = build_stochastic_instance(instance, [scenario])
-        return column_generation_algorithm(stochastic_inst)
-    end
+    return AnticipativeSolver()
+end
+
+"""
+$TYPEDSIGNATURES
+
+Return the parametric anticipative solver: a callable `(θ, scenario::VSPScenario; instance, kwargs...) -> y`.
+"""
+function Utils.generate_parametric_anticipative_solver(
+    ::StochasticVehicleSchedulingBenchmark
+)
+    return AnticipativeSolver()
 end
 
 """

src/StochasticVehicleScheduling/solution/algorithms/anticipative_solver.jl

@@ -0,0 +1,17 @@
+@kwdef struct AnticipativeSolver{A}
+    single_scenario_algorithm::A = compact_mip
+end
+
+function Base.show(io::IO, ::AnticipativeSolver)
+    return print(io, "Anticipative solver for StochasticVehicleSchedulingBenchmark")
+end
+
+function (solver::AnticipativeSolver)(scenario; instance::Instance, kwargs...)
+    stochastic_inst = build_stochastic_instance(instance, [scenario])
+    return solver.single_scenario_algorithm(stochastic_inst)
+end
+
+function (solver::AnticipativeSolver)(θ, scenario; instance::Instance, kwargs...)
+    stochastic_inst = build_stochastic_instance(instance, [scenario])
+    return solver.single_scenario_algorithm(stochastic_inst, θ)
+end

src/StochasticVehicleScheduling/solution/algorithms/mip.jl

@@ -6,7 +6,11 @@ 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
+    instance::Instance,
+    θ=nothing;
+    scenario_range=nothing,
+    model_builder=scip_model,
+    silent=true,
 )
     (; graph, slacks, intrinsic_delays, vehicle_cost, delay_cost) = instance
     nb_nodes = nv(graph)
@@ -28,13 +32,17 @@ function compact_linearized_mip(
     @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,
+    obj = (
         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
+        +
+        vehicle_cost * sum(y[1, v] for v in job_indices) # nb_vehicles
     )
+    if !isnothing(θ)
+        @assert length(θ) == ne(graph)
+        obj += sum(θ[a] * y[src(edge), dst(edge)] for (a, edge) in enumerate(edges(graph)))
+    end
+
+    @objective(model, Min, obj)
 
     # Flow contraints
     @constraint(
@@ -103,7 +111,11 @@ Note: If you have Gurobi, use `grb_model` as `model_builder` instead of `highs_m
     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
+    instance::Instance,
+    θ=nothing;
+    scenario_range=nothing,
+    model_builder=scip_model,
+    silent=true,
 )
     (; graph, slacks, intrinsic_delays, vehicle_cost, delay_cost) = instance
     nb_nodes = nv(graph)
@@ -124,13 +136,17 @@ function compact_mip(
     @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,
+    obj = (
         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
+        +
+        vehicle_cost * sum(y[1, v] for v in job_indices) # nb_vehicles
     )
+    if !isnothing(θ)
+        @assert length(θ) == ne(graph)
+        obj += sum(θ[a] * y[src(edge), dst(edge)] for (a, edge) in enumerate(edges(graph)))
+    end
+
+    @objective(model, Min, obj)
 
     # Flow contraints
     @constraint(
@@ -170,6 +186,8 @@ $TYPEDSIGNATURES
 SAA variant: build stochastic instance from `scenarios` then solve via
 [`compact_mip`](@ref).
 """
-function compact_mip(instance::Instance, scenarios::Vector{VSPScenario}; kwargs...)
-    return compact_mip(build_stochastic_instance(instance, scenarios); kwargs...)
+function compact_mip(
+    instance::Instance, scenarios::Vector{VSPScenario}, θ=nothing; kwargs...
+)
+    return compact_mip(build_stochastic_instance(instance, scenarios), θ; kwargs...)
 end