Skip to content

Commit

Permalink
Merge pull request #17 from psrenergy/vj/fix-cvar-ub
Browse files Browse the repository at this point in the history
fix cvar ub and add cvar tests
  • Loading branch information
viniciusjusten authored Dec 11, 2024
2 parents 0f53546 + c1a7281 commit 5cfe12c
Show file tree
Hide file tree
Showing 10 changed files with 138 additions and 50 deletions.
2 changes: 1 addition & 1 deletion src/cut_strategies/cuts_base.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
18 changes: 15 additions & 3 deletions src/cut_strategies/multi_cut.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
13 changes: 13 additions & 0 deletions src/policy.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
24 changes: 13 additions & 11 deletions src/progress_logs/benders_training_iterations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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]"],
Expand Down Expand Up @@ -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
9 changes: 4 additions & 5 deletions src/progress_logs/deterministic_equivalent.jl
Original file line number Diff line number Diff line change
@@ -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
4 changes: 2 additions & 2 deletions src/simulation_strategies/benders_serial.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand Down
11 changes: 4 additions & 7 deletions src/training_strategies/benders_job_queue.jl
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,6 @@ mutable struct SecondStageMessage
end

mutable struct SecondStageAnswer
UB::Float64
coefs
rhs
obj
Expand Down Expand Up @@ -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,
Expand All @@ -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
Expand Down Expand Up @@ -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,
Expand Down
9 changes: 4 additions & 5 deletions src/training_strategies/benders_serial.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
55 changes: 48 additions & 7 deletions test/script_newsvendor_benders_job_queue.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand All @@ -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(;
Expand All @@ -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,
Expand All @@ -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

Expand Down
43 changes: 34 additions & 9 deletions test/test_newsvendor_benders.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand All @@ -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(;
Expand All @@ -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,
Expand All @@ -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()
Expand All @@ -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()
Expand Down

0 comments on commit 5cfe12c

Please sign in to comment.