From 7e8a2ee026587756903683909156faf7efd074c9 Mon Sep 17 00:00:00 2001 From: oyurdakul Date: Wed, 22 Feb 2023 12:44:46 -0600 Subject: [PATCH] stochastic extension --- deps/formatter/format.jl | 1 + src/format.jl | 2 + src/instance/read.jl | 101 ++- src/instance/structs.jl | 23 +- src/model/build.jl | 29 +- src/model/formulations/ArrCon2000/ramp.jl | 5 +- src/model/formulations/CarArr2006/pwlcosts.jl | 10 +- .../formulations/DamKucRajAta2016/ramp.jl | 5 +- src/model/formulations/Gar1962/prod.jl | 11 +- src/model/formulations/Gar1962/pwlcosts.jl | 13 +- .../formulations/KnuOstWat2018/pwlcosts.jl | 10 +- src/model/formulations/MorLatRam2013/ramp.jl | 8 +- src/model/formulations/PanGua2016/ramp.jl | 21 +- src/model/formulations/WanHob2016/ramp.jl | 57 +- src/model/formulations/base/bus.jl | 12 +- src/model/formulations/base/line.jl | 20 +- src/model/formulations/base/psload.jl | 9 +- src/model/formulations/base/system.jl | 95 +-- src/model/formulations/base/unit.jl | 154 +++-- src/solution/fix.jl | 46 +- .../methods/XavQiuWanThi2019/enforce.jl | 2 +- src/solution/methods/XavQiuWanThi2019/find.jl | 16 +- .../methods/XavQiuWanThi2019/optimize.jl | 8 +- src/solution/solution.jl | 72 ++- src/transform/initcond.jl | 12 +- src/transform/randomize/XavQiuAhm2021.jl | 50 +- src/transform/slice.jl | 50 +- src/validation/validate.jl | 611 +++++++++--------- test/Project.toml | 1 + test/instance/migrate_test.jl | 14 +- test/instance/read_test.jl | 191 +++--- .../methods/XavQiuWanThi19/filter_test.jl | 143 ++-- .../methods/XavQiuWanThi19/find_test.jl | 51 +- .../XavQiuWanThi19/sensitivity_test.jl | 246 +++---- test/transform/initcond_test.jl | 25 +- .../transform/randomize/XavQiuAhm2021_test.jl | 96 +-- test/transform/slice_test.jl | 46 +- test/usage.jl | 9 +- test/validation/repair_test.jl | 14 +- 39 files changed, 1197 insertions(+), 1092 deletions(-) create mode 100644 src/format.jl diff --git a/deps/formatter/format.jl b/deps/formatter/format.jl index 59ebd54..ecceaf1 100644 --- a/deps/formatter/format.jl +++ b/deps/formatter/format.jl @@ -1,4 +1,5 @@ using JuliaFormatter +print(pwd()) format( [ "../../src", diff --git a/src/format.jl b/src/format.jl new file mode 100644 index 0000000..8b5e7a6 --- /dev/null +++ b/src/format.jl @@ -0,0 +1,2 @@ +using JuliaFormatter +format(["../../src", "../../test", "../../benchmark/run.jl"], verbose = true) diff --git a/src/instance/read.jl b/src/instance/read.jl index 81e7723..935d774 100644 --- a/src/instance/read.jl +++ b/src/instance/read.jl @@ -44,25 +44,6 @@ function read_benchmark( return UnitCommitment.read(filename) end -function read_scenarios(path::AbstractString)::UnitCommitmentInstance - scenario_paths = glob("*.json", path) - scenarios = Vector{UnitCommitmentScenario}() - total_number_of_scenarios = length(scenario_paths) - for (scenario_index, scenario_path) in enumerate(scenario_paths) - scenario = read_scenario(scenario_path, - total_number_of_scenarios = total_number_of_scenarios, - scenario_index = scenario_index) - push!(scenarios, scenario) - end - instance = UnitCommitmentInstance( - time = scenarios[1].time, - scenarios = scenarios - ) - abs(sum(scenario.probability for scenario in instance.scenarios) - 1.0) <= 0.01 || - error("scenario probabilities do not add up to one") - return instance -end - """ read(path::AbstractString)::UnitCommitmentInstance @@ -74,33 +55,54 @@ Read instance from a file. The file may be gzipped. instance = UnitCommitment.read("/path/to/input.json.gz") ``` """ -function read_scenario(path::AbstractString; total_number_of_scenarios = 1, scenario_index = 1)::UnitCommitmentScenario - if endswith(path, ".gz") - return _read(gzopen(path),total_number_of_scenarios, scenario_index) - else - return _read(open(path), total_number_of_scenarios, scenario_index) - end + +function _repair_scenario_name_and_probability( + sc::UnitCommitmentScenario, + path::String, + number_of_scenarios::Int, +)::UnitCommitmentScenario + sc.name !== nothing || (sc.name = first(split(last(split(path, "/")), "."))) + sc.probability !== nothing || (sc.probability = (1 / number_of_scenarios)) + return sc end -function read(path::AbstractString; total_number_of_scenarios = 1, scenario_index = 1)::UnitCommitmentInstance - if endswith(path, ".gz") - scenario = _read(gzopen(path),total_number_of_scenarios, scenario_index) +function read(path)::UnitCommitmentInstance + scenarios = Vector{UnitCommitmentScenario}() + if (endswith(path, ".gz") || endswith(path, ".json")) + endswith(path, ".gz") ? (scenario = _read(gzopen(path))) : + (scenario = _read(open(path))) + scenario = _repair_scenario_name_and_probability(scenario, "s1", 1) + scenarios = [scenario] + elseif typeof(path) == Vector{String} + number_of_scenarios = length(paths) + for scenario_path in path + if endswith(scenario_path, ".gz") + scenario = _read(gzopen(scenario_path)) + elseif endswith(scenario_path, ".json") + scenario = _read(open(scenario_path)) + else + error("Unsupported input format") + end + scenario = _repair_scenario_name_and_probability( + scenario, + scenario_path, + number_of_scenarios, + ) + push!(scenarios, scenario) + end else - scenario = _read(open(path), total_number_of_scenarios, scenario_index) + error("Unsupported input format") end - instance = UnitCommitmentInstance( - time = scenario.time, - scenarios = [scenario] - ) + + instance = + UnitCommitmentInstance(time = scenarios[1].time, scenarios = scenarios) return instance end -function _read(file::IO, total_number_of_scenarios::Int, scenario_index::Int)::UnitCommitmentScenario +function _read(file::IO)::UnitCommitmentScenario return _from_json( JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing)), - total_number_of_scenarios, - scenario_index - ) + ) end function _read_json(path::String)::OrderedDict @@ -112,7 +114,7 @@ function _read_json(path::String)::OrderedDict return JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing)) end -function _from_json(json, total_number_of_scenarios::Int, scenario_index::Int; repair = true)::UnitCommitmentScenario +function _from_json(json; repair = true)::UnitCommitmentScenario _migrate(json) units = Unit[] buses = Bus[] @@ -136,18 +138,8 @@ function _from_json(json, total_number_of_scenarios::Int, scenario_index::Int; r error("Time step $time_step is not a divisor of 60") time_multiplier = 60 ÷ time_step T = time_horizon * time_multiplier - ##### probability = json["Parameters"]["Scenario probability"] - if probability === nothing - probability = (1 / total_number_of_scenarios) - end scenario_name = json["Parameters"]["Scenario name"] - if scenario_name === nothing - scenario_name = "s$(scenario_index)" - end - - ###### - name_to_bus = Dict{String,Bus}() name_to_line = Dict{String,TransmissionLine}() name_to_unit = Dict{String,Unit}() @@ -164,15 +156,6 @@ function _from_json(json, total_number_of_scenarios::Int, scenario_index::Int; r json["Parameters"]["Power balance penalty (\$/MW)"], default = [1000.0 for t in 1:T], ) - # Penalty price for shortage in meeting system-wide flexiramp requirements - flexiramp_shortfall_penalty = timeseries( - json["Parameters"]["Flexiramp penalty (\$/MW)"], - default = [500.0 for t in 1:T], - ) - shortfall_penalty = timeseries( - json["Parameters"]["Reserve shortfall penalty (\$/MW)"], - default = [-1.0 for t in 1:T], - ) # Read buses for (bus_name, dict) in json["Buses"] @@ -368,13 +351,11 @@ function _from_json(json, total_number_of_scenarios::Int, scenario_index::Int; r price_sensitive_loads = loads, reserves = reserves, reserves_by_name = name_to_reserve, - # shortfall_penalty = shortfall_penalty, - # flexiramp_shortfall_penalty = flexiramp_shortfall_penalty, time = T, units_by_name = Dict(g.name => g for g in units), units = units, isf = spzeros(Float64, length(lines), length(buses) - 1), - lodf = spzeros(Float64, length(lines), length(lines)) + lodf = spzeros(Float64, length(lines), length(lines)), ) if repair UnitCommitment.repair!(scenario) diff --git a/src/instance/structs.jl b/src/instance/structs.jl index 56019b1..445a9c3 100644 --- a/src/instance/structs.jl +++ b/src/instance/structs.jl @@ -74,8 +74,8 @@ mutable struct PriceSensitiveLoad end Base.@kwdef mutable struct UnitCommitmentScenario - name::String - probability::Float64 + name::Any + probability::Any buses_by_name::Dict{AbstractString,Bus} buses::Vector{Bus} contingencies_by_name::Dict{AbstractString,Contingency} @@ -87,8 +87,6 @@ Base.@kwdef mutable struct UnitCommitmentScenario price_sensitive_loads::Vector{PriceSensitiveLoad} reserves::Vector{Reserve} reserves_by_name::Dict{AbstractString,Reserve} - # shortfall_penalty::Vector{Float64} - # flexiramp_shortfall_penalty::Vector{Float64} units_by_name::Dict{AbstractString,Unit} units::Vector{Unit} time::Int @@ -104,16 +102,13 @@ end function Base.show(io::IO, instance::UnitCommitmentInstance) print(io, "UnitCommitmentInstance(") print(io, "$(length(instance.scenarios)) scenarios: ") - for scenario in instance.scenarios - print(io, "Scenario $(scenario.name): ") - print(io, "$(length(scenario.units)) units, ") - print(io, "$(length(scenario.buses)) buses, ") - print(io, "$(length(scenario.lines)) lines, ") - print(io, "$(length(scenario.contingencies)) contingencies, ") - print( - io, - "$(length(scenario.price_sensitive_loads)) price sensitive loads, ", - ) + for sc in instance.scenarios + print(io, "Scenario $(sc.name): ") + print(io, "$(length(sc.units)) units, ") + print(io, "$(length(sc.buses)) buses, ") + print(io, "$(length(sc.lines)) lines, ") + print(io, "$(length(sc.contingencies)) contingencies, ") + print(io, "$(length(sc.price_sensitive_loads)) price sensitive loads, ") end print(io, "$(instance.time) time steps") print(io, ")") diff --git a/src/model/build.jl b/src/model/build.jl index 84b10ee..d78e75f 100644 --- a/src/model/build.jl +++ b/src/model/build.jl @@ -78,26 +78,25 @@ function build_model(; model[:obj] = AffExpr() model[:instance] = instance for g in instance.scenarios[1].units - _add_unit_first_stage!(model, g, formulation) + _add_unit_commitment!(model, g, formulation) end - for scenario in instance.scenarios - @info "Building scenario $(scenario.name) with" * - "probability $(scenario.probability)" - _setup_transmission(model, formulation.transmission, scenario) - for l in scenario.lines - _add_transmission_line!(model, l, formulation.transmission, - scenario) + for sc in instance.scenarios + @info "Building scenario $(sc.name) with" * + "probability $(sc.probability)" + _setup_transmission(formulation.transmission, sc) + for l in sc.lines + _add_transmission_line!(model, l, formulation.transmission, sc) end - for b in scenario.buses - _add_bus!(model, b, scenario) + for b in sc.buses + _add_bus!(model, b, sc) end - for ps in scenario.price_sensitive_loads - _add_price_sensitive_load!(model, ps, scenario) + for ps in sc.price_sensitive_loads + _add_price_sensitive_load!(model, ps, sc) end - for g in scenario.units - _add_unit_second_stage!(model, g, formulation, scenario) + for g in sc.units + _add_unit_dispatch!(model, g, formulation, sc) end - _add_system_wide_eqs!(model, scenario) + _add_system_wide_eqs!(model, sc) end @objective(model, Min, model[:obj]) end diff --git a/src/model/formulations/ArrCon2000/ramp.jl b/src/model/formulations/ArrCon2000/ramp.jl index 0af66b6..a043311 100644 --- a/src/model/formulations/ArrCon2000/ramp.jl +++ b/src/model/formulations/ArrCon2000/ramp.jl @@ -8,7 +8,7 @@ function _add_ramp_eqs!( formulation_prod_vars::Gar1962.ProdVars, formulation_ramping::ArrCon2000.Ramping, formulation_status_vars::Gar1962.StatusVars, - sc::UnitCommitmentScenario + sc::UnitCommitmentScenario, )::Nothing # TODO: Move upper case constants to model[:instance] RESERVES_WHEN_START_UP = true @@ -74,7 +74,8 @@ function _add_ramp_eqs!( # but the constraint below will force the unit to produce power eq_ramp_down[sc.name, gn, t] = @constraint( model, - g.initial_power - (g.min_power[t] + prod_above[sc.name, gn, t]) <= RD + g.initial_power - + (g.min_power[t] + prod_above[sc.name, gn, t]) <= RD ) end else diff --git a/src/model/formulations/CarArr2006/pwlcosts.jl b/src/model/formulations/CarArr2006/pwlcosts.jl index 2bd2de6..9b236f7 100644 --- a/src/model/formulations/CarArr2006/pwlcosts.jl +++ b/src/model/formulations/CarArr2006/pwlcosts.jl @@ -8,7 +8,7 @@ function _add_production_piecewise_linear_eqs!( formulation_prod_vars::Gar1962.ProdVars, formulation_pwl_costs::CarArr2006.PwlCosts, formulation_status_vars::StatusVarsFormulation, - sc::UnitCommitmentScenario + sc::UnitCommitmentScenario, )::Nothing eq_prod_above_def = _init(model, :eq_prod_above_def) eq_segprod_limit = _init(model, :eq_segprod_limit) @@ -34,13 +34,17 @@ function _add_production_piecewise_linear_eqs!( # Also add this as an explicit upper bound on segprod to make the # solver's work a bit easier - set_upper_bound(segprod[sc.name, gn, t, k], g.cost_segments[k].mw[t]) + set_upper_bound( + segprod[sc.name, gn, t, k], + g.cost_segments[k].mw[t], + ) # Definition of production # Equation (43) in Kneuven et al. (2020) eq_prod_above_def[sc.name, gn, t] = @constraint( model, - prod_above[sc.name, gn, t] == sum(segprod[sc.name, gn, t, k] for k in 1:K) + prod_above[sc.name, gn, t] == + sum(segprod[sc.name, gn, t, k] for k in 1:K) ) # Objective function diff --git a/src/model/formulations/DamKucRajAta2016/ramp.jl b/src/model/formulations/DamKucRajAta2016/ramp.jl index ad73538..fa380c3 100644 --- a/src/model/formulations/DamKucRajAta2016/ramp.jl +++ b/src/model/formulations/DamKucRajAta2016/ramp.jl @@ -8,7 +8,7 @@ function _add_ramp_eqs!( formulation_prod_vars::Gar1962.ProdVars, formulation_ramping::DamKucRajAta2016.Ramping, formulation_status_vars::Gar1962.StatusVars, - sc::UnitCommitmentScenario + sc::UnitCommitmentScenario, )::Nothing # TODO: Move upper case constants to model[:instance] RESERVES_WHEN_START_UP = true @@ -66,7 +66,8 @@ function _add_ramp_eqs!( elseif (t == 1 && is_initially_on) || (t > 1 && !time_invariant) if t > 1 min_prod_last_period = - prod_above[sc.name, gn, t-1] + g.min_power[t-1] * is_on[gn, t-1] + prod_above[sc.name, gn, t-1] + + g.min_power[t-1] * is_on[gn, t-1] else min_prod_last_period = max(g.initial_power, 0.0) end diff --git a/src/model/formulations/Gar1962/prod.jl b/src/model/formulations/Gar1962/prod.jl index 5571ae1..3e70a04 100644 --- a/src/model/formulations/Gar1962/prod.jl +++ b/src/model/formulations/Gar1962/prod.jl @@ -6,7 +6,7 @@ function _add_production_vars!( model::JuMP.Model, g::Unit, formulation_prod_vars::Gar1962.ProdVars, - sc::UnitCommitmentScenario + sc::UnitCommitmentScenario, )::Nothing prod_above = _init(model, :prod_above) segprod = _init(model, :segprod) @@ -23,7 +23,7 @@ function _add_production_limit_eqs!( model::JuMP.Model, g::Unit, formulation_prod_vars::Gar1962.ProdVars, - sc::UnitCommitmentScenario + sc::UnitCommitmentScenario, )::Nothing eq_prod_limit = _init(model, :eq_prod_limit) is_on = model[:is_on] @@ -34,10 +34,6 @@ function _add_production_limit_eqs!( # Objective function terms for production costs # Part of (69) of Kneuven et al. (2020) as C^R_g * u_g(t) term - ### Moving this term to another function - # add_to_expression!(model[:obj], is_on[gn, t], g.min_power_cost[t]) - ### - # Production limit # Equation (18) in Kneuven et al. (2020) # as \bar{p}_g(t) \le \bar{P}_g u_g(t) @@ -49,7 +45,8 @@ function _add_production_limit_eqs!( end eq_prod_limit[sc.name, gn, t] = @constraint( model, - prod_above[sc.name, gn, t] + reserve[t] <= power_diff * is_on[gn, t] + prod_above[sc.name, gn, t] + reserve[t] <= + power_diff * is_on[gn, t] ) end end diff --git a/src/model/formulations/Gar1962/pwlcosts.jl b/src/model/formulations/Gar1962/pwlcosts.jl index 1566ac7..62d0b5c 100644 --- a/src/model/formulations/Gar1962/pwlcosts.jl +++ b/src/model/formulations/Gar1962/pwlcosts.jl @@ -8,7 +8,7 @@ function _add_production_piecewise_linear_eqs!( formulation_prod_vars::Gar1962.ProdVars, formulation_pwl_costs::Gar1962.PwlCosts, formulation_status_vars::Gar1962.StatusVars, - sc::UnitCommitmentScenario + sc::UnitCommitmentScenario, )::Nothing eq_prod_above_def = _init(model, :eq_prod_above_def) eq_segprod_limit = _init(model, :eq_segprod_limit) @@ -27,7 +27,8 @@ function _add_production_piecewise_linear_eqs!( # Equation (43) in Kneuven et al. (2020) eq_prod_above_def[sc.name, gn, t] = @constraint( model, - prod_above[sc.name, gn, t] == sum(segprod[sc.name, gn, t, k] for k in 1:K) + prod_above[sc.name, gn, t] == + sum(segprod[sc.name, gn, t, k] for k in 1:K) ) for k in 1:K @@ -40,12 +41,16 @@ function _add_production_piecewise_linear_eqs!( # that segment* eq_segprod_limit[sc.name, gn, t, k] = @constraint( model, - segprod[sc.name, gn, t, k] <= g.cost_segments[k].mw[t] * is_on[gn, t] + segprod[sc.name, gn, t, k] <= + g.cost_segments[k].mw[t] * is_on[gn, t] ) # Also add this as an explicit upper bound on segprod to make the # solver's work a bit easier - set_upper_bound(segprod[sc.name, gn, t, k], g.cost_segments[k].mw[t]) + set_upper_bound( + segprod[sc.name, gn, t, k], + g.cost_segments[k].mw[t], + ) # Objective function # Equation (44) in Kneuven et al. (2020) diff --git a/src/model/formulations/KnuOstWat2018/pwlcosts.jl b/src/model/formulations/KnuOstWat2018/pwlcosts.jl index ab95e05..0da8217 100644 --- a/src/model/formulations/KnuOstWat2018/pwlcosts.jl +++ b/src/model/formulations/KnuOstWat2018/pwlcosts.jl @@ -8,7 +8,7 @@ function _add_production_piecewise_linear_eqs!( formulation_prod_vars::Gar1962.ProdVars, formulation_pwl_costs::KnuOstWat2018.PwlCosts, formulation_status_vars::Gar1962.StatusVars, - sc::UnitCommitmentScenario + sc::UnitCommitmentScenario, )::Nothing eq_prod_above_def = _init(model, :eq_prod_above_def) eq_segprod_limit_a = _init(model, :eq_segprod_limit_a) @@ -90,7 +90,8 @@ function _add_production_piecewise_linear_eqs!( # Equation (43) in Kneuven et al. (2020) eq_prod_above_def[sc.name, gn, t] = @constraint( model, - prod_above[sc.name, gn, t] == sum(segprod[sc.name, gn, t, k] for k in 1:K) + prod_above[sc.name, gn, t] == + sum(segprod[sc.name, gn, t, k] for k in 1:K) ) # Objective function @@ -103,7 +104,10 @@ function _add_production_piecewise_linear_eqs!( # Also add an explicit upper bound on segprod to make the solver's # work a bit easier - set_upper_bound(segprod[sc.name, gn, t, k], g.cost_segments[k].mw[t]) + set_upper_bound( + segprod[sc.name, gn, t, k], + g.cost_segments[k].mw[t], + ) end end end diff --git a/src/model/formulations/MorLatRam2013/ramp.jl b/src/model/formulations/MorLatRam2013/ramp.jl index 0672011..90afd71 100644 --- a/src/model/formulations/MorLatRam2013/ramp.jl +++ b/src/model/formulations/MorLatRam2013/ramp.jl @@ -8,7 +8,7 @@ function _add_ramp_eqs!( formulation_prod_vars::Gar1962.ProdVars, formulation_ramping::MorLatRam2013.Ramping, formulation_status_vars::Gar1962.StatusVars, - sc::UnitCommitmentScenario + sc::UnitCommitmentScenario, )::Nothing # TODO: Move upper case constants to model[:instance] RESERVES_WHEN_START_UP = true @@ -65,7 +65,8 @@ function _add_ramp_eqs!( reserve[t] : 0.0 ) min_prod_last_period = - g.min_power[t-1] * is_on[gn, t-1] + prod_above[sc.name, gn, t-1] + g.min_power[t-1] * is_on[gn, t-1] + + prod_above[sc.name, gn, t-1] eq_ramp_up[gn, t] = @constraint( model, max_prod_this_period - min_prod_last_period <= @@ -93,7 +94,8 @@ function _add_ramp_eqs!( # but the constraint below will force the unit to produce power eq_ramp_down[sc.name, gn, t] = @constraint( model, - g.initial_power - (g.min_power[t] + prod_above[sc.name, gn, t]) <= RD + g.initial_power - + (g.min_power[t] + prod_above[sc.name, gn, t]) <= RD ) end else diff --git a/src/model/formulations/PanGua2016/ramp.jl b/src/model/formulations/PanGua2016/ramp.jl index eb85ff5..62a83ba 100644 --- a/src/model/formulations/PanGua2016/ramp.jl +++ b/src/model/formulations/PanGua2016/ramp.jl @@ -8,7 +8,7 @@ function _add_ramp_eqs!( formulation_prod_vars::Gar1962.ProdVars, formulation_ramping::PanGua2016.Ramping, formulation_status_vars::Gar1962.StatusVars, - sc::UnitCommitmentScenario + sc::UnitCommitmentScenario, )::Nothing # TODO: Move upper case constants to model[:instance] RESERVES_WHEN_SHUT_DOWN = true @@ -68,16 +68,17 @@ function _add_ramp_eqs!( if UT - 2 < TRU # Equation (40) in Kneuven et al. (2020) # Covers an additional time period of the ramp-up trajectory, compared to (38) - eq_prod_limit_ramp_up_extra_period[sc.name, gn, t] = @constraint( - model, - prod_above[sc.name, gn, t] + - g.min_power[t] * is_on[gn, t] + - reserve[t] <= - Pbar * is_on[gn, t] - sum( - (Pbar - (SU + i * RU)) * switch_on[gn, t-i] for - i in 0:min(UT - 1, TRU, t - 1) + eq_prod_limit_ramp_up_extra_period[sc.name, gn, t] = + @constraint( + model, + prod_above[sc.name, gn, t] + + g.min_power[t] * is_on[gn, t] + + reserve[t] <= + Pbar * is_on[gn, t] - sum( + (Pbar - (SU + i * RU)) * switch_on[gn, t-i] for + i in 0:min(UT - 1, TRU, t - 1) + ) ) - ) end # Add in shutdown trajectory if KSD >= 0 (else this is dominated by (38)) diff --git a/src/model/formulations/WanHob2016/ramp.jl b/src/model/formulations/WanHob2016/ramp.jl index cc92ae9..417d2a3 100644 --- a/src/model/formulations/WanHob2016/ramp.jl +++ b/src/model/formulations/WanHob2016/ramp.jl @@ -8,7 +8,7 @@ function _add_ramp_eqs!( ::Gar1962.ProdVars, ::WanHob2016.Ramping, ::Gar1962.StatusVars, - sc::UnitCommitmentScenario + sc::UnitCommitmentScenario, )::Nothing is_initially_on = (g.initial_status > 0) SU = g.startup_limit @@ -30,7 +30,7 @@ function _add_ramp_eqs!( error("Each generator may only provide one flexiramp reserve") end for r in g.reserves - if r.type !== "up-frp" && r.type !== "down-frp" + if r.type !== "flexiramp" error( "This formulation only supports flexiramp reserves, not $(r.type)", ) @@ -39,21 +39,23 @@ function _add_ramp_eqs!( for t in 1:model[:instance].time @constraint( model, - prod_above[sc.name, gn, t] + (is_on[gn, t] * minp[t]) <= mfg[sc.name, rn, gn, t] + prod_above[sc.name, gn, t] + (is_on[gn, t] * minp[t]) <= + mfg[sc.name, gn, t] ) # Eq. (19) in Wang & Hobbs (2016) - @constraint(model, mfg[sc.name, rn, gn, t] <= is_on[gn, t] * maxp[t]) # Eq. (22) in Wang & Hobbs (2016) + @constraint(model, mfg[sc.name, gn, t] <= is_on[gn, t] * maxp[t]) # Eq. (22) in Wang & Hobbs (2016) if t != model[:instance].time @constraint( model, minp[t] * (is_on[gn, t+1] + is_on[gn, t] - 1) <= - prod_above[sc.name, gn, t] - dwflexiramp[sc.name, rn, gn, t] + - (is_on[gn, t] * minp[t]) + prod_above[sc.name, gn, t] - + dwflexiramp[sc.name, rn, gn, t] + (is_on[gn, t] * minp[t]) ) # first inequality of Eq. (20) in Wang & Hobbs (2016) @constraint( model, - prod_above[sc.name, gn, t] - dwflexiramp[sc.name, rn, gn, t] + + prod_above[sc.name, gn, t] - + dwflexiramp[sc.name, rn, gn, t] + (is_on[gn, t] * minp[t]) <= - mfg[sc.name, rn, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1])) + mfg[sc.name, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1])) ) # second inequality of Eq. (20) in Wang & Hobbs (2016) @constraint( model, @@ -67,12 +69,12 @@ function _add_ramp_eqs!( prod_above[sc.name, gn, t] + upflexiramp[sc.name, rn, gn, t] + (is_on[gn, t] * minp[t]) <= - mfg[sc.name, rn, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1])) + mfg[sc.name, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1])) ) # second inequality of Eq. (21) in Wang & Hobbs (2016) if t != 1 @constraint( model, - mfg[sc.name, rn, gn, t] <= + mfg[sc.name, gn, t] <= prod_above[sc.name, gn, t-1] + (is_on[gn, t-1] * minp[t]) + (RU * is_on[gn, t-1]) + @@ -81,8 +83,13 @@ function _add_ramp_eqs!( ) # Eq. (23) in Wang & Hobbs (2016) @constraint( model, - (prod_above[sc.name, gn, t-1] + (is_on[gn, t-1] * minp[t])) - - (prod_above[sc.name, gn, t] + (is_on[gn, t] * minp[t])) <= + ( + prod_above[sc.name, gn, t-1] + + (is_on[gn, t-1] * minp[t]) + ) - ( + prod_above[sc.name, gn, t] + + (is_on[gn, t] * minp[t]) + ) <= RD * is_on[gn, t] + SD * (is_on[gn, t-1] - is_on[gn, t]) + maxp[t] * (1 - is_on[gn, t-1]) @@ -90,7 +97,7 @@ function _add_ramp_eqs!( else @constraint( model, - mfg[sc.name, rn, gn, t] <= + mfg[sc.name, gn, t] <= initial_power + (RU * is_initially_on) + (SU * (is_on[gn, t] - is_initially_on)) + @@ -98,8 +105,10 @@ function _add_ramp_eqs!( ) # Eq. (23) in Wang & Hobbs (2016) for the first time period @constraint( model, - initial_power - - (prod_above[sc.name, gn, t] + (is_on[gn, t] * minp[t])) <= + initial_power - ( + prod_above[sc.name, gn, t] + + (is_on[gn, t] * minp[t]) + ) <= RD * is_on[gn, t] + SD * (is_initially_on - is_on[gn, t]) + maxp[t] * (1 - is_initially_on) @@ -107,7 +116,7 @@ function _add_ramp_eqs!( end @constraint( model, - mfg[sc.name, rn, gn, t] <= + mfg[sc.name, gn, t] <= (SD * (is_on[gn, t] - is_on[gn, t+1])) + (maxp[t] * is_on[gn, t+1]) ) # Eq. (24) in Wang & Hobbs (2016) @@ -115,7 +124,8 @@ function _add_ramp_eqs!( model, -RD * is_on[gn, t+1] - SD * (is_on[gn, t] - is_on[gn, t+1]) - - maxp[t] * (1 - is_on[gn, t]) <= upflexiramp[sc.name, rn, gn, t] + maxp[t] * (1 - is_on[gn, t]) <= + upflexiramp[sc.name, rn, gn, t] ) # first inequality of Eq. (26) in Wang & Hobbs (2016) @constraint( model, @@ -127,7 +137,8 @@ function _add_ramp_eqs!( @constraint( model, -RU * is_on[gn, t] - SU * (is_on[gn, t+1] - is_on[gn, t]) - - maxp[t] * (1 - is_on[gn, t+1]) <= dwflexiramp[sc.name, rn, gn, t] + maxp[t] * (1 - is_on[gn, t+1]) <= + dwflexiramp[sc.name, rn, gn, t] ) # first inequality of Eq. (27) in Wang & Hobbs (2016) @constraint( model, @@ -147,7 +158,8 @@ function _add_ramp_eqs!( ) # second inequality of Eq. (28) in Wang & Hobbs (2016) @constraint( model, - -maxp[t] * is_on[gn, t+1] <= dwflexiramp[sc.name, rn, gn, t] + -maxp[t] * is_on[gn, t+1] <= + dwflexiramp[sc.name, rn, gn, t] ) # first inequality of Eq. (29) in Wang & Hobbs (2016) @constraint( model, @@ -157,7 +169,7 @@ function _add_ramp_eqs!( else @constraint( model, - mfg[sc.name, rn, gn, t] <= + mfg[sc.name, gn, t] <= prod_above[sc.name, gn, t-1] + (is_on[gn, t-1] * minp[t]) + (RU * is_on[gn, t-1]) + @@ -166,7 +178,10 @@ function _add_ramp_eqs!( ) # Eq. (23) in Wang & Hobbs (2016) for the last time period @constraint( model, - (prod_above[sc.name, gn, t-1] + (is_on[gn, t-1] * minp[t])) - + ( + prod_above[sc.name, gn, t-1] + + (is_on[gn, t-1] * minp[t]) + ) - (prod_above[sc.name, gn, t] + (is_on[gn, t] * minp[t])) <= RD * is_on[gn, t] + SD * (is_on[gn, t-1] - is_on[gn, t]) + diff --git a/src/model/formulations/base/bus.jl b/src/model/formulations/base/bus.jl index 926bc2c..b28d3aa 100644 --- a/src/model/formulations/base/bus.jl +++ b/src/model/formulations/base/bus.jl @@ -2,7 +2,11 @@ # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. # Released under the modified BSD license. See COPYING.md for more details. -function _add_bus!(model::JuMP.Model, b::Bus, sc::UnitCommitmentScenario)::Nothing +function _add_bus!( + model::JuMP.Model, + b::Bus, + sc::UnitCommitmentScenario, +)::Nothing net_injection = _init(model, :expr_net_injection) curtail = _init(model, :curtail) for t in 1:model[:instance].time @@ -13,7 +17,11 @@ function _add_bus!(model::JuMP.Model, b::Bus, sc::UnitCommitmentScenario)::Nothi curtail[sc.name, b.name, t] = @variable(model, lower_bound = 0, upper_bound = b.load[t]) - add_to_expression!(net_injection[sc.name, b.name, t], curtail[sc.name, b.name, t], 1.0) + add_to_expression!( + net_injection[sc.name, b.name, t], + curtail[sc.name, b.name, t], + 1.0, + ) add_to_expression!( model[:obj], curtail[sc.name, b.name, t], diff --git a/src/model/formulations/base/line.jl b/src/model/formulations/base/line.jl index 0a4a586..98d5f2c 100644 --- a/src/model/formulations/base/line.jl +++ b/src/model/formulations/base/line.jl @@ -6,7 +6,7 @@ function _add_transmission_line!( model::JuMP.Model, lm::TransmissionLine, f::ShiftFactorsFormulation, - sc::UnitCommitmentScenario + sc::UnitCommitmentScenario, )::Nothing overflow = _init(model, :overflow) for t in 1:model[:instance].time @@ -21,30 +21,28 @@ function _add_transmission_line!( end function _setup_transmission( - model::JuMP.Model, formulation::ShiftFactorsFormulation, - scenario::UnitCommitmentScenario + sc::UnitCommitmentScenario, )::Nothing - instance = model[:instance] isf = formulation.precomputed_isf lodf = formulation.precomputed_lodf - if length(scenario.buses) == 1 + if length(sc.buses) == 1 isf = zeros(0, 0) lodf = zeros(0, 0) elseif isf === nothing @info "Computing injection shift factors..." time_isf = @elapsed begin isf = UnitCommitment._injection_shift_factors( - lines = scenario.lines, - buses = scenario.buses, + buses = sc.buses, + lines = sc.lines, ) end @info @sprintf("Computed ISF in %.2f seconds", time_isf) @info "Computing line outage factors..." time_lodf = @elapsed begin lodf = UnitCommitment._line_outage_factors( - lines = scenario.lines, - buses = scenario.buses, + buses = sc.buses, + lines = sc.lines, isf = isf, ) end @@ -57,7 +55,7 @@ function _setup_transmission( isf[abs.(isf).= r.amount[t] + sum( + model[:reserve][sc.name, r.name, g.name, t] for + g in r.units + ) + model[:reserve_shortfall][sc.name, r.name, t] >= + r.amount[t] ) # Account for shortfall contribution to objective @@ -57,7 +71,10 @@ function _add_spinning_reserve_eqs!(model::JuMP.Model, sc::UnitCommitmentScenari return end -function _add_flexiramp_reserve_eqs!(model::JuMP.Model, sc::UnitCommitmentScenario)::Nothing +function _add_flexiramp_reserve_eqs!( + model::JuMP.Model, + sc::UnitCommitmentScenario, +)::Nothing # Note: The flexpramp requirements in Wang & Hobbs (2016) are imposed as hard constraints # through Eq. (17) and Eq. (18). The constraints eq_min_upflexiramp and eq_min_dwflexiramp # provided below are modified versions of Eq. (17) and Eq. (18), respectively, in that @@ -67,39 +84,37 @@ function _add_flexiramp_reserve_eqs!(model::JuMP.Model, sc::UnitCommitmentScenar eq_min_dwflexiramp = _init(model, :eq_min_dwflexiramp) T = model[:instance].time for r in sc.reserves - if r.type == "up-frp" - for t in 1:T - # Eq. (17) in Wang & Hobbs (2016) - eq_min_upflexiramp[sc.name, r.name, t] = @constraint( - model, - sum(model[:upflexiramp][sc.name, r.name, g.name, t] for g in r.units) + - model[:upflexiramp_shortfall][sc.name, r.name, t] >= r.amount[t] - ) - # Account for flexiramp shortfall contribution to objective - if r.shortfall_penalty >= 0 - add_to_expression!( - model[:obj], - r.shortfall_penalty * sc.probability, - model[:upflexiramp_shortfall][sc.name, r.name, t] - ) - end - end - elseif r.type == "down-frp" - for t in 1:T - # Eq. (18) in Wang & Hobbs (2016) - eq_min_dwflexiramp[sc.name, r.name, t] = @constraint( - model, - sum(model[:dwflexiramp][sc.name, r.name, g.name, t] for g in r.units) + - model[:dwflexiramp_shortfall][sc.name, r.name, t] >= r.amount[t] - ) - # Account for flexiramp shortfall contribution to objective - if r.shortfall_penalty >= 0 - add_to_expression!( - model[:obj], - r.shortfall_penalty * sc.probability, + r.type == "flexiramp" || continue + for t in 1:T + # Eq. (17) in Wang & Hobbs (2016) + eq_min_upflexiramp[sc.name, r.name, t] = @constraint( + model, + sum( + model[:upflexiramp][sc.name, r.name, g.name, t] for + g in r.units + ) + model[:upflexiramp_shortfall][sc.name, r.name, t] >= + r.amount[t] + ) + # Eq. (18) in Wang & Hobbs (2016) + eq_min_dwflexiramp[sc.name, r.name, t] = @constraint( + model, + sum( + model[:dwflexiramp][sc.name, r.name, g.name, t] for + g in r.units + ) + model[:dwflexiramp_shortfall][sc.name, r.name, t] >= + r.amount[t] + ) + + # Account for flexiramp shortfall contribution to objective + if r.shortfall_penalty >= 0 + add_to_expression!( + model[:obj], + r.shortfall_penalty * sc.probability, + ( + model[:upflexiramp_shortfall][sc.name, r.name, t] + model[:dwflexiramp_shortfall][sc.name, r.name, t] - ) - end + ), + ) end end end diff --git a/src/model/formulations/base/unit.jl b/src/model/formulations/base/unit.jl index 788ac0a..27cbba6 100644 --- a/src/model/formulations/base/unit.jl +++ b/src/model/formulations/base/unit.jl @@ -2,7 +2,13 @@ # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. # Released under the modified BSD license. See COPYING.md for more details. -function _add_unit_first_stage!(model::JuMP.Model, g::Unit, formulation::Formulation) +# Function for adding variables, constraints, and objective function terms +# related to the binary commitment, startup and shutdown decisions of units +function _add_unit_commitment!( + model::JuMP.Model, + g::Unit, + formulation::Formulation, +) if !all(g.must_run) && any(g.must_run) error("Partially must-run units are not currently supported") end @@ -21,24 +27,30 @@ function _add_unit_first_stage!(model::JuMP.Model, g::Unit, formulation::Formula return end -function _add_unit_second_stage!(model::JuMP.Model, g::Unit, formulation::Formulation, - scenario::UnitCommitmentScenario) +# Function for adding variables, constraints, and objective function terms +# related to the continuous dispatch decisions of units +function _add_unit_dispatch!( + model::JuMP.Model, + g::Unit, + formulation::Formulation, + sc::UnitCommitmentScenario, +) # Variables - _add_production_vars!(model, g, formulation.prod_vars, scenario) - _add_spinning_reserve_vars!(model, g, scenario) - _add_flexiramp_reserve_vars!(model, g, scenario) + _add_production_vars!(model, g, formulation.prod_vars, sc) + _add_spinning_reserve_vars!(model, g, sc) + _add_flexiramp_reserve_vars!(model, g, sc) # Constraints and objective function - _add_net_injection_eqs!(model, g, scenario) - _add_production_limit_eqs!(model, g, formulation.prod_vars, scenario) + _add_net_injection_eqs!(model, g, sc) + _add_production_limit_eqs!(model, g, formulation.prod_vars, sc) _add_production_piecewise_linear_eqs!( model, g, formulation.prod_vars, formulation.pwl_costs, formulation.status_vars, - scenario + sc, ) _add_ramp_eqs!( model, @@ -46,62 +58,29 @@ function _add_unit_second_stage!(model::JuMP.Model, g::Unit, formulation::Formul formulation.prod_vars, formulation.ramping, formulation.status_vars, - scenario + sc, ) - _add_startup_shutdown_limit_eqs!(model, g, scenario) + _add_startup_shutdown_limit_eqs!(model, g, sc) return end -# function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation) -# if !all(g.must_run) && any(g.must_run) -# error("Partially must-run units are not currently supported") -# end -# if g.initial_power === nothing || g.initial_status === nothing -# error("Initial conditions for $(g.name) must be provided") -# end - -# # Variables -# _add_production_vars!(model, g, formulation.prod_vars) -# _add_spinning_reserve_vars!(model, g) -# _add_flexiramp_reserve_vars!(model, g) -# _add_startup_shutdown_vars!(model, g) -# _add_status_vars!(model, g, formulation.status_vars) - -# # Constraints and objective function -# _add_min_uptime_downtime_eqs!(model, g) -# _add_net_injection_eqs!(model, g) -# _add_production_limit_eqs!(model, g, formulation.prod_vars) -# _add_production_piecewise_linear_eqs!( -# model, -# g, -# formulation.prod_vars, -# formulation.pwl_costs, -# formulation.status_vars, -# ) -# _add_ramp_eqs!( -# model, -# g, -# formulation.prod_vars, -# formulation.ramping, -# formulation.status_vars, -# ) -# _add_startup_cost_eqs!(model, g, formulation.startup_costs) -# _add_startup_shutdown_limit_eqs!(model, g) -# _add_status_eqs!(model, g, formulation.status_vars) -# return -# end - _is_initially_on(g::Unit)::Float64 = (g.initial_status > 0 ? 1.0 : 0.0) -function _add_spinning_reserve_vars!(model::JuMP.Model, g::Unit, sc::UnitCommitmentScenario)::Nothing +function _add_spinning_reserve_vars!( + model::JuMP.Model, + g::Unit, + sc::UnitCommitmentScenario, +)::Nothing reserve = _init(model, :reserve) reserve_shortfall = _init(model, :reserve_shortfall) for r in g.reserves r.type == "spinning" || continue for t in 1:model[:instance].time - reserve[sc.name, r.name, g.name, t] = @variable(model, lower_bound = 0) + reserve[sc.name, r.name, g.name, t] = + @variable(model, lower_bound = 0) if (sc.name, r.name, t) ∉ keys(reserve_shortfall) - reserve_shortfall[sc.name, r.name, t] = @variable(model, lower_bound = 0) + reserve_shortfall[sc.name, r.name, t] = + @variable(model, lower_bound = 0) if r.shortfall_penalty < 0 set_upper_bound(reserve_shortfall[sc.name, r.name, t], 0.0) end @@ -111,35 +90,37 @@ function _add_spinning_reserve_vars!(model::JuMP.Model, g::Unit, sc::UnitCommitm return end -function _add_flexiramp_reserve_vars!(model::JuMP.Model, g::Unit, sc::UnitCommitmentScenario)::Nothing +function _add_flexiramp_reserve_vars!( + model::JuMP.Model, + g::Unit, + sc::UnitCommitmentScenario, +)::Nothing upflexiramp = _init(model, :upflexiramp) upflexiramp_shortfall = _init(model, :upflexiramp_shortfall) mfg = _init(model, :mfg) dwflexiramp = _init(model, :dwflexiramp) dwflexiramp_shortfall = _init(model, :dwflexiramp_shortfall) - for r in g.reserves - if r.type == "up-frp" - for t in 1:model[:instance].time - # maximum feasible generation, \bar{g_{its}} in Wang & Hobbs (2016) - mfg[sc.name, r.name, g.name, t] = @variable(model, lower_bound = 0) - upflexiramp[sc.name, r.name, g.name, t] = @variable(model) # up-flexiramp, ur_{it} in Wang & Hobbs (2016) - if (sc.name, r.name, t) ∉ keys(upflexiramp_shortfall) - upflexiramp_shortfall[sc.name, r.name, t] = - @variable(model, lower_bound = 0) - if r.shortfall_penalty < 0 - set_upper_bound(upflexiramp_shortfall[sc.name, r.name, t], 0.0) - end - end - end - elseif r.type == "down-frp" - for t in 1:model[:instance].time - dwflexiramp[sc.name, r.name, g.name, t] = @variable(model) # down-flexiramp, dr_{it} in Wang & Hobbs (2016) - if (sc.name, r.name, t) ∉ keys(dwflexiramp_shortfall) - dwflexiramp_shortfall[sc.name, r.name, t] = - @variable(model, lower_bound = 0) - if r.shortfall_penalty < 0 - set_upper_bound(dwflexiramp_shortfall[sc.name, r.name, t], 0.0) - end + for t in 1:model[:instance].time + # maximum feasible generation, \bar{g_{its}} in Wang & Hobbs (2016) + mfg[sc.name, g.name, t] = @variable(model, lower_bound = 0) + for r in g.reserves + r.type == "flexiramp" || continue + upflexiramp[sc.name, r.name, g.name, t] = @variable(model) # up-flexiramp, ur_{it} in Wang & Hobbs (2016) + dwflexiramp[sc.name, r.name, g.name, t] = @variable(model) # down-flexiramp, dr_{it} in Wang & Hobbs (2016) + if (sc.name, r.name, t) ∉ keys(upflexiramp_shortfall) + upflexiramp_shortfall[sc.name, r.name, t] = + @variable(model, lower_bound = 0) + dwflexiramp_shortfall[sc.name, r.name, t] = + @variable(model, lower_bound = 0) + if r.shortfall_penalty < 0 + set_upper_bound( + upflexiramp_shortfall[sc.name, r.name, t], + 0.0, + ) + set_upper_bound( + dwflexiramp_shortfall[sc.name, r.name, t], + 0.0, + ) end end end @@ -157,7 +138,11 @@ function _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing return end -function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit, sc::UnitCommitmentScenario)::Nothing +function _add_startup_shutdown_limit_eqs!( + model::JuMP.Model, + g::Unit, + sc::UnitCommitmentScenario, +)::Nothing eq_shutdown_limit = _init(model, :eq_shutdown_limit) eq_startup_limit = _init(model, :eq_startup_limit) is_on = model[:is_on] @@ -196,7 +181,7 @@ function _add_ramp_eqs!( model::JuMP.Model, g::Unit, formulation::RampingFormulation, - sc::UnitCommitmentScenario + sc::UnitCommitmentScenario, )::Nothing prod_above = model[:prod_above] reserve = _total_reserves(model, g, sc) @@ -282,7 +267,11 @@ function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing end end -function _add_net_injection_eqs!(model::JuMP.Model, g::Unit, sc::UnitCommitmentScenario)::Nothing +function _add_net_injection_eqs!( + model::JuMP.Model, + g::Unit, + sc::UnitCommitmentScenario, +)::Nothing expr_net_injection = model[:expr_net_injection] for t in 1:model[:instance].time # Add to net injection expression @@ -305,7 +294,10 @@ function _total_reserves(model, g, sc)::Vector spinning_reserves = [r for r in g.reserves if r.type == "spinning"] if !isempty(spinning_reserves) reserve += [ - sum(model[:reserve][sc.name, r.name, g.name, t] for r in spinning_reserves) for t in 1:model[:instance].time + sum( + model[:reserve][sc.name, r.name, g.name, t] for + r in spinning_reserves + ) for t in 1:model[:instance].time ] end return reserve diff --git a/src/solution/fix.jl b/src/solution/fix.jl index b2d37d8..ca7703d 100644 --- a/src/solution/fix.jl +++ b/src/solution/fix.jl @@ -10,37 +10,43 @@ solution. Useful for computing LMPs. """ function fix!(model::JuMP.Model, solution::AbstractDict)::Nothing instance, T = model[:instance], model[:instance].time + "Production (MW)" ∈ keys(solution) ? solution = Dict("s1" => solution) : + nothing is_on = model[:is_on] prod_above = model[:prod_above] reserve = model[:reserve] - for g in instance.units - for t in 1:T - is_on_value = round(solution["Is on"][g.name][t]) - prod_value = - round(solution["Production (MW)"][g.name][t], digits = 5) - JuMP.fix(is_on[g.name, t], is_on_value, force = true) - JuMP.fix( - prod_above[g.name, t], - prod_value - is_on_value * g.min_power[t], - force = true, - ) - end - end - for r in instance.reserves - r.type == "spinning" || continue - for g in r.units + for sc in instance.scenarios + for g in sc.units for t in 1:T - reserve_value = round( - solution["Spinning reserve (MW)"][r.name][g.name][t], + is_on_value = round(solution[sc.name]["Is on"][g.name][t]) + prod_value = round( + solution[sc.name]["Production (MW)"][g.name][t], digits = 5, ) + JuMP.fix(is_on[g.name, t], is_on_value, force = true) JuMP.fix( - reserve[r.name, g.name, t], - reserve_value, + prod_above[sc.name, g.name, t], + prod_value - is_on_value * g.min_power[t], force = true, ) end end + for r in sc.reserves + r.type == "spinning" || continue + for g in r.units + for t in 1:T + reserve_value = round( + solution[sc.name]["Spinning reserve (MW)"][r.name][g.name][t], + digits = 5, + ) + JuMP.fix( + reserve[sc.name, r.name, g.name, t], + reserve_value, + force = true, + ) + end + end + end end return end diff --git a/src/solution/methods/XavQiuWanThi2019/enforce.jl b/src/solution/methods/XavQiuWanThi2019/enforce.jl index a4e6ce7..8c5f3b0 100644 --- a/src/solution/methods/XavQiuWanThi2019/enforce.jl +++ b/src/solution/methods/XavQiuWanThi2019/enforce.jl @@ -5,7 +5,7 @@ function _enforce_transmission( model::JuMP.Model, violations::Vector{_Violation}, - sc::UnitCommitmentScenario + sc::UnitCommitmentScenario, )::Nothing for v in violations _enforce_transmission( diff --git a/src/solution/methods/XavQiuWanThi2019/find.jl b/src/solution/methods/XavQiuWanThi2019/find.jl index 43e3ad3..b9d97e6 100644 --- a/src/solution/methods/XavQiuWanThi2019/find.jl +++ b/src/solution/methods/XavQiuWanThi2019/find.jl @@ -19,8 +19,8 @@ function _find_violations( time_screening = @elapsed begin non_slack_buses = [b for b in sc.buses if b.offset > 0] net_injection_values = [ - value(net_injection[sc.name, b.name, t]) for b in non_slack_buses, - t in 1:instance.time + value(net_injection[sc.name, b.name, t]) for + b in non_slack_buses, t in 1:instance.time ] overflow_values = [ value(overflow[sc.name, lm.name, t]) for lm in sc.lines, @@ -72,7 +72,7 @@ function _find_violations(; isf::Array{Float64,2}, lodf::Array{Float64,2}, max_per_line::Int, - max_per_period::Int + max_per_period::Int, )::Array{_Violation,1} B = length(sc.buses) - 1 L = length(sc.lines) @@ -96,13 +96,13 @@ function _find_violations(; post_v::Array{Float64} = zeros(L, L, K) # post_v[lm, lc, thread] normal_limits::Array{Float64,2} = [ - l.normal_flow_limit[t] + overflow[l.offset, t] for - l in sc.lines, t in 1:T + l.normal_flow_limit[t] + overflow[l.offset, t] for l in sc.lines, + t in 1:T ] emergency_limits::Array{Float64,2} = [ - l.emergency_flow_limit[t] + overflow[l.offset, t] for - l in sc.lines, t in 1:T + l.emergency_flow_limit[t] + overflow[l.offset, t] for l in sc.lines, + t in 1:T ] is_vulnerable::Array{Bool} = zeros(Bool, L) @@ -114,7 +114,7 @@ function _find_violations(; k = threadid() # Pre-contingency flows - pre_flow[:, k] = isf * net_injections[ :, t] + pre_flow[:, k] = isf * net_injections[:, t] # Post-contingency flows for lc in 1:L, lm in 1:L diff --git a/src/solution/methods/XavQiuWanThi2019/optimize.jl b/src/solution/methods/XavQiuWanThi2019/optimize.jl index c7b3a7e..14a34d2 100644 --- a/src/solution/methods/XavQiuWanThi2019/optimize.jl +++ b/src/solution/methods/XavQiuWanThi2019/optimize.jl @@ -10,9 +10,9 @@ function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing JuMP.set_optimizer_attribute(model, "MIPGap", gap) @info @sprintf("MIP gap tolerance set to %f", gap) end - for scenario in model[:instance].scenarios + for sc in model[:instance].scenarios large_gap = false - has_transmission = (length(scenario.isf) > 0) + has_transmission = (length(sc.isf) > 0) if has_transmission && method.two_phase_gap set_gap(1e-2) large_gap = true @@ -36,7 +36,7 @@ function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing has_transmission || break violations = _find_violations( model, - scenario, + sc, max_per_line = method.max_violations_per_line, max_per_period = method.max_violations_per_period, ) @@ -49,7 +49,7 @@ function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing break end else - _enforce_transmission(model, violations, scenario) + _enforce_transmission(model, violations, sc) end end end diff --git a/src/solution/solution.jl b/src/solution/solution.jl index f02fb5b..b170f30 100644 --- a/src/solution/solution.jl +++ b/src/solution/solution.jl @@ -16,17 +16,21 @@ solution = UnitCommitment.solution(model) """ function solution(model::JuMP.Model)::OrderedDict instance, T = model[:instance], model[:instance].time - function timeseries_first_stage(vars, collection) - return OrderedDict( - b.name => [round(value(vars[b.name, t]), digits = 5) for t in 1:T] - for b in collection - ) - end - function timeseries_second_stage(vars, collection, sc) - return OrderedDict( - b.name => [round(value(vars[sc.name, b.name, t]), digits = 5) for t in 1:T] - for b in collection - ) + function timeseries(vars, collection; sc = nothing) + if sc === nothing + return OrderedDict( + b.name => + [round(value(vars[b.name, t]), digits = 5) for t in 1:T] for + b in collection + ) + else + return OrderedDict( + b.name => [ + round(value(vars[sc.name, b.name, t]), digits = 5) for + t in 1:T + ] for b in collection + ) + end end function production_cost(g, sc) return [ @@ -67,24 +71,25 @@ function solution(model::JuMP.Model)::OrderedDict OrderedDict(g.name => production_cost(g, sc) for g in sc.units) sol[sc.name]["Startup cost (\$)"] = OrderedDict(g.name => startup_cost(g, sc) for g in sc.units) - sol[sc.name]["Is on"] = timeseries_first_stage(model[:is_on], sc.units) - sol[sc.name]["Switch on"] = timeseries_first_stage(model[:switch_on], sc.units) - sol[sc.name]["Switch off"] = timeseries_first_stage(model[:switch_off], sc.units) + sol[sc.name]["Is on"] = timeseries(model[:is_on], sc.units) + sol[sc.name]["Switch on"] = timeseries(model[:switch_on], sc.units) + sol[sc.name]["Switch off"] = timeseries(model[:switch_off], sc.units) sol[sc.name]["Net injection (MW)"] = - timeseries_second_stage(model[:net_injection], sc.buses, sc) - sol[sc.name]["Load curtail (MW)"] = timeseries_second_stage(model[:curtail], sc.buses, sc) + timeseries(model[:net_injection], sc.buses, sc = sc) + sol[sc.name]["Load curtail (MW)"] = + timeseries(model[:curtail], sc.buses, sc = sc) if !isempty(sc.lines) - sol[sc.name]["Line overflow (MW)"] = timeseries_second_stage(model[:overflow], sc.lines, sc) + sol[sc.name]["Line overflow (MW)"] = + timeseries(model[:overflow], sc.lines, sc = sc) end if !isempty(sc.price_sensitive_loads) sol[sc.name]["Price-sensitive loads (MW)"] = - timeseries_second_stage(model[:loads], sc.price_sensitive_loads, sc) + timeseries(model[:loads], sc.price_sensitive_loads, sc = sc) end sol[sc.name]["Spinning reserve (MW)"] = OrderedDict( r.name => OrderedDict( g.name => [ - value(model[:reserve][sc.name, r.name, g.name, t]) for - t in 1:instance.time + value(model[:reserve][sc.name, r.name, g.name, t]) for t in 1:instance.time ] for g in r.units ) for r in sc.reserves if r.type == "spinning" ) @@ -97,32 +102,31 @@ function solution(model::JuMP.Model)::OrderedDict sol[sc.name]["Up-flexiramp (MW)"] = OrderedDict( r.name => OrderedDict( g.name => [ - value(model[:upflexiramp][sc.name, r.name, g.name, t]) for - t in 1:instance.time + value(model[:upflexiramp][sc.name, r.name, g.name, t]) for t in 1:instance.time ] for g in r.units - ) for r in sc.reserves if r.type == "up-frp" + ) for r in sc.reserves if r.type == "flexiramp" ) sol[sc.name]["Up-flexiramp shortfall (MW)"] = OrderedDict( r.name => [ - value(model[:upflexiramp_shortfall][sc.name, r.name, t]) for - t in 1:instance.time - ] for r in sc.reserves if r.type == "up-frp" + value(model[:upflexiramp_shortfall][sc.name, r.name, t]) for t in 1:instance.time + ] for r in sc.reserves if r.type == "flexiramp" ) sol[sc.name]["Down-flexiramp (MW)"] = OrderedDict( r.name => OrderedDict( g.name => [ - value(model[:dwflexiramp][sc.name, r.name, g.name, t]) for - t in 1:instance.time + value(model[:dwflexiramp][sc.name, r.name, g.name, t]) for t in 1:instance.time ] for g in r.units - ) for r in sc.reserves if r.type == "down-frp" + ) for r in sc.reserves if r.type == "flexiramp" ) sol[sc.name]["Down-flexiramp shortfall (MW)"] = OrderedDict( r.name => [ - value(model[:dwflexiramp_shortfall][sc.name, r.name, t]) for - t in 1:instance.time - ] for r in sc.reserves if r.type == "down-frp" + value(model[:dwflexiramp_shortfall][sc.name, r.name, t]) for t in 1:instance.time + ] for r in sc.reserves if r.type == "flexiramp" ) end - length(keys(sol)) > 1 ? nothing : sol = Dict(sol_key => sol_val for scen_key in keys(sol) for (sol_key, sol_val) in sol[scen_key]) - return sol + if length(instance.scenarios) == 1 + return first(values(sol)) + else + return sol + end end diff --git a/src/transform/initcond.jl b/src/transform/initcond.jl index 6bbf0ea..cfbce02 100644 --- a/src/transform/initcond.jl +++ b/src/transform/initcond.jl @@ -5,18 +5,18 @@ using JuMP """ - generate_initial_conditions!(instance, optimizer) + generate_initial_conditions!(sc, optimizer) -Generates feasible initial conditions for the given instance, by constructing +Generates feasible initial conditions for the given scenario, by constructing and solving a single-period mixed-integer optimization problem, using the given -optimizer. The instance is modified in-place. +optimizer. The scenario is modified in-place. """ function generate_initial_conditions!( - instance::UnitCommitmentInstance, + sc::UnitCommitmentScenario, optimizer, )::Nothing - G = instance.units - B = instance.buses + G = sc.units + B = sc.buses t = 1 mip = JuMP.Model(optimizer) diff --git a/src/transform/randomize/XavQiuAhm2021.jl b/src/transform/randomize/XavQiuAhm2021.jl index 27bade1..a8ed08e 100644 --- a/src/transform/randomize/XavQiuAhm2021.jl +++ b/src/transform/randomize/XavQiuAhm2021.jl @@ -6,6 +6,7 @@ module XavQiuAhm2021 using Distributions import ..UnitCommitmentInstance +import ..UnitCommitmentScenario """ struct Randomization @@ -119,10 +120,10 @@ end function _randomize_costs( rng, - instance::UnitCommitmentInstance, + sc::UnitCommitmentScenario, distribution, )::Nothing - for unit in instance.units + for unit in sc.units α = rand(rng, distribution) unit.min_power_cost *= α for k in unit.cost_segments @@ -137,17 +138,15 @@ end function _randomize_load_share( rng, - instance::UnitCommitmentInstance, + sc::UnitCommitmentScenario, distribution, )::Nothing - α = rand(rng, distribution, length(instance.buses)) - for t in 1:instance.time - total = sum(bus.load[t] for bus in instance.buses) - den = sum( - bus.load[t] / total * α[i] for - (i, bus) in enumerate(instance.buses) - ) - for (i, bus) in enumerate(instance.buses) + α = rand(rng, distribution, length(sc.buses)) + for t in 1:sc.time + total = sum(bus.load[t] for bus in sc.buses) + den = + sum(bus.load[t] / total * α[i] for (i, bus) in enumerate(sc.buses)) + for (i, bus) in enumerate(sc.buses) bus.load[t] *= α[i] / den end end @@ -156,12 +155,12 @@ end function _randomize_load_profile( rng, - instance::UnitCommitmentInstance, + sc::UnitCommitmentScenario, params::Randomization, )::Nothing # Generate new system load system_load = [1.0] - for t in 2:instance.time + for t in 2:sc.time idx = (t - 1) % length(params.load_profile_mu) + 1 gamma = rand( rng, @@ -169,14 +168,14 @@ function _randomize_load_profile( ) push!(system_load, system_load[t-1] * gamma) end - capacity = sum(maximum(u.max_power) for u in instance.units) + capacity = sum(maximum(u.max_power) for u in sc.units) peak_load = rand(rng, params.peak_load) * capacity system_load = system_load ./ maximum(system_load) .* peak_load # Scale bus loads to match the new system load - prev_system_load = sum(b.load for b in instance.buses) - for b in instance.buses - for t in 1:instance.time + prev_system_load = sum(b.load for b in sc.buses) + for b in sc.buses + for t in 1:sc.time b.load[t] *= system_load[t] / prev_system_load[t] end end @@ -199,15 +198,26 @@ function randomize!( instance::UnitCommitment.UnitCommitmentInstance, method::XavQiuAhm2021.Randomization; rng = MersenneTwister(), +)::Nothing + for sc in instance.scenarios + randomize!(sc; method, rng) + end + return +end + +function randomize!( + sc::UnitCommitment.UnitCommitmentScenario; + method::XavQiuAhm2021.Randomization, + rng = MersenneTwister(), )::Nothing if method.randomize_costs - XavQiuAhm2021._randomize_costs(rng, instance, method.cost) + XavQiuAhm2021._randomize_costs(rng, sc, method.cost) end if method.randomize_load_share - XavQiuAhm2021._randomize_load_share(rng, instance, method.load_share) + XavQiuAhm2021._randomize_load_share(rng, sc, method.load_share) end if method.randomize_load_profile - XavQiuAhm2021._randomize_load_profile(rng, instance, method) + XavQiuAhm2021._randomize_load_profile(rng, sc, method) end return end diff --git a/src/transform/slice.jl b/src/transform/slice.jl index 0ae861f..e1c6799 100644 --- a/src/transform/slice.jl +++ b/src/transform/slice.jl @@ -24,31 +24,33 @@ function slice( )::UnitCommitmentInstance modified = deepcopy(instance) modified.time = length(range) - modified.power_balance_penalty = modified.power_balance_penalty[range] - for r in modified.reserves - r.amount = r.amount[range] - end - for u in modified.units - u.max_power = u.max_power[range] - u.min_power = u.min_power[range] - u.must_run = u.must_run[range] - u.min_power_cost = u.min_power_cost[range] - for s in u.cost_segments - s.mw = s.mw[range] - s.cost = s.cost[range] + for sc in modified.scenarios + sc.power_balance_penalty = sc.power_balance_penalty[range] + for r in sc.reserves + r.amount = r.amount[range] + end + for u in sc.units + u.max_power = u.max_power[range] + u.min_power = u.min_power[range] + u.must_run = u.must_run[range] + u.min_power_cost = u.min_power_cost[range] + for s in u.cost_segments + s.mw = s.mw[range] + s.cost = s.cost[range] + end + end + for b in sc.buses + b.load = b.load[range] + end + for l in sc.lines + l.normal_flow_limit = l.normal_flow_limit[range] + l.emergency_flow_limit = l.emergency_flow_limit[range] + l.flow_limit_penalty = l.flow_limit_penalty[range] + end + for ps in sc.price_sensitive_loads + ps.demand = ps.demand[range] + ps.revenue = ps.revenue[range] end - end - for b in modified.buses - b.load = b.load[range] - end - for l in modified.lines - l.normal_flow_limit = l.normal_flow_limit[range] - l.emergency_flow_limit = l.emergency_flow_limit[range] - l.flow_limit_penalty = l.flow_limit_penalty[range] - end - for ps in modified.price_sensitive_loads - ps.demand = ps.demand[range] - ps.revenue = ps.revenue[range] end return modified end diff --git a/src/validation/validate.jl b/src/validation/validate.jl index 9677661..be8234d 100644 --- a/src/validation/validate.jl +++ b/src/validation/validate.jl @@ -28,6 +28,8 @@ function validate( instance::UnitCommitmentInstance, solution::Union{Dict,OrderedDict}, )::Bool + "Production (MW)" ∈ keys(solution) ? solution = Dict("s1" => solution) : + nothing err_count = 0 err_count += _validate_units(instance, solution) err_count += _validate_reserve_and_demand(instance, solution) @@ -42,358 +44,369 @@ end function _validate_units(instance::UnitCommitmentInstance, solution; tol = 0.01) err_count = 0 - - for unit in instance.units - production = solution["Production (MW)"][unit.name] - reserve = [0.0 for _ in 1:instance.time] - spinning_reserves = [r for r in unit.reserves if r.type == "spinning"] - if !isempty(spinning_reserves) - reserve += sum( - solution["Spinning reserve (MW)"][r.name][unit.name] for - r in spinning_reserves - ) - end - actual_production_cost = solution["Production cost (\$)"][unit.name] - actual_startup_cost = solution["Startup cost (\$)"][unit.name] - is_on = bin(solution["Is on"][unit.name]) - - for t in 1:instance.time - # Auxiliary variables - if t == 1 - is_starting_up = (unit.initial_status < 0) && is_on[t] - is_shutting_down = (unit.initial_status > 0) && !is_on[t] - ramp_up = - max(0, production[t] + reserve[t] - unit.initial_power) - ramp_down = max(0, unit.initial_power - production[t]) - else - is_starting_up = !is_on[t-1] && is_on[t] - is_shutting_down = is_on[t-1] && !is_on[t] - ramp_up = max(0, production[t] + reserve[t] - production[t-1]) - ramp_down = max(0, production[t-1] - production[t]) + for sc in instance.scenarios + for unit in sc.units + production = solution[sc.name]["Production (MW)"][unit.name] + reserve = [0.0 for _ in 1:instance.time] + spinning_reserves = + [r for r in unit.reserves if r.type == "spinning"] + if !isempty(spinning_reserves) + reserve += sum( + solution[sc.name]["Spinning reserve (MW)"][r.name][unit.name] + for r in spinning_reserves + ) end + actual_production_cost = + solution[sc.name]["Production cost (\$)"][unit.name] + actual_startup_cost = + solution[sc.name]["Startup cost (\$)"][unit.name] + is_on = bin(solution[sc.name]["Is on"][unit.name]) + + for t in 1:instance.time + # Auxiliary variables + if t == 1 + is_starting_up = (unit.initial_status < 0) && is_on[t] + is_shutting_down = (unit.initial_status > 0) && !is_on[t] + ramp_up = + max(0, production[t] + reserve[t] - unit.initial_power) + ramp_down = max(0, unit.initial_power - production[t]) + else + is_starting_up = !is_on[t-1] && is_on[t] + is_shutting_down = is_on[t-1] && !is_on[t] + ramp_up = + max(0, production[t] + reserve[t] - production[t-1]) + ramp_down = max(0, production[t-1] - production[t]) + end - # Compute production costs - production_cost, startup_cost = 0, 0 - if is_on[t] - production_cost += unit.min_power_cost[t] - residual = max(0, production[t] - unit.min_power[t]) - for s in unit.cost_segments - cleared = min(residual, s.mw[t]) - production_cost += cleared * s.cost[t] - residual = max(0, residual - s.mw[t]) + # Compute production costs + production_cost, startup_cost = 0, 0 + if is_on[t] + production_cost += unit.min_power_cost[t] + residual = max(0, production[t] - unit.min_power[t]) + for s in unit.cost_segments + cleared = min(residual, s.mw[t]) + production_cost += cleared * s.cost[t] + residual = max(0, residual - s.mw[t]) + end end - end - # Production should be non-negative - if production[t] < -tol - @error @sprintf( - "Unit %s produces negative amount of power at time %d (%.2f)", - unit.name, - t, - production[t] - ) - err_count += 1 - end + # Production should be non-negative + if production[t] < -tol + @error @sprintf( + "Unit %s produces negative amount of power at time %d (%.2f)", + unit.name, + t, + production[t] + ) + err_count += 1 + end - # Verify must-run - if !is_on[t] && unit.must_run[t] - @error @sprintf( - "Must-run unit %s is offline at time %d", - unit.name, - t - ) - err_count += 1 - end + # Verify must-run + if !is_on[t] && unit.must_run[t] + @error @sprintf( + "Must-run unit %s is offline at time %d", + unit.name, + t + ) + err_count += 1 + end - # Verify reserve eligibility - for r in instance.reserves - if r.type == "spinning" - if unit ∉ r.units && - (unit in keys(solution["Spinning reserve (MW)"][r.name])) - @error @sprintf( - "Unit %s is not eligible to provide reserve %s", - unit.name, - r.name, + # Verify reserve eligibility + for r in sc.reserves + if r.type == "spinning" + if unit ∉ r.units && ( + unit in keys( + solution[sc.name]["Spinning reserve (MW)"][r.name], + ) ) - err_count += 1 + @error @sprintf( + "Unit %s is not eligible to provide reserve %s", + unit.name, + r.name, + ) + err_count += 1 + end end end - end - - # If unit is on, must produce at least its minimum power - if is_on[t] && (production[t] < unit.min_power[t] - tol) - @error @sprintf( - "Unit %s produces below its minimum limit at time %d (%.2f < %.2f)", - unit.name, - t, - production[t], - unit.min_power[t] - ) - err_count += 1 - end - - # If unit is on, must produce at most its maximum power - if is_on[t] && - (production[t] + reserve[t] > unit.max_power[t] + tol) - @error @sprintf( - "Unit %s produces above its maximum limit at time %d (%.2f + %.2f> %.2f)", - unit.name, - t, - production[t], - reserve[t], - unit.max_power[t] - ) - err_count += 1 - end - # If unit is off, must produce zero - if !is_on[t] && production[t] + reserve[t] > tol - @error @sprintf( - "Unit %s produces power at time %d while off (%.2f + %.2f > 0)", - unit.name, - t, - production[t], - reserve[t], - ) - err_count += 1 - end - - # Startup limit - if is_starting_up && (ramp_up > unit.startup_limit + tol) - @error @sprintf( - "Unit %s exceeds startup limit at time %d (%.2f > %.2f)", - unit.name, - t, - ramp_up, - unit.startup_limit - ) - err_count += 1 - end - - # Shutdown limit - if is_shutting_down && (ramp_down > unit.shutdown_limit + tol) - @error @sprintf( - "Unit %s exceeds shutdown limit at time %d (%.2f > %.2f)", - unit.name, - t, - ramp_down, - unit.shutdown_limit - ) - err_count += 1 - end - - # Ramp-up limit - if !is_starting_up && - !is_shutting_down && - (ramp_up > unit.ramp_up_limit + tol) - @error @sprintf( - "Unit %s exceeds ramp up limit at time %d (%.2f > %.2f)", - unit.name, - t, - ramp_up, - unit.ramp_up_limit - ) - err_count += 1 - end + # If unit is on, must produce at least its minimum power + if is_on[t] && (production[t] < unit.min_power[t] - tol) + @error @sprintf( + "Unit %s produces below its minimum limit at time %d (%.2f < %.2f)", + unit.name, + t, + production[t], + unit.min_power[t] + ) + err_count += 1 + end - # Ramp-down limit - if !is_starting_up && - !is_shutting_down && - (ramp_down > unit.ramp_down_limit + tol) - @error @sprintf( - "Unit %s exceeds ramp down limit at time %d (%.2f > %.2f)", - unit.name, - t, - ramp_down, - unit.ramp_down_limit - ) - err_count += 1 - end + # If unit is on, must produce at most its maximum power + if is_on[t] && + (production[t] + reserve[t] > unit.max_power[t] + tol) + @error @sprintf( + "Unit %s produces above its maximum limit at time %d (%.2f + %.2f> %.2f)", + unit.name, + t, + production[t], + reserve[t], + unit.max_power[t] + ) + err_count += 1 + end - # Verify startup costs & minimum downtime - if is_starting_up + # If unit is off, must produce zero + if !is_on[t] && production[t] + reserve[t] > tol + @error @sprintf( + "Unit %s produces power at time %d while off (%.2f + %.2f > 0)", + unit.name, + t, + production[t], + reserve[t], + ) + err_count += 1 + end - # Calculate how much time the unit has been offline - time_down = 0 - for k in 1:(t-1) - if !is_on[t-k] - time_down += 1 - else - break - end + # Startup limit + if is_starting_up && (ramp_up > unit.startup_limit + tol) + @error @sprintf( + "Unit %s exceeds startup limit at time %d (%.2f > %.2f)", + unit.name, + t, + ramp_up, + unit.startup_limit + ) + err_count += 1 end - if (t == time_down + 1) && (unit.initial_status < 0) - time_down -= unit.initial_status + + # Shutdown limit + if is_shutting_down && (ramp_down > unit.shutdown_limit + tol) + @error @sprintf( + "Unit %s exceeds shutdown limit at time %d (%.2f > %.2f)", + unit.name, + t, + ramp_down, + unit.shutdown_limit + ) + err_count += 1 end - # Calculate startup costs - for c in unit.startup_categories - if time_down >= c.delay - startup_cost = c.cost - end + # Ramp-up limit + if !is_starting_up && + !is_shutting_down && + (ramp_up > unit.ramp_up_limit + tol) + @error @sprintf( + "Unit %s exceeds ramp up limit at time %d (%.2f > %.2f)", + unit.name, + t, + ramp_up, + unit.ramp_up_limit + ) + err_count += 1 end - # Check minimum downtime - if time_down < unit.min_downtime + # Ramp-down limit + if !is_starting_up && + !is_shutting_down && + (ramp_down > unit.ramp_down_limit + tol) @error @sprintf( - "Unit %s violates minimum downtime at time %d", + "Unit %s exceeds ramp down limit at time %d (%.2f > %.2f)", unit.name, - t + t, + ramp_down, + unit.ramp_down_limit ) err_count += 1 end - end - # Verify minimum uptime - if is_shutting_down + # Verify startup costs & minimum downtime + if is_starting_up + + # Calculate how much time the unit has been offline + time_down = 0 + for k in 1:(t-1) + if !is_on[t-k] + time_down += 1 + else + break + end + end + if (t == time_down + 1) && (unit.initial_status < 0) + time_down -= unit.initial_status + end + + # Calculate startup costs + for c in unit.startup_categories + if time_down >= c.delay + startup_cost = c.cost + end + end - # Calculate how much time the unit has been online - time_up = 0 - for k in 1:(t-1) - if is_on[t-k] - time_up += 1 - else - break + # Check minimum downtime + if time_down < unit.min_downtime + @error @sprintf( + "Unit %s violates minimum downtime at time %d", + unit.name, + t + ) + err_count += 1 end end - if (t == time_up + 1) && (unit.initial_status > 0) - time_up += unit.initial_status + + # Verify minimum uptime + if is_shutting_down + + # Calculate how much time the unit has been online + time_up = 0 + for k in 1:(t-1) + if is_on[t-k] + time_up += 1 + else + break + end + end + if (t == time_up + 1) && (unit.initial_status > 0) + time_up += unit.initial_status + end + + # Check minimum uptime + if time_up < unit.min_uptime + @error @sprintf( + "Unit %s violates minimum uptime at time %d", + unit.name, + t + ) + err_count += 1 + end end - # Check minimum uptime - if time_up < unit.min_uptime + # Verify production costs + if abs(actual_production_cost[t] - production_cost) > 1.00 @error @sprintf( - "Unit %s violates minimum uptime at time %d", + "Unit %s has unexpected production cost at time %d (%.2f should be %.2f)", unit.name, - t + t, + actual_production_cost[t], + production_cost ) err_count += 1 end - end - - # Verify production costs - if abs(actual_production_cost[t] - production_cost) > 1.00 - @error @sprintf( - "Unit %s has unexpected production cost at time %d (%.2f should be %.2f)", - unit.name, - t, - actual_production_cost[t], - production_cost - ) - err_count += 1 - end - # Verify startup costs - if abs(actual_startup_cost[t] - startup_cost) > 1.00 - @error @sprintf( - "Unit %s has unexpected startup cost at time %d (%.2f should be %.2f)", - unit.name, - t, - actual_startup_cost[t], - startup_cost - ) - err_count += 1 + # Verify startup costs + if abs(actual_startup_cost[t] - startup_cost) > 1.00 + @error @sprintf( + "Unit %s has unexpected startup cost at time %d (%.2f should be %.2f)", + unit.name, + t, + actual_startup_cost[t], + startup_cost + ) + err_count += 1 + end end end end - return err_count end function _validate_reserve_and_demand(instance, solution, tol = 0.01) err_count = 0 - for t in 1:instance.time - load_curtail = 0 - fixed_load = sum(b.load[t] for b in instance.buses) - ps_load = 0 - if length(instance.price_sensitive_loads) > 0 - ps_load = sum( - solution["Price-sensitive loads (MW)"][ps.name][t] for - ps in instance.price_sensitive_loads - ) - end - production = - sum(solution["Production (MW)"][g.name][t] for g in instance.units) - if "Load curtail (MW)" in keys(solution) - load_curtail = sum( - solution["Load curtail (MW)"][b.name][t] for - b in instance.buses - ) - end - balance = fixed_load - load_curtail - production + ps_load - - # Verify that production equals demand - if abs(balance) > tol - @error @sprintf( - "Non-zero power balance at time %d (%.2f + %.2f - %.2f - %.2f != 0)", - t, - fixed_load, - ps_load, - load_curtail, - production, + for sc in instance.scenarios + for t in 1:instance.time + load_curtail = 0 + fixed_load = sum(b.load[t] for b in sc.buses) + ps_load = 0 + if length(sc.price_sensitive_loads) > 0 + ps_load = sum( + solution[sc.name]["Price-sensitive loads (MW)"][ps.name][t] + for ps in sc.price_sensitive_loads + ) + end + production = sum( + solution[sc.name]["Production (MW)"][g.name][t] for + g in sc.units ) - err_count += 1 - end + if "Load curtail (MW)" in keys(solution) + load_curtail = sum( + solution[sc.name]["Load curtail (MW)"][b.name][t] for + b in sc.buses + ) + end + balance = fixed_load - load_curtail - production + ps_load - # Verify reserves - for r in instance.reserves - if r.type == "spinning" - provided = sum( - solution["Spinning reserve (MW)"][r.name][g.name][t] for - g in r.units + # Verify that production equals demand + if abs(balance) > tol + @error @sprintf( + "Non-zero power balance at time %d (%.2f + %.2f - %.2f - %.2f != 0)", + t, + fixed_load, + ps_load, + load_curtail, + production, ) - shortfall = - solution["Spinning reserve shortfall (MW)"][r.name][t] - required = r.amount[t] + err_count += 1 + end - if provided + shortfall < required - tol - @error @sprintf( - "Insufficient reserve %s at time %d (%.2f + %.2f < %.2f)", - r.name, - t, - provided, - shortfall, - required, + # Verify reserves + for r in sc.reserves + if r.type == "spinning" + provided = sum( + solution[sc.name]["Spinning reserve (MW)"][r.name][g.name][t] + for g in r.units ) - end - elseif r.type == "up-frp" - upflexiramp = sum( - solution["Up-flexiramp (MW)"][r.name][g.name][t] for - g in r.units - ) - upflexiramp_shortfall = - solution["Up-flexiramp shortfall (MW)"][r.name][t] + shortfall = + solution[sc.name]["Spinning reserve shortfall (MW)"][r.name][t] + required = r.amount[t] - if upflexiramp + upflexiramp_shortfall < r.amount[t] - tol - @error @sprintf( - "Insufficient up-flexiramp at time %d (%.2f + %.2f < %.2f)", - t, - upflexiramp, - upflexiramp_shortfall, - r.amount[t], + if provided + shortfall < required - tol + @error @sprintf( + "Insufficient reserve %s at time %d (%.2f + %.2f < %.2f)", + r.name, + t, + provided, + shortfall, + required, + ) + end + elseif r.type == "flexiramp" + upflexiramp = sum( + solution[sc.name]["Up-flexiramp (MW)"][r.name][g.name][t] + for g in r.units ) - err_count += 1 - end - elseif r.type == "down-frp" - dwflexiramp = sum( - solution["Down-flexiramp (MW)"][r.name][g.name][t] for - g in r.units - ) - dwflexiramp_shortfall = - solution["Down-flexiramp shortfall (MW)"][r.name][t] + upflexiramp_shortfall = + solution[sc.name]["Up-flexiramp shortfall (MW)"][r.name][t] - if dwflexiramp + dwflexiramp_shortfall < r.amount[t] - tol - @error @sprintf( - "Insufficient down-flexiramp at time %d (%.2f + %.2f < %.2f)", - t, - dwflexiramp, - dwflexiramp_shortfall, - r.amount[t], + if upflexiramp + upflexiramp_shortfall < r.amount[t] - tol + @error @sprintf( + "Insufficient up-flexiramp at time %d (%.2f + %.2f < %.2f)", + t, + upflexiramp, + upflexiramp_shortfall, + r.amount[t], + ) + err_count += 1 + end + + dwflexiramp = sum( + solution[sc.name]["Down-flexiramp (MW)"][r.name][g.name][t] + for g in r.units ) - err_count += 1 + dwflexiramp_shortfall = + solution[sc.name]["Down-flexiramp shortfall (MW)"][r.name][t] + + if dwflexiramp + dwflexiramp_shortfall < r.amount[t] - tol + @error @sprintf( + "Insufficient down-flexiramp at time %d (%.2f + %.2f < %.2f)", + t, + dwflexiramp, + dwflexiramp_shortfall, + r.amount[t], + ) + err_count += 1 + end + else + error("Unknown reserve type: $(r.type)") end - else - error("Unknown reserve type: $(r.type)") end end end diff --git a/test/Project.toml b/test/Project.toml index b87293d..e65d966 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -2,6 +2,7 @@ Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" diff --git a/test/instance/migrate_test.jl b/test/instance/migrate_test.jl index 2ea6d46..ab8fb59 100644 --- a/test/instance/migrate_test.jl +++ b/test/instance/migrate_test.jl @@ -6,13 +6,17 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip @testset "read v0.2" begin instance = UnitCommitment.read("$FIXTURES/ucjl-0.2.json.gz") - @test length(instance.reserves_by_name["r1"].amount) == 4 - @test instance.units_by_name["g2"].reserves[1].name == "r1" + for sc in instance.scenarios + @test length(sc.reserves_by_name["r1"].amount) == 4 + @test sc.units_by_name["g2"].reserves[1].name == "r1" + end end @testset "read v0.3" begin instance = UnitCommitment.read("$FIXTURES/ucjl-0.3.json.gz") - @test length(instance.units) == 6 - @test length(instance.buses) == 14 - @test length(instance.lines) == 20 + for sc in instance.scenarios + @test length(sc.units) == 6 + @test length(sc.buses) == 14 + @test length(sc.lines) == 20 + end end diff --git a/test/instance/read_test.jl b/test/instance/read_test.jl index cc9e332..64003a2 100644 --- a/test/instance/read_test.jl +++ b/test/instance/read_test.jl @@ -6,115 +6,116 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip @testset "read_benchmark" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") + for sc in instance.scenarios + @test length(sc.lines) == 20 + @test length(sc.buses) == 14 + @test length(sc.units) == 6 + @test length(sc.contingencies) == 19 + @test length(sc.price_sensitive_loads) == 1 + @test instance.time == 4 - @test length(instance.lines) == 20 - @test length(instance.buses) == 14 - @test length(instance.units) == 6 - @test length(instance.contingencies) == 19 - @test length(instance.price_sensitive_loads) == 1 - @test instance.time == 4 + @test sc.lines[5].name == "l5" + @test sc.lines[5].source.name == "b2" + @test sc.lines[5].target.name == "b5" + @test sc.lines[5].reactance ≈ 0.17388 + @test sc.lines[5].susceptance ≈ 10.037550333 + @test sc.lines[5].normal_flow_limit == [1e8 for t in 1:4] + @test sc.lines[5].emergency_flow_limit == [1e8 for t in 1:4] + @test sc.lines[5].flow_limit_penalty == [5e3 for t in 1:4] + @test sc.lines_by_name["l5"].name == "l5" - @test instance.lines[5].name == "l5" - @test instance.lines[5].source.name == "b2" - @test instance.lines[5].target.name == "b5" - @test instance.lines[5].reactance ≈ 0.17388 - @test instance.lines[5].susceptance ≈ 10.037550333 - @test instance.lines[5].normal_flow_limit == [1e8 for t in 1:4] - @test instance.lines[5].emergency_flow_limit == [1e8 for t in 1:4] - @test instance.lines[5].flow_limit_penalty == [5e3 for t in 1:4] - @test instance.lines_by_name["l5"].name == "l5" + @test sc.lines[1].name == "l1" + @test sc.lines[1].source.name == "b1" + @test sc.lines[1].target.name == "b2" + @test sc.lines[1].reactance ≈ 0.059170 + @test sc.lines[1].susceptance ≈ 29.496860773945 + @test sc.lines[1].normal_flow_limit == [300.0 for t in 1:4] + @test sc.lines[1].emergency_flow_limit == [400.0 for t in 1:4] + @test sc.lines[1].flow_limit_penalty == [1e3 for t in 1:4] - @test instance.lines[1].name == "l1" - @test instance.lines[1].source.name == "b1" - @test instance.lines[1].target.name == "b2" - @test instance.lines[1].reactance ≈ 0.059170 - @test instance.lines[1].susceptance ≈ 29.496860773945 - @test instance.lines[1].normal_flow_limit == [300.0 for t in 1:4] - @test instance.lines[1].emergency_flow_limit == [400.0 for t in 1:4] - @test instance.lines[1].flow_limit_penalty == [1e3 for t in 1:4] + @test sc.buses[9].name == "b9" + @test sc.buses[9].load == [35.36638, 33.25495, 31.67138, 31.14353] + @test sc.buses_by_name["b9"].name == "b9" - @test instance.buses[9].name == "b9" - @test instance.buses[9].load == [35.36638, 33.25495, 31.67138, 31.14353] - @test instance.buses_by_name["b9"].name == "b9" + @test sc.reserves[1].name == "r1" + @test sc.reserves[1].type == "spinning" + @test sc.reserves[1].amount == [100.0, 100.0, 100.0, 100.0] + @test sc.reserves_by_name["r1"].name == "r1" - @test instance.reserves[1].name == "r1" - @test instance.reserves[1].type == "spinning" - @test instance.reserves[1].amount == [100.0, 100.0, 100.0, 100.0] - @test instance.reserves_by_name["r1"].name == "r1" + unit = sc.units[1] + @test unit.name == "g1" + @test unit.bus.name == "b1" + @test unit.ramp_up_limit == 1e6 + @test unit.ramp_down_limit == 1e6 + @test unit.startup_limit == 1e6 + @test unit.shutdown_limit == 1e6 + @test unit.must_run == [false for t in 1:4] + @test unit.min_power_cost == [1400.0 for t in 1:4] + @test unit.min_uptime == 1 + @test unit.min_downtime == 1 + for t in 1:1 + @test unit.cost_segments[1].mw[t] == 10.0 + @test unit.cost_segments[2].mw[t] == 20.0 + @test unit.cost_segments[3].mw[t] == 5.0 + @test unit.cost_segments[1].cost[t] ≈ 20.0 + @test unit.cost_segments[2].cost[t] ≈ 30.0 + @test unit.cost_segments[3].cost[t] ≈ 40.0 + end + @test length(unit.startup_categories) == 3 + @test unit.startup_categories[1].delay == 1 + @test unit.startup_categories[2].delay == 2 + @test unit.startup_categories[3].delay == 3 + @test unit.startup_categories[1].cost == 1000.0 + @test unit.startup_categories[2].cost == 1500.0 + @test unit.startup_categories[3].cost == 2000.0 + @test length(unit.reserves) == 0 + @test sc.units_by_name["g1"].name == "g1" - unit = instance.units[1] - @test unit.name == "g1" - @test unit.bus.name == "b1" - @test unit.ramp_up_limit == 1e6 - @test unit.ramp_down_limit == 1e6 - @test unit.startup_limit == 1e6 - @test unit.shutdown_limit == 1e6 - @test unit.must_run == [false for t in 1:4] - @test unit.min_power_cost == [1400.0 for t in 1:4] - @test unit.min_uptime == 1 - @test unit.min_downtime == 1 - for t in 1:1 - @test unit.cost_segments[1].mw[t] == 10.0 - @test unit.cost_segments[2].mw[t] == 20.0 - @test unit.cost_segments[3].mw[t] == 5.0 - @test unit.cost_segments[1].cost[t] ≈ 20.0 - @test unit.cost_segments[2].cost[t] ≈ 30.0 - @test unit.cost_segments[3].cost[t] ≈ 40.0 - end - @test length(unit.startup_categories) == 3 - @test unit.startup_categories[1].delay == 1 - @test unit.startup_categories[2].delay == 2 - @test unit.startup_categories[3].delay == 3 - @test unit.startup_categories[1].cost == 1000.0 - @test unit.startup_categories[2].cost == 1500.0 - @test unit.startup_categories[3].cost == 2000.0 - @test length(unit.reserves) == 0 - @test instance.units_by_name["g1"].name == "g1" + unit = sc.units[2] + @test unit.name == "g2" + @test unit.must_run == [false for t in 1:4] + @test length(unit.reserves) == 1 - unit = instance.units[2] - @test unit.name == "g2" - @test unit.must_run == [false for t in 1:4] - @test length(unit.reserves) == 1 + unit = sc.units[3] + @test unit.name == "g3" + @test unit.bus.name == "b3" + @test unit.ramp_up_limit == 70.0 + @test unit.ramp_down_limit == 70.0 + @test unit.startup_limit == 70.0 + @test unit.shutdown_limit == 70.0 + @test unit.must_run == [true for t in 1:4] + @test unit.min_power_cost == [0.0 for t in 1:4] + @test unit.min_uptime == 1 + @test unit.min_downtime == 1 + for t in 1:4 + @test unit.cost_segments[1].mw[t] ≈ 33 + @test unit.cost_segments[2].mw[t] ≈ 33 + @test unit.cost_segments[3].mw[t] ≈ 34 + @test unit.cost_segments[1].cost[t] ≈ 33.75 + @test unit.cost_segments[2].cost[t] ≈ 38.04 + @test unit.cost_segments[3].cost[t] ≈ 44.77853 + end + @test length(unit.reserves) == 1 + @test unit.reserves[1].name == "r1" - unit = instance.units[3] - @test unit.name == "g3" - @test unit.bus.name == "b3" - @test unit.ramp_up_limit == 70.0 - @test unit.ramp_down_limit == 70.0 - @test unit.startup_limit == 70.0 - @test unit.shutdown_limit == 70.0 - @test unit.must_run == [true for t in 1:4] - @test unit.min_power_cost == [0.0 for t in 1:4] - @test unit.min_uptime == 1 - @test unit.min_downtime == 1 - for t in 1:4 - @test unit.cost_segments[1].mw[t] ≈ 33 - @test unit.cost_segments[2].mw[t] ≈ 33 - @test unit.cost_segments[3].mw[t] ≈ 34 - @test unit.cost_segments[1].cost[t] ≈ 33.75 - @test unit.cost_segments[2].cost[t] ≈ 38.04 - @test unit.cost_segments[3].cost[t] ≈ 44.77853 - end - @test length(unit.reserves) == 1 - @test unit.reserves[1].name == "r1" - - @test instance.contingencies[1].lines == [instance.lines[1]] - @test instance.contingencies[1].units == [] - @test instance.contingencies[1].name == "c1" - @test instance.contingencies_by_name["c1"].name == "c1" + @test sc.contingencies[1].lines == [sc.lines[1]] + @test sc.contingencies[1].units == [] + @test sc.contingencies[1].name == "c1" + @test sc.contingencies_by_name["c1"].name == "c1" - load = instance.price_sensitive_loads[1] - @test load.name == "ps1" - @test load.bus.name == "b3" - @test load.revenue == [100.0 for t in 1:4] - @test load.demand == [50.0 for t in 1:4] - @test instance.price_sensitive_loads_by_name["ps1"].name == "ps1" + load = sc.price_sensitive_loads[1] + @test load.name == "ps1" + @test load.bus.name == "b3" + @test load.revenue == [100.0 for t in 1:4] + @test load.demand == [50.0 for t in 1:4] + @test sc.price_sensitive_loads_by_name["ps1"].name == "ps1" + end end @testset "read_benchmark sub-hourly" begin instance = UnitCommitment.read("$FIXTURES/case14-sub-hourly.json.gz") @test instance.time == 4 - unit = instance.units[1] + unit = instance.scenarios[1].units[1] @test unit.name == "g1" @test unit.min_uptime == 2 @test unit.min_downtime == 2 diff --git a/test/solution/methods/XavQiuWanThi19/filter_test.jl b/test/solution/methods/XavQiuWanThi19/filter_test.jl index e4a8b48..fa964d3 100644 --- a/test/solution/methods/XavQiuWanThi19/filter_test.jl +++ b/test/solution/methods/XavQiuWanThi19/filter_test.jl @@ -8,76 +8,77 @@ import UnitCommitment: _Violation, _offer, _query @testset "_ViolationFilter" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") filter = UnitCommitment._ViolationFilter(max_per_line = 1, max_total = 2) + for sc in instance.scenarios + _offer( + filter, + _Violation( + time = 1, + monitored_line = sc.lines[1], + outage_line = nothing, + amount = 100.0, + ), + ) + _offer( + filter, + _Violation( + time = 1, + monitored_line = sc.lines[1], + outage_line = sc.lines[1], + amount = 300.0, + ), + ) + _offer( + filter, + _Violation( + time = 1, + monitored_line = sc.lines[1], + outage_line = sc.lines[5], + amount = 500.0, + ), + ) + _offer( + filter, + _Violation( + time = 1, + monitored_line = sc.lines[1], + outage_line = sc.lines[4], + amount = 400.0, + ), + ) + _offer( + filter, + _Violation( + time = 1, + monitored_line = sc.lines[2], + outage_line = sc.lines[1], + amount = 200.0, + ), + ) + _offer( + filter, + _Violation( + time = 1, + monitored_line = sc.lines[2], + outage_line = sc.lines[8], + amount = 100.0, + ), + ) - _offer( - filter, - _Violation( - time = 1, - monitored_line = instance.lines[1], - outage_line = nothing, - amount = 100.0, - ), - ) - _offer( - filter, - _Violation( - time = 1, - monitored_line = instance.lines[1], - outage_line = instance.lines[1], - amount = 300.0, - ), - ) - _offer( - filter, - _Violation( - time = 1, - monitored_line = instance.lines[1], - outage_line = instance.lines[5], - amount = 500.0, - ), - ) - _offer( - filter, - _Violation( - time = 1, - monitored_line = instance.lines[1], - outage_line = instance.lines[4], - amount = 400.0, - ), - ) - _offer( - filter, - _Violation( - time = 1, - monitored_line = instance.lines[2], - outage_line = instance.lines[1], - amount = 200.0, - ), - ) - _offer( - filter, - _Violation( - time = 1, - monitored_line = instance.lines[2], - outage_line = instance.lines[8], - amount = 100.0, - ), - ) - - actual = _query(filter) - expected = [ - _Violation( - time = 1, - monitored_line = instance.lines[2], - outage_line = instance.lines[1], - amount = 200.0, - ), - _Violation( - time = 1, - monitored_line = instance.lines[1], - outage_line = instance.lines[5], - amount = 500.0, - ), - ] - @test actual == expected + actual = _query(filter) + expected = [ + _Violation( + time = 1, + monitored_line = sc.lines[2], + outage_line = sc.lines[1], + amount = 200.0, + ), + _Violation( + time = 1, + monitored_line = sc.lines[1], + outage_line = sc.lines[5], + amount = 500.0, + ), + ] + @test actual == expected + end end diff --git a/test/solution/methods/XavQiuWanThi19/find_test.jl b/test/solution/methods/XavQiuWanThi19/find_test.jl index 27110cb..e4be83f 100644 --- a/test/solution/methods/XavQiuWanThi19/find_test.jl +++ b/test/solution/methods/XavQiuWanThi19/find_test.jl @@ -7,29 +7,32 @@ import UnitCommitment: _Violation, _offer, _query @testset "find_violations" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") - for line in instance.lines, t in 1:instance.time - line.normal_flow_limit[t] = 1.0 - line.emergency_flow_limit[t] = 1.0 + for sc in instance.scenarios + for line in sc.lines, t in 1:instance.time + line.normal_flow_limit[t] = 1.0 + line.emergency_flow_limit[t] = 1.0 + end + isf = UnitCommitment._injection_shift_factors( + lines = sc.lines, + buses = sc.buses, + ) + lodf = UnitCommitment._line_outage_factors( + lines = sc.lines, + buses = sc.buses, + isf = isf, + ) + inj = [1000.0 for b in 1:13, t in 1:instance.time] + overflow = [0.0 for l in sc.lines, t in 1:instance.time] + violations = UnitCommitment._find_violations( + instance = instance, + sc = sc, + net_injections = inj, + overflow = overflow, + isf = isf, + lodf = lodf, + max_per_line = 1, + max_per_period = 5, + ) + @test length(violations) == 20 end - isf = UnitCommitment._injection_shift_factors( - lines = instance.lines, - buses = instance.buses, - ) - lodf = UnitCommitment._line_outage_factors( - lines = instance.lines, - buses = instance.buses, - isf = isf, - ) - inj = [1000.0 for b in 1:13, t in 1:instance.time] - overflow = [0.0 for l in instance.lines, t in 1:instance.time] - violations = UnitCommitment._find_violations( - instance = instance, - net_injections = inj, - overflow = overflow, - isf = isf, - lodf = lodf, - max_per_line = 1, - max_per_period = 5, - ) - @test length(violations) == 20 end diff --git a/test/solution/methods/XavQiuWanThi19/sensitivity_test.jl b/test/solution/methods/XavQiuWanThi19/sensitivity_test.jl index 5d64a15..8b3d93f 100644 --- a/test/solution/methods/XavQiuWanThi19/sensitivity_test.jl +++ b/test/solution/methods/XavQiuWanThi19/sensitivity_test.jl @@ -6,137 +6,145 @@ using UnitCommitment, Test, LinearAlgebra @testset "_susceptance_matrix" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") - actual = UnitCommitment._susceptance_matrix(instance.lines) - @test size(actual) == (20, 20) - expected = Diagonal([ - 29.5, - 7.83, - 8.82, - 9.9, - 10.04, - 10.2, - 41.45, - 8.35, - 3.14, - 6.93, - 8.77, - 6.82, - 13.4, - 9.91, - 15.87, - 20.65, - 6.46, - 9.09, - 8.73, - 5.02, - ]) - @test round.(actual, digits = 2) == expected + for sc in instance.scenarios + actual = UnitCommitment._susceptance_matrix(sc.lines) + @test size(actual) == (20, 20) + expected = Diagonal([ + 29.5, + 7.83, + 8.82, + 9.9, + 10.04, + 10.2, + 41.45, + 8.35, + 3.14, + 6.93, + 8.77, + 6.82, + 13.4, + 9.91, + 15.87, + 20.65, + 6.46, + 9.09, + 8.73, + 5.02, + ]) + @test round.(actual, digits = 2) == expected + end end @testset "_reduced_incidence_matrix" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") - actual = UnitCommitment._reduced_incidence_matrix( - lines = instance.lines, - buses = instance.buses, - ) - @test size(actual) == (20, 13) - @test actual[1, 1] == -1.0 - @test actual[3, 1] == 1.0 - @test actual[4, 1] == 1.0 - @test actual[5, 1] == 1.0 - @test actual[3, 2] == -1.0 - @test actual[6, 2] == 1.0 - @test actual[4, 3] == -1.0 - @test actual[6, 3] == -1.0 - @test actual[7, 3] == 1.0 - @test actual[8, 3] == 1.0 - @test actual[9, 3] == 1.0 - @test actual[2, 4] == -1.0 - @test actual[5, 4] == -1.0 - @test actual[7, 4] == -1.0 - @test actual[10, 4] == 1.0 - @test actual[10, 5] == -1.0 - @test actual[11, 5] == 1.0 - @test actual[12, 5] == 1.0 - @test actual[13, 5] == 1.0 - @test actual[8, 6] == -1.0 - @test actual[14, 6] == 1.0 - @test actual[15, 6] == 1.0 - @test actual[14, 7] == -1.0 - @test actual[9, 8] == -1.0 - @test actual[15, 8] == -1.0 - @test actual[16, 8] == 1.0 - @test actual[17, 8] == 1.0 - @test actual[16, 9] == -1.0 - @test actual[18, 9] == 1.0 - @test actual[11, 10] == -1.0 - @test actual[18, 10] == -1.0 - @test actual[12, 11] == -1.0 - @test actual[19, 11] == 1.0 - @test actual[13, 12] == -1.0 - @test actual[19, 12] == -1.0 - @test actual[20, 12] == 1.0 - @test actual[17, 13] == -1.0 - @test actual[20, 13] == -1.0 + for sc in instance.scenarios + actual = UnitCommitment._reduced_incidence_matrix( + lines = sc.lines, + buses = sc.buses, + ) + @test size(actual) == (20, 13) + @test actual[1, 1] == -1.0 + @test actual[3, 1] == 1.0 + @test actual[4, 1] == 1.0 + @test actual[5, 1] == 1.0 + @test actual[3, 2] == -1.0 + @test actual[6, 2] == 1.0 + @test actual[4, 3] == -1.0 + @test actual[6, 3] == -1.0 + @test actual[7, 3] == 1.0 + @test actual[8, 3] == 1.0 + @test actual[9, 3] == 1.0 + @test actual[2, 4] == -1.0 + @test actual[5, 4] == -1.0 + @test actual[7, 4] == -1.0 + @test actual[10, 4] == 1.0 + @test actual[10, 5] == -1.0 + @test actual[11, 5] == 1.0 + @test actual[12, 5] == 1.0 + @test actual[13, 5] == 1.0 + @test actual[8, 6] == -1.0 + @test actual[14, 6] == 1.0 + @test actual[15, 6] == 1.0 + @test actual[14, 7] == -1.0 + @test actual[9, 8] == -1.0 + @test actual[15, 8] == -1.0 + @test actual[16, 8] == 1.0 + @test actual[17, 8] == 1.0 + @test actual[16, 9] == -1.0 + @test actual[18, 9] == 1.0 + @test actual[11, 10] == -1.0 + @test actual[18, 10] == -1.0 + @test actual[12, 11] == -1.0 + @test actual[19, 11] == 1.0 + @test actual[13, 12] == -1.0 + @test actual[19, 12] == -1.0 + @test actual[20, 12] == 1.0 + @test actual[17, 13] == -1.0 + @test actual[20, 13] == -1.0 + end end @testset "_injection_shift_factors" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") - actual = UnitCommitment._injection_shift_factors( - lines = instance.lines, - buses = instance.buses, - ) - @test size(actual) == (20, 13) - @test round.(actual, digits = 2) == [ - -0.84 -0.75 -0.67 -0.61 -0.63 -0.66 -0.66 -0.65 -0.65 -0.64 -0.63 -0.63 -0.64 - -0.16 -0.25 -0.33 -0.39 -0.37 -0.34 -0.34 -0.35 -0.35 -0.36 -0.37 -0.37 -0.36 - 0.03 -0.53 -0.15 -0.1 -0.12 -0.14 -0.14 -0.14 -0.13 -0.13 -0.12 -0.12 -0.13 - 0.06 -0.14 -0.32 -0.22 -0.25 -0.3 -0.3 -0.29 -0.28 -0.27 -0.25 -0.26 -0.27 - 0.08 -0.07 -0.2 -0.29 -0.26 -0.22 -0.22 -0.22 -0.23 -0.25 -0.26 -0.26 -0.24 - 0.03 0.47 -0.15 -0.1 -0.12 -0.14 -0.14 -0.14 -0.13 -0.13 -0.12 -0.12 -0.13 - 0.08 0.31 0.5 -0.3 -0.03 0.36 0.36 0.28 0.23 0.1 -0.0 0.02 0.17 - 0.0 0.01 0.02 -0.01 -0.22 -0.63 -0.63 -0.45 -0.41 -0.32 -0.24 -0.25 -0.36 - 0.0 0.01 0.01 -0.01 -0.12 -0.17 -0.17 -0.26 -0.24 -0.18 -0.14 -0.14 -0.21 - -0.0 -0.02 -0.03 0.02 -0.66 -0.2 -0.2 -0.29 -0.36 -0.5 -0.63 -0.61 -0.43 - -0.0 -0.01 -0.02 0.01 0.21 -0.12 -0.12 -0.17 -0.28 -0.53 0.18 0.15 -0.03 - -0.0 -0.0 -0.0 0.0 0.03 -0.02 -0.02 -0.03 -0.02 0.01 -0.52 -0.17 -0.09 - -0.0 -0.01 -0.01 0.01 0.11 -0.06 -0.06 -0.09 -0.05 0.02 -0.28 -0.59 -0.31 - -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 - 0.0 0.01 0.02 -0.01 -0.22 0.37 0.37 -0.45 -0.41 -0.32 -0.24 -0.25 -0.36 - 0.0 0.01 0.02 -0.01 -0.21 0.12 0.12 0.17 -0.72 -0.47 -0.18 -0.15 0.03 - 0.0 0.01 0.01 -0.01 -0.14 0.08 0.08 0.12 0.07 -0.03 -0.2 -0.24 -0.6 - 0.0 0.01 0.02 -0.01 -0.21 0.12 0.12 0.17 0.28 -0.47 -0.18 -0.15 0.03 - -0.0 -0.0 -0.0 0.0 0.03 -0.02 -0.02 -0.03 -0.02 0.01 0.48 -0.17 -0.09 - -0.0 -0.01 -0.01 0.01 0.14 -0.08 -0.08 -0.12 -0.07 0.03 0.2 0.24 -0.4 - ] + for sc in instance.scenarios + actual = UnitCommitment._injection_shift_factors( + lines = sc.lines, + buses = sc.buses, + ) + @test size(actual) == (20, 13) + @test round.(actual, digits = 2) == [ + -0.84 -0.75 -0.67 -0.61 -0.63 -0.66 -0.66 -0.65 -0.65 -0.64 -0.63 -0.63 -0.64 + -0.16 -0.25 -0.33 -0.39 -0.37 -0.34 -0.34 -0.35 -0.35 -0.36 -0.37 -0.37 -0.36 + 0.03 -0.53 -0.15 -0.1 -0.12 -0.14 -0.14 -0.14 -0.13 -0.13 -0.12 -0.12 -0.13 + 0.06 -0.14 -0.32 -0.22 -0.25 -0.3 -0.3 -0.29 -0.28 -0.27 -0.25 -0.26 -0.27 + 0.08 -0.07 -0.2 -0.29 -0.26 -0.22 -0.22 -0.22 -0.23 -0.25 -0.26 -0.26 -0.24 + 0.03 0.47 -0.15 -0.1 -0.12 -0.14 -0.14 -0.14 -0.13 -0.13 -0.12 -0.12 -0.13 + 0.08 0.31 0.5 -0.3 -0.03 0.36 0.36 0.28 0.23 0.1 -0.0 0.02 0.17 + 0.0 0.01 0.02 -0.01 -0.22 -0.63 -0.63 -0.45 -0.41 -0.32 -0.24 -0.25 -0.36 + 0.0 0.01 0.01 -0.01 -0.12 -0.17 -0.17 -0.26 -0.24 -0.18 -0.14 -0.14 -0.21 + -0.0 -0.02 -0.03 0.02 -0.66 -0.2 -0.2 -0.29 -0.36 -0.5 -0.63 -0.61 -0.43 + -0.0 -0.01 -0.02 0.01 0.21 -0.12 -0.12 -0.17 -0.28 -0.53 0.18 0.15 -0.03 + -0.0 -0.0 -0.0 0.0 0.03 -0.02 -0.02 -0.03 -0.02 0.01 -0.52 -0.17 -0.09 + -0.0 -0.01 -0.01 0.01 0.11 -0.06 -0.06 -0.09 -0.05 0.02 -0.28 -0.59 -0.31 + -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 + 0.0 0.01 0.02 -0.01 -0.22 0.37 0.37 -0.45 -0.41 -0.32 -0.24 -0.25 -0.36 + 0.0 0.01 0.02 -0.01 -0.21 0.12 0.12 0.17 -0.72 -0.47 -0.18 -0.15 0.03 + 0.0 0.01 0.01 -0.01 -0.14 0.08 0.08 0.12 0.07 -0.03 -0.2 -0.24 -0.6 + 0.0 0.01 0.02 -0.01 -0.21 0.12 0.12 0.17 0.28 -0.47 -0.18 -0.15 0.03 + -0.0 -0.0 -0.0 0.0 0.03 -0.02 -0.02 -0.03 -0.02 0.01 0.48 -0.17 -0.09 + -0.0 -0.01 -0.01 0.01 0.14 -0.08 -0.08 -0.12 -0.07 0.03 0.2 0.24 -0.4 + ] + end end @testset "_line_outage_factors" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") - isf_before = UnitCommitment._injection_shift_factors( - lines = instance.lines, - buses = instance.buses, - ) - lodf = UnitCommitment._line_outage_factors( - lines = instance.lines, - buses = instance.buses, - isf = isf_before, - ) - for contingency in instance.contingencies - for lc in contingency.lines - prev_susceptance = lc.susceptance - lc.susceptance = 0.0 - isf_after = UnitCommitment._injection_shift_factors( - lines = instance.lines, - buses = instance.buses, - ) - lc.susceptance = prev_susceptance - for lm in instance.lines - expected = isf_after[lm.offset, :] - actual = - isf_before[lm.offset, :] + - lodf[lm.offset, lc.offset] * isf_before[lc.offset, :] - @test norm(expected - actual) < 1e-6 + for sc in instance.scenarios + isf_before = UnitCommitment._injection_shift_factors( + lines = sc.lines, + buses = sc.buses, + ) + lodf = UnitCommitment._line_outage_factors( + lines = sc.lines, + buses = sc.buses, + isf = isf_before, + ) + for contingency in sc.contingencies + for lc in contingency.lines + prev_susceptance = lc.susceptance + lc.susceptance = 0.0 + isf_after = UnitCommitment._injection_shift_factors( + lines = sc.lines, + buses = sc.buses, + ) + lc.susceptance = prev_susceptance + for lm in sc.lines + expected = isf_after[lm.offset, :] + actual = + isf_before[lm.offset, :] + + lodf[lm.offset, lc.offset] * isf_before[lc.offset, :] + @test norm(expected - actual) < 1e-6 + end end end end diff --git a/test/transform/initcond_test.jl b/test/transform/initcond_test.jl index 56c9831..d01e311 100644 --- a/test/transform/initcond_test.jl +++ b/test/transform/initcond_test.jl @@ -8,20 +8,21 @@ using UnitCommitment, Cbc, JuMP # Load instance instance = UnitCommitment.read("$FIXTURES/case118-initcond.json.gz") optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0) + for sc in instance.scenarios + # All units should have unknown initial conditions + for g in sc.units + @test g.initial_power === nothing + @test g.initial_status === nothing + end - # All units should have unknown initial conditions - for g in instance.units - @test g.initial_power === nothing - @test g.initial_status === nothing - end - - # Generate initial conditions - UnitCommitment.generate_initial_conditions!(instance, optimizer) + # Generate initial conditions + UnitCommitment.generate_initial_conditions!(sc, optimizer) - # All units should now have known initial conditions - for g in instance.units - @test g.initial_power !== nothing - @test g.initial_status !== nothing + # All units should now have known initial conditions + for g in sc.units + @test g.initial_power !== nothing + @test g.initial_status !== nothing + end end # TODO: Check that initial conditions are feasible diff --git a/test/transform/randomize/XavQiuAhm2021_test.jl b/test/transform/randomize/XavQiuAhm2021_test.jl index 6d35829..1717e21 100644 --- a/test/transform/randomize/XavQiuAhm2021_test.jl +++ b/test/transform/randomize/XavQiuAhm2021_test.jl @@ -10,60 +10,78 @@ using Random using UnitCommitment, Cbc, JuMP get_instance() = UnitCommitment.read_benchmark("matpower/case118/2017-02-01") -system_load(instance) = sum(b.load for b in instance.buses) +system_load(sc) = sum(b.load for b in sc.buses) test_approx(x, y) = @test isapprox(x, y, atol = 1e-3) @testset "XavQiuAhm2021" begin @testset "cost and load share" begin instance = get_instance() + for sc in instance.scenarios + # Check original costs + unit = sc.units[10] + test_approx(unit.min_power_cost[1], 825.023) + test_approx(unit.cost_segments[1].cost[1], 36.659) + test_approx(unit.startup_categories[1].cost[1], 7570.42) - # Check original costs - unit = instance.units[10] - test_approx(unit.min_power_cost[1], 825.023) - test_approx(unit.cost_segments[1].cost[1], 36.659) - test_approx(unit.startup_categories[1].cost[1], 7570.42) + # Check original load share + bus = sc.buses[1] + prev_system_load = system_load(sc) + test_approx(bus.load[1] / prev_system_load[1], 0.012) - # Check original load share - bus = instance.buses[1] - prev_system_load = system_load(instance) - test_approx(bus.load[1] / prev_system_load[1], 0.012) + randomize!( + sc, + method = XavQiuAhm2021.Randomization( + randomize_load_profile = false, + ), + rng = MersenneTwister(42), + ) - randomize!( - instance, - method = XavQiuAhm2021.Randomization( - randomize_load_profile = false, - ), - rng = MersenneTwister(42), - ) + # Check randomized costs + test_approx(unit.min_power_cost[1], 831.977) + test_approx(unit.cost_segments[1].cost[1], 36.968) + test_approx(unit.startup_categories[1].cost[1], 7634.226) - # Check randomized costs - test_approx(unit.min_power_cost[1], 831.977) - test_approx(unit.cost_segments[1].cost[1], 36.968) - test_approx(unit.startup_categories[1].cost[1], 7634.226) + # Check randomized load share + curr_system_load = system_load(sc) + test_approx(bus.load[1] / curr_system_load[1], 0.013) - # Check randomized load share - curr_system_load = system_load(instance) - test_approx(bus.load[1] / curr_system_load[1], 0.013) - - # System load should not change - @test prev_system_load ≈ curr_system_load + # System load should not change + @test prev_system_load ≈ curr_system_load + end end @testset "load profile" begin instance = get_instance() + for sc in instance.scenarios + # Check original load profile + @test round.(system_load(sc), digits = 1)[1:8] ≈ [ + 3059.5, + 2983.2, + 2937.5, + 2953.9, + 3073.1, + 3356.4, + 4068.5, + 4018.8, + ] - # Check original load profile - @test round.(system_load(instance), digits = 1)[1:8] ≈ - [3059.5, 2983.2, 2937.5, 2953.9, 3073.1, 3356.4, 4068.5, 4018.8] - - randomize!( - instance, - XavQiuAhm2021.Randomization(), - rng = MersenneTwister(42), - ) + randomize!( + sc; + method = XavQiuAhm2021.Randomization(), + rng = MersenneTwister(42), + ) - # Check randomized load profile - @test round.(system_load(instance), digits = 1)[1:8] ≈ - [4854.7, 4849.2, 4732.7, 4848.2, 4948.4, 5231.1, 5874.8, 5934.8] + # Check randomized load profile + @test round.(system_load(sc), digits = 1)[1:8] ≈ [ + 4854.7, + 4849.2, + 4732.7, + 4848.2, + 4948.4, + 5231.1, + 5874.8, + 5934.8, + ] + end end end diff --git a/test/transform/slice_test.jl b/test/transform/slice_test.jl index 9985144..3e75339 100644 --- a/test/transform/slice_test.jl +++ b/test/transform/slice_test.jl @@ -10,29 +10,31 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip # Should update all time-dependent fields @test modified.time == 2 - @test length(modified.power_balance_penalty) == 2 - @test length(modified.reserves_by_name["r1"].amount) == 2 - for u in modified.units - @test length(u.max_power) == 2 - @test length(u.min_power) == 2 - @test length(u.must_run) == 2 - @test length(u.min_power_cost) == 2 - for s in u.cost_segments - @test length(s.mw) == 2 - @test length(s.cost) == 2 + for sc in modified.scenarios + @test length(sc.power_balance_penalty) == 2 + @test length(sc.reserves_by_name["r1"].amount) == 2 + for u in sc.units + @test length(u.max_power) == 2 + @test length(u.min_power) == 2 + @test length(u.must_run) == 2 + @test length(u.min_power_cost) == 2 + for s in u.cost_segments + @test length(s.mw) == 2 + @test length(s.cost) == 2 + end + end + for b in sc.buses + @test length(b.load) == 2 + end + for l in sc.lines + @test length(l.normal_flow_limit) == 2 + @test length(l.emergency_flow_limit) == 2 + @test length(l.flow_limit_penalty) == 2 + end + for ps in sc.price_sensitive_loads + @test length(ps.demand) == 2 + @test length(ps.revenue) == 2 end - end - for b in modified.buses - @test length(b.load) == 2 - end - for l in modified.lines - @test length(l.normal_flow_limit) == 2 - @test length(l.emergency_flow_limit) == 2 - @test length(l.flow_limit_penalty) == 2 - end - for ps in modified.price_sensitive_loads - @test length(ps.demand) == 2 - @test length(ps.revenue) == 2 end # Should be able to build model without errors optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0) diff --git a/test/usage.jl b/test/usage.jl index 4b941f4..a699c1d 100644 --- a/test/usage.jl +++ b/test/usage.jl @@ -6,8 +6,13 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON @testset "usage" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") - for line in instance.lines, t in 1:4 - line.normal_flow_limit[t] = 10.0 + for sc in instance.scenarios + sc.power_balance_penalty = [100000 for _ in 1:instance.time] + end + for sc in instance.scenarios + for line in sc.lines, t in 1:4 + line.normal_flow_limit[t] = 10.0 + end end optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0) model = UnitCommitment.build_model( diff --git a/test/validation/repair_test.jl b/test/validation/repair_test.jl index d7e5663..30a96ae 100644 --- a/test/validation/repair_test.jl +++ b/test/validation/repair_test.jl @@ -16,8 +16,10 @@ end json = parse_case14() json["Generators"]["g1"]["Production cost curve (MW)"] = [100, 150, 200] json["Generators"]["g1"]["Production cost curve (\$)"] = [10, 25, 30] - instance = UnitCommitment._from_json(json, repair = false) - @test UnitCommitment.repair!(instance) == 4 + sc = UnitCommitment._from_json(json, repair = false) + sc.name = "s1" + sc.probability = 1.0 + @test UnitCommitment.repair!(sc) == 4 end @testset "Startup limit must be greater than Pmin" begin @@ -25,15 +27,15 @@ end json["Generators"]["g1"]["Production cost curve (MW)"] = [100, 150] json["Generators"]["g1"]["Production cost curve (\$)"] = [100, 150] json["Generators"]["g1"]["Startup limit (MW)"] = 80 - instance = UnitCommitment._from_json(json, repair = false) - @test UnitCommitment.repair!(instance) == 1 + sc = UnitCommitment._from_json(json, repair = false) + @test UnitCommitment.repair!(sc) == 1 end @testset "Startup costs and delays must be increasing" begin json = parse_case14() json["Generators"]["g1"]["Startup costs (\$)"] = [300, 200, 100] json["Generators"]["g1"]["Startup delays (h)"] = [8, 4, 2] - instance = UnitCommitment._from_json(json, repair = false) - @test UnitCommitment.repair!(instance) == 4 + sc = UnitCommitment._from_json(json, repair = false) + @test UnitCommitment.repair!(sc) == 4 end end