Fix bug in validation script; create large tests

CHANGELOG.md

@@ -11,6 +11,10 @@ All notable changes to this project will be documented in this file.
 [semver]: https://semver.org/spec/v2.0.0.html
 [pkjjl]: https://pkgdocs.julialang.org/v1/compatibility/#compat-pre-1.0
 
+## [0.2.2] - Unreleased
+### Fixed
+- Fix small bug in validation scripts related to startup costs
+
 ## [0.2.1] - 2021-06-02
 ### Added
 - Add multiple ramping formulations (ArrCon2000, MorLatRam2013, DamKucRajAta2016, PanGua2016)

Project.toml

@@ -2,10 +2,11 @@ name = "UnitCommitment"
 uuid = "64606440-39ea-11e9-0f29-3303a1d3d877"
 authors = ["Santos Xavier, Alinson <axavier@anl.gov>"]
 repo = "https://github.com/ANL-CEEESA/UnitCommitment.jl"
-version = "0.2.1"
+version = "0.2.2"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
@@ -19,6 +20,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 [compat]
 Cbc = "0.7"
 DataStructures = "0.18"
+Distributions = "0.25"
 GZip = "0.5"
 JSON = "0.21"
 JuMP = "0.21"
@@ -29,6 +31,7 @@ julia = "1"
 [extras]
 Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b"
 
 [targets]
-test = ["Cbc", "Test"]
+test = ["Cbc", "Test", "Gurobi"]

README.md

@@ -20,7 +20,7 @@
 
 * **Data Format:** The package proposes an extensible and fully-documented JSON-based data format for SCUC, developed in collaboration with Independent System Operators (ISOs), which describes the most important aspects of the problem. The format supports the most common generator characteristics (including ramping, piecewise-linear production cost curves and time-dependent startup costs), as well as operating reserves, price-sensitive loads, transmission networks and contingencies.
 * **Benchmark Instances:** The package provides a diverse collection of large-scale benchmark instances collected from the literature, converted into a common data format, and extended using data-driven methods to make them more challenging and realistic.
-* **Model Implementation**: The package provides a Julia/JuMP implementations of state-of-the-art formulations and solution methods for SCUC, including multiple ramping formulations ([ArrCon2000][ArrCon2000], [MorLatRam2013][MorLatRam2013], [DamKucRajAta2016][DamKucRajAta2016], [PanGua2016][PanGua2016]), multiple piecewise-linear costs formulations ([Gar1962][Gar1962], [CarArr2006][CarArr2006], [KnuOstWat2018][KnuOstWat2018]) and contingency screening methods ([XavQiuWanThi2019][XavQiuWanThi2019]). Our goal is to keep these implementations up-to-date as new methods are proposed in the literature.
+* **Model Implementation**: The package provides Julia/JuMP implementations of state-of-the-art formulations and solution methods for SCUC, including multiple ramping formulations ([ArrCon2000][ArrCon2000], [MorLatRam2013][MorLatRam2013], [DamKucRajAta2016][DamKucRajAta2016], [PanGua2016][PanGua2016]), multiple piecewise-linear costs formulations ([Gar1962][Gar1962], [CarArr2006][CarArr2006], [KnuOstWat2018][KnuOstWat2018]) and contingency screening methods ([XavQiuWanThi2019][XavQiuWanThi2019]). Our goal is to keep these implementations up-to-date as new methods are proposed in the literature.
 * **Benchmark Tools:** The package provides automated benchmark scripts to accurately evaluate the performance impact of proposed code changes.
 
 [ArrCon2000]: https://doi.org/10.1109/59.871739

benchmark/scripts/table.py

@@ -6,6 +6,9 @@ from pathlib import Path
 import pandas as pd
 import re
 from tabulate import tabulate
+from colorama import init, Fore, Back, Style
+
+init()
 
 
 def process_all_log_files():
