Reformat source code

.gitignore

@@ -1,21 +1,38 @@
 *.bak
 *.gz
-*.lastrun
-*.so
-*.mps
 *.ipynb
+*.lastrun
+*.mps
+*.so
+*/Manifest.toml
+.AppleDB
+.AppleDesktop
+.AppleDouble
+.DS_Store
+.DocumentRevisions-V100
+.LSOverride
+.Spotlight-V100
+.TemporaryItems
+.Trashes
+.VolumeIcon.icns
+._*
+.apdisk
+.com.apple.timemachine.donotpresent
+.fseventsd
 .ipy*
+.vscode
+Icon
+Manifest.toml
+Network Trash Folder
+TODO.md
+Temporary Items
 benchmark/results
 benchmark/runs
 benchmark/tables
 benchmark/tmp.json
 build
+docs/_build
 instances/**/*.json
 instances/_source
 local
 notebooks
-TODO.md
-docs/_build
-.vscode
-Manifest.toml
-*/Manifest.toml

CHANGELOG.md

@@ -11,6 +11,17 @@ 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
 
+## [Unreleased]
+### Added
+- Add multiple reserve products
+
+### Changed
+- To support multiple reserve products, the input data format has been modified as follows:
+    - In `Generators`, replace `Provides spinning reserves?` by `Reserve eligibility`
+    - In `Parameters`, remove `Reserve shortfall penalty`
+    - Revise `Reserves` section
+
+
 ## [0.2.2] - 2021-07-21
 ### Fixed
 - Fix small bug in validation scripts related to startup costs

Makefile

@@ -5,7 +5,7 @@
 VERSION := 0.2
 
 clean:
-	rm -rfv build
+	rm -rfv build Manifest.toml test/Manifest.toml deps/formatter/build deps/formatter/Manifest.toml
 
 docs:
 	cd docs; make clean; make dirhtml

Project.toml

@@ -2,7 +2,7 @@ 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.2"
+version = "0.3.0"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
@@ -17,6 +17,7 @@ MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
@@ -24,7 +25,7 @@ DataStructures = "0.18"
 Distributions = "0.25"
 GZip = "0.5"
 JSON = "0.21"
-JuMP = "0.21"
-MathOptInterface = "0.9"
+JuMP = "1"
+MathOptInterface = "1"
 PackageCompiler = "1"
 julia = "1"

README.md

@@ -95,6 +95,7 @@ UnitCommitment.write("/tmp/output.json", solution)
 ## Authors
 * **Alinson S. Xavier** (Argonne National Laboratory)
 * **Aleksandr M. Kazachkov** (University of Florida)
+* **Ogün Yurdakul** (Technische Universität Berlin)
 * **Feng Qiu** (Argonne National Laboratory)
 
 ## Acknowledgments

docs/format.md

@@ -28,14 +28,13 @@ Each section is described in detail below. For a complete example, see [case14](
 
 ### Parameters
 
-This section describes system-wide parameters, such as power balance and reserve shortfall penalties, and optimization parameters, such as the length of the planning horizon and the time.
+This section describes system-wide parameters, such as power balance penalty, and optimization parameters, such as the length of the planning horizon and the time.
 
 | Key                            | Description                                       | Default  | Time series?
 | :----------------------------- | :------------------------------------------------ | :------: | :------------:
 | `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
@@ -44,7 +43,6 @@ This section describes system-wide parameters, such as power balance and reserve
     "Parameters": {
         "Time horizon (h)": 4,
         "Power balance penalty ($/MW)": 1000.0,
-        "Reserve shortfall penalty ($/MW)": -1.0
     }
 }
 ```
@@ -96,7 +94,7 @@ This section describes all generators in the system, including thermal units, re
 | `Initial status (h)`  | If set to a positive number, indicates the amount of time (in hours) the generator has been on at the beginning of the simulation, and if set to a negative number, the amount of time the generator has been off. For example, if `Initial status (h)` is `-2`, this means that the generator was off since `-02:00` (h:min). The simulation starts at time `00:00`. If `Initial status (h)` is `3`, this means that the generator was on since `-03:00`. A value of zero is not acceptable. | Required | N
 | `Initial power (MW)`  | Amount of power the generator at time step `-1`, immediately before the planning horizon starts. | Required | N
 | `Must run?`               | If `true`, the generator should be committed, even if that is not economical (Boolean). | `false` | Y
-| `Provides spinning reserves?`    | If `true`, this generator may provide spinning reserves (Boolean). | `true` | Y
+| `Reserve eligibility` | List of reserve products this generator is eligibe to provide. By default, the generator is not eligible to provide any reserves. | `[]` | N
 
 #### Production costs and limits
 
@@ -135,13 +133,13 @@ Note that this curve also specifies the production limits. Specifically, the fir
             "Minimum uptime (h)": 4,
             "Initial status (h)": 12,
             "Must run?": false,
-            "Provides spinning reserves?": true,
+            "Reserve eligibility": ["r1"],
         },
         "gen2": {
             "Bus": "b5",
             "Production cost curve (MW)": [0.0, [10.0, 8.0, 0.0, 3.0]],
             "Production cost curve ($)": [0.0, 0.0],
-            "Provides spinning reserves?": true,
+            "Reserve eligibility": ["r1", "r2"],
         }
     }
 }
@@ -206,24 +204,39 @@ This section describes the characteristics of transmission system, such as its t
 
 ### Reserves
 
-This section describes the hourly amount of operating reserves required.
+This section describes the hourly amount of reserves required.
 
 
 | Key                   | Description                                        | Default   |  Time series?
 | :-------------------- | :------------------------------------------------- | --------- |  :----:
-| `Spinning (MW)`       | Minimum amount of system-wide spinning reserves (in MW). Only generators which are online may provide this reserve. | `0.0` | Y
+| `Type` | Type of reserve product. Must be either "spinning" or "flexiramp". | Required | N
+| `Amount (MW)` | Amount of reserves required. | Required | Y
+| `Shortfall penalty ($/MW)` | Penalty for shortage in meeting the reserve requirements (in $/MW). This is charged per time step. Negative value implies reserve constraints must always be satisfied. | `-1` | Y
 
-#### Example
+#### Example 1
 
 ```json
 {
     "Reserves": {
-        "Spinning (MW)": [
-            57.30552,
-            53.88429,
-            51.31838,
-            50.46307
-        ]
+        "r1": {
+            "Type": "spinning",
+            "Amount (MW)": [
+                57.30552,
+                53.88429,
+                51.31838,
+                50.46307
+            ],
+            "Shortfall penalty ($/MW)": 5.0
+        },
+        "r2": {
+            "Type": "flexiramp",
+            "Amount (MW)": [
+                20.31042,
+                23.65273,
+                27.41784,
+                25.34057
+            ],
+        }
     }
 }
 ```
@@ -286,9 +299,7 @@ The output data format is also JSON-based, but it is not currently documented si
 Current limitations
 -------------------
 
-* All reserves are system-wide. Zonal reserves are not currently supported.
 * Network topology remains the same for all time periods
 * Only N-1 transmission contingencies are supported. Generator contingencies are not currently supported.
 * Time-varying minimum production amounts are not currently compatible with ramp/startup/shutdown limits.
-
-
+* Flexible ramping products can only be acquired under the `WanHob2016` formulation, which does not support spinning reserves. 

docs/model.md

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

juliaw

