From c95b01dadf5bad36f973c1a289cf7067b4a03e6a Mon Sep 17 00:00:00 2001 From: oyurdakul Date: Wed, 8 Feb 2023 23:46:10 -0600 Subject: [PATCH 01/11] stochastic extension w/ scenarios --- Project.toml | 1 + src/UnitCommitment.jl | 1 + src/instance/read.jl | 72 ++++++-- src/instance/structs.jl | 37 ++-- src/model/build.jl | 32 ++-- src/model/formulations/ArrCon2000/ramp.jl | 23 +-- src/model/formulations/CarArr2006/pwlcosts.jl | 15 +- .../formulations/DamKucRajAta2016/ramp.jl | 19 +- src/model/formulations/Gar1962/prod.jl | 17 +- src/model/formulations/Gar1962/pwlcosts.jl | 15 +- src/model/formulations/Gar1962/status.jl | 1 + .../formulations/KnuOstWat2018/pwlcosts.jl | 21 +-- src/model/formulations/MorLatRam2013/ramp.jl | 31 ++-- src/model/formulations/PanGua2016/ramp.jl | 15 +- src/model/formulations/WanHob2016/ramp.jl | 61 +++---- src/model/formulations/base/bus.jl | 12 +- src/model/formulations/base/line.jl | 22 +-- src/model/formulations/base/psload.jl | 10 +- src/model/formulations/base/sensitivity.jl | 2 +- src/model/formulations/base/system.jl | 100 ++++++----- src/model/formulations/base/unit.jl | 165 ++++++++++++------ .../methods/XavQiuWanThi2019/enforce.jl | 17 +- src/solution/methods/XavQiuWanThi2019/find.jl | 37 ++-- .../methods/XavQiuWanThi2019/optimize.jl | 70 ++++---- src/solution/solution.jl | 138 ++++++++------- src/validation/repair.jl | 10 +- src/validation/validate.jl | 4 +- 27 files changed, 564 insertions(+), 384 deletions(-) diff --git a/Project.toml b/Project.toml index 9c3832e..1f3fe09 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63" +Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/UnitCommitment.jl b/src/UnitCommitment.jl index 57e795b..87e45d5 100644 --- a/src/UnitCommitment.jl +++ b/src/UnitCommitment.jl @@ -4,6 +4,7 @@ module UnitCommitment +using Base: String include("instance/structs.jl") include("model/formulations/base/structs.jl") include("solution/structs.jl") diff --git a/src/instance/read.jl b/src/instance/read.jl index 1586e63..81e7723 100644 --- a/src/instance/read.jl +++ b/src/instance/read.jl @@ -7,6 +7,7 @@ using JSON using DataStructures using GZip import Base: getindex, time +using Glob const INSTANCES_URL = "https://axavier.org/UnitCommitment.jl/0.3/instances" @@ -43,6 +44,25 @@ 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 @@ -54,18 +74,33 @@ Read instance from a file. The file may be gzipped. instance = UnitCommitment.read("/path/to/input.json.gz") ``` """ -function read(path::AbstractString)::UnitCommitmentInstance +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 +end + +function read(path::AbstractString; total_number_of_scenarios = 1, scenario_index = 1)::UnitCommitmentInstance if endswith(path, ".gz") - return _read(gzopen(path)) + scenario = _read(gzopen(path),total_number_of_scenarios, scenario_index) else - return _read(open(path)) + scenario = _read(open(path), total_number_of_scenarios, scenario_index) end + instance = UnitCommitmentInstance( + time = scenario.time, + scenarios = [scenario] + ) + return instance end -function _read(file::IO)::UnitCommitmentInstance +function _read(file::IO, total_number_of_scenarios::Int, scenario_index::Int)::UnitCommitmentScenario return _from_json( JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing)), - ) + total_number_of_scenarios, + scenario_index + ) end function _read_json(path::String)::OrderedDict @@ -77,7 +112,7 @@ function _read_json(path::String)::OrderedDict return JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing)) end -function _from_json(json; repair = true) +function _from_json(json, total_number_of_scenarios::Int, scenario_index::Int; repair = true)::UnitCommitmentScenario _migrate(json) units = Unit[] buses = Bus[] @@ -101,6 +136,17 @@ function _from_json(json; repair = true) 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}() @@ -308,7 +354,9 @@ function _from_json(json; repair = true) end end - instance = UnitCommitmentInstance( + scenario = UnitCommitmentScenario( + name = scenario_name, + probability = probability, buses_by_name = Dict(b.name => b for b in buses), buses = buses, contingencies_by_name = Dict(c.name => c for c in contingencies), @@ -320,14 +368,16 @@ function _from_json(json; repair = true) price_sensitive_loads = loads, reserves = reserves, reserves_by_name = name_to_reserve, - shortfall_penalty = shortfall_penalty, - flexiramp_shortfall_penalty = flexiramp_shortfall_penalty, + # 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)) ) if repair - UnitCommitment.repair!(instance) + UnitCommitment.repair!(scenario) end - return instance + return scenario end diff --git a/src/instance/structs.jl b/src/instance/structs.jl index f07f694..56019b1 100644 --- a/src/instance/structs.jl +++ b/src/instance/structs.jl @@ -73,7 +73,9 @@ mutable struct PriceSensitiveLoad revenue::Vector{Float64} end -Base.@kwdef mutable struct UnitCommitmentInstance +Base.@kwdef mutable struct UnitCommitmentScenario + name::String + probability::Float64 buses_by_name::Dict{AbstractString,Bus} buses::Vector{Bus} contingencies_by_name::Dict{AbstractString,Contingency} @@ -85,23 +87,34 @@ Base.@kwdef mutable struct UnitCommitmentInstance price_sensitive_loads::Vector{PriceSensitiveLoad} reserves::Vector{Reserve} reserves_by_name::Dict{AbstractString,Reserve} - shortfall_penalty::Vector{Float64} - flexiramp_shortfall_penalty::Vector{Float64} - time::Int + # shortfall_penalty::Vector{Float64} + # flexiramp_shortfall_penalty::Vector{Float64} units_by_name::Dict{AbstractString,Unit} units::Vector{Unit} + time::Int + isf::Array{Float64,2} + lodf::Array{Float64,2} +end + +Base.@kwdef mutable struct UnitCommitmentInstance + time::Int + scenarios::Vector{UnitCommitmentScenario} end function Base.show(io::IO, instance::UnitCommitmentInstance) print(io, "UnitCommitmentInstance(") - print(io, "$(length(instance.units)) units, ") - print(io, "$(length(instance.buses)) buses, ") - print(io, "$(length(instance.lines)) lines, ") - print(io, "$(length(instance.contingencies)) contingencies, ") - print( - io, - "$(length(instance.price_sensitive_loads)) price sensitive loads, ", - ) + 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, ", + ) + end print(io, "$(instance.time) time steps") print(io, ")") return diff --git a/src/model/build.jl b/src/model/build.jl index 8ee7235..84b10ee 100644 --- a/src/model/build.jl +++ b/src/model/build.jl @@ -77,20 +77,28 @@ function build_model(; end model[:obj] = AffExpr() model[:instance] = instance - _setup_transmission(model, formulation.transmission) - for l in instance.lines - _add_transmission_line!(model, l, formulation.transmission) + for g in instance.scenarios[1].units + _add_unit_first_stage!(model, g, formulation) end - for b in instance.buses - _add_bus!(model, b) + 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) + end + for b in scenario.buses + _add_bus!(model, b, scenario) + end + for ps in scenario.price_sensitive_loads + _add_price_sensitive_load!(model, ps, scenario) + end + for g in scenario.units + _add_unit_second_stage!(model, g, formulation, scenario) + end + _add_system_wide_eqs!(model, scenario) end - for g in instance.units - _add_unit!(model, g, formulation) - end - for ps in instance.price_sensitive_loads - _add_price_sensitive_load!(model, ps) - end - _add_system_wide_eqs!(model) @objective(model, Min, model[:obj]) end @info @sprintf("Built model in %.2f seconds", time_model) diff --git a/src/model/formulations/ArrCon2000/ramp.jl b/src/model/formulations/ArrCon2000/ramp.jl index 080abe6..0af66b6 100644 --- a/src/model/formulations/ArrCon2000/ramp.jl +++ b/src/model/formulations/ArrCon2000/ramp.jl @@ -8,6 +8,7 @@ function _add_ramp_eqs!( formulation_prod_vars::Gar1962.ProdVars, formulation_ramping::ArrCon2000.Ramping, formulation_status_vars::Gar1962.StatusVars, + sc::UnitCommitmentScenario )::Nothing # TODO: Move upper case constants to model[:instance] RESERVES_WHEN_START_UP = true @@ -22,7 +23,7 @@ function _add_ramp_eqs!( eq_ramp_down = _init(model, :eq_ramp_down) eq_ramp_up = _init(model, :eq_ramp_up) is_initially_on = (g.initial_status > 0) - reserve = _total_reserves(model, g) + reserve = _total_reserves(model, g, sc) # Gar1962.ProdVars prod_above = model[:prod_above] @@ -37,10 +38,10 @@ function _add_ramp_eqs!( if t == 1 if is_initially_on # min power is _not_ multiplied by is_on because if !is_on, then ramp up is irrelevant - eq_ramp_up[gn, t] = @constraint( + eq_ramp_up[sc.name, gn, t] = @constraint( model, g.min_power[t] + - prod_above[gn, t] + + prod_above[sc.name, gn, t] + (RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) <= g.initial_power + RU ) @@ -48,16 +49,16 @@ function _add_ramp_eqs!( else max_prod_this_period = g.min_power[t] * is_on[gn, t] + - prod_above[gn, t] + + prod_above[sc.name, gn, t] + ( RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0 ) min_prod_last_period = - g.min_power[t-1] * is_on[gn, t-1] + prod_above[gn, t-1] + g.min_power[t-1] * is_on[gn, t-1] + prod_above[sc.name, gn, t-1] # Equation (24) in Kneuven et al. (2020) - eq_ramp_up[gn, t] = @constraint( + eq_ramp_up[sc.name, gn, t] = @constraint( model, max_prod_this_period - min_prod_last_period <= RU * is_on[gn, t-1] + SU * switch_on[gn, t] @@ -71,24 +72,24 @@ function _add_ramp_eqs!( # min_power + RD < initial_power < SD # then the generator should be able to shut down at time t = 1, # but the constraint below will force the unit to produce power - eq_ramp_down[gn, t] = @constraint( + eq_ramp_down[sc.name, gn, t] = @constraint( model, - g.initial_power - (g.min_power[t] + prod_above[gn, t]) <= RD + g.initial_power - (g.min_power[t] + prod_above[sc.name, gn, t]) <= RD ) end else max_prod_last_period = g.min_power[t-1] * is_on[gn, t-1] + - prod_above[gn, t-1] + + prod_above[sc.name, gn, t-1] + ( RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ? reserve[t-1] : 0.0 ) min_prod_this_period = - g.min_power[t] * is_on[gn, t] + prod_above[gn, t] + g.min_power[t] * is_on[gn, t] + prod_above[sc.name, gn, t] # Equation (25) in Kneuven et al. (2020) - eq_ramp_down[gn, t] = @constraint( + eq_ramp_down[sc.name, gn, t] = @constraint( model, max_prod_last_period - min_prod_this_period <= RD * is_on[gn, t] + SD * switch_off[gn, t] diff --git a/src/model/formulations/CarArr2006/pwlcosts.jl b/src/model/formulations/CarArr2006/pwlcosts.jl index 2f13e9a..2bd2de6 100644 --- a/src/model/formulations/CarArr2006/pwlcosts.jl +++ b/src/model/formulations/CarArr2006/pwlcosts.jl @@ -8,6 +8,7 @@ function _add_production_piecewise_linear_eqs!( formulation_prod_vars::Gar1962.ProdVars, formulation_pwl_costs::CarArr2006.PwlCosts, formulation_status_vars::StatusVarsFormulation, + sc::UnitCommitmentScenario )::Nothing eq_prod_above_def = _init(model, :eq_prod_above_def) eq_segprod_limit = _init(model, :eq_segprod_limit) @@ -26,28 +27,28 @@ function _add_production_piecewise_linear_eqs!( # difference between max power for segments k and k-1 so the # value of cost_segments[k].mw[t] is the max production *for # that segment* - eq_segprod_limit[gn, t, k] = @constraint( + eq_segprod_limit[sc.name, gn, t, k] = @constraint( model, - segprod[gn, t, k] <= g.cost_segments[k].mw[t] + segprod[sc.name, gn, t, k] <= g.cost_segments[k].mw[t] ) # Also add this as an explicit upper bound on segprod to make the # solver's work a bit easier - set_upper_bound(segprod[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[gn, t] = @constraint( + eq_prod_above_def[sc.name, gn, t] = @constraint( model, - prod_above[gn, t] == sum(segprod[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 # Equation (44) in Kneuven et al. (2020) add_to_expression!( model[:obj], - segprod[gn, t, k], - g.cost_segments[k].cost[t], + segprod[sc.name, gn, t, k], + sc.probability * g.cost_segments[k].cost[t], ) end end diff --git a/src/model/formulations/DamKucRajAta2016/ramp.jl b/src/model/formulations/DamKucRajAta2016/ramp.jl index cee2db9..ad73538 100644 --- a/src/model/formulations/DamKucRajAta2016/ramp.jl +++ b/src/model/formulations/DamKucRajAta2016/ramp.jl @@ -8,6 +8,7 @@ function _add_ramp_eqs!( formulation_prod_vars::Gar1962.ProdVars, formulation_ramping::DamKucRajAta2016.Ramping, formulation_status_vars::Gar1962.StatusVars, + sc::UnitCommitmentScenario )::Nothing # TODO: Move upper case constants to model[:instance] RESERVES_WHEN_START_UP = true @@ -23,7 +24,7 @@ function _add_ramp_eqs!( gn = g.name eq_str_ramp_down = _init(model, :eq_str_ramp_down) eq_str_ramp_up = _init(model, :eq_str_ramp_up) - reserve = _total_reserves(model, g) + reserve = _total_reserves(model, g, sc) # Gar1962.ProdVars prod_above = model[:prod_above] @@ -48,15 +49,15 @@ function _add_ramp_eqs!( # end max_prod_this_period = - prod_above[gn, t] + + prod_above[sc.name, gn, t] + (RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) min_prod_last_period = 0.0 if t > 1 && time_invariant - min_prod_last_period = prod_above[gn, t-1] + min_prod_last_period = prod_above[sc.name, gn, t-1] # Equation (35) in Kneuven et al. (2020) # Sparser version of (24) - eq_str_ramp_up[gn, t] = @constraint( + eq_str_ramp_up[sc.name, gn, t] = @constraint( model, max_prod_this_period - min_prod_last_period <= (SU - g.min_power[t] - RU) * switch_on[gn, t] + @@ -65,7 +66,7 @@ function _add_ramp_eqs!( elseif (t == 1 && is_initially_on) || (t > 1 && !time_invariant) if t > 1 min_prod_last_period = - prod_above[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 @@ -76,7 +77,7 @@ function _add_ramp_eqs!( # Modified version of equation (35) in Kneuven et al. (2020) # Equivalent to (24) - eq_str_ramp_up[gn, t] = @constraint( + eq_str_ramp_up[sc.name, gn, t] = @constraint( model, max_prod_this_period - min_prod_last_period <= (SU - RU) * switch_on[gn, t] + RU * is_on[gn, t] @@ -88,7 +89,7 @@ function _add_ramp_eqs!( t > 1 && (RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN) ? reserve[t-1] : 0.0 ) - min_prod_this_period = prod_above[gn, t] + min_prod_this_period = prod_above[sc.name, gn, t] on_last_period = 0.0 if t > 1 on_last_period = is_on[gn, t-1] @@ -98,7 +99,7 @@ function _add_ramp_eqs!( if t > 1 && time_invariant # Equation (36) in Kneuven et al. (2020) - eq_str_ramp_down[gn, t] = @constraint( + eq_str_ramp_down[sc.name, gn, t] = @constraint( model, max_prod_last_period - min_prod_this_period <= (SD - g.min_power[t] - RD) * switch_off[gn, t] + @@ -110,7 +111,7 @@ function _add_ramp_eqs!( # Modified version of equation (36) in Kneuven et al. (2020) # Equivalent to (25) - eq_str_ramp_down[gn, t] = @constraint( + eq_str_ramp_down[sc.name, gn, t] = @constraint( model, max_prod_last_period - min_prod_this_period <= (SD - RD) * switch_off[gn, t] + RD * on_last_period diff --git a/src/model/formulations/Gar1962/prod.jl b/src/model/formulations/Gar1962/prod.jl index fa83e79..5571ae1 100644 --- a/src/model/formulations/Gar1962/prod.jl +++ b/src/model/formulations/Gar1962/prod.jl @@ -6,14 +6,15 @@ function _add_production_vars!( model::JuMP.Model, g::Unit, formulation_prod_vars::Gar1962.ProdVars, + sc::UnitCommitmentScenario )::Nothing prod_above = _init(model, :prod_above) segprod = _init(model, :segprod) for t in 1:model[:instance].time for k in 1:length(g.cost_segments) - segprod[g.name, t, k] = @variable(model, lower_bound = 0) + segprod[sc.name, g.name, t, k] = @variable(model, lower_bound = 0) end - prod_above[g.name, t] = @variable(model, lower_bound = 0) + prod_above[sc.name, g.name, t] = @variable(model, lower_bound = 0) end return end @@ -22,16 +23,20 @@ function _add_production_limit_eqs!( model::JuMP.Model, g::Unit, formulation_prod_vars::Gar1962.ProdVars, + sc::UnitCommitmentScenario )::Nothing eq_prod_limit = _init(model, :eq_prod_limit) is_on = model[:is_on] prod_above = model[:prod_above] - reserve = _total_reserves(model, g) + reserve = _total_reserves(model, g, sc) gn = g.name for t in 1:model[:instance].time # Objective function terms for production costs # Part of (69) of Kneuven et al. (2020) as C^R_g * u_g(t) term - add_to_expression!(model[:obj], is_on[gn, t], g.min_power_cost[t]) + + ### 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) @@ -42,9 +47,9 @@ function _add_production_limit_eqs!( if power_diff < 1e-7 power_diff = 0.0 end - eq_prod_limit[gn, t] = @constraint( + eq_prod_limit[sc.name, gn, t] = @constraint( model, - prod_above[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 3ac4871..1566ac7 100644 --- a/src/model/formulations/Gar1962/pwlcosts.jl +++ b/src/model/formulations/Gar1962/pwlcosts.jl @@ -8,6 +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 )::Nothing eq_prod_above_def = _init(model, :eq_prod_above_def) eq_segprod_limit = _init(model, :eq_segprod_limit) @@ -24,9 +25,9 @@ function _add_production_piecewise_linear_eqs!( for t in 1:model[:instance].time # Definition of production # Equation (43) in Kneuven et al. (2020) - eq_prod_above_def[gn, t] = @constraint( + eq_prod_above_def[sc.name, gn, t] = @constraint( model, - prod_above[gn, t] == sum(segprod[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 @@ -37,21 +38,21 @@ function _add_production_piecewise_linear_eqs!( # difference between max power for segments k and k-1 so the # value of cost_segments[k].mw[t] is the max production *for # that segment* - eq_segprod_limit[gn, t, k] = @constraint( + eq_segprod_limit[sc.name, gn, t, k] = @constraint( model, - segprod[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[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) add_to_expression!( model[:obj], - segprod[gn, t, k], - g.cost_segments[k].cost[t], + segprod[sc.name, gn, t, k], + sc.probability * g.cost_segments[k].cost[t], ) end end diff --git a/src/model/formulations/Gar1962/status.jl b/src/model/formulations/Gar1962/status.jl index 14c055f..2a4a911 100644 --- a/src/model/formulations/Gar1962/status.jl +++ b/src/model/formulations/Gar1962/status.jl @@ -20,6 +20,7 @@ function _add_status_vars!( switch_on[g.name, t] = @variable(model, binary = true) switch_off[g.name, t] = @variable(model, binary = true) end + add_to_expression!(model[:obj], is_on[g.name, t], g.min_power_cost[t]) end return end diff --git a/src/model/formulations/KnuOstWat2018/pwlcosts.jl b/src/model/formulations/KnuOstWat2018/pwlcosts.jl index 85afa1e..ab95e05 100644 --- a/src/model/formulations/KnuOstWat2018/pwlcosts.jl +++ b/src/model/formulations/KnuOstWat2018/pwlcosts.jl @@ -8,6 +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 )::Nothing eq_prod_above_def = _init(model, :eq_prod_above_def) eq_segprod_limit_a = _init(model, :eq_segprod_limit_a) @@ -58,27 +59,27 @@ function _add_production_piecewise_linear_eqs!( if g.min_uptime > 1 # Equation (46) in Kneuven et al. (2020) - eq_segprod_limit_a[gn, t, k] = @constraint( + eq_segprod_limit_a[sc.name, gn, t, k] = @constraint( model, - segprod[gn, t, k] <= + segprod[sc.name, gn, t, k] <= g.cost_segments[k].mw[t] * is_on[gn, t] - Cv * switch_on[gn, t] - (t < T ? Cw * switch_off[gn, t+1] : 0.0) ) else # Equation (47a)/(48a) in Kneuven et al. (2020) - eq_segprod_limit_b[gn, t, k] = @constraint( + eq_segprod_limit_b[sc.name, gn, t, k] = @constraint( model, - segprod[gn, t, k] <= + segprod[sc.name, gn, t, k] <= g.cost_segments[k].mw[t] * is_on[gn, t] - Cv * switch_on[gn, t] - (t < T ? max(0, Cv - Cw) * switch_off[gn, t+1] : 0.0) ) # Equation (47b)/(48b) in Kneuven et al. (2020) - eq_segprod_limit_c[gn, t, k] = @constraint( + eq_segprod_limit_c[sc.name, gn, t, k] = @constraint( model, - segprod[gn, t, k] <= + segprod[sc.name, gn, t, k] <= g.cost_segments[k].mw[t] * is_on[gn, t] - max(0, Cw - Cv) * switch_on[gn, t] - (t < T ? Cw * switch_off[gn, t+1] : 0.0) @@ -87,22 +88,22 @@ function _add_production_piecewise_linear_eqs!( # Definition of production # Equation (43) in Kneuven et al. (2020) - eq_prod_above_def[gn, t] = @constraint( + eq_prod_above_def[sc.name, gn, t] = @constraint( model, - prod_above[gn, t] == sum(segprod[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 # Equation (44) in Kneuven et al. (2020) add_to_expression!( model[:obj], - segprod[gn, t, k], + segprod[sc.name, gn, t, k], g.cost_segments[k].cost[t], ) # Also add an explicit upper bound on segprod to make the solver's # work a bit easier - set_upper_bound(segprod[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 5eaa178..0672011 100644 --- a/src/model/formulations/MorLatRam2013/ramp.jl +++ b/src/model/formulations/MorLatRam2013/ramp.jl @@ -8,6 +8,7 @@ function _add_ramp_eqs!( formulation_prod_vars::Gar1962.ProdVars, formulation_ramping::MorLatRam2013.Ramping, formulation_status_vars::Gar1962.StatusVars, + sc::UnitCommitmentScenario )::Nothing # TODO: Move upper case constants to model[:instance] RESERVES_WHEN_START_UP = true @@ -22,7 +23,7 @@ function _add_ramp_eqs!( gn = g.name eq_ramp_down = _init(model, :eq_ramp_down) eq_ramp_up = _init(model, :eq_str_ramp_up) - reserve = _total_reserves(model, g) + reserve = _total_reserves(model, g, sc) # Gar1962.ProdVars prod_above = model[:prod_above] @@ -39,10 +40,10 @@ function _add_ramp_eqs!( # Ramp up limit if t == 1 if is_initially_on - eq_ramp_up[gn, t] = @constraint( + eq_ramp_up[sc.name, gn, t] = @constraint( model, g.min_power[t] + - prod_above[gn, t] + + prod_above[sc.name, gn, t] + (RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) <= g.initial_power + RU ) @@ -58,13 +59,13 @@ function _add_ramp_eqs!( SU = g.startup_limit max_prod_this_period = g.min_power[t] * is_on[gn, t] + - prod_above[gn, t] + + prod_above[sc.name, gn, t] + ( RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0 ) min_prod_last_period = - g.min_power[t-1] * is_on[gn, t-1] + prod_above[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 <= @@ -74,11 +75,11 @@ function _add_ramp_eqs!( # Equation (26) in Kneuven et al. (2020) # TODO: what if RU < SU? places too stringent upper bound # prod_above[gn, t] when starting up, and creates diff with (24). - eq_ramp_up[gn, t] = @constraint( + eq_ramp_up[sc.name, gn, t] = @constraint( model, - prod_above[gn, t] + + prod_above[sc.name, gn, t] + (RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) - - prod_above[gn, t-1] <= RU + prod_above[sc.name, gn, t-1] <= RU ) end end @@ -90,9 +91,9 @@ function _add_ramp_eqs!( # min_power + RD < initial_power < SD # then the generator should be able to shut down at time t = 1, # but the constraint below will force the unit to produce power - eq_ramp_down[gn, t] = @constraint( + eq_ramp_down[sc.name, gn, t] = @constraint( model, - g.initial_power - (g.min_power[t] + prod_above[gn, t]) <= RD + g.initial_power - (g.min_power[t] + prod_above[sc.name, gn, t]) <= RD ) end else @@ -102,13 +103,13 @@ function _add_ramp_eqs!( SD = g.shutdown_limit max_prod_last_period = g.min_power[t-1] * is_on[gn, t-1] + - prod_above[gn, t-1] + + prod_above[sc.name, gn, t-1] + ( RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ? reserve[t-1] : 0.0 ) min_prod_this_period = - g.min_power[t] * is_on[gn, t] + prod_above[gn, t] + g.min_power[t] * is_on[gn, t] + prod_above[sc.name, gn, t] eq_ramp_down[gn, t] = @constraint( model, max_prod_last_period - min_prod_this_period <= @@ -118,11 +119,11 @@ function _add_ramp_eqs!( # Equation (27) in Kneuven et al. (2020) # TODO: Similar to above, what to do if shutting down in time t # and RD < SD? There is a difference with (25). - eq_ramp_down[gn, t] = @constraint( + eq_ramp_down[sc.name, gn, t] = @constraint( model, - prod_above[gn, t-1] + + prod_above[sc.name, gn, t-1] + (RESERVES_WHEN_RAMP_DOWN ? reserve[t-1] : 0.0) - - prod_above[gn, t] <= RD + prod_above[sc.name, gn, t] <= RD ) end end diff --git a/src/model/formulations/PanGua2016/ramp.jl b/src/model/formulations/PanGua2016/ramp.jl index fe5d334..eb85ff5 100644 --- a/src/model/formulations/PanGua2016/ramp.jl +++ b/src/model/formulations/PanGua2016/ramp.jl @@ -8,11 +8,12 @@ function _add_ramp_eqs!( formulation_prod_vars::Gar1962.ProdVars, formulation_ramping::PanGua2016.Ramping, formulation_status_vars::Gar1962.StatusVars, + sc::UnitCommitmentScenario )::Nothing # TODO: Move upper case constants to model[:instance] RESERVES_WHEN_SHUT_DOWN = true gn = g.name - reserve = _total_reserves(model, g) + reserve = _total_reserves(model, g, sc) eq_str_prod_limit = _init(model, :eq_str_prod_limit) eq_prod_limit_ramp_up_extra_period = _init(model, :eq_prod_limit_ramp_up_extra_period) @@ -52,9 +53,9 @@ function _add_ramp_eqs!( # Generalization of (20) # Necessary that if any of the switch_on = 1 in the sum, # then switch_off[gn, t+1] = 0 - eq_str_prod_limit[gn, t] = @constraint( + eq_str_prod_limit[sc.name, gn, t] = @constraint( model, - prod_above[gn, t] + + prod_above[sc.name, gn, t] + g.min_power[t] * is_on[gn, t] + reserve[t] <= Pbar * is_on[gn, t] - @@ -67,9 +68,9 @@ 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[gn, t] = @constraint( + eq_prod_limit_ramp_up_extra_period[sc.name, gn, t] = @constraint( model, - prod_above[gn, t] + + prod_above[sc.name, gn, t] + g.min_power[t] * is_on[gn, t] + reserve[t] <= Pbar * is_on[gn, t] - sum( @@ -84,9 +85,9 @@ function _add_ramp_eqs!( if KSD > 0 KSU = min(TRU, UT - 2 - KSD, t - 1) # Equation (41) in Kneuven et al. (2020) - eq_prod_limit_shutdown_trajectory[gn, t] = @constraint( + eq_prod_limit_shutdown_trajectory[sc.name, gn, t] = @constraint( model, - prod_above[gn, t] + + prod_above[sc.name, gn, t] + g.min_power[t] * is_on[gn, t] + (RESERVES_WHEN_SHUT_DOWN ? reserve[t] : 0.0) <= Pbar * is_on[gn, t] - sum( diff --git a/src/model/formulations/WanHob2016/ramp.jl b/src/model/formulations/WanHob2016/ramp.jl index e36e4d2..cc92ae9 100644 --- a/src/model/formulations/WanHob2016/ramp.jl +++ b/src/model/formulations/WanHob2016/ramp.jl @@ -8,6 +8,7 @@ function _add_ramp_eqs!( ::Gar1962.ProdVars, ::WanHob2016.Ramping, ::Gar1962.StatusVars, + sc::UnitCommitmentScenario )::Nothing is_initially_on = (g.initial_status > 0) SU = g.startup_limit @@ -29,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 !== "flexiramp" + if r.type !== "up-frp" && r.type !== "down-frp" error( "This formulation only supports flexiramp reserves, not $(r.type)", ) @@ -38,41 +39,41 @@ function _add_ramp_eqs!( for t in 1:model[:instance].time @constraint( model, - prod_above[gn, t] + (is_on[gn, t] * minp[t]) <= mfg[rn, gn, t] + prod_above[sc.name, gn, t] + (is_on[gn, t] * minp[t]) <= mfg[sc.name, rn, gn, t] ) # Eq. (19) in Wang & Hobbs (2016) - @constraint(model, mfg[rn, gn, t] <= is_on[gn, t] * maxp[t]) # Eq. (22) in Wang & Hobbs (2016) + @constraint(model, mfg[sc.name, rn, 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[gn, t] - dwflexiramp[rn, gn, 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[gn, t] - dwflexiramp[rn, gn, t] + + prod_above[sc.name, gn, t] - dwflexiramp[sc.name, rn, gn, t] + (is_on[gn, t] * minp[t]) <= - mfg[rn, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1])) + mfg[sc.name, rn, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1])) ) # second inequality of Eq. (20) in Wang & Hobbs (2016) @constraint( model, minp[t] * (is_on[gn, t+1] + is_on[gn, t] - 1) <= - prod_above[gn, t] + - upflexiramp[rn, gn, t] + + prod_above[sc.name, gn, t] + + upflexiramp[sc.name, rn, gn, t] + (is_on[gn, t] * minp[t]) ) # first inequality of Eq. (21) in Wang & Hobbs (2016) @constraint( model, - prod_above[gn, t] + - upflexiramp[rn, gn, t] + + prod_above[sc.name, gn, t] + + upflexiramp[sc.name, rn, gn, t] + (is_on[gn, t] * minp[t]) <= - mfg[rn, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1])) + mfg[sc.name, rn, 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[rn, gn, t] <= - prod_above[gn, t-1] + + mfg[sc.name, rn, gn, t] <= + prod_above[sc.name, gn, t-1] + (is_on[gn, t-1] * minp[t]) + (RU * is_on[gn, t-1]) + (SU * (is_on[gn, t] - is_on[gn, t-1])) + @@ -80,8 +81,8 @@ function _add_ramp_eqs!( ) # Eq. (23) in Wang & Hobbs (2016) @constraint( model, - (prod_above[gn, t-1] + (is_on[gn, t-1] * minp[t])) - - (prod_above[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]) @@ -89,7 +90,7 @@ function _add_ramp_eqs!( else @constraint( model, - mfg[rn, gn, t] <= + mfg[sc.name, rn, gn, t] <= initial_power + (RU * is_initially_on) + (SU * (is_on[gn, t] - is_initially_on)) + @@ -98,7 +99,7 @@ function _add_ramp_eqs!( @constraint( model, initial_power - - (prod_above[gn, t] + (is_on[gn, t] * minp[t])) <= + (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) @@ -106,7 +107,7 @@ function _add_ramp_eqs!( end @constraint( model, - mfg[rn, gn, t] <= + mfg[sc.name, rn, 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) @@ -114,11 +115,11 @@ 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[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, - upflexiramp[rn, gn, t] <= + upflexiramp[sc.name, rn, gn, t] <= RU * is_on[gn, t] + SU * (is_on[gn, t+1] - is_on[gn, t]) + maxp[t] * (1 - is_on[gn, t+1]) @@ -126,11 +127,11 @@ 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[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, - dwflexiramp[rn, gn, t] <= + dwflexiramp[sc.name, rn, gn, t] <= RD * is_on[gn, t+1] + SD * (is_on[gn, t] - is_on[gn, t+1]) + maxp[t] * (1 - is_on[gn, t]) @@ -138,26 +139,26 @@ function _add_ramp_eqs!( @constraint( model, -maxp[t] * is_on[gn, t] + minp[t] * is_on[gn, t+1] <= - upflexiramp[rn, gn, t] + upflexiramp[sc.name, rn, gn, t] ) # first inequality of Eq. (28) in Wang & Hobbs (2016) @constraint( model, - upflexiramp[rn, gn, t] <= maxp[t] * is_on[gn, t+1] + upflexiramp[sc.name, rn, gn, t] <= maxp[t] * is_on[gn, t+1] ) # second inequality of Eq. (28) in Wang & Hobbs (2016) @constraint( model, - -maxp[t] * is_on[gn, t+1] <= dwflexiramp[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, - dwflexiramp[rn, gn, t] <= + dwflexiramp[sc.name, rn, gn, t] <= (maxp[t] * is_on[gn, t]) - (minp[t] * is_on[gn, t+1]) ) # second inequality of Eq. (29) in Wang & Hobbs (2016) else @constraint( model, - mfg[rn, gn, t] <= - prod_above[gn, t-1] + + mfg[sc.name, rn, gn, t] <= + prod_above[sc.name, gn, t-1] + (is_on[gn, t-1] * minp[t]) + (RU * is_on[gn, t-1]) + (SU * (is_on[gn, t] - is_on[gn, t-1])) + @@ -165,8 +166,8 @@ function _add_ramp_eqs!( ) # Eq. (23) in Wang & Hobbs (2016) for the last time period @constraint( model, - (prod_above[gn, t-1] + (is_on[gn, t-1] * minp[t])) - - (prod_above[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]) diff --git a/src/model/formulations/base/bus.jl b/src/model/formulations/base/bus.jl index d881fc6..926bc2c 100644 --- a/src/model/formulations/base/bus.jl +++ b/src/model/formulations/base/bus.jl @@ -2,22 +2,22 @@ # 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)::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 # Fixed load - net_injection[b.name, t] = AffExpr(-b.load[t]) + net_injection[sc.name, b.name, t] = AffExpr(-b.load[t]) # Load curtailment - curtail[b.name, t] = + curtail[sc.name, b.name, t] = @variable(model, lower_bound = 0, upper_bound = b.load[t]) - add_to_expression!(net_injection[b.name, t], curtail[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[b.name, t], - model[:instance].power_balance_penalty[t], + curtail[sc.name, b.name, t], + sc.power_balance_penalty[t] * sc.probability, ) end return diff --git a/src/model/formulations/base/line.jl b/src/model/formulations/base/line.jl index 6398c8d..0a4a586 100644 --- a/src/model/formulations/base/line.jl +++ b/src/model/formulations/base/line.jl @@ -6,14 +6,15 @@ function _add_transmission_line!( model::JuMP.Model, lm::TransmissionLine, f::ShiftFactorsFormulation, + sc::UnitCommitmentScenario )::Nothing overflow = _init(model, :overflow) for t in 1:model[:instance].time - overflow[lm.name, t] = @variable(model, lower_bound = 0) + overflow[sc.name, lm.name, t] = @variable(model, lower_bound = 0) add_to_expression!( model[:obj], - overflow[lm.name, t], - lm.flow_limit_penalty[t], + overflow[sc.name, lm.name, t], + lm.flow_limit_penalty[t] * sc.probability, ) end return @@ -22,27 +23,28 @@ end function _setup_transmission( model::JuMP.Model, formulation::ShiftFactorsFormulation, + scenario::UnitCommitmentScenario )::Nothing instance = model[:instance] isf = formulation.precomputed_isf lodf = formulation.precomputed_lodf - if length(instance.buses) == 1 + if length(scenario.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 = instance.lines, - buses = instance.buses, + lines = scenario.lines, + buses = scenario.buses, ) 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 = instance.lines, - buses = instance.buses, + lines = scenario.lines, + buses = scenario.buses, isf = isf, ) end @@ -55,7 +57,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 if r.shortfall_penalty >= 0 add_to_expression!( model[:obj], - r.shortfall_penalty, - model[:reserve_shortfall][r.name, t], + r.shortfall_penalty * sc.probability, + model[:reserve_shortfall][sc.name, r.name, t], ) end end @@ -57,7 +57,7 @@ function _add_spinning_reserve_eqs!(model::JuMP.Model)::Nothing return end -function _add_flexiramp_reserve_eqs!(model::JuMP.Model)::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 @@ -65,31 +65,41 @@ function _add_flexiramp_reserve_eqs!(model::JuMP.Model)::Nothing # objective function. eq_min_upflexiramp = _init(model, :eq_min_upflexiramp) eq_min_dwflexiramp = _init(model, :eq_min_dwflexiramp) - instance = model[:instance] - for r in instance.reserves - r.type == "flexiramp" || continue - for t in 1:instance.time - # Eq. (17) in Wang & Hobbs (2016) - eq_min_upflexiramp[r.name, t] = @constraint( - model, - sum(model[:upflexiramp][r.name, g.name, t] for g in r.units) + model[:upflexiramp_shortfall][r.name, t] >= r.amount[t] - ) - # Eq. (18) in Wang & Hobbs (2016) - eq_min_dwflexiramp[r.name, t] = @constraint( - model, - sum(model[:dwflexiramp][r.name, g.name, t] for g in r.units) + model[:dwflexiramp_shortfall][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, - ( - model[:upflexiramp_shortfall][r.name, t] + - model[:dwflexiramp_shortfall][r.name, t] - ), + 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, + 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 bcccf16..788ac0a 100644 --- a/src/model/formulations/base/unit.jl +++ b/src/model/formulations/base/unit.jl @@ -2,7 +2,7 @@ # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. # Released under the modified BSD license. See COPYING.md for more details. -function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation) +function _add_unit_first_stage!(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 @@ -11,22 +11,34 @@ function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation) 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_startup_cost_eqs!(model, g, formulation.startup_costs) + _add_status_eqs!(model, g, formulation.status_vars) + return +end + +function _add_unit_second_stage!(model::JuMP.Model, g::Unit, formulation::Formulation, + scenario::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) + + # Constraints and objective function + _add_net_injection_eqs!(model, g, scenario) + _add_production_limit_eqs!(model, g, formulation.prod_vars, scenario) _add_production_piecewise_linear_eqs!( model, g, formulation.prod_vars, formulation.pwl_costs, formulation.status_vars, + scenario ) _add_ramp_eqs!( model, @@ -34,26 +46,64 @@ function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation) formulation.prod_vars, formulation.ramping, formulation.status_vars, + scenario ) - _add_startup_cost_eqs!(model, g, formulation.startup_costs) - _add_startup_shutdown_limit_eqs!(model, g) - _add_status_eqs!(model, g, formulation.status_vars) + _add_startup_shutdown_limit_eqs!(model, g, scenario) 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)::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[r.name, g.name, t] = @variable(model, lower_bound = 0) - if (r.name, t) ∉ keys(reserve_shortfall) - reserve_shortfall[r.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) if r.shortfall_penalty < 0 - set_upper_bound(reserve_shortfall[r.name, t], 0.0) + set_upper_bound(reserve_shortfall[sc.name, r.name, t], 0.0) end end end @@ -61,27 +111,35 @@ function _add_spinning_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing return end -function _add_flexiramp_reserve_vars!(model::JuMP.Model, g::Unit)::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 - r.type == "flexiramp" || continue - for t in 1:model[:instance].time - # maximum feasible generation, \bar{g_{its}} in Wang & Hobbs (2016) - mfg[r.name, g.name, t] = @variable(model, lower_bound = 0) - upflexiramp[r.name, g.name, t] = @variable(model) # up-flexiramp, ur_{it} in Wang & Hobbs (2016) - dwflexiramp[r.name, g.name, t] = @variable(model) # down-flexiramp, dr_{it} in Wang & Hobbs (2016) - if (r.name, t) ∉ keys(upflexiramp_shortfall) - upflexiramp_shortfall[r.name, t] = - @variable(model, lower_bound = 0) - dwflexiramp_shortfall[r.name, t] = - @variable(model, lower_bound = 0) - if r.shortfall_penalty < 0 - set_upper_bound(upflexiramp_shortfall[r.name, t], 0.0) - set_upper_bound(dwflexiramp_shortfall[r.name, t], 0.0) + 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 end end end @@ -99,32 +157,32 @@ function _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing return end -function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::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] prod_above = model[:prod_above] - reserve = _total_reserves(model, g) + reserve = _total_reserves(model, g, sc) switch_off = model[:switch_off] switch_on = model[:switch_on] T = model[:instance].time for t in 1:T # Startup limit - eq_startup_limit[g.name, t] = @constraint( + eq_startup_limit[sc.name, g.name, t] = @constraint( model, - prod_above[g.name, t] + reserve[t] <= + prod_above[sc.name, g.name, t] + reserve[t] <= (g.max_power[t] - g.min_power[t]) * is_on[g.name, t] - max(0, g.max_power[t] - g.startup_limit) * switch_on[g.name, t] ) # Shutdown limit if g.initial_power > g.shutdown_limit - eq_shutdown_limit[g.name, 0] = + eq_shutdown_limit[sc.name, g.name, 0] = @constraint(model, switch_off[g.name, 1] <= 0) end if t < T - eq_shutdown_limit[g.name, t] = @constraint( + eq_shutdown_limit[sc.name, g.name, t] = @constraint( model, - prod_above[g.name, t] <= + prod_above[sc.name, g.name, t] <= (g.max_power[t] - g.min_power[t]) * is_on[g.name, t] - max(0, g.max_power[t] - g.shutdown_limit) * switch_off[g.name, t+1] @@ -138,43 +196,44 @@ function _add_ramp_eqs!( model::JuMP.Model, g::Unit, formulation::RampingFormulation, + sc::UnitCommitmentScenario )::Nothing prod_above = model[:prod_above] - reserve = _total_reserves(model, g) + reserve = _total_reserves(model, g, sc) eq_ramp_up = _init(model, :eq_ramp_up) eq_ramp_down = _init(model, :eq_ramp_down) for t in 1:model[:instance].time # Ramp up limit if t == 1 if _is_initially_on(g) == 1 - eq_ramp_up[g.name, t] = @constraint( + eq_ramp_up[sc.name, g.name, t] = @constraint( model, - prod_above[g.name, t] + reserve[t] <= + prod_above[sc.name, g.name, t] + reserve[t] <= (g.initial_power - g.min_power[t]) + g.ramp_up_limit ) end else - eq_ramp_up[g.name, t] = @constraint( + eq_ramp_up[sc.name, g.name, t] = @constraint( model, - prod_above[g.name, t] + reserve[t] <= - prod_above[g.name, t-1] + g.ramp_up_limit + prod_above[sc.name, g.name, t] + reserve[t] <= + prod_above[sc.name, g.name, t-1] + g.ramp_up_limit ) end # Ramp down limit if t == 1 if _is_initially_on(g) == 1 - eq_ramp_down[g.name, t] = @constraint( + eq_ramp_down[sc.name, g.name, t] = @constraint( model, - prod_above[g.name, t] >= + prod_above[sc.name, g.name, t] >= (g.initial_power - g.min_power[t]) - g.ramp_down_limit ) end else - eq_ramp_down[g.name, t] = @constraint( + eq_ramp_down[sc.name, g.name, t] = @constraint( model, - prod_above[g.name, t] >= - prod_above[g.name, t-1] - g.ramp_down_limit + prod_above[sc.name, g.name, t] >= + prod_above[sc.name, g.name, t-1] - g.ramp_down_limit ) end end @@ -223,30 +282,30 @@ function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing end end -function _add_net_injection_eqs!(model::JuMP.Model, g::Unit)::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 add_to_expression!( - expr_net_injection[g.bus.name, t], - model[:prod_above][g.name, t], + expr_net_injection[sc.name, g.bus.name, t], + model[:prod_above][sc.name, g.name, t], 1.0, ) add_to_expression!( - expr_net_injection[g.bus.name, t], + expr_net_injection[sc.name, g.bus.name, t], model[:is_on][g.name, t], g.min_power[t], ) end end -function _total_reserves(model, g)::Vector +function _total_reserves(model, g, sc)::Vector T = model[:instance].time reserve = [0.0 for _ in 1:T] spinning_reserves = [r for r in g.reserves if r.type == "spinning"] if !isempty(spinning_reserves) reserve += [ - sum(model[:reserve][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/methods/XavQiuWanThi2019/enforce.jl b/src/solution/methods/XavQiuWanThi2019/enforce.jl index 526274f..a4e6ce7 100644 --- a/src/solution/methods/XavQiuWanThi2019/enforce.jl +++ b/src/solution/methods/XavQiuWanThi2019/enforce.jl @@ -5,13 +5,15 @@ function _enforce_transmission( model::JuMP.Model, violations::Vector{_Violation}, + sc::UnitCommitmentScenario )::Nothing for v in violations _enforce_transmission( model = model, + sc = sc, violation = v, - isf = model[:isf], - lodf = model[:lodf], + isf = sc.isf, + lodf = sc.lodf, ) end return @@ -19,6 +21,7 @@ end function _enforce_transmission(; model::JuMP.Model, + sc::UnitCommitmentScenario, violation::_Violation, isf::Matrix{Float64}, lodf::Matrix{Float64}, @@ -51,7 +54,7 @@ function _enforce_transmission(; t = violation.time flow = @variable(model, base_name = "flow[$fm,$t]") - v = overflow[violation.monitored_line.name, violation.time] + v = overflow[sc.name, violation.monitored_line.name, violation.time] @constraint(model, flow <= limit + v) @constraint(model, -flow <= limit + v) @@ -59,23 +62,23 @@ function _enforce_transmission(; @constraint( model, flow == sum( - net_injection[b.name, violation.time] * + net_injection[sc.name, b.name, violation.time] * isf[violation.monitored_line.offset, b.offset] for - b in instance.buses if b.offset > 0 + b in sc.buses if b.offset > 0 ) ) else @constraint( model, flow == sum( - net_injection[b.name, violation.time] * ( + net_injection[sc.name, b.name, violation.time] * ( isf[violation.monitored_line.offset, b.offset] + ( lodf[ violation.monitored_line.offset, violation.outage_line.offset, ] * isf[violation.outage_line.offset, b.offset] ) - ) for b in instance.buses if b.offset > 0 + ) for b in sc.buses if b.offset > 0 ) ) end diff --git a/src/solution/methods/XavQiuWanThi2019/find.jl b/src/solution/methods/XavQiuWanThi2019/find.jl index 3cf9094..43e3ad3 100644 --- a/src/solution/methods/XavQiuWanThi2019/find.jl +++ b/src/solution/methods/XavQiuWanThi2019/find.jl @@ -5,32 +5,34 @@ import Base.Threads: @threads function _find_violations( - model::JuMP.Model; + model::JuMP.Model, + sc::UnitCommitmentScenario; max_per_line::Int, max_per_period::Int, ) instance = model[:instance] net_injection = model[:net_injection] overflow = model[:overflow] - length(instance.buses) > 1 || return [] + length(sc.buses) > 1 || return [] violations = [] @info "Verifying transmission limits..." time_screening = @elapsed begin - non_slack_buses = [b for b in instance.buses if b.offset > 0] + non_slack_buses = [b for b in sc.buses if b.offset > 0] net_injection_values = [ - value(net_injection[b.name, t]) for b in non_slack_buses, + value(net_injection[sc.name, b.name, t]) for b in non_slack_buses, t in 1:instance.time ] overflow_values = [ - value(overflow[lm.name, t]) for lm in instance.lines, + value(overflow[sc.name, lm.name, t]) for lm in sc.lines, t in 1:instance.time ] violations = UnitCommitment._find_violations( instance = instance, + sc = sc, net_injections = net_injection_values, overflow = overflow_values, - isf = model[:isf], - lodf = model[:lodf], + isf = sc.isf, + lodf = sc.lodf, max_per_line = max_per_line, max_per_period = max_per_period, ) @@ -64,15 +66,16 @@ matrix, where L is the number of transmission lines. """ function _find_violations(; instance::UnitCommitmentInstance, + sc::UnitCommitmentScenario, net_injections::Array{Float64,2}, overflow::Array{Float64,2}, 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(instance.buses) - 1 - L = length(instance.lines) + B = length(sc.buses) - 1 + L = length(sc.lines) T = instance.time K = nthreads() @@ -94,16 +97,16 @@ function _find_violations(; normal_limits::Array{Float64,2} = [ l.normal_flow_limit[t] + overflow[l.offset, t] for - l in instance.lines, t in 1:T + 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 instance.lines, t in 1:T + l in sc.lines, t in 1:T ] is_vulnerable::Array{Bool} = zeros(Bool, L) - for c in instance.contingencies + for c in sc.contingencies is_vulnerable[c.lines[1].offset] = true end @@ -111,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 @@ -144,7 +147,7 @@ function _find_violations(; filters[t], _Violation( time = t, - monitored_line = instance.lines[lm], + monitored_line = sc.lines[lm], outage_line = nothing, amount = pre_v[lm, k], ), @@ -159,8 +162,8 @@ function _find_violations(; filters[t], _Violation( time = t, - monitored_line = instance.lines[lm], - outage_line = instance.lines[lc], + monitored_line = sc.lines[lm], + outage_line = sc.lines[lc], amount = post_v[lm, lc, k], ), ) diff --git a/src/solution/methods/XavQiuWanThi2019/optimize.jl b/src/solution/methods/XavQiuWanThi2019/optimize.jl index 78c8381..c7b3a7e 100644 --- a/src/solution/methods/XavQiuWanThi2019/optimize.jl +++ b/src/solution/methods/XavQiuWanThi2019/optimize.jl @@ -10,43 +10,47 @@ 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 - initial_time = time() - large_gap = false - has_transmission = (length(model[:isf]) > 0) - if has_transmission && method.two_phase_gap - set_gap(1e-2) - large_gap = true - end - while true - time_elapsed = time() - initial_time - time_remaining = method.time_limit - time_elapsed - if time_remaining < 0 - @info "Time limit exceeded" - break + for scenario in model[:instance].scenarios + large_gap = false + has_transmission = (length(scenario.isf) > 0) + if has_transmission && method.two_phase_gap + set_gap(1e-2) + large_gap = true end - @info @sprintf( - "Setting MILP time limit to %.2f seconds", - time_remaining - ) - JuMP.set_time_limit_sec(model, time_remaining) - @info "Solving MILP..." JuMP.optimize!(model) - has_transmission || break - violations = _find_violations( - model, - max_per_line = method.max_violations_per_line, - max_per_period = method.max_violations_per_period, - ) - if isempty(violations) - @info "No violations found" - if large_gap - large_gap = false - set_gap(method.gap_limit) - else + while true + initial_time = time() + time_elapsed = time() - initial_time + time_remaining = method.time_limit - time_elapsed + if time_remaining < 0 + @info "Time limit exceeded" break end - else - _enforce_transmission(model, violations) + @info @sprintf( + "Setting MILP time limit to %.2f seconds", + time_remaining + ) + JuMP.set_time_limit_sec(model, time_remaining) + @info "Solving MILP..." + JuMP.optimize!(model) + has_transmission || break + violations = _find_violations( + model, + scenario, + max_per_line = method.max_violations_per_line, + max_per_period = method.max_violations_per_period, + ) + if isempty(violations) + @info "No violations found" + if large_gap + large_gap = false + set_gap(method.gap_limit) + else + break + end + else + _enforce_transmission(model, violations, scenario) + end end end return diff --git a/src/solution/solution.jl b/src/solution/solution.jl index 0b76d3a..f02fb5b 100644 --- a/src/solution/solution.jl +++ b/src/solution/solution.jl @@ -16,34 +16,40 @@ solution = UnitCommitment.solution(model) """ function solution(model::JuMP.Model)::OrderedDict instance, T = model[:instance], model[:instance].time - function timeseries(vars, collection) + 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 production_cost(g) + 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 + ) + end + function production_cost(g, sc) return [ value(model[:is_on][g.name, t]) * g.min_power_cost[t] + sum( Float64[ - value(model[:segprod][g.name, t, k]) * + value(model[:segprod][sc.name, g.name, t, k]) * g.cost_segments[k].cost[t] for k in 1:length(g.cost_segments) ], ) for t in 1:T ] end - function production(g) + function production(g, sc) return [ value(model[:is_on][g.name, t]) * g.min_power[t] + sum( Float64[ - value(model[:segprod][g.name, t, k]) for + value(model[:segprod][sc.name, g.name, t, k]) for k in 1:length(g.cost_segments) ], ) for t in 1:T ] end - function startup_cost(g) + function startup_cost(g, sc) S = length(g.startup_categories) return [ sum( @@ -53,66 +59,70 @@ function solution(model::JuMP.Model)::OrderedDict ] end sol = OrderedDict() - sol["Production (MW)"] = - OrderedDict(g.name => production(g) for g in instance.units) - sol["Production cost (\$)"] = - OrderedDict(g.name => production_cost(g) for g in instance.units) - sol["Startup cost (\$)"] = - OrderedDict(g.name => startup_cost(g) for g in instance.units) - sol["Is on"] = timeseries(model[:is_on], instance.units) - sol["Switch on"] = timeseries(model[:switch_on], instance.units) - sol["Switch off"] = timeseries(model[:switch_off], instance.units) - sol["Net injection (MW)"] = - timeseries(model[:net_injection], instance.buses) - sol["Load curtail (MW)"] = timeseries(model[:curtail], instance.buses) - if !isempty(instance.lines) - sol["Line overflow (MW)"] = timeseries(model[:overflow], instance.lines) - end - if !isempty(instance.price_sensitive_loads) - sol["Price-sensitive loads (MW)"] = - timeseries(model[:loads], instance.price_sensitive_loads) - end - sol["Spinning reserve (MW)"] = OrderedDict( - r.name => OrderedDict( - g.name => [ - value(model[:reserve][r.name, g.name, t]) for + for sc in instance.scenarios + sol[sc.name] = OrderedDict() + sol[sc.name]["Production (MW)"] = + OrderedDict(g.name => production(g, sc) for g in sc.units) + sol[sc.name]["Production cost (\$)"] = + 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]["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) + if !isempty(sc.lines) + sol[sc.name]["Line overflow (MW)"] = timeseries_second_stage(model[:overflow], sc.lines, 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) + 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 + ] for g in r.units + ) for r in sc.reserves if r.type == "spinning" + ) + sol[sc.name]["Spinning reserve shortfall (MW)"] = OrderedDict( + r.name => [ + value(model[:reserve_shortfall][sc.name, r.name, t]) for t in 1:instance.time - ] for g in r.units - ) for r in instance.reserves if r.type == "spinning" - ) - sol["Spinning reserve shortfall (MW)"] = OrderedDict( - r.name => [ - value(model[:reserve_shortfall][r.name, t]) for - t in 1:instance.time - ] for r in instance.reserves if r.type == "spinning" - ) - sol["Up-flexiramp (MW)"] = OrderedDict( - r.name => OrderedDict( - g.name => [ - value(model[:upflexiramp][r.name, g.name, t]) for + ] for r in sc.reserves if r.type == "spinning" + ) + 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 + ] for g in r.units + ) for r in sc.reserves if r.type == "up-frp" + ) + 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 g in r.units - ) for r in instance.reserves if r.type == "flexiramp" - ) - sol["Up-flexiramp shortfall (MW)"] = OrderedDict( - r.name => [ - value(model[:upflexiramp_shortfall][r.name, t]) for - t in 1:instance.time - ] for r in instance.reserves if r.type == "flexiramp" - ) - sol["Down-flexiramp (MW)"] = OrderedDict( - r.name => OrderedDict( - g.name => [ - value(model[:dwflexiramp][r.name, g.name, t]) for + ] for r in sc.reserves if r.type == "up-frp" + ) + 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 + ] for g in r.units + ) for r in sc.reserves if r.type == "down-frp" + ) + 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 g in r.units - ) for r in instance.reserves if r.type == "flexiramp" - ) - sol["Down-flexiramp shortfall (MW)"] = OrderedDict( - r.name => [ - value(model[:upflexiramp_shortfall][r.name, t]) for - t in 1:instance.time - ] for r in instance.reserves if r.type == "flexiramp" - ) + ] for r in sc.reserves if r.type == "down-frp" + ) + 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 end diff --git a/src/validation/repair.jl b/src/validation/repair.jl index e1dd964..50ee614 100644 --- a/src/validation/repair.jl +++ b/src/validation/repair.jl @@ -3,19 +3,19 @@ # Released under the modified BSD license. See COPYING.md for more details. """ - repair!(instance) + repair!(sc) -Verifies that the given unit commitment instance is valid and automatically +Verifies that the given unit commitment scenario is valid and automatically fixes some validation errors if possible, issuing a warning for each error found. If a validation error cannot be automatically fixed, issues an exception. Returns the number of validation errors found. """ -function repair!(instance::UnitCommitmentInstance)::Int +function repair!(sc::UnitCommitmentScenario)::Int n_errors = 0 - for g in instance.units + for g in sc.units # Startup costs and delays must be increasing for s in 2:length(g.startup_categories) @@ -38,7 +38,7 @@ function repair!(instance::UnitCommitmentInstance)::Int end end - for t in 1:instance.time + for t in 1:sc.time # Production cost curve should be convex for k in 2:length(g.cost_segments) cost = g.cost_segments[k].cost[t] diff --git a/src/validation/validate.jl b/src/validation/validate.jl index 66a580f..9677661 100644 --- a/src/validation/validate.jl +++ b/src/validation/validate.jl @@ -356,7 +356,7 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01) required, ) end - elseif r.type == "flexiramp" + elseif r.type == "up-frp" upflexiramp = sum( solution["Up-flexiramp (MW)"][r.name][g.name][t] for g in r.units @@ -374,7 +374,7 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01) ) 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 From 7e8a2ee026587756903683909156faf7efd074c9 Mon Sep 17 00:00:00 2001 From: oyurdakul Date: Wed, 22 Feb 2023 12:44:46 -0600 Subject: [PATCH 02/11] 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 From 481f5a904ce8b612d3e9b968956a3b242f6817a1 Mon Sep 17 00:00:00 2001 From: oyurdakul Date: Mon, 6 Mar 2023 17:03:34 -0600 Subject: [PATCH 03/11] read and repair scenario --- src/UnitCommitment.jl | 4 + src/instance/read.jl | 85 ++++++++++++------- src/instance/structs.jl | 8 +- src/model/build.jl | 2 +- src/transform/randomize/XavQiuAhm2021.jl | 6 +- .../transform/randomize/XavQiuAhm2021_test.jl | 6 +- 6 files changed, 67 insertions(+), 44 deletions(-) diff --git a/src/UnitCommitment.jl b/src/UnitCommitment.jl index 87e45d5..81989ee 100644 --- a/src/UnitCommitment.jl +++ b/src/UnitCommitment.jl @@ -4,6 +4,10 @@ module UnitCommitment +function hello() + @show "Hello 2" +end + using Base: String include("instance/structs.jl") include("model/formulations/base/structs.jl") diff --git a/src/instance/read.jl b/src/instance/read.jl index 935d774..7a7e947 100644 --- a/src/instance/read.jl +++ b/src/instance/read.jl @@ -57,48 +57,65 @@ instance = UnitCommitment.read("/path/to/input.json.gz") """ 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)::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") + scenarios::Vector{UnitCommitmentScenario}, + path::Vector{String} +)::Vector{UnitCommitmentScenario} + number_of_scenarios = length(scenarios) + probs = [sc.probability for sc in scenarios] + total_weight = number_of_scenarios + if Float64 in typeof.(probs) + try + total_weight = sum(probs) + catch e + if isa(e, MethodError) + error("If any of the scenarios is assigned a weight, then all scenarios must be assigned weights.") end - scenario = _repair_scenario_name_and_probability( - scenario, - scenario_path, - number_of_scenarios, - ) - push!(scenarios, scenario) end else - error("Unsupported input format") + [sc.probability = 1 for sc in scenarios] end + for (sc_path, sc) in zip(path, scenarios) + sc.name !== "" || (sc.name = first(split(last(split(sc_path, "/")), "."))) + sc.probability = (sc.probability / total_weight) + end + return scenarios +end + +function read(path::String)::UnitCommitmentInstance + scenarios = Vector{UnitCommitmentScenario}() + scenario = _read_scenario(path) + scenario.name = "s1" + scenario.probability = 1.0 + scenarios = [scenario] + instance = + UnitCommitmentInstance(time = scenario.time, scenarios = scenarios) + return instance +end + +function read(path::Vector{String})::UnitCommitmentInstance + scenarios = Vector{UnitCommitmentScenario}() + for scenario_path in path + scenario = _read_scenario(scenario_path) + push!(scenarios, scenario) + end + scenarios = _repair_scenario_name_and_probability(scenarios, path) instance = UnitCommitmentInstance(time = scenarios[1].time, scenarios = scenarios) return instance end +function _read_scenario(path::String)::UnitCommitmentScenario + if endswith(path, ".gz") + scenario = _read(gzopen(path)) + elseif endswith(path, ".json") + scenario = _read(open(path)) + else + error("Unsupported input format") + end + return scenario +end + function _read(file::IO)::UnitCommitmentScenario return _from_json( JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing)), @@ -138,8 +155,10 @@ function _from_json(json; repair = true)::UnitCommitmentScenario 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"] + probability = nothing + probability = json["Parameters"]["Scenario weight"] scenario_name = json["Parameters"]["Scenario name"] + scenario_name !== nothing || (scenario_name = "") name_to_bus = Dict{String,Bus}() name_to_line = Dict{String,TransmissionLine}() name_to_unit = Dict{String,Unit}() diff --git a/src/instance/structs.jl b/src/instance/structs.jl index 445a9c3..1760114 100644 --- a/src/instance/structs.jl +++ b/src/instance/structs.jl @@ -30,9 +30,9 @@ end mutable struct Unit name::String - bus::Bus + bus::Bus max_power::Vector{Float64} - min_power::Vector{Float64} + min_power::Vector{Float64} must_run::Vector{Bool} min_power_cost::Vector{Float64} cost_segments::Vector{CostSegment} @@ -74,8 +74,8 @@ mutable struct PriceSensitiveLoad end Base.@kwdef mutable struct UnitCommitmentScenario - name::Any - probability::Any + name::String + probability::Union{Float64,Nothing} buses_by_name::Dict{AbstractString,Bus} buses::Vector{Bus} contingencies_by_name::Dict{AbstractString,Contingency} diff --git a/src/model/build.jl b/src/model/build.jl index d78e75f..fdfc27e 100644 --- a/src/model/build.jl +++ b/src/model/build.jl @@ -81,7 +81,7 @@ function build_model(; _add_unit_commitment!(model, g, formulation) end for sc in instance.scenarios - @info "Building scenario $(sc.name) with" * + @info "Building scenario $(sc.name) with " * "probability $(sc.probability)" _setup_transmission(formulation.transmission, sc) for l in sc.lines diff --git a/src/transform/randomize/XavQiuAhm2021.jl b/src/transform/randomize/XavQiuAhm2021.jl index a8ed08e..b60cafd 100644 --- a/src/transform/randomize/XavQiuAhm2021.jl +++ b/src/transform/randomize/XavQiuAhm2021.jl @@ -200,14 +200,14 @@ function randomize!( rng = MersenneTwister(), )::Nothing for sc in instance.scenarios - randomize!(sc; method, rng) + randomize!(sc, method; rng) end return end function randomize!( - sc::UnitCommitment.UnitCommitmentScenario; - method::XavQiuAhm2021.Randomization, + sc::UnitCommitment.UnitCommitmentScenario, + method::XavQiuAhm2021.Randomization; rng = MersenneTwister(), )::Nothing if method.randomize_costs diff --git a/test/transform/randomize/XavQiuAhm2021_test.jl b/test/transform/randomize/XavQiuAhm2021_test.jl index 1717e21..1322b19 100644 --- a/test/transform/randomize/XavQiuAhm2021_test.jl +++ b/test/transform/randomize/XavQiuAhm2021_test.jl @@ -30,7 +30,7 @@ test_approx(x, y) = @test isapprox(x, y, atol = 1e-3) randomize!( sc, - method = XavQiuAhm2021.Randomization( + XavQiuAhm2021.Randomization( randomize_load_profile = false, ), rng = MersenneTwister(42), @@ -66,8 +66,8 @@ test_approx(x, y) = @test isapprox(x, y, atol = 1e-3) ] randomize!( - sc; - method = XavQiuAhm2021.Randomization(), + sc, + XavQiuAhm2021.Randomization(); rng = MersenneTwister(42), ) From ad4a754d63490044403df2667dfb7600486dd5d8 Mon Sep 17 00:00:00 2001 From: oyurdakul Date: Mon, 6 Mar 2023 17:07:54 -0600 Subject: [PATCH 04/11] read and repair scenario --- src/instance/read.jl | 13 ++++++++----- src/instance/structs.jl | 4 ++-- test/transform/randomize/XavQiuAhm2021_test.jl | 4 +--- 3 files changed, 11 insertions(+), 10 deletions(-) diff --git a/src/instance/read.jl b/src/instance/read.jl index 7a7e947..ee4de2c 100644 --- a/src/instance/read.jl +++ b/src/instance/read.jl @@ -58,17 +58,19 @@ instance = UnitCommitment.read("/path/to/input.json.gz") function _repair_scenario_name_and_probability( scenarios::Vector{UnitCommitmentScenario}, - path::Vector{String} + path::Vector{String}, )::Vector{UnitCommitmentScenario} number_of_scenarios = length(scenarios) probs = [sc.probability for sc in scenarios] total_weight = number_of_scenarios if Float64 in typeof.(probs) - try + try total_weight = sum(probs) catch e if isa(e, MethodError) - error("If any of the scenarios is assigned a weight, then all scenarios must be assigned weights.") + error( + "If any of the scenarios is assigned a weight, then all scenarios must be assigned weights.", + ) end end else @@ -76,7 +78,8 @@ function _repair_scenario_name_and_probability( end for (sc_path, sc) in zip(path, scenarios) - sc.name !== "" || (sc.name = first(split(last(split(sc_path, "/")), "."))) + sc.name !== "" || + (sc.name = first(split(last(split(sc_path, "/")), "."))) sc.probability = (sc.probability / total_weight) end return scenarios @@ -86,7 +89,7 @@ function read(path::String)::UnitCommitmentInstance scenarios = Vector{UnitCommitmentScenario}() scenario = _read_scenario(path) scenario.name = "s1" - scenario.probability = 1.0 + scenario.probability = 1.0 scenarios = [scenario] instance = UnitCommitmentInstance(time = scenario.time, scenarios = scenarios) diff --git a/src/instance/structs.jl b/src/instance/structs.jl index 1760114..9d565af 100644 --- a/src/instance/structs.jl +++ b/src/instance/structs.jl @@ -30,9 +30,9 @@ end mutable struct Unit name::String - bus::Bus + bus::Bus max_power::Vector{Float64} - min_power::Vector{Float64} + min_power::Vector{Float64} must_run::Vector{Bool} min_power_cost::Vector{Float64} cost_segments::Vector{CostSegment} diff --git a/test/transform/randomize/XavQiuAhm2021_test.jl b/test/transform/randomize/XavQiuAhm2021_test.jl index 1322b19..85b52e0 100644 --- a/test/transform/randomize/XavQiuAhm2021_test.jl +++ b/test/transform/randomize/XavQiuAhm2021_test.jl @@ -30,9 +30,7 @@ test_approx(x, y) = @test isapprox(x, y, atol = 1e-3) randomize!( sc, - XavQiuAhm2021.Randomization( - randomize_load_profile = false, - ), + XavQiuAhm2021.Randomization(randomize_load_profile = false), rng = MersenneTwister(42), ) From 31e06131345a127501b704aad7ed69d10ae2326a Mon Sep 17 00:00:00 2001 From: "Alinson S. Xavier" Date: Thu, 16 Mar 2023 10:09:01 -0500 Subject: [PATCH 05/11] Remove unused dependency & debug statements --- deps/formatter/format.jl | 1 - src/UnitCommitment.jl | 5 +---- src/instance/read.jl | 1 - test/Project.toml | 1 - 4 files changed, 1 insertion(+), 7 deletions(-) diff --git a/deps/formatter/format.jl b/deps/formatter/format.jl index ecceaf1..59ebd54 100644 --- a/deps/formatter/format.jl +++ b/deps/formatter/format.jl @@ -1,5 +1,4 @@ using JuliaFormatter -print(pwd()) format( [ "../../src", diff --git a/src/UnitCommitment.jl b/src/UnitCommitment.jl index 81989ee..b7a88e8 100644 --- a/src/UnitCommitment.jl +++ b/src/UnitCommitment.jl @@ -4,11 +4,8 @@ module UnitCommitment -function hello() - @show "Hello 2" -end - using Base: String + include("instance/structs.jl") include("model/formulations/base/structs.jl") include("solution/structs.jl") diff --git a/src/instance/read.jl b/src/instance/read.jl index ee4de2c..cfd07a3 100644 --- a/src/instance/read.jl +++ b/src/instance/read.jl @@ -7,7 +7,6 @@ using JSON using DataStructures using GZip import Base: getindex, time -using Glob const INSTANCES_URL = "https://axavier.org/UnitCommitment.jl/0.3/instances" diff --git a/test/Project.toml b/test/Project.toml index e65d966..b87293d 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -2,7 +2,6 @@ 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" From cb9334c0a3deec041e9a50969359caae69ef1c56 Mon Sep 17 00:00:00 2001 From: "Alinson S. Xavier" Date: Thu, 16 Mar 2023 10:21:31 -0500 Subject: [PATCH 06/11] Minor changes to tests --- test/instance/migrate_test.jl | 18 +- test/instance/read_test.jl | 190 ++++++------- .../methods/XavQiuWanThi19/filter_test.jl | 143 +++++----- .../methods/XavQiuWanThi19/find_test.jl | 53 ++-- .../XavQiuWanThi19/sensitivity_test.jl | 250 +++++++++--------- test/transform/initcond_test.jl | 25 +- .../transform/randomize/XavQiuAhm2021_test.jl | 94 +++---- test/transform/slice_test.jl | 48 ++-- test/usage.jl | 4 +- test/validation/repair_test.jl | 2 - 10 files changed, 398 insertions(+), 429 deletions(-) diff --git a/test/instance/migrate_test.jl b/test/instance/migrate_test.jl index ab8fb59..3fe1e09 100644 --- a/test/instance/migrate_test.jl +++ b/test/instance/migrate_test.jl @@ -6,17 +6,17 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip @testset "read v0.2" begin instance = UnitCommitment.read("$FIXTURES/ucjl-0.2.json.gz") - 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 + @test length(instance.scenarios) == 1 + sc = instance.scenarios[1] + @test length(sc.reserves_by_name["r1"].amount) == 4 + @test sc.units_by_name["g2"].reserves[1].name == "r1" end @testset "read v0.3" begin instance = UnitCommitment.read("$FIXTURES/ucjl-0.3.json.gz") - for sc in instance.scenarios - @test length(sc.units) == 6 - @test length(sc.buses) == 14 - @test length(sc.lines) == 20 - end + @test length(instance.scenarios) == 1 + sc = instance.scenarios[1] + @test length(sc.units) == 6 + @test length(sc.buses) == 14 + @test length(sc.lines) == 20 end diff --git a/test/instance/read_test.jl b/test/instance/read_test.jl index 64003a2..8f88fe1 100644 --- a/test/instance/read_test.jl +++ b/test/instance/read_test.jl @@ -6,110 +6,110 @@ 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.scenarios) == 1 + sc = instance.scenarios[1] + @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 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 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 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 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 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 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 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 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" - 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 = 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 = sc.units[2] - @test unit.name == "g2" - @test unit.must_run == [false for t in 1:4] - @test length(unit.reserves) == 1 + unit = sc.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 = 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" - @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" + @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 = 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 + 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 @testset "read_benchmark sub-hourly" begin diff --git a/test/solution/methods/XavQiuWanThi19/filter_test.jl b/test/solution/methods/XavQiuWanThi19/filter_test.jl index fa964d3..4ed69e7 100644 --- a/test/solution/methods/XavQiuWanThi19/filter_test.jl +++ b/test/solution/methods/XavQiuWanThi19/filter_test.jl @@ -7,78 +7,77 @@ import UnitCommitment: _Violation, _offer, _query @testset "_ViolationFilter" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") + sc = instance.scenarios[1] 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 = 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, + ), + ) - 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 + 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 diff --git a/test/solution/methods/XavQiuWanThi19/find_test.jl b/test/solution/methods/XavQiuWanThi19/find_test.jl index e4be83f..d3d3452 100644 --- a/test/solution/methods/XavQiuWanThi19/find_test.jl +++ b/test/solution/methods/XavQiuWanThi19/find_test.jl @@ -7,32 +7,31 @@ import UnitCommitment: _Violation, _offer, _query @testset "find_violations" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") - 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 + sc = instance.scenarios[1] + 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 diff --git a/test/solution/methods/XavQiuWanThi19/sensitivity_test.jl b/test/solution/methods/XavQiuWanThi19/sensitivity_test.jl index 8b3d93f..1ab5950 100644 --- a/test/solution/methods/XavQiuWanThi19/sensitivity_test.jl +++ b/test/solution/methods/XavQiuWanThi19/sensitivity_test.jl @@ -6,145 +6,141 @@ using UnitCommitment, Test, LinearAlgebra @testset "_susceptance_matrix" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") - 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 + sc = instance.scenarios[1] + 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 @testset "_reduced_incidence_matrix" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") - 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 + sc = instance.scenarios[1] + 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 @testset "_injection_shift_factors" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") - 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 + sc = instance.scenarios[1] + 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 @testset "_line_outage_factors" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") - 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 + sc = instance.scenarios[1] + 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 diff --git a/test/transform/initcond_test.jl b/test/transform/initcond_test.jl index d01e311..11304bb 100644 --- a/test/transform/initcond_test.jl +++ b/test/transform/initcond_test.jl @@ -8,21 +8,20 @@ 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 + sc = instance.scenarios[1] + # All units should have unknown initial conditions + for g in sc.units + @test g.initial_power === nothing + @test g.initial_status === nothing + end - # Generate initial conditions - UnitCommitment.generate_initial_conditions!(sc, optimizer) + # Generate initial conditions + UnitCommitment.generate_initial_conditions!(sc, optimizer) - # All units should now have known initial conditions - for g in sc.units - @test g.initial_power !== nothing - @test g.initial_status !== nothing - end + # All units should now have known initial conditions + for g in sc.units + @test g.initial_power !== nothing + @test g.initial_status !== nothing 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 85b52e0..226db47 100644 --- a/test/transform/randomize/XavQiuAhm2021_test.jl +++ b/test/transform/randomize/XavQiuAhm2021_test.jl @@ -9,77 +9,57 @@ using Distributions using Random using UnitCommitment, Cbc, JuMP -get_instance() = UnitCommitment.read_benchmark("matpower/case118/2017-02-01") +function get_scenario() + return UnitCommitment.read_benchmark( + "matpower/case118/2017-02-01", + ).scenarios[1] +end 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) + sc = get_scenario() + # 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 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 = sc.buses[1] + prev_system_load = system_load(sc) + test_approx(bus.load[1] / prev_system_load[1], 0.012) - randomize!( - sc, - XavQiuAhm2021.Randomization(randomize_load_profile = false), - rng = MersenneTwister(42), - ) + randomize!( + sc, + 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(sc) + test_approx(bus.load[1] / curr_system_load[1], 0.013) - # System load should not change - @test prev_system_load ≈ curr_system_load - end + # System load should not change + @test prev_system_load ≈ curr_system_load 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, - ] + sc = get_scenario() + # 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] - randomize!( - sc, - XavQiuAhm2021.Randomization(); - rng = MersenneTwister(42), - ) + randomize!(sc, XavQiuAhm2021.Randomization(); rng = MersenneTwister(42)) - # 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 + # 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 diff --git a/test/transform/slice_test.jl b/test/transform/slice_test.jl index 3e75339..173ebbb 100644 --- a/test/transform/slice_test.jl +++ b/test/transform/slice_test.jl @@ -7,35 +7,35 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip @testset "slice" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") modified = UnitCommitment.slice(instance, 1:2) + sc = modified.scenarios[1] # Should update all time-dependent fields @test modified.time == 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 + @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 + # Should be able to build model without errors optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0) model = UnitCommitment.build_model( diff --git a/test/usage.jl b/test/usage.jl index a699c1d..74e0719 100644 --- a/test/usage.jl +++ b/test/usage.jl @@ -7,9 +7,7 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON @testset "usage" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") for sc in instance.scenarios - sc.power_balance_penalty = [100000 for _ in 1:instance.time] - end - for sc in instance.scenarios + sc.power_balance_penalty = [100_000 for _ in 1:instance.time] for line in sc.lines, t in 1:4 line.normal_flow_limit[t] = 10.0 end diff --git a/test/validation/repair_test.jl b/test/validation/repair_test.jl index 30a96ae..060ca9f 100644 --- a/test/validation/repair_test.jl +++ b/test/validation/repair_test.jl @@ -17,8 +17,6 @@ end json["Generators"]["g1"]["Production cost curve (MW)"] = [100, 150, 200] json["Generators"]["g1"]["Production cost curve (\$)"] = [10, 25, 30] sc = UnitCommitment._from_json(json, repair = false) - sc.name = "s1" - sc.probability = 1.0 @test UnitCommitment.repair!(sc) == 4 end From 204c5d900f0f687b337488e2c96faa455d7e9012 Mon Sep 17 00:00:00 2001 From: "Alinson S. Xavier" Date: Thu, 16 Mar 2023 10:23:40 -0500 Subject: [PATCH 07/11] Remove unused dependency --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index 1f3fe09..9c3832e 100644 --- a/Project.toml +++ b/Project.toml @@ -9,7 +9,6 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63" -Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" From 3b6d810884b0340ac1b3db7dc79425bd36206796 Mon Sep 17 00:00:00 2001 From: "Alinson S. Xavier" Date: Thu, 16 Mar 2023 10:24:31 -0500 Subject: [PATCH 08/11] Remove duplicate format.jl file --- src/format.jl | 2 -- 1 file changed, 2 deletions(-) delete mode 100644 src/format.jl diff --git a/src/format.jl b/src/format.jl deleted file mode 100644 index 8b5e7a6..0000000 --- a/src/format.jl +++ /dev/null @@ -1,2 +0,0 @@ -using JuliaFormatter -format(["../../src", "../../test", "../../benchmark/run.jl"], verbose = true) From d8741f04a009c94d1de8dcd6fc50601c178bd953 Mon Sep 17 00:00:00 2001 From: "Alinson S. Xavier" Date: Thu, 16 Mar 2023 10:38:08 -0500 Subject: [PATCH 09/11] Minor edits to instance/read.jl --- src/instance/read.jl | 74 +++++++++++++++++++++----------------------- 1 file changed, 35 insertions(+), 39 deletions(-) diff --git a/src/instance/read.jl b/src/instance/read.jl index cfd07a3..b823dc1 100644 --- a/src/instance/read.jl +++ b/src/instance/read.jl @@ -43,47 +43,30 @@ function read_benchmark( return UnitCommitment.read(filename) end -""" - read(path::AbstractString)::UnitCommitmentInstance - -Read instance from a file. The file may be gzipped. - -# Example - -```julia -instance = UnitCommitment.read("/path/to/input.json.gz") -``` -""" - -function _repair_scenario_name_and_probability( +function _repair_scenario_names_and_probabilities!( scenarios::Vector{UnitCommitmentScenario}, path::Vector{String}, -)::Vector{UnitCommitmentScenario} - number_of_scenarios = length(scenarios) - probs = [sc.probability for sc in scenarios] - total_weight = number_of_scenarios - if Float64 in typeof.(probs) - try - total_weight = sum(probs) - catch e - if isa(e, MethodError) - error( - "If any of the scenarios is assigned a weight, then all scenarios must be assigned weights.", - ) - end - end - else - [sc.probability = 1 for sc in scenarios] - end - +)::Nothing + total_weight = sum([sc.probability for sc in scenarios]) for (sc_path, sc) in zip(path, scenarios) sc.name !== "" || (sc.name = first(split(last(split(sc_path, "/")), "."))) sc.probability = (sc.probability / total_weight) end - return scenarios + return end +""" + read(path::AbstractString)::UnitCommitmentInstance + +Read a deterministic test case from the given file. The file may be gzipped. + +# Example + +```julia +instance = UnitCommitment.read("s1.json.gz") +``` +""" function read(path::String)::UnitCommitmentInstance scenarios = Vector{UnitCommitmentScenario}() scenario = _read_scenario(path) @@ -95,13 +78,24 @@ function read(path::String)::UnitCommitmentInstance return instance end -function read(path::Vector{String})::UnitCommitmentInstance - scenarios = Vector{UnitCommitmentScenario}() - for scenario_path in path - scenario = _read_scenario(scenario_path) - push!(scenarios, scenario) +""" + read(path::Vector{String})::UnitCommitmentInstance + +Read a stochastic unit commitment instance from the given files. Each file +describes a scenario. The files may be gzipped. + +# Example + +```julia +instance = UnitCommitment.read(["s1.json.gz", "s2.json.gz"]) +``` +""" +function read(paths::Vector{String})::UnitCommitmentInstance + scenarios = UnitCommitmentScenario[] + for p in paths + push!(scenarios, _read_scenario(p)) end - scenarios = _repair_scenario_name_and_probability(scenarios, path) + _repair_scenario_names_and_probabilities!(scenarios, paths) instance = UnitCommitmentInstance(time = scenarios[1].time, scenarios = scenarios) return instance @@ -157,10 +151,12 @@ function _from_json(json; repair = true)::UnitCommitmentScenario error("Time step $time_step is not a divisor of 60") time_multiplier = 60 ÷ time_step T = time_horizon * time_multiplier - probability = nothing + probability = json["Parameters"]["Scenario weight"] + probability !== nothing || (probability = 1) scenario_name = json["Parameters"]["Scenario name"] scenario_name !== nothing || (scenario_name = "") + name_to_bus = Dict{String,Bus}() name_to_line = Dict{String,TransmissionLine}() name_to_unit = Dict{String,Unit}() From 20939dc4b7c41518e898dfc92285ef0c46c1cef4 Mon Sep 17 00:00:00 2001 From: "Alinson S. Xavier" Date: Thu, 16 Mar 2023 10:43:30 -0500 Subject: [PATCH 10/11] Minor edits to instance/structs.jl --- src/instance/structs.jl | 28 +++++++++++++--------------- test/instance/read_test.jl | 6 ++++++ 2 files changed, 19 insertions(+), 15 deletions(-) diff --git a/src/instance/structs.jl b/src/instance/structs.jl index 9d565af..0c72a23 100644 --- a/src/instance/structs.jl +++ b/src/instance/structs.jl @@ -74,24 +74,24 @@ mutable struct PriceSensitiveLoad end Base.@kwdef mutable struct UnitCommitmentScenario - name::String - probability::Union{Float64,Nothing} buses_by_name::Dict{AbstractString,Bus} buses::Vector{Bus} contingencies_by_name::Dict{AbstractString,Contingency} contingencies::Vector{Contingency} + isf::Array{Float64,2} lines_by_name::Dict{AbstractString,TransmissionLine} lines::Vector{TransmissionLine} + lodf::Array{Float64,2} + name::String power_balance_penalty::Vector{Float64} price_sensitive_loads_by_name::Dict{AbstractString,PriceSensitiveLoad} price_sensitive_loads::Vector{PriceSensitiveLoad} - reserves::Vector{Reserve} + probability::Float64 reserves_by_name::Dict{AbstractString,Reserve} + reserves::Vector{Reserve} + time::Int units_by_name::Dict{AbstractString,Unit} units::Vector{Unit} - time::Int - isf::Array{Float64,2} - lodf::Array{Float64,2} end Base.@kwdef mutable struct UnitCommitmentInstance @@ -100,16 +100,14 @@ Base.@kwdef mutable struct UnitCommitmentInstance end function Base.show(io::IO, instance::UnitCommitmentInstance) + sc = instance.scenarios[1] print(io, "UnitCommitmentInstance(") - print(io, "$(length(instance.scenarios)) scenarios: ") - 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, "$(length(instance.scenarios)) scenarios, ") + 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, ") print(io, "$(instance.time) time steps") print(io, ")") return diff --git a/test/instance/read_test.jl b/test/instance/read_test.jl index 8f88fe1..5fa4a93 100644 --- a/test/instance/read_test.jl +++ b/test/instance/read_test.jl @@ -6,6 +6,12 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip @testset "read_benchmark" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") + + @test repr(instance) == ( + "UnitCommitmentInstance(1 scenarios, 6 units, 14 buses, " * + "20 lines, 19 contingencies, 1 price sensitive loads, 4 time steps)" + ) + @test length(instance.scenarios) == 1 sc = instance.scenarios[1] @test length(sc.lines) == 20 From 414128cc0b034e20886f66eb009b6fcab545daf1 Mon Sep 17 00:00:00 2001 From: "Alinson S. Xavier" Date: Thu, 16 Mar 2023 12:03:40 -0500 Subject: [PATCH 11/11] Correct optimize!, add stochastic test case --- .../methods/XavQiuWanThi2019/enforce.jl | 6 +- src/solution/methods/XavQiuWanThi2019/find.jl | 44 ++++----- .../methods/XavQiuWanThi2019/optimize.jl | 91 ++++++++++++------- test/usage.jl | 76 ++++++++++------ 4 files changed, 131 insertions(+), 86 deletions(-) diff --git a/src/solution/methods/XavQiuWanThi2019/enforce.jl b/src/solution/methods/XavQiuWanThi2019/enforce.jl index 8c5f3b0..87b39f0 100644 --- a/src/solution/methods/XavQiuWanThi2019/enforce.jl +++ b/src/solution/methods/XavQiuWanThi2019/enforce.jl @@ -34,19 +34,21 @@ function _enforce_transmission(; if violation.outage_line === nothing limit = violation.monitored_line.normal_flow_limit[violation.time] @info @sprintf( - " %8.3f MW overflow in %-5s time %3d (pre-contingency)", + " %8.3f MW overflow in %-5s time %3d (pre-contingency, scenario %s)", violation.amount, violation.monitored_line.name, violation.time, + sc.name, ) else limit = violation.monitored_line.emergency_flow_limit[violation.time] @info @sprintf( - " %8.3f MW overflow in %-5s time %3d (outage: line %s)", + " %8.3f MW overflow in %-5s time %3d (outage: line %s, scenario %s)", violation.amount, violation.monitored_line.name, violation.time, violation.outage_line.name, + sc.name, ) end diff --git a/src/solution/methods/XavQiuWanThi2019/find.jl b/src/solution/methods/XavQiuWanThi2019/find.jl index b9d97e6..406989f 100644 --- a/src/solution/methods/XavQiuWanThi2019/find.jl +++ b/src/solution/methods/XavQiuWanThi2019/find.jl @@ -15,31 +15,25 @@ function _find_violations( overflow = model[:overflow] length(sc.buses) > 1 || return [] violations = [] - @info "Verifying transmission limits..." - 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 - ] - overflow_values = [ - value(overflow[sc.name, lm.name, t]) for lm in sc.lines, - t in 1:instance.time - ] - violations = UnitCommitment._find_violations( - instance = instance, - sc = sc, - net_injections = net_injection_values, - overflow = overflow_values, - isf = sc.isf, - lodf = sc.lodf, - max_per_line = max_per_line, - max_per_period = max_per_period, - ) - end - @info @sprintf( - "Verified transmission limits in %.2f seconds", - time_screening + + 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 + ] + overflow_values = [ + value(overflow[sc.name, lm.name, t]) for lm in sc.lines, + t in 1:instance.time + ] + violations = UnitCommitment._find_violations( + instance = instance, + sc = sc, + net_injections = net_injection_values, + overflow = overflow_values, + isf = sc.isf, + lodf = sc.lodf, + max_per_line = max_per_line, + max_per_period = max_per_period, ) return violations end diff --git a/src/solution/methods/XavQiuWanThi2019/optimize.jl b/src/solution/methods/XavQiuWanThi2019/optimize.jl index 14a34d2..57a202e 100644 --- a/src/solution/methods/XavQiuWanThi2019/optimize.jl +++ b/src/solution/methods/XavQiuWanThi2019/optimize.jl @@ -10,46 +10,73 @@ 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 + initial_time = time() + large_gap = false + has_transmission = false for sc in model[:instance].scenarios - large_gap = false - has_transmission = (length(sc.isf) > 0) + if length(sc.isf) > 0 + has_transmission = true + end if has_transmission && method.two_phase_gap set_gap(1e-2) large_gap = true end + end + while true + time_elapsed = time() - initial_time + time_remaining = method.time_limit - time_elapsed + if time_remaining < 0 + @info "Time limit exceeded" + break + end + @info @sprintf( + "Setting MILP time limit to %.2f seconds", + time_remaining + ) + JuMP.set_time_limit_sec(model, time_remaining) + @info "Solving MILP..." JuMP.optimize!(model) - while true - initial_time = time() - time_elapsed = time() - initial_time - time_remaining = method.time_limit - time_elapsed - if time_remaining < 0 - @info "Time limit exceeded" - break + + has_transmission || break + + @info "Verifying transmission limits..." + time_screening = @elapsed begin + violations = [] + for sc in model[:instance].scenarios + push!( + violations, + _find_violations( + model, + sc, + max_per_line = method.max_violations_per_line, + max_per_period = method.max_violations_per_period, + ), + ) end - @info @sprintf( - "Setting MILP time limit to %.2f seconds", - time_remaining - ) - JuMP.set_time_limit_sec(model, time_remaining) - @info "Solving MILP..." - JuMP.optimize!(model) - has_transmission || break - violations = _find_violations( - model, - sc, - max_per_line = method.max_violations_per_line, - max_per_period = method.max_violations_per_period, - ) - if isempty(violations) - @info "No violations found" - if large_gap - large_gap = false - set_gap(method.gap_limit) - else - break - end + end + @info @sprintf( + "Verified transmission limits in %.2f seconds", + time_screening + ) + + violations_found = false + for v in violations + if !isempty(v) + violations_found = true + end + end + + if violations_found + for (i, v) in enumerate(violations) + _enforce_transmission(model, v, model[:instance].scenarios[i]) + end + else + @info "No violations found" + if large_gap + large_gap = false + set_gap(method.gap_limit) else - _enforce_transmission(model, violations, sc) + break end end end diff --git a/test/usage.jl b/test/usage.jl index 74e0719..63b3f36 100644 --- a/test/usage.jl +++ b/test/usage.jl @@ -4,37 +4,59 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON -@testset "usage" begin - instance = UnitCommitment.read("$FIXTURES/case14.json.gz") +function _set_flow_limits!(instance) for sc in instance.scenarios sc.power_balance_penalty = [100_000 for _ in 1:instance.time] 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( - instance = instance, - optimizer = optimizer, - variable_names = true, - ) - @test name(model[:is_on]["g1", 1]) == "is_on[g1,1]" - - # Optimize and retrieve solution - UnitCommitment.optimize!(model) - solution = UnitCommitment.solution(model) - - # Write solution to a file - filename = tempname() - UnitCommitment.write(filename, solution) - loaded = JSON.parsefile(filename) - @test length(loaded["Is on"]) == 6 - - # Verify solution - @test UnitCommitment.validate(instance, solution) - - # Reoptimize with fixed solution - UnitCommitment.fix!(model, solution) - UnitCommitment.optimize!(model) - @test UnitCommitment.validate(instance, solution) +end + +@testset "usage" begin + @testset "deterministic" begin + instance = UnitCommitment.read("$FIXTURES/case14.json.gz") + _set_flow_limits!(instance) + optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0) + model = UnitCommitment.build_model( + instance = instance, + optimizer = optimizer, + variable_names = true, + ) + @test name(model[:is_on]["g1", 1]) == "is_on[g1,1]" + + # Optimize and retrieve solution + UnitCommitment.optimize!(model) + solution = UnitCommitment.solution(model) + + # Write solution to a file + filename = tempname() + UnitCommitment.write(filename, solution) + loaded = JSON.parsefile(filename) + @test length(loaded["Is on"]) == 6 + + # Verify solution + @test UnitCommitment.validate(instance, solution) + + # Reoptimize with fixed solution + UnitCommitment.fix!(model, solution) + UnitCommitment.optimize!(model) + @test UnitCommitment.validate(instance, solution) + end + + @testset "stochastic" begin + instance = UnitCommitment.read([ + "$FIXTURES/case14.json.gz", + "$FIXTURES/case14.json.gz", + ]) + _set_flow_limits!(instance) + @test length(instance.scenarios) == 2 + optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0) + model = UnitCommitment.build_model( + instance = instance, + optimizer = optimizer, + ) + UnitCommitment.optimize!(model) + solution = UnitCommitment.solution(model) + end end