fix cvar ub and add cvar tests

Project.toml

@@ -5,6 +5,7 @@ version = "0.2.0"
 
 [deps]
 EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
+HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 JobQueueMPI = "32d208e1-246e-420c-b6ff-18b71b410923"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/cut_strategies/cuts_base.jl

@@ -37,7 +37,7 @@ function get_future_cost(model::JuMP.Model, policy_training_options)::Float64
     if policy_training_options.cut_strategy == CutStrategy.SingleCut
         return LightBenders.get_single_cut_future_cost(model)
     else
-        return LightBenders.get_multi_cut_future_cost(model)
+        return LightBenders.get_multi_cut_future_cost(model, policy_training_options)
     end
 end
 

src/cut_strategies/multi_cut.jl

@@ -115,10 +115,22 @@ function add_all_cuts!(model::JuMP.Model, pool::CutPoolMultiCut, policy_training
     return nothing
 end
 
-function get_multi_cut_future_cost(model::JuMP.Model)::Float64
+function get_multi_cut_future_cost(model::JuMP.Model, policy_training_options)::Float64
     if !haskey(model, :epi_multi_cut)
         return 0.0
     end
-    alphas = model[:epi_multi_cut]::Vector{JuMP.VariableRef}
-    return mean(JuMP.value.(alphas))
+    alphas = JuMP.value.(model[:epi_multi_cut])
+    if policy_training_options.risk_measure isa RiskNeutral
+        return mean(alphas)
+    elseif policy_training_options.risk_measure isa CVaR
+        discount_rate_multiplier = (1.0 - policy_training_options.discount_rate)
+        z_explicit_cvar = JuMP.value(model[:z_explicit_cvar])
+        delta_explicit_cvar = JuMP.value.(model[:delta_explicit_cvar])
+        fcf = z_explicit_cvar * discount_rate_multiplier * (policy_training_options.risk_measure.lambda)
+        for scen in 1:policy_training_options.num_scenarios
+            fcf += alphas[scen] * discount_rate_multiplier * (1 - policy_training_options.risk_measure.lambda) / policy_training_options.num_scenarios
+            fcf += delta_explicit_cvar[scen] * discount_rate_multiplier * (policy_training_options.risk_measure.lambda) / ((1 - policy_training_options.risk_measure.alpha) * policy_training_options.num_scenarios)
+        end
+        return fcf
+    end
 end

src/policy.jl

@@ -12,3 +12,16 @@ end
 function upper_bound(policy::Policy)
     return last_upper_bound(policy.progress)
 end
+
+function second_stage_upper_bound_contribution(policy_training_options::PolicyTrainingOptions, objectives::Vector{Float64})
+    num_scenarios = policy_training_options.num_scenarios
+    expected_value = sum(objectives[s] / num_scenarios for s in 1:num_scenarios)
+    if policy_training_options.risk_measure isa RiskNeutral
+        return expected_value
+    elseif policy_training_options.risk_measure isa CVaR
+        alpha = policy_training_options.risk_measure.alpha
+        lambda = policy_training_options.risk_measure.lambda
+        weights = build_cvar_weights(objectives, alpha, lambda)
+        return dot(weights, objectives)
+    end
+end

src/progress_logs/benders_training_iterations.jl

@@ -11,7 +11,6 @@ function BendersTrainingIterationsLog(policy_training_options::PolicyTrainingOpt
     println("Benders Training")
     println(" ")
     println("Training options:")
-    println("Number of stages: 2")
     println("Number of scenarios: ", policy_training_options.num_scenarios)
     println("Cut strategy: ", policy_training_options.cut_strategy)
     println("Risk measure: ", policy_training_options.risk_measure)

src/progress_logs/deterministic_equivalent.jl

@@ -2,7 +2,6 @@ function DeterministicEquivalentLog(num_scenarios::Int)
     println(" ")
     println("Deterministic Equivalent")
     println(" ")
-    println("Number of stages: 2")
     println("Number of scenarios: ", num_scenarios)
     return nothing
 end

src/training_strategies/benders_job_queue.jl

@@ -16,7 +16,6 @@ mutable struct SecondStageMessage
 end
 
 mutable struct SecondStageAnswer
-    UB::Float64
     coefs
     rhs
     obj
@@ -93,7 +92,6 @@ function job_queue_benders_train(;
                 if !isnothing(job_answer)
                     message = JQM.get_message(job_answer)
                     if message isa SecondStageAnswer
-                        progress.UB[progress.current_iteration] += message.UB
                         store_cut!(
                             local_pools,
                             message.coefs,
@@ -111,6 +109,9 @@ function job_queue_benders_train(;
         # Cuts here can be following the single cut strategy or 
         # the multi cut strategy
         store_cut!(pool, local_pools, state, policy_training_options, t)
+        progress.UB[progress.current_iteration] += second_stage_upper_bound_contribution(
+            policy_training_options, local_pools.obj
+        )
         report_current_bounds(progress)
         convergence_result =
             convergence_test(progress, policy_training_options.stopping_rule)
@@ -224,10 +225,7 @@ function worker_second_stage(
     optimize_with_retry(second_stage_model)
     treat_termination_status(second_stage_model, policy_training_options, t, scenario, iteration)
     coefs, rhs, obj = get_cut(second_stage_model, state)
-    future_cost = get_future_cost(second_stage_model, policy_training_options)
-    UB = (JuMP.objective_value(second_stage_model) - future_cost) / policy_training_options.num_scenarios
     return SecondStageAnswer(
-        UB,
         coefs,
         rhs,
         obj,

src/training_strategies/benders_serial.jl

@@ -40,10 +40,10 @@ function serial_benders_train(;
             coefs, rhs, obj = get_cut(second_stage_model, state)
             # Store the opening cut in a temporary cut pool
             store_cut!(local_pools, coefs, state, rhs, obj)
-            future_cost = get_future_cost(second_stage_model, policy_training_options)
-            progress.UB[progress.current_iteration] +=
-                (JuMP.objective_value(second_stage_model) - future_cost) / policy_training_options.num_scenarios
         end
+        progress.UB[progress.current_iteration] += second_stage_upper_bound_contribution(
+            policy_training_options, local_pools.obj
+        )
         # Store the (stage, scenario) cut(s) in a persitent pool.
         # Cuts here can be following the single cut strategy or 
         # the multi cut strategy

test/script_newsvendor_benders_job_queue.jl

@@ -49,7 +49,10 @@ function second_stage_modifier(sp, inputs, s)
     return nothing
 end
 
-function newsvendor_benders(; cut_strategy = LightBenders.CutStrategy.MultiCut)
+function newsvendor_benders(;
+        cut_strategy = LightBenders.CutStrategy.MultiCut,
+        risk_measure = LightBenders.RiskNeutral(),
+    )
     inputs = Inputs(5, 10, 1, 100, [10, 20, 30])
     num_scenarios = length(inputs.demand)
 
@@ -62,6 +65,7 @@ function newsvendor_benders(; cut_strategy = LightBenders.CutStrategy.MultiCut)
             min_iterations = 2,
         ),
         cut_strategy = cut_strategy,
+        risk_measure = risk_measure,
     )
 
     policy = LightBenders.train(;
@@ -77,9 +81,6 @@ function newsvendor_benders(; cut_strategy = LightBenders.CutStrategy.MultiCut)
         return nothing
     end
 
-    @test LightBenders.lower_bound(policy) ≈ -70
-    @test LightBenders.upper_bound(policy) ≈ -70
-
     results = LightBenders.simulate(;
         state_variables_builder,
         first_stage_builder,
@@ -93,15 +94,55 @@ function newsvendor_benders(; cut_strategy = LightBenders.CutStrategy.MultiCut)
         ),
     )
 
-    @test results["objective", 0] ≈ -70 atol = 1e-2
+    return policy, results
 end
 
 function test_newsvendor_benders()
     @testset "[Job Queue] Benders Newsvendor single cut" begin
-        newsvendor_benders(; cut_strategy = LightBenders.CutStrategy.SingleCut)
+        v = newsvendor_benders(;
+            cut_strategy = LightBenders.CutStrategy.SingleCut
+        )
+        if v !== nothing
+            policy, results = v
+            @test LightBenders.lower_bound(policy) ≈ -70
+            @test LightBenders.upper_bound(policy) ≈ -70
+            @test results["objective", 0] ≈ -70 atol = 1e-2
+        end
     end
     @testset "[Job Queue] Benders Newsvendor multi cut" begin
-        newsvendor_benders(; cut_strategy = LightBenders.CutStrategy.MultiCut)
+        v = newsvendor_benders(;
+            cut_strategy = LightBenders.CutStrategy.MultiCut
+        )
+        if v !== nothing
+            policy, results = v
+            @test LightBenders.lower_bound(policy) ≈ -70
+            @test LightBenders.upper_bound(policy) ≈ -70
+            @test results["objective", 0] ≈ -70 atol = 1e-2
+        end
+    end
+    @testset "[Job Queue] Benders Newsvendor single cut CVaR" begin
+        v = newsvendor_benders(;
+            cut_strategy = LightBenders.CutStrategy.SingleCut,
+            risk_measure = LightBenders.CVaR(alpha = 0.9, lambda = 0.5)
+        )
+        if v !== nothing
+            policy, results = v
+            @test LightBenders.lower_bound(policy) ≈ -50
+            @test LightBenders.upper_bound(policy) ≈ -50
+            @test results["objective", 0] ≈ -50 atol = 1e-2
+        end
+    end
+    @testset "[Job Queue] Benders Newsvendor multi cut CVaR" begin
+        v = newsvendor_benders(;
+            cut_strategy = LightBenders.CutStrategy.MultiCut,
+            risk_measure = LightBenders.CVaR(alpha = 0.9, lambda = 0.5)
+        )
+        if v !== nothing
+            policy, results = v
+            @test LightBenders.lower_bound(policy) ≈ -50
+            @test LightBenders.upper_bound(policy) ≈ -50
+            @test results["objective", 0] ≈ -50 atol = 1e-2
+        end
     end
 end
 

test/test_newsvendor_benders.jl

@@ -49,7 +49,10 @@ function second_stage_modifier(sp, inputs, s)
     return nothing
 end
 
-function newsvendor_benders(; cut_strategy = LightBenders.CutStrategy.MultiCut)
+function newsvendor_benders(;
+    cut_strategy = LightBenders.CutStrategy.MultiCut,
+    risk_measure = LightBenders.RiskNeutral(),
+)
     inputs = Inputs(5, 10, 1, 100, [10, 20, 30])
     num_scenarios = length(inputs.demand)
 
@@ -62,6 +65,7 @@ function newsvendor_benders(; cut_strategy = LightBenders.CutStrategy.MultiCut)
             min_iterations = 2,
         ),
         cut_strategy = cut_strategy,
+        risk_measure = risk_measure,
     )
 
     policy = LightBenders.train(;
@@ -73,9 +77,6 @@ function newsvendor_benders(; cut_strategy = LightBenders.CutStrategy.MultiCut)
         policy_training_options,
     )
 
-    @test LightBenders.lower_bound(policy) ≈ -70
-    @test LightBenders.upper_bound(policy) ≈ -70
-
     results = LightBenders.simulate(;
         state_variables_builder,
         first_stage_builder,
@@ -89,7 +90,7 @@ function newsvendor_benders(; cut_strategy = LightBenders.CutStrategy.MultiCut)
         ),
     )
 
-    @test results["objective", 0] ≈ -70 atol = 1e-2
+    return policy, results
 end
 
 function newsvendor_deterministic()
@@ -111,11 +112,35 @@ function newsvendor_deterministic()
 end
 
 function test_newsvendor_benders()
-    @testset "Benders Newsvendor single cut" begin
-        newsvendor_benders(; cut_strategy = LightBenders.CutStrategy.SingleCut)
+    @testset "Benders Newsvendor single cut risk neutral" begin
+        policy, results = newsvendor_benders(; cut_strategy = LightBenders.CutStrategy.SingleCut)
+        @test LightBenders.lower_bound(policy) ≈ -70
+        @test LightBenders.upper_bound(policy) ≈ -70
+        @test results["objective", 0] ≈ -70 atol = 1e-2
+    end
+    @testset "Benders Newsvendor multi cut risk neutral" begin
+        policy, results = newsvendor_benders(; cut_strategy = LightBenders.CutStrategy.MultiCut)
+        @test LightBenders.lower_bound(policy) ≈ -70
+        @test LightBenders.upper_bound(policy) ≈ -70
+        @test results["objective", 0] ≈ -70 atol = 1e-2
+    end
+    @testset "Benders Newsvendor single cut CVaR" begin
+        policy, results = newsvendor_benders(;
+            cut_strategy = LightBenders.CutStrategy.SingleCut,
+            risk_measure = LightBenders.CVaR(alpha = 0.9, lambda = 0.5),
+        )
+        @test LightBenders.lower_bound(policy) ≈ -50
+        @test LightBenders.upper_bound(policy) ≈ -50
+        @test results["objective", 0] ≈ -50 atol = 1e-2
     end
-    @testset "Benders Newsvendor multi cut" begin
-        newsvendor_benders(; cut_strategy = LightBenders.CutStrategy.MultiCut)
+    @testset "Benders Newsvendor multi cut CVaR" begin
+        policy, results = newsvendor_benders(;
+            cut_strategy = LightBenders.CutStrategy.MultiCut,
+            risk_measure = LightBenders.CVaR(alpha = 0.9, lambda = 0.5),
+        )
+        @test LightBenders.lower_bound(policy) ≈ -50
+        @test LightBenders.upper_bound(policy) ≈ -50
+        @test results["objective", 0] ≈ -50 atol = 1e-2
     end
     @testset "Deterministic equivalent Newsvendor" begin
         newsvendor_deterministic()