GitHub Actions: Remove older non-LTS Julia versions

.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
-*.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

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

@@ -24,7 +24,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

@@ -36,6 +36,7 @@ This section describes system-wide parameters, such as power balance and reserve
 | `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
+| `Flexiramp penalty ($/MW)` | Penalty for system-wide shortage in meeting flexible ramping product requirements (in $/MW). This is charged per time step. | `500` | Y
 
 
 #### Example
@@ -44,7 +45,8 @@ 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
+        "Reserve shortfall penalty ($/MW)": -1.0,
+        "Flexiramp penalty ($/MW)": 100.0
     }
 }
 ```
@@ -97,6 +99,8 @@ This section describes all generators in the system, including thermal units, re
 | `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
+| `Provides flexible capacity?`    | If `true`, this generator may provide flexible ramping product (Boolean). | `true` | Y
+
 
 #### Production costs and limits
 
@@ -136,6 +140,7 @@ Note that this curve also specifies the production limits. Specifically, the fir
             "Initial status (h)": 12,
             "Must run?": false,
             "Provides spinning reserves?": true,
+            "Provides flexible capacity?": false,
         },
         "gen2": {
             "Bus": "b5",
@@ -206,14 +211,16 @@ 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
+| `Up-flexiramp (MW)`       | Minimum amount of system-wide upward flexible ramping product (in MW). Only generators which are online may provide this reserve. | `0.0` | Y
+| `Down-flexiramp (MW)`       | Minimum amount of system-wide downward flexible ramping product (in MW). Only generators which are online may provide this reserve. | `0.0` | Y
 
-#### Example
+#### Example 1
 
 ```json
 {
@@ -228,6 +235,27 @@ This section describes the hourly amount of operating reserves required.
 }
 ```
 
+#### Example 2
+
+```json
+{
+    "Reserves": {
+        "Up-flexiramp (MW)": [
+            20.31042,
+            23.65273,
+            27.41784,
+            25.34057
+        ],
+         "Down-flexiramp (MW)": [
+            19.41546,
+            21.45377,
+            23.53402,
+            24.80973
+        ]
+    }
+}
+```
+
 ### Contingencies
 
 This section describes credible contingency scenarios in the optimization, such as the loss of a transmission line or generator.
@@ -287,6 +315,7 @@ Current limitations
 -------------------
 
 * All reserves are system-wide. Zonal reserves are not currently supported.
+* Upward and downward flexible ramping products can only be acquired under the WanHob2016 formulation, which does not support spinning reserves. 
 * 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.

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

@@ -117,6 +117,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],
@@ -209,6 +214,10 @@ function _from_json(json; repair = true)
                 dict["Provides spinning reserves?"],
                 default = [true for t in 1:T],
             ),
+            timeseries(
+                dict["Provides flexible capacity?"],
+                default = [true for t in 1:T],
+            ),
             startup_categories,
         )
         push!(bus.units, unit)
@@ -216,11 +225,19 @@ function _from_json(json; repair = true)
         push!(units, unit)
     end
 
-    # Read reserves
+    # Read spinning, up-flexiramp, and down-flexiramp reserve requirements
-    reserves = Reserves(zeros(T))
+    reserves = Reserves(zeros(T), zeros(T), zeros(T))
     if "Reserves" in keys(json)
         reserves.spinning =
             timeseries(json["Reserves"]["Spinning (MW)"], default = zeros(T))
+        reserves.upflexiramp = timeseries(
+            json["Reserves"]["Up-flexiramp (MW)"],
+            default = zeros(T),
+        )
+        reserves.dwflexiramp = timeseries(
+            json["Reserves"]["Down-flexiramp (MW)"],
+            default = zeros(T),
+        )
     end
 
     # Read transmission lines
@@ -296,6 +313,7 @@ function _from_json(json; repair = true)
         price_sensitive_loads = loads,
         reserves = reserves,
         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

@@ -37,6 +37,7 @@ mutable struct Unit
     initial_status::Union{Int,Nothing}
     initial_power::Union{Float64,Nothing}
     provides_spinning_reserves::Vector{Bool}
+    provides_flexiramp_reserves::Vector{Bool}
     startup_categories::Vector{StartupCategory}
 end
 
@@ -54,6 +55,8 @@ end
 
 mutable struct Reserves
     spinning::Vector{Float64}
+    upflexiramp::Vector{Float64}
+    dwflexiramp::Vector{Float64}
 end
 
 mutable struct Contingency
@@ -81,6 +84,7 @@ Base.@kwdef mutable struct UnitCommitmentInstance
     price_sensitive_loads::Vector{PriceSensitiveLoad}
     reserves::Reserves
     shortfall_penalty::Vector{Float64}
+    flexiramp_shortfall_penalty::Vector{Float64}
     time::Int
     units_by_name::Dict{AbstractString,Unit}
     units::Vector{Unit}

src/model/build.jl

@@ -32,6 +32,22 @@ function build_model(;
     formulation = Formulation(),
     variable_names::Bool = false,
 )::JuMP.Model
+    if formulation.ramping == WanHob2016.Ramping() &&
+       instance.reserves.spinning >= ones(instance.time) .* 1e-6
+        error(
+            "Spinning reserves are not supported by the WanHob2016 ramping formulation",
+        )
+    end
+
+    if formulation.ramping !== WanHob2016.Ramping() && (
+        instance.reserves.upflexiramp >= ones(instance.time) .* 1e-6 ||
+        instance.reserves.dwflexiramp >= ones(instance.time) .* 1e-6
+    )
+        error(
+            "Flexiramp is supported only by the WanHob2016 ramping formulation",
+        )
+    end
+
     @info "Building model..."
     time_model = @elapsed begin
         model = Model()

src/model/formulations/WanHob2016/ramp.jl

@@ -0,0 +1,185 @@
+# 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_flexiramp_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 t in 1:model[:instance].time
+        # maximum feasible generation, \bar{g_{its}} in Wang & Hobbs (2016)
+        mfg[g.name, t] = @variable(model, lower_bound = 0)
+        if g.provides_flexiramp_reserves[t]
+            upflexiramp[g.name, t] = @variable(model) # up-flexiramp, ur_{it} in Wang & Hobbs (2016)
+            dwflexiramp[g.name, t] = @variable(model) # down-flexiramp, dr_{it} in Wang & Hobbs (2016)
+        else
+            upflexiramp[g.name, t] = 0.0
+            dwflexiramp[g.name, t] = 0.0
+        end
+        upflexiramp_shortfall[t] =
+            (model[:instance].flexiramp_shortfall_penalty[t] >= 0) ?
+            @variable(model, lower_bound = 0) : 0.0
+        dwflexiramp_shortfall[t] =
+            (model[:instance].flexiramp_shortfall_penalty[t] >= 0) ?
+            @variable(model, lower_bound = 0) : 0.0
+    end
+    return
+end
+
+function _add_ramp_eqs!(
+    model::JuMP.Model,
+    g::Unit,
+    formulation_prod_vars::Gar1962.ProdVars,
+    formulation_ramping::WanHob2016.Ramping,
+    formulation_status_vars::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]
+
+    for t in 1:model[:instance].time
+        @constraint(
+            model,
+            prod_above[gn, t] + (is_on[gn, t] * minp[t]) <= mfg[gn, t]
+        ) # Eq. (19) in Wang & Hobbs (2016)
+        @constraint(model, mfg[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[gn, t] +
+                (is_on[gn, t] * minp[t])
+            ) # first inequality of Eq. (20) in Wang & Hobbs (2016)
+            @constraint(
+                model,
+                prod_above[gn, t] - dwflexiramp[gn, t] +
+                (is_on[gn, t] * minp[t]) <=
+                mfg[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[gn, t] +
+                (is_on[gn, t] * minp[t])
+            ) # first inequality of Eq. (21) in Wang & Hobbs (2016)
+            @constraint(
+                model,
+                prod_above[gn, t] +
+                upflexiramp[gn, t] +
+                (is_on[gn, t] * minp[t]) <=
+                mfg[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[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[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[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[gn, t]
+            ) # first inequality of Eq. (26) in Wang & Hobbs (2016)
+            @constraint(
+                model,
+                upflexiramp[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[gn, t]
+            ) # first inequality of Eq. (27) in Wang & Hobbs (2016)
+            @constraint(
+                model,
+                dwflexiramp[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[gn, t]
+            ) # first inequality of Eq. (28) in Wang & Hobbs (2016)
+            @constraint(model, upflexiramp[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[gn, t]) # first inequality of Eq. (29) in Wang & Hobbs (2016)
+            @constraint(
+                model,
+                dwflexiramp[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[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

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

@@ -5,6 +5,7 @@
 function _add_system_wide_eqs!(model::JuMP.Model)::Nothing
     _add_net_injection_eqs!(model)
     _add_reserve_eqs!(model)
+    _add_flexiramp_eqs!(model)
     return
 end
 
@@ -54,3 +55,48 @@ function _add_reserve_eqs!(model::JuMP.Model)::Nothing
     end
     return
 end
+
+function _add_flexiramp_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[t] and eq_min_dwflexiramp[t] 
+    #       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 t in 1:instance.time
+        flexiramp_shortfall_penalty = instance.flexiramp_shortfall_penalty[t]
+        # Eq. (17) in Wang & Hobbs (2016)
+        eq_min_upflexiramp[t] = @constraint(
+            model,
+            sum(model[:upflexiramp][g.name, t] for g in instance.units) +
+            (
+                flexiramp_shortfall_penalty >= 0 ?
+                model[:upflexiramp_shortfall][t] : 0.0
+            ) >= instance.reserves.upflexiramp[t]
+        )
+        # Eq. (18) in Wang & Hobbs (2016)
+        eq_min_dwflexiramp[t] = @constraint(
+            model,
+            sum(model[:dwflexiramp][g.name, t] for g in instance.units) +
+            (
+                flexiramp_shortfall_penalty >= 0 ?
+                model[:dwflexiramp_shortfall][t] : 0.0
+            ) >= instance.reserves.dwflexiramp[t]
+        )
+
+        # Account for flexiramp shortfall contribution to objective
+        if flexiramp_shortfall_penalty >= 0
+            add_to_expression!(
+                model[:obj],
+                flexiramp_shortfall_penalty,
+                (
+                    model[:upflexiramp_shortfall][t] +
+                    model[:dwflexiramp_shortfall][t]
+                ),
+            )
+        end
+    end
+    return
+end

src/model/formulations/base/unit.jl

@@ -13,6 +13,7 @@ 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_flexiramp_vars!(model, g)
     _add_startup_shutdown_vars!(model, g)
     _add_status_vars!(model, g, formulation.status_vars)
 

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

@@ -50,13 +50,37 @@ 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)
+    if instance.reserves.upflexiramp != zeros(T) ||
-    sol["Reserve shortfall (MW)"] = OrderedDict(
+       instance.reserves.dwflexiramp != zeros(T)
-        t =>
+        # Report flexiramp solutions only if either of the up-flexiramp and  
-            (instance.shortfall_penalty[t] >= 0) ?
+        # down-flexiramp requirements is not a default array of zeros
-            round(value(model[:reserve_shortfall][t]), digits = 5) : 0.0 for
+        sol["Up-flexiramp (MW)"] =
-        t in 1:instance.time
+            timeseries(model[:upflexiramp], instance.units)
-    )
+        sol["Up-flexiramp shortfall (MW)"] = OrderedDict(
+            t =>
+                (instance.flexiramp_shortfall_penalty[t] >= 0) ?
+                round(value(model[:upflexiramp_shortfall][t]), digits = 5) :
+                0.0 for t in 1:instance.time
+        )
+        sol["Down-flexiramp (MW)"] =
+            timeseries(model[:dwflexiramp], instance.units)
+        sol["Down-flexiramp shortfall (MW)"] = OrderedDict(
+            t =>
+                (instance.flexiramp_shortfall_penalty[t] >= 0) ?
+                round(value(model[:dwflexiramp_shortfall][t]), digits = 5) :
+                0.0 for t in 1:instance.time
+        )
+    else
+        # Report spinning reserve solutions only if both up-flexiramp and  
+        # down-flexiramp requirements are arrays of zeros.
+        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
+        )
+    end
     sol["Net injection (MW)"] =
         timeseries(model[:net_injection], instance.buses)
     sol["Load curtail (MW)"] = timeseries(model[:curtail], instance.buses)

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

@@ -321,22 +321,68 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01)
             err_count += 1
         end
 
-        # Verify spinning reserves
+        # Verify flexiramp solutions only if either of the up-flexiramp and  
-        reserve =
+        # down-flexiramp requirements is not a default array of zeros
-            sum(solution["Reserve (MW)"][g.name][t] for g in instance.units)
+        if instance.reserves.upflexiramp != zeros(instance.time) ||
-        reserve_shortfall =
+           instance.reserves.dwflexiramp != zeros(instance.time)
-            (instance.shortfall_penalty[t] >= 0) ?
+            upflexiramp = sum(
-            solution["Reserve shortfall (MW)"][t] : 0
+                solution["Up-flexiramp (MW)"][g.name][t] for
-
+                g in instance.units
-        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
+            upflexiramp_shortfall =
+                (instance.flexiramp_shortfall_penalty[t] >= 0) ?
+                solution["Up-flexiramp shortfall (MW)"][t] : 0
+
+            if upflexiramp + upflexiramp_shortfall <
+               instance.reserves.upflexiramp[t] - tol
+                @error @sprintf(
+                    "Insufficient up-flexiramp at time %d (%.2f + %.2f should be %.2f)",
+                    t,
+                    upflexiramp,
+                    upflexiramp_shortfall,
+                    instance.reserves.upflexiramp[t],
+                )
+                err_count += 1
+            end
+
+            dwflexiramp = sum(
+                solution["Down-flexiramp (MW)"][g.name][t] for
+                g in instance.units
+            )
+            dwflexiramp_shortfall =
+                (instance.flexiramp_shortfall_penalty[t] >= 0) ?
+                solution["Down-flexiramp shortfall (MW)"][t] : 0
+
+            if dwflexiramp + dwflexiramp_shortfall <
+               instance.reserves.dwflexiramp[t] - tol
+                @error @sprintf(
+                    "Insufficient down-flexiramp at time %d (%.2f + %.2f should be %.2f)",
+                    t,
+                    dwflexiramp,
+                    dwflexiramp_shortfall,
+                    instance.reserves.dwflexiramp[t],
+                )
+                err_count += 1
+            end
+            # Verify spinning reserve solutions only if both up-flexiramp and  
+            # down-flexiramp requirements are arrays of zeros.
+        else
+            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
+
+            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
+            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

@@ -4,6 +4,7 @@
 
 using UnitCommitment
 using JuMP
+using Cbc
 import UnitCommitment:
     ArrCon2000,
     CarArr2006,
@@ -13,32 +14,20 @@ import UnitCommitment:
     KnuOstWat2018,
     MorLatRam2013,
     PanGua2016,
-    XavQiuWanThi2019
+    XavQiuWanThi2019,
+    WanHob2016
 
-if ENABLE_LARGE_TESTS
+function _test(
-    using Gurobi
+    formulation::Formulation;
-end
+    instances::Array{String} = ["test/case14"],
-
+    optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0),
-function _small_test(formulation::Formulation)::Nothing
+)::Nothing
-    instances = ["matpower/case118/2017-02-01", "test/case14"]
+    for instance_name in instances
-    for instance in instances
+        instance = UnitCommitment.read_benchmark(instance_name)
-        # Should not crash
-        UnitCommitment.build_model(
-            instance = UnitCommitment.read_benchmark(instance),
-            formulation = formulation,
-        )
-    end
-    return
-end
-
-function _large_test(formulation::Formulation)::Nothing
-    instances = ["pglib-uc/ca/Scenario400_reserves_1"]
-    for instance in instances
-        instance = UnitCommitment.read_benchmark(instance)
         model = UnitCommitment.build_model(
             instance = instance,
             formulation = formulation,
-            optimizer = Gurobi.Optimizer,
+            optimizer = optimizer,
         )
         UnitCommitment.optimize!(
             model,
@@ -47,19 +36,12 @@ function _large_test(formulation::Formulation)::Nothing
         solution = UnitCommitment.solution(model)
         @test UnitCommitment.validate(instance, solution)
     end
-    return
-end
-
-function _test(formulation::Formulation)::Nothing
-    _small_test(formulation)
-    if ENABLE_LARGE_TESTS
-        _large_test(formulation)
-    end
 end
 
 @testset "formulations" begin
     _test(Formulation())
     _test(Formulation(ramping = ArrCon2000.Ramping()))
+
     # _test(Formulation(ramping = DamKucRajAta2016.Ramping()))
     _test(
         Formulation(
@@ -71,4 +53,8 @@ end
     _test(Formulation(pwl_costs = Gar1962.PwlCosts()))
     _test(Formulation(pwl_costs = CarArr2006.PwlCosts()))
     _test(Formulation(pwl_costs = KnuOstWat2018.PwlCosts()))
+    _test(
+        Formulation(ramping = WanHob2016.Ramping()),
+        instances = ["test/case14-flex"],
+    )
 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

@@ -35,7 +35,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(