Merge remote-tracking branch 'upstream/dev' into dev

src/UnitCommitment.jl

@@ -4,6 +4,8 @@
 
 module UnitCommitment
 
+using Base: String
+
 include("instance/structs.jl")
 include("model/formulations/base/structs.jl")
 include("solution/structs.jl")

src/instance/read.jl

@@ -43,26 +43,76 @@ function read_benchmark(
     return UnitCommitment.read(filename)
 end
 
+function _repair_scenario_names_and_probabilities!(
+    scenarios::Vector{UnitCommitmentScenario},
+    path::Vector{String},
+)::Nothing
+    total_weight = sum([sc.probability for sc in scenarios])
+    for (sc_path, sc) in zip(path, scenarios)
+        sc.name !== "" ||
+            (sc.name = first(split(last(split(sc_path, "/")), ".")))
+        sc.probability = (sc.probability / total_weight)
+    end
+    return
+end
+
 """
     read(path::AbstractString)::UnitCommitmentInstance
 
-Read instance from a file. The file may be gzipped.
+Read a deterministic test case from the given file. The file may be gzipped.
 
 # Example
 
 ```julia
-instance = UnitCommitment.read("/path/to/input.json.gz")
+instance = UnitCommitment.read("s1.json.gz")
 ```
 """
-function read(path::AbstractString)::UnitCommitmentInstance
+function read(path::String)::UnitCommitmentInstance
-    if endswith(path, ".gz")
+    scenarios = Vector{UnitCommitmentScenario}()
-        return _read(gzopen(path))
+    scenario = _read_scenario(path)
-    else
+    scenario.name = "s1"
-        return _read(open(path))
+    scenario.probability = 1.0
-    end
+    scenarios = [scenario]
+    instance =
+        UnitCommitmentInstance(time = scenario.time, scenarios = scenarios)
+    return instance
 end
 
-function _read(file::IO)::UnitCommitmentInstance
+"""
+    read(path::Vector{String})::UnitCommitmentInstance
+
+Read a stochastic unit commitment instance from the given files. Each file
+describes a scenario. The files may be gzipped.
+
+# Example
+
+```julia
+instance = UnitCommitment.read(["s1.json.gz", "s2.json.gz"])
+```
+"""
+function read(paths::Vector{String})::UnitCommitmentInstance
+    scenarios = UnitCommitmentScenario[]
+    for p in paths
+        push!(scenarios, _read_scenario(p))
+    end
+    _repair_scenario_names_and_probabilities!(scenarios, paths)
+    instance =
+        UnitCommitmentInstance(time = scenarios[1].time, scenarios = scenarios)
+    return instance
+end
+
+function _read_scenario(path::String)::UnitCommitmentScenario
+    if endswith(path, ".gz")
+        scenario = _read(gzopen(path))
+    elseif endswith(path, ".json")
+        scenario = _read(open(path))
+    else
+        error("Unsupported input format")
+    end
+    return scenario
+end
+
+function _read(file::IO)::UnitCommitmentScenario
     return _from_json(
         JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing)),
     )
@@ -77,7 +127,7 @@ function _read_json(path::String)::OrderedDict
     return JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing))
 end
 
-function _from_json(json; repair = true)
+function _from_json(json; repair = true)::UnitCommitmentScenario
     _migrate(json)
     units = Unit[]
     buses = Bus[]
@@ -102,6 +152,11 @@ function _from_json(json; repair = true)
     time_multiplier = 60 ÷ time_step
     T = time_horizon * time_multiplier
 
+    probability = json["Parameters"]["Scenario weight"]
+    probability !== nothing || (probability = 1)
+    scenario_name = json["Parameters"]["Scenario name"]
+    scenario_name !== nothing || (scenario_name = "")
+
     name_to_bus = Dict{String,Bus}()
     name_to_line = Dict{String,TransmissionLine}()
     name_to_unit = Dict{String,Unit}()
@@ -118,15 +173,6 @@ function _from_json(json; repair = true)
         json["Parameters"]["Power balance penalty (\$/MW)"],
         default = [1000.0 for t in 1:T],
     )
-    # Penalty price for shortage in meeting system-wide flexiramp requirements
-    flexiramp_shortfall_penalty = timeseries(
-        json["Parameters"]["Flexiramp penalty (\$/MW)"],
-        default = [500.0 for t in 1:T],
-    )
-    shortfall_penalty = timeseries(
-        json["Parameters"]["Reserve shortfall penalty (\$/MW)"],
-        default = [-1.0 for t in 1:T],
-    )
 
     # Read buses
     for (bus_name, dict) in json["Buses"]
@@ -308,7 +354,9 @@ function _from_json(json; repair = true)
         end
     end
 
