diff --git a/src/UnitCommitment.jl b/src/UnitCommitment.jl index 484e5be..5d0bf24 100644 --- a/src/UnitCommitment.jl +++ b/src/UnitCommitment.jl @@ -4,6 +4,8 @@ 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..b823dc1 100644 --- a/src/instance/read.jl +++ b/src/instance/read.jl @@ -43,26 +43,76 @@ function read_benchmark( return UnitCommitment.read(filename) end +function _repair_scenario_names_and_probabilities!( + scenarios::Vector{UnitCommitmentScenario}, + path::Vector{String}, +)::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 +end + """ read(path::AbstractString)::UnitCommitmentInstance -Read instance from a file. The file may be gzipped. +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) + scenario.name = "s1" + scenario.probability = 1.0 + scenarios = [scenario] + instance = + UnitCommitmentInstance(time = scenario.time, scenarios = scenarios) + return instance +end + +""" + 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("/path/to/input.json.gz") +instance = UnitCommitment.read(["s1.json.gz", "s2.json.gz"]) ``` """ -function read(path::AbstractString)::UnitCommitmentInstance +function read(paths::Vector{String})::UnitCommitmentInstance + scenarios = UnitCommitmentScenario[] + for p in paths + push!(scenarios, _read_scenario(p)) + end + _repair_scenario_names_and_probabilities!(scenarios, paths) + instance = + UnitCommitmentInstance(time = scenarios[1].time, scenarios = scenarios) + return instance +end + +function _read_scenario(path::String)::UnitCommitmentScenario if endswith(path, ".gz") - return _read(gzopen(path)) + scenario = _read(gzopen(path)) + elseif endswith(path, ".json") + scenario = _read(open(path)) else - return _read(open(path)) + error("Unsupported input format") end + return scenario end -function _read(file::IO)::UnitCommitmentInstance +function _read(file::IO)::UnitCommitmentScenario return _from_json( JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing)), ) @@ -77,7 +127,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; repair = true)::UnitCommitmentScenario _migrate(json) units = Unit[] buses = Bus[] @@ -102,6 +152,11 @@ function _from_json(json; repair = true) time_multiplier = 60 ÷ time_step T = time_horizon * time_multiplier + 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}() @@ -118,15 +173,6 @@ function _from_json(json; repair = true) 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"] @@ -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,14 @@ 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, 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..0c72a23 100644 --- a/src/instance/structs.jl +++ b/src/instance/structs.jl @@ -73,35 +73,41 @@ mutable struct PriceSensitiveLoad revenue::Vector{Float64} end -Base.@kwdef mutable struct UnitCommitmentInstance +Base.@kwdef mutable struct UnitCommitmentScenario 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} - shortfall_penalty::Vector{Float64} - flexiramp_shortfall_penalty::Vector{Float64} + reserves::Vector{Reserve} time::Int units_by_name::Dict{AbstractString,Unit} units::Vector{Unit} end +Base.@kwdef mutable struct UnitCommitmentInstance + time::Int + scenarios::Vector{UnitCommitmentScenario} +end + function Base.show(io::IO, instance::UnitCommitmentInstance) + sc = instance.scenarios[1] 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, ") + 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/src/model/build.jl b/src/model/build.jl index 8ee7235..fdfc27e 100644 --- a/src/model/build.jl +++ b/src/model/build.jl @@ -77,20 +77,27 @@ 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_commitment!(model, g, formulation) end - for b in instance.buses - _add_bus!(model, b) + 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 sc.buses + _add_bus!(model, b, sc) + end + for ps in sc.price_sensitive_loads + _add_price_sensitive_load!(model, ps, sc) + end + for g in sc.units + _add_unit_dispatch!(model, g, formulation, sc) + end + _add_system_wide_eqs!(model, sc) 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..a043311 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,25 @@ 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..9b236f7 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,32 @@ 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..fa380c3 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,8 @@ 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 +78,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 +90,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 +100,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 +112,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..3e70a04 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,16 @@ 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]) # Production limit # Equation (18) in Kneuven et al. (2020) @@ -42,9 +43,10 @@ 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..62d0b5c 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,10 @@ 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 +39,25 @@ 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..0da8217 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,26 @@ 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..90afd71 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,14 @@ 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 +76,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 +92,10 @@ 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 +105,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 +121,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..62a83ba 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,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[gn, t] = @constraint( - model, - prod_above[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)) @@ -84,9 +86,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..417d2a3 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 @@ -38,41 +39,43 @@ 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, 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, 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] + - (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[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, 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, 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, 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 +83,13 @@ 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 +97,7 @@ function _add_ramp_eqs!( else @constraint( model, - mfg[rn, gn, t] <= + mfg[sc.name, gn, t] <= initial_power + (RU * is_initially_on) + (SU * (is_on[gn, t] - is_initially_on)) + @@ -97,8 +105,10 @@ function _add_ramp_eqs!( ) # Eq. (23) in Wang & Hobbs (2016) for the first time period @constraint( model, - initial_power - - (prod_above[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) @@ -106,7 +116,7 @@ function _add_ramp_eqs!( end @constraint( model, - mfg[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) @@ -114,11 +124,12 @@ 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 +137,12 @@ 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 +150,27 @@ 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, 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 +178,11 @@ 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..b28d3aa 100644 --- a/src/model/formulations/base/bus.jl +++ b/src/model/formulations/base/bus.jl @@ -2,22 +2,30 @@ # 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..98d5f2c 100644 --- a/src/model/formulations/base/line.jl +++ b/src/model/formulations/base/line.jl @@ -6,43 +6,43 @@ 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 end function _setup_transmission( - model::JuMP.Model, formulation::ShiftFactorsFormulation, + sc::UnitCommitmentScenario, )::Nothing - instance = model[:instance] isf = formulation.precomputed_isf lodf = formulation.precomputed_lodf - if length(instance.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 = instance.lines, - buses = instance.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 = instance.lines, - buses = instance.buses, + buses = sc.buses, + lines = sc.lines, isf = isf, ) end @@ -55,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 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 +71,10 @@ 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,29 +82,37 @@ 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 + T = model[:instance].time + for r in sc.reserves r.type == "flexiramp" || continue - for t in 1:instance.time + for t in 1:T # Eq. (17) in Wang & Hobbs (2016) - eq_min_upflexiramp[r.name, t] = @constraint( + eq_min_upflexiramp[sc.name, 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] + 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[r.name, t] = @constraint( + eq_min_dwflexiramp[sc.name, 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] + 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, + r.shortfall_penalty * sc.probability, ( - model[:upflexiramp_shortfall][r.name, t] + - model[:dwflexiramp_shortfall][r.name, t] + model[:upflexiramp_shortfall][sc.name, r.name, t] + + model[:dwflexiramp_shortfall][sc.name, r.name, t] ), ) end diff --git a/src/model/formulations/base/unit.jl b/src/model/formulations/base/unit.jl index bcccf16..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!(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 @@ -11,22 +17,40 @@ 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 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, 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, 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, + sc, ) _add_ramp_eqs!( model, @@ -34,26 +58,31 @@ function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation) formulation.prod_vars, formulation.ramping, formulation.status_vars, + sc, ) - _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, sc) 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 +90,37 @@ 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] = + 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[r.name, t] = + dwflexiramp_shortfall[sc.name, 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) + 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 @@ -99,32 +138,36 @@ 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 +181,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 +267,37 @@ 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/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 526274f..87b39f0 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}, @@ -31,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 @@ -51,7 +56,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 +64,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..406989f 100644 --- a/src/solution/methods/XavQiuWanThi2019/find.jl +++ b/src/solution/methods/XavQiuWanThi2019/find.jl @@ -5,39 +5,35 @@ 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] - net_injection_values = [ - value(net_injection[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, - t in 1:instance.time - ] - violations = UnitCommitment._find_violations( - instance = instance, - net_injections = net_injection_values, - overflow = overflow_values, - isf = model[:isf], - lodf = model[: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 @@ -64,6 +60,7 @@ 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}, @@ -71,8 +68,8 @@ function _find_violations(; max_per_line::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() @@ -93,17 +90,17 @@ 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 instance.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 instance.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) - for c in instance.contingencies + for c in sc.contingencies is_vulnerable[c.lines[1].offset] = true end @@ -144,7 +141,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 +156,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..57a202e 100644 --- a/src/solution/methods/XavQiuWanThi2019/optimize.jl +++ b/src/solution/methods/XavQiuWanThi2019/optimize.jl @@ -12,10 +12,15 @@ function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing 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 + has_transmission = false + for sc in model[:instance].scenarios + 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 @@ -31,13 +36,41 @@ function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing 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, + + @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 + end + @info @sprintf( + "Verified transmission limits in %.2f seconds", + time_screening ) - if isempty(violations) + + 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 @@ -45,8 +78,6 @@ function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing else break end - else - _enforce_transmission(model, violations) end end return diff --git a/src/solution/solution.jl b/src/solution/solution.jl index 0b76d3a..b170f30 100644 --- a/src/solution/solution.jl +++ b/src/solution/solution.jl @@ -16,34 +16,44 @@ solution = UnitCommitment.solution(model) """ function solution(model::JuMP.Model)::OrderedDict instance, T = model[:instance], model[:instance].time - function timeseries(vars, collection) - return OrderedDict( - b.name => [round(value(vars[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) + 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 +63,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) + 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(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(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(model[:overflow], sc.lines, sc = sc) + end + if !isempty(sc.price_sensitive_loads) + sol[sc.name]["Price-sensitive loads (MW)"] = + 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 + ] 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 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 == "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 == "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 + ] for g in r.units + ) 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 == "flexiramp" + ) end - if !isempty(instance.price_sensitive_loads) - sol["Price-sensitive loads (MW)"] = - timeseries(model[:loads], instance.price_sensitive_loads) + if length(instance.scenarios) == 1 + return first(values(sol)) + else + return sol end - sol["Spinning reserve (MW)"] = OrderedDict( - r.name => OrderedDict( - g.name => [ - value(model[:reserve][r.name, g.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 - 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 - 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" - ) - return sol 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..b60cafd 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/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..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 == "flexiramp" - 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 + upflexiramp_shortfall = + solution[sc.name]["Up-flexiramp shortfall (MW)"][r.name][t] - 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] + 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 - 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], + 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/instance/migrate_test.jl b/test/instance/migrate_test.jl index 2ea6d46..3fe1e09 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" + @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") - @test length(instance.units) == 6 - @test length(instance.buses) == 14 - @test length(instance.lines) == 20 + @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 cc9e332..5fa4a93 100644 --- a/test/instance/read_test.jl +++ b/test/instance/read_test.jl @@ -7,42 +7,49 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip @testset "read_benchmark" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") - @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 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 + @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 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[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[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.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.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.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.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" + @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 = instance.units[1] + unit = sc.units[1] @test unit.name == "g1" @test unit.bus.name == "b1" @test unit.ramp_up_limit == 1e6 @@ -69,14 +76,14 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip @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" + @test sc.units_by_name["g1"].name == "g1" - unit = instance.units[2] + 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[3] + unit = sc.units[3] @test unit.name == "g3" @test unit.bus.name == "b3" @test unit.ramp_up_limit == 70.0 @@ -98,23 +105,23 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip @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] + 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 instance.price_sensitive_loads_by_name["ps1"].name == "ps1" + @test sc.price_sensitive_loads_by_name["ps1"].name == "ps1" 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..4ed69e7 100644 --- a/test/solution/methods/XavQiuWanThi19/filter_test.jl +++ b/test/solution/methods/XavQiuWanThi19/filter_test.jl @@ -7,13 +7,13 @@ 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) - _offer( filter, _Violation( time = 1, - monitored_line = instance.lines[1], + monitored_line = sc.lines[1], outage_line = nothing, amount = 100.0, ), @@ -22,8 +22,8 @@ import UnitCommitment: _Violation, _offer, _query filter, _Violation( time = 1, - monitored_line = instance.lines[1], - outage_line = instance.lines[1], + monitored_line = sc.lines[1], + outage_line = sc.lines[1], amount = 300.0, ), ) @@ -31,8 +31,8 @@ import UnitCommitment: _Violation, _offer, _query filter, _Violation( time = 1, - monitored_line = instance.lines[1], - outage_line = instance.lines[5], + monitored_line = sc.lines[1], + outage_line = sc.lines[5], amount = 500.0, ), ) @@ -40,8 +40,8 @@ import UnitCommitment: _Violation, _offer, _query filter, _Violation( time = 1, - monitored_line = instance.lines[1], - outage_line = instance.lines[4], + monitored_line = sc.lines[1], + outage_line = sc.lines[4], amount = 400.0, ), ) @@ -49,8 +49,8 @@ import UnitCommitment: _Violation, _offer, _query filter, _Violation( time = 1, - monitored_line = instance.lines[2], - outage_line = instance.lines[1], + monitored_line = sc.lines[2], + outage_line = sc.lines[1], amount = 200.0, ), ) @@ -58,8 +58,8 @@ import UnitCommitment: _Violation, _offer, _query filter, _Violation( time = 1, - monitored_line = instance.lines[2], - outage_line = instance.lines[8], + monitored_line = sc.lines[2], + outage_line = sc.lines[8], amount = 100.0, ), ) @@ -68,14 +68,14 @@ import UnitCommitment: _Violation, _offer, _query expected = [ _Violation( time = 1, - monitored_line = instance.lines[2], - outage_line = instance.lines[1], + monitored_line = sc.lines[2], + outage_line = sc.lines[1], amount = 200.0, ), _Violation( time = 1, - monitored_line = instance.lines[1], - outage_line = instance.lines[5], + monitored_line = sc.lines[1], + outage_line = sc.lines[5], amount = 500.0, ), ] diff --git a/test/solution/methods/XavQiuWanThi19/find_test.jl b/test/solution/methods/XavQiuWanThi19/find_test.jl index 27110cb..d3d3452 100644 --- a/test/solution/methods/XavQiuWanThi19/find_test.jl +++ b/test/solution/methods/XavQiuWanThi19/find_test.jl @@ -7,23 +7,25 @@ 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 + 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 = instance.lines, - buses = instance.buses, + lines = sc.lines, + buses = sc.buses, ) lodf = UnitCommitment._line_outage_factors( - lines = instance.lines, - buses = instance.buses, + 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 instance.lines, 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, diff --git a/test/solution/methods/XavQiuWanThi19/sensitivity_test.jl b/test/solution/methods/XavQiuWanThi19/sensitivity_test.jl index 5d64a15..1ab5950 100644 --- a/test/solution/methods/XavQiuWanThi19/sensitivity_test.jl +++ b/test/solution/methods/XavQiuWanThi19/sensitivity_test.jl @@ -6,7 +6,8 @@ using UnitCommitment, Test, LinearAlgebra @testset "_susceptance_matrix" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") - actual = UnitCommitment._susceptance_matrix(instance.lines) + sc = instance.scenarios[1] + actual = UnitCommitment._susceptance_matrix(sc.lines) @test size(actual) == (20, 20) expected = Diagonal([ 29.5, @@ -35,9 +36,10 @@ end @testset "_reduced_incidence_matrix" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") + sc = instance.scenarios[1] actual = UnitCommitment._reduced_incidence_matrix( - lines = instance.lines, - buses = instance.buses, + lines = sc.lines, + buses = sc.buses, ) @test size(actual) == (20, 13) @test actual[1, 1] == -1.0 @@ -82,9 +84,10 @@ end @testset "_injection_shift_factors" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") + sc = instance.scenarios[1] actual = UnitCommitment._injection_shift_factors( - lines = instance.lines, - buses = instance.buses, + lines = sc.lines, + buses = sc.buses, ) @test size(actual) == (20, 13) @test round.(actual, digits = 2) == [ @@ -113,25 +116,26 @@ end @testset "_line_outage_factors" begin instance = UnitCommitment.read("$FIXTURES/case14.json.gz") + sc = instance.scenarios[1] isf_before = UnitCommitment._injection_shift_factors( - lines = instance.lines, - buses = instance.buses, + lines = sc.lines, + buses = sc.buses, ) lodf = UnitCommitment._line_outage_factors( - lines = instance.lines, - buses = instance.buses, + lines = sc.lines, + buses = sc.buses, isf = isf_before, ) - for contingency in instance.contingencies + 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 = instance.lines, - buses = instance.buses, + lines = sc.lines, + buses = sc.buses, ) lc.susceptance = prev_susceptance - for lm in instance.lines + for lm in sc.lines expected = isf_after[lm.offset, :] actual = isf_before[lm.offset, :] + diff --git a/test/transform/initcond_test.jl b/test/transform/initcond_test.jl index 56c9831..11304bb 100644 --- a/test/transform/initcond_test.jl +++ b/test/transform/initcond_test.jl @@ -8,18 +8,18 @@ using UnitCommitment, Cbc, JuMP # Load instance instance = UnitCommitment.read("$FIXTURES/case118-initcond.json.gz") optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0) - + sc = instance.scenarios[1] # All units should have unknown initial conditions - for g in instance.units + for g in sc.units @test g.initial_power === nothing @test g.initial_status === nothing end # Generate initial conditions - UnitCommitment.generate_initial_conditions!(instance, optimizer) + UnitCommitment.generate_initial_conditions!(sc, optimizer) # All units should now have known initial conditions - for g in instance.units + for g in sc.units @test g.initial_power !== nothing @test g.initial_status !== nothing end diff --git a/test/transform/randomize/XavQiuAhm2021_test.jl b/test/transform/randomize/XavQiuAhm2021_test.jl index 6d35829..226db47 100644 --- a/test/transform/randomize/XavQiuAhm2021_test.jl +++ b/test/transform/randomize/XavQiuAhm2021_test.jl @@ -9,30 +9,31 @@ using Distributions 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) +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() - + sc = get_scenario() # Check original costs - unit = instance.units[10] + 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 = instance.buses[1] - prev_system_load = system_load(instance) + bus = sc.buses[1] + prev_system_load = system_load(sc) test_approx(bus.load[1] / prev_system_load[1], 0.012) randomize!( - instance, - method = XavQiuAhm2021.Randomization( - randomize_load_profile = false, - ), + sc, + XavQiuAhm2021.Randomization(randomize_load_profile = false), rng = MersenneTwister(42), ) @@ -42,7 +43,7 @@ test_approx(x, y) = @test isapprox(x, y, atol = 1e-3) test_approx(unit.startup_categories[1].cost[1], 7634.226) # Check randomized load share - curr_system_load = system_load(instance) + curr_system_load = system_load(sc) test_approx(bus.load[1] / curr_system_load[1], 0.013) # System load should not change @@ -50,20 +51,15 @@ test_approx(x, y) = @test isapprox(x, y, atol = 1e-3) end @testset "load profile" begin - instance = get_instance() - + sc = get_scenario() # Check original load profile - @test round.(system_load(instance), digits = 1)[1:8] ≈ + @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!( - instance, - XavQiuAhm2021.Randomization(), - rng = MersenneTwister(42), - ) + randomize!(sc, XavQiuAhm2021.Randomization(); rng = MersenneTwister(42)) # Check randomized load profile - @test round.(system_load(instance), digits = 1)[1:8] ≈ + @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 9985144..173ebbb 100644 --- a/test/transform/slice_test.jl +++ b/test/transform/slice_test.jl @@ -7,12 +7,13 @@ 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 - @test length(modified.power_balance_penalty) == 2 - @test length(modified.reserves_by_name["r1"].amount) == 2 - for u in modified.units + @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 @@ -22,18 +23,19 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip @test length(s.cost) == 2 end end - for b in modified.buses + for b in sc.buses @test length(b.load) == 2 end - for l in modified.lines + 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 modified.price_sensitive_loads + 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 4b941f4..63b3f36 100644 --- a/test/usage.jl +++ b/test/usage.jl @@ -4,34 +4,59 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON +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 +end + @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 + @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 - 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 diff --git a/test/validation/repair_test.jl b/test/validation/repair_test.jl index d7e5663..060ca9f 100644 --- a/test/validation/repair_test.jl +++ b/test/validation/repair_test.jl @@ -16,8 +16,8 @@ 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) + @test UnitCommitment.repair!(sc) == 4 end @testset "Startup limit must be greater than Pmin" begin @@ -25,15 +25,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