Merge branch 'dev' into feature/reserves

.github/workflows/test.yml

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

.gitignore

@@ -1,21 +1,38 @@
 *.bak
 *.gz
+*.ipynb
 *.lastrun
-*.so
 *.mps
-*.ipynb
+*.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

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

@@ -25,7 +25,7 @@ DataStructures = "0.18"
 Distributions = "0.25"
 GZip = "0.5"
 JSON = "0.21"
-JuMP = "0.21"
+JuMP = "1"
-MathOptInterface = "0.9"
+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

@@ -209,17 +209,17 @@ This section describes the hourly amount of reserves required.
 
 | Key                   | Description                                        | Default   |  Time series?
 | :-------------------- | :------------------------------------------------- | --------- |  :----:
-| `Type` | Type of reserve product. Currently, only `Spinning` is supported. | Required | N
+| `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": {
         "r1": {
-            "Type": "Spinning",
+            "Type": "spinning",
             "Amount (MW)": [
                 57.30552,
                 53.88429,
@@ -227,6 +227,15 @@ This section describes the hourly amount of reserves required.
                 50.46307
             ],
             "Shortfall penalty ($/MW)": 5.0
+        },
+        "r2": {
+            "Type": "flexiramp",
+            "Amount (MW)": [
+                20.31042,
+                23.65273,
+                27.41784,
+                25.34057
+            ],
         }
     }
 }
@@ -293,3 +302,4 @@ Current limitations
 * 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. 

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

@@ -119,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],
@@ -317,6 +322,7 @@ function _from_json(json; repair = true)
         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

@@ -86,6 +86,7 @@ Base.@kwdef mutable struct UnitCommitmentInstance
     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/WanHob2016/ramp.jl

@@ -0,0 +1,168 @@
+# 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,16 +28,17 @@ function _add_net_injection_eqs!(model::JuMP.Model)::Nothing
     return
 end
 
-function _add_reserve_eqs!(model::JuMP.Model)::Nothing
+function _add_spinning_reserve_eqs!(model::JuMP.Model)::Nothing
     instance = model[:instance]
-    eq_min_reserve = _init(model, :eq_min_reserve)
+    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_reserve[r.name, t] = @constraint(
+            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]
@@ -54,3 +56,44 @@ function _add_reserve_eqs!(model::JuMP.Model)::Nothing
     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
+end

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,10 +43,11 @@ 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 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)
@@ -59,6 +61,32 @@ function _add_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
     return
 end
 
+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
+
 function _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing
     startup = _init(model, :startup)
     for t in 1:model[:instance].time
@@ -176,16 +204,14 @@ function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing
                 eq_min_uptime[g.name, 0] = @constraint(
                     model,
                     sum(
-                        switch_off[g.name, i] for
+                        switch_off[g.name, i] for i in 1:(g.min_uptime-g.initial_status) if i <= T
-                        i in 1:(g.min_uptime-g.initial_status) if i <= T
                     ) == 0
                 )
             else
                 eq_min_downtime[g.name, 0] = @constraint(
                     model,
                     sum(
-                        switch_on[g.name, i] for
+                        switch_on[g.name, i] for i in 1:(g.min_downtime+g.initial_status) if i <= T
-                        i in 1:(g.min_downtime+g.initial_status) if i <= T
                     ) == 0
                 )
             end
@@ -213,10 +239,10 @@ end
 function _total_reserves(model, g)::Vector
     T = model[:instance].time
     reserve = [0.0 for _ in 1:T]
-    if !isempty(g.reserves)
+    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 g.reserves) for
+            sum(model[:reserve][r.name, g.name, t] for r in spinning_reserves) for t in 1:model[:instance].time
-            t in 1:model[:instance].time
         ]
     end
     return reserve

src/solution/fix.jl

@@ -27,10 +27,11 @@ function fix!(model::JuMP.Model, solution::AbstractDict)::Nothing
         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["Reserve (MW)"][r.name][g.name][t],
+                    solution["Spinning reserve (MW)"][r.name][g.name][t],
                     digits = 5,
                 )
                 JuMP.fix(

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)
-            JuMP.set_optimizer_attribute(model, "MIPGap", gap)
+        @info @sprintf("MIP gap tolerance set to %f", gap)
-            @info @sprintf("MIP gap tolerance set to %f", gap)
-        catch
-            @warn "Could not change MIP gap tolerance"
-        end
     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

@@ -15,8 +15,7 @@ function solution(model::JuMP.Model)::OrderedDict
             value(model[:is_on][g.name, t]) * g.min_power_cost[t] + sum(
                 Float64[
                     value(model[:segprod][g.name, t, k]) *
-                    g.cost_segments[k].cost[t] for
+                    g.cost_segments[k].cost[t] for k in 1:length(g.cost_segments)
-                    k in 1:length(g.cost_segments)
                 ],
             ) for t in 1:T
         ]
@@ -25,8 +24,7 @@ function solution(model::JuMP.Model)::OrderedDict
         return [
             value(model[:is_on][g.name, t]) * g.min_power[t] + sum(
                 Float64[
-                    value(model[:segprod][g.name, t, k]) for
+                    value(model[:segprod][g.name, t, k]) for k in 1:length(g.cost_segments)
-                    k in 1:length(g.cost_segments)
                 ],
             ) for t in 1:T
         ]
@@ -60,19 +58,45 @@ function solution(model::JuMP.Model)::OrderedDict
         sol["Price-sensitive loads (MW)"] =
             timeseries(model[:loads], instance.price_sensitive_loads)
     end
-    sol["Reserve (MW)"] = OrderedDict(
+    sol["Spinning reserve (MW)"] = OrderedDict(
         r.name => OrderedDict(
             g.name => [
-                value(model[:reserve][r.name, g.name, t]) for
+                value(model[:reserve][r.name, g.name, t]) for t in 1:instance.time
-                t in 1:instance.time
             ] for g in r.units
-        ) for r in instance.reserves
+        ) for r in instance.reserves if r.type == "spinning"
     )
-    sol["Reserve shortfall (MW)"] = OrderedDict(
+    sol["Spinning reserve shortfall (MW)"] = OrderedDict(
         r.name => [
-            value(model[:reserve_shortfall][r.name, t]) for
+            value(model[:reserve_shortfall][r.name, t]) for t in 1:instance.time
-            t in 1:instance.time
+        ] for r in instance.reserves if r.type == "spinning"
-        ] for r in instance.reserves
+    )
+    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/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}
+    file::Union{Nothing,IOStream} = nothing
-    screen_log_level::Any
+    screen_log_level::Any = CoreLogging.Info
-    io_log_level::Any
+    io_log_level::Any = CoreLogging.Info
-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)
 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

@@ -46,10 +46,11 @@ function _validate_units(instance::UnitCommitmentInstance, solution; tol = 0.01)
     for unit in instance.units
         production = solution["Production (MW)"][unit.name]
         reserve = [0.0 for _ in 1:instance.time]
-        if !isempty(unit.reserves)
+        spinning_reserves = [r for r in unit.reserves if r.type == "spinning"]
+        if !isempty(spinning_reserves)
             reserve += sum(
-                solution["Reserve (MW)"][r.name][unit.name] for
+                solution["Spinning reserve (MW)"][r.name][unit.name]
-                r in unit.reserves
+                for r in spinning_reserves
             )
         end
         actual_production_cost = solution["Production cost (\$)"][unit.name]
@@ -106,14 +107,16 @@ function _validate_units(instance::UnitCommitmentInstance, solution; tol = 0.01)
 
             # Verify reserve eligibility
             for r in instance.reserves
-                if unit ∉ r.units &&
+                if r.type == "spinning"
-                   (unit in keys(solution["Reserve (MW)"][r.name]))
+                    if unit ∉ r.units &&
-                    @error @sprintf(
+                    (unit in keys(solution["Spinning reserve (MW)"][r.name]))
-                        "Unit %s is not eligible to provide reserve %s",
+                        @error @sprintf(
-                        unit.name,
+                            "Unit %s is not eligible to provide reserve %s",
-                        r.name,
+                            unit.name,
-                    )
+                            r.name,
-                    err_count += 1
+                        )
+                        err_count += 1
+                    end
                 end
             end
 
@@ -305,16 +308,14 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01)
         ps_load = 0
         if length(instance.price_sensitive_loads) > 0
             ps_load = sum(
-                solution["Price-sensitive loads (MW)"][ps.name][t] for
+                solution["Price-sensitive loads (MW)"][ps.name][t] for ps in instance.price_sensitive_loads
-                ps in instance.price_sensitive_loads
             )
         end
         production =
             sum(solution["Production (MW)"][g.name][t] for g in instance.units)
         if "Load curtail (MW)" in keys(solution)
             load_curtail = sum(
-                solution["Load curtail (MW)"][b.name][t] for
+                solution["Load curtail (MW)"][b.name][t] for b in instance.buses
-                b in instance.buses
             )
         end
         balance = fixed_load - load_curtail - production + ps_load
@@ -334,22 +335,58 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01)
 
         # Verify reserves
         for r in instance.reserves
-            provided = sum(
+            if r.type == "spinning"
-                solution["Reserve (MW)"][r.name][g.name][t] for g in r.units
+                provided = sum(
-            )
+                    solution["Spinning reserve (MW)"][r.name][g.name][t] for g in r.units
-            shortfall = solution["Reserve shortfall (MW)"][r.name][t]
+                )
-            required = r.amount[t]
+                shortfall = solution["Spinning reserve shortfall (MW)"][r.name][t]
+                required = r.amount[t]
 
-            if provided + shortfall < required - tol
+                if provided + shortfall < required - tol
-                @error @sprintf(
+                    @error @sprintf(
-                    "Insufficient reserve %s at time %d (%.2f + %.2f < %.2f)",
+                        "Insufficient reserve %s at time %d (%.2f + %.2f < %.2f)",
-                    r.name,
+                        r.name,
-                    t,
+                        t,
-                    provided,
+                        provided,
-                    shortfall,
+                        shortfall,
-                    required,
+                        required,
+                    )
+                end
+            elseif r.type == "flexiramp"
+                upflexiramp = sum(
+                    solution["Up-flexiramp (MW)"][r.name][g.name][t] for g in r.units
                 )
-                err_count += 1
+                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"
+JuMP = "1"
-MathOptInterface = "0.9"
+MathOptInterface = "1"
 PackageCompiler = "1"
 julia = "1"

test/model/formulations_test.jl

@@ -15,56 +15,36 @@ import UnitCommitment:
     KnuOstWat2018,
     MorLatRam2013,
     PanGua2016,
-    XavQiuWanThi2019
+    XavQiuWanThi2019,
+    WanHob2016
 
-if ENABLE_LARGE_TESTS
+function _test(
-    using Gurobi
+    formulation::Formulation;
-end
+    instances=["test/case14"],
-
+    dump::Bool = false,
-function _small_test(formulation::Formulation; dump::Bool = false)::Nothing
+)::Nothing
-    instance = UnitCommitment.read_benchmark("test/case14")
+    for instance_name in instances
-    model = UnitCommitment.build_model(
+        instance = UnitCommitment.read_benchmark(instance_name)
-        instance = instance,
+        model = UnitCommitment.build_model(
-        formulation = formulation,
+            instance = instance,
-        optimizer = Cbc.Optimizer,
+            formulation = formulation,
-        variable_names = true,
+            optimizer = Cbc.Optimizer,
-    )
+            variable_names = true,
-    UnitCommitment.optimize!(model)
+        )
-    solution = UnitCommitment.solution(model)
+        set_silent(model)
-    if dump
+        UnitCommitment.optimize!(model)
-        open("/tmp/ucjl.json", "w") do f
+        solution = UnitCommitment.solution(model)
-            return write(f, JSON.json(solution, 2))
+        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
-        write_to_file(model, "/tmp/ucjl.lp")
+        @test UnitCommitment.validate(instance, solution)
     end
-    @test UnitCommitment.validate(instance, solution)
-    return
-end
-
-function _large_test(formulation::Formulation)::Nothing
-    instance =
-        UnitCommitment.read_benchmark("pglib-uc/ca/Scenario400_reserves_1")
-    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)
     return
 end
 
-function _test(formulation::Formulation; dump::Bool = false)::Nothing
-    _small_test(formulation; dump)
-    if ENABLE_LARGE_TESTS
-        _large_test(formulation)
-    end
-end
-
 @testset "formulations" begin
     @testset "default" begin
         _test(Formulation())
@@ -95,4 +75,10 @@ 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

test/runtests.jl

@@ -6,9 +6,7 @@ using Test
 using UnitCommitment
 
 push!(Base.LOAD_PATH, @__DIR__)
-UnitCommitment._setup_logger()
+UnitCommitment._setup_logger(level = Base.CoreLogging.Error)
-
-const ENABLE_LARGE_TESTS = ("UCJL_LARGE_TESTS" in keys(ENV))
 
 @testset "UnitCommitment" begin
     include("usage.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!(
-        randomize!(instance, XavQiuAhm2021.Randomization())
+            instance,
+            XavQiuAhm2021.Randomization(),
+            rng = MersenneTwister(42),
+        )
 
         # Check randomized load profile
         @test round.(system_load(instance), digits = 1)[1:8] ≈

test/transform/slice_test.jl

@@ -34,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(