diff --git a/src/cut_strategies/cuts_base.jl b/src/cut_strategies/cuts_base.jl index 8bc9d8a..3eb56c2 100644 --- a/src/cut_strategies/cuts_base.jl +++ b/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 diff --git a/src/cut_strategies/multi_cut.jl b/src/cut_strategies/multi_cut.jl index 2cac30f..e44c784 100644 --- a/src/cut_strategies/multi_cut.jl +++ b/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 diff --git a/src/policy.jl b/src/policy.jl index 9c9e850..5d2a532 100644 --- a/src/policy.jl +++ b/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 diff --git a/src/progress_logs/benders_training_iterations.jl b/src/progress_logs/benders_training_iterations.jl index 1de9bcd..fbb04e8 100644 --- a/src/progress_logs/benders_training_iterations.jl +++ b/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 diff --git a/src/progress_logs/deterministic_equivalent.jl b/src/progress_logs/deterministic_equivalent.jl index 4d065e3..ddc2778 100644 --- a/src/progress_logs/deterministic_equivalent.jl +++ b/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 diff --git a/src/simulation_strategies/benders_serial.jl b/src/simulation_strategies/benders_serial.jl index 1c47647..311da74 100644 --- a/src/simulation_strategies/benders_serial.jl +++ b/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) diff --git a/src/training_strategies/benders_job_queue.jl b/src/training_strategies/benders_job_queue.jl index d7ffb92..22173b8 100644 --- a/src/training_strategies/benders_job_queue.jl +++ b/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, diff --git a/src/training_strategies/benders_serial.jl b/src/training_strategies/benders_serial.jl index 1a7b934..b4b5739 100644 --- a/src/training_strategies/benders_serial.jl +++ b/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 diff --git a/test/script_newsvendor_benders_job_queue.jl b/test/script_newsvendor_benders_job_queue.jl index 0ed5574..03d664b 100644 --- a/test/script_newsvendor_benders_job_queue.jl +++ b/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 diff --git a/test/test_newsvendor_benders.jl b/test/test_newsvendor_benders.jl index e9374ec..fefd004 100644 --- a/test/test_newsvendor_benders.jl +++ b/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()