Set up multi-formulation architecture; start merging akazachk's code

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])

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).<formulation.isf_cutoff] .= 0
+        lodf[abs.(lodf).<formulation.lodf_cutoff] .= 0
+    end
+    model[:isf] = isf
+    model[:lodf] = lodf
+    return
+end

src/model/formulations/base/psload.jl

@@ -6,8 +6,8 @@ function _add_price_sensitive_load!(
     model::JuMP.Model,
     ps::PriceSensitiveLoad,
 )::Nothing
-    loads = _get(model, :loads)
-    net_injection = _get(model, :expr_net_injection)
+    loads = _init(model, :loads)
+    net_injection = _init(model, :expr_net_injection)
     for t in 1:model[:instance].time
         # Decision variable
         loads[ps.name, t] =

src/model/formulations/base/sensitivity.jl

@@ -0,0 +1,83 @@
+# 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.
+
+using SparseArrays, Base.Threads, LinearAlgebra, JuMP
+
+"""
+    _injection_shift_factors(; buses, lines)
+
+Returns a (B-1)xL matrix M, where B is the number of buses and L is the number
+of transmission lines. For a given bus b and transmission line l, the entry
+M[l.offset, b.offset] indicates the amount of power (in MW) that flows through
+transmission line l when 1 MW of power is injected at the slack bus (the bus
+that has offset zero) and withdrawn from b.
+"""
+function _injection_shift_factors(;
+    buses::Array{Bus},
+    lines::Array{TransmissionLine},
+)
+    susceptance = _susceptance_matrix(lines)
+    incidence = _reduced_incidence_matrix(lines = lines, buses = buses)
+    laplacian = transpose(incidence) * susceptance * incidence
+    isf = susceptance * incidence * inv(Array(laplacian))
+    return isf
+end
+
+"""
+    _reduced_incidence_matrix(; buses::Array{Bus}, lines::Array{TransmissionLine})
+
+Returns the incidence matrix for the network, with the column corresponding to
+the slack bus is removed. More precisely, returns a (B-1) x L matrix, where B
+is the number of buses and L is the number of lines. For each row, there is a 1
+element and a -1 element, indicating the source and target buses, respectively,
+for that line.
+"""
+function _reduced_incidence_matrix(;
+    buses::Array{Bus},
+    lines::Array{TransmissionLine},
+)
+    matrix = spzeros(Float64, length(lines), length(buses) - 1)
+    for line in lines
+        if line.source.offset > 0
+            matrix[line.offset, line.source.offset] = 1
+        end
+        if line.target.offset > 0
+            matrix[line.offset, line.target.offset] = -1
+        end
+    end
+    return matrix
+end
+
+"""
+    _susceptance_matrix(lines::Array{TransmissionLine})
+
+Returns a LxL diagonal matrix, where each diagonal entry is the susceptance of
+the corresponding transmission line.
+"""
+function _susceptance_matrix(lines::Array{TransmissionLine})
+    return Diagonal([l.susceptance for l in lines])
+end
+
+"""
+
+    _line_outage_factors(; buses, lines, isf)
+
+Returns a LxL matrix containing the Line Outage Distribution Factors (LODFs)
+for the given network. This matrix how does the pre-contingency flow change
+when each individual transmission line is removed.
+"""
+function _line_outage_factors(;
+    buses::Array{Bus,1},
+    lines::Array{TransmissionLine,1},
+    isf::Array{Float64,2},
+)::Array{Float64,2}
+    incidence = Array(_reduced_incidence_matrix(lines = lines, buses = buses))
+    lodf::Array{Float64,2} = isf * transpose(incidence)
+    _, n = size(lodf)
+    for i in 1:n
+        lodf[:, i] *= 1.0 / (1.0 - lodf[i, i])
+        lodf[i, i] = -1
+    end
+    return lodf
+end

src/model/formulations/base/structs.jl

