diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 4ee6045..a217247 100644 --- a/.github/workflows/test.yml +++ b/.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: diff --git a/.gitignore b/.gitignore index bfcf178..4ce337b 100644 --- a/.gitignore +++ b/.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 diff --git a/Makefile b/Makefile index a331d19..0d2446f 100644 --- a/Makefile +++ b/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 diff --git a/Project.toml b/Project.toml index ce2d6a6..d7deab8 100644 --- a/Project.toml +++ b/Project.toml @@ -25,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" diff --git a/README.md b/README.md index a2427d7..39248fa 100755 --- a/README.md +++ b/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 diff --git a/docs/format.md b/docs/format.md index 02ef89f..979bcbe 100644 --- a/docs/format.md +++ b/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. diff --git a/instances/test/case14-flex.json.gz b/instances/test/case14-flex.json.gz new file mode 100644 index 0000000..8e88af3 Binary files /dev/null and b/instances/test/case14-flex.json.gz differ diff --git a/juliaw b/juliaw old mode 100644 new mode 100755 index ad5826e..b78bc72 --- a/juliaw +++ b/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)") diff --git a/src/UnitCommitment.jl b/src/UnitCommitment.jl index 9ada14c..89049a1 100644 --- a/src/UnitCommitment.jl +++ b/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") diff --git a/src/instance/read.jl b/src/instance/read.jl index d3c6a5b..90cb169 100644 --- a/src/instance/read.jl +++ b/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, diff --git a/src/instance/structs.jl b/src/instance/structs.jl index 91e87bb..f07f694 100644 --- a/src/instance/structs.jl +++ b/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} diff --git a/src/model/formulations/WanHob2016/ramp.jl b/src/model/formulations/WanHob2016/ramp.jl new file mode 100644 index 0000000..8321273 --- /dev/null +++ b/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 diff --git a/src/model/formulations/WanHob2016/structs.jl b/src/model/formulations/WanHob2016/structs.jl new file mode 100644 index 0000000..52f621b --- /dev/null +++ b/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 diff --git a/src/model/formulations/base/system.jl b/src/model/formulations/base/system.jl index a265b2e..79832ea 100644 --- a/src/model/formulations/base/system.jl +++ b/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 diff --git a/src/model/formulations/base/unit.jl b/src/model/formulations/base/unit.jl index 1ed6571..ea39cb9 100644 --- a/src/model/formulations/base/unit.jl +++ b/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 - i in 1:(g.min_uptime-g.initial_status) if i <= T + switch_off[g.name, i] for 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 - i in 1:(g.min_downtime+g.initial_status) if i <= T + switch_on[g.name, i] for 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 - t in 1:model[:instance].time + sum(model[:reserve][r.name, g.name, t] for r in spinning_reserves) for t in 1:model[:instance].time ] end return reserve diff --git a/src/solution/fix.jl b/src/solution/fix.jl index 114e3ce..b2d37d8 100644 --- a/src/solution/fix.jl +++ b/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( diff --git a/src/solution/methods/XavQiuWanThi2019/optimize.jl b/src/solution/methods/XavQiuWanThi2019/optimize.jl index b513bd1..78c8381 100644 --- a/src/solution/methods/XavQiuWanThi2019/optimize.jl +++ b/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 diff --git a/src/solution/methods/XavQiuWanThi2019/structs.jl b/src/solution/methods/XavQiuWanThi2019/structs.jl index fc16383..4125563 100644 --- a/src/solution/methods/XavQiuWanThi2019/structs.jl +++ b/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 diff --git a/src/solution/solution.jl b/src/solution/solution.jl index 57d6ef7..bc5f0a6 100644 --- a/src/solution/solution.jl +++ b/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 - k in 1:length(g.cost_segments) + g.cost_segments[k].cost[t] for 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 - k in 1:length(g.cost_segments) + value(model[:segprod][g.name, t, k]) for 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 - t in 1:instance.time + value(model[:reserve][r.name, g.name, t]) for 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 - t in 1:instance.time - ] for r in instance.reserves + 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 diff --git a/src/transform/randomize/XavQiuAhm2021.jl b/src/transform/randomize/XavQiuAhm2021.jl index adedcd2..a3716c3 100644 --- a/src/transform/randomize/XavQiuAhm2021.jl +++ b/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 diff --git a/src/utils/log.jl b/src/utils/log.jl index b69a7d5..6b5f7c5 100644 --- a/src/utils/log.jl +++ b/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 diff --git a/src/validation/validate.jl b/src/validation/validate.jl index 27638b1..8da2904 100644 --- a/src/validation/validate.jl +++ b/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 - r in unit.reserves + solution["Spinning reserve (MW)"][r.name][unit.name] + 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 && - (unit in keys(solution["Reserve (MW)"][r.name])) - @error @sprintf( - "Unit %s is not eligible to provide reserve %s", - unit.name, - r.name, - ) - err_count += 1 + 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 @@ -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 - ps in instance.price_sensitive_loads + solution["Price-sensitive loads (MW)"][ps.name][t] for 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 - b in instance.buses + solution["Load curtail (MW)"][b.name][t] for 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( - solution["Reserve (MW)"][r.name][g.name][t] for g in r.units - ) - shortfall = solution["Reserve shortfall (MW)"][r.name][t] - required = r.amount[t] + 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 provided + shortfall < required - tol - @error @sprintf( - "Insufficient reserve %s at time %d (%.2f + %.2f < %.2f)", - r.name, - t, - provided, - shortfall, - required, + 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 ) - 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 diff --git a/test/Project.toml b/test/Project.toml index 7d70fe1..b87293d 100644 --- a/test/Project.toml +++ b/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" diff --git a/test/model/formulations_test.jl b/test/model/formulations_test.jl index eba99e0..978b8e8 100644 --- a/test/model/formulations_test.jl +++ b/test/model/formulations_test.jl @@ -15,56 +15,36 @@ import UnitCommitment: KnuOstWat2018, MorLatRam2013, PanGua2016, - XavQiuWanThi2019 + XavQiuWanThi2019, + WanHob2016 -if ENABLE_LARGE_TESTS - using Gurobi -end - -function _small_test(formulation::Formulation; dump::Bool = false)::Nothing - instance = UnitCommitment.read_benchmark("test/case14") - model = UnitCommitment.build_model( - instance = instance, - formulation = formulation, - optimizer = Cbc.Optimizer, - variable_names = true, - ) - UnitCommitment.optimize!(model) - solution = UnitCommitment.solution(model) - if dump - open("/tmp/ucjl.json", "w") do f - return write(f, JSON.json(solution, 2)) +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 = 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 - 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 diff --git a/test/runtests.jl b/test/runtests.jl index 7f1a1be..8d43ca4 100644 --- a/test/runtests.jl +++ b/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") diff --git a/test/transform/randomize/XavQiuAhm2021_test.jl b/test/transform/randomize/XavQiuAhm2021_test.jl index a3107c4..667fcf1 100644 --- a/test/transform/randomize/XavQiuAhm2021_test.jl +++ b/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] ≈ diff --git a/test/transform/slice_test.jl b/test/transform/slice_test.jl index 84b0f84..b736d16 100644 --- a/test/transform/slice_test.jl +++ b/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(