@@ -47,7 +47,14 @@ project = TOML.parsefile("Project.toml")
 manifest = TOML.parsefile("Manifest.toml")
 deps = Symbol[]
 for dep in keys(project["deps"])
-    if "path" in keys(manifest[dep][1])
+    if dep in keys(manifest)
+	# Up to Julia 1.6
+        dep_entry = manifest[dep][1]
+    else
+        # Julia 1.7+
+        dep_entry = manifest["deps"][dep][1]
+    end
+    if "path" in keys(dep_entry)
         println("  - \$(dep) [skip]")
     else
         println("  - \$(dep)")

src/UnitCommitment.jl

@@ -16,6 +16,7 @@ include("model/formulations/KnuOstWat2018/structs.jl")
 include("model/formulations/MorLatRam2013/structs.jl")
 include("model/formulations/PanGua2016/structs.jl")
 include("solution/methods/XavQiuWanThi2019/structs.jl")
+include("model/formulations/WanHob2016/structs.jl")
 
 include("import/egret.jl")
 include("instance/read.jl")
@@ -36,6 +37,7 @@ include("model/formulations/KnuOstWat2018/pwlcosts.jl")
 include("model/formulations/MorLatRam2013/ramp.jl")
 include("model/formulations/MorLatRam2013/scosts.jl")
 include("model/formulations/PanGua2016/ramp.jl")
+include("model/formulations/WanHob2016/ramp.jl")
 include("model/jumpext.jl")
 include("solution/fix.jl")
 include("solution/methods/XavQiuWanThi2019/enforce.jl")

src/instance/read.jl

@@ -8,7 +8,7 @@ using DataStructures
 using GZip
 import Base: getindex, time
 
