From c95b01dadf5bad36f973c1a289cf7067b4a03e6a Mon Sep 17 00:00:00 2001 From: oyurdakul Date: Wed, 8 Feb 2023 23:46:10 -0600 Subject: [PATCH] stochastic extension w/ scenarios --- Project.toml | 1 + src/UnitCommitment.jl | 1 + src/instance/read.jl | 72 ++++++-- src/instance/structs.jl | 37 ++-- src/model/build.jl | 32 ++-- src/model/formulations/ArrCon2000/ramp.jl | 23 +-- src/model/formulations/CarArr2006/pwlcosts.jl | 15 +- .../formulations/DamKucRajAta2016/ramp.jl | 19 +- src/model/formulations/Gar1962/prod.jl | 17 +- src/model/formulations/Gar1962/pwlcosts.jl | 15 +- src/model/formulations/Gar1962/status.jl | 1 + .../formulations/KnuOstWat2018/pwlcosts.jl | 21 +-- src/model/formulations/MorLatRam2013/ramp.jl | 31 ++-- src/model/formulations/PanGua2016/ramp.jl | 15 +- src/model/formulations/WanHob2016/ramp.jl | 61 +++---- src/model/formulations/base/bus.jl | 12 +- src/model/formulations/base/line.jl | 22 +-- src/model/formulations/base/psload.jl | 10 +- src/model/formulations/base/sensitivity.jl | 2 +- src/model/formulations/base/system.jl | 100 ++++++----- src/model/formulations/base/unit.jl | 165 ++++++++++++------ .../methods/XavQiuWanThi2019/enforce.jl | 17 +- src/solution/methods/XavQiuWanThi2019/find.jl | 37 ++-- .../methods/XavQiuWanThi2019/optimize.jl | 70 ++++---- src/solution/solution.jl | 138 ++++++++------- src/validation/repair.jl | 10 +- src/validation/validate.jl | 4 +- 27 files changed, 564 insertions(+), 384 deletions(-) diff --git a/Project.toml b/Project.toml index 9c3832e..1f3fe09 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63" +Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/UnitCommitment.jl b/src/UnitCommitment.jl index 57e795b..87e45d5 100644 --- a/src/UnitCommitment.jl +++ b/src/UnitCommitment.jl @@ -4,6 +4,7 @@ module UnitCommitment +using Base: String include("instance/structs.jl") include("model/formulations/base/structs.jl") include("solution/structs.jl") diff --git a/src/instance/read.jl b/src/instance/read.jl index 1586e63..81e7723 100644 --- a/src/instance/read.jl +++ b/src/instance/read.jl @@ -7,6 +7,7 @@ using JSON using DataStructures using GZip import Base: getindex, time +using Glob const INSTANCES_URL = "https://axavier.org/UnitCommitment.jl/0.3/instances" @@ -43,6 +44,25 @@ function read_benchmark( return UnitCommitment.read(filename) end +function read_scenarios(path::AbstractString)::UnitCommitmentInstance + scenario_paths = glob("*.json", path) + scenarios = Vector{UnitCommitmentScenario}() + total_number_of_scenarios = length(scenario_paths) + for (scenario_index, scenario_path) in enumerate(scenario_paths) + scenario = read_scenario(scenario_path, + total_number_of_scenarios = total_number_of_scenarios, + scenario_index = scenario_index) + push!(scenarios, scenario) + end + instance = UnitCommitmentInstance( + time = scenarios[1].time, + scenarios = scenarios + ) + abs(sum(scenario.probability for scenario in instance.scenarios) - 1.0) <= 0.01 || + error("scenario probabilities do not add up to one") + return instance +end + """ read(path::AbstractString)::UnitCommitmentInstance @@ -54,18 +74,33 @@ Read instance from a file. The file may be gzipped. instance = UnitCommitment.read("/path/to/input.json.gz") ``` """ -function read(path::AbstractString)::UnitCommitmentInstance +function read_scenario(path::AbstractString; total_number_of_scenarios = 1, scenario_index = 1)::UnitCommitmentScenario + if endswith(path, ".gz") + return _read(gzopen(path),total_number_of_scenarios, scenario_index) + else + return _read(open(path), total_number_of_scenarios, scenario_index) + end +end + +function read(path::AbstractString; total_number_of_scenarios = 1, scenario_index = 1)::UnitCommitmentInstance if endswith(path, ".gz") - return _read(gzopen(path)) + scenario = _read(gzopen(path),total_number_of_scenarios, scenario_index) else - return _read(open(path)) + scenario = _read(open(path), total_number_of_scenarios, scenario_index) end + instance = UnitCommitmentInstance( + time = scenario.time, + scenarios = [scenario] + ) + return instance end -function _read(file::IO)::UnitCommitmentInstance +function _read(file::IO, total_number_of_scenarios::Int, scenario_index::Int)::UnitCommitmentScenario return _from_json( JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing)), - ) + total_number_of_scenarios, + scenario_index + ) end function _read_json(path::String)::OrderedDict @@ -77,7 +112,7 @@ function _read_json(path::String)::OrderedDict return JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing)) end -function _from_json(json; repair = true) +function _from_json(json, total_number_of_scenarios::Int, scenario_index::Int; repair = true)::UnitCommitmentScenario _migrate(json) units = Unit[] buses = Bus[] @@ -101,6 +136,17 @@ function _from_json(json; repair = true) error("Time step $time_step is not a divisor of 60") time_multiplier = 60 ÷ time_step T = time_horizon * time_multiplier + ##### + probability = json["Parameters"]["Scenario probability"] + if probability === nothing + probability = (1 / total_number_of_scenarios) + end + scenario_name = json["Parameters"]["Scenario name"] + if scenario_name === nothing + scenario_name = "s$(scenario_index)" + end + + ###### name_to_bus = Dict{String,Bus}() name_to_line = Dict{String,TransmissionLine}() @@ -308,7 +354,9 @@ function _from_json(json; repair = true) end end - instance = UnitCommitmentInstance( + scenario = UnitCommitmentScenario( + name = scenario_name, + probability = probability, buses_by_name = Dict(b.name => b for b in buses), buses = buses, contingencies_by_name = Dict(c.name => c for c in contingencies), @@ -320,14 +368,16 @@ function _from_json(json; repair = true) price_sensitive_loads = loads, reserves = reserves, reserves_by_name = name_to_reserve, - shortfall_penalty = shortfall_penalty, - flexiramp_shortfall_penalty = flexiramp_shortfall_penalty, + # shortfall_penalty = shortfall_penalty, + # flexiramp_shortfall_penalty = flexiramp_shortfall_penalty, time = T, units_by_name = Dict(g.name => g for g in units), units = units, + isf = spzeros(Float64, length(lines), length(buses) - 1), + lodf = spzeros(Float64, length(lines), length(lines)) ) if repair - UnitCommitment.repair!(instance) + UnitCommitment.repair!(scenario) end - return instance + return scenario end diff --git a/src/instance/structs.jl b/src/instance/structs.jl index f07f694..56019b1 100644 --- a/src/instance/structs.jl +++ b/src/instance/structs.jl @@ -73,7 +73,9 @@ mutable struct PriceSensitiveLoad revenue::Vector{Float64} end -Base.@kwdef mutable struct UnitCommitmentInstance +Base.@kwdef mutable struct UnitCommitmentScenario + name::String + probability::Float64 buses_by_name::Dict{AbstractString,Bus} buses::Vector{Bus} contingencies_by_name::Dict{AbstractString,Contingency} @@ -85,23 +87,34 @@ Base.@kwdef mutable struct UnitCommitmentInstance price_sensitive_loads::Vector{PriceSensitiveLoad} reserves::Vector{Reserve} reserves_by_name::Dict{AbstractString,Reserve} - shortfall_penalty::Vector{Float64} - flexiramp_shortfall_penalty::Vector{Float64} - time::Int + # shortfall_penalty::Vector{Float64} + # flexiramp_shortfall_penalty::Vector{Float64} units_by_name::Dict{AbstractString,Unit} units::Vector{Unit} + time::Int + isf::Array{Float64,2} + lodf::Array{Float64,2} +end + +Base.@kwdef mutable struct UnitCommitmentInstance + time::Int + scenarios::Vector{UnitCommitmentScenario} end function Base.show(io::IO, instance::UnitCommitmentInstance) print(io, "UnitCommitmentInstance(") - print(io, "$(length(instance.units)) units, ") - print(io, "$(length(instance.buses)) buses, ") - print(io, "$(length(instance.lines)) lines, ") - print(io, "$(length(instance.contingencies)) contingencies, ") - print( - io, - "$(length(instance.price_sensitive_loads)) price sensitive loads, ", - ) + print(io, "$(length(instance.scenarios)) scenarios: ") + for scenario in instance.scenarios + print(io, "Scenario $(scenario.name): ") + print(io, "$(length(scenario.units)) units, ") + print(io, "$(length(scenario.buses)) buses, ") + print(io, "$(length(scenario.lines)) lines, ") + print(io, "$(length(scenario.contingencies)) contingencies, ") + print( + io, + "$(length(scenario.price_sensitive_loads)) price sensitive loads, ", + ) + end print(io, "$(instance.time) time steps") print(io, ")") return diff --git a/src/model/build.jl b/src/model/build.jl index 8ee7235..84b10ee 100644 --- a/src/model/build.jl +++ b/src/model/build.jl @@ -77,20 +77,28 @@ function build_model(; end model[:obj] = AffExpr() model[:instance] = instance - _setup_transmission(model, formulation.transmission) - for l in instance.lines - _add_transmission_line!(model, l, formulation.transmission) + for g in instance.scenarios[1].units + _add_unit_first_stage!(model, g, formulation) end - for b in instance.buses - _add_bus!(model, b) + for scenario in instance.scenarios + @info "Building scenario $(scenario.name) with" * + "probability $(scenario.probability)" + _setup_transmission(model, formulation.transmission, scenario) + for l in scenario.lines + _add_transmission_line!(model, l, formulation.transmission, + scenario) + end + for b in scenario.buses + _add_bus!(model, b, scenario) + end + for ps in scenario.price_sensitive_loads + _add_price_sensitive_load!(model, ps, scenario) + end + for g in scenario.units + _add_unit_second_stage!(model, g, formulation, scenario) + end + _add_system_wide_eqs!(model, scenario) end - for g in instance.units - _add_unit!(model, g, formulation) - end - for ps in instance.price_sensitive_loads - _add_price_sensitive_load!(model, ps) - end - _add_system_wide_eqs!(model) @objective(model, Min, model[:obj]) end @info @sprintf("Built model in %.2f seconds", time_model) diff --git a/src/model/formulations/ArrCon2000/ramp.jl b/src/model/formulations/ArrCon2000/ramp.jl index 080abe6..0af66b6 100644 --- a/src/model/formulations/ArrCon2000/ramp.jl +++ b/src/model/formulations/ArrCon2000/ramp.jl @@ -8,6 +8,7 @@ function _add_ramp_eqs!( formulation_prod_vars::Gar1962.ProdVars, formulation_ramping::ArrCon2000.Ramping, formulation_status_vars::Gar1962.StatusVars, + sc::UnitCommitmentScenario )::Nothing # TODO: Move upper case constants to model[:instance] RESERVES_WHEN_START_UP = true @@ -22,7 +23,7 @@ function _add_ramp_eqs!( eq_ramp_down = _init(model, :eq_ramp_down) eq_ramp_up = _init(model, :eq_ramp_up) is_initially_on = (g.initial_status > 0) - reserve = _total_reserves(model, g) + reserve = _total_reserves(model, g, sc) # Gar1962.ProdVars prod_above = model[:prod_above] @@ -37,10 +38,10 @@ function _add_ramp_eqs!( if t == 1 if is_initially_on # min power is _not_ multiplied by is_on because if !is_on, then ramp up is irrelevant - eq_ramp_up[gn, t] = @constraint( + eq_ramp_up[sc.name, gn, t] = @constraint( model, g.min_power[t] + - prod_above[gn, t] + + prod_above[sc.name, gn, t] + (RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) <= g.initial_power + RU ) @@ -48,16 +49,16 @@ function _add_ramp_eqs!( else max_prod_this_period = g.min_power[t] * is_on[gn, t] + - prod_above[gn, t] + + prod_above[sc.name, gn, t] + ( RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0 ) min_prod_last_period = - g.min_power[t-1] * is_on[gn, t-1] + prod_above[gn, t-1] + g.min_power[t-1] * is_on[gn, t-1] + prod_above[sc.name, gn, t-1] # Equation (24) in Kneuven et al. (2020) - eq_ramp_up[gn, t] = @constraint( + eq_ramp_up[sc.name, gn, t] = @constraint( model, max_prod_this_period - min_prod_last_period <= RU * is_on[gn, t-1] + SU * switch_on[gn, t] @@ -71,24 +72,24 @@ function _add_ramp_eqs!( # min_power + RD < initial_power < SD # then the generator should be able to shut down at time t = 1, # but the constraint below will force the unit to produce power - eq_ramp_down[gn, t] = @constraint( + eq_ramp_down[sc.name, gn, t] = @constraint( model, - g.initial_power - (g.min_power[t] + prod_above[gn, t]) <= RD + g.initial_power - (g.min_power[t] + prod_above[sc.name, gn, t]) <= RD ) end else max_prod_last_period = g.min_power[t-1] * is_on[gn, t-1] + - prod_above[gn, t-1] + + prod_above[sc.name, gn, t-1] + ( RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ? reserve[t-1] : 0.0 ) min_prod_this_period = - g.min_power[t] * is_on[gn, t] + prod_above[gn, t] + g.min_power[t] * is_on[gn, t] + prod_above[sc.name, gn, t] # Equation (25) in Kneuven et al. (2020) - eq_ramp_down[gn, t] = @constraint( + eq_ramp_down[sc.name, gn, t] = @constraint( model, max_prod_last_period - min_prod_this_period <= RD * is_on[gn, t] + SD * switch_off[gn, t] diff --git a/src/model/formulations/CarArr2006/pwlcosts.jl b/src/model/formulations/CarArr2006/pwlcosts.jl index 2f13e9a..2bd2de6 100644 --- a/src/model/formulations/CarArr2006/pwlcosts.jl +++ b/src/model/formulations/CarArr2006/pwlcosts.jl @@ -8,6 +8,7 @@ function _add_production_piecewise_linear_eqs!( formulation_prod_vars::Gar1962.ProdVars, formulation_pwl_costs::CarArr2006.PwlCosts, formulation_status_vars::StatusVarsFormulation, + sc::UnitCommitmentScenario )::Nothing eq_prod_above_def = _init(model, :eq_prod_above_def) eq_segprod_limit = _init(model, :eq_segprod_limit) @@ -26,28 +27,28 @@ function _add_production_piecewise_linear_eqs!( # difference between max power for segments k and k-1 so the # value of cost_segments[k].mw[t] is the max production *for # that segment* - eq_segprod_limit[gn, t, k] = @constraint( + eq_segprod_limit[sc.name, gn, t, k] = @constraint( model, - segprod[gn, t, k] <= g.cost_segments[k].mw[t] + segprod[sc.name, gn, t, k] <= g.cost_segments[k].mw[t] ) # Also add this as an explicit upper bound on segprod to make the # solver's work a bit easier - set_upper_bound(segprod[gn, t, k], g.cost_segments[k].mw[t]) + set_upper_bound(segprod[sc.name, gn, t, k], g.cost_segments[k].mw[t]) # Definition of production # Equation (43) in Kneuven et al. (2020) - eq_prod_above_def[gn, t] = @constraint( + eq_prod_above_def[sc.name, gn, t] = @constraint( model, - prod_above[gn, t] == sum(segprod[gn, t, k] for k in 1:K) + prod_above[sc.name, gn, t] == sum(segprod[sc.name, gn, t, k] for k in 1:K) ) # Objective function # Equation (44) in Kneuven et al. (2020) add_to_expression!( model[:obj], - segprod[gn, t, k], - g.cost_segments[k].cost[t], + segprod[sc.name, gn, t, k], + sc.probability * g.cost_segments[k].cost[t], ) end end diff --git a/src/model/formulations/DamKucRajAta2016/ramp.jl b/src/model/formulations/DamKucRajAta2016/ramp.jl index cee2db9..ad73538 100644 --- a/src/model/formulations/DamKucRajAta2016/ramp.jl +++ b/src/model/formulations/DamKucRajAta2016/ramp.jl @@ -8,6 +8,7 @@ function _add_ramp_eqs!( formulation_prod_vars::Gar1962.ProdVars, formulation_ramping::DamKucRajAta2016.Ramping, formulation_status_vars::Gar1962.StatusVars, + sc::UnitCommitmentScenario )::Nothing # TODO: Move upper case constants to model[:instance] RESERVES_WHEN_START_UP = true @@ -23,7 +24,7 @@ function _add_ramp_eqs!( gn = g.name eq_str_ramp_down = _init(model, :eq_str_ramp_down) eq_str_ramp_up = _init(model, :eq_str_ramp_up) - reserve = _total_reserves(model, g) + reserve = _total_reserves(model, g, sc) # Gar1962.ProdVars prod_above = model[:prod_above] @@ -48,15 +49,15 @@ function _add_ramp_eqs!( # end max_prod_this_period = - prod_above[gn, t] + + prod_above[sc.name, gn, t] + (RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) min_prod_last_period = 0.0 if t > 1 && time_invariant - min_prod_last_period = prod_above[gn, t-1] + min_prod_last_period = prod_above[sc.name, gn, t-1] # Equation (35) in Kneuven et al. (2020) # Sparser version of (24) - eq_str_ramp_up[gn, t] = @constraint( + eq_str_ramp_up[sc.name, gn, t] = @constraint( model, max_prod_this_period - min_prod_last_period <= (SU - g.min_power[t] - RU) * switch_on[gn, t] + @@ -65,7 +66,7 @@ function _add_ramp_eqs!( elseif (t == 1 && is_initially_on) || (t > 1 && !time_invariant) if t > 1 min_prod_last_period = - prod_above[gn, t-1] + g.min_power[t-1] * is_on[gn, t-1] + prod_above[sc.name, gn, t-1] + g.min_power[t-1] * is_on[gn, t-1] else min_prod_last_period = max(g.initial_power, 0.0) end @@ -76,7 +77,7 @@ function _add_ramp_eqs!( # Modified version of equation (35) in Kneuven et al. (2020) # Equivalent to (24) - eq_str_ramp_up[gn, t] = @constraint( + eq_str_ramp_up[sc.name, gn, t] = @constraint( model, max_prod_this_period - min_prod_last_period <= (SU - RU) * switch_on[gn, t] + RU * is_on[gn, t] @@ -88,7 +89,7 @@ function _add_ramp_eqs!( t > 1 && (RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN) ? reserve[t-1] : 0.0 ) - min_prod_this_period = prod_above[gn, t] + min_prod_this_period = prod_above[sc.name, gn, t] on_last_period = 0.0 if t > 1 on_last_period = is_on[gn, t-1] @@ -98,7 +99,7 @@ function _add_ramp_eqs!( if t > 1 && time_invariant # Equation (36) in Kneuven et al. (2020) - eq_str_ramp_down[gn, t] = @constraint( + eq_str_ramp_down[sc.name, gn, t] = @constraint( model, max_prod_last_period - min_prod_this_period <= (SD - g.min_power[t] - RD) * switch_off[gn, t] + @@ -110,7 +111,7 @@ function _add_ramp_eqs!( # Modified version of equation (36) in Kneuven et al. (2020) # Equivalent to (25) - eq_str_ramp_down[gn, t] = @constraint( + eq_str_ramp_down[sc.name, gn, t] = @constraint( model, max_prod_last_period - min_prod_this_period <= (SD - RD) * switch_off[gn, t] + RD * on_last_period diff --git a/src/model/formulations/Gar1962/prod.jl b/src/model/formulations/Gar1962/prod.jl index fa83e79..5571ae1 100644 --- a/src/model/formulations/Gar1962/prod.jl +++ b/src/model/formulations/Gar1962/prod.jl @@ -6,14 +6,15 @@ function _add_production_vars!( model::JuMP.Model, g::Unit, formulation_prod_vars::Gar1962.ProdVars, + sc::UnitCommitmentScenario )::Nothing prod_above = _init(model, :prod_above) segprod = _init(model, :segprod) for t in 1:model[:instance].time for k in 1:length(g.cost_segments) - segprod[g.name, t, k] = @variable(model, lower_bound = 0) + segprod[sc.name, g.name, t, k] = @variable(model, lower_bound = 0) end - prod_above[g.name, t] = @variable(model, lower_bound = 0) + prod_above[sc.name, g.name, t] = @variable(model, lower_bound = 0) end return end @@ -22,16 +23,20 @@ function _add_production_limit_eqs!( model::JuMP.Model, g::Unit, formulation_prod_vars::Gar1962.ProdVars, + sc::UnitCommitmentScenario )::Nothing eq_prod_limit = _init(model, :eq_prod_limit) is_on = model[:is_on] prod_above = model[:prod_above] - reserve = _total_reserves(model, g) + reserve = _total_reserves(model, g, sc) gn = g.name for t in 1:model[:instance].time # Objective function terms for production costs # Part of (69) of Kneuven et al. (2020) as C^R_g * u_g(t) term - add_to_expression!(model[:obj], is_on[gn, t], g.min_power_cost[t]) + + ### Moving this term to another function + # add_to_expression!(model[:obj], is_on[gn, t], g.min_power_cost[t]) + ### # Production limit # Equation (18) in Kneuven et al. (2020) @@ -42,9 +47,9 @@ function _add_production_limit_eqs!( if power_diff < 1e-7 power_diff = 0.0 end - eq_prod_limit[gn, t] = @constraint( + eq_prod_limit[sc.name, gn, t] = @constraint( model, - prod_above[gn, t] + reserve[t] <= power_diff * is_on[gn, t] + prod_above[sc.name, gn, t] + reserve[t] <= power_diff * is_on[gn, t] ) end end diff --git a/src/model/formulations/Gar1962/pwlcosts.jl b/src/model/formulations/Gar1962/pwlcosts.jl index 3ac4871..1566ac7 100644 --- a/src/model/formulations/Gar1962/pwlcosts.jl +++ b/src/model/formulations/Gar1962/pwlcosts.jl @@ -8,6 +8,7 @@ function _add_production_piecewise_linear_eqs!( formulation_prod_vars::Gar1962.ProdVars, formulation_pwl_costs::Gar1962.PwlCosts, formulation_status_vars::Gar1962.StatusVars, + sc::UnitCommitmentScenario )::Nothing eq_prod_above_def = _init(model, :eq_prod_above_def) eq_segprod_limit = _init(model, :eq_segprod_limit) @@ -24,9 +25,9 @@ function _add_production_piecewise_linear_eqs!( for t in 1:model[:instance].time # Definition of production # Equation (43) in Kneuven et al. (2020) - eq_prod_above_def[gn, t] = @constraint( + eq_prod_above_def[sc.name, gn, t] = @constraint( model, - prod_above[gn, t] == sum(segprod[gn, t, k] for k in 1:K) + prod_above[sc.name, gn, t] == sum(segprod[sc.name, gn, t, k] for k in 1:K) ) for k in 1:K @@ -37,21 +38,21 @@ function _add_production_piecewise_linear_eqs!( # difference between max power for segments k and k-1 so the # value of cost_segments[k].mw[t] is the max production *for # that segment* - eq_segprod_limit[gn, t, k] = @constraint( + eq_segprod_limit[sc.name, gn, t, k] = @constraint( model, - segprod[gn, t, k] <= g.cost_segments[k].mw[t] * is_on[gn, t] + segprod[sc.name, gn, t, k] <= g.cost_segments[k].mw[t] * is_on[gn, t] ) # Also add this as an explicit upper bound on segprod to make the # solver's work a bit easier - set_upper_bound(segprod[gn, t, k], g.cost_segments[k].mw[t]) + set_upper_bound(segprod[sc.name, gn, t, k], g.cost_segments[k].mw[t]) # Objective function # Equation (44) in Kneuven et al. (2020) add_to_expression!( model[:obj], - segprod[gn, t, k], - g.cost_segments[k].cost[t], + segprod[sc.name, gn, t, k], + sc.probability * g.cost_segments[k].cost[t], ) end end diff --git a/src/model/formulations/Gar1962/status.jl b/src/model/formulations/Gar1962/status.jl index 14c055f..2a4a911 100644 --- a/src/model/formulations/Gar1962/status.jl +++ b/src/model/formulations/Gar1962/status.jl @@ -20,6 +20,7 @@ function _add_status_vars!( switch_on[g.name, t] = @variable(model, binary = true) switch_off[g.name, t] = @variable(model, binary = true) end + add_to_expression!(model[:obj], is_on[g.name, t], g.min_power_cost[t]) end return end diff --git a/src/model/formulations/KnuOstWat2018/pwlcosts.jl b/src/model/formulations/KnuOstWat2018/pwlcosts.jl index 85afa1e..ab95e05 100644 --- a/src/model/formulations/KnuOstWat2018/pwlcosts.jl +++ b/src/model/formulations/KnuOstWat2018/pwlcosts.jl @@ -8,6 +8,7 @@ function _add_production_piecewise_linear_eqs!( formulation_prod_vars::Gar1962.ProdVars, formulation_pwl_costs::KnuOstWat2018.PwlCosts, formulation_status_vars::Gar1962.StatusVars, + sc::UnitCommitmentScenario )::Nothing eq_prod_above_def = _init(model, :eq_prod_above_def) eq_segprod_limit_a = _init(model, :eq_segprod_limit_a) @@ -58,27 +59,27 @@ function _add_production_piecewise_linear_eqs!( if g.min_uptime > 1 # Equation (46) in Kneuven et al. (2020) - eq_segprod_limit_a[gn, t, k] = @constraint( + eq_segprod_limit_a[sc.name, gn, t, k] = @constraint( model, - segprod[gn, t, k] <= + segprod[sc.name, gn, t, k] <= g.cost_segments[k].mw[t] * is_on[gn, t] - Cv * switch_on[gn, t] - (t < T ? Cw * switch_off[gn, t+1] : 0.0) ) else # Equation (47a)/(48a) in Kneuven et al. (2020) - eq_segprod_limit_b[gn, t, k] = @constraint( + eq_segprod_limit_b[sc.name, gn, t, k] = @constraint( model, - segprod[gn, t, k] <= + segprod[sc.name, gn, t, k] <= g.cost_segments[k].mw[t] * is_on[gn, t] - Cv * switch_on[gn, t] - (t < T ? max(0, Cv - Cw) * switch_off[gn, t+1] : 0.0) ) # Equation (47b)/(48b) in Kneuven et al. (2020) - eq_segprod_limit_c[gn, t, k] = @constraint( + eq_segprod_limit_c[sc.name, gn, t, k] = @constraint( model, - segprod[gn, t, k] <= + segprod[sc.name, gn, t, k] <= g.cost_segments[k].mw[t] * is_on[gn, t] - max(0, Cw - Cv) * switch_on[gn, t] - (t < T ? Cw * switch_off[gn, t+1] : 0.0) @@ -87,22 +88,22 @@ function _add_production_piecewise_linear_eqs!( # Definition of production # Equation (43) in Kneuven et al. (2020) - eq_prod_above_def[gn, t] = @constraint( + eq_prod_above_def[sc.name, gn, t] = @constraint( model, - prod_above[gn, t] == sum(segprod[gn, t, k] for k in 1:K) + prod_above[sc.name, gn, t] == sum(segprod[sc.name, gn, t, k] for k in 1:K) ) # Objective function # Equation (44) in Kneuven et al. (2020) add_to_expression!( model[:obj], - segprod[gn, t, k], + segprod[sc.name, gn, t, k], g.cost_segments[k].cost[t], ) # Also add an explicit upper bound on segprod to make the solver's # work a bit easier - set_upper_bound(segprod[gn, t, k], g.cost_segments[k].mw[t]) + set_upper_bound(segprod[sc.name, gn, t, k], g.cost_segments[k].mw[t]) end end end diff --git a/src/model/formulations/MorLatRam2013/ramp.jl b/src/model/formulations/MorLatRam2013/ramp.jl index 5eaa178..0672011 100644 --- a/src/model/formulations/MorLatRam2013/ramp.jl +++ b/src/model/formulations/MorLatRam2013/ramp.jl @@ -8,6 +8,7 @@ function _add_ramp_eqs!( formulation_prod_vars::Gar1962.ProdVars, formulation_ramping::MorLatRam2013.Ramping, formulation_status_vars::Gar1962.StatusVars, + sc::UnitCommitmentScenario )::Nothing # TODO: Move upper case constants to model[:instance] RESERVES_WHEN_START_UP = true @@ -22,7 +23,7 @@ function _add_ramp_eqs!( gn = g.name eq_ramp_down = _init(model, :eq_ramp_down) eq_ramp_up = _init(model, :eq_str_ramp_up) - reserve = _total_reserves(model, g) + reserve = _total_reserves(model, g, sc) # Gar1962.ProdVars prod_above = model[:prod_above] @@ -39,10 +40,10 @@ function _add_ramp_eqs!( # Ramp up limit if t == 1 if is_initially_on - eq_ramp_up[gn, t] = @constraint( + eq_ramp_up[sc.name, gn, t] = @constraint( model, g.min_power[t] + - prod_above[gn, t] + + prod_above[sc.name, gn, t] + (RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) <= g.initial_power + RU ) @@ -58,13 +59,13 @@ function _add_ramp_eqs!( SU = g.startup_limit max_prod_this_period = g.min_power[t] * is_on[gn, t] + - prod_above[gn, t] + + prod_above[sc.name, gn, t] + ( RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0 ) min_prod_last_period = - g.min_power[t-1] * is_on[gn, t-1] + prod_above[gn, t-1] + g.min_power[t-1] * is_on[gn, t-1] + prod_above[sc.name, gn, t-1] eq_ramp_up[gn, t] = @constraint( model, max_prod_this_period - min_prod_last_period <= @@ -74,11 +75,11 @@ function _add_ramp_eqs!( # Equation (26) in Kneuven et al. (2020) # TODO: what if RU < SU? places too stringent upper bound # prod_above[gn, t] when starting up, and creates diff with (24). - eq_ramp_up[gn, t] = @constraint( + eq_ramp_up[sc.name, gn, t] = @constraint( model, - prod_above[gn, t] + + prod_above[sc.name, gn, t] + (RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) - - prod_above[gn, t-1] <= RU + prod_above[sc.name, gn, t-1] <= RU ) end end @@ -90,9 +91,9 @@ function _add_ramp_eqs!( # min_power + RD < initial_power < SD # then the generator should be able to shut down at time t = 1, # but the constraint below will force the unit to produce power - eq_ramp_down[gn, t] = @constraint( + eq_ramp_down[sc.name, gn, t] = @constraint( model, - g.initial_power - (g.min_power[t] + prod_above[gn, t]) <= RD + g.initial_power - (g.min_power[t] + prod_above[sc.name, gn, t]) <= RD ) end else @@ -102,13 +103,13 @@ function _add_ramp_eqs!( SD = g.shutdown_limit max_prod_last_period = g.min_power[t-1] * is_on[gn, t-1] + - prod_above[gn, t-1] + + prod_above[sc.name, gn, t-1] + ( RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ? reserve[t-1] : 0.0 ) min_prod_this_period = - g.min_power[t] * is_on[gn, t] + prod_above[gn, t] + g.min_power[t] * is_on[gn, t] + prod_above[sc.name, gn, t] eq_ramp_down[gn, t] = @constraint( model, max_prod_last_period - min_prod_this_period <= @@ -118,11 +119,11 @@ function _add_ramp_eqs!( # Equation (27) in Kneuven et al. (2020) # TODO: Similar to above, what to do if shutting down in time t # and RD < SD? There is a difference with (25). - eq_ramp_down[gn, t] = @constraint( + eq_ramp_down[sc.name, gn, t] = @constraint( model, - prod_above[gn, t-1] + + prod_above[sc.name, gn, t-1] + (RESERVES_WHEN_RAMP_DOWN ? reserve[t-1] : 0.0) - - prod_above[gn, t] <= RD + prod_above[sc.name, gn, t] <= RD ) end end diff --git a/src/model/formulations/PanGua2016/ramp.jl b/src/model/formulations/PanGua2016/ramp.jl index fe5d334..eb85ff5 100644 --- a/src/model/formulations/PanGua2016/ramp.jl +++ b/src/model/formulations/PanGua2016/ramp.jl @@ -8,11 +8,12 @@ function _add_ramp_eqs!( formulation_prod_vars::Gar1962.ProdVars, formulation_ramping::PanGua2016.Ramping, formulation_status_vars::Gar1962.StatusVars, + sc::UnitCommitmentScenario )::Nothing # TODO: Move upper case constants to model[:instance] RESERVES_WHEN_SHUT_DOWN = true gn = g.name - reserve = _total_reserves(model, g) + reserve = _total_reserves(model, g, sc) eq_str_prod_limit = _init(model, :eq_str_prod_limit) eq_prod_limit_ramp_up_extra_period = _init(model, :eq_prod_limit_ramp_up_extra_period) @@ -52,9 +53,9 @@ function _add_ramp_eqs!( # Generalization of (20) # Necessary that if any of the switch_on = 1 in the sum, # then switch_off[gn, t+1] = 0 - eq_str_prod_limit[gn, t] = @constraint( + eq_str_prod_limit[sc.name, gn, t] = @constraint( model, - prod_above[gn, t] + + prod_above[sc.name, gn, t] + g.min_power[t] * is_on[gn, t] + reserve[t] <= Pbar * is_on[gn, t] - @@ -67,9 +68,9 @@ function _add_ramp_eqs!( if UT - 2 < TRU # Equation (40) in Kneuven et al. (2020) # Covers an additional time period of the ramp-up trajectory, compared to (38) - eq_prod_limit_ramp_up_extra_period[gn, t] = @constraint( + eq_prod_limit_ramp_up_extra_period[sc.name, gn, t] = @constraint( model, - prod_above[gn, t] + + prod_above[sc.name, gn, t] + g.min_power[t] * is_on[gn, t] + reserve[t] <= Pbar * is_on[gn, t] - sum( @@ -84,9 +85,9 @@ function _add_ramp_eqs!( if KSD > 0 KSU = min(TRU, UT - 2 - KSD, t - 1) # Equation (41) in Kneuven et al. (2020) - eq_prod_limit_shutdown_trajectory[gn, t] = @constraint( + eq_prod_limit_shutdown_trajectory[sc.name, gn, t] = @constraint( model, - prod_above[gn, t] + + prod_above[sc.name, gn, t] + g.min_power[t] * is_on[gn, t] + (RESERVES_WHEN_SHUT_DOWN ? reserve[t] : 0.0) <= Pbar * is_on[gn, t] - sum( diff --git a/src/model/formulations/WanHob2016/ramp.jl b/src/model/formulations/WanHob2016/ramp.jl index e36e4d2..cc92ae9 100644 --- a/src/model/formulations/WanHob2016/ramp.jl +++ b/src/model/formulations/WanHob2016/ramp.jl @@ -8,6 +8,7 @@ function _add_ramp_eqs!( ::Gar1962.ProdVars, ::WanHob2016.Ramping, ::Gar1962.StatusVars, + sc::UnitCommitmentScenario )::Nothing is_initially_on = (g.initial_status > 0) SU = g.startup_limit @@ -29,7 +30,7 @@ function _add_ramp_eqs!( error("Each generator may only provide one flexiramp reserve") end for r in g.reserves - if r.type !== "flexiramp" + if r.type !== "up-frp" && r.type !== "down-frp" error( "This formulation only supports flexiramp reserves, not $(r.type)", ) @@ -38,41 +39,41 @@ function _add_ramp_eqs!( for t in 1:model[:instance].time @constraint( model, - prod_above[gn, t] + (is_on[gn, t] * minp[t]) <= mfg[rn, gn, t] + prod_above[sc.name, gn, t] + (is_on[gn, t] * minp[t]) <= mfg[sc.name, rn, gn, t] ) # Eq. (19) in Wang & Hobbs (2016) - @constraint(model, mfg[rn, gn, t] <= is_on[gn, t] * maxp[t]) # Eq. (22) in Wang & Hobbs (2016) + @constraint(model, mfg[sc.name, rn, gn, t] <= is_on[gn, t] * maxp[t]) # Eq. (22) in Wang & Hobbs (2016) if t != model[:instance].time @constraint( model, minp[t] * (is_on[gn, t+1] + is_on[gn, t] - 1) <= - prod_above[gn, t] - dwflexiramp[rn, gn, t] + + prod_above[sc.name, gn, t] - dwflexiramp[sc.name, rn, gn, t] + (is_on[gn, t] * minp[t]) ) # first inequality of Eq. (20) in Wang & Hobbs (2016) @constraint( model, - prod_above[gn, t] - dwflexiramp[rn, gn, t] + + prod_above[sc.name, gn, t] - dwflexiramp[sc.name, rn, gn, t] + (is_on[gn, t] * minp[t]) <= - mfg[rn, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1])) + mfg[sc.name, rn, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1])) ) # second inequality of Eq. (20) in Wang & Hobbs (2016) @constraint( model, minp[t] * (is_on[gn, t+1] + is_on[gn, t] - 1) <= - prod_above[gn, t] + - upflexiramp[rn, gn, t] + + prod_above[sc.name, gn, t] + + upflexiramp[sc.name, rn, gn, t] + (is_on[gn, t] * minp[t]) ) # first inequality of Eq. (21) in Wang & Hobbs (2016) @constraint( model, - prod_above[gn, t] + - upflexiramp[rn, gn, t] + + prod_above[sc.name, gn, t] + + upflexiramp[sc.name, rn, gn, t] + (is_on[gn, t] * minp[t]) <= - mfg[rn, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1])) + mfg[sc.name, rn, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1])) ) # second inequality of Eq. (21) in Wang & Hobbs (2016) if t != 1 @constraint( model, - mfg[rn, gn, t] <= - prod_above[gn, t-1] + + mfg[sc.name, rn, gn, t] <= + prod_above[sc.name, gn, t-1] + (is_on[gn, t-1] * minp[t]) + (RU * is_on[gn, t-1]) + (SU * (is_on[gn, t] - is_on[gn, t-1])) + @@ -80,8 +81,8 @@ function _add_ramp_eqs!( ) # Eq. (23) in Wang & Hobbs (2016) @constraint( model, - (prod_above[gn, t-1] + (is_on[gn, t-1] * minp[t])) - - (prod_above[gn, t] + (is_on[gn, t] * minp[t])) <= + (prod_above[sc.name, gn, t-1] + (is_on[gn, t-1] * minp[t])) - + (prod_above[sc.name, gn, t] + (is_on[gn, t] * minp[t])) <= RD * is_on[gn, t] + SD * (is_on[gn, t-1] - is_on[gn, t]) + maxp[t] * (1 - is_on[gn, t-1]) @@ -89,7 +90,7 @@ function _add_ramp_eqs!( else @constraint( model, - mfg[rn, gn, t] <= + mfg[sc.name, rn, gn, t] <= initial_power + (RU * is_initially_on) + (SU * (is_on[gn, t] - is_initially_on)) + @@ -98,7 +99,7 @@ function _add_ramp_eqs!( @constraint( model, initial_power - - (prod_above[gn, t] + (is_on[gn, t] * minp[t])) <= + (prod_above[sc.name, gn, t] + (is_on[gn, t] * minp[t])) <= RD * is_on[gn, t] + SD * (is_initially_on - is_on[gn, t]) + maxp[t] * (1 - is_initially_on) @@ -106,7 +107,7 @@ function _add_ramp_eqs!( end @constraint( model, - mfg[rn, gn, t] <= + mfg[sc.name, rn, gn, t] <= (SD * (is_on[gn, t] - is_on[gn, t+1])) + (maxp[t] * is_on[gn, t+1]) ) # Eq. (24) in Wang & Hobbs (2016) @@ -114,11 +115,11 @@ function _add_ramp_eqs!( model, -RD * is_on[gn, t+1] - SD * (is_on[gn, t] - is_on[gn, t+1]) - - maxp[t] * (1 - is_on[gn, t]) <= upflexiramp[rn, gn, t] + maxp[t] * (1 - is_on[gn, t]) <= upflexiramp[sc.name, rn, gn, t] ) # first inequality of Eq. (26) in Wang & Hobbs (2016) @constraint( model, - upflexiramp[rn, gn, t] <= + upflexiramp[sc.name, rn, gn, t] <= RU * is_on[gn, t] + SU * (is_on[gn, t+1] - is_on[gn, t]) + maxp[t] * (1 - is_on[gn, t+1]) @@ -126,11 +127,11 @@ function _add_ramp_eqs!( @constraint( model, -RU * is_on[gn, t] - SU * (is_on[gn, t+1] - is_on[gn, t]) - - maxp[t] * (1 - is_on[gn, t+1]) <= dwflexiramp[rn, gn, t] + maxp[t] * (1 - is_on[gn, t+1]) <= dwflexiramp[sc.name, rn, gn, t] ) # first inequality of Eq. (27) in Wang & Hobbs (2016) @constraint( model, - dwflexiramp[rn, gn, t] <= + dwflexiramp[sc.name, rn, gn, t] <= RD * is_on[gn, t+1] + SD * (is_on[gn, t] - is_on[gn, t+1]) + maxp[t] * (1 - is_on[gn, t]) @@ -138,26 +139,26 @@ function _add_ramp_eqs!( @constraint( model, -maxp[t] * is_on[gn, t] + minp[t] * is_on[gn, t+1] <= - upflexiramp[rn, gn, t] + upflexiramp[sc.name, rn, gn, t] ) # first inequality of Eq. (28) in Wang & Hobbs (2016) @constraint( model, - upflexiramp[rn, gn, t] <= maxp[t] * is_on[gn, t+1] + upflexiramp[sc.name, rn, gn, t] <= maxp[t] * is_on[gn, t+1] ) # second inequality of Eq. (28) in Wang & Hobbs (2016) @constraint( model, - -maxp[t] * is_on[gn, t+1] <= dwflexiramp[rn, gn, t] + -maxp[t] * is_on[gn, t+1] <= dwflexiramp[sc.name, rn, gn, t] ) # first inequality of Eq. (29) in Wang & Hobbs (2016) @constraint( model, - dwflexiramp[rn, gn, t] <= + dwflexiramp[sc.name, rn, gn, t] <= (maxp[t] * is_on[gn, t]) - (minp[t] * is_on[gn, t+1]) ) # second inequality of Eq. (29) in Wang & Hobbs (2016) else @constraint( model, - mfg[rn, gn, t] <= - prod_above[gn, t-1] + + mfg[sc.name, rn, gn, t] <= + prod_above[sc.name, gn, t-1] + (is_on[gn, t-1] * minp[t]) + (RU * is_on[gn, t-1]) + (SU * (is_on[gn, t] - is_on[gn, t-1])) + @@ -165,8 +166,8 @@ function _add_ramp_eqs!( ) # Eq. (23) in Wang & Hobbs (2016) for the last time period @constraint( model, - (prod_above[gn, t-1] + (is_on[gn, t-1] * minp[t])) - - (prod_above[gn, t] + (is_on[gn, t] * minp[t])) <= + (prod_above[sc.name, gn, t-1] + (is_on[gn, t-1] * minp[t])) - + (prod_above[sc.name, gn, t] + (is_on[gn, t] * minp[t])) <= RD * is_on[gn, t] + SD * (is_on[gn, t-1] - is_on[gn, t]) + maxp[t] * (1 - is_on[gn, t-1]) diff --git a/src/model/formulations/base/bus.jl b/src/model/formulations/base/bus.jl index d881fc6..926bc2c 100644 --- a/src/model/formulations/base/bus.jl +++ b/src/model/formulations/base/bus.jl @@ -2,22 +2,22 @@ # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. # Released under the modified BSD license. See COPYING.md for more details. -function _add_bus!(model::JuMP.Model, b::Bus)::Nothing +function _add_bus!(model::JuMP.Model, b::Bus, sc::UnitCommitmentScenario)::Nothing net_injection = _init(model, :expr_net_injection) curtail = _init(model, :curtail) for t in 1:model[:instance].time # Fixed load - net_injection[b.name, t] = AffExpr(-b.load[t]) + net_injection[sc.name, b.name, t] = AffExpr(-b.load[t]) # Load curtailment - curtail[b.name, t] = + curtail[sc.name, b.name, t] = @variable(model, lower_bound = 0, upper_bound = b.load[t]) - add_to_expression!(net_injection[b.name, t], curtail[b.name, t], 1.0) + add_to_expression!(net_injection[sc.name, b.name, t], curtail[sc.name, b.name, t], 1.0) add_to_expression!( model[:obj], - curtail[b.name, t], - model[:instance].power_balance_penalty[t], + curtail[sc.name, b.name, t], + sc.power_balance_penalty[t] * sc.probability, ) end return diff --git a/src/model/formulations/base/line.jl b/src/model/formulations/base/line.jl index 6398c8d..0a4a586 100644 --- a/src/model/formulations/base/line.jl +++ b/src/model/formulations/base/line.jl @@ -6,14 +6,15 @@ function _add_transmission_line!( model::JuMP.Model, lm::TransmissionLine, f::ShiftFactorsFormulation, + sc::UnitCommitmentScenario )::Nothing overflow = _init(model, :overflow) for t in 1:model[:instance].time - overflow[lm.name, t] = @variable(model, lower_bound = 0) + overflow[sc.name, lm.name, t] = @variable(model, lower_bound = 0) add_to_expression!( model[:obj], - overflow[lm.name, t], - lm.flow_limit_penalty[t], + overflow[sc.name, lm.name, t], + lm.flow_limit_penalty[t] * sc.probability, ) end return @@ -22,27 +23,28 @@ end function _setup_transmission( model::JuMP.Model, formulation::ShiftFactorsFormulation, + scenario::UnitCommitmentScenario )::Nothing instance = model[:instance] isf = formulation.precomputed_isf lodf = formulation.precomputed_lodf - if length(instance.buses) == 1 + if length(scenario.buses) == 1 isf = zeros(0, 0) lodf = zeros(0, 0) elseif isf === nothing @info "Computing injection shift factors..." time_isf = @elapsed begin isf = UnitCommitment._injection_shift_factors( - lines = instance.lines, - buses = instance.buses, + lines = scenario.lines, + buses = scenario.buses, ) end @info @sprintf("Computed ISF in %.2f seconds", time_isf) @info "Computing line outage factors..." time_lodf = @elapsed begin lodf = UnitCommitment._line_outage_factors( - lines = instance.lines, - buses = instance.buses, + lines = scenario.lines, + buses = scenario.buses, isf = isf, ) end @@ -55,7 +57,7 @@ function _setup_transmission( isf[abs.(isf).= r.amount[t] + sum(model[:reserve][sc.name, r.name, g.name, t] for g in r.units) + + model[:reserve_shortfall][sc.name, r.name, t] >= r.amount[t] ) # Account for shortfall contribution to objective if r.shortfall_penalty >= 0 add_to_expression!( model[:obj], - r.shortfall_penalty, - model[:reserve_shortfall][r.name, t], + r.shortfall_penalty * sc.probability, + model[:reserve_shortfall][sc.name, r.name, t], ) end end @@ -57,7 +57,7 @@ function _add_spinning_reserve_eqs!(model::JuMP.Model)::Nothing return end -function _add_flexiramp_reserve_eqs!(model::JuMP.Model)::Nothing +function _add_flexiramp_reserve_eqs!(model::JuMP.Model, sc::UnitCommitmentScenario)::Nothing # Note: The flexpramp requirements in Wang & Hobbs (2016) are imposed as hard constraints # through Eq. (17) and Eq. (18). The constraints eq_min_upflexiramp and eq_min_dwflexiramp # provided below are modified versions of Eq. (17) and Eq. (18), respectively, in that @@ -65,31 +65,41 @@ function _add_flexiramp_reserve_eqs!(model::JuMP.Model)::Nothing # objective function. eq_min_upflexiramp = _init(model, :eq_min_upflexiramp) eq_min_dwflexiramp = _init(model, :eq_min_dwflexiramp) - instance = model[:instance] - for r in instance.reserves - r.type == "flexiramp" || continue - for t in 1:instance.time - # Eq. (17) in Wang & Hobbs (2016) - eq_min_upflexiramp[r.name, t] = @constraint( - model, - sum(model[:upflexiramp][r.name, g.name, t] for g in r.units) + model[:upflexiramp_shortfall][r.name, t] >= r.amount[t] - ) - # Eq. (18) in Wang & Hobbs (2016) - eq_min_dwflexiramp[r.name, t] = @constraint( - model, - sum(model[:dwflexiramp][r.name, g.name, t] for g in r.units) + model[:dwflexiramp_shortfall][r.name, t] >= r.amount[t] - ) - - # Account for flexiramp shortfall contribution to objective - if r.shortfall_penalty >= 0 - add_to_expression!( - model[:obj], - r.shortfall_penalty, - ( - model[:upflexiramp_shortfall][r.name, t] + - model[:dwflexiramp_shortfall][r.name, t] - ), + T = model[:instance].time + for r in sc.reserves + if r.type == "up-frp" + for t in 1:T + # Eq. (17) in Wang & Hobbs (2016) + eq_min_upflexiramp[sc.name, r.name, t] = @constraint( + model, + sum(model[:upflexiramp][sc.name, r.name, g.name, t] for g in r.units) + + model[:upflexiramp_shortfall][sc.name, r.name, t] >= r.amount[t] + ) + # Account for flexiramp shortfall contribution to objective + if r.shortfall_penalty >= 0 + add_to_expression!( + model[:obj], + r.shortfall_penalty * sc.probability, + model[:upflexiramp_shortfall][sc.name, r.name, t] + ) + end + end + elseif r.type == "down-frp" + for t in 1:T + # Eq. (18) in Wang & Hobbs (2016) + eq_min_dwflexiramp[sc.name, r.name, t] = @constraint( + model, + sum(model[:dwflexiramp][sc.name, r.name, g.name, t] for g in r.units) + + model[:dwflexiramp_shortfall][sc.name, r.name, t] >= r.amount[t] ) + # Account for flexiramp shortfall contribution to objective + if r.shortfall_penalty >= 0 + add_to_expression!( + model[:obj], + r.shortfall_penalty * sc.probability, + model[:dwflexiramp_shortfall][sc.name, r.name, t] + ) + end end end end diff --git a/src/model/formulations/base/unit.jl b/src/model/formulations/base/unit.jl index bcccf16..788ac0a 100644 --- a/src/model/formulations/base/unit.jl +++ b/src/model/formulations/base/unit.jl @@ -2,7 +2,7 @@ # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. # Released under the modified BSD license. See COPYING.md for more details. -function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation) +function _add_unit_first_stage!(model::JuMP.Model, g::Unit, formulation::Formulation) if !all(g.must_run) && any(g.must_run) error("Partially must-run units are not currently supported") end @@ -11,22 +11,34 @@ function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation) end # Variables - _add_production_vars!(model, g, formulation.prod_vars) - _add_spinning_reserve_vars!(model, g) - _add_flexiramp_reserve_vars!(model, g) _add_startup_shutdown_vars!(model, g) _add_status_vars!(model, g, formulation.status_vars) # Constraints and objective function _add_min_uptime_downtime_eqs!(model, g) - _add_net_injection_eqs!(model, g) - _add_production_limit_eqs!(model, g, formulation.prod_vars) + _add_startup_cost_eqs!(model, g, formulation.startup_costs) + _add_status_eqs!(model, g, formulation.status_vars) + return +end + +function _add_unit_second_stage!(model::JuMP.Model, g::Unit, formulation::Formulation, + scenario::UnitCommitmentScenario) + + # Variables + _add_production_vars!(model, g, formulation.prod_vars, scenario) + _add_spinning_reserve_vars!(model, g, scenario) + _add_flexiramp_reserve_vars!(model, g, scenario) + + # Constraints and objective function + _add_net_injection_eqs!(model, g, scenario) + _add_production_limit_eqs!(model, g, formulation.prod_vars, scenario) _add_production_piecewise_linear_eqs!( model, g, formulation.prod_vars, formulation.pwl_costs, formulation.status_vars, + scenario ) _add_ramp_eqs!( model, @@ -34,26 +46,64 @@ function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation) formulation.prod_vars, formulation.ramping, formulation.status_vars, + scenario ) - _add_startup_cost_eqs!(model, g, formulation.startup_costs) - _add_startup_shutdown_limit_eqs!(model, g) - _add_status_eqs!(model, g, formulation.status_vars) + _add_startup_shutdown_limit_eqs!(model, g, scenario) return end +# function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation) +# if !all(g.must_run) && any(g.must_run) +# error("Partially must-run units are not currently supported") +# end +# if g.initial_power === nothing || g.initial_status === nothing +# error("Initial conditions for $(g.name) must be provided") +# end + +# # Variables +# _add_production_vars!(model, g, formulation.prod_vars) +# _add_spinning_reserve_vars!(model, g) +# _add_flexiramp_reserve_vars!(model, g) +# _add_startup_shutdown_vars!(model, g) +# _add_status_vars!(model, g, formulation.status_vars) + +# # Constraints and objective function +# _add_min_uptime_downtime_eqs!(model, g) +# _add_net_injection_eqs!(model, g) +# _add_production_limit_eqs!(model, g, formulation.prod_vars) +# _add_production_piecewise_linear_eqs!( +# model, +# g, +# formulation.prod_vars, +# formulation.pwl_costs, +# formulation.status_vars, +# ) +# _add_ramp_eqs!( +# model, +# g, +# formulation.prod_vars, +# formulation.ramping, +# formulation.status_vars, +# ) +# _add_startup_cost_eqs!(model, g, formulation.startup_costs) +# _add_startup_shutdown_limit_eqs!(model, g) +# _add_status_eqs!(model, g, formulation.status_vars) +# return +# end + _is_initially_on(g::Unit)::Float64 = (g.initial_status > 0 ? 1.0 : 0.0) -function _add_spinning_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing +function _add_spinning_reserve_vars!(model::JuMP.Model, g::Unit, sc::UnitCommitmentScenario)::Nothing reserve = _init(model, :reserve) reserve_shortfall = _init(model, :reserve_shortfall) for r in g.reserves r.type == "spinning" || continue for t in 1:model[:instance].time - reserve[r.name, g.name, t] = @variable(model, lower_bound = 0) - if (r.name, t) ∉ keys(reserve_shortfall) - reserve_shortfall[r.name, t] = @variable(model, lower_bound = 0) + reserve[sc.name, r.name, g.name, t] = @variable(model, lower_bound = 0) + if (sc.name, r.name, t) ∉ keys(reserve_shortfall) + reserve_shortfall[sc.name, r.name, t] = @variable(model, lower_bound = 0) if r.shortfall_penalty < 0 - set_upper_bound(reserve_shortfall[r.name, t], 0.0) + set_upper_bound(reserve_shortfall[sc.name, r.name, t], 0.0) end end end @@ -61,27 +111,35 @@ function _add_spinning_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing return end -function _add_flexiramp_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing +function _add_flexiramp_reserve_vars!(model::JuMP.Model, g::Unit, sc::UnitCommitmentScenario)::Nothing upflexiramp = _init(model, :upflexiramp) upflexiramp_shortfall = _init(model, :upflexiramp_shortfall) mfg = _init(model, :mfg) dwflexiramp = _init(model, :dwflexiramp) dwflexiramp_shortfall = _init(model, :dwflexiramp_shortfall) for r in g.reserves - r.type == "flexiramp" || continue - for t in 1:model[:instance].time - # maximum feasible generation, \bar{g_{its}} in Wang & Hobbs (2016) - mfg[r.name, g.name, t] = @variable(model, lower_bound = 0) - upflexiramp[r.name, g.name, t] = @variable(model) # up-flexiramp, ur_{it} in Wang & Hobbs (2016) - dwflexiramp[r.name, g.name, t] = @variable(model) # down-flexiramp, dr_{it} in Wang & Hobbs (2016) - if (r.name, t) ∉ keys(upflexiramp_shortfall) - upflexiramp_shortfall[r.name, t] = - @variable(model, lower_bound = 0) - dwflexiramp_shortfall[r.name, t] = - @variable(model, lower_bound = 0) - if r.shortfall_penalty < 0 - set_upper_bound(upflexiramp_shortfall[r.name, t], 0.0) - set_upper_bound(dwflexiramp_shortfall[r.name, t], 0.0) + if r.type == "up-frp" + for t in 1:model[:instance].time + # maximum feasible generation, \bar{g_{its}} in Wang & Hobbs (2016) + mfg[sc.name, r.name, g.name, t] = @variable(model, lower_bound = 0) + upflexiramp[sc.name, r.name, g.name, t] = @variable(model) # up-flexiramp, ur_{it} in Wang & Hobbs (2016) + if (sc.name, r.name, t) ∉ keys(upflexiramp_shortfall) + upflexiramp_shortfall[sc.name, r.name, t] = + @variable(model, lower_bound = 0) + if r.shortfall_penalty < 0 + set_upper_bound(upflexiramp_shortfall[sc.name, r.name, t], 0.0) + end + end + end + elseif r.type == "down-frp" + for t in 1:model[:instance].time + dwflexiramp[sc.name, r.name, g.name, t] = @variable(model) # down-flexiramp, dr_{it} in Wang & Hobbs (2016) + if (sc.name, r.name, t) ∉ keys(dwflexiramp_shortfall) + dwflexiramp_shortfall[sc.name, r.name, t] = + @variable(model, lower_bound = 0) + if r.shortfall_penalty < 0 + set_upper_bound(dwflexiramp_shortfall[sc.name, r.name, t], 0.0) + end end end end @@ -99,32 +157,32 @@ function _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing return end -function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing +function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit, sc::UnitCommitmentScenario)::Nothing eq_shutdown_limit = _init(model, :eq_shutdown_limit) eq_startup_limit = _init(model, :eq_startup_limit) is_on = model[:is_on] prod_above = model[:prod_above] - reserve = _total_reserves(model, g) + reserve = _total_reserves(model, g, sc) switch_off = model[:switch_off] switch_on = model[:switch_on] T = model[:instance].time for t in 1:T # Startup limit - eq_startup_limit[g.name, t] = @constraint( + eq_startup_limit[sc.name, g.name, t] = @constraint( model, - prod_above[g.name, t] + reserve[t] <= + prod_above[sc.name, g.name, t] + reserve[t] <= (g.max_power[t] - g.min_power[t]) * is_on[g.name, t] - max(0, g.max_power[t] - g.startup_limit) * switch_on[g.name, t] ) # Shutdown limit if g.initial_power > g.shutdown_limit - eq_shutdown_limit[g.name, 0] = + eq_shutdown_limit[sc.name, g.name, 0] = @constraint(model, switch_off[g.name, 1] <= 0) end if t < T - eq_shutdown_limit[g.name, t] = @constraint( + eq_shutdown_limit[sc.name, g.name, t] = @constraint( model, - prod_above[g.name, t] <= + prod_above[sc.name, g.name, t] <= (g.max_power[t] - g.min_power[t]) * is_on[g.name, t] - max(0, g.max_power[t] - g.shutdown_limit) * switch_off[g.name, t+1] @@ -138,43 +196,44 @@ function _add_ramp_eqs!( model::JuMP.Model, g::Unit, formulation::RampingFormulation, + sc::UnitCommitmentScenario )::Nothing prod_above = model[:prod_above] - reserve = _total_reserves(model, g) + reserve = _total_reserves(model, g, sc) eq_ramp_up = _init(model, :eq_ramp_up) eq_ramp_down = _init(model, :eq_ramp_down) for t in 1:model[:instance].time # Ramp up limit if t == 1 if _is_initially_on(g) == 1 - eq_ramp_up[g.name, t] = @constraint( + eq_ramp_up[sc.name, g.name, t] = @constraint( model, - prod_above[g.name, t] + reserve[t] <= + prod_above[sc.name, g.name, t] + reserve[t] <= (g.initial_power - g.min_power[t]) + g.ramp_up_limit ) end else - eq_ramp_up[g.name, t] = @constraint( + eq_ramp_up[sc.name, g.name, t] = @constraint( model, - prod_above[g.name, t] + reserve[t] <= - prod_above[g.name, t-1] + g.ramp_up_limit + prod_above[sc.name, g.name, t] + reserve[t] <= + prod_above[sc.name, g.name, t-1] + g.ramp_up_limit ) end # Ramp down limit if t == 1 if _is_initially_on(g) == 1 - eq_ramp_down[g.name, t] = @constraint( + eq_ramp_down[sc.name, g.name, t] = @constraint( model, - prod_above[g.name, t] >= + prod_above[sc.name, g.name, t] >= (g.initial_power - g.min_power[t]) - g.ramp_down_limit ) end else - eq_ramp_down[g.name, t] = @constraint( + eq_ramp_down[sc.name, g.name, t] = @constraint( model, - prod_above[g.name, t] >= - prod_above[g.name, t-1] - g.ramp_down_limit + prod_above[sc.name, g.name, t] >= + prod_above[sc.name, g.name, t-1] - g.ramp_down_limit ) end end @@ -223,30 +282,30 @@ function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing end end -function _add_net_injection_eqs!(model::JuMP.Model, g::Unit)::Nothing +function _add_net_injection_eqs!(model::JuMP.Model, g::Unit, sc::UnitCommitmentScenario)::Nothing expr_net_injection = model[:expr_net_injection] for t in 1:model[:instance].time # Add to net injection expression add_to_expression!( - expr_net_injection[g.bus.name, t], - model[:prod_above][g.name, t], + expr_net_injection[sc.name, g.bus.name, t], + model[:prod_above][sc.name, g.name, t], 1.0, ) add_to_expression!( - expr_net_injection[g.bus.name, t], + expr_net_injection[sc.name, g.bus.name, t], model[:is_on][g.name, t], g.min_power[t], ) end end -function _total_reserves(model, g)::Vector +function _total_reserves(model, g, sc)::Vector T = model[:instance].time reserve = [0.0 for _ in 1:T] spinning_reserves = [r for r in g.reserves if r.type == "spinning"] if !isempty(spinning_reserves) reserve += [ - sum(model[:reserve][r.name, g.name, t] for r in spinning_reserves) for t in 1:model[:instance].time + sum(model[:reserve][sc.name, r.name, g.name, t] for r in spinning_reserves) for t in 1:model[:instance].time ] end return reserve diff --git a/src/solution/methods/XavQiuWanThi2019/enforce.jl b/src/solution/methods/XavQiuWanThi2019/enforce.jl index 526274f..a4e6ce7 100644 --- a/src/solution/methods/XavQiuWanThi2019/enforce.jl +++ b/src/solution/methods/XavQiuWanThi2019/enforce.jl @@ -5,13 +5,15 @@ function _enforce_transmission( model::JuMP.Model, violations::Vector{_Violation}, + sc::UnitCommitmentScenario )::Nothing for v in violations _enforce_transmission( model = model, + sc = sc, violation = v, - isf = model[:isf], - lodf = model[:lodf], + isf = sc.isf, + lodf = sc.lodf, ) end return @@ -19,6 +21,7 @@ end function _enforce_transmission(; model::JuMP.Model, + sc::UnitCommitmentScenario, violation::_Violation, isf::Matrix{Float64}, lodf::Matrix{Float64}, @@ -51,7 +54,7 @@ function _enforce_transmission(; t = violation.time flow = @variable(model, base_name = "flow[$fm,$t]") - v = overflow[violation.monitored_line.name, violation.time] + v = overflow[sc.name, violation.monitored_line.name, violation.time] @constraint(model, flow <= limit + v) @constraint(model, -flow <= limit + v) @@ -59,23 +62,23 @@ function _enforce_transmission(; @constraint( model, flow == sum( - net_injection[b.name, violation.time] * + net_injection[sc.name, b.name, violation.time] * isf[violation.monitored_line.offset, b.offset] for - b in instance.buses if b.offset > 0 + b in sc.buses if b.offset > 0 ) ) else @constraint( model, flow == sum( - net_injection[b.name, violation.time] * ( + net_injection[sc.name, b.name, violation.time] * ( isf[violation.monitored_line.offset, b.offset] + ( lodf[ violation.monitored_line.offset, violation.outage_line.offset, ] * isf[violation.outage_line.offset, b.offset] ) - ) for b in instance.buses if b.offset > 0 + ) for b in sc.buses if b.offset > 0 ) ) end diff --git a/src/solution/methods/XavQiuWanThi2019/find.jl b/src/solution/methods/XavQiuWanThi2019/find.jl index 3cf9094..43e3ad3 100644 --- a/src/solution/methods/XavQiuWanThi2019/find.jl +++ b/src/solution/methods/XavQiuWanThi2019/find.jl @@ -5,32 +5,34 @@ import Base.Threads: @threads function _find_violations( - model::JuMP.Model; + model::JuMP.Model, + sc::UnitCommitmentScenario; max_per_line::Int, max_per_period::Int, ) instance = model[:instance] net_injection = model[:net_injection] overflow = model[:overflow] - length(instance.buses) > 1 || return [] + length(sc.buses) > 1 || return [] violations = [] @info "Verifying transmission limits..." time_screening = @elapsed begin - non_slack_buses = [b for b in instance.buses if b.offset > 0] + non_slack_buses = [b for b in sc.buses if b.offset > 0] net_injection_values = [ - value(net_injection[b.name, t]) for b in non_slack_buses, + value(net_injection[sc.name, b.name, t]) for b in non_slack_buses, t in 1:instance.time ] overflow_values = [ - value(overflow[lm.name, t]) for lm in instance.lines, + value(overflow[sc.name, lm.name, t]) for lm in sc.lines, t in 1:instance.time ] violations = UnitCommitment._find_violations( instance = instance, + sc = sc, net_injections = net_injection_values, overflow = overflow_values, - isf = model[:isf], - lodf = model[:lodf], + isf = sc.isf, + lodf = sc.lodf, max_per_line = max_per_line, max_per_period = max_per_period, ) @@ -64,15 +66,16 @@ matrix, where L is the number of transmission lines. """ function _find_violations(; instance::UnitCommitmentInstance, + sc::UnitCommitmentScenario, net_injections::Array{Float64,2}, overflow::Array{Float64,2}, isf::Array{Float64,2}, lodf::Array{Float64,2}, max_per_line::Int, - max_per_period::Int, + max_per_period::Int )::Array{_Violation,1} - B = length(instance.buses) - 1 - L = length(instance.lines) + B = length(sc.buses) - 1 + L = length(sc.lines) T = instance.time K = nthreads() @@ -94,16 +97,16 @@ function _find_violations(; normal_limits::Array{Float64,2} = [ l.normal_flow_limit[t] + overflow[l.offset, t] for - l in instance.lines, t in 1:T + l in sc.lines, t in 1:T ] emergency_limits::Array{Float64,2} = [ l.emergency_flow_limit[t] + overflow[l.offset, t] for - l in instance.lines, t in 1:T + l in sc.lines, t in 1:T ] is_vulnerable::Array{Bool} = zeros(Bool, L) - for c in instance.contingencies + for c in sc.contingencies is_vulnerable[c.lines[1].offset] = true end @@ -111,7 +114,7 @@ function _find_violations(; k = threadid() # Pre-contingency flows - pre_flow[:, k] = isf * net_injections[:, t] + pre_flow[:, k] = isf * net_injections[ :, t] # Post-contingency flows for lc in 1:L, lm in 1:L @@ -144,7 +147,7 @@ function _find_violations(; filters[t], _Violation( time = t, - monitored_line = instance.lines[lm], + monitored_line = sc.lines[lm], outage_line = nothing, amount = pre_v[lm, k], ), @@ -159,8 +162,8 @@ function _find_violations(; filters[t], _Violation( time = t, - monitored_line = instance.lines[lm], - outage_line = instance.lines[lc], + monitored_line = sc.lines[lm], + outage_line = sc.lines[lc], amount = post_v[lm, lc, k], ), ) diff --git a/src/solution/methods/XavQiuWanThi2019/optimize.jl b/src/solution/methods/XavQiuWanThi2019/optimize.jl index 78c8381..c7b3a7e 100644 --- a/src/solution/methods/XavQiuWanThi2019/optimize.jl +++ b/src/solution/methods/XavQiuWanThi2019/optimize.jl @@ -10,43 +10,47 @@ function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing JuMP.set_optimizer_attribute(model, "MIPGap", gap) @info @sprintf("MIP gap tolerance set to %f", gap) end - initial_time = time() - large_gap = false - has_transmission = (length(model[:isf]) > 0) - if has_transmission && method.two_phase_gap - set_gap(1e-2) - large_gap = true - end - while true - time_elapsed = time() - initial_time - time_remaining = method.time_limit - time_elapsed - if time_remaining < 0 - @info "Time limit exceeded" - break + for scenario in model[:instance].scenarios + large_gap = false + has_transmission = (length(scenario.isf) > 0) + if has_transmission && method.two_phase_gap + set_gap(1e-2) + large_gap = true end - @info @sprintf( - "Setting MILP time limit to %.2f seconds", - time_remaining - ) - JuMP.set_time_limit_sec(model, time_remaining) - @info "Solving MILP..." JuMP.optimize!(model) - has_transmission || break - violations = _find_violations( - model, - max_per_line = method.max_violations_per_line, - max_per_period = method.max_violations_per_period, - ) - if isempty(violations) - @info "No violations found" - if large_gap - large_gap = false - set_gap(method.gap_limit) - else + while true + initial_time = time() + time_elapsed = time() - initial_time + time_remaining = method.time_limit - time_elapsed + if time_remaining < 0 + @info "Time limit exceeded" break end - else - _enforce_transmission(model, violations) + @info @sprintf( + "Setting MILP time limit to %.2f seconds", + time_remaining + ) + JuMP.set_time_limit_sec(model, time_remaining) + @info "Solving MILP..." + JuMP.optimize!(model) + has_transmission || break + violations = _find_violations( + model, + scenario, + max_per_line = method.max_violations_per_line, + max_per_period = method.max_violations_per_period, + ) + if isempty(violations) + @info "No violations found" + if large_gap + large_gap = false + set_gap(method.gap_limit) + else + break + end + else + _enforce_transmission(model, violations, scenario) + end end end return diff --git a/src/solution/solution.jl b/src/solution/solution.jl index 0b76d3a..f02fb5b 100644 --- a/src/solution/solution.jl +++ b/src/solution/solution.jl @@ -16,34 +16,40 @@ solution = UnitCommitment.solution(model) """ function solution(model::JuMP.Model)::OrderedDict instance, T = model[:instance], model[:instance].time - function timeseries(vars, collection) + function timeseries_first_stage(vars, collection) return OrderedDict( b.name => [round(value(vars[b.name, t]), digits = 5) for t in 1:T] for b in collection ) end - function production_cost(g) + function timeseries_second_stage(vars, collection, sc) + return OrderedDict( + b.name => [round(value(vars[sc.name, b.name, t]), digits = 5) for t in 1:T] + for b in collection + ) + end + function production_cost(g, sc) return [ value(model[:is_on][g.name, t]) * g.min_power_cost[t] + sum( Float64[ - value(model[:segprod][g.name, t, k]) * + value(model[:segprod][sc.name, g.name, t, k]) * g.cost_segments[k].cost[t] for k in 1:length(g.cost_segments) ], ) for t in 1:T ] end - function production(g) + function production(g, sc) return [ value(model[:is_on][g.name, t]) * g.min_power[t] + sum( Float64[ - value(model[:segprod][g.name, t, k]) for + value(model[:segprod][sc.name, g.name, t, k]) for k in 1:length(g.cost_segments) ], ) for t in 1:T ] end - function startup_cost(g) + function startup_cost(g, sc) S = length(g.startup_categories) return [ sum( @@ -53,66 +59,70 @@ function solution(model::JuMP.Model)::OrderedDict ] end sol = OrderedDict() - sol["Production (MW)"] = - OrderedDict(g.name => production(g) for g in instance.units) - sol["Production cost (\$)"] = - OrderedDict(g.name => production_cost(g) for g in instance.units) - sol["Startup cost (\$)"] = - OrderedDict(g.name => startup_cost(g) for g in instance.units) - sol["Is on"] = timeseries(model[:is_on], instance.units) - sol["Switch on"] = timeseries(model[:switch_on], instance.units) - sol["Switch off"] = timeseries(model[:switch_off], instance.units) - sol["Net injection (MW)"] = - timeseries(model[:net_injection], instance.buses) - sol["Load curtail (MW)"] = timeseries(model[:curtail], instance.buses) - if !isempty(instance.lines) - sol["Line overflow (MW)"] = timeseries(model[:overflow], instance.lines) - end - if !isempty(instance.price_sensitive_loads) - sol["Price-sensitive loads (MW)"] = - timeseries(model[:loads], instance.price_sensitive_loads) - end - sol["Spinning reserve (MW)"] = OrderedDict( - r.name => OrderedDict( - g.name => [ - value(model[:reserve][r.name, g.name, t]) for + for sc in instance.scenarios + sol[sc.name] = OrderedDict() + sol[sc.name]["Production (MW)"] = + OrderedDict(g.name => production(g, sc) for g in sc.units) + sol[sc.name]["Production cost (\$)"] = + OrderedDict(g.name => production_cost(g, sc) for g in sc.units) + sol[sc.name]["Startup cost (\$)"] = + OrderedDict(g.name => startup_cost(g, sc) for g in sc.units) + sol[sc.name]["Is on"] = timeseries_first_stage(model[:is_on], sc.units) + sol[sc.name]["Switch on"] = timeseries_first_stage(model[:switch_on], sc.units) + sol[sc.name]["Switch off"] = timeseries_first_stage(model[:switch_off], sc.units) + sol[sc.name]["Net injection (MW)"] = + timeseries_second_stage(model[:net_injection], sc.buses, sc) + sol[sc.name]["Load curtail (MW)"] = timeseries_second_stage(model[:curtail], sc.buses, sc) + if !isempty(sc.lines) + sol[sc.name]["Line overflow (MW)"] = timeseries_second_stage(model[:overflow], sc.lines, sc) + end + if !isempty(sc.price_sensitive_loads) + sol[sc.name]["Price-sensitive loads (MW)"] = + timeseries_second_stage(model[:loads], sc.price_sensitive_loads, sc) + end + sol[sc.name]["Spinning reserve (MW)"] = OrderedDict( + r.name => OrderedDict( + g.name => [ + value(model[:reserve][sc.name, r.name, g.name, t]) for + t in 1:instance.time + ] for g in r.units + ) for r in sc.reserves if r.type == "spinning" + ) + sol[sc.name]["Spinning reserve shortfall (MW)"] = OrderedDict( + r.name => [ + value(model[:reserve_shortfall][sc.name, r.name, t]) for t in 1:instance.time - ] for g in r.units - ) for r in instance.reserves if r.type == "spinning" - ) - sol["Spinning reserve shortfall (MW)"] = OrderedDict( - r.name => [ - value(model[:reserve_shortfall][r.name, t]) for - t in 1:instance.time - ] for r in instance.reserves if r.type == "spinning" - ) - sol["Up-flexiramp (MW)"] = OrderedDict( - r.name => OrderedDict( - g.name => [ - value(model[:upflexiramp][r.name, g.name, t]) for + ] for r in sc.reserves if r.type == "spinning" + ) + sol[sc.name]["Up-flexiramp (MW)"] = OrderedDict( + r.name => OrderedDict( + g.name => [ + value(model[:upflexiramp][sc.name, r.name, g.name, t]) for + t in 1:instance.time + ] for g in r.units + ) for r in sc.reserves if r.type == "up-frp" + ) + sol[sc.name]["Up-flexiramp shortfall (MW)"] = OrderedDict( + r.name => [ + value(model[:upflexiramp_shortfall][sc.name, r.name, t]) for t in 1:instance.time - ] for g in r.units - ) for r in instance.reserves if r.type == "flexiramp" - ) - sol["Up-flexiramp shortfall (MW)"] = OrderedDict( - r.name => [ - value(model[:upflexiramp_shortfall][r.name, t]) for - t in 1:instance.time - ] for r in instance.reserves if r.type == "flexiramp" - ) - sol["Down-flexiramp (MW)"] = OrderedDict( - r.name => OrderedDict( - g.name => [ - value(model[:dwflexiramp][r.name, g.name, t]) for + ] for r in sc.reserves if r.type == "up-frp" + ) + sol[sc.name]["Down-flexiramp (MW)"] = OrderedDict( + r.name => OrderedDict( + g.name => [ + value(model[:dwflexiramp][sc.name, r.name, g.name, t]) for + t in 1:instance.time + ] for g in r.units + ) for r in sc.reserves if r.type == "down-frp" + ) + sol[sc.name]["Down-flexiramp shortfall (MW)"] = OrderedDict( + r.name => [ + value(model[:dwflexiramp_shortfall][sc.name, r.name, t]) for t in 1:instance.time - ] for g in r.units - ) for r in instance.reserves if r.type == "flexiramp" - ) - sol["Down-flexiramp shortfall (MW)"] = OrderedDict( - r.name => [ - value(model[:upflexiramp_shortfall][r.name, t]) for - t in 1:instance.time - ] for r in instance.reserves if r.type == "flexiramp" - ) + ] for r in sc.reserves if r.type == "down-frp" + ) + end + length(keys(sol)) > 1 ? nothing : sol = Dict(sol_key => sol_val for scen_key in keys(sol) for (sol_key, sol_val) in sol[scen_key]) return sol end diff --git a/src/validation/repair.jl b/src/validation/repair.jl index e1dd964..50ee614 100644 --- a/src/validation/repair.jl +++ b/src/validation/repair.jl @@ -3,19 +3,19 @@ # Released under the modified BSD license. See COPYING.md for more details. """ - repair!(instance) + repair!(sc) -Verifies that the given unit commitment instance is valid and automatically +Verifies that the given unit commitment scenario is valid and automatically fixes some validation errors if possible, issuing a warning for each error found. If a validation error cannot be automatically fixed, issues an exception. Returns the number of validation errors found. """ -function repair!(instance::UnitCommitmentInstance)::Int +function repair!(sc::UnitCommitmentScenario)::Int n_errors = 0 - for g in instance.units + for g in sc.units # Startup costs and delays must be increasing for s in 2:length(g.startup_categories) @@ -38,7 +38,7 @@ function repair!(instance::UnitCommitmentInstance)::Int end end - for t in 1:instance.time + for t in 1:sc.time # Production cost curve should be convex for k in 2:length(g.cost_segments) cost = g.cost_segments[k].cost[t] diff --git a/src/validation/validate.jl b/src/validation/validate.jl index 66a580f..9677661 100644 --- a/src/validation/validate.jl +++ b/src/validation/validate.jl @@ -356,7 +356,7 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01) required, ) end - elseif r.type == "flexiramp" + elseif r.type == "up-frp" upflexiramp = sum( solution["Up-flexiramp (MW)"][r.name][g.name][t] for g in r.units @@ -374,7 +374,7 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01) ) err_count += 1 end - + elseif r.type == "down-frp" dwflexiramp = sum( solution["Down-flexiramp (MW)"][r.name][g.name][t] for g in r.units