Ran JuliaFormatter

Fix duplicated startup constraint

.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-arch: [x64, x86]
+        julia-version: ['1.4', '1.5', '1.6']
+        julia-arch: [x64]
         os: [ubuntu-latest, windows-latest, macOS-latest]
         exclude:
           - os: macOS-latest

CHANGELOG.md

@@ -11,6 +11,11 @@ All notable changes to this project will be documented in this file.
 [semver]: https://semver.org/spec/v2.0.0.html
 [pkjjl]: https://pkgdocs.julialang.org/v1/compatibility/#compat-pre-1.0
 
+## [0.2.2] - 2021-07-21
+### Fixed
+- Fix small bug in validation scripts related to startup costs
+- Fix duplicated startup constraints (@mtanneau, #12)
+
 ## [0.2.1] - 2021-06-02
 ### Added
 - Add multiple ramping formulations (ArrCon2000, MorLatRam2013, DamKucRajAta2016, PanGua2016)

Project.toml

@@ -2,10 +2,11 @@ name = "UnitCommitment"
 uuid = "64606440-39ea-11e9-0f29-3303a1d3d877"
 authors = ["Santos Xavier, Alinson <axavier@anl.gov>"]
 repo = "https://github.com/ANL-CEEESA/UnitCommitment.jl"
-version = "0.2.1"
+version = "0.2.2"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
@@ -19,6 +20,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 [compat]
 Cbc = "0.7"
 DataStructures = "0.18"
+Distributions = "0.25"
 GZip = "0.5"
 JSON = "0.21"
 JuMP = "0.21"
@@ -29,6 +31,7 @@ julia = "1"
 [extras]
 Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b"
 
 [targets]
-test = ["Cbc", "Test"]
+test = ["Cbc", "Test", "Gurobi"]

README.md

@@ -20,7 +20,7 @@
 
 * **Data Format:** The package proposes an extensible and fully-documented JSON-based data format for SCUC, developed in collaboration with Independent System Operators (ISOs), which describes the most important aspects of the problem. The format supports the most common generator characteristics (including ramping, piecewise-linear production cost curves and time-dependent startup costs), as well as operating reserves, price-sensitive loads, transmission networks and contingencies.
 * **Benchmark Instances:** The package provides a diverse collection of large-scale benchmark instances collected from the literature, converted into a common data format, and extended using data-driven methods to make them more challenging and realistic.
-* **Model Implementation**: The package provides a Julia/JuMP implementations of state-of-the-art formulations and solution methods for SCUC, including multiple ramping formulations ([ArrCon2000][ArrCon2000], [MorLatRam2013][MorLatRam2013], [DamKucRajAta2016][DamKucRajAta2016], [PanGua2016][PanGua2016]), multiple piecewise-linear costs formulations ([Gar1962][Gar1962], [CarArr2006][CarArr2006], [KnuOstWat2018][KnuOstWat2018]) and contingency screening methods ([XavQiuWanThi2019][XavQiuWanThi2019]). Our goal is to keep these implementations up-to-date as new methods are proposed in the literature.
+* **Model Implementation**: The package provides Julia/JuMP implementations of state-of-the-art formulations and solution methods for SCUC, including multiple ramping formulations ([ArrCon2000][ArrCon2000], [MorLatRam2013][MorLatRam2013], [DamKucRajAta2016][DamKucRajAta2016], [PanGua2016][PanGua2016]), multiple piecewise-linear costs formulations ([Gar1962][Gar1962], [CarArr2006][CarArr2006], [KnuOstWat2018][KnuOstWat2018]) and contingency screening methods ([XavQiuWanThi2019][XavQiuWanThi2019]). Our goal is to keep these implementations up-to-date as new methods are proposed in the literature.
 * **Benchmark Tools:** The package provides automated benchmark scripts to accurately evaluate the performance impact of proposed code changes.
 
 [ArrCon2000]: https://doi.org/10.1109/59.871739

benchmark/scripts/table.py

@@ -6,6 +6,9 @@ from pathlib import Path
 import pandas as pd
 import re
 from tabulate import tabulate
+from colorama import init, Fore, Back, Style
+
+init()
 
 
 def process_all_log_files():
@@ -48,6 +51,7 @@ def process(filename):
     # m = re.search("case([0-9]*)", instance_name)
     # n_buses = int(m.group(1))
     n_buses = 0
+    validation_errors = 0
 
     with open(filename) as file:
         for line in file.readlines():
@@ -137,6 +141,14 @@ def process(filename):
             if m is not None:
                 transmission_count += 1
 
+            m = re.search(r".*Found ([0-9]*) validation errors", line)
+            if m is not None:
+                validation_errors += int(m.group(1))
+                print(
+                    f"{Fore.YELLOW}{Style.BRIGHT}Warning:{Style.RESET_ALL} {validation_errors:8d} "
+                    f"{Style.DIM}validation errors in {Style.RESET_ALL}{group_name}/{instance_name}/{sample_name}"
+                )
+
     return {
         "Group": group_name,
         "Instance": instance_name,
@@ -168,6 +180,7 @@ def process(filename):
         "Transmission screening constraints": transmission_count,
         "Transmission screening time": transmission_time,
         "Transmission screening calls": transmission_calls,
+        "Validation errors": validation_errors,
     }
 
 

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 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
     }
 }
 ```

docs/model.md

@@ -148,7 +148,7 @@ for g in instance.units
 end
 ```
 