-const INSTANCES_URL = "https://axavier.org/UnitCommitment.jl/0.2/instances"
+const INSTANCES_URL = "https://axavier.org/UnitCommitment.jl/0.3/instances"
 
 """
     read_benchmark(name::AbstractString)::UnitCommitmentInstance
@@ -85,6 +85,7 @@ function _from_json(json; repair = true)
     contingencies = Contingency[]
     lines = TransmissionLine[]
     loads = PriceSensitiveLoad[]
+    reserves = Reserve[]
 
     function scalar(x; default = nothing)
         x !== nothing || return default
@@ -105,6 +106,7 @@ function _from_json(json; repair = true)
     name_to_bus = Dict{String,Bus}()
     name_to_line = Dict{String,TransmissionLine}()
     name_to_unit = Dict{String,Unit}()
+    name_to_reserve = Dict{String,Reserve}()
 
     function timeseries(x; default = nothing)
         x !== nothing || return default
@@ -117,6 +119,11 @@ function _from_json(json; repair = true)
         json["Parameters"]["Power balance penalty (\$/MW)"],
         default = [1000.0 for t in 1:T],
     )
+    # Penalty price for shortage in meeting system-wide flexiramp requirements
+    flexiramp_shortfall_penalty = timeseries(
+        json["Parameters"]["Flexiramp penalty (\$/MW)"],
+        default = [500.0 for t in 1:T],
+    )
     shortfall_penalty = timeseries(
         json["Parameters"]["Reserve shortfall penalty (\$/MW)"],
         default = [-1.0 for t in 1:T],
@@ -135,6 +142,24 @@ function _from_json(json; repair = true)
         push!(buses, bus)
     end
 
+    # Read reserves
+    if "Reserves" in keys(json)
+        for (reserve_name, dict) in json["Reserves"]
+            r = Reserve(
+                name = reserve_name,
+                type = lowercase(dict["Type"]),
+                amount = timeseries(dict["Amount (MW)"]),
+                units = [],
+                shortfall_penalty = scalar(
+                    dict["Shortfall penalty (\$/MW)"],
+                    default = -1,
+                ),
+            )
+            name_to_reserve[reserve_name] = r
+            push!(reserves, r)
+        end
+    end
+
     # Read units
     for (unit_name, dict) in json["Generators"]
         bus = name_to_bus[dict["Bus"]]
@@ -172,6 +197,13 @@ function _from_json(json; repair = true)
             )
         end
 
+        # Read reserve eligibility
+        unit_reserves = Reserve[]
+        if "Reserve eligibility" in keys(dict)
+            unit_reserves =
+                [name_to_reserve[n] for n in dict["Reserve eligibility"]]
+        end
+
         # Read and validate initial conditions
         initial_power = scalar(dict["Initial power (MW)"], default = nothing)
         initial_status = scalar(dict["Initial status (h)"], default = nothing)
@@ -205,24 +237,17 @@ function _from_json(json; repair = true)
             scalar(dict["Shutdown limit (MW)"], default = 1e6),
             initial_status,
             initial_power,
-            timeseries(
-                dict["Provides spinning reserves?"],
-                default = [true for t in 1:T],
-            ),
             startup_categories,
+            unit_reserves,
         )
         push!(bus.units, unit)
+        for r in unit_reserves
+            push!(r.units, unit)
+        end
         name_to_unit[unit_name] = unit
         push!(units, unit)
     end
 
-    # Read reserves
-    reserves = Reserves(zeros(T))
-    if "Reserves" in keys(json)
-        reserves.spinning =
-            timeseries(json["Reserves"]["Spinning (MW)"], default = zeros(T))
-    end
-
     # Read transmission lines
     if "Transmission lines" in keys(json)
         for (line_name, dict) in json["Transmission lines"]
@@ -295,7 +320,9 @@ function _from_json(json; repair = true)
         price_sensitive_loads_by_name = Dict(ps.name => ps for ps in loads),
         price_sensitive_loads = loads,
         reserves = reserves,
+        reserves_by_name = name_to_reserve,
         shortfall_penalty = shortfall_penalty,
+        flexiramp_shortfall_penalty = flexiramp_shortfall_penalty,
         time = T,
         units_by_name = Dict(g.name => g for g in units),
         units = units,

src/instance/structs.jl

@@ -20,6 +20,14 @@ mutable struct StartupCategory
     cost::Float64
 end
 
+Base.@kwdef mutable struct Reserve
+    name::String
+    type::String
+    amount::Vector{Float64}
+    units::Vector
+    shortfall_penalty::Float64
+end
+
 mutable struct Unit
     name::String
     bus::Bus
@@ -36,8 +44,8 @@ mutable struct Unit
     shutdown_limit::Float64
     initial_status::Union{Int,Nothing}
     initial_power::Union{Float64,Nothing}
-    provides_spinning_reserves::Vector{Bool}
     startup_categories::Vector{StartupCategory}
+    reserves::Vector{Reserve}
 end
 
 mutable struct TransmissionLine
@@ -52,10 +60,6 @@ mutable struct TransmissionLine
     flow_limit_penalty::Vector{Float64}
 end
 
-mutable struct Reserves
-    spinning::Vector{Float64}
-end
-
 mutable struct Contingency
     name::String
     lines::Vector{TransmissionLine}
@@ -79,8 +83,10 @@ Base.@kwdef mutable struct UnitCommitmentInstance
     power_balance_penalty::Vector{Float64}
     price_sensitive_loads_by_name::Dict{AbstractString,PriceSensitiveLoad}
     price_sensitive_loads::Vector{PriceSensitiveLoad}
-    reserves::Reserves
+    reserves::Vector{Reserve}
+    reserves_by_name::Dict{AbstractString,Reserve}
     shortfall_penalty::Vector{Float64}
+    flexiramp_shortfall_penalty::Vector{Float64}
     time::Int
     units_by_name::Dict{AbstractString,Unit}
     units::Vector{Unit}

src/model/formulations/ArrCon2000/ramp.jl

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

src/model/formulations/DamKucRajAta2016/ramp.jl

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

src/model/formulations/Gar1962/prod.jl

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

src/model/formulations/MorLatRam2013/ramp.jl

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

src/model/formulations/PanGua2016/ramp.jl

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

src/model/formulations/WanHob2016/ramp.jl

@@ -0,0 +1,177 @@
+# UnitCommitmentFL.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.
+
+function _add_ramp_eqs!(
+    model::JuMP.Model,
+    g::Unit,
+    ::Gar1962.ProdVars,
+    ::WanHob2016.Ramping,
+    ::Gar1962.StatusVars,
+)::Nothing
+    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
+    gn = g.name
+    minp = g.min_power
+    maxp = g.max_power
+    initial_power = g.initial_power
+
+    is_on = model[:is_on]
+    prod_above = model[:prod_above]
+    upflexiramp = model[:upflexiramp]
+    dwflexiramp = model[:dwflexiramp]
+    mfg = model[:mfg]
+
+    if length(g.reserves) > 1
+        error("Each generator may only provide one flexiramp reserve")
+    end
+    for r in g.reserves
+        if r.type !== "flexiramp"
+            error(
+                "This formulation only supports flexiramp reserves, not $(r.type)",
+            )
+        end
+        rn = r.name
+        for t in 1:model[:instance].time
+            @constraint(
+                model,
+                prod_above[gn, t] + (is_on[gn, t] * minp[t]) <= mfg[rn, gn, t]
+            ) # Eq. (19) in Wang & Hobbs (2016)
+            @constraint(model, mfg[rn, gn, t] <= is_on[gn, t] * maxp[t]) # Eq. (22) in Wang & Hobbs (2016)
+            if t != model[:instance].time
+                @constraint(
+                    model,
+                    minp[t] * (is_on[gn, t+1] + is_on[gn, t] - 1) <=
+                    prod_above[gn, t] - dwflexiramp[rn, gn, t] +
+                    (is_on[gn, t] * minp[t])
+                ) # first inequality of Eq. (20) in Wang & Hobbs (2016)
+                @constraint(
+                    model,
+                    prod_above[gn, t] - dwflexiramp[rn, gn, t] +
+                    (is_on[gn, t] * minp[t]) <=
+                    mfg[rn, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1]))
+                ) # second inequality of Eq. (20) in Wang & Hobbs (2016)
+                @constraint(
+                    model,
+                    minp[t] * (is_on[gn, t+1] + is_on[gn, t] - 1) <=
+                    prod_above[gn, t] +
+                    upflexiramp[rn, gn, t] +
+                    (is_on[gn, t] * minp[t])
+                ) # first inequality of Eq. (21) in Wang & Hobbs (2016)
+                @constraint(
+                    model,
+                    prod_above[gn, t] +
+                    upflexiramp[rn, gn, t] +
+                    (is_on[gn, t] * minp[t]) <=
+                    mfg[rn, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1]))
+                ) # second inequality of Eq. (21) in Wang & Hobbs (2016)
+                if t != 1
+                    @constraint(
+                        model,
+                        mfg[rn, gn, t] <=
+                        prod_above[gn, t-1] +
+                        (is_on[gn, t-1] * minp[t]) +
+                        (RU * is_on[gn, t-1]) +
+                        (SU * (is_on[gn, t] - is_on[gn, t-1])) +
+                        maxp[t] * (1 - is_on[gn, t])
+                    ) # Eq. (23) in Wang & Hobbs (2016)
+                    @constraint(
+                        model,
+                        (prod_above[gn, t-1] + (is_on[gn, t-1] * minp[t])) -
+                        (prod_above[gn, t] + (is_on[gn, t] * minp[t])) <=
+                        RD * is_on[gn, t] +
+                        SD * (is_on[gn, t-1] - is_on[gn, t]) +
+                        maxp[t] * (1 - is_on[gn, t-1])
+                    ) # Eq. (25) in Wang & Hobbs (2016)
+                else
+                    @constraint(
+                        model,
+                        mfg[rn, gn, t] <=
+                        initial_power +
+                        (RU * is_initially_on) +
+                        (SU * (is_on[gn, t] - is_initially_on)) +
+                        maxp[t] * (1 - is_on[gn, t])
+                    ) # Eq. (23) in Wang & Hobbs (2016) for the first time period
+                    @constraint(
+                        model,
+                        initial_power -
+                        (prod_above[gn, t] + (is_on[gn, t] * minp[t])) <=
+                        RD * is_on[gn, t] +
+                        SD * (is_initially_on - is_on[gn, t]) +
+                        maxp[t] * (1 - is_initially_on)
+                    ) # Eq. (25) in Wang & Hobbs (2016) for the first time period
+                end
+                @constraint(
+                    model,
+                    mfg[rn, gn, t] <=
+                    (SD * (is_on[gn, t] - is_on[gn, t+1])) +
+                    (maxp[t] * is_on[gn, t+1])
+                ) # Eq. (24) in Wang & Hobbs (2016)
+                @constraint(
+                    model,
+                    -RD * is_on[gn, t+1] -
+                    SD * (is_on[gn, t] - is_on[gn, t+1]) -
+                    maxp[t] * (1 - is_on[gn, t]) <= upflexiramp[rn, gn, t]
+                ) # first inequality of Eq. (26) in Wang & Hobbs (2016)
+                @constraint(
+                    model,
+                    upflexiramp[rn, gn, t] <=
+                    RU * is_on[gn, t] +
+                    SU * (is_on[gn, t+1] - is_on[gn, t]) +
+                    maxp[t] * (1 - is_on[gn, t+1])
+                ) # second inequality of Eq. (26) in Wang & Hobbs (2016)
+                @constraint(
+                    model,
+                    -RU * is_on[gn, t] - SU * (is_on[gn, t+1] - is_on[gn, t]) -
+                    maxp[t] * (1 - is_on[gn, t+1]) <= dwflexiramp[rn, gn, t]
+                ) # first inequality of Eq. (27) in Wang & Hobbs (2016)
+                @constraint(
+                    model,
+                    dwflexiramp[rn, gn, t] <=
+                    RD * is_on[gn, t+1] +
+                    SD * (is_on[gn, t] - is_on[gn, t+1]) +
+                    maxp[t] * (1 - is_on[gn, t])
+                ) # second inequality of Eq. (27) in Wang & Hobbs (2016)
+                @constraint(
+                    model,
+                    -maxp[t] * is_on[gn, t] + minp[t] * is_on[gn, t+1] <=
+                    upflexiramp[rn, gn, t]
+                ) # first inequality of Eq. (28) in Wang & Hobbs (2016)
+                @constraint(
+                    model,
+                    upflexiramp[rn, gn, t] <= maxp[t] * is_on[gn, t+1]
+                ) # second inequality of Eq. (28) in Wang & Hobbs (2016)
+                @constraint(
+                    model,
+                    -maxp[t] * is_on[gn, t+1] <= dwflexiramp[rn, gn, t]
+                ) # first inequality of Eq. (29) in Wang & Hobbs (2016)
+                @constraint(
+                    model,
+                    dwflexiramp[rn, gn, t] <=
+                    (maxp[t] * is_on[gn, t]) - (minp[t] * is_on[gn, t+1])
+                ) # second inequality of Eq. (29) in Wang & Hobbs (2016)
+            else
+                @constraint(
+                    model,
+                    mfg[rn, gn, t] <=
+                    prod_above[gn, t-1] +
+                    (is_on[gn, t-1] * minp[t]) +
+                    (RU * is_on[gn, t-1]) +
+                    (SU * (is_on[gn, t] - is_on[gn, t-1])) +
+                    maxp[t] * (1 - is_on[gn, t])
+                ) # Eq. (23) in Wang & Hobbs (2016) for the last time period
+                @constraint(
+                    model,
+                    (prod_above[gn, t-1] + (is_on[gn, t-1] * minp[t])) -
+                    (prod_above[gn, t] + (is_on[gn, t] * minp[t])) <=
+                    RD * is_on[gn, t] +
+                    SD * (is_on[gn, t-1] - is_on[gn, t]) +
+                    maxp[t] * (1 - is_on[gn, t-1])
+                ) # Eq. (25) in Wang & Hobbs (2016) for the last time period
+            end
+        end
+    end
+end

src/model/formulations/WanHob2016/structs.jl

@@ -0,0 +1,17 @@
+# UnitCommitmentFL.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.
+
+"""
+Formulation described in:
+    B. Wang and B. F. Hobbs, "Real-Time Markets for Flexiramp: A Stochastic 
+    Unit Commitment-Based Analysis," in IEEE Transactions on Power Systems, 
+    vol. 31, no. 2, pp. 846-860, March 2016, doi: 10.1109/TPWRS.2015.2411268.
+"""
+module WanHob2016
+
+import ..RampingFormulation
+
+struct Ramping <: RampingFormulation end
+
+end

src/model/formulations/base/system.jl

@@ -4,7 +4,8 @@
 
 function _add_system_wide_eqs!(model::JuMP.Model)::Nothing
     _add_net_injection_eqs!(model)
-    _add_reserve_eqs!(model)
+    _add_spinning_reserve_eqs!(model)
+    _add_flexiramp_reserve_eqs!(model)
     return
 end
 
@@ -27,29 +28,69 @@ function _add_net_injection_eqs!(model::JuMP.Model)::Nothing
     return
 end
 
-function _add_reserve_eqs!(model::JuMP.Model)::Nothing
-    eq_min_reserve = _init(model, :eq_min_reserve)
+function _add_spinning_reserve_eqs!(model::JuMP.Model)::Nothing
     instance = model[:instance]
-    for t in 1:instance.time
-        # Equation (68) in Kneuven et al. (2020)
-        # As in Morales-España et al. (2013a)
-        # Akin to the alternative formulation with max_power_avail
-        # from Carrión and Arroyo (2006) and Ostrowski et al. (2012)
-        shortfall_penalty = instance.shortfall_penalty[t]
-        eq_min_reserve[t] = @constraint(
-            model,
-            sum(model[: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],
+    eq_min_spinning_reserve = _init(model, :eq_min_spinning_reserve)
+    for r in instance.reserves
+        r.type == "spinning" || continue
+        for t in 1:instance.time
+            # Equation (68) in Kneuven et al. (2020)
+            # As in Morales-España et al. (2013a)
+            # Akin to the alternative formulation with max_power_avail
+            # from Carrión and Arroyo (2006) and Ostrowski et al. (2012)
+            eq_min_spinning_reserve[r.name, t] = @constraint(
+                model,
+                sum(model[:reserve][r.name, g.name, t] for g in r.units) +
+                model[:reserve_shortfall][r.name, t] >= r.amount[t]
             )
+
+            # Account for shortfall contribution to objective
+            if r.shortfall_penalty >= 0
+                add_to_expression!(
+                    model[:obj],
+                    r.shortfall_penalty,
+                    model[:reserve_shortfall][r.name, t],
+                )
+            end
+        end
+    end
+    return
+end
+
+function _add_flexiramp_reserve_eqs!(model::JuMP.Model)::Nothing
+    # Note: The flexpramp requirements in Wang & Hobbs (2016) are imposed as hard constraints 
+    #       through Eq. (17) and Eq. (18). The constraints eq_min_upflexiramp and eq_min_dwflexiramp 
+    #       provided below are modified versions of Eq. (17) and Eq. (18), respectively, in that   
+    #       they include slack variables for flexiramp shortfall, which are penalized in the
+    #       objective function.
+    eq_min_upflexiramp = _init(model, :eq_min_upflexiramp)
+    eq_min_dwflexiramp = _init(model, :eq_min_dwflexiramp)
+    instance = model[:instance]
+    for r in instance.reserves
+        r.type == "flexiramp" || continue
+        for t in 1:instance.time
+            # Eq. (17) in Wang & Hobbs (2016)
+            eq_min_upflexiramp[r.name, t] = @constraint(
+                model,
+                sum(model[:upflexiramp][r.name, g.name, t] for g in r.units) + model[:upflexiramp_shortfall][r.name, t] >= r.amount[t]
+            )
+            # Eq. (18) in Wang & Hobbs (2016)
+            eq_min_dwflexiramp[r.name, t] = @constraint(
+                model,
+                sum(model[:dwflexiramp][r.name, g.name, t] for g in r.units) + model[:dwflexiramp_shortfall][r.name, t] >= r.amount[t]
+            )
+
+            # Account for flexiramp shortfall contribution to objective
+            if r.shortfall_penalty >= 0
+                add_to_expression!(
+                    model[:obj],
+                    r.shortfall_penalty,
+                    (
+                        model[:upflexiramp_shortfall][r.name, t] +
+                        model[:dwflexiramp_shortfall][r.name, t]
+                    ),
+                )
+            end
         end
     end
     return

src/model/formulations/base/unit.jl

@@ -12,7 +12,8 @@ function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation)
 
     # Variables
     _add_production_vars!(model, g, formulation.prod_vars)
-    _add_reserve_vars!(model, g)
+    _add_spinning_reserve_vars!(model, g)
+    _add_flexiramp_reserve_vars!(model, g)
     _add_startup_shutdown_vars!(model, g)
     _add_status_vars!(model, g, formulation.status_vars)
 
@@ -42,26 +43,48 @@ end
 
 _is_initially_on(g::Unit)::Float64 = (g.initial_status > 0 ? 1.0 : 0.0)
 
-function _add_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
+function _add_spinning_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
     reserve = _init(model, :reserve)
     reserve_shortfall = _init(model, :reserve_shortfall)
-    for t in 1:model[:instance].time
-        if g.provides_spinning_reserves[t]
-            reserve[g.name, t] = @variable(model, lower_bound = 0)
-        else
-            reserve[g.name, t] = 0.0
+    for r in g.reserves
+        r.type == "spinning" || continue
+        for t in 1:model[:instance].time
+            reserve[r.name, g.name, t] = @variable(model, lower_bound = 0)
+            if (r.name, t) ∉ keys(reserve_shortfall)
+                reserve_shortfall[r.name, t] = @variable(model, lower_bound = 0)
+                if r.shortfall_penalty < 0
+                    set_upper_bound(reserve_shortfall[r.name, t], 0.0)
+                end
+            end
         end
-        reserve_shortfall[t] =
-            (model[:instance].shortfall_penalty[t] >= 0) ?
-            @variable(model, lower_bound = 0) : 0.0
     end
     return
 end
 
-function _add_reserve_eqs!(model::JuMP.Model, g::Unit)::Nothing
-    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)
+function _add_flexiramp_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
+    upflexiramp = _init(model, :upflexiramp)
+    upflexiramp_shortfall = _init(model, :upflexiramp_shortfall)
+    mfg = _init(model, :mfg)
+    dwflexiramp = _init(model, :dwflexiramp)
+    dwflexiramp_shortfall = _init(model, :dwflexiramp_shortfall)
+    for r in g.reserves
+        r.type == "flexiramp" || continue
+        for t in 1:model[:instance].time
+            # maximum feasible generation, \bar{g_{its}} in Wang & Hobbs (2016)
+            mfg[r.name, g.name, t] = @variable(model, lower_bound = 0)
+            upflexiramp[r.name, g.name, t] = @variable(model) # up-flexiramp, ur_{it} in Wang & Hobbs (2016)
+            dwflexiramp[r.name, g.name, t] = @variable(model) # down-flexiramp, dr_{it} in Wang & Hobbs (2016)
+            if (r.name, t) ∉ keys(upflexiramp_shortfall)
+                upflexiramp_shortfall[r.name, t] =
+                    @variable(model, lower_bound = 0)
+                dwflexiramp_shortfall[r.name, t] =
+                    @variable(model, lower_bound = 0)
+                if r.shortfall_penalty < 0
+                    set_upper_bound(upflexiramp_shortfall[r.name, t], 0.0)
+                    set_upper_bound(dwflexiramp_shortfall[r.name, t], 0.0)
+                end
+            end
+        end
     end
     return
 end
@@ -81,7 +104,7 @@ function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
     eq_startup_limit = _init(model, :eq_startup_limit)
     is_on = model[:is_on]
     prod_above = model[:prod_above]
-    reserve = model[:reserve]
+    reserve = _total_reserves(model, g)
     switch_off = model[:switch_off]
     switch_on = model[:switch_on]
     T = model[:instance].time
@@ -89,7 +112,7 @@ function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
         # Startup limit
         eq_startup_limit[g.name, t] = @constraint(
             model,
-            prod_above[g.name, t] + reserve[g.name, t] <=
+            prod_above[g.name, t] + reserve[t] <=
             (g.max_power[t] - g.min_power[t]) * is_on[g.name, t] -
             max(0, g.max_power[t] - g.startup_limit) * switch_on[g.name, t]
         )
@@ -117,7 +140,7 @@ function _add_ramp_eqs!(
     formulation::RampingFormulation,
 )::Nothing
     prod_above = model[:prod_above]
-    reserve = model[:reserve]
+    reserve = _total_reserves(model, g)
     eq_ramp_up = _init(model, :eq_ramp_up)
     eq_ramp_down = _init(model, :eq_ramp_down)
     for t in 1:model[:instance].time
@@ -126,14 +149,14 @@ function _add_ramp_eqs!(
             if _is_initially_on(g) == 1
                 eq_ramp_up[g.name, t] = @constraint(
                     model,
-                    prod_above[g.name, t] + reserve[g.name, t] <=
+                    prod_above[g.name, t] + reserve[t] <=
                     (g.initial_power - g.min_power[t]) + g.ramp_up_limit
                 )
             end
         else
             eq_ramp_up[g.name, t] = @constraint(
                 model,
-                prod_above[g.name, t] + reserve[g.name, t] <=
+                prod_above[g.name, t] + reserve[t] <=
                 prod_above[g.name, t-1] + g.ramp_up_limit
             )
         end
@@ -216,3 +239,15 @@ function _add_net_injection_eqs!(model::JuMP.Model, g::Unit)::Nothing
         )
     end
 end
+
+function _total_reserves(model, g)::Vector
+    T = model[:instance].time
+    reserve = [0.0 for _ in 1:T]
+    spinning_reserves = [r for r in g.reserves if r.type == "spinning"]
+    if !isempty(spinning_reserves)
+        reserve += [
+            sum(model[:reserve][r.name, g.name, t] for r in spinning_reserves) for t in 1:model[:instance].time
+        ]
+    end
+    return reserve
+end

src/solution/fix.jl

@@ -18,15 +18,28 @@ function fix!(model::JuMP.Model, solution::AbstractDict)::Nothing
             is_on_value = round(solution["Is on"][g.name][t])
             prod_value =
                 round(solution["Production (MW)"][g.name][t], digits = 5)
-            reserve_value =
-                round(solution["Reserve (MW)"][g.name][t], digits = 5)
             JuMP.fix(is_on[g.name, t], is_on_value, force = true)
             JuMP.fix(
                 prod_above[g.name, t],
                 prod_value - is_on_value * g.min_power[t],
                 force = true,
             )
-            JuMP.fix(reserve[g.name, t], reserve_value, force = true)
+        end
+    end
+    for r in instance.reserves
+        r.type == "spinning" || continue
+        for g in r.units
+            for t in 1:T
+                reserve_value = round(
+                    solution["Spinning reserve (MW)"][r.name][g.name][t],
+                    digits = 5,
+                )
+                JuMP.fix(
+                    reserve[r.name, g.name, t],
+                    reserve_value,
+                    force = true,
+                )
+            end
         end
     end
     return

src/solution/methods/XavQiuWanThi2019/optimize.jl

@@ -3,13 +3,12 @@
 # Released under the modified BSD license. See COPYING.md for more details.
 
 function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing
+    if !occursin("Gurobi", JuMP.solver_name(model))
+        method.two_phase_gap = false
+    end
     function set_gap(gap)
-        try
-            JuMP.set_optimizer_attribute(model, "MIPGap", gap)
-            @info @sprintf("MIP gap tolerance set to %f", gap)
-        catch
-            @warn "Could not change MIP gap tolerance"
-        end
+        JuMP.set_optimizer_attribute(model, "MIPGap", gap)
+        @info @sprintf("MIP gap tolerance set to %f", gap)
     end
     initial_time = time()
     large_gap = false
@@ -17,8 +16,6 @@ function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing
     if has_transmission && method.two_phase_gap
         set_gap(1e-2)
         large_gap = true
-    else
-        set_gap(method.gap_limit)
     end
     while true
         time_elapsed = time() - initial_time

src/solution/methods/XavQiuWanThi2019/structs.jl

@@ -13,7 +13,7 @@ Lazy constraint solution method described in:
 module XavQiuWanThi2019
 import ..SolutionMethod
 """