@@ -48,6 +51,7 @@ def process(filename):
     # m = re.search("case([0-9]*)", instance_name)
     # n_buses = int(m.group(1))
     n_buses = 0
+    validation_errors = 0
 
     with open(filename) as file:
         for line in file.readlines():
@@ -137,6 +141,14 @@ def process(filename):
             if m is not None:
                 transmission_count += 1
 
+            m = re.search(r".*Found ([0-9]*) validation errors", line)
+            if m is not None:
+                validation_errors += int(m.group(1))
+                print(
+                    f"{Fore.YELLOW}{Style.BRIGHT}Warning:{Style.RESET_ALL} {validation_errors:8d} "
+                    f"{Style.DIM}validation errors in {Style.RESET_ALL}{group_name}/{instance_name}/{sample_name}"
+                )
+
     return {
         "Group": group_name,
         "Instance": instance_name,
@@ -168,6 +180,7 @@ def process(filename):
         "Transmission screening constraints": transmission_count,
         "Transmission screening time": transmission_time,
         "Transmission screening calls": transmission_calls,
+        "Validation errors": validation_errors,
     }
 
 

docs/usage.md

@@ -71,7 +71,7 @@ Advanced usage
 
 ### Customizing the formulation
 
-By default, `build_model` uses a formulation that combines modeling components from different publications, and that has been carefully tested, using our own benchmark scripts, to provide good performance across a wide variety of instances. This default formulation is expected to change over time, as new methods are proposed in the literature. You can, however, construct your own formulation, based the modeling components that you choose, as shown in the next example.
+By default, `build_model` uses a formulation that combines modeling components from different publications, and that has been carefully tested, using our own benchmark scripts, to provide good performance across a wide variety of instances. This default formulation is expected to change over time, as new methods are proposed in the literature. You can, however, construct your own formulation, based on the modeling components that you choose, as shown in the next example.
 
 ```julia
 using Cbc
@@ -119,7 +119,10 @@ instance = UnitCommitment.read("instance.json")
 UnitCommitment.generate_initial_conditions!(instance, Cbc.Optimizer)
 
 # Construct and solve optimization model
-model = UnitCommitment.build_model(instance, Cbc.Optimizer)
+model = UnitCommitment.build_model(
+    instance=instance,
+    optimizer=Cbc.Optimizer,
+)
 UnitCommitment.optimize!(model)
 ```
 

src/UnitCommitment.jl

@@ -30,6 +30,8 @@ include("model/formulations/base/unit.jl")
 include("model/formulations/CarArr2006/pwlcosts.jl")
 include("model/formulations/DamKucRajAta2016/ramp.jl")
 include("model/formulations/Gar1962/pwlcosts.jl")
+include("model/formulations/Gar1962/status.jl")
+include("model/formulations/Gar1962/prod.jl")
 include("model/formulations/KnuOstWat2018/pwlcosts.jl")
 include("model/formulations/MorLatRam2013/ramp.jl")
 include("model/formulations/MorLatRam2013/scosts.jl")
@@ -46,6 +48,7 @@ include("solution/warmstart.jl")
 include("solution/write.jl")
 include("transform/initcond.jl")
 include("transform/slice.jl")
+include("transform/randomize.jl")
 include("utils/log.jl")
 include("validation/repair.jl")
 include("validation/validate.jl")

src/model/formulations/ArrCon2000/ramp.jl

@@ -5,7 +5,9 @@
 function _add_ramp_eqs!(
     model::JuMP.Model,
     g::Unit,
-    formulation::ArrCon2000.Ramping,
+    formulation_prod_vars::Gar1962.ProdVars,
+    formulation_ramping::ArrCon2000.Ramping,
+    formulation_status_vars::Gar1962.StatusVars,
 )::Nothing
     # TODO: Move upper case constants to model[:instance]
     RESERVES_WHEN_START_UP = true
