Implement new reserves

docs/format.md

@@ -28,7 +28,7 @@ Each section is described in detail below. For a complete example, see [case14](
 
 ### Parameters
 
-This section describes system-wide parameters, such as power balance and reserve shortfall penalties, and optimization parameters, such as the length of the planning horizon and the time.
+This section describes system-wide parameters, such as power balance penalty, and optimization parameters, such as the length of the planning horizon and the time.
 
 | Key                            | Description                                       | Default  | Time series?
 | :----------------------------- | :------------------------------------------------ | :------: | :------------:

docs/model.md

@@ -23,7 +23,7 @@ Name | Symbol | Description | Unit
 `switch_off[g,t]` | $w_{g}(t)$ | True if generator `g` switches off at time `t`. | Binary
 `prod_above[g,t]` |$p'_{g}(t)$ | Amount of power produced by generator `g` above its minimum power output at time `t`. For example, if the minimum power of generator `g` is 100 MW and `g` is producing 115 MW of power at time `t`, then `prod_above[g,t]` equals `15.0`. | MW
 `segprod[g,t,k]` | $p^k_g(t)$ | Amount of power from piecewise linear segment `k` produced by generator `g` at time `t`. For example, if cost curve for generator `g` is defined by the points `(100, 1400)`, `(110, 1600)`, `(130, 2200)` and `(135, 2400)`, and if the generator is producing 115 MW of power at time `t`, then `segprod[g,t,:]` equals `[10.0, 5.0, 0.0]`.| MW
-`reserve[g,t]` | $r_g(t)$ | Amount of reserves provided by generator `g` at time `t`. | MW
+`reserve[r,g,t]` | $r_g(t)$ | Amount of reserve `r` provided by unit `g` at time `t`. | MW
 `startup[g,t,s]` | $\delta^s_g(t)$ | True if generator `g` switches on at time `t` incurring start-up costs from start-up category `s`. | Binary
 
 

src/instance/read.jl

@@ -125,6 +125,11 @@ function _from_json(json; repair = true)
                 name = reserve_name,
                 type = lowercase(dict["Type"]),
                 amount = timeseries(dict["Amount (MW)"]),
+                units = [],
+                shortfall_penalty = scalar(
+                    dict["Shortfall penalty (\$/MW)"],
+                    default = -1,
+                ),
             )
             name_to_reserve[reserve_name] = reserve
             push!(reserves2, reserve)
@@ -171,7 +176,8 @@ function _from_json(json; repair = true)
         # Read reserves
         unit_reserves = Reserve[]
         if "Reserve eligibility" in keys(dict)
-            unit_reserves = [name_to_reserve[n] for n in dict["Reserve eligibility"]]
+            unit_reserves =
+                [name_to_reserve[n] for n in dict["Reserve eligibility"]]
         end
 
         # Read and validate initial conditions
@@ -215,6 +221,9 @@ function _from_json(json; repair = true)
             unit_reserves,
         )
         push!(bus.units, unit)
+        for r in unit_reserves
+            push!(r.units, unit)
+        end
         name_to_unit[unit_name] = unit
         push!(units, unit)
     end

src/instance/structs.jl

@@ -24,6 +24,8 @@ Base.@kwdef mutable struct Reserve
     name::String
     type::String
     amount::Vector{Float64}
+    units::Vector
+    shortfall_penalty::Float64
 end
 
 mutable struct Unit

src/model/formulations/ArrCon2000/ramp.jl

@@ -19,10 +19,10 @@ function _add_ramp_eqs!(
     RD = g.ramp_down_limit
     SU = g.startup_limit
     SD = g.shutdown_limit
-    reserve = model[:reserve]
     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)
 
     # Gar1962.ProdVars
     prod_above = model[:prod_above]
@@ -41,7 +41,7 @@ function _add_ramp_eqs!(
                     model,
                     g.min_power[t] +
                     prod_above[gn, t] +
-                    (RESERVES_WHEN_RAMP_UP ? reserve[gn, t] : 0.0) <=
+                    (RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) <=
                     g.initial_power + RU
                 )
             end
