Implement randomization method from XavQiuAhm2021

.github/workflows/test.yml

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

Project.toml

@@ -15,6 +15,7 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
@@ -30,8 +31,8 @@ julia = "1"
 
 [extras]
 Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["Cbc", "Test", "Gurobi"]

docs/format.md

@@ -28,13 +28,14 @@ Each section is described in detail below. For a complete example, see [case14](
 
 ### Parameters
 
-This section describes system-wide parameters, such as power balance penalties,  optimization parameters, such as the length of the planning horizon and the time.
+This section describes system-wide parameters, such as power balance and reserve shortfall penalties, and optimization parameters, such as the length of the planning horizon and the time.
 
 | Key                            | Description                                       | Default  | Time series?
 | :----------------------------- | :------------------------------------------------ | :------: | :------------:
-| `Time horizon (h)`                     | Length of the planning horizon (in hours). | Required | N
+| `Time horizon (h)` | Length of the planning horizon (in hours). | Required | N
 | `Time step (min)` | Length of each time step (in minutes). Must be a divisor of 60 (e.g. 60, 30, 20, 15, etc). | `60` | N
 | `Power balance penalty ($/MW)` | Penalty for system-wide shortage or surplus in production (in $/MW). This is charged per time step. For example, if there is a shortage of 1 MW for three time steps, three times this amount will be charged. | `1000.0` | Y
+| `Reserve shortfall penalty ($/MW)` | Penalty for system-wide shortage in meeting reserve requirements (in $/MW). This is charged per time step. Negative value implies reserve constraints must always be satisfied. | `-1` | Y
 
 
 #### Example
@@ -42,7 +43,8 @@ This section describes system-wide parameters, such as power balance penalties,
 {
     "Parameters": {
         "Time horizon (h)": 4,
-        "Power balance penalty ($/MW)": 1000.0
+        "Power balance penalty ($/MW)": 1000.0,
+        "Reserve shortfall penalty ($/MW)": -1.0
     }
 }
 ```

src/UnitCommitment.jl

@@ -48,7 +48,7 @@ include("solution/warmstart.jl")
 include("solution/write.jl")
 include("transform/initcond.jl")
 include("transform/slice.jl")
-include("transform/randomize.jl")
+include("transform/randomize/XavQiuAhm2021.jl")
 include("utils/log.jl")
 include("validation/repair.jl")
 include("validation/validate.jl")

src/instance/read.jl

@@ -98,6 +98,10 @@ function _from_json(json; repair = true)
         json["Parameters"]["Power balance penalty (\$/MW)"],
         default = [1000.0 for t in 1:T],
     )
+    shortfall_penalty = timeseries(
+        json["Parameters"]["Reserve shortfall penalty (\$/MW)"],
+        default = [-1.0 for t in 1:T],
+    )
 
     # Read buses
     for (bus_name, dict) in json["Buses"]
@@ -264,6 +268,7 @@ function _from_json(json; repair = true)
     instance = UnitCommitmentInstance(
         T,
         power_balance_penalty,
+        shortfall_penalty,
         units,
         buses,
         lines,

src/instance/structs.jl

@@ -72,6 +72,8 @@ end
 mutable struct UnitCommitmentInstance
     time::Int
     power_balance_penalty::Vector{Float64}
+    "Penalty for failing to meet reserve requirement."
+    shortfall_penalty::Vector{Float64}
     units::Vector{Unit}
     buses::Vector{Bus}
     lines::Vector{TransmissionLine}

src/model/formulations/base/bus.jl

@@ -4,15 +4,11 @@
 
 function _add_bus!(model::JuMP.Model, b::Bus)::Nothing
     net_injection = _init(model, :expr_net_injection)
-    reserve = _init(model, :expr_reserve)
     curtail = _init(model, :curtail)
     for t in 1:model[:instance].time
         # Fixed load
         net_injection[b.name, t] = AffExpr(-b.load[t])
 
-        # Reserves
-        reserve[b.name, t] = AffExpr()
-
         # Load curtailment
         curtail[b.name, t] =
             @variable(model, lower_bound = 0, upper_bound = b.load[t])

src/model/formulations/base/system.jl

@@ -29,13 +29,28 @@ end
 
 function _add_reserve_eqs!(model::JuMP.Model)::Nothing
     eq_min_reserve = _init(model, :eq_min_reserve)
-    for t in 1:model[:instance].time
+    instance = model[:instance]
+    for t in 1:instance.time
+        # Equation (68) in Kneuven et al. (2020)
+        # As in Morales-España et al. (2013a)
+        # Akin to the alternative formulation with max_power_avail
+        # from Carrión and Arroyo (2006) and Ostrowski et al. (2012)
+        shortfall_penalty = instance.shortfall_penalty[t]
         eq_min_reserve[t] = @constraint(
             model,
-            sum(
+            sum(model[:reserve][g.name, t] for g in instance.units) +
-                model[:expr_reserve][b.name, t] for b in model[:instance].buses
+            (shortfall_penalty >= 0 ? model[:reserve_shortfall][t] : 0.0) >=
-            ) >= model[:instance].reserves.spinning[t]
+            instance.reserves.spinning[t]
         )
+
+        # Account for shortfall contribution to objective
+        if shortfall_penalty >= 0
+            add_to_expression!(
+                model[:obj],
+                shortfall_penalty,
+                model[:reserve_shortfall][t],
+            )
+        end
     end
     return
 end

src/model/formulations/base/unit.jl

@@ -44,12 +44,16 @@ _is_initially_on(g::Unit)::Float64 = (g.initial_status > 0 ? 1.0 : 0.0)
 
 function _add_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
     reserve = _init(model, :reserve)
+    reserve_shortfall = _init(model, :reserve_shortfall)
     for t in 1:model[:instance].time
         if g.provides_spinning_reserves[t]
             reserve[g.name, t] = @variable(model, lower_bound = 0)
         else
             reserve[g.name, t] = 0.0
         end
+        reserve_shortfall[t] =
+            (model[:instance].shortfall_penalty[t] >= 0) ?
+            @variable(model, lower_bound = 0) : 0.0
     end
     return
 end
@@ -210,11 +214,5 @@ function _add_net_injection_eqs!(model::JuMP.Model, g::Unit)::Nothing
             model[:is_on][g.name, t],
             g.min_power[t],
         )
-        # Add to reserves expression
-        add_to_expression!(
-            model[:expr_reserve][g.bus.name, t],
-            model[:reserve][g.name, t],
-            1.0,
-        )
     end
 end

src/solution/solution.jl

@@ -51,6 +51,12 @@ function solution(model::JuMP.Model)::OrderedDict
     sol["Switch on"] = timeseries(model[:switch_on], instance.units)
     sol["Switch off"] = timeseries(model[:switch_off], instance.units)
     sol["Reserve (MW)"] = timeseries(model[:reserve], instance.units)
+    sol["Reserve shortfall (MW)"] = OrderedDict(
+        t =>
+            (instance.shortfall_penalty[t] >= 0) ?
+            round(value(model[:reserve_shortfall][t]), digits = 5) : 0.0 for
+        t in 1:instance.time
+    )
     sol["Net injection (MW)"] =
         timeseries(model[:net_injection], instance.buses)
     sol["Load curtail (MW)"] = timeseries(model[:curtail], instance.buses)

src/transform/randomize.jl

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

src/transform/randomize/XavQiuAhm2021.jl

@@ -0,0 +1,209 @@
+# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
+# Copyright (C) 2020-2021, UChicago Argonne, LLC. All rights reserved.
+# Released under the modified BSD license. See COPYING.md for more details.
+
+"""
+Methods described in:
+
+    Xavier, Álinson S., Feng Qiu, and Shabbir Ahmed. "Learning to solve
+    large-scale security-constrained unit commitment problems." INFORMS
+    Journal on Computing 33.2 (2021): 739-756. DOI: 10.1287/ijoc.2020.0976
+"""
+module XavQiuAhm2021
+
+using Distributions
+import ..UnitCommitmentInstance
+
+"""
+    struct Randomization
+        cost = Uniform(0.95, 1.05)
+        load_profile_mu = [...]
+        load_profile_sigma = [...]
+        load_share = Uniform(0.90, 1.10)
+        peak_load = Uniform(0.6 * 0.925, 0.6 * 1.075)
+        randomize_costs = true
+        randomize_load_profile = true
+        randomize_load_share = true
+    end
+
+Randomization method that changes: (1) production and startup costs, (2)
+share of load coming from each bus, (3) peak system load, and (4) temporal
+load profile, as follows:
+
+1. **Production and startup costs:**
+    For each unit `u`, the vectors `u.min_power_cost` and `u.cost_segments`
+    are multiplied by a constant `α[u]` sampled from the provided `cost`
+    distribution. If `randomize_costs` is false, skips this step.
+
+2. **Load share:**
+   For each bus `b` and time `t`, the value `b.load[t]` is multiplied by
+   `(β[b] * b.load[t]) / sum(β[b2] * b2.load[t] for b2 in buses)`, where
+   `β[b]` is sampled from the provided `load_share` distribution. If
+   `randomize_load_share` is false, skips this step.
+
+3. **Peak system load and temporal load profile:**
+    Sets the peak load to `ρ * C`, where `ρ` is sampled from `peak_load` and `C`
+    is the maximum system capacity, at any time. Also scales the loads of all
+    buses, so that `system_load[t+1]` becomes equal to `system_load[t] * γ[t]`,
+    where `γ[t]` is sampled from `Normal(load_profile_mu[t], load_profile_sigma[t])`.
+    
+    The system load for the first time period is set so that the peak load
+    matches `ρ * C`. If `load_profile_sigma` and `load_profile_mu` have fewer
+    elements than `instance.time`, wraps around. If `randomize_load_profile`
+    is false, skips this step.
+
+The default parameters were obtained based on an analysis of publicly available
+bid and hourly data from PJM, corresponding to the month of January, 2017. For
+more details, see Section 4.2 of the paper.
+"""
+Base.@kwdef struct Randomization
+    cost = Uniform(0.95, 1.05)
+    load_profile_mu::Vector{Float64} = [
+        1.0,
+        0.978,
+        0.98,
+        1.004,
+        1.02,
+        1.078,
+        1.132,
+        1.018,
+        0.999,
+        1.006,
+        0.999,
+        0.987,
+        0.975,
+        0.984,
+        0.995,
+        1.005,
+        1.045,
+        1.106,
+        0.981,
+        0.981,
+        0.978,
+        0.948,
+        0.928,
+        0.953,
+    ]
+    load_profile_sigma::Vector{Float64} = [
+        0.0,
+        0.011,
+        0.015,
+        0.01,
+        0.012,
+        0.029,
+        0.055,
+        0.027,
+        0.026,
+        0.023,
+        0.013,
+        0.012,
+        0.014,
+        0.011,
+        0.008,
+        0.008,
+        0.02,
+        0.02,
+        0.016,
+        0.012,
+        0.014,
+        0.015,
+        0.017,
+        0.024,
+    ]
+    load_share = Uniform(0.90, 1.10)
+    peak_load = Uniform(0.6 * 0.925, 0.6 * 1.075)
+    randomize_load_profile::Bool = true
+    randomize_costs::Bool = true
+    randomize_load_share::Bool = true
+end
+
+function _randomize_costs(
+    instance::UnitCommitmentInstance,
+    distribution,
+)::Nothing
+    for unit in instance.units
+        α = rand(distribution)
+        unit.min_power_cost *= α
+        for k in unit.cost_segments
+            k.cost *= α
+        end
+        for s in unit.startup_categories
+            s.cost *= α
+        end
+    end
+    return
+end
+
+function _randomize_load_share(
+    instance::UnitCommitmentInstance,
+    distribution,
+)::Nothing
+    α = rand(distribution, length(instance.buses))
+    for t in 1:instance.time
+        total = sum(bus.load[t] for bus in instance.buses)
+        den = sum(
+            bus.load[t] / total * α[i] for
+            (i, bus) in enumerate(instance.buses)
+        )
+        for (i, bus) in enumerate(instance.buses)
+            bus.load[t] *= α[i] / den
+        end
+    end
+    return
+end
+
+function _randomize_load_profile(
+    instance::UnitCommitmentInstance,
+    params::Randomization,
+)::Nothing
+    # Generate new system load
+    system_load = [1.0]
+    for t in 2:instance.time
+        idx = (t - 1) % length(params.load_profile_mu) + 1
+        gamma = rand(
+            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
+    system_load = system_load ./ maximum(system_load) .* peak_load
+
+    # Scale bus loads to match the new system load
+    prev_system_load = sum(b.load for b in instance.buses)
+    for b in instance.buses
+        for t in 1:instance.time
+            b.load[t] *= system_load[t] / prev_system_load[t]
+        end
+    end
+
+    return
+end
+
+end
+
+"""
+    function randomize!(
+        instance::UnitCommitment.UnitCommitmentInstance,
+        method::XavQiuAhm2021.Randomization,
+    )::Nothing
+
+Randomize costs and loads based on the method described in XavQiuAhm2021.
+"""
+function randomize!(
+    instance::UnitCommitment.UnitCommitmentInstance,
+    method::XavQiuAhm2021.Randomization,
+)::Nothing
+    if method.randomize_costs
+        XavQiuAhm2021._randomize_costs(instance, method.cost)
+    end
+    if method.randomize_load_share
+        XavQiuAhm2021._randomize_load_share(instance, method.load_share)
+    end
+    if method.randomize_load_profile
+        XavQiuAhm2021._randomize_load_profile(instance, method)
+    end
+    return
+end
+
+export randomize!

src/validation/validate.jl

@@ -324,11 +324,16 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01)
         # Verify spinning reserves
         reserve =
             sum(solution["Reserve (MW)"][g.name][t] for g in instance.units)
-        if reserve < instance.reserves.spinning[t] - tol
+        reserve_shortfall =
+            (instance.shortfall_penalty[t] >= 0) ?
+            solution["Reserve shortfall (MW)"][t] : 0
+
+        if reserve + reserve_shortfall < instance.reserves.spinning[t] - tol
             @error @sprintf(
-                "Insufficient spinning reserves at time %d (%.2f should be %.2f)",
+                "Insufficient spinning reserves at time %d (%.2f + %.2f should be %.2f)",
                 t,
                 reserve,
+                reserve_shortfall,
                 instance.reserves.spinning[t],
             )
             err_count += 1

test/model/formulations_test.jl

@@ -20,8 +20,14 @@ if ENABLE_LARGE_TESTS
 end
 
 function _small_test(formulation::Formulation)::Nothing
-    instance = UnitCommitment.read_benchmark("matpower/case118/2017-02-01")
+    instances = ["matpower/case118/2017-02-01", "test/case14"]
-    UnitCommitment.build_model(instance = instance, formulation = formulation)  # should not crash
+    for instance in instances
+        # Should not crash
+        UnitCommitment.build_model(
+            instance = UnitCommitment.read_benchmark(instance),
+            formulation = formulation,
+        )
+    end
     return
 end
 

test/runtests.jl

@@ -28,7 +28,9 @@ const ENABLE_LARGE_TESTS = ("UCJL_LARGE_TESTS" in keys(ENV))
     @testset "transform" begin
         include("transform/initcond_test.jl")
         include("transform/slice_test.jl")
-        include("transform/randomize_test.jl")
+        @testset "randomize" begin
+            include("transform/randomize/XavQiuAhm2021_test.jl")
+        end
     end
     @testset "validation" begin
         include("validation/repair_test.jl")

test/transform/randomize/XavQiuAhm2021_test.jl

@@ -0,0 +1,63 @@
+# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
+# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+# Released under the modified BSD license. See COPYING.md for more details.
+
+import Random
+import UnitCommitment: XavQiuAhm2021
+
+using Distributions
+using UnitCommitment, Cbc, JuMP
+
+get_instance() = UnitCommitment.read_benchmark("matpower/case118/2017-02-01")
+system_load(instance) = sum(b.load for b in instance.buses)
+test_approx(x, y) = @test isapprox(x, y, atol = 1e-3)
+
+@testset "XavQiuAhm2021" begin
+    @testset "cost and load share" begin
+        instance = get_instance()
+
+        # Check original costs
+        unit = instance.units[10]
+        test_approx(unit.min_power_cost[1], 825.023)
+        test_approx(unit.cost_segments[1].cost[1], 36.659)
+        test_approx(unit.startup_categories[1].cost[1], 7570.42)
+
+        # Check original load share
+        bus = instance.buses[1]
+        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),
+        )
+
+        # Check randomized costs
+        test_approx(unit.min_power_cost[1], 831.977)
+        test_approx(unit.cost_segments[1].cost[1], 36.968)
+        test_approx(unit.startup_categories[1].cost[1], 7634.226)
+
+        # Check randomized load share
+        curr_system_load = system_load(instance)
+        test_approx(bus.load[1] / curr_system_load[1], 0.013)
+
+        # System load should not change
+        @test prev_system_load ≈ curr_system_load
+    end
+
+    @testset "load profile" begin
+        instance = get_instance()
+
+        # Check original load profile
+        @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())
+
+        # Check randomized load profile
+        @test round.(system_load(instance), digits = 1)[1:8] ≈
+              [4854.7, 4849.2, 4732.7, 4848.2, 4948.4, 5231.1, 5874.8, 5934.8]
+    end
+end

test/transform/randomize_test.jl

@@ -1,43 +0,0 @@
-# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
-# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
-# Released under the modified BSD license. See COPYING.md for more details.
-
-using UnitCommitment, Cbc, JuMP
-
-_get_instance() = UnitCommitment.read_benchmark("matpower/case118/2017-02-01")
-_total_load(instance) = sum(b.load[1] for b in instance.buses)
-
-@testset "randomize_unit_costs!" begin
-    instance = _get_instance()
-    unit = instance.units[10]
-    prev_min_power_cost = unit.min_power_cost
-    prev_prod_cost = unit.cost_segments[1].cost
-    prev_startup_cost = unit.startup_categories[1].cost
-    randomize_unit_costs!(instance)
-    @test prev_min_power_cost != unit.min_power_cost
-    @test prev_prod_cost != unit.cost_segments[1].cost
-    @test prev_startup_cost != unit.startup_categories[1].cost
-end
-
-@testset "randomize_load_distribution!" begin
-    instance = _get_instance()
-    bus = instance.buses[1]
-    prev_load = instance.buses[1].load[1]
-    prev_total_load = _total_load(instance)
-    randomize_load_distribution!(instance)
-    curr_total_load = _total_load(instance)
-    @test prev_load != instance.buses[1].load[1]
-    @test abs(prev_total_load - curr_total_load) < 1e-3
-end
-
-@testset "randomize_peak_load!" begin
-    instance = _get_instance()
-    bus = instance.buses[1]
-    prev_total_load = _total_load(instance)
-    prev_share = bus.load[1] / prev_total_load
-    randomize_peak_load!(instance)
-    curr_total_load = _total_load(instance)
-    curr_share = bus.load[1] / prev_total_load
-    @test curr_total_load != prev_total_load
-    @test abs(curr_share - prev_share) < 1e-3
-end