Ran JuliaFormatter

.github/workflows/test.yml

@@ -9,8 +9,8 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        julia-version: ['1.3', '1.4', '1.5', '1.6']
+        julia-version: ['1.4', '1.5', '1.6']
-        julia-arch: [x64, x86]
+        julia-arch: [x64]
         os: [ubuntu-latest, windows-latest, macOS-latest]
         exclude:
           - os: macOS-latest

docs/format.md

@@ -28,13 +28,14 @@ Each section is described in detail below. For a complete example, see [case14](
 
 ### Parameters
 
-This section describes system-wide parameters, such as power balance penalties,  optimization parameters, such as the length of the planning horizon and the time.
+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.
 
 | Key                            | Description                                       | Default  | Time series?
 | :----------------------------- | :------------------------------------------------ | :------: | :------------:
 | `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. Negative value implies reserve constraints must always be satisfied. | `-1` | Y
 
 
 #### Example
@@ -42,7 +43,8 @@ This section describes system-wide parameters, such as power balance penalties,
 {
     "Parameters": {
         "Time horizon (h)": 4,
-        "Power balance penalty ($/MW)": 1000.0
+        "Power balance penalty ($/MW)": 1000.0,
+        "Reserve shortfall penalty ($/MW)": -1.0
     }
 }
 ```

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 = [-1.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,

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}

src/model/build.jl

@@ -21,6 +21,8 @@ Arguments
 - `optimizer`:
     the optimizer factory that should be attached to this model (e.g. Cbc.Optimizer).
     If not provided, no optimizer will be attached.
+- `formulation`:
+    the details of which constraints, variables, etc. to use
 - `variable_names`: 
     If true, set variable and constraint names. Important if the model is going
     to be exported to an MPS file. For large models, this can take significant

src/model/formulations/ArrCon2000/ramp.jl

@@ -2,6 +2,15 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
+"""
+    _add_ramp_eqs!(model, unit, formulation_prod_vars, formulation_ramping, formulation_status_vars)::Nothing
+
+Ensure constraints on ramping are met.
+Based on Arroyo and Conejo (2000).
+Eqns. (24), (25) in Knueven et al. (2020).
+
+Adds constraints identified by `ArrCon200.Ramping` to `model` using variables `Gar1962.ProdVars` and `is_on` from `Gar1962.StatusVars`.
+"""
 function _add_ramp_eqs!(
     model::JuMP.Model,
     g::Unit,
@@ -22,7 +31,6 @@ function _add_ramp_eqs!(
     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)
 
     # Gar1962.ProdVars
     prod_above = model[:prod_above]
@@ -32,6 +40,8 @@ function _add_ramp_eqs!(
     switch_off = model[:switch_off]
     switch_on = model[:switch_on]
 
+    is_initially_on = _is_initially_on(g) > 0
+
     for t in 1:model[:instance].time
         # Ramp up limit
         if t == 1
@@ -56,7 +66,7 @@ function _add_ramp_eqs!(
             min_prod_last_period =
                 g.min_power[t-1] * is_on[gn, t-1] + prod_above[gn, t-1]
 
-            # Equation (24) in Kneuven et al. (2020)
+            # Equation (24) in Knueven et al. (2020)
             eq_ramp_up[gn, t] = @constraint(
                 model,
                 max_prod_this_period - min_prod_last_period <=
@@ -87,7 +97,7 @@ function _add_ramp_eqs!(
             min_prod_this_period =
                 g.min_power[t] * is_on[gn, t] + prod_above[gn, t]
 
-            # Equation (25) in Kneuven et al. (2020)
+            # Equation (25) in Knueven et al. (2020)
             eq_ramp_down[gn, t] = @constraint(
                 model,
                 max_prod_last_period - min_prod_this_period <=

src/model/formulations/ArrCon2000/structs.jl

@@ -13,6 +13,12 @@ module ArrCon2000
 
 import ..RampingFormulation
 
+"""
+Constraints
+---
+* `eq_ramp_up`: Equation (24) in Knueven et al. (2020)
+* `eq_ramp_down`: Equation (25) in Knueven et al. (2020)
+"""
 struct Ramping <: RampingFormulation end
 
 end

src/model/formulations/CarArr2006/pwlcosts.jl

@@ -2,6 +2,12 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
+"""
+    _add_production_piecewise_linear_eqs!
+
+Ensure respect of production limits along each segment.
+Creates constraints `CarArr2006.PwlCosts` using variables `Gar1962.StatusVars`
+"""
 function _add_production_piecewise_linear_eqs!(
     model::JuMP.Model,
     g::Unit,
@@ -21,7 +27,7 @@ function _add_production_piecewise_linear_eqs!(
     for t in 1:model[:instance].time
         gn = g.name
         for k in 1:K
-            # Equation (45) in Kneuven et al. (2020)
+            # Equation (45) in Knueven et al. (2020)
             # NB: when reading instance, UnitCommitment.jl already calculates
             #     difference between max power for segments k and k-1 so the
             #     value of cost_segments[k].mw[t] is the max production *for
@@ -36,14 +42,14 @@ function _add_production_piecewise_linear_eqs!(
             set_upper_bound(segprod[gn, t, k], g.cost_segments[k].mw[t])
 
             # Definition of production
-            # Equation (43) in Kneuven et al. (2020)
+            # Equation (43) in Knueven et al. (2020)
             eq_prod_above_def[gn, t] = @constraint(
                 model,
                 prod_above[gn, t] == sum(segprod[gn, t, k] for k in 1:K)
             )
 
             # Objective function
-            # Equation (44) in Kneuven et al. (2020)
+            # Equation (44) in Knueven et al. (2020)
             add_to_expression!(
                 model[:obj],
                 segprod[gn, t, k],

src/model/formulations/CarArr2006/structs.jl

@@ -14,6 +14,16 @@ module CarArr2006
 
 import ..PiecewiseLinearCostsFormulation
 
+"""
+Based on Garver (1962) and Carrión and Arryo (2006),
+which replaces (42) in Knueven et al. (2020) with a weaker version missing the on/off variable.
+Equations (45), (43), (44) in Knueven et al. (2020).
+
+Constraints
+---
+* `eq_prod_above_def`: Equation (43) in Knueven et al. (2020)
+* `eq_segprod_limit`: Equation (45) in Knueven et al. (2020)
+"""
 struct PwlCosts <: PiecewiseLinearCostsFormulation end
 
 end

src/model/formulations/DamKucRajAta2016/ramp.jl

@@ -2,6 +2,26 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
+"""
+    _add_ramp_eqs!
+
+Ensure constraints on ramping are met.
+Based on Damcı-Kurt et al. (2016).
+Eqns. (35), (36) in Knueven et al. (2020).
+
+Variables
+---
+* :prod_above
+* :reserve
+* :is_on
+* :switch_on
+* :switch_off],
+
+Constraints
+---
+* :eq_str_ramp_up
+* :eq_str_ramp_down
+"""
 function _add_ramp_eqs!(
     model::JuMP.Model,
     g::Unit,
@@ -14,12 +34,14 @@ function _add_ramp_eqs!(
     RESERVES_WHEN_RAMP_UP = true
     RESERVES_WHEN_RAMP_DOWN = true
     RESERVES_WHEN_SHUT_DOWN = true
-    known_initial_conditions = true
+    is_initially_on = _is_initially_on(g)
-    is_initially_on = (g.initial_status > 0)
+
-    SU = g.startup_limit
+    # The following are the same for generator g across all time periods
-    SD = g.shutdown_limit
+    SU = g.startup_limit   # startup rate
-    RU = g.ramp_up_limit
+    SD = g.shutdown_limit  # shutdown rate
-    RD = g.ramp_down_limit
+    RU = g.ramp_up_limit   # ramp up rate
+    RD = g.ramp_down_limit # ramp down rate
+
     gn = g.name
     eq_str_ramp_down = _init(model, :eq_str_ramp_down)
     eq_str_ramp_up = _init(model, :eq_str_ramp_up)
@@ -56,7 +78,7 @@ function _add_ramp_eqs!(
         if t > 1 && time_invariant
             min_prod_last_period = prod_above[gn, t-1]
 
-            # Equation (35) in Kneuven et al. (2020)
+            # Equation (35) in Knueven et al. (2020)
             # Sparser version of (24)
             eq_str_ramp_up[gn, t] = @constraint(
                 model,
@@ -76,7 +98,7 @@ function _add_ramp_eqs!(
             # (instead of using the amount above minimum, as min prod for t < 1 is unknown)
             max_prod_this_period += g.min_power[t] * is_on[gn, t]
 
-            # Modified version of equation (35) in Kneuven et al. (2020)
+            # Modified version of equation (35) in Knueven et al. (2020)
             # Equivalent to (24)
             eq_str_ramp_up[gn, t] = @constraint(
                 model,
@@ -94,12 +116,12 @@ function _add_ramp_eqs!(
         on_last_period = 0.0
         if t > 1
             on_last_period = is_on[gn, t-1]
-        elseif (known_initial_conditions && g.initial_status > 0)
+        elseif is_initially_on
             on_last_period = 1.0
         end
 
         if t > 1 && time_invariant
-            # Equation (36) in Kneuven et al. (2020)
+            # Equation (36) in Knueven et al. (2020)
             eq_str_ramp_down[gn, t] = @constraint(
                 model,
                 max_prod_last_period - min_prod_this_period <=
@@ -110,7 +132,7 @@ function _add_ramp_eqs!(
             # Add back in min power
             min_prod_this_period += g.min_power[t] * is_on[gn, t]
 
-            # Modified version of equation (36) in Kneuven et al. (2020)
+            # Modified version of equation (36) in Knueven et al. (2020)
             # Equivalent to (25)
             eq_str_ramp_down[gn, t] = @constraint(
                 model,

src/model/formulations/Gar1962/prod.jl

@@ -2,6 +2,11 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
+"""
+    _add_production_vars!(model, unit, formulation_prod_vars)
+
+Adds symbols identified by `Gar1962.ProdVars` to `model`.
+"""
 function _add_production_vars!(
     model::JuMP.Model,
     g::Unit,
@@ -18,6 +23,23 @@ function _add_production_vars!(
     return
 end
 
+"""
+    _add_production_limit_eqs!(model, unit, formulation_prod_vars)
+
+Ensure production limit constraints are met.
+Based on Garver (1962) and Morales-España et al. (2013).
+Eqns. (18), part of (69) in Knueven et al. (2020).
+
+===
+Variables
+* :is_on
+* :prod_above
+* :reserve
+
+===
+Constraints
+* :eq_prod_limit
+"""
 function _add_production_limit_eqs!(
     model::JuMP.Model,
     g::Unit,
@@ -30,13 +52,13 @@ function _add_production_limit_eqs!(
     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
+        # Part of (69) of Knueven 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)
+        # Equation (18) in Knueven et al. (2020)
         #   as \bar{p}_g(t) \le \bar{P}_g u_g(t)
-        # amk: this is a weaker version of (20) and (21) in Kneuven et al. (2020)
+        # amk: this is a weaker version of (20) and (21) in Knueven et al. (2020)
         #      but keeping it here in case those are not present
         power_diff = max(g.max_power[t], 0.0) - max(g.min_power[t], 0.0)
         if power_diff < 1e-7

src/model/formulations/Gar1962/pwlcosts.jl

@@ -2,6 +2,28 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
+"""
+    _add_production_piecewise_linear_eqs!
+
+Ensure respect of production limits along each segment.
+Based on Garver (1962).
+Equations (42), (43), (44) in Knueven et al. (2020).
+NB: when reading instance, UnitCommitment.jl already calculates 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*.
+
+Added to `model`: `:eq_prod_above_def` and `:eq_segprod_limit`.
+
+===
+Variables
+* :segprod
+* :is_on
+* :prod_above
+
+===
+Constraints
+* :eq_prod_above_def
+* :eq_segprod_limit
+"""
 function _add_production_piecewise_linear_eqs!(
     model::JuMP.Model,
     g::Unit,
@@ -23,14 +45,14 @@ function _add_production_piecewise_linear_eqs!(
     K = length(g.cost_segments)
     for t in 1:model[:instance].time
         # Definition of production
-        # Equation (43) in Kneuven et al. (2020)
+        # Equation (43) in Knueven et al. (2020)
         eq_prod_above_def[gn, t] = @constraint(
             model,
             prod_above[gn, t] == sum(segprod[gn, t, k] for k in 1:K)
         )
 
         for k in 1:K
-            # Equation (42) in Kneuven et al. (2020)
+            # Equation (42) in Knueven et al. (2020)
             # Without this, solvers will add a lot of implied bound cuts to
             # have this same effect.
             # NB: when reading instance, UnitCommitment.jl already calculates
@@ -47,7 +69,7 @@ function _add_production_piecewise_linear_eqs!(
             set_upper_bound(segprod[gn, t, k], g.cost_segments[k].mw[t])
 
             # Objective function
-            # Equation (44) in Kneuven et al. (2020)
+            # Equation (44) in Knueven et al. (2020)
             add_to_expression!(
                 model[:obj],
                 segprod[gn, t, k],

src/model/formulations/Gar1962/status.jl

@@ -2,6 +2,12 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
+"""
+    _add_status_vars!
+
+Adds symbols identified by `Gar1962.StatusVars` to `model`.
+Fix variables if a certain generator _must_ run or based on initial conditions.
+"""
 function _add_status_vars!(
     model::JuMP.Model,
     g::Unit,
@@ -10,20 +16,108 @@ function _add_status_vars!(
     is_on = _init(model, :is_on)
     switch_on = _init(model, :switch_on)
     switch_off = _init(model, :switch_off)
+    FIX_VARS = !formulation_status_vars.fix_vars_via_constraint
+    is_initially_on = _is_initially_on(g) > 0
     for t in 1:model[:instance].time
-        if g.must_run[t]
-            is_on[g.name, t] = 1.0
-            switch_on[g.name, t] = (t == 1 ? 1.0 - _is_initially_on(g) : 0.0)
-            switch_off[g.name, t] = 0.0
-        else
         is_on[g.name, t] = @variable(model, binary = true)
         switch_on[g.name, t] = @variable(model, binary = true)
         switch_off[g.name, t] = @variable(model, binary = true)
+
+        # Use initial conditions and whether a unit must run to fix variables
+        if FIX_VARS
+            # Fix variables using fix function
+            if g.must_run[t]
+                # If the generator _must_ run, then it is obviously on and cannot be switched off
+                # In the first time period, force unit to switch on if was off before
+                # Otherwise, unit is on, and will never turn off, so will never need to turn on
+                fix(is_on[g.name, t], 1.0; force = true)
+                fix(
+                    switch_on[g.name, t],
+                    (t == 1 ? 1.0 - _is_initially_on(g) : 0.0);
+                    force = true,
+                )
+                fix(switch_off[g.name, t], 0.0; force = true)
+            elseif t == 1
+                if is_initially_on
+                    # Generator was on (for g.initial_status time periods),
+                    # so cannot be more switched on until the period after the first time it can be turned off
+                    fix(switch_on[g.name, 1], 0.0; force = true)
+                else
+                    # Generator is initially off (for -g.initial_status time periods)
+                    # Cannot be switched off more
+                    fix(switch_off[g.name, 1], 0.0; force = true)
+                end
+            end
+        else
+            # Add explicit constraint if !FIX_VARS
+            if g.must_run[t]
+                is_on[g.name, t] = 1.0
+                switch_on[g.name, t] =
+                    (t == 1 ? 1.0 - _is_initially_on(g) : 0.0)
+                switch_off[g.name, t] = 0.0
+            elseif t == 1
+                if is_initially_on
+                    switch_on[g.name, t] = 0.0
+                else
+                    switch_off[g.name, t] = 0.0
+                end
+            end
+        end
+
+        # Use initial conditions and whether a unit must run to fix variables
+        if FIX_VARS
+            # Fix variables using fix function
+            if g.must_run[t]
+                # If the generator _must_ run, then it is obviously on and cannot be switched off
+                # In the first time period, force unit to switch on if was off before
+                # Otherwise, unit is on, and will never turn off, so will never need to turn on
+                fix(is_on[g.name, t], 1.0; force = true)
+                fix(
+                    switch_on[g.name, t],
+                    (t == 1 ? 1.0 - _is_initially_on(g) : 0.0);
+                    force = true,
+                )
+                fix(switch_off[g.name, t], 0.0; force = true)
+            elseif t == 1
+                if is_initially_on
+                    # Generator was on (for g.initial_status time periods),
+                    # so cannot be more switched on until the period after the first time it can be turned off
+                    fix(switch_on[g.name, 1], 0.0; force = true)
+                else
+                    # Generator is initially off (for -g.initial_status time periods)
+                    # Cannot be switched off more
+                    fix(switch_off[g.name, 1], 0.0; force = true)
+                end
+            end
+        else
+            # Add explicit constraint if !FIX_VARS
+            if g.must_run[t]
+                is_on[g.name, t] = 1.0
+                switch_on[g.name, t] =
+                    (t == 1 ? 1.0 - _is_initially_on(g) : 0.0)
+                switch_off[g.name, t] = 0.0
+            elseif t == 1
+                if is_initially_on
+                    switch_on[g.name, t] = 0.0
+                else
+                    switch_off[g.name, t] = 0.0
+                end
+            end
         end
     end
     return
 end
 
+"""
+    _add_status_eqs!
+
+Creates constraints `eq_binary_link` and `eq_switch_on_off` using variables in `Gar1962.StatusVars`.
+
+Constraints
+---
+* `eq_binary_link`
+* `eq_switch_on_off`
+"""
 function _add_status_eqs!(
     model::JuMP.Model,
     g::Unit,
@@ -35,8 +129,12 @@ function _add_status_eqs!(
     switch_off = model[:switch_off]
     switch_on = model[:switch_on]
     for t in 1:model[:instance].time
-        if !g.must_run[t]
+        if g.must_run[t]
+            continue
+        end
+
         # Link binary variables
+        # Equation (2) in Knueven et al. (2020), originally from Garver (1962)
         if t == 1
             eq_binary_link[g.name, t] = @constraint(
                 model,
@@ -50,12 +148,13 @@ function _add_status_eqs!(
                 switch_on[g.name, t] - switch_off[g.name, t]
             )
         end
+
         # Cannot switch on and off at the same time
+        # amk: I am not sure this is in Knueven et al. (2020)
         eq_switch_on_off[g.name, t] = @constraint(
             model,
             switch_on[g.name, t] + switch_off[g.name, t] <= 1
         )
     end
-    end
     return
 end

src/model/formulations/Gar1962/structs.jl

@@ -17,8 +17,53 @@ import ..PiecewiseLinearCostsFormulation
 import ..ProductionVarsFormulation
 import ..StatusVarsFormulation
 
+"""
+Variables
+---
+* `prod_above`:
+        [gen, t];
+        *production above minimum required level*;
+        lb: 0, ub: Inf.
+        KnuOstWat2020: `p'_g(t)`
+* `segprod`:
+        [gen, segment, t];
+        *how much generator produces on cost segment in time t*;
+        lb: 0, ub: Inf.
+        KnuOstWat2020: `p_g^l(t)`
+"""
 struct ProdVars <: ProductionVarsFormulation end
+
 struct PwlCosts <: PiecewiseLinearCostsFormulation end
-struct StatusVars <: StatusVarsFormulation end
+
+"""
+Variables
+---
+* `is_on`:
+        [gen, t];
+        *is generator on at time t?*
+        lb: 0, ub: 1, binary.
+        KnuOstWat2020: `u_g(t)`
+* `switch_on`:
+        [gen, t];
+        *indicator that generator will be turned on at t*;
+        lb: 0, ub: 1, binary.
+        KnuOstWat2020: `v_g(t)`
+* `switch_off`: binary;
+        [gen, t];
+        *indicator that generator will be turned off at t*;
+        lb: 0, ub: 1, binary.
+        KnuOstWat2020: `w_g(t)`
+
+Arguments
+---
+* `fix_vars_via_constraint`:
+        indicator for whether to set vars to a constant using `fix` or by adding an explicit constraint
+        (particulary useful for debugging purposes).
+"""
+struct StatusVars <: StatusVarsFormulation
+    fix_vars_via_constraint::Bool
+
+    StatusVars() = new(false)
+end
 
 end

src/model/formulations/GenMorRam2017/startstop.jl

@@ -0,0 +1,100 @@
+# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
+# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+# Released under the modified BSD license. See COPYING.md for more details.
+
+"""
+    _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
+
+Startup and shutdown limits from Gentile et al. (2017).
+Eqns. (20), (23a), and (23b) in Knueven et al. (2020).
+
+Creates constraints `eq_startstop_limit`,  `eq_startup_limit`, and `eq_shutdown_limit`
+using variables `Gar1962.StatusVars`, `prod_above` from `Gar1962.ProdVars`, and `reserve`.
+
+Constraints
+---
+* `eq_startstop_limit`
+* `eq_startup_limit`
+* `eq_shutdown_limit`
+"""
+function _add_startup_shutdown_limit_eqs!(
+    model::JuMP.Model,
+    g::Unit,
+    formulation_status_vars::Gar1962.StatusVars,
+)::Nothing
+    # TODO: Move upper case constants to model[:instance]
+    RESERVES_WHEN_START_UP = true
+    RESERVES_WHEN_RAMP_UP = true
+    RESERVES_WHEN_RAMP_DOWN = true
+    RESERVES_WHEN_SHUT_DOWN = true
+
+    eq_startstop_limit = _init(model, :eq_startstop_limit)
+    eq_shutdown_limit = _init(model, :eq_shutdown_limit)
+    eq_startup_limit = _init(model, :eq_startup_limit)
+
+    is_on = model[:is_on]
+    prod_above = model[:prod_above]
+    reserve = model[:reserve]
+    switch_off = model[:switch_off]
+    switch_on = model[:switch_on]
+
+    T = model[:instance].time
+    gi = g.name
+
+    if g.initial_power > g.shutdown_limit
+        #eqs.shutdown_limit[gi, 0] = @constraint(mip, vars.switch_off[gi, 1] <= 0)
+        if formulation_status_vars.always_create_vars
+            fix(switch_off[gi, 1], 0.0; force = true)
+            @constraint(mip, vars.switch_off[gi, 1] <= 0)
+        else
+            switch_off[gi, 1] = 0.0
+        end
+    end
+
+    for t in 1:T
+        ## 2020-10-09 amk: added eqn (20) and check of g.min_uptime
+        # Not present in (23) in Kneueven et al.
+        if g.min_uptime > 1
+            # Equation (20) in Knueven et al. (2020)
+            eqs.startstop_limit[gi, t] = @constraint(
+                model,
+                prod_above[gi, t] + reserve[gi, t] <=
+                (g.max_power[t] - g.min_power[t]) * is_on[gi, t] -
+                max(0, g.max_power[t] - g.startup_limit) * switch_on[gi, t] - (
+                    t < T ?
+                    max(0, g.max_power[t] - g.shutdown_limit) *
+                    switch_off[gi, t+1] : 0.0
+                )
+            )
+        else
+            ## Startup limits
+            # Equation (23a) in Knueven et al. (2020)
+            eqs.startup_limit[gi, t] = @constraint(
+                model,
+                prod_above[gi, t] + reserve[gi, t] <=
+                (g.max_power[t] - g.min_power[t]) * is_on[gi, t] -
+                max(0, g.max_power[t] - g.startup_limit) * switch_on[gi, t] - (
+                    t < T ?
+                    max(0, g.startup_limit - g.shutdown_limit) *
+                    switch_off[gi, t+1] : 0.0
+                )
+            )
+
+            ## Shutdown limits
+            if t < T
+                # Equation (23b) in Knueven et al. (2020)
+                eqs.shutdown_limit[gi, t] = @constraint(
+                    model,
+                    prod_above[gi, t] + reserve[gi, t] <=
+                    (g.max_power[t] - g.min_power[t]) * xis_on[gi, t] - (
+                        t < T ?
+                        max(0, g.max_power[t] - g.shutdown_limit) *
+                        switch_off[gi, t+1] : 0.0
+                    ) -
+                    max(0, g.shutdown_limit - g.startup_limit) *
+                    switch_on[gi, t]
+                )
+            end
+        end
+    end
+end

src/model/formulations/KnuOstWat2018/pwlcosts.jl

@@ -2,6 +2,30 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
+"""
+    _add_production_piecewise_linear_eqs!
+
+Ensure respect of production limits along each segment.
+Based on Knueven et al. (2018b).
+Eqns. (43), (44), (46), (48) in Knueven et al. (2020).
+NB: when reading instance, UnitCommitment.jl already calculates 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*.
+
+===
+Variables
+* :segprod
+* :is_on
+* :switch_on
+* :switch_off
+* :prod_above
+
+===
+Constraints
+* :eq_prod_above_def
+* :eq_segprod_limit_a
+* :eq_segprod_limit_b
+* :eq_segprod_limit_c
+"""
 function _add_production_piecewise_linear_eqs!(
     model::JuMP.Model,
     g::Unit,
@@ -57,7 +81,7 @@ function _add_production_piecewise_linear_eqs!(
             end
 
             if g.min_uptime > 1
-                # Equation (46) in Kneuven et al. (2020)
+                # Equation (46) in Knueven et al. (2020)
                 eq_segprod_limit_a[gn, t, k] = @constraint(
                     model,
                     segprod[gn, t, k] <=
@@ -66,7 +90,7 @@ function _add_production_piecewise_linear_eqs!(
                     (t < T ? Cw * switch_off[gn, t+1] : 0.0)
                 )
             else
-                # Equation (47a)/(48a) in Kneuven et al. (2020)
+                # Equation (47a)/(48a) in Knueven et al. (2020)
                 eq_segprod_limit_b[gn, t, k] = @constraint(
                     model,
                     segprod[gn, t, k] <=
@@ -75,7 +99,7 @@ function _add_production_piecewise_linear_eqs!(
                     (t < T ? max(0, Cv - Cw) * switch_off[gn, t+1] : 0.0)
                 )
 
-                # Equation (47b)/(48b) in Kneuven et al. (2020)
+                # Equation (47b)/(48b) in Knueven et al. (2020)
                 eq_segprod_limit_c[gn, t, k] = @constraint(
                     model,
                     segprod[gn, t, k] <=
@@ -86,14 +110,14 @@ function _add_production_piecewise_linear_eqs!(
             end
 
             # Definition of production
-            # Equation (43) in Kneuven et al. (2020)
+            # Equation (43) in Knueven et al. (2020)
             eq_prod_above_def[gn, t] = @constraint(
                 model,
                 prod_above[gn, t] == sum(segprod[gn, t, k] for k in 1:K)
             )
 
             # Objective function
-            # Equation (44) in Kneuven et al. (2020)
+            # Equation (44) in Knueven et al. (2020)
             add_to_expression!(
                 model[:obj],
                 segprod[gn, t, k],

src/model/formulations/KnuOstWat2018/scosts.jl

@@ -0,0 +1,123 @@
+# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
+# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+# Released under the modified BSD license. See COPYING.md for more details.
+
+"""
+    _add_startup_cost_eqs!
+
+Extended formulation of startup costs using indicator variables
+based on Knueven, Ostrowski, and Watson, 2020
+--- equations (59), (60), (61).
+
+Variables
+---
+* switch_on
+* switch_off
+* downtime_arc
+
+Constraints
+---
+* eq_startup_at_t
+* eq_shutdown_at_t
+"""
+function _add_startup_cost_eqs!(
+    model::JuMP.Model,
+    g::Unit,
+    formulation::MorLatRam2013.StartupCosts,
+)::Nothing
+    S = length(g.startup_categories)
+    if S == 0
+        return
+    end
+    gn = g.name
+
+    _init(model, eq_startup_at_t)
+    _init(model, eq_shutdown_at_t)
+
+    switch_on = model[:switch_on]
+    switch_off = model[:switch_off]
+    downtime_arc = model[:downtime_arc]
+
+    DT = g.min_downtime # minimum time offline
+    TC = g.startup_categories[S].delay # time offline until totally cold
+
+    # If initial_status < 0, then this is the amount of time the generator has been off
+    initial_time_shutdown = (g.initial_status < 0 ? -g.initial_status : 0)
+
+    for t in 1:model[:instance].time
+        # Fix to zero values of downtime_arc outside the feasible time pairs
+        # Specifically, x(t,t') = 0 if t' does not belong to 𝒢 = [t+DT, t+TC-1]
+        # This is because DT is the minimum downtime, so there is no way x(t,t')=1 for t'<t+DT
+        # and TC is the "time until cold" => if the generator starts afterwards, always has max cost
+        #start_time = min(t + DT, T)
+        #end_time = min(t + TC - 1, T)
+        #for tmp_t in t+1:start_time
+        #  fix(vars.downtime_arc[gn, t, tmp_t], 0.; force = true)
+        #end
+        #for tmp_t in end_time+1:T
+        #  fix(vars.downtime_arc[gn, t, tmp_t], 0.; force = true)
+        #end
+
+        # Equation (59) in Knueven et al. (2020)
+        # Relate downtime_arc with switch_on
+        # "switch_on[g,t] >= x_g(t',t) for all t' \in [t-TC+1, t-DT]"
+        eq_startup_at_t[gn, t] = @constraint(
+            model,
+            switch_on[gn, t] >= sum(
+                downtime_arc[gn, tmp_t, t] for
+                tmp_t in t-TC+1:t-DT if tmp_t >= 1
+            )
+        )
+
+        # Equation (60) in Knueven et al. (2020)
+        # "switch_off[g,t] >= x_g(t,t') for all t' \in [t+DT, t+TC-1]"
+        eqs.shutdown_at_t[gn, t] = @constraint(
+            model,
+            switch_off[gn, t] >= sum(
+                downtime_arc[gn, t, tmp_t] for
+                tmp_t in t+DT:t+TC-1 if tmp_t <= T
+            )
+        )
+
+        # Objective function terms for start-up costs
+        # Equation (61) in Knueven et al. (2020)
+        default_category = S
+        if initial_time_shutdown > 0 && t + initial_time_shutdown - 1 < TC
+            for s in 1:S-1
+                # If off for x periods before, then belongs to category s
+                # if -x+1 in [t-delay[s+1]+1,t-delay[s]]
+                # or, equivalently, if total time off in [delay[s], delay[s+1]-1]
+                # where total time off = t - 1 + initial_time_shutdown
+                # (the -1 because not off for current time period)
+                if t + initial_time_shutdown - 1 <
+                   g.startup_categories[s+1].delay
+                    default_category = s
+                    break # does not go into next category
+                end
+            end
+        end
+        add_to_expression!(
+            model[:obj],
+            switch_on[gn, t],
+            g.startup_categories[default_category].cost,
+        )
+
+        for s in 1:S-1
+            # Objective function terms for start-up costs
+            # Equation (61) in Knueven et al. (2020)
+            # Says to replace the cost of last category with cost of category s
+            start_range = max((t - g.startup_categories[s+1].delay + 1), 1)
+            end_range = min((t - g.startup_categories[s].delay), T - 1)
+            for tmp_t in start_range:end_range
+                if (t < tmp_t + DT) || (t >= tmp_t + TC) # the second clause should never be true for s < S
+                    continue
+                end
+                add_to_expression!(
+                    model[:obj],
+                    downtime_arc[gn, tmp_t, t],
+                    g.startup_categories[s].cost - g.startup_categories[S].cost,
+                )
+            end
+        end # iterate over startup categories
+    end # iterate over time
+end # add_startup_costs_KneOstWat20

src/model/formulations/MorLatRam2013/ramp.jl

@@ -2,6 +2,27 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
+"""
+    _add_ramp_eqs!
+
+Ensure constraints on ramping are met.
+Needs to be used in combination with shutdown rate constraints, e.g., (21b) in Knueven et al. (2020).
+Based on Morales-España, Latorre, and Ramos, 2013.
+Eqns. (26)+(27) [replaced by (24)+(25) if time-varying min demand] in Knueven et al. (2020).
+
+Variables
+---
+* :is_on
+* :switch_off
+* :switch_on
+* :prod_above
+* :reserve
+
+Constraints
+---
+* :eq_ramp_up
+* :eq_ramp_down
+"""
 function _add_ramp_eqs!(
     model::JuMP.Model,
     g::Unit,
@@ -15,10 +36,10 @@ function _add_ramp_eqs!(
     RESERVES_WHEN_RAMP_DOWN = true
     RESERVES_WHEN_SHUT_DOWN = true
     is_initially_on = (g.initial_status > 0)
-    SU = g.startup_limit
+    SU = g.startup_limit   # startup rate
-    SD = g.shutdown_limit
+    SD = g.shutdown_limit  # shutdown rate
-    RU = g.ramp_up_limit
+    RU = g.ramp_up_limit   # ramp up rate
-    RD = g.ramp_down_limit
+    RD = g.ramp_down_limit # ramp down rate
     gn = g.name
     eq_ramp_down = _init(model, :eq_ramp_down)
     eq_ramp_up = _init(model, :eq_str_ramp_up)
@@ -71,7 +92,7 @@ function _add_ramp_eqs!(
                     RU * is_on[gn, t-1] + SU * switch_on[gn, t]
                 )
             else
-                # Equation (26) in Kneuven et al. (2020)
+                # Equation (26) in Knueven 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(
@@ -115,7 +136,7 @@ function _add_ramp_eqs!(
                     RD * is_on[gn, t] + SD * switch_off[gn, t]
                 )
             else
-                # Equation (27) in Kneuven et al. (2020)
+                # Equation (27) in Knueven 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(

src/model/formulations/MorLatRam2013/scosts.jl

@@ -2,21 +2,56 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
+"""
+    _add_startup_cost_eqs!
+
+Extended formulation of startup costs using indicator variables
+based on Muckstadt and Wilson, 1968;
+this version by Morales-España, Latorre, and Ramos, 2013.
+Eqns. (54), (55), and (56) in Knueven et al. (2020).
+Note that the last 'constraint' is actually setting the objective.
+
+\tstartup[gi,s,t] ≤ sum_{i=s.delay}^{(s+1).delay-1} switch_off[gi,t-i]
+\tswitch_on[gi,t] = sum_{s=1}^{length(startup_categories)} startup[gi,s,t]
+\tstartup_cost[gi,t] = sum_{s=1}^{length(startup_categories)} cost_segments[s].cost * startup[gi,s,t]
+
+Variables
+---
+* startup
+* switch_on
+* switch_off
+
+Constraints
+---
+* eq_startup_choose
+* eq_startup_restrict
+"""
 function _add_startup_cost_eqs!(
     model::JuMP.Model,
     g::Unit,
     formulation::MorLatRam2013.StartupCosts,
 )::Nothing
+    S = length(g.startup_categories)
+    if S == 0
+        return
+    end
+
+    # Constraints created
     eq_startup_choose = _init(model, :eq_startup_choose)
     eq_startup_restrict = _init(model, :eq_startup_restrict)
-    S = length(g.startup_categories)
+
+    # Variables needed
     startup = model[:startup]
+    switch_on = model[:switch_on]
+    switch_off = model[:switch_off]
+
+    gn = g.name
     for t in 1:model[:instance].time
         # If unit is switching on, we must choose a startup category
-        eq_startup_choose[g.name, t] = @constraint(
+        # Equation (55) in Knueven et al. (2020)
+        eq_startup_choose[gn, t] = @constraint(
             model,
-            model[:switch_on][g.name, t] ==
+            switch_on[gn, t] == sum(startup[gn, t, s] for s in 1:S)
-            sum(startup[g.name, t, s] for s in 1:S)
         )
 
         for s in 1:S
@@ -26,25 +61,28 @@ function _add_startup_cost_eqs!(
                 range_start = t - g.startup_categories[s+1].delay + 1
                 range_end = t - g.startup_categories[s].delay
                 range = (range_start:range_end)
+                # If initial_status < 0, then this is the amount of time the generator has been off
                 initial_sum = (
                     g.initial_status < 0 && (g.initial_status + 1 in range) ? 1.0 : 0.0
                 )
-                eq_startup_restrict[g.name, t, s] = @constraint(
+                # Change of index version of equation (54) in Knueven et al. (2020):
+                #   startup[gi,s,t] ≤ sum_{i=s.delay}^{(s+1).delay-1} switch_off[gi,t-i]
+                eq_startup_restrict[gn, t, s] = @constraint(
                     model,
-                    startup[g.name, t, s] <=
+                    startup[gn, t, s] <=
-                    initial_sum + sum(
+                    initial_sum +
-                        model[:switch_off][g.name, i] for i in range if i >= 1
+                    sum(switch_off[gn, i] for i in range if i >= 1)
                 )
-                )
+            end # if s < S (not the last category)
-            end
 
             # Objective function terms for start-up costs
+            # Equation (56) in Knueven et al. (2020)
             add_to_expression!(
                 model[:obj],
-                startup[g.name, t, s],
+                startup[gn, t, s],
                 g.startup_categories[s].cost,
             )
-        end
+        end # iterate over startup categories
-    end
+    end # iterate over time
     return
 end

src/model/formulations/MorLatRam2013/startstop.jl

@@ -0,0 +1,89 @@
+# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
+# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+# Released under the modified BSD license. See COPYING.md for more details.
+
+"""
+    _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
+
+Startup and shutdown limits from Morales-España et al. (2013a).
+Eqns. (20), (21a), and (21b) in Knueven et al. (2020).
+
+Variables
+---
+* :is_on
+* :prod_above
+* :reserve
+* :switch_on
+* :switch_off
+
+Constraints
+---
+* :eq_startstop_limit
+* :eq_startup_limit
+* :eq_shutdown_limit
+"""
+function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
+    # TODO: Move upper case constants to model[:instance]
+    RESERVES_WHEN_START_UP = true
+    RESERVES_WHEN_RAMP_UP = true
+    RESERVES_WHEN_RAMP_DOWN = true
+    RESERVES_WHEN_SHUT_DOWN = true
+
+    eq_startstop_limit = _init(model, :eq_startstop_limit)
+    eq_shutdown_limit = _init(model, :eq_shutdown_limit)
+    eq_startup_limit = _init(model, :eq_startup_limit)
+
+    is_on = model[:is_on]
+    prod_above = model[:prod_above]
+    reserve = model[:reserve]
+    switch_off = model[:switch_off]
+    switch_on = model[:switch_on]
+
+    T = model[:instance].time
+    gi = g.name
+    for t in 1:T
+        ## 2020-10-09 amk: added eqn (20) and check of g.min_uptime
+        if g.min_uptime > 1 && t < T
+            # Equation (20) in Knueven et al. (2020)
+            # UT > 1 required, to guarantee that vars.switch_on[gi, t] and vars.switch_off[gi, t+1] are not both = 1 at the same time
+            eq_startstop_limit[gi, t] = @constraint(
+                model,
+                prod_above[gi, t] + reserve[gi, t] <=
+                (g.max_power[t] - g.min_power[t]) * is_on[gi, t] -
+                max(0, g.max_power[t] - g.startup_limit) * switch_on[gi, t] -
+                max(0, g.max_power[t] - g.shutdown_limit) * switch_off[gi, t+1]
+            )
+        else
+            ## Startup limits
+            # Equation (21a) in Knueven et al. (2020)
+            # Proposed by Morales-España et al. (2013a)
+            eqs_startup_limit[gi, t] = @constraint(
+                model,
+                prod_above[gi, t] + reserve[gi, t] <=
+                (g.max_power[t] - g.min_power[t]) * is_on[gi, t] -
+                max(0, g.max_power[t] - g.startup_limit) * switch_on[gi, t]
+            )
+
+            ## Shutdown limits
+            if t < T
+                # Equation (21b) in Knueven et al. (2020)
+                # TODO different from what was in previous model, due to reserve variable
+                # ax: ideally should have reserve_up and reserve_down variables
+                #     i.e., the generator should be able to increase/decrease production as specified
+                #     (this is a heuristic for a "robust" solution,
+                #     in case there is an outage or a surge, and flow has to be redirected)
+                # amk: if shutdown_limit is the max prod of generator in time period before shutting down,
+                #      then it makes sense to count reserves, because otherwise, if reserves ≠ 0,
+                #      then the generator will actually produce more than the limit
+                eqs.shutdown_limit[gi, t] = @constraint(
+                    model,
+                    prod_above[gi, t] +
+                    (RESERVES_WHEN_SHUT_DOWN ? reserve[gi, t] : 0.0) <=
+                    (g.max_power[t] - g.min_power[t]) * is_on[gi, t] -
+                    max(0, g.max_power[t] - g.shutdown_limit) *
+                    switch_off[gi, t+1]
+                )
+            end
+        end # check if g.min_uptime > 1
+    end # loop over time
+end # _add_startup_shutdown_limit_eqs!

src/model/formulations/NowRom2000/scosts.jl

@@ -0,0 +1,20 @@
+# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
+# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+# Released under the modified BSD license. See COPYING.md for more details.
+
+"""
+    _add_startup_cost_eqs!
+
+Based on Nowak and Römisch, 2000.
+Introduces auxiliary startup cost variable, c_g^SU(t) for each time period,
+and uses startup status variable, u_g(t);
+there are exponentially many facets in this space,
+but there is a linear-time separation algorithm (Brandenburg et al., 2017).
+"""
+function _add_startup_cost_eqs!(
+    model::JuMP.Model,
+    g::Unit,
+    formulation::MorLatRam2013.StartupCosts,
+)::Nothing
+    return error("Not implemented.")
+end

src/model/formulations/OstAnjVan2012/ramp.jl

@@ -0,0 +1,85 @@
+# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
+# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+# Released under the modified BSD license. See COPYING.md for more details.
+
+"""
+    _add_ramp_eqs!
+
+Ensure constraints on ramping are met.
+Based on Ostrowski, Anjos, Vannelli (2012).
+Eqn (37) in Knueven et al. (2020).
+
+Variables
+---
+* :is_on
+* :prod_above
+* :reserve
+
+Constraints
+---
+* :eq_str_prod_limit
+"""
+function _add_ramp_eqs!(
+    model::JuMP.Model,
+    g::Unit,
+    formulation_prod_vars::Gar1962.ProdVars,
+    formulation_ramping::MorLatRam2013.Ramping,
+    formulation_status_vars::Gar1962.StatusVars,
+)::Nothing
+    # TODO: Move upper case constants to model[:instance]
+    RESERVES_WHEN_START_UP = true
+    RESERVES_WHEN_RAMP_UP = true
+    RESERVES_WHEN_RAMP_DOWN = true
+    RESERVES_WHEN_SHUT_DOWN = true
+    is_initially_on = _is_initially_on(g)
+
+    gn = g.name
+    eq_str_prod_limit = _init(model, :eq_str_prod_limit)
+
+    # Variables that we need
+    reserve = model[:reserve]
+
+    # Gar1962.ProdVars
+    prod_above = model[:prod_above]
+
+    # Gar1962.StatusVars
+    is_on = model[:is_on]
+    switch_off = model[:switch_off]
+
+    # The following are the same for generator g across all time periods
+    UT = g.min_uptime
+
+    SU = g.startup_limit   # startup rate
+    SD = g.shutdown_limit  # shutdown rate
+    RU = g.ramp_up_limit   # ramp up rate
+    RD = g.ramp_down_limit # ramp down rate
+
+    # TODO check initial conditions, but maybe okay as long as (35) and (36) are also used
+    for t in 1:model[:instance].time
+        Pbar = g.max_power[t]
+
+        #TRD = floor((Pbar - SU)/RD)
+        # TODO check amk changed TRD wrt Knueven et al.
+        TRD = ceil((Pbar - SD) / RD) # ramp down time
+
+        if Pbar < 1e-7
+            # Skip this time period if max power = 0
+            continue
+        end
+
+        if UT >= 1
+            # Equation (37) in Knueven et al. (2020)
+            KSD = min(TRD, UT - 1, T - t - 1)
+            eq_str_prod_limit[gn, t] = @constraint(
+                model,
+                prod_above[gn, t] +
+                g.min_power[t] * is_on[gn, t] +
+                (RESERVES_WHEN_RAMP_DOWN ? reserve[gn, t] : 0.0) <=
+                Pbar * is_on[gi, t] - sum(
+                    (Pbar - (SD + i * RD)) * switch_off[gi, t+1+i] for
+                    i in 0:KSD
+                )
+            )
+        end # check UT >= 1
+    end # loop over time
+end

src/model/formulations/PanGua2016/ramp.jl

@@ -2,6 +2,28 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
+"""
+    _add_ramp_eqs!
+
+Add tighter upper bounds on production based on ramp-down trajectory.
+Based on (28) in Pan and Guan (2016).
+But there is an extra time period covered using (40) of Knueven et al. (2020).
+Eqns. (38), (40), (41) in Knueven et al. (2020).
+
+Variables
+---
+* :prod_above
+* :reserve
+* :is_on
+* :switch_on
+* :switch_off
+
+Constraints
+---
+* :str_prod_limit
+* :prod_limit_ramp_up_extra_period
+* :prod_limit_shutdown_trajectory
+"""
 function _add_ramp_eqs!(
     model::JuMP.Model,
     g::Unit,
@@ -41,14 +63,14 @@ function _add_ramp_eqs!(
         end
 
         #TRD = floor((Pbar - SU) / RD) # ramp down time
-        # TODO check amk changed TRD wrt Kneuven et al.
+        # TODO check amk changed TRD wrt Knueven et al.
         TRD = ceil((Pbar - SD) / RD)  # ramp down time
         TRU = floor((Pbar - SU) / RU) # ramp up time, can be negative if Pbar < SU
 
         # TODO check initial time periods: what if generator has been running for x periods?
         # But maybe ok as long as (35) and (36) are also used...
         if UT > 1
-            # Equation (38) in Kneuven et al. (2020)
+            # Equation (38) in Knueven et al. (2020)
             # Generalization of (20)
             # Necessary that if any of the switch_on = 1 in the sum,
             # then switch_off[gn, t+1] = 0
@@ -65,7 +87,7 @@ function _add_ramp_eqs!(
             )
 
             if UT - 2 < TRU
-                # Equation (40) in Kneuven et al. (2020)
+                # Equation (40) in Knueven 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(
                     model,
@@ -83,7 +105,7 @@ function _add_ramp_eqs!(
             KSD = min(TRD, UT - 1, T - t - 1)
             if KSD > 0
                 KSU = min(TRU, UT - 2 - KSD, t - 1)
-                # Equation (41) in Kneuven et al. (2020)
+                # Equation (41) in Knueven et al. (2020)
                 eq_prod_limit_shutdown_trajectory[gn, t] = @constraint(
                     model,
                     prod_above[gn, t] +

src/model/formulations/base/bus.jl

@@ -2,17 +2,18 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
+"""
+    _add_bus!(model::JuMP.Model, b::Bus)::Nothing
+
+Creates `expr_net_injection` and adds `curtail` variable to `model`.
+"""
 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])

src/model/formulations/base/psload.jl

@@ -2,6 +2,9 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
+"""
+    _add_price_sensitive_load!(model::JuMP.Model, ps::PriceSensitiveLoad)
+"""
 function _add_price_sensitive_load!(
     model::JuMP.Model,
     ps::PriceSensitiveLoad,

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/structs.jl

@@ -9,6 +9,15 @@ abstract type StartupCostsFormulation end
 abstract type StatusVarsFormulation end
 abstract type ProductionVarsFormulation end
 
+"""
+    struct Formulation
+
+Every formulation has to specify components for setting production variables and limits,
+startup costs and piecewise-linear costs, ramping variables and constraints,
+status variables (on/off, starting up, shutting down), and transmission constraints.
+
+Some of these components are allowed to be empty, as long as overall validity of the formulation is maintained.
+"""
 struct Formulation
     prod_vars::ProductionVarsFormulation
     pwl_costs::PiecewiseLinearCostsFormulation

src/model/formulations/base/system.jl

@@ -2,12 +2,32 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
+"""
+    _add_system_wide_eqs!(model::JuMP.Model)::Nothing
+
+Adds constraints that apply to the whole system, such as relating to net injection and reserves.
+"""
 function _add_system_wide_eqs!(model::JuMP.Model)::Nothing
     _add_net_injection_eqs!(model)
     _add_reserve_eqs!(model)
     return
 end
 
+"""
+    _add_net_injection_eqs!(model::JuMP.Model)::Nothing
+
+Adds `net_injection`, `eq_net_injection_def`, and `eq_power_balance` identifiers into `model`.
+
+Variables
+---
+* `expr_net_injection`
+* `net_injection`
+
+Constraints
+---
+* `eq_net_injection_def`
+* `eq_power_balance`
+"""
 function _add_net_injection_eqs!(model::JuMP.Model)::Nothing
     T = model[:instance].time
     net_injection = _init(model, :net_injection)
@@ -27,15 +47,49 @@ function _add_net_injection_eqs!(model::JuMP.Model)::Nothing
     return
 end
 
+"""
+    _add_reserve_eqs!(model::JuMP.Model)::Nothing
+
+Ensure constraints on reserves are met.
+Based on Morales-España et al. (2013a).
+Eqn. (68) from Knueven et al. (2020).
+
+Adds `eq_min_reserve` identifier to `model`, and corresponding constraint.
+
+Variables
+---
+* `reserve`
+* `reserve_shortfall`
+
+Constraints
+---
+* `eq_min_reserve`
+"""
 function _add_reserve_eqs!(model::JuMP.Model)::Nothing
+    instance = model[:instance]
     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 Knueven 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(
+            sum(model[:reserve][g.name, t] for g in instance.units) +
-                model[:expr_reserve][b.name, t] for b in model[:instance].buses
+            (shortfall_penalty >= 0 ? model[:reserve_shortfall][t] : 0.0) >=
-            ) >= model[:instance].reserves.spinning[t]
+            instance.reserves.spinning[t]
         )
+
+        # Account for shortfall contribution to objective
+        if shortfall_penalty >= 0
+            add_to_expression!(
+                model[:obj],
+                shortfall_penalty,
+                model[:reserve_shortfall][t],
+            )
+        end
     end
     return
 end

src/model/formulations/base/unit.jl

@@ -2,6 +2,15 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
+"""
+    _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation)
+
+Add production, reserve, startup, shutdown, and status variables,
+and constraints for min uptime/downtime, net injection, production, ramping, startup, shutdown, and status.
+
+Fix variables if a certain generator _must_ run or if a generator provides spinning reserves.
+Also, add overflow penalty to objective for each transmission line.
+"""
 function _add_unit!(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")
@@ -35,33 +44,62 @@ function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation)
         formulation.status_vars,
     )
     _add_startup_cost_eqs!(model, g, formulation.startup_costs)
-    _add_startup_shutdown_limit_eqs!(model, g)
+    _add_shutdown_cost_eqs!(model, g)
+    _add_startup_shutdown_limit_eqs!(
+        model,
+        g,
+        formulation.status_vars,
+        formulation.prod_vars,
+    )
     _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_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
+"""
+    _add_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
+
+Add `:reserve` variable to `model`, fixed to zero if no spinning reserves specified.
+"""
+function _add_reserve_vars!(
+    model::JuMP.Model,
+    g::Unit,
+    ALWAYS_CREATE_VARS = false,
+)::Nothing
     reserve = _init(model, :reserve)
+    reserve_shortfall = _init(model, :reserve_shortfall) # for accounting for shortfall penalty in the objective
     for t in 1:model[:instance].time
         if g.provides_spinning_reserves[t]
             reserve[g.name, t] = @variable(model, lower_bound = 0)
+        else
+            if ALWAYS_CREATE_VARS
+                reserve[g.name, t] = @variable(model, lower_bound = 0)
+                fix(reserve[g.name, t], 0.0; force = true)
             else
                 reserve[g.name, t] = 0.0
             end
         end
-    return
+        reserve_shortfall[t] =
-end
+            (model[:instance].shortfall_penalty[t] >= 0) ?
-
+            @variable(model, lower_bound = 0) : 0.0
-function _add_reserve_eqs!(model::JuMP.Model, g::Unit)::Nothing
-    reserve = model[:reserve]
-    for t in 1:model[:instance].time
-        add_to_expression!(expr_reserve[g.bus.name, t], reserve[g.name, t], 1.0)
     end
     return
 end
 
+"""
+    _add_reserve_eqs!(model::JuMP.Model, g::Unit)::Nothing
+"""
+function _add_reserve_eqs!(model::JuMP.Model, g::Unit)::Nothing
+    # nothing to do here
+    return
+end
+
+"""
+    _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing
+
+Add `startup` to model.
+"""
 function _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing
     startup = _init(model, :startup)
     for t in 1:model[:instance].time
@@ -72,14 +110,31 @@ 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
+"""
+    _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
+
+Creates startup/shutdown limit constraints below based on variables `Gar1962.StatusVars`, `prod_above` from `Gar1962.ProdVars`, and `reserve`.
+
+Constraints
+---
+* :eq_startup_limit
+* :eq_shutdown_limit
+"""
+function _add_startup_shutdown_limit_eqs!(
+    model::JuMP.Model,
+    g::Unit,
+    formulation_status_vars::Gar1962.StatusVars,
+    formulation_prod_vars::Gar1962.ProdVars,
+)::Nothing
     eq_shutdown_limit = _init(model, :eq_shutdown_limit)
     eq_startup_limit = _init(model, :eq_startup_limit)
+
     is_on = model[:is_on]
-    prod_above = model[:prod_above]
-    reserve = model[:reserve]
     switch_off = model[:switch_off]
     switch_on = model[:switch_on]
+    prod_above = model[:prod_above]
+    reserve = model[:reserve]
+
     T = model[:instance].time
     for t in 1:T
         # Startup limit
@@ -91,8 +146,15 @@ function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
         )
         # Shutdown limit
         if g.initial_power > g.shutdown_limit
+            # TODO check what happens with these variables when exporting the model
+            # Generator producing too much to be turned off in the first time period
+            # (can a binary variable have bounds x = 0?)
+            if formulation_status_vars.fix_vars_via_constraint
                 eq_shutdown_limit[g.name, 0] =
-                @constraint(model, switch_off[g.name, 1] <= 0)
+                    @constraint(model, model[:switch_off][g.name, 1] <= 0.0)
+            else
+                fix(model[:switch_off][g.name, 1], 0.0; force = true)
+            end
         end
         if t < T
             eq_shutdown_limit[g.name, t] = @constraint(
@@ -107,6 +169,32 @@ function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
     return
 end
 
+"""
+    _add_shutdown_cost_eqs!(model::JuMP.Model, g::Unit)::Nothing
+
+Variables
+---
+* `switch_off`
+"""
+function _add_shutdown_cost_eqs!(model::JuMP.Model, g::Unit)::Nothing
+    T = model[:instance].time
+    gi = g.name
+    for t in 1:T
+        shutdown_cost = 0.0
+        if shutdown_cost > 1e-7
+            # Equation (62) in Knueven et al. (2020)
+            add_to_expression!(
+                model[:obj],
+                model[:switch_off][gi, t],
+                shutdown_cost,
+            )
+        end
+    end # loop over time
+end
+
+"""
+    _add_ramp_eqs!(model, unit, formulation)
+"""
 function _add_ramp_eqs!(
     model::JuMP.Model,
     g::Unit,
@@ -153,6 +241,24 @@ function _add_ramp_eqs!(
     end
 end
 
+"""
+    _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing
+
+Ensure constraints on up/down time are met.
+Based on Garver (1962), Malkin (2003), and Rajan and Takritti (2005).
+Eqns. (3), (4), (5) in Knueven et al. (2020).
+
+Variables
+---
+* `is_on`
+* `switch_off`
+* `switch_on`
+
+Constraints
+---
+* `eq_min_uptime`
+* `eq_min_downtime`
+"""
 function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing
     is_on = model[:is_on]
     switch_off = model[:switch_off]
@@ -162,18 +268,24 @@ function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing
     T = model[:instance].time
     for t in 1:T
         # Minimum up-time
+        # Equation (4) in Knueven et al. (2020)
         eq_min_uptime[g.name, t] = @constraint(
             model,
             sum(switch_on[g.name, i] for i in (t-g.min_uptime+1):t if i >= 1) <= is_on[g.name, t]
         )
+
         # Minimum down-time
+        # Equation (5) in Knueven et al. (2020)
         eq_min_downtime[g.name, t] = @constraint(
             model,
             sum(
                 switch_off[g.name, i] for i in (t-g.min_downtime+1):t if i >= 1
             ) <= 1 - is_on[g.name, t]
         )
+
         # Minimum up/down-time for initial periods
+        # Equations (3a) and (3b) in Knueven et al. (2020)
+        # (using :switch_on and :switch_off instead of :is_on)
         if t == 1
             if g.initial_status > 0
                 eq_min_uptime[g.name, 0] = @constraint(
@@ -196,6 +308,9 @@ function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing
     end
 end
 
+"""
+    _add_net_injection_eqs!(model::JuMP.Model, g::Unit)::Nothing
+"""
 function _add_net_injection_eqs!(model::JuMP.Model, g::Unit)::Nothing
     expr_net_injection = model[:expr_net_injection]
     for t in 1:model[:instance].time
@@ -210,11 +325,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

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