-    struct Method
+    mutable struct Method
         time_limit::Float64
         gap_limit::Float64
         two_phase_gap::Bool
@@ -27,7 +27,7 @@ Fields
 - `time_limit`:
     the time limit over the entire optimization procedure.
 - `gap_limit`: 
-    the desired relative optimality gap.
+    the desired relative optimality gap. Only used when `two_phase_gap=true`.
 - `two_phase_gap`: 
     if true, solve the problem with large gap tolerance first, then reduce
     the gap tolerance when no further violated constraints are found.
@@ -39,7 +39,7 @@ Fields
     formulation per time period.
 
 """
-struct Method <: SolutionMethod
+mutable struct Method <: SolutionMethod
     time_limit::Float64
     gap_limit::Float64
     two_phase_gap::Bool

src/solution/solution.jl

@@ -50,13 +50,6 @@ function solution(model::JuMP.Model)::OrderedDict
     sol["Is on"] = timeseries(model[:is_on], instance.units)
     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)
@@ -67,5 +60,47 @@ function solution(model::JuMP.Model)::OrderedDict
         sol["Price-sensitive loads (MW)"] =
             timeseries(model[:loads], instance.price_sensitive_loads)
     end
+    sol["Spinning reserve (MW)"] = OrderedDict(
+        r.name => OrderedDict(
+            g.name => [
+                value(model[:reserve][r.name, g.name, t]) for
+                t in 1:instance.time
+            ] for g in r.units
+        ) for r in instance.reserves if r.type == "spinning"
+    )
+    sol["Spinning reserve shortfall (MW)"] = OrderedDict(
+        r.name => [
+            value(model[:reserve_shortfall][r.name, t]) for
+            t in 1:instance.time
+        ] for r in instance.reserves if r.type == "spinning"
+    )
+    sol["Up-flexiramp (MW)"] = OrderedDict(
+        r.name => OrderedDict(
+            g.name => [
+                value(model[:upflexiramp][r.name, g.name, t]) for
+                t in 1:instance.time
+            ] for g in r.units
+        ) for r in instance.reserves if r.type == "flexiramp"
+    )
+    sol["Up-flexiramp shortfall (MW)"] = OrderedDict(
+        r.name => [
+            value(model[:upflexiramp_shortfall][r.name, t]) for
+            t in 1:instance.time
+        ] for r in instance.reserves if r.type == "flexiramp"
+    )
+    sol["Down-flexiramp (MW)"] = OrderedDict(
+        r.name => OrderedDict(
+            g.name => [
+                value(model[:dwflexiramp][r.name, g.name, t]) for
+                t in 1:instance.time
+            ] for g in r.units
+        ) for r in instance.reserves if r.type == "flexiramp"
+    )
+    sol["Down-flexiramp shortfall (MW)"] = OrderedDict(
+        r.name => [
+            value(model[:upflexiramp_shortfall][r.name, t]) for
+            t in 1:instance.time
+        ] for r in instance.reserves if r.type == "flexiramp"
+    )
     return sol
 end

src/transform/randomize/XavQiuAhm2021.jl

@@ -118,11 +118,12 @@ Base.@kwdef struct Randomization
 end
 
 function _randomize_costs(
+    rng,
     instance::UnitCommitmentInstance,
     distribution,
 )::Nothing
     for unit in instance.units
-        α = rand(distribution)
+        α = rand(rng, distribution)
         unit.min_power_cost *= α
         for k in unit.cost_segments
             k.cost *= α
@@ -135,10 +136,11 @@ function _randomize_costs(
 end
 
 function _randomize_load_share(
+    rng,
     instance::UnitCommitmentInstance,
     distribution,
 )::Nothing
-    α = rand(distribution, length(instance.buses))
+    α = rand(rng, distribution, length(instance.buses))
     for t in 1:instance.time
         total = sum(bus.load[t] for bus in instance.buses)
         den = sum(
@@ -153,6 +155,7 @@ function _randomize_load_share(
 end
 
 function _randomize_load_profile(
+    rng,
     instance::UnitCommitmentInstance,
     params::Randomization,
 )::Nothing
@@ -161,12 +164,13 @@ function _randomize_load_profile(
     for t in 2:instance.time
         idx = (t - 1) % length(params.load_profile_mu) + 1
         gamma = rand(
+            rng,
             Normal(params.load_profile_mu[idx], params.load_profile_sigma[idx]),
         )
         push!(system_load, system_load[t-1] * gamma)
     end
     capacity = sum(maximum(u.max_power) for u in instance.units)
-    peak_load = rand(params.peak_load) * capacity
+    peak_load = rand(rng, params.peak_load) * capacity
     system_load = system_load ./ maximum(system_load) .* peak_load
 
     # Scale bus loads to match the new system load
@@ -186,22 +190,24 @@ end
     function randomize!(
         instance::UnitCommitment.UnitCommitmentInstance,
         method::XavQiuAhm2021.Randomization,
+        rng = MersenneTwister(),
     )::Nothing
 
 Randomize costs and loads based on the method described in XavQiuAhm2021.
 """
 function randomize!(
     instance::UnitCommitment.UnitCommitmentInstance,
-    method::XavQiuAhm2021.Randomization,
+    method::XavQiuAhm2021.Randomization;
+    rng = MersenneTwister(),
 )::Nothing
     if method.randomize_costs
