Add reserve shortfall penalty

docs/format.md

@@ -35,7 +35,7 @@ This section describes system-wide parameters, such as power balance and reserve
 | `Time horizon (h)` | Length of the planning horizon (in hours). | Required | N
 | `Time step (min)` | Length of each time step (in minutes). Must be a divisor of 60 (e.g. 60, 30, 20, 15, etc). | `60` | N
 | `Power balance penalty ($/MW)` | Penalty for system-wide shortage or surplus in production (in $/MW). This is charged per time step. For example, if there is a shortage of 1 MW for three time steps, three times this amount will be charged. | `1000.0` | Y
-| `Reserve shortfall penalty (\$/MW)` | Penalty for system-wide shortage in meeting reserve requirements (in $/MW). This is charged per time step. | `0` | Y
+| `Reserve shortfall penalty ($/MW)` | Penalty for system-wide shortage in meeting reserve requirements (in $/MW). This is charged per time step. Negative value implies reserve constraints must always be satisfied. | `-1` | Y
 
 
 #### Example
@@ -43,8 +43,8 @@ This section describes system-wide parameters, such as power balance and reserve
 {
     "Parameters": {
         "Time horizon (h)": 4,
-        "Power balance penalty ($/MW)": 1000.0
+        "Power balance penalty ($/MW)": 1000.0,
-        "Reserve shortfall penalty ($/MW)": 0.0
+        "Reserve shortfall penalty ($/MW)": -1.0
     }
 }
 ```

src/instance/read.jl

@@ -100,7 +100,7 @@ function _from_json(json; repair = true)
     )
     shortfall_penalty = timeseries(
         json["Parameters"]["Reserve shortfall penalty (\$/MW)"],
-        default = [0.0 for t in 1:T],
+        default = [-1.0 for t in 1:T],
     )
 
     # Read buses

src/model/formulations/base/system.jl

@@ -68,30 +68,28 @@ Constraints
 function _add_reserve_eqs!(model::JuMP.Model)::Nothing
     instance = model[:instance]
     eq_min_reserve = _init(model, :eq_min_reserve)
+    instance = model[:instance]
     for t in 1:instance.time
-        # Equation (68) in Knueven et al. (2020)
+        # 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)
         shortfall_penalty = instance.shortfall_penalty[t]
         eq_min_reserve[t] = @constraint(
             model,
-            sum(model[:reserve][g.name, t] for g in instance.units) + (
+            sum(model[:reserve][g.name, t] for g in instance.units) +
-                shortfall_penalty > 1e-7 ? model[:reserve_shortfall][t] : 0.0
+            (shortfall_penalty >= 0 ? model[:reserve_shortfall][t] : 0.0) >=
-            ) >= instance.reserves.spinning[t]
+            instance.reserves.spinning[t]
         )
 
         # Account for shortfall contribution to objective
-        if shortfall_penalty > 1e-7
+        if shortfall_penalty >= 0
             add_to_expression!(
                 model[:obj],
                 shortfall_penalty,
                 model[:reserve_shortfall][t],
             )
-        else
-            # Not added to the model at all
-            #fix(model.vars.reserve_shortfall[t], 0.; force=true)
         end
-    end # loop over time
+    end
     return
 end

src/model/formulations/base/unit.jl

@@ -75,6 +75,9 @@ function _add_reserve_vars!(
                 reserve[g.name, t] = 0.0
             end
         end
+        reserve_shortfall[t] =
+            (model[:instance].shortfall_penalty[t] >= 0) ?
+            @variable(model, lower_bound = 0) : 0.0
     end
     return
 end

src/solution/solution.jl

@@ -51,6 +51,12 @@ function solution(model::JuMP.Model)::OrderedDict
     sol["Switch on"] = timeseries(model[:switch_on], instance.units)
     sol["Switch off"] = timeseries(model[:switch_off], instance.units)
     sol["Reserve (MW)"] = timeseries(model[:reserve], instance.units)
+    sol["Reserve shortfall (MW)"] = OrderedDict(
+        t =>
+            (instance.shortfall_penalty[t] >= 0) ?
+            round(value(model[:reserve_shortfall][t]), digits = 5) : 0.0 for
+        t in 1:instance.time
+    )
     sol["Net injection (MW)"] =
         timeseries(model[:net_injection], instance.buses)
     sol["Load curtail (MW)"] = timeseries(model[:curtail], instance.buses)

src/validation/validate.jl

@@ -324,11 +324,16 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01)
         # Verify spinning reserves
         reserve =
             sum(solution["Reserve (MW)"][g.name][t] for g in instance.units)
-        if reserve < instance.reserves.spinning[t] - tol
+        reserve_shortfall =
+            (instance.shortfall_penalty[t] >= 0) ?
+            solution["Reserve shortfall (MW)"][t] : 0
+
+        if reserve + reserve_shortfall < instance.reserves.spinning[t] - tol
             @error @sprintf(
-                "Insufficient spinning reserves at time %d (%.2f should be %.2f)",
+                "Insufficient spinning reserves at time %d (%.2f + %.2f should be %.2f)",
                 t,
                 reserve,
+                reserve_shortfall,
                 instance.reserves.spinning[t],
             )
             err_count += 1

test/model/formulations_test.jl

@@ -20,8 +20,14 @@ if ENABLE_LARGE_TESTS
 end
 
 function _small_test(formulation::Formulation)::Nothing
-    instance = UnitCommitment.read_benchmark("matpower/case118/2017-02-01")
+    instances = ["matpower/case118/2017-02-01", "test/case14"]
-    UnitCommitment.build_model(instance = instance, formulation = formulation)  # should not crash
+    for instance in instances
+        # Should not crash
+        UnitCommitment.build_model(
+            instance = UnitCommitment.read_benchmark(instance),
+            formulation = formulation,
+        )
+    end
     return
 end