-    instance = UnitCommitmentInstance(
+    scenario = UnitCommitmentScenario(
+        name = scenario_name,
+        probability = probability,
         buses_by_name = Dict(b.name => b for b in buses),
         buses = buses,
         contingencies_by_name = Dict(c.name => c for c in contingencies),
@@ -320,14 +368,14 @@ function _from_json(json; repair = true)
         price_sensitive_loads = loads,
         reserves = reserves,
         reserves_by_name = name_to_reserve,
-        shortfall_penalty = shortfall_penalty,
-        flexiramp_shortfall_penalty = flexiramp_shortfall_penalty,
         time = T,
         units_by_name = Dict(g.name => g for g in units),
         units = units,
+        isf = spzeros(Float64, length(lines), length(buses) - 1),
+        lodf = spzeros(Float64, length(lines), length(lines)),
     )
     if repair
-        UnitCommitment.repair!(instance)
+        UnitCommitment.repair!(scenario)
     end
-    return instance
+    return scenario
 end

src/instance/structs.jl

@@ -73,35 +73,41 @@ mutable struct PriceSensitiveLoad
     revenue::Vector{Float64}
 end
 
-Base.@kwdef mutable struct UnitCommitmentInstance
+Base.@kwdef mutable struct UnitCommitmentScenario
     buses_by_name::Dict{AbstractString,Bus}
     buses::Vector{Bus}
     contingencies_by_name::Dict{AbstractString,Contingency}
     contingencies::Vector{Contingency}
+    isf::Array{Float64,2}
     lines_by_name::Dict{AbstractString,TransmissionLine}
     lines::Vector{TransmissionLine}
+    lodf::Array{Float64,2}
+    name::String
     power_balance_penalty::Vector{Float64}
     price_sensitive_loads_by_name::Dict{AbstractString,PriceSensitiveLoad}
     price_sensitive_loads::Vector{PriceSensitiveLoad}
-    reserves::Vector{Reserve}
+    probability::Float64
     reserves_by_name::Dict{AbstractString,Reserve}
-    shortfall_penalty::Vector{Float64}
+    reserves::Vector{Reserve}
-    flexiramp_shortfall_penalty::Vector{Float64}
     time::Int
     units_by_name::Dict{AbstractString,Unit}
     units::Vector{Unit}
 end
 
+Base.@kwdef mutable struct UnitCommitmentInstance
+    time::Int
+    scenarios::Vector{UnitCommitmentScenario}
+end
+
 function Base.show(io::IO, instance::UnitCommitmentInstance)
+    sc = instance.scenarios[1]
     print(io, "UnitCommitmentInstance(")
-    print(io, "$(length(instance.units)) units, ")
+    print(io, "$(length(instance.scenarios)) scenarios, ")
-    print(io, "$(length(instance.buses)) buses, ")
+    print(io, "$(length(sc.units)) units, ")
-    print(io, "$(length(instance.lines)) lines, ")
+    print(io, "$(length(sc.buses)) buses, ")
-    print(io, "$(length(instance.contingencies)) contingencies, ")
+    print(io, "$(length(sc.lines)) lines, ")
-    print(
+    print(io, "$(length(sc.contingencies)) contingencies, ")
-        io,
+    print(io, "$(length(sc.price_sensitive_loads)) price sensitive loads, ")
-        "$(length(instance.price_sensitive_loads)) price sensitive loads, ",
-    )
     print(io, "$(instance.time) time steps")
     print(io, ")")
     return

src/model/build.jl

@@ -77,20 +77,27 @@ function build_model(;
         end
         model[:obj] = AffExpr()
         model[:instance] = instance
-        _setup_transmission(model, formulation.transmission)
+        for g in instance.scenarios[1].units
-        for l in instance.lines
+            _add_unit_commitment!(model, g, formulation)
-            _add_transmission_line!(model, l, formulation.transmission)
         end
-        for b in instance.buses
+        for sc in instance.scenarios
-            _add_bus!(model, b)
+            @info "Building scenario $(sc.name) with " *
+                  "probability $(sc.probability)"
+            _setup_transmission(formulation.transmission, sc)
+            for l in sc.lines
+                _add_transmission_line!(model, l, formulation.transmission, sc)
+            end
+            for b in sc.buses
+                _add_bus!(model, b, sc)
+            end
+            for ps in sc.price_sensitive_loads
+                _add_price_sensitive_load!(model, ps, sc)
+            end
+            for g in sc.units
+                _add_unit_dispatch!(model, g, formulation, sc)
+            end
+            _add_system_wide_eqs!(model, sc)
         end
-        for g in instance.units
-            _add_unit!(model, g, formulation)
-        end
-        for ps in instance.price_sensitive_loads
-            _add_price_sensitive_load!(model, ps)
-        end
-        _add_system_wide_eqs!(model)
         @objective(model, Min, model[:obj])
     end
     @info @sprintf("Built model in %.2f seconds", time_model)

src/model/formulations/ArrCon2000/ramp.jl

@@ -8,6 +8,7 @@ function _add_ramp_eqs!(
     formulation_prod_vars::Gar1962.ProdVars,
     formulation_ramping::ArrCon2000.Ramping,
     formulation_status_vars::Gar1962.StatusVars,
+    sc::UnitCommitmentScenario,
 )::Nothing
     # TODO: Move upper case constants to model[:instance]
     RESERVES_WHEN_START_UP = true
@@ -22,7 +23,7 @@ function _add_ramp_eqs!(
     eq_ramp_down = _init(model, :eq_ramp_down)
     eq_ramp_up = _init(model, :eq_ramp_up)
     is_initially_on = (g.initial_status > 0)
-    reserve = _total_reserves(model, g)
+    reserve = _total_reserves(model, g, sc)
 
     # Gar1962.ProdVars
     prod_above = model[:prod_above]
@@ -37,10 +38,10 @@ function _add_ramp_eqs!(
         if t == 1
             if is_initially_on
                 # min power is _not_ multiplied by is_on because if !is_on, then ramp up is irrelevant
-                eq_ramp_up[gn, t] = @constraint(
+                eq_ramp_up[sc.name, gn, t] = @constraint(
                     model,
                     g.min_power[t] +
-                    prod_above[gn, t] +
+                    prod_above[sc.name, gn, t] +
                     (RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) <=
                     g.initial_power + RU
                 )
@@ -48,16 +49,16 @@ function _add_ramp_eqs!(
         else
             max_prod_this_period =
                 g.min_power[t] * is_on[gn, t] +
-                prod_above[gn, t] +
+                prod_above[sc.name, gn, t] +
                 (
                     RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ?
                     reserve[t] : 0.0
                 )
             min_prod_last_period =
-                g.min_power[t-1] * is_on[gn, t-1] + prod_above[gn, t-1]
+                g.min_power[t-1] * is_on[gn, t-1] + prod_above[sc.name, gn, t-1]
 
             # Equation (24) in Kneuven et al. (2020)
-            eq_ramp_up[gn, t] = @constraint(
+            eq_ramp_up[sc.name, gn, t] = @constraint(
                 model,
                 max_prod_this_period - min_prod_last_period <=
                 RU * is_on[gn, t-1] + SU * switch_on[gn, t]
@@ -71,24 +72,25 @@ function _add_ramp_eqs!(
                 #      min_power + RD < initial_power < SD
                 #      then the generator should be able to shut down at time t = 1,
                 #      but the constraint below will force the unit to produce power
-                eq_ramp_down[gn, t] = @constraint(
+                eq_ramp_down[sc.name, gn, t] = @constraint(
                     model,
-                    g.initial_power - (g.min_power[t] + prod_above[gn, t]) <= RD
+                    g.initial_power -
+                    (g.min_power[t] + prod_above[sc.name, gn, t]) <= RD
                 )
             end
         else
             max_prod_last_period =
                 g.min_power[t-1] * is_on[gn, t-1] +
-                prod_above[gn, t-1] +
+                prod_above[sc.name, gn, t-1] +
                 (
                     RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ?
                     reserve[t-1] : 0.0
                 )
             min_prod_this_period =
-                g.min_power[t] * is_on[gn, t] + prod_above[gn, t]
+                g.min_power[t] * is_on[gn, t] + prod_above[sc.name, gn, t]
 
             # Equation (25) in Kneuven et al. (2020)
-            eq_ramp_down[gn, t] = @constraint(
+            eq_ramp_down[sc.name, gn, t] = @constraint(
                 model,
                 max_prod_last_period - min_prod_this_period <=
                 RD * is_on[gn, t] + SD * switch_off[gn, t]

src/model/formulations/CarArr2006/pwlcosts.jl

@@ -8,6 +8,7 @@ function _add_production_piecewise_linear_eqs!(
     formulation_prod_vars::Gar1962.ProdVars,
     formulation_pwl_costs::CarArr2006.PwlCosts,
     formulation_status_vars::StatusVarsFormulation,
+    sc::UnitCommitmentScenario,
 )::Nothing
     eq_prod_above_def = _init(model, :eq_prod_above_def)
     eq_segprod_limit = _init(model, :eq_segprod_limit)
@@ -26,28 +27,32 @@ function _add_production_piecewise_linear_eqs!(
             #     difference between max power for segments k and k-1 so the
             #     value of cost_segments[k].mw[t] is the max production *for
             #     that segment*
-            eq_segprod_limit[gn, t, k] = @constraint(
+            eq_segprod_limit[sc.name, gn, t, k] = @constraint(
                 model,
-                segprod[gn, t, k] <= g.cost_segments[k].mw[t]
+                segprod[sc.name, gn, t, k] <= g.cost_segments[k].mw[t]
             )
 
             # Also add this as an explicit upper bound on segprod to make the
             # solver's work a bit easier
-            set_upper_bound(segprod[gn, t, k], g.cost_segments[k].mw[t])
+            set_upper_bound(
+                segprod[sc.name, gn, t, k],
+                g.cost_segments[k].mw[t],
+            )
 
             # Definition of production
             # Equation (43) in Kneuven et al. (2020)
-            eq_prod_above_def[gn, t] = @constraint(
+            eq_prod_above_def[sc.name, gn, t] = @constraint(
                 model,
-                prod_above[gn, t] == sum(segprod[gn, t, k] for k in 1:K)
+                prod_above[sc.name, gn, t] ==
+                sum(segprod[sc.name, gn, t, k] for k in 1:K)
             )
 
             # Objective function
             # Equation (44) in Kneuven et al. (2020)
             add_to_expression!(
                 model[:obj],
-                segprod[gn, t, k],
+                segprod[sc.name, gn, t, k],
-                g.cost_segments[k].cost[t],
+                sc.probability * g.cost_segments[k].cost[t],
             )
         end
     end

src/model/formulations/DamKucRajAta2016/ramp.jl

@@ -8,6 +8,7 @@ function _add_ramp_eqs!(
     formulation_prod_vars::Gar1962.ProdVars,
     formulation_ramping::DamKucRajAta2016.Ramping,
     formulation_status_vars::Gar1962.StatusVars,
+    sc::UnitCommitmentScenario,
 )::Nothing
     # TODO: Move upper case constants to model[:instance]
     RESERVES_WHEN_START_UP = true
@@ -23,7 +24,7 @@ function _add_ramp_eqs!(
     gn = g.name
     eq_str_ramp_down = _init(model, :eq_str_ramp_down)
     eq_str_ramp_up = _init(model, :eq_str_ramp_up)
-    reserve = _total_reserves(model, g)
+    reserve = _total_reserves(model, g, sc)
 
     # Gar1962.ProdVars
     prod_above = model[:prod_above]
@@ -48,15 +49,15 @@ function _add_ramp_eqs!(
         # end
 
         max_prod_this_period =
-            prod_above[gn, t] +
+            prod_above[sc.name, gn, t] +
             (RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0)
         min_prod_last_period = 0.0
         if t > 1 && time_invariant
-            min_prod_last_period = prod_above[gn, t-1]
+            min_prod_last_period = prod_above[sc.name, gn, t-1]
 
             # Equation (35) in Kneuven et al. (2020)
             # Sparser version of (24)
-            eq_str_ramp_up[gn, t] = @constraint(
+            eq_str_ramp_up[sc.name, gn, t] = @constraint(
                 model,
                 max_prod_this_period - min_prod_last_period <=
                 (SU - g.min_power[t] - RU) * switch_on[gn, t] +
@@ -65,7 +66,8 @@ function _add_ramp_eqs!(
         elseif (t == 1 && is_initially_on) || (t > 1 && !time_invariant)
             if t > 1
                 min_prod_last_period =
-                    prod_above[gn, t-1] + g.min_power[t-1] * is_on[gn, t-1]
+                    prod_above[sc.name, gn, t-1] +
+                    g.min_power[t-1] * is_on[gn, t-1]
             else
                 min_prod_last_period = max(g.initial_power, 0.0)
             end
@@ -76,7 +78,7 @@ function _add_ramp_eqs!(
 
             # Modified version of equation (35) in Kneuven et al. (2020)
             # Equivalent to (24)
-            eq_str_ramp_up[gn, t] = @constraint(
+            eq_str_ramp_up[sc.name, gn, t] = @constraint(
                 model,
                 max_prod_this_period - min_prod_last_period <=
                 (SU - RU) * switch_on[gn, t] + RU * is_on[gn, t]
@@ -88,7 +90,7 @@ function _add_ramp_eqs!(
                 t > 1 && (RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN) ?
                 reserve[t-1] : 0.0
             )
-        min_prod_this_period = prod_above[gn, t]
+        min_prod_this_period = prod_above[sc.name, gn, t]
         on_last_period = 0.0
         if t > 1
             on_last_period = is_on[gn, t-1]
@@ -98,7 +100,7 @@ function _add_ramp_eqs!(
 
         if t > 1 && time_invariant
             # Equation (36) in Kneuven et al. (2020)
-            eq_str_ramp_down[gn, t] = @constraint(
+            eq_str_ramp_down[sc.name, gn, t] = @constraint(
                 model,
                 max_prod_last_period - min_prod_this_period <=
                 (SD - g.min_power[t] - RD) * switch_off[gn, t] +
@@ -110,7 +112,7 @@ function _add_ramp_eqs!(
 
             # Modified version of equation (36) in Kneuven et al. (2020)
             # Equivalent to (25)
-            eq_str_ramp_down[gn, t] = @constraint(
+            eq_str_ramp_down[sc.name, gn, t] = @constraint(
                 model,
                 max_prod_last_period - min_prod_this_period <=
                 (SD - RD) * switch_off[gn, t] + RD * on_last_period

src/model/formulations/Gar1962/prod.jl

@@ -6,14 +6,15 @@ function _add_production_vars!(
     model::JuMP.Model,
     g::Unit,
     formulation_prod_vars::Gar1962.ProdVars,
+    sc::UnitCommitmentScenario,
 )::Nothing
     prod_above = _init(model, :prod_above)
     segprod = _init(model, :segprod)
     for t in 1:model[:instance].time
         for k in 1:length(g.cost_segments)
-            segprod[g.name, t, k] = @variable(model, lower_bound = 0)
+            segprod[sc.name, g.name, t, k] = @variable(model, lower_bound = 0)
         end
-        prod_above[g.name, t] = @variable(model, lower_bound = 0)
+        prod_above[sc.name, g.name, t] = @variable(model, lower_bound = 0)
     end
     return
 end
@@ -22,16 +23,16 @@ function _add_production_limit_eqs!(
     model::JuMP.Model,
     g::Unit,
     formulation_prod_vars::Gar1962.ProdVars,
+    sc::UnitCommitmentScenario,
 )::Nothing
     eq_prod_limit = _init(model, :eq_prod_limit)
     is_on = model[:is_on]
     prod_above = model[:prod_above]
-    reserve = _total_reserves(model, g)
+    reserve = _total_reserves(model, g, sc)
     gn = g.name
     for t in 1:model[:instance].time
         # Objective function terms for production costs
         # Part of (69) of Kneuven et al. (2020) as C^R_g * u_g(t) term
-        add_to_expression!(model[:obj], is_on[gn, t], g.min_power_cost[t])
 
         # Production limit
         # Equation (18) in Kneuven et al. (2020)
@@ -42,9 +43,10 @@ function _add_production_limit_eqs!(
         if power_diff < 1e-7
             power_diff = 0.0
         end
-        eq_prod_limit[gn, t] = @constraint(
+        eq_prod_limit[sc.name, gn, t] = @constraint(
             model,
-            prod_above[gn, t] + reserve[t] <= power_diff * is_on[gn, t]
+            prod_above[sc.name, gn, t] + reserve[t] <=
+            power_diff * is_on[gn, t]
         )
     end
 end

src/model/formulations/Gar1962/pwlcosts.jl

@@ -8,6 +8,7 @@ function _add_production_piecewise_linear_eqs!(
     formulation_prod_vars::Gar1962.ProdVars,
     formulation_pwl_costs::Gar1962.PwlCosts,
     formulation_status_vars::Gar1962.StatusVars,
+    sc::UnitCommitmentScenario,
 )::Nothing
     eq_prod_above_def = _init(model, :eq_prod_above_def)
     eq_segprod_limit = _init(model, :eq_segprod_limit)
@@ -24,9 +25,10 @@ function _add_production_piecewise_linear_eqs!(
     for t in 1:model[:instance].time
         # Definition of production
         # Equation (43) in Kneuven et al. (2020)
-        eq_prod_above_def[gn, t] = @constraint(
+        eq_prod_above_def[sc.name, gn, t] = @constraint(
             model,
-            prod_above[gn, t] == sum(segprod[gn, t, k] for k in 1:K)
+            prod_above[sc.name, gn, t] ==
+            sum(segprod[sc.name, gn, t, k] for k in 1:K)
         )
 
         for k in 1:K
@@ -37,21 +39,25 @@ function _add_production_piecewise_linear_eqs!(
             #     difference between max power for segments k and k-1 so the
             #     value of cost_segments[k].mw[t] is the max production *for
             #     that segment*
-            eq_segprod_limit[gn, t, k] = @constraint(
+            eq_segprod_limit[sc.name, gn, t, k] = @constraint(
                 model,
-                segprod[gn, t, k] <= g.cost_segments[k].mw[t] * is_on[gn, t]
+                segprod[sc.name, gn, t, k] <=
+                g.cost_segments[k].mw[t] * is_on[gn, t]
             )
 
             # Also add this as an explicit upper bound on segprod to make the
             # solver's work a bit easier
-            set_upper_bound(segprod[gn, t, k], g.cost_segments[k].mw[t])
+            set_upper_bound(
+                segprod[sc.name, gn, t, k],
+                g.cost_segments[k].mw[t],
+            )
 
             # Objective function
             # Equation (44) in Kneuven et al. (2020)
             add_to_expression!(
                 model[:obj],
-                segprod[gn, t, k],
+                segprod[sc.name, gn, t, k],
-                g.cost_segments[k].cost[t],
+                sc.probability * g.cost_segments[k].cost[t],
             )
         end
     end

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

src/model/formulations/KnuOstWat2018/pwlcosts.jl

@@ -8,6 +8,7 @@ function _add_production_piecewise_linear_eqs!(
     formulation_prod_vars::Gar1962.ProdVars,
     formulation_pwl_costs::KnuOstWat2018.PwlCosts,
     formulation_status_vars::Gar1962.StatusVars,
+    sc::UnitCommitmentScenario,
 )::Nothing
     eq_prod_above_def = _init(model, :eq_prod_above_def)
     eq_segprod_limit_a = _init(model, :eq_segprod_limit_a)
@@ -58,27 +59,27 @@ function _add_production_piecewise_linear_eqs!(
 
             if g.min_uptime > 1
                 # Equation (46) in Kneuven et al. (2020)
-                eq_segprod_limit_a[gn, t, k] = @constraint(
+                eq_segprod_limit_a[sc.name, gn, t, k] = @constraint(
                     model,
-                    segprod[gn, t, k] <=
+                    segprod[sc.name, gn, t, k] <=
                     g.cost_segments[k].mw[t] * is_on[gn, t] -
                     Cv * switch_on[gn, t] -
                     (t < T ? Cw * switch_off[gn, t+1] : 0.0)
                 )
             else
                 # Equation (47a)/(48a) in Kneuven et al. (2020)
-                eq_segprod_limit_b[gn, t, k] = @constraint(
+                eq_segprod_limit_b[sc.name, gn, t, k] = @constraint(
                     model,
-                    segprod[gn, t, k] <=
+                    segprod[sc.name, gn, t, k] <=
                     g.cost_segments[k].mw[t] * is_on[gn, t] -
                     Cv * switch_on[gn, t] -
                     (t < T ? max(0, Cv - Cw) * switch_off[gn, t+1] : 0.0)
                 )
 
                 # Equation (47b)/(48b) in Kneuven et al. (2020)
-                eq_segprod_limit_c[gn, t, k] = @constraint(
+                eq_segprod_limit_c[sc.name, gn, t, k] = @constraint(
                     model,
-                    segprod[gn, t, k] <=
+                    segprod[sc.name, gn, t, k] <=
                     g.cost_segments[k].mw[t] * is_on[gn, t] -
                     max(0, Cw - Cv) * switch_on[gn, t] -
                     (t < T ? Cw * switch_off[gn, t+1] : 0.0)
@@ -87,22 +88,26 @@ function _add_production_piecewise_linear_eqs!(
 
             # Definition of production
             # Equation (43) in Kneuven et al. (2020)
-            eq_prod_above_def[gn, t] = @constraint(
+            eq_prod_above_def[sc.name, gn, t] = @constraint(
                 model,
-                prod_above[gn, t] == sum(segprod[gn, t, k] for k in 1:K)
+                prod_above[sc.name, gn, t] ==
+                sum(segprod[sc.name, gn, t, k] for k in 1:K)
             )
 
             # Objective function
             # Equation (44) in Kneuven et al. (2020)
             add_to_expression!(
                 model[:obj],
-                segprod[gn, t, k],
+                segprod[sc.name, gn, t, k],
                 g.cost_segments[k].cost[t],
             )
 
             # Also add an explicit upper bound on segprod to make the solver's
             # work a bit easier
-            set_upper_bound(segprod[gn, t, k], g.cost_segments[k].mw[t])
+            set_upper_bound(
+                segprod[sc.name, gn, t, k],
+                g.cost_segments[k].mw[t],
+            )
         end
     end
 end

src/model/formulations/MorLatRam2013/ramp.jl

@@ -8,6 +8,7 @@ function _add_ramp_eqs!(
     formulation_prod_vars::Gar1962.ProdVars,
     formulation_ramping::MorLatRam2013.Ramping,
     formulation_status_vars::Gar1962.StatusVars,
+    sc::UnitCommitmentScenario,
 )::Nothing
     # TODO: Move upper case constants to model[:instance]
     RESERVES_WHEN_START_UP = true
@@ -22,7 +23,7 @@ function _add_ramp_eqs!(
     gn = g.name
     eq_ramp_down = _init(model, :eq_ramp_down)
     eq_ramp_up = _init(model, :eq_str_ramp_up)
-    reserve = _total_reserves(model, g)
+    reserve = _total_reserves(model, g, sc)
 
     # Gar1962.ProdVars
     prod_above = model[:prod_above]
@@ -39,10 +40,10 @@ function _add_ramp_eqs!(
         # Ramp up limit
         if t == 1
             if is_initially_on
-                eq_ramp_up[gn, t] = @constraint(
+                eq_ramp_up[sc.name, gn, t] = @constraint(
                     model,
                     g.min_power[t] +
-                    prod_above[gn, t] +
+                    prod_above[sc.name, gn, t] +
                     (RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) <=
                     g.initial_power + RU
                 )
@@ -58,13 +59,14 @@ function _add_ramp_eqs!(
                 SU = g.startup_limit
                 max_prod_this_period =
                     g.min_power[t] * is_on[gn, t] +
-                    prod_above[gn, t] +
+                    prod_above[sc.name, gn, t] +
                     (
                         RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ?
                         reserve[t] : 0.0
                     )
                 min_prod_last_period =
-                    g.min_power[t-1] * is_on[gn, t-1] + prod_above[gn, t-1]
+                    g.min_power[t-1] * is_on[gn, t-1] +
+                    prod_above[sc.name, gn, t-1]
                 eq_ramp_up[gn, t] = @constraint(
                     model,
                     max_prod_this_period - min_prod_last_period <=
@@ -74,11 +76,11 @@ function _add_ramp_eqs!(
                 # Equation (26) in Kneuven et al. (2020)
                 # TODO: what if RU < SU? places too stringent upper bound
                 # prod_above[gn, t] when starting up, and creates diff with (24).
-                eq_ramp_up[gn, t] = @constraint(
+                eq_ramp_up[sc.name, gn, t] = @constraint(
                     model,
-                    prod_above[gn, t] +
+                    prod_above[sc.name, gn, t] +
                     (RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) -
-                    prod_above[gn, t-1] <= RU
+                    prod_above[sc.name, gn, t-1] <= RU
                 )
             end
         end
@@ -90,9 +92,10 @@ function _add_ramp_eqs!(
                 #        min_power + RD < initial_power < SD
                 #      then the generator should be able to shut down at time t = 1,
                 #      but the constraint below will force the unit to produce power
-                eq_ramp_down[gn, t] = @constraint(
+                eq_ramp_down[sc.name, gn, t] = @constraint(
                     model,
-                    g.initial_power - (g.min_power[t] + prod_above[gn, t]) <= RD
+                    g.initial_power -
+                    (g.min_power[t] + prod_above[sc.name, gn, t]) <= RD
                 )
             end
         else
@@ -102,13 +105,13 @@ function _add_ramp_eqs!(
                 SD = g.shutdown_limit
                 max_prod_last_period =
                     g.min_power[t-1] * is_on[gn, t-1] +
-                    prod_above[gn, t-1] +
+                    prod_above[sc.name, gn, t-1] +
                     (
                         RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ?
                         reserve[t-1] : 0.0
                     )
                 min_prod_this_period =
-                    g.min_power[t] * is_on[gn, t] + prod_above[gn, t]
+                    g.min_power[t] * is_on[gn, t] + prod_above[sc.name, gn, t]
                 eq_ramp_down[gn, t] = @constraint(
                     model,
                     max_prod_last_period - min_prod_this_period <=
@@ -118,11 +121,11 @@ function _add_ramp_eqs!(
                 # Equation (27) in Kneuven et al. (2020)
                 # TODO: Similar to above, what to do if shutting down in time t
                 # and RD < SD? There is a difference with (25).
-                eq_ramp_down[gn, t] = @constraint(
+                eq_ramp_down[sc.name, gn, t] = @constraint(
                     model,
-                    prod_above[gn, t-1] +
+                    prod_above[sc.name, gn, t-1] +
                     (RESERVES_WHEN_RAMP_DOWN ? reserve[t-1] : 0.0) -
-                    prod_above[gn, t] <= RD
+                    prod_above[sc.name, gn, t] <= RD
                 )
             end
         end

src/model/formulations/PanGua2016/ramp.jl

@@ -8,11 +8,12 @@ function _add_ramp_eqs!(
     formulation_prod_vars::Gar1962.ProdVars,
     formulation_ramping::PanGua2016.Ramping,
     formulation_status_vars::Gar1962.StatusVars,
+    sc::UnitCommitmentScenario,
 )::Nothing
     # TODO: Move upper case constants to model[:instance]
     RESERVES_WHEN_SHUT_DOWN = true
     gn = g.name
-    reserve = _total_reserves(model, g)
+    reserve = _total_reserves(model, g, sc)
     eq_str_prod_limit = _init(model, :eq_str_prod_limit)
     eq_prod_limit_ramp_up_extra_period =
         _init(model, :eq_prod_limit_ramp_up_extra_period)
@@ -52,9 +53,9 @@ function _add_ramp_eqs!(
             # Generalization of (20)
             # Necessary that if any of the switch_on = 1 in the sum,
             # then switch_off[gn, t+1] = 0
-            eq_str_prod_limit[gn, t] = @constraint(
+            eq_str_prod_limit[sc.name, gn, t] = @constraint(
                 model,
-                prod_above[gn, t] +
+                prod_above[sc.name, gn, t] +
                 g.min_power[t] * is_on[gn, t] +
                 reserve[t] <=
                 Pbar * is_on[gn, t] -
@@ -67,16 +68,17 @@ function _add_ramp_eqs!(
             if UT - 2 < TRU
                 # Equation (40) in Kneuven et al. (2020)
                 # Covers an additional time period of the ramp-up trajectory, compared to (38)
-                eq_prod_limit_ramp_up_extra_period[gn, t] = @constraint(
+                eq_prod_limit_ramp_up_extra_period[sc.name, gn, t] =
-                    model,
+                    @constraint(
-                    prod_above[gn, t] +
+                        model,
-                    g.min_power[t] * is_on[gn, t] +
+                        prod_above[sc.name, gn, t] +
-                    reserve[t] <=
+                        g.min_power[t] * is_on[gn, t] +
-                    Pbar * is_on[gn, t] - sum(
+                        reserve[t] <=
-                        (Pbar - (SU + i * RU)) * switch_on[gn, t-i] for
+                        Pbar * is_on[gn, t] - sum(
-                        i in 0:min(UT - 1, TRU, t - 1)
+                            (Pbar - (SU + i * RU)) * switch_on[gn, t-i] for
+                            i in 0:min(UT - 1, TRU, t - 1)
+                        )
                     )
-                )
             end
 
             # Add in shutdown trajectory if KSD >= 0 (else this is dominated by (38))
@@ -84,9 +86,9 @@ function _add_ramp_eqs!(
             if KSD > 0
                 KSU = min(TRU, UT - 2 - KSD, t - 1)
                 # Equation (41) in Kneuven et al. (2020)
-                eq_prod_limit_shutdown_trajectory[gn, t] = @constraint(
+                eq_prod_limit_shutdown_trajectory[sc.name, gn, t] = @constraint(
                     model,
-                    prod_above[gn, t] +
+                    prod_above[sc.name, gn, t] +
                     g.min_power[t] * is_on[gn, t] +
                     (RESERVES_WHEN_SHUT_DOWN ? reserve[t] : 0.0) <=
                     Pbar * is_on[gn, t] - sum(

src/model/formulations/WanHob2016/ramp.jl

@@ -8,6 +8,7 @@ function _add_ramp_eqs!(
     ::Gar1962.ProdVars,
     ::WanHob2016.Ramping,
     ::Gar1962.StatusVars,
+    sc::UnitCommitmentScenario,
 )::Nothing
     is_initially_on = (g.initial_status > 0)
     SU = g.startup_limit
@@ -38,41 +39,43 @@ function _add_ramp_eqs!(
         for t in 1:model[:instance].time
             @constraint(
                 model,
-                prod_above[gn, t] + (is_on[gn, t] * minp[t]) <= mfg[rn, gn, t]
+                prod_above[sc.name, gn, t] + (is_on[gn, t] * minp[t]) <=
+                mfg[sc.name, gn, t]
             ) # Eq. (19) in Wang & Hobbs (2016)
-            @constraint(model, mfg[rn, gn, t] <= is_on[gn, t] * maxp[t]) # Eq. (22) in Wang & Hobbs (2016)
+            @constraint(model, mfg[sc.name, gn, t] <= is_on[gn, t] * maxp[t]) # Eq. (22) in Wang & Hobbs (2016)
             if t != model[:instance].time
                 @constraint(
                     model,
                     minp[t] * (is_on[gn, t+1] + is_on[gn, t] - 1) <=
-                    prod_above[gn, t] - dwflexiramp[rn, gn, t] +
+                    prod_above[sc.name, gn, t] -
-                    (is_on[gn, t] * minp[t])
+                    dwflexiramp[sc.name, rn, gn, t] + (is_on[gn, t] * minp[t])
                 ) # first inequality of Eq. (20) in Wang & Hobbs (2016)
                 @constraint(
                     model,
-                    prod_above[gn, t] - dwflexiramp[rn, gn, t] +
+                    prod_above[sc.name, gn, t] -
+                    dwflexiramp[sc.name, rn, gn, t] +
                     (is_on[gn, t] * minp[t]) <=
-                    mfg[rn, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1]))
+                    mfg[sc.name, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1]))
                 ) # second inequality of Eq. (20) in Wang & Hobbs (2016)
                 @constraint(
                     model,
                     minp[t] * (is_on[gn, t+1] + is_on[gn, t] - 1) <=
-                    prod_above[gn, t] +
+                    prod_above[sc.name, gn, t] +
-                    upflexiramp[rn, 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] +
+                    prod_above[sc.name, gn, t] +
-                    upflexiramp[rn, gn, t] +
+                    upflexiramp[sc.name, rn, gn, t] +
                     (is_on[gn, t] * minp[t]) <=
-                    mfg[rn, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1]))
+                    mfg[sc.name, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1]))
                 ) # second inequality of Eq. (21) in Wang & Hobbs (2016)
                 if t != 1
                     @constraint(
                         model,
-                        mfg[rn, gn, t] <=
+                        mfg[sc.name, gn, t] <=
-                        prod_above[gn, t-1] +
+                        prod_above[sc.name, gn, t-1] +
                         (is_on[gn, t-1] * minp[t]) +
                         (RU * is_on[gn, t-1]) +
                         (SU * (is_on[gn, t] - is_on[gn, t-1])) +
@@ -80,8 +83,13 @@ function _add_ramp_eqs!(
                     ) # Eq. (23) in Wang & Hobbs (2016)
                     @constraint(
                         model,
-                        (prod_above[gn, t-1] + (is_on[gn, t-1] * minp[t])) -
+                        (
-                        (prod_above[gn, t] + (is_on[gn, t] * minp[t])) <=
+                            prod_above[sc.name, gn, t-1] +
+                            (is_on[gn, t-1] * minp[t])
+                        ) - (
+                            prod_above[sc.name, gn, t] +
+                            (is_on[gn, t] * minp[t])
+                        ) <=
                         RD * is_on[gn, t] +
                         SD * (is_on[gn, t-1] - is_on[gn, t]) +
                         maxp[t] * (1 - is_on[gn, t-1])
@@ -89,7 +97,7 @@ function _add_ramp_eqs!(
                 else
                     @constraint(
                         model,
-                        mfg[rn, gn, t] <=
+                        mfg[sc.name, gn, t] <=
                         initial_power +
                         (RU * is_initially_on) +
                         (SU * (is_on[gn, t] - is_initially_on)) +
@@ -97,8 +105,10 @@ function _add_ramp_eqs!(
                     ) # Eq. (23) in Wang & Hobbs (2016) for the first time period
                     @constraint(
                         model,
-                        initial_power -
+                        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 +116,7 @@ function _add_ramp_eqs!(
                 end
                 @constraint(
                     model,
-                    mfg[rn, gn, t] <=
+                    mfg[sc.name, gn, t] <=
                     (SD * (is_on[gn, t] - is_on[gn, t+1])) +
                     (maxp[t] * is_on[gn, t+1])
                 ) # Eq. (24) in Wang & Hobbs (2016)
@@ -114,11 +124,12 @@ function _add_ramp_eqs!(
                     model,
                     -RD * is_on[gn, t+1] -
                     SD * (is_on[gn, t] - is_on[gn, t+1]) -
-                    maxp[t] * (1 - is_on[gn, t]) <= upflexiramp[rn, gn, t]
+                    maxp[t] * (1 - is_on[gn, t]) <=
+                    upflexiramp[sc.name, rn, gn, t]
                 ) # first inequality of Eq. (26) in Wang & Hobbs (2016)
                 @constraint(
                     model,
-                    upflexiramp[rn, gn, t] <=
+                    upflexiramp[sc.name, rn, gn, t] <=
                     RU * is_on[gn, t] +
                     SU * (is_on[gn, t+1] - is_on[gn, t]) +
                     maxp[t] * (1 - is_on[gn, t+1])
@@ -126,11 +137,12 @@ function _add_ramp_eqs!(
                 @constraint(
                     model,
                     -RU * is_on[gn, t] - SU * (is_on[gn, t+1] - is_on[gn, t]) -
-                    maxp[t] * (1 - is_on[gn, t+1]) <= dwflexiramp[rn, gn, t]
+                    maxp[t] * (1 - is_on[gn, t+1]) <=
+                    dwflexiramp[sc.name, rn, gn, t]
                 ) # first inequality of Eq. (27) in Wang & Hobbs (2016)
                 @constraint(
                     model,
-                    dwflexiramp[rn, gn, t] <=
+                    dwflexiramp[sc.name, rn, gn, t] <=
                     RD * is_on[gn, t+1] +
                     SD * (is_on[gn, t] - is_on[gn, t+1]) +
                     maxp[t] * (1 - is_on[gn, t])
@@ -138,26 +150,27 @@ function _add_ramp_eqs!(
                 @constraint(
                     model,
                     -maxp[t] * is_on[gn, t] + minp[t] * is_on[gn, t+1] <=
-                    upflexiramp[rn, gn, t]
+                    upflexiramp[sc.name, rn, gn, t]
                 ) # first inequality of Eq. (28) in Wang & Hobbs (2016)
                 @constraint(
                     model,
-                    upflexiramp[rn, gn, t] <= maxp[t] * is_on[gn, t+1]
+                    upflexiramp[sc.name, rn, gn, t] <= maxp[t] * is_on[gn, t+1]
                 ) # second inequality of Eq. (28) in Wang & Hobbs (2016)
                 @constraint(
                     model,
-                    -maxp[t] * is_on[gn, t+1] <= dwflexiramp[rn, gn, t]
+                    -maxp[t] * is_on[gn, t+1] <=
+                    dwflexiramp[sc.name, rn, gn, t]
                 ) # first inequality of Eq. (29) in Wang & Hobbs (2016)
                 @constraint(
                     model,
-                    dwflexiramp[rn, gn, t] <=
+                    dwflexiramp[sc.name, rn, gn, t] <=
                     (maxp[t] * is_on[gn, t]) - (minp[t] * is_on[gn, t+1])
                 ) # second inequality of Eq. (29) in Wang & Hobbs (2016)
             else
                 @constraint(
                     model,
-                    mfg[rn, gn, t] <=
+                    mfg[sc.name, gn, t] <=
-                    prod_above[gn, t-1] +
+                    prod_above[sc.name, gn, t-1] +
                     (is_on[gn, t-1] * minp[t]) +
                     (RU * is_on[gn, t-1]) +
                     (SU * (is_on[gn, t] - is_on[gn, t-1])) +
@@ -165,8 +178,11 @@ function _add_ramp_eqs!(
                 ) # Eq. (23) in Wang & Hobbs (2016) for the last time period
                 @constraint(
                     model,
-                    (prod_above[gn, t-1] + (is_on[gn, t-1] * minp[t])) -
+                    (
-                    (prod_above[gn, t] + (is_on[gn, t] * minp[t])) <=
+                        prod_above[sc.name, gn, t-1] +
+                        (is_on[gn, t-1] * minp[t])
+                    ) -
+                    (prod_above[sc.name, gn, t] + (is_on[gn, t] * minp[t])) <=
                     RD * is_on[gn, t] +
                     SD * (is_on[gn, t-1] - is_on[gn, t]) +
                     maxp[t] * (1 - is_on[gn, t-1])

src/model/formulations/base/bus.jl

@@ -2,22 +2,30 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
-function _add_bus!(model::JuMP.Model, b::Bus)::Nothing
+function _add_bus!(
+    model::JuMP.Model,
+    b::Bus,
+    sc::UnitCommitmentScenario,
+)::Nothing
     net_injection = _init(model, :expr_net_injection)
     curtail = _init(model, :curtail)
     for t in 1:model[:instance].time
         # Fixed load
-        net_injection[b.name, t] = AffExpr(-b.load[t])
+        net_injection[sc.name, b.name, t] = AffExpr(-b.load[t])
 
         # Load curtailment
-        curtail[b.name, t] =
+        curtail[sc.name, b.name, t] =
             @variable(model, lower_bound = 0, upper_bound = b.load[t])
 
-        add_to_expression!(net_injection[b.name, t], curtail[b.name, t], 1.0)
+        add_to_expression!(
+            net_injection[sc.name, b.name, t],
+            curtail[sc.name, b.name, t],
+            1.0,
+        )
         add_to_expression!(
             model[:obj],
-            curtail[b.name, t],
+            curtail[sc.name, b.name, t],
-            model[:instance].power_balance_penalty[t],
+            sc.power_balance_penalty[t] * sc.probability,
         )
     end
     return

src/model/formulations/base/line.jl

@@ -6,43 +6,43 @@ function _add_transmission_line!(
     model::JuMP.Model,
     lm::TransmissionLine,
     f::ShiftFactorsFormulation,
+    sc::UnitCommitmentScenario,
 )::Nothing
     overflow = _init(model, :overflow)
     for t in 1:model[:instance].time
-        overflow[lm.name, t] = @variable(model, lower_bound = 0)
+        overflow[sc.name, lm.name, t] = @variable(model, lower_bound = 0)
         add_to_expression!(
             model[:obj],
-            overflow[lm.name, t],
+            overflow[sc.name, lm.name, t],
-            lm.flow_limit_penalty[t],
+            lm.flow_limit_penalty[t] * sc.probability,
         )
     end
     return
 end
 
 function _setup_transmission(
-    model::JuMP.Model,
     formulation::ShiftFactorsFormulation,
+    sc::UnitCommitmentScenario,
 )::Nothing
-    instance = model[:instance]
     isf = formulation.precomputed_isf
     lodf = formulation.precomputed_lodf
-    if length(instance.buses) == 1
+    if length(sc.buses) == 1
         isf = zeros(0, 0)
         lodf = zeros(0, 0)
     elseif isf === nothing
         @info "Computing injection shift factors..."
         time_isf = @elapsed begin
             isf = UnitCommitment._injection_shift_factors(
-                lines = instance.lines,
+                buses = sc.buses,
-                buses = instance.buses,
+                lines = sc.lines,
             )
         end
         @info @sprintf("Computed ISF in %.2f seconds", time_isf)
         @info "Computing line outage factors..."
         time_lodf = @elapsed begin
             lodf = UnitCommitment._line_outage_factors(
-                lines = instance.lines,
+                buses = sc.buses,
-                buses = instance.buses,
+                lines = sc.lines,
                 isf = isf,
             )
         end
@@ -55,7 +55,7 @@ function _setup_transmission(
         isf[abs.(isf).<formulation.isf_cutoff] .= 0
         lodf[abs.(lodf).<formulation.lodf_cutoff] .= 0
     end
-    model[:isf] = isf
+    sc.isf = isf
-    model[:lodf] = lodf
+    sc.lodf = lodf
     return
 end

src/model/formulations/base/psload.jl

@@ -5,21 +5,26 @@
 function _add_price_sensitive_load!(
     model::JuMP.Model,
     ps::PriceSensitiveLoad,
+    sc::UnitCommitmentScenario,
 )::Nothing
     loads = _init(model, :loads)
     net_injection = _init(model, :expr_net_injection)
     for t in 1:model[:instance].time
         # Decision variable
-        loads[ps.name, t] =
+        loads[sc.name, ps.name, t] =
             @variable(model, lower_bound = 0, upper_bound = ps.demand[t])
 
         # Objective function terms
-        add_to_expression!(model[:obj], loads[ps.name, t], -ps.revenue[t])
+        add_to_expression!(
+            model[:obj],
+            loads[sc.name, ps.name, t],
+            -ps.revenue[t] * sc.probability,
+        )
 
         # Net injection
         add_to_expression!(
-            net_injection[ps.bus.name, t],
+            net_injection[sc.name, ps.bus.name, t],
-            loads[ps.name, t],
+            loads[sc.name, ps.name, t],
             -1.0,
         )
     end

src/model/formulations/base/sensitivity.jl

@@ -18,7 +18,7 @@ function _injection_shift_factors(;
     lines::Array{TransmissionLine},
 )
     susceptance = _susceptance_matrix(lines)
-    incidence = _reduced_incidence_matrix(lines = lines, buses = buses)
+    incidence = _reduced_incidence_matrix(buses = buses, lines = lines)
     laplacian = transpose(incidence) * susceptance * incidence
     isf = susceptance * incidence * inv(Array(laplacian))
     return isf

src/model/formulations/base/system.jl

@@ -2,54 +2,68 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
-function _add_system_wide_eqs!(model::JuMP.Model)::Nothing
+function _add_system_wide_eqs!(
-    _add_net_injection_eqs!(model)
+    model::JuMP.Model,
-    _add_spinning_reserve_eqs!(model)
+    sc::UnitCommitmentScenario,
-    _add_flexiramp_reserve_eqs!(model)
+)::Nothing
+    _add_net_injection_eqs!(model, sc)
+    _add_spinning_reserve_eqs!(model, sc)
+    _add_flexiramp_reserve_eqs!(model, sc)
     return
 end
 
-function _add_net_injection_eqs!(model::JuMP.Model)::Nothing
+function _add_net_injection_eqs!(
+    model::JuMP.Model,
+    sc::UnitCommitmentScenario,
+)::Nothing
     T = model[:instance].time
     net_injection = _init(model, :net_injection)
     eq_net_injection = _init(model, :eq_net_injection)
     eq_power_balance = _init(model, :eq_power_balance)
-    for t in 1:T, b in model[:instance].buses
+    for t in 1:T, b in sc.buses
-        n = net_injection[b.name, t] = @variable(model)
+        n = net_injection[sc.name, b.name, t] = @variable(model)
-        eq_net_injection[b.name, t] =
+        eq_net_injection[sc.name, b.name, t] = @constraint(
-            @constraint(model, -n + model[:expr_net_injection][b.name, t] == 0)
+            model,
+            -n + model[:expr_net_injection][sc.name, b.name, t] == 0
+        )
     end
     for t in 1:T
-        eq_power_balance[t] = @constraint(
+        eq_power_balance[sc.name, t] = @constraint(
             model,
-            sum(net_injection[b.name, t] for b in model[:instance].buses) == 0
+            sum(net_injection[sc.name, b.name, t] for b in sc.buses) == 0
         )
     end
     return
 end
 
-function _add_spinning_reserve_eqs!(model::JuMP.Model)::Nothing
+function _add_spinning_reserve_eqs!(
-    instance = model[:instance]
+    model::JuMP.Model,
+    sc::UnitCommitmentScenario,
+)::Nothing
+    T = model[:instance].time
     eq_min_spinning_reserve = _init(model, :eq_min_spinning_reserve)
-    for r in instance.reserves
+    for r in sc.reserves
         r.type == "spinning" || continue
-        for t in 1:instance.time
+        for t in 1:T
             # Equation (68) in Kneuven et al. (2020)
             # As in Morales-España et al. (2013a)
             # Akin to the alternative formulation with max_power_avail
             # from Carrión and Arroyo (2006) and Ostrowski et al. (2012)
-            eq_min_spinning_reserve[r.name, t] = @constraint(
+            eq_min_spinning_reserve[sc.name, r.name, t] = @constraint(
                 model,
-                sum(model[:reserve][r.name, g.name, t] for g in r.units) +
+                sum(
-                model[:reserve_shortfall][r.name, t] >= r.amount[t]
+                    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,
+                    r.shortfall_penalty * sc.probability,
-                    model[:reserve_shortfall][r.name, t],
+                    model[:reserve_shortfall][sc.name, r.name, t],
                 )
             end
         end
@@ -57,7 +71,10 @@ function _add_spinning_reserve_eqs!(model::JuMP.Model)::Nothing
     return
 end
 
-function _add_flexiramp_reserve_eqs!(model::JuMP.Model)::Nothing
+function _add_flexiramp_reserve_eqs!(
+    model::JuMP.Model,
+    sc::UnitCommitmentScenario,
+)::Nothing
     # Note: The flexpramp requirements in Wang & Hobbs (2016) are imposed as hard constraints 
     #       through Eq. (17) and Eq. (18). The constraints eq_min_upflexiramp and eq_min_dwflexiramp 
     #       provided below are modified versions of Eq. (17) and Eq. (18), respectively, in that   
@@ -65,29 +82,37 @@ function _add_flexiramp_reserve_eqs!(model::JuMP.Model)::Nothing
     #       objective function.
     eq_min_upflexiramp = _init(model, :eq_min_upflexiramp)
     eq_min_dwflexiramp = _init(model, :eq_min_dwflexiramp)
-    instance = model[:instance]
+    T = model[:instance].time
-    for r in instance.reserves
+    for r in sc.reserves
         r.type == "flexiramp" || continue
-        for t in 1:instance.time
+        for t in 1:T
             # Eq. (17) in Wang & Hobbs (2016)
-            eq_min_upflexiramp[r.name, t] = @constraint(
+            eq_min_upflexiramp[sc.name, r.name, t] = @constraint(
                 model,
-                sum(model[:upflexiramp][r.name, g.name, t] for g in r.units) + model[:upflexiramp_shortfall][r.name, t] >= r.amount[t]
+                sum(
+                    model[:upflexiramp][sc.name, r.name, g.name, t] for
+                    g in r.units
+                ) + model[:upflexiramp_shortfall][sc.name, r.name, t] >=
+                r.amount[t]
             )
             # Eq. (18) in Wang & Hobbs (2016)
-            eq_min_dwflexiramp[r.name, t] = @constraint(
+            eq_min_dwflexiramp[sc.name, r.name, t] = @constraint(
                 model,
-                sum(model[:dwflexiramp][r.name, g.name, t] for g in r.units) + model[:dwflexiramp_shortfall][r.name, t] >= r.amount[t]
+                sum(
+                    model[:dwflexiramp][sc.name, r.name, g.name, t] for
+                    g in r.units
+                ) + model[:dwflexiramp_shortfall][sc.name, r.name, t] >=
+                r.amount[t]
             )
 
             # Account for flexiramp shortfall contribution to objective
             if r.shortfall_penalty >= 0
                 add_to_expression!(
                     model[:obj],
-                    r.shortfall_penalty,
+                    r.shortfall_penalty * sc.probability,
                     (
-                        model[:upflexiramp_shortfall][r.name, t] +
+                        model[:upflexiramp_shortfall][sc.name, r.name, t] +
-                        model[:dwflexiramp_shortfall][r.name, t]
+                        model[:dwflexiramp_shortfall][sc.name, r.name, t]
                     ),
                 )
             end

src/model/formulations/base/unit.jl

@@ -2,7 +2,13 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
-function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation)
+# Function for adding variables, constraints, and objective function terms
+# related to the binary commitment, startup and shutdown decisions of units
+function _add_unit_commitment!(
+    model::JuMP.Model,
+    g::Unit,
+    formulation::Formulation,
+)
     if !all(g.must_run) && any(g.must_run)
         error("Partially must-run units are not currently supported")
     end
@@ -11,22 +17,40 @@ function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation)
     end
 
     # Variables
-    _add_production_vars!(model, g, formulation.prod_vars)
-    _add_spinning_reserve_vars!(model, g)
-    _add_flexiramp_reserve_vars!(model, g)
     _add_startup_shutdown_vars!(model, g)
     _add_status_vars!(model, g, formulation.status_vars)
 
     # Constraints and objective function
     _add_min_uptime_downtime_eqs!(model, g)
-    _add_net_injection_eqs!(model, g)
+    _add_startup_cost_eqs!(model, g, formulation.startup_costs)
-    _add_production_limit_eqs!(model, g, formulation.prod_vars)
+    _add_status_eqs!(model, g, formulation.status_vars)
+    return
+end
+
+# Function for adding variables, constraints, and objective function terms
+# related to the continuous dispatch decisions of units
+function _add_unit_dispatch!(
+    model::JuMP.Model,
+    g::Unit,
+    formulation::Formulation,
+    sc::UnitCommitmentScenario,
+)
+
+    # Variables
+    _add_production_vars!(model, g, formulation.prod_vars, sc)
+    _add_spinning_reserve_vars!(model, g, sc)
+    _add_flexiramp_reserve_vars!(model, g, sc)
+
+    # Constraints and objective function
+    _add_net_injection_eqs!(model, g, sc)
+    _add_production_limit_eqs!(model, g, formulation.prod_vars, sc)
     _add_production_piecewise_linear_eqs!(
         model,
         g,
         formulation.prod_vars,
         formulation.pwl_costs,
         formulation.status_vars,
+        sc,
     )
     _add_ramp_eqs!(
         model,
@@ -34,26 +58,31 @@ function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation)
         formulation.prod_vars,
         formulation.ramping,
         formulation.status_vars,
+        sc,
     )
-    _add_startup_cost_eqs!(model, g, formulation.startup_costs)
+    _add_startup_shutdown_limit_eqs!(model, g, sc)
-    _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)
+            reserve[sc.name, r.name, g.name, t] =
-            if (r.name, t) ∉ keys(reserve_shortfall)
+                @variable(model, lower_bound = 0)
-                reserve_shortfall[r.name, t] = @variable(model, lower_bound = 0)
+            if (sc.name, r.name, t) ∉ keys(reserve_shortfall)
+                reserve_shortfall[sc.name, r.name, t] =
+                    @variable(model, lower_bound = 0)
                 if r.shortfall_penalty < 0
-                    set_upper_bound(reserve_shortfall[r.name, t], 0.0)
+                    set_upper_bound(reserve_shortfall[sc.name, r.name, t], 0.0)
                 end
             end
         end
@@ -61,27 +90,37 @@ function _add_spinning_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
     return
 end
 
-function _add_flexiramp_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
+function _add_flexiramp_reserve_vars!(
+    model::JuMP.Model,
+    g::Unit,
+    sc::UnitCommitmentScenario,
+)::Nothing
     upflexiramp = _init(model, :upflexiramp)
     upflexiramp_shortfall = _init(model, :upflexiramp_shortfall)
     mfg = _init(model, :mfg)
     dwflexiramp = _init(model, :dwflexiramp)
     dwflexiramp_shortfall = _init(model, :dwflexiramp_shortfall)
-    for r in g.reserves
+    for t in 1:model[:instance].time
-        r.type == "flexiramp" || continue
+        # maximum feasible generation, \bar{g_{its}} in Wang & Hobbs (2016)
-        for t in 1:model[:instance].time
+        mfg[sc.name, g.name, t] = @variable(model, lower_bound = 0)
-            # maximum feasible generation, \bar{g_{its}} in Wang & Hobbs (2016)
+        for r in g.reserves
-            mfg[r.name, g.name, t] = @variable(model, lower_bound = 0)
+            r.type == "flexiramp" || continue
-            upflexiramp[r.name, g.name, t] = @variable(model) # up-flexiramp, ur_{it} in Wang & Hobbs (2016)
+            upflexiramp[sc.name, 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)
+            dwflexiramp[sc.name, r.name, g.name, t] = @variable(model) # down-flexiramp, dr_{it} in Wang & Hobbs (2016)
-            if (r.name, t) ∉ keys(upflexiramp_shortfall)
+            if (sc.name, r.name, t) ∉ keys(upflexiramp_shortfall)
-                upflexiramp_shortfall[r.name, t] =
+                upflexiramp_shortfall[sc.name, r.name, t] =
                     @variable(model, lower_bound = 0)
-                dwflexiramp_shortfall[r.name, t] =
+                dwflexiramp_shortfall[sc.name, r.name, t] =
                     @variable(model, lower_bound = 0)
                 if r.shortfall_penalty < 0
-                    set_upper_bound(upflexiramp_shortfall[r.name, t], 0.0)
+                    set_upper_bound(
-                    set_upper_bound(dwflexiramp_shortfall[r.name, t], 0.0)
+                        upflexiramp_shortfall[sc.name, r.name, t],
+                        0.0,
+                    )
+                    set_upper_bound(
+                        dwflexiramp_shortfall[sc.name, r.name, t],
+                        0.0,
+                    )
                 end
             end
         end
@@ -99,32 +138,36 @@ function _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing
     return
 end
 
-function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
+function _add_startup_shutdown_limit_eqs!(
+    model::JuMP.Model,
+    g::Unit,
+    sc::UnitCommitmentScenario,
+)::Nothing
     eq_shutdown_limit = _init(model, :eq_shutdown_limit)
     eq_startup_limit = _init(model, :eq_startup_limit)
     is_on = model[:is_on]
     prod_above = model[:prod_above]
-    reserve = _total_reserves(model, g)
+    reserve = _total_reserves(model, g, sc)
     switch_off = model[:switch_off]
     switch_on = model[:switch_on]
     T = model[:instance].time
     for t in 1:T
         # Startup limit
-        eq_startup_limit[g.name, t] = @constraint(
+        eq_startup_limit[sc.name, g.name, t] = @constraint(
             model,
-            prod_above[g.name, t] + reserve[t] <=
+            prod_above[sc.name, g.name, t] + reserve[t] <=
             (g.max_power[t] - g.min_power[t]) * is_on[g.name, t] -
             max(0, g.max_power[t] - g.startup_limit) * switch_on[g.name, t]
         )
         # Shutdown limit
         if g.initial_power > g.shutdown_limit
-            eq_shutdown_limit[g.name, 0] =
+            eq_shutdown_limit[sc.name, g.name, 0] =
                 @constraint(model, switch_off[g.name, 1] <= 0)
         end
         if t < T
-            eq_shutdown_limit[g.name, t] = @constraint(
+            eq_shutdown_limit[sc.name, g.name, t] = @constraint(
                 model,
-                prod_above[g.name, t] <=
+                prod_above[sc.name, g.name, t] <=
                 (g.max_power[t] - g.min_power[t]) * is_on[g.name, t] -
                 max(0, g.max_power[t] - g.shutdown_limit) *
                 switch_off[g.name, t+1]
@@ -138,43 +181,44 @@ function _add_ramp_eqs!(
     model::JuMP.Model,
     g::Unit,
     formulation::RampingFormulation,
+    sc::UnitCommitmentScenario,
 )::Nothing
     prod_above = model[:prod_above]
-    reserve = _total_reserves(model, g)
+    reserve = _total_reserves(model, g, sc)
     eq_ramp_up = _init(model, :eq_ramp_up)
     eq_ramp_down = _init(model, :eq_ramp_down)
     for t in 1:model[:instance].time
         # Ramp up limit 
         if t == 1
             if _is_initially_on(g) == 1
-                eq_ramp_up[g.name, t] = @constraint(
+                eq_ramp_up[sc.name, g.name, t] = @constraint(
                     model,
-                    prod_above[g.name, t] + reserve[t] <=
+                    prod_above[sc.name, g.name, t] + reserve[t] <=
                     (g.initial_power - g.min_power[t]) + g.ramp_up_limit
                 )
             end
         else
-            eq_ramp_up[g.name, t] = @constraint(
+            eq_ramp_up[sc.name, g.name, t] = @constraint(
                 model,
-                prod_above[g.name, t] + reserve[t] <=
+                prod_above[sc.name, g.name, t] + reserve[t] <=
-                prod_above[g.name, t-1] + g.ramp_up_limit
+                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[sc.name, g.name, t] >=
-                prod_above[g.name, t-1] - g.ramp_down_limit
+                prod_above[sc.name, g.name, t-1] - g.ramp_down_limit
             )
         end
     end
@@ -223,30 +267,37 @@ function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing
     end
 end
 
-function _add_net_injection_eqs!(model::JuMP.Model, g::Unit)::Nothing
+function _add_net_injection_eqs!(
+    model::JuMP.Model,
+    g::Unit,
+    sc::UnitCommitmentScenario,
+)::Nothing
     expr_net_injection = model[:expr_net_injection]
     for t in 1:model[:instance].time
         # Add to net injection expression
         add_to_expression!(
-            expr_net_injection[g.bus.name, t],
+            expr_net_injection[sc.name, g.bus.name, t],
-            model[:prod_above][g.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

src/solution/fix.jl

@@ -10,37 +10,43 @@ solution. Useful for computing LMPs.
 """
 function fix!(model::JuMP.Model, solution::AbstractDict)::Nothing
     instance, T = model[:instance], model[:instance].time
+    "Production (MW)" ∈ keys(solution) ? solution = Dict("s1" => solution) :
+    nothing
     is_on = model[:is_on]
     prod_above = model[:prod_above]
     reserve = model[:reserve]
-    for g in instance.units
+    for sc in instance.scenarios
-        for t in 1:T
+        for g in sc.units
-            is_on_value = round(solution["Is on"][g.name][t])
-            prod_value =
-                round(solution["Production (MW)"][g.name][t], digits = 5)
-            JuMP.fix(is_on[g.name, t], is_on_value, force = true)
-            JuMP.fix(
-                prod_above[g.name, t],
-                prod_value - is_on_value * g.min_power[t],
-                force = true,
-            )
-        end
-    end
-    for r in instance.reserves
-        r.type == "spinning" || continue
-        for g in r.units
             for t in 1:T
-                reserve_value = round(
+                is_on_value = round(solution[sc.name]["Is on"][g.name][t])
-                    solution["Spinning reserve (MW)"][r.name][g.name][t],
+                prod_value = round(
+                    solution[sc.name]["Production (MW)"][g.name][t],
                     digits = 5,
                 )
+                JuMP.fix(is_on[g.name, t], is_on_value, force = true)
                 JuMP.fix(
-                    reserve[r.name, g.name, t],
+                    prod_above[sc.name, g.name, t],
-                    reserve_value,
+                    prod_value - is_on_value * g.min_power[t],
                     force = true,
                 )
             end
         end
+        for r in sc.reserves
+            r.type == "spinning" || continue
+            for g in r.units
+                for t in 1:T
+                    reserve_value = round(
+                        solution[sc.name]["Spinning reserve (MW)"][r.name][g.name][t],
+                        digits = 5,
+                    )
+                    JuMP.fix(
+                        reserve[sc.name, r.name, g.name, t],
+                        reserve_value,
+                        force = true,
+                    )
+                end
+            end
+        end
     end
     return
 end

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],
+            isf = sc.isf,
-            lodf = model[:lodf],
+            lodf = sc.lodf,
         )
     end
     return
@@ -19,6 +21,7 @@ end
 
 function _enforce_transmission(;
     model::JuMP.Model,
+    sc::UnitCommitmentScenario,
     violation::_Violation,
     isf::Matrix{Float64},
     lodf::Matrix{Float64},
@@ -31,19 +34,21 @@ function _enforce_transmission(;
     if violation.outage_line === nothing
         limit = violation.monitored_line.normal_flow_limit[violation.time]
         @info @sprintf(
-            "    %8.3f MW overflow in %-5s time %3d (pre-contingency)",
+            "    %8.3f MW overflow in %-5s time %3d (pre-contingency, scenario %s)",
             violation.amount,
             violation.monitored_line.name,
             violation.time,
+            sc.name,
         )
     else
         limit = violation.monitored_line.emergency_flow_limit[violation.time]
         @info @sprintf(
-            "    %8.3f MW overflow in %-5s time %3d (outage: line %s)",
+            "    %8.3f MW overflow in %-5s time %3d (outage: line %s, scenario %s)",
             violation.amount,
             violation.monitored_line.name,
             violation.time,
             violation.outage_line.name,
+            sc.name,
         )
     end
 
@@ -51,7 +56,7 @@ function _enforce_transmission(;
     t = violation.time
     flow = @variable(model, base_name = "flow[$fm,$t]")
 
-    v = overflow[violation.monitored_line.name, violation.time]
+    v = overflow[sc.name, violation.monitored_line.name, violation.time]
     @constraint(model, flow <= limit + v)
     @constraint(model, -flow <= limit + v)
 
@@ -59,23 +64,23 @@ function _enforce_transmission(;
         @constraint(
             model,
             flow == sum(
-                net_injection[b.name, violation.time] *
+                net_injection[sc.name, b.name, violation.time] *
                 isf[violation.monitored_line.offset, b.offset] for
-                b in instance.buses if b.offset > 0
+                b in sc.buses if b.offset > 0
             )
         )
     else
         @constraint(
             model,
             flow == sum(
-                net_injection[b.name, violation.time] * (
+                net_injection[sc.name, b.name, violation.time] * (
                     isf[violation.monitored_line.offset, b.offset] + (
                         lodf[
                             violation.monitored_line.offset,
                             violation.outage_line.offset,
                         ] * isf[violation.outage_line.offset, b.offset]
                     )
-                ) for b in instance.buses if b.offset > 0
+                ) for b in sc.buses if b.offset > 0
             )
         )
     end

src/solution/methods/XavQiuWanThi2019/find.jl

@@ -5,39 +5,35 @@
 import Base.Threads: @threads
 
 function _find_violations(
-    model::JuMP.Model;
+    model::JuMP.Model,
+    sc::UnitCommitmentScenario;
     max_per_line::Int,
     max_per_period::Int,
 )
     instance = model[:instance]
     net_injection = model[:net_injection]
     overflow = model[:overflow]
-    length(instance.buses) > 1 || return []
+    length(sc.buses) > 1 || return []
     violations = []
-    @info "Verifying transmission limits..."
+
-    time_screening = @elapsed begin
+    non_slack_buses = [b for b in sc.buses if b.offset > 0]
-        non_slack_buses = [b for b in instance.buses if b.offset > 0]
+    net_injection_values = [
-        net_injection_values = [
+        value(net_injection[sc.name, b.name, t]) for b in non_slack_buses,
-            value(net_injection[b.name, t]) for b in non_slack_buses,
+        t in 1:instance.time
-            t in 1:instance.time
+    ]
-        ]
+    overflow_values = [
-        overflow_values = [
+        value(overflow[sc.name, lm.name, t]) for lm in sc.lines,
-            value(overflow[lm.name, t]) for lm in instance.lines,
+        t in 1:instance.time
-            t in 1:instance.time
+    ]
-        ]
+    violations = UnitCommitment._find_violations(
-        violations = UnitCommitment._find_violations(
+        instance = instance,
-            instance = instance,
+        sc = sc,
-            net_injections = net_injection_values,
+        net_injections = net_injection_values,
-            overflow = overflow_values,
+        overflow = overflow_values,
-            isf = model[:isf],
+        isf = sc.isf,
-            lodf = model[:lodf],
+        lodf = sc.lodf,
-            max_per_line = max_per_line,
+        max_per_line = max_per_line,
-            max_per_period = max_per_period,
+        max_per_period = max_per_period,
-        )
-    end
-    @info @sprintf(
-        "Verified transmission limits in %.2f seconds",
-        time_screening
     )
     return violations
 end
@@ -64,6 +60,7 @@ matrix, where L is the number of transmission lines.
 """
 function _find_violations(;
     instance::UnitCommitmentInstance,
+    sc::UnitCommitmentScenario,
     net_injections::Array{Float64,2},
     overflow::Array{Float64,2},
     isf::Array{Float64,2},
@@ -71,8 +68,8 @@ function _find_violations(;
     max_per_line::Int,
     max_per_period::Int,
 )::Array{_Violation,1}
-    B = length(instance.buses) - 1
+    B = length(sc.buses) - 1
-    L = length(instance.lines)
+    L = length(sc.lines)
     T = instance.time
     K = nthreads()
 
@@ -93,17 +90,17 @@ function _find_violations(;
     post_v::Array{Float64} = zeros(L, L, K)          # post_v[lm, lc, thread]
 
     normal_limits::Array{Float64,2} = [
-        l.normal_flow_limit[t] + overflow[l.offset, t] for
+        l.normal_flow_limit[t] + overflow[l.offset, t] for l in sc.lines,
-        l in instance.lines, t in 1:T
+        t in 1:T
     ]
 
     emergency_limits::Array{Float64,2} = [
-        l.emergency_flow_limit[t] + overflow[l.offset, t] for
+        l.emergency_flow_limit[t] + overflow[l.offset, t] for l in sc.lines,
-        l in instance.lines, t in 1:T
+        t in 1:T
     ]
 
     is_vulnerable::Array{Bool} = zeros(Bool, L)
-    for c in instance.contingencies
+    for c in sc.contingencies
         is_vulnerable[c.lines[1].offset] = true
     end
 
@@ -144,7 +141,7 @@ function _find_violations(;
                     filters[t],
                     _Violation(
                         time = t,
-                        monitored_line = instance.lines[lm],
+                        monitored_line = sc.lines[lm],
                         outage_line = nothing,
                         amount = pre_v[lm, k],
                     ),
@@ -159,8 +156,8 @@ function _find_violations(;
                     filters[t],
                     _Violation(
                         time = t,
-                        monitored_line = instance.lines[lm],
+                        monitored_line = sc.lines[lm],
-                        outage_line = instance.lines[lc],
+                        outage_line = sc.lines[lc],
                         amount = post_v[lm, lc, k],
                     ),
                 )

src/solution/methods/XavQiuWanThi2019/optimize.jl

@@ -12,10 +12,15 @@ function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing
     end
     initial_time = time()
     large_gap = false
-    has_transmission = (length(model[:isf]) > 0)
+    has_transmission = false
-    if has_transmission && method.two_phase_gap
+    for sc in model[:instance].scenarios
-        set_gap(1e-2)
+        if length(sc.isf) > 0
-        large_gap = true
+            has_transmission = true
+        end
+        if has_transmission && method.two_phase_gap
+            set_gap(1e-2)
+            large_gap = true
+        end
     end
     while true
         time_elapsed = time() - initial_time
@@ -31,13 +36,41 @@ function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing
         JuMP.set_time_limit_sec(model, time_remaining)
         @info "Solving MILP..."
         JuMP.optimize!(model)
+
         has_transmission || break
-        violations = _find_violations(
+
-            model,
+        @info "Verifying transmission limits..."
-            max_per_line = method.max_violations_per_line,
+        time_screening = @elapsed begin
-            max_per_period = method.max_violations_per_period,
+            violations = []
+            for sc in model[:instance].scenarios
+                push!(
+                    violations,
+                    _find_violations(
+                        model,
+                        sc,
+                        max_per_line = method.max_violations_per_line,
+                        max_per_period = method.max_violations_per_period,
+                    ),
+                )
+            end
+        end
+        @info @sprintf(
+            "Verified transmission limits in %.2f seconds",
+            time_screening
         )
-        if isempty(violations)
+
+        violations_found = false
+        for v in violations
+            if !isempty(v)
+                violations_found = true
+            end
+        end
+
+        if violations_found
+            for (i, v) in enumerate(violations)
+                _enforce_transmission(model, v, model[:instance].scenarios[i])
+            end
+        else
             @info "No violations found"
             if large_gap
                 large_gap = false
@@ -45,8 +78,6 @@ function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing
             else
                 break
             end
-        else
-            _enforce_transmission(model, violations)
         end
     end
     return

src/solution/solution.jl

@@ -16,34 +16,44 @@ solution = UnitCommitment.solution(model)
 """
 function solution(model::JuMP.Model)::OrderedDict
     instance, T = model[:instance], model[:instance].time
-    function timeseries(vars, collection)
+    function timeseries(vars, collection; sc = nothing)
-        return OrderedDict(
+        if sc === nothing
-            b.name => [round(value(vars[b.name, t]), digits = 5) for t in 1:T]
+            return OrderedDict(
-            for b in collection
+                b.name =>
-        )
+                    [round(value(vars[b.name, t]), digits = 5) for t in 1:T] for
+                b in collection
+            )
+        else
+            return OrderedDict(
+                b.name => [
+                    round(value(vars[sc.name, b.name, t]), digits = 5) for
+                    t in 1:T
+                ] for b in collection
+            )
+        end
     end
-    function production_cost(g)
+    function production_cost(g, sc)
         return [
             value(model[:is_on][g.name, t]) * g.min_power_cost[t] + sum(
                 Float64[
-                    value(model[:segprod][g.name, t, k]) *
+                    value(model[:segprod][sc.name, g.name, t, k]) *
                     g.cost_segments[k].cost[t] for
                     k in 1:length(g.cost_segments)
                 ],
             ) for t in 1:T
         ]
     end
-    function production(g)
+    function production(g, sc)
         return [
             value(model[:is_on][g.name, t]) * g.min_power[t] + sum(
                 Float64[
-                    value(model[:segprod][g.name, t, k]) for
+                    value(model[:segprod][sc.name, g.name, t, k]) for
                     k in 1:length(g.cost_segments)
                 ],
             ) for t in 1:T
         ]
     end
-    function startup_cost(g)
+    function startup_cost(g, sc)
         S = length(g.startup_categories)
         return [
             sum(
@@ -53,66 +63,70 @@ function solution(model::JuMP.Model)::OrderedDict
         ]
     end
     sol = OrderedDict()
-    sol["Production (MW)"] =
+    for sc in instance.scenarios
-        OrderedDict(g.name => production(g) for g in instance.units)
+        sol[sc.name] = OrderedDict()
-    sol["Production cost (\$)"] =
+        sol[sc.name]["Production (MW)"] =
-        OrderedDict(g.name => production_cost(g) for g in instance.units)
+            OrderedDict(g.name => production(g, sc) for g in sc.units)
-    sol["Startup cost (\$)"] =
+        sol[sc.name]["Production cost (\$)"] =
-        OrderedDict(g.name => startup_cost(g) for g in instance.units)
+            OrderedDict(g.name => production_cost(g, sc) for g in sc.units)
-    sol["Is on"] = timeseries(model[:is_on], instance.units)
+        sol[sc.name]["Startup cost (\$)"] =
-    sol["Switch on"] = timeseries(model[:switch_on], instance.units)
+            OrderedDict(g.name => startup_cost(g, sc) for g in sc.units)
-    sol["Switch off"] = timeseries(model[:switch_off], instance.units)
+        sol[sc.name]["Is on"] = timeseries(model[:is_on], sc.units)
-    sol["Net injection (MW)"] =
+        sol[sc.name]["Switch on"] = timeseries(model[:switch_on], sc.units)
-        timeseries(model[:net_injection], instance.buses)
+        sol[sc.name]["Switch off"] = timeseries(model[:switch_off], sc.units)
-    sol["Load curtail (MW)"] = timeseries(model[:curtail], instance.buses)
+        sol[sc.name]["Net injection (MW)"] =
-    if !isempty(instance.lines)
+            timeseries(model[:net_injection], sc.buses, sc = sc)
-        sol["Line overflow (MW)"] = timeseries(model[:overflow], instance.lines)
+        sol[sc.name]["Load curtail (MW)"] =
+            timeseries(model[:curtail], sc.buses, sc = sc)
+        if !isempty(sc.lines)
+            sol[sc.name]["Line overflow (MW)"] =
+                timeseries(model[:overflow], sc.lines, sc = sc)
+        end
+        if !isempty(sc.price_sensitive_loads)
+            sol[sc.name]["Price-sensitive loads (MW)"] =
+                timeseries(model[:loads], sc.price_sensitive_loads, sc = sc)
+        end
+        sol[sc.name]["Spinning reserve (MW)"] = OrderedDict(
+            r.name => OrderedDict(
+                g.name => [
+                    value(model[:reserve][sc.name, r.name, g.name, t]) for t in 1:instance.time
+                ] for g in r.units
+            ) for r in sc.reserves if r.type == "spinning"
+        )
+        sol[sc.name]["Spinning reserve shortfall (MW)"] = OrderedDict(
+            r.name => [
+                value(model[:reserve_shortfall][sc.name, r.name, t]) for
+                t in 1:instance.time
+            ] for r in sc.reserves if r.type == "spinning"
+        )
+        sol[sc.name]["Up-flexiramp (MW)"] = OrderedDict(
+            r.name => OrderedDict(
+                g.name => [
+                    value(model[:upflexiramp][sc.name, r.name, g.name, t]) for t in 1:instance.time
+                ] for g in r.units
+            ) for r in sc.reserves if r.type == "flexiramp"
+        )
+        sol[sc.name]["Up-flexiramp shortfall (MW)"] = OrderedDict(
+            r.name => [
+                value(model[:upflexiramp_shortfall][sc.name, r.name, t]) for t in 1:instance.time
+            ] for r in sc.reserves if r.type == "flexiramp"
+        )
+        sol[sc.name]["Down-flexiramp (MW)"] = OrderedDict(
+            r.name => OrderedDict(
+                g.name => [
+                    value(model[:dwflexiramp][sc.name, r.name, g.name, t]) for t in 1:instance.time
+                ] for g in r.units
+            ) for r in sc.reserves if r.type == "flexiramp"
+        )
+        sol[sc.name]["Down-flexiramp shortfall (MW)"] = OrderedDict(
+            r.name => [
+                value(model[:dwflexiramp_shortfall][sc.name, r.name, t]) for t in 1:instance.time
+            ] for r in sc.reserves if r.type == "flexiramp"
+        )
     end
-    if !isempty(instance.price_sensitive_loads)
+    if length(instance.scenarios) == 1
-        sol["Price-sensitive loads (MW)"] =
+        return first(values(sol))
-            timeseries(model[:loads], instance.price_sensitive_loads)
+    else
+        return sol
     end
-    sol["Spinning reserve (MW)"] = OrderedDict(
-        r.name => OrderedDict(
-            g.name => [
-                value(model[:reserve][r.name, g.name, t]) for
-                t in 1:instance.time
-            ] for g in r.units
-        ) for r in instance.reserves if r.type == "spinning"
-    )
-    sol["Spinning reserve shortfall (MW)"] = OrderedDict(
-        r.name => [
-            value(model[:reserve_shortfall][r.name, t]) for
-            t in 1:instance.time
-        ] for r in instance.reserves if r.type == "spinning"
-    )
-    sol["Up-flexiramp (MW)"] = OrderedDict(
-        r.name => OrderedDict(
-            g.name => [
-                value(model[:upflexiramp][r.name, g.name, t]) for
-                t in 1:instance.time
-            ] for g in r.units
-        ) for r in instance.reserves if r.type == "flexiramp"
-    )
-    sol["Up-flexiramp shortfall (MW)"] = OrderedDict(
-        r.name => [
-            value(model[:upflexiramp_shortfall][r.name, t]) for
-            t in 1:instance.time
-        ] for r in instance.reserves if r.type == "flexiramp"
-    )
-    sol["Down-flexiramp (MW)"] = OrderedDict(
-        r.name => OrderedDict(
-            g.name => [
-                value(model[:dwflexiramp][r.name, g.name, t]) for
-                t in 1:instance.time
-            ] for g in r.units
-        ) for r in instance.reserves if r.type == "flexiramp"
-    )
-    sol["Down-flexiramp shortfall (MW)"] = OrderedDict(
-        r.name => [
-            value(model[:upflexiramp_shortfall][r.name, t]) for
-            t in 1:instance.time
-        ] for r in instance.reserves if r.type == "flexiramp"
-    )
-    return sol
 end

src/transform/initcond.jl

@@ -5,18 +5,18 @@
 using JuMP
 
 """
-    generate_initial_conditions!(instance, optimizer)
+    generate_initial_conditions!(sc, optimizer)
 
-Generates feasible initial conditions for the given instance, by constructing
+Generates feasible initial conditions for the given scenario, by constructing
 and solving a single-period mixed-integer optimization problem, using the given
-optimizer. The instance is modified in-place.
+optimizer. The scenario is modified in-place.
 """
 function generate_initial_conditions!(
-    instance::UnitCommitmentInstance,
+    sc::UnitCommitmentScenario,
     optimizer,
 )::Nothing
-    G = instance.units
+    G = sc.units
-    B = instance.buses
+    B = sc.buses
     t = 1
     mip = JuMP.Model(optimizer)
 

src/transform/randomize/XavQiuAhm2021.jl

@@ -6,6 +6,7 @@ module XavQiuAhm2021
 
 using Distributions
 import ..UnitCommitmentInstance
+import ..UnitCommitmentScenario
 
 """
     struct Randomization
@@ -119,10 +120,10 @@ end
 
 function _randomize_costs(
     rng,
-    instance::UnitCommitmentInstance,
+    sc::UnitCommitmentScenario,
     distribution,
 )::Nothing
-    for unit in instance.units
+    for unit in sc.units
         α = rand(rng, distribution)
         unit.min_power_cost *= α
         for k in unit.cost_segments
@@ -137,17 +138,15 @@ end
 
 function _randomize_load_share(
     rng,
-    instance::UnitCommitmentInstance,
+    sc::UnitCommitmentScenario,
     distribution,
 )::Nothing
-    α = rand(rng, distribution, length(instance.buses))
+    α = rand(rng, distribution, length(sc.buses))
-    for t in 1:instance.time
+    for t in 1:sc.time
-        total = sum(bus.load[t] for bus in instance.buses)
+        total = sum(bus.load[t] for bus in sc.buses)
-        den = sum(
+        den =
-            bus.load[t] / total * α[i] for
+            sum(bus.load[t] / total * α[i] for (i, bus) in enumerate(sc.buses))
-            (i, bus) in enumerate(instance.buses)
+        for (i, bus) in enumerate(sc.buses)
-        )
-        for (i, bus) in enumerate(instance.buses)
             bus.load[t] *= α[i] / den
         end
     end
@@ -156,12 +155,12 @@ end
 
 function _randomize_load_profile(
     rng,
-    instance::UnitCommitmentInstance,
+    sc::UnitCommitmentScenario,
     params::Randomization,
 )::Nothing
     # Generate new system load
     system_load = [1.0]
-    for t in 2:instance.time
+    for t in 2:sc.time
         idx = (t - 1) % length(params.load_profile_mu) + 1
         gamma = rand(
             rng,
@@ -169,14 +168,14 @@ function _randomize_load_profile(
         )
         push!(system_load, system_load[t-1] * gamma)
     end
-    capacity = sum(maximum(u.max_power) for u in instance.units)
+    capacity = sum(maximum(u.max_power) for u in sc.units)
     peak_load = rand(rng, params.peak_load) * capacity
     system_load = system_load ./ maximum(system_load) .* peak_load
 
     # Scale bus loads to match the new system load
-    prev_system_load = sum(b.load for b in instance.buses)
+    prev_system_load = sum(b.load for b in sc.buses)
-    for b in instance.buses
+    for b in sc.buses
-        for t in 1:instance.time
+        for t in 1:sc.time
             b.load[t] *= system_load[t] / prev_system_load[t]
         end
     end
@@ -199,15 +198,26 @@ function randomize!(
     instance::UnitCommitment.UnitCommitmentInstance,
     method::XavQiuAhm2021.Randomization;
     rng = MersenneTwister(),
+)::Nothing
+    for sc in instance.scenarios
+        randomize!(sc, method; rng)
+    end
+    return
+end
+
+function randomize!(
+    sc::UnitCommitment.UnitCommitmentScenario,
+    method::XavQiuAhm2021.Randomization;
+    rng = MersenneTwister(),
 )::Nothing
     if method.randomize_costs
-        XavQiuAhm2021._randomize_costs(rng, instance, method.cost)
+        XavQiuAhm2021._randomize_costs(rng, sc, method.cost)
     end
     if method.randomize_load_share
-        XavQiuAhm2021._randomize_load_share(rng, instance, method.load_share)
+        XavQiuAhm2021._randomize_load_share(rng, sc, method.load_share)
     end
     if method.randomize_load_profile
-        XavQiuAhm2021._randomize_load_profile(rng, instance, method)
+        XavQiuAhm2021._randomize_load_profile(rng, sc, method)
     end
     return
 end

src/transform/slice.jl

@@ -24,31 +24,33 @@ function slice(
 )::UnitCommitmentInstance
     modified = deepcopy(instance)
     modified.time = length(range)
-    modified.power_balance_penalty = modified.power_balance_penalty[range]
+    for sc in modified.scenarios
-    for r in modified.reserves
+        sc.power_balance_penalty = sc.power_balance_penalty[range]
-        r.amount = r.amount[range]
+        for r in sc.reserves
-    end
+            r.amount = r.amount[range]
-    for u in modified.units
+        end
-        u.max_power = u.max_power[range]
+        for u in sc.units
-        u.min_power = u.min_power[range]
+            u.max_power = u.max_power[range]
-        u.must_run = u.must_run[range]
+            u.min_power = u.min_power[range]
-        u.min_power_cost = u.min_power_cost[range]
+            u.must_run = u.must_run[range]
-        for s in u.cost_segments
+            u.min_power_cost = u.min_power_cost[range]
-            s.mw = s.mw[range]
+            for s in u.cost_segments
-            s.cost = s.cost[range]
+                s.mw = s.mw[range]
+                s.cost = s.cost[range]
+            end
+        end
+        for b in sc.buses
+            b.load = b.load[range]
+        end
+        for l in sc.lines
+            l.normal_flow_limit = l.normal_flow_limit[range]
+            l.emergency_flow_limit = l.emergency_flow_limit[range]
+            l.flow_limit_penalty = l.flow_limit_penalty[range]
+        end
+        for ps in sc.price_sensitive_loads
+            ps.demand = ps.demand[range]
+            ps.revenue = ps.revenue[range]
         end
-    end
-    for b in modified.buses
-        b.load = b.load[range]
-    end
-    for l in modified.lines
-        l.normal_flow_limit = l.normal_flow_limit[range]
-        l.emergency_flow_limit = l.emergency_flow_limit[range]
-        l.flow_limit_penalty = l.flow_limit_penalty[range]
-    end
-    for ps in modified.price_sensitive_loads
-        ps.demand = ps.demand[range]
-        ps.revenue = ps.revenue[range]
     end
     return modified
 end

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]

src/validation/validate.jl

@@ -28,6 +28,8 @@ function validate(
     instance::UnitCommitmentInstance,
     solution::Union{Dict,OrderedDict},
 )::Bool
+    "Production (MW)" ∈ keys(solution) ? solution = Dict("s1" => solution) :
+    nothing
     err_count = 0
     err_count += _validate_units(instance, solution)
     err_count += _validate_reserve_and_demand(instance, solution)
@@ -42,358 +44,369 @@ end
 
 function _validate_units(instance::UnitCommitmentInstance, solution; tol = 0.01)
     err_count = 0
-
+    for sc in instance.scenarios
-    for unit in instance.units
+        for unit in sc.units
-        production = solution["Production (MW)"][unit.name]
+            production = solution[sc.name]["Production (MW)"][unit.name]
-        reserve = [0.0 for _ in 1:instance.time]
+            reserve = [0.0 for _ in 1:instance.time]
-        spinning_reserves = [r for r in unit.reserves if r.type == "spinning"]
+            spinning_reserves =
-        if !isempty(spinning_reserves)
+                [r for r in unit.reserves if r.type == "spinning"]
-            reserve += sum(
+            if !isempty(spinning_reserves)
-                solution["Spinning reserve (MW)"][r.name][unit.name] for
+                reserve += sum(
-                r in spinning_reserves
+                    solution[sc.name]["Spinning reserve (MW)"][r.name][unit.name]
-            )
+                    for r in spinning_reserves
-        end
+                )
-        actual_production_cost = solution["Production cost (\$)"][unit.name]
-        actual_startup_cost = solution["Startup cost (\$)"][unit.name]
-        is_on = bin(solution["Is on"][unit.name])
-
-        for t in 1:instance.time
-            # Auxiliary variables
-            if t == 1
-                is_starting_up = (unit.initial_status < 0) && is_on[t]
-                is_shutting_down = (unit.initial_status > 0) && !is_on[t]
-                ramp_up =
-                    max(0, production[t] + reserve[t] - unit.initial_power)
-                ramp_down = max(0, unit.initial_power - production[t])
-            else
-                is_starting_up = !is_on[t-1] && is_on[t]
-                is_shutting_down = is_on[t-1] && !is_on[t]
-                ramp_up = max(0, production[t] + reserve[t] - production[t-1])
-                ramp_down = max(0, production[t-1] - production[t])
             end
+            actual_production_cost =
+                solution[sc.name]["Production cost (\$)"][unit.name]
+            actual_startup_cost =
+                solution[sc.name]["Startup cost (\$)"][unit.name]
+            is_on = bin(solution[sc.name]["Is on"][unit.name])
 
-            # Compute production costs
+            for t in 1:instance.time
-            production_cost, startup_cost = 0, 0
+                # Auxiliary variables
-            if is_on[t]
+                if t == 1
-                production_cost += unit.min_power_cost[t]
+                    is_starting_up = (unit.initial_status < 0) && is_on[t]
-                residual = max(0, production[t] - unit.min_power[t])
+                    is_shutting_down = (unit.initial_status > 0) && !is_on[t]
-                for s in unit.cost_segments
+                    ramp_up =
-                    cleared = min(residual, s.mw[t])
+                        max(0, production[t] + reserve[t] - unit.initial_power)
-                    production_cost += cleared * s.cost[t]
+                    ramp_down = max(0, unit.initial_power - production[t])
-                    residual = max(0, residual - s.mw[t])
+                else
+                    is_starting_up = !is_on[t-1] && is_on[t]
+                    is_shutting_down = is_on[t-1] && !is_on[t]
+                    ramp_up =
+                        max(0, production[t] + reserve[t] - production[t-1])
+                    ramp_down = max(0, production[t-1] - production[t])
                 end
-            end
 
-            # Production should be non-negative
+                # Compute production costs
-            if production[t] < -tol
+                production_cost, startup_cost = 0, 0
-                @error @sprintf(
+                if is_on[t]
-                    "Unit %s produces negative amount of power at time %d (%.2f)",
+                    production_cost += unit.min_power_cost[t]
-                    unit.name,
+                    residual = max(0, production[t] - unit.min_power[t])
-                    t,
+                    for s in unit.cost_segments
-                    production[t]
+                        cleared = min(residual, s.mw[t])
-                )
+                        production_cost += cleared * s.cost[t]
-                err_count += 1
+                        residual = max(0, residual - s.mw[t])
-            end
+                    end
+                end
 
-            # Verify must-run
+                # Production should be non-negative
-            if !is_on[t] && unit.must_run[t]
+                if production[t] < -tol
-                @error @sprintf(
+                    @error @sprintf(
-                    "Must-run unit %s is offline at time %d",
+                        "Unit %s produces negative amount of power at time %d (%.2f)",
-                    unit.name,
+                        unit.name,
-                    t
+                        t,
-                )
+                        production[t]
-                err_count += 1
+                    )
-            end
+                    err_count += 1
+                end
 
-            # Verify reserve eligibility
+                # Verify must-run
-            for r in instance.reserves
+                if !is_on[t] && unit.must_run[t]
-                if r.type == "spinning"
+                    @error @sprintf(
-                    if unit ∉ r.units &&
+                        "Must-run unit %s is offline at time %d",
-                       (unit in keys(solution["Spinning reserve (MW)"][r.name]))
+                        unit.name,
+                        t
+                    )
+                    err_count += 1
+                end
+
+                # Verify reserve eligibility
+                for r in sc.reserves
+                    if r.type == "spinning"
+                        if unit ∉ r.units && (
+                            unit in keys(
+                                solution[sc.name]["Spinning reserve (MW)"][r.name],
+                            )
+                        )
+                            @error @sprintf(
+                                "Unit %s is not eligible to provide reserve %s",
+                                unit.name,
+                                r.name,
+                            )
+                            err_count += 1
+                        end
+                    end
+                end
+
+                # If unit is on, must produce at least its minimum power
+                if is_on[t] && (production[t] < unit.min_power[t] - tol)
+                    @error @sprintf(
+                        "Unit %s produces below its minimum limit at time %d (%.2f < %.2f)",
+                        unit.name,
+                        t,
+                        production[t],
+                        unit.min_power[t]
+                    )
+                    err_count += 1
+                end
+
+                # If unit is on, must produce at most its maximum power
+                if is_on[t] &&
+                   (production[t] + reserve[t] > unit.max_power[t] + tol)
+                    @error @sprintf(
+                        "Unit %s produces above its maximum limit at time %d (%.2f + %.2f> %.2f)",
+                        unit.name,
+                        t,
+                        production[t],
+                        reserve[t],
+                        unit.max_power[t]
+                    )
+                    err_count += 1
+                end
+
+                # If unit is off, must produce zero
+                if !is_on[t] && production[t] + reserve[t] > tol
+                    @error @sprintf(
+                        "Unit %s produces power at time %d while off (%.2f + %.2f > 0)",
+                        unit.name,
+                        t,
+                        production[t],
+                        reserve[t],
+                    )
+                    err_count += 1
+                end
+
+                # Startup limit
+                if is_starting_up && (ramp_up > unit.startup_limit + tol)
+                    @error @sprintf(
+                        "Unit %s exceeds startup limit at time %d (%.2f > %.2f)",
+                        unit.name,
+                        t,
+                        ramp_up,
+                        unit.startup_limit
+                    )
+                    err_count += 1
+                end
+
+                # Shutdown limit
+                if is_shutting_down && (ramp_down > unit.shutdown_limit + tol)
+                    @error @sprintf(
+                        "Unit %s exceeds shutdown limit at time %d (%.2f > %.2f)",
+                        unit.name,
+                        t,
+                        ramp_down,
+                        unit.shutdown_limit
+                    )
+                    err_count += 1
+                end
+
+                # Ramp-up limit
+                if !is_starting_up &&
+                   !is_shutting_down &&
+                   (ramp_up > unit.ramp_up_limit + tol)
+                    @error @sprintf(
+                        "Unit %s exceeds ramp up limit at time %d (%.2f > %.2f)",
+                        unit.name,
+                        t,
+                        ramp_up,
+                        unit.ramp_up_limit
+                    )
+                    err_count += 1
+                end
+
+                # Ramp-down limit
+                if !is_starting_up &&
+                   !is_shutting_down &&
+                   (ramp_down > unit.ramp_down_limit + tol)
+                    @error @sprintf(
+                        "Unit %s exceeds ramp down limit at time %d (%.2f > %.2f)",
+                        unit.name,
+                        t,
+                        ramp_down,
+                        unit.ramp_down_limit
+                    )
+                    err_count += 1
+                end
+
+                # Verify startup costs & minimum downtime
+                if is_starting_up
+
+                    # Calculate how much time the unit has been offline
+                    time_down = 0
+                    for k in 1:(t-1)
+                        if !is_on[t-k]
+                            time_down += 1
+                        else
+                            break
+                        end
+                    end
+                    if (t == time_down + 1) && (unit.initial_status < 0)
+                        time_down -= unit.initial_status
+                    end
+
+                    # Calculate startup costs
+                    for c in unit.startup_categories
+                        if time_down >= c.delay
+                            startup_cost = c.cost
+                        end
+                    end
+
+                    # Check minimum downtime
+                    if time_down < unit.min_downtime
                         @error @sprintf(
-                            "Unit %s is not eligible to provide reserve %s",
+                            "Unit %s violates minimum downtime at time %d",
                             unit.name,
-                            r.name,
+                            t
                         )
                         err_count += 1
                     end
                 end
-            end
 
-            # If unit is on, must produce at least its minimum power
+                # Verify minimum uptime
-            if is_on[t] && (production[t] < unit.min_power[t] - tol)
+                if is_shutting_down
-                @error @sprintf(
-                    "Unit %s produces below its minimum limit at time %d (%.2f < %.2f)",
-                    unit.name,
-                    t,
-                    production[t],
-                    unit.min_power[t]
-                )
-                err_count += 1
-            end
 
-            # If unit is on, must produce at most its maximum power
+                    # Calculate how much time the unit has been online
-            if is_on[t] &&
+                    time_up = 0
-               (production[t] + reserve[t] > unit.max_power[t] + tol)
+                    for k in 1:(t-1)
-                @error @sprintf(
+                        if is_on[t-k]
-                    "Unit %s produces above its maximum limit at time %d (%.2f + %.2f> %.2f)",
+                            time_up += 1
-                    unit.name,
+                        else
-                    t,
+                            break
-                    production[t],
+                        end
-                    reserve[t],
+                    end
-                    unit.max_power[t]
+                    if (t == time_up + 1) && (unit.initial_status > 0)
-                )
+                        time_up += unit.initial_status
-                err_count += 1
-            end
-
-            # If unit is off, must produce zero
-            if !is_on[t] && production[t] + reserve[t] > tol
-                @error @sprintf(
-                    "Unit %s produces power at time %d while off (%.2f + %.2f > 0)",
-                    unit.name,
-                    t,
-                    production[t],
-                    reserve[t],
-                )
-                err_count += 1
-            end
-
-            # Startup limit
-            if is_starting_up && (ramp_up > unit.startup_limit + tol)
-                @error @sprintf(
-                    "Unit %s exceeds startup limit at time %d (%.2f > %.2f)",
-                    unit.name,
-                    t,
-                    ramp_up,
-                    unit.startup_limit
-                )
-                err_count += 1
-            end
-
-            # Shutdown limit
-            if is_shutting_down && (ramp_down > unit.shutdown_limit + tol)
-                @error @sprintf(
-                    "Unit %s exceeds shutdown limit at time %d (%.2f > %.2f)",
-                    unit.name,
-                    t,
-                    ramp_down,
-                    unit.shutdown_limit
-                )
-                err_count += 1
-            end
-
-            # Ramp-up limit
-            if !is_starting_up &&
-               !is_shutting_down &&
-               (ramp_up > unit.ramp_up_limit + tol)
-                @error @sprintf(
-                    "Unit %s exceeds ramp up limit at time %d (%.2f > %.2f)",
-                    unit.name,
-                    t,
-                    ramp_up,
-                    unit.ramp_up_limit
-                )
-                err_count += 1
-            end
-
-            # Ramp-down limit
-            if !is_starting_up &&
-               !is_shutting_down &&
-               (ramp_down > unit.ramp_down_limit + tol)
-                @error @sprintf(
-                    "Unit %s exceeds ramp down limit at time %d (%.2f > %.2f)",
-                    unit.name,
-                    t,
-                    ramp_down,
-                    unit.ramp_down_limit
-                )
-                err_count += 1
-            end
-
-            # Verify startup costs & minimum downtime
-            if is_starting_up
-
-                # Calculate how much time the unit has been offline
-                time_down = 0
-                for k in 1:(t-1)
-                    if !is_on[t-k]
-                        time_down += 1
-                    else
-                        break
                     end
-                end
-                if (t == time_down + 1) && (unit.initial_status < 0)
-                    time_down -= unit.initial_status
-                end
 
-                # Calculate startup costs
+                    # Check minimum uptime
-                for c in unit.startup_categories
+                    if time_up < unit.min_uptime
-                    if time_down >= c.delay
+                        @error @sprintf(
-                        startup_cost = c.cost
+                            "Unit %s violates minimum uptime at time %d",
+                            unit.name,
+                            t
+                        )
+                        err_count += 1
                     end
                 end
 
-                # Check minimum downtime
+                # Verify production costs
-                if time_down < unit.min_downtime
+                if abs(actual_production_cost[t] - production_cost) > 1.00
                     @error @sprintf(
-                        "Unit %s violates minimum downtime at time %d",
+                        "Unit %s has unexpected production cost at time %d (%.2f should be %.2f)",
                         unit.name,
-                        t
+                        t,
+                        actual_production_cost[t],
+                        production_cost
                     )
                     err_count += 1
                 end
-            end
 
-            # Verify minimum uptime
+                # Verify startup costs
-            if is_shutting_down
+                if abs(actual_startup_cost[t] - startup_cost) > 1.00
-
-                # Calculate how much time the unit has been online
-                time_up = 0
-                for k in 1:(t-1)
-                    if is_on[t-k]
-                        time_up += 1
-                    else
-                        break
-                    end
-                end
-                if (t == time_up + 1) && (unit.initial_status > 0)
-                    time_up += unit.initial_status
-                end
-
-                # Check minimum uptime
-                if time_up < unit.min_uptime
                     @error @sprintf(
-                        "Unit %s violates minimum uptime at time %d",
+                        "Unit %s has unexpected startup cost at time %d (%.2f should be %.2f)",
                         unit.name,
-                        t
+                        t,
+                        actual_startup_cost[t],
+                        startup_cost
                     )
                     err_count += 1
                 end
             end
-
-            # Verify production costs
-            if abs(actual_production_cost[t] - production_cost) > 1.00
-                @error @sprintf(
-                    "Unit %s has unexpected production cost at time %d (%.2f should be %.2f)",
-                    unit.name,
-                    t,
-                    actual_production_cost[t],
-                    production_cost
-                )
-                err_count += 1
-            end
-
-            # Verify startup costs
-            if abs(actual_startup_cost[t] - startup_cost) > 1.00
-                @error @sprintf(
-                    "Unit %s has unexpected startup cost at time %d (%.2f should be %.2f)",
-                    unit.name,
-                    t,
-                    actual_startup_cost[t],
-                    startup_cost
-                )
-                err_count += 1
-            end
         end
     end
-
     return err_count
 end
 
 function _validate_reserve_and_demand(instance, solution, tol = 0.01)
     err_count = 0
-    for t in 1:instance.time
+    for sc in instance.scenarios
-        load_curtail = 0
+        for t in 1:instance.time
-        fixed_load = sum(b.load[t] for b in instance.buses)
+            load_curtail = 0
-        ps_load = 0
+            fixed_load = sum(b.load[t] for b in sc.buses)
-        if length(instance.price_sensitive_loads) > 0
+            ps_load = 0
-            ps_load = sum(
+            if length(sc.price_sensitive_loads) > 0
-                solution["Price-sensitive loads (MW)"][ps.name][t] for
+                ps_load = sum(
-                ps in instance.price_sensitive_loads
+                    solution[sc.name]["Price-sensitive loads (MW)"][ps.name][t]
-            )
+                    for ps in sc.price_sensitive_loads
-        end
-        production =
-            sum(solution["Production (MW)"][g.name][t] for g in instance.units)
-        if "Load curtail (MW)" in keys(solution)
-            load_curtail = sum(
-                solution["Load curtail (MW)"][b.name][t] for
-                b in instance.buses
-            )
-        end
-        balance = fixed_load - load_curtail - production + ps_load
-
-        # Verify that production equals demand
-        if abs(balance) > tol
-            @error @sprintf(
-                "Non-zero power balance at time %d (%.2f + %.2f - %.2f - %.2f != 0)",
-                t,
-                fixed_load,
-                ps_load,
-                load_curtail,
-                production,
-            )
-            err_count += 1
-        end
-
-        # Verify reserves
-        for r in instance.reserves
-            if r.type == "spinning"
-                provided = sum(
-                    solution["Spinning reserve (MW)"][r.name][g.name][t] for
-                    g in r.units
                 )
-                shortfall =
+            end
-                    solution["Spinning reserve shortfall (MW)"][r.name][t]
+            production = sum(
-                required = r.amount[t]
+                solution[sc.name]["Production (MW)"][g.name][t] for
-
+                g in sc.units
-                if provided + shortfall < required - tol
+            )
-                    @error @sprintf(
+            if "Load curtail (MW)" in keys(solution)
-                        "Insufficient reserve %s at time %d (%.2f + %.2f < %.2f)",
+                load_curtail = sum(
-                        r.name,
+                    solution[sc.name]["Load curtail (MW)"][b.name][t] for
-                        t,
+                    b in sc.buses
-                        provided,
-                        shortfall,
-                        required,
-                    )
-                end
-            elseif r.type == "flexiramp"
-                upflexiramp = sum(
-                    solution["Up-flexiramp (MW)"][r.name][g.name][t] for
-                    g in r.units
                 )
-                upflexiramp_shortfall =
+            end
-                    solution["Up-flexiramp shortfall (MW)"][r.name][t]
+            balance = fixed_load - load_curtail - production + ps_load
 
-                if upflexiramp + upflexiramp_shortfall < r.amount[t] - tol
+            # Verify that production equals demand
-                    @error @sprintf(
+            if abs(balance) > tol
-                        "Insufficient up-flexiramp at time %d (%.2f + %.2f < %.2f)",
+                @error @sprintf(
-                        t,
+                    "Non-zero power balance at time %d (%.2f + %.2f - %.2f - %.2f != 0)",
-                        upflexiramp,
+                    t,
-                        upflexiramp_shortfall,
+                    fixed_load,
-                        r.amount[t],
+                    ps_load,
-                    )
+                    load_curtail,
-                    err_count += 1
+                    production,
-                end
-
-                dwflexiramp = sum(
-                    solution["Down-flexiramp (MW)"][r.name][g.name][t] for
-                    g in r.units
                 )
-                dwflexiramp_shortfall =
+                err_count += 1
-                    solution["Down-flexiramp shortfall (MW)"][r.name][t]
+            end
 
-                if dwflexiramp + dwflexiramp_shortfall < r.amount[t] - tol
+            # Verify reserves
-                    @error @sprintf(
+            for r in sc.reserves
-                        "Insufficient down-flexiramp at time %d (%.2f + %.2f < %.2f)",
+                if r.type == "spinning"
-                        t,
+                    provided = sum(
-                        dwflexiramp,
+                        solution[sc.name]["Spinning reserve (MW)"][r.name][g.name][t]
-                        dwflexiramp_shortfall,
+                        for g in r.units
-                        r.amount[t],
                     )
-                    err_count += 1
+                    shortfall =
+                        solution[sc.name]["Spinning reserve shortfall (MW)"][r.name][t]
+                    required = r.amount[t]
+
+                    if provided + shortfall < required - tol
+                        @error @sprintf(
+                            "Insufficient reserve %s at time %d (%.2f + %.2f < %.2f)",
+                            r.name,
+                            t,
+                            provided,
+                            shortfall,
+                            required,
+                        )
+                    end
+                elseif r.type == "flexiramp"
+                    upflexiramp = sum(
+                        solution[sc.name]["Up-flexiramp (MW)"][r.name][g.name][t]
+                        for g in r.units
+                    )
+                    upflexiramp_shortfall =
+                        solution[sc.name]["Up-flexiramp shortfall (MW)"][r.name][t]
+
+                    if upflexiramp + upflexiramp_shortfall < r.amount[t] - tol
+                        @error @sprintf(
+                            "Insufficient up-flexiramp at time %d (%.2f + %.2f < %.2f)",
+                            t,
+                            upflexiramp,
+                            upflexiramp_shortfall,
+                            r.amount[t],
+                        )
+                        err_count += 1
+                    end
+
+                    dwflexiramp = sum(
+                        solution[sc.name]["Down-flexiramp (MW)"][r.name][g.name][t]
+                        for g in r.units
+                    )
+                    dwflexiramp_shortfall =
+                        solution[sc.name]["Down-flexiramp shortfall (MW)"][r.name][t]
+
+                    if dwflexiramp + dwflexiramp_shortfall < r.amount[t] - tol
+                        @error @sprintf(
+                            "Insufficient down-flexiramp at time %d (%.2f + %.2f < %.2f)",
+                            t,
+                            dwflexiramp,
+                            dwflexiramp_shortfall,
+                            r.amount[t],
+                        )
+                        err_count += 1
+                    end
+                else
+                    error("Unknown reserve type: $(r.type)")
                 end
-            else
-                error("Unknown reserve type: $(r.type)")
             end
         end
     end

test/instance/migrate_test.jl

@@ -6,13 +6,17 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
 
 @testset "read v0.2" begin
     instance = UnitCommitment.read("$FIXTURES/ucjl-0.2.json.gz")
-    @test length(instance.reserves_by_name["r1"].amount) == 4
+    @test length(instance.scenarios) == 1
-    @test instance.units_by_name["g2"].reserves[1].name == "r1"
+    sc = instance.scenarios[1]
+    @test length(sc.reserves_by_name["r1"].amount) == 4
+    @test sc.units_by_name["g2"].reserves[1].name == "r1"
 end
 
 @testset "read v0.3" begin
     instance = UnitCommitment.read("$FIXTURES/ucjl-0.3.json.gz")
-    @test length(instance.units) == 6
+    @test length(instance.scenarios) == 1
-    @test length(instance.buses) == 14
+    sc = instance.scenarios[1]
-    @test length(instance.lines) == 20
+    @test length(sc.units) == 6
+    @test length(sc.buses) == 14
+    @test length(sc.lines) == 20
 end

test/instance/read_test.jl

@@ -7,42 +7,49 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
 @testset "read_benchmark" begin
     instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
 
-    @test length(instance.lines) == 20
+    @test repr(instance) == (
-    @test length(instance.buses) == 14
+        "UnitCommitmentInstance(1 scenarios, 6 units, 14 buses, " *
-    @test length(instance.units) == 6
+        "20 lines, 19 contingencies, 1 price sensitive loads, 4 time steps)"
-    @test length(instance.contingencies) == 19
+    )
-    @test length(instance.price_sensitive_loads) == 1
+
+    @test length(instance.scenarios) == 1
+    sc = instance.scenarios[1]
+    @test length(sc.lines) == 20
+    @test length(sc.buses) == 14
+    @test length(sc.units) == 6
+    @test length(sc.contingencies) == 19
+    @test length(sc.price_sensitive_loads) == 1
     @test instance.time == 4
 
-    @test instance.lines[5].name == "l5"
+    @test sc.lines[5].name == "l5"
-    @test instance.lines[5].source.name == "b2"
+    @test sc.lines[5].source.name == "b2"
-    @test instance.lines[5].target.name == "b5"
+    @test sc.lines[5].target.name == "b5"
-    @test instance.lines[5].reactance ≈ 0.17388
+    @test sc.lines[5].reactance ≈ 0.17388
-    @test instance.lines[5].susceptance ≈ 10.037550333
+    @test sc.lines[5].susceptance ≈ 10.037550333
-    @test instance.lines[5].normal_flow_limit == [1e8 for t in 1:4]
+    @test sc.lines[5].normal_flow_limit == [1e8 for t in 1:4]
-    @test instance.lines[5].emergency_flow_limit == [1e8 for t in 1:4]
+    @test sc.lines[5].emergency_flow_limit == [1e8 for t in 1:4]
-    @test instance.lines[5].flow_limit_penalty == [5e3 for t in 1:4]
+    @test sc.lines[5].flow_limit_penalty == [5e3 for t in 1:4]
-    @test instance.lines_by_name["l5"].name == "l5"
+    @test sc.lines_by_name["l5"].name == "l5"
 
-    @test instance.lines[1].name == "l1"
+    @test sc.lines[1].name == "l1"
-    @test instance.lines[1].source.name == "b1"
+    @test sc.lines[1].source.name == "b1"
-    @test instance.lines[1].target.name == "b2"
+    @test sc.lines[1].target.name == "b2"
-    @test instance.lines[1].reactance ≈ 0.059170
+    @test sc.lines[1].reactance ≈ 0.059170
-    @test instance.lines[1].susceptance ≈ 29.496860773945
+    @test sc.lines[1].susceptance ≈ 29.496860773945
-    @test instance.lines[1].normal_flow_limit == [300.0 for t in 1:4]
+    @test sc.lines[1].normal_flow_limit == [300.0 for t in 1:4]
-    @test instance.lines[1].emergency_flow_limit == [400.0 for t in 1:4]
+    @test sc.lines[1].emergency_flow_limit == [400.0 for t in 1:4]
-    @test instance.lines[1].flow_limit_penalty == [1e3 for t in 1:4]
+    @test sc.lines[1].flow_limit_penalty == [1e3 for t in 1:4]
 
-    @test instance.buses[9].name == "b9"
+    @test sc.buses[9].name == "b9"
-    @test instance.buses[9].load == [35.36638, 33.25495, 31.67138, 31.14353]
+    @test sc.buses[9].load == [35.36638, 33.25495, 31.67138, 31.14353]
-    @test instance.buses_by_name["b9"].name == "b9"
+    @test sc.buses_by_name["b9"].name == "b9"
 
-    @test instance.reserves[1].name == "r1"
+    @test sc.reserves[1].name == "r1"
-    @test instance.reserves[1].type == "spinning"
+    @test sc.reserves[1].type == "spinning"
-    @test instance.reserves[1].amount == [100.0, 100.0, 100.0, 100.0]
+    @test sc.reserves[1].amount == [100.0, 100.0, 100.0, 100.0]
-    @test instance.reserves_by_name["r1"].name == "r1"
+    @test sc.reserves_by_name["r1"].name == "r1"
 
-    unit = instance.units[1]
+    unit = sc.units[1]
     @test unit.name == "g1"
     @test unit.bus.name == "b1"
     @test unit.ramp_up_limit == 1e6
@@ -69,14 +76,14 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
     @test unit.startup_categories[2].cost == 1500.0
     @test unit.startup_categories[3].cost == 2000.0
     @test length(unit.reserves) == 0
-    @test instance.units_by_name["g1"].name == "g1"
+    @test sc.units_by_name["g1"].name == "g1"
 
-    unit = instance.units[2]
+    unit = sc.units[2]
     @test unit.name == "g2"
     @test unit.must_run == [false for t in 1:4]
     @test length(unit.reserves) == 1
 
-    unit = instance.units[3]
+    unit = sc.units[3]
     @test unit.name == "g3"
     @test unit.bus.name == "b3"
     @test unit.ramp_up_limit == 70.0
@@ -98,23 +105,23 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
     @test length(unit.reserves) == 1
     @test unit.reserves[1].name == "r1"
 
-    @test instance.contingencies[1].lines == [instance.lines[1]]
+    @test sc.contingencies[1].lines == [sc.lines[1]]
-    @test instance.contingencies[1].units == []
+    @test sc.contingencies[1].units == []
-    @test instance.contingencies[1].name == "c1"
+    @test sc.contingencies[1].name == "c1"
-    @test instance.contingencies_by_name["c1"].name == "c1"
+    @test sc.contingencies_by_name["c1"].name == "c1"
 
-    load = instance.price_sensitive_loads[1]
+    load = sc.price_sensitive_loads[1]
     @test load.name == "ps1"
     @test load.bus.name == "b3"
     @test load.revenue == [100.0 for t in 1:4]
     @test load.demand == [50.0 for t in 1:4]
-    @test instance.price_sensitive_loads_by_name["ps1"].name == "ps1"
+    @test sc.price_sensitive_loads_by_name["ps1"].name == "ps1"
 end
 
 @testset "read_benchmark sub-hourly" begin
     instance = UnitCommitment.read("$FIXTURES/case14-sub-hourly.json.gz")
     @test instance.time == 4
-    unit = instance.units[1]
+    unit = instance.scenarios[1].units[1]
     @test unit.name == "g1"
     @test unit.min_uptime == 2
     @test unit.min_downtime == 2

test/solution/methods/XavQiuWanThi19/filter_test.jl

@@ -7,13 +7,13 @@ import UnitCommitment: _Violation, _offer, _query
 
 @testset "_ViolationFilter" begin
     instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
+    sc = instance.scenarios[1]
     filter = UnitCommitment._ViolationFilter(max_per_line = 1, max_total = 2)
-
     _offer(
         filter,
         _Violation(
             time = 1,
-            monitored_line = instance.lines[1],
+            monitored_line = sc.lines[1],
             outage_line = nothing,
             amount = 100.0,
         ),
@@ -22,8 +22,8 @@ import UnitCommitment: _Violation, _offer, _query
         filter,
         _Violation(
             time = 1,
-            monitored_line = instance.lines[1],
+            monitored_line = sc.lines[1],
-            outage_line = instance.lines[1],
+            outage_line = sc.lines[1],
             amount = 300.0,
         ),
     )
@@ -31,8 +31,8 @@ import UnitCommitment: _Violation, _offer, _query
         filter,
         _Violation(
             time = 1,
-            monitored_line = instance.lines[1],
+            monitored_line = sc.lines[1],
-            outage_line = instance.lines[5],
+            outage_line = sc.lines[5],
             amount = 500.0,
         ),
     )
@@ -40,8 +40,8 @@ import UnitCommitment: _Violation, _offer, _query
         filter,
         _Violation(
             time = 1,
-            monitored_line = instance.lines[1],
+            monitored_line = sc.lines[1],
-            outage_line = instance.lines[4],
+            outage_line = sc.lines[4],
             amount = 400.0,
         ),
     )
@@ -49,8 +49,8 @@ import UnitCommitment: _Violation, _offer, _query
         filter,
         _Violation(
             time = 1,
-            monitored_line = instance.lines[2],
+            monitored_line = sc.lines[2],
-            outage_line = instance.lines[1],
+            outage_line = sc.lines[1],
             amount = 200.0,
         ),
     )
@@ -58,8 +58,8 @@ import UnitCommitment: _Violation, _offer, _query
         filter,
         _Violation(
             time = 1,
-            monitored_line = instance.lines[2],
+            monitored_line = sc.lines[2],
-            outage_line = instance.lines[8],
+            outage_line = sc.lines[8],
             amount = 100.0,
         ),
     )
@@ -68,14 +68,14 @@ import UnitCommitment: _Violation, _offer, _query
     expected = [
         _Violation(
             time = 1,
-            monitored_line = instance.lines[2],
+            monitored_line = sc.lines[2],
-            outage_line = instance.lines[1],
+            outage_line = sc.lines[1],
             amount = 200.0,
         ),
         _Violation(
             time = 1,
-            monitored_line = instance.lines[1],
+            monitored_line = sc.lines[1],
-            outage_line = instance.lines[5],
+            outage_line = sc.lines[5],
             amount = 500.0,
         ),
     ]

test/solution/methods/XavQiuWanThi19/find_test.jl

@@ -7,23 +7,25 @@ import UnitCommitment: _Violation, _offer, _query
 
 @testset "find_violations" begin
     instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
-    for line in instance.lines, t in 1:instance.time
+    sc = instance.scenarios[1]
+    for line in sc.lines, t in 1:instance.time
         line.normal_flow_limit[t] = 1.0
         line.emergency_flow_limit[t] = 1.0
     end
     isf = UnitCommitment._injection_shift_factors(
-        lines = instance.lines,
+        lines = sc.lines,
-        buses = instance.buses,
+        buses = sc.buses,
     )
     lodf = UnitCommitment._line_outage_factors(
-        lines = instance.lines,
+        lines = sc.lines,
-        buses = instance.buses,
+        buses = sc.buses,
         isf = isf,
     )
     inj = [1000.0 for b in 1:13, t in 1:instance.time]
-    overflow = [0.0 for l in instance.lines, t in 1:instance.time]
+    overflow = [0.0 for l in sc.lines, t in 1:instance.time]
     violations = UnitCommitment._find_violations(
         instance = instance,
+        sc = sc,
         net_injections = inj,
         overflow = overflow,
         isf = isf,

test/solution/methods/XavQiuWanThi19/sensitivity_test.jl

@@ -6,7 +6,8 @@ using UnitCommitment, Test, LinearAlgebra
 
 @testset "_susceptance_matrix" begin
     instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
-    actual = UnitCommitment._susceptance_matrix(instance.lines)
+    sc = instance.scenarios[1]
+    actual = UnitCommitment._susceptance_matrix(sc.lines)
     @test size(actual) == (20, 20)
     expected = Diagonal([
         29.5,
@@ -35,9 +36,10 @@ end
 
 @testset "_reduced_incidence_matrix" begin
     instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
+    sc = instance.scenarios[1]
     actual = UnitCommitment._reduced_incidence_matrix(
-        lines = instance.lines,
+        lines = sc.lines,
-        buses = instance.buses,
+        buses = sc.buses,
     )
     @test size(actual) == (20, 13)
     @test actual[1, 1] == -1.0
@@ -82,9 +84,10 @@ end
 
 @testset "_injection_shift_factors" begin
     instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
+    sc = instance.scenarios[1]
     actual = UnitCommitment._injection_shift_factors(
-        lines = instance.lines,
+        lines = sc.lines,
-        buses = instance.buses,
+        buses = sc.buses,
     )
     @test size(actual) == (20, 13)
     @test round.(actual, digits = 2) == [
@@ -113,25 +116,26 @@ end
 
 @testset "_line_outage_factors" begin
     instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
+    sc = instance.scenarios[1]
     isf_before = UnitCommitment._injection_shift_factors(
-        lines = instance.lines,
+        lines = sc.lines,
-        buses = instance.buses,
+        buses = sc.buses,
     )
     lodf = UnitCommitment._line_outage_factors(
-        lines = instance.lines,
+        lines = sc.lines,
-        buses = instance.buses,
+        buses = sc.buses,
         isf = isf_before,
     )
-    for contingency in instance.contingencies
+    for contingency in sc.contingencies
         for lc in contingency.lines
             prev_susceptance = lc.susceptance
             lc.susceptance = 0.0
             isf_after = UnitCommitment._injection_shift_factors(
-                lines = instance.lines,
+                lines = sc.lines,
-                buses = instance.buses,
+                buses = sc.buses,
             )
             lc.susceptance = prev_susceptance
-            for lm in instance.lines
+            for lm in sc.lines
                 expected = isf_after[lm.offset, :]
                 actual =
                     isf_before[lm.offset, :] +

test/transform/initcond_test.jl

@@ -8,18 +8,18 @@ using UnitCommitment, Cbc, JuMP
     # Load instance
     instance = UnitCommitment.read("$FIXTURES/case118-initcond.json.gz")
     optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
-
+    sc = instance.scenarios[1]
     # All units should have unknown initial conditions
-    for g in instance.units
+    for g in sc.units
         @test g.initial_power === nothing
         @test g.initial_status === nothing
     end
 
     # Generate initial conditions
-    UnitCommitment.generate_initial_conditions!(instance, optimizer)
+    UnitCommitment.generate_initial_conditions!(sc, optimizer)
 
     # All units should now have known initial conditions
-    for g in instance.units
+    for g in sc.units
         @test g.initial_power !== nothing
         @test g.initial_status !== nothing
     end

test/transform/randomize/XavQiuAhm2021_test.jl

@@ -9,30 +9,31 @@ using Distributions
 using Random
 using UnitCommitment, Cbc, JuMP
 
-get_instance() = UnitCommitment.read_benchmark("matpower/case118/2017-02-01")
+function get_scenario()
-system_load(instance) = sum(b.load for b in instance.buses)
+    return UnitCommitment.read_benchmark(
+        "matpower/case118/2017-02-01",
+    ).scenarios[1]
+end
+system_load(sc) = sum(b.load for b in sc.buses)
 test_approx(x, y) = @test isapprox(x, y, atol = 1e-3)
 
 @testset "XavQiuAhm2021" begin
     @testset "cost and load share" begin
-        instance = get_instance()
+        sc = get_scenario()
-
         # Check original costs
-        unit = instance.units[10]
+        unit = sc.units[10]
         test_approx(unit.min_power_cost[1], 825.023)
         test_approx(unit.cost_segments[1].cost[1], 36.659)
         test_approx(unit.startup_categories[1].cost[1], 7570.42)
 
         # Check original load share
-        bus = instance.buses[1]
+        bus = sc.buses[1]
-        prev_system_load = system_load(instance)
+        prev_system_load = system_load(sc)
         test_approx(bus.load[1] / prev_system_load[1], 0.012)
 
         randomize!(
-            instance,
+            sc,
-            method = XavQiuAhm2021.Randomization(
+            XavQiuAhm2021.Randomization(randomize_load_profile = false),
-                randomize_load_profile = false,
-            ),
             rng = MersenneTwister(42),
         )
 
@@ -42,7 +43,7 @@ test_approx(x, y) = @test isapprox(x, y, atol = 1e-3)
         test_approx(unit.startup_categories[1].cost[1], 7634.226)
 
         # Check randomized load share
-        curr_system_load = system_load(instance)
+        curr_system_load = system_load(sc)
         test_approx(bus.load[1] / curr_system_load[1], 0.013)
 
         # System load should not change
@@ -50,20 +51,15 @@ test_approx(x, y) = @test isapprox(x, y, atol = 1e-3)
     end
 
     @testset "load profile" begin
-        instance = get_instance()
+        sc = get_scenario()
-
         # Check original load profile
-        @test round.(system_load(instance), digits = 1)[1:8] ≈
+        @test round.(system_load(sc), digits = 1)[1:8] ≈
               [3059.5, 2983.2, 2937.5, 2953.9, 3073.1, 3356.4, 4068.5, 4018.8]
 
-        randomize!(
+        randomize!(sc, XavQiuAhm2021.Randomization(); rng = MersenneTwister(42))
-            instance,
-            XavQiuAhm2021.Randomization(),
-            rng = MersenneTwister(42),
-        )
 
         # Check randomized load profile
-        @test round.(system_load(instance), digits = 1)[1:8] ≈
+        @test round.(system_load(sc), digits = 1)[1:8] ≈
               [4854.7, 4849.2, 4732.7, 4848.2, 4948.4, 5231.1, 5874.8, 5934.8]
     end
 end

test/transform/slice_test.jl

@@ -7,12 +7,13 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
 @testset "slice" begin
     instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
     modified = UnitCommitment.slice(instance, 1:2)
+    sc = modified.scenarios[1]
 
     # Should update all time-dependent fields
     @test modified.time == 2
-    @test length(modified.power_balance_penalty) == 2
+    @test length(sc.power_balance_penalty) == 2
-    @test length(modified.reserves_by_name["r1"].amount) == 2
+    @test length(sc.reserves_by_name["r1"].amount) == 2
-    for u in modified.units
+    for u in sc.units
         @test length(u.max_power) == 2
         @test length(u.min_power) == 2
         @test length(u.must_run) == 2
@@ -22,18 +23,19 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
             @test length(s.cost) == 2
         end
     end
-    for b in modified.buses
+    for b in sc.buses
         @test length(b.load) == 2
     end
-    for l in modified.lines
+    for l in sc.lines
         @test length(l.normal_flow_limit) == 2
         @test length(l.emergency_flow_limit) == 2
         @test length(l.flow_limit_penalty) == 2
     end
-    for ps in modified.price_sensitive_loads
+    for ps in sc.price_sensitive_loads
         @test length(ps.demand) == 2
         @test length(ps.revenue) == 2
     end
+
     # Should be able to build model without errors
     optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
     model = UnitCommitment.build_model(

test/usage.jl

@@ -4,34 +4,59 @@
 
 using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON
 
-@testset "usage" begin
+function _set_flow_limits!(instance)
-    instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
+    for sc in instance.scenarios
-    for line in instance.lines, t in 1:4
+        sc.power_balance_penalty = [100_000 for _ in 1:instance.time]
-        line.normal_flow_limit[t] = 10.0
+        for line in sc.lines, t in 1:4
+            line.normal_flow_limit[t] = 10.0
+        end
+    end
+end
+
+@testset "usage" begin
+    @testset "deterministic" begin
+        instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
+        _set_flow_limits!(instance)
+        optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
+        model = UnitCommitment.build_model(
+            instance = instance,
+            optimizer = optimizer,
+            variable_names = true,
+        )
+        @test name(model[:is_on]["g1", 1]) == "is_on[g1,1]"
+
+        # Optimize and retrieve solution
+        UnitCommitment.optimize!(model)
+        solution = UnitCommitment.solution(model)
+
+        # Write solution to a file
+        filename = tempname()
+        UnitCommitment.write(filename, solution)
+        loaded = JSON.parsefile(filename)
+        @test length(loaded["Is on"]) == 6
+
+        # Verify solution
+        @test UnitCommitment.validate(instance, solution)
+
+        # Reoptimize with fixed solution
+        UnitCommitment.fix!(model, solution)
+        UnitCommitment.optimize!(model)
+        @test UnitCommitment.validate(instance, solution)
+    end
+
+    @testset "stochastic" begin
+        instance = UnitCommitment.read([
+            "$FIXTURES/case14.json.gz",
+            "$FIXTURES/case14.json.gz",
+        ])
+        _set_flow_limits!(instance)
+        @test length(instance.scenarios) == 2
+        optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
+        model = UnitCommitment.build_model(
+            instance = instance,
+            optimizer = optimizer,
+        )
+        UnitCommitment.optimize!(model)
+        solution = UnitCommitment.solution(model)
     end
-    optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
-    model = UnitCommitment.build_model(
-        instance = instance,
-        optimizer = optimizer,
-        variable_names = true,
-    )
-    @test name(model[:is_on]["g1", 1]) == "is_on[g1,1]"
-
-    # Optimize and retrieve solution
-    UnitCommitment.optimize!(model)
-    solution = UnitCommitment.solution(model)
-
-    # Write solution to a file
-    filename = tempname()
-    UnitCommitment.write(filename, solution)
-    loaded = JSON.parsefile(filename)
-    @test length(loaded["Is on"]) == 6
-
-    # Verify solution
-    @test UnitCommitment.validate(instance, solution)
-
-    # Reoptimize with fixed solution
-    UnitCommitment.fix!(model, solution)
-    UnitCommitment.optimize!(model)
-    @test UnitCommitment.validate(instance, solution)
 end

test/validation/repair_test.jl

@@ -16,8 +16,8 @@ end
         json = parse_case14()
         json["Generators"]["g1"]["Production cost curve (MW)"] = [100, 150, 200]
         json["Generators"]["g1"]["Production cost curve (\$)"] = [10, 25, 30]
-        instance = UnitCommitment._from_json(json, repair = false)
+        sc = UnitCommitment._from_json(json, repair = false)
-        @test UnitCommitment.repair!(instance) == 4
+        @test UnitCommitment.repair!(sc) == 4
     end
 
     @testset "Startup limit must be greater than Pmin" begin
@@ -25,15 +25,15 @@ end
         json["Generators"]["g1"]["Production cost curve (MW)"] = [100, 150]
         json["Generators"]["g1"]["Production cost curve (\$)"] = [100, 150]
         json["Generators"]["g1"]["Startup limit (MW)"] = 80
-        instance = UnitCommitment._from_json(json, repair = false)
+        sc = UnitCommitment._from_json(json, repair = false)
-        @test UnitCommitment.repair!(instance) == 1
+        @test UnitCommitment.repair!(sc) == 1
     end
 
     @testset "Startup costs and delays must be increasing" begin
         json = parse_case14()
         json["Generators"]["g1"]["Startup costs (\$)"] = [300, 200, 100]
         json["Generators"]["g1"]["Startup delays (h)"] = [8, 4, 2]
-        instance = UnitCommitment._from_json(json, repair = false)
+        sc = UnitCommitment._from_json(json, repair = false)
-        @test UnitCommitment.repair!(instance) == 4
+        @test UnitCommitment.repair!(sc) == 4
     end
 end