-        XavQiuAhm2021._randomize_costs(instance, method.cost)
+        XavQiuAhm2021._randomize_costs(rng, instance, method.cost)
     end
     if method.randomize_load_share
-        XavQiuAhm2021._randomize_load_share(instance, method.load_share)
+        XavQiuAhm2021._randomize_load_share(rng, instance, method.load_share)
     end
     if method.randomize_load_profile
-        XavQiuAhm2021._randomize_load_profile(instance, method)
+        XavQiuAhm2021._randomize_load_profile(rng, instance, method)
     end
     return
 end

src/transform/slice.jl

@@ -24,13 +24,14 @@ function slice(
     modified = deepcopy(instance)
     modified.time = length(range)
     modified.power_balance_penalty = modified.power_balance_penalty[range]
-    modified.reserves.spinning = modified.reserves.spinning[range]
+    for r in modified.reserves
+        r.amount = r.amount[range]
+    end
     for u in modified.units
         u.max_power = u.max_power[range]
         u.min_power = u.min_power[range]
         u.must_run = u.must_run[range]
         u.min_power_cost = u.min_power_cost[range]
-        u.provides_spinning_reserves = u.provides_spinning_reserves[range]
         for s in u.cost_segments
             s.mw = s.mw[range]
             s.cost = s.cost[range]

src/utils/log.jl

@@ -5,20 +5,11 @@
 import Logging: min_enabled_level, shouldlog, handle_message
 using Base.CoreLogging, Logging, Printf
 
-struct TimeLogger <: AbstractLogger
+Base.@kwdef struct TimeLogger <: AbstractLogger
     initial_time::Float64
-    file::Union{Nothing,IOStream}
-    screen_log_level::Any
-    io_log_level::Any
-end
-
-function TimeLogger(;
-    initial_time::Float64,
-    file::Union{Nothing,IOStream} = nothing,
-    screen_log_level = CoreLogging.Info,
-    io_log_level = CoreLogging.Info,
-)::TimeLogger
-    return TimeLogger(initial_time, file, screen_log_level, io_log_level)
+    file::Union{Nothing,IOStream} = nothing
+    screen_log_level::Any = CoreLogging.Info
+    io_log_level::Any = CoreLogging.Info
 end
 
 min_enabled_level(logger::TimeLogger) = logger.io_log_level
@@ -61,7 +52,9 @@ function handle_message(
     end
 end
 
-function _setup_logger()
+function _setup_logger(; level = CoreLogging.Info)
     initial_time = time()
-    return global_logger(TimeLogger(initial_time = initial_time))
+    return global_logger(
+        TimeLogger(initial_time = initial_time, screen_log_level = level),
+    )
 end

src/validation/validate.jl

@@ -40,12 +40,19 @@ function validate(
     return true
 end
 
-function _validate_units(instance, solution; tol = 0.01)
+function _validate_units(instance::UnitCommitmentInstance, solution; tol = 0.01)
     err_count = 0
 
     for unit in instance.units
         production = solution["Production (MW)"][unit.name]
-        reserve = solution["Reserve (MW)"][unit.name]
+        reserve = [0.0 for _ in 1:instance.time]
+        spinning_reserves = [r for r in unit.reserves if r.type == "spinning"]
+        if !isempty(spinning_reserves)
+            reserve += sum(
+                solution["Spinning reserve (MW)"][r.name][unit.name] for
+                r in spinning_reserves
+            )
+        end
         actual_production_cost = solution["Production cost (\$)"][unit.name]
         actual_startup_cost = solution["Startup cost (\$)"][unit.name]
         is_on = bin(solution["Is on"][unit.name])
@@ -99,13 +106,18 @@ function _validate_units(instance, solution; tol = 0.01)
             end
 
             # Verify reserve eligibility
-            if !unit.provides_spinning_reserves[t] && reserve[t] > tol
-                @error @sprintf(
-                    "Unit %s is not eligible to provide spinning reserves at time %d",
-                    unit.name,
-                    t
-                )
-                err_count += 1
+            for r in instance.reserves
+                if r.type == "spinning"
+                    if unit ∉ r.units &&
+                       (unit in keys(solution["Spinning reserve (MW)"][r.name]))
+                        @error @sprintf(
+                            "Unit %s is not eligible to provide reserve %s",
+                            unit.name,
+                            r.name,
+                        )
+                        err_count += 1
+                    end
+                end
             end
 
             # If unit is on, must produce at least its minimum power
@@ -137,9 +149,11 @@ function _validate_units(instance, solution; tol = 0.01)
             # If unit is off, must produce zero
             if !is_on[t] && production[t] + reserve[t] > tol
                 @error @sprintf(
-                    "Unit %s produces power at time %d while off",
+                    "Unit %s produces power at time %d while off (%.2f + %.2f > 0)",
                     unit.name,
-                    t
+                    t,
+                    production[t],
+                    reserve[t],
                 )
                 err_count += 1
             end
@@ -321,22 +335,66 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01)
             err_count += 1
         end
 
-        # Verify spinning reserves
-        reserve =
-            sum(solution["Reserve (MW)"][g.name][t] for g in instance.units)
-        reserve_shortfall =
-            (instance.shortfall_penalty[t] >= 0) ?
-            solution["Reserve shortfall (MW)"][t] : 0
+        # Verify reserves
+        for r in instance.reserves
+            if r.type == "spinning"
+                provided = sum(
+                    solution["Spinning reserve (MW)"][r.name][g.name][t] for
+                    g in r.units
+                )
+                shortfall =
+                    solution["Spinning reserve shortfall (MW)"][r.name][t]
+                required = r.amount[t]
 
-        if reserve + reserve_shortfall < instance.reserves.spinning[t] - tol
-            @error @sprintf(
-                "Insufficient spinning reserves at time %d (%.2f + %.2f should be %.2f)",
-                t,
-                reserve,
-                reserve_shortfall,
-                instance.reserves.spinning[t],
-            )
-            err_count += 1
+                if provided + shortfall < required - tol
+                    @error @sprintf(
+                        "Insufficient reserve %s at time %d (%.2f + %.2f < %.2f)",
+                        r.name,
+                        t,
+                        provided,
+                        shortfall,
+                        required,
+                    )
+                end
+            elseif r.type == "flexiramp"
+                upflexiramp = sum(
+                    solution["Up-flexiramp (MW)"][r.name][g.name][t] for
+                    g in r.units
+                )
+                upflexiramp_shortfall =
+                    solution["Up-flexiramp shortfall (MW)"][r.name][t]
+
+                if upflexiramp + upflexiramp_shortfall < r.amount[t] - tol
+                    @error @sprintf(
+                        "Insufficient up-flexiramp at time %d (%.2f + %.2f < %.2f)",
+                        t,
+                        upflexiramp,
+                        upflexiramp_shortfall,
+                        r.amount[t],
+                    )
+                    err_count += 1
+                end
+
+                dwflexiramp = sum(
+                    solution["Down-flexiramp (MW)"][r.name][g.name][t] for
+                    g in r.units
+                )
+                dwflexiramp_shortfall =
+                    solution["Down-flexiramp shortfall (MW)"][r.name][t]
+
+                if dwflexiramp + dwflexiramp_shortfall < r.amount[t] - tol
+                    @error @sprintf(
+                        "Insufficient down-flexiramp at time %d (%.2f + %.2f < %.2f)",
+                        t,
+                        dwflexiramp,
+                        dwflexiramp_shortfall,
+                        r.amount[t],
+                    )
+                    err_count += 1
+                end
+            else
+                error("Unknown reserve type: $(r.type)")
+            end
         end
     end
 

test/Project.toml

@@ -3,7 +3,6 @@ Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63"
-Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -20,7 +19,7 @@ DataStructures = "0.18"
 Distributions = "0.25"
 GZip = "0.5"
 JSON = "0.21"
-JuMP = "0.21"
-MathOptInterface = "0.9"
+JuMP = "1"
+MathOptInterface = "1"
 PackageCompiler = "1"
 julia = "1"

test/instance/read_test.jl

@@ -37,6 +37,11 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
     @test instance.buses[9].load == [35.36638, 33.25495, 31.67138, 31.14353]
     @test instance.buses_by_name["b9"].name == "b9"
 
+    @test instance.reserves[1].name == "r1"
+    @test instance.reserves[1].type == "spinning"
+    @test instance.reserves[1].amount == [100.0, 100.0, 100.0, 100.0]
+    @test instance.reserves_by_name["r1"].name == "r1"
+
     unit = instance.units[1]
     @test unit.name == "g1"
     @test unit.bus.name == "b1"
@@ -48,7 +53,6 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
     @test unit.min_power_cost == [1400.0 for t in 1:4]
     @test unit.min_uptime == 1
     @test unit.min_downtime == 1
-    @test unit.provides_spinning_reserves == [true for t in 1:4]
     for t in 1:1
         @test unit.cost_segments[1].mw[t] == 10.0
         @test unit.cost_segments[2].mw[t] == 20.0
@@ -64,11 +68,13 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
     @test unit.startup_categories[1].cost == 1000.0
     @test unit.startup_categories[2].cost == 1500.0
     @test unit.startup_categories[3].cost == 2000.0
+    @test length(unit.reserves) == 0
     @test instance.units_by_name["g1"].name == "g1"
 
     unit = instance.units[2]
     @test unit.name == "g2"
     @test unit.must_run == [false for t in 1:4]
+    @test length(unit.reserves) == 1
 
     unit = instance.units[3]
     @test unit.name == "g3"
@@ -81,7 +87,6 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
     @test unit.min_power_cost == [0.0 for t in 1:4]
     @test unit.min_uptime == 1
     @test unit.min_downtime == 1
-    @test unit.provides_spinning_reserves == [true for t in 1:4]
     for t in 1:4
         @test unit.cost_segments[1].mw[t] ≈ 33
         @test unit.cost_segments[2].mw[t] ≈ 33
@@ -90,8 +95,8 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
         @test unit.cost_segments[2].cost[t] ≈ 38.04
         @test unit.cost_segments[3].cost[t] ≈ 44.77853
     end
-
-    @test instance.reserves.spinning == zeros(4)
+    @test length(unit.reserves) == 1
+    @test unit.reserves[1].name == "r1"
 
     @test instance.contingencies[1].lines == [instance.lines[1]]
     @test instance.contingencies[1].units == []

test/model/formulations_test.jl

@@ -4,6 +4,8 @@
 
 using UnitCommitment
 using JuMP
+using Cbc
+using JSON
 import UnitCommitment:
     ArrCon2000,
     CarArr2006,
@@ -13,62 +15,70 @@ import UnitCommitment:
     KnuOstWat2018,
     MorLatRam2013,
     PanGua2016,
-    XavQiuWanThi2019
+    XavQiuWanThi2019,
+    WanHob2016
 
-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)
+function _test(
+    formulation::Formulation;
+    instances = ["test/case14"],
+    dump::Bool = false,
+)::Nothing
+    for instance_name in instances
+        instance = UnitCommitment.read_benchmark(instance_name)
         model = UnitCommitment.build_model(
             instance = instance,
             formulation = formulation,
-            optimizer = Gurobi.Optimizer,
-        )
-        UnitCommitment.optimize!(
-            model,
-            XavQiuWanThi2019.Method(two_phase_gap = false, gap_limit = 0.1),
+            optimizer = Cbc.Optimizer,
+            variable_names = true,
         )
+        set_silent(model)
+        UnitCommitment.optimize!(model)
         solution = UnitCommitment.solution(model)
+        if dump
+            open("/tmp/ucjl.json", "w") do f
+                return write(f, JSON.json(solution, 2))
+            end
+            write_to_file(model, "/tmp/ucjl.lp")
+        end
         @test UnitCommitment.validate(instance, solution)
     end
     return
 end
 
-function _test(formulation::Formulation)::Nothing
-    _small_test(formulation)
-    if ENABLE_LARGE_TESTS
-        _large_test(formulation)
+@testset "formulations" begin
+    @testset "default" begin
+        _test(Formulation())
+    end
+    @testset "ArrCon2000" begin
+        _test(Formulation(ramping = ArrCon2000.Ramping()))
+    end
+    @testset "DamKucRajAta2016" begin
+        _test(Formulation(ramping = DamKucRajAta2016.Ramping()))
+    end
+    @testset "MorLatRam2013" begin
+        _test(
+            Formulation(
+                ramping = MorLatRam2013.Ramping(),
+                startup_costs = MorLatRam2013.StartupCosts(),
+            ),
+        )
+    end
+    @testset "PanGua2016" begin
+        _test(Formulation(ramping = PanGua2016.Ramping()))
+    end
+    @testset "Gar1962" begin
+        _test(Formulation(pwl_costs = Gar1962.PwlCosts()))
+    end
+    @testset "CarArr2006" begin
+        _test(Formulation(pwl_costs = CarArr2006.PwlCosts()))
+    end
+    @testset "KnuOstWat2018" begin
+        _test(Formulation(pwl_costs = KnuOstWat2018.PwlCosts()))
+    end
+    @testset "WanHob2016" begin
+        _test(
+            Formulation(ramping = WanHob2016.Ramping()),
+            instances = ["test/case14-flex"],
+        )
     end
 end
-
-@testset "formulations" begin
-    _test(Formulation())
-    _test(Formulation(ramping = ArrCon2000.Ramping()))
-    # _test(Formulation(ramping = DamKucRajAta2016.Ramping()))
-    _test(
-        Formulation(
-            ramping = MorLatRam2013.Ramping(),
-            startup_costs = MorLatRam2013.StartupCosts(),
-        ),
-    )
-    _test(Formulation(ramping = PanGua2016.Ramping()))
-    _test(Formulation(pwl_costs = Gar1962.PwlCosts()))
-    _test(Formulation(pwl_costs = CarArr2006.PwlCosts()))
-    _test(Formulation(pwl_costs = KnuOstWat2018.PwlCosts()))
-end

test/runtests.jl

@@ -6,9 +6,7 @@ using Test
 using UnitCommitment
 
 push!(Base.LOAD_PATH, @__DIR__)
-UnitCommitment._setup_logger()
-
-const ENABLE_LARGE_TESTS = ("UCJL_LARGE_TESTS" in keys(ENV))
+UnitCommitment._setup_logger(level = Base.CoreLogging.Error)
 
 @testset "UnitCommitment" begin
     include("usage.jl")
@@ -21,10 +19,12 @@ const ENABLE_LARGE_TESTS = ("UCJL_LARGE_TESTS" in keys(ENV))
     @testset "model" begin
         include("model/formulations_test.jl")
     end
-    @testset "XavQiuWanThi19" begin
-        include("solution/methods/XavQiuWanThi19/filter_test.jl")
-        include("solution/methods/XavQiuWanThi19/find_test.jl")
-        include("solution/methods/XavQiuWanThi19/sensitivity_test.jl")
+    @testset "solution" begin
+        @testset "XavQiuWanThi19" begin
+            include("solution/methods/XavQiuWanThi19/filter_test.jl")
+            include("solution/methods/XavQiuWanThi19/find_test.jl")
+            include("solution/methods/XavQiuWanThi19/sensitivity_test.jl")
+        end
     end
     @testset "transform" begin
         include("transform/initcond_test.jl")

test/transform/randomize/XavQiuAhm2021_test.jl

@@ -6,6 +6,7 @@ import Random
 import UnitCommitment: XavQiuAhm2021
 
 using Distributions
+using Random
 using UnitCommitment, Cbc, JuMP
 
 get_instance() = UnitCommitment.read_benchmark("matpower/case118/2017-02-01")
@@ -27,10 +28,10 @@ test_approx(x, y) = @test isapprox(x, y, atol = 1e-3)
         prev_system_load = system_load(instance)
         test_approx(bus.load[1] / prev_system_load[1], 0.012)
 
-        Random.seed!(42)
         randomize!(
             instance,
             XavQiuAhm2021.Randomization(randomize_load_profile = false),
+            rng = MersenneTwister(42),
         )
 
         # Check randomized costs
@@ -53,8 +54,11 @@ test_approx(x, y) = @test isapprox(x, y, atol = 1e-3)
         @test round.(system_load(instance), digits = 1)[1:8] ≈
               [3059.5, 2983.2, 2937.5, 2953.9, 3073.1, 3356.4, 4068.5, 4018.8]
 
-        Random.seed!(42)
-        randomize!(instance, XavQiuAhm2021.Randomization())
+        randomize!(
+            instance,
+            XavQiuAhm2021.Randomization(),
+            rng = MersenneTwister(42),
+        )
 
         # Check randomized load profile
         @test round.(system_load(instance), digits = 1)[1:8] ≈

test/transform/slice_test.jl

@@ -11,13 +11,12 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
     # Should update all time-dependent fields
     @test modified.time == 2
     @test length(modified.power_balance_penalty) == 2
-    @test length(modified.reserves.spinning) == 2
+    @test length(modified.reserves_by_name["r1"].amount) == 2
     for u in modified.units
         @test length(u.max_power) == 2
         @test length(u.min_power) == 2
         @test length(u.must_run) == 2
         @test length(u.min_power_cost) == 2
-        @test length(u.provides_spinning_reserves) == 2
         for s in u.cost_segments
             @test length(s.mw) == 2
             @test length(s.cost) == 2
@@ -35,7 +34,6 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
         @test length(ps.demand) == 2
         @test length(ps.revenue) == 2
     end
-
     # Should be able to build model without errors
     optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
     model = UnitCommitment.build_model(

test/usage.jl

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