From bf6d19343e27aa9e6ab7e7daea11f22fb6b64c84 Mon Sep 17 00:00:00 2001 From: Alinson S Xavier Date: Sun, 30 May 2021 07:14:28 -0500 Subject: [PATCH] Set up multi-formulation architecture; start merging akazachk's code --- src/UnitCommitment.jl | 7 +- src/model/build.jl | 69 ++++---- src/model/formulations/ArrCon00/structs.jl | 12 ++ src/model/formulations/ArrCon00/unit.jl | 92 +++++++++++ .../formulations/DamKucRajAta16/structs.jl | 11 ++ src/model/formulations/DamKucRajAta16/unit.jl | 116 ++++++++++++++ src/model/formulations/base/bus.jl | 6 +- src/model/formulations/base/line.jl | 57 ++++++- src/model/formulations/base/psload.jl | 4 +- .../formulations/base}/sensitivity.jl | 3 +- src/model/formulations/base/structs.jl | 57 +++++++ src/model/formulations/base/system.jl | 8 +- src/model/formulations/base/unit.jl | 151 ++++++++++++------ src/model/jumpext.jl | 2 +- src/model/structs.jl | 5 - 15 files changed, 486 insertions(+), 114 deletions(-) create mode 100644 src/model/formulations/ArrCon00/structs.jl create mode 100644 src/model/formulations/ArrCon00/unit.jl create mode 100644 src/model/formulations/DamKucRajAta16/structs.jl create mode 100644 src/model/formulations/DamKucRajAta16/unit.jl rename src/{transmission => model/formulations/base}/sensitivity.jl (98%) create mode 100644 src/model/formulations/base/structs.jl delete mode 100644 src/model/structs.jl diff --git a/src/UnitCommitment.jl b/src/UnitCommitment.jl index e7692ae..c8e6765 100644 --- a/src/UnitCommitment.jl +++ b/src/UnitCommitment.jl @@ -6,6 +6,9 @@ module UnitCommitment include("instance/structs.jl") include("solution/structs.jl") +include("model/formulations/base/structs.jl") +include("model/formulations/ArrCon00/structs.jl") +include("model/formulations/DamKucRajAta16/structs.jl") include("solution/methods/XaQiWaTh19/structs.jl") include("import/egret.jl") @@ -16,6 +19,9 @@ include("model/formulations/base/line.jl") include("model/formulations/base/psload.jl") include("model/formulations/base/system.jl") include("model/formulations/base/unit.jl") +include("model/formulations/base/sensitivity.jl") +include("model/formulations/DamKucRajAta16/unit.jl") +include("model/formulations/ArrCon00/unit.jl") include("model/jumpext.jl") include("solution/fix.jl") include("solution/methods/XaQiWaTh19/enforce.jl") @@ -28,7 +34,6 @@ include("solution/warmstart.jl") include("solution/write.jl") include("transform/initcond.jl") include("transform/slice.jl") -include("transmission/sensitivity.jl") include("utils/log.jl") include("validation/repair.jl") include("validation/validate.jl") diff --git a/src/model/build.jl b/src/model/build.jl index 9ff8fe4..7e7259c 100644 --- a/src/model/build.jl +++ b/src/model/build.jl @@ -55,45 +55,25 @@ julia> model = UnitCommitment.build_model( """ function build_model(; instance::UnitCommitmentInstance, - isf::Union{Matrix{Float64},Nothing} = nothing, - lodf::Union{Matrix{Float64},Nothing} = nothing, - isf_cutoff::Float64 = 0.005, - lodf_cutoff::Float64 = 0.001, optimizer = nothing, variable_names::Bool = false, )::JuMP.Model - if length(instance.buses) == 1 - isf = zeros(0, 0) - lodf = zeros(0, 0) - else - if isf === nothing - @info "Computing injection shift factors..." - time_isf = @elapsed begin - isf = UnitCommitment._injection_shift_factors( - lines = instance.lines, - buses = instance.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, - isf = isf, - ) - end - @info @sprintf("Computed LODF in %.2f seconds", time_lodf) + return _build_model( + instance, + _GeneratorFormulation(), + _ShiftFactorsFormulation(), + optimizer = optimizer, + variable_names = variable_names, + ) +end - @info @sprintf( - "Applying PTDF and LODF cutoffs (%.5f, %.5f)", - isf_cutoff, - lodf_cutoff - ) - isf[abs.(isf). 0) + + for t in 1:model[:instance].time + # Ramp up limit + 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( + mip, + g.min_power[t] + + prod_above[gn, t] + + (RESERVES_WHEN_RAMP_UP ? reserve[gn, t] : 0.0) <= + g.initial_power + RU + ) + end + else + max_prod_this_period = + g.min_power[t] * is_on[gn, t] + + prod_above[gn, t] + + ( + RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? + reserve[gn, t] : 0.0 + ) + min_prod_last_period = + g.min_power[t-1] * is_on[gn, t-1] + prod_above[gn, t-1] + + # Equation (24) in Kneuven et al. (2020) + eq_ramp_up[gn, t] = @constraint( + mip, + max_prod_this_period - min_prod_last_period <= + RU * is_on[gn, t-1] + SU * switch_on[gn, t] + ) + end + + # Ramp down limit + if t == 1 + if is_initially_on + # TODO If RD < SD, or more specifically if + # 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( + mip, + g.initial_power - (g.min_power[t] + prod_above[gn, t]) <= RD + ) + end + else + max_prod_last_period = + g.min_power[t-1] * is_on[gn, t-1] + + prod_above[gn, t-1] + + ( + RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ? + reserve[gn, t-1] : 0.0 + ) + min_prod_this_period = + g.min_power[t] * is_on[gn, t] + prod_above[gn, t] + + # Equation (25) in Kneuven et al. (2020) + eq_ramp_down[gn, t] = @constraint( + mip, + max_prod_last_period - min_prod_this_period <= + RD * is_on[gn, t] + SD * switch_off[gn, t] + ) + end + end +end diff --git a/src/model/formulations/DamKucRajAta16/structs.jl b/src/model/formulations/DamKucRajAta16/structs.jl new file mode 100644 index 0000000..1777691 --- /dev/null +++ b/src/model/formulations/DamKucRajAta16/structs.jl @@ -0,0 +1,11 @@ +# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +""" +Formulation described in: + + Damcı-Kurt, P., Küçükyavuz, S., Rajan, D., & Atamtürk, A. (2016). A polyhedral + study of production ramping. Mathematical Programming, 158(1), 175-205. +""" +mutable struct _DamKucRajAta16 <: _RampingFormulation end diff --git a/src/model/formulations/DamKucRajAta16/unit.jl b/src/model/formulations/DamKucRajAta16/unit.jl new file mode 100644 index 0000000..64a851e --- /dev/null +++ b/src/model/formulations/DamKucRajAta16/unit.jl @@ -0,0 +1,116 @@ +# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +function _add_ramp_eqs!( + model::JuMP.Model, + g::Unit, + formulation::_DamKucRajAta16, +)::Nothing + # TODO: Move upper case constants to model[:instance] + RESERVES_WHEN_START_UP = true + RESERVES_WHEN_RAMP_UP = true + RESERVES_WHEN_RAMP_DOWN = true + RESERVES_WHEN_SHUT_DOWN = true + known_initial_conditions = true + is_initially_on = (g.initial_status > 0) + SU = g.startup_limit + SD = g.shutdown_limit + RU = g.ramp_up_limit + RD = g.ramp_down_limit + gn = g.name + eq_str_ramp_down = _init(model, :eq_str_ramp_down) + eq_str_ramp_up = _init(model, :eq_str_ramp_up) + is_on = model[:is_on] + prod_above = model[:prod_above] + reserve = model[:reserve] + switch_off = model[:switch_off] + switch_on = model[:switch_on] + + for t in 1:model[:instance].time + time_invariant = + (t > 1) ? (abs(g.min_power[t] - g.min_power[t-1]) < 1e-7) : true + + if t > 1 && !time_invariant + warning( + "Ramping according to Damcı-Kurt et al. (2016) requires " * + "time-invariant minimum power. This does not hold for " * + "generator $(gn): min_power[$t] = $(g.min_power[t]); " * + "min_power[$(t-1)] = $(g.min_power[t-1]). Reverting to " * + "Arroyo and Conejo (2000) formulation for this generator.", + ) + end + + max_prod_this_period = + prod_above[gn, t] + ( + RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? + reserve[gn, t] : 0.0 + ) + min_prod_last_period = 0.0 + if t > 1 && time_invariant + min_prod_last_period = prod_above[gn, t-1] + + # Equation (35) in Kneuven et al. (2020) + # Sparser version of (24) + eq_str_ramp_up[gn, t] = @constraint( + model, + max_prod_this_period - min_prod_last_period <= + (SU - g.min_power[t] - RU) * switch_on[gn, t] + + RU * is_on[gn, t] + ) + 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] + else + min_prod_last_period = max(g.initial_power, 0.0) + end + + # Add the min prod at time t back in to max_prod_this_period to get _total_ production + # (instead of using the amount above minimum, as min prod for t < 1 is unknown) + max_prod_this_period += g.min_power[t] * is_on[gn, t] + + # Modified version of equation (35) in Kneuven et al. (2020) + # Equivalent to (24) + eq_str_ramp_up[gn, t] = @constraint( + model, + max_prod_this_period - min_prod_last_period <= + (SU - RU) * switch_on[gn, t] + RU * is_on[gn, t] + ) + end + + max_prod_last_period = + min_prod_last_period + ( + t > 1 && (RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN) ? + reserve[gn, t-1] : 0.0 + ) + min_prod_this_period = prod_above[gn, t] + on_last_period = 0.0 + if t > 1 + on_last_period = is_on[gn, t-1] + elseif (known_initial_conditions && g.initial_status > 0) + on_last_period = 1.0 + end + + if t > 1 && time_invariant + # Equation (36) in Kneuven et al. (2020) + eq_str_ramp_down[gn, t] = @constraint( + model, + max_prod_last_period - min_prod_this_period <= + (SD - g.min_power[t] - RD) * switch_off[gn, t] + + RD * on_last_period + ) + elseif (t == 1 && is_initially_on) || (t > 1 && !time_invariant) + # Add back in min power + min_prod_this_period += g.min_power[t] * is_on[gn, t] + + # Modified version of equation (36) in Kneuven et al. (2020) + # Equivalent to (25) + eq_str_ramp_down[gn, t] = @constraint( + model, + max_prod_last_period - min_prod_this_period <= + (SD - RD) * switch_off[gn, t] + RD * on_last_period + ) + end + end +end diff --git a/src/model/formulations/base/bus.jl b/src/model/formulations/base/bus.jl index cc176ca..00cacba 100644 --- a/src/model/formulations/base/bus.jl +++ b/src/model/formulations/base/bus.jl @@ -3,9 +3,9 @@ # Released under the modified BSD license. See COPYING.md for more details. function _add_bus!(model::JuMP.Model, b::Bus)::Nothing - net_injection = _get(model, :expr_net_injection) - reserve = _get(model, :expr_reserve) - curtail = _get(model, :curtail) + net_injection = _init(model, :expr_net_injection) + reserve = _init(model, :expr_reserve) + curtail = _init(model, :curtail) for t in 1:model[:instance].time # Fixed load net_injection[b.name, t] = AffExpr(-b.load[t]) diff --git a/src/model/formulations/base/line.jl b/src/model/formulations/base/line.jl index 8370336..6501982 100644 --- a/src/model/formulations/base/line.jl +++ b/src/model/formulations/base/line.jl @@ -2,11 +2,60 @@ # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. # Released under the modified BSD license. See COPYING.md for more details. -function _add_transmission_line!(model, lm)::Nothing - overflow = _get(model, :overflow) +function _add_transmission_line!( + model::JuMP.Model, + lm::TransmissionLine, + f::_TransmissionFormulation, +)::Nothing + overflow = _init(model, :overflow) for t in 1:model[:instance].time - v = overflow[lm.name, t] = @variable(model, lower_bound = 0) - add_to_expression!(model[:obj], v, lm.flow_limit_penalty[t]) + overflow[lm.name, t] = @variable(model, lower_bound = 0) + add_to_expression!( + model[:obj], + overflow[lm.name, t], + lm.flow_limit_penalty[t], + ) end return end + +function _setup_transmission( + model::JuMP.Model, + formulation::_TransmissionFormulation, +)::Nothing + instance = model[:instance] + isf = formulation.precomputed_isf + lodf = formulation.precomputed_lodf + if length(instance.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, + ) + 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, + isf = isf, + ) + end + @info @sprintf("Computed LODF in %.2f seconds", time_lodf) + @info @sprintf( + "Applying PTDF and LODF cutoffs (%.5f, %.5f)", + formulation.isf_cutoff, + formulation.lodf_cutoff + ) + isf[abs.(isf). 0 ? 1.0 : 0.0) -function _add_production_vars!(model::JuMP.Model, g::Unit)::Nothing - prod_above = _get(model, :prod_above) - segprod = _get(model, :segprod) +function _add_production_vars!( + model::JuMP.Model, + g::Unit, + formulation::_GeneratorFormulation, +)::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) @@ -40,10 +45,14 @@ function _add_production_vars!(model::JuMP.Model, g::Unit)::Nothing return end -function _add_production_eqs!(model::JuMP.Model, g::Unit)::Nothing - eq_prod_above_def = _get(model, :eq_prod_above_def) - eq_prod_limit = _get(model, :eq_prod_limit) - eq_segprod_limit = _get(model, :eq_segprod_limit) +function _add_production_eqs!( + model::JuMP.Model, + g::Unit, + formulation::_GeneratorFormulation, +)::Nothing + eq_prod_above_def = _init(model, :eq_prod_above_def) + eq_prod_limit = _init(model, :eq_prod_limit) + eq_segprod_limit = _init(model, :eq_segprod_limit) is_on = model[:is_on] K = length(g.cost_segments) prod_above = model[:prod_above] @@ -82,8 +91,12 @@ function _add_production_eqs!(model::JuMP.Model, g::Unit)::Nothing return end -function _add_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing - reserve = _get(model, :reserve) +function _add_reserve_vars!( + model::JuMP.Model, + g::Unit, + formulation::_GeneratorFormulation, +)::Nothing + reserve = _init(model, :reserve) for t in 1:model[:instance].time if g.provides_spinning_reserves[t] reserve[g.name, t] = @variable(model, lower_bound = 0) @@ -94,7 +107,11 @@ function _add_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing return end -function _add_reserve_eqs!(model::JuMP.Model, g::Unit)::Nothing +function _add_reserve_eqs!( + model::JuMP.Model, + g::Unit, + formulation::_GeneratorFormulation, +)::Nothing reserve = model[:reserve] for t in 1:model[:instance].time add_to_expression!(expr_reserve[g.bus.name, t], reserve[g.name, t], 1.0) @@ -102,8 +119,12 @@ function _add_reserve_eqs!(model::JuMP.Model, g::Unit)::Nothing return end -function _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing - startup = _get(model, :startup) +function _add_startup_shutdown_vars!( + model::JuMP.Model, + g::Unit, + formulation::_GeneratorFormulation, +)::Nothing + startup = _init(model, :startup) for t in 1:model[:instance].time for s in 1:length(g.startup_categories) startup[g.name, t, s] = @variable(model, binary = true) @@ -112,9 +133,13 @@ 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 - eq_shutdown_limit = _get(model, :eq_shutdown_limit) - eq_startup_limit = _get(model, :eq_startup_limit) +function _add_startup_shutdown_limit_eqs!( + model::JuMP.Model, + g::Unit, + formulation::_GeneratorFormulation, +)::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 = model[:reserve] @@ -147,9 +172,13 @@ function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing return end -function _add_startup_shutdown_costs_eqs!(model::JuMP.Model, g::Unit)::Nothing - eq_startup_choose = _get(model, :eq_startup_choose) - eq_startup_restrict = _get(model, :eq_startup_restrict) +function _add_startup_shutdown_costs_eqs!( + model::JuMP.Model, + g::Unit, + formulation::_GeneratorFormulation, +)::Nothing + eq_startup_choose = _init(model, :eq_startup_choose) + eq_startup_restrict = _init(model, :eq_startup_restrict) S = length(g.startup_categories) startup = model[:startup] for t in 1:model[:instance].time @@ -190,10 +219,14 @@ function _add_startup_shutdown_costs_eqs!(model::JuMP.Model, g::Unit)::Nothing return end -function _add_status_vars!(model::JuMP.Model, g::Unit)::Nothing - is_on = _get(model, :is_on) - switch_on = _get(model, :switch_on) - switch_off = _get(model, :switch_off) +function _add_status_vars!( + model::JuMP.Model, + g::Unit, + formulation::_GeneratorFormulation, +)::Nothing + is_on = _init(model, :is_on) + switch_on = _init(model, :switch_on) + switch_off = _init(model, :switch_off) for t in 1:model[:instance].time if g.must_run[t] is_on[g.name, t] = 1.0 @@ -208,9 +241,13 @@ function _add_status_vars!(model::JuMP.Model, g::Unit)::Nothing return end -function _add_status_eqs!(model::JuMP.Model, g::Unit)::Nothing - eq_binary_link = _get(model, :eq_binary_link) - eq_switch_on_off = _get(model, :eq_switch_on_off) +function _add_status_eqs!( + model::JuMP.Model, + g::Unit, + formulation::_GeneratorFormulation, +)::Nothing + eq_binary_link = _init(model, :eq_binary_link) + eq_switch_on_off = _init(model, :eq_switch_on_off) is_on = model[:is_on] switch_off = model[:switch_off] switch_on = model[:switch_on] @@ -240,11 +277,15 @@ function _add_status_eqs!(model::JuMP.Model, g::Unit)::Nothing return end -function _add_ramp_eqs!(model::JuMP.Model, g::Unit)::Nothing +function _add_ramp_eqs!( + model::JuMP.Model, + g::Unit, + formulation::_GeneratorFormulation, +)::Nothing prod_above = model[:prod_above] reserve = model[:reserve] - eq_ramp_up = _get(model, :eq_ramp_up) - eq_ramp_down = _get(model, :eq_ramp_down) + 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 @@ -282,12 +323,16 @@ function _add_ramp_eqs!(model::JuMP.Model, g::Unit)::Nothing end end -function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing +function _add_min_uptime_downtime_eqs!( + model::JuMP.Model, + g::Unit, + formulation::_GeneratorFormulation, +)::Nothing is_on = model[:is_on] switch_off = model[:switch_off] switch_on = model[:switch_on] - eq_min_uptime = _get(model, :eq_min_uptime) - eq_min_downtime = _get(model, :eq_min_downtime) + eq_min_uptime = _init(model, :eq_min_uptime) + eq_min_downtime = _init(model, :eq_min_downtime) T = model[:instance].time for t in 1:T # Minimum up-time @@ -325,25 +370,29 @@ 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, + formulation::_GeneratorFormulation, +)::Nothing expr_net_injection = model[:expr_net_injection] - expr_reserve = model[:expr_reserve] - is_on = model[:is_on] - prod_above = model[:prod_above] - reserve = model[:reserve] for t in 1:model[:instance].time # Add to net injection expression add_to_expression!( expr_net_injection[g.bus.name, t], - prod_above[g.name, t], + model[:prod_above][g.name, t], 1.0, ) add_to_expression!( expr_net_injection[g.bus.name, t], - is_on[g.name, t], + model[:is_on][g.name, t], g.min_power[t], ) # Add to reserves expression - add_to_expression!(expr_reserve[g.bus.name, t], reserve[g.name, t], 1.0) + add_to_expression!( + model[:expr_reserve][g.bus.name, t], + model[:reserve][g.name, t], + 1.0, + ) end end diff --git a/src/model/jumpext.jl b/src/model/jumpext.jl index 1586e92..be73692 100644 --- a/src/model/jumpext.jl +++ b/src/model/jumpext.jl @@ -19,7 +19,7 @@ function set_name(x::Float64, n::String) # nop end -function _get(model::JuMP.Model, key::Symbol)::OrderedDict +function _init(model::JuMP.Model, key::Symbol)::OrderedDict if !(key in keys(object_dictionary(model))) model[key] = OrderedDict() end diff --git a/src/model/structs.jl b/src/model/structs.jl deleted file mode 100644 index 567bec2..0000000 --- a/src/model/structs.jl +++ /dev/null @@ -1,5 +0,0 @@ -# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment -# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. -# Released under the modified BSD license. See COPYING.md for more details. - -abstract type Formulation end