-### Modifying the model
+### Fixing variables, modifying objective function and adding constraints
 
 Since we now have a direct reference to the JuMP decision variables, it is possible to fix variables, change the coefficients in the objective function, or even add new constraints to the model before solving it. The script below shows how can this be accomplished. For more information on modifying an existing model, [see the JuMP documentation](https://jump.dev/JuMP.jl/stable/manual/variables/).
 
@@ -190,6 +190,54 @@ JuMP.set_objective_coefficient(
 UnitCommitment.optimize!(model)
 ```
 
+### Adding new component to a bus
+
+The following snippet shows how to add a new grid component to a particular bus. For each time step, we create decision variables for the new grid component, add these variables to the objective function, then attach the component to a particular bus by modifying some existing model constraints. 
+
+```julia
+using Cbc
+using JuMP
+using UnitCommitment
+
+# Load instance and build base model
+instance = UnitCommitment.read_benchmark("matpower/case118/2017-02-01")
+model = UnitCommitment.build_model(
+    instance=instance,
+    optimizer=Cbc.Optimizer,
+)
+
+# Get the number of time steps in the original instance
+T = instance.time
+
+# Create decision variables for the new grid component.
+# In this example, we assume that the new component can
+# inject up to 10 MW of power at each time step, so we
+# create new continuous variables 0 ≤ x[t] ≤ 10.
+@variable(model, x[1:T], lower_bound=0.0, upper_bound=10.0)
+
+# For each time step
+for t in 1:T
+
+    # Add production costs to the objective function.
+    # In this example, we assume a cost of $5/MW.
+    set_objective_coefficient(model, x[t], 5.0)
+
+    # Attach the new component to bus b1, by modifying the
+    # constraint `eq_net_injection`.
+    set_normalized_coefficient(
+        model[:eq_net_injection]["b1", t],
+        x[t],
+        1.0,
+    )
+end
+
+# Solve the model
+UnitCommitment.optimize!(model)
+
+# Show optimal values for the x variables
+@show value.(x)
+```
+
 References
 ----------
 * [KnOsWa20] **Bernard Knueven, James Ostrowski and Jean-Paul Watson.** "On Mixed-Integer Programming Formulations for the Unit Commitment Problem". INFORMS Journal on Computing (2020). [DOI: 10.1287/ijoc.2019.0944](https://doi.org/10.1287/ijoc.2019.0944)

docs/usage.md

@@ -71,7 +71,7 @@ Advanced usage
 
 ### Customizing the formulation
 
-By default, `build_model` uses a formulation that combines modeling components from different publications, and that has been carefully tested, using our own benchmark scripts, to provide good performance across a wide variety of instances. This default formulation is expected to change over time, as new methods are proposed in the literature. You can, however, construct your own formulation, based the modeling components that you choose, as shown in the next example.
+By default, `build_model` uses a formulation that combines modeling components from different publications, and that has been carefully tested, using our own benchmark scripts, to provide good performance across a wide variety of instances. This default formulation is expected to change over time, as new methods are proposed in the literature. You can, however, construct your own formulation, based on the modeling components that you choose, as shown in the next example.
 
 ```julia
 using Cbc
@@ -119,7 +119,10 @@ instance = UnitCommitment.read("instance.json")
 UnitCommitment.generate_initial_conditions!(instance, Cbc.Optimizer)
 
 # Construct and solve optimization model
-model = UnitCommitment.build_model(instance, Cbc.Optimizer)
+model = UnitCommitment.build_model(
+    instance=instance,
+    optimizer=Cbc.Optimizer,
+)
 UnitCommitment.optimize!(model)
 ```
 

src/UnitCommitment.jl

@@ -30,6 +30,8 @@ include("model/formulations/base/unit.jl")
 include("model/formulations/CarArr2006/pwlcosts.jl")
 include("model/formulations/DamKucRajAta2016/ramp.jl")
 include("model/formulations/Gar1962/pwlcosts.jl")
+include("model/formulations/Gar1962/status.jl")
+include("model/formulations/Gar1962/prod.jl")
 include("model/formulations/KnuOstWat2018/pwlcosts.jl")
 include("model/formulations/MorLatRam2013/ramp.jl")
 include("model/formulations/MorLatRam2013/scosts.jl")
@@ -46,6 +48,7 @@ include("solution/warmstart.jl")
 include("solution/write.jl")
 include("transform/initcond.jl")
 include("transform/slice.jl")
+include("transform/randomize.jl")
 include("utils/log.jl")
 include("validation/repair.jl")
 include("validation/validate.jl")

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,10 +2,21 @@
 # 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,
-    formulation::ArrCon2000.Ramping,
+    formulation_prod_vars::Gar1962.ProdVars,
+    formulation_ramping::ArrCon2000.Ramping,
+    formulation_status_vars::Gar1962.StatusVars,
 )::Nothing
     # TODO: Move upper case constants to model[:instance]
     RESERVES_WHEN_START_UP = true
@@ -17,14 +28,19 @@ function _add_ramp_eqs!(
     RD = g.ramp_down_limit
     SU = g.startup_limit
     SD = g.shutdown_limit
-    is_on = model[:is_on]
-    prod_above = model[:prod_above]
     reserve = model[:reserve]
-    switch_off = model[:switch_off]
-    switch_on = model[:switch_on]
     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]
+
+    # Gar1962.StatusVars
+    is_on = model[:is_on]
+    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
@@ -50,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 <=
@@ -81,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,21 +2,32 @@
 # 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,
-    formulation::CarArr2006.PwlCosts,
+    formulation_prod_vars::Gar1962.ProdVars,
+    formulation_pwl_costs::CarArr2006.PwlCosts,
+    formulation_status_vars::StatusVarsFormulation,
 )::Nothing
     eq_prod_above_def = _init(model, :eq_prod_above_def)
     eq_segprod_limit = _init(model, :eq_segprod_limit)
-    prod_above = model[:prod_above]
     segprod = model[:segprod]
     gn = g.name
+
+    # Gar1962.ProdVars
+    prod_above = model[:prod_above]
+
     K = length(g.cost_segments)
     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
@@ -31,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,28 +2,56 @@
 # 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,
-    formulation::DamKucRajAta2016.Ramping,
+    formulation_prod_vars::Gar1962.ProdVars,
+    formulation_ramping::DamKucRajAta2016.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
-    known_initial_conditions = true
-    is_initially_on = (g.initial_status > 0)
-    SU = g.startup_limit
-    SD = g.shutdown_limit
-    RU = g.ramp_up_limit
-    RD = g.ramp_down_limit
+    is_initially_on = _is_initially_on(g)
+
+    # The following are the same for generator g across all time periods
+    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
+
     gn = g.name
     eq_str_ramp_down = _init(model, :eq_str_ramp_down)
     eq_str_ramp_up = _init(model, :eq_str_ramp_up)
-    is_on = model[:is_on]
-    prod_above = model[:prod_above]
     reserve = model[:reserve]
+
+    # Gar1962.ProdVars
+    prod_above = model[:prod_above]
+
+    # Gar1962.StatusVars
+    is_on = model[:is_on]
     switch_off = model[:switch_off]
     switch_on = model[:switch_on]
 
@@ -50,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,
@@ -70,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,
@@ -88,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 <=
@@ -104,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

@@ -0,0 +1,72 @@
+# 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_production_vars!(model, unit, formulation_prod_vars)
+
+Adds symbols identified by `Gar1962.ProdVars` to `model`.
+"""
+function _add_production_vars!(
+    model::JuMP.Model,
+    g::Unit,
+    formulation_prod_vars::Gar1962.ProdVars,
+)::Nothing
+    prod_above = _init(model, :prod_above)
+    segprod = _init(model, :segprod)
+    for t in 1:model[:instance].time
+        for k in 1:length(g.cost_segments)
+            segprod[g.name, t, k] = @variable(model, lower_bound = 0)
+        end
+        prod_above[g.name, t] = @variable(model, lower_bound = 0)
+    end
+    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,
+    formulation_prod_vars::Gar1962.ProdVars,
+)::Nothing
+    eq_prod_limit = _init(model, :eq_prod_limit)
+    is_on = model[:is_on]
+    prod_above = model[:prod_above]
+    reserve = model[:reserve]
+    gn = g.name
+    for t in 1:model[:instance].time
+        # Objective function terms for production costs
+        # 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 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 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
+            power_diff = 0.0
+        end
+        eq_prod_limit[gn, t] = @constraint(
+            model,
+            prod_above[gn, t] + reserve[gn, t] <= power_diff * is_on[gn, t]
+        )
+    end
+end

src/model/formulations/Gar1962/pwlcosts.jl

@@ -2,28 +2,57 @@
 # 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,
-    formulation::Gar1962.PwlCosts,
+    formulation_prod_vars::Gar1962.ProdVars,
+    formulation_pwl_costs::Gar1962.PwlCosts,
+    formulation_status_vars::Gar1962.StatusVars,
 )::Nothing
     eq_prod_above_def = _init(model, :eq_prod_above_def)
     eq_segprod_limit = _init(model, :eq_segprod_limit)
-    is_on = model[:is_on]
-    prod_above = model[:prod_above]
     segprod = model[:segprod]
     gn = g.name
+
+    # Gar1962.ProdVars
+    prod_above = model[:prod_above]
+
+    # Gar1962.StatusVars
+    is_on = model[:is_on]
+
     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
@@ -40,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

@@ -0,0 +1,160 @@
+# 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_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,
+    formulation_status_vars::Gar1962.StatusVars,
+)::Nothing
+    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
+        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,
+    formulation_status_vars::Gar1962.StatusVars,
+)::Nothing
+    eq_binary_link = _init(model, :eq_binary_link)
+    eq_switch_on_off = _init(model, :eq_switch_on_off)
+    is_on = model[:is_on]
+    switch_off = model[:switch_off]
+    switch_on = model[:switch_on]
+    for t in 1:model[:instance].time
+        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,
+                is_on[g.name, t] - _is_initially_on(g) ==
+                switch_on[g.name, t] - switch_off[g.name, t]
+            )
+        else
+            eq_binary_link[g.name, t] = @constraint(
+                model,
+                is_on[g.name, t] - is_on[g.name, t-1] ==
+                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
+    return
+end

src/model/formulations/Gar1962/structs.jl

@@ -14,7 +14,56 @@ Formulation described in:
 module Gar1962
 
 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
 
+"""
+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,24 +2,54 @@
 # 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,
-    formulation::KnuOstWat2018.PwlCosts,
+    formulation_prod_vars::Gar1962.ProdVars,
+    formulation_pwl_costs::KnuOstWat2018.PwlCosts,
+    formulation_status_vars::Gar1962.StatusVars,
 )::Nothing
     eq_prod_above_def = _init(model, :eq_prod_above_def)
     eq_segprod_limit_a = _init(model, :eq_segprod_limit_a)
     eq_segprod_limit_b = _init(model, :eq_segprod_limit_b)
     eq_segprod_limit_c = _init(model, :eq_segprod_limit_c)
-    prod_above = model[:prod_above]
     segprod = model[:segprod]
-    is_on = model[:is_on]
-    switch_on = model[:switch_on]
-    switch_off = model[:switch_off]
     gn = g.name
     K = length(g.cost_segments)
     T = model[:instance].time
 
+    # Gar1962.ProdVars
+    prod_above = model[:prod_above]
+
+    # Gar1962.StatusVars
+    is_on = model[:is_on]
+    switch_on = model[:switch_on]
+    switch_off = model[:switch_off]
+
     for t in 1:T
         for k in 1:K
             # Pbar^{k-1)
@@ -51,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] <=
@@ -60,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] <=
@@ -69,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] <=
@@ -80,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,10 +2,33 @@
 # 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,
-    formulation::MorLatRam2013.Ramping,
+    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
@@ -13,16 +36,20 @@ 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
-    SD = g.shutdown_limit
-    RU = g.ramp_up_limit
-    RD = g.ramp_down_limit
+    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
     gn = g.name
     eq_ramp_down = _init(model, :eq_ramp_down)
     eq_ramp_up = _init(model, :eq_str_ramp_up)
-    is_on = model[:is_on]
-    prod_above = model[:prod_above]
     reserve = model[:reserve]
+
+    # Gar1962.ProdVars
+    prod_above = model[:prod_above]
+
+    # Gar1962.StatusVars
+    is_on = model[:is_on]
     switch_off = model[:switch_off]
     switch_on = model[:switch_on]
 
@@ -65,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(
@@ -109,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,49 +2,87 @@
 # 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)
-    startup = model[:startup]
-    for t in 1:model[:instance].time
-        for s in 1:S
-            # If unit is switching on, we must choose a startup category
-            eq_startup_choose[g.name, t, s] = @constraint(
-                model,
-                model[:switch_on][g.name, t] ==
-                sum(startup[g.name, t, s] for s in 1:S)
-            )
 
+    # 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
+        # Equation (55) in Knueven et al. (2020)
+        eq_startup_choose[gn, t] = @constraint(
+            model,
+            switch_on[gn, t] == sum(startup[gn, t, s] for s in 1:S)
+        )
+
+        for s in 1:S
             # If unit has not switched off in the last `delay` time periods, startup category is forbidden.
             # The last startup category is always allowed.
             if s < S
                 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] <=
-                    initial_sum + sum(
-                        model[:switch_off][g.name, i] for i in range if i >= 1
-                    )
+                    startup[gn, t, s] <=
+                    initial_sum +
+                    sum(switch_off[gn, i] for i in range if i >= 1)
                 )
-            end
+            end # if s < S (not the last category)
 
             # 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
+        end # iterate over startup categories
+    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

@@ -1,16 +1,40 @@
+# 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!
+
+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,
-    formulation::PanGua2016.Ramping,
+    formulation_prod_vars::Gar1962.ProdVars,
+    formulation_ramping::PanGua2016.Ramping,
+    formulation_status_vars::Gar1962.StatusVars,
 )::Nothing
     # TODO: Move upper case constants to model[:instance]
     RESERVES_WHEN_SHUT_DOWN = true
     gn = g.name
-    is_on = model[:is_on]
-    prod_above = model[:prod_above]
     reserve = model[:reserve]
-    switch_off = model[:switch_off]
-    switch_on = model[:switch_on]
     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)
@@ -23,6 +47,14 @@ function _add_ramp_eqs!(
     RD = g.ramp_down_limit # ramp down rate
     T = model[:instance].time
 
+    # Gar1962.ProdVars
+    prod_above = model[:prod_above]
+
+    # Gar1962.StatusVars
+    is_on = model[:is_on]
+    switch_off = model[:switch_off]
+    switch_on = model[:switch_on]
+
     for t in 1:T
         Pbar = g.max_power[t]
         if Pbar < 1e-7
@@ -31,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
@@ -55,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,
@@ -73,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

@@ -6,20 +6,42 @@ abstract type TransmissionFormulation end
 abstract type RampingFormulation end
 abstract type PiecewiseLinearCostsFormulation end
 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
     ramping::RampingFormulation
     startup_costs::StartupCostsFormulation
+    status_vars::StatusVarsFormulation
     transmission::TransmissionFormulation
 
     function Formulation(;
+        prod_vars::ProductionVarsFormulation = Gar1962.ProdVars(),
         pwl_costs::PiecewiseLinearCostsFormulation = KnuOstWat2018.PwlCosts(),
         ramping::RampingFormulation = MorLatRam2013.Ramping(),
         startup_costs::StartupCostsFormulation = MorLatRam2013.StartupCosts(),
+        status_vars::StatusVarsFormulation = Gar1962.StatusVars(),
         transmission::TransmissionFormulation = ShiftFactorsFormulation(),
     )
-        return new(pwl_costs, ramping, startup_costs, transmission)
+        return new(
+            prod_vars,
+            pwl_costs,
+            ramping,
+            startup_costs,
+            status_vars,
+            transmission,
+        )
     end
 end
 

src/model/formulations/base/system.jl

@@ -2,21 +2,41 @@
 # 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)
-    eq_net_injection_def = _init(model, :eq_net_injection_def)
+    eq_net_injection = _init(model, :eq_net_injection)
     eq_power_balance = _init(model, :eq_power_balance)
     for t in 1:T, b in model[:instance].buses
         n = net_injection[b.name, t] = @variable(model)
-        eq_net_injection_def[t, b.name] =
-            @constraint(model, n == model[:expr_net_injection][b.name, t])
+        eq_net_injection[b.name, t] =
+            @constraint(model, -n + model[:expr_net_injection][b.name, t] == 0)
     end
     for t in 1:T
         eq_power_balance[t] = @constraint(
@@ -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(
-                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 >= 0 ? model[:reserve_shortfall][t] : 0.0) >=
+            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,7 +2,16 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
-function _add_unit!(model::JuMP.Model, g::Unit, f::Formulation)
+"""
+    _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")
     end
@@ -11,84 +20,86 @@ function _add_unit!(model::JuMP.Model, g::Unit, f::Formulation)
     end
 
     # Variables
-    _add_production_vars!(model, g)
+    _add_production_vars!(model, g, formulation.prod_vars)
     _add_reserve_vars!(model, g)
     _add_startup_shutdown_vars!(model, g)
-    _add_status_vars!(model, g)
+    _add_status_vars!(model, g, formulation.status_vars)
 
     # Constraints and objective function
     _add_min_uptime_downtime_eqs!(model, g)
     _add_net_injection_eqs!(model, g)
-    _add_production_limit_eqs!(model, g)
-    _add_production_piecewise_linear_eqs!(model, g, f.pwl_costs)
-    _add_ramp_eqs!(model, g, f.ramping)
-    _add_startup_cost_eqs!(model, g, f.startup_costs)
-    _add_startup_shutdown_limit_eqs!(model, g)
-    _add_status_eqs!(model, g)
+    _add_production_limit_eqs!(model, g, formulation.prod_vars)
+    _add_production_piecewise_linear_eqs!(
+        model,
+        g,
+        formulation.prod_vars,
+        formulation.pwl_costs,
+        formulation.status_vars,
+    )
+    _add_ramp_eqs!(
+        model,
+        g,
+        formulation.prod_vars,
+        formulation.ramping,
+        formulation.status_vars,
+    )
+    _add_startup_cost_eqs!(model, g, formulation.startup_costs)
+    _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_production_vars!(model::JuMP.Model, g::Unit)::Nothing
-    prod_above = _init(model, :prod_above)
-    segprod = _init(model, :segprod)
-    for t in 1:model[:instance].time
-        for k in 1:length(g.cost_segments)
-            segprod[g.name, t, k] = @variable(model, lower_bound = 0)
-        end
-        prod_above[g.name, t] = @variable(model, lower_bound = 0)
-    end
-    return
-end
+"""
+    _add_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
 
-function _add_production_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
-    eq_prod_limit = _init(model, :eq_prod_limit)
-    is_on = model[:is_on]
-    prod_above = model[:prod_above]
-    reserve = model[:reserve]
-    gn = g.name
-    for t in 1:model[:instance].time
-        # Objective function terms for production costs
-        # Part of (69) of Kneuven et al. (2020) as C^R_g * u_g(t) term
-        add_to_expression!(model[:obj], is_on[gn, t], g.min_power_cost[t])
-
-        # Production limit
-        # Equation (18) in Kneuven et al. (2020)
-        #   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)
-        #      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
-            power_diff = 0.0
-        end
-        eq_prod_limit[gn, t] = @constraint(
-            model,
-            prod_above[gn, t] + reserve[gn, t] <= power_diff * is_on[gn, t]
-        )
-    end
-end
-
-function _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
-            reserve[g.name, t] = 0.0
+            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
+        reserve_shortfall[t] =
+            (model[:instance].shortfall_penalty[t] >= 0) ?
+            @variable(model, lower_bound = 0) : 0.0
     end
     return
 end
 
+"""
+    _add_reserve_eqs!(model::JuMP.Model, g::Unit)::Nothing
+"""
 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
+    # 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
@@ -99,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
@@ -118,8 +146,15 @@ function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
         )
         # Shutdown limit
         if g.initial_power > g.shutdown_limit
-            eq_shutdown_limit[g.name, 0] =
-                @constraint(model, switch_off[g.name, 1] <= 0)
+            # 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, 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(
@@ -134,56 +169,32 @@ function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
     return
 end
 
-function _add_status_vars!(model::JuMP.Model, g::Unit)::Nothing
-    is_on = _init(model, :is_on)
-    switch_on = _init(model, :switch_on)
-    switch_off = _init(model, :switch_off)
-    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)
-        end
-    end
-    return
-end
+"""
+    _add_shutdown_cost_eqs!(model::JuMP.Model, g::Unit)::Nothing
 
-function _add_status_eqs!(model::JuMP.Model, g::Unit)::Nothing
-    eq_binary_link = _init(model, :eq_binary_link)
-    eq_switch_on_off = _init(model, :eq_switch_on_off)
-    is_on = model[:is_on]
-    switch_off = model[:switch_off]
-    switch_on = model[:switch_on]
-    for t in 1:model[:instance].time
-        if !g.must_run[t]
-            # Link binary variables
-            if t == 1
-                eq_binary_link[g.name, t] = @constraint(
-                    model,
-                    is_on[g.name, t] - _is_initially_on(g) ==
-                    switch_on[g.name, t] - switch_off[g.name, t]
-                )
-            else
-                eq_binary_link[g.name, t] = @constraint(
-                    model,
-                    is_on[g.name, t] - is_on[g.name, t-1] ==
-                    switch_on[g.name, t] - switch_off[g.name, t]
-                )
-            end
-            # Cannot switch on and off at the same time
-            eq_switch_on_off[g.name, t] = @constraint(
-                model,
-                switch_on[g.name, t] + switch_off[g.name, t] <= 1
+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
-    return
+    end # loop over time
 end
 
+"""
+    _add_ramp_eqs!(model, unit, formulation)
+"""
 function _add_ramp_eqs!(
     model::JuMP.Model,
     g::Unit,
@@ -230,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]
@@ -239,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(
@@ -273,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
@@ -287,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/transform/randomize.jl

@@ -0,0 +1,53 @@
+# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
+# Copyright (C) 2020-2021, UChicago Argonne, LLC. All rights reserved.
+# Released under the modified BSD license. See COPYING.md for more details.
+
+using Distributions
+
+function randomize_unit_costs!(
+    instance::UnitCommitmentInstance;
+    distribution = Uniform(0.95, 1.05),
+)::Nothing
+    for unit in instance.units
+        α = rand(distribution)
+        unit.min_power_cost *= α
+        for k in unit.cost_segments
+            k.cost *= α
+        end
+        for s in unit.startup_categories
+            s.cost *= α
+        end
+    end
+    return
+end
+
+function randomize_load_distribution!(
+    instance::UnitCommitmentInstance;
+    distribution = Uniform(0.90, 1.10),
+)::Nothing
+    α = rand(distribution, length(instance.buses))
+    for t in 1:instance.time
+        total = sum(bus.load[t] for bus in instance.buses)
+        den = sum(
+            bus.load[t] / total * α[i] for
+            (i, bus) in enumerate(instance.buses)
+        )
+        for (i, bus) in enumerate(instance.buses)
+            bus.load[t] *= α[i] / den
+        end
+    end
+    return
+end
+
+function randomize_peak_load!(
+    instance::UnitCommitmentInstance;
+    distribution = Uniform(0.925, 1.075),
+)::Nothing
+    α = rand(distribution)
+    for bus in instance.buses
+        bus.load *= α
+    end
+    return
+end
+
+export randomize_unit_costs!, randomize_load_distribution!, randomize_peak_load!

src/utils/sysimage.jl

@@ -5,12 +5,20 @@
 using PackageCompiler
 
 using DataStructures
+using Distributions
 using JSON
 using JuMP
 using MathOptInterface
 using SparseArrays
 
-pkg = [:DataStructures, :JSON, :JuMP, :MathOptInterface, :SparseArrays]
+pkg = [
+    :DataStructures,
+    :Distributions,
+    :JSON,
+    :JuMP,
+    :MathOptInterface,
+    :SparseArrays,
+]
 
 @info "Building system image..."
 create_sysimage(

src/validation/validate.jl

@@ -208,12 +208,8 @@ function _validate_units(instance, solution; tol = 0.01)
                         break
                     end
                 end
-                if t == time_down + 1
-                    initial_down = unit.min_downtime
-                    if unit.initial_status < 0
-                        initial_down = -unit.initial_status
-                    end
-                    time_down += initial_down
+                if (t == time_down + 1) && (unit.initial_status < 0)
+                    time_down -= unit.initial_status
                 end
 
                 # Calculate startup costs
@@ -246,14 +242,6 @@ function _validate_units(instance, solution; tol = 0.01)
                         break
                     end
                 end
-                if t == time_up + 1
-                    initial_up = unit.min_uptime
-                    if unit.initial_status > 0
-                        initial_up = unit.initial_status
-                    end
-                    time_up += initial_up
-                end
-
                 if (t == time_up + 1) && (unit.initial_status > 0)
                     time_up += unit.initial_status
                 end
@@ -336,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

@@ -3,6 +3,7 @@
 # Released under the modified BSD license. See COPYING.md for more details.
 
 using UnitCommitment
+using JuMP
 import UnitCommitment:
     ArrCon2000,
     CarArr2006,
@@ -11,17 +12,55 @@ import UnitCommitment:
     Gar1962,
     KnuOstWat2018,
     MorLatRam2013,
-    PanGua2016
+    PanGua2016,
+    XavQiuWanThi2019
 
-function _test(formulation::Formulation)::Nothing
-    instance = UnitCommitment.read_benchmark("matpower/case118/2017-02-01")
-    UnitCommitment.build_model(instance = instance, formulation = formulation)  # should not crash
+if ENABLE_LARGE_TESTS
+    using Gurobi
+end
+
+function _small_test(formulation::Formulation)::Nothing
+    instances = ["matpower/case118/2017-02-01", "test/case14"]
+    for instance in instances
+        # Should not crash
+        UnitCommitment.build_model(
+            instance = UnitCommitment.read_benchmark(instance),
+            formulation = formulation,
+        )
+    end
     return
 end
 
+function _large_test(formulation::Formulation)::Nothing
+    instances = ["pglib-uc/ca/Scenario400_reserves_1"]
+    for instance in instances
+        instance = UnitCommitment.read_benchmark(instance)
+        model = UnitCommitment.build_model(
+            instance = instance,
+            formulation = formulation,
+            optimizer = Gurobi.Optimizer,
+        )
+        UnitCommitment.optimize!(
+            model,
+            XavQiuWanThi2019.Method(two_phase_gap = false, gap_limit = 0.1),
+        )
+        solution = UnitCommitment.solution(model)
+        @test UnitCommitment.validate(instance, solution)
+    end
+    return
+end
+
+function _test(formulation::Formulation)::Nothing
+    _small_test(formulation)
+    if ENABLE_LARGE_TESTS
+        _large_test(formulation)
+    end
+end
+
 @testset "formulations" begin
+    _test(Formulation())
     _test(Formulation(ramping = ArrCon2000.Ramping()))
-    _test(Formulation(ramping = DamKucRajAta2016.Ramping()))
+    # _test(Formulation(ramping = DamKucRajAta2016.Ramping()))
     _test(
         Formulation(
             ramping = MorLatRam2013.Ramping(),

test/runtests.jl

@@ -7,6 +7,8 @@ using UnitCommitment
 
 UnitCommitment._setup_logger()
 
+const ENABLE_LARGE_TESTS = ("UCJL_LARGE_TESTS" in keys(ENV))
+
 @testset "UnitCommitment" begin
     include("usage.jl")
     @testset "import" begin
@@ -26,6 +28,7 @@ UnitCommitment._setup_logger()
     @testset "transform" begin
         include("transform/initcond_test.jl")
         include("transform/slice_test.jl")
+        include("transform/randomize_test.jl")
     end
     @testset "validation" begin
         include("validation/repair_test.jl")

test/transform/randomize_test.jl

@@ -0,0 +1,43 @@
+# 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.
+
+using UnitCommitment, Cbc, JuMP
+
+_get_instance() = UnitCommitment.read_benchmark("matpower/case118/2017-02-01")
+_total_load(instance) = sum(b.load[1] for b in instance.buses)
+
+@testset "randomize_unit_costs!" begin
+    instance = _get_instance()
+    unit = instance.units[10]
+    prev_min_power_cost = unit.min_power_cost
+    prev_prod_cost = unit.cost_segments[1].cost
+    prev_startup_cost = unit.startup_categories[1].cost
+    randomize_unit_costs!(instance)
+    @test prev_min_power_cost != unit.min_power_cost
+    @test prev_prod_cost != unit.cost_segments[1].cost
+    @test prev_startup_cost != unit.startup_categories[1].cost
+end
+
+@testset "randomize_load_distribution!" begin
+    instance = _get_instance()
+    bus = instance.buses[1]
+    prev_load = instance.buses[1].load[1]
+    prev_total_load = _total_load(instance)
+    randomize_load_distribution!(instance)
+    curr_total_load = _total_load(instance)
+    @test prev_load != instance.buses[1].load[1]
+    @test abs(prev_total_load - curr_total_load) < 1e-3
+end
+
+@testset "randomize_peak_load!" begin
+    instance = _get_instance()
+    bus = instance.buses[1]
+    prev_total_load = _total_load(instance)
+    prev_share = bus.load[1] / prev_total_load
+    randomize_peak_load!(instance)
+    curr_total_load = _total_load(instance)
+    curr_share = bus.load[1] / prev_total_load
+    @test curr_total_load != prev_total_load
+    @test abs(curr_share - prev_share) < 1e-3
+end