@@ -51,7 +51,7 @@ function _add_ramp_eqs!(
                 prod_above[gn, t] +
                 (
                     RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ?
-                    reserve[gn, t] : 0.0
+                    reserve[t] : 0.0
                 )
             min_prod_last_period =
                 g.min_power[t-1] * is_on[gn, t-1] + prod_above[gn, t-1]
@@ -82,7 +82,7 @@ function _add_ramp_eqs!(
                 prod_above[gn, t-1] +
                 (
                     RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ?
-                    reserve[gn, t-1] : 0.0
+                    reserve[t-1] : 0.0
                 )
             min_prod_this_period =
                 g.min_power[t] * is_on[gn, t] + prod_above[gn, t]

src/model/formulations/DamKucRajAta2016/ramp.jl

@@ -23,7 +23,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 = model[:reserve]
+    reserve = _total_reserves(model, g)
 
     # Gar1962.ProdVars
     prod_above = model[:prod_above]
@@ -48,10 +48,8 @@ function _add_ramp_eqs!(
         # end
 
         max_prod_this_period =
-            prod_above[gn, t] + (
+            prod_above[gn, t] +
-                RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ?
+            (RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0)
-                reserve[gn, t] : 0.0
-            )
         min_prod_last_period = 0.0
         if t > 1 && time_invariant
             min_prod_last_period = prod_above[gn, t-1]
@@ -88,7 +86,7 @@ function _add_ramp_eqs!(
         max_prod_last_period =
             min_prod_last_period + (
                 t > 1 && (RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN) ?
-                reserve[gn, t-1] : 0.0
+                reserve[t-1] : 0.0
             )
         min_prod_this_period = prod_above[gn, t]
         on_last_period = 0.0

src/model/formulations/Gar1962/prod.jl

@@ -26,7 +26,7 @@ function _add_production_limit_eqs!(
     eq_prod_limit = _init(model, :eq_prod_limit)
     is_on = model[:is_on]
     prod_above = model[:prod_above]
-    reserve = model[:reserve]
+    reserve = _total_reserves(model, g)
     gn = g.name
     for t in 1:model[:instance].time
         # Objective function terms for production costs
@@ -44,7 +44,7 @@ function _add_production_limit_eqs!(
         end
         eq_prod_limit[gn, t] = @constraint(
             model,
-            prod_above[gn, t] + reserve[gn, t] <= power_diff * is_on[gn, t]
+            prod_above[gn, t] + reserve[t] <= power_diff * is_on[gn, t]
         )
     end
 end

src/model/formulations/MorLatRam2013/ramp.jl

@@ -22,7 +22,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 = model[:reserve]
+    reserve = _total_reserves(model, g)
 
     # Gar1962.ProdVars
     prod_above = model[:prod_above]
@@ -43,7 +43,7 @@ function _add_ramp_eqs!(
                     model,
                     g.min_power[t] +
                     prod_above[gn, t] +
-                    (RESERVES_WHEN_RAMP_UP ? reserve[gn, t] : 0.0) <=
+                    (RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) <=
                     g.initial_power + RU
                 )
             end
@@ -61,7 +61,7 @@ function _add_ramp_eqs!(
                     prod_above[gn, t] +
                     (
                         RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ?
-                        reserve[gn, t] : 0.0
+                        reserve[t] : 0.0
                     )
                 min_prod_last_period =
                     g.min_power[t-1] * is_on[gn, t-1] + prod_above[gn, t-1]
@@ -77,7 +77,7 @@ function _add_ramp_eqs!(
                 eq_ramp_up[gn, t] = @constraint(
                     model,
                     prod_above[gn, t] +
-                    (RESERVES_WHEN_RAMP_UP ? reserve[gn, t] : 0.0) -
+                    (RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) -
                     prod_above[gn, t-1] <= RU
                 )
             end
@@ -105,7 +105,7 @@ function _add_ramp_eqs!(
                     prod_above[gn, t-1] +
                     (
                         RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ?
-                        reserve[gn, t-1] : 0.0
+                        reserve[t-1] : 0.0
                     )
                 min_prod_this_period =
                     g.min_power[t] * is_on[gn, t] + prod_above[gn, t]
@@ -121,7 +121,7 @@ function _add_ramp_eqs!(
                 eq_ramp_down[gn, t] = @constraint(
                     model,
                     prod_above[gn, t-1] +
-                    (RESERVES_WHEN_RAMP_DOWN ? reserve[gn, t-1] : 0.0) -
+                    (RESERVES_WHEN_RAMP_DOWN ? reserve[t-1] : 0.0) -
                     prod_above[gn, t] <= RD
                 )
             end

src/model/formulations/PanGua2016/ramp.jl

@@ -12,7 +12,7 @@ function _add_ramp_eqs!(
     # TODO: Move upper case constants to model[:instance]
     RESERVES_WHEN_SHUT_DOWN = true
     gn = g.name
-    reserve = model[:reserve]
+    reserve = _total_reserves(model, g)
     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)
@@ -56,7 +56,7 @@ function _add_ramp_eqs!(
                 model,
                 prod_above[gn, t] +
                 g.min_power[t] * is_on[gn, t] +
-                reserve[gn, t] <=
+                reserve[t] <=
                 Pbar * is_on[gn, t] -
                 (t < T ? (Pbar - SD) * switch_off[gn, t+1] : 0.0) - sum(
                     (Pbar - (SU + i * RU)) * switch_on[gn, t-i] for
@@ -71,7 +71,7 @@ function _add_ramp_eqs!(
                     model,
                     prod_above[gn, t] +
                     g.min_power[t] * is_on[gn, t] +
-                    reserve[gn, t] <=
+                    reserve[t] <=
                     Pbar * is_on[gn, t] - sum(
                         (Pbar - (SU + i * RU)) * switch_on[gn, t-i] for
                         i in 0:min(UT - 1, TRU, t - 1)
@@ -88,7 +88,7 @@ function _add_ramp_eqs!(
                     model,
                     prod_above[gn, t] +
                     g.min_power[t] * is_on[gn, t] +
-                    (RESERVES_WHEN_SHUT_DOWN ? reserve[gn, t] : 0.0) <=
+                    (RESERVES_WHEN_SHUT_DOWN ? reserve[t] : 0.0) <=
                     Pbar * is_on[gn, t] - sum(
                         (Pbar - (SD + i * RD)) * switch_off[gn, t+1+i] for
                         i in 0:KSD

src/model/formulations/base/system.jl

@@ -52,5 +52,29 @@ function _add_reserve_eqs!(model::JuMP.Model)::Nothing
             )
         end
     end
+
+    eq_min_reserve2 = _init(model, :eq_min_reserve2)
+    for r in instance.reserves2
+        for t in 1:instance.time
+            # 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_reserve2[r.name, t] = @constraint(
+                model,
+                sum(model[:reserve2][r.name, g.name, t] for g in r.units) +
+                model[:reserve_shortfall2][r.name, t] >= r.amount[t]
+            )
+
+            # Account for shortfall contribution to objective
+            if r.shortfall_penalty >= 0
+                add_to_expression!(
+                    model[:obj],
+                    r.shortfall_penalty,
+                    model[:reserve_shortfall2][r.name, t],
+                )
+            end
+        end
+    end
     return
 end

src/model/formulations/base/unit.jl

@@ -55,13 +55,19 @@ function _add_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
             (model[:instance].shortfall_penalty[t] >= 0) ?
             @variable(model, lower_bound = 0) : 0.0
     end
-    return
-end
 
-function _add_reserve_eqs!(model::JuMP.Model, g::Unit)::Nothing
+    reserve2 = _init(model, :reserve2)
-    reserve = model[:reserve]
+    reserve_shortfall2 = _init(model, :reserve_shortfall2)
+    for r in g.reserves
         for t in 1:model[:instance].time
-        add_to_expression!(expr_reserve[g.bus.name, t], reserve[g.name, t], 1.0)
+            reserve2[r.name, g.name, t] = @variable(model, lower_bound = 0)
+            if (r.name, t) ∉ keys(reserve_shortfall2)
+                reserve_shortfall2[r.name, t] = @variable(model, lower_bound = 0)
+                if r.shortfall_penalty < 0
+                    set_upper_bound(reserve_shortfall2[r.name, t], 0.0)
+                end
+            end
+        end
     end
     return
 end
@@ -81,7 +87,7 @@ function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
     eq_startup_limit = _init(model, :eq_startup_limit)
     is_on = model[:is_on]
     prod_above = model[:prod_above]
-    reserve = model[:reserve]
+    reserve = _total_reserves(model, g)
     switch_off = model[:switch_off]
     switch_on = model[:switch_on]
     T = model[:instance].time
@@ -89,7 +95,7 @@ function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
         # Startup limit
         eq_startup_limit[g.name, t] = @constraint(
             model,
-            prod_above[g.name, t] + reserve[g.name, t] <=
+            prod_above[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]
         )
@@ -117,7 +123,7 @@ function _add_ramp_eqs!(
     formulation::RampingFormulation,
 )::Nothing
     prod_above = model[:prod_above]
-    reserve = model[:reserve]
+    reserve = _total_reserves(model, g)
     eq_ramp_up = _init(model, :eq_ramp_up)
     eq_ramp_down = _init(model, :eq_ramp_down)
     for t in 1:model[:instance].time
@@ -126,14 +132,14 @@ function _add_ramp_eqs!(
             if _is_initially_on(g) == 1
                 eq_ramp_up[g.name, t] = @constraint(
                     model,
-                    prod_above[g.name, t] + reserve[g.name, t] <=
+                    prod_above[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(
                 model,
-                prod_above[g.name, t] + reserve[g.name, t] <=
+                prod_above[g.name, t] + reserve[t] <=
                 prod_above[g.name, t-1] + g.ramp_up_limit
             )
         end
@@ -216,3 +222,15 @@ function _add_net_injection_eqs!(model::JuMP.Model, g::Unit)::Nothing
         )
     end
 end
+
+function _total_reserves(model, g)::Vector
+    T = model[:instance].time
+    reserve = [model[:reserve][g.name, t] for t in 1:T]
+    if !isempty(g.reserves)
+        reserve += [
+            sum(model[:reserve2][r.name, g.name, t] for r in g.reserves) for
+            t in 1:model[:instance].time
+        ]
+    end
+    return reserve
+end

src/solution/solution.jl

@@ -67,5 +67,19 @@ function solution(model::JuMP.Model)::OrderedDict
         sol["Price-sensitive loads (MW)"] =
             timeseries(model[:loads], instance.price_sensitive_loads)
     end
+    sol["Reserve 2 (MW)"] = OrderedDict(
+        r.name => OrderedDict(
+            g.name => [
+                value(model[:reserve2][r.name, g.name, t]) for
+                t in 1:instance.time
+            ] for g in r.units
+        ) for r in instance.reserves2
+    )
+    sol["Reserve shortfall 2 (MW)"] = OrderedDict(
+        r.name => [
+            value(model[:reserve_shortfall2][r.name, t]) for
+            t in 1:instance.time
+        ] for r in instance.reserves2
+    )
     return sol
 end

src/validation/validate.jl

@@ -46,6 +46,12 @@ function _validate_units(instance, solution; tol = 0.01)
     for unit in instance.units
         production = solution["Production (MW)"][unit.name]
         reserve = solution["Reserve (MW)"][unit.name]
+        if !isempty(unit.reserves)
+            reserve += sum(
+                solution["Reserve 2 (MW)"][r.name][unit.name] for
+                r in unit.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])
@@ -137,9 +143,11 @@ function _validate_units(instance, solution; tol = 0.01)
             # 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",
+                    "Unit %s produces power at time %d while off (%.2f + %.2f > 0)",
                     unit.name,
-                    t
+                    t,
+                    production[t],
+                    reserve[t],
                 )
                 err_count += 1
             end
@@ -338,6 +346,27 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01)
             )
             err_count += 1
         end
+
+        # Verify reserves
+        for r in instance.reserves2
+            provided = sum(
+                solution["Reserve 2 (MW)"][r.name][g.name][t] for g in r.units
+            )
+            shortfall = solution["Reserve shortfall 2 (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,
+                )
+                err_count += 1
+            end
+        end
     end
 
     return err_count

test/model/formulations_test.jl

@@ -4,6 +4,8 @@
 
 using UnitCommitment
 using JuMP
+using Cbc
+using JSON
 import UnitCommitment:
     ArrCon2000,
     CarArr2006,
@@ -19,22 +21,29 @@ if ENABLE_LARGE_TESTS
     using Gurobi
 end
 
-function _small_test(formulation::Formulation)::Nothing
+function _small_test(formulation::Formulation; dump::Bool = false)::Nothing
-    instances = ["matpower/case118/2017-02-01", "test/case14"]
+    instance = UnitCommitment.read_benchmark("test/case14")
-    for instance in instances
+    model = UnitCommitment.build_model(
-        # Should not crash
+        instance = instance,
-        UnitCommitment.build_model(
-            instance = UnitCommitment.read_benchmark(instance),
         formulation = formulation,
+        optimizer = Cbc.Optimizer,
+        variable_names = true,
     )
+    UnitCommitment.optimize!(model)
+    solution = UnitCommitment.solution(model)
+    if dump
+        open("/tmp/ucjl.json", "w") do f
+            return write(f, JSON.json(solution, 2))
+        end
+        write_to_file(model, "/tmp/ucjl.lp")
     end
+    @test UnitCommitment.validate(instance, solution)
     return
 end
 
 function _large_test(formulation::Formulation)::Nothing
-    instances = ["pglib-uc/ca/Scenario400_reserves_1"]
+    instance =
-    for instance in instances
+        UnitCommitment.read_benchmark("pglib-uc/ca/Scenario400_reserves_1")
-        instance = UnitCommitment.read_benchmark(instance)
     model = UnitCommitment.build_model(
         instance = instance,
         formulation = formulation,
@@ -46,29 +55,44 @@ function _large_test(formulation::Formulation)::Nothing
     )
     solution = UnitCommitment.solution(model)
     @test UnitCommitment.validate(instance, solution)
-    end
     return
 end
 
-function _test(formulation::Formulation)::Nothing
+function _test(formulation::Formulation; dump::Bool = false)::Nothing
-    _small_test(formulation)
+    _small_test(formulation; dump)
     if ENABLE_LARGE_TESTS
         _large_test(formulation)
     end
 end
 
 @testset "formulations" begin
+    @testset "default" begin
         _test(Formulation())
+    end
+    @testset "ArrCon2000" begin
         _test(Formulation(ramping = ArrCon2000.Ramping()))
-    # _test(Formulation(ramping = DamKucRajAta2016.Ramping()))
+    end
+    @testset "DamKucRajAta2016" begin
+        _test(Formulation(ramping = DamKucRajAta2016.Ramping()))
+    end
+    @testset "MorLatRam2013" begin
         _test(
             Formulation(
                 ramping = MorLatRam2013.Ramping(),
                 startup_costs = MorLatRam2013.StartupCosts(),
             ),
         )
+    end
+    @testset "PanGua2016" begin
         _test(Formulation(ramping = PanGua2016.Ramping()))
+    end
+    @testset "Gar1962" begin
         _test(Formulation(pwl_costs = Gar1962.PwlCosts()))
+    end
+    @testset "CarArr2006" begin
         _test(Formulation(pwl_costs = CarArr2006.PwlCosts()))
+    end
+    @testset "KnuOstWat2018" begin
         _test(Formulation(pwl_costs = KnuOstWat2018.PwlCosts()))
     end
+end

test/runtests.jl

@@ -21,11 +21,13 @@ const ENABLE_LARGE_TESTS = ("UCJL_LARGE_TESTS" in keys(ENV))
     @testset "model" begin
         include("model/formulations_test.jl")
     end
+    @testset "solution" begin
         @testset "XavQiuWanThi19" begin
             include("solution/methods/XavQiuWanThi19/filter_test.jl")
             include("solution/methods/XavQiuWanThi19/find_test.jl")
             include("solution/methods/XavQiuWanThi19/sensitivity_test.jl")
         end
+    end
     @testset "transform" begin
         include("transform/initcond_test.jl")
         include("transform/slice_test.jl")

test/usage.jl

@@ -4,7 +4,7 @@
 
 using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON
 
-@testset "build_model" begin
+@testset "usage" begin
     instance = UnitCommitment.read_benchmark("test/case14")
     for line in instance.lines, t in 1:4
         line.normal_flow_limit[t] = 10.0