fix cvar ub and add cvar tests

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

@@ -7,16 +7,14 @@ Base.@kwdef mutable struct BendersTrainingIterationsLog <: AbstractProgressLog
 end
 
 function BendersTrainingIterationsLog(policy_training_options::PolicyTrainingOptions)
-    println(" ")
-    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)
-    println("Stopping rule: ", policy_training_options.stopping_rule)
-    # TODO add more prints
+    @info(" ")
+    @info("Benders Training")
+    @info(" ")
+    @info("Training options:")
+    @info("Number of scenarios: ", policy_training_options.num_scenarios)
+    @info("Cut strategy: ", policy_training_options.cut_strategy)
+    @info("Risk measure: ", policy_training_options.risk_measure)
+    @info("Stopping rule: ", policy_training_options.stopping_rule)
 
     progress_table = ProgressTable(
         header = ["Iteration", "Lower bound", "Upper bound", "Gap", "Time [s]"],
@@ -71,7 +69,11 @@ function report_current_bounds(progress::BendersTrainingIterationsLog)
     return nothing
 end
 
-function finish_training!(progress::BendersTrainingIterationsLog)
+function finish_training!(
+    progress::BendersTrainingIterationsLog,
+    convergence_result::ConvergenceResult,    
+)
     finalize(progress.progress_table)
+    @info(results_message(convergence_result))
     return nothing
 end

src/progress_logs/deterministic_equivalent.jl

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

src/simulation_strategies/benders_serial.jl

@@ -14,7 +14,7 @@ function serial_benders_simulate(;
     results = Dict{Tuple{String, Int}, Any}() # (variable_name, scenario) => value
 
     # first stage
-    println("Simulating first stage...")
+    @info("Simulating first stage...")
 
     state_variables_model = state_variables_builder(inputs)
     model = first_stage_builder(state_variables_model, inputs)
@@ -31,7 +31,7 @@ function serial_benders_simulate(;
     end
 
     # second stage
-    println("Simulating second stage...")
+    @info("Simulating second stage...")
 
     state = if simulation_options.state_handling == SimulationStateHandling.StatesRecalculatedInSimulation
         get_state(model)

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,12 +109,14 @@ 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)
         if has_converged(convergence_result)
-            println(results_message(convergence_result))
-            finish_training!(progress)
+            finish_training!(progress, convergence_result)
             JQM.send_termination_message()
             break
         end
@@ -224,10 +224,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
@@ -54,8 +54,7 @@ function serial_benders_train(;
         convergence_result =
             convergence_test(progress, policy_training_options.stopping_rule)
         if has_converged(convergence_result)
-            finish_training!(progress)
-            println(results_message(convergence_result))
+            finish_training!(progress, convergence_result)
             break
         end
     end

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()