From 7a03f4bbb09874486a3e7fa1022a54960af9beb6 Mon Sep 17 00:00:00 2001
From: Aleksandr Kazachkov <akazachk@users.noreply.github.com>
Date: Fri, 23 Jul 2021 11:23:16 -0500
Subject: [PATCH] Add reserve shortfall penalty

---
 src/instance/read.jl                  |  5 +++++
 src/instance/structs.jl               |  2 ++
 src/model/formulations/base/bus.jl    |  4 ----
 src/model/formulations/base/system.jl | 28 ++++++++++++++++++++++-----
 src/model/formulations/base/unit.jl   |  7 +------
 5 files changed, 31 insertions(+), 15 deletions(-)

diff --git a/src/instance/read.jl b/src/instance/read.jl
index 06ebfad..84e76c2 100644
--- a/src/instance/read.jl
+++ b/src/instance/read.jl
@@ -98,6 +98,10 @@ function _from_json(json; repair = true)
         json["Parameters"]["Power balance penalty (\$/MW)"],
         default = [1000.0 for t in 1:T],
     )
+    shortfall_penalty = timeseries(
+        json["Parameters"]["Reserve shortfall penalty (\$/MW)"],
+        default = [0.0 for t in 1:T],
+    )
 
     # Read buses
     for (bus_name, dict) in json["Buses"]
@@ -264,6 +268,7 @@ function _from_json(json; repair = true)
     instance = UnitCommitmentInstance(
         T,
         power_balance_penalty,
+        shortfall_penalty,
         units,
         buses,
         lines,
diff --git a/src/instance/structs.jl b/src/instance/structs.jl
index d75fba9..bf7360c 100644
--- a/src/instance/structs.jl
+++ b/src/instance/structs.jl
@@ -72,6 +72,8 @@ end
 mutable struct UnitCommitmentInstance
     time::Int
     power_balance_penalty::Vector{Float64}
+    "Penalty for failing to meet reserve requirement."
+    shortfall_penalty::Vector{Float64}
     units::Vector{Unit}
     buses::Vector{Bus}
     lines::Vector{TransmissionLine}
diff --git a/src/model/formulations/base/bus.jl b/src/model/formulations/base/bus.jl
index 00cacba..d881fc6 100644
--- a/src/model/formulations/base/bus.jl
+++ b/src/model/formulations/base/bus.jl
@@ -4,15 +4,11 @@
 
 function _add_bus!(model::JuMP.Model, b::Bus)::Nothing
     net_injection = _init(model, :expr_net_injection)
-    reserve = _init(model, :expr_reserve)
     curtail = _init(model, :curtail)
     for t in 1:model[:instance].time
         # Fixed load
         net_injection[b.name, t] = AffExpr(-b.load[t])
 
-        # Reserves
-        reserve[b.name, t] = AffExpr()
-
         # Load curtailment
         curtail[b.name, t] =
             @variable(model, lower_bound = 0, upper_bound = b.load[t])
diff --git a/src/model/formulations/base/system.jl b/src/model/formulations/base/system.jl
index 496ea2b..88bc78a 100644
--- a/src/model/formulations/base/system.jl
+++ b/src/model/formulations/base/system.jl
@@ -29,13 +29,31 @@ end
 
 function _add_reserve_eqs!(model::JuMP.Model)::Nothing
     eq_min_reserve = _init(model, :eq_min_reserve)
-    for t in 1:model[:instance].time
+    instance = model[:instance]
+    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)
+        shortfall_penalty = instance.shortfall_penalty[t]
         eq_min_reserve[t] = @constraint(
             model,
-            sum(
-                model[:expr_reserve][b.name, t] for b in model[:instance].buses
-            ) >= model[:instance].reserves.spinning[t]
+            sum(model[:reserve][g.name, t] for g in instance.units) + (
+                shortfall_penalty > 1e-7 ? model[:reserve_shortfall][t] : 0.0
+            ) >= instance.reserves.spinning[t]
         )
-    end
+
+        # Account for shortfall contribution to objective
+        if shortfall_penalty > 1e-7
+            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
     return
 end
diff --git a/src/model/formulations/base/unit.jl b/src/model/formulations/base/unit.jl
index ad00d44..2666ced 100644
--- a/src/model/formulations/base/unit.jl
+++ b/src/model/formulations/base/unit.jl
@@ -44,6 +44,7 @@ _is_initially_on(g::Unit)::Float64 = (g.initial_status > 0 ? 1.0 : 0.0)
 
 function _add_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
     reserve = _init(model, :reserve)
+    reserve_shortfall = _init(model, :reserve_shortfall)
     for t in 1:model[:instance].time
         if g.provides_spinning_reserves[t]
             reserve[g.name, t] = @variable(model, lower_bound = 0)
@@ -210,11 +211,5 @@ function _add_net_injection_eqs!(model::JuMP.Model, g::Unit)::Nothing
             model[:is_on][g.name, t],
             g.min_power[t],
         )
-        # Add to reserves expression
-        add_to_expression!(
-            model[:expr_reserve][g.bus.name, t],
-            model[:reserve][g.name, t],
-            1.0,
-        )
     end
 end