@@ -17,15 +19,19 @@ function _add_ramp_eqs!(
     RD = g.ramp_down_limit
     SU = g.startup_limit
     SD = g.shutdown_limit
-    is_on = model[:is_on]
-    prod_above = model[:prod_above]
     reserve = model[:reserve]
-    switch_off = model[:switch_off]
-    switch_on = model[:switch_on]
     eq_ramp_down = _init(model, :eq_ramp_down)
     eq_ramp_up = _init(model, :eq_ramp_up)
     is_initially_on = (g.initial_status > 0)
 
+    # Gar1962.ProdVars
+    prod_above = model[:prod_above]
+
+    # Gar1962.StatusVars
+    is_on = model[:is_on]
+    switch_off = model[:switch_off]
+    switch_on = model[:switch_on]
+
     for t in 1:model[:instance].time
         # Ramp up limit
         if t == 1

src/model/formulations/CarArr2006/pwlcosts.jl

@@ -5,13 +5,18 @@
 function _add_production_piecewise_linear_eqs!(
     model::JuMP.Model,
     g::Unit,
-    formulation::CarArr2006.PwlCosts,
+    formulation_prod_vars::Gar1962.ProdVars,
+    formulation_pwl_costs::CarArr2006.PwlCosts,
+    formulation_status_vars::StatusVarsFormulation,
 )::Nothing
     eq_prod_above_def = _init(model, :eq_prod_above_def)
     eq_segprod_limit = _init(model, :eq_segprod_limit)
-    prod_above = model[:prod_above]
     segprod = model[:segprod]
     gn = g.name
+
+    # Gar1962.ProdVars
+    prod_above = model[:prod_above]
+
     K = length(g.cost_segments)
     for t in 1:model[:instance].time
         gn = g.name

src/model/formulations/DamKucRajAta2016/ramp.jl

@@ -5,7 +5,9 @@
 function _add_ramp_eqs!(
     model::JuMP.Model,
     g::Unit,
-    formulation::DamKucRajAta2016.Ramping,
+    formulation_prod_vars::Gar1962.ProdVars,
+    formulation_ramping::DamKucRajAta2016.Ramping,
+    formulation_status_vars::Gar1962.StatusVars,
 )::Nothing
     # TODO: Move upper case constants to model[:instance]
     RESERVES_WHEN_START_UP = true
@@ -21,9 +23,13 @@ function _add_ramp_eqs!(
     gn = g.name
     eq_str_ramp_down = _init(model, :eq_str_ramp_down)
     eq_str_ramp_up = _init(model, :eq_str_ramp_up)
-    is_on = model[:is_on]
-    prod_above = model[:prod_above]
     reserve = model[:reserve]
+
+    # Gar1962.ProdVars
+    prod_above = model[:prod_above]
+
+    # Gar1962.StatusVars
+    is_on = model[:is_on]
     switch_off = model[:switch_off]
     switch_on = model[:switch_on]
 

src/model/formulations/Gar1962/prod.jl

@@ -0,0 +1,50 @@
+# 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.
+
+function _add_production_vars!(
+    model::JuMP.Model,
+    g::Unit,
+    formulation_prod_vars::Gar1962.ProdVars,
+)::Nothing
+    prod_above = _init(model, :prod_above)
+    segprod = _init(model, :segprod)
+    for t in 1:model[:instance].time
+        for k in 1:length(g.cost_segments)
+            segprod[g.name, t, k] = @variable(model, lower_bound = 0)
+        end
+        prod_above[g.name, t] = @variable(model, lower_bound = 0)
+    end
+    return
+end
+
+function _add_production_limit_eqs!(
+    model::JuMP.Model,
+    g::Unit,
+    formulation_prod_vars::Gar1962.ProdVars,
+)::Nothing
+    eq_prod_limit = _init(model, :eq_prod_limit)
+    is_on = model[:is_on]
+    prod_above = model[:prod_above]
+    reserve = model[:reserve]
+    gn = g.name
+    for t in 1:model[:instance].time
+        # Objective function terms for production costs
+        # Part of (69) of Kneuven et al. (2020) as C^R_g * u_g(t) term
+        add_to_expression!(model[:obj], is_on[gn, t], g.min_power_cost[t])
+
+        # Production limit
+        # Equation (18) in Kneuven et al. (2020)
+        #   as \bar{p}_g(t) \le \bar{P}_g u_g(t)
+        # amk: this is a weaker version of (20) and (21) in Kneuven et al. (2020)
+        #      but keeping it here in case those are not present
+        power_diff = max(g.max_power[t], 0.0) - max(g.min_power[t], 0.0)
+        if power_diff < 1e-7
+            power_diff = 0.0
+        end
+        eq_prod_limit[gn, t] = @constraint(
+            model,
+            prod_above[gn, t] + reserve[gn, t] <= power_diff * is_on[gn, t]
+        )
+    end
+end

src/model/formulations/Gar1962/pwlcosts.jl

@@ -5,14 +5,21 @@
 function _add_production_piecewise_linear_eqs!(
     model::JuMP.Model,
     g::Unit,
-    formulation::Gar1962.PwlCosts,
+    formulation_prod_vars::Gar1962.ProdVars,
+    formulation_pwl_costs::Gar1962.PwlCosts,
+    formulation_status_vars::Gar1962.StatusVars,
 )::Nothing
     eq_prod_above_def = _init(model, :eq_prod_above_def)
     eq_segprod_limit = _init(model, :eq_segprod_limit)
-    is_on = model[:is_on]
-    prod_above = model[:prod_above]
     segprod = model[:segprod]
     gn = g.name
+
+    # Gar1962.ProdVars
+    prod_above = model[:prod_above]
+
+    # Gar1962.StatusVars
+    is_on = model[:is_on]
+
     K = length(g.cost_segments)
     for t in 1:model[:instance].time
         # Definition of production

src/model/formulations/Gar1962/status.jl

@@ -0,0 +1,61 @@
+# 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.
+
+function _add_status_vars!(
+    model::JuMP.Model,
+    g::Unit,
+    formulation_status_vars::Gar1962.StatusVars,
+)::Nothing
+    is_on = _init(model, :is_on)
+    switch_on = _init(model, :switch_on)
+    switch_off = _init(model, :switch_off)
+    for t in 1:model[:instance].time
+        if g.must_run[t]
+            is_on[g.name, t] = 1.0
+            switch_on[g.name, t] = (t == 1 ? 1.0 - _is_initially_on(g) : 0.0)
+            switch_off[g.name, t] = 0.0
+        else
+            is_on[g.name, t] = @variable(model, binary = true)
+            switch_on[g.name, t] = @variable(model, binary = true)
+            switch_off[g.name, t] = @variable(model, binary = true)
+        end
+    end
+    return
+end
+
+function _add_status_eqs!(
+    model::JuMP.Model,
+    g::Unit,
+    formulation_status_vars::Gar1962.StatusVars,
+)::Nothing
+    eq_binary_link = _init(model, :eq_binary_link)
+    eq_switch_on_off = _init(model, :eq_switch_on_off)
+    is_on = model[:is_on]
+    switch_off = model[:switch_off]
+    switch_on = model[:switch_on]
+    for t in 1:model[:instance].time
+        if !g.must_run[t]
+            # Link binary variables
+            if t == 1
+                eq_binary_link[g.name, t] = @constraint(
+                    model,
+                    is_on[g.name, t] - _is_initially_on(g) ==
+                    switch_on[g.name, t] - switch_off[g.name, t]
+                )
+            else
+                eq_binary_link[g.name, t] = @constraint(
+                    model,
+                    is_on[g.name, t] - is_on[g.name, t-1] ==
+                    switch_on[g.name, t] - switch_off[g.name, t]
+                )
+            end
+            # Cannot switch on and off at the same time
+            eq_switch_on_off[g.name, t] = @constraint(
+                model,
+                switch_on[g.name, t] + switch_off[g.name, t] <= 1
+            )
+        end
+    end
+    return
+end

src/model/formulations/Gar1962/structs.jl

@@ -14,7 +14,11 @@ Formulation described in:
 module Gar1962
 
 import ..PiecewiseLinearCostsFormulation
+import ..ProductionVarsFormulation
+import ..StatusVarsFormulation
 
+struct ProdVars <: ProductionVarsFormulation end
 struct PwlCosts <: PiecewiseLinearCostsFormulation end
+struct StatusVars <: StatusVarsFormulation end
 
 end

src/model/formulations/KnuOstWat2018/pwlcosts.jl

@@ -5,21 +5,27 @@
 function _add_production_piecewise_linear_eqs!(
     model::JuMP.Model,
     g::Unit,
-    formulation::KnuOstWat2018.PwlCosts,
+    formulation_prod_vars::Gar1962.ProdVars,
+    formulation_pwl_costs::KnuOstWat2018.PwlCosts,
+    formulation_status_vars::Gar1962.StatusVars,
 )::Nothing
     eq_prod_above_def = _init(model, :eq_prod_above_def)
     eq_segprod_limit_a = _init(model, :eq_segprod_limit_a)
     eq_segprod_limit_b = _init(model, :eq_segprod_limit_b)
     eq_segprod_limit_c = _init(model, :eq_segprod_limit_c)
-    prod_above = model[:prod_above]
     segprod = model[:segprod]
-    is_on = model[:is_on]
-    switch_on = model[:switch_on]
-    switch_off = model[:switch_off]
     gn = g.name
     K = length(g.cost_segments)
     T = model[:instance].time
 
+    # Gar1962.ProdVars
+    prod_above = model[:prod_above]
+
+    # Gar1962.StatusVars
+    is_on = model[:is_on]
+    switch_on = model[:switch_on]
+    switch_off = model[:switch_off]
+
     for t in 1:T
         for k in 1:K
             # Pbar^{k-1)

src/model/formulations/MorLatRam2013/ramp.jl

@@ -5,7 +5,9 @@
 function _add_ramp_eqs!(
     model::JuMP.Model,
     g::Unit,
-    formulation::MorLatRam2013.Ramping,
+    formulation_prod_vars::Gar1962.ProdVars,
+    formulation_ramping::MorLatRam2013.Ramping,
+    formulation_status_vars::Gar1962.StatusVars,
 )::Nothing
     # TODO: Move upper case constants to model[:instance]
     RESERVES_WHEN_START_UP = true
@@ -20,9 +22,13 @@ function _add_ramp_eqs!(
     gn = g.name
     eq_ramp_down = _init(model, :eq_ramp_down)
     eq_ramp_up = _init(model, :eq_str_ramp_up)
-    is_on = model[:is_on]
-    prod_above = model[:prod_above]
     reserve = model[:reserve]
+
+    # Gar1962.ProdVars
+    prod_above = model[:prod_above]
+
+    # Gar1962.StatusVars
+    is_on = model[:is_on]
     switch_off = model[:switch_off]
     switch_on = model[:switch_on]
 

src/model/formulations/PanGua2016/ramp.jl

@@ -1,16 +1,18 @@
+# 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.
+
 function _add_ramp_eqs!(
     model::JuMP.Model,
     g::Unit,
-    formulation::PanGua2016.Ramping,
+    formulation_prod_vars::Gar1962.ProdVars,
+    formulation_ramping::PanGua2016.Ramping,
+    formulation_status_vars::Gar1962.StatusVars,
 )::Nothing
     # TODO: Move upper case constants to model[:instance]
     RESERVES_WHEN_SHUT_DOWN = true
     gn = g.name
-    is_on = model[:is_on]
-    prod_above = model[:prod_above]
     reserve = model[:reserve]
-    switch_off = model[:switch_off]
-    switch_on = model[:switch_on]
     eq_str_prod_limit = _init(model, :eq_str_prod_limit)
     eq_prod_limit_ramp_up_extra_period =
         _init(model, :eq_prod_limit_ramp_up_extra_period)
@@ -23,6 +25,14 @@ function _add_ramp_eqs!(
     RD = g.ramp_down_limit # ramp down rate
     T = model[:instance].time
 
+    # Gar1962.ProdVars
+    prod_above = model[:prod_above]
+
+    # Gar1962.StatusVars
+    is_on = model[:is_on]
+    switch_off = model[:switch_off]
+    switch_on = model[:switch_on]
+
     for t in 1:T
         Pbar = g.max_power[t]
         if Pbar < 1e-7

src/model/formulations/base/structs.jl

@@ -6,20 +6,33 @@ abstract type TransmissionFormulation end
 abstract type RampingFormulation end
 abstract type PiecewiseLinearCostsFormulation end
 abstract type StartupCostsFormulation end
+abstract type StatusVarsFormulation end
+abstract type ProductionVarsFormulation end
 
 struct Formulation
+    prod_vars::ProductionVarsFormulation
     pwl_costs::PiecewiseLinearCostsFormulation
     ramping::RampingFormulation
     startup_costs::StartupCostsFormulation
+    status_vars::StatusVarsFormulation
     transmission::TransmissionFormulation
 
     function Formulation(;
+        prod_vars::ProductionVarsFormulation = Gar1962.ProdVars(),
         pwl_costs::PiecewiseLinearCostsFormulation = KnuOstWat2018.PwlCosts(),
         ramping::RampingFormulation = MorLatRam2013.Ramping(),
         startup_costs::StartupCostsFormulation = MorLatRam2013.StartupCosts(),
+        status_vars::StatusVarsFormulation = Gar1962.StatusVars(),
         transmission::TransmissionFormulation = ShiftFactorsFormulation(),
     )
-        return new(pwl_costs, ramping, startup_costs, transmission)
+        return new(
+            prod_vars,
+            pwl_costs,
+            ramping,
+            startup_costs,
+            status_vars,
+            transmission,
+        )
     end
 end
 

src/model/formulations/base/unit.jl

@@ -2,7 +2,7 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
-function _add_unit!(model::JuMP.Model, g::Unit, f::Formulation)
+function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation)
     if !all(g.must_run) && any(g.must_run)
         error("Partially must-run units are not currently supported")
     end
@@ -11,64 +11,37 @@ function _add_unit!(model::JuMP.Model, g::Unit, f::Formulation)
     end
 
     # Variables
-    _add_production_vars!(model, g)
+    _add_production_vars!(model, g, formulation.prod_vars)
     _add_reserve_vars!(model, g)
     _add_startup_shutdown_vars!(model, g)
-    _add_status_vars!(model, g)
+    _add_status_vars!(model, g, formulation.status_vars)
 
     # Constraints and objective function
     _add_min_uptime_downtime_eqs!(model, g)
     _add_net_injection_eqs!(model, g)
-    _add_production_limit_eqs!(model, g)
-    _add_production_piecewise_linear_eqs!(model, g, f.pwl_costs)
-    _add_ramp_eqs!(model, g, f.ramping)
-    _add_startup_cost_eqs!(model, g, f.startup_costs)
+    _add_production_limit_eqs!(model, g, formulation.prod_vars)
+    _add_production_piecewise_linear_eqs!(
+        model,
+        g,
+        formulation.prod_vars,
+        formulation.pwl_costs,
+        formulation.status_vars,
+    )
+    _add_ramp_eqs!(
+        model,
+        g,
+        formulation.prod_vars,
+        formulation.ramping,
+        formulation.status_vars,
+    )
+    _add_startup_cost_eqs!(model, g, formulation.startup_costs)
     _add_startup_shutdown_limit_eqs!(model, g)
-    _add_status_eqs!(model, g)
+    _add_status_eqs!(model, g, formulation.status_vars)
     return
 end
 
 _is_initially_on(g::Unit)::Float64 = (g.initial_status > 0 ? 1.0 : 0.0)
 
-function _add_production_vars!(model::JuMP.Model, g::Unit)::Nothing
-    prod_above = _init(model, :prod_above)
-    segprod = _init(model, :segprod)
-    for t in 1:model[:instance].time
-        for k in 1:length(g.cost_segments)
-            segprod[g.name, t, k] = @variable(model, lower_bound = 0)
-        end
-        prod_above[g.name, t] = @variable(model, lower_bound = 0)
-    end
-    return
-end
-
-function _add_production_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
-    eq_prod_limit = _init(model, :eq_prod_limit)
-    is_on = model[:is_on]
-    prod_above = model[:prod_above]
-    reserve = model[:reserve]
-    gn = g.name
-    for t in 1:model[:instance].time
-        # Objective function terms for production costs
-        # Part of (69) of Kneuven et al. (2020) as C^R_g * u_g(t) term
-        add_to_expression!(model[:obj], is_on[gn, t], g.min_power_cost[t])
-
-        # Production limit
-        # Equation (18) in Kneuven et al. (2020)
-        #   as \bar{p}_g(t) \le \bar{P}_g u_g(t)
-        # amk: this is a weaker version of (20) and (21) in Kneuven et al. (2020)
-        #      but keeping it here in case those are not present
-        power_diff = max(g.max_power[t], 0.0) - max(g.min_power[t], 0.0)
-        if power_diff < 1e-7
-            power_diff = 0.0
-        end
-        eq_prod_limit[gn, t] = @constraint(
-            model,
-            prod_above[gn, t] + reserve[gn, t] <= power_diff * is_on[gn, t]
-        )
-    end
-end
-
 function _add_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
     reserve = _init(model, :reserve)
     for t in 1:model[:instance].time
@@ -134,56 +107,6 @@ function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
     return
 end
 
-function _add_status_vars!(model::JuMP.Model, g::Unit)::Nothing
-    is_on = _init(model, :is_on)
-    switch_on = _init(model, :switch_on)
-    switch_off = _init(model, :switch_off)
-    for t in 1:model[:instance].time
-        if g.must_run[t]
-            is_on[g.name, t] = 1.0
-            switch_on[g.name, t] = (t == 1 ? 1.0 - _is_initially_on(g) : 0.0)
-            switch_off[g.name, t] = 0.0
-        else
-            is_on[g.name, t] = @variable(model, binary = true)
-            switch_on[g.name, t] = @variable(model, binary = true)
-            switch_off[g.name, t] = @variable(model, binary = true)
-        end
-    end
-    return
-end
-
-function _add_status_eqs!(model::JuMP.Model, g::Unit)::Nothing
-    eq_binary_link = _init(model, :eq_binary_link)
-    eq_switch_on_off = _init(model, :eq_switch_on_off)
-    is_on = model[:is_on]
-    switch_off = model[:switch_off]
-    switch_on = model[:switch_on]
-    for t in 1:model[:instance].time
-        if !g.must_run[t]
-            # Link binary variables
-            if t == 1
-                eq_binary_link[g.name, t] = @constraint(
-                    model,
-                    is_on[g.name, t] - _is_initially_on(g) ==
-                    switch_on[g.name, t] - switch_off[g.name, t]
-                )
-            else
-                eq_binary_link[g.name, t] = @constraint(
-                    model,
-                    is_on[g.name, t] - is_on[g.name, t-1] ==
-                    switch_on[g.name, t] - switch_off[g.name, t]
-                )
-            end
-            # Cannot switch on and off at the same time
-            eq_switch_on_off[g.name, t] = @constraint(
-                model,
-                switch_on[g.name, t] + switch_off[g.name, t] <= 1
-            )
-        end
-    end
-    return
-end
-
 function _add_ramp_eqs!(
     model::JuMP.Model,
     g::Unit,
@@ -194,7 +117,7 @@ function _add_ramp_eqs!(
     eq_ramp_up = _init(model, :eq_ramp_up)
     eq_ramp_down = _init(model, :eq_ramp_down)
     for t in 1:model[:instance].time
-        # Ramp up limit
+        # Ramp up limit 
         if t == 1
             if _is_initially_on(g) == 1
                 eq_ramp_up[g.name, t] = @constraint(

src/transform/randomize.jl

@@ -0,0 +1,53 @@
+# 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/utils/sysimage.jl

@@ -5,12 +5,20 @@
 using PackageCompiler
 
 using DataStructures
+using Distributions
 using JSON
 using JuMP
 using MathOptInterface
 using SparseArrays
 
-pkg = [:DataStructures, :JSON, :JuMP, :MathOptInterface, :SparseArrays]
+pkg = [
+    :DataStructures,
+    :Distributions,
+    :JSON,
+    :JuMP,
+    :MathOptInterface,
+    :SparseArrays,
+]
 
 @info "Building system image..."
 create_sysimage(

src/validation/validate.jl

@@ -208,12 +208,8 @@ function _validate_units(instance, solution; tol = 0.01)
                         break
                     end
                 end
-                if t == time_down + 1
-                    initial_down = unit.min_downtime
-                    if unit.initial_status < 0
-                        initial_down = -unit.initial_status
-                    end
-                    time_down += initial_down
+                if (t == time_down + 1) && (unit.initial_status < 0)
+                    time_down -= unit.initial_status
                 end
 
                 # Calculate startup costs
@@ -246,14 +242,6 @@ function _validate_units(instance, solution; tol = 0.01)
                         break
                     end
                 end
-                if t == time_up + 1
-                    initial_up = unit.min_uptime
-                    if unit.initial_status > 0
-                        initial_up = unit.initial_status
-                    end
-                    time_up += initial_up
-                end
-
                 if (t == time_up + 1) && (unit.initial_status > 0)
                     time_up += unit.initial_status
                 end

test/model/formulations_test.jl

@@ -3,6 +3,7 @@
 # Released under the modified BSD license. See COPYING.md for more details.
 
 using UnitCommitment
+using JuMP
 import UnitCommitment:
     ArrCon2000,
     CarArr2006,
@@ -11,17 +12,49 @@ import UnitCommitment:
     Gar1962,
     KnuOstWat2018,
     MorLatRam2013,
-    PanGua2016
+    PanGua2016,
+    XavQiuWanThi2019
 
-function _test(formulation::Formulation)::Nothing
+if ENABLE_LARGE_TESTS
+    using Gurobi
+end
+
+function _small_test(formulation::Formulation)::Nothing
     instance = UnitCommitment.read_benchmark("matpower/case118/2017-02-01")
     UnitCommitment.build_model(instance = instance, formulation = formulation)  # should not crash
     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,
+        )
+        UnitCommitment.optimize!(
+            model,
+            XavQiuWanThi2019.Method(two_phase_gap = false, gap_limit = 0.1),
+        )
+        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(ramping = DamKucRajAta2016.Ramping()))
     _test(
         Formulation(
             ramping = MorLatRam2013.Ramping(),

test/runtests.jl

@@ -7,6 +7,8 @@ using UnitCommitment
 
 UnitCommitment._setup_logger()
 
+const ENABLE_LARGE_TESTS = ("UCJL_LARGE_TESTS" in keys(ENV))
+
 @testset "UnitCommitment" begin
     include("usage.jl")
     @testset "import" begin
@@ -26,6 +28,7 @@ UnitCommitment._setup_logger()
     @testset "transform" begin
         include("transform/initcond_test.jl")
         include("transform/slice_test.jl")
+        include("transform/randomize_test.jl")
     end
     @testset "validation" begin
         include("validation/repair_test.jl")

test/transform/randomize_test.jl

@@ -0,0 +1,43 @@
+# 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