@@ -0,0 +1,57 @@
+# 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 _RampingFormulation end
+abstract type _TransmissionFormulation end
+
+struct _GeneratorFormulation
+    ramping::_RampingFormulation
+
+    function _GeneratorFormulation(
+        ramping::_RampingFormulation = _DamKucRajAta16(),
+    )
+        return new(ramping)
+    end
+end
+
+"""
+    mutable struct _ShiftFactorsFormulation <: _TransmissionFormulation
+        isf_cutoff::Float64
+        lodf_cutoff::Float64
+        precomputed_isf::Union{Nothing,Matrix{Float64}}
+        precomputed_lodf::Union{Nothing,Matrix{Float64}}
+    end
+
+Transmission formulation based on Injection Shift Factors (ISF) and Line
+Outage Distribution Factors (LODF). Constraints are enforced in a lazy way.
+
+Arguments
+---------
+- `precomputed_isf::Union{Matrix{Float64},Nothing} = nothing`:
+    the injection shift factors matrix. If not provided, it will be computed.
+- `precomputed_lodf::Union{Matrix{Float64},Nothing} = nothing`: 
+    the line outage distribution factors matrix. If not provided, it will be
+    computed.
+- `isf_cutoff::Float64 = 0.005`: 
+    the cutoff that should be applied to the ISF matrix. Entries with magnitude
+    smaller than this value will be set to zero.
+- `lodf_cutoff::Float64 = 0.001`: 
+    the cutoff that should be applied to the LODF matrix. Entries with magnitude
+    smaller than this value will be set to zero.
+"""
+mutable struct _ShiftFactorsFormulation <: _TransmissionFormulation
+    isf_cutoff::Float64
+    lodf_cutoff::Float64
+    precomputed_isf::Union{Nothing,Matrix{Float64}}
+    precomputed_lodf::Union{Nothing,Matrix{Float64}}
+
+    function _ShiftFactorsFormulation(;
+        isf_cutoff = 0.005,
+        lodf_cutoff = 0.001,
+        precomputed_isf = nothing,
+        precomputed_lodf = nothing,
+    )
+        return new(isf_cutoff, lodf_cutoff, precomputed_isf, precomputed_lodf)
+    end
+end

src/model/formulations/base/system.jl

@@ -10,9 +10,9 @@ end
 
 function _add_net_injection_eqs!(model::JuMP.Model)::Nothing
     T = model[:instance].time
-    net_injection = _get(model, :net_injection)
-    eq_net_injection_def = _get(model, :eq_net_injection_def)
-    eq_power_balance = _get(model, :eq_power_balance)
+    net_injection = _init(model, :net_injection)
+    eq_net_injection_def = _init(model, :eq_net_injection_def)
+    eq_power_balance = _init(model, :eq_power_balance)
     for t in 1:T, b in model[:instance].buses
         n = net_injection[b.name, t] = @variable(model)
         eq_net_injection_def[t, b.name] =
@@ -28,7 +28,7 @@ function _add_net_injection_eqs!(model::JuMP.Model)::Nothing
 end
 
 function _add_reserve_eqs!(model::JuMP.Model)::Nothing
-    eq_min_reserve = _get(model, :eq_min_reserve)
+    eq_min_reserve = _init(model, :eq_min_reserve)
     for t in 1:model[:instance].time
         eq_min_reserve[t] = @constraint(
             model,

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)
+function _add_unit!(model::JuMP.Model, g::Unit, f::_GeneratorFormulation)
     if !all(g.must_run) && any(g.must_run)
         error("Partially must-run units are not currently supported")
     end
@@ -11,26 +11,31 @@ function _add_unit!(model::JuMP.Model, g::Unit)
     end
 
     # Variables
-    _add_production_vars!(model, g)
-    _add_reserve_vars!(model, g)
-    _add_startup_shutdown_vars!(model, g)
-    _add_status_vars!(model, g)
+    _add_production_vars!(model, g, f)
+    _add_reserve_vars!(model, g, f)
+    _add_startup_shutdown_vars!(model, g, f)
+    _add_status_vars!(model, g, f)
 
     # Constraints and objective function
-    _add_min_uptime_downtime_eqs!(model, g)
-    _add_net_injection_eqs!(model, g)
-    _add_production_eqs!(model, g)
-    _add_ramp_eqs!(model, g)
-    _add_startup_shutdown_costs_eqs!(model, g)
-    _add_startup_shutdown_limit_eqs!(model, g)
-    return _add_status_eqs!(model, g)
+    _add_min_uptime_downtime_eqs!(model, g, f)
+    _add_net_injection_eqs!(model, g, f)
+    _add_production_eqs!(model, g, f)
+    _add_ramp_eqs!(model, g, f.ramping)
+    _add_startup_shutdown_costs_eqs!(model, g, f)
+    _add_startup_shutdown_limit_eqs!(model, g, f)
+    _add_status_eqs!(model, g, f)
+    return
 end
 
 _is_initially_on(g::Unit)::Float64 = (g